tMassEnergyBudget.cc - 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
---
tMassEnergyBudget.cc (8259B)
---
1 #include <iostream>
2 #include <pism/icebin/MassEnergyBudget.hh>
3 #include <pism/util/error_handling.hh>
4
5 namespace pism{
6 namespace icebin{
7
8
9 void MassEnthVec2S::create(pism::IceGrid::ConstPtr my_grid, const std::string &my_name,
10 pism::IceModelVecKind ghostedp, int width)
11 {
12 mass.create(my_grid, my_name + ".mass", ghostedp, width);
13 enth.create(my_grid, my_name + ".enth", ghostedp, width);
14 }
15
16 void MassEnthVec2S::set_attrs(
17 const std::string &my_pism_intent,
18 const std::string &my_long_name,
19 const std::string &my_units,
20 const std::string &my_standard_name)
21 {
22 auto mass_units = "kg " + my_units;
23 auto energy_units = "J " + my_units;
24 mass.set_attrs(my_pism_intent, my_long_name + " (mass portion)",
25 mass_units, mass_units, my_standard_name + ".mass", 0);
26 enth.set_attrs(my_pism_intent, my_long_name + " (enthalpy portion)",
27 energy_units, energy_units, my_standard_name + ".enth", 0);
28 }
29
30 MassEnergyBudget::MassEnergyBudget()
31 {
32 }
33
34 void MassEnergyBudget::create(pism::IceGrid::ConstPtr grid, std::string const &prefix,
35 pism::IceModelVecKind ghostedp, unsigned int width)
36 {
37 if (all_vecs.size() != 0) throw RuntimeError(PISM_ERROR_LOCATION,
38 "MassEnergyBudget::create() cannot be called twice, fix your code!");
39 printf("MassEnergyBudget(%p)::create()\n", (void*)this);
40
41 // ----------- Mass and Enthalpy State of the Ice Sheet
42 total.create(grid, prefix+"total",
43 ghostedp, width);
44 total.set_attrs("diagnostic",
45 "State of the ice sheet (NOT a difference between timetseps)",
46 "m-2", "total");
47 add_massenth(total, TOTAL, "", "");
48
49 // ----------- Heat generation of flows [vertical]
50 // Postive means heat is flowing INTO the ice sheet.
51 basal_frictional_heating.create(grid, prefix+"basal_frictional_heating",
52 ghostedp, width);
53 basal_frictional_heating.set_attrs("internal",
54 "Basal frictional heating",
55 "W m-2", "W m-2", "", 0);
56 add_enth(basal_frictional_heating, DELTA, "basal_frictional_heating");
57
58 strain_heating.create(grid, prefix+"strain_heating",
59 ghostedp, width);
60 strain_heating.set_attrs("internal",
61 "Strain heating",
62 "W m-2", "W m-2", "", 0);
63 add_enth(strain_heating, DELTA, "strain_heating"); //!< Total amount of strain heating [J/m^2]
64
65 geothermal_flux.create(grid, prefix+"geothermal_flux",
66 ghostedp, width);
67 geothermal_flux.set_attrs("internal",
68 "Geothermal energy through (compare to upward_geothermal_flux?)",
69 "W m-2", "W m-2", "", 0);
70 add_enth(geothermal_flux, 0, "geothermal_flux"); //!< Total amount of geothermal energy [J/m^2]
71
72 upward_geothermal_flux.create(grid, prefix+"upward_geothermal_flux",
73 ghostedp, width);
74 upward_geothermal_flux.set_attrs("internal",
75 "Geothermal energy through (compare to geothermal_flux?)",
76 "W m-2", "W m-2", "", 0);
77 add_enth(upward_geothermal_flux, DELTA, "upward_geothermal_flux"); //!< Total amount of geothermal energy [J/m^2]
78
79 // ----------- Mass advection, with accompanying enthalpy change
80 // Postive means mass/enthalpy is flowing INTO the ice sheet.
81 std::string name;
82
83 calving.create(grid, prefix+"calving",
84 ghostedp, width);
85 calving.set_attrs("diagnostic",
86 "Mass/Enthalpy gain from calving. Should be negative.",
87 "m-2 s-1", "calving");
88 add_massenth(calving, DELTA, "calving.mass", "calving.enth");
89
90 pism_smb.create(grid, prefix+"pism_smb",
91 ghostedp, width);
92 pism_smb.set_attrs("diagnostic",
93 "pism_smb",
94 "m-2 s-1", "pism_smb");
95 // No DELTA< does not participate in epsilon computation
96 add_massenth(pism_smb, 0, "pism_smb.mass", "pism_smb.enth");
97
98 icebin_xfer.create(grid, prefix+"icebin_xfer",
99 ghostedp, width);
100 icebin_xfer.set_attrs("diagnostic",
101 "icebin_xfer",
102 "m-2 s-1", "icebin_xfer");
103 add_massenth(icebin_xfer, DELTA, "icebin_xfer.mass", "icebin_xfer.enth");
104
105 icebin_deltah.create(grid, prefix+"icebin_deltah",
106 ghostedp, width);
107 icebin_deltah.set_attrs("diagnostic",
108 "icebin_deltah",
109 "J m-2 s-1", "J m-2 s-1", "icebin_deltah", 0);
110 add_enth(icebin_deltah, DELTA, "");
111
112 href_to_h.create(grid, prefix+"href_to_h",
113 ghostedp, width);
114 href_to_h.set_attrs("diagnostic",
115 "href_to_h",
116 "kg m-2 s-1", "kg m-2 s-1", "href_to_h", 0);
117 add_mass(href_to_h, 0, "");
118
119 nonneg_rule.create(grid, prefix+"nonneg_rule",
120 ghostedp, width);
121 nonneg_rule.set_attrs("diagnostic",
122 "nonneg_rule",
123 "kg m-2 s-1", "kg m-2 s-1", "nonneg_rule", 0);
124 add_mass(nonneg_rule, 0, "");
125
126
127
128 melt_grounded.create(grid, prefix+"melt_grounded",
129 ghostedp, width);
130 melt_grounded.set_attrs("diagnostic",
131 "Basal melting of grounded ice (negative)",
132 "m-2 s-1", "melt_grounded");
133 add_massenth(melt_grounded, DELTA, "melt_grounded.mass", "melt_grounded.enth");
134
135 melt_floating.create(grid, prefix+"melt_floating",
136 ghostedp, width);
137 melt_floating.set_attrs("diagnostic",
138 "Sub-shelf melting (negative)",
139 "m-2 s-1", "melt_floating");
140 add_massenth(melt_floating, DELTA, "melt_floating.mass", "melt_floating.enth");
141
142 // ----------- Advection WITHIN the ice sheet
143 internal_advection.create(grid, prefix+"internal_advection",
144 ghostedp, width);
145 internal_advection.set_attrs("diagnostic",
146 "Advection within the ice sheet",
147 "m-2 s-1", "internal_advection");
148 add_massenth(internal_advection, DELTA, "internal_advection.mass", "internal_advection.enth");
149
150 // ----------- Balance the Budget
151 epsilon.create(grid, prefix+"epsilon",
152 ghostedp, width);
153 epsilon.set_attrs("diagnostic",
154 "Unaccounted-for changes",
155 "m-2 s-1", "epsilon");
156 add_massenth(epsilon, EPSILON, "epsilon.mass", "epsilon.enth");
157 }
158
159 std::ostream &MassEnergyBudget::print_formulas(std::ostream &out)
160 {
161 // MASS
162 out << "epsilon.mass = total.mass -" << std::endl;
163 out << " (";
164 for (auto ii=all_vecs.begin(); ii != all_vecs.end(); ++ii) {
165 if ((ii->flags & (DELTA | MASS)) != (DELTA | MASS)) continue;
166 char str[20];
167 sprintf(str, "%p", (void*)&ii->vec);
168 out << ii->vec.get_name() << " + ";
169 }
170 out << ")" << std::endl;
171
172 // Energy
173 out << "epsilon.enth = total.enth -" << std::endl;
174 out << " (";
175 for (auto ii=all_vecs.begin(); ii != all_vecs.end(); ++ii) {
176 if ((ii->flags & (DELTA | ENTH)) != (DELTA | ENTH)) continue;
177 char str[20];
178 sprintf(str, "%p", (void*)&ii->vec);
179 out << ii->vec.get_name() << " + ";
180 }
181 out << ")" << std::endl;
182
183 return out;
184 }
185
186
187 void MassEnergyBudget::set_epsilon(pism::IceGrid::ConstPtr grid)
188 {
189 // ==> epsilon = (sum of fluxes) - total
190 printf("BEGIN MassEnergyBudget::set_epsilon()\n");
191
192 // -------- Mass
193 epsilon.mass.begin_access();
194 total.mass.begin_access();
195 for (int i = grid->xs(); i < grid->xs() + grid->xm(); ++i) {
196 for (int j = grid->ys(); j < grid->ys() + grid->ym(); ++j) {
197 epsilon.mass(i,j) = total.mass(i,j);
198 }}
199 total.mass.end_access();
200
201 for (auto &ii : all_vecs) {
202 if ((ii.flags & (DELTA | MASS)) != (DELTA | MASS)) continue;
203
204 printf("epsilon.mass: %s\n", ii.vec.get_name().c_str());
205
206 ii.vec.begin_access();
207 for (int i = grid->xs(); i < grid->xs() + grid->xm(); ++i) {
208 for (int j = grid->ys(); j < grid->ys() + grid->ym(); ++j) {
209 epsilon.mass(i,j) -= ii.vec(i,j);
210 }}
211 ii.vec.end_access();
212 }
213 epsilon.mass.end_access();
214
215 // -------- Energy
216 epsilon.enth.begin_access();
217 total.enth.begin_access();
218 for (int i = grid->xs(); i < grid->xs() + grid->xm(); ++i) {
219 for (int j = grid->ys(); j < grid->ys() + grid->ym(); ++j) {
220 epsilon.enth(i,j) = total.enth(i,j);
221 }}
222 total.enth.end_access();
223
224 #if 1
225 for (auto &ii : all_vecs) {
226 if ((ii.flags & (DELTA | ENTH)) != (DELTA | ENTH)) continue;
227
228 printf("epsilon.energy: %s\n", ii.vec.get_name().c_str());
229
230 ii.vec.begin_access();
231 for (int i = grid->xs(); i < grid->xs() + grid->xm(); ++i) {
232 for (int j = grid->ys(); j < grid->ys() + grid->ym(); ++j) {
233 epsilon.enth(i,j) -= ii.vec(i,j);
234 }}
235 ii.vec.end_access();
236 }
237 #endif
238 epsilon.enth.end_access();
239
240 printf("END MassEnergyBudget::set_epsilon()\n");
241 }
242
243 }}