ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Morrison.H
Go to the documentation of this file.
1 /*
2  * Implementation of a 2-moment microphysics model
3  * This model is based on the Morrison model as implemented in WRF (in Fortran)
4  */
5 #ifndef ERF_Morrison_H
6 #define ERF_Morrison_H
7 
8 #include <string>
9 #include <vector>
10 #include <memory>
11 
12 #include <AMReX_FArrayBox.H>
13 #include <AMReX_Geometry.H>
14 #include <AMReX_TableData.H>
15 #include <AMReX_MultiFabUtil.H>
16 
17 #include "ERF_Constants.H"
18 #include "ERF_MicrophysicsUtils.H"
19 #include "ERF_IndexDefines.H"
20 #include "ERF_DataStruct.H"
21 #include "ERF_NullMoist.H"
23 
24 namespace MicVar_Morr {
25  enum {
26  // independent variables
27  rho=0, // density
28  theta, // liquid/ice water potential temperature
29  tabs, // temperature
30  pres, // pressure
31  // non-precipitating vars
32  qt, // total cloud
33  qn, // cloud condesnate (liquid + ice)
34  qv, // cloud vapor
35  qcl, // cloud water
36  qci, // cloud ice
37  // precipitating vars
38  qp, // total precip
39  qpr, // precip rain
40  qps, // precip ice
41  qpg, // graupel
42  // The "numbers" and the "accums" are written/read at checkpoint/restart
43  // number concentrations
44  nc, // cloud droplet number
45  nr, // rain number
46  ni, // cloud ice number
47  ns, // snow number
48  ng, // graupel number
49  // derived vars
54  NumVars
55  };
56 }
57 
58 class Morrison : public NullMoist {
59 
60  using FabPtr = std::shared_ptr<amrex::MultiFab>;
61 
62 public:
63  // constructor
64  Morrison () {}
65 
66  // destructor
67  virtual ~Morrison () = default;
68 
69  // Set up for first time
70  void
71  Define (SolverChoice& sc) override
72  {
73  m_fac_cond = lcond / sc.c_p;
74  m_fac_fus = lfus / sc.c_p;
75  m_fac_sub = lsub / sc.c_p;
76  m_gOcp = CONST_GRAV / sc.c_p;
77  m_axis = sc.ave_plane;
78  m_rdOcp = sc.rdOcp;
79  m_do_cond = (!sc.use_shoc);
80  }
81 
82  // init
83  void
84  Init (const amrex::MultiFab& cons_in,
85  const amrex::BoxArray& grids,
86  const amrex::Geometry& geom,
87  const amrex::Real& dt_advance,
88  std::unique_ptr<amrex::MultiFab>& z_phys_nd,
89  std::unique_ptr<amrex::MultiFab>& detJ_cc) override;
90 
91  // Copy state into micro vars
92  void
93  Copy_State_to_Micro (const amrex::MultiFab& cons_in) override;
94 
95  // Copy state into micro vars
96  void
97  Copy_Micro_to_State (amrex::MultiFab& cons_in) override;
98 
99  void
100  Update_Micro_Vars (amrex::MultiFab& cons_in) override
101  {
102  this->Copy_State_to_Micro(cons_in);
103  this->Compute_Coefficients();
104  }
105 
106  void
107  Update_State_Vars (amrex::MultiFab& cons_in) override
108  {
109  this->Copy_Micro_to_State(cons_in);
110  }
111 
112  // wrapper to do all the updating
113  void
114  Advance (const amrex::Real& dt_advance,
115  const SolverChoice& sc) override;
116 
117  amrex::MultiFab*
118  Qmoist_Ptr (const int& varIdx) override
119  {
120  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
121  return mic_fab_vars[MicVarMap[varIdx]].get();
122  }
123 
124  void
126 
127  int
128  Qmoist_Size () override { return Morrison::m_qmoist_size; }
129 
130  int
132 
133  AMREX_GPU_HOST_DEVICE
134  AMREX_FORCE_INLINE
135  static amrex::Real
136  NewtonIterSat (int& i, int& j, int& k,
137  const amrex::Real& fac_cond,
138  const amrex::Real& fac_fus,
139  const amrex::Real& fac_sub,
140  const amrex::Real& an,
141  const amrex::Real& bn,
142  const amrex::Array4<amrex::Real>& tabs_array,
143  const amrex::Array4<amrex::Real>& pres_array,
144  const amrex::Array4<amrex::Real>& qv_array,
145  const amrex::Array4<amrex::Real>& qc_array,
146  const amrex::Array4<amrex::Real>& qi_array,
147  const amrex::Array4<amrex::Real>& qn_array,
148  const amrex::Array4<amrex::Real>& qt_array)
149  {
150  // Solution tolerance
151  amrex::Real tol = 1.0e-4;
152 
153  // Saturation moisture fractions
154  amrex::Real omn, domn;
155  amrex::Real qsat, dqsat;
156  amrex::Real qsatw, dqsatw;
157  amrex::Real qsati, dqsati;
158 
159  // Newton iteration vars
160  int niter;
161  amrex::Real fff, dfff, dtabs;
162  amrex::Real lstar, dlstar;
163  amrex::Real lstarw, lstari;
164  amrex::Real delta_qv, delta_qc, delta_qi;
165 
166  // Initial guess for temperature & pressure
167  amrex::Real tabs = tabs_array(i,j,k);
168  amrex::Real pres = pres_array(i,j,k);
169 
170  niter = 0;
171  dtabs = 1;
172  //==================================================
173  // Newton iteration to qv=qsat (cloud phase only)
174  //==================================================
175  do {
176  // Latent heats and their derivatives wrt to T
177  lstarw = fac_cond;
178  lstari = fac_fus;
179  domn = 0.0;
180 
181  // Saturation moisture fractions
182  erf_qsatw(tabs, pres, qsatw);
183  erf_qsati(tabs, pres, qsati);
184  erf_dtqsatw(tabs, pres, dqsatw);
185  erf_dtqsati(tabs, pres, dqsati);
186 
187  // Cloud ice not permitted (condensation & fusion)
188  if(tabs >= tbgmax) {
189  omn = 1.0;
190  }
191  // Cloud water not permitted (sublimation & fusion)
192  else if(tabs <= tbgmin) {
193  omn = 0.0;
194  lstarw = fac_sub;
195  }
196  // Mixed cloud phase (condensation & fusion)
197  else {
198  omn = an*tabs-bn;
199  domn = an;
200  }
201 
202  // Linear combination of each component
203  qsat = omn * qsatw + (1.0-omn) * qsati;
204  dqsat = omn * dqsatw + (1.0-omn) * dqsati
205  + domn * qsatw - domn * qsati;
206  lstar = omn * lstarw + (1.0-omn) * lstari;
207  dlstar = domn * lstarw - domn * lstari;
208 
209  // Function for root finding:
210  // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
211  fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat);
212 
213  // Derivative of function (T_new iterated on)
214  dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;
215 
216  // Update the temperature
217  dtabs = -fff/dfff;
218  tabs += dtabs;
219 
220  // Update iteration
221  niter = niter+1;
222  } while(std::abs(dtabs) > tol && niter < 20);
223 
224  // Update qsat from last iteration (dq = dq/dt * dt)
225  qsat += dqsat*dtabs;
226 
227  // Changes in each component
228  delta_qv = qv_array(i,j,k) - qsat;
229  delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn);
230  delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn));
231 
232  // Partition the change in non-precipitating q
233  qv_array(i,j,k) = qsat;
234  qc_array(i,j,k) += delta_qc;
235  qi_array(i,j,k) += delta_qi;
236  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
237  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
238 
239  // Return to temperature
240  return tabs;
241  }
242 
243  void
245  std::vector<int>& a_idx,
246  std::vector<std::string>& a_names) const override
247  {
248  a_idx.clear();
249  a_names.clear();
250 
251  // The ordering here needs to match that in
252  // MicVarMap = {MicVar_Morr::nc, MicVar_Morr::nr, MicVar_Morr::ni, MicVar_Morr::ns, MicVar_Morr::ng,
253  // MicVar_Morr::rain_accum, MicVar_Morr::snow_accum, MicVar_Morr::graup_accum};
254  //
255  a_idx.push_back(0); a_names.push_back("Nc");
256  a_idx.push_back(1); a_names.push_back("Nr");
257  a_idx.push_back(2); a_names.push_back("Ni");
258  a_idx.push_back(3); a_names.push_back("Ns");
259  a_idx.push_back(4); a_names.push_back("Ng");
260  a_idx.push_back(5); a_names.push_back("RainAccum");
261  a_idx.push_back(6); a_names.push_back("SnowAccum");
262  a_idx.push_back(7); a_names.push_back("GraupAccum");
263  }
264 
265 private:
266  // Number of qmoist variables (Ncloud, Nrain, Nice, Nsnow, Ngraup, rain_accum, snow_accum, graup_accum)
267  int m_qmoist_size = 8;
268 
269  // Number of qstate variables
271 
272  // CFL MAX for vertical advection
273  static constexpr amrex::Real CFL_MAX = 0.5;
274 
275  // MicVar map (Qmoist indices -> MicVar enum)
276  amrex::Vector<int> MicVarMap;
277 
278  // geometry
279  amrex::Geometry m_geom;
280 
281  // valid boxes on which to evolve the solution
282  amrex::BoxArray m_gtoe;
283 
284  // number of vertical levels
285  int nlev, zlo, zhi;
286 
287  // plane average axis
288  int m_axis;
289 
290  // constants
296  bool m_do_cond;
297 
298  // Pointer to terrain data
299  amrex::MultiFab* m_z_phys_nd;
300  amrex::MultiFab* m_detJ_cc;
301 
302  // independent variables
303  amrex::Array<FabPtr, MicVar_Morr::NumVars> mic_fab_vars;
304 
305  // microphysics parameters/coefficients
306  amrex::TableData<amrex::Real, 1> accrrc;
307  amrex::TableData<amrex::Real, 1> accrsi;
308  amrex::TableData<amrex::Real, 1> accrsc;
309  amrex::TableData<amrex::Real, 1> coefice;
310  amrex::TableData<amrex::Real, 1> evaps1;
311  amrex::TableData<amrex::Real, 1> evaps2;
312  amrex::TableData<amrex::Real, 1> accrgi;
313  amrex::TableData<amrex::Real, 1> accrgc;
314  amrex::TableData<amrex::Real, 1> evapg1;
315  amrex::TableData<amrex::Real, 1> evapg2;
316  amrex::TableData<amrex::Real, 1> evapr1;
317  amrex::TableData<amrex::Real, 1> evapr2;
318 
319  // vertical plane average data
320  amrex::TableData<amrex::Real, 1> rho1d;
321  amrex::TableData<amrex::Real, 1> pres1d;
322  amrex::TableData<amrex::Real, 1> tabs1d;
323  amrex::TableData<amrex::Real, 1> qt1d;
324  amrex::TableData<amrex::Real, 1> qv1d;
325  amrex::TableData<amrex::Real, 1> qn1d;
326 
327  amrex::TableData<amrex::Real, 1> gamaz;
328  amrex::TableData<amrex::Real, 1> zmid; // mid value of vertical coordinate in physical domain
329 };
330 #endif
constexpr amrex::Real lfus
Definition: ERF_Constants.H:67
constexpr amrex::Real tbgmax
Definition: ERF_Constants.H:32
constexpr amrex::Real lsub
Definition: ERF_Constants.H:68
constexpr amrex::Real lcond
Definition: ERF_Constants.H:66
constexpr amrex::Real tbgmin
Definition: ERF_Constants.H:31
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:178
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsati(amrex::Real t, amrex::Real p, amrex::Real &dtqsati)
Definition: ERF_MicrophysicsUtils.H:170
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:163
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsati(amrex::Real t, amrex::Real p, amrex::Real &qsati)
Definition: ERF_MicrophysicsUtils.H:152
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Definition: ERF_Morrison.H:58
virtual ~Morrison()=default
amrex::Geometry m_geom
Definition: ERF_Morrison.H:279
Morrison()
Definition: ERF_Morrison.H:64
void Advance(const amrex::Real &dt_advance, const SolverChoice &sc) override
Definition: ERF_AdvanceMorrison.cpp:503
amrex::TableData< amrex::Real, 1 > evapg2
Definition: ERF_Morrison.H:315
void Qmoist_Restart_Vars(const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
Definition: ERF_Morrison.H:244
amrex::TableData< amrex::Real, 1 > qn1d
Definition: ERF_Morrison.H:325
int Qstate_Moist_Size() override
Definition: ERF_Morrison.H:131
int n_qstate_moist_size
Definition: ERF_Morrison.H:270
void Update_Micro_Vars(amrex::MultiFab &cons_in) override
Definition: ERF_Morrison.H:100
amrex::TableData< amrex::Real, 1 > pres1d
Definition: ERF_Morrison.H:321
amrex::TableData< amrex::Real, 1 > zmid
Definition: ERF_Morrison.H:328
amrex::Real m_gOcp
Definition: ERF_Morrison.H:294
amrex::TableData< amrex::Real, 1 > accrgi
Definition: ERF_Morrison.H:312
static constexpr amrex::Real CFL_MAX
Definition: ERF_Morrison.H:273
amrex::TableData< amrex::Real, 1 > accrgc
Definition: ERF_Morrison.H:313
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitMorrison.cpp:98
amrex::TableData< amrex::Real, 1 > tabs1d
Definition: ERF_Morrison.H:322
amrex::Vector< int > MicVarMap
Definition: ERF_Morrison.H:276
amrex::MultiFab * Qmoist_Ptr(const int &varIdx) override
Definition: ERF_Morrison.H:118
amrex::MultiFab * m_z_phys_nd
Definition: ERF_Morrison.H:299
amrex::TableData< amrex::Real, 1 > qt1d
Definition: ERF_Morrison.H:323
int nlev
Definition: ERF_Morrison.H:285
void Compute_Coefficients()
Definition: ERF_InitMorrison.cpp:152
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateMorrison.cpp:17
void Define(SolverChoice &sc) override
Definition: ERF_Morrison.H:71
amrex::MultiFab * m_detJ_cc
Definition: ERF_Morrison.H:300
amrex::TableData< amrex::Real, 1 > gamaz
Definition: ERF_Morrison.H:327
amrex::TableData< amrex::Real, 1 > qv1d
Definition: ERF_Morrison.H:324
amrex::Real m_fac_cond
Definition: ERF_Morrison.H:291
amrex::Real m_fac_sub
Definition: ERF_Morrison.H:293
amrex::TableData< amrex::Real, 1 > coefice
Definition: ERF_Morrison.H:309
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat(int &i, int &j, int &k, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &fac_sub, const amrex::Real &an, const amrex::Real &bn, const amrex::Array4< amrex::Real > &tabs_array, const amrex::Array4< amrex::Real > &pres_array, const amrex::Array4< amrex::Real > &qv_array, const amrex::Array4< amrex::Real > &qc_array, const amrex::Array4< amrex::Real > &qi_array, const amrex::Array4< amrex::Real > &qn_array, const amrex::Array4< amrex::Real > &qt_array)
Definition: ERF_Morrison.H:136
int m_axis
Definition: ERF_Morrison.H:288
amrex::TableData< amrex::Real, 1 > accrsi
Definition: ERF_Morrison.H:307
void Init(const amrex::MultiFab &cons_in, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &z_phys_nd, std::unique_ptr< amrex::MultiFab > &detJ_cc) override
Definition: ERF_InitMorrison.cpp:22
int Qmoist_Size() override
Definition: ERF_Morrison.H:128
amrex::TableData< amrex::Real, 1 > accrrc
Definition: ERF_Morrison.H:306
amrex::TableData< amrex::Real, 1 > evapg1
Definition: ERF_Morrison.H:314
int m_qmoist_size
Definition: ERF_Morrison.H:267
amrex::TableData< amrex::Real, 1 > evapr1
Definition: ERF_Morrison.H:316
int zlo
Definition: ERF_Morrison.H:285
amrex::TableData< amrex::Real, 1 > evaps1
Definition: ERF_Morrison.H:310
amrex::BoxArray m_gtoe
Definition: ERF_Morrison.H:282
amrex::TableData< amrex::Real, 1 > evaps2
Definition: ERF_Morrison.H:311
amrex::Real m_rdOcp
Definition: ERF_Morrison.H:295
amrex::Array< FabPtr, MicVar_Morr::NumVars > mic_fab_vars
Definition: ERF_Morrison.H:303
amrex::TableData< amrex::Real, 1 > rho1d
Definition: ERF_Morrison.H:320
std::shared_ptr< amrex::MultiFab > FabPtr
Definition: ERF_Morrison.H:60
int zhi
Definition: ERF_Morrison.H:285
amrex::TableData< amrex::Real, 1 > evapr2
Definition: ERF_Morrison.H:317
void Update_State_Vars(amrex::MultiFab &cons_in) override
Definition: ERF_Morrison.H:107
amrex::TableData< amrex::Real, 1 > accrsc
Definition: ERF_Morrison.H:308
bool m_do_cond
Definition: ERF_Morrison.H:296
amrex::Real m_fac_fus
Definition: ERF_Morrison.H:292
Definition: ERF_NullMoist.H:8
Definition: ERF_Morrison.H:24
@ qv
Definition: ERF_Morrison.H:34
@ ng
Definition: ERF_Morrison.H:48
@ nc
Definition: ERF_Morrison.H:44
@ qpg
Definition: ERF_Morrison.H:41
@ pres
Definition: ERF_Morrison.H:30
@ nr
Definition: ERF_Morrison.H:45
@ qcl
Definition: ERF_Morrison.H:35
@ tabs
Definition: ERF_Morrison.H:29
@ theta
Definition: ERF_Morrison.H:28
@ qp
Definition: ERF_Morrison.H:38
@ ni
Definition: ERF_Morrison.H:46
@ ns
Definition: ERF_Morrison.H:47
@ omega
Definition: ERF_Morrison.H:53
@ qps
Definition: ERF_Morrison.H:40
@ qn
Definition: ERF_Morrison.H:33
@ graup_accum
Definition: ERF_Morrison.H:52
@ rho
Definition: ERF_Morrison.H:27
@ qpr
Definition: ERF_Morrison.H:39
@ qci
Definition: ERF_Morrison.H:36
@ NumVars
Definition: ERF_Morrison.H:54
@ rain_accum
Definition: ERF_Morrison.H:50
@ snow_accum
Definition: ERF_Morrison.H:51
@ qt
Definition: ERF_Morrison.H:32
Definition: ERF_DataStruct.H:123
amrex::Real rdOcp
Definition: ERF_DataStruct.H:866
amrex::Real c_p
Definition: ERF_DataStruct.H:865
bool use_shoc
Definition: ERF_DataStruct.H:895
int ave_plane
Definition: ERF_DataStruct.H:926