ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Morrison Class Reference

#include <ERF_Morrison.H>

Inheritance diagram for Morrison:
Collaboration diagram for Morrison:

Public Member Functions

 Morrison ()
 
virtual ~Morrison ()=default
 
void Define (SolverChoice &sc) override
 
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
 
void Set_dzmin (const amrex::Real dz_min) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &sc) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
int Qstate_Moist_NumConc_Size () override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const override
 Populate a vector with names of all available Morrison plot variables. More...
 
void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf) const override
 
void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const override
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void SetCurrentLevel (const int &)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_RealWidth (const int)
 

Private Types

using FabPtr = std::shared_ptr< amrex::MultiFab >
 

Private Attributes

int m_qmoist_size = 3
 
int n_qstate_moist_size = 11
 
int n_qstate_moist_numconc_size = 5
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::Real m_rdOcp
 
bool m_do_cond
 
amrex::Real m_dzmin
 
amrex::MultiFab * m_z_phys_nd
 
amrex::MultiFab * m_detJ_cc
 
amrex::Array< FabPtr, MicVar_Morr::NumVarsmic_fab_vars
 

Member Typedef Documentation

◆ FabPtr

using Morrison::FabPtr = std::shared_ptr<amrex::MultiFab>
private

Constructor & Destructor Documentation

◆ Morrison()

Morrison::Morrison ( )
inline
64 {}

◆ ~Morrison()

virtual Morrison::~Morrison ( )
virtualdefault

Member Function Documentation

◆ Advance()

void Morrison::Advance ( const amrex::Real dt_advance,
const SolverChoice sc 
)
overridevirtual

Reimplemented from NullMoist.

