tMask.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tMask.hh (4782B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Constantine Khroulev and David Maxwell
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #ifndef _MASK_H_
20 #define _MASK_H_
21
22 // the following three includes are needed here because of inlined code
23 #include "iceModelVec.hh"
24 #include "ConfigInterface.hh"
25 #include "error_handling.hh"
26
27 namespace pism {
28
29 enum MaskValue {
30 MASK_UNKNOWN = -1,
31 MASK_ICE_FREE_BEDROCK = 0,
32 MASK_GROUNDED = 2,
33 MASK_FLOATING = 3,
34 MASK_ICE_FREE_OCEAN = 4
35 };
36
37 namespace mask {
38 //! \brief An ocean cell (floating ice or ice-free).
39 inline bool ocean(int M) {
40 return M >= MASK_FLOATING;
41 }
42 //! \brief Grounded cell (grounded ice or ice-free).
43 inline bool grounded(int M) {
44 return not ocean(M);
45 }
46 //! \brief Ice-filled cell (grounded or floating).
47 inline bool icy(int M) {
48 return (M == MASK_GROUNDED) || (M == MASK_FLOATING);
49 }
50 inline bool grounded_ice(int M) {
51 return icy(M) && grounded(M);
52 }
53 inline bool floating_ice(int M) {
54 return icy(M) && ocean(M);
55 }
56 //! \brief Ice-free cell (grounded or ocean).
57 inline bool ice_free(int M) {
58 return not icy(M);
59 }
60 inline bool ice_free_ocean(int M) {
61 return ocean(M) && ice_free(M);
62 }
63 inline bool ice_free_land(int M) {
64 return grounded(M) && ice_free(M);
65 }
66 }
67
68 class GeometryCalculator {
69 public:
70 GeometryCalculator(const Config &config) {
71 m_alpha = 1 - config.get_number("constants.ice.density") / config.get_number("constants.sea_water.density");
72 m_is_dry_simulation = config.get_flag("ocean.always_grounded");
73 m_icefree_thickness = config.get_number("geometry.ice_free_thickness_standard");
74 if (m_icefree_thickness < 0.0) {
75 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
76 "invalid ice-free thickness threshold: %f", m_icefree_thickness);
77 }
78 }
79
80 void set_icefree_thickness(double threshold) {
81 if (threshold < 0.0) {
82 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid ice-free thickness threshold: %f", threshold);
83 }
84 m_icefree_thickness = threshold;
85 }
86
87 void compute(const IceModelVec2S &sea_level, const IceModelVec2S &bed, const IceModelVec2S &thickness,
88 IceModelVec2Int &out_mask, IceModelVec2S &out_surface) const;
89
90 void compute_mask(const IceModelVec2S& sea_level, const IceModelVec2S& bed,
91 const IceModelVec2S& thickness, IceModelVec2Int& result) const;
92
93 void compute_surface(const IceModelVec2S& sea_level, const IceModelVec2S& bed,
94 const IceModelVec2S& thickness, IceModelVec2S& result) const;
95
96 inline void compute(double sea_level, double bed, double thickness,
97 int *out_mask, double *out_surface) const {
98 const double hgrounded = bed + thickness; // FIXME issue #15
99 const double hfloating = sea_level + m_alpha*thickness;
100
101 const bool
102 is_floating = (hfloating > hgrounded),
103 ice_free = (thickness <= m_icefree_thickness);
104
105 int mask_result;
106 double surface_result;
107
108 if (is_floating && (not m_is_dry_simulation)) {
109 surface_result = hfloating;
110
111 if (ice_free) {
112 mask_result = MASK_ICE_FREE_OCEAN;
113 } else {
114 mask_result = MASK_FLOATING;
115 }
116 } else { // Grounded
117 surface_result = hgrounded;
118
119 if (ice_free) {
120 mask_result = MASK_ICE_FREE_BEDROCK;
121 } else {
122 mask_result = MASK_GROUNDED;
123 }
124 }
125
126 if (out_surface != NULL) {
127 *out_surface = surface_result;
128 }
129
130 if (out_mask != NULL) {
131 *out_mask = mask_result;
132 }
133 }
134
135 inline int mask(double sea_level, double bed, double thickness) const {
136 int result;
137 compute(sea_level, bed, thickness, &result, NULL);
138 return result;
139 }
140
141 inline double surface(double sea_level, double bed, double thickness) const {
142 double result;
143 compute(sea_level, bed, thickness, NULL, &result);
144 return result;
145 }
146
147 protected:
148 double m_alpha;
149 double m_icefree_thickness;
150 bool m_is_dry_simulation;
151 };
152
153 } // end of namespace pism
154
155 #endif /* _MASK_H_ */