ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SuperDropletPCMassChange.H
Go to the documentation of this file.
1 #ifndef SUPERDROPLET_PC_MASSCHANGE_H_
2 #define SUPERDROPLET_PC_MASSCHANGE_H_
3 
4 #include <cmath>
5 #include <limits>
6 #include <AMReX_GpuPrint.H>
7 #include <AMReX_Math.H>
8 #include "ERF_NumericalConstants.H"
9 #include "ERF_MicrophysicsConstants.H"
10 
11 #ifdef ERF_USE_PARTICLES
12 
13 /*! \brief Namespace with classes and functions for condensation/evaporation
14  * (liquid <--> vapour) */
15 namespace SDMassChangeUtils_LV
16 {
17  /*! \brief Ventilation factor as a function of droplet radius (Bayley et al., 2025, Eq. 3) */
18  template <typename RT /*!< real-type */ >
19  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
20  RT ventilationFactor ( const RT a_R, /*!< radius */
21  const RT a_alpha1, /*!< fit coefficient alpha_1 */
22  const RT a_beta1, /*!< fit exponent beta_1 */
23  const RT a_alpha2, /*!< fit coefficient alpha_2 */
24  const RT a_beta2, /*!< fit exponent beta_2 */
25  const RT a_fcap /*!< maximum ventilation factor */ ) noexcept
26  {
27  RT x1 = a_alpha1 * std::pow(a_R, a_beta1);
28  RT x2 = a_alpha2 * std::pow(a_R, a_beta2);
29  RT fv = RT(1.0) + (x1*x2) / (x1 + x2);
30  return std::min(fv, a_fcap);
31  }
32 
33 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
34  /*! \brief Phase change equation (in terms of R) */
35  template <typename RT /*!< real-type */ >
36  struct dRdt
37  {
38  RT L; /*!< latent heat of vaporization (condensate) */
39  RT K; /*!< thermal conductivity (condensate) */
40  RT Rv; /*!< gas constant of air with vapour */
41  RT rho_l; /*!< density of condensate */
42 
43  /*! \brief Right-hand-side of the phase change ODE */
44  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
45  RT operator () ( const RT a_R, /*!< radius */
46  const RT a_S, /*!< saturation ratio */
47  const RT a_T, /*!< temperature */
48  const RT a_e_s, /*!< saturation pressure */
49  const RT a_D, /*!< mol. diffusion coeff */
50  const RT a_a, /*!< curvature coeff */
51  const RT a_b, /*!< solute coeff */
52  const RT a_N_s, /*!< total solute moles */
53  const RT a_f_v /*!< ventilation factor */ ) const noexcept
54  {
55  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
56  RT Kn = lambda_v/a_R;
57  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
58 
59  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
60  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
61 
62  RT R_inv = RT(1.0)/a_R;
63  RT R_inv_cubed = R_inv*R_inv*R_inv;
64 
65  RT alpha = (a_S-RT(1.0)) / (F_k + F_d);
66  RT retval = alpha*R_inv;
67 
68  RT beta = -(a_a/a_T) / (F_k + F_d);
69  retval += beta*R_inv*R_inv;
70 
71  RT gamma = a_b * a_N_s / (F_k + F_d);
72  retval += gamma*R_inv_cubed*R_inv;
73 
74  return a_f_v * retval;
75  }
76 
77  };
78 #endif
79 
80  /*! \brief Phase change equation (in terms of R^2) */
81  template <typename RT /*!< real-type */ >
82  struct dRsqdt
83  {
84  RT L; /*!< latent heat of vaporization (condensate) */
85  RT K; /*!< thermal conductivity (condensate) */
86  RT Rv; /*!< gas constant of air with vapour */
87  RT rho_l; /*!< density of condensate */
88 
89  /*! \brief Right-hand-side of the phase change ODE */
90  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
91  RT rhs_func ( const RT a_R_sq, /*!< radius squared */
92  const RT a_S, /*!< saturation ratio */
93  const RT a_T, /*!< temperature */
94  const RT a_e_s, /*!< saturation pressure */
95  const RT a_D, /*!< mol. diffusion coeff */
96  const RT a_a, /*!< curvature coeff */
97  const RT a_b, /*!< solute coeff */
98  const RT a_N_s, /*!< total solute moles */
99  const RT a_f_v /*!< ventilation factor */ ) const noexcept
100  {
101  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
102  RT Kn = lambda_v/std::sqrt(a_R_sq);
103  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
104 
105  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
106  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
107 
108  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
109  RT R_inv_cubed = R_inv*R_inv*R_inv;
110 
111  RT alpha = RT(2.0) * (a_S-RT(1.0)) / (F_k + F_d);
112  RT retval = alpha;
113 
114  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
115  retval += beta*R_inv;
116 
117  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
118  retval += gamma*R_inv_cubed;
119 
120  return a_f_v * retval;
121  }
122 
123  /*! \brief Jacobian of right-hand-side of the phase change ODE */
124  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
125  RT rhs_jac ( const RT a_R_sq, /*!< radius squared */
126  const RT a_T, /*!< temperature */
127  const RT a_e_s, /*!< saturation pressure */
128  const RT a_D, /*!< mol. diffusion coeff */
129  const RT a_a, /*!< curvature coeff */
130  const RT a_b, /*!< solute coeff */
131  const RT a_N_s, /*!< total solute moles */
132  const RT a_f_v /*!< ventilation factor */ ) const noexcept
133  {
134  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
135  RT Kn = lambda_v/std::sqrt(a_R_sq);
136  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
137 
138  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
139  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
140 
141  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
142  RT R_inv_3 = R_inv*R_inv*R_inv;
143  RT R_inv_5 = R_inv_3*R_inv*R_inv;
144 
145  RT retval = RT(0.0);
146 
147  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
148  retval -= myhalf * beta*R_inv_3;
149 
150  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
151  retval -= myhalf * three*gamma*R_inv_5;
152 
153  return a_f_v * retval;
154  }
155 
156  };
157 
158  /*! \brief Scalar Newton solver for phase change equation
159  *
160  * Solves the following nonlinear equation:
161  * mu * u - F(u) - R = 0,
162  * where:
163  * u: solution variable
164  * mu: constant
165  * R: right-hand-side (constant)
166  * F(u): function
167  */
168  template<typename NE /*!< Nonlinear equation */, typename RT /*!< real-type */>
169  struct NewtonSolver
170  {
171  const NE m_ne; /*!< nonlinear equation */
172 
173  RT m_rtol; /*!< relative tolerance */
174  RT m_atol; /*!< absolute tolerance */
175  RT m_stol; /*!< step size tolerance */
176  int m_maxits; /*!< max number of iterations */
177 
178  /*! \brief solve the nonlinear equation */
179  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
180  void operator() ( RT& a_u, /*!< solution variable */
181  RT& a_r, /*!< right-hand-side */
182  const RT& a_mu, /*!< mu */
183  const RT& a_S, /*!< saturation ratio */
184  const RT& a_T, /*!< temperature */
185  const RT& a_e_s, /*!< saturation pressure */
186  const RT& a_D, /*!< mol. diff. coeff */
187  const RT& a_a, /*!< curvature coeff */
188  const RT& a_b, /*!< solute coeff */
189  const RT& a_N_s, /*!< total solute moles */
190  const RT& a_f_v, /*!< ventilation factor */
191  RT& a_res_norm_a, /*!< absolute norm at exit */
192  RT& a_res_norm_r /*!< relative norm at exit */,
193  bool& a_converged /*!< convergence status at exit */
194  ) const
195  {
196  a_converged = false;
197  RT res_norm0 = RT(0.0);
198 
199  for (int k = 0; k < m_maxits; k++) {
200  RT residual = a_mu * a_u
201  - ( a_r
202  + m_ne.rhs_func( a_u, a_S, a_T, a_e_s, a_D, a_a, a_b, a_N_s, a_f_v ) );
203  a_res_norm_a = std::sqrt(residual*residual);
204 
205  if (k == 0) {
206  if (a_res_norm_a > 0) {
207  res_norm0 = a_res_norm_a;
208  } else {
209  res_norm0 = RT(1.0);
210  }
211  }
212  a_res_norm_r = a_res_norm_a / res_norm0;
213 
214  if (a_res_norm_a <= m_atol) {
215  a_converged = true;
216  break;
217  }
218  if (a_res_norm_r <= m_rtol) {
219  a_converged = true;
220  break;
221  }
222  if (!amrex::Math::isfinite(a_res_norm_a)) {
223  a_converged = false;
224  break;
225  }
226 
227  RT slope = a_mu - m_ne.rhs_jac( a_u, a_T, a_e_s, a_D, a_a, a_b, a_N_s, a_f_v );
228  RT du = RT(0.0);
229  du = - residual / slope;
230 
231  RT du_norm = std::sqrt(du*du);
232  RT u_norm = std::sqrt(a_u*a_u);
233  if (du_norm/u_norm <= m_stol) {
234  a_converged = true;
235  break;
236  }
237 
238  a_u += du;
239  if (a_u <= RT(0.0)) {
240  a_converged = false;
241  break;
242  }
243  }
244  }
245  };
246 
247  /*! \brief Implicit and explicit time integrators for the phase change equation */
248  template< typename ODE /*!< ODE */,
249  typename NewtonSolver /*!< Newton solver */,
250  /**
251  * @brief Implicit and explicit time integrators for the phase change equation.
252  *
253  * @tparam ODE ODE type.
254  * @tparam NewtonSolver Newton solver type.
255  * @tparam RT Real-type.
256  */
257  typename RT /*!< real-type */ >
258  struct TI
259  {
260  const ODE m_ode; /*!< ODE */
261  const NewtonSolver m_newton; /*!< Newton solver */
262 
263  RT m_t_final; /*!< final time */
264  RT m_max_steps; /*!< max number of timesteps */
265  RT m_S; /*!< saturation ratio */
266  RT m_T; /*!< temperature */
267  RT m_e_s; /*!< saturation pressure */
268  RT m_D; /*!< mol. diff. coeff */
269  RT m_a; /*!< coefficient of curvature */
270  RT m_b; /*!< coefficient of solute */
271  RT m_N_s; /*!< total solute moles */
272 
273  RT m_cfl; /*!< CFL */
274  RT m_atol; /*!< absolute tolerance (for adaptive dt) */
275  RT m_rtol; /*!< absolute tolerance (for adaptive dt) */
276  RT m_stol; /*!< solution update tolerance for exit due to steady state */
277 
278  bool m_adapt_dt; /*!< use error-based adaptive dt? */
279  bool m_verbose; /*!< verbosity */
280 
281  bool m_ventilation; /*!< include ventilation factor? */
282  RT m_vent_alpha1; /*!< ventilation factor fit coefficient alpha_1 */
283  RT m_vent_beta1; /*!< ventilation factor fit exponent beta_1 */
284  RT m_vent_alpha2; /*!< ventilation factor fit coefficient alpha_2 */
285  RT m_vent_beta2; /*!< ventilation factor fit exponent beta_2 */
286  RT m_vent_fcap; /*!< maximum ventilation factor */
287 
288  /*! \brief Ventilation factor at the given radius-squared, frozen over a
289  * step to keep the implicit Jacobian exact (unity if disabled) */
290  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
291  RT ventFactor (const RT& a_u) const
292  {
293  if (!m_ventilation) { return RT(1.0); }
294  return ventilationFactor( std::sqrt(a_u),
295  m_vent_alpha1, m_vent_beta1,
296  m_vent_alpha2, m_vent_beta2,
297  m_vent_fcap );
298  }
299 
300  /*! \brief Compute timestep from stiffness estimate */
301  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
302  RT computeTimestep (const RT& a_u) const
303  {
304  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
305  return m_cfl / std::sqrt(tau*tau);
306  }
307 
308  /*! \brief Compute stiffness estimate (tau) */
309  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
310  RT computeTau (const RT& a_u) const
311  {
312  return m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
313  }
314 
315  /*! \brief Limit timestep to not exceed final time */
316  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
317  RT limitTimestep (RT dt, RT cur_time) const
318  {
319  if ((cur_time + dt) > m_t_final) {
320  dt = m_t_final - cur_time;
321  }
322  return dt;
323  }
324 
325  /*! \brief Check if timestep is too small to continue */
326  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
327  bool isTimestepTooSmall (RT dt, RT tau) const
328  {
329  return (dt < (RT(1.0e-12) * m_cfl / std::sqrt(tau*tau)))
330  && (dt < (RT(1.0e-12) * m_t_final));
331  }
332 
333  /*! \brief Evaluate ODE right-hand side */
334  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
335  RT evalRHS (const RT& a_u) const
336  {
337  return m_ode.rhs_func(a_u, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
338  }
339 
340  /*! \brief Print verbose step information */
341  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
342  void printStepInfo (RT cur_time, RT dt, RT tau, RT radius, RT snorm) const
343  {
344  if (m_verbose) {
345  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
346  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
347  }
348  }
349 
350  /*! \brief Print verbose step info with Newton solver details */
351  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
352  void printStepInfoNewton (RT cur_time, RT dt, RT tau, RT radius, RT snorm,
353  RT res_norm_a, RT res_norm_r, bool converged) const
354  {
355  if (m_verbose) {
356  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
357  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
358  AMREX_DEVICE_PRINTF(" norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
359  res_norm_a, res_norm_r, (converged ? "yes" : "no"));
360  }
361  }
362 
363  /*! \brief 3rd-order, 4-stage Bogacki-Shampine explicit RK method
364  * with 2nd order embedded method */
365  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
366  void rk3bs ( RT& a_u, /*!< solution */
367  bool& a_success /*!< success/failure flag */ ) const
368  {
369  RT cur_time = RT(0.0);
370  a_success = true;
371 
372  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
373  RT dt = m_cfl / std::sqrt(tau*tau);
374 
375  RT dt_new = dt;
376  RT a_u_old = a_u;
377 
378  int n_step = 0;
379  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
380 
381  if (!m_adapt_dt) {
382  tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
383  dt = m_cfl / std::sqrt(tau*tau);
384  } else {
385  dt = dt_new;
386  }
387 
388  if ((cur_time + dt) > m_t_final) {
389  dt = m_t_final - cur_time;
390  }
391  if (!amrex::Math::isfinite(dt)) {
392  a_success = false;
393  break;
394  }
395 
396  RT u_new = RT(0.0);
397  bool step_success = false;
398  while (!step_success) {
399 
400  if ( (dt < (RT(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
401  && (dt < (RT(1.0e-12)*m_t_final)) ) {
402  break;
403  }
404 
405  RT u1 = a_u;
406  if (u1 <= 0) {
407  dt *= myhalf;
408  continue;
409  }
410  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
411 
412  RT u2 = a_u + myhalf*dt*f1;
413  if (u2 <= 0) {
414  dt *= myhalf;
415  continue;
416  }
417  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
418 
419  RT u3 = a_u + RT(0.75)*dt*f2;
420  if (u3 <= 0) {
421  dt *= myhalf;
422  continue;
423  }
424  RT f3 = m_ode.rhs_func(u3, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
425 
426  RT u4 = a_u + (RT(1.0)/RT(9.0))*dt * (RT(2.0)*f1 + three*f2 + RT(4.0)*f3);
427  if (u4 <= 0) {
428  dt *= myhalf;
429  continue;
430  }
431  RT f4 = m_ode.rhs_func(u4, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
432 
433  u_new = u4;
434 
435  if (m_adapt_dt) {
436  RT u_embed = a_u + (RT(1.0)/RT(24.0))*dt * (RT(7.0)*f1 + RT(6.0)*f2 + RT(8.0)*f3 + three*f4);
437  RT err = std::sqrt((u_new-u_embed)*(u_new-u_embed));
438  RT tol = m_atol + m_rtol * std::max(a_u, a_u_old);
439  RT E = err / tol;
440  dt_new = dt / std::cbrt(E);
441  }
442 
443  if (amrex::Math::isfinite(u_new)) {
444  if (u_new > 0) {
445  step_success = true;
446  break;
447  }
448  }
449  dt *= myhalf;
450  }
451 
452  if (step_success) {
453 
454  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
455  a_u_old = a_u;
456  a_u = u_new;
457  cur_time += dt;
458 
459  if (m_verbose) {
460  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
461  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
462  }
463  if (snorm < m_stol) {
464  break;
465  }
466 
467  } else {
468 
469  a_success = false;
470  break;
471 
472  }
473 
474  n_step++;
475  }
476 
477  return;
478  }
479 
480  /*! \brief 4th-order, 4-stage explicit Runge-Kutta method */
481  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
482  void rk4 ( RT& a_u, /*!< solution */
483  bool& a_success /*!< success/failure flag */ ) const
484  {
485  RT cur_time = RT(0.0);
486  a_success = true;
487 
488  int n_step = 0;
489  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
490 
491  RT tau = computeTau(a_u);
492  RT dt = m_cfl / std::sqrt(tau*tau);
493  dt = limitTimestep(dt, cur_time);
494  if (!amrex::Math::isfinite(dt)) {
495  a_success = false;
496  break;
497  }
498 
499  RT u_new = RT(0.0);
500  bool step_success = false;
501  while (!step_success) {
502 
503  if (isTimestepTooSmall(dt, tau)) { break; }
504 
505  RT u1 = a_u;
506  if (u1 <= 0) { dt *= myhalf; continue; }
507  RT f1 = evalRHS(u1);
508 
509  RT u2 = a_u + myhalf*dt*f1;
510  if (u2 <= 0) { dt *= myhalf; continue; }
511  RT f2 = evalRHS(u2);
512 
513  RT u3 = a_u + myhalf*dt*f2;
514  if (u3 <= 0) { dt *= myhalf; continue; }
515  RT f3 = evalRHS(u3);
516 
517  RT u4 = a_u + dt*f3;
518  if (u4 <= 0) { dt *= myhalf; continue; }
519  RT f4 = evalRHS(u4);
520 
521  u_new = a_u + dt*(f1 + RT(2.0)*f2 + RT(2.0)*f3 + f4)/RT(6.0);
522 
523  if (amrex::Math::isfinite(u_new) && (u_new > 0)) {
524  step_success = true;
525  break;
526  }
527  dt *= myhalf;
528  }
529 
530  if (step_success) {
531  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
532  a_u = u_new;
533  cur_time += dt;
534  printStepInfo(cur_time, dt, tau, std::sqrt(a_u), snorm);
535  if (snorm < m_stol) { break; }
536  } else {
537  a_success = false;
538  break;
539  }
540 
541  n_step++;
542  }
543 
544  return;
545  }
546 
547  /*! \brief 1st order implicit backward Euler method */
548  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
549  void be ( RT& a_u, /*!< solution */
550  bool& a_success /*!< success/failure flag */ ) const
551  {
552  RT cur_time = RT(0.0);
553  a_success = true;
554 
555  int n_step = 0;
556  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
557 
558  RT tau = computeTau(a_u);
559  RT dt = m_cfl / std::sqrt(tau*tau);
560  dt = limitTimestep(dt, cur_time);
561  if (!amrex::Math::isfinite(dt)) {
562  a_success = false;
563  break;
564  }
565 
566 #ifdef AMREX_USE_FLOAT
567  RT res_norm_a = std::numeric_limits<float>::max();
568  RT res_norm_r = std::numeric_limits<float>::max();
569 #else
570  RT res_norm_a = std::numeric_limits<double>::max();
571  RT res_norm_r = std::numeric_limits<double>::max();
572 #endif
573  bool converged = false;
574 
575  RT u_new = RT(0.0);
576  bool step_success = false;
577  while (!step_success) {
578 
579  if (isTimestepTooSmall(dt, tau)) { break; }
580 
581  RT mu = RT(1.0) / dt;
582  RT rhs = mu * a_u;
583  u_new = a_u;
584  m_newton(u_new, rhs, mu,
585  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
586  res_norm_a, res_norm_r, converged);
587 
588  if (converged && amrex::Math::isfinite(u_new) && (u_new > 0)) {
589  step_success = true;
590  break;
591  }
592  dt *= myhalf;
593  }
594 
595  if (step_success) {
596  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
597  a_u = u_new;
598  cur_time += dt;
599  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
600  res_norm_a, res_norm_r, converged);
601  if (snorm < m_stol) { break; }
602  } else {
603  a_success = false;
604  break;
605  }
606 
607  n_step++;
608  }
609 
610  return;
611  }
612 
613  /*! \brief 2nd-order implicit Crank-Nicolson method */
614  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
615  void cn ( RT& a_u, /*!< solution */
616  bool& a_success /*!< success/failure flag */ ) const
617  {
618  RT cur_time = RT(0.0);
619  a_success = true;
620 
621  int n_step = 0;
622  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
623 
624  RT tau = computeTau(a_u);
625  RT dt = m_cfl / std::sqrt(tau*tau);
626  dt = limitTimestep(dt, cur_time);
627  if (!amrex::Math::isfinite(dt)) {
628  a_success = false;
629  break;
630  }
631 
632  RT res_norm_a = std::numeric_limits<RT>::max();
633  RT res_norm_r = std::numeric_limits<RT>::max();
634  bool converged = false;
635 
636  RT u_new = RT(0.0);
637  bool step_success = false;
638  while (!step_success) {
639 
640  if (isTimestepTooSmall(dt, tau)) { break; }
641 
642  RT mu = RT(1.0) / (myhalf*dt);
643  RT u1 = a_u;
644  RT f1 = evalRHS(u1);
645 
646  RT u2 = u1;
647  RT rhs = mu * (a_u + myhalf*dt*f1);
648  m_newton(u2, rhs, mu,
649  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
650  res_norm_a, res_norm_r, converged);
651  if (u2 <= 0) { dt *= myhalf; continue; }
652  RT f2 = evalRHS(u2);
653 
654  u_new = a_u + myhalf * dt * (f1 + f2);
655 
656  if (converged && amrex::Math::isfinite(u_new) && (u_new > 0)) {
657  step_success = true;
658  break;
659  }
660  dt *= myhalf;
661  }
662 
663  if (step_success) {
664  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
665  a_u = u_new;
666  cur_time += dt;
667  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
668  res_norm_a, res_norm_r, converged);
669  if (snorm < m_stol) { break; }
670  } else {
671  a_success = false;
672  break;
673  }
674 
675  n_step++;
676  }
677 
678  return;
679  }
680 
681  /*! \brief 2nd-order, 2-stage diagonally-implicit RK method */
682  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
683  void dirk212 ( RT& a_u, /*!< solution */
684  bool& a_success /*!< success/failure flag */ ) const
685  {
686  RT cur_time = RT(0.0);
687  a_success = true;
688 
689  int n_step = 0;
690  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
691 
692  RT tau = computeTau(a_u);
693  RT dt = m_cfl / std::sqrt(tau*tau);
694  dt = limitTimestep(dt, cur_time);
695  if (!amrex::Math::isfinite(dt)) {
696  a_success = false;
697  break;
698  }
699 
700  bool converged = false;
701  RT res_norm_a = std::numeric_limits<RT>::max();
702  RT res_norm_r = std::numeric_limits<RT>::max();
703 
704  RT u_new = RT(0.0);
705  bool step_success = false;
706  while (!step_success) {
707 
708  if (isTimestepTooSmall(dt, tau)) {
709  break;
710  }
711  RT mu = RT(1.0) / dt;
712 
713  converged = true;
714  res_norm_a = RT(0.0);
715  res_norm_r = RT(0.0);
716 
717  RT u1 = a_u;
718  {
719  RT rhs = mu * a_u;
720  RT res_norm_a_i = std::numeric_limits<RT>::max();
721  RT res_norm_r_i = std::numeric_limits<RT>::max();
722  bool converged_i = false;
723  m_newton( u1, rhs, mu,
724  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
725  res_norm_a_i, res_norm_r_i, converged_i );
726  converged = converged && converged_i;
727  res_norm_a = std::max(res_norm_a, res_norm_a_i);
728  res_norm_r = std::max(res_norm_r, res_norm_r_i);
729  }
730  if (u1 <= 0) {
731  dt *= myhalf;
732  continue;
733  }
734  RT f1 = evalRHS(u1);
735 
736  RT u2 = u1;
737  {
738  RT rhs = mu * (a_u - dt*f1);
739  RT res_norm_a_i = std::numeric_limits<RT>::max();
740  RT res_norm_r_i = std::numeric_limits<RT>::max();
741  bool converged_i = false;
742  m_newton( u2, rhs, mu,
743  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
744  res_norm_a_i, res_norm_r_i, converged_i );
745  converged = converged && converged_i;
746  res_norm_a = std::max(res_norm_a, res_norm_a_i);
747  res_norm_r = std::max(res_norm_r, res_norm_r_i);
748  }
749  if (u2 <= 0) {
750  dt *= myhalf;
751  continue;
752  }
753  RT f2 = evalRHS(u2);
754 
755  u_new = a_u + myhalf * dt * (f1 + f2);
756 
757  if (converged) {
758  if (amrex::Math::isfinite(u_new)) {
759  if (u_new > 0) {
760  step_success = true;
761  break;
762  }
763  }
764  }
765  dt *= myhalf;
766  }
767 
768  if (step_success) {
769 
770  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
771  a_u = u_new;
772  cur_time += dt;
773 
774  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
775  res_norm_a, res_norm_r, converged);
776  if (snorm < m_stol) {
777  break;
778  }
779 
780  } else {
781 
782  a_success = false;
783  break;
784 
785  }
786 
787  n_step++;
788  }
789 
790  return;
791  }
792 
793  };
794 }
795 
796 /*! \brief Namespace with classes and functions for deposition/sublimation
797  * (solid <--> vapour) */
798 namespace SDMassChangeUtils_SV
799 {
800  /*! \brief Phase change equation (in terms of mass) */
801  template <typename RT /*!< real-type */ >
802  struct dMdt
803  {
804  RT L; /*!< latent heat of vaporization (condensate) */
805  RT K; /*!< thermal conductivity */
806  RT Rv; /*!< gas constant of air with vapour */
807  RT rho_ice; /*!< True density of ice */
808 
809  /*! \brief Compute growth ratio */
810  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
811  RT growthRatio ( const RT a_T /*!< temperature */ ) const noexcept
812  {
813  // Table of the inherent growth ratio (Chen and Lamb 1994)
814  // T(n) = (1.375-n)/4. Here, T is the temperature in [degreeC], and n is the index.
815  static constexpr amrex::Real tb_habit[121] = {
816  amrex::Real(1.000000e+00), amrex::Real(9.799412e-01), amrex::Real(9.600636e-01), amrex::Real(9.399397e-01), amrex::Real(9.200258e-01),
817  amrex::Real(9.001191e-01), amrex::Real(8.800351e-01), amrex::Real(8.570378e-01), amrex::Real(8.340652e-01), amrex::Real(8.109611e-01),
818  amrex::Real(7.830689e-01), amrex::Real(7.550922e-01), amrex::Real(7.030723e-01), amrex::Real(5.370318e-01), amrex::Real(4.677351e-01),
819  amrex::Real(5.248075e-01), amrex::Real(6.309573e-01), amrex::Real(8.128305e-01), amrex::Real(1.096478e+00), amrex::Real(1.479108e+00),
820  amrex::Real(1.905461e+00), amrex::Real(2.089296e+00), amrex::Real(2.290868e+00), amrex::Real(2.398833e+00), amrex::Real(2.454709e+00),
821  amrex::Real(2.426610e+00), amrex::Real(2.371374e+00), amrex::Real(2.290868e+00), amrex::Real(2.137962e+00), amrex::Real(1.995262e+00),
822  amrex::Real(1.862087e+00), amrex::Real(1.737801e+00), amrex::Real(1.621810e+00), amrex::Real(1.513561e+00), amrex::Real(1.396368e+00),
823  amrex::Real(1.288250e+00), amrex::Real(1.188502e+00), amrex::Real(1.096478e+00), amrex::Real(1.000000e+00), amrex::Real(9.225714e-01),
824  amrex::Real(8.511380e-01), amrex::Real(7.852356e-01), amrex::Real(7.244360e-01), amrex::Real(6.683439e-01), amrex::Real(6.165950e-01),
825  amrex::Real(5.754399e-01), amrex::Real(5.370318e-01), amrex::Real(5.011872e-01), amrex::Real(4.677351e-01), amrex::Real(4.365158e-01),
826  amrex::Real(4.073803e-01), amrex::Real(3.801894e-01), amrex::Real(3.548134e-01), amrex::Real(3.311311e-01), amrex::Real(3.162278e-01),
827  amrex::Real(3.019952e-01), amrex::Real(2.917427e-01), amrex::Real(2.851018e-01), amrex::Real(2.818383e-01), amrex::Real(2.786121e-01),
828  amrex::Real(2.754229e-01), amrex::Real(2.786121e-01), amrex::Real(2.818383e-01), amrex::Real(2.851018e-01), amrex::Real(2.917427e-01),
829  amrex::Real(2.985383e-01), amrex::Real(3.090295e-01), amrex::Real(3.198895e-01), amrex::Real(3.311311e-01), amrex::Real(3.467369e-01),
830  amrex::Real(3.672823e-01), amrex::Real(3.935501e-01), amrex::Real(4.265795e-01), amrex::Real(4.570882e-01), amrex::Real(4.897788e-01),
831  amrex::Real(5.248075e-01), amrex::Real(5.623413e-01), amrex::Real(6.095369e-01), amrex::Real(6.606934e-01), amrex::Real(7.161434e-01),
832  amrex::Real(7.852356e-01), amrex::Real(8.609938e-01), amrex::Real(9.549926e-01), amrex::Real(1.047129e+00), amrex::Real(1.148154e+00),
833  amrex::Real(1.258925e+00), amrex::Real(1.380384e+00), amrex::Real(1.496236e+00), amrex::Real(1.603245e+00), amrex::Real(1.698244e+00),
834  amrex::Real(1.778279e+00), amrex::Real(1.840772e+00), amrex::Real(1.883649e+00), amrex::Real(1.905461e+00), amrex::Real(1.905461e+00),
835  amrex::Real(1.883649e+00), amrex::Real(1.862087e+00), amrex::Real(1.840772e+00), amrex::Real(1.798871e+00), amrex::Real(1.737801e+00),
836  amrex::Real(1.698244e+00), amrex::Real(1.640590e+00), amrex::Real(1.584893e+00), amrex::Real(1.549173e+00), amrex::Real(1.513910e+00),
837  amrex::Real(1.476046e+00), amrex::Real(1.452112e+00), amrex::Real(1.428894e+00), amrex::Real(1.412863e+00), amrex::Real(1.393157e+00),
838  amrex::Real(1.376892e+00), amrex::Real(1.361131e+00), amrex::Real(1.348963e+00), amrex::Real(1.336903e+00), amrex::Real(1.327089e+00),
839  amrex::Real(1.317953e+00), amrex::Real(1.308881e+00), amrex::Real(1.302867e+00), amrex::Real(1.293898e+00), amrex::Real(1.287953e+00),
840  amrex::Real(1.279087e+00)
841  };
842  auto T_degC = a_T - RT(273.15); // [K] => [degC]
843  int n = std::min(std::max(static_cast<int>(std::round(RT(1.375)-RT(4.0)*T_degC)),1),121);
844  return static_cast<RT>(tb_habit[n-1]);
845 
846  }
847 
848  /*! \brief Compute F_k + F_d */
849  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
850  RT Fk_plus_Fd ( const RT a_T, /*!< temperature */
851  const RT a_e_s, /*!< saturation pressure */
852  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
853  {
854  auto F_k = (L/(Rv*a_T) - RT(1.0)) * (L / (K*a_T));
855  auto F_d = (Rv*a_T) / (a_D*a_e_s);
856  return (F_k + F_d);
857  }
858 
859  /*! \brief Compute dynamic viscosity */
860  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
861  RT viscCoeff ( const RT a_T /*!< temperature */ ) const noexcept
862  {
863  auto T_degC = a_T - RT(273.15); // [K] => [degC]
864  RT visc_coeff = RT(0.0);
865  if( T_degC >= RT(0.0) ) {
866  visc_coeff = ( RT(1.7180) + RT(4.9E-3)*T_degC ) * RT(1.E-5);
867  } else {
868  visc_coeff = ( RT(1.7180) + RT(4.9E-3)*T_degC -RT(1.2E-5)*T_degC*T_degC ) * RT(1.E-5);
869  }
870  return visc_coeff;
871  }
872 
873  /*! \brief Compute capacitance */
874  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
875  RT capacitance( const RT a_a, /*!< equatorial radius */
876  const RT a_c /*!< polar radius */ ) const noexcept
877  {
878  auto phi = a_c / a_a;
879  RT capaci = RT(0.0);
880  if(phi >= RT(1.001)) {
881  auto ecent = std::sqrt(RT(1.0) - RT(1.0)/(phi*phi));
882  capaci = ecent*a_c / std::log((RT(1.0)+ecent)*phi);
883  } else if (phi <= RT(0.999)) {
884  auto ecent = std::sqrt(RT(1.0) - phi*phi);
885  capaci = ecent*a_a / std::asin(ecent);
886  } else {
887  capaci = (RT(2.0)*a_a + a_c)/three;
888  }
889  return capaci;
890  }
891 
892  /*! \brief Compute ventilation effect */
893  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
894  RT ventX ( const RT a_a, /*!< equatorial radius */
895  const RT a_c, /*!< polar radius */
896  const RT a_vterm,/*!< terminal velocity */
897  const RT a_T, /*!< temperature */
898  const RT a_rhom, /*!< moist density */
899  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
900  {
901  auto visc_coeff = viscCoeff(a_T);
902  auto maxdim = RT(2.0) * std::max(a_a, a_c);
903  auto n_re = a_rhom * maxdim * a_vterm / visc_coeff;
904  auto n_sc = visc_coeff / (a_rhom*a_D);
905  auto X_vent = std::cbrt(n_sc) * std::sqrt(n_re);
906  return X_vent;
907  }
908 
909  /*! \brief Radius change */
910  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
911  RT growthRatioStar ( const RT a_dmass, /*!< delta mass */
912  const RT a_a, /*!< equatorial radius */
913  const RT a_c, /*!< polar radius */
914  const RT a_vterm, /*!< terminal velocity */
915  const RT a_rhom, /*!< moist air density */
916  const RT a_T, /*!< temperature */
917  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
918  {
919  auto capaci = capacitance(a_a, a_c);
920  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
921  RT fratio = RT(0.0);
922  if(X_vent <= RT(1.0)) {
923  fratio = (RT(1.0) + RT(0.14)*X_vent*X_vent*std::sqrt(a_c/capaci))
924  / (RT(1.0) + RT(0.14)*X_vent*X_vent*std::sqrt(a_a/capaci));
925  } else {
926  fratio = (RT(0.86) + RT(0.28)*X_vent*std::sqrt(a_c/capaci))
927  / (RT(0.86) + RT(0.28)*X_vent*std::sqrt(a_a/capaci));
928  }
929 
930  auto gr = growthRatio(a_T);
931  if ((RT(2.0)*std::max(a_a,a_c) < RT(1.0e-5)) || (a_dmass <= RT(0.0))) { gr = RT(1.0); }
932  auto gr_star = gr * fratio;
933  if((a_dmass > RT(0.0)) && (a_c/a_a > RT(40.0))) { gr_star = std::min(gr_star,RT(1.0)); }
934 
935  return gr_star;
936  }
937 
938  /*! \brief Volume change */
939  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
940  RT dVolume ( const RT a_dmass, /*!< delta mass */
941  const RT a_a, /*!< equatorial radius */
942  const RT a_c, /*!< polar radius */
943  const RT a_rho, /*!< ice density of particle */
944  const RT a_S, /*!< saturation ratio */
945  const RT a_T, /*!< temperature */
946  const RT a_e_s, /*!< saturation pressure */
947  const RT a_es_ratio_wi, /*!< ratio of saturation pressures of water to ice */
948  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
949  {
950  RT rho_i_ds = a_rho;
951  auto gr = growthRatio(a_T);
952  if (a_dmass > RT(0.0)) {
953  if ((gr < RT(1.0)) && (a_a < RT(1.0e-4))) {
954  rho_i_ds = rho_ice;
955  } else {
956  if (RT(2.0)*std::max(a_a,a_c) < RT(1.0e-5)) { gr = RT(1.0); }
957  auto ex_vap_dens = RT(1000.0) * (std::min(a_S,a_es_ratio_wi) - RT(1.0)) / (Fk_plus_Fd(a_T,a_e_s,a_D) * a_D);
958  rho_i_ds = rho_ice * std::exp(-three*std::max(ex_vap_dens-RT(0.05),RT(0.0))/gr);
959  }
960  }
961 
962  auto dvol = a_dmass / rho_i_ds;
963  return dvol;
964  }
965 
966  /*! \brief Radius change */
967  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
968  RT dLogRadius ( const RT a_grstar, /*!< growth ratio */
969  const RT a_dlogvol /*!< delta log volume */ ) const noexcept
970  {
971  return a_dlogvol / (a_grstar+RT(2.0));
972  }
973 
974  /*! \brief Right-hand-side of the phase change ODE */
975  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
976  RT rhs_func ( const RT a_a, /*!< equatorial radius */
977  const RT a_c, /*!< polar radius */
978  const RT a_vterm,/*!< terminal velocity */
979  const RT a_S, /*!< saturation ratio */
980  const RT a_rhom, /*!< moist density */
981  const RT a_T, /*!< temperature */
982  const RT a_e_s, /*!< saturation pressure */
983  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
984  {
985  auto capaci = capacitance(a_a, a_c);
986  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
987 
988  RT f_vent = RT(0.0);
989  if (X_vent <= RT(1.0)) { f_vent = RT(1.0) + RT(0.14)*X_vent*X_vent; }
990  else { f_vent = RT(0.86) + RT(0.28)*X_vent; }
991 
992  auto retval = RT(4.0)*PI * capaci * (a_S-RT(1.0)) * f_vent / Fk_plus_Fd(a_T, a_e_s, a_D);
993  return retval;
994  }
995 
996  /*! \brief Ice melt rate dm_ice/dt (Pruppacher-Klett; Seifert-Beheng 2006 Eq. 72).
997  * Negative when melting. Heat conducted from warm air plus latent heat of the
998  * vapour flux to the 0 degC wet surface, divided by the latent heat of fusion.
999  * Reuses capacitance and ventilation. For this functor built with the water
1000  * material, L is the latent heat of vaporization, K the air thermal
1001  * conductivity, and Rv the vapour gas constant. */
1002  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1003  RT meltRate ( const RT a_a, /*!< equatorial radius */
1004  const RT a_c, /*!< polar radius */
1005  const RT a_vterm, /*!< terminal velocity */
1006  const RT a_T, /*!< temperature */
1007  const RT a_e_inf, /*!< ambient vapour pressure */
1008  const RT a_e_sat_wT0, /*!< sat. pressure over water at melting point */
1009  const RT a_rhom, /*!< moist density */
1010  const RT a_D, /*!< vapour diffusion coeff */
1011  const RT a_L_f /*!< latent heat of fusion */ ) const noexcept
1012  {
1013  auto capaci = capacitance(a_a, a_c);
1014  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
1015 
1016  RT f_vent = RT(0.0);
1017  if (X_vent <= RT(1.0)) { f_vent = RT(1.0) + RT(0.14)*X_vent*X_vent; }
1018  else { f_vent = RT(0.86) + RT(0.28)*X_vent; }
1019 
1020  RT T0 = RT(tmelt);
1021  RT heat = K * f_vent * (a_T - T0)
1022  + (L*a_D/Rv) * (a_e_inf/a_T - a_e_sat_wT0/T0) * f_vent;
1023  return -(RT(4.0)*PI*capaci/a_L_f) * heat;
1024  }
1025 
1026  };
1027 
1028  /*! \brief Implicit and explicit time integrators for the phase change equation */
1029  template< typename ODE /*!< ODE */,
1030  typename RT /*!< real-type */ >
1031  struct TI
1032  {
1033  const ODE m_ode; /*!< ODE */
1034 
1035  RT m_t_final; /*!< final time */
1036  mutable RT m_dt;/*!< dt */
1037  RT m_max_steps; /*!< max number of timesteps */
1038 
1039  RT m_a; /*!< equatorial radius */
1040  RT m_c; /*!< polar radius */
1041  RT m_vt; /*!< terminal velocity */
1042  RT m_S; /*!< saturation ratio */
1043  RT m_rhom; /*!< moist density */
1044  RT m_T; /*!< temperature */
1045  RT m_e_s; /*!< saturation pressure */
1046  RT m_D; /*!< mol. diff. coeff */
1047 
1048  bool m_verbose; /*!< verbosity */
1049 
1050  /*! \brief 1st-order, forward Euler method */
1051  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1052  void fe ( RT& a_u, /*!< solution */
1053  bool& a_success /*!< success/failure flag */ ) const
1054  {
1055  RT cur_time = RT(0.0);
1056  a_success = true;
1057 
1058  int n_step = 0;
1059  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
1060 
1061  if ((cur_time + m_dt) > m_t_final) {
1062  m_dt = m_t_final - cur_time;
1063  }
1064 
1065  RT rhs = m_ode.rhs_func(m_a, m_c, m_vt, m_S, m_rhom, m_T, m_e_s, m_D);
1066  RT u_new = a_u + m_dt*rhs;
1067  if (!amrex::Math::isfinite(u_new)) {
1068  a_success = false;
1069  break;
1070  }
1071 
1072  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
1073  a_u = u_new;
1074  cur_time += m_dt;
1075 
1076  if (m_verbose) {
1077  AMREX_DEVICE_PRINTF( "Time %1.2e, m_dt = %1.2e, mass = %1.4e, snorm = %1.1e\n",
1078  cur_time, m_dt, a_u, snorm);
1079  }
1080 
1081  n_step++;
1082  }
1083 
1084  return;
1085  }
1086 
1087  };
1088 }
1089 
1090 #endif
1091 #endif
gamma
amrex::Real gamma
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:9
beta
amrex::Real beta
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:10
ERF_MicrophysicsConstants.H
Physical constants and tuning parameters used only by the moisture and cloud-physics code.
tmelt
constexpr amrex::Real tmelt
Definition: ERF_MicrophysicsConstants.H:122
ERF_NumericalConstants.H
Dimensionless numeric literals and pure mathematical constants.
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:32
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
PI
constexpr amrex::Real PI
Definition: ERF_NumericalConstants.H:39
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MORRInd::mu
@ mu
Definition: ERF_AdvanceMorrison.cpp:92
mp_wdm6::alpha
real(kind=kind_phys), parameter, private alpha
Definition: ERF_module_mp_wdm6.F90:62