166  {
167  // Expose for GPU
168  bool do_cond = m_do_cond;
169 
170  // Store timestep
171  Real dt = dt_advance;
172 
173  // Check if CPP or FORT answer is used
174  ParmParse pp("erf");
175  bool use_morr_cpp_answer = true;
176  pp.query("use_morr_cpp_answer", use_morr_cpp_answer);
177 
178  // Ensure that only one of these is true
179  bool run_morr_cpp = use_morr_cpp_answer;
180  bool run_morr_fort = !run_morr_cpp;
181 
182  std::string filename = std::string("output_cpp") + std::to_string(use_morr_cpp_answer) + ".txt";
183 
184  // Allow user to override constant droplet concentration from inputs file
185  // Constant droplet concentration (if INUM = 1)
186  Real m_ndcnst = Real(250.0); // Droplet number concentration (cm^-3)
187  pp.query("morrison_ndcnst", m_ndcnst);
188 
189  // Loop through the grids
190  for (MFIter mfi(*mic_fab_vars[MicVar_Morr::qcl],TileNoZ()); mfi.isValid(); ++mfi)
191  {
192  auto box = mfi.tilebox();
193 
194  if (!box.ok()) { // Avoid going farther if the box is inverted (i.e., ilo > ihi or jlo > jhi).
195  continue;
196  }
197 
198  // Get array data from class member variables
199  auto const& theta_arr = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
200 
201  auto const& qv_arr = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
202 
203  auto const& qcl_arr = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
204  auto const& qpr_arr = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
205  auto const& qci_arr = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
206  auto const& qps_arr = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
207  auto const& qpg_arr = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
208 
209  auto const& nc_arr = mic_fab_vars[MicVar_Morr::nc]->array(mfi);
210  auto const& ni_arr = mic_fab_vars[MicVar_Morr::ni]->array(mfi);
211  auto const& nr_arr = mic_fab_vars[MicVar_Morr::nr]->array(mfi);
212  auto const& ns_arr = mic_fab_vars[MicVar_Morr::ns]->array(mfi);
213  auto const& ng_arr = mic_fab_vars[MicVar_Morr::ng]->array(mfi);
214 
215 #ifdef ERF_USE_MORR_FORT
216  auto const& rho_arr = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
217 #endif
218  auto const& pres_arr = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
219  auto const& rain_accum_arr = mic_fab_vars[MicVar_Morr::rain_accum]->array(mfi);
220  auto const& snow_accum_arr = mic_fab_vars[MicVar_Morr::snow_accum]->array(mfi);
221  auto const& graup_accum_arr = mic_fab_vars[MicVar_Morr::graup_accum]->array(mfi);
222  auto const& w_arr = mic_fab_vars[MicVar_Morr::omega]->array(mfi);
223 
224  // Get radar reflectivity array if radar diagnostics enabled
225  // auto const& refl_arr = m_do_radar_ref ? m_radar->array(mfi) : nullptr;
226  // auto const& refl_arr = m_radar->array(mfi);
227 
228  // Extract box dimensions
229  const int ilo = box.loVect()[0];
230  const int ihi = box.hiVect()[0];
231  const int jlo = box.loVect()[1];
232  const int jhi = box.hiVect()[1];
233  const int klo = box.loVect()[2];
234  const int khi = box.hiVect()[2];
235 
236  Box grown_box(box); grown_box.grow(3);
237 
238 #if defined(ERF_USE_MORR_FORT) && defined(AMREX_USE_GPU)
239  Arena* Arena_Used = The_Pinned_Arena();
240 #else
241  Arena* Arena_Used = The_Async_Arena();
242 #endif
243 
244  // Calculate Exner function (PII) to convert potential temperature to temperature
245  // PII = (P/P0)^(R/cp)
246  FArrayBox pii_fab(grown_box, 1, Arena_Used);
247  auto const& pii_arr = pii_fab.array();
248 
249  const Real p0 = Real(100000.0); // Reference pressure (Pa)
250 
251  const Real rdcp = m_rdOcp; // R/cp ratio
252 
253  // Calculate Exner function
254  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
255  // NOTE: the Morrison Fortran version uses Pa not hPa so we didn't divide p by 100
256  // so we don't need to multiply by 100 here
257  pii_arr(i,j,k) = std::pow((pres_arr(i,j,k)) / p0, rdcp);
258  });
259 
260  // Create arrays for height differences (dz)
261  FArrayBox dz_fab(grown_box, 1, Arena_Used);
262  auto const& dz_arr = dz_fab.array();
263 
264  // Calculate height differences
265  const Real dz_val = m_geom.CellSize(2);
266  const Array4<const Real> z_arr = (m_z_phys_nd) ? m_z_phys_nd->const_array(mfi) : Array4<const Real> {};
267  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
268  dz_arr(i,j,k) = (z_arr) ? Real(0.25) * ( (z_arr(i ,j ,k+1) - z_arr(i ,j ,k))
269  + (z_arr(i+1,j ,k+1) - z_arr(i+1,j ,k))
270  + (z_arr(i ,j+1,k+1) - z_arr(i ,j+1,k))
271  + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz_val;
272  });
273 
274  Box grown_boxD(grown_box); grown_boxD.makeSlab(2,0);
275 
276  // Arrays to store precipitation rates
277  FArrayBox rainncv_fab(grown_boxD, 1, Arena_Used);
278  FArrayBox sr_fab(grown_boxD, 1, Arena_Used); // Ratio of snow to total precipitation
279  FArrayBox snowncv_fab(grown_boxD, 1, Arena_Used);
280  FArrayBox graupelncv_fab(grown_boxD, 1, Arena_Used);
281 
282  auto const& rainncv_arr = rainncv_fab.array();
283  auto const& sr_arr = sr_fab.array();
284  auto const& snowncv_arr = snowncv_fab.array();
285  auto const& graupelncv_arr = graupelncv_fab.array();
286 
287  // Initialize precipitation rate arrays to Real(0)
288  ParallelFor(grown_boxD, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
289  rainncv_arr(i,j,k) = Real(0);
290  sr_arr(i,j,k) = Real(0);
291  snowncv_arr(i,j,k) = Real(0);
292  graupelncv_arr(i,j,k) = Real(0);
293  });
294 
295  // Create terrain height array (not actually used by Morrison scheme)
296  FArrayBox ht_fab(Box(IntVect(ilo, jlo, 0), IntVect(ihi, jhi, 0)), 1, Arena_Used);
297  [[maybe_unused]] auto const& ht_arr = ht_fab.array();
298  // ParallelFor(Box(IntVect(ilo, jlo, 0), IntVect(ihi, jhi, 0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
299  // ht_arr(i,j,k) = (z_arr) ? Real(0.25) * ( z_arr(i ,j ,k) + z_arr(i+1,j ,k)
300  // + z_arr(i ,j+1,k) + z_arr(i+1,j+1,k) ) : Real(0.); // Not used by Morrison scheme
301  //});
302 
303  // Microphysics options/switches
304  int m_iact = 2; // CCN activation option (1:std::power-law, 2: lognormal aerosol)
305  int m_inum = 1; // Droplet number option (0: predict, 1: constant)
306 
307  int m_iliq = 0; // Liquid-only option (0: include ice, 1: liquid only)
308  int m_inuc = 0; // Ice nucleation option (0: mid-latitude, 1: arctic)
309  // int m_ibase = 2; // Cloud base activation option
310  // int m_isub = 0; // Sub-grid vertical velocity option
311  int m_igraup = 0; // Graupel option (0: include graupel, 1: no graupel)
312  int m_ihail = 0; // Graupel/hail option (0: graupel, 1: hail)
313 
314  if(sc.moisture_type == MoistureType::Morrison_NoIce) {
315  m_iliq = 1; // Liquid-only option (0: include ice, 1: liquid only)
316  m_inuc = 0; // Ice nucleation option (0: mid-latitude, 1: arctic)
317  // m_ibase = 2; // Cloud base activation option
318  // m_isub = 0; // Sub-grid vertical velocity option
319  m_igraup = 1; // Graupel option (0: include graupel, 1: no graupel)
320  m_ihail = 0; // Graupel/hail option (0: graupel, 1: hail)
321  }
322  // bool m_do_radar_ref = false; // Radar reflectivity calculation flag
323 
324  // Physical constants
325  Real m_pi; // Pi constant
326  Real m_R; // Gas constant for dry air (J/kg/K)
327  Real m_Rd; // Gas constant for dry air (J/kg/K)
328  Real m_Rv; // Gas constant for water vapor (J/kg/K)
329  // Real m_cp; // Specific heat at constant pressure (J/kg/K)
330  Real m_g; // Gravitational acceleration (m/s^2)
331  Real m_ep_2; // Molecular weight ratio (Rd/Rv)
332 
333  // Reference density and species densities
334  Real m_rhosu; // Standard air density at 850 mb (kg/m^3)
335  Real m_rhow; // Density of liquid water (kg/m^3)
336  Real m_rhoi; // Bulk density of cloud ice (kg/m^3)
337  Real m_rhosn; // Bulk density of snow (kg/m^3)
338  Real m_rhog; // Bulk density of graupel/hail (kg/m^3)
339 
340  // Fall speed parameters (V=AD^B)
341  Real m_ai, m_bi; // Cloud ice fall speed parameters
342  // Real m_ac; // Cloud droplet fall speed parameters
343  Real m_bc; // Cloud droplet fall speed parameters
344  Real m_as, m_bs; // Snow fall speed parameters
345  Real m_ar, m_br; // Rain fall speed parameters
346  Real m_ag, m_bg; // Graupel/hail fall speed parameters
347 
348  // Microphysical parameters
349  Real m_aimm; // Parameter in Bigg immersion freezing
350  Real m_bimm; // Parameter in Bigg immersion freezing
351  Real m_ecr; // Collection efficiency between droplets/rain and snow/rain
352  Real m_dcs; // Threshold size for cloud ice autoconversion (m)
353  Real m_mi0; // Initial mass of nucleated ice crystal (kg)
354  Real m_mg0; // Mass of embryo graupel (kg)
355  Real m_f1s; // Ventilation parameter for snow
356  Real m_f2s; // Ventilation parameter for snow
357  Real m_f1r; // Ventilation parameter for rain
358  Real m_f2r; // Ventilation parameter for rain
359  Real m_qsmall; // Smallest allowed hydrometeor mixing ratio
360  Real m_eii; // Collection efficiency, ice-ice collisions
361  Real m_eci; // Collection efficiency, ice-droplet collisions
362  Real m_cpw; // Specific heat of liquid water (J/kg/K)
363  Real m_rin; // Radius of contact nuclei (m)
364  Real m_mmult; // Mass of splintered ice particle (kg)
365 
366  // Size distribution parameters
367  Real m_ci, m_di; // Cloud ice size distribution parameters
368  Real m_cs, m_ds; // Snow size distribution parameters
369  Real m_cg, m_dg; // Graupel size distribution parameters
370 
371  // Lambda limits for size distributions
372  Real m_lammaxi, m_lammini; // Cloud ice lambda limits
373  Real m_lammaxr, m_lamminr; // Rain lambda limits
374  Real m_lammaxs, m_lammins; // Snow lambda limits
375  Real m_lammaxg, m_lamming; // Graupel lambda limits
376 
377  // CCN spectra parameters (for IACT = 1)
378  // Real m_k1; // Exponent in CCN activation formula
379  // Real m_c1; // Coefficient in CCN activation formula (cm^-3)
380 
381  // Aerosol activation parameters (for IACT = 2)
382  // Real m_mw; // Molecular weight water (kg/mol)
383  // Real m_osm; // Osmotic coefficient
384  // Real m_vi; // Number of ions dissociated in solution
385  // Real m_epsm; // Aerosol soluble fraction
386  // Real m_rhoa; // Aerosol bulk density (kg/m^3)
387  // Real m_map; // Molecular weight aerosol (kg/mol)
388  // Real m_ma; // Molecular weight of air (kg/mol)
389  // Real m_rr; // Universal gas constant (J/mol/K)
390  // Real m_bact; // Activation parameter
391  // Real m_rm1; // Geometric mean radius, mode 1 (m)
392  // Real m_rm2; // Geometric mean radius, mode 2 (m)
393  Real m_nanew1; // Total aerosol concentration, mode 1 (m^-3)
394  Real m_nanew2; // Total aerosol concentration, mode 2 (m^-3)
395  // Real m_sig1; // Standard deviation of aerosol dist, mode 1
396  // Real m_sig2; // Standard deviation of aerosol dist, mode 2
397  // Real m_f11; // Correction factor for activation, mode 1
398  // Real m_f12; // Correction factor for activation, mode 1
399  // Real m_f21; // Correction factor for activation, mode 2
400  // Real m_f22; // Correction factor for activation, mode 2
401 
402  // Precomputed constants for efficiency
403  Real m_cons1, m_cons2, m_cons3, m_cons4, m_cons5;
404  Real m_cons6, m_cons7, m_cons8, m_cons9, m_cons10;
405  Real m_cons11, m_cons12, m_cons13, m_cons14, m_cons15;
406  Real m_cons16, m_cons17, m_cons18, m_cons19, m_cons20;
407  Real m_cons21, m_cons22, m_cons23, m_cons24, m_cons25;
408  Real m_cons26, m_cons27, m_cons28, m_cons29;
409  Real m_cons31, m_cons32, m_cons34, m_cons35;
410  Real m_cons36, m_cons37, m_cons38, m_cons39, m_cons40;
411  Real m_cons41;
412 
413  // Set microphysics control parameters
414  m_inum = 1; // Use constant droplet number concentration
415  m_ndcnst = Real(250.0); // Droplet number concentration (cm^-3)
416  // Mathematical constants
417  m_pi = Real(3.1415926535897932384626434);
418 
419  m_R = Real(287.0); // Gas constant for dry air (J/kg/K)
420  m_Rd = Real(287.0); // Gas constant for dry air (J/kg/K)
421  m_Rv = Real(461.6); // Gas constant for water vapor (J/kg/K)
422  // m_cp = Real(7.0)*Real(287.0)/Real(2); // Specific heat at constant pressure (J/kg/K)
423  m_g = Real(9.81); // Gravitational acceleration (m/s^2)
424  m_ep_2 = m_Rd / m_Rv; // Molecular weight ratio (Rd/Rv)
425 
426  // Reference density
427  m_rhosu = Real(85000.0)/(Real(287.15)*Real(273.15)); // Standard air density at 850 mb (kg/m^3)
428 
429  // Densities for different hydrometeor species
430  m_rhow = Real(997.0); // Density of liquid water (kg/m^3)
431  m_rhoi = Real(500.0); // Bulk density of cloud ice (kg/m^3)
432  m_rhosn = Real(100.0); // Bulk density of snow (kg/m^3)
433 
434  // Set density for graupel or hail based on configuration
435  if (m_ihail == 0) {
436  m_rhog = Real(400.0); // Bulk density of graupel (kg/m^3)
437  } else {
438  m_rhog = Real(900.0); // Bulk density of hail (kg/m^3)
439  }
440 
441  // Fall speed parameters (V=AD^B) for different hydrometeors
442  // Cloud ice
443  m_ai = Real(700.0);
444  m_bi = one;
445 
446  // Cloud droplets
447  // m_ac = Real(3.0E7);
448  m_bc = Real(2);
449 
450  // Snow
451  m_as = Real(11.72);
452  m_bs = Real(0.41);
453 
454  // Rain
455  m_ar = Real(841.99667);
456  m_br = Real(0.8);
457 
458  // Graupel/hail (dependent on configuration)
459  if (m_ihail == 0) {
460  // Graupel parameters
461  m_ag = Real(19.3);
462  m_bg = Real(0.37);
463  } else {
464  // Hail parameters (Matsun and Huggins 1980)
465  m_ag = Real(114.5);
466  m_bg = myhalf;
467  }
468 
469  // Microphysical parameters
470  m_aimm = Real(0.66); // Parameter in Bigg immersion freezing
471  m_bimm = Real(100.0); // Parameter in Bigg immersion freezing
472  m_ecr = one; // Collection efficiency between rain and snow/graupel
473  m_dcs = Real(125.0E-6); // Threshold size for cloud ice autoconversion (m)
474  m_mi0 = Real(4.0)/three*m_pi*m_rhoi*amrex::Math::powi<3>(Real(10.0E-6)); // Initial mass of nucleated ice crystal (kg)
475  m_mg0 = Real(1.6E-10); // Mass of embryo graupel (kg)
476 
477  // Ventilation parameters
478  m_f1s = Real(0.86); // Ventilation parameter for snow
479  m_f2s = Real(0.28); // Ventilation parameter for snow
480  m_f1r = Real(0.78); // Ventilation parameter for rain
481  m_f2r = Real(0.308); // Ventilation parameter for rain
482 
483  // Smallest allowed hydrometeor mixing ratio
484  m_qsmall = Real(1.0E-14);
485 
486  // Collection efficiencies
487  m_eii = Real(0.1); // Ice-ice collision efficiency
488  m_eci = Real(0.7); // Ice-droplet collision efficiency
489 
490  // Specific heat of liquid water (J/kg/K)
491  m_cpw = Real(4187.0);
492 
493  // Size distribution parameters
494  m_ci = m_rhoi * m_pi / Real(6.0);
495  m_di = three;
496  m_cs = m_rhosn * m_pi / Real(6.0);
497  m_ds = three;
498  m_cg = m_rhog * m_pi / Real(6.0);
499  m_dg = three;
500 
501  // Radius of contact nuclei (m)
502  m_rin = Real(0.1E-6);
503 
504  // Mass of splintered ice particle (kg)
505  m_mmult = Real(4.0)/three*m_pi*m_rhoi*amrex::Math::powi<3>(Real(5.0E-6));
506 
507  // Set lambda limits for size distributions
508  // Maximum and minimum values for lambda parameter in size distributions
509  m_lammaxi = one/Real(1.0E-6);
510  m_lammini = one/(Real(2)*m_dcs + Real(100.0E-6));
511  m_lammaxr = one/Real(20.0E-6);
512  m_lamminr = one/Real(2800.0E-6);
513  m_lammaxs = one/Real(10.0E-6);
514  m_lammins = one/Real(2000.0E-6);
515  m_lammaxg = one/Real(20.0E-6);
516  m_lamming = one/Real(2000.0E-6);
517 
518  // Set CCN parameters for different environments
519  if (m_iact == 1) {
520  // Maritime CCN spectrum parameters (modified from Rasmussen et al. 2002)
521  // NCCN = C*S^K, where S is supersaturation in %
522  // m_k1 = Real(0.4); // Exponent in CCN activation formula
523  // m_c1 = Real(120.0); // Coefficient in CCN activation formula (cm^-3)
524  }
525 
526  // Initialize aerosol activation parameters for lognormal distribution
527  if (m_iact == 2) {
528  // Parameters for ammonium sulfate
529  // m_mw = Real(0.018); // Molecular weight of water (kg/mol)
530  // m_osm = one; // Osmotic coefficient
531  // m_vi = three; // Number of ions dissociated in solution
532  // m_epsm = Real(0.7); // Aerosol soluble fraction
533  // m_rhoa = Real(1777.0); // Aerosol bulk density (kg/m^3)
534  // m_map = Real(0.132); // Molecular weight of aerosol (kg/mol)
535  // m_ma = Real(0.0284); // Molecular weight of air (kg/mol)
536  // m_rr = Real(8.3145); // Universal gas constant (J/mol/K)
537  // m_bact = m_vi * m_osm * m_epsm * m_mw * m_rhoa / (m_map * m_rhow);
538  // m_a_w = two * m_mw * Real(0.0761) / (m_rhow * m_r_v * Real(293.15)); // "A" parameter
539 
540  // Aerosol size distribution parameters for MPACE (Morrison et al. 2007, JGR)
541  // Mode 1
542  // m_rm1 = Real(0.052E-6); // Geometric mean radius, mode 1 (m)
543  // m_sig1 = Real(2.04); // Standard deviation of aerosol size distribution, mode 1
544  m_nanew1 = Real(72.2E6); // Total aerosol concentration, mode 1 (m^-3)
545  // m_f11 = myhalf * std::exp(Real(2.5) * amrex::Math::powi<2>(std::log(m_sig1)));
546  // m_f21 = one + fourth * std::log(m_sig1);
547 
548  // Mode 2
549  // m_rm2 = Real(1.3E-6); // Geometric mean radius, mode 2 (m)
550  // m_sig2 = Real(2.5); // Standard deviation of aerosol size distribution, mode 2
551  m_nanew2 = Real(1.8E6); // Total aerosol concentration, mode 2 (m^-3)
552  // m_f12 = myhalf * std::exp(Real(2.5) * amrex::Math::powi<2>(std::log(m_sig2)));
553  // m_f22 = one + fourth * std::log(m_sig2);
554  }
555 
556  // Precompute constants for efficiency
557  m_cons1 = gamma_function(one + m_ds) * m_cs;
558  m_cons2 = gamma_function(one + m_dg) * m_cg;
559  m_cons3 = gamma_function(Real(4.0) + m_bs) / Real(6.0);
560  m_cons4 = gamma_function(Real(4.0) + m_br) / Real(6.0);
561  m_cons5 = gamma_function(one + m_bs);
562  m_cons6 = gamma_function(one + m_br);
563  m_cons7 = gamma_function(Real(4.0) + m_bg) / Real(6.0);
564  m_cons8 = gamma_function(one + m_bg);
565  m_cons9 = gamma_function(Real(5.0)/Real(2) + m_br/Real(2));
566  m_cons10 = gamma_function(Real(5.0)/Real(2) + m_bs/Real(2));
567  m_cons11 = gamma_function(Real(5.0)/Real(2) + m_bg/Real(2));
568  m_cons12 = gamma_function(one + m_di) * m_ci;
569  m_cons13 = gamma_function(m_bs + three) * m_pi / Real(4.0) * m_eci;
570  m_cons14 = gamma_function(m_bg + three) * m_pi / Real(4.0) * m_eci;
571  m_cons15 = -Real(1108.0) * m_eii * std::pow(m_pi, (one-m_bs)/three) *
572  std::pow(m_rhosn, (-Real(2)-m_bs)/three) / (Real(4.0)*Real(720.0));
573  m_cons16 = gamma_function(m_bi + three) * m_pi / Real(4.0) * m_eci;
574  m_cons17 = Real(4.0) * Real(2) * three * m_rhosu * m_pi * m_eci * m_eci *
575  gamma_function(Real(2)*m_bs + Real(2)) / (Real(8.0)*(m_rhog-m_rhosn));
576  m_cons18 = m_rhosn * m_rhosn;
577  m_cons19 = m_rhow * m_rhow;
578  m_cons20 = Real(20.0) * m_pi * m_pi * m_rhow * m_bimm;
579  m_cons21 = Real(4.0) / (m_dcs * m_rhoi);
580  m_cons22 = m_pi * m_rhoi * amrex::Math::powi<3>(m_dcs) / Real(6.0);
581  m_cons23 = m_pi / Real(4.0) * m_eii * gamma_function(m_bs + three);
582  m_cons24 = m_pi / Real(4.0) * m_ecr * gamma_function(m_br + three);
583  m_cons25 = m_pi * m_pi / Real(24.0) * m_rhow * m_ecr * gamma_function(m_br + Real(6.0));
584  m_cons26 = m_pi / Real(6.0) * m_rhow;
585  m_cons27 = gamma_function(one + m_bi);
586  m_cons28 = gamma_function(Real(4.0) + m_bi) / Real(6.0);
587  m_cons29 = Real(4.0)/three * m_pi * m_rhow * amrex::Math::powi<3>(Real(25.0E-6));
588  m_cons31 = m_pi * m_pi * m_ecr * m_rhosn;
589  m_cons32 = m_pi / Real(2) * m_ecr;
590  m_cons34 = Real(5.0)/Real(2) + m_br/Real(2);
591  m_cons35 = Real(5.0)/Real(2) + m_bs/Real(2);
592  m_cons36 = Real(5.0)/Real(2) + m_bg/Real(2);
593  m_cons37 = Real(4.0) * m_pi * Real(1.38E-23) / (Real(6.0) * m_pi * m_rin);
594  m_cons38 = m_pi * m_pi / three * m_rhow;
595  m_cons39 = m_pi * m_pi / Real(36.0) * m_rhow * m_bimm;
596  m_cons40 = m_pi / Real(6.0) * m_bimm;
597  m_cons41 = m_pi * m_pi * m_ecr * m_rhow;
598 
599  // Set CCN parameters for different environments
600  if (m_iact == 1) {
601  // Maritime CCN spectrum parameters (modified from Rasmussen et al. 2002)
602  // NCCN = C*S^K, where S is supersaturation in %
603  // m_k1 = Real(0.4); // Exponent in CCN activation formula
604  // m_c1 = Real(120.0); // Coefficient in CCN activation formula (cm^-3)
605  }
606 
607  // Initialize aerosol activation parameters for IACT=2
608  if (m_iact == 2) {
609  // Parameters for ammonium sulfate
610  // m_mw = Real(0.018); // Molecular weight of water (kg/mol)
611  // m_osm = one; // Osmotic coefficient
612  // m_vi = three; // Number of ions dissociated in solution
613  // m_epsm = Real(0.7); // Aerosol soluble fraction
614  // m_rhoa = Real(1777.0); // Aerosol bulk density (kg/m^3)
615  // m_map = Real(0.132); // Molecular weight of aerosol (kg/mol)
616  // m_ma = Real(0.0284); // Molecular weight of air (kg/mol)
617  // m_rr = Real(8.3145); // Universal gas constant (J/mol/K)
618  // m_bact = m_vi * m_osm * m_epsm * m_mw * m_rhoa / (m_map * m_rhow);
619 
620  // Aerosol size distribution parameters for MPACE (Morrison et al. 2007, JGR)
621  // Mode 1
622  // m_rm1 = Real(0.052E-6); // Geometric mean radius, mode 1 (m)
623  // m_sig1 = Real(2.04); // Standard deviation of aerosol size distribution, mode 1
624  m_nanew1 = Real(72.2E6); // Total aerosol concentration, mode 1 (m^-3)
625  // m_f11 = myhalf * std::exp(Real(2.5) * amrex::Math::powi<2>(std::log(m_sig1)));
626  // m_f21 = one + fourth * std::log(m_sig1);
627 
628  // Mode 2
629  // m_rm2 = Real(1.3E-6); // Geometric mean radius, mode 2 (m)
630  // m_sig2 = Real(2.5); // Standard deviation of aerosol size distribution, mode 2
631  m_nanew2 = Real(1.8E6); // Total aerosol concentration, mode 2 (m^-3)
632  // m_f12 = myhalf * std::exp(Real(2.5) * amrex::Math::powi<2>(std::log(m_sig2)));
633  // m_f22 = one + fourth * std::log(m_sig2);
634  }
635  // Set microphysics control parameters
636  m_iact = 2; // Lognormal aerosol activation
637  m_inuc = 0; // Mid-latitude ice nucleation (Cooper)
638  if (sc.moisture_type == MoistureType::Morrison_NoIce) {
639  m_iliq = 1; // Include ice processes
640  m_igraup = 1; // Include graupel processes
641  } else {
642  m_iliq = 0; // Include ice processes
643  m_igraup = 0; // Include graupel processes
644  }
645  m_ihail = 0; // Use graupel (0) instead of hail (1)
646  // m_isub = 0; // Sub-grid vertical velocity option
647  // m_do_radar_ref = false; // Disable radar reflectivity by default
648  Box boxD(box); boxD.makeSlab(2,0);
649 
650  if(run_morr_cpp) {
651 
652  // One FAB to rule them all
653  FArrayBox morr_fab(grown_box, MORRInd::NumInds, Arena_Used);
654  morr_fab.template setVal<RunOn::Device>(0);
655  auto const& morr_arr = morr_fab.array();
656 
657  ////////////////////////////////////////////////////////////
658  // ParallelFor for testing partial C++ implementation
659  // NOTE: Currently all Array4 values are copied to locals
660  // This means we're not updating or outputting anything
661  ////////////////////////////////////////////////////////////
662  ParallelFor( box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
663  {
664  // Tendencies and mixing ratios
665  morr_arr(i,j,k,MORRInd::qc3d) = qcl_arr(i,j,k); // CLOUD WATER MIXING RATIO
666  morr_arr(i,j,k,MORRInd::qi3d) = qci_arr(i,j,k); // CLOUD ICE MIXING RATIO
667  morr_arr(i,j,k,MORRInd::qni3d) = qps_arr(i,j,k); // SNOW MIXING RATIO
668  morr_arr(i,j,k,MORRInd::qr3d) = qpr_arr(i,j,k); // RAIN MIXING RATIO
669  morr_arr(i,j,k,MORRInd::qg3d) = qpg_arr(i,j,k); // GRAUPEL MIXING RATIO
670 
671  morr_arr(i,j,k,MORRInd::nc3d) = nc_arr(i,j,k); // CLOUD WATER NUMBER CONCENTRATION
672  morr_arr(i,j,k,MORRInd::ni3d) = ni_arr(i,j,k); // CLOUD ICE NUMBER CONCENTRATION
673  morr_arr(i,j,k,MORRInd::ns3d) = ns_arr(i,j,k); // SNOW NUMBER CONCENTRATION
674  morr_arr(i,j,k,MORRInd::nr3d) = nr_arr(i,j,k); // RAIN NUMBER CONCENTRATION
675  morr_arr(i,j,k,MORRInd::ng3d) = ng_arr(i,j,k); // GRAUPEL NUMBER CONCENTRATION
676 
677  morr_arr(i,j,k,MORRInd::t3d) = theta_arr(i,j,k) * pii_arr(i,j,k); // TEMPERATURE
678  morr_arr(i,j,k,MORRInd::qv3d) = qv_arr(i,j,k); // WATER VAPOR MIXING RATIO
679  morr_arr(i,j,k,MORRInd::pres) = pres_arr(i,j,k); // ATMOSPHERIC PRESSURE
680  morr_arr(i,j,k,MORRInd::dzq) = dz_arr(i,j,k); // DIFFERENCE IN HEIGHT ACROSS LEVEL
681  morr_arr(i,j,k,MORRInd::w3d) = w_arr(i,j,k); // GRID-SCALE VERTICAL VELOCITY
682 
683  // NOTE: There are no cumulus tendencies passed to Morrison
684  // and the FORTRAN version zeros these out.
685  morr_arr(i,j,k,MORRInd::qrcu1d) = Real(0); //morr_arr(i,j,k,MORRInd::qrcuten_arr); // RAIN FROM CUMULUS PARAMETERIZATION
686  morr_arr(i,j,k,MORRInd::qscu1d) = Real(0); //morr_arr(i,j,k,MORRInd::qscuten_arr); // SNOW FROM CUMULUS PARAMETERIZATION
687  morr_arr(i,j,k,MORRInd::qicu1d) = Real(0); //morr_arr(i,j,k,MORRInd::qicuten_arr); // ICE FROM CUMULUS PARAMETERIZATION
688  });
689 
690  ParallelFor( boxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
691  {
692  int ltrue=0; // LTRUE: SWITCH = 0: NO HYDROMETEORS IN COLUMN, = 1: HYDROMETEORS IN COLUMN
693  int nstep; // NSTEP: Timestep counter
694  int iinum=m_inum; // iinum: Integer control variable
695 
696  for (int k=klo; k<=khi; k++) {
697  // Microphysical processes
698  // Real nsubc; // NSUBC: Loss of NC during evaporation
699  Real nsubi; // NSUBI: Loss of NI during sublimation
700  Real nsubs; // NSUBS: Loss of NS during sublimation
701  Real nsubr; // NSUBR: Loss of NR during evaporation
702  Real prd; // PRD: Deposition cloud ice
703  Real pre; // PRE: Evaporation of rain
704  Real prds; // PRDS: Deposition snow
705  Real nnuccc; // NNUCCC: Change N due to contact freezing droplets
706  Real mnuccc; // MNUCCC: Change Q due to contact freezing droplets
707  Real pra; // PRA: Accretion droplets by rain
708  Real prc; // PRC: Autoconversion droplets
709  Real pcc; // PCC: Condensation/evaporation droplets
710  Real nnuccd; // NNUCCD: Change N freezing aerosol (primary ice nucleation)
711  Real mnuccd; // MNUCCD: Change Q freezing aerosol (primary ice nucleation)
712  Real mnuccr; // MNUCCR: Change Q due to contact freezing rain
713  Real nnuccr; // NNUCCR: Change N due to contact freezing rain
714  Real npra; // NPRA: Change N due to droplet accretion by rain
715  Real nragg; // NRAGG: Self-collection/breakup of rain
716  Real nsagg; // NSAGG: Self-collection of snow
717  Real nprc; // NPRC: Change NC autoconversion droplets
718  Real nprc1; // NPRC1: Change NR autoconversion droplets
719  Real prai; // PRAI: Change Q accretion cloud ice by snow
720  Real prci; // PRCI: Change Q autoconversion cloud ice to snow
721  Real psacws; // PSACWS: Change Q droplet accretion by snow
722  Real npsacws; // NPSACWS: Change N droplet accretion by snow
723  Real psacwi; // PSACWI: Change Q droplet accretion by cloud ice
724  Real npsacwi; // NPSACWI: Change N droplet accretion by cloud ice
725  Real nprci; // NPRCI: Change N autoconversion cloud ice by snow
726  Real nprai; // NPRAI: Change N accretion cloud ice
727  Real nmults; // NMULTS: Ice multiplication due to riming droplets by snow
728  Real nmultr; // NMULTR: Ice multiplication due to riming rain by snow
729  Real qmults; // QMULTS: Change Q due to ice multiplication droplets/snow
730  Real qmultr; // QMULTR: Change Q due to ice multiplication rain/snow
731  Real pracs; // PRACS: Change Q rain-snow collection
732  Real npracs; // NPRACS: Change N rain-snow collection
733  // Real pccn; // PCCN: Change Q droplet activation
734  Real psmlt; // PSMLT: Change Q melting snow to rain
735  Real evpms; // EVPMS: Change Q melting snow evaporating
736  Real nsmlts; // NSMLTS: Change N melting snow
737  Real nsmltr; // NSMLTR: Change N melting snow to rain
738  Real piacr; // PIACR: Change QR, ice-rain collection
739  Real niacr; // NIACR: Change N, ice-rain collection
740  Real praci; // PRACI: Change QI, ice-rain collection
741  Real piacrs; // PIACRS: Change QR, ice rain collision, added to snow
742  Real niacrs; // NIACRS: Change N, ice rain collision, added to snow
743  Real pracis; // PRACIS: Change QI, ice rain collision, added to snow
744  Real eprd; // EPRD: Sublimation cloud ice
745  Real eprds; // EPRDS: Sublimation snow
746 
747  // Graupel processes
748  Real pracg; // PRACG: Change in Q collection rain by graupel
749  Real psacwg; // PSACWG: Change in Q collection droplets by graupel
750  Real pgsacw; // PGSACW: Conversion Q to graupel due to collection droplets by snow
751  Real pgracs; // PGRACS: Conversion Q to graupel due to collection rain by snow
752  Real prdg; // PRDG: Deposition of graupel
753  Real eprdg; // EPRDG: Sublimation of graupel
754  Real evpmg; // EVPMG: Change Q melting of graupel and evaporation
755  Real pgmlt; // PGMLT: Change Q melting of graupel
756  Real npracg; // NPRACG: Change N collection rain by graupel
757  Real npsacwg; // NPSACWG: Change N collection droplets by graupel
758  Real nscng; // NSCNG: Change N conversion to graupel due to collection droplets by snow
759  Real ngracs; // NGRACS: Change N conversion to graupel due to collection rain by snow
760  Real ngmltg; // NGMLTG: Change N melting graupel
761  Real ngmltr; // NGMLTR: Change N melting graupel to rain
762  Real nsubg; // NSUBG: Change N sublimation/deposition of graupel
763  Real psacr; // PSACR: Conversion due to collection of snow by rain
764  Real nmultg; // NMULTG: Ice multiplication due to accretion droplets by graupel
765  Real nmultrg; // NMULTRG: Ice multiplication due to accretion rain by graupel
766  Real qmultg; // QMULTG: Change Q due to ice multiplication droplets/graupel
767  Real qmultrg; // QMULTRG: Change Q due to ice multiplication rain/graupel
768 
769  // Time-varying atmospheric parameters
770  Real kap; // KAP: Thermal conductivity of air
771  Real evs; // EVS: Saturation vapor pressure
772  Real eis; // EIS: Ice saturation vapor pressure
773  Real qvs; // QVS: Saturation mixing ratio
774  Real qvi; // QVI: Ice saturation mixing ratio
775  Real qvqvs; // QVQVS: Saturation ratio
776  Real qvqvsi; // QVQVSI: Ice saturation ratio
777  Real dv; // DV: Diffusivity of water vapor in air
778  Real sc_schmidt; // SC: Schmidt number
779  Real ab; // AB: Correction to condensation rate due to latent heating
780  Real abi; // ABI: Correction to deposition rate due to latent heating
781 
782  // Dummy variables
783  Real dum; // DUM: General dummy variable
784  Real dum1; // DUM1: General dummy variable
785  Real dumt; // DUMT: Dummy variable for temperature
786  Real dumqv; // DUMQV: Dummy variable for water vapor
787  Real dumqss; // DUMQSS: Dummy saturation mixing ratio
788  Real dums; // DUMS: General dummy variable
789 
790  // Prognostic supersaturation
791  Real dqsdt; // DQSDT: Change of saturation mixing ratio with temperature
792  Real dqsidt; // DQSIDT: Change in ice saturation mixing ratio with temperature
793 
794  Real epsi; // EPSI: 1/phase relaxation time (see M2005), ice
795  Real epss; // EPSS: 1/phase relaxation time (see M2005), snow
796  Real epsr; // EPSR: 1/phase relaxation time (see M2005), rain
797  Real epsg; // EPSG: 1/phase relaxation time (see M2005), graupel
798  Real kc2; // KC2: Total ice nucleation rate
799  Real di0; // DC0: Characteristic diameter for ice
800  Real ds0; // DS0: Characteristic diameter for snow
801  Real dg0; // DG0: Characteristic diameter for graupel
802  Real dumqc; // DUMQC: Dummy variable for cloud water mixing ratio
803  Real ratio; // RATIO: General ratio variable
804  Real sum_dep; // SUM_DEP: Sum of deposition/sublimation
805  Real fudgef; // FUDGEF: Adjustment factor
806 
807  // Real dum2; // DUM2: General dummy variable
808  // Real dumqsi; // DUMQSI: Dummy ice saturation mixing ratio
809  // Real dc0; // DC0: Characteristic diameter for cloud droplets
810  // Real dumqr; // DUMQR: Dummy variable for rain mixing ratio
811 
812  // For WRF-CHEM
813  // Real c2prec; // C2PREC: Cloud to precipitation conversion
814  // Real csed; // CSED: Cloud sedimentation
815  // Real ised; // ISED: Ice sedimentation
816  // Real ssed; // SSED: Snow sedimentation
817  // Real gsed; // GSED: Graupel sedimentation
818  // Real rsed; // RSED: Rain sedimentation
819  // Real tqimelt; // tqimelt: Melting of cloud ice (tendency)
820 
821  // NC3DTEN LOCAL ARRAY INITIALIZED
822  morr_arr(i,j,k,MORRInd::nc3dten) = Real(0);
823 
824  // INITIALIZE VARIABLES FOR WRF-CHEM OUTPUT TO ZERO
825  // c2prec = Real(0);
826  // csed = Real(0);
827  // ised = Real(0);
828  // ssed = Real(0);
829  // gsed = Real(0);
830  // rsed = Real(0);
831 
832  // LATENT HEAT OF VAPORIZATION
833  morr_arr(i,j,k,MORRInd::xxlv) = Real(3.1484E6) - Real(2370.0) * morr_arr(i,j,k,MORRInd::t3d);
834  // LATENT HEAT OF SUBLIMATION
835  morr_arr(i,j,k,MORRInd::xxls) = Real(3.15E6) - Real(2370.0) * morr_arr(i,j,k,MORRInd::t3d) + Real(0.3337E6);
836 
837  // Assuming CP is a constant defined elsewhere (specific heat of dry air at constant pressure)
838  const Real CP = Real(1004.5); // J/kg/K
839  morr_arr(i,j,k,MORRInd::cpm) = CP * (one + Real(0.887) * morr_arr(i,j,k,MORRInd::qv3d));
840 
841  // SATURATION VAPOR PRESSURE AND MIXING RATIO
842  // hm, add fix for low pressure, 5/12/10
843  // Assuming POLYSVP is defined elsewhere
844  evs = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(morr_arr(i,j,k,MORRInd::t3d), 0)); // PA
845  eis = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(morr_arr(i,j,k,MORRInd::t3d), 1)); // PA
846  // MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
847  if (eis > evs) {
848  eis = evs; // temporary update: adjust ice saturation pressure
849  }
850 
851  // SATURATION MIXING RATIOS
852  qvs = m_ep_2 * evs / (morr_arr(i,j,k,MORRInd::pres) - evs); // budget equation: calculate water saturation mixing ratio
853  qvi = m_ep_2 * eis / (morr_arr(i,j,k,MORRInd::pres) - eis); // budget equation: calculate ice saturation mixing ratio
854 
855  // SATURATION RATIOS
856  qvqvs = morr_arr(i,j,k,MORRInd::qv3d) / qvs; // budget equation: calculate water saturation ratio
857  qvqvsi = morr_arr(i,j,k,MORRInd::qv3d) / qvi; // budget equation: calculate ice saturation ratio
858 
859  // AIR DENSITY
860  morr_arr(i,j,k,MORRInd::rho) = morr_arr(i,j,k,MORRInd::pres) / (m_R * morr_arr(i,j,k,MORRInd::t3d)); // budget equation: calculate air density
861 
862  ds0 = three; // Size distribution parameter for snow
863  di0 = three; // Size distribution parameter for cloud ice
864  dg0 = three; // Size distribution parameter for graupel
865 
866  // ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
867  // ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
868  // ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
869  // FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
870  if (morr_arr(i,j,k,MORRInd::qrcu1d) >= Real(1.0e-10)) {
871  dum = Real(1.8e5) * std::pow(morr_arr(i,j,k,MORRInd::qrcu1d) * dt / (m_pi * m_rhow * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::rho))), fourth); // rate equation: calculate rain number concentration from cumulus
872  morr_arr(i,j,k,MORRInd::nr3d) += dum; // budget equation: update rain number concentration
873  }
874  if (morr_arr(i,j,k,MORRInd::qscu1d) >= Real(1.0e-10)) {
875  dum = Real(3.e5) * std::pow(morr_arr(i,j,k,MORRInd::qscu1d) * dt / (m_cons1 * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::rho))), one / (ds0 + one)); // rate equation: calculate snow number concentration from cumulus
876  morr_arr(i,j,k,MORRInd::ns3d) += dum; // budget equation: update snow number concentration
877  }
878  if (morr_arr(i,j,k,MORRInd::qicu1d) >= Real(1.0e-10)) {
879  dum = morr_arr(i,j,k,MORRInd::qicu1d) * dt / (m_ci * std::pow(Real(80.0e-6), di0)); // rate equation: calculate cloud ice number concentration from cumulus
880  morr_arr(i,j,k,MORRInd::ni3d) += dum; // budget equation: update cloud ice number concentration
881  }
882 
883  // AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
884  // hm modify 7/0/09 change limit to Real(1.e-8)
885  if (qvqvs < Real(0.9)) {
886  if (morr_arr(i,j,k,MORRInd::qr3d) < Real(1.0e-8)) {
887  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qr3d); // budget equation: transfer rain to vapor
888  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: adjust temperature
889  morr_arr(i,j,k,MORRInd::qr3d) = Real(0); // temporary update: set rain to Real(0)
890  }
891  if (morr_arr(i,j,k,MORRInd::qc3d) < Real(1.0e-8)) {
892  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qc3d); // budget equation: transfer cloud water to vapor
893  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: adjust temperature
894  morr_arr(i,j,k,MORRInd::qc3d) = Real(0); // temporary update: set cloud water to Real(0)
895  }
896  }
897  if (qvqvsi < Real(0.9)) {
898  if (morr_arr(i,j,k,MORRInd::qi3d) < Real(1.0e-8)) {
899  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qi3d); // budget equation: transfer cloud ice to vapor
900  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: adjust temperature
901  morr_arr(i,j,k,MORRInd::qi3d) = Real(0); // temporary update: set cloud ice to Real(0)
902  }
903  if (morr_arr(i,j,k,MORRInd::qni3d) < Real(1.0e-8)) {
904  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qni3d); // budget equation: transfer snow to vapor
905  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qni3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: adjust temperature
906  morr_arr(i,j,k,MORRInd::qni3d) = Real(0); // temporary update: set snow to Real(0)
907  }
908  if (morr_arr(i,j,k,MORRInd::qg3d) < Real(1.0e-8)) {
909  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qg3d); // budget equation: transfer graupel to vapor
910  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qg3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: adjust temperature
911  morr_arr(i,j,k,MORRInd::qg3d) = Real(0); // temporary update: set graupel to Real(0)
912  }
913  }
914  // HEAT OF FUSION
915  morr_arr(i,j,k,MORRInd::xlf) = morr_arr(i,j,k,MORRInd::xxls) - morr_arr(i,j,k,MORRInd::xxlv);
916 
917  // IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
918  // Note: QSMALL is not defined in the variable list, so I'll define it
919  const Real QSMALL = m_qsmall;
920 
921  if (morr_arr(i,j,k,MORRInd::qc3d) < QSMALL) {
922  morr_arr(i,j,k,MORRInd::qc3d) = Real(0);
923  morr_arr(i,j,k,MORRInd::nc3d) = Real(0);
924  morr_arr(i,j,k,MORRInd::effc) = Real(0);
925  }
926  if (morr_arr(i,j,k,MORRInd::qr3d) < QSMALL) {
927  morr_arr(i,j,k,MORRInd::qr3d) = Real(0);
928  morr_arr(i,j,k,MORRInd::nr3d) = Real(0);
929  morr_arr(i,j,k,MORRInd::effr) = Real(0);
930  }
931  if (morr_arr(i,j,k,MORRInd::qi3d) < QSMALL) {
932  morr_arr(i,j,k,MORRInd::qi3d) = Real(0);
933  morr_arr(i,j,k,MORRInd::ni3d) = Real(0);
934  morr_arr(i,j,k,MORRInd::effi) = Real(0);
935  }
936  if (morr_arr(i,j,k,MORRInd::qni3d) < QSMALL) {
937  morr_arr(i,j,k,MORRInd::qni3d) = Real(0);
938  morr_arr(i,j,k,MORRInd::ns3d) = Real(0);
939  morr_arr(i,j,k,MORRInd::effs) = Real(0);
940  }
941  if (morr_arr(i,j,k,MORRInd::qg3d) < QSMALL) {
942  morr_arr(i,j,k,MORRInd::qg3d) = Real(0);
943  morr_arr(i,j,k,MORRInd::ng3d) = Real(0);
944  morr_arr(i,j,k,MORRInd::effg) = Real(0);
945  }
946  // INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
947  morr_arr(i,j,k,MORRInd::qrsten) = Real(0); // temporary update: initialize QRSTEN
948  morr_arr(i,j,k,MORRInd::qisten) = Real(0); // temporary update: initialize QISTEN
949  morr_arr(i,j,k,MORRInd::qnisten) = Real(0); // temporary update: initialize QNISTEN
950  morr_arr(i,j,k,MORRInd::qcsten) = Real(0); // temporary update: initialize QCSTEN
951  morr_arr(i,j,k,MORRInd::qgsten) = Real(0); // temporary update: initialize QGSTEN
952 
953  // MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
954  morr_arr(i,j,k,MORRInd::mu) = Real(1.496e-6) * std::pow(morr_arr(i,j,k,MORRInd::t3d), Real(1.5)) / (morr_arr(i,j,k,MORRInd::t3d) + Real(120.0)); // budget equation: calculate air viscosity
955 
956  // Fall speed with density correction (Heymsfield and Benssemer 2006)
957  dum = std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.54)); // temporary update: calculate density correction factor
958 
959  // AA revision 4/1/11: Ikawa and Saito 1991 air-density correction
960  morr_arr(i,j,k,MORRInd::ain) = std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.35)) * m_ai; // budget equation: calculate ice fall speed parameter
961  morr_arr(i,j,k,MORRInd::arn) = dum * m_ar; // budget equation: calculate rain fall speed parameter
962  morr_arr(i,j,k,MORRInd::asn) = dum * m_as; // budget equation: calculate snow fall speed parameter
963 
964  // AA revision 4/1/11: temperature-dependent Stokes fall speed
965  morr_arr(i,j,k,MORRInd::acn) = m_g * m_rhow / (Real(18.0) * morr_arr(i,j,k,MORRInd::mu)); // budget equation: calculate cloud droplet fall speed parameter
966 
967  // HM ADD GRAUPEL 8/28/06
968  morr_arr(i,j,k,MORRInd::agn) = dum * m_ag; // budget equation: calculate graupel fall speed parameter
969  // hm 4/7/09 bug fix, initialize morr_arr(i,j,k,MORRInd::lami) to prevent later division by Real(0)
970  morr_arr(i,j,k,MORRInd::lami) = Real(0); // temporary update: initialize LAMI
971 
972  // If there is no cloud/precip water, and if subsaturated, then skip microphysics for this level
973  bool skipMicrophysics = false;
974  bool skipConcentrations = false;
975  if (morr_arr(i,j,k,MORRInd::qc3d) < QSMALL && morr_arr(i,j,k,MORRInd::qi3d) < QSMALL && morr_arr(i,j,k,MORRInd::qni3d) < QSMALL && morr_arr(i,j,k,MORRInd::qr3d) < QSMALL && morr_arr(i,j,k,MORRInd::qg3d) < QSMALL) {
976  if ((morr_arr(i,j,k,MORRInd::t3d) < Real(273.15) && qvqvsi < Real(0.999)) || (morr_arr(i,j,k,MORRInd::t3d) >= Real(273.15) && qvqvs < Real(0.999))) {
977  skipMicrophysics = true;// goto label_200;
978  }
979  }
980 
981  if(!skipMicrophysics) {
982 
983  // Thermal conductivity for air
984  kap = Real(1.414e3) * morr_arr(i,j,k,MORRInd::mu); // budget equation: calculate thermal conductivity
985 
986  // Diffusivity of water vapor
987  dv = Real(8.794e-5) * std::pow(morr_arr(i,j,k,MORRInd::t3d), Real(1.81)) / morr_arr(i,j,k,MORRInd::pres); // budget equation: calculate vapor diffusivity
988 
989  // Schmidt number
990  sc_schmidt = morr_arr(i,j,k,MORRInd::mu) / (morr_arr(i,j,k,MORRInd::rho) * dv); // budget equation: calculate Schmidt number
991 
992  // Psychometric corrections
993  // Rate of change sat. mix. ratio with temperature
994  dum = (m_Rv * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::t3d))); // temporary update: calculate temperature factor
995  dqsdt = morr_arr(i,j,k,MORRInd::xxlv) * qvs / dum; // budget equation: calculate DQSDT
996  dqsidt = morr_arr(i,j,k,MORRInd::xxls) * qvi / dum; // budget equation: calculate DQSIDT
997  abi = one + dqsidt * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: calculate ABI
998  ab = one + dqsdt * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm); // budget equation: calculate AB
999 
1000  // CASE FOR TEMPERATURE ABOVE FREEZING
1001  if (morr_arr(i,j,k,MORRInd::t3d) >= Real(273.15)) {
1002  //......................................................................
1003  // ALLOW FOR CONSTANT DROPLET NUMBER
1004  // INUM = 0, PREDICT DROPLET NUMBER
1005  // INUM = 1, SET CONSTANT DROPLET NUMBER
1006 
1007  if (m_inum == 1) {
1008  // CONVERT NDCNST FROM CM-3 TO KG-1
1009  // Note: NDCNST constant would need to be defined elsewhere
1010  morr_arr(i,j,k,MORRInd::nc3d) = m_ndcnst * Real(1.0e6) / morr_arr(i,j,k,MORRInd::rho); // Set cloud droplet number concentration
1011  }
1012 
1013  // GET SIZE DISTRIBUTION PARAMETERS
1014  // MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1015  if (morr_arr(i,j,k,MORRInd::qni3d) < Real(1.0e-6)) {
1016  morr_arr(i,j,k,MORRInd::qr3d) = morr_arr(i,j,k,MORRInd::qr3d) + morr_arr(i,j,k,MORRInd::qni3d); // Transfer snow to rain
1017  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::nr3d) + morr_arr(i,j,k,MORRInd::ns3d); // Transfer snow number to rain
1018  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) - morr_arr(i,j,k,MORRInd::qni3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm); // Adjust temperature
1019  morr_arr(i,j,k,MORRInd::qni3d) = Real(0); // Set snow to Real(0)
1020  morr_arr(i,j,k,MORRInd::ns3d) = Real(0); // Set snow number to Real(0)
1021  }
1022 
1023  if (morr_arr(i,j,k,MORRInd::qg3d) < Real(1.0e-6)) {
1024  morr_arr(i,j,k,MORRInd::qr3d) = morr_arr(i,j,k,MORRInd::qr3d) + morr_arr(i,j,k,MORRInd::qg3d); // Transfer graupel to rain
1025  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::nr3d) + morr_arr(i,j,k,MORRInd::ng3d); // Transfer graupel number to rain
1026  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) - morr_arr(i,j,k,MORRInd::qg3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm); // Adjust temperature
1027  morr_arr(i,j,k,MORRInd::qg3d) = Real(0); // Set graupel to Real(0)
1028  morr_arr(i,j,k,MORRInd::ng3d) = Real(0); // Set graupel number to Real(0)
1029  }
1030  // Skip to label 300 if concentrations are below thresholds
1031  if (morr_arr(i,j,k,MORRInd::qc3d) < m_qsmall && morr_arr(i,j,k,MORRInd::qni3d) < Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qr3d) < m_qsmall && morr_arr(i,j,k,MORRInd::qg3d) < Real(1.0e-8)) {
1032  skipConcentrations=true;// goto label_300;
1033  }
1034  if(!skipConcentrations) {
1035  morr_arr(i,j,k,MORRInd::ns3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::ns3d));
1036  morr_arr(i,j,k,MORRInd::nc3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::nc3d));
1037  morr_arr(i,j,k,MORRInd::nr3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::nr3d));
1038  morr_arr(i,j,k,MORRInd::ng3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::ng3d));
1039 
1040  // ========================================================================
1041  // USING WRF APPROACH FOR SIZE DISTRIBUTION PARAMETERS
1042  // ========================================================================
1043  // Rain
1044  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
1045  // Calculate lambda parameter using cons26 (pi*rhow/6)
1046  morr_arr(i,j,k,MORRInd::lamr) = std::pow(m_pi * m_rhow * morr_arr(i,j,k,MORRInd::nr3d) / morr_arr(i,j,k,MORRInd::qr3d), one/three);
1047  morr_arr(i,j,k,MORRInd::n0r) = morr_arr(i,j,k,MORRInd::nr3d)*morr_arr(i,j,k,MORRInd::lamr);
1048 
1049  // Check for slope and adjust vars
1050  if (morr_arr(i,j,k,MORRInd::lamr) < m_lamminr) {
1051  morr_arr(i,j,k,MORRInd::lamr) = m_lamminr;
1052  morr_arr(i,j,k,MORRInd::n0r) = std::pow(morr_arr(i,j,k,MORRInd::lamr), Real(4.0)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
1053  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr); // Update number concentration
1054  } else if (morr_arr(i,j,k,MORRInd::lamr) > m_lammaxr) {
1055  morr_arr(i,j,k,MORRInd::lamr) = m_lammaxr;
1056  morr_arr(i,j,k,MORRInd::n0r) = std::pow(morr_arr(i,j,k,MORRInd::lamr), Real(4.0)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
1057  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr); // Update number concentration
1058  }
1059  }
1060 
1061  // Cloud droplets
1062  if (morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1063  // Calculate air density factor (moist air density)
1064  dum = morr_arr(i,j,k,MORRInd::pres)/(Real(287.15)*morr_arr(i,j,k,MORRInd::t3d));
1065 
1066  // MARTIN ET AL. (1994) FORMULA FOR PGAM (WRF implementation)
1067  morr_arr(i,j,k,MORRInd::pgam) = Real(0.0005714)*(morr_arr(i,j,k,MORRInd::nc3d)/Real(1.0e6)*dum) + Real(0.2714);
1068  morr_arr(i,j,k,MORRInd::pgam) = one/(morr_arr(i,j,k,MORRInd::pgam)*morr_arr(i,j,k,MORRInd::pgam)) - one;
1069  morr_arr(i,j,k,MORRInd::pgam) = amrex::max(morr_arr(i,j,k,MORRInd::pgam), Real(2));
1070  morr_arr(i,j,k,MORRInd::pgam) = amrex::min(morr_arr(i,j,k,MORRInd::pgam), Real(10.0));
1071 
1072  // Calculate gamma function values
1073  Real gamma_pgam_plus_1 = gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one);
1074  Real gamma_pgam_plus_4 = gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0));
1075 
1076  // Calculate lambda parameter
1077  morr_arr(i,j,k,MORRInd::lamc) = std::pow((m_cons26 * morr_arr(i,j,k,MORRInd::nc3d) * gamma_pgam_plus_4) / (morr_arr(i,j,k,MORRInd::qc3d) * gamma_pgam_plus_1), one/three);
1078 
1079  // Lambda bounds from WRF - 60 micron max diameter, 1 micron min diameter
1080  Real lambda_min = (morr_arr(i,j,k,MORRInd::pgam) + one)/Real(60.0e-6);
1081  Real lambda_max = (morr_arr(i,j,k,MORRInd::pgam) + one)/Real(1.0e-6);
1082 
1083  // Check bounds and update number concentration if needed
1084  if (morr_arr(i,j,k,MORRInd::lamc) < lambda_min) {
1085  morr_arr(i,j,k,MORRInd::lamc) = lambda_min;
1086  // Update cloud droplet number using the same formula as in WRF
1087  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three*std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
1088  std::log(gamma_pgam_plus_1) - std::log(gamma_pgam_plus_4))/ m_cons26;
1089  } else if (morr_arr(i,j,k,MORRInd::lamc) > lambda_max) {
1090  morr_arr(i,j,k,MORRInd::lamc) = lambda_max;
1091  // Update cloud droplet number using the same formula as in WRF
1092  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three*std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
1093  std::log(gamma_pgam_plus_1) - std::log(gamma_pgam_plus_4))/ m_cons26;
1094  }
1095 
1096  // Calculate intercept parameter
1097  morr_arr(i,j,k,MORRInd::cdist1) = morr_arr(i,j,k,MORRInd::nc3d) / gamma_pgam_plus_1;
1098  }
1099 
1100  // Snow
1101  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
1102  // Calculate lambda parameter
1103  morr_arr(i,j,k,MORRInd::lams) = std::pow(m_cons1 * morr_arr(i,j,k,MORRInd::ns3d) / morr_arr(i,j,k,MORRInd::qni3d), one/ds0);
1104 
1105  // Calculate intercept parameter
1106  morr_arr(i,j,k,MORRInd::n0s) = morr_arr(i,j,k,MORRInd::ns3d) * morr_arr(i,j,k,MORRInd::lams);
1107 
1108  // Check for slope and adjust vars
1109  if (morr_arr(i,j,k,MORRInd::lams) < m_lammins) {
1110  morr_arr(i,j,k,MORRInd::lams) = m_lammins;
1111  morr_arr(i,j,k,MORRInd::n0s) = std::pow(morr_arr(i,j,k,MORRInd::lams), Real(4.0)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
1112  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams); // Update number concentration
1113  } else if (morr_arr(i,j,k,MORRInd::lams) > m_lammaxs) {
1114  morr_arr(i,j,k,MORRInd::lams) = m_lammaxs;
1115  morr_arr(i,j,k,MORRInd::n0s) = std::pow(morr_arr(i,j,k,MORRInd::lams), Real(4.0)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
1116  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams); // Update number concentration
1117  }
1118  }
1119 
1120  // Graupel
1121  if (morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
1122  // Calculate lambda parameter
1123  morr_arr(i,j,k,MORRInd::lamg) = std::pow(m_cons2 * morr_arr(i,j,k,MORRInd::ng3d) / morr_arr(i,j,k,MORRInd::qg3d), one/dg0);
1124 
1125  // Calculate intercept parameter
1126  morr_arr(i,j,k,MORRInd::n0g) = morr_arr(i,j,k,MORRInd::ng3d) * morr_arr(i,j,k,MORRInd::lamg);
1127 
1128  // Check for slope and adjust vars
1129  if (morr_arr(i,j,k,MORRInd::lamg) < m_lamming) {
1130  morr_arr(i,j,k,MORRInd::lamg) = m_lamming;
1131  morr_arr(i,j,k,MORRInd::n0g) = std::pow(morr_arr(i,j,k,MORRInd::lamg), Real(4.0)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
1132  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg); // Update number concentration
1133  } else if (morr_arr(i,j,k,MORRInd::lamg) > m_lammaxg) {
1134  morr_arr(i,j,k,MORRInd::lamg) = m_lammaxg;
1135  morr_arr(i,j,k,MORRInd::n0g) = std::pow(morr_arr(i,j,k,MORRInd::lamg), Real(4.0)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
1136  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg); // Update number concentration
1137  }
1138  }
1139  ////////////////////// First instance of ZERO OUT PROCESS RATES
1140  // Zero out process rates
1141  prc = Real(0); // Cloud water to rain conversion rate (PRC)
1142  nprc = Real(0); // Change in cloud droplet number due to autoconversion (NPRC)
1143  nprc1 = Real(0); // Change in rain number due to autoconversion (NPRC1)
1144  pra = Real(0); // Accretion of cloud water by rain (PRA)
1145  npra = Real(0); // Change in cloud droplet number due to accretion by rain (NPRA)
1146  nragg = Real(0); // Self-collection/breakup of rain (NRAGG)
1147  nsmlts = Real(0); // Loss of snow number during melting (NSMLTS)
1148  nsmltr = Real(0); // Change in rain number due to snow melting (NSMLTR)
1149  evpms = Real(0); // Melting snow evaporation rate (EVPMS)
1150  pcc = Real(0); // Condensation/evaporation of cloud water (PCC)
1151  pre = Real(0); // Evaporation of rain (PRE)
1152  // nsubc = Real(0); // Loss of cloud droplet number during evaporation (NSUBC)
1153  nsubr = Real(0); // Loss of rain number during evaporation (NSUBR)
1154  pracg = Real(0); // Collection of rain by graupel (PRACG)
1155  npracg = Real(0); // Change in number due to collection of rain by graupel (NPRACG)
1156  psmlt = Real(0); // Melting of snow (PSMLT)
1157  pgmlt = Real(0); // Melting of graupel (PGMLT)
1158  evpmg = Real(0); // Evaporation of melting graupel (EVPMG)
1159  pracs = Real(0); // Collection of snow by rain (PRACS)
1160  npracs = Real(0); // Change in number due to collection of snow by rain (NPRACS)
1161  ngmltg = Real(0); // Loss of graupel number during melting (NGMLTG)
1162  ngmltr = Real(0); // Change in rain number due to graupel melting (NGMLTR)
1163 
1164  // CALCULATION OF MICROPHYSICAL PROCESS RATES, T > Real(273.15) K
1165 
1166  // AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1167  // FORMULA FROM BEHENG (1994)
1168  // USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1169  // AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1170  // AS A GAMMA DISTRIBUTION
1171 
1172  // USE MINIMUM VALUE OF Real(1.E-6) TO PREVENT FLOATING POINT ERROR
1173 
1174  if (morr_arr(i,j,k,MORRInd::qc3d) >= Real(1.0e-6)) {
1175  // HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1176  // FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1177  prc = Real(1350.0) * std::pow(morr_arr(i,j,k,MORRInd::qc3d), Real(2.47)) *
1178  std::pow((morr_arr(i,j,k,MORRInd::nc3d)/Real(1.0e6)*morr_arr(i,j,k,MORRInd::rho)), -Real(1.79));
1179 
1180  // note: nprc1 is change in Nr,
1181  // nprc is change in Nc
1182  nprc1 = prc / m_cons29;
1183  nprc = prc / (morr_arr(i,j,k,MORRInd::qc3d) / morr_arr(i,j,k,MORRInd::nc3d));
1184 
1185  // hm bug fix 3/20/12
1186  nprc = std::min(nprc, morr_arr(i,j,k,MORRInd::nc3d) / dt);
1187  nprc1 = std::min(nprc1, nprc);
1188  }
1189 
1190  // HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1191  // FORMULA FROM IKAWA AND SAITO (1991)
1192 
1193  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8)) {
1194  Real ums_local = morr_arr(i,j,k,MORRInd::asn) * m_cons3 / std::pow(morr_arr(i,j,k,MORRInd::lams), m_bs);
1195  Real umr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons4 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1196  Real uns_local = morr_arr(i,j,k,MORRInd::asn) * m_cons5 / std::pow(morr_arr(i,j,k,MORRInd::lams), m_bs);
1197  Real unr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons6 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1198 
1199  // SET REALISTIC LIMITS ON FALLSPEEDS
1200  // bug fix, 10/08/09
1201  dum = std::pow(m_rhosu/morr_arr(i,j,k,MORRInd::rho), Real(0.54));
1202  ums_local = std::min(ums_local, Real(1.2)*dum);
1203  uns_local = std::min(uns_local, Real(1.2)*dum);
1204  umr_local = std::min(umr_local, Real(9.1)*dum);
1205  unr_local = std::min(unr_local, Real(9.1)*dum);
1206 
1207 
1208  // hm fix, 2/12/13
1209  // for above freezing conditions to get accelerated melting of snow,
1210  // we need collection of rain by snow (following Lin et al. 1983)
1211  ////////////////////////Might needstd::pow expanding
1212  pracs = m_cons41 * (std::sqrt(amrex::Math::powi<2>(Real(1.2)*umr_local-Real(0.95)*ums_local) +
1213  Real(0.08)*ums_local*umr_local) * morr_arr(i,j,k,MORRInd::rho) *
1214  morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0s) / amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) *
1215  (Real(5.0)/(amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lams)) +
1216  Real(2)/(amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lams))) +
1217  myhalf/(morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lams)))));
1218  }
1219  // ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
1220  // ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
1221  // ASSUME SHED DROPS ARE 1 MM IN SIZE
1222 
1223  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qg3d) >= Real(1.0e-8)) {
1224 
1225  Real umg_local = morr_arr(i,j,k,MORRInd::agn) * m_cons7 / std::pow(morr_arr(i,j,k,MORRInd::lamg), m_bg);
1226  Real umr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons4 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1227  Real ung_local = morr_arr(i,j,k,MORRInd::agn) * m_cons8 / std::pow(morr_arr(i,j,k,MORRInd::lamg), m_bg);
1228  Real unr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons6 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1229 
1230  // SET REALISTIC LIMITS ON FALLSPEEDS
1231  // bug fix, 10/08/09
1232  dum = std::pow(m_rhosu/morr_arr(i,j,k,MORRInd::rho), Real(0.54));
1233  umg_local = std::min(umg_local, Real(20.0)*dum);
1234  ung_local = std::min(ung_local, Real(20.0)*dum);
1235  umr_local = std::min(umr_local, Real(9.1)*dum);
1236  unr_local = std::min(unr_local, Real(9.1)*dum);
1237 
1238  // PRACG IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
1239  pracg = m_cons41 * (std::sqrt(amrex::Math::powi<2>(Real(1.2)*umr_local-Real(0.95)*umg_local) +
1240  Real(0.08)*umg_local*umr_local) * morr_arr(i,j,k,MORRInd::rho) *
1241  morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0g) / amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) *
1242  (Real(5.0)/(amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lamg)) +
1243  Real(2)/(amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamg))) +
1244  myhalf/(morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamg)))));
1245 
1246  // ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1247  dum = pracg/Real(5.2e-7);
1248 
1249  npracg = m_cons32 * morr_arr(i,j,k,MORRInd::rho) * (std::sqrt(Real(1.7)*amrex::Math::powi<2>(unr_local-ung_local) +
1250  Real(0.3)*unr_local*ung_local) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0g) *
1251  (one/(amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lamg)) +
1252  one/(amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamg))) +
1253  one/(morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamg)))));
1254  // hm 7/15/13, remove limit so that the number of collected drops can smaller than
1255  // number of shed drops
1256  npracg = npracg - dum;
1257  }
1258  // ACCRETION OF CLOUD LIQUID WATER BY RAIN
1259  // CONTINUOUS COLLECTION EQUATION WITH
1260  // GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1261 
1262  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qc3d) >= Real(1.0e-8)) {
1263  // 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1264  // KHAIROUTDINOV AND KOGAN 2000, MWR
1265  dum = morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::qr3d);
1266  pra = Real(67.0) * std::pow(dum, Real(1.15));
1267  npra = pra / (morr_arr(i,j,k,MORRInd::qc3d) / morr_arr(i,j,k,MORRInd::nc3d));
1268  }
1269 
1270  // SELF-COLLECTION OF RAIN DROPS
1271  // FROM BEHENG(1994)
1272  // FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
1273  // AS DESCRIBED ABOVE FOR AUTOCONVERSION
1274 
1275  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8)) {
1276  // include breakup add 10/09/09
1277  dum1 = Real(300.0e-6);
1278  if (one/morr_arr(i,j,k,MORRInd::lamr) < dum1) {
1279  dum = one;
1280  } else {
1281  dum = Real(2) - std::exp(Real(2300.0) * (one/morr_arr(i,j,k,MORRInd::lamr) - dum1));
1282  }
1283  nragg = -Real(5.78) * dum * morr_arr(i,j,k,MORRInd::nr3d) * morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::rho);
1284  }
1285  // CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1286  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
1287  epsr = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::rho) * dv *
1288  (m_f1r/(morr_arr(i,j,k,MORRInd::lamr)*morr_arr(i,j,k,MORRInd::lamr)) +
1289  m_f2r * std::sqrt(morr_arr(i,j,k,MORRInd::arn)*morr_arr(i,j,k,MORRInd::rho)/morr_arr(i,j,k,MORRInd::mu)) *
1290  std::pow(sc_schmidt, one/three) * m_cons9 /
1291  std::pow(morr_arr(i,j,k,MORRInd::lamr), m_cons34));
1292  } else {
1293  epsr = Real(0);
1294  }
1295  // NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1296  if (morr_arr(i,j,k,MORRInd::qv3d) < qvs) {
1297  pre = epsr * (morr_arr(i,j,k,MORRInd::qv3d) - qvs) / ab;
1298  pre = std::min(pre, Real(0));
1299  } else {
1300  pre = Real(0);
1301  }
1302  // MELTING OF SNOW
1303  // SNOW MAY PERSIST ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1304  // IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1305 
1306  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8)) {
1307  // fix 053011
1308  // HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1309  dum = -m_cpw/morr_arr(i,j,k,MORRInd::xlf) * (morr_arr(i,j,k,MORRInd::t3d) - Real(273.15)) * pracs;
1310 
1311  // hm fix 1/20/15
1312  psmlt = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0s) * kap * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d)) /
1313  morr_arr(i,j,k,MORRInd::xlf) * (m_f1s/(morr_arr(i,j,k,MORRInd::lams)*morr_arr(i,j,k,MORRInd::lams)) +
1314  m_f2s * std::sqrt(morr_arr(i,j,k,MORRInd::asn)*morr_arr(i,j,k,MORRInd::rho)/morr_arr(i,j,k,MORRInd::mu)) *
1315  std::pow(sc_schmidt, one/three) * m_cons10 /
1316  std::pow(morr_arr(i,j,k,MORRInd::lams), m_cons35)) + dum;
1317 
1318  // IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1319  if (qvqvs < one) {
1320  epss = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0s) * morr_arr(i,j,k,MORRInd::rho) * dv *
1321  (m_f1s/(morr_arr(i,j,k,MORRInd::lams)*morr_arr(i,j,k,MORRInd::lams)) +
1322  m_f2s * std::sqrt(morr_arr(i,j,k,MORRInd::asn)*morr_arr(i,j,k,MORRInd::rho)/morr_arr(i,j,k,MORRInd::mu)) *
1323  std::pow(sc_schmidt, one/three) * m_cons10 /
1324  std::pow(morr_arr(i,j,k,MORRInd::lams), m_cons35));
1325 
1326  // hm fix 8/4/08
1327  evpms = (morr_arr(i,j,k,MORRInd::qv3d) - qvs) * epss / ab;
1328  evpms = std::max(evpms, psmlt);
1329  psmlt = psmlt - evpms;
1330  }
1331  }
1332  // MELTING OF GRAUPEL
1333  // GRAUPEL MAY PERSIST ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1334  // IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
1335 
1336  if (morr_arr(i,j,k,MORRInd::qg3d) >= Real(1.0e-8)) {
1337  // fix 053011
1338  // HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1339 
1340  dum = -m_cpw/morr_arr(i,j,k,MORRInd::xlf) * (morr_arr(i,j,k,MORRInd::t3d) - Real(273.15)) * pracg;
1341 
1342  // hm fix 1/20/15
1343  pgmlt = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0g) * kap * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d)) /
1344  morr_arr(i,j,k,MORRInd::xlf) * (m_f1s/(morr_arr(i,j,k,MORRInd::lamg)*morr_arr(i,j,k,MORRInd::lamg)) +
1345  m_f2s * std::sqrt(morr_arr(i,j,k,MORRInd::agn)*morr_arr(i,j,k,MORRInd::rho)/morr_arr(i,j,k,MORRInd::mu)) *
1346  std::pow(sc_schmidt, one/three) * m_cons11 /
1347  std::pow(morr_arr(i,j,k,MORRInd::lamg), m_cons36)) + dum;
1348 
1349  // IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
1350  if (qvqvs < one) {
1351  epsg = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0g) * morr_arr(i,j,k,MORRInd::rho) * dv *
1352  (m_f1s/(morr_arr(i,j,k,MORRInd::lamg)*morr_arr(i,j,k,MORRInd::lamg)) +
1353  m_f2s * std::sqrt(morr_arr(i,j,k,MORRInd::agn)*morr_arr(i,j,k,MORRInd::rho)/morr_arr(i,j,k,MORRInd::mu)) *
1354  std::pow(sc_schmidt, one/three) * m_cons11 /
1355  std::pow(morr_arr(i,j,k,MORRInd::lamg), m_cons36));
1356 
1357  // hm fix 8/4/08
1358  evpmg = (morr_arr(i,j,k,MORRInd::qv3d) - qvs) * epsg / ab;
1359  evpmg = std::max(evpmg, pgmlt);
1360  pgmlt = pgmlt - evpmg;
1361  }
1362  }
1363  // HM, V3.2
1364  // RESET PRACG AND PRACS TO ZERO, THIS IS DONE BECAUSE THERE IS NO
1365  // TRANSFER OF MASS FROM SNOW AND GRAUPEL TO RAIN DIRECTLY FROM COLLECTION
1366  // ABOVE FREEZING, IT IS ONLY USED FOR ENHANCEMENT OF MELTING AND SHEDDING
1367 
1368  pracg = Real(0);
1369  pracs = Real(0);
1370  // CONSERVATION OF QC
1371  dum = (prc + pra) * dt;
1372 
1373  if (dum > morr_arr(i,j,k,MORRInd::qc3d) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1374  ratio = morr_arr(i,j,k,MORRInd::qc3d) / dum;
1375  prc = prc * ratio;
1376  pra = pra * ratio;
1377  }
1378 
1379  // CONSERVATION OF SNOW
1380  dum = (-psmlt - evpms + pracs) * dt;
1381 
1382  if (dum > morr_arr(i,j,k,MORRInd::qni3d) && morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
1383  // NO SOURCE TERMS FOR SNOW AT T > FREEZING
1384  ratio = morr_arr(i,j,k,MORRInd::qni3d) / dum;
1385  psmlt = psmlt * ratio;
1386  evpms = evpms * ratio;
1387  pracs = pracs * ratio;
1388  }
1389 
1390  // CONSERVATION OF GRAUPEL
1391  dum = (-pgmlt - evpmg + pracg) * dt;
1392 
1393  if (dum > morr_arr(i,j,k,MORRInd::qg3d) && morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
1394  // NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
1395  ratio = morr_arr(i,j,k,MORRInd::qg3d) / dum;
1396  pgmlt = pgmlt * ratio;
1397  evpmg = evpmg * ratio;
1398  pracg = pracg * ratio;
1399  }
1400 
1401  // CONSERVATION OF QR
1402  // HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
1403 
1404  dum = (-pracs - pracg - pre - pra - prc + psmlt + pgmlt) * dt;
1405 
1406  if (dum > morr_arr(i,j,k,MORRInd::qr3d) && morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
1407  ratio = (morr_arr(i,j,k,MORRInd::qr3d)/dt + pracs + pracg + pra + prc - psmlt - pgmlt) / (-pre);
1408  pre = pre * ratio;
1409  }
1410  // Update tendencies
1411  morr_arr(i,j,k,MORRInd::qv3dten) = morr_arr(i,j,k,MORRInd::qv3dten) + (-pre - evpms - evpmg);
1412 
1413  morr_arr(i,j,k,MORRInd::t3dten) = morr_arr(i,j,k,MORRInd::t3dten) + (pre * morr_arr(i,j,k,MORRInd::xxlv) +
1414  (evpms + evpmg) * morr_arr(i,j,k,MORRInd::xxls) +
1415  (psmlt + pgmlt - pracs - pracg) * morr_arr(i,j,k,MORRInd::xlf)) / morr_arr(i,j,k,MORRInd::cpm);
1416 
1417  morr_arr(i,j,k,MORRInd::qc3dten) = morr_arr(i,j,k,MORRInd::qc3dten) + (-pra - prc);
1418  morr_arr(i,j,k,MORRInd::qr3dten) = morr_arr(i,j,k,MORRInd::qr3dten) + (pre + pra + prc - psmlt - pgmlt + pracs + pracg);
1419  morr_arr(i,j,k,MORRInd::qni3dten) = morr_arr(i,j,k,MORRInd::qni3dten) + (psmlt + evpms - pracs);
1420  morr_arr(i,j,k,MORRInd::qg3dten) = morr_arr(i,j,k,MORRInd::qg3dten) + (pgmlt + evpmg - pracg);
1421 
1422  // fix 053011
1423  // HM, bug fix 5/12/08, npracg is subtracted from nr not ng
1424  morr_arr(i,j,k,MORRInd::nc3dten) = morr_arr(i,j,k,MORRInd::nc3dten) + (-npra - nprc);
1425  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) + (nprc1 + nragg - npracg);
1426 
1427  // HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
1428  // c2prec = pra + prc;
1429 
1430  if (pre < Real(0)) {
1431  dum = pre * dt / morr_arr(i,j,k,MORRInd::qr3d);
1432  dum = std::max(-one, dum);
1433  nsubr = dum * morr_arr(i,j,k,MORRInd::nr3d) / dt;
1434  }
1435 
1436  if (evpms + psmlt < Real(0)) {
1437  dum = (evpms + psmlt) * dt / morr_arr(i,j,k,MORRInd::qni3d);
1438  dum = std::max(-one, dum);
1439  nsmlts = dum * morr_arr(i,j,k,MORRInd::ns3d) / dt;
1440  }
1441 
1442  if (psmlt < Real(0)) {
1443  dum = psmlt * dt / morr_arr(i,j,k,MORRInd::qni3d);
1444  dum = std::max(-one, dum);
1445  nsmltr = dum * morr_arr(i,j,k,MORRInd::ns3d) / dt;
1446  }
1447 
1448  if (evpmg + pgmlt < Real(0)) {
1449  dum = (evpmg + pgmlt) * dt / morr_arr(i,j,k,MORRInd::qg3d);
1450  dum = std::max(-one, dum);
1451  ngmltg = dum * morr_arr(i,j,k,MORRInd::ng3d) / dt;
1452  }
1453 
1454  if (pgmlt < Real(0)) {
1455  dum = pgmlt * dt / morr_arr(i,j,k,MORRInd::qg3d);
1456  dum = std::max(-one, dum);
1457  ngmltr = dum * morr_arr(i,j,k,MORRInd::ng3d) / dt;
1458  }
1459 
1460  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + nsmlts;
1461  morr_arr(i,j,k,MORRInd::ng3dten) = morr_arr(i,j,k,MORRInd::ng3dten) + ngmltg;
1462  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) + (nsubr - nsmltr - ngmltr);
1463 
1464  }
1465  //Right after 300 CONTINUE
1466 // label_300:
1467  // Calculate saturation adjustment to condense extra vapor above water saturation
1468  dumt = morr_arr(i,j,k,MORRInd::t3d) + dt * morr_arr(i,j,k,MORRInd::t3dten);
1469  dumqv = morr_arr(i,j,k,MORRInd::qv3d) + dt * morr_arr(i,j,k,MORRInd::qv3dten);
1470 
1471  // Fix for low pressure (added 5/12/10)
1472  dum = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(dumt, 0));
1473  dumqss = m_ep_2 * dum / (morr_arr(i,j,k,MORRInd::pres) - dum);
1474  dumqc = morr_arr(i,j,k,MORRInd::qc3d) + dt * morr_arr(i,j,k,MORRInd::qc3dten);
1475  dumqc = std::max(dumqc, Real(0));
1476 
1477  // Saturation adjustment for liquid
1478  dums = dumqv - dumqss;
1479  pcc = dums / (one + amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::xxlv)) * dumqss / (morr_arr(i,j,k,MORRInd::cpm) * m_Rv * amrex::Math::powi<2>(dumt))) / dt;
1480  if (pcc * dt + dumqc < Real(0)) {
1481  pcc = -dumqc / dt;
1482  }
1483 
1484  if (!do_cond) { pcc = Real(0); }
1485 
1486  // Update tendencies
1487  morr_arr(i,j,k,MORRInd::qv3dten) -= pcc;
1488  morr_arr(i,j,k,MORRInd::t3dten) += pcc * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm);
1489  morr_arr(i,j,k,MORRInd::qc3dten) += pcc;
1490  } else { //cold
1491  //......................................................................
1492  // ALLOW FOR CONSTANT DROPLET NUMBER
1493  // INUM = 0, PREDICT DROPLET NUMBER
1494  // INUM = 1, SET CONSTANT DROPLET NUMBER
1495 
1496  if (m_inum == 1) {
1497  // CONVERT NDCNST FROM CM-3 TO KG-1
1498  // Note: NDCNST constant would need to be defined elsewhere
1499  morr_arr(i,j,k,MORRInd::nc3d) = m_ndcnst * Real(1.0e6) / morr_arr(i,j,k,MORRInd::rho); // Set cloud droplet number concentration
1500  }
1501 
1502  morr_arr(i,j,k,MORRInd::ni3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::ni3d));
1503  morr_arr(i,j,k,MORRInd::ns3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::ns3d));
1504  morr_arr(i,j,k,MORRInd::nc3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::nc3d));
1505  morr_arr(i,j,k,MORRInd::nr3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::nr3d));
1506  morr_arr(i,j,k,MORRInd::ng3d) = amrex::max(Real(0),morr_arr(i,j,k,MORRInd::ng3d));
1507 
1508  // ========================================================================
1509  // USING WRF APPROACH FOR SIZE DISTRIBUTION PARAMETERS
1510  // ========================================================================
1511  // Rain
1512  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
1513  // Calculate lambda parameter using cons26 (pi*rhow/6)
1514  morr_arr(i,j,k,MORRInd::lamr) = std::pow(m_pi * m_rhow * morr_arr(i,j,k,MORRInd::nr3d) / morr_arr(i,j,k,MORRInd::qr3d), one/three);
1515 
1516  // Check for slope and adjust vars
1517  if (morr_arr(i,j,k,MORRInd::lamr) < m_lamminr) {
1518  morr_arr(i,j,k,MORRInd::lamr) = m_lamminr;
1519  morr_arr(i,j,k,MORRInd::n0r) = std::pow(morr_arr(i,j,k,MORRInd::lamr), Real(4.0)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
1520  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr); // Update number concentration
1521  } else if (morr_arr(i,j,k,MORRInd::lamr) > m_lammaxr) {
1522  morr_arr(i,j,k,MORRInd::lamr) = m_lammaxr;
1523  morr_arr(i,j,k,MORRInd::n0r) = std::pow(morr_arr(i,j,k,MORRInd::lamr), Real(4.0)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
1524  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr); // Update number concentration
1525  } else {
1526  // Calculate intercept parameter using WRF formula
1527  morr_arr(i,j,k,MORRInd::n0r) = std::pow(morr_arr(i,j,k,MORRInd::lamr), Real(4.0)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
1528  }
1529  }
1530 
1531 
1532  // Cloud droplets
1533  if (morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1534  // Calculate air density factor (moist air density)
1535  dum = morr_arr(i,j,k,MORRInd::pres)/(Real(287.15)*morr_arr(i,j,k,MORRInd::t3d));
1536 
1537  // MARTIN ET AL. (1994) FORMULA FOR PGAM (WRF implementation)
1538  morr_arr(i,j,k,MORRInd::pgam) = Real(0.0005714)*(morr_arr(i,j,k,MORRInd::nc3d)/Real(1.0e6)*dum) + Real(0.2714);
1539  morr_arr(i,j,k,MORRInd::pgam) = one/(amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::pgam))) - one;
1540  morr_arr(i,j,k,MORRInd::pgam) = amrex::max(morr_arr(i,j,k,MORRInd::pgam), Real(2));
1541  morr_arr(i,j,k,MORRInd::pgam) = amrex::min(morr_arr(i,j,k,MORRInd::pgam), Real(10.0));
1542 
1543  // Calculate gamma function values
1544  Real gamma_pgam_plus_1 = gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one);
1545  Real gamma_pgam_plus_4 = gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0));
1546 
1547  // Calculate lambda parameter
1548  morr_arr(i,j,k,MORRInd::lamc) = std::pow((m_cons26 * morr_arr(i,j,k,MORRInd::nc3d) * gamma_pgam_plus_4) / (morr_arr(i,j,k,MORRInd::qc3d) * gamma_pgam_plus_1), one/three);
1549 
1550  // Lambda bounds from WRF - 60 micron max diameter, 1 micron min diameter
1551  Real lambda_min = (morr_arr(i,j,k,MORRInd::pgam) + one)/Real(60.0e-6);
1552  Real lambda_max = (morr_arr(i,j,k,MORRInd::pgam) + one)/Real(1.0e-6);
1553 
1554  // Check bounds and update number concentration if needed
1555  if (morr_arr(i,j,k,MORRInd::lamc) < lambda_min) {
1556  morr_arr(i,j,k,MORRInd::lamc) = lambda_min;
1557  // Update cloud droplet number using the same formula as in WRF
1558  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three*std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
1559  std::log(gamma_pgam_plus_1) - std::log(gamma_pgam_plus_4))/ m_cons26;
1560  } else if (morr_arr(i,j,k,MORRInd::lamc) > lambda_max) {
1561  morr_arr(i,j,k,MORRInd::lamc) = lambda_max;
1562  // Update cloud droplet number using the same formula as in WRF
1563  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three*std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
1564  std::log(gamma_pgam_plus_1) - std::log(gamma_pgam_plus_4))/ m_cons26;
1565  }
1566 
1567  // Calculate intercept parameter
1568  morr_arr(i,j,k,MORRInd::cdist1) = morr_arr(i,j,k,MORRInd::nc3d) / gamma_pgam_plus_1;
1569  }
1570 
1571  // Snow
1572  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
1573  // Calculate lambda parameter
1574  morr_arr(i,j,k,MORRInd::lams) = std::pow(m_cons1 * morr_arr(i,j,k,MORRInd::ns3d) / morr_arr(i,j,k,MORRInd::qni3d), one/ds0);
1575 
1576  // Calculate intercept parameter
1577  morr_arr(i,j,k,MORRInd::n0s) = morr_arr(i,j,k,MORRInd::ns3d) * morr_arr(i,j,k,MORRInd::lams);
1578 
1579  // Check for slope and adjust vars
1580  if (morr_arr(i,j,k,MORRInd::lams) < m_lammins) {
1581  morr_arr(i,j,k,MORRInd::lams) = m_lammins;
1582  morr_arr(i,j,k,MORRInd::n0s) = std::pow(morr_arr(i,j,k,MORRInd::lams), Real(4.0)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
1583  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams); // Update number concentration
1584  } else if (morr_arr(i,j,k,MORRInd::lams) > m_lammaxs) {
1585  morr_arr(i,j,k,MORRInd::lams) = m_lammaxs;
1586  morr_arr(i,j,k,MORRInd::n0s) = std::pow(morr_arr(i,j,k,MORRInd::lams), Real(4.0)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
1587  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams); // Update number concentration
1588  }
1589  }
1590 
1591  // Cloud ice
1592  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
1593  // Calculate lambda parameter
1594  morr_arr(i,j,k,MORRInd::lami) = std::pow(m_cons12 * morr_arr(i,j,k,MORRInd::ni3d) / morr_arr(i,j,k,MORRInd::qi3d), one/three);
1595 
1596  // Calculate intercept parameter (initial calculation)
1597  morr_arr(i,j,k,MORRInd::n0i) = morr_arr(i,j,k,MORRInd::ni3d) * morr_arr(i,j,k,MORRInd::lami);
1598 
1599  // Check for slope (apply bounds)
1600  if (morr_arr(i,j,k,MORRInd::lami) < m_lammini) {
1601  morr_arr(i,j,k,MORRInd::lami) = m_lammini;
1602  // Recalculate morr_arr(i,j,k,MORRInd::n0i) when lambda is adjusted
1603  morr_arr(i,j,k,MORRInd::n0i) = std::pow(morr_arr(i,j,k,MORRInd::lami), Real(4.0)) * morr_arr(i,j,k,MORRInd::qi3d) / m_cons12;
1604  // Update ni3d when lambda is adjusted
1605  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::n0i) / morr_arr(i,j,k,MORRInd::lami);
1606  } else if (morr_arr(i,j,k,MORRInd::lami) > m_lammaxi) {
1607  morr_arr(i,j,k,MORRInd::lami) = m_lammaxi;
1608  // Recalculate morr_arr(i,j,k,MORRInd::n0i) when lambda is adjusted
1609  morr_arr(i,j,k,MORRInd::n0i) = std::pow(morr_arr(i,j,k,MORRInd::lami), Real(4.0)) * morr_arr(i,j,k,MORRInd::qi3d) / m_cons12;
1610  // Update ni3d when lambda is adjusted
1611  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::n0i) / morr_arr(i,j,k,MORRInd::lami);
1612  }
1613  }
1614  // Graupel
1615  if (morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
1616  // Calculate lambda parameter
1617  morr_arr(i,j,k,MORRInd::lamg) = std::pow(m_cons2 * morr_arr(i,j,k,MORRInd::ng3d) / morr_arr(i,j,k,MORRInd::qg3d), one/dg0);
1618 
1619  // Calculate intercept parameter
1620  morr_arr(i,j,k,MORRInd::n0g) = morr_arr(i,j,k,MORRInd::ng3d) * morr_arr(i,j,k,MORRInd::lamg);
1621 
1622  // Check for slope and adjust vars
1623  if (morr_arr(i,j,k,MORRInd::lamg) < m_lamming) {
1624  morr_arr(i,j,k,MORRInd::lamg) = m_lamming;
1625  morr_arr(i,j,k,MORRInd::n0g) = std::pow(morr_arr(i,j,k,MORRInd::lamg), Real(4.0)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
1626  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg); // Update number concentration
1627  } else if (morr_arr(i,j,k,MORRInd::lamg) > m_lammaxg) {
1628  morr_arr(i,j,k,MORRInd::lamg) = m_lammaxg;
1629  morr_arr(i,j,k,MORRInd::n0g) = std::pow(morr_arr(i,j,k,MORRInd::lamg), Real(4.0)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
1630  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg); // Update number concentration
1631  }
1632  }
1633  ////////////////////// Second instance of ZERO OUT PROCESS RATES
1634  // Zero out process rates
1635  mnuccc = Real(0); // Change Q due to contact freezing droplets (MNUCCC)
1636  nnuccc = Real(0); // Change N due to contact freezing droplets (NNUCCC)
1637  prc = Real(0); // Autoconversion droplets (PRC)
1638  nprc = Real(0); // Change NC autoconversion droplets (NPRC)
1639  nprc1 = Real(0); // Change NR autoconversion droplets (NPRC1)
1640  nsagg = Real(0); // Self-collection of snow (NSAGG)
1641  psacws = Real(0); // Change Q droplet accretion by snow (PSACWS)
1642  npsacws = Real(0); // Change N droplet accretion by snow (NPSACWS)
1643  psacwi = Real(0); // Change Q droplet accretion by cloud ice (PSACWI)
1644  npsacwi = Real(0); // Change N droplet accretion by cloud ice (NPSACWI)
1645  pracs = Real(0); // Change Q rain-snow collection (PRACS)
1646  npracs = Real(0); // Change N rain-snow collection (NPRACS)
1647  nmults = Real(0); // Ice multiplication due to riming droplets by snow (NMULTS)
1648  qmults = Real(0); // Change Q due to ice multiplication droplets/snow (QMULTS)
1649  nmultr = Real(0); // Ice multiplication due to riming rain by snow (NMULTR)
1650  qmultr = Real(0); // Change Q due to ice multiplication rain/snow (QMULTR)
1651  nmultg = Real(0); // Ice multiplication due to accretion droplets by graupel (NMULTG)
1652  qmultg = Real(0); // Change Q due to ice multiplication droplets/graupel (QMULTG)
1653  nmultrg = Real(0); // Ice multiplication due to accretion rain by graupel (NMULTRG)
1654  qmultrg = Real(0); // Change Q due to ice multiplication rain/graupel (QMULTRG)
1655  mnuccr = Real(0); // Change Q due to contact freezing rain (MNUCCR)
1656  nnuccr = Real(0); // Change N due to contact freezing rain (NNUCCR)
1657  pra = Real(0); // Accretion droplets by rain (PRA)
1658  npra = Real(0); // Change N due to droplet accretion by rain (NPRA)
1659  nragg = Real(0); // Self-collection/breakup of rain (NRAGG)
1660  prci = Real(0); // Change Q autoconversion cloud ice to snow (PRCI)
1661  nprci = Real(0); // Change N autoconversion cloud ice by snow (NPRCI)
1662  prai = Real(0); // Change Q accretion cloud ice by snow (PRAI)
1663  nprai = Real(0); // Change N accretion cloud ice (NPRAI)
1664  nnuccd = Real(0); // Change N freezing aerosol (primary ice nucleation) (NNUCCD)
1665  mnuccd = Real(0); // Change Q freezing aerosol (primary ice nucleation) (MNUCCD)
1666  pcc = Real(0); // Condensation/evaporation droplets (PCC)
1667  pre = Real(0); // Evaporation of rain (PRE)
1668  prd = Real(0); // Deposition cloud ice (PRD)
1669  prds = Real(0); // Deposition snow (PRDS)
1670  eprd = Real(0); // Sublimation cloud ice (EPRD)
1671  eprds = Real(0); // Sublimation snow (EPRDS)
1672  // nsubc = Real(0); // Loss of NC during evaporation (NSUBC)
1673  nsubi = Real(0); // Loss of NI during sublimation (NSUBI)
1674  nsubs = Real(0); // Loss of NS during sublimation (NSUBS)
1675  nsubr = Real(0); // Loss of NR during evaporation (NSUBR)
1676  piacr = Real(0); // Change QR, ice-rain collection (PIACR)
1677  niacr = Real(0); // Change N, ice-rain collection (NIACR)
1678  praci = Real(0); // Change QI, ice-rain collection (PRACI)
1679  piacrs = Real(0); // Change QR, ice rain collision, added to snow (PIACRS)
1680  niacrs = Real(0); // Change N, ice rain collision, added to snow (NIACRS)
1681  pracis = Real(0); // Change QI, ice rain collision, added to snow (PRACIS)
1682 
1683  // Graupel processes
1684  pracg = Real(0); // Change in Q collection rain by graupel (PRACG)
1685  psacr = Real(0); // Conversion due to collection of snow by rain (PSACR)
1686  psacwg = Real(0); // Change in Q collection droplets by graupel (PSACWG)
1687  pgsacw = Real(0); // Conversion Q to graupel due to collection droplets by snow (PGSACW)
1688  pgracs = Real(0); // Conversion Q to graupel due to collection rain by snow (PGRACS)
1689  prdg = Real(0); // Deposition of graupel (PRDG)
1690  eprdg = Real(0); // Sublimation of graupel (EPRDG)
1691  npracg = Real(0); // Change N collection rain by graupel (NPRACG)
1692  npsacwg = Real(0); // Change N collection droplets by graupel (NPSACWG)
1693  nscng = Real(0); // Change N conversion to graupel due to collection droplets by snow (NSCNG)
1694  ngracs = Real(0); // Change N conversion to graupel due to collection rain by snow (NGRACS)
1695  nsubg = Real(0); // Change N sublimation/deposition of graupel (NSUBG)
1696 
1697  ////////////////////// CALCULATION OF MICROPHYSICAL PROCESS RATES
1698  // FREEZING OF CLOUD DROPLETS - ONLY ALLOWED BELOW -4C
1699  if (morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall && morr_arr(i,j,k,MORRInd::t3d) < Real(269.15)) {
1700  // NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
1701  // FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
1702  // MEYERS CURVE
1703  Real nacnt = std::exp(-Real(2.80) + Real(0.262) * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) * Real(1000.0);
1704 
1705  // MEAN FREE PATH
1706  dum = Real(7.37) * morr_arr(i,j,k,MORRInd::t3d) / (Real(288.0) * Real(10.0) * morr_arr(i,j,k,MORRInd::pres)) / Real(100.0);
1707 
1708  // EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
1709  // BASED ON BROWNIAN DIFFUSION
1710  Real dap = m_cons37 * morr_arr(i,j,k,MORRInd::t3d) * (one + dum / m_rin) / morr_arr(i,j,k,MORRInd::mu);
1711 
1712  // CONTACT FREEZING
1713  mnuccc = m_cons38 * dap * nacnt * std::exp(std::log(morr_arr(i,j,k,MORRInd::cdist1)) +
1714  std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(5.0))) - Real(4.0) * std::log(morr_arr(i,j,k,MORRInd::lamc)));
1715  nnuccc = Real(2) * m_pi * dap * nacnt * morr_arr(i,j,k,MORRInd::cdist1) *
1716  gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(2)) / morr_arr(i,j,k,MORRInd::lamc);
1717 
1718  // IMMERSION FREEZING (BIGG 1953)
1719  // hm 7/15/13 fix for consistency w/ original formula
1720  mnuccc = mnuccc + m_cons39 *
1721  std::exp(std::log(morr_arr(i,j,k,MORRInd::cdist1)) + std::log(gamma_function(Real(7.0) + morr_arr(i,j,k,MORRInd::pgam))) - Real(6.0) * std::log(morr_arr(i,j,k,MORRInd::lamc))) *
1722  (std::exp(m_aimm * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) - one);
1723 
1724  nnuccc = nnuccc +
1725  m_cons40 * std::exp(std::log(morr_arr(i,j,k,MORRInd::cdist1)) + std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0))) - three * std::log(morr_arr(i,j,k,MORRInd::lamc))) *
1726  (std::exp(m_aimm * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) - one);
1727 
1728  // PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
1729  // MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
1730  nnuccc = std::min(nnuccc, morr_arr(i,j,k,MORRInd::nc3d) / dt);
1731  }
1732 
1733  // AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1734  // FORMULA FROM BEHENG (1994)
1735  // USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1736  // AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1737  // AS A GAMMA DISTRIBUTION
1738 
1739  // USE MINIMUM VALUE OF Real(1.E-6) TO PREVENT FLOATING POINT ERROR
1740  if (morr_arr(i,j,k,MORRInd::qc3d) >= Real(1.0e-6)) {
1741  // hm add 12/13/06, replace with newer formula
1742  // from khairoutdinov and kogan 2000, mwr
1743  prc = Real(1350.0) * std::pow(morr_arr(i,j,k,MORRInd::qc3d), Real(2.47)) *
1744  std::pow((morr_arr(i,j,k,MORRInd::nc3d) / Real(1.0e6) * morr_arr(i,j,k,MORRInd::rho)), -Real(1.79));
1745 
1746  // note: nprc1 is change in nr,
1747  // nprc is change in nc
1748  nprc1 = prc / m_cons29;
1749  nprc = prc / (morr_arr(i,j,k,MORRInd::qc3d) / morr_arr(i,j,k,MORRInd::nc3d));
1750 
1751  // hm bug fix 3/20/12
1752  nprc = std::min(nprc, morr_arr(i,j,k,MORRInd::nc3d) / dt);
1753  nprc1 = std::min(nprc1, nprc);
1754  }
1755  // SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
1756  // THIS IS HARD-WIRED FOR BS = Real(0.4) FOR NOW
1757  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8)) {
1758  nsagg = m_cons15 * morr_arr(i,j,k,MORRInd::asn) * std::pow(morr_arr(i,j,k,MORRInd::rho), ((Real(2) + m_bs) / three)) *
1759  std::pow(morr_arr(i,j,k,MORRInd::qni3d), ((Real(2) + m_bs) / three)) *
1760  std::pow((morr_arr(i,j,k,MORRInd::ns3d) * morr_arr(i,j,k,MORRInd::rho)), ((Real(4.0) - m_bs) / three)) / morr_arr(i,j,k,MORRInd::rho);
1761  }
1762 
1763  // ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
1764  // HERE USE CONTINUOUS COLLECTION EQUATION WITH
1765  // SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
1766 
1767  // SNOW
1768  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1769  psacws = m_cons13 * morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::rho) *
1770  morr_arr(i,j,k,MORRInd::n0s) / std::pow(morr_arr(i,j,k,MORRInd::lams), (m_bs + three));
1771 
1772  npsacws = m_cons13 * morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::nc3d) * morr_arr(i,j,k,MORRInd::rho) *
1773  morr_arr(i,j,k,MORRInd::n0s) / std::pow(morr_arr(i,j,k,MORRInd::lams), (m_bs + three));
1774  }
1775 
1776  // COLLECTION OF CLOUD WATER BY GRAUPEL
1777  if (morr_arr(i,j,k,MORRInd::qg3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1778  psacwg = m_cons14 * morr_arr(i,j,k,MORRInd::agn) * morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::rho) *
1779  morr_arr(i,j,k,MORRInd::n0g) / std::pow(morr_arr(i,j,k,MORRInd::lamg), (m_bg + three));
1780 
1781  npsacwg = m_cons14 * morr_arr(i,j,k,MORRInd::agn) * morr_arr(i,j,k,MORRInd::nc3d) * morr_arr(i,j,k,MORRInd::rho) *
1782  morr_arr(i,j,k,MORRInd::n0g) / std::pow(morr_arr(i,j,k,MORRInd::lamg), (m_bg + three));
1783  }
1784  // hm, add 12/13/06
1785  // CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
1786  // BEFORE RIMING CAN OCCUR
1787  // ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
1788  // TO HALLET-MOSSOP SPLINTERING
1789  if (morr_arr(i,j,k,MORRInd::qi3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
1790  // PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
1791  // FROM THOMPSON ET AL. 2004, MWR
1792  if (one / morr_arr(i,j,k,MORRInd::lami) >= Real(100.0e-6)) {
1793  psacwi = m_cons16 * morr_arr(i,j,k,MORRInd::ain) * morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::rho) *
1794  morr_arr(i,j,k,MORRInd::n0i) / std::pow(morr_arr(i,j,k,MORRInd::lami), (m_bi + three));
1795 
1796  npsacwi = m_cons16 * morr_arr(i,j,k,MORRInd::ain) * morr_arr(i,j,k,MORRInd::nc3d) * morr_arr(i,j,k,MORRInd::rho) *
1797  morr_arr(i,j,k,MORRInd::n0i) / std::pow(morr_arr(i,j,k,MORRInd::lami), (m_bi + three));
1798  }
1799  }
1800 
1801  // ACCRETION OF RAIN WATER BY SNOW
1802  // FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
1803  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8)) {
1804  Real ums_local = morr_arr(i,j,k,MORRInd::asn) * m_cons3 / std::pow(morr_arr(i,j,k,MORRInd::lams), m_bs);
1805  Real umr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons4 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1806  Real uns_local = morr_arr(i,j,k,MORRInd::asn) * m_cons5 / std::pow(morr_arr(i,j,k,MORRInd::lams), m_bs);
1807  Real unr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons6 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1808 
1809  // SET REASLISTIC LIMITS ON FALLSPEEDS
1810  // bug fix, 10/08/09
1811  dum = std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.54));
1812  ums_local = std::min(ums_local, Real(1.2) * dum);
1813  uns_local = std::min(uns_local, Real(1.2) * dum);
1814  umr_local = std::min(umr_local, Real(9.1) * dum);
1815  unr_local = std::min(unr_local, Real(9.1) * dum);
1816 
1817  pracs = m_cons41 * (std::sqrt(amrex::Math::powi<2>(Real(1.2) * umr_local - Real(0.95) * ums_local) +
1818  Real(0.08) * ums_local * umr_local) * morr_arr(i,j,k,MORRInd::rho) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0s) /
1819  amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * (Real(5.0) / (amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lams)) +
1820  Real(2) / (amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lams))) +
1821  myhalf / (morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lams)))));
1822 
1823  npracs = m_cons32 * morr_arr(i,j,k,MORRInd::rho) * std::sqrt(Real(1.7) * amrex::Math::powi<2>(unr_local - uns_local) +
1824  Real(0.3) * unr_local * uns_local) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0s) *
1825  (one / (amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lams)) +
1826  one / (amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lams))) +
1827  one / (morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lams))));
1828 
1829  // MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
1830  // AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
1831  // RIME-SPLINTERING
1832  pracs = std::min(pracs, morr_arr(i,j,k,MORRInd::qr3d) / dt);
1833 
1834  // COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
1835  // ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED Real(0.1) G/KG
1836  // hm modify for wrfv3.1
1837  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(0.1e-3) && morr_arr(i,j,k,MORRInd::qr3d) >= Real(0.1e-3)) {
1838  psacr = m_cons31 * (std::sqrt(amrex::Math::powi<2>(Real(1.2) * umr_local - Real(0.95) * ums_local) +
1839  Real(0.08) * ums_local * umr_local) * morr_arr(i,j,k,MORRInd::rho) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0s) /
1840  amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lams)) * (Real(5.0) / (amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lams)) * morr_arr(i,j,k,MORRInd::lamr)) +
1841  Real(2) / (amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lams)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr))) +
1842  myhalf / (morr_arr(i,j,k,MORRInd::lams) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)))));
1843  }
1844  }
1845 
1846  // COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
1847  // USED BY REISNER ET AL 1998
1848  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qg3d) >= Real(1.0e-8)) {
1849  Real umg_local = morr_arr(i,j,k,MORRInd::agn) * m_cons7 / std::pow(morr_arr(i,j,k,MORRInd::lamg), m_bg);
1850  Real umr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons4 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1851  Real ung_local = morr_arr(i,j,k,MORRInd::agn) * m_cons8 / std::pow(morr_arr(i,j,k,MORRInd::lamg), m_bg);
1852  Real unr_local = morr_arr(i,j,k,MORRInd::arn) * m_cons6 / std::pow(morr_arr(i,j,k,MORRInd::lamr), m_br);
1853 
1854  // SET REASLISTIC LIMITS ON FALLSPEEDS
1855  // bug fix, 10/08/09
1856  dum = std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.54));
1857  umg_local = std::min(umg_local, Real(20.0) * dum);
1858  ung_local = std::min(ung_local, Real(20.0) * dum);
1859  umr_local = std::min(umr_local, Real(9.1) * dum);
1860  unr_local = std::min(unr_local, Real(9.1) * dum);
1861 
1862  pracg = m_cons41 * (std::sqrt(amrex::Math::powi<2>(Real(1.2) * umr_local - Real(0.95) * umg_local) +
1863  Real(0.08) * umg_local * umr_local) * morr_arr(i,j,k,MORRInd::rho) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0g) /
1864  amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * (Real(5.0) / (amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lamg)) +
1865  Real(2) / (amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamg))) +
1866  myhalf / (morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamg)))));
1867 
1868  npracg = m_cons32 * morr_arr(i,j,k,MORRInd::rho) * std::sqrt(Real(1.7) * amrex::Math::powi<2>(unr_local - ung_local) +
1869  Real(0.3) * unr_local * ung_local) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::n0g) *
1870  (one / (amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::lamg)) +
1871  one / (amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::lamg))) +
1872  one / (morr_arr(i,j,k,MORRInd::lamr) * amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamg))));
1873 
1874  // MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
1875  // AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
1876  // RIME-SPLINTERING
1877  pracg = std::min(pracg, morr_arr(i,j,k,MORRInd::qr3d) / dt);
1878  }
1879 
1880  // RIME-SPLINTERING - SNOW
1881  // HALLET-MOSSOP (1974)
1882  // NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
1883  // hm add threshold snow and droplet mixing ratio for rime-splintering
1884  // to limit rime-splintering in stratiform clouds
1885  // these thresholds correspond with graupel thresholds in rh 1984
1886  //v1.4
1887  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(0.1e-3)) {
1888  if (morr_arr(i,j,k,MORRInd::qc3d) >= Real(0.5e-3) || morr_arr(i,j,k,MORRInd::qr3d) >= Real(0.1e-3)) {
1889  if (psacws > Real(0) || pracs > Real(0)) {
1890  if (morr_arr(i,j,k,MORRInd::t3d) < Real(270.16) && morr_arr(i,j,k,MORRInd::t3d) > Real(265.16)) {
1891  Real fmult = Real(0);
1892 
1893  if (morr_arr(i,j,k,MORRInd::t3d) > Real(270.16)) {
1894  fmult = Real(0);
1895  } else if (morr_arr(i,j,k,MORRInd::t3d) <= Real(270.16) && morr_arr(i,j,k,MORRInd::t3d) > Real(268.16)) {
1896  fmult = (Real(270.16) - morr_arr(i,j,k,MORRInd::t3d)) / Real(2);
1897  } else if (morr_arr(i,j,k,MORRInd::t3d) >= Real(265.16) && morr_arr(i,j,k,MORRInd::t3d) <= Real(268.16)) {
1898  fmult = (morr_arr(i,j,k,MORRInd::t3d) - Real(265.16)) / three;
1899  } else if (morr_arr(i,j,k,MORRInd::t3d) < Real(265.16)) {
1900  fmult = Real(0);
1901  }
1902 
1903  // 1000 IS TO CONVERT FROM KG TO G
1904  // SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
1905  if (psacws > Real(0)) {
1906  nmults = Real(35.0e4) * psacws * fmult * Real(1000.0);
1907  qmults = nmults * m_mmult;
1908 
1909  // CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
1910  // THAN WAS RIMED ONTO SNOW
1911  qmults = std::min(qmults, psacws);
1912  psacws = psacws - qmults;
1913  }
1914 
1915  // RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
1916  if (pracs > Real(0)) {
1917  nmultr = Real(35.0e4) * pracs * fmult * Real(1000.0);
1918  qmultr = nmultr * m_mmult;
1919 
1920  // CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
1921  // THAN WAS RIMED ONTO SNOW
1922  qmultr = std::min(qmultr, pracs);
1923  pracs = pracs - qmultr;
1924  }
1925  }
1926  }
1927  }
1928  }
1929 
1930  // RIME-SPLINTERING - GRAUPEL
1931  // HALLET-MOSSOP (1974)
1932  // NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
1933  // hm add threshold snow mixing ratio for rime-splintering
1934  // to limit rime-splintering in stratiform clouds
1935  // v1.4
1936  if (morr_arr(i,j,k,MORRInd::qg3d) >= Real(0.1e-3)) {
1937  if (morr_arr(i,j,k,MORRInd::qc3d) >= Real(0.5e-3) || morr_arr(i,j,k,MORRInd::qr3d) >= Real(0.1e-3)) {
1938  if (psacwg > Real(0) || pracg > Real(0)) {
1939  if (morr_arr(i,j,k,MORRInd::t3d) < Real(270.16) && morr_arr(i,j,k,MORRInd::t3d) > Real(265.16)) {
1940  Real fmult = Real(0);
1941 
1942  if (morr_arr(i,j,k,MORRInd::t3d) > Real(270.16)) {
1943  fmult = Real(0);
1944  } else if (morr_arr(i,j,k,MORRInd::t3d) <= Real(270.16) && morr_arr(i,j,k,MORRInd::t3d) > Real(268.16)) {
1945  fmult = (Real(270.16) - morr_arr(i,j,k,MORRInd::t3d)) / Real(2);
1946  } else if (morr_arr(i,j,k,MORRInd::t3d) >= Real(265.16) && morr_arr(i,j,k,MORRInd::t3d) <= Real(268.16)) {
1947  fmult = (morr_arr(i,j,k,MORRInd::t3d) - Real(265.16)) / three;
1948  } else if (morr_arr(i,j,k,MORRInd::t3d) < Real(265.16)) {
1949  fmult = Real(0);
1950  }
1951 
1952  // 1000 IS TO CONVERT FROM KG TO G
1953  // SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
1954  if (psacwg > Real(0)) {
1955  nmultg = Real(35.0e4) * psacwg * fmult * Real(1000.0);
1956  qmultg = nmultg * m_mmult;
1957 
1958  // CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
1959  // THAN WAS RIMED ONTO GRAUPEL
1960  qmultg = std::min(qmultg, psacwg);
1961  psacwg = psacwg - qmultg;
1962  }
1963 
1964  // RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
1965  if (pracg > Real(0)) {
1966  nmultrg = Real(35.0e4) * pracg * fmult * Real(1000.0);
1967  qmultrg = nmultrg * m_mmult;
1968 
1969  // CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
1970  // THAN WAS RIMED ONTO GRAUPEL
1971  qmultrg = std::min(qmultrg, pracg);
1972  pracg = pracg - qmultrg;
1973  }
1974  }
1975  }
1976  }
1977  }
1978 
1979  // CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL/HAIL
1980  if (psacws > Real(0)) {
1981  // ONLY ALLOW CONVERSION IF QNI > Real(0.1) AND QC > myhalf G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
1982  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(0.1e-3) && morr_arr(i,j,k,MORRInd::qc3d) >= Real(0.5e-3)) {
1983  // PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
1984  pgsacw = std::min(psacws, m_cons17 * dt * morr_arr(i,j,k,MORRInd::n0s) * morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::qc3d) *
1985  morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::asn) /
1986  (morr_arr(i,j,k,MORRInd::rho) * std::pow(morr_arr(i,j,k,MORRInd::lams), (Real(2) * m_bs + Real(2)))));
1987 
1988  // MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
1989  dum = std::max(m_rhosn / (m_rhog - m_rhosn) * pgsacw, Real(0));
1990 
1991  // NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
1992  nscng = dum / m_mg0 * morr_arr(i,j,k,MORRInd::rho);
1993  // LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
1994  nscng = std::min(nscng, morr_arr(i,j,k,MORRInd::ns3d) / dt);
1995 
1996  // PORTION OF RIMING LEFT FOR SNOW
1997  psacws = psacws - pgsacw;
1998  }
1999  }
2000 
2001  // CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2002  if (pracs > Real(0)) {
2003  // ONLY ALLOW CONVERSION IF QNI > Real(0.1) AND QR > Real(0.1) G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2004  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(0.1e-3) && morr_arr(i,j,k,MORRInd::qr3d) >= Real(0.1e-3)) {
2005  // PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2006  dum = m_cons18 * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lams)) * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lams)) /
2007  (m_cons18 * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lams)) * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lams)) +
2008  m_cons19 * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lamr)) * amrex::Math::powi<3>(Real(4.0) / morr_arr(i,j,k,MORRInd::lamr)));
2009  dum = std::min(dum, Real(1));
2010  dum = std::max(dum, Real(0));
2011 
2012  pgracs = (one - dum) * pracs;
2013  ngracs = (one - dum) * npracs;
2014 
2015  // LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2016  ngracs = std::min(ngracs, morr_arr(i,j,k,MORRInd::nr3d) / dt);
2017  ngracs = std::min(ngracs, morr_arr(i,j,k,MORRInd::ns3d) / dt);
2018 
2019  // AMOUNT LEFT FOR SNOW PRODUCTION
2020  pracs = pracs - pgracs;
2021  npracs = npracs - ngracs;
2022 
2023  // CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2024  psacr = psacr * (one - dum);
2025  }
2026  }
2027 
2028  // FREEZING OF RAIN DROPS
2029  // FREEZING ALLOWED BELOW -4 C
2030  if (morr_arr(i,j,k,MORRInd::t3d) < Real(269.15) && morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
2031  // IMMERSION FREEZING (BIGG 1953)
2032  // hm fix 7/15/13 for consistency w/ original formula
2033  mnuccr = m_cons20 * morr_arr(i,j,k,MORRInd::nr3d) * (std::exp(m_aimm * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) - one) /
2034  amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) / amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr));
2035 
2036  nnuccr = m_pi * morr_arr(i,j,k,MORRInd::nr3d) * m_bimm * (std::exp(m_aimm * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) - one) /
2037  amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr));
2038 
2039  // PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2040  nnuccr = std::min(nnuccr, morr_arr(i,j,k,MORRInd::nr3d) / dt);
2041  }
2042 
2043  // ACCRETION OF CLOUD LIQUID WATER BY RAIN
2044  // CONTINUOUS COLLECTION EQUATION WITH
2045  // GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2046  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qc3d) >= Real(1.0e-8)) {
2047  // 12/13/06 hm add, replace with newer formula from
2048  // khairoutdinov and kogan 2000, mwr
2049  dum = morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::qr3d);
2050  pra = Real(67.0) * std::pow(dum, Real(1.15));
2051  npra = pra / (morr_arr(i,j,k,MORRInd::qc3d) / morr_arr(i,j,k,MORRInd::nc3d));
2052  }
2053 
2054  // SELF-COLLECTION OF RAIN DROPS
2055  // FROM BEHENG(1994)
2056  // FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2057  // AS DESCRINED ABOVE FOR AUTOCONVERSION
2058  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8)) {
2059  // include breakup add 10/09/09
2060  dum1 = Real(300.0e-6);
2061  if (one / morr_arr(i,j,k,MORRInd::lamr) < dum1) {
2062  dum = one;
2063  } else if (one / morr_arr(i,j,k,MORRInd::lamr) >= dum1) {
2064  dum = Real(2) - std::exp(Real(2300.0) * (one / morr_arr(i,j,k,MORRInd::lamr) - dum1));
2065  }
2066  nragg = -Real(5.78) * dum * morr_arr(i,j,k,MORRInd::nr3d) * morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::rho);
2067  }
2068 
2069  // AUTOCONVERSION OF CLOUD ICE TO SNOW
2070  // FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2071  // HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2072  // ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2073  if (morr_arr(i,j,k,MORRInd::qi3d) >= Real(1.0e-8) && qvqvsi >= one) {
2074  nprci = m_cons21 * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) * morr_arr(i,j,k,MORRInd::rho) *
2075  morr_arr(i,j,k,MORRInd::n0i) * std::exp(-morr_arr(i,j,k,MORRInd::lami) * m_dcs) * dv / abi;
2076  prci = m_cons22 * nprci;
2077  nprci = std::min(nprci, morr_arr(i,j,k,MORRInd::ni3d) / dt);
2078  }
2079 
2080  // ACCRETION OF CLOUD ICE BY SNOW
2081  // FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2082  // AND DS >> DI FOR CONTINUOUS COLLECTION
2083  if (morr_arr(i,j,k,MORRInd::qni3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
2084  prai = m_cons23 * morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::rho) * morr_arr(i,j,k,MORRInd::n0s) /
2085  std::pow(morr_arr(i,j,k,MORRInd::lams), (m_bs + three));
2086  nprai = m_cons23 * morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::ni3d) *
2087  morr_arr(i,j,k,MORRInd::rho) * morr_arr(i,j,k,MORRInd::n0s) /
2088  std::pow(morr_arr(i,j,k,MORRInd::lams), (m_bs + three));
2089  nprai = std::min(nprai, morr_arr(i,j,k,MORRInd::ni3d) / dt);
2090  }
2091 
2092  // hm, add 12/13/06, collision of rain and ice to produce snow or graupel
2093  // follows reisner et al. 1998
2094  // assumed fallspeed and size of ice crystal << than for rain
2095  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::qi3d) >= Real(1.0e-8) && morr_arr(i,j,k,MORRInd::t3d) <= Real(273.15)) {
2096  // allow graupel formation from rain-ice collisions only if rain mixing ratio > Real(0.1) g/kg,
2097  // otherwise add to snow
2098  if (morr_arr(i,j,k,MORRInd::qr3d) >= Real(0.1e-3)) {
2099  niacr = m_cons24 * morr_arr(i,j,k,MORRInd::ni3d) * morr_arr(i,j,k,MORRInd::n0r)* morr_arr(i,j,k,MORRInd::arn) /
2100  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) * morr_arr(i,j,k,MORRInd::rho);
2101  piacr = m_cons25 * morr_arr(i,j,k,MORRInd::ni3d) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::arn) /
2102  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) / amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::rho);
2103  praci = m_cons24 * morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::arn) /
2104  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) * morr_arr(i,j,k,MORRInd::rho);
2105  niacr = std::min(niacr, morr_arr(i,j,k,MORRInd::nr3d) / dt);
2106  niacr = std::min(niacr, morr_arr(i,j,k,MORRInd::ni3d) / dt);
2107  } else {
2108  niacrs = m_cons24 * morr_arr(i,j,k,MORRInd::ni3d) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::arn) /
2109  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) * morr_arr(i,j,k,MORRInd::rho);
2110  piacrs = m_cons25 * morr_arr(i,j,k,MORRInd::ni3d) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::arn) /
2111  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) / amrex::Math::powi<3>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::rho);
2112  pracis = m_cons24 * morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::arn) /
2113  std::pow(morr_arr(i,j,k,MORRInd::lamr), (m_br + three)) * morr_arr(i,j,k,MORRInd::rho);
2114  niacrs = std::min(niacrs, morr_arr(i,j,k,MORRInd::nr3d) / dt);
2115  niacrs = std::min(niacrs, morr_arr(i,j,k,MORRInd::ni3d) / dt);
2116  }
2117  }
2118 
2119  // NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2120  if (m_inuc == 0) {
2121  // ADD THRESHOLD ACCORDING TO GREG THOMSPON
2122  if ((qvqvs >= Real(0.999) && morr_arr(i,j,k,MORRInd::t3d) <= Real(265.15)) || qvqvsi >= Real(1.08)) {
2123  // hm, modify dec. 5, 2006, replace with cooper curve
2124  kc2 = Real(0.005) * std::exp(Real(0.304) * (Real(273.15) - morr_arr(i,j,k,MORRInd::t3d))) * Real(1000.0); // CONVERT FROM L-1 TO M-3
2125  // LIMIT TO 500 L-1
2126  kc2 = std::min(kc2, Real(500.0e3));
2127  kc2 = std::max(kc2 / morr_arr(i,j,k,MORRInd::rho), Real(0)); // CONVERT TO KG-1
2128 
2129  if (kc2 > morr_arr(i,j,k,MORRInd::ni3d) + morr_arr(i,j,k,MORRInd::ns3d) + morr_arr(i,j,k,MORRInd::ng3d)) {
2130  nnuccd = (kc2 - morr_arr(i,j,k,MORRInd::ni3d) - morr_arr(i,j,k,MORRInd::ns3d) - morr_arr(i,j,k,MORRInd::ng3d)) / dt;
2131  mnuccd = nnuccd * m_mi0;
2132  }
2133  }
2134  } else if (m_inuc == 1) {
2135  if (morr_arr(i,j,k,MORRInd::t3d) < Real(273.15) && qvqvsi > one) {
2136  kc2 = Real(0.16) * Real(1000.0) / morr_arr(i,j,k,MORRInd::rho); // CONVERT FROM L-1 TO KG-1
2137  if (kc2 > morr_arr(i,j,k,MORRInd::ni3d) + morr_arr(i,j,k,MORRInd::ns3d) + morr_arr(i,j,k,MORRInd::ng3d)) {
2138  nnuccd = (kc2 - morr_arr(i,j,k,MORRInd::ni3d) - morr_arr(i,j,k,MORRInd::ns3d) - morr_arr(i,j,k,MORRInd::ng3d)) / dt;
2139  mnuccd = nnuccd * m_mi0;
2140  }
2141  }
2142  }
2143 
2144  // CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2145  // NO VENTILATION FOR CLOUD ICE
2146  epsi = Real(0);
2147  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
2148  epsi = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0i) * morr_arr(i,j,k,MORRInd::rho) * dv / (morr_arr(i,j,k,MORRInd::lami) * morr_arr(i,j,k,MORRInd::lami));
2149  }
2150 
2151  // VENTILATION FOR SNOW
2152  epss = Real(0);
2153  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
2154  epss = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0s) * morr_arr(i,j,k,MORRInd::rho) * dv *
2155  (m_f1s / (morr_arr(i,j,k,MORRInd::lams) * morr_arr(i,j,k,MORRInd::lams)) +
2156  m_f2s * std::pow(morr_arr(i,j,k,MORRInd::asn) * morr_arr(i,j,k,MORRInd::rho) / morr_arr(i,j,k,MORRInd::mu), myhalf) *
2157  std::pow(sc_schmidt, (one / three)) * m_cons10 /
2158  std::pow(morr_arr(i,j,k,MORRInd::lams), m_cons35));
2159  }
2160 
2161  // Ventilation for graupel
2162  epsg = Real(0);
2163  if (morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
2164  epsg = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0g) * morr_arr(i,j,k,MORRInd::rho) * dv *
2165  (m_f1s / (morr_arr(i,j,k,MORRInd::lamg) * morr_arr(i,j,k,MORRInd::lamg)) +
2166  m_f2s * std::pow(morr_arr(i,j,k,MORRInd::agn) * morr_arr(i,j,k,MORRInd::rho) / morr_arr(i,j,k,MORRInd::mu), myhalf) *
2167  std::pow(sc_schmidt, (one / three)) * m_cons11 /
2168  std::pow(morr_arr(i,j,k,MORRInd::lamg), m_cons36));
2169  }
2170 
2171  // VENTILATION FOR RAIN
2172  epsr = Real(0);
2173  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
2174  epsr = Real(2) * m_pi * morr_arr(i,j,k,MORRInd::n0r) * morr_arr(i,j,k,MORRInd::rho) * dv *
2175  (m_f1r / (morr_arr(i,j,k,MORRInd::lamr) * morr_arr(i,j,k,MORRInd::lamr)) +
2176  m_f2r * std::pow(morr_arr(i,j,k,MORRInd::arn) * morr_arr(i,j,k,MORRInd::rho) / morr_arr(i,j,k,MORRInd::mu), myhalf) *
2177  std::pow(sc_schmidt, (one / three)) * m_cons9 /
2178  std::pow(morr_arr(i,j,k,MORRInd::lamr), m_cons34));
2179  }
2180 
2181  // ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2182  // DUM IS FRACTION OF D*N(D) < DCS
2183  // LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2184  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
2185  dum = (one - std::exp(-morr_arr(i,j,k,MORRInd::lami) * m_dcs) * (one + morr_arr(i,j,k,MORRInd::lami) * m_dcs));
2186  prd = epsi * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / abi * dum;
2187  } else {
2188  dum = Real(0);
2189  }
2190 
2191  // ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2192  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
2193  prds = epss * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / abi +
2194  epsi * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / abi * (one - dum);
2195  } else {
2196  // OTHERWISE ADD TO CLOUD ICE
2197  prd = prd + epsi * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / abi * (one - dum);
2198  }
2199 
2200  // VAPOR DPEOSITION ON GRAUPEL
2201  prdg = epsg * (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / abi;
2202 
2203  // NO CONDENSATION ONTO RAIN, ONLY EVAP
2204  if (morr_arr(i,j,k,MORRInd::qv3d) < qvs) {
2205  pre = epsr * (morr_arr(i,j,k,MORRInd::qv3d) - qvs) / ab;
2206  pre = std::min(pre, Real(0));
2207  } else {
2208  pre = Real(0);
2209  }
2210 
2211  // MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2212  // FORMULA FROM REISNER 2 SCHEME
2213  dum = (morr_arr(i,j,k,MORRInd::qv3d) - qvi) / dt;
2214 
2215  fudgef = Real(0.9999);
2216  sum_dep = prd + prds + mnuccd + prdg;
2217 
2218  if ((dum > Real(0) && sum_dep > dum * fudgef) ||
2219  (dum < Real(0) && sum_dep < dum * fudgef)) {
2220  mnuccd = fudgef * mnuccd * dum / sum_dep;
2221  prd = fudgef * prd * dum / sum_dep;
2222  prds = fudgef * prds * dum / sum_dep;
2223  prdg = fudgef * prdg * dum / sum_dep;
2224  }
2225 
2226  // IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2227  if (prd < Real(0)) {
2228  eprd = prd;
2229  prd = Real(0);
2230  }
2231  if (prds < Real(0)) {
2232  eprds = prds;
2233  prds = Real(0);
2234  }
2235  if (prdg < Real(0)) {
2236  eprdg = prdg;
2237  prdg = Real(0);
2238  }
2239  // CONSERVATION OF WATER
2240  // THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2241  // ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2242 
2243  // IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2244  // THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2245 
2246  // NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2247  // BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2248  // FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2249  // TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2250  // STEP
2251 
2252  // ****SENSITIVITY - NO ICE
2253  if (m_iliq == 1) {
2254  mnuccc = Real(0);
2255  nnuccc = Real(0);
2256  mnuccr = Real(0);
2257  nnuccr = Real(0);
2258  mnuccd = Real(0);
2259  nnuccd = Real(0);
2260  }
2261 
2262  // ****SENSITIVITY - NO GRAUPEL
2263  if (m_igraup == 1) {
2264  pracg = Real(0);
2265  psacr = Real(0);
2266  psacwg = Real(0);
2267  prdg = Real(0);
2268  eprdg = Real(0);
2269  evpmg = Real(0);
2270  pgmlt = Real(0);
2271  npracg = Real(0);
2272  npsacwg = Real(0);
2273  nscng = Real(0);
2274  ngracs = Real(0);
2275  nsubg = Real(0);
2276  ngmltg = Real(0);
2277  ngmltr = Real(0);
2278 
2279  // fix 053011
2280  piacrs = piacrs + piacr;
2281  piacr = Real(0);
2282 
2283  // fix 070713
2284  pracis = pracis + praci;
2285  praci = Real(0);
2286  psacws = psacws + pgsacw;
2287  pgsacw = Real(0);
2288  pracs = pracs + pgracs;
2289  pgracs = Real(0);
2290  }
2291 
2292  // CONSERVATION OF QC
2293  dum = (prc + pra + mnuccc + psacws + psacwi + qmults + psacwg + pgsacw + qmultg) * dt;
2294 
2295  if (dum > morr_arr(i,j,k,MORRInd::qc3d) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
2296  ratio = morr_arr(i,j,k,MORRInd::qc3d) / dum;
2297 
2298  prc = prc * ratio;
2299  pra = pra * ratio;
2300  mnuccc = mnuccc * ratio;
2301  psacws = psacws * ratio;
2302  psacwi = psacwi * ratio;
2303  qmults = qmults * ratio;
2304  qmultg = qmultg * ratio;
2305  psacwg = psacwg * ratio;
2306  pgsacw = pgsacw * ratio;
2307  }
2308 
2309  // CONSERVATION OF QI
2310  dum = (-prd - mnuccc + prci + prai - qmults - qmultg - qmultr - qmultrg
2311  - mnuccd + praci + pracis - eprd - psacwi) * dt;
2312 
2313  if (dum > morr_arr(i,j,k,MORRInd::qi3d) && morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
2314  ratio = (morr_arr(i,j,k,MORRInd::qi3d) / dt + prd + mnuccc + qmults + qmultg + qmultr + qmultrg +
2315  mnuccd + psacwi) /
2316  (prci + prai + praci + pracis - eprd);
2317 
2318  prci = prci * ratio;
2319  prai = prai * ratio;
2320  praci = praci * ratio;
2321  pracis = pracis * ratio;
2322  eprd = eprd * ratio;
2323  }
2324 
2325  // CONSERVATION OF QR
2326  dum = ((pracs - pre) + (qmultr + qmultrg - prc) + (mnuccr - pra) +
2327  piacr + piacrs + pgracs + pracg) * dt;
2328 
2329  if (dum > morr_arr(i,j,k,MORRInd::qr3d) && morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
2330  ratio = (morr_arr(i,j,k,MORRInd::qr3d) / dt + prc + pra) /
2331  (-pre + qmultr + qmultrg + pracs + mnuccr + piacr + piacrs + pgracs + pracg);
2332 
2333  pre = pre * ratio;
2334  pracs = pracs * ratio;
2335  qmultr = qmultr * ratio;
2336  qmultrg = qmultrg * ratio;
2337  mnuccr = mnuccr * ratio;
2338  piacr = piacr * ratio;
2339  piacrs = piacrs * ratio;
2340  pgracs = pgracs * ratio;
2341  pracg = pracg * ratio;
2342  }
2343 
2344  // CONSERVATION OF QNI
2345  if (m_igraup == 0) {
2346  dum = (-prds - psacws - prai - prci - pracs - eprds + psacr - piacrs - pracis) * dt;
2347 
2348  if (dum > morr_arr(i,j,k,MORRInd::qni3d) && morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
2349  ratio = (morr_arr(i,j,k,MORRInd::qni3d) / dt + prds + psacws + prai + prci + pracs + piacrs + pracis) /
2350  (-eprds + psacr);
2351 
2352  eprds = eprds * ratio;
2353  psacr = psacr * ratio;
2354  }
2355  } else if (m_igraup == 1) {
2356  // FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2357  dum = (-prds - psacws - prai - prci - pracs - eprds + psacr - piacrs - pracis - mnuccr) * dt;
2358 
2359  if (dum > morr_arr(i,j,k,MORRInd::qni3d) && morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
2360  ratio = (morr_arr(i,j,k,MORRInd::qni3d) / dt + prds + psacws + prai + prci + pracs + piacrs + pracis + mnuccr) /
2361  (-eprds + psacr);
2362 
2363  eprds = eprds * ratio;
2364  psacr = psacr * ratio;
2365  }
2366  }
2367 
2368  // CONSERVATION OF QG
2369  dum = (-psacwg - pracg - pgsacw - pgracs - prdg - mnuccr - eprdg - piacr - praci - psacr) * dt;
2370 
2371  if (dum > morr_arr(i,j,k,MORRInd::qg3d) && morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
2372  ratio = (morr_arr(i,j,k,MORRInd::qg3d) / dt + psacwg + pracg + pgsacw + pgracs + prdg + mnuccr + psacr +
2373  piacr + praci) / (-eprdg);
2374 
2375  eprdg = eprdg * ratio;
2376  }
2377 
2378  // TENDENCIES
2379  morr_arr(i,j,k,MORRInd::qv3dten) = morr_arr(i,j,k,MORRInd::qv3dten) + (-pre - prd - prds - mnuccd - eprd - eprds - prdg - eprdg);
2380 
2381  // bug fix hm, 3/1/11, include piacr and piacrs
2382  morr_arr(i,j,k,MORRInd::t3dten) = morr_arr(i,j,k,MORRInd::t3dten) + (pre * morr_arr(i,j,k,MORRInd::xxlv) +
2383  (prd + prds + mnuccd + eprd + eprds + prdg + eprdg) * morr_arr(i,j,k,MORRInd::xxls) +
2384  (psacws + psacwi + mnuccc + mnuccr + qmults + qmultg + qmultr + qmultrg + pracs +
2385  psacwg + pracg + pgsacw + pgracs + piacr + piacrs) * morr_arr(i,j,k,MORRInd::xlf)) / morr_arr(i,j,k,MORRInd::cpm);
2386 
2387  morr_arr(i,j,k,MORRInd::qc3dten) = morr_arr(i,j,k,MORRInd::qc3dten) +
2388  (-pra - prc - mnuccc + pcc -
2389  psacws - psacwi - qmults - qmultg - psacwg - pgsacw);
2390 
2391  morr_arr(i,j,k,MORRInd::qi3dten) = morr_arr(i,j,k,MORRInd::qi3dten) +
2392  (prd + eprd + psacwi + mnuccc - prci -
2393  prai + qmults + qmultg + qmultr + qmultrg + mnuccd - praci - pracis);
2394 
2395  morr_arr(i,j,k,MORRInd::qr3dten) = morr_arr(i,j,k,MORRInd::qr3dten) +
2396  (pre + pra + prc - pracs - mnuccr - qmultr - qmultrg -
2397  piacr - piacrs - pracg - pgracs);
2398  if (m_igraup == 0) {
2399  morr_arr(i,j,k,MORRInd::qni3dten) = morr_arr(i,j,k,MORRInd::qni3dten) +
2400  (prai + psacws + prds + pracs + prci + eprds - psacr + piacrs + pracis);
2401 
2402  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + (nsagg + nprci - nscng - ngracs + niacrs);
2403 
2404  morr_arr(i,j,k,MORRInd::qg3dten) = morr_arr(i,j,k,MORRInd::qg3dten) + (pracg + psacwg + pgsacw + pgracs +
2405  prdg + eprdg + mnuccr + piacr + praci + psacr);
2406 
2407  morr_arr(i,j,k,MORRInd::ng3dten) = morr_arr(i,j,k,MORRInd::ng3dten) + (nscng + ngracs + nnuccr + niacr);
2408  } else if (m_igraup == 1) {
2409  // FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
2410  morr_arr(i,j,k,MORRInd::qni3dten) = morr_arr(i,j,k,MORRInd::qni3dten) +
2411  (prai + psacws + prds + pracs + prci + eprds - psacr + piacrs + pracis + mnuccr);
2412 
2413  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + (nsagg + nprci - nscng - ngracs + niacrs + nnuccr);
2414  }
2415 
2416  morr_arr(i,j,k,MORRInd::nc3dten) = morr_arr(i,j,k,MORRInd::nc3dten) + (-nnuccc - npsacws -
2417  npra - nprc - npsacwi - npsacwg);
2418 
2419  morr_arr(i,j,k,MORRInd::ni3dten) = morr_arr(i,j,k,MORRInd::ni3dten) +
2420  (nnuccc - nprci - nprai + nmults + nmultg + nmultr + nmultrg +
2421  nnuccd - niacr - niacrs);
2422 
2423  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) + (nprc1 - npracs - nnuccr +
2424  nragg - niacr - niacrs - npracg - ngracs);
2425 
2426  // hm add, wrf-chem, add tendencies for c2prec
2427  // c2prec = pra + prc + psacws + qmults + qmultg + psacwg + pgsacw + mnuccc + psacwi;
2428 
2429  // CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
2430  // WATER SATURATION
2431  dumt = morr_arr(i,j,k,MORRInd::t3d) + dt * morr_arr(i,j,k,MORRInd::t3dten);
2432  dumqv = morr_arr(i,j,k,MORRInd::qv3d) + dt * morr_arr(i,j,k,MORRInd::qv3dten);
2433 
2434  // hm, add fix for low pressure, 5/12/10
2435  dum = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(dumt, 0));
2436  dumqss = m_ep_2 * dum / (morr_arr(i,j,k,MORRInd::pres) - dum);
2437 
2438  dumqc = morr_arr(i,j,k,MORRInd::qc3d) + dt * morr_arr(i,j,k,MORRInd::qc3dten);
2439  dumqc = std::max(dumqc, Real(0));
2440 
2441  // SATURATION ADJUSTMENT FOR LIQUID
2442  dums = dumqv - dumqss;
2443 
2444  pcc = dums / (one + amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::xxlv)) * dumqss / (morr_arr(i,j,k,MORRInd::cpm) * m_Rv * amrex::Math::powi<2>(dumt))) / dt;
2445 
2446  if (pcc * dt + dumqc < Real(0)) {
2447  pcc = -dumqc / dt;
2448  }
2449 
2450  morr_arr(i,j,k,MORRInd::qv3dten) = morr_arr(i,j,k,MORRInd::qv3dten) - pcc;
2451  morr_arr(i,j,k,MORRInd::t3dten) = morr_arr(i,j,k,MORRInd::t3dten) + pcc * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm);
2452  morr_arr(i,j,k,MORRInd::qc3dten) = morr_arr(i,j,k,MORRInd::qc3dten) + pcc;
2453  // SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
2454  // THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
2455  // LOSS OF NUMBER CONCENTRATION
2456  if (eprd < Real(0)) {
2457  dum = eprd * dt / morr_arr(i,j,k,MORRInd::qi3d);
2458  dum = std::max(-one, dum);
2459  nsubi = dum * morr_arr(i,j,k,MORRInd::ni3d) / dt;
2460  }
2461 
2462  if (eprds < Real(0)) {
2463  dum = eprds * dt / morr_arr(i,j,k,MORRInd::qni3d);
2464  dum = std::max(-one, dum);
2465  nsubs = dum * morr_arr(i,j,k,MORRInd::ns3d) / dt;
2466  }
2467 
2468  if (pre < Real(0)) {
2469  dum = pre * dt / morr_arr(i,j,k,MORRInd::qr3d);
2470  dum = std::max(-one, dum);
2471  nsubr = dum * morr_arr(i,j,k,MORRInd::nr3d) / dt;
2472  }
2473 
2474  if (eprdg < Real(0)) {
2475  dum = eprdg * dt / morr_arr(i,j,k,MORRInd::qg3d);
2476  dum = std::max(-one, dum);
2477  nsubg = dum * morr_arr(i,j,k,MORRInd::ng3d) / dt;
2478  }
2479 
2480  // UPDATE TENDENCIES
2481  morr_arr(i,j,k,MORRInd::ni3dten) = morr_arr(i,j,k,MORRInd::ni3dten) + nsubi;
2482  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + nsubs;
2483  morr_arr(i,j,k,MORRInd::ng3dten) = morr_arr(i,j,k,MORRInd::ng3dten) + nsubg;
2484  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) + nsubr;
2485  }
2486  ltrue = 1;
2487  }
2488  // label_200:
2489  } // k
2490 
2491  for(int k=klo; k<=khi; k++) {
2492  // INITIALIZE PRECIP AND SNOW RATES
2493  morr_arr(i,j,k,MORRInd::precrt) = Real(0);
2494  morr_arr(i,j,k,MORRInd::snowrt) = Real(0);
2495  // hm added 7/13/13
2496  morr_arr(i,j,k,MORRInd::snowprt) = Real(0);
2497  morr_arr(i,j,k,MORRInd::grplprt) = Real(0);
2498  } // k
2499 
2500  nstep = 1;
2501 
2502  if (ltrue != 0) {
2503  //goto 400
2504  // CALCULATE SEDIMENTATION
2505  // THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
2506  // FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
2507  // STABILITY, I.E. COURANT# < 1
2508  // Loop from top to bottom (KTE to KTS)
2509  for(int k=khi; k>=klo; k--) {
2510 
2511  Real dum; // DUM: General dummy variable
2512 
2513  Real di0; // DI0: Characteristic diameter for ice
2514  Real ds0; // DS0: Characteristic diameter for snow
2515  Real dg0; // DG0: Characteristic diameter for graupel
2516  Real lammax; // LAMMAX: Maximum value for slope parameter
2517  Real lammin; // LAMMIN: Minimum value for slope parameter
2518 
2519  ds0 = three; // Size distribution parameter for snow
2520  di0 = three; // Size distribution parameter for cloud ice
2521  dg0 = three; // Size distribution parameter for graupel
2522 
2523  // Update prognostic variables with tendencies
2524  morr_arr(i,j,k,MORRInd::dumi) = morr_arr(i,j,k,MORRInd::qi3d) + morr_arr(i,j,k,MORRInd::qi3dten) * dt;
2525  morr_arr(i,j,k,MORRInd::dumqs) = morr_arr(i,j,k,MORRInd::qni3d) + morr_arr(i,j,k,MORRInd::qni3dten) * dt;
2526  morr_arr(i,j,k,MORRInd::dumr) = morr_arr(i,j,k,MORRInd::qr3d) + morr_arr(i,j,k,MORRInd::qr3dten) * dt;
2527  morr_arr(i,j,k,MORRInd::dumfni) = morr_arr(i,j,k,MORRInd::ni3d) + morr_arr(i,j,k,MORRInd::ni3dten) * dt;
2528  morr_arr(i,j,k,MORRInd::dumfns) = morr_arr(i,j,k,MORRInd::ns3d) + morr_arr(i,j,k,MORRInd::ns3dten) * dt;
2529  morr_arr(i,j,k,MORRInd::dumfnr) = morr_arr(i,j,k,MORRInd::nr3d) + morr_arr(i,j,k,MORRInd::nr3dten) * dt;
2530  morr_arr(i,j,k,MORRInd::dumc) = morr_arr(i,j,k,MORRInd::qc3d) + morr_arr(i,j,k,MORRInd::qc3dten) * dt;
2531  morr_arr(i,j,k,MORRInd::dumfnc) = morr_arr(i,j,k,MORRInd::nc3d) + morr_arr(i,j,k,MORRInd::nc3dten) * dt;
2532  morr_arr(i,j,k,MORRInd::dumg) = morr_arr(i,j,k,MORRInd::qg3d) + morr_arr(i,j,k,MORRInd::qg3dten) * dt;
2533  morr_arr(i,j,k,MORRInd::dumfng) = morr_arr(i,j,k,MORRInd::ng3d) + morr_arr(i,j,k,MORRInd::ng3dten) * dt;
2534 
2535  // SWITCH FOR CONSTANT DROPLET NUMBER
2536  if (iinum == 1) {
2537  morr_arr(i,j,k,MORRInd::dumfnc) = morr_arr(i,j,k,MORRInd::nc3d);
2538  }
2539 
2540  // MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
2541  morr_arr(i,j,k,MORRInd::dumfni) = amrex::max(Real(0), morr_arr(i,j,k,MORRInd::dumfni));
2542  morr_arr(i,j,k,MORRInd::dumfns) = amrex::max(Real(0), morr_arr(i,j,k,MORRInd::dumfns));
2543  morr_arr(i,j,k,MORRInd::dumfnc) = amrex::max(Real(0), morr_arr(i,j,k,MORRInd::dumfnc));
2544  morr_arr(i,j,k,MORRInd::dumfnr) = amrex::max(Real(0), morr_arr(i,j,k,MORRInd::dumfnr));
2545  morr_arr(i,j,k,MORRInd::dumfng) = amrex::max(Real(0), morr_arr(i,j,k,MORRInd::dumfng));
2546 
2547  // CLOUD ICE
2548  if (morr_arr(i,j,k,MORRInd::dumi) >= m_qsmall) {
2549  morr_arr(i,j,k,MORRInd::dlami) = std::pow(m_cons12 * morr_arr(i,j,k,MORRInd::dumfni) / morr_arr(i,j,k,MORRInd::dumi), one/di0);
2550  morr_arr(i,j,k,MORRInd::dlami) = amrex::max(morr_arr(i,j,k,MORRInd::dlami), m_lammini);
2551  morr_arr(i,j,k,MORRInd::dlami) = amrex::min(morr_arr(i,j,k,MORRInd::dlami), m_lammaxi);
2552  }
2553 
2554  // RAIN
2555  if (morr_arr(i,j,k,MORRInd::dumr) >= m_qsmall) {
2556  morr_arr(i,j,k,MORRInd::dlamr) = std::pow(m_pi * m_rhow * morr_arr(i,j,k,MORRInd::dumfnr) / morr_arr(i,j,k,MORRInd::dumr), one/three);
2557  morr_arr(i,j,k,MORRInd::dlamr) = amrex::max(morr_arr(i,j,k,MORRInd::dlamr), m_lamminr);
2558  morr_arr(i,j,k,MORRInd::dlamr) = amrex::min(morr_arr(i,j,k,MORRInd::dlamr), m_lammaxr);
2559  }
2560 
2561  // CLOUD DROPLETS
2562  if (morr_arr(i,j,k,MORRInd::dumc) >= m_qsmall) {
2563  dum = morr_arr(i,j,k,MORRInd::pres) / (Real(287.15) * morr_arr(i,j,k,MORRInd::t3d));
2564  morr_arr(i,j,k,MORRInd::pgam) = Real(0.0005714) * (morr_arr(i,j,k,MORRInd::nc3d) / Real(1.0e6) * dum) + Real(0.2714);
2565  morr_arr(i,j,k,MORRInd::pgam) = one / (morr_arr(i,j,k,MORRInd::pgam) * morr_arr(i,j,k,MORRInd::pgam)) - one;
2566  morr_arr(i,j,k,MORRInd::pgam) = amrex::max(morr_arr(i,j,k,MORRInd::pgam), Real(2));
2567  morr_arr(i,j,k,MORRInd::pgam) = amrex::min(morr_arr(i,j,k,MORRInd::pgam), Real(10.0));
2568 
2569  morr_arr(i,j,k,MORRInd::dlamc) = std::pow(m_cons26 * morr_arr(i,j,k,MORRInd::dumfnc) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0)) /
2570  (morr_arr(i,j,k,MORRInd::dumc) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one)), one/three);
2571  lammin = (morr_arr(i,j,k,MORRInd::pgam) + one) / Real(60.0e-6);
2572  lammax = (morr_arr(i,j,k,MORRInd::pgam) + one) / Real(1.0e-6);
2573  morr_arr(i,j,k,MORRInd::dlamc) = amrex::max(morr_arr(i,j,k,MORRInd::dlamc), lammin);
2574  morr_arr(i,j,k,MORRInd::dlamc) = amrex::min(morr_arr(i,j,k,MORRInd::dlamc), lammax);
2575  }
2576 
2577  // SNOW
2578  if (morr_arr(i,j,k,MORRInd::dumqs) >= m_qsmall) {
2579  morr_arr(i,j,k,MORRInd::dlams) = std::pow(m_cons1 * morr_arr(i,j,k,MORRInd::dumfns) / morr_arr(i,j,k,MORRInd::dumqs), one/ds0);
2580  morr_arr(i,j,k,MORRInd::dlams) = amrex::max(morr_arr(i,j,k,MORRInd::dlams), m_lammins);
2581  morr_arr(i,j,k,MORRInd::dlams) = amrex::min(morr_arr(i,j,k,MORRInd::dlams), m_lammaxs);
2582  }
2583 
2584  // GRAUPEL
2585  if (morr_arr(i,j,k,MORRInd::dumg) >= m_qsmall) {
2586  morr_arr(i,j,k,MORRInd::dlamg) = std::pow(m_cons2 * morr_arr(i,j,k,MORRInd::dumfng) / morr_arr(i,j,k,MORRInd::dumg), one/dg0);
2587  morr_arr(i,j,k,MORRInd::dlamg) = amrex::max(morr_arr(i,j,k,MORRInd::dlamg), m_lamming);
2588  morr_arr(i,j,k,MORRInd::dlamg) = amrex::min(morr_arr(i,j,k,MORRInd::dlamg), m_lammaxg);
2589  }
2590 
2591  // Calculate number-weighted and mass-weighted terminal fall speeds
2592  // CLOUD WATER
2593  if (morr_arr(i,j,k,MORRInd::dumc) >= m_qsmall) {
2594  morr_arr(i,j,k,MORRInd::unc) = morr_arr(i,j,k,MORRInd::acn) * gamma_function(one + m_bc + morr_arr(i,j,k,MORRInd::pgam)) /
2595  (std::pow(morr_arr(i,j,k,MORRInd::dlamc), m_bc) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one));
2596  morr_arr(i,j,k,MORRInd::umc) = morr_arr(i,j,k,MORRInd::acn) * gamma_function(Real(4.) + m_bc + morr_arr(i,j,k,MORRInd::pgam)) /
2597  (std::pow(morr_arr(i,j,k,MORRInd::dlamc), m_bc) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.)));
2598  } else {
2599  morr_arr(i,j,k,MORRInd::umc) = Real(0);
2600  morr_arr(i,j,k,MORRInd::unc) = Real(0);
2601  }
2602 
2603  // CLOUD ICE
2604  if (morr_arr(i,j,k,MORRInd::dumi) >= m_qsmall) {
2605  morr_arr(i,j,k,MORRInd::uni) = morr_arr(i,j,k,MORRInd::ain) * m_cons27 / std::pow(morr_arr(i,j,k,MORRInd::dlami), m_bi);
2606  morr_arr(i,j,k,MORRInd::umi) = morr_arr(i,j,k,MORRInd::ain) * m_cons28 / std::pow(morr_arr(i,j,k,MORRInd::dlami), m_bi);
2607  } else {
2608  morr_arr(i,j,k,MORRInd::umi) = Real(0);
2609  morr_arr(i,j,k,MORRInd::uni) = Real(0);
2610  }
2611 
2612  // RAIN
2613  if (morr_arr(i,j,k,MORRInd::dumr) >= m_qsmall) {
2614  morr_arr(i,j,k,MORRInd::unr) = morr_arr(i,j,k,MORRInd::arn) * m_cons6 / std::pow(morr_arr(i,j,k,MORRInd::dlamr), m_br);
2615  morr_arr(i,j,k,MORRInd::umr) = morr_arr(i,j,k,MORRInd::arn) * m_cons4 / std::pow(morr_arr(i,j,k,MORRInd::dlamr), m_br);
2616  } else {
2617  morr_arr(i,j,k,MORRInd::umr) = Real(0);
2618  morr_arr(i,j,k,MORRInd::unr) = Real(0);
2619  }
2620 
2621  // SNOW
2622  if (morr_arr(i,j,k,MORRInd::dumqs) >= m_qsmall) {
2623  morr_arr(i,j,k,MORRInd::ums) = morr_arr(i,j,k,MORRInd::asn) * m_cons3 / std::pow(morr_arr(i,j,k,MORRInd::dlams), m_bs);
2624  morr_arr(i,j,k,MORRInd::uns) = morr_arr(i,j,k,MORRInd::asn) * m_cons5 / std::pow(morr_arr(i,j,k,MORRInd::dlams), m_bs);
2625  } else {
2626  morr_arr(i,j,k,MORRInd::ums) = Real(0);
2627  morr_arr(i,j,k,MORRInd::uns) = Real(0);
2628  }
2629 
2630  // GRAUPEL
2631  if (morr_arr(i,j,k,MORRInd::dumg) >= m_qsmall) {
2632  morr_arr(i,j,k,MORRInd::umg) = morr_arr(i,j,k,MORRInd::agn) * m_cons7 / std::pow(morr_arr(i,j,k,MORRInd::dlamg), m_bg);
2633  morr_arr(i,j,k,MORRInd::ung) = morr_arr(i,j,k,MORRInd::agn) * m_cons8 / std::pow(morr_arr(i,j,k,MORRInd::dlamg), m_bg);
2634  } else {
2635  morr_arr(i,j,k,MORRInd::umg) = Real(0);
2636  morr_arr(i,j,k,MORRInd::ung) = Real(0);
2637  }
2638 
2639  // SET REALISTIC LIMITS ON FALLSPEED
2640  // Bug fix, 10/08/09
2641  dum = std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.54));
2642  morr_arr(i,j,k,MORRInd::ums) = std::min(morr_arr(i,j,k,MORRInd::ums), Real(1.2) * dum);
2643  morr_arr(i,j,k,MORRInd::uns) = std::min(morr_arr(i,j,k,MORRInd::uns), Real(1.2) * dum);
2644 
2645  // Fix 053011
2646  // Fix for correction by AA 4/6/11
2647  morr_arr(i,j,k,MORRInd::umi) = std::min(morr_arr(i,j,k,MORRInd::umi), Real(1.2) * std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.35)));
2648  morr_arr(i,j,k,MORRInd::uni) = std::min(morr_arr(i,j,k,MORRInd::uni), Real(1.2) * std::pow(m_rhosu / morr_arr(i,j,k,MORRInd::rho), Real(0.35)));
2649  morr_arr(i,j,k,MORRInd::umr) = std::min(morr_arr(i,j,k,MORRInd::umr), Real(9.1) * dum);
2650  morr_arr(i,j,k,MORRInd::unr) = std::min(morr_arr(i,j,k,MORRInd::unr), Real(9.1) * dum);
2651  morr_arr(i,j,k,MORRInd::umg) = std::min(morr_arr(i,j,k,MORRInd::umg), Real(20.) * dum);
2652  morr_arr(i,j,k,MORRInd::ung) = std::min(morr_arr(i,j,k,MORRInd::ung), Real(20.) * dum);
2653 
2654  // Set fall speed values
2655  morr_arr(i,j,k,MORRInd::fr) = morr_arr(i,j,k,MORRInd::umr); // RAIN FALL SPEED
2656  morr_arr(i,j,k,MORRInd::fi) = morr_arr(i,j,k,MORRInd::umi); // CLOUD ICE FALL SPEED
2657  morr_arr(i,j,k,MORRInd::fni) = morr_arr(i,j,k,MORRInd::uni); // CLOUD ICE NUMBER FALL SPEED
2658  morr_arr(i,j,k,MORRInd::fs) = morr_arr(i,j,k,MORRInd::ums); // SNOW FALL SPEED
2659  morr_arr(i,j,k,MORRInd::fns) = morr_arr(i,j,k,MORRInd::uns); // SNOW NUMBER FALL SPEED
2660  morr_arr(i,j,k,MORRInd::fnr) = morr_arr(i,j,k,MORRInd::unr); // RAIN NUMBER FALL SPEED
2661  morr_arr(i,j,k,MORRInd::fc) = morr_arr(i,j,k,MORRInd::umc); // CLOUD WATER FALL SPEED
2662  morr_arr(i,j,k,MORRInd::fnc) = morr_arr(i,j,k,MORRInd::unc); // CLOUD NUMBER FALL SPEED
2663  morr_arr(i,j,k,MORRInd::fg) = morr_arr(i,j,k,MORRInd::umg); // GRAUPEL FALL SPEED
2664  morr_arr(i,j,k,MORRInd::fng) = morr_arr(i,j,k,MORRInd::ung); // GRAUPEL NUMBER FALL SPEED
2665 
2666  // V3.3 MODIFY FALLSPEED BELOW LEVEL OF PRECIP
2667  if (morr_arr(i,j,k,MORRInd::fr) < Real(1.e-10)) {
2668  morr_arr(i,j,k,MORRInd::fr) = morr_arr(i,j,k+1,MORRInd::fr);
2669  }
2670  if (morr_arr(i,j,k,MORRInd::fi) < Real(1.e-10)) {
2671  morr_arr(i,j,k,MORRInd::fi) = morr_arr(i,j,k+1,MORRInd::fi);
2672  }
2673  if (morr_arr(i,j,k,MORRInd::fni) < Real(1.e-10)) {
2674  morr_arr(i,j,k,MORRInd::fni) = morr_arr(i,j,k+1,MORRInd::fni);
2675  }
2676  if (morr_arr(i,j,k,MORRInd::fs) < Real(1.e-10)) {
2677  morr_arr(i,j,k,MORRInd::fs) = morr_arr(i,j,k+1,MORRInd::fs);
2678  }
2679  if (morr_arr(i,j,k,MORRInd::fns) < Real(1.e-10)) {
2680  morr_arr(i,j,k,MORRInd::fns) = morr_arr(i,j,k+1,MORRInd::fns);
2681  }
2682  if (morr_arr(i,j,k,MORRInd::fnr) < Real(1.e-10)) {
2683  morr_arr(i,j,k,MORRInd::fnr) = morr_arr(i,j,k+1,MORRInd::fnr);
2684  }
2685  if (morr_arr(i,j,k,MORRInd::fc) < Real(1.e-10)) {
2686  morr_arr(i,j,k,MORRInd::fc) = morr_arr(i,j,k+1,MORRInd::fc);
2687  }
2688  if (morr_arr(i,j,k,MORRInd::fnc) < Real(1.e-10)) {
2689  morr_arr(i,j,k,MORRInd::fnc) = morr_arr(i,j,k+1,MORRInd::fnc);
2690  }
2691  if (morr_arr(i,j,k,MORRInd::fg) < Real(1.e-10)) {
2692  morr_arr(i,j,k,MORRInd::fg) = morr_arr(i,j,k+1,MORRInd::fg);
2693  }
2694  if (morr_arr(i,j,k,MORRInd::fng) < Real(1.e-10)) {
2695  morr_arr(i,j,k,MORRInd::fng) = morr_arr(i,j,k+1,MORRInd::fng);
2696  }
2697 
2698  // CALCULATE NUMBER OF SPLIT TIME STEPS
2699  // Find maximum fall speed at this point
2700  morr_arr(i,j,k,MORRInd::rgvm) = std::max({morr_arr(i,j,k,MORRInd::fr), morr_arr(i,j,k,MORRInd::fi), morr_arr(i,j,k,MORRInd::fs), morr_arr(i,j,k,MORRInd::fc),
2701  morr_arr(i,j,k,MORRInd::fni), morr_arr(i,j,k,MORRInd::fnr), morr_arr(i,j,k,MORRInd::fns), morr_arr(i,j,k,MORRInd::fnc),
2702  morr_arr(i,j,k,MORRInd::fg), morr_arr(i,j,k,MORRInd::fng)});
2703 
2704  // Calculate number of steps (dt and nstep would need to be defined elsewhere)
2705  nstep = std::max(static_cast<int>(morr_arr(i,j,k,MORRInd::rgvm) * dt / morr_arr(i,j,k,MORRInd::dzq) + one), nstep);
2706  // MULTIPLY VARIABLES BY RHO
2707  morr_arr(i,j,k,MORRInd::dumr) = morr_arr(i,j,k,MORRInd::dumr) * morr_arr(i,j,k,MORRInd::rho); // Rain water content * density
2708  morr_arr(i,j,k,MORRInd::dumi) = morr_arr(i,j,k,MORRInd::dumi) * morr_arr(i,j,k,MORRInd::rho); // Cloud ice content * density
2709  morr_arr(i,j,k,MORRInd::dumfni) = morr_arr(i,j,k,MORRInd::dumfni) * morr_arr(i,j,k,MORRInd::rho); // Cloud ice number * density
2710  morr_arr(i,j,k,MORRInd::dumqs) = morr_arr(i,j,k,MORRInd::dumqs) * morr_arr(i,j,k,MORRInd::rho); // Snow content * density
2711  morr_arr(i,j,k,MORRInd::dumfns) = morr_arr(i,j,k,MORRInd::dumfns) * morr_arr(i,j,k,MORRInd::rho); // Snow number * density
2712  morr_arr(i,j,k,MORRInd::dumfnr) = morr_arr(i,j,k,MORRInd::dumfnr) * morr_arr(i,j,k,MORRInd::rho); // Rain number * density
2713  morr_arr(i,j,k,MORRInd::dumc) = morr_arr(i,j,k,MORRInd::dumc) * morr_arr(i,j,k,MORRInd::rho); // Cloud water content * density
2714  morr_arr(i,j,k,MORRInd::dumfnc) = morr_arr(i,j,k,MORRInd::dumfnc) * morr_arr(i,j,k,MORRInd::rho); // Cloud droplet number * density
2715  morr_arr(i,j,k,MORRInd::dumg) = morr_arr(i,j,k,MORRInd::dumg) * morr_arr(i,j,k,MORRInd::rho); // Graupel content * density
2716  morr_arr(i,j,k,MORRInd::dumfng) = morr_arr(i,j,k,MORRInd::dumfng) * morr_arr(i,j,k,MORRInd::rho); // Graupel number * density
2717  } // k
2718 
2719  // Main time stepping loop for sedimentation
2720  for (int n = 1; n <= nstep; n++) {
2721  // Calculate initial fallout for each hydrometeor type for all levels
2722  for (int k = klo; k <= khi; k++) {
2723  morr_arr(i,j,k,MORRInd::faloutr) = morr_arr(i,j,k,MORRInd::fr) * morr_arr(i,j,k,MORRInd::dumr);
2724  morr_arr(i,j,k,MORRInd::falouti) = morr_arr(i,j,k,MORRInd::fi) * morr_arr(i,j,k,MORRInd::dumi);
2725  morr_arr(i,j,k,MORRInd::faloutni) = morr_arr(i,j,k,MORRInd::fni) * morr_arr(i,j,k,MORRInd::dumfni);
2726  morr_arr(i,j,k,MORRInd::falouts) = morr_arr(i,j,k,MORRInd::fs) * morr_arr(i,j,k,MORRInd::dumqs);
2727  morr_arr(i,j,k,MORRInd::faloutns) = morr_arr(i,j,k,MORRInd::fns) * morr_arr(i,j,k,MORRInd::dumfns);
2728  morr_arr(i,j,k,MORRInd::faloutnr) = morr_arr(i,j,k,MORRInd::fnr) * morr_arr(i,j,k,MORRInd::dumfnr);
2729  morr_arr(i,j,k,MORRInd::faloutc) = morr_arr(i,j,k,MORRInd::fc) * morr_arr(i,j,k,MORRInd::dumc);
2730  morr_arr(i,j,k,MORRInd::faloutnc) = morr_arr(i,j,k,MORRInd::fnc) * morr_arr(i,j,k,MORRInd::dumfnc);
2731  morr_arr(i,j,k,MORRInd::faloutg) = morr_arr(i,j,k,MORRInd::fg) * morr_arr(i,j,k,MORRInd::dumg);
2732  morr_arr(i,j,k,MORRInd::faloutng) = morr_arr(i,j,k,MORRInd::fng) * morr_arr(i,j,k,MORRInd::dumfng);
2733  } //k
2734 
2735  // Process top of model level
2736  int k = khi;
2737 
2738  // Calculate tendencies at top level
2739  morr_arr(i,j,k,MORRInd::faltndr) = morr_arr(i,j,k,MORRInd::faloutr) / morr_arr(i,j,k,MORRInd::dzq);
2740  morr_arr(i,j,k,MORRInd::faltndi) = morr_arr(i,j,k,MORRInd::falouti) / morr_arr(i,j,k,MORRInd::dzq);
2741  morr_arr(i,j,k,MORRInd::faltndni) = morr_arr(i,j,k,MORRInd::faloutni) / morr_arr(i,j,k,MORRInd::dzq);
2742  morr_arr(i,j,k,MORRInd::faltnds) = morr_arr(i,j,k,MORRInd::falouts) / morr_arr(i,j,k,MORRInd::dzq);
2743  morr_arr(i,j,k,MORRInd::faltndns) = morr_arr(i,j,k,MORRInd::faloutns) / morr_arr(i,j,k,MORRInd::dzq);
2744  morr_arr(i,j,k,MORRInd::faltndnr) = morr_arr(i,j,k,MORRInd::faloutnr) / morr_arr(i,j,k,MORRInd::dzq);
2745  morr_arr(i,j,k,MORRInd::faltndc) = morr_arr(i,j,k,MORRInd::faloutc) / morr_arr(i,j,k,MORRInd::dzq);
2746  morr_arr(i,j,k,MORRInd::faltndnc) = morr_arr(i,j,k,MORRInd::faloutnc) / morr_arr(i,j,k,MORRInd::dzq);
2747  morr_arr(i,j,k,MORRInd::faltndg) = morr_arr(i,j,k,MORRInd::faloutg) / morr_arr(i,j,k,MORRInd::dzq);
2748  morr_arr(i,j,k,MORRInd::faltndng) = morr_arr(i,j,k,MORRInd::faloutng) / morr_arr(i,j,k,MORRInd::dzq);
2749 
2750  // Add fallout terms to Eulerian tendencies (scaled by time step and density)
2751  morr_arr(i,j,k,MORRInd::qrsten) = morr_arr(i,j,k,MORRInd::qrsten) - morr_arr(i,j,k,MORRInd::faltndr) / nstep / morr_arr(i,j,k,MORRInd::rho);
2752  morr_arr(i,j,k,MORRInd::qisten) = morr_arr(i,j,k,MORRInd::qisten) - morr_arr(i,j,k,MORRInd::faltndi) / nstep / morr_arr(i,j,k,MORRInd::rho);
2753  morr_arr(i,j,k,MORRInd::ni3dten) = morr_arr(i,j,k,MORRInd::ni3dten) - morr_arr(i,j,k,MORRInd::faltndni) / nstep / morr_arr(i,j,k,MORRInd::rho);
2754  morr_arr(i,j,k,MORRInd::qnisten) = morr_arr(i,j,k,MORRInd::qnisten) - morr_arr(i,j,k,MORRInd::faltnds) / nstep / morr_arr(i,j,k,MORRInd::rho);
2755  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) - morr_arr(i,j,k,MORRInd::faltndns) / nstep / morr_arr(i,j,k,MORRInd::rho);
2756  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) - morr_arr(i,j,k,MORRInd::faltndnr) / nstep / morr_arr(i,j,k,MORRInd::rho);
2757  morr_arr(i,j,k,MORRInd::qcsten) = morr_arr(i,j,k,MORRInd::qcsten) - morr_arr(i,j,k,MORRInd::faltndc) / nstep / morr_arr(i,j,k,MORRInd::rho);
2758  morr_arr(i,j,k,MORRInd::nc3dten) = morr_arr(i,j,k,MORRInd::nc3dten) - morr_arr(i,j,k,MORRInd::faltndnc) / nstep / morr_arr(i,j,k,MORRInd::rho);
2759  morr_arr(i,j,k,MORRInd::qgsten) = morr_arr(i,j,k,MORRInd::qgsten) - morr_arr(i,j,k,MORRInd::faltndg) / nstep / morr_arr(i,j,k,MORRInd::rho);
2760  morr_arr(i,j,k,MORRInd::ng3dten) = morr_arr(i,j,k,MORRInd::ng3dten) - morr_arr(i,j,k,MORRInd::faltndng) / nstep / morr_arr(i,j,k,MORRInd::rho);
2761 
2762  // Update temporary working variables
2763  morr_arr(i,j,k,MORRInd::dumr) = morr_arr(i,j,k,MORRInd::dumr) - morr_arr(i,j,k,MORRInd::faltndr) * dt / nstep;
2764  morr_arr(i,j,k,MORRInd::dumi) = morr_arr(i,j,k,MORRInd::dumi) - morr_arr(i,j,k,MORRInd::faltndi) * dt / nstep;
2765  morr_arr(i,j,k,MORRInd::dumfni) = morr_arr(i,j,k,MORRInd::dumfni) - morr_arr(i,j,k,MORRInd::faltndni) * dt / nstep;
2766  morr_arr(i,j,k,MORRInd::dumqs) = morr_arr(i,j,k,MORRInd::dumqs) - morr_arr(i,j,k,MORRInd::faltnds) * dt / nstep;
2767  morr_arr(i,j,k,MORRInd::dumfns) = morr_arr(i,j,k,MORRInd::dumfns) - morr_arr(i,j,k,MORRInd::faltndns) * dt / nstep;
2768  morr_arr(i,j,k,MORRInd::dumfnr) = morr_arr(i,j,k,MORRInd::dumfnr) - morr_arr(i,j,k,MORRInd::faltndnr) * dt / nstep;
2769  morr_arr(i,j,k,MORRInd::dumc) = morr_arr(i,j,k,MORRInd::dumc) - morr_arr(i,j,k,MORRInd::faltndc) * dt / nstep;
2770  morr_arr(i,j,k,MORRInd::dumfnc) = morr_arr(i,j,k,MORRInd::dumfnc) - morr_arr(i,j,k,MORRInd::faltndnc) * dt / nstep;
2771  morr_arr(i,j,k,MORRInd::dumg) = morr_arr(i,j,k,MORRInd::dumg) - morr_arr(i,j,k,MORRInd::faltndg) * dt / nstep;
2772  morr_arr(i,j,k,MORRInd::dumfng) = morr_arr(i,j,k,MORRInd::dumfng) - morr_arr(i,j,k,MORRInd::faltndng) * dt / nstep;
2773 
2774  // Process remaining levels from top to bottom
2775  for (k = khi-1; k >= klo; k--) {
2776  // Calculate tendencies based on difference between levels
2777  morr_arr(i,j,k,MORRInd::faltndr) = (morr_arr(i,j,k+1,MORRInd::faloutr) - morr_arr(i,j,k,MORRInd::faloutr)) / morr_arr(i,j,k,MORRInd::dzq);
2778  morr_arr(i,j,k,MORRInd::faltndi) = (morr_arr(i,j,k+1,MORRInd::falouti) - morr_arr(i,j,k,MORRInd::falouti)) / morr_arr(i,j,k,MORRInd::dzq);
2779  morr_arr(i,j,k,MORRInd::faltndni) = (morr_arr(i,j,k+1,MORRInd::faloutni) - morr_arr(i,j,k,MORRInd::faloutni)) / morr_arr(i,j,k,MORRInd::dzq);
2780  morr_arr(i,j,k,MORRInd::faltnds) = (morr_arr(i,j,k+1,MORRInd::falouts) - morr_arr(i,j,k,MORRInd::falouts)) / morr_arr(i,j,k,MORRInd::dzq);
2781  morr_arr(i,j,k,MORRInd::faltndns) = (morr_arr(i,j,k+1,MORRInd::faloutns) - morr_arr(i,j,k,MORRInd::faloutns)) / morr_arr(i,j,k,MORRInd::dzq);
2782  morr_arr(i,j,k,MORRInd::faltndnr) = (morr_arr(i,j,k+1,MORRInd::faloutnr) - morr_arr(i,j,k,MORRInd::faloutnr)) / morr_arr(i,j,k,MORRInd::dzq);
2783  morr_arr(i,j,k,MORRInd::faltndc) = (morr_arr(i,j,k+1,MORRInd::faloutc) - morr_arr(i,j,k,MORRInd::faloutc)) / morr_arr(i,j,k,MORRInd::dzq);
2784  morr_arr(i,j,k,MORRInd::faltndnc) = (morr_arr(i,j,k+1,MORRInd::faloutnc) - morr_arr(i,j,k,MORRInd::faloutnc)) / morr_arr(i,j,k,MORRInd::dzq);
2785  morr_arr(i,j,k,MORRInd::faltndg) = (morr_arr(i,j,k+1,MORRInd::faloutg) - morr_arr(i,j,k,MORRInd::faloutg)) / morr_arr(i,j,k,MORRInd::dzq);
2786  morr_arr(i,j,k,MORRInd::faltndng) = (morr_arr(i,j,k+1,MORRInd::faloutng) - morr_arr(i,j,k,MORRInd::faloutng)) / morr_arr(i,j,k,MORRInd::dzq);
2787 
2788  // Add fallout terms to Eulerian tendencies (positive here, as mass flows in from above)
2789  morr_arr(i,j,k,MORRInd::qrsten) = morr_arr(i,j,k,MORRInd::qrsten) + morr_arr(i,j,k,MORRInd::faltndr) / nstep / morr_arr(i,j,k,MORRInd::rho);
2790  morr_arr(i,j,k,MORRInd::qisten) = morr_arr(i,j,k,MORRInd::qisten) + morr_arr(i,j,k,MORRInd::faltndi) / nstep / morr_arr(i,j,k,MORRInd::rho);
2791  morr_arr(i,j,k,MORRInd::ni3dten) = morr_arr(i,j,k,MORRInd::ni3dten) + morr_arr(i,j,k,MORRInd::faltndni) / nstep / morr_arr(i,j,k,MORRInd::rho);
2792  morr_arr(i,j,k,MORRInd::qnisten) = morr_arr(i,j,k,MORRInd::qnisten) + morr_arr(i,j,k,MORRInd::faltnds) / nstep / morr_arr(i,j,k,MORRInd::rho);
2793  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + morr_arr(i,j,k,MORRInd::faltndns) / nstep / morr_arr(i,j,k,MORRInd::rho);
2794  morr_arr(i,j,k,MORRInd::nr3dten) = morr_arr(i,j,k,MORRInd::nr3dten) + morr_arr(i,j,k,MORRInd::faltndnr) / nstep / morr_arr(i,j,k,MORRInd::rho);
2795  morr_arr(i,j,k,MORRInd::qcsten) = morr_arr(i,j,k,MORRInd::qcsten) + morr_arr(i,j,k,MORRInd::faltndc) / nstep / morr_arr(i,j,k,MORRInd::rho);
2796  morr_arr(i,j,k,MORRInd::nc3dten) = morr_arr(i,j,k,MORRInd::nc3dten) + morr_arr(i,j,k,MORRInd::faltndnc) / nstep / morr_arr(i,j,k,MORRInd::rho);
2797  morr_arr(i,j,k,MORRInd::qgsten) = morr_arr(i,j,k,MORRInd::qgsten) + morr_arr(i,j,k,MORRInd::faltndg) / nstep / morr_arr(i,j,k,MORRInd::rho);
2798  morr_arr(i,j,k,MORRInd::ng3dten) = morr_arr(i,j,k,MORRInd::ng3dten) + morr_arr(i,j,k,MORRInd::faltndng) / nstep / morr_arr(i,j,k,MORRInd::rho);
2799  // Update temporary working variables
2800  morr_arr(i,j,k,MORRInd::dumr) = morr_arr(i,j,k,MORRInd::dumr) + morr_arr(i,j,k,MORRInd::faltndr) * dt / nstep;
2801  morr_arr(i,j,k,MORRInd::dumi) = morr_arr(i,j,k,MORRInd::dumi) + morr_arr(i,j,k,MORRInd::faltndi) * dt / nstep;
2802  morr_arr(i,j,k,MORRInd::dumfni) = morr_arr(i,j,k,MORRInd::dumfni) + morr_arr(i,j,k,MORRInd::faltndni) * dt / nstep;
2803  morr_arr(i,j,k,MORRInd::dumqs) = morr_arr(i,j,k,MORRInd::dumqs) + morr_arr(i,j,k,MORRInd::faltnds) * dt / nstep;
2804  morr_arr(i,j,k,MORRInd::dumfns) = morr_arr(i,j,k,MORRInd::dumfns) + morr_arr(i,j,k,MORRInd::faltndns) * dt / nstep;
2805  morr_arr(i,j,k,MORRInd::dumfnr) = morr_arr(i,j,k,MORRInd::dumfnr) + morr_arr(i,j,k,MORRInd::faltndnr) * dt / nstep;
2806  morr_arr(i,j,k,MORRInd::dumc) = morr_arr(i,j,k,MORRInd::dumc) + morr_arr(i,j,k,MORRInd::faltndc) * dt / nstep;
2807  morr_arr(i,j,k,MORRInd::dumfnc) = morr_arr(i,j,k,MORRInd::dumfnc) + morr_arr(i,j,k,MORRInd::faltndnc) * dt / nstep;
2808  morr_arr(i,j,k,MORRInd::dumg) = morr_arr(i,j,k,MORRInd::dumg) + morr_arr(i,j,k,MORRInd::faltndg) * dt / nstep;
2809  morr_arr(i,j,k,MORRInd::dumfng) = morr_arr(i,j,k,MORRInd::dumfng) + morr_arr(i,j,k,MORRInd::faltndng) * dt / nstep;
2810  }
2811  // Get precipitation and snowfall accumulation during the time step
2812  // Factor of 1000 converts from m to mm, but division by density
2813  // of liquid water cancels this factor of 1000
2814  int kts=klo;
2815  morr_arr(i,j,klo,MORRInd::precrt) += (morr_arr(i,j,kts,MORRInd::faloutr) + morr_arr(i,j,kts,MORRInd::faloutc) + morr_arr(i,j,kts,MORRInd::falouts) +
2816  morr_arr(i,j,kts,MORRInd::falouti) + morr_arr(i,j,kts,MORRInd::faloutg)) * dt / nstep;
2817  morr_arr(i,j,klo,MORRInd::snowrt) += (morr_arr(i,j,kts,MORRInd::falouts) + morr_arr(i,j,kts,MORRInd::falouti) + morr_arr(i,j,kts,MORRInd::faloutg)) * dt / nstep;
2818 
2819  // Added 7/13/13
2820  morr_arr(i,j,klo,MORRInd::snowprt) += (morr_arr(i,j,kts,MORRInd::falouti) + morr_arr(i,j,kts,MORRInd::falouts)) * dt / nstep;
2821  morr_arr(i,j,klo,MORRInd::grplprt) += morr_arr(i,j,kts,MORRInd::faloutg) * dt / nstep;
2822  }
2823 
2824  for(int k=klo; k<=khi; k++) {
2825  Real evs; // EVS: Saturation vapor pressure
2826  Real eis; // EIS: Ice saturation vapor pressure
2827  Real qvs; // QVS: Saturation mixing ratio
2828  Real qvi; // QVI: Ice saturation mixing ratio
2829  Real qvqvs; // QVQVS: Saturation ratio
2830  Real qvqvsi; // QVQVSI: Ice saturation ratio
2831 
2832  // ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
2833  morr_arr(i,j,k,MORRInd::qr3dten) = morr_arr(i,j,k,MORRInd::qr3dten) + morr_arr(i,j,k,MORRInd::qrsten);
2834  morr_arr(i,j,k,MORRInd::qi3dten) = morr_arr(i,j,k,MORRInd::qi3dten) + morr_arr(i,j,k,MORRInd::qisten);
2835  morr_arr(i,j,k,MORRInd::qc3dten) = morr_arr(i,j,k,MORRInd::qc3dten) + morr_arr(i,j,k,MORRInd::qcsten);
2836  morr_arr(i,j,k,MORRInd::qg3dten) = morr_arr(i,j,k,MORRInd::qg3dten) + morr_arr(i,j,k,MORRInd::qgsten);
2837  morr_arr(i,j,k,MORRInd::qni3dten) = morr_arr(i,j,k,MORRInd::qni3dten) + morr_arr(i,j,k,MORRInd::qnisten);
2838  // PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
2839  // bug fix
2840  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall && morr_arr(i,j,k,MORRInd::t3d) < Real(273.15) && morr_arr(i,j,k,MORRInd::lami) >= Real(1.e-10)) {
2841  if (one/morr_arr(i,j,k,MORRInd::lami) >= Real(2)*m_dcs) {
2842  morr_arr(i,j,k,MORRInd::qni3dten) = morr_arr(i,j,k,MORRInd::qni3dten) + morr_arr(i,j,k,MORRInd::qi3d)/dt + morr_arr(i,j,k,MORRInd::qi3dten);
2843  morr_arr(i,j,k,MORRInd::ns3dten) = morr_arr(i,j,k,MORRInd::ns3dten) + morr_arr(i,j,k,MORRInd::ni3d)/dt + morr_arr(i,j,k,MORRInd::ni3dten);
2844  morr_arr(i,j,k,MORRInd::qi3dten) = -morr_arr(i,j,k,MORRInd::qi3d)/dt;
2845  morr_arr(i,j,k,MORRInd::ni3dten) = -morr_arr(i,j,k,MORRInd::ni3d)/dt;
2846  }
2847  }
2848 
2849  // Add tendencies to ensure consistency between mixing ratio and number concentration
2850  morr_arr(i,j,k,MORRInd::qc3d) = morr_arr(i,j,k,MORRInd::qc3d) + morr_arr(i,j,k,MORRInd::qc3dten)*dt;
2851  morr_arr(i,j,k,MORRInd::qi3d) = morr_arr(i,j,k,MORRInd::qi3d) + morr_arr(i,j,k,MORRInd::qi3dten)*dt;
2852  morr_arr(i,j,k,MORRInd::qni3d) = morr_arr(i,j,k,MORRInd::qni3d) + morr_arr(i,j,k,MORRInd::qni3dten)*dt;
2853  morr_arr(i,j,k,MORRInd::qr3d) = morr_arr(i,j,k,MORRInd::qr3d) + morr_arr(i,j,k,MORRInd::qr3dten)*dt;
2854  morr_arr(i,j,k,MORRInd::nc3d) = morr_arr(i,j,k,MORRInd::nc3d) + morr_arr(i,j,k,MORRInd::nc3dten)*dt;
2855  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::ni3d) + morr_arr(i,j,k,MORRInd::ni3dten)*dt;
2856  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::ns3d) + morr_arr(i,j,k,MORRInd::ns3dten)*dt;
2857  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::nr3d) + morr_arr(i,j,k,MORRInd::nr3dten)*dt;
2858  if (m_igraup == 0) {
2859  morr_arr(i,j,k,MORRInd::qg3d) = morr_arr(i,j,k,MORRInd::qg3d) + morr_arr(i,j,k,MORRInd::qg3dten)*dt;
2860  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::ng3d) + morr_arr(i,j,k,MORRInd::ng3dten)*dt;
2861  }
2862 
2863  // ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
2864  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) + morr_arr(i,j,k,MORRInd::t3dten)*dt;
2865  morr_arr(i,j,k,MORRInd::qv3d) = morr_arr(i,j,k,MORRInd::qv3d) + morr_arr(i,j,k,MORRInd::qv3dten)*dt;
2866  // SATURATION VAPOR PRESSURE AND MIXING RATIO
2867  // hm, add fix for low pressure, 5/12/10
2868  // Assuming POLYSVP is defined elsewhere
2869  evs = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(morr_arr(i,j,k,MORRInd::t3d), 0)); // PA
2870  eis = std::min(Real(0.99) * morr_arr(i,j,k,MORRInd::pres), calc_saturation_vapor_pressure(morr_arr(i,j,k,MORRInd::t3d), 1)); // PA
2871 
2872  // MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
2873  if (eis > evs) {
2874  eis = evs; // temporary update: adjust ice saturation pressure
2875  }
2876 
2877  // SATURATION MIXING RATIOS
2878  qvs = m_ep_2 * evs / (morr_arr(i,j,k,MORRInd::pres) - evs); // budget equation: calculate water saturation mixing ratio
2879  qvi = m_ep_2 * eis / (morr_arr(i,j,k,MORRInd::pres) - eis); // budget equation: calculate ice saturation mixing ratio
2880 
2881  // SATURATION RATIOS
2882  qvqvs = morr_arr(i,j,k,MORRInd::qv3d) / qvs; // budget equation: calculate water saturation ratio
2883  qvqvsi = morr_arr(i,j,k,MORRInd::qv3d) / qvi; // budget equation: calculate ice saturation ratio
2884  // AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
2885  if (qvqvs < Real(0.9)) {
2886  if (morr_arr(i,j,k,MORRInd::qr3d) < Real(1.0e-8)) {
2887  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qr3d);
2888  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm);
2889  morr_arr(i,j,k,MORRInd::qr3d) = Real(0);
2890  }
2891  if (morr_arr(i,j,k,MORRInd::qc3d) < Real(1.0e-8)) {
2892  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qc3d);
2893  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::xxlv) / morr_arr(i,j,k,MORRInd::cpm);
2894  morr_arr(i,j,k,MORRInd::qc3d) = Real(0);
2895  }
2896  }
2897  if (qvqvsi < Real(0.9)) {
2898  if (morr_arr(i,j,k,MORRInd::qi3d) < Real(1.0e-8)) {
2899  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qi3d);
2900  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm);
2901  morr_arr(i,j,k,MORRInd::qi3d) = Real(0);
2902  }
2903  if (morr_arr(i,j,k,MORRInd::qni3d) < Real(1.0e-8)) {
2904  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qni3d);
2905  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qni3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm);
2906  morr_arr(i,j,k,MORRInd::qni3d) = Real(0);
2907  }
2908  if (morr_arr(i,j,k,MORRInd::qg3d) < Real(1.0e-8)) {
2909  morr_arr(i,j,k,MORRInd::qv3d) += morr_arr(i,j,k,MORRInd::qg3d);
2910  morr_arr(i,j,k,MORRInd::t3d) -= morr_arr(i,j,k,MORRInd::qg3d) * morr_arr(i,j,k,MORRInd::xxls) / morr_arr(i,j,k,MORRInd::cpm);
2911  morr_arr(i,j,k,MORRInd::qg3d) = Real(0);
2912  }
2913  }
2914  // IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
2915  if (morr_arr(i,j,k,MORRInd::qc3d) < m_qsmall) {
2916  morr_arr(i,j,k,MORRInd::qc3d) = Real(0);
2917  morr_arr(i,j,k,MORRInd::nc3d) = Real(0);
2918  morr_arr(i,j,k,MORRInd::effc) = Real(0);
2919  }
2920  if (morr_arr(i,j,k,MORRInd::qr3d) < m_qsmall) {
2921  morr_arr(i,j,k,MORRInd::qr3d) = Real(0);
2922  morr_arr(i,j,k,MORRInd::nr3d) = Real(0);
2923  morr_arr(i,j,k,MORRInd::effr) = Real(0);
2924  }
2925  if (morr_arr(i,j,k,MORRInd::qi3d) < m_qsmall) {
2926  morr_arr(i,j,k,MORRInd::qi3d) = Real(0);
2927  morr_arr(i,j,k,MORRInd::ni3d) = Real(0);
2928  morr_arr(i,j,k,MORRInd::effi) = Real(0);
2929  }
2930  if (morr_arr(i,j,k,MORRInd::qni3d) < m_qsmall) {
2931  morr_arr(i,j,k,MORRInd::qni3d) = Real(0);
2932  morr_arr(i,j,k,MORRInd::ns3d) = Real(0);
2933  morr_arr(i,j,k,MORRInd::effs) = Real(0);
2934  }
2935  if (morr_arr(i,j,k,MORRInd::qg3d) < m_qsmall) {
2936  morr_arr(i,j,k,MORRInd::qg3d) = Real(0);
2937  morr_arr(i,j,k,MORRInd::ng3d) = Real(0);
2938  morr_arr(i,j,k,MORRInd::effg) = Real(0);
2939  }
2940  /*
2941  // Skip calculations if there is no cloud/precipitation water
2942  if ((morr_arr(i,j,k,MORRInd::qc3d) < m_qsmall && // CLOUD WATER MIXING RATIO (KG/KG)
2943  morr_arr(i,j,k,MORRInd::qi3d) < m_qsmall && // CLOUD ICE MIXING RATIO (KG/KG)
2944  morr_arr(i,j,k,MORRInd::qni3d) < m_qsmall && // SNOW MIXING RATIO (KG/KG)
2945  morr_arr(i,j,k,MORRInd::qr3d) < m_qsmall && // RAIN MIXING RATIO (KG/KG)
2946  morr_arr(i,j,k,MORRInd::qg3d) < m_qsmall)) { // GRAUPEL MIX RATIO (KG/KG)
2947  goto label_500;
2948  } else {*/
2949  if (!(morr_arr(i,j,k,MORRInd::qc3d) < m_qsmall && // CLOUD WATER MIXING RATIO (KG/KG)
2950  morr_arr(i,j,k,MORRInd::qi3d) < m_qsmall && // CLOUD ICE MIXING RATIO (KG/KG)
2951  morr_arr(i,j,k,MORRInd::qni3d) < m_qsmall && // SNOW MIXING RATIO (KG/KG)
2952  morr_arr(i,j,k,MORRInd::qr3d) < m_qsmall && // RAIN MIXING RATIO (KG/KG)
2953  morr_arr(i,j,k,MORRInd::qg3d) < m_qsmall)) { // GRAUPEL MIX RATIO (KG/KG)
2954  // CALCULATE INSTANTANEOUS PROCESSES
2955 
2956  // ADD MELTING OF CLOUD ICE TO FORM RAIN
2957  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall && morr_arr(i,j,k,MORRInd::t3d) >= Real(273.15)) {
2958  morr_arr(i,j,k,MORRInd::qr3d) = morr_arr(i,j,k,MORRInd::qr3d) + morr_arr(i,j,k,MORRInd::qi3d);
2959  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) - morr_arr(i,j,k,MORRInd::qi3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm);
2960  morr_arr(i,j,k,MORRInd::qi3d) = Real(0);
2961  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::nr3d) + morr_arr(i,j,k,MORRInd::ni3d);
2962  morr_arr(i,j,k,MORRInd::ni3d) = Real(0);
2963  }
2964  // ****SENSITIVITY - NO ICE
2965  if ((m_iliq != 1)) {
2966 
2967  // HOMOGENEOUS FREEZING OF CLOUD WATER
2968  if (morr_arr(i,j,k,MORRInd::t3d) <= Real(233.15) && morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
2969  morr_arr(i,j,k,MORRInd::qi3d) = morr_arr(i,j,k,MORRInd::qi3d) + morr_arr(i,j,k,MORRInd::qc3d);
2970  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) + morr_arr(i,j,k,MORRInd::qc3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm);
2971  morr_arr(i,j,k,MORRInd::qc3d) = Real(0);
2972  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::ni3d) + morr_arr(i,j,k,MORRInd::nc3d);
2973  morr_arr(i,j,k,MORRInd::nc3d) = Real(0);
2974  }
2975  // HOMOGENEOUS FREEZING OF RAIN
2976  if (m_igraup == 0) {
2977  if (morr_arr(i,j,k,MORRInd::t3d) <= Real(233.15) && morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
2978  morr_arr(i,j,k,MORRInd::qg3d) = morr_arr(i,j,k,MORRInd::qg3d) + morr_arr(i,j,k,MORRInd::qr3d);
2979  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) + morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm);
2980  morr_arr(i,j,k,MORRInd::qr3d) = Real(0);
2981  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::ng3d) + morr_arr(i,j,k,MORRInd::nr3d);
2982  morr_arr(i,j,k,MORRInd::nr3d) = Real(0);
2983  }
2984  } else if (m_igraup == 1) {
2985  if (morr_arr(i,j,k,MORRInd::t3d) <= Real(233.15) && morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
2986  morr_arr(i,j,k,MORRInd::qni3d) = morr_arr(i,j,k,MORRInd::qni3d) + morr_arr(i,j,k,MORRInd::qr3d);
2987  morr_arr(i,j,k,MORRInd::t3d) = morr_arr(i,j,k,MORRInd::t3d) + morr_arr(i,j,k,MORRInd::qr3d) * morr_arr(i,j,k,MORRInd::xlf) / morr_arr(i,j,k,MORRInd::cpm);
2988  morr_arr(i,j,k,MORRInd::qr3d) = Real(0);
2989  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::ns3d) + morr_arr(i,j,k,MORRInd::nr3d);
2990  morr_arr(i,j,k,MORRInd::nr3d) = Real(0);
2991  }
2992  }
2993 
2994  }/* else {
2995  Real dontdoanything=m_iliq;//printf("m_iliq: %d\n",m_iliq);// goto label_778;
2996  }*/
2997 
2998 // label_778:
2999  // MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3000  morr_arr(i,j,k,MORRInd::ni3d) = std::max(Real(0), morr_arr(i,j,k,MORRInd::ni3d));
3001  morr_arr(i,j,k,MORRInd::ns3d) = std::max(Real(0), morr_arr(i,j,k,MORRInd::ns3d));
3002  morr_arr(i,j,k,MORRInd::nc3d) = std::max(Real(0), morr_arr(i,j,k,MORRInd::nc3d));
3003  morr_arr(i,j,k,MORRInd::nr3d) = std::max(Real(0), morr_arr(i,j,k,MORRInd::nr3d));
3004  morr_arr(i,j,k,MORRInd::ng3d) = std::max(Real(0), morr_arr(i,j,k,MORRInd::ng3d));
3005 
3006  // CLOUD ICE
3007  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
3008  morr_arr(i,j,k,MORRInd::lami) = std::pow(m_cons12 * morr_arr(i,j,k,MORRInd::ni3d) / morr_arr(i,j,k,MORRInd::qi3d), one/m_di);
3009  // CHECK FOR SLOPE
3010  // ADJUST VARS
3011  if (morr_arr(i,j,k,MORRInd::lami) < m_lammini) {
3012  morr_arr(i,j,k,MORRInd::lami) = m_lammini;
3013  morr_arr(i,j,k,MORRInd::n0i) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lami)) * morr_arr(i,j,k,MORRInd::qi3d) / m_cons12;
3014  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::n0i) / morr_arr(i,j,k,MORRInd::lami);
3015  } else if (morr_arr(i,j,k,MORRInd::lami) > m_lammaxi) {
3016  morr_arr(i,j,k,MORRInd::lami) = m_lammaxi;
3017  morr_arr(i,j,k,MORRInd::n0i) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lami)) * morr_arr(i,j,k,MORRInd::qi3d) / m_cons12;
3018  morr_arr(i,j,k,MORRInd::ni3d) = morr_arr(i,j,k,MORRInd::n0i) / morr_arr(i,j,k,MORRInd::lami);
3019  }
3020  }
3021 
3022  // RAIN
3023  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
3024  morr_arr(i,j,k,MORRInd::lamr) = std::pow(m_pi * m_rhow * morr_arr(i,j,k,MORRInd::nr3d) / morr_arr(i,j,k,MORRInd::qr3d), one/three);
3025 
3026  // CHECK FOR SLOPE
3027  // ADJUST VARS
3028  if (morr_arr(i,j,k,MORRInd::lamr) < m_lamminr) {
3029  morr_arr(i,j,k,MORRInd::lamr) = m_lamminr;
3030  morr_arr(i,j,k,MORRInd::n0r) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
3031  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr);
3032  } else if (morr_arr(i,j,k,MORRInd::lamr) > m_lammaxr) {
3033  morr_arr(i,j,k,MORRInd::lamr) = m_lammaxr;
3034  morr_arr(i,j,k,MORRInd::n0r) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lamr)) * morr_arr(i,j,k,MORRInd::qr3d) / (m_pi * m_rhow);
3035  morr_arr(i,j,k,MORRInd::nr3d) = morr_arr(i,j,k,MORRInd::n0r) / morr_arr(i,j,k,MORRInd::lamr);
3036  }
3037  }
3038 
3039  // CLOUD DROPLETS
3040  // MARTIN ET AL. (1994) FORMULA FOR PGAM
3041  if (morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
3042  Real dum = morr_arr(i,j,k,MORRInd::pres) / (Real(287.15) * morr_arr(i,j,k,MORRInd::t3d));
3043  morr_arr(i,j,k,MORRInd::pgam) = Real(0.0005714) * (morr_arr(i,j,k,MORRInd::nc3d) / Real(1.0e6) * dum) + Real(0.2714);
3044  morr_arr(i,j,k,MORRInd::pgam) = one/(amrex::Math::powi<2>(morr_arr(i,j,k,MORRInd::pgam))) - one;
3045  morr_arr(i,j,k,MORRInd::pgam) = std::max(morr_arr(i,j,k,MORRInd::pgam), Real(2));
3046  morr_arr(i,j,k,MORRInd::pgam) = std::min(morr_arr(i,j,k,MORRInd::pgam), Real(10.0));
3047 
3048  // CALCULATE LAMC
3049  morr_arr(i,j,k,MORRInd::lamc) = std::pow(m_cons26 * morr_arr(i,j,k,MORRInd::nc3d) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0)) /
3050  (morr_arr(i,j,k,MORRInd::qc3d) * gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one)), one/three);
3051 
3052  // LAMMIN, 60 MICRON DIAMETER
3053  // LAMMAX, 1 MICRON
3054  Real lammin = (morr_arr(i,j,k,MORRInd::pgam) + one) / Real(60.0e-6);
3055  Real lammax = (morr_arr(i,j,k,MORRInd::pgam) + one) / Real(1.0e-6);
3056 
3057  if (morr_arr(i,j,k,MORRInd::lamc) < lammin) {
3058  morr_arr(i,j,k,MORRInd::lamc) = lammin;
3059  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three * std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
3060  std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one)) - std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0)))) / m_cons26;
3061  } else if (morr_arr(i,j,k,MORRInd::lamc) > lammax) {
3062  morr_arr(i,j,k,MORRInd::lamc) = lammax;
3063  morr_arr(i,j,k,MORRInd::nc3d) = std::exp(three * std::log(morr_arr(i,j,k,MORRInd::lamc)) + std::log(morr_arr(i,j,k,MORRInd::qc3d)) +
3064  std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + one)) - std::log(gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0)))) / m_cons26;
3065  }
3066  }
3067 
3068  // SNOW
3069  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
3070  morr_arr(i,j,k,MORRInd::lams) = std::pow(m_cons1 * morr_arr(i,j,k,MORRInd::ns3d) / morr_arr(i,j,k,MORRInd::qni3d), one/m_ds);
3071 
3072  // CHECK FOR SLOPE
3073  // ADJUST VARS
3074  if (morr_arr(i,j,k,MORRInd::lams) < m_lammins) {
3075  morr_arr(i,j,k,MORRInd::lams) = m_lammins;
3076  morr_arr(i,j,k,MORRInd::n0s) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lams)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
3077  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams);
3078  } else if (morr_arr(i,j,k,MORRInd::lams) > m_lammaxs) {
3079  morr_arr(i,j,k,MORRInd::lams) = m_lammaxs;
3080  morr_arr(i,j,k,MORRInd::n0s) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lams)) * morr_arr(i,j,k,MORRInd::qni3d) / m_cons1;
3081  morr_arr(i,j,k,MORRInd::ns3d) = morr_arr(i,j,k,MORRInd::n0s) / morr_arr(i,j,k,MORRInd::lams);
3082  }
3083  }
3084 
3085  // GRAUPEL
3086  if (morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
3087  morr_arr(i,j,k,MORRInd::lamg) = std::pow(m_cons2 * morr_arr(i,j,k,MORRInd::ng3d) / morr_arr(i,j,k,MORRInd::qg3d), one/m_dg);
3088 
3089  // CHECK FOR SLOPE
3090  // ADJUST VARS
3091  if (morr_arr(i,j,k,MORRInd::lamg) < m_lamming) {
3092  morr_arr(i,j,k,MORRInd::lamg) = m_lamming;
3093  morr_arr(i,j,k,MORRInd::n0g) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lamg)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
3094  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg);
3095  } else if (morr_arr(i,j,k,MORRInd::lamg) > m_lammaxg) {
3096  morr_arr(i,j,k,MORRInd::lamg) = m_lammaxg;
3097  morr_arr(i,j,k,MORRInd::n0g) = amrex::Math::powi<4>(morr_arr(i,j,k,MORRInd::lamg)) * morr_arr(i,j,k,MORRInd::qg3d) / m_cons2;
3098  morr_arr(i,j,k,MORRInd::ng3d) = morr_arr(i,j,k,MORRInd::n0g) / morr_arr(i,j,k,MORRInd::lamg);
3099  }
3100  }
3101  }
3102 
3103 // label_500:
3104  // CALCULATE EFFECTIVE RADIUS
3105  if (morr_arr(i,j,k,MORRInd::qi3d) >= m_qsmall) {
3106  morr_arr(i,j,k,MORRInd::effi) = three / morr_arr(i,j,k,MORRInd::lami) / Real(2) * Real(1.0e6);
3107  } else {
3108  morr_arr(i,j,k,MORRInd::effi) = Real(25.0);
3109  }
3110 
3111  if (morr_arr(i,j,k,MORRInd::qni3d) >= m_qsmall) {
3112  morr_arr(i,j,k,MORRInd::effs) = three / morr_arr(i,j,k,MORRInd::lams) / Real(2) * Real(1.0e6);
3113  } else {
3114  morr_arr(i,j,k,MORRInd::effs) = Real(25.0);
3115  }
3116 
3117  if (morr_arr(i,j,k,MORRInd::qr3d) >= m_qsmall) {
3118  morr_arr(i,j,k,MORRInd::effr) = three / morr_arr(i,j,k,MORRInd::lamr) / Real(2) * Real(1.0e6);
3119  } else {
3120  morr_arr(i,j,k,MORRInd::effr) = Real(25.0);
3121  }
3122 
3123  if (morr_arr(i,j,k,MORRInd::qc3d) >= m_qsmall) {
3124  morr_arr(i,j,k,MORRInd::effc) = gamma_function(morr_arr(i,j,k,MORRInd::pgam) + Real(4.0)) / gamma_function(morr_arr(i,j,k,MORRInd::pgam) + three) / morr_arr(i,j,k,MORRInd::lamc) / Real(2) * Real(1.0e6);
3125  } else {
3126  morr_arr(i,j,k,MORRInd::effc) = Real(25.0);
3127  }
3128 
3129  if (morr_arr(i,j,k,MORRInd::qg3d) >= m_qsmall) {
3130  morr_arr(i,j,k,MORRInd::effg) = three / morr_arr(i,j,k,MORRInd::lamg) / Real(2) * Real(1.0e6);
3131  } else {
3132  morr_arr(i,j,k,MORRInd::effg) = Real(25.0);
3133  }
3134 
3135  // HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3136  // TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3137  // OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3138  // HM, 12/28/12, LOWER MAXIMUM ICE CONCENTRATION TO ADDRESS PROBLEM
3139  // OF EXCESSIVE AND PERSISTENT ANVIL
3140  // NOTE: THIS MAY CHANGE/REDUCE SENSITIVITY TO AEROSOL/CCN CONCENTRATION
3141  morr_arr(i,j,k,MORRInd::ni3d) = std::min(morr_arr(i,j,k,MORRInd::ni3d), Real(0.3e6) / morr_arr(i,j,k,MORRInd::rho));
3142 
3143  // ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3144  if (iinum == 0 && m_iact == 2) {
3145  morr_arr(i,j,k,MORRInd::nc3d) = std::min(morr_arr(i,j,k,MORRInd::nc3d), (m_nanew1 + m_nanew2) / morr_arr(i,j,k,MORRInd::rho));
3146  }
3147 
3148  // SWITCH FOR CONSTANT DROPLET NUMBER
3149  if (iinum == 1) {
3150  // CHANGE NDCNST FROM CM-3 TO KG-1
3151  morr_arr(i,j,k,MORRInd::nc3d) = m_ndcnst * Real(1.0e6) / morr_arr(i,j,k,MORRInd::rho);
3152  }
3153  }
3154 
3155  }/* else {
3156  goto label_400;
3157  }
3158  label_400:*/
3159  //End of _micro
3160 
3161  if(use_morr_cpp_answer) {
3162  for (int k=klo; k<=khi; k++) {
3163  // Transfer 1D variables back to 3D arrays
3164  qcl_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qc3d);
3165  qci_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qi3d);
3166  qps_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qni3d);
3167  qpr_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qr3d);
3168  ni_arr(i,j,k) = morr_arr(i,j,k,MORRInd::ni3d);
3169  ns_arr(i,j,k) = morr_arr(i,j,k,MORRInd::ns3d);
3170  nr_arr(i,j,k) = morr_arr(i,j,k,MORRInd::nr3d);
3171  qpg_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qg3d);
3172  ng_arr(i,j,k) = morr_arr(i,j,k,MORRInd::ng3d);
3173 
3174  // Temperature and potential temperature conversion
3175  theta_arr(i,j,k) = morr_arr(i,j,k,MORRInd::t3d) / pii_arr(i,j,k); // Convert temp back to potential temp
3176  qv_arr(i,j,k) = morr_arr(i,j,k,MORRInd::qv3d);
3177 
3178  //Deleted wrf-check, effc, and precr type data as not used by ERF
3179  /*
3180  // NEED gpu-compatible summation for rain_accum, check SAM or Kessler for better example
3181  rain_accum_arr(i,j,k) = rain_accum_arr(i,j,k) + morr_arr(i,j,k,MORRInd::precrt);
3182  snow_accum_arr(i,j,k) = snow_accum_arr(i,j,k) + morr_arr(i,j,k,MORRInd::snowprt);
3183  graup_accum_arr(i,j,k) = graup_accum_arr(i,j,k) + morr_arr(i,j,k,MORRInd::grplprt);
3184  */
3185 
3186  rainncv_arr(i,j,0) = morr_arr(i,j,klo,MORRInd::precrt);
3187  snowncv_arr(i,j,0) = morr_arr(i,j,klo,MORRInd::snowprt);
3188  graupelncv_arr(i,j,0) = morr_arr(i,j,klo,MORRInd::grplprt);
3189  sr_arr(i,j,0) = morr_arr(i,j,klo,MORRInd::snowrt) / (morr_arr(i,j,klo,MORRInd::precrt) + Real(1.e-12));
3190  } // k
3191 
3192  // Update precipitation accumulation variables
3193  // These are outside the k-loop in the original code
3194  rain_accum_arr(i,j,klo) = rain_accum_arr(i,j,klo) + morr_arr(i,j,klo,MORRInd::precrt);
3195  snow_accum_arr(i,j,klo) = snow_accum_arr(i,j,klo) + morr_arr(i,j,klo,MORRInd::snowprt);
3196  graup_accum_arr(i,j,klo) = graup_accum_arr(i,j,klo) + morr_arr(i,j,klo,MORRInd::grplprt);
3197 
3198  } // cpp
3199  });
3200 
3201  }
3202 
3203  if (run_morr_fort) {
3204 #ifdef ERF_USE_MORR_FORT
3205 #include "ERF_Morrison_Advance_F.H"
3206 #endif
3207  } // run_morr_fort
3208 
3209  }
3210  }
constexpr amrex::Real three
Definition: ERF_Constants.H:9
constexpr amrex::Real one
Definition: ERF_Constants.H:7
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
ParmParse pp("prob")
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
auto qv_arr
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:210
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real gamma_function(Real x)
Definition: ERF_MorrisonGammaFunction.H:300
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_saturation_vapor_pressure(const amrex::Real T, const int type)
Definition: ERF_MorrisonVaporPressure.H:36
Arena * Arena_Used
Definition: ERF_Morrison_Advance_F.H:17
ParallelFor(grown_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
auto rho_arr
Definition: ERF_UpdateWSubsidence_SineMassFlux.H:3
amrex::Geometry m_geom
Definition: ERF_Morrison.H:180
amrex::MultiFab * m_z_phys_nd
Definition: ERF_Morrison.H:190
amrex::Real m_rdOcp
Definition: ERF_Morrison.H:183
amrex::Array< FabPtr, MicVar_Morr::NumVars > mic_fab_vars
Definition: ERF_Morrison.H:194
bool m_do_cond
Definition: ERF_Morrison.H:184
@ qisten
Definition: ERF_AdvanceMorrison.cpp:76
@ dumfnc
Definition: ERF_AdvanceMorrison.cpp:134
@ qi3dten
Definition: ERF_AdvanceMorrison.cpp:48
@ agn
Definition: ERF_AdvanceMorrison.cpp:97
@ pres
Definition: ERF_AdvanceMorrison.cpp:65
@ precrt
Definition: ERF_AdvanceMorrison.cpp:82
@ dumc
Definition: ERF_AdvanceMorrison.cpp:133
@ falouti
Definition: ERF_AdvanceMorrison.cpp:113
@ ni3d
Definition: ERF_AdvanceMorrison.cpp:58
@ dumi
Definition: ERF_AdvanceMorrison.cpp:98
@ faloutnc
Definition: ERF_AdvanceMorrison.cpp:141
@ snowprt
Definition: ERF_AdvanceMorrison.cpp:84
@ lamg
Definition: ERF_AdvanceMorrison.cpp:40
@ xxls
Definition: ERF_AdvanceMorrison.cpp:154
@ n0s
Definition: ERF_AdvanceMorrison.cpp:43
@ fns
Definition: ERF_AdvanceMorrison.cpp:123
@ grplprt
Definition: ERF_AdvanceMorrison.cpp:85
@ dumfns
Definition: ERF_AdvanceMorrison.cpp:119
@ qni3d
Definition: ERF_AdvanceMorrison.cpp:56
@ ung
Definition: ERF_AdvanceMorrison.cpp:137
@ faloutng
Definition: ERF_AdvanceMorrison.cpp:127
@ qc3dten
Definition: ERF_AdvanceMorrison.cpp:47
@ w3d
Definition: ERF_AdvanceMorrison.cpp:67
@ rgvm
Definition: ERF_AdvanceMorrison.cpp:111
@ arn
Definition: ERF_AdvanceMorrison.cpp:94
@ lami
Definition: ERF_AdvanceMorrison.cpp:37
@ qscu1d
Definition: ERF_AdvanceMorrison.cpp:80
@ qcsten
Definition: ERF_AdvanceMorrison.cpp:78
@ fc
Definition: ERF_AdvanceMorrison.cpp:139
@ fnr
Definition: ERF_AdvanceMorrison.cpp:148
@ pgam
Definition: ERF_AdvanceMorrison.cpp:46
@ qrcu1d
Definition: ERF_AdvanceMorrison.cpp:79
@ fng
Definition: ERF_AdvanceMorrison.cpp:110
@ ng3dten
Definition: ERF_AdvanceMorrison.cpp:71
@ dlami
Definition: ERF_AdvanceMorrison.cpp:151
@ faloutni
Definition: ERF_AdvanceMorrison.cpp:114
@ faloutnr
Definition: ERF_AdvanceMorrison.cpp:146
@ NumInds
Definition: ERF_AdvanceMorrison.cpp:158
@ dumfni
Definition: ERF_AdvanceMorrison.cpp:100
@ n0g
Definition: ERF_AdvanceMorrison.cpp:45
@ dlams
Definition: ERF_AdvanceMorrison.cpp:149
@ n0i
Definition: ERF_AdvanceMorrison.cpp:42
@ effs
Definition: ERF_AdvanceMorrison.cpp:88
@ faltndg
Definition: ERF_AdvanceMorrison.cpp:131
@ cpm
Definition: ERF_AdvanceMorrison.cpp:156
@ qr3dten
Definition: ERF_AdvanceMorrison.cpp:50
@ t3d
Definition: ERF_AdvanceMorrison.cpp:63
@ dumqs
Definition: ERF_AdvanceMorrison.cpp:118
@ qg3d
Definition: ERF_AdvanceMorrison.cpp:72
@ lamr
Definition: ERF_AdvanceMorrison.cpp:39
@ qr3d
Definition: ERF_AdvanceMorrison.cpp:57
@ nc3d
Definition: ERF_AdvanceMorrison.cpp:68
@ qg3dten
Definition: ERF_AdvanceMorrison.cpp:70
@ nr3dten
Definition: ERF_AdvanceMorrison.cpp:53
@ dzq
Definition: ERF_AdvanceMorrison.cpp:66
@ dumfnr
Definition: ERF_AdvanceMorrison.cpp:145
@ uns
Definition: ERF_AdvanceMorrison.cpp:121
@ faltndng
Definition: ERF_AdvanceMorrison.cpp:132
@ effc
Definition: ERF_AdvanceMorrison.cpp:86
@ qnisten
Definition: ERF_AdvanceMorrison.cpp:77
@ faloutg
Definition: ERF_AdvanceMorrison.cpp:126
@ t3dten
Definition: ERF_AdvanceMorrison.cpp:61
@ qv3d
Definition: ERF_AdvanceMorrison.cpp:64
@ ni3dten
Definition: ERF_AdvanceMorrison.cpp:51
@ uni
Definition: ERF_AdvanceMorrison.cpp:103
@ umi
Definition: ERF_AdvanceMorrison.cpp:104
@ qni3dten
Definition: ERF_AdvanceMorrison.cpp:49
@ faloutr
Definition: ERF_AdvanceMorrison.cpp:112
@ dumr
Definition: ERF_AdvanceMorrison.cpp:99
@ faloutns
Definition: ERF_AdvanceMorrison.cpp:125
@ effi
Definition: ERF_AdvanceMorrison.cpp:87
@ faltndni
Definition: ERF_AdvanceMorrison.cpp:117
@ unc
Definition: ERF_AdvanceMorrison.cpp:135
@ umc
Definition: ERF_AdvanceMorrison.cpp:136
@ qv3dten
Definition: ERF_AdvanceMorrison.cpp:62
@ faltndns
Definition: ERF_AdvanceMorrison.cpp:129
@ nc3dten
Definition: ERF_AdvanceMorrison.cpp:69
@ dumg
Definition: ERF_AdvanceMorrison.cpp:101
@ dlamc
Definition: ERF_AdvanceMorrison.cpp:152
@ rho
Definition: ERF_AdvanceMorrison.cpp:91
@ effr
Definition: ERF_AdvanceMorrison.cpp:89
@ faltndnc
Definition: ERF_AdvanceMorrison.cpp:143
@ xxlv
Definition: ERF_AdvanceMorrison.cpp:155
@ faltndc
Definition: ERF_AdvanceMorrison.cpp:142
@ faltndnr
Definition: ERF_AdvanceMorrison.cpp:147
@ fni
Definition: ERF_AdvanceMorrison.cpp:108
@ umr
Definition: ERF_AdvanceMorrison.cpp:105
@ faloutc
Definition: ERF_AdvanceMorrison.cpp:140
@ ain
Definition: ERF_AdvanceMorrison.cpp:93
@ effg
Definition: ERF_AdvanceMorrison.cpp:90
@ faltnds
Definition: ERF_AdvanceMorrison.cpp:128
@ xlf
Definition: ERF_AdvanceMorrison.cpp:157
@ asn
Definition: ERF_AdvanceMorrison.cpp:95
@ fr
Definition: ERF_AdvanceMorrison.cpp:106
@ fnc
Definition: ERF_AdvanceMorrison.cpp:144
@ fi
Definition: ERF_AdvanceMorrison.cpp:107
@ dumfng
Definition: ERF_AdvanceMorrison.cpp:102
@ faltndr
Definition: ERF_AdvanceMorrison.cpp:115
@ lams
Definition: ERF_AdvanceMorrison.cpp:38
@ fs
Definition: ERF_AdvanceMorrison.cpp:122
@ qicu1d
Definition: ERF_AdvanceMorrison.cpp:81
@ qrsten
Definition: ERF_AdvanceMorrison.cpp:75
@ qgsten
Definition: ERF_AdvanceMorrison.cpp:74
@ cdist1
Definition: ERF_AdvanceMorrison.cpp:41
@ dlamg
Definition: ERF_AdvanceMorrison.cpp:153
@ mu
Definition: ERF_AdvanceMorrison.cpp:92
@ lamc
Definition: ERF_AdvanceMorrison.cpp:36
@ ns3dten
Definition: ERF_AdvanceMorrison.cpp:52
@ ns3d
Definition: ERF_AdvanceMorrison.cpp:59
@ dlamr
Definition: ERF_AdvanceMorrison.cpp:150
@ ng3d
Definition: ERF_AdvanceMorrison.cpp:73
@ qi3d
Definition: ERF_AdvanceMorrison.cpp:55
@ acn
Definition: ERF_AdvanceMorrison.cpp:96
@ falouts
Definition: ERF_AdvanceMorrison.cpp:124
@ faltndi
Definition: ERF_AdvanceMorrison.cpp:116
@ fg
Definition: ERF_AdvanceMorrison.cpp:109
@ qc3d
Definition: ERF_AdvanceMorrison.cpp:54
@ umg
Definition: ERF_AdvanceMorrison.cpp:138
@ unr
Definition: ERF_AdvanceMorrison.cpp:130
@ ums
Definition: ERF_AdvanceMorrison.cpp:120
@ n0r
Definition: ERF_AdvanceMorrison.cpp:44
@ snowrt
Definition: ERF_AdvanceMorrison.cpp:83
@ nr3d
Definition: ERF_AdvanceMorrison.cpp:60
@ 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
@ theta
Definition: ERF_Morrison.H:28
@ ni
Definition: ERF_Morrison.H:46
@ ns
Definition: ERF_Morrison.H:47
@ omega
Definition: ERF_Morrison.H:53
@ qps
Definition: ERF_Morrison.H:40
@ 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
@ rain_accum
Definition: ERF_Morrison.H:50
@ snow_accum
Definition: ERF_Morrison.H:51
@ psacr
Definition: ERF_WSM6.H:207
@ praci
Definition: ERF_WSM6.H:203
@ psmlt
Definition: ERF_WSM6.H:213
@ pgmlt
Definition: ERF_WSM6.H:214
@ piacr
Definition: ERF_WSM6.H:201
@ pracs
Definition: ERF_WSM6.H:204
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
MoistureType moisture_type
Definition: ERF_DataStruct.H:1299
Here is the call graph for this function:

◆ Copy_Micro_to_State()

void Morrison::Copy_Micro_to_State ( amrex::MultiFab &  cons_in)
overridevirtual

Updates conserved and microphysics variables in the provided MultiFabs from the internal MultiFabs that store Microphysics module data.

Parameters
[out]consConserved variables
[out]qmoistqv, qc, qi, qr, qs, qg

Reimplemented from NullMoist.

18 {
19  // Get the temperature, density, theta, qt and qp from input
20  for ( MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
21  const auto& box3d = mfi.tilebox();
22 
23  auto states_arr = cons.array(mfi);
24 
25  auto rho_arr = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
26  auto theta_arr = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
27 
28  auto qv_arr = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
29  auto qc_arr = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
30  auto qi_arr = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
31 
32  auto qpr_arr = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
33  auto qps_arr = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
34  auto qpg_arr = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
35 
36  auto nc_arr = mic_fab_vars[MicVar_Morr::nc]->array(mfi);
37  auto ni_arr = mic_fab_vars[MicVar_Morr::ni]->array(mfi);
38  auto nr_arr = mic_fab_vars[MicVar_Morr::nr]->array(mfi);
39  auto ns_arr = mic_fab_vars[MicVar_Morr::ns]->array(mfi);
40  auto ng_arr = mic_fab_vars[MicVar_Morr::ng]->array(mfi);
41 
42  // get potential total density, temperature, qt, qp
43  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
44  {
45  states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
46 
47  states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*std::max(Real(0),qv_arr(i,j,k));
48  states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*std::max(Real(0),qc_arr(i,j,k));
49  states_arr(i,j,k,RhoQ3_comp) = rho_arr(i,j,k)*std::max(Real(0),qi_arr(i,j,k));
50 
51  states_arr(i,j,k,RhoQ4_comp) = rho_arr(i,j,k)*std::max(Real(0),qpr_arr(i,j,k));
52  states_arr(i,j,k,RhoQ5_comp) = rho_arr(i,j,k)*std::max(Real(0),qps_arr(i,j,k));
53  states_arr(i,j,k,RhoQ6_comp) = rho_arr(i,j,k)*std::max(Real(0),qpg_arr(i,j,k));
54 
55  states_arr(i,j,k,RhoQ7_comp) = rho_arr(i,j,k)*std::max(Real(0),nc_arr(i,j,k));
56  states_arr(i,j,k,RhoQ8_comp) = rho_arr(i,j,k)*std::max(Real(0),ni_arr(i,j,k));
57  states_arr(i,j,k,RhoQ9_comp) = rho_arr(i,j,k)*std::max(Real(0),nr_arr(i,j,k));
58  states_arr(i,j,k,RhoQ10_comp) = rho_arr(i,j,k)*std::max(Real(0),ns_arr(i,j,k));
59  states_arr(i,j,k,RhoQ11_comp) = rho_arr(i,j,k)*std::max(Real(0),ng_arr(i,j,k));
60  });
61  }
62 
63  // Fill interior ghost cells and periodic boundaries
64  cons.FillBoundary(m_geom.periodicity());
65 }
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:44
#define RhoQ11_comp
Definition: ERF_IndexDefines.H:52
#define RhoQ9_comp
Definition: ERF_IndexDefines.H:50
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:47
#define RhoQ8_comp
Definition: ERF_IndexDefines.H:49
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:46
#define RhoQ7_comp
Definition: ERF_IndexDefines.H:48
#define RhoQ10_comp
Definition: ERF_IndexDefines.H:51
@ cons
Definition: ERF_IndexDefines.H:174

