tMassEnergyBudget.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
---
tMassEnergyBudget.hh (6559B)
---
1 #pragma once
2
3 // --------------------------------
4 // PISM Includes... want to be included first
5 #include <petsc.h>
6 #include <pism/util/IceGrid.hh>
7 #include <pism/util/iceModelVec.hh>
8 // --------------------------------
9
10 namespace pism {
11 namespace icebin {
12
13 /** Encapsulates mass and enthalpy together. Used to tabulate total
14 enthalpy of a bunch of advected H2O based on its mass and specific
15 enthalpy. This allows us to have only one C++ variable per advected
16 quanity, instead of two. */
17 struct MassEnthVec2S : public pism::PetscAccessible
18 {
19 pism::IceModelVec2S mass;
20 pism::IceModelVec2S enth;
21
22 ~MassEnthVec2S() {}
23
24 void create(pism::IceGrid::ConstPtr my_grid, const std::string &my_name,
25 pism::IceModelVecKind ghostedp, int width = 1);
26
27
28 void set_attrs(
29 const std::string &my_pism_intent,
30 const std::string &my_long_name,
31 const std::string &my_units,
32 const std::string &my_standard_name);
33
34 virtual void begin_access() const
35 {
36 mass.begin_access();
37 enth.begin_access();
38 }
39
40 virtual void end_access() const
41 {
42 mass.end_access();
43 enth.end_access();
44 }
45
46 #if 0
47 /** @param mass [kg m-2]
48 @param specific_enthalpy [J kg-1] */
49 void add_mass(int i, int j, double mass, double specific_enthalpy) {
50 mass(i,j) += mass;
51 enthalpy(i,j) += mass * specific_enthalpy;
52 }
53
54 void clear() {
55 mass.set(0);
56 enth.set(0);
57 }
58 #endif
59 };
60
61 struct VecWithFlags {
62 pism::IceModelVec2S &vec;
63 int flags;
64
65 /** IF this variable is used directly as a contractual ice model
66 output, the name of that contract entry. */
67 std::string contract_name;
68
69 VecWithFlags(pism::IceModelVec2S &_vec, int _flags, std::string const &_contract_name) :
70 vec(_vec), flags(_flags), contract_name(_contract_name) {}
71 };
72
73 class MassEnergyBudget {
74 public:
75 // ============================================================
76 // Total State
77
78 // ------------ Enthalpy State
79 MassEnthVec2S total; // Total mass [kg m-2] and enthalpy [J m-2] of the ice sheet
80
81 // =============================================================
82 // Cumulative Fluxes
83
84 // ======================= Variables to accumulate PISM output
85 // These are accumulated as [kg m-2] or [J m-2]
86 // They are accumulated for the life of the simulation, and never zeroed.
87 // Other instances of MassEnergyBudget are used to compute differences.
88
89 // ----------- Heat generation of flows [vertical]
90 pism::IceModelVec2S basal_frictional_heating; //!< Total amount of basal friction heating [J/m^2]
91 pism::IceModelVec2S strain_heating; //!< Total amount of strain heating [J/m^2]
92 pism::IceModelVec2S geothermal_flux; //!< Total amount of geothermal energy [J/m^2]
93 pism::IceModelVec2S upward_geothermal_flux; //!< Total amount of geothermal energy [J/m^2]
94
95 // ----------- Mass advection, with accompanying enthalpy change
96 // The enthalpy reported for these mass fluxes are the enthalpies
97 // AS REPORTED TO ICEBIN! That is not necessarily the same as the enthalpy
98 // that PISM sees internally.
99 MassEnthVec2S calving; //!< Equal to IceModel::discharge_flux_2D_cumulative
100
101 // Let's not use basal_runoff (which would be derived from the variables
102 // in NullTransportHydrology).
103 //
104 // 1. It is derived directly from
105 // bmelt/basal_meltrate/melt_grounded/meltrate_grounded, which is
106 // already computed here.
107 //
108 // 2. bmelt plays directly into removal of mass from the ice sheet (see
109 // accumulateFluxes_massContExplicitStep() in iMgeometry.cc).
110 // Including basal_runoff would be double-counting it.
111 // MassEnthVec2S basal_runoff; //!< Enthalpy here is predictable, since runoff is 0C 100% water fraction.
112
113 MassEnthVec2S icebin_xfer; //!< accumulation / ablation, as provided by Icebin
114 pism::IceModelVec2S icebin_deltah; //!< Change in enthalpy of top layer
115 MassEnthVec2S pism_smb; //! SMB as seen by PISM in iMgeometry.cc massContExplicitSte(). Used to check icebin_xfer.mass, but does not figure into contract.
116 pism::IceModelVec2S href_to_h;
117 pism::IceModelVec2S nonneg_rule;
118 MassEnthVec2S melt_grounded; //!< basal melt (grounded) (from summing meltrate_grounded)
119 MassEnthVec2S melt_floating; //!< sub-shelf melt (from summing meltrate_floating)
120
121 // ----------- Mass advection WITHIN the ice sheet
122 MassEnthVec2S internal_advection;
123 // MassEnthVec2S divQ_SIA;
124 // MassEnthVec2S divQ_SSA; //!< flux divergence
125
126 // ======================= Balance the Budget
127 // At each step, we set epsilon as follows:
128 // total - (sum of fluxes) + epsilon = 0
129 // ==> epsilon = (sum of fluxes) - total
130 MassEnthVec2S epsilon;
131
132 // ======================== Different sets (flags)
133 static const int MASS = 1;
134 static const int ENTH = 2;
135
136 static const int TOTAL = 4; // To be differenced at the end.
137 static const int DELTA = 8;
138 static const int EPSILON = 16; // To be differenced at the end.
139
140 // ======================== Summary of above variables
141 // This makes it easy to difference two MassEnergyBudget instances.
142 std::vector<VecWithFlags> all_vecs;
143
144
145 // =====================================================================
146 std::ostream &print_formulas(std::ostream &out);
147
148 protected:
149 void add_mass(pism::IceModelVec2S &vec, int flags,
150 std::string const &contract_name)
151 {
152 all_vecs.push_back(VecWithFlags(vec, MASS | flags, contract_name));
153 }
154
155 /** @param contract_name The name of this variable in the ice model's output contract. */
156 void add_enth(pism::IceModelVec2S &vec, int flags,
157 std::string const &contract_name)
158 {
159 all_vecs.push_back(VecWithFlags(vec, ENTH | flags, contract_name));
160 }
161
162 void add_massenth(MassEnthVec2S &massenth, int flags,
163 std::string const &contract_name_mass,
164 std::string const &contract_name_enth)
165 {
166 all_vecs.push_back(VecWithFlags(massenth.mass, MASS | flags,
167 contract_name_mass));
168 all_vecs.push_back(VecWithFlags(massenth.enth, ENTH | flags,
169 contract_name_enth));
170 }
171
172 public:
173
174
175 MassEnergyBudget();
176
177 void create(pism::IceGrid::ConstPtr grid, std::string const &prefix,
178 pism::IceModelVecKind ghostedp, unsigned int width = 1);
179
180 void set_epsilon(pism::IceGrid::ConstPtr grid);
181 };
182
183 } // end of namespace icebin
184 } // end of namespace pism