tIBIceModel.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
---
tIBIceModel.hh (5884B)
---
1 #pragma once
2
3 // --------------------------------
4 // PISM Includes... want to be included first
5 #include <petsc.h>
6
7 #include <pism/icemodel/IceModel.hh>
8 #include <pism/util/IceGrid.hh>
9
10 #include <pism/util/pism_options.hh>
11 #include <pism/coupler/atmosphere/Factory.hh>
12 #include <pism/coupler/ocean/Factory.hh>
13 #include <pism/coupler/surface/Factory.hh>
14 #include <pism/hydrology/NullTransport.hh>
15
16 #include <pism/util/Time.hh>
17 // --------------------------------
18 #include <pism/icebin/IBSurfaceModel.hh>
19 #include <pism/icebin/MassEnergyBudget.hh>
20
21 // Stuff defined in the icebin library
22 // (NOT a dependency of ours)
23 namespace icebin {
24 namespace gpism {
25 class IceModel_PISM;
26 }
27 }
28
29 namespace pism {
30 namespace icebin {
31
32
33 /** This is the IceBin customized version of PISM's pism::IceModel class.
34
35 See https://github.com/pism/pism/issues/219
36
37 Here's a short term solution, though: create a new class
38 IBIceModel derived from IceModel and re-implement
39 IceModel::allocate_couplers(). In it, set
40 IceModel::external_surface_model and IceModel::external_ocean_model as
41 you see fit (IceModel will not de-allocate a surface (ocean) model if
42 it is set to true) and allocate PSConstantICEBIN. You might also want
43 to add IBIceModel::get_surface_model() which returns
44 IceModel::surface to get access to PSConstantICEBIN from outside of
45 IBIceModel.
46 */
47
48 class IBIceModel : public pism::IceModel {
49 friend class ::icebin::gpism::IceModel_PISM;
50
51 public:
52 typedef pism::IceModel super;
53 struct Params {
54 double time_start_s;
55 std::string output_dir;
56 };
57 Params const params;
58
59 protected:
60 MassEnergyBudget base; // Cumulative totals at start of this timestep
61 MassEnergyBudget cur; // Cumulative totals now
62 MassEnergyBudget rate; // At end of coupling timestep, set to (cur - base) / dt
63
64 // Output variables prepared for return to GCM
65 // (relevant ice model state to be exported)
66 pism::IceModelVec2S M1, M2;
67 pism::IceModelVec2S H1, H2;
68 pism::IceModelVec2S V1, V2;
69
70 protected:
71 // see iceModel.cc
72 virtual void allocate_storage();
73
74 public:
75 virtual void massContExplicitStep(double dt,
76 const IceModelVec2Stag &diffusive_flux,
77 const IceModelVec2V &advective_velocity);
78 virtual void accumulateFluxes_massContExplicitStep(int i, int j,
79 double surface_mass_balance, // [m s-1] ice equivalent (from PISM)
80 double meltrate_grounded, // [m s-1] ice equivalent
81 double meltrate_floating, // [m s-1] ice equivalent
82 double divQ_SIA, // [m s-1] ice equivalent
83 double divQ_SSA, // [m s-1] ice equivalent
84 double Href_to_H_flux, // [m s-1] ice equivalent
85 double nonneg_rule_flux); // [m s-1] ice equivalent
86 private:
87 // Temporary variables inside massContExplicitStep()
88 double _ice_density; // From config
89 double _meter_per_s_to_kg_per_m2; // Conversion factor computed from _ice_density
90
91
92 private:
93 // Utility function
94 void prepare_nc(std::string const &fname, std::unique_ptr<pism::File> &nc);
95
96 public:
97 /** @param t0 Time of last time we coupled. */
98 void set_rate(double dt);
99
100 void reset_rate();
101
102
103 std::unique_ptr<pism::File> pre_mass_nc; //!< Write variables every time massContPostHook() is called.
104 std::unique_ptr<pism::File> post_mass_nc;
105 std::unique_ptr<pism::File> pre_energy_nc;
106 std::unique_ptr<pism::File> post_energy_nc;
107
108 // see iceModel.cc for implementation of constructor and destructor:
109 /** @param gcm_params Pointer to IceModel::gcm_params. Lives at least as long as this object. */
110 IBIceModel(IceGrid::Ptr g, Context::Ptr context, IBIceModel::Params const &_params);
111 virtual ~IBIceModel(); // must be virtual merely because some members are virtual
112
113 virtual void allocate_subglacial_hydrology();
114 virtual void allocate_couplers();
115 virtual void time_setup();
116 virtual void misc_setup();
117
118 void compute_enth2(pism::IceModelVec2S &enth2, pism::IceModelVec2S &mass2);
119
120 /** @return Our instance of IBSurfaceModel */
121 pism::icebin::IBSurfaceModel *ib_surface_model() {
122 return dynamic_cast<IBSurfaceModel *>(m_surface.get());
123 }
124 pism::hydrology::NullTransport* null_hydrology() {
125 return dynamic_cast<hydrology::NullTransport *>(pism::IceModel::m_subglacial_hydrology.get());
126 }
127
128
129 /** @return Current time for mass timestepping */
130 double mass_t() const {
131 return m_time->current();
132 }
133 /** @return Current time for enthalpy timestepping */
134 double enthalpy_t() const {
135 return t_TempAge;
136 }
137
138 void energy_step();
139
140 void prepare_outputs(double time_s);
141
142 /** Read things out of the ice model that will be sent back BEFORE
143 the first coupling timestep (eg, ice surface enthalpy) */
144 void prepare_initial_outputs();
145
146 /** Merges surface temperature derived from the energy balance model into any NaN values
147 in the vector provided.
148 @param deltah IN: Input from Icebin (change in enthalpy of each grid
149 cell over the timestep) [W m-2].
150 @param default_val: The value that deltah(i,j) will have if no value
151 is listed for that grid cell
152 @param timestep_s: Length of the current coupling timestep [s]
153 @param surface_temp OUT: Resulting surface temperature to use as the Dirichlet B.C.
154 */
155 void construct_surface_temp(pism::IceModelVec2S &deltah, // IN: Input from Icebin
156 double default_val,
157 double timestep_s, // Length of this coupling interval [s]
158 pism::IceModelVec2S &surface_temp);
159 };
160 }
161 }