Referenced by Update_State_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Copy_State_to_Micro()

void Morrison::Copy_State_to_Micro ( const amrex::MultiFab &  cons_in)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input

Reimplemented from NullMoist.

74 {
75  // Get the temperature, density, theta, qt and qp from input
76  for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
77  const auto& box3d = mfi.growntilebox();
78 
79  auto states_array = cons_in.array(mfi);
80 
81  // Non-precipitating
82  auto qv_array = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
83  auto qc_array = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
84  auto qi_array = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
85  auto qn_array = mic_fab_vars[MicVar_Morr::qn]->array(mfi);
86  auto qt_array = mic_fab_vars[MicVar_Morr::qt]->array(mfi);
87 
88  // Precipitating
89  auto qpr_array = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
90  auto qps_array = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
91  auto qpg_array = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
92  auto qp_array = mic_fab_vars[MicVar_Morr::qp]->array(mfi);
93 
94  auto nc_array = mic_fab_vars[MicVar_Morr::nc]->array(mfi);
95  auto ni_array = mic_fab_vars[MicVar_Morr::ni]->array(mfi);
96  auto nr_array = mic_fab_vars[MicVar_Morr::nr]->array(mfi);
97  auto ns_array = mic_fab_vars[MicVar_Morr::ns]->array(mfi);
98  auto ng_array = mic_fab_vars[MicVar_Morr::ng]->array(mfi);
99 
100  auto rho_array = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
101  auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
102  auto tabs_array = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
103  auto pres_array = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
104 
105  // Get pressure, theta, temperature, density, and qt, qp
106  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
107  {
108  rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
109  theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
110 
111  qv_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp));
112  qc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp));
113  qi_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp));
114  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
115  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
116 
117  qpr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ4_comp)/states_array(i,j,k,Rho_comp));
118  qps_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ5_comp)/states_array(i,j,k,Rho_comp));
119  qpg_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ6_comp)/states_array(i,j,k,Rho_comp));
120 
121  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
122 
123  nc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ7_comp) /states_array(i,j,k,Rho_comp));
124  ni_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ8_comp) /states_array(i,j,k,Rho_comp));
125  nr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ9_comp) /states_array(i,j,k,Rho_comp));
126  ns_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ10_comp)/states_array(i,j,k,Rho_comp));
127  ng_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ11_comp)/states_array(i,j,k,Rho_comp));
128 
129  tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
130  states_array(i,j,k,RhoTheta_comp),
131  qv_array(i,j,k));
132 
133  // NOTE: the Morrison Fortran version uses Pa not hPa so we don't divideby 100!
134  pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); // * Real(0.01);
135  });
136  }
137 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:46
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
#define Rho_comp
Definition: ERF_IndexDefines.H:36
@ tabs
Definition: ERF_Morrison.H:29
@ qp
Definition: ERF_Morrison.H:38
@ qn
Definition: ERF_Morrison.H:33
@ qt
Definition: ERF_Morrison.H:32

