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