Referenced by Update_Micro_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Define()

void Morrison::Define ( SolverChoice sc)
inlineoverridevirtual

Reimplemented from NullMoist.

72  {
73  m_rdOcp = sc.rdOcp;
74  m_do_cond = (!sc.use_shoc);
75  }
amrex::Real rdOcp
Definition: ERF_DataStruct.H:1222
bool use_shoc
Definition: ERF_DataStruct.H:1253

◆ GetPlotVar() [1/2]

void Morrison::GetPlotVar ( const std::string &  a_name,
amrex::MultiFab &  a_mf 
) const
overridevirtual

Reimplemented from NullMoist.

Referenced by GetPlotVar().

Here is the caller graph for this function:

◆ GetPlotVar() [2/2]

void Morrison::GetPlotVar ( const std::string &  a_name,
amrex::MultiFab &  a_mf,
const int   
) const
inlineoverridevirtual

Reimplemented from NullMoist.

162  {
163  GetPlotVar(a_name, a_mf);
164  }
void GetPlotVar(const std::string &a_name, amrex::MultiFab &a_mf) const override
Here is the call graph for this function:

◆ GetPlotVarNames()

void Morrison::GetPlotVarNames ( amrex::Vector< std::string > &  a_vec) const
overridevirtual

Populate a vector with names of all available Morrison plot variables.

This function returns the names of all Morrison microphysics variables that can be written to plotfiles. These names correspond to the diagnostic fields defined in the plot_entries table.

Parameters
[out]a_vecVector to be populated with plot variable names
Note
The returned vector will contain 19 variable names, including both prognostic quantities (mixing ratios, number concentrations) and diagnostic quantities (temperature, pressure, derived moisture totals).

Reimplemented from NullMoist.

102 {
103  a_vec.clear();
104  a_vec.reserve(sizeof(plot_entries) / sizeof(plot_entries[0]));
105  for (const auto& entry : plot_entries) {
106  a_vec.emplace_back(entry.name);
107  }
108 }

◆ Init()

void Morrison::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 
)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input
[in]qc_inCloud variables input
[in,out]qv_inVapor variables input
[in]qi_inIce variables input
[in]gridsThe boxes on which we will evolve the solution
[in]geomGeometry associated with these MultiFabs and grids
[in]dt_advanceTimestep for the advance

Reimplemented from NullMoist.

27 {
28  [[maybe_unused]] amrex::Real dt = dt_advance;
29  m_geom = geom;
30 
31  m_z_phys_nd = z_phys_nd.get();
32  m_detJ_cc = detJ_cc.get();
33 
34  MicVarMap.resize(m_qmoist_size);
36 
37 #if defined(ERF_USE_MORR_FORT) && defined(AMREX_USE_GPU)
38  Arena* Arena_Used = The_Managed_Arena();
39 #else
40  Arena* Arena_Used = The_Arena();
41 #endif
42 
43  // initialize microphysics variables
44  for (auto ivar = 0; ivar < MicVar_Morr::NumVars; ++ivar) {
45  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
46  1, cons_in.nGrowVect(),
47  MFInfo().SetArena(Arena_Used));
48  mic_fab_vars[ivar]->setVal(0.);
49  }
50 
51 #ifdef ERF_USE_MORR_FORT
52  bool use_cpp;
53  amrex::ParmParse pp("erf");
54  pp.query("use_morr_cpp_answer", use_cpp);
55 
56  if (!use_cpp) {
57  MoistureType moisture_type;
58  pp.query_enum_case_insensitive("moisture_model",moisture_type);
59  int morr_noice = (moisture_type == MoistureType::Morrison_NoIce);
60  int morr_rimed_ice = 0; // This is used to set something called "ihail"
61  morr_two_moment_init_c(morr_rimed_ice, morr_noice);
62  }
63 #endif
64 }
void morr_two_moment_init_c(int morr_rimed_ice, int morr_noice)
amrex::Vector< int > MicVarMap
Definition: ERF_Morrison.H:177
amrex::MultiFab * m_detJ_cc
Definition: ERF_Morrison.H:191
int m_qmoist_size
Definition: ERF_Morrison.H:168
@ NumVars
Definition: ERF_Morrison.H:54
Here is the call graph for this function:

◆ Qmoist_Ptr()

amrex::MultiFab* Morrison::Qmoist_Ptr ( const int &  varIdx)
inlineoverridevirtual

Reimplemented from NullMoist.

121  {
123  return mic_fab_vars[MicVarMap[varIdx]].get();
124  }
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
Here is the call graph for this function:

◆ Qmoist_Restart_Vars()

void Morrison::Qmoist_Restart_Vars ( const SolverChoice ,
std::vector< int > &  a_idx,
std::vector< std::string > &  a_names 
) const
inlineoverridevirtual

Reimplemented from NullMoist.

139  {
140  a_idx.clear();
141  a_names.clear();
142 
143  // The ordering here needs to match that in
144  // MicVarMap = {MicVar_Morr::rain_accum, MicVar_Morr::snow_accum, MicVar_Morr::graup_accum};
145  //
146  a_idx.push_back(0); a_names.push_back("RainAccum");
147  a_idx.push_back(1); a_names.push_back("SnowAccum");
148  a_idx.push_back(2); a_names.push_back("GraupAccum");
149  }

◆ Qmoist_Size()

int Morrison::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

127 { return Morrison::m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int Morrison::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

int n_qstate_moist_numconc_size
Definition: ERF_Morrison.H:174

◆ Qstate_Moist_Size()

int Morrison::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

int n_qstate_moist_size
Definition: ERF_Morrison.H:171

◆ Set_dzmin()

void Morrison::Set_dzmin ( const amrex::Real  dz_min)
inlineoverridevirtual

Reimplemented from NullMoist.

89  {
90  m_dzmin = dz_min;
91  }
amrex::Real m_dzmin
Definition: ERF_Morrison.H:187

◆ Update_Micro_Vars()

void Morrison::Update_Micro_Vars ( amrex::MultiFab &  cons_in)
inlineoverridevirtual

Reimplemented from NullMoist.

103  {
104  this->Copy_State_to_Micro(cons_in);
105  }
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitMorrison.cpp:73
Here is the call graph for this function:

◆ Update_State_Vars()

void Morrison::Update_State_Vars ( amrex::MultiFab &  cons_in,
const amrex::MultiFab &   
)
inlineoverridevirtual

Reimplemented from NullMoist.

110  {
111  this->Copy_Micro_to_State(cons_in);
112  }
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateMorrison.cpp:17
Here is the call graph for this function:

Member Data Documentation

◆ m_detJ_cc

amrex::MultiFab* Morrison::m_detJ_cc
private

◆ m_do_cond

bool Morrison::m_do_cond
private

Referenced by Define().

◆ m_dzmin

amrex::Real Morrison::m_dzmin
private

Referenced by Set_dzmin().

◆ m_geom

amrex::Geometry Morrison::m_geom
private

◆ m_qmoist_size

int Morrison::m_qmoist_size = 3
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

amrex::Real Morrison::m_rdOcp
private

Referenced by Define().

◆ m_z_phys_nd

amrex::MultiFab* Morrison::m_z_phys_nd
private

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar_Morr::NumVars> Morrison::mic_fab_vars
private

Referenced by Qmoist_Ptr().

◆ MicVarMap

amrex::Vector<int> Morrison::MicVarMap
private

Referenced by Qmoist_Ptr().

◆ n_qstate_moist_numconc_size

int Morrison::n_qstate_moist_numconc_size = 5
private

◆ n_qstate_moist_size

int Morrison::n_qstate_moist_size = 11
private

Referenced by Qstate_Moist_Size().


The documentation for this class was generated from the following files: