ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Public Member Functions | Static Public Member Functions | Private Types | Private Attributes | Static Private Attributes | List of all members
Morrison Class Reference

#include <ERF_Morrison.H>

Inheritance diagram for Morrison:
Inheritance graph
[legend]
Collaboration diagram for Morrison:
Collaboration graph
[legend]

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 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) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &sc) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
void Compute_Coefficients ()
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 

Static Public Member Functions

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat (int &i, int &j, int &k, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &fac_sub, const amrex::Real &an, const amrex::Real &bn, const amrex::Array4< amrex::Real > &tabs_array, const amrex::Array4< amrex::Real > &pres_array, const amrex::Array4< amrex::Real > &qv_array, const amrex::Array4< amrex::Real > &qc_array, const amrex::Array4< amrex::Real > &qi_array, const amrex::Array4< amrex::Real > &qn_array, const amrex::Array4< amrex::Real > &qt_array)
 

Private Types

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

Private Attributes

int m_qmoist_size = 8
 
int n_qstate_moist_size = 6
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::BoxArray m_gtoe
 
int nlev
 
int zlo
 
int zhi
 
int m_axis
 
amrex::Real m_fac_cond
 
amrex::Real m_fac_fus
 
amrex::Real m_fac_sub
 
amrex::Real m_gOcp
 
amrex::Real m_rdOcp
 
amrex::MultiFab * m_z_phys_nd
 
amrex::MultiFab * m_detJ_cc
 
amrex::Array< FabPtr, MicVar_Morr::NumVars > mic_fab_vars
 
amrex::TableData< amrex::Real, 1 > accrrc
 
amrex::TableData< amrex::Real, 1 > accrsi
 
amrex::TableData< amrex::Real, 1 > accrsc
 
amrex::TableData< amrex::Real, 1 > coefice
 
amrex::TableData< amrex::Real, 1 > evaps1
 
amrex::TableData< amrex::Real, 1 > evaps2
 
amrex::TableData< amrex::Real, 1 > accrgi
 
amrex::TableData< amrex::Real, 1 > accrgc
 
amrex::TableData< amrex::Real, 1 > evapg1
 
amrex::TableData< amrex::Real, 1 > evapg2
 
amrex::TableData< amrex::Real, 1 > evapr1
 
amrex::TableData< amrex::Real, 1 > evapr2
 
amrex::TableData< amrex::Real, 1 > rho1d
 
amrex::TableData< amrex::Real, 1 > pres1d
 
amrex::TableData< amrex::Real, 1 > tabs1d
 
amrex::TableData< amrex::Real, 1 > qt1d
 
amrex::TableData< amrex::Real, 1 > qv1d
 
amrex::TableData< amrex::Real, 1 > qn1d
 
amrex::TableData< amrex::Real, 1 > gamaz
 
amrex::TableData< amrex::Real, 1 > zmid
 

Static Private Attributes

static constexpr amrex::Real CFL_MAX = 0.5
 

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.

374  {
375  // Store timestep
376  amrex::Real dt = dt_advance;
377 
378  // Loop through the grids
379  for (amrex::MFIter mfi(*mic_fab_vars[MicVar_Morr::qcl],TileNoZ()); mfi.isValid(); ++mfi)
380  {
381  const amrex::Box& box = mfi.tilebox();
382 
383  // Get array data from class member variables
384  auto const& theta_arr = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
385  auto const& qv_arr = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
386  auto const& qcl_arr = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
387  auto const& qpr_arr = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
388  auto const& qci_arr = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
389  auto const& qps_arr = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
390  auto const& qpg_arr = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
391  auto const& ni_arr = mic_fab_vars[MicVar_Morr::ni]->array(mfi);
392  [[maybe_unused]] auto const& nc_arr = mic_fab_vars[MicVar_Morr::nc]->array(mfi);
393  auto const& ns_arr = mic_fab_vars[MicVar_Morr::ns]->array(mfi);
394  auto const& nr_arr = mic_fab_vars[MicVar_Morr::nr]->array(mfi);
395  auto const& ng_arr = mic_fab_vars[MicVar_Morr::ng]->array(mfi);
396  [[maybe_unused]] auto const& rho_arr = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
397  auto const& pres_arr = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
398  [[maybe_unused]] auto const& tabs_arr = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
399  auto const& rain_accum_arr = mic_fab_vars[MicVar_Morr::rain_accum]->array(mfi);
400  auto const& snow_accum_arr = mic_fab_vars[MicVar_Morr::snow_accum]->array(mfi);
401  auto const& graup_accum_arr = mic_fab_vars[MicVar_Morr::graup_accum]->array(mfi);
402  auto const& w_arr = mic_fab_vars[MicVar_Morr::omega]->array(mfi);
403 
404  // Get radar reflectivity array if radar diagnostics enabled
405  // auto const& refl_arr = m_do_radar_ref ? m_radar->array(mfi) : nullptr;
406  // auto const& refl_arr = m_radar->array(mfi);
407 
408  // Extract box dimensions
409  const int ilo = box.loVect()[0];
410  const int ihi = box.hiVect()[0];
411  const int jlo = box.loVect()[1];
412  const int jhi = box.hiVect()[1];
413  const int klo = box.loVect()[2];
414  const int khi = box.hiVect()[2];
415 
416  amrex::Box grown_box(box); grown_box.grow(3);
417 #ifdef ERF_USE_MORR_FORT
418  const int ilom = grown_box.loVect()[0];
419  const int ihim = grown_box.hiVect()[0];
420  const int jlom = grown_box.loVect()[1];
421  const int jhim = grown_box.hiVect()[1];
422  const int klom = grown_box.loVect()[2];
423  const int khim = grown_box.hiVect()[2];
424 #endif
425  // Calculate Exner function (PII) to convert potential temperature to temperature
426  // PII = (P/P0)^(R/cp)
427  FArrayBox pii_fab(grown_box, 1);
428  auto const& pii_arr = pii_fab.array();
429 
430  const amrex::Real p0 = 100000.0; // Reference pressure (Pa)
431 
432  const amrex::Real rdcp = m_rdOcp; // R/cp ratio
433 
434  // Calculate Exner function
435  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
436  // NOTE: the Morrison Fortran version uses Pa not hPa so we didn't divide p by 100
437  // so we don't need to multiply by 100 here
438  pii_arr(i,j,k) = std::pow((pres_arr(i,j,k)) / p0, rdcp);
439  });
440 
441  // Create arrays for height differences (dz)
442  FArrayBox dz_fab(grown_box, 1);
443  auto const& dz_arr = dz_fab.array();
444 
445  // Calculate height differences
446  const amrex::Real dz_val = m_geom.CellSize(m_axis);
447  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
448  dz_arr(i,j,k) = dz_val;
449  });
450  amrex::Box grown_boxD(grown_box); grown_boxD.makeSlab(2,0);
451 
452  // Arrays to store precipitation rates
453  FArrayBox rainncv_fab(grown_boxD, 1);
454  FArrayBox sr_fab(grown_boxD, 1); // Ratio of snow to total precipitation
455  FArrayBox snowncv_fab(grown_boxD, 1);
456  FArrayBox graupelncv_fab(grown_boxD, 1);
457 
458  auto const& rainncv_arr = rainncv_fab.array();
459  auto const& sr_arr = sr_fab.array();
460  auto const& snowncv_arr = snowncv_fab.array();
461  auto const& graupelncv_arr = graupelncv_fab.array();
462 
463  // Initialize precipitation rate arrays to zero
464  ParallelFor(grown_boxD, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
465  rainncv_arr(i,j,k) = 0.0;
466  sr_arr(i,j,k) = 0.0;
467  snowncv_arr(i,j,k) = 0.0;
468  graupelncv_arr(i,j,k) = 0.0;
469  });
470 
471  // Create terrain height array (not actually used by Morrison scheme)
472  FArrayBox ht_fab(amrex::Box(amrex::IntVect(ilo, jlo, 0), amrex::IntVect(ihi, jhi, 0)), 1);
473  [[maybe_unused]] auto const& ht_arr = ht_fab.array();
474  ParallelFor(amrex::Box(amrex::IntVect(ilo, jlo, 0), amrex::IntVect(ihi, jhi, 0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
475  ht_arr(i,j,k) = 0.0; // Not used by Morrison scheme
476  });
477 
478  // Create dummy arrays for cumulus tendencies (if needed)
479  FArrayBox qrcuten_fab(grown_box, 1);
480  FArrayBox qscuten_fab(grown_box, 1);
481  FArrayBox qicuten_fab(grown_box, 1);
482  auto const& qrcuten_arr = qrcuten_fab.array();
483  auto const& qscuten_arr = qscuten_fab.array();
484  auto const& qicuten_arr = qicuten_fab.array();
485 
486  // Initialize tendencies to zero (no cumulus parameterization in this example)
487  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
488  qrcuten_arr(i,j,k) = 0.0;
489  qscuten_arr(i,j,k) = 0.0;
490  qicuten_arr(i,j,k) = 0.0;
491  });
492 
493 #ifdef ERF_USE_MORR_FORT
494  // WRF-Chem related variables (optional)
495  bool flag_qndrop = false; // Flag to indicate droplet number prediction
496 
497  // Now create arrays for other optional variables
498  FArrayBox rainprod_fab(grown_box, 1);
499  FArrayBox evapprod_fab(grown_box, 1);
500  FArrayBox qlsink_fab(grown_box, 1);
501  FArrayBox precr_fab(grown_box, 1);
502  FArrayBox preci_fab(grown_box, 1);
503  FArrayBox precs_fab(grown_box, 1);
504  FArrayBox precg_fab(grown_box, 1);
505 
506  auto const& rainprod_arr = rainprod_fab.array();
507  auto const& evapprod_arr = evapprod_fab.array();
508  auto const& qlsink_arr = qlsink_fab.array();
509  auto const& precr_arr = precr_fab.array();
510  auto const& preci_arr = preci_fab.array();
511  auto const& precs_arr = precs_fab.array();
512  auto const& precg_arr = precg_fab.array();
513 
514  // Initialize WRF-Chem arrays to zero
515  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
516  rainprod_arr(i,j,k) = 0.0;
517  evapprod_arr(i,j,k) = 0.0;
518  qlsink_arr(i,j,k) = 0.0;
519  precr_arr(i,j,k) = 0.0;
520  preci_arr(i,j,k) = 0.0;
521  precs_arr(i,j,k) = 0.0;
522  precg_arr(i,j,k) = 0.0;
523  });
524 #endif
525 
526  // Create FArrayBox for slope parameters and PSD variables
527  FArrayBox lamc_fab(grown_box, 1);
528  FArrayBox lami_fab(grown_box, 1);
529  FArrayBox lams_fab(grown_box, 1);
530  FArrayBox lamr_fab(grown_box, 1);
531  FArrayBox lamg_fab(grown_box, 1);
532  FArrayBox cdist1_fab(grown_box, 1);
533  FArrayBox n0i_fab(grown_box, 1);
534  FArrayBox n0s_fab(grown_box, 1);
535  FArrayBox n0r_fab(grown_box, 1);
536  FArrayBox n0g_fab(grown_box, 1);
537  FArrayBox pgam_fab(grown_box, 1);
538 
539  // Get Array4 objects for each parameter
540  auto const& lamc = lamc_fab.array();
541  auto const& lami = lami_fab.array();
542  auto const& lams = lams_fab.array();
543  auto const& lamr = lamr_fab.array();
544  auto const& lamg = lamg_fab.array();
545  auto const& cdist1 = cdist1_fab.array();
546  auto const& n0i = n0i_fab.array();
547  auto const& n0s = n0s_fab.array();
548  auto const& n0r = n0r_fab.array();
549  auto const& n0g = n0g_fab.array();
550  auto const& pgam = pgam_fab.array();
551 
552  // Initialize all values to zero
553  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
554  lamc(i,j,k) = 0.0;
555  lami(i,j,k) = 0.0;
556  lams(i,j,k) = 0.0;
557  lamr(i,j,k) = 0.0;
558  lamg(i,j,k) = 0.0;
559  cdist1(i,j,k) = 0.0;
560  n0i(i,j,k) = 0.0;
561  n0s(i,j,k) = 0.0;
562  n0r(i,j,k) = 0.0;
563  n0g(i,j,k) = 0.0;
564  pgam(i,j,k) = 0.0;
565  });
566 
567 #ifdef ERF_USE_MORR_FORT
568  // Prepare data pointers for Fortran call
569  // These would be passed directly to the Fortran interface
570  double dummy_reflectivity = 0.0;
571  double* dummy_reflectivity_ptr = &dummy_reflectivity;
572 #endif
573  // Example call (pseudo-code - actual interface would depend on your Fortran interop setup)
574 
575  // Microphysics options/switches
576  int m_iact = 2; // CCN activation option (1: power-law, 2: lognormal aerosol)
577  int m_inum = 1; // Droplet number option (0: predict, 1: constant)
578 
579  int m_iliq = 0; // Liquid-only option (0: include ice, 1: liquid only)
580  int m_inuc = 0; // Ice nucleation option (0: mid-latitude, 1: arctic)
581  [[maybe_unused]] int m_ibase = 2; // Cloud base activation option
582  [[maybe_unused]] int m_isub = 0; // Sub-grid vertical velocity option
583  int m_igraup = 0; // Graupel option (0: include graupel, 1: no graupel)
584  int m_ihail = 0; // Graupel/hail option (0: graupel, 1: hail)
585 
586  if(sc.moisture_type == MoistureType::Morrison_NoIce) {
587  m_iliq = 1; // Liquid-only option (0: include ice, 1: liquid only)
588  m_inuc = 0; // Ice nucleation option (0: mid-latitude, 1: arctic)
589  m_ibase = 2; // Cloud base activation option
590  m_isub = 0; // Sub-grid vertical velocity option
591  m_igraup = 1; // Graupel option (0: include graupel, 1: no graupel)
592  m_ihail = 0; // Graupel/hail option (0: graupel, 1: hail)
593  }
594  [[maybe_unused]] bool m_do_radar_ref = false; // Radar reflectivity calculation flag
595 
596  // Physical constants
597  amrex::Real m_pi; // Pi constant
598  amrex::Real m_R; // Gas constant for dry air (J/kg/K)
599  amrex::Real m_Rd; // Gas constant for dry air (J/kg/K)
600  amrex::Real m_Rv; // Gas constant for water vapor (J/kg/K)
601  [[maybe_unused]] amrex::Real m_cp; // Specific heat at constant pressure (J/kg/K)
602  amrex::Real m_g; // Gravitational acceleration (m/s^2)
603  amrex::Real m_ep_2; // Molecular weight ratio (Rd/Rv)
604 
605  // Reference density and species densities
606  amrex::Real m_rhosu; // Standard air density at 850 mb (kg/m^3)
607  amrex::Real m_rhow; // Density of liquid water (kg/m^3)
608  amrex::Real m_rhoi; // Bulk density of cloud ice (kg/m^3)
609  amrex::Real m_rhosn; // Bulk density of snow (kg/m^3)
610  amrex::Real m_rhog; // Bulk density of graupel/hail (kg/m^3)
611 
612  // Fall speed parameters (V=AD^B)
613  amrex::Real m_ai, m_bi; // Cloud ice fall speed parameters
614  [[maybe_unused]] amrex::Real m_ac, m_bc; // Cloud droplet fall speed parameters
615  amrex::Real m_as, m_bs; // Snow fall speed parameters
616  amrex::Real m_ar, m_br; // Rain fall speed parameters
617  amrex::Real m_ag, m_bg; // Graupel/hail fall speed parameters
618 
619  // Microphysical parameters
620  amrex::Real m_aimm; // Parameter in Bigg immersion freezing
621  amrex::Real m_bimm; // Parameter in Bigg immersion freezing
622  amrex::Real m_ecr; // Collection efficiency between droplets/rain and snow/rain
623  amrex::Real m_dcs; // Threshold size for cloud ice autoconversion (m)
624  amrex::Real m_mi0; // Initial mass of nucleated ice crystal (kg)
625  amrex::Real m_mg0; // Mass of embryo graupel (kg)
626  amrex::Real m_f1s; // Ventilation parameter for snow
627  amrex::Real m_f2s; // Ventilation parameter for snow
628  amrex::Real m_f1r; // Ventilation parameter for rain
629  amrex::Real m_f2r; // Ventilation parameter for rain
630  amrex::Real m_qsmall; // Smallest allowed hydrometeor mixing ratio
631  amrex::Real m_eii; // Collection efficiency, ice-ice collisions
632  amrex::Real m_eci; // Collection efficiency, ice-droplet collisions
633  amrex::Real m_cpw; // Specific heat of liquid water (J/kg/K)
634  amrex::Real m_rin; // Radius of contact nuclei (m)
635  amrex::Real m_mmult; // Mass of splintered ice particle (kg)
636 
637  // Size distribution parameters
638  amrex::Real m_ci, m_di; // Cloud ice size distribution parameters
639  amrex::Real m_cs, m_ds; // Snow size distribution parameters
640  amrex::Real m_cg, m_dg; // Graupel size distribution parameters
641 
642  // Lambda limits for size distributions
643  amrex::Real m_lammaxi, m_lammini; // Cloud ice lambda limits
644  amrex::Real m_lammaxr, m_lamminr; // Rain lambda limits
645  amrex::Real m_lammaxs, m_lammins; // Snow lambda limits
646  amrex::Real m_lammaxg, m_lamming; // Graupel lambda limits
647 
648  // Constant droplet concentration (if INUM = 1)
649  amrex::Real m_ndcnst = 250.0; // Droplet number concentration (cm^-3)
650 
651  // CCN spectra parameters (for IACT = 1)
652  [[maybe_unused]] amrex::Real m_k1; // Exponent in CCN activation formula
653  [[maybe_unused]] amrex::Real m_c1; // Coefficient in CCN activation formula (cm^-3)
654 
655  // Aerosol activation parameters (for IACT = 2)
656  [[maybe_unused]] amrex::Real m_mw; // Molecular weight water (kg/mol)
657  [[maybe_unused]] amrex::Real m_osm; // Osmotic coefficient
658  [[maybe_unused]] amrex::Real m_vi; // Number of ions dissociated in solution
659  [[maybe_unused]] amrex::Real m_epsm; // Aerosol soluble fraction
660  [[maybe_unused]] amrex::Real m_rhoa; // Aerosol bulk density (kg/m^3)
661  [[maybe_unused]] amrex::Real m_map; // Molecular weight aerosol (kg/mol)
662  [[maybe_unused]] amrex::Real m_ma; // Molecular weight of air (kg/mol)
663  [[maybe_unused]] amrex::Real m_rr; // Universal gas constant (J/mol/K)
664  [[maybe_unused]] amrex::Real m_bact; // Activation parameter
665  [[maybe_unused]] amrex::Real m_rm1; // Geometric mean radius, mode 1 (m)
666  [[maybe_unused]] amrex::Real m_rm2; // Geometric mean radius, mode 2 (m)
667  amrex::Real m_nanew1; // Total aerosol concentration, mode 1 (m^-3)
668  amrex::Real m_nanew2; // Total aerosol concentration, mode 2 (m^-3)
669  [[maybe_unused]] amrex::Real m_sig1; // Standard deviation of aerosol dist, mode 1
670  [[maybe_unused]] amrex::Real m_sig2; // Standard deviation of aerosol dist, mode 2
671  [[maybe_unused]] amrex::Real m_f11; // Correction factor for activation, mode 1
672  [[maybe_unused]] amrex::Real m_f12; // Correction factor for activation, mode 1
673  [[maybe_unused]] amrex::Real m_f21; // Correction factor for activation, mode 2
674  [[maybe_unused]] amrex::Real m_f22; // Correction factor for activation, mode 2
675 
676  // Precomputed constants for efficiency
677  amrex::Real m_cons1, m_cons2, m_cons3, m_cons4, m_cons5;
678  amrex::Real m_cons6, m_cons7, m_cons8, m_cons9, m_cons10;
679  amrex::Real m_cons11, m_cons12, m_cons13, m_cons14, m_cons15;
680  amrex::Real m_cons16, m_cons17, m_cons18, m_cons19, m_cons20;
681  amrex::Real m_cons21, m_cons22, m_cons23, m_cons24, m_cons25;
682  amrex::Real m_cons26, m_cons27, m_cons28, m_cons29; [[maybe_unused]] amrex::Real m_cons30;
683  amrex::Real m_cons31, m_cons32, m_cons34, m_cons35; [[maybe_unused]] amrex::Real m_cons33;
684  amrex::Real m_cons36, m_cons37, m_cons38, m_cons39, m_cons40;
685  amrex::Real m_cons41;
686  // Set microphysics control parameters
687  m_inum = 1; // Use constant droplet number concentration
688  m_ndcnst = 250.0; // Droplet number concentration (cm^-3)
689  // Mathematical constants
690  m_pi = 3.1415926535897932384626434;
691 
692  m_R = 287.0; // Gas constant for dry air (J/kg/K)
693  m_Rd = 287.0; // Gas constant for dry air (J/kg/K)
694  m_Rv = 461.6; // Gas constant for water vapor (J/kg/K)
695  m_cp = 7.0*287.0/2.0; // Specific heat at constant pressure (J/kg/K)
696  m_g = 9.81; // Gravitational acceleration (m/s^2)
697  m_ep_2 = m_Rd / m_Rv; // Molecular weight ratio (Rd/Rv)
698 
699  // Reference density
700  m_rhosu = 85000.0/(287.15*273.15); // Standard air density at 850 mb (kg/m^3)
701 
702  // Densities for different hydrometeor species
703  m_rhow = 997.0; // Density of liquid water (kg/m^3)
704  m_rhoi = 500.0; // Bulk density of cloud ice (kg/m^3)
705  m_rhosn = 100.0; // Bulk density of snow (kg/m^3)
706 
707  // Set density for graupel or hail based on configuration
708  if (m_ihail == 0) {
709  m_rhog = 400.0; // Bulk density of graupel (kg/m^3)
710  } else {
711  m_rhog = 900.0; // Bulk density of hail (kg/m^3)
712  }
713 
714  // Fall speed parameters (V=AD^B) for different hydrometeors
715  // Cloud ice
716  m_ai = 700.0;
717  m_bi = 1.0;
718 
719  // Cloud droplets
720  m_ac = 3.0E7;
721  m_bc = 2.0;
722 
723  // Snow
724  m_as = 11.72;
725  m_bs = 0.41;
726 
727  // Rain
728  m_ar = 841.99667;
729  m_br = 0.8;
730 
731  // Graupel/hail (dependent on configuration)
732  if (m_ihail == 0) {
733  // Graupel parameters
734  m_ag = 19.3;
735  m_bg = 0.37;
736  } else {
737  // Hail parameters (Matsun and Huggins 1980)
738  m_ag = 114.5;
739  m_bg = 0.5;
740  }
741 
742  // Microphysical parameters
743  m_aimm = 0.66; // Parameter in Bigg immersion freezing
744  m_bimm = 100.0; // Parameter in Bigg immersion freezing
745  m_ecr = 1.0; // Collection efficiency between rain and snow/graupel
746  m_dcs = 125.0E-6; // Threshold size for cloud ice autoconversion (m)
747  m_mi0 = 4.0/3.0*m_pi*m_rhoi*std::pow(10.0E-6, 3); // Initial mass of nucleated ice crystal (kg)
748  m_mg0 = 1.6E-10; // Mass of embryo graupel (kg)
749 
750  // Ventilation parameters
751  m_f1s = 0.86; // Ventilation parameter for snow
752  m_f2s = 0.28; // Ventilation parameter for snow
753  m_f1r = 0.78; // Ventilation parameter for rain
754  m_f2r = 0.308; // Ventilation parameter for rain
755 
756  // Smallest allowed hydrometeor mixing ratio
757  m_qsmall = 1.0E-14;
758 
759  // Collection efficiencies
760  m_eii = 0.1; // Ice-ice collision efficiency
761  m_eci = 0.7; // Ice-droplet collision efficiency
762 
763  // Specific heat of liquid water (J/kg/K)
764  m_cpw = 4187.0;
765 
766  // Size distribution parameters
767  m_ci = m_rhoi * m_pi / 6.0;
768  m_di = 3.0;
769  m_cs = m_rhosn * m_pi / 6.0;
770  m_ds = 3.0;
771  m_cg = m_rhog * m_pi / 6.0;
772  m_dg = 3.0;
773 
774  // Radius of contact nuclei (m)
775  m_rin = 0.1E-6;
776 
777  // Mass of splintered ice particle (kg)
778  m_mmult = 4.0/3.0*m_pi*m_rhoi*std::pow(5.0E-6, 3);
779 
780  // Set lambda limits for size distributions
781  // Maximum and minimum values for lambda parameter in size distributions
782  m_lammaxi = 1.0/1.0E-6;
783  m_lammini = 1.0/(2.0*m_dcs + 100.0E-6);
784  m_lammaxr = 1.0/20.0E-6;
785  m_lamminr = 1.0/2800.0E-6;
786  m_lammaxs = 1.0/10.0E-6;
787  m_lammins = 1.0/2000.0E-6;
788  m_lammaxg = 1.0/20.0E-6;
789  m_lamming = 1.0/2000.0E-6;
790 
791  // Set CCN parameters for different environments
792  if (m_iact == 1) {
793  // Maritime CCN spectrum parameters (modified from Rasmussen et al. 2002)
794  // NCCN = C*S^K, where S is supersaturation in %
795  m_k1 = 0.4; // Exponent in CCN activation formula
796  m_c1 = 120.0; // Coefficient in CCN activation formula (cm^-3)
797  }
798 
799  // Initialize aerosol activation parameters for lognormal distribution
800  if (m_iact == 2) {
801  // Parameters for ammonium sulfate
802  m_mw = 0.018; // Molecular weight of water (kg/mol)
803  m_osm = 1.0; // Osmotic coefficient
804  m_vi = 3.0; // Number of ions dissociated in solution
805  m_epsm = 0.7; // Aerosol soluble fraction
806  m_rhoa = 1777.0; // Aerosol bulk density (kg/m^3)
807  m_map = 0.132; // Molecular weight of aerosol (kg/mol)
808  m_ma = 0.0284; // Molecular weight of air (kg/mol)
809  m_rr = 8.3145; // Universal gas constant (J/mol/K)
810  m_bact = m_vi * m_osm * m_epsm * m_mw * m_rhoa / (m_map * m_rhow);
811  // m_a_w = 2.0 * m_mw * 0.0761 / (m_rhow * m_r_v * 293.15); // "A" parameter
812 
813  // Aerosol size distribution parameters for MPACE (Morrison et al. 2007, JGR)
814  // Mode 1
815  m_rm1 = 0.052E-6; // Geometric mean radius, mode 1 (m)
816  m_sig1 = 2.04; // Standard deviation of aerosol size distribution, mode 1
817  m_nanew1 = 72.2E6; // Total aerosol concentration, mode 1 (m^-3)
818  m_f11 = 0.5 * std::exp(2.5 * std::pow(std::log(m_sig1), 2));
819  m_f21 = 1.0 + 0.25 * std::log(m_sig1);
820 
821  // Mode 2
822  m_rm2 = 1.3E-6; // Geometric mean radius, mode 2 (m)
823  m_sig2 = 2.5; // Standard deviation of aerosol size distribution, mode 2
824  m_nanew2 = 1.8E6; // Total aerosol concentration, mode 2 (m^-3)
825  m_f12 = 0.5 * std::exp(2.5 * std::pow(std::log(m_sig2), 2));
826  m_f22 = 1.0 + 0.25 * std::log(m_sig2);
827  }
828 
829  // Precompute constants for efficiency
830  m_cons1 = gamma_function(1.0 + m_ds) * m_cs;
831  m_cons2 = gamma_function(1.0 + m_dg) * m_cg;
832  m_cons3 = gamma_function(4.0 + m_bs) / 6.0;
833  m_cons4 = gamma_function(4.0 + m_br) / 6.0;
834  m_cons5 = gamma_function(1.0 + m_bs);
835  m_cons6 = gamma_function(1.0 + m_br);
836  m_cons7 = gamma_function(4.0 + m_bg) / 6.0;
837  m_cons8 = gamma_function(1.0 + m_bg);
838  m_cons9 = gamma_function(5.0/2.0 + m_br/2.0);
839  m_cons10 = gamma_function(5.0/2.0 + m_bs/2.0);
840  m_cons11 = gamma_function(5.0/2.0 + m_bg/2.0);
841  m_cons12 = gamma_function(1.0 + m_di) * m_ci;
842  m_cons13 = gamma_function(m_bs + 3.0) * m_pi / 4.0 * m_eci;
843  m_cons14 = gamma_function(m_bg + 3.0) * m_pi / 4.0 * m_eci;
844  m_cons15 = -1108.0 * m_eii * std::pow(m_pi, (1.0-m_bs)/3.0) *
845  std::pow(m_rhosn, (-2.0-m_bs)/3.0) / (4.0*720.0);
846  m_cons16 = gamma_function(m_bi + 3.0) * m_pi / 4.0 * m_eci;
847  m_cons17 = 4.0 * 2.0 * 3.0 * m_rhosu * m_pi * m_eci * m_eci *
848  gamma_function(2.0*m_bs + 2.0) / (8.0*(m_rhog-m_rhosn));
849  m_cons18 = m_rhosn * m_rhosn;
850  m_cons19 = m_rhow * m_rhow;
851  m_cons20 = 20.0 * m_pi * m_pi * m_rhow * m_bimm;
852  m_cons21 = 4.0 / (m_dcs * m_rhoi);
853  m_cons22 = m_pi * m_rhoi * std::pow(m_dcs, 3) / 6.0;
854  m_cons23 = m_pi / 4.0 * m_eii * gamma_function(m_bs + 3.0);
855  m_cons24 = m_pi / 4.0 * m_ecr * gamma_function(m_br + 3.0);
856  m_cons25 = m_pi * m_pi / 24.0 * m_rhow * m_ecr * gamma_function(m_br + 6.0);
857  m_cons26 = m_pi / 6.0 * m_rhow;
858  m_cons27 = gamma_function(1.0 + m_bi);
859  m_cons28 = gamma_function(4.0 + m_bi) / 6.0;
860  m_cons29 = 4.0/3.0 * m_pi * m_rhow * std::pow(25.0E-6, 3);
861  m_cons30 = 4.0/3.0 * m_pi * m_rhow;
862  m_cons31 = m_pi * m_pi * m_ecr * m_rhosn;
863  m_cons32 = m_pi / 2.0 * m_ecr;
864  m_cons33 = m_pi * m_pi * m_ecr * m_rhog;
865  m_cons34 = 5.0/2.0 + m_br/2.0;
866  m_cons35 = 5.0/2.0 + m_bs/2.0;
867  m_cons36 = 5.0/2.0 + m_bg/2.0;
868  m_cons37 = 4.0 * m_pi * 1.38E-23 / (6.0 * m_pi * m_rin);
869  m_cons38 = m_pi * m_pi / 3.0 * m_rhow;
870  m_cons39 = m_pi * m_pi / 36.0 * m_rhow * m_bimm;
871  m_cons40 = m_pi / 6.0 * m_bimm;
872  m_cons41 = m_pi * m_pi * m_ecr * m_rhow;
873 
874  // Set CCN parameters for different environments
875  if (m_iact == 1) {
876  // Maritime CCN spectrum parameters (modified from Rasmussen et al. 2002)
877  // NCCN = C*S^K, where S is supersaturation in %
878  m_k1 = 0.4; // Exponent in CCN activation formula
879  m_c1 = 120.0; // Coefficient in CCN activation formula (cm^-3)
880  }
881 
882  // Initialize aerosol activation parameters for IACT=2
883  if (m_iact == 2) {
884  // Parameters for ammonium sulfate
885  m_mw = 0.018; // Molecular weight of water (kg/mol)
886  m_osm = 1.0; // Osmotic coefficient
887  m_vi = 3.0; // Number of ions dissociated in solution
888  m_epsm = 0.7; // Aerosol soluble fraction
889  m_rhoa = 1777.0; // Aerosol bulk density (kg/m^3)
890  m_map = 0.132; // Molecular weight of aerosol (kg/mol)
891  m_ma = 0.0284; // Molecular weight of air (kg/mol)
892  m_rr = 8.3145; // Universal gas constant (J/mol/K)
893  m_bact = m_vi * m_osm * m_epsm * m_mw * m_rhoa / (m_map * m_rhow);
894 
895  // Aerosol size distribution parameters for MPACE (Morrison et al. 2007, JGR)
896  // Mode 1
897  m_rm1 = 0.052E-6; // Geometric mean radius, mode 1 (m)
898  m_sig1 = 2.04; // Standard deviation of aerosol size distribution, mode 1
899  m_nanew1 = 72.2E6; // Total aerosol concentration, mode 1 (m^-3)
900  m_f11 = 0.5 * std::exp(2.5 * std::pow(std::log(m_sig1), 2));
901  m_f21 = 1.0 + 0.25 * std::log(m_sig1);
902 
903  // Mode 2
904  m_rm2 = 1.3E-6; // Geometric mean radius, mode 2 (m)
905  m_sig2 = 2.5; // Standard deviation of aerosol size distribution, mode 2
906  m_nanew2 = 1.8E6; // Total aerosol concentration, mode 2 (m^-3)
907  m_f12 = 0.5 * std::exp(2.5 * std::pow(std::log(m_sig2), 2));
908  m_f22 = 1.0 + 0.25 * std::log(m_sig2);
909  }
910  // Set microphysics control parameters
911  m_iact = 2; // Lognormal aerosol activation
912  m_inuc = 0; // Mid-latitude ice nucleation (Cooper)
913  if (sc.moisture_type == MoistureType::Morrison_NoIce) {
914  m_iliq = 1; // Include ice processes
915  m_igraup = 1; // Include graupel processes
916  } else {
917  m_iliq = 0; // Include ice processes
918  m_igraup = 0; // Include graupel processes
919  }
920  m_ihail = 0; // Use graupel (0) instead of hail (1)
921  m_isub = 0; // Sub-grid vertical velocity option
922  m_do_radar_ref = false; // Disable radar reflectivity by default
923  amrex::Box boxD(box); boxD.makeSlab(2,0);
924 
925  ParmParse pp("erf");
926 
927  bool run_morr_cpp = true;
928 
929  bool use_morr_cpp_answer = true;
930  pp.query("use_morr_cpp_answer", use_morr_cpp_answer);
931  Print() << "use_morr_cpp_answer" << use_morr_cpp_answer <<std::endl;
932 
933  bool run_morr_fort = !use_morr_cpp_answer;
934 
935  std::string filename = std::string("output_cpp") + std::to_string(use_morr_cpp_answer) + ".txt";
936 
937  if(run_morr_cpp) {
938 
939  // Create dummy arrays for tendencies
940  FArrayBox qc3dten_fab(grown_box, 1); // CLOUD WATER MIXING RATIO TENDENCY
941  FArrayBox qi3dten_fab(grown_box, 1); // CLOUD ICE MIXING RATIO TENDENCY
942  FArrayBox qni3dten_fab(grown_box, 1); // SNOW MIXING RATIO TENDENCY
943  FArrayBox qr3dten_fab(grown_box, 1); // RAIN MIXING RATIO TENDENCY
944  FArrayBox ni3dten_fab(grown_box, 1); // CLOUD ICE NUMBER CONCENTRATION
945  FArrayBox ns3dten_fab(grown_box, 1); // SNOW NUMBER CONCENTRATION
946  FArrayBox nr3dten_fab(grown_box, 1); // RAIN NUMBER CONCENTRATION
947 
948  // Get array references
949  auto const& qc3dten = qc3dten_fab.array(); // CLOUD WATER MIXING RATIO TENDENCY (KG/KG/S)
950  auto const& qi3dten = qi3dten_fab.array(); // CLOUD ICE MIXING RATIO TENDENCY (KG/KG/S)
951  auto const& qni3dten = qni3dten_fab.array(); // SNOW MIXING RATIO TENDENCY (KG/KG/S)
952  auto const& qr3dten = qr3dten_fab.array(); // RAIN MIXING RATIO TENDENCY (KG/KG/S)
953  auto const& ni3dten = ni3dten_fab.array(); // CLOUD ICE NUMBER CONCENTRATION (1/KG/S)
954  auto const& ns3dten = ns3dten_fab.array(); // SNOW NUMBER CONCENTRATION (1/KG/S)
955  auto const& nr3dten = nr3dten_fab.array(); // RAIN NUMBER CONCENTRATION (1/KG/S)
956 
957  // Initialize tendencies to zero (no precipitation implemented yet)
958  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
959  qc3dten(i,j,k) = 0.0;
960  qi3dten(i,j,k) = 0.0;
961  qni3dten(i,j,k) = 0.0;
962  qr3dten(i,j,k) = 0.0;
963  ni3dten(i,j,k) = 0.0;
964  ns3dten(i,j,k) = 0.0;
965  nr3dten(i,j,k) = 0.0;
966  });
967 
968  // Create arrays for mixing ratios and number concentrations
969  FArrayBox qc3d_fab(grown_box, 1); // CLOUD WATER MIXING RATIO
970  FArrayBox qi3d_fab(grown_box, 1); // CLOUD ICE MIXING RATIO
971  FArrayBox qni3d_fab(grown_box, 1); // SNOW MIXING RATIO
972  FArrayBox qr3d_fab(grown_box, 1); // RAIN MIXING RATIO
973  FArrayBox ni3d_fab(grown_box, 1); // CLOUD ICE NUMBER CONCENTRATION
974  FArrayBox ns3d_fab(grown_box, 1); // SNOW NUMBER CONCENTRATION
975  FArrayBox nr3d_fab(grown_box, 1); // RAIN NUMBER CONCENTRATION
976 
977  // Get array references
978  auto const& qc3d = qc3d_fab.array(); // CLOUD WATER MIXING RATIO (KG/KG)
979  auto const& qi3d = qi3d_fab.array(); // CLOUD ICE MIXING RATIO (KG/KG)
980  auto const& qni3d = qni3d_fab.array(); // SNOW MIXING RATIO (KG/KG)
981  auto const& qr3d = qr3d_fab.array(); // RAIN MIXING RATIO (KG/KG)
982  auto const& ni3d = ni3d_fab.array(); // CLOUD ICE NUMBER CONCENTRATION (1/KG)
983  auto const& ns3d = ns3d_fab.array(); // SNOW NUMBER CONCENTRATION (1/KG)
984  auto const& nr3d = nr3d_fab.array(); // RAIN NUMBER CONCENTRATION (1/KG)
985 
986  // Initialize mixing ratios and number concentrations to zero
987  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
988  qc3d(i,j,k) = 0.0;
989  qi3d(i,j,k) = 0.0;
990  qni3d(i,j,k) = 0.0;
991  qr3d(i,j,k) = 0.0;
992  ni3d(i,j,k) = 0.0;
993  ns3d(i,j,k) = 0.0;
994  nr3d(i,j,k) = 0.0;
995  });
996 
997  // Create arrays for temperature, vapor, and pressure variables
998  FArrayBox t3dten_fab(grown_box, 1); // TEMPERATURE TENDENCY
999  FArrayBox qv3dten_fab(grown_box, 1); // WATER VAPOR MIXING RATIO TENDENCY
1000  FArrayBox t3d_fab(grown_box, 1); // TEMPERATURE
1001  FArrayBox qv3d_fab(grown_box, 1); // WATER VAPOR MIXING RATIO
1002  FArrayBox pres_fab(grown_box, 1); // ATMOSPHERIC PRESSURE
1003  FArrayBox dzq_fab(grown_box, 1); // DIFFERENCE IN HEIGHT ACROSS LEVEL
1004  FArrayBox w3d_fab(grown_box, 1); // GRID-SCALE VERTICAL VELOCITY
1005 
1006  // WRF-chem variables
1007  FArrayBox nc3d_fab(grown_box, 1); // CLOUD DROPLET NUMBER CONCENTRATION
1008  FArrayBox nc3dten_fab(grown_box, 1); // CLOUD DROPLET NUMBER CONCENTRATION TENDENCY
1009 
1010  // Graupel variables
1011  FArrayBox qg3dten_fab(grown_box, 1); // GRAUPEL MIX RATIO TENDENCY
1012  FArrayBox ng3dten_fab(grown_box, 1); // GRAUPEL NUMB CONC TENDENCY
1013  FArrayBox qg3d_fab(grown_box, 1); // GRAUPEL MIX RATIO
1014  FArrayBox ng3d_fab(grown_box, 1); // GRAUPEL NUMBER CONC
1015 
1016  // Sedimentation tendencies
1017  FArrayBox qgsten_fab(grown_box, 1); // GRAUPEL SED TEND
1018  FArrayBox qrsten_fab(grown_box, 1); // RAIN SED TEND
1019  FArrayBox qisten_fab(grown_box, 1); // CLOUD ICE SED TEND
1020  FArrayBox qnisten_fab(grown_box, 1); // SNOW SED TEND
1021  FArrayBox qcsten_fab(grown_box, 1); // CLOUD WAT SED TEND
1022 
1023  // Cumulus tendencies
1024  FArrayBox qrcu1d_fab(grown_box, 1); // RAIN FROM CUMULUS PARAMETERIZATION
1025  FArrayBox qscu1d_fab(grown_box, 1); // SNOW FROM CUMULUS PARAMETERIZATION
1026  FArrayBox qicu1d_fab(grown_box, 1); // ICE FROM CUMULUS PARAMETERIZATION
1027 
1028  // Get array references
1029  auto const& t3dten = t3dten_fab.array(); // TEMPERATURE TENDENCY (K/S)
1030  auto const& qv3dten = qv3dten_fab.array(); // WATER VAPOR MIXING RATIO TENDENCY (KG/KG/S)
1031  auto const& t3d = t3d_fab.array(); // TEMPERATURE (K)
1032  auto const& qv3d = qv3d_fab.array(); // WATER VAPOR MIXING RATIO (KG/KG)
1033  auto const& pres = pres_fab.array(); // ATMOSPHERIC PRESSURE (PA)
1034  auto const& dzq = dzq_fab.array(); // DIFFERENCE IN HEIGHT ACROSS LEVEL (m)
1035  auto const& w3d = w3d_fab.array(); // GRID-SCALE VERTICAL VELOCITY (M/S)
1036 
1037  // WRF-chem variables
1038  auto const& nc3d = nc3d_fab.array(); // CLOUD DROPLET NUMBER CONCENTRATION
1039  auto const& nc3dten = nc3dten_fab.array(); // CLOUD DROPLET NUMBER CONCENTRATION TENDENCY
1040 
1041  // Graupel variables
1042  auto const& qg3dten = qg3dten_fab.array(); // GRAUPEL MIX RATIO TENDENCY (KG/KG/S)
1043  auto const& ng3dten = ng3dten_fab.array(); // GRAUPEL NUMB CONC TENDENCY (1/KG/S)
1044  auto const& qg3d = qg3d_fab.array(); // GRAUPEL MIX RATIO (KG/KG)
1045  auto const& ng3d = ng3d_fab.array(); // GRAUPEL NUMBER CONC (1/KG)
1046 
1047  // Sedimentation tendencies
1048  auto const& qgsten = qgsten_fab.array(); // GRAUPEL SED TEND (KG/KG/S)
1049  auto const& qrsten = qrsten_fab.array(); // RAIN SED TEND (KG/KG/S)
1050  auto const& qisten = qisten_fab.array(); // CLOUD ICE SED TEND (KG/KG/S)
1051  auto const& qnisten = qnisten_fab.array(); // SNOW SED TEND (KG/KG/S)
1052  auto const& qcsten = qcsten_fab.array(); // CLOUD WAT SED TEND (KG/KG/S)
1053 
1054  // Cumulus tendencies
1055  auto const& qrcu1d = qrcu1d_fab.array(); // RAIN FROM CUMULUS PARAMETERIZATION
1056  auto const& qscu1d = qscu1d_fab.array(); // SNOW FROM CUMULUS PARAMETERIZATION
1057  auto const& qicu1d = qicu1d_fab.array(); // ICE FROM CUMULUS PARAMETERIZATION
1058 
1059  // Initialize tendency arrays to zero
1060  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1061  t3dten(i,j,k) = 0.0;
1062  qv3dten(i,j,k) = 0.0;
1063  nc3dten(i,j,k) = 0.0;
1064  qg3dten(i,j,k) = 0.0;
1065  ng3dten(i,j,k) = 0.0;
1066  qgsten(i,j,k) = 0.0;
1067  qrsten(i,j,k) = 0.0;
1068  qisten(i,j,k) = 0.0;
1069  qnisten(i,j,k) = 0.0;
1070  qcsten(i,j,k) = 0.0;
1071  });
1072 
1073  // Create arrays for precipitation rates
1074  FArrayBox precrt_fab(grown_box, 1); // TOTAL PRECIP PER TIME STEP
1075  FArrayBox snowrt_fab(grown_box, 1); // SNOW PER TIME STEP
1076  FArrayBox snowprt_fab(grown_box, 1); // TOTAL CLOUD ICE PLUS SNOW PER TIME STEP
1077  FArrayBox grplprt_fab(grown_box, 1); // TOTAL GRAUPEL PER TIME STEP
1078 
1079  // Create arrays for effective radii
1080  FArrayBox effc_fab(grown_box, 1); // DROPLET EFFECTIVE RADIUS
1081  FArrayBox effi_fab(grown_box, 1); // CLOUD ICE EFFECTIVE RADIUS
1082  FArrayBox effs_fab(grown_box, 1); // SNOW EFFECTIVE RADIUS
1083  FArrayBox effr_fab(grown_box, 1); // RAIN EFFECTIVE RADIUS
1084  FArrayBox effg_fab(grown_box, 1); // GRAUPEL EFFECTIVE RADIUS
1085 
1086  // Get array references for precipitation rates
1087  auto const& precrt = precrt_fab.array(); // TOTAL PRECIP PER TIME STEP (mm)
1088  auto const& snowrt = snowrt_fab.array(); // SNOW PER TIME STEP (mm)
1089  auto const& snowprt = snowprt_fab.array(); // TOTAL CLOUD ICE PLUS SNOW PER TIME STEP (mm)
1090  auto const& grplprt = grplprt_fab.array(); // TOTAL GRAUPEL PER TIME STEP (mm)
1091 
1092  // Get array references for effective radii
1093  auto const& effc = effc_fab.array(); // DROPLET EFFECTIVE RADIUS (MICRON)
1094  auto const& effi = effi_fab.array(); // CLOUD ICE EFFECTIVE RADIUS (MICRON)
1095  auto const& effs = effs_fab.array(); // SNOW EFFECTIVE RADIUS (MICRON)
1096  auto const& effr = effr_fab.array(); // RAIN EFFECTIVE RADIUS (MICRON)
1097  auto const& effg = effg_fab.array(); // GRAUPEL EFFECTIVE RADIUS (MICRON)
1098 
1099  // Initialize these arrays to zero (they will be computed later)
1100  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1101  precrt(i,j,k) = 0.0;
1102  snowrt(i,j,k) = 0.0;
1103  snowprt(i,j,k) = 0.0;
1104  grplprt(i,j,k) = 0.0;
1105  effc(i,j,k) = 0.0;
1106  effi(i,j,k) = 0.0;
1107  effs(i,j,k) = 0.0;
1108  effr(i,j,k) = 0.0;
1109  effg(i,j,k) = 0.0;
1110  });
1111 
1112  // Create FArrayBoxes for scalar variables
1113  FArrayBox rho_fab(grown_box, 1);
1114  FArrayBox mu_fab(grown_box, 1);
1115  FArrayBox ain_fab(grown_box, 1);
1116  FArrayBox arn_fab(grown_box, 1);
1117  FArrayBox asn_fab(grown_box, 1);
1118  FArrayBox acn_fab(grown_box, 1);
1119  FArrayBox agn_fab(grown_box, 1);
1120 
1121  // Get Array4 views
1122  auto const& rho = rho_fab.array();
1123  auto const& mu = mu_fab.array();
1124  auto const& ain = ain_fab.array();
1125  auto const& arn = arn_fab.array();
1126  auto const& asn = asn_fab.array();
1127  auto const& acn = acn_fab.array();
1128  auto const& agn = agn_fab.array();
1129 
1130  // Initialize all values to zero
1131  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1132  rho(i,j,k) = 0.0;
1133  mu(i,j,k) = 0.0;
1134  ain(i,j,k) = 0.0;
1135  arn(i,j,k) = 0.0;
1136  asn(i,j,k) = 0.0;
1137  acn(i,j,k) = 0.0;
1138  agn(i,j,k) = 0.0;
1139  });
1140 
1141  // Create FArrayBoxes for fall speed working variables
1142  FArrayBox dumi_fab(grown_box, 1);
1143  FArrayBox dumr_fab(grown_box, 1);
1144  FArrayBox dumfni_fab(grown_box, 1);
1145  FArrayBox dumg_fab(grown_box, 1);
1146  FArrayBox dumfng_fab(grown_box, 1);
1147  FArrayBox uni_fab(grown_box, 1);
1148  FArrayBox umi_fab(grown_box, 1);
1149  FArrayBox umr_fab(grown_box, 1);
1150  FArrayBox fr_fab(grown_box, 1);
1151  FArrayBox fi_fab(grown_box, 1);
1152  FArrayBox fni_fab(grown_box, 1);
1153  FArrayBox fg_fab(grown_box, 1);
1154  FArrayBox fng_fab(grown_box, 1);
1155  FArrayBox rgvm_fab(grown_box, 1);
1156  FArrayBox faloutr_fab(grown_box, 1);
1157  FArrayBox falouti_fab(grown_box, 1);
1158  FArrayBox faloutni_fab(grown_box, 1);
1159  FArrayBox faltndr_fab(grown_box, 1);
1160  FArrayBox faltndi_fab(grown_box, 1);
1161  FArrayBox faltndni_fab(grown_box, 1);
1162  FArrayBox dumqs_fab(grown_box, 1);
1163  FArrayBox dumfns_fab(grown_box, 1);
1164  FArrayBox ums_fab(grown_box, 1);
1165  FArrayBox uns_fab(grown_box, 1);
1166  FArrayBox fs_fab(grown_box, 1);
1167  FArrayBox fns_fab(grown_box, 1);
1168  FArrayBox falouts_fab(grown_box, 1);
1169  FArrayBox faloutns_fab(grown_box, 1);
1170  FArrayBox faloutg_fab(grown_box, 1);
1171  FArrayBox faloutng_fab(grown_box, 1);
1172  FArrayBox faltnds_fab(grown_box, 1);
1173  FArrayBox faltndns_fab(grown_box, 1);
1174  FArrayBox unr_fab(grown_box, 1);
1175  FArrayBox faltndg_fab(grown_box, 1);
1176  FArrayBox faltndng_fab(grown_box, 1);
1177  FArrayBox dumc_fab(grown_box, 1);
1178  FArrayBox dumfnc_fab(grown_box, 1);
1179  FArrayBox unc_fab(grown_box, 1);
1180  FArrayBox umc_fab(grown_box, 1);
1181  FArrayBox ung_fab(grown_box, 1);
1182  FArrayBox umg_fab(grown_box, 1);
1183  FArrayBox fc_fab(grown_box, 1);
1184  FArrayBox faloutc_fab(grown_box, 1);
1185  FArrayBox faloutnc_fab(grown_box, 1);
1186  FArrayBox faltndc_fab(grown_box, 1);
1187  FArrayBox faltndnc_fab(grown_box, 1);
1188  FArrayBox fnc_fab(grown_box, 1);
1189  FArrayBox dumfnr_fab(grown_box, 1);
1190  FArrayBox faloutnr_fab(grown_box, 1);
1191  FArrayBox faltndnr_fab(grown_box, 1);
1192  FArrayBox fnr_fab(grown_box, 1);
1193  FArrayBox dlams_fab(grown_box, 1);
1194  FArrayBox dlamr_fab(grown_box, 1);
1195  FArrayBox dlami_fab(grown_box, 1);
1196  FArrayBox dlamc_fab(grown_box, 1);
1197  FArrayBox dlamg_fab(grown_box, 1);
1198 
1199  // Create Array4 references
1200  auto const& dumi = dumi_fab.array();
1201  auto const& dumr = dumr_fab.array();
1202  auto const& dumfni = dumfni_fab.array();
1203  auto const& dumg = dumg_fab.array();
1204  auto const& dumfng = dumfng_fab.array();
1205  auto const& uni = uni_fab.array();
1206  auto const& umi = umi_fab.array();
1207  auto const& umr = umr_fab.array();
1208  auto const& fr = fr_fab.array();
1209  auto const& fi = fi_fab.array();
1210  auto const& fni = fni_fab.array();
1211  auto const& fg = fg_fab.array();
1212  auto const& fng = fng_fab.array();
1213  auto const& rgvm = rgvm_fab.array();
1214  auto const& faloutr = faloutr_fab.array();
1215  auto const& falouti = falouti_fab.array();
1216  auto const& faloutni = faloutni_fab.array();
1217  auto const& faltndr = faltndr_fab.array();
1218  auto const& faltndi = faltndi_fab.array();
1219  auto const& faltndni = faltndni_fab.array();
1220  auto const& dumqs = dumqs_fab.array();
1221  auto const& dumfns = dumfns_fab.array();
1222  auto const& ums = ums_fab.array();
1223  auto const& uns = uns_fab.array();
1224  auto const& fs = fs_fab.array();
1225  auto const& fns = fns_fab.array();
1226  auto const& falouts = falouts_fab.array();
1227  auto const& faloutns = faloutns_fab.array();
1228  auto const& faloutg = faloutg_fab.array();
1229  auto const& faloutng = faloutng_fab.array();
1230  auto const& faltnds = faltnds_fab.array();
1231  auto const& faltndns = faltndns_fab.array();
1232  auto const& unr = unr_fab.array();
1233  auto const& faltndg = faltndg_fab.array();
1234  auto const& faltndng = faltndng_fab.array();
1235  auto const& dumc = dumc_fab.array();
1236  auto const& dumfnc = dumfnc_fab.array();
1237  auto const& unc = unc_fab.array();
1238  auto const& umc = umc_fab.array();
1239  auto const& ung = ung_fab.array();
1240  auto const& umg = umg_fab.array();
1241  auto const& fc = fc_fab.array();
1242  auto const& faloutc = faloutc_fab.array();
1243  auto const& faloutnc = faloutnc_fab.array();
1244  auto const& faltndc = faltndc_fab.array();
1245  auto const& faltndnc = faltndnc_fab.array();
1246  auto const& fnc = fnc_fab.array();
1247  auto const& dumfnr = dumfnr_fab.array();
1248  auto const& faloutnr = faloutnr_fab.array();
1249  auto const& faltndnr = faltndnr_fab.array();
1250  auto const& fnr = fnr_fab.array();
1251  auto const& dlams = dlams_fab.array();
1252  auto const& dlamr = dlamr_fab.array();
1253  auto const& dlami = dlami_fab.array();
1254  auto const& dlamc = dlamc_fab.array();
1255  auto const& dlamg = dlamg_fab.array();
1256 
1257  // Initialize arrays to 0
1258  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1259  dlams(i,j,k) = 0.0;
1260  dlamr(i,j,k) = 0.0;
1261  dlami(i,j,k) = 0.0;
1262  dlamc(i,j,k) = 0.0;
1263  dlamg(i,j,k) = 0.0;
1264  dumi(i,j,k) = 0.0;
1265  dumr(i,j,k) = 0.0;
1266  dumfni(i,j,k) = 0.0;
1267  dumg(i,j,k) = 0.0;
1268  dumfng(i,j,k) = 0.0;
1269  uni(i,j,k) = 0.0;
1270  umi(i,j,k) = 0.0;
1271  umr(i,j,k) = 0.0;
1272  fr(i,j,k) = 0.0;
1273  fi(i,j,k) = 0.0;
1274  fni(i,j,k) = 0.0;
1275  fg(i,j,k) = 0.0;
1276  fng(i,j,k) = 0.0;
1277  rgvm(i,j,k) = 0.0;
1278  faloutr(i,j,k) = 0.0;
1279  falouti(i,j,k) = 0.0;
1280  faloutni(i,j,k) = 0.0;
1281  faltndr(i,j,k) = 0.0;
1282  faltndi(i,j,k) = 0.0;
1283  faltndni(i,j,k) = 0.0;
1284  dumqs(i,j,k) = 0.0;
1285  dumfns(i,j,k) = 0.0;
1286  ums(i,j,k) = 0.0;
1287  uns(i,j,k) = 0.0;
1288  fs(i,j,k) = 0.0;
1289  fns(i,j,k) = 0.0;
1290  falouts(i,j,k) = 0.0;
1291  faloutns(i,j,k) = 0.0;
1292  faloutg(i,j,k) = 0.0;
1293  faloutng(i,j,k) = 0.0;
1294  faltnds(i,j,k) = 0.0;
1295  faltndns(i,j,k) = 0.0;
1296  unr(i,j,k) = 0.0;
1297  faltndg(i,j,k) = 0.0;
1298  faltndng(i,j,k) = 0.0;
1299  dumc(i,j,k) = 0.0;
1300  dumfnc(i,j,k) = 0.0;
1301  unc(i,j,k) = 0.0;
1302  umc(i,j,k) = 0.0;
1303  ung(i,j,k) = 0.0;
1304  umg(i,j,k) = 0.0;
1305  fc(i,j,k) = 0.0;
1306  faloutc(i,j,k) = 0.0;
1307  faloutnc(i,j,k) = 0.0;
1308  faltndc(i,j,k) = 0.0;
1309  faltndnc(i,j,k) = 0.0;
1310  fnc(i,j,k) = 0.0;
1311  dumfnr(i,j,k) = 0.0;
1312  faloutnr(i,j,k) = 0.0;
1313  faltndnr(i,j,k) = 0.0;
1314  fnr(i,j,k) = 0.0;
1315  });
1316 
1317  // Create FArrayBoxes for thermodynamic variables
1318  FArrayBox xxls_fab(grown_box, 1); // Latent heat of sublimation
1319  FArrayBox xxlv_fab(grown_box, 1); // Latent heat of vaporization
1320  FArrayBox cpm_fab(grown_box, 1); // Specific heat at constant pressure for moist air
1321  FArrayBox xlf_fab(grown_box, 1); // Latent heat of freezing
1322 
1323  // Get Array4 references
1324  auto const& xxls = xxls_fab.array(); // XXLS: Latent heat of sublimation
1325  auto const& xxlv = xxlv_fab.array(); // XXLV: Latent heat of vaporization
1326  auto const& cpm = cpm_fab.array(); // CPM: Specific heat at constant pressure for moist air
1327  auto const& xlf = xlf_fab.array(); // XLF: Latent heat of freezing
1328 
1329  // Initialize values to zero
1330  ParallelFor(grown_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1331  xxls(i,j,k) = 0.0;
1332  xxlv(i,j,k) = 0.0;
1333  cpm(i,j,k) = 0.0;
1334  xlf(i,j,k) = 0.0;
1335  });
1336 
1337  ////////////////////////////////////////////////////////////
1338  // ParallelFor for testing partial C++ implementation
1339  // NOTE: Currently all Array4 values are copied to locals
1340  // This means we're not updating or outputting anything
1341  ////////////////////////////////////////////////////////////
1342  ParallelFor( box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
1343  {
1344  // Tendencies and mixing ratios
1345  qc3d(i,j,k) = qcl_arr(i,j,k); // CLOUD WATER MIXING RATIO
1346  qi3d(i,j,k) = qci_arr(i,j,k); // CLOUD ICE MIXING RATIO
1347  qni3d(i,j,k) = qps_arr(i,j,k); // SNOW MIXING RATIO
1348  qr3d(i,j,k) = qpr_arr(i,j,k); // RAIN MIXING RATIO
1349  ni3d(i,j,k) = ni_arr(i,j,k); // CLOUD ICE NUMBER CONCENTRATION
1350  ns3d(i,j,k) = ns_arr(i,j,k); // SNOW NUMBER CONCENTRATION
1351  nr3d(i,j,k) = nr_arr(i,j,k); // RAIN NUMBER CONCENTRATION
1352  nc3d(i,j,k) = nc_arr(i,j,k); // RAIN NUMBER CONCENTRATION
1353 
1354  t3d(i,j,k) = theta_arr(i,j,k) * pii_arr(i,j,k); // TEMPERATURE
1355  qv3d(i,j,k) = qv_arr(i,j,k); // WATER VAPOR MIXING RATIO
1356  pres(i,j,k) = pres_arr(i,j,k); // ATMOSPHERIC PRESSURE
1357  dzq(i,j,k) = dz_arr(i,j,k); // DIFFERENCE IN HEIGHT ACROSS LEVEL
1358  w3d(i,j,k) = w_arr(i,j,k); // GRID-SCALE VERTICAL VELOCITY
1359  qg3d(i,j,k) = qpg_arr(i,j,k); // GRAUPEL MIX RATIO
1360  ng3d(i,j,k) = ng_arr(i,j,k); // GRAUPEL NUMBER CONC
1361  qrcu1d(i,j,k) = qrcuten_arr(i,j,k); // RAIN FROM CUMULUS PARAMETERIZATION
1362  qscu1d(i,j,k) = qscuten_arr(i,j,k); // SNOW FROM CUMULUS PARAMETERIZATION
1363  qicu1d(i,j,k) = qicuten_arr(i,j,k); // ICE FROM CUMULUS PARAMETERIZATION
1364  });
1365  ParallelFor( boxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
1366  {
1367  int ltrue=0; // LTRUE: SWITCH = 0: NO HYDROMETEORS IN COLUMN, = 1: HYDROMETEORS IN COLUMN
1368  int nstep; // NSTEP: Timestep counter
1369  int iinum=m_inum; // iinum: Integer control variable
1370 
1371  for(int k=klo; k<=khi; k++) {
1372  // Model input parameters
1373  //amrex::Real dt; // DT: MODEL TIME STEP (SEC)
1374  //amrex::Real lami(i,j,k); // LAMI: Slope parameter for cloud ice (m^-1)
1375 
1376  // Microphysical processes
1377  [[maybe_unused]] amrex::Real nsubc; // NSUBC: Loss of NC during evaporation
1378  amrex::Real nsubi; // NSUBI: Loss of NI during sublimation
1379  amrex::Real nsubs; // NSUBS: Loss of NS during sublimation
1380  amrex::Real nsubr; // NSUBR: Loss of NR during evaporation
1381  amrex::Real prd; // PRD: Deposition cloud ice
1382  amrex::Real pre; // PRE: Evaporation of rain
1383  amrex::Real prds; // PRDS: Deposition snow
1384  amrex::Real nnuccc; // NNUCCC: Change N due to contact freezing droplets
1385  amrex::Real mnuccc; // MNUCCC: Change Q due to contact freezing droplets
1386  amrex::Real pra; // PRA: Accretion droplets by rain
1387  amrex::Real prc; // PRC: Autoconversion droplets
1388  amrex::Real pcc; // PCC: Condensation/evaporation droplets
1389  amrex::Real nnuccd; // NNUCCD: Change N freezing aerosol (primary ice nucleation)
1390  amrex::Real mnuccd; // MNUCCD: Change Q freezing aerosol (primary ice nucleation)
1391  amrex::Real mnuccr; // MNUCCR: Change Q due to contact freezing rain
1392  amrex::Real nnuccr; // NNUCCR: Change N due to contact freezing rain
1393  amrex::Real npra; // NPRA: Change N due to droplet accretion by rain
1394  amrex::Real nragg; // NRAGG: Self-collection/breakup of rain
1395  amrex::Real nsagg; // NSAGG: Self-collection of snow
1396  amrex::Real nprc; // NPRC: Change NC autoconversion droplets
1397  amrex::Real nprc1; // NPRC1: Change NR autoconversion droplets
1398  amrex::Real prai; // PRAI: Change Q accretion cloud ice by snow
1399  amrex::Real prci; // PRCI: Change Q autoconversion cloud ice to snow
1400  amrex::Real psacws; // PSACWS: Change Q droplet accretion by snow
1401  amrex::Real npsacws; // NPSACWS: Change N droplet accretion by snow
1402  amrex::Real psacwi; // PSACWI: Change Q droplet accretion by cloud ice
1403  amrex::Real npsacwi; // NPSACWI: Change N droplet accretion by cloud ice
1404  amrex::Real nprci; // NPRCI: Change N autoconversion cloud ice by snow
1405  amrex::Real nprai; // NPRAI: Change N accretion cloud ice
1406  amrex::Real nmults; // NMULTS: Ice multiplication due to riming droplets by snow
1407  amrex::Real nmultr; // NMULTR: Ice multiplication due to riming rain by snow
1408  amrex::Real qmults; // QMULTS: Change Q due to ice multiplication droplets/snow
1409  amrex::Real qmultr; // QMULTR: Change Q due to ice multiplication rain/snow
1410  amrex::Real pracs; // PRACS: Change Q rain-snow collection
1411  amrex::Real npracs; // NPRACS: Change N rain-snow collection
1412  [[maybe_unused]] amrex::Real pccn; // PCCN: Change Q droplet activation
1413  amrex::Real psmlt; // PSMLT: Change Q melting snow to rain
1414  amrex::Real evpms; // EVPMS: Change Q melting snow evaporating
1415  amrex::Real nsmlts; // NSMLTS: Change N melting snow
1416  amrex::Real nsmltr; // NSMLTR: Change N melting snow to rain
1417  amrex::Real piacr; // PIACR: Change QR, ice-rain collection
1418  amrex::Real niacr; // NIACR: Change N, ice-rain collection
1419  amrex::Real praci; // PRACI: Change QI, ice-rain collection
1420  amrex::Real piacrs; // PIACRS: Change QR, ice rain collision, added to snow
1421  amrex::Real niacrs; // NIACRS: Change N, ice rain collision, added to snow
1422  amrex::Real pracis; // PRACIS: Change QI, ice rain collision, added to snow
1423  amrex::Real eprd; // EPRD: Sublimation cloud ice
1424  amrex::Real eprds; // EPRDS: Sublimation snow
1425 
1426  // Graupel processes
1427  amrex::Real pracg; // PRACG: Change in Q collection rain by graupel
1428  amrex::Real psacwg; // PSACWG: Change in Q collection droplets by graupel
1429  amrex::Real pgsacw; // PGSACW: Conversion Q to graupel due to collection droplets by snow
1430  amrex::Real pgracs; // PGRACS: Conversion Q to graupel due to collection rain by snow
1431  amrex::Real prdg; // PRDG: Deposition of graupel
1432  amrex::Real eprdg; // EPRDG: Sublimation of graupel
1433  amrex::Real evpmg; // EVPMG: Change Q melting of graupel and evaporation
1434  amrex::Real pgmlt; // PGMLT: Change Q melting of graupel
1435  amrex::Real npracg; // NPRACG: Change N collection rain by graupel
1436  amrex::Real npsacwg; // NPSACWG: Change N collection droplets by graupel
1437  amrex::Real nscng; // NSCNG: Change N conversion to graupel due to collection droplets by snow
1438  amrex::Real ngracs; // NGRACS: Change N conversion to graupel due to collection rain by snow
1439  amrex::Real ngmltg; // NGMLTG: Change N melting graupel
1440  amrex::Real ngmltr; // NGMLTR: Change N melting graupel to rain
1441  amrex::Real nsubg; // NSUBG: Change N sublimation/deposition of graupel
1442  amrex::Real psacr; // PSACR: Conversion due to collection of snow by rain
1443  amrex::Real nmultg; // NMULTG: Ice multiplication due to accretion droplets by graupel
1444  amrex::Real nmultrg; // NMULTRG: Ice multiplication due to accretion rain by graupel
1445  amrex::Real qmultg; // QMULTG: Change Q due to ice multiplication droplets/graupel
1446  amrex::Real qmultrg; // QMULTRG: Change Q due to ice multiplication rain/graupel
1447 
1448  // Time-varying atmospheric parameters
1449  amrex::Real kap; // KAP: Thermal conductivity of air
1450  amrex::Real evs; // EVS: Saturation vapor pressure
1451  amrex::Real eis; // EIS: Ice saturation vapor pressure
1452  amrex::Real qvs; // QVS: Saturation mixing ratio
1453  amrex::Real qvi; // QVI: Ice saturation mixing ratio
1454  amrex::Real qvqvs; // QVQVS: Saturation ratio
1455  amrex::Real qvqvsi; // QVQVSI: Ice saturation ratio
1456  amrex::Real dv; // DV: Diffusivity of water vapor in air
1457  amrex::Real sc_schmidt; // SC: Schmidt number
1458  amrex::Real ab; // AB: Correction to condensation rate due to latent heating
1459  amrex::Real abi; // ABI: Correction to deposition rate due to latent heating
1460 
1461  // Dummy variables
1462  amrex::Real dum; // DUM: General dummy variable
1463  amrex::Real dum1; // DUM1: General dummy variable
1464  [[maybe_unused]] amrex::Real dum2; // DUM2: General dummy variable
1465  amrex::Real dumt; // DUMT: Dummy variable for temperature
1466  amrex::Real dumqv; // DUMQV: Dummy variable for water vapor
1467  amrex::Real dumqss; // DUMQSS: Dummy saturation mixing ratio
1468  [[maybe_unused]] amrex::Real dumqsi; // DUMQSI: Dummy ice saturation mixing ratio
1469  amrex::Real dums; // DUMS: General dummy variable
1470 
1471  // Prognostic supersaturation
1472  amrex::Real dqsdt; // DQSDT: Change of saturation mixing ratio with temperature
1473  amrex::Real dqsidt; // DQSIDT: Change in ice saturation mixing ratio with temperature
1474 
1475  amrex::Real epsi; // EPSI: 1/phase relaxation time (see M2005), ice
1476  amrex::Real epss; // EPSS: 1/phase relaxation time (see M2005), snow
1477  amrex::Real epsr; // EPSR: 1/phase relaxation time (see M2005), rain
1478  amrex::Real epsg; // EPSG: 1/phase relaxation time (see M2005), graupel
1479  amrex::Real kc2; // KC2: Total ice nucleation rate
1480  amrex::Real di0; // DC0: Characteristic diameter for ice
1481  [[maybe_unused]] amrex::Real dc0; // DC0: Characteristic diameter for cloud droplets
1482  amrex::Real ds0; // DS0: Characteristic diameter for snow
1483  amrex::Real dg0; // DG0: Characteristic diameter for graupel
1484  amrex::Real dumqc; // DUMQC: Dummy variable for cloud water mixing ratio
1485  [[maybe_unused]] amrex::Real dumqr; // DUMQR: Dummy variable for rain mixing ratio
1486  amrex::Real ratio; // RATIO: General ratio variable
1487  amrex::Real sum_dep; // SUM_DEP: Sum of deposition/sublimation
1488  amrex::Real fudgef; // FUDGEF: Adjustment factor
1489  // For WRF-CHEM
1490  [[maybe_unused]] amrex::Real c2prec; // C2PREC: Cloud to precipitation conversion
1491  [[maybe_unused]] amrex::Real csed; // CSED: Cloud sedimentation
1492  [[maybe_unused]] amrex::Real ised; // ISED: Ice sedimentation
1493  [[maybe_unused]] amrex::Real ssed; // SSED: Snow sedimentation
1494  [[maybe_unused]] amrex::Real gsed; // GSED: Graupel sedimentation
1495  [[maybe_unused]] amrex::Real rsed; // RSED: Rain sedimentation
1496  [[maybe_unused]] amrex::Real tqimelt; // tqimelt: Melting of cloud ice (tendency)
1497 
1498  // NC3DTEN LOCAL ARRAY INITIALIZED
1499  nc3dten(i,j,k) = 0.0;
1500 
1501  // INITIALIZE VARIABLES FOR WRF-CHEM OUTPUT TO ZERO
1502  c2prec = 0.0;
1503  csed = 0.0;
1504  ised = 0.0;
1505  ssed = 0.0;
1506  gsed = 0.0;
1507  rsed = 0.0;
1508 
1509  // LATENT HEAT OF VAPORIZATION
1510  xxlv(i,j,k) = 3.1484E6 - 2370.0 * t3d(i,j,k);
1511  // LATENT HEAT OF SUBLIMATION
1512  xxls(i,j,k) = 3.15E6 - 2370.0 * t3d(i,j,k) + 0.3337E6;
1513 
1514  // Assuming CP is a constant defined elsewhere (specific heat of dry air at constant pressure)
1515  const amrex::Real CP = 1004.5; // J/kg/K
1516  cpm(i,j,k) = CP * (1.0 + 0.887 * qv3d(i,j,k));
1517 
1518  // SATURATION VAPOR PRESSURE AND MIXING RATIO
1519  // hm, add fix for low pressure, 5/12/10
1520  // Assuming POLYSVP is defined elsewhere
1521  evs = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(t3d(i,j,k), 0)); // PA
1522  eis = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(t3d(i,j,k), 1)); // PA
1523  // MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
1524  if (eis > evs) {
1525  eis = evs; // temporary update: adjust ice saturation pressure
1526  }
1527 
1528  // SATURATION MIXING RATIOS
1529  qvs = m_ep_2 * evs / (pres(i,j,k) - evs); // budget equation: calculate water saturation mixing ratio
1530  qvi = m_ep_2 * eis / (pres(i,j,k) - eis); // budget equation: calculate ice saturation mixing ratio
1531 
1532  // SATURATION RATIOS
1533  qvqvs = qv3d(i,j,k) / qvs; // budget equation: calculate water saturation ratio
1534  qvqvsi = qv3d(i,j,k) / qvi; // budget equation: calculate ice saturation ratio
1535 
1536  // AIR DENSITY
1537  rho(i,j,k) = pres(i,j,k) / (m_R * t3d(i,j,k)); // budget equation: calculate air density
1538 
1539  ds0 = 3.0; // Size distribution parameter for snow
1540  di0 = 3.0; // Size distribution parameter for cloud ice
1541  dg0 = 3.0; // Size distribution parameter for graupel
1542  const double CI = 800.0; // Mass-diameter relationship parameter for cloud ice
1543  // ADD NUMBER CONCENTRATION DUE TO CUMULUS TENDENCY
1544  // ASSUME N0 ASSOCIATED WITH CUMULUS PARAM RAIN IS 10^7 M^-4
1545  // ASSUME N0 ASSOCIATED WITH CUMULUS PARAM SNOW IS 2 X 10^7 M^-4
1546  // FOR DETRAINED CLOUD ICE, ASSUME MEAN VOLUME DIAM OF 80 MICRON
1547  if (qrcu1d(i,j,k) >= 1.0e-10) {
1548  dum = 1.8e5 * std::pow(qrcu1d(i,j,k) * dt / (m_pi * m_rhow * std::pow(rho(i,j,k), 3)), 0.25); // rate equation: calculate rain number concentration from cumulus
1549  nr3d(i,j,k) += dum; // budget equation: update rain number concentration
1550  }
1551  if (qscu1d(i,j,k) >= 1.0e-10) {
1552  dum = 3.e5 * std::pow(qscu1d(i,j,k) * dt / (m_cons1 * std::pow(rho(i,j,k), 3)), 1.0 / (ds0 + 1.0)); // rate equation: calculate snow number concentration from cumulus
1553  ns3d(i,j,k) += dum; // budget equation: update snow number concentration
1554  }
1555  if (qicu1d(i,j,k) >= 1.0e-10) {
1556  dum = qicu1d(i,j,k) * dt / (CI * std::pow(80.0e-6, di0)); // rate equation: calculate cloud ice number concentration from cumulus
1557  ni3d(i,j,k) += dum; // budget equation: update cloud ice number concentration
1558  }
1559 
1560  // AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
1561  // hm modify 7/0/09 change limit to 1.e-8
1562  if (qvqvs < 0.9) {
1563  if (qr3d(i,j,k) < 1.0e-8) {
1564  qv3d(i,j,k) += qr3d(i,j,k); // budget equation: transfer rain to vapor
1565  t3d(i,j,k) -= qr3d(i,j,k) * xxlv(i,j,k) / cpm(i,j,k); // budget equation: adjust temperature
1566  qr3d(i,j,k) = 0.0; // temporary update: set rain to zero
1567  }
1568  if (qc3d(i,j,k) < 1.0e-8) {
1569  qv3d(i,j,k) += qc3d(i,j,k); // budget equation: transfer cloud water to vapor
1570  t3d(i,j,k) -= qc3d(i,j,k) * xxlv(i,j,k) / cpm(i,j,k); // budget equation: adjust temperature
1571  qc3d(i,j,k) = 0.0; // temporary update: set cloud water to zero
1572  }
1573  }
1574  if (qvqvsi < 0.9) {
1575  if (qi3d(i,j,k) < 1.0e-8) {
1576  qv3d(i,j,k) += qi3d(i,j,k); // budget equation: transfer cloud ice to vapor
1577  t3d(i,j,k) -= qi3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k); // budget equation: adjust temperature
1578  qi3d(i,j,k) = 0.0; // temporary update: set cloud ice to zero
1579  }
1580  if (qni3d(i,j,k) < 1.0e-8) {
1581  qv3d(i,j,k) += qni3d(i,j,k); // budget equation: transfer snow to vapor
1582  t3d(i,j,k) -= qni3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k); // budget equation: adjust temperature
1583  qni3d(i,j,k) = 0.0; // temporary update: set snow to zero
1584  }
1585  if (qg3d(i,j,k) < 1.0e-8) {
1586  qv3d(i,j,k) += qg3d(i,j,k); // budget equation: transfer graupel to vapor
1587  t3d(i,j,k) -= qg3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k); // budget equation: adjust temperature
1588  qg3d(i,j,k) = 0.0; // temporary update: set graupel to zero
1589  }
1590  }
1591  // HEAT OF FUSION
1592  xlf(i,j,k) = xxls(i,j,k) - xxlv(i,j,k);
1593 
1594  // IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
1595  // Note: QSMALL is not defined in the variable list, so I'll define it
1596  const amrex::Real QSMALL = m_qsmall;
1597 
1598  if (qc3d(i,j,k) < QSMALL) {
1599  qc3d(i,j,k) = 0.0;
1600  nc3d(i,j,k) = 0.0;
1601  effc(i,j,k) = 0.0;
1602  }
1603  if (qr3d(i,j,k) < QSMALL) {
1604  qr3d(i,j,k) = 0.0;
1605  nr3d(i,j,k) = 0.0;
1606  effr(i,j,k) = 0.0;
1607  }
1608  if (qi3d(i,j,k) < QSMALL) {
1609  qi3d(i,j,k) = 0.0;
1610  ni3d(i,j,k) = 0.0;
1611  effi(i,j,k) = 0.0;
1612  }
1613  if (qni3d(i,j,k) < QSMALL) {
1614  qni3d(i,j,k) = 0.0;
1615  ns3d(i,j,k) = 0.0;
1616  effs(i,j,k) = 0.0;
1617  }
1618  if (qg3d(i,j,k) < QSMALL) {
1619  qg3d(i,j,k) = 0.0;
1620  ng3d(i,j,k) = 0.0;
1621  effg(i,j,k) = 0.0;
1622  }
1623  // INITIALIZE SEDIMENTATION TENDENCIES FOR MIXING RATIO
1624  qrsten(i,j,k) = 0.0; // temporary update: initialize QRSTEN
1625  qisten(i,j,k) = 0.0; // temporary update: initialize QISTEN
1626  qnisten(i,j,k) = 0.0; // temporary update: initialize QNISTEN
1627  qcsten(i,j,k) = 0.0; // temporary update: initialize QCSTEN
1628  qgsten(i,j,k) = 0.0; // temporary update: initialize QGSTEN
1629 
1630  // MICROPHYSICS PARAMETERS VARYING IN TIME/HEIGHT
1631  mu(i,j,k) = 1.496e-6 * std::pow(t3d(i,j,k), 1.5) / (t3d(i,j,k) + 120.0); // budget equation: calculate air viscosity
1632 
1633  // Fall speed with density correction (Heymsfield and Benssemer 2006)
1634  dum = std::pow(m_rhosu / rho(i,j,k), 0.54); // temporary update: calculate density correction factor
1635 
1636  // AA revision 4/1/11: Ikawa and Saito 1991 air-density correction
1637  ain(i,j,k) = std::pow(m_rhosu / rho(i,j,k), 0.35) * m_ai; // budget equation: calculate ice fall speed parameter
1638  arn(i,j,k) = dum * m_ar; // budget equation: calculate rain fall speed parameter
1639  asn(i,j,k) = dum * m_as; // budget equation: calculate snow fall speed parameter
1640 
1641  // AA revision 4/1/11: temperature-dependent Stokes fall speed
1642  acn(i,j,k) = m_g * m_rhow / (18.0 * mu(i,j,k)); // budget equation: calculate cloud droplet fall speed parameter
1643 
1644  // HM ADD GRAUPEL 8/28/06
1645  agn(i,j,k) = dum * m_ag; // budget equation: calculate graupel fall speed parameter
1646  // hm 4/7/09 bug fix, initialize lami(i,j,k) to prevent later division by zero
1647  lami(i,j,k) = 0.0; // temporary update: initialize LAMI
1648 
1649  // If there is no cloud/precip water, and if subsaturated, then skip microphysics for this level
1650  bool skipMicrophysics = false;
1651  bool skipConcentrations = false;
1652  if (qc3d(i,j,k) < QSMALL && qi3d(i,j,k) < QSMALL && qni3d(i,j,k) < QSMALL && qr3d(i,j,k) < QSMALL && qg3d(i,j,k) < QSMALL) {
1653  if ((t3d(i,j,k) < 273.15 && qvqvsi < 0.999) || (t3d(i,j,k) >= 273.15 && qvqvs < 0.999)) {
1654  skipMicrophysics = true;// goto label_200;
1655  }
1656  }
1657 
1658  if(!skipMicrophysics) {
1659 
1660  // Thermal conductivity for air
1661  kap = 1.414e3 * mu(i,j,k); // budget equation: calculate thermal conductivity
1662 
1663  // Diffusivity of water vapor
1664  dv = 8.794e-5 * std::pow(t3d(i,j,k), 1.81) / pres(i,j,k); // budget equation: calculate vapor diffusivity
1665 
1666  // Schmidt number
1667  sc_schmidt = mu(i,j,k) / (rho(i,j,k) * dv); // budget equation: calculate Schmidt number
1668 
1669  // Psychometric corrections
1670  // Rate of change sat. mix. ratio with temperature
1671  dum = (m_Rv * std::pow(t3d(i,j,k),2)); // temporary update: calculate temperature factor
1672  dqsdt = xxlv(i,j,k) * qvs / dum; // budget equation: calculate DQSDT
1673  dqsidt = xxls(i,j,k) * qvi / dum; // budget equation: calculate DQSIDT
1674  abi = 1.0 + dqsidt * xxls(i,j,k) / cpm(i,j,k); // budget equation: calculate ABI
1675  ab = 1.0 + dqsdt * xxlv(i,j,k) / cpm(i,j,k); // budget equation: calculate AB
1676 
1677  // CASE FOR TEMPERATURE ABOVE FREEZING
1678  if (t3d(i,j,k) >= 273.15) {
1679  //......................................................................
1680  // ALLOW FOR CONSTANT DROPLET NUMBER
1681  // INUM = 0, PREDICT DROPLET NUMBER
1682  // INUM = 1, SET CONSTANT DROPLET NUMBER
1683 
1684  if (m_inum == 1) {
1685  // CONVERT NDCNST FROM CM-3 TO KG-1
1686  // Note: NDCNST constant would need to be defined elsewhere
1687  nc3d(i,j,k) = m_ndcnst * 1.0e6 / rho(i,j,k); // Set cloud droplet number concentration
1688  }
1689 
1690  // GET SIZE DISTRIBUTION PARAMETERS
1691  // MELT VERY SMALL SNOW AND GRAUPEL MIXING RATIOS, ADD TO RAIN
1692  if (qni3d(i,j,k) < 1.0e-6) {
1693  qr3d(i,j,k) = qr3d(i,j,k) + qni3d(i,j,k); // Transfer snow to rain
1694  nr3d(i,j,k) = nr3d(i,j,k) + ns3d(i,j,k); // Transfer snow number to rain
1695  t3d(i,j,k) = t3d(i,j,k) - qni3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k); // Adjust temperature
1696  qni3d(i,j,k) = 0.0; // Set snow to zero
1697  ns3d(i,j,k) = 0.0; // Set snow number to zero
1698  }
1699 
1700  if (qg3d(i,j,k) < 1.0e-6) {
1701  qr3d(i,j,k) = qr3d(i,j,k) + qg3d(i,j,k); // Transfer graupel to rain
1702  nr3d(i,j,k) = nr3d(i,j,k) + ng3d(i,j,k); // Transfer graupel number to rain
1703  t3d(i,j,k) = t3d(i,j,k) - qg3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k); // Adjust temperature
1704  qg3d(i,j,k) = 0.0; // Set graupel to zero
1705  ng3d(i,j,k) = 0.0; // Set graupel number to zero
1706  }
1707  // Skip to label 300 if concentrations are below thresholds
1708  if (qc3d(i,j,k) < m_qsmall && qni3d(i,j,k) < 1.0e-8 && qr3d(i,j,k) < m_qsmall && qg3d(i,j,k) < 1.0e-8) {
1709  skipConcentrations=true;// goto label_300;
1710  }
1711  if(!skipConcentrations) {
1712  ns3d(i,j,k) = amrex::max(0.0,ns3d(i,j,k));
1713  nc3d(i,j,k) = amrex::max(0.0,nc3d(i,j,k));
1714  nr3d(i,j,k) = amrex::max(0.0,nr3d(i,j,k));
1715  ng3d(i,j,k) = amrex::max(0.0,ng3d(i,j,k));
1716 
1717  // ========================================================================
1718  // USING WRF APPROACH FOR SIZE DISTRIBUTION PARAMETERS
1719  // ========================================================================
1720  // Rain
1721  if (qr3d(i,j,k) >= m_qsmall) {
1722  // Calculate lambda parameter using cons26 (pi*rhow/6)
1723  lamr(i,j,k) = pow(m_pi * m_rhow * nr3d(i,j,k) / qr3d(i,j,k), 1.0/3.0);
1724  n0r(i,j,k) = nr3d(i,j,k)*lamr(i,j,k);
1725 
1726  // Check for slope and adjust vars
1727  if (lamr(i,j,k) < m_lamminr) {
1728  lamr(i,j,k) = m_lamminr;
1729  n0r(i,j,k) = pow(lamr(i,j,k), 4.0) * qr3d(i,j,k) / (m_pi * m_rhow);
1730  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k); // Update number concentration
1731  } else if (lamr(i,j,k) > m_lammaxr) {
1732  lamr(i,j,k) = m_lammaxr;
1733  n0r(i,j,k) = pow(lamr(i,j,k), 4.0) * qr3d(i,j,k) / (m_pi * m_rhow);
1734  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k); // Update number concentration
1735  }
1736  }
1737 
1738  // Cloud droplets
1739  if (qc3d(i,j,k) >= m_qsmall) {
1740  // Calculate air density factor (moist air density)
1741  dum = pres(i,j,k)/(287.15*t3d(i,j,k));
1742 
1743  // MARTIN ET AL. (1994) FORMULA FOR PGAM (WRF implementation)
1744  pgam(i,j,k) = 0.0005714*(nc3d(i,j,k)/1.0e6*dum) + 0.2714;
1745  pgam(i,j,k) = 1.0/(pgam(i,j,k)*pgam(i,j,k)) - 1.0;
1746  pgam(i,j,k) = amrex::max(pgam(i,j,k), 2.0);
1747  pgam(i,j,k) = amrex::min(pgam(i,j,k), 10.0);
1748 
1749  // Calculate gamma function values
1750  amrex::Real gamma_pgam_plus_1 = gamma_function(pgam(i,j,k) + 1.0);
1751  amrex::Real gamma_pgam_plus_4 = gamma_function(pgam(i,j,k) + 4.0);
1752 
1753  // Calculate lambda parameter
1754  lamc(i,j,k) = pow((m_cons26 * nc3d(i,j,k) * gamma_pgam_plus_4) / (qc3d(i,j,k) * gamma_pgam_plus_1), 1.0/3.0);
1755 
1756  // Lambda bounds from WRF - 60 micron max diameter, 1 micron min diameter
1757  amrex::Real lambda_min = (pgam(i,j,k) + 1.0)/60.0e-6;
1758  amrex::Real lambda_max = (pgam(i,j,k) + 1.0)/1.0e-6;
1759 
1760  // Check bounds and update number concentration if needed
1761  if (lamc(i,j,k) < lambda_min) {
1762  lamc(i,j,k) = lambda_min;
1763  // Update cloud droplet number using the same formula as in WRF
1764  nc3d(i,j,k) = exp(3.0*log(lamc(i,j,k)) + log(qc3d(i,j,k)) +
1765  log(gamma_pgam_plus_1) - log(gamma_pgam_plus_4))/ m_cons26;
1766  } else if (lamc(i,j,k) > lambda_max) {
1767  lamc(i,j,k) = lambda_max;
1768  // Update cloud droplet number using the same formula as in WRF
1769  nc3d(i,j,k) = exp(3.0*log(lamc(i,j,k)) + log(qc3d(i,j,k)) +
1770  log(gamma_pgam_plus_1) - log(gamma_pgam_plus_4))/ m_cons26;
1771  }
1772 
1773  // Calculate intercept parameter
1774  cdist1(i,j,k) = nc3d(i,j,k) * pow(lamc(i,j,k), pgam(i,j,k)+1) / gamma_pgam_plus_1;
1775  }
1776 
1777  // Snow
1778  if (qni3d(i,j,k) >= m_qsmall) {
1779  // Calculate lambda parameter
1780  lams(i,j,k) = pow(m_cons1 * ns3d(i,j,k) / qni3d(i,j,k), 1.0/ds0);
1781 
1782  // Calculate intercept parameter
1783  n0s(i,j,k) = ns3d(i,j,k) * lams(i,j,k);
1784 
1785  // Check for slope and adjust vars
1786  if (lams(i,j,k) < m_lammins) {
1787  lams(i,j,k) = m_lammins;
1788  n0s(i,j,k) = pow(lams(i,j,k), 4.0) * qni3d(i,j,k) / m_cons1;
1789  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k); // Update number concentration
1790  } else if (lams(i,j,k) > m_lammaxs) {
1791  lams(i,j,k) = m_lammaxs;
1792  n0s(i,j,k) = pow(lams(i,j,k), 4.0) * qni3d(i,j,k) / m_cons1;
1793  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k); // Update number concentration
1794  }
1795  }
1796 
1797  // Graupel
1798  if (qg3d(i,j,k) >= m_qsmall) {
1799  // Calculate lambda parameter
1800  lamg(i,j,k) = pow(m_cons2 * ng3d(i,j,k) / qg3d(i,j,k), 1.0/dg0);
1801 
1802  // Calculate intercept parameter
1803  n0g(i,j,k) = ng3d(i,j,k) * lamg(i,j,k);
1804 
1805  // Check for slope and adjust vars
1806  if (lamg(i,j,k) < m_lamming) {
1807  lamg(i,j,k) = m_lamming;
1808  n0g(i,j,k) = pow(lamg(i,j,k), 4.0) * qg3d(i,j,k) / m_cons2;
1809  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k); // Update number concentration
1810  } else if (lamg(i,j,k) > m_lammaxg) {
1811  lamg(i,j,k) = m_lammaxg;
1812  n0g(i,j,k) = pow(lamg(i,j,k), 4.0) * qg3d(i,j,k) / m_cons2;
1813  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k); // Update number concentration
1814  }
1815  }
1816  ////////////////////// First instance of ZERO OUT PROCESS RATES
1817  // Zero out process rates
1818  prc = 0.0; // Cloud water to rain conversion rate (PRC)
1819  nprc = 0.0; // Change in cloud droplet number due to autoconversion (NPRC)
1820  nprc1 = 0.0; // Change in rain number due to autoconversion (NPRC1)
1821  pra = 0.0; // Accretion of cloud water by rain (PRA)
1822  npra = 0.0; // Change in cloud droplet number due to accretion by rain (NPRA)
1823  nragg = 0.0; // Self-collection/breakup of rain (NRAGG)
1824  nsmlts = 0.0; // Loss of snow number during melting (NSMLTS)
1825  nsmltr = 0.0; // Change in rain number due to snow melting (NSMLTR)
1826  evpms = 0.0; // Melting snow evaporation rate (EVPMS)
1827  pcc = 0.0; // Condensation/evaporation of cloud water (PCC)
1828  pre = 0.0; // Evaporation of rain (PRE)
1829  nsubc = 0.0; // Loss of cloud droplet number during evaporation (NSUBC)
1830  nsubr = 0.0; // Loss of rain number during evaporation (NSUBR)
1831  pracg = 0.0; // Collection of rain by graupel (PRACG)
1832  npracg = 0.0; // Change in number due to collection of rain by graupel (NPRACG)
1833  psmlt = 0.0; // Melting of snow (PSMLT)
1834  pgmlt = 0.0; // Melting of graupel (PGMLT)
1835  evpmg = 0.0; // Evaporation of melting graupel (EVPMG)
1836  pracs = 0.0; // Collection of snow by rain (PRACS)
1837  npracs = 0.0; // Change in number due to collection of snow by rain (NPRACS)
1838  ngmltg = 0.0; // Loss of graupel number during melting (NGMLTG)
1839  ngmltr = 0.0; // Change in rain number due to graupel melting (NGMLTR)
1840 
1841  // CALCULATION OF MICROPHYSICAL PROCESS RATES, T > 273.15 K
1842 
1843  // AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
1844  // FORMULA FROM BEHENG (1994)
1845  // USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
1846  // AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
1847  // AS A GAMMA DISTRIBUTION
1848 
1849  // USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
1850 
1851  if (qc3d(i,j,k) >= 1.0e-6) {
1852  // HM ADD 12/13/06, REPLACE WITH NEWER FORMULA
1853  // FROM KHAIROUTDINOV AND KOGAN 2000, MWR
1854  prc = 1350.0 * std::pow(qc3d(i,j,k), 2.47) *
1855  std::pow((nc3d(i,j,k)/1.0e6*rho(i,j,k)), -1.79);
1856 
1857  // note: nprc1 is change in Nr,
1858  // nprc is change in Nc
1859  nprc1 = prc / m_cons29;
1860  nprc = prc / (qc3d(i,j,k) / nc3d(i,j,k));
1861 
1862  // hm bug fix 3/20/12
1863  nprc = std::min(nprc, nc3d(i,j,k) / dt);
1864  nprc1 = std::min(nprc1, nprc);
1865  }
1866 
1867  // HM ADD 12/13/06, COLLECTION OF SNOW BY RAIN ABOVE FREEZING
1868  // FORMULA FROM IKAWA AND SAITO (1991)
1869 
1870  if (qr3d(i,j,k) >= 1.0e-8 && qni3d(i,j,k) >= 1.0e-8) {
1871  amrex::Real ums_local = asn(i,j,k) * m_cons3 / std::pow(lams(i,j,k), m_bs);
1872  amrex::Real umr_local = arn(i,j,k) * m_cons4 / std::pow(lamr(i,j,k), m_br);
1873  amrex::Real uns_local = asn(i,j,k) * m_cons5 / std::pow(lams(i,j,k), m_bs);
1874  amrex::Real unr_local = arn(i,j,k) * m_cons6 / std::pow(lamr(i,j,k), m_br);
1875 
1876  // SET REALISTIC LIMITS ON FALLSPEEDS
1877  // bug fix, 10/08/09
1878  dum = std::pow(m_rhosu/rho(i,j,k), 0.54);
1879  ums_local = std::min(ums_local, 1.2*dum);
1880  uns_local = std::min(uns_local, 1.2*dum);
1881  umr_local = std::min(umr_local, 9.1*dum);
1882  unr_local = std::min(unr_local, 9.1*dum);
1883 
1884 
1885  // hm fix, 2/12/13
1886  // for above freezing conditions to get accelerated melting of snow,
1887  // we need collection of rain by snow (following Lin et al. 1983)
1888  ////////////////////////Might need pow expanding
1889  pracs = m_cons41 * (std::sqrt(std::pow(1.2*umr_local-0.95*ums_local, 2) +
1890  0.08*ums_local*umr_local) * rho(i,j,k) *
1891  n0r(i,j,k) * n0s(i,j,k) / std::pow(lamr(i,j,k), 3) *
1892  (5.0/(std::pow(lamr(i,j,k), 3) * lams(i,j,k)) +
1893  2.0/(std::pow(lamr(i,j,k), 2) * std::pow(lams(i,j,k), 2)) +
1894  0.5/(lamr(i,j,k) * std::pow(lams(i,j,k), 3))));
1895  }
1896  // ADD COLLECTION OF GRAUPEL BY RAIN ABOVE FREEZING
1897  // ASSUME ALL RAIN COLLECTION BY GRAUPEL ABOVE FREEZING IS SHED
1898  // ASSUME SHED DROPS ARE 1 MM IN SIZE
1899 
1900  if (qr3d(i,j,k) >= 1.0e-8 && qg3d(i,j,k) >= 1.0e-8) {
1901 
1902  amrex::Real umg_local = agn(i,j,k) * m_cons7 / std::pow(lamg(i,j,k), m_bg);
1903  amrex::Real umr_local = arn(i,j,k) * m_cons4 / std::pow(lamr(i,j,k), m_br);
1904  amrex::Real ung_local = agn(i,j,k) * m_cons8 / std::pow(lamg(i,j,k), m_bg);
1905  amrex::Real unr_local = arn(i,j,k) * m_cons6 / std::pow(lamr(i,j,k), m_br);
1906 
1907  // SET REALISTIC LIMITS ON FALLSPEEDS
1908  // bug fix, 10/08/09
1909  dum = std::pow(m_rhosu/rho(i,j,k), 0.54);
1910  umg_local = std::min(umg_local, 20.0*dum);
1911  ung_local = std::min(ung_local, 20.0*dum);
1912  umr_local = std::min(umr_local, 9.1*dum);
1913  unr_local = std::min(unr_local, 9.1*dum);
1914 
1915  // PRACG IS MIXING RATIO OF RAIN PER SEC COLLECTED BY GRAUPEL/HAIL
1916  pracg = m_cons41 * (std::sqrt(std::pow(1.2*umr_local-0.95*umg_local, 2) +
1917  0.08*umg_local*umr_local) * rho(i,j,k) *
1918  n0r(i,j,k) * n0g(i,j,k) / std::pow(lamr(i,j,k), 3) *
1919  (5.0/(std::pow(lamr(i,j,k), 3) * lamg(i,j,k)) +
1920  2.0/(std::pow(lamr(i,j,k), 2) * std::pow(lamg(i,j,k), 2)) +
1921  0.5/(lamr(i,j,k) * std::pow(lamg(i,j,k), 3))));
1922 
1923  // ASSUME 1 MM DROPS ARE SHED, GET NUMBER SHED PER SEC
1924  dum = pracg/5.2e-7;
1925 
1926  npracg = m_cons32 * rho(i,j,k) * (std::sqrt(1.7*std::pow(unr_local-ung_local, 2) +
1927  0.3*unr_local*ung_local) * n0r(i,j,k) * n0g(i,j,k) *
1928  (1.0/(std::pow(lamr(i,j,k), 3) * lamg(i,j,k)) +
1929  1.0/(std::pow(lamr(i,j,k), 2) * std::pow(lamg(i,j,k), 2)) +
1930  1.0/(lamr(i,j,k) * std::pow(lamg(i,j,k), 3))));
1931  // hm 7/15/13, remove limit so that the number of collected drops can smaller than
1932  // number of shed drops
1933  npracg = npracg - dum;
1934  }
1935  // ACCRETION OF CLOUD LIQUID WATER BY RAIN
1936  // CONTINUOUS COLLECTION EQUATION WITH
1937  // GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
1938 
1939  if (qr3d(i,j,k) >= 1.0e-8 && qc3d(i,j,k) >= 1.0e-8) {
1940  // 12/13/06 HM ADD, REPLACE WITH NEWER FORMULA FROM
1941  // KHAIROUTDINOV AND KOGAN 2000, MWR
1942  dum = qc3d(i,j,k) * qr3d(i,j,k);
1943  pra = 67.0 * std::pow(dum, 1.15);
1944  npra = pra / (qc3d(i,j,k) / nc3d(i,j,k));
1945  }
1946 
1947  // SELF-COLLECTION OF RAIN DROPS
1948  // FROM BEHENG(1994)
1949  // FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
1950  // AS DESCRIBED ABOVE FOR AUTOCONVERSION
1951 
1952  if (qr3d(i,j,k) >= 1.0e-8) {
1953  // include breakup add 10/09/09
1954  dum1 = 300.0e-6;
1955  if (1.0/lamr(i,j,k) < dum1) {
1956  dum = 1.0;
1957  } else {
1958  dum = 2.0 - std::exp(2300.0 * (1.0/lamr(i,j,k) - dum1));
1959  }
1960  nragg = -5.78 * dum * nr3d(i,j,k) * qr3d(i,j,k) * rho(i,j,k);
1961  }
1962  // CALCULATE EVAP OF RAIN (RUTLEDGE AND HOBBS 1983)
1963  if (qr3d(i,j,k) >= m_qsmall) {
1964  epsr = 2.0 * m_pi * n0r(i,j,k) * rho(i,j,k) * dv *
1965  (m_f1r/(lamr(i,j,k)*lamr(i,j,k)) +
1966  m_f2r * std::sqrt(arn(i,j,k)*rho(i,j,k)/mu(i,j,k)) *
1967  std::pow(sc_schmidt, 1.0/3.0) * m_cons9 /
1968  std::pow(lamr(i,j,k), m_cons34));
1969  } else {
1970  epsr = 0.0;
1971  }
1972  // NO CONDENSATION ONTO RAIN, ONLY EVAP ALLOWED
1973  if (qv3d(i,j,k) < qvs) {
1974  pre = epsr * (qv3d(i,j,k) - qvs) / ab;
1975  pre = std::min(pre, 0.0);
1976  } else {
1977  pre = 0.0;
1978  }
1979  // MELTING OF SNOW
1980  // SNOW MAY PERSIST ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
1981  // IF WATER SUPERSATURATION, SNOW MELTS TO FORM RAIN
1982 
1983  if (qni3d(i,j,k) >= 1.0e-8) {
1984  // fix 053011
1985  // HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
1986  dum = -m_cpw/xlf(i,j,k) * (t3d(i,j,k) - 273.15) * pracs;
1987 
1988  // hm fix 1/20/15
1989  psmlt = 2.0 * m_pi * n0s(i,j,k) * kap * (273.15 - t3d(i,j,k)) /
1990  xlf(i,j,k) * (m_f1s/(lams(i,j,k)*lams(i,j,k)) +
1991  m_f2s * std::sqrt(asn(i,j,k)*rho(i,j,k)/mu(i,j,k)) *
1992  std::pow(sc_schmidt, 1.0/3.0) * m_cons10 /
1993  std::pow(lams(i,j,k), m_cons35)) + dum;
1994 
1995  // IN WATER SUBSATURATION, SNOW MELTS AND EVAPORATES
1996  if (qvqvs < 1.0) {
1997  epss = 2.0 * m_pi * n0s(i,j,k) * rho(i,j,k) * dv *
1998  (m_f1s/(lams(i,j,k)*lams(i,j,k)) +
1999  m_f2s * std::sqrt(asn(i,j,k)*rho(i,j,k)/mu(i,j,k)) *
2000  std::pow(sc_schmidt, 1.0/3.0) * m_cons10 /
2001  std::pow(lams(i,j,k), m_cons35));
2002 
2003  // hm fix 8/4/08
2004  evpms = (qv3d(i,j,k) - qvs) * epss / ab;
2005  evpms = std::max(evpms, psmlt);
2006  psmlt = psmlt - evpms;
2007  }
2008  }
2009  // MELTING OF GRAUPEL
2010  // GRAUPEL MAY PERSIST ABOVE FREEZING, FORMULA FROM RUTLEDGE AND HOBBS, 1984
2011  // IF WATER SUPERSATURATION, GRAUPEL MELTS TO FORM RAIN
2012 
2013  if (qg3d(i,j,k) >= 1.0e-8) {
2014  // fix 053011
2015  // HM, MODIFY FOR V3.2, ADD ACCELERATED MELTING DUE TO COLLISION WITH RAIN
2016 
2017  dum = -m_cpw/xlf(i,j,k) * (t3d(i,j,k) - 273.15) * pracg;
2018 
2019  // hm fix 1/20/15
2020  pgmlt = 2.0 * m_pi * n0g(i,j,k) * kap * (273.15 - t3d(i,j,k)) /
2021  xlf(i,j,k) * (m_f1s/(lamg(i,j,k)*lamg(i,j,k)) +
2022  m_f2s * std::sqrt(agn(i,j,k)*rho(i,j,k)/mu(i,j,k)) *
2023  std::pow(sc_schmidt, 1.0/3.0) * m_cons11 /
2024  std::pow(lamg(i,j,k), m_cons36)) + dum;
2025 
2026  // IN WATER SUBSATURATION, GRAUPEL MELTS AND EVAPORATES
2027  if (qvqvs < 1.0) {
2028  epsg = 2.0 * m_pi * n0g(i,j,k) * rho(i,j,k) * dv *
2029  (m_f1s/(lamg(i,j,k)*lamg(i,j,k)) +
2030  m_f2s * std::sqrt(agn(i,j,k)*rho(i,j,k)/mu(i,j,k)) *
2031  std::pow(sc_schmidt, 1.0/3.0) * m_cons11 /
2032  std::pow(lamg(i,j,k), m_cons36));
2033 
2034  // hm fix 8/4/08
2035  evpmg = (qv3d(i,j,k) - qvs) * epsg / ab;
2036  evpmg = std::max(evpmg, pgmlt);
2037  pgmlt = pgmlt - evpmg;
2038  }
2039  }
2040  // HM, V3.2
2041  // RESET PRACG AND PRACS TO ZERO, THIS IS DONE BECAUSE THERE IS NO
2042  // TRANSFER OF MASS FROM SNOW AND GRAUPEL TO RAIN DIRECTLY FROM COLLECTION
2043  // ABOVE FREEZING, IT IS ONLY USED FOR ENHANCEMENT OF MELTING AND SHEDDING
2044 
2045  pracg = 0.0;
2046  pracs = 0.0;
2047  // CONSERVATION OF QC
2048  dum = (prc + pra) * dt;
2049 
2050  if (dum > qc3d(i,j,k) && qc3d(i,j,k) >= m_qsmall) {
2051  ratio = qc3d(i,j,k) / dum;
2052  prc = prc * ratio;
2053  pra = pra * ratio;
2054  }
2055 
2056  // CONSERVATION OF SNOW
2057  dum = (-psmlt - evpms + pracs) * dt;
2058 
2059  if (dum > qni3d(i,j,k) && qni3d(i,j,k) >= m_qsmall) {
2060  // NO SOURCE TERMS FOR SNOW AT T > FREEZING
2061  ratio = qni3d(i,j,k) / dum;
2062  psmlt = psmlt * ratio;
2063  evpms = evpms * ratio;
2064  pracs = pracs * ratio;
2065  }
2066 
2067  // CONSERVATION OF GRAUPEL
2068  dum = (-pgmlt - evpmg + pracg) * dt;
2069 
2070  if (dum > qg3d(i,j,k) && qg3d(i,j,k) >= m_qsmall) {
2071  // NO SOURCE TERM FOR GRAUPEL ABOVE FREEZING
2072  ratio = qg3d(i,j,k) / dum;
2073  pgmlt = pgmlt * ratio;
2074  evpmg = evpmg * ratio;
2075  pracg = pracg * ratio;
2076  }
2077 
2078  // CONSERVATION OF QR
2079  // HM 12/13/06, ADDED CONSERVATION OF RAIN SINCE PRE IS NEGATIVE
2080 
2081  dum = (-pracs - pracg - pre - pra - prc + psmlt + pgmlt) * dt;
2082 
2083  if (dum > qr3d(i,j,k) && qr3d(i,j,k) >= m_qsmall) {
2084  ratio = (qr3d(i,j,k)/dt + pracs + pracg + pra + prc - psmlt - pgmlt) / (-pre);
2085  pre = pre * ratio;
2086  }
2087  // Update tendencies
2088  qv3dten(i,j,k) = qv3dten(i,j,k) + (-pre - evpms - evpmg);
2089 
2090  t3dten(i,j,k) = t3dten(i,j,k) + (pre * xxlv(i,j,k) +
2091  (evpms + evpmg) * xxls(i,j,k) +
2092  (psmlt + pgmlt - pracs - pracg) * xlf(i,j,k)) / cpm(i,j,k);
2093 
2094  qc3dten(i,j,k) = qc3dten(i,j,k) + (-pra - prc);
2095  qr3dten(i,j,k) = qr3dten(i,j,k) + (pre + pra + prc - psmlt - pgmlt + pracs + pracg);
2096  qni3dten(i,j,k) = qni3dten(i,j,k) + (psmlt + evpms - pracs);
2097  qg3dten(i,j,k) = qg3dten(i,j,k) + (pgmlt + evpmg - pracg);
2098 
2099  // fix 053011
2100  // HM, bug fix 5/12/08, npracg is subtracted from nr not ng
2101  nc3dten(i,j,k) = nc3dten(i,j,k) + (-npra - nprc);
2102  nr3dten(i,j,k) = nr3dten(i,j,k) + (nprc1 + nragg - npracg);
2103 
2104  // HM ADD, WRF-CHEM, ADD TENDENCIES FOR C2PREC
2105  c2prec = pra + prc;
2106 
2107  if (pre < 0.0) {
2108  dum = pre * dt / qr3d(i,j,k);
2109  dum = std::max(-1.0, dum);
2110  nsubr = dum * nr3d(i,j,k) / dt;
2111  }
2112 
2113  if (evpms + psmlt < 0.0) {
2114  dum = (evpms + psmlt) * dt / qni3d(i,j,k);
2115  dum = std::max(-1.0, dum);
2116  nsmlts = dum * ns3d(i,j,k) / dt;
2117  }
2118 
2119  if (psmlt < 0.0) {
2120  dum = psmlt * dt / qni3d(i,j,k);
2121  dum = std::max(-1.0, dum);
2122  nsmltr = dum * ns3d(i,j,k) / dt;
2123  }
2124 
2125  if (evpmg + pgmlt < 0.0) {
2126  dum = (evpmg + pgmlt) * dt / qg3d(i,j,k);
2127  dum = std::max(-1.0, dum);
2128  ngmltg = dum * ng3d(i,j,k) / dt;
2129  }
2130 
2131  if (pgmlt < 0.0) {
2132  dum = pgmlt * dt / qg3d(i,j,k);
2133  dum = std::max(-1.0, dum);
2134  ngmltr = dum * ng3d(i,j,k) / dt;
2135  }
2136 
2137  ns3dten(i,j,k) = ns3dten(i,j,k) + nsmlts;
2138  ng3dten(i,j,k) = ng3dten(i,j,k) + ngmltg;
2139  nr3dten(i,j,k) = nr3dten(i,j,k) + (nsubr - nsmltr - ngmltr);
2140 
2141  }
2142  //Right after 300 CONTINUE
2143 // label_300:
2144  // Calculate saturation adjustment to condense extra vapor above water saturation
2145  dumt = t3d(i,j,k) + dt * t3dten(i,j,k);
2146  dumqv = qv3d(i,j,k) + dt * qv3dten(i,j,k);
2147 
2148  // Fix for low pressure (added 5/12/10)
2149  dum = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(dumt, 0));
2150  dumqss = m_ep_2 * dum / (pres(i,j,k) - dum);
2151  dumqc = qc3d(i,j,k) + dt * qc3dten(i,j,k);
2152  dumqc = std::max(dumqc, 0.0);
2153 
2154  // Saturation adjustment for liquid
2155  dums = dumqv - dumqss;
2156  pcc = dums / (1.0 + std::pow(xxlv(i,j,k), 2) * dumqss / (cpm(i,j,k) * m_Rv * std::pow(dumt, 2))) / dt;
2157  if (pcc * dt + dumqc < 0.0) {
2158  pcc = -dumqc / dt;
2159  }
2160 
2161  // Update tendencies
2162  qv3dten(i,j,k) -= pcc;
2163  t3dten(i,j,k) += pcc * xxlv(i,j,k) / cpm(i,j,k);
2164  qc3dten(i,j,k) += pcc;
2165  } else { //cold
2166  //......................................................................
2167  // ALLOW FOR CONSTANT DROPLET NUMBER
2168  // INUM = 0, PREDICT DROPLET NUMBER
2169  // INUM = 1, SET CONSTANT DROPLET NUMBER
2170 
2171  if (m_inum == 1) {
2172  // CONVERT NDCNST FROM CM-3 TO KG-1
2173  // Note: NDCNST constant would need to be defined elsewhere
2174  nc3d(i,j,k) = m_ndcnst * 1.0e6 / rho(i,j,k); // Set cloud droplet number concentration
2175  }
2176 
2177  ni3d(i,j,k) = amrex::max(0.0,ni3d(i,j,k));
2178  ns3d(i,j,k) = amrex::max(0.0,ns3d(i,j,k));
2179  nc3d(i,j,k) = amrex::max(0.0,nc3d(i,j,k));
2180  nr3d(i,j,k) = amrex::max(0.0,nr3d(i,j,k));
2181  ng3d(i,j,k) = amrex::max(0.0,ng3d(i,j,k));
2182 
2183  // ========================================================================
2184  // USING WRF APPROACH FOR SIZE DISTRIBUTION PARAMETERS
2185  // ========================================================================
2186  // Rain
2187  if (qr3d(i,j,k) >= m_qsmall) {
2188  // Calculate lambda parameter using cons26 (pi*rhow/6)
2189  lamr(i,j,k) = pow(m_pi * m_rhow * nr3d(i,j,k) / qr3d(i,j,k), 1.0/3.0);
2190 
2191  // Check for slope and adjust vars
2192  if (lamr(i,j,k) < m_lamminr) {
2193  lamr(i,j,k) = m_lamminr;
2194  n0r(i,j,k) = pow(lamr(i,j,k), 4.0) * qr3d(i,j,k) / (m_pi * m_rhow);
2195  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k); // Update number concentration
2196  } else if (lamr(i,j,k) > m_lammaxr) {
2197  lamr(i,j,k) = m_lammaxr;
2198  n0r(i,j,k) = pow(lamr(i,j,k), 4.0) * qr3d(i,j,k) / (m_pi * m_rhow);
2199  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k); // Update number concentration
2200  } else {
2201  // Calculate intercept parameter using WRF formula
2202  n0r(i,j,k) = pow(lamr(i,j,k), 4.0) * qr3d(i,j,k) / (m_pi * m_rhow);
2203  }
2204  }
2205 
2206 
2207  // Cloud droplets
2208  if (qc3d(i,j,k) >= m_qsmall) {
2209  // Calculate air density factor (moist air density)
2210  dum = pres(i,j,k)/(287.15*t3d(i,j,k));
2211 
2212  // MARTIN ET AL. (1994) FORMULA FOR PGAM (WRF implementation)
2213  pgam(i,j,k) = 0.0005714*(nc3d(i,j,k)/1.0e6*dum) + 0.2714;
2214  pgam(i,j,k) = 1.0/(std::pow(pgam(i,j,k), 2)) - 1.0;
2215  pgam(i,j,k) = amrex::max(pgam(i,j,k), 2.0);
2216  pgam(i,j,k) = amrex::min(pgam(i,j,k), 10.0);
2217 
2218  // Calculate gamma function values
2219  amrex::Real gamma_pgam_plus_1 = gamma_function(pgam(i,j,k) + 1.0);
2220  amrex::Real gamma_pgam_plus_4 = gamma_function(pgam(i,j,k) + 4.0);
2221 
2222  // Calculate lambda parameter
2223  lamc(i,j,k) = pow((m_cons26 * nc3d(i,j,k) * gamma_pgam_plus_4) / (qc3d(i,j,k) * gamma_pgam_plus_1), 1.0/3.0);
2224 
2225  // Lambda bounds from WRF - 60 micron max diameter, 1 micron min diameter
2226  amrex::Real lambda_min = (pgam(i,j,k) + 1.0)/60.0e-6;
2227  amrex::Real lambda_max = (pgam(i,j,k) + 1.0)/1.0e-6;
2228 
2229  // Check bounds and update number concentration if needed
2230  if (lamc(i,j,k) < lambda_min) {
2231  lamc(i,j,k) = lambda_min;
2232  // Update cloud droplet number using the same formula as in WRF
2233  nc3d(i,j,k) = exp(3.0*log(lamc(i,j,k)) + log(qc3d(i,j,k)) +
2234  log(gamma_pgam_plus_1) - log(gamma_pgam_plus_4))/ m_cons26;
2235  } else if (lamc(i,j,k) > lambda_max) {
2236  lamc(i,j,k) = lambda_max;
2237  // Update cloud droplet number using the same formula as in WRF
2238  nc3d(i,j,k) = exp(3.0*log(lamc(i,j,k)) + log(qc3d(i,j,k)) +
2239  log(gamma_pgam_plus_1) - log(gamma_pgam_plus_4))/ m_cons26;
2240  }
2241 
2242  // Calculate intercept parameter
2243  cdist1(i,j,k) = nc3d(i,j,k) / gamma_pgam_plus_1;
2244  }
2245 
2246  // Snow
2247  if (qni3d(i,j,k) >= m_qsmall) {
2248  // Calculate lambda parameter
2249  lams(i,j,k) = pow(m_cons1 * ns3d(i,j,k) / qni3d(i,j,k), 1.0/ds0);
2250 
2251  // Calculate intercept parameter
2252  n0s(i,j,k) = ns3d(i,j,k) * lams(i,j,k);
2253 
2254  // Check for slope and adjust vars
2255  if (lams(i,j,k) < m_lammins) {
2256  lams(i,j,k) = m_lammins;
2257  n0s(i,j,k) = pow(lams(i,j,k), 4.0) * qni3d(i,j,k) / m_cons1;
2258  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k); // Update number concentration
2259  } else if (lams(i,j,k) > m_lammaxs) {
2260  lams(i,j,k) = m_lammaxs;
2261  n0s(i,j,k) = pow(lams(i,j,k), 4.0) * qni3d(i,j,k) / m_cons1;
2262  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k); // Update number concentration
2263  }
2264  }
2265 
2266  // Cloud ice
2267  if (qi3d(i,j,k) >= m_qsmall) {
2268  // Calculate lambda parameter
2269  lami(i,j,k) = pow(m_cons12 * ni3d(i,j,k) / qi3d(i,j,k), 1.0/3.0);
2270 
2271  // Calculate intercept parameter (initial calculation)
2272  n0i(i,j,k) = ni3d(i,j,k) * lami(i,j,k);
2273 
2274  // Check for slope (apply bounds)
2275  if (lami(i,j,k) < m_lammini) {
2276  lami(i,j,k) = m_lammini;
2277  // Recalculate n0i(i,j,k) when lambda is adjusted
2278  n0i(i,j,k) = pow(lami(i,j,k), 4.0) * qi3d(i,j,k) / m_cons12;
2279  // Update ni3d when lambda is adjusted
2280  ni3d(i,j,k) = n0i(i,j,k) / lami(i,j,k);
2281  } else if (lami(i,j,k) > m_lammaxi) {
2282  lami(i,j,k) = m_lammaxi;
2283  // Recalculate n0i(i,j,k) when lambda is adjusted
2284  n0i(i,j,k) = pow(lami(i,j,k), 4.0) * qi3d(i,j,k) / m_cons12;
2285  // Update ni3d when lambda is adjusted
2286  ni3d(i,j,k) = n0i(i,j,k) / lami(i,j,k);
2287  }
2288  }
2289  // Graupel
2290  if (qg3d(i,j,k) >= m_qsmall) {
2291  // Calculate lambda parameter
2292  lamg(i,j,k) = pow(m_cons2 * ng3d(i,j,k) / qg3d(i,j,k), 1.0/dg0);
2293 
2294  // Calculate intercept parameter
2295  n0g(i,j,k) = ng3d(i,j,k) * lamg(i,j,k);
2296 
2297  // Check for slope and adjust vars
2298  if (lamg(i,j,k) < m_lamming) {
2299  lamg(i,j,k) = m_lamming;
2300  n0g(i,j,k) = pow(lamg(i,j,k), 4.0) * qg3d(i,j,k) / m_cons2;
2301  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k); // Update number concentration
2302  } else if (lamg(i,j,k) > m_lammaxg) {
2303  lamg(i,j,k) = m_lammaxg;
2304  n0g(i,j,k) = pow(lamg(i,j,k), 4.0) * qg3d(i,j,k) / m_cons2;
2305  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k); // Update number concentration
2306  }
2307  }
2308  ////////////////////// Second instance of ZERO OUT PROCESS RATES
2309  // Zero out process rates
2310  mnuccc = 0.0; // Change Q due to contact freezing droplets (MNUCCC)
2311  nnuccc = 0.0; // Change N due to contact freezing droplets (NNUCCC)
2312  prc = 0.0; // Autoconversion droplets (PRC)
2313  nprc = 0.0; // Change NC autoconversion droplets (NPRC)
2314  nprc1 = 0.0; // Change NR autoconversion droplets (NPRC1)
2315  nsagg = 0.0; // Self-collection of snow (NSAGG)
2316  psacws = 0.0; // Change Q droplet accretion by snow (PSACWS)
2317  npsacws = 0.0; // Change N droplet accretion by snow (NPSACWS)
2318  psacwi = 0.0; // Change Q droplet accretion by cloud ice (PSACWI)
2319  npsacwi = 0.0; // Change N droplet accretion by cloud ice (NPSACWI)
2320  pracs = 0.0; // Change Q rain-snow collection (PRACS)
2321  npracs = 0.0; // Change N rain-snow collection (NPRACS)
2322  nmults = 0.0; // Ice multiplication due to riming droplets by snow (NMULTS)
2323  qmults = 0.0; // Change Q due to ice multiplication droplets/snow (QMULTS)
2324  nmultr = 0.0; // Ice multiplication due to riming rain by snow (NMULTR)
2325  qmultr = 0.0; // Change Q due to ice multiplication rain/snow (QMULTR)
2326  nmultg = 0.0; // Ice multiplication due to accretion droplets by graupel (NMULTG)
2327  qmultg = 0.0; // Change Q due to ice multiplication droplets/graupel (QMULTG)
2328  nmultrg = 0.0; // Ice multiplication due to accretion rain by graupel (NMULTRG)
2329  qmultrg = 0.0; // Change Q due to ice multiplication rain/graupel (QMULTRG)
2330  mnuccr = 0.0; // Change Q due to contact freezing rain (MNUCCR)
2331  nnuccr = 0.0; // Change N due to contact freezing rain (NNUCCR)
2332  pra = 0.0; // Accretion droplets by rain (PRA)
2333  npra = 0.0; // Change N due to droplet accretion by rain (NPRA)
2334  nragg = 0.0; // Self-collection/breakup of rain (NRAGG)
2335  prci = 0.0; // Change Q autoconversion cloud ice to snow (PRCI)
2336  nprci = 0.0; // Change N autoconversion cloud ice by snow (NPRCI)
2337  prai = 0.0; // Change Q accretion cloud ice by snow (PRAI)
2338  nprai = 0.0; // Change N accretion cloud ice (NPRAI)
2339  nnuccd = 0.0; // Change N freezing aerosol (primary ice nucleation) (NNUCCD)
2340  mnuccd = 0.0; // Change Q freezing aerosol (primary ice nucleation) (MNUCCD)
2341  pcc = 0.0; // Condensation/evaporation droplets (PCC)
2342  pre = 0.0; // Evaporation of rain (PRE)
2343  prd = 0.0; // Deposition cloud ice (PRD)
2344  prds = 0.0; // Deposition snow (PRDS)
2345  eprd = 0.0; // Sublimation cloud ice (EPRD)
2346  eprds = 0.0; // Sublimation snow (EPRDS)
2347  nsubc = 0.0; // Loss of NC during evaporation (NSUBC)
2348  nsubi = 0.0; // Loss of NI during sublimation (NSUBI)
2349  nsubs = 0.0; // Loss of NS during sublimation (NSUBS)
2350  nsubr = 0.0; // Loss of NR during evaporation (NSUBR)
2351  piacr = 0.0; // Change QR, ice-rain collection (PIACR)
2352  niacr = 0.0; // Change N, ice-rain collection (NIACR)
2353  praci = 0.0; // Change QI, ice-rain collection (PRACI)
2354  piacrs = 0.0; // Change QR, ice rain collision, added to snow (PIACRS)
2355  niacrs = 0.0; // Change N, ice rain collision, added to snow (NIACRS)
2356  pracis = 0.0; // Change QI, ice rain collision, added to snow (PRACIS)
2357 
2358  // Graupel processes
2359  pracg = 0.0; // Change in Q collection rain by graupel (PRACG)
2360  psacr = 0.0; // Conversion due to collection of snow by rain (PSACR)
2361  psacwg = 0.0; // Change in Q collection droplets by graupel (PSACWG)
2362  pgsacw = 0.0; // Conversion Q to graupel due to collection droplets by snow (PGSACW)
2363  pgracs = 0.0; // Conversion Q to graupel due to collection rain by snow (PGRACS)
2364  prdg = 0.0; // Deposition of graupel (PRDG)
2365  eprdg = 0.0; // Sublimation of graupel (EPRDG)
2366  npracg = 0.0; // Change N collection rain by graupel (NPRACG)
2367  npsacwg = 0.0; // Change N collection droplets by graupel (NPSACWG)
2368  nscng = 0.0; // Change N conversion to graupel due to collection droplets by snow (NSCNG)
2369  ngracs = 0.0; // Change N conversion to graupel due to collection rain by snow (NGRACS)
2370  nsubg = 0.0; // Change N sublimation/deposition of graupel (NSUBG)
2371 
2372  ////////////////////// CALCULATION OF MICROPHYSICAL PROCESS RATES
2373  // FREEZING OF CLOUD DROPLETS - ONLY ALLOWED BELOW -4C
2374  if (qc3d(i,j,k) >= m_qsmall && t3d(i,j,k) < 269.15) {
2375  // NUMBER OF CONTACT NUCLEI (M^-3) FROM MEYERS ET AL., 1992
2376  // FACTOR OF 1000 IS TO CONVERT FROM L^-1 TO M^-3
2377  // MEYERS CURVE
2378  amrex::Real nacnt = std::exp(-2.80 + 0.262 * (273.15 - t3d(i,j,k))) * 1000.0;
2379 
2380  // MEAN FREE PATH
2381  dum = 7.37 * t3d(i,j,k) / (288.0 * 10.0 * pres(i,j,k)) / 100.0;
2382 
2383  // EFFECTIVE DIFFUSIVITY OF CONTACT NUCLEI
2384  // BASED ON BROWNIAN DIFFUSION
2385  amrex::Real dap = m_cons37 * t3d(i,j,k) * (1.0 + dum / m_rin) / mu(i,j,k);
2386 
2387  // CONTACT FREEZING
2388  mnuccc = m_cons38 * dap * nacnt * std::exp(std::log(cdist1(i,j,k)) +
2389  std::log(gamma_function(pgam(i,j,k) + 5.0)) - 4.0 * std::log(lamc(i,j,k)));
2390  nnuccc = 2.0 * m_pi * dap * nacnt * cdist1(i,j,k) *
2391  gamma_function(pgam(i,j,k) + 2.0) / lamc(i,j,k);
2392 
2393  // IMMERSION FREEZING (BIGG 1953)
2394  // hm 7/15/13 fix for consistency w/ original formula
2395  mnuccc = mnuccc + m_cons39 *
2396  std::exp(std::log(cdist1(i,j,k)) + std::log(gamma_function(7.0 + pgam(i,j,k))) - 6.0 * std::log(lamc(i,j,k))) *
2397  (std::exp(m_aimm * (273.15 - t3d(i,j,k))) - 1.0);
2398 
2399  nnuccc = nnuccc +
2400  m_cons40 * std::exp(std::log(cdist1(i,j,k)) + std::log(gamma_function(pgam(i,j,k) + 4.0)) - 3.0 * std::log(lamc(i,j,k))) *
2401  (std::exp(m_aimm * (273.15 - t3d(i,j,k))) - 1.0);
2402 
2403  // PUT IN A CATCH HERE TO PREVENT DIVERGENCE BETWEEN NUMBER CONC. AND
2404  // MIXING RATIO, SINCE STRICT CONSERVATION NOT CHECKED FOR NUMBER CONC
2405  nnuccc = std::min(nnuccc, nc3d(i,j,k) / dt);
2406  }
2407 
2408  // AUTOCONVERSION OF CLOUD LIQUID WATER TO RAIN
2409  // FORMULA FROM BEHENG (1994)
2410  // USING NUMERICAL SIMULATION OF STOCHASTIC COLLECTION EQUATION
2411  // AND INITIAL CLOUD DROPLET SIZE DISTRIBUTION SPECIFIED
2412  // AS A GAMMA DISTRIBUTION
2413 
2414  // USE MINIMUM VALUE OF 1.E-6 TO PREVENT FLOATING POINT ERROR
2415  if (qc3d(i,j,k) >= 1.0e-6) {
2416  // hm add 12/13/06, replace with newer formula
2417  // from khairoutdinov and kogan 2000, mwr
2418  prc = 1350.0 * std::pow(qc3d(i,j,k), 2.47) *
2419  std::pow((nc3d(i,j,k) / 1.0e6 * rho(i,j,k)), -1.79);
2420 
2421  // note: nprc1 is change in nr,
2422  // nprc is change in nc
2423  nprc1 = prc / m_cons29;
2424  nprc = prc / (qc3d(i,j,k) / nc3d(i,j,k));
2425 
2426  // hm bug fix 3/20/12
2427  nprc = std::min(nprc, nc3d(i,j,k) / dt);
2428  nprc1 = std::min(nprc1, nprc);
2429  }
2430  // SNOW AGGREGATION FROM PASSARELLI, 1978, USED BY REISNER, 1998
2431  // THIS IS HARD-WIRED FOR BS = 0.4 FOR NOW
2432  if (qni3d(i,j,k) >= 1.0e-8) {
2433  nsagg = m_cons15 * asn(i,j,k) * std::pow(rho(i,j,k), ((2.0 + m_bs) / 3.0)) *
2434  std::pow(qni3d(i,j,k), ((2.0 + m_bs) / 3.0)) *
2435  std::pow((ns3d(i,j,k) * rho(i,j,k)), ((4.0 - m_bs) / 3.0)) / rho(i,j,k);
2436  }
2437 
2438  // ACCRETION OF CLOUD DROPLETS ONTO SNOW/GRAUPEL
2439  // HERE USE CONTINUOUS COLLECTION EQUATION WITH
2440  // SIMPLE GRAVITATIONAL COLLECTION KERNEL IGNORING
2441 
2442  // SNOW
2443  if (qni3d(i,j,k) >= 1.0e-8 && qc3d(i,j,k) >= m_qsmall) {
2444  psacws = m_cons13 * asn(i,j,k) * qc3d(i,j,k) * rho(i,j,k) *
2445  n0s(i,j,k) / std::pow(lams(i,j,k), (m_bs + 3.0));
2446 
2447  npsacws = m_cons13 * asn(i,j,k) * nc3d(i,j,k) * rho(i,j,k) *
2448  n0s(i,j,k) / std::pow(lams(i,j,k), (m_bs + 3.0));
2449  }
2450 
2451  // COLLECTION OF CLOUD WATER BY GRAUPEL
2452  if (qg3d(i,j,k) >= 1.0e-8 && qc3d(i,j,k) >= m_qsmall) {
2453  psacwg = m_cons14 * agn(i,j,k) * qc3d(i,j,k) * rho(i,j,k) *
2454  n0g(i,j,k) / std::pow(lamg(i,j,k), (m_bg + 3.0));
2455 
2456  npsacwg = m_cons14 * agn(i,j,k) * nc3d(i,j,k) * rho(i,j,k) *
2457  n0g(i,j,k) / std::pow(lamg(i,j,k), (m_bg + 3.0));
2458  }
2459  // hm, add 12/13/06
2460  // CLOUD ICE COLLECTING DROPLETS, ASSUME THAT CLOUD ICE MEAN DIAM > 100 MICRON
2461  // BEFORE RIMING CAN OCCUR
2462  // ASSUME THAT RIME COLLECTED ON CLOUD ICE DOES NOT LEAD
2463  // TO HALLET-MOSSOP SPLINTERING
2464  if (qi3d(i,j,k) >= 1.0e-8 && qc3d(i,j,k) >= m_qsmall) {
2465  // PUT IN SIZE DEPENDENT COLLECTION EFFICIENCY BASED ON STOKES LAW
2466  // FROM THOMPSON ET AL. 2004, MWR
2467  if (1.0 / lami(i,j,k) >= 100.0e-6) {
2468  psacwi = m_cons16 * ain(i,j,k) * qc3d(i,j,k) * rho(i,j,k) *
2469  n0i(i,j,k) / std::pow(lami(i,j,k), (m_bi + 3.0));
2470 
2471  npsacwi = m_cons16 * ain(i,j,k) * nc3d(i,j,k) * rho(i,j,k) *
2472  n0i(i,j,k) / std::pow(lami(i,j,k), (m_bi + 3.0));
2473  }
2474  }
2475 
2476  // ACCRETION OF RAIN WATER BY SNOW
2477  // FORMULA FROM IKAWA AND SAITO, 1991, USED BY REISNER ET AL, 1998
2478  if (qr3d(i,j,k) >= 1.0e-8 && qni3d(i,j,k) >= 1.0e-8) {
2479  amrex::Real ums_local = asn(i,j,k) * m_cons3 / std::pow(lams(i,j,k), m_bs);
2480  amrex::Real umr_local = arn(i,j,k) * m_cons4 / std::pow(lamr(i,j,k), m_br);
2481  amrex::Real uns_local = asn(i,j,k) * m_cons5 / std::pow(lams(i,j,k), m_bs);
2482  amrex::Real unr_local = arn(i,j,k) * m_cons6 / std::pow(lamr(i,j,k), m_br);
2483 
2484  // SET REASLISTIC LIMITS ON FALLSPEEDS
2485  // bug fix, 10/08/09
2486  dum = std::pow(m_rhosu / rho(i,j,k), 0.54);
2487  ums_local = std::min(ums_local, 1.2 * dum);
2488  uns_local = std::min(uns_local, 1.2 * dum);
2489  umr_local = std::min(umr_local, 9.1 * dum);
2490  unr_local = std::min(unr_local, 9.1 * dum);
2491 
2492  pracs = m_cons41 * (std::sqrt(std::pow(1.2 * umr_local - 0.95 * ums_local, 2) +
2493  0.08 * ums_local * umr_local) * rho(i,j,k) * n0r(i,j,k) * n0s(i,j,k) /
2494  std::pow(lamr(i,j,k), 3) * (5.0 / (std::pow(lamr(i,j,k), 3) * lams(i,j,k)) +
2495  2.0 / (std::pow(lamr(i,j,k), 2) * std::pow(lams(i,j,k), 2)) +
2496  0.5 / (lamr(i,j,k) * std::pow(lams(i,j,k), 3))));
2497 
2498  npracs = m_cons32 * rho(i,j,k) * std::sqrt(1.7 * std::pow(unr_local - uns_local, 2) +
2499  0.3 * unr_local * uns_local) * n0r(i,j,k) * n0s(i,j,k) *
2500  (1.0 / (std::pow(lamr(i,j,k), 3) * lams(i,j,k)) +
2501  1.0 / (std::pow(lamr(i,j,k), 2) * std::pow(lams(i,j,k), 2)) +
2502  1.0 / (lamr(i,j,k) * std::pow(lams(i,j,k), 3)));
2503 
2504  // MAKE SURE PRACS DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2505  // AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2506  // RIME-SPLINTERING
2507  pracs = std::min(pracs, qr3d(i,j,k) / dt);
2508 
2509  // COLLECTION OF SNOW BY RAIN - NEEDED FOR GRAUPEL CONVERSION CALCULATIONS
2510  // ONLY CALCULATE IF SNOW AND RAIN MIXING RATIOS EXCEED 0.1 G/KG
2511  // hm modify for wrfv3.1
2512  if (qni3d(i,j,k) >= 0.1e-3 && qr3d(i,j,k) >= 0.1e-3) {
2513  psacr = m_cons31 * (std::sqrt(std::pow(1.2 * umr_local - 0.95 * ums_local, 2) +
2514  0.08 * ums_local * umr_local) * rho(i,j,k) * n0r(i,j,k) * n0s(i,j,k) /
2515  std::pow(lams(i,j,k), 3) * (5.0 / (std::pow(lams(i,j,k), 3) * lamr(i,j,k)) +
2516  2.0 / (std::pow(lams(i,j,k), 2) * std::pow(lamr(i,j,k), 2)) +
2517  0.5 / (lams(i,j,k) * std::pow(lamr(i,j,k), 3))));
2518  }
2519  }
2520 
2521  // COLLECTION OF RAINWATER BY GRAUPEL, FROM IKAWA AND SAITO 1990,
2522  // USED BY REISNER ET AL 1998
2523  if (qr3d(i,j,k) >= 1.0e-8 && qg3d(i,j,k) >= 1.0e-8) {
2524  amrex::Real umg_local = agn(i,j,k) * m_cons7 / std::pow(lamg(i,j,k), m_bg);
2525  amrex::Real umr_local = arn(i,j,k) * m_cons4 / std::pow(lamr(i,j,k), m_br);
2526  amrex::Real ung_local = agn(i,j,k) * m_cons8 / std::pow(lamg(i,j,k), m_bg);
2527  amrex::Real unr_local = arn(i,j,k) * m_cons6 / std::pow(lamr(i,j,k), m_br);
2528 
2529  // SET REASLISTIC LIMITS ON FALLSPEEDS
2530  // bug fix, 10/08/09
2531  dum = std::pow(m_rhosu / rho(i,j,k), 0.54);
2532  umg_local = std::min(umg_local, 20.0 * dum);
2533  ung_local = std::min(ung_local, 20.0 * dum);
2534  umr_local = std::min(umr_local, 9.1 * dum);
2535  unr_local = std::min(unr_local, 9.1 * dum);
2536 
2537  pracg = m_cons41 * (std::sqrt(std::pow(1.2 * umr_local - 0.95 * umg_local, 2) +
2538  0.08 * umg_local * umr_local) * rho(i,j,k) * n0r(i,j,k) * n0g(i,j,k) /
2539  std::pow(lamr(i,j,k), 3) * (5.0 / (std::pow(lamr(i,j,k), 3) * lamg(i,j,k)) +
2540  2.0 / (std::pow(lamr(i,j,k), 2) * std::pow(lamg(i,j,k), 2)) +
2541  0.5 / (lamr(i,j,k) * std::pow(lamg(i,j,k), 3))));
2542 
2543  npracg = m_cons32 * rho(i,j,k) * std::sqrt(1.7 * std::pow(unr_local - ung_local, 2) +
2544  0.3 * unr_local * ung_local) * n0r(i,j,k) * n0g(i,j,k) *
2545  (1.0 / (std::pow(lamr(i,j,k), 3) * lamg(i,j,k)) +
2546  1.0 / (std::pow(lamr(i,j,k), 2) * std::pow(lamg(i,j,k), 2)) +
2547  1.0 / (lamr(i,j,k) * std::pow(lamg(i,j,k), 3)));
2548 
2549  // MAKE SURE PRACG DOESN'T EXCEED TOTAL RAIN MIXING RATIO
2550  // AS THIS MAY OTHERWISE RESULT IN TOO MUCH TRANSFER OF WATER DURING
2551  // RIME-SPLINTERING
2552  pracg = std::min(pracg, qr3d(i,j,k) / dt);
2553  }
2554 
2555  // RIME-SPLINTERING - SNOW
2556  // HALLET-MOSSOP (1974)
2557  // NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2558  // hm add threshold snow and droplet mixing ratio for rime-splintering
2559  // to limit rime-splintering in stratiform clouds
2560  // these thresholds correspond with graupel thresholds in rh 1984
2561  //v1.4
2562  if (qni3d(i,j,k) >= 0.1e-3) {
2563  if (qc3d(i,j,k) >= 0.5e-3 || qr3d(i,j,k) >= 0.1e-3) {
2564  if (psacws > 0.0 || pracs > 0.0) {
2565  if (t3d(i,j,k) < 270.16 && t3d(i,j,k) > 265.16) {
2566  amrex::Real fmult = 0.0;
2567 
2568  if (t3d(i,j,k) > 270.16) {
2569  fmult = 0.0;
2570  } else if (t3d(i,j,k) <= 270.16 && t3d(i,j,k) > 268.16) {
2571  fmult = (270.16 - t3d(i,j,k)) / 2.0;
2572  } else if (t3d(i,j,k) >= 265.16 && t3d(i,j,k) <= 268.16) {
2573  fmult = (t3d(i,j,k) - 265.16) / 3.0;
2574  } else if (t3d(i,j,k) < 265.16) {
2575  fmult = 0.0;
2576  }
2577 
2578  // 1000 IS TO CONVERT FROM KG TO G
2579  // SPLINTERING FROM DROPLETS ACCRETED ONTO SNOW
2580  if (psacws > 0.0) {
2581  nmults = 35.0e4 * psacws * fmult * 1000.0;
2582  qmults = nmults * m_mmult;
2583 
2584  // CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2585  // THAN WAS RIMED ONTO SNOW
2586  qmults = std::min(qmults, psacws);
2587  psacws = psacws - qmults;
2588  }
2589 
2590  // RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2591  if (pracs > 0.0) {
2592  nmultr = 35.0e4 * pracs * fmult * 1000.0;
2593  qmultr = nmultr * m_mmult;
2594 
2595  // CONSTRAIN SO THAT TRANSFER OF MASS FROM SNOW TO ICE CANNOT BE MORE MASS
2596  // THAN WAS RIMED ONTO SNOW
2597  qmultr = std::min(qmultr, pracs);
2598  pracs = pracs - qmultr;
2599  }
2600  }
2601  }
2602  }
2603  }
2604 
2605  // RIME-SPLINTERING - GRAUPEL
2606  // HALLET-MOSSOP (1974)
2607  // NUMBER OF SPLINTERS FORMED IS BASED ON MASS OF RIMED WATER
2608  // hm add threshold snow mixing ratio for rime-splintering
2609  // to limit rime-splintering in stratiform clouds
2610  // v1.4
2611  if (qg3d(i,j,k) >= 0.1e-3) {
2612  if (qc3d(i,j,k) >= 0.5e-3 || qr3d(i,j,k) >= 0.1e-3) {
2613  if (psacwg > 0.0 || pracg > 0.0) {
2614  if (t3d(i,j,k) < 270.16 && t3d(i,j,k) > 265.16) {
2615  amrex::Real fmult = 0.0;
2616 
2617  if (t3d(i,j,k) > 270.16) {
2618  fmult = 0.0;
2619  } else if (t3d(i,j,k) <= 270.16 && t3d(i,j,k) > 268.16) {
2620  fmult = (270.16 - t3d(i,j,k)) / 2.0;
2621  } else if (t3d(i,j,k) >= 265.16 && t3d(i,j,k) <= 268.16) {
2622  fmult = (t3d(i,j,k) - 265.16) / 3.0;
2623  } else if (t3d(i,j,k) < 265.16) {
2624  fmult = 0.0;
2625  }
2626 
2627  // 1000 IS TO CONVERT FROM KG TO G
2628  // SPLINTERING FROM DROPLETS ACCRETED ONTO GRAUPEL
2629  if (psacwg > 0.0) {
2630  nmultg = 35.0e4 * psacwg * fmult * 1000.0;
2631  qmultg = nmultg * m_mmult;
2632 
2633  // CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2634  // THAN WAS RIMED ONTO GRAUPEL
2635  qmultg = std::min(qmultg, psacwg);
2636  psacwg = psacwg - qmultg;
2637  }
2638 
2639  // RIMING AND SPLINTERING FROM ACCRETED RAINDROPS
2640  if (pracg > 0.0) {
2641  nmultrg = 35.0e4 * pracg * fmult * 1000.0;
2642  qmultrg = nmultrg * m_mmult;
2643 
2644  // CONSTRAIN SO THAT TRANSFER OF MASS FROM GRAUPEL TO ICE CANNOT BE MORE MASS
2645  // THAN WAS RIMED ONTO GRAUPEL
2646  qmultrg = std::min(qmultrg, pracg);
2647  pracg = pracg - qmultrg;
2648  }
2649  }
2650  }
2651  }
2652  }
2653 
2654  // CONVERSION OF RIMED CLOUD WATER ONTO SNOW TO GRAUPEL/HAIL
2655  if (psacws > 0.0) {
2656  // ONLY ALLOW CONVERSION IF QNI > 0.1 AND QC > 0.5 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2657  if (qni3d(i,j,k) >= 0.1e-3 && qc3d(i,j,k) >= 0.5e-3) {
2658  // PORTION OF RIMING CONVERTED TO GRAUPEL (REISNER ET AL. 1998, ORIGINALLY IS1991)
2659  pgsacw = std::min(psacws, m_cons17 * dt * n0s(i,j,k) * qc3d(i,j,k) * qc3d(i,j,k) *
2660  asn(i,j,k) * asn(i,j,k) /
2661  (rho(i,j,k) * std::pow(lams(i,j,k), (2.0 * m_bs + 2.0))));
2662 
2663  // MIX RAT CONVERTED INTO GRAUPEL AS EMBRYO (REISNER ET AL. 1998, ORIG M1990)
2664  dum = std::max(m_rhosn / (m_rhog - m_rhosn) * pgsacw, 0.0);
2665 
2666  // NUMBER CONCENTRAITON OF EMBRYO GRAUPEL FROM RIMING OF SNOW
2667  nscng = dum / m_mg0 * rho(i,j,k);
2668  // LIMIT MAX NUMBER CONVERTED TO SNOW NUMBER
2669  nscng = std::min(nscng, ns3d(i,j,k) / dt);
2670 
2671  // PORTION OF RIMING LEFT FOR SNOW
2672  psacws = psacws - pgsacw;
2673  }
2674  }
2675 
2676  // CONVERSION OF RIMED RAINWATER ONTO SNOW CONVERTED TO GRAUPEL
2677  if (pracs > 0.0) {
2678  // ONLY ALLOW CONVERSION IF QNI > 0.1 AND QR > 0.1 G/KG FOLLOWING RUTLEDGE AND HOBBS (1984)
2679  if (qni3d(i,j,k) >= 0.1e-3 && qr3d(i,j,k) >= 0.1e-3) {
2680  // PORTION OF COLLECTED RAINWATER CONVERTED TO GRAUPEL (REISNER ET AL. 1998)
2681  dum = m_cons18 * std::pow(4.0 / lams(i,j,k), 3) * std::pow(4.0 / lams(i,j,k), 3) /
2682  (m_cons18 * std::pow(4.0 / lams(i,j,k), 3) * std::pow(4.0 / lams(i,j,k), 3) +
2683  m_cons19 * std::pow(4.0 / lamr(i,j,k), 3) * std::pow(4.0 / lamr(i,j,k), 3));
2684  dum = std::min(dum, 1.0);
2685  dum = std::max(dum, 0.0);
2686 
2687  pgracs = (1.0 - dum) * pracs;
2688  ngracs = (1.0 - dum) * npracs;
2689 
2690  // LIMIT MAX NUMBER CONVERTED TO MIN OF EITHER RAIN OR SNOW NUMBER CONCENTRATION
2691  ngracs = std::min(ngracs, nr3d(i,j,k) / dt);
2692  ngracs = std::min(ngracs, ns3d(i,j,k) / dt);
2693 
2694  // AMOUNT LEFT FOR SNOW PRODUCTION
2695  pracs = pracs - pgracs;
2696  npracs = npracs - ngracs;
2697 
2698  // CONVERSION TO GRAUPEL DUE TO COLLECTION OF SNOW BY RAIN
2699  psacr = psacr * (1.0 - dum);
2700  }
2701  }
2702 
2703  // FREEZING OF RAIN DROPS
2704  // FREEZING ALLOWED BELOW -4 C
2705  if (t3d(i,j,k) < 269.15 && qr3d(i,j,k) >= m_qsmall) {
2706  // IMMERSION FREEZING (BIGG 1953)
2707  // hm fix 7/15/13 for consistency w/ original formula
2708  mnuccr = m_cons20 * nr3d(i,j,k) * (std::exp(m_aimm * (273.15 - t3d(i,j,k))) - 1.0) /
2709  std::pow(lamr(i,j,k), 3) / std::pow(lamr(i,j,k), 3);
2710 
2711  nnuccr = m_pi * nr3d(i,j,k) * m_bimm * (std::exp(m_aimm * (273.15 - t3d(i,j,k))) - 1.0) /
2712  std::pow(lamr(i,j,k), 3);
2713 
2714  // PREVENT DIVERGENCE BETWEEN MIXING RATIO AND NUMBER CONC
2715  nnuccr = std::min(nnuccr, nr3d(i,j,k) / dt);
2716  }
2717 
2718  // ACCRETION OF CLOUD LIQUID WATER BY RAIN
2719  // CONTINUOUS COLLECTION EQUATION WITH
2720  // GRAVITATIONAL COLLECTION KERNEL, DROPLET FALL SPEED NEGLECTED
2721  if (qr3d(i,j,k) >= 1.0e-8 && qc3d(i,j,k) >= 1.0e-8) {
2722  // 12/13/06 hm add, replace with newer formula from
2723  // khairoutdinov and kogan 2000, mwr
2724  dum = qc3d(i,j,k) * qr3d(i,j,k);
2725  pra = 67.0 * std::pow(dum, 1.15);
2726  npra = pra / (qc3d(i,j,k) / nc3d(i,j,k));
2727  }
2728 
2729  // SELF-COLLECTION OF RAIN DROPS
2730  // FROM BEHENG(1994)
2731  // FROM NUMERICAL SIMULATION OF THE STOCHASTIC COLLECTION EQUATION
2732  // AS DESCRINED ABOVE FOR AUTOCONVERSION
2733  if (qr3d(i,j,k) >= 1.0e-8) {
2734  // include breakup add 10/09/09
2735  dum1 = 300.0e-6;
2736  if (1.0 / lamr(i,j,k) < dum1) {
2737  dum = 1.0;
2738  } else if (1.0 / lamr(i,j,k) >= dum1) {
2739  dum = 2.0 - std::exp(2300.0 * (1.0 / lamr(i,j,k) - dum1));
2740  }
2741  nragg = -5.78 * dum * nr3d(i,j,k) * qr3d(i,j,k) * rho(i,j,k);
2742  }
2743 
2744  // AUTOCONVERSION OF CLOUD ICE TO SNOW
2745  // FOLLOWING HARRINGTON ET AL. (1995) WITH MODIFICATION
2746  // HERE IT IS ASSUMED THAT AUTOCONVERSION CAN ONLY OCCUR WHEN THE
2747  // ICE IS GROWING, I.E. IN CONDITIONS OF ICE SUPERSATURATION
2748  if (qi3d(i,j,k) >= 1.0e-8 && qvqvsi >= 1.0) {
2749  nprci = m_cons21 * (qv3d(i,j,k) - qvi) * rho(i,j,k) *
2750  n0i(i,j,k) * std::exp(-lami(i,j,k) * m_dcs) * dv / abi;
2751  prci = m_cons22 * nprci;
2752  nprci = std::min(nprci, ni3d(i,j,k) / dt);
2753  }
2754 
2755  // ACCRETION OF CLOUD ICE BY SNOW
2756  // FOR THIS CALCULATION, IT IS ASSUMED THAT THE VS >> VI
2757  // AND DS >> DI FOR CONTINUOUS COLLECTION
2758  if (qni3d(i,j,k) >= 1.0e-8 && qi3d(i,j,k) >= m_qsmall) {
2759  prai = m_cons23 * asn(i,j,k) * qi3d(i,j,k) * rho(i,j,k) * n0s(i,j,k) /
2760  std::pow(lams(i,j,k), (m_bs + 3.0));
2761  nprai = m_cons23 * asn(i,j,k) * ni3d(i,j,k) *
2762  rho(i,j,k) * n0s(i,j,k) /
2763  std::pow(lams(i,j,k), (m_bs + 3.0));
2764  nprai = std::min(nprai, ni3d(i,j,k) / dt);
2765  }
2766 
2767  // hm, add 12/13/06, collision of rain and ice to produce snow or graupel
2768  // follows reisner et al. 1998
2769  // assumed fallspeed and size of ice crystal << than for rain
2770  if (qr3d(i,j,k) >= 1.0e-8 && qi3d(i,j,k) >= 1.0e-8 && t3d(i,j,k) <= 273.15) {
2771  // allow graupel formation from rain-ice collisions only if rain mixing ratio > 0.1 g/kg,
2772  // otherwise add to snow
2773  if (qr3d(i,j,k) >= 0.1e-3) {
2774  niacr = m_cons24 * ni3d(i,j,k) * n0r(i,j,k)* arn(i,j,k) /
2775  std::pow(lamr(i,j,k), (m_br + 3.0)) * rho(i,j,k);
2776  piacr = m_cons25 * ni3d(i,j,k) * n0r(i,j,k) * arn(i,j,k) /
2777  std::pow(lamr(i,j,k), (m_br + 3.0)) / std::pow(lamr(i,j,k), 3) * rho(i,j,k);
2778  praci = m_cons24 * qi3d(i,j,k) * n0r(i,j,k) * arn(i,j,k) /
2779  std::pow(lamr(i,j,k), (m_br + 3.0)) * rho(i,j,k);
2780  niacr = std::min(niacr, nr3d(i,j,k) / dt);
2781  niacr = std::min(niacr, ni3d(i,j,k) / dt);
2782  } else {
2783  niacrs = m_cons24 * ni3d(i,j,k) * n0r(i,j,k) * arn(i,j,k) /
2784  std::pow(lamr(i,j,k), (m_br + 3.0)) * rho(i,j,k);
2785  piacrs = m_cons25 * ni3d(i,j,k) * n0r(i,j,k) * arn(i,j,k) /
2786  std::pow(lamr(i,j,k), (m_br + 3.0)) / std::pow(lamr(i,j,k), 3) * rho(i,j,k);
2787  pracis = m_cons24 * qi3d(i,j,k) * n0r(i,j,k) * arn(i,j,k) /
2788  std::pow(lamr(i,j,k), (m_br + 3.0)) * rho(i,j,k);
2789  niacrs = std::min(niacrs, nr3d(i,j,k) / dt);
2790  niacrs = std::min(niacrs, ni3d(i,j,k) / dt);
2791  }
2792  }
2793 
2794  // NUCLEATION OF CLOUD ICE FROM HOMOGENEOUS AND HETEROGENEOUS FREEZING ON AEROSOL
2795  if (m_inuc == 0) {
2796  // ADD THRESHOLD ACCORDING TO GREG THOMSPON
2797  if ((qvqvs >= 0.999 && t3d(i,j,k) <= 265.15) || qvqvsi >= 1.08) {
2798  // hm, modify dec. 5, 2006, replace with cooper curve
2799  kc2 = 0.005 * std::exp(0.304 * (273.15 - t3d(i,j,k))) * 1000.0; // CONVERT FROM L-1 TO M-3
2800  // LIMIT TO 500 L-1
2801  kc2 = std::min(kc2, 500.0e3);
2802  kc2 = std::max(kc2 / rho(i,j,k), 0.0); // CONVERT TO KG-1
2803 
2804  if (kc2 > ni3d(i,j,k) + ns3d(i,j,k) + ng3d(i,j,k)) {
2805  nnuccd = (kc2 - ni3d(i,j,k) - ns3d(i,j,k) - ng3d(i,j,k)) / dt;
2806  mnuccd = nnuccd * m_mi0;
2807  }
2808  }
2809  } else if (m_inuc == 1) {
2810  if (t3d(i,j,k) < 273.15 && qvqvsi > 1.0) {
2811  kc2 = 0.16 * 1000.0 / rho(i,j,k); // CONVERT FROM L-1 TO KG-1
2812  if (kc2 > ni3d(i,j,k) + ns3d(i,j,k) + ng3d(i,j,k)) {
2813  nnuccd = (kc2 - ni3d(i,j,k) - ns3d(i,j,k) - ng3d(i,j,k)) / dt;
2814  mnuccd = nnuccd * m_mi0;
2815  }
2816  }
2817  }
2818 
2819  // CALCULATE EVAP/SUB/DEP TERMS FOR QI,QNI,QR
2820  // NO VENTILATION FOR CLOUD ICE
2821  epsi = 0.0;
2822  if (qi3d(i,j,k) >= m_qsmall) {
2823  epsi = 2.0 * m_pi * n0i(i,j,k) * rho(i,j,k) * dv / (lami(i,j,k) * lami(i,j,k));
2824  }
2825 
2826  // VENTILATION FOR SNOW
2827  epss = 0.0;
2828  if (qni3d(i,j,k) >= m_qsmall) {
2829  epss = 2.0 * m_pi * n0s(i,j,k) * rho(i,j,k) * dv *
2830  (m_f1s / (lams(i,j,k) * lams(i,j,k)) +
2831  m_f2s * std::pow(asn(i,j,k) * rho(i,j,k) / mu(i,j,k), 0.5) *
2832  std::pow(sc_schmidt, (1.0 / 3.0)) * m_cons10 /
2833  std::pow(lams(i,j,k), m_cons35));
2834  }
2835 
2836  // Ventilation for graupel
2837  epsg = 0.0;
2838  if (qg3d(i,j,k) >= m_qsmall) {
2839  epsg = 2.0 * m_pi * n0g(i,j,k) * rho(i,j,k) * dv *
2840  (m_f1s / (lamg(i,j,k) * lamg(i,j,k)) +
2841  m_f2s * std::pow(agn(i,j,k) * rho(i,j,k) / mu(i,j,k), 0.5) *
2842  std::pow(sc_schmidt, (1.0 / 3.0)) * m_cons11 /
2843  std::pow(lamg(i,j,k), m_cons36));
2844  }
2845 
2846  // VENTILATION FOR RAIN
2847  epsr = 0.0;
2848  if (qr3d(i,j,k) >= m_qsmall) {
2849  epsr = 2.0 * m_pi * n0r(i,j,k) * rho(i,j,k) * dv *
2850  (m_f1r / (lamr(i,j,k) * lamr(i,j,k)) +
2851  m_f2r * std::pow(arn(i,j,k) * rho(i,j,k) / mu(i,j,k), 0.5) *
2852  std::pow(sc_schmidt, (1.0 / 3.0)) * m_cons9 /
2853  std::pow(lamr(i,j,k), m_cons34));
2854  }
2855 
2856  // ONLY INCLUDE REGION OF ICE SIZE DIST < DCS
2857  // DUM IS FRACTION OF D*N(D) < DCS
2858  // LOGIC BELOW FOLLOWS THAT OF HARRINGTON ET AL. 1995 (JAS)
2859  if (qi3d(i,j,k) >= m_qsmall) {
2860  dum = (1.0 - std::exp(-lami(i,j,k) * m_dcs) * (1.0 + lami(i,j,k) * m_dcs));
2861  prd = epsi * (qv3d(i,j,k) - qvi) / abi * dum;
2862  } else {
2863  dum = 0.0;
2864  }
2865 
2866  // ADD DEPOSITION IN TAIL OF ICE SIZE DIST TO SNOW IF SNOW IS PRESENT
2867  if (qni3d(i,j,k) >= m_qsmall) {
2868  prds = epss * (qv3d(i,j,k) - qvi) / abi +
2869  epsi * (qv3d(i,j,k) - qvi) / abi * (1.0 - dum);
2870  } else {
2871  // OTHERWISE ADD TO CLOUD ICE
2872  prd = prd + epsi * (qv3d(i,j,k) - qvi) / abi * (1.0 - dum);
2873  }
2874 
2875  // VAPOR DPEOSITION ON GRAUPEL
2876  prdg = epsg * (qv3d(i,j,k) - qvi) / abi;
2877 
2878  // NO CONDENSATION ONTO RAIN, ONLY EVAP
2879  if (qv3d(i,j,k) < qvs) {
2880  pre = epsr * (qv3d(i,j,k) - qvs) / ab;
2881  pre = std::min(pre, 0.0);
2882  } else {
2883  pre = 0.0;
2884  }
2885 
2886  // MAKE SURE NOT PUSHED INTO ICE SUPERSAT/SUBSAT
2887  // FORMULA FROM REISNER 2 SCHEME
2888  dum = (qv3d(i,j,k) - qvi) / dt;
2889 
2890  fudgef = 0.9999;
2891  sum_dep = prd + prds + mnuccd + prdg;
2892 
2893  if ((dum > 0.0 && sum_dep > dum * fudgef) ||
2894  (dum < 0.0 && sum_dep < dum * fudgef)) {
2895  mnuccd = fudgef * mnuccd * dum / sum_dep;
2896  prd = fudgef * prd * dum / sum_dep;
2897  prds = fudgef * prds * dum / sum_dep;
2898  prdg = fudgef * prdg * dum / sum_dep;
2899  }
2900 
2901  // IF CLOUD ICE/SNOW/GRAUPEL VAP DEPOSITION IS NEG, THEN ASSIGN TO SUBLIMATION PROCESSES
2902  if (prd < 0.0) {
2903  eprd = prd;
2904  prd = 0.0;
2905  }
2906  if (prds < 0.0) {
2907  eprds = prds;
2908  prds = 0.0;
2909  }
2910  if (prdg < 0.0) {
2911  eprdg = prdg;
2912  prdg = 0.0;
2913  }
2914  // CONSERVATION OF WATER
2915  // THIS IS ADOPTED LOOSELY FROM MM5 RESINER CODE. HOWEVER, HERE WE
2916  // ONLY ADJUST PROCESSES THAT ARE NEGATIVE, RATHER THAN ALL PROCESSES.
2917 
2918  // IF MIXING RATIOS LESS THAN QSMALL, THEN NO DEPLETION OF WATER
2919  // THROUGH MICROPHYSICAL PROCESSES, SKIP CONSERVATION
2920 
2921  // NOTE: CONSERVATION CHECK NOT APPLIED TO NUMBER CONCENTRATION SPECIES. ADDITIONAL CATCH
2922  // BELOW WILL PREVENT NEGATIVE NUMBER CONCENTRATION
2923  // FOR EACH MICROPHYSICAL PROCESS WHICH PROVIDES A SOURCE FOR NUMBER, THERE IS A CHECK
2924  // TO MAKE SURE THAT CAN'T EXCEED TOTAL NUMBER OF DEPLETED SPECIES WITH THE TIME
2925  // STEP
2926 
2927  // ****SENSITIVITY - NO ICE
2928  if (m_iliq == 1) {
2929  mnuccc = 0.0;
2930  nnuccc = 0.0;
2931  mnuccr = 0.0;
2932  nnuccr = 0.0;
2933  mnuccd = 0.0;
2934  nnuccd = 0.0;
2935  }
2936 
2937  // ****SENSITIVITY - NO GRAUPEL
2938  if (m_igraup == 1) {
2939  pracg = 0.0;
2940  psacr = 0.0;
2941  psacwg = 0.0;
2942  prdg = 0.0;
2943  eprdg = 0.0;
2944  evpmg = 0.0;
2945  pgmlt = 0.0;
2946  npracg = 0.0;
2947  npsacwg = 0.0;
2948  nscng = 0.0;
2949  ngracs = 0.0;
2950  nsubg = 0.0;
2951  ngmltg = 0.0;
2952  ngmltr = 0.0;
2953 
2954  // fix 053011
2955  piacrs = piacrs + piacr;
2956  piacr = 0.0;
2957 
2958  // fix 070713
2959  pracis = pracis + praci;
2960  praci = 0.0;
2961  psacws = psacws + pgsacw;
2962  pgsacw = 0.0;
2963  pracs = pracs + pgracs;
2964  pgracs = 0.0;
2965  }
2966 
2967  // CONSERVATION OF QC
2968  dum = (prc + pra + mnuccc + psacws + psacwi + qmults + psacwg + pgsacw + qmultg) * dt;
2969 
2970  if (dum > qc3d(i,j,k) && qc3d(i,j,k) >= m_qsmall) {
2971  ratio = qc3d(i,j,k) / dum;
2972 
2973  prc = prc * ratio;
2974  pra = pra * ratio;
2975  mnuccc = mnuccc * ratio;
2976  psacws = psacws * ratio;
2977  psacwi = psacwi * ratio;
2978  qmults = qmults * ratio;
2979  qmultg = qmultg * ratio;
2980  psacwg = psacwg * ratio;
2981  pgsacw = pgsacw * ratio;
2982  }
2983 
2984  // CONSERVATION OF QI
2985  dum = (-prd - mnuccc + prci + prai - qmults - qmultg - qmultr - qmultrg
2986  - mnuccd + praci + pracis - eprd - psacwi) * dt;
2987 
2988  if (dum > qi3d(i,j,k) && qi3d(i,j,k) >= m_qsmall) {
2989  ratio = (qi3d(i,j,k) / dt + prd + mnuccc + qmults + qmultg + qmultr + qmultrg +
2990  mnuccd + psacwi) /
2991  (prci + prai + praci + pracis - eprd);
2992 
2993  prci = prci * ratio;
2994  prai = prai * ratio;
2995  praci = praci * ratio;
2996  pracis = pracis * ratio;
2997  eprd = eprd * ratio;
2998  }
2999 
3000  // CONSERVATION OF QR
3001  dum = ((pracs - pre) + (qmultr + qmultrg - prc) + (mnuccr - pra) +
3002  piacr + piacrs + pgracs + pracg) * dt;
3003 
3004  if (dum > qr3d(i,j,k) && qr3d(i,j,k) >= m_qsmall) {
3005  ratio = (qr3d(i,j,k) / dt + prc + pra) /
3006  (-pre + qmultr + qmultrg + pracs + mnuccr + piacr + piacrs + pgracs + pracg);
3007 
3008  pre = pre * ratio;
3009  pracs = pracs * ratio;
3010  qmultr = qmultr * ratio;
3011  qmultrg = qmultrg * ratio;
3012  mnuccr = mnuccr * ratio;
3013  piacr = piacr * ratio;
3014  piacrs = piacrs * ratio;
3015  pgracs = pgracs * ratio;
3016  pracg = pracg * ratio;
3017  }
3018 
3019  // CONSERVATION OF QNI
3020  if (m_igraup == 0) {
3021  dum = (-prds - psacws - prai - prci - pracs - eprds + psacr - piacrs - pracis) * dt;
3022 
3023  if (dum > qni3d(i,j,k) && qni3d(i,j,k) >= m_qsmall) {
3024  ratio = (qni3d(i,j,k) / dt + prds + psacws + prai + prci + pracs + piacrs + pracis) /
3025  (-eprds + psacr);
3026 
3027  eprds = eprds * ratio;
3028  psacr = psacr * ratio;
3029  }
3030  } else if (m_igraup == 1) {
3031  // FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3032  dum = (-prds - psacws - prai - prci - pracs - eprds + psacr - piacrs - pracis - mnuccr) * dt;
3033 
3034  if (dum > qni3d(i,j,k) && qni3d(i,j,k) >= m_qsmall) {
3035  ratio = (qni3d(i,j,k) / dt + prds + psacws + prai + prci + pracs + piacrs + pracis + mnuccr) /
3036  (-eprds + psacr);
3037 
3038  eprds = eprds * ratio;
3039  psacr = psacr * ratio;
3040  }
3041  }
3042 
3043  // CONSERVATION OF QG
3044  dum = (-psacwg - pracg - pgsacw - pgracs - prdg - mnuccr - eprdg - piacr - praci - psacr) * dt;
3045 
3046  if (dum > qg3d(i,j,k) && qg3d(i,j,k) >= m_qsmall) {
3047  ratio = (qg3d(i,j,k) / dt + psacwg + pracg + pgsacw + pgracs + prdg + mnuccr + psacr +
3048  piacr + praci) / (-eprdg);
3049 
3050  eprdg = eprdg * ratio;
3051  }
3052 
3053  // TENDENCIES
3054  qv3dten(i,j,k) = qv3dten(i,j,k) + (-pre - prd - prds - mnuccd - eprd - eprds - prdg - eprdg);
3055 
3056  // bug fix hm, 3/1/11, include piacr and piacrs
3057  t3dten(i,j,k) = t3dten(i,j,k) + (pre * xxlv(i,j,k) +
3058  (prd + prds + mnuccd + eprd + eprds + prdg + eprdg) * xxls(i,j,k) +
3059  (psacws + psacwi + mnuccc + mnuccr + qmults + qmultg + qmultr + qmultrg + pracs +
3060  psacwg + pracg + pgsacw + pgracs + piacr + piacrs) * xlf(i,j,k)) / cpm(i,j,k);
3061 
3062  qc3dten(i,j,k) = qc3dten(i,j,k) +
3063  (-pra - prc - mnuccc + pcc -
3064  psacws - psacwi - qmults - qmultg - psacwg - pgsacw);
3065 
3066  qi3dten(i,j,k) = qi3dten(i,j,k) +
3067  (prd + eprd + psacwi + mnuccc - prci -
3068  prai + qmults + qmultg + qmultr + qmultrg + mnuccd - praci - pracis);
3069 
3070  qr3dten(i,j,k) = qr3dten(i,j,k) +
3071  (pre + pra + prc - pracs - mnuccr - qmultr - qmultrg -
3072  piacr - piacrs - pracg - pgracs);
3073  if (m_igraup == 0) {
3074  qni3dten(i,j,k) = qni3dten(i,j,k) +
3075  (prai + psacws + prds + pracs + prci + eprds - psacr + piacrs + pracis);
3076 
3077  ns3dten(i,j,k) = ns3dten(i,j,k) + (nsagg + nprci - nscng - ngracs + niacrs);
3078 
3079  qg3dten(i,j,k) = qg3dten(i,j,k) + (pracg + psacwg + pgsacw + pgracs +
3080  prdg + eprdg + mnuccr + piacr + praci + psacr);
3081 
3082  ng3dten(i,j,k) = ng3dten(i,j,k) + (nscng + ngracs + nnuccr + niacr);
3083  } else if (m_igraup == 1) {
3084  // FOR NO GRAUPEL, NEED TO INCLUDE FREEZING OF RAIN FOR SNOW
3085  qni3dten(i,j,k) = qni3dten(i,j,k) +
3086  (prai + psacws + prds + pracs + prci + eprds - psacr + piacrs + pracis + mnuccr);
3087 
3088  ns3dten(i,j,k) = ns3dten(i,j,k) + (nsagg + nprci - nscng - ngracs + niacrs + nnuccr);
3089  }
3090 
3091  nc3dten(i,j,k) = nc3dten(i,j,k) + (-nnuccc - npsacws -
3092  npra - nprc - npsacwi - npsacwg);
3093 
3094  ni3dten(i,j,k) = ni3dten(i,j,k) +
3095  (nnuccc - nprci - nprai + nmults + nmultg + nmultr + nmultrg +
3096  nnuccd - niacr - niacrs);
3097 
3098  nr3dten(i,j,k) = nr3dten(i,j,k) + (nprc1 - npracs - nnuccr +
3099  nragg - niacr - niacrs - npracg - ngracs);
3100 
3101  // hm add, wrf-chem, add tendencies for c2prec
3102  c2prec = pra + prc + psacws + qmults + qmultg + psacwg +
3103  pgsacw + mnuccc + psacwi;
3104 
3105  // CALCULATE SATURATION ADJUSTMENT TO CONDENSE EXTRA VAPOR ABOVE
3106  // WATER SATURATION
3107  dumt = t3d(i,j,k) + dt * t3dten(i,j,k);
3108  dumqv = qv3d(i,j,k) + dt * qv3dten(i,j,k);
3109 
3110  // hm, add fix for low pressure, 5/12/10
3111  dum = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(dumt, 0));
3112  dumqss = m_ep_2 * dum / (pres(i,j,k) - dum);
3113 
3114  dumqc = qc3d(i,j,k) + dt * qc3dten(i,j,k);
3115  dumqc = std::max(dumqc, 0.0);
3116 
3117  // SATURATION ADJUSTMENT FOR LIQUID
3118  dums = dumqv - dumqss;
3119 
3120  pcc = dums / (1.0 + std::pow(xxlv(i,j,k), 2) * dumqss / (cpm(i,j,k) * m_Rv * std::pow(dumt, 2))) / dt;
3121 
3122  if (pcc * dt + dumqc < 0.0) {
3123  pcc = -dumqc / dt;
3124  }
3125 
3126  qv3dten(i,j,k) = qv3dten(i,j,k) - pcc;
3127  t3dten(i,j,k) = t3dten(i,j,k) + pcc * xxlv(i,j,k) / cpm(i,j,k);
3128  qc3dten(i,j,k) = qc3dten(i,j,k) + pcc;
3129  // SUBLIMATE, MELT, OR EVAPORATE NUMBER CONCENTRATION
3130  // THIS FORMULATION ASSUMES 1:1 RATIO BETWEEN MASS LOSS AND
3131  // LOSS OF NUMBER CONCENTRATION
3132  if (eprd < 0.0) {
3133  dum = eprd * dt / qi3d(i,j,k);
3134  dum = std::max(-1.0, dum);
3135  nsubi = dum * ni3d(i,j,k) / dt;
3136  }
3137 
3138  if (eprds < 0.0) {
3139  dum = eprds * dt / qni3d(i,j,k);
3140  dum = std::max(-1.0, dum);
3141  nsubs = dum * ns3d(i,j,k) / dt;
3142  }
3143 
3144  if (pre < 0.0) {
3145  dum = pre * dt / qr3d(i,j,k);
3146  dum = std::max(-1.0, dum);
3147  nsubr = dum * nr3d(i,j,k) / dt;
3148  }
3149 
3150  if (eprdg < 0.0) {
3151  dum = eprdg * dt / qg3d(i,j,k);
3152  dum = std::max(-1.0, dum);
3153  nsubg = dum * ng3d(i,j,k) / dt;
3154  }
3155 
3156  // UPDATE TENDENCIES
3157  ni3dten(i,j,k) = ni3dten(i,j,k) + nsubi;
3158  ns3dten(i,j,k) = ns3dten(i,j,k) + nsubs;
3159  ng3dten(i,j,k) = ng3dten(i,j,k) + nsubg;
3160  nr3dten(i,j,k) = nr3dten(i,j,k) + nsubr;
3161  }
3162  ltrue = 1;
3163  }
3164  // label_200:
3165 
3166  }
3167  for(int k=klo; k<=khi; k++) {
3168  // INITIALIZE PRECIP AND SNOW RATES
3169  precrt(i,j,k) = 0.0;
3170  snowrt(i,j,k) = 0.0;
3171  // hm added 7/13/13
3172  snowprt(i,j,k) = 0.0;
3173  grplprt(i,j,k) = 0.0;
3174  }
3175  nstep = 1;
3176  if(ltrue != 0) {
3177  //goto 400
3178  // CALCULATE SEDIMENTATION
3179  // THE NUMERICS HERE FOLLOW FROM REISNER ET AL. (1998)
3180  // FALLOUT TERMS ARE CALCULATED ON SPLIT TIME STEPS TO ENSURE NUMERICAL
3181  // STABILITY, I.E. COURANT# < 1
3182  // Loop from top to bottom (KTE to KTS)
3183  for(int k=khi; k>=klo; k--) {
3184 
3185  amrex::Real dum; // DUM: General dummy variable
3186 
3187  amrex::Real di0; // DI0: Characteristic diameter for ice
3188  amrex::Real ds0; // DS0: Characteristic diameter for snow
3189  amrex::Real dg0; // DG0: Characteristic diameter for graupel
3190  amrex::Real lammax; // LAMMAX: Maximum value for slope parameter
3191  amrex::Real lammin; // LAMMIN: Minimum value for slope parameter
3192 
3193  ds0 = 3.0; // Size distribution parameter for snow
3194  di0 = 3.0; // Size distribution parameter for cloud ice
3195  dg0 = 3.0; // Size distribution parameter for graupel
3196 
3197  // Update prognostic variables with tendencies
3198  dumi(i,j,k) = qi3d(i,j,k) + qi3dten(i,j,k) * dt;
3199  dumqs(i,j,k) = qni3d(i,j,k) + qni3dten(i,j,k) * dt;
3200  dumr(i,j,k) = qr3d(i,j,k) + qr3dten(i,j,k) * dt;
3201  dumfni(i,j,k) = ni3d(i,j,k) + ni3dten(i,j,k) * dt;
3202  dumfns(i,j,k) = ns3d(i,j,k) + ns3dten(i,j,k) * dt;
3203  dumfnr(i,j,k) = nr3d(i,j,k) + nr3dten(i,j,k) * dt;
3204  dumc(i,j,k) = qc3d(i,j,k) + qc3dten(i,j,k) * dt;
3205  dumfnc(i,j,k) = nc3d(i,j,k) + nc3dten(i,j,k) * dt;
3206  dumg(i,j,k) = qg3d(i,j,k) + qg3dten(i,j,k) * dt;
3207  dumfng(i,j,k) = ng3d(i,j,k) + ng3dten(i,j,k) * dt;
3208 
3209  // SWITCH FOR CONSTANT DROPLET NUMBER
3210  if (iinum == 1) {
3211  dumfnc(i,j,k) = nc3d(i,j,k);
3212  }
3213 
3214  // MAKE SURE NUMBER CONCENTRATIONS ARE POSITIVE
3215  dumfni(i,j,k) = amrex::max(0., dumfni(i,j,k));
3216  dumfns(i,j,k) = amrex::max(0., dumfns(i,j,k));
3217  dumfnc(i,j,k) = amrex::max(0., dumfnc(i,j,k));
3218  dumfnr(i,j,k) = amrex::max(0., dumfnr(i,j,k));
3219  dumfng(i,j,k) = amrex::max(0., dumfng(i,j,k));
3220 
3221  // CLOUD ICE
3222  if (dumi(i,j,k) >= m_qsmall) {
3223  dlami(i,j,k) = std::pow(m_cons12 * dumfni(i,j,k) / dumi(i,j,k), 1.0/di0);
3224  dlami(i,j,k) = amrex::max(dlami(i,j,k), m_lammini);
3225  dlami(i,j,k) = amrex::min(dlami(i,j,k), m_lammaxi);
3226  }
3227 
3228  // RAIN
3229  if (dumr(i,j,k) >= m_qsmall) {
3230  dlamr(i,j,k) = std::pow(m_pi * m_rhow * dumfnr(i,j,k) / dumr(i,j,k), 1.0/3.0);
3231  dlamr(i,j,k) = amrex::max(dlamr(i,j,k), m_lamminr);
3232  dlamr(i,j,k) = amrex::min(dlamr(i,j,k), m_lammaxr);
3233  }
3234 
3235  // CLOUD DROPLETS
3236  if (dumc(i,j,k) >= m_qsmall) {
3237  dum = pres(i,j,k) / (287.15 * t3d(i,j,k));
3238  pgam(i,j,k) = 0.0005714 * (nc3d(i,j,k) / 1.0e6 * dum) + 0.2714;
3239  pgam(i,j,k) = 1.0 / (pgam(i,j,k) * pgam(i,j,k)) - 1.0;
3240  pgam(i,j,k) = amrex::max(pgam(i,j,k), 2.0);
3241  pgam(i,j,k) = amrex::min(pgam(i,j,k), 10.0);
3242 
3243  dlamc(i,j,k) = std::pow(m_cons26 * dumfnc(i,j,k) * gamma_function(pgam(i,j,k) + 4.0) /
3244  (dumc(i,j,k) * gamma_function(pgam(i,j,k) + 1.0)), 1.0/3.0);
3245  lammin = (pgam(i,j,k) + 1.0) / 60.0e-6;
3246  lammax = (pgam(i,j,k) + 1.0) / 1.0e-6;
3247  dlamc(i,j,k) = amrex::max(dlamc(i,j,k), lammin);
3248  dlamc(i,j,k) = amrex::min(dlamc(i,j,k), lammax);
3249  }
3250 
3251  // SNOW
3252  if (dumqs(i,j,k) >= m_qsmall) {
3253  dlams(i,j,k) = std::pow(m_cons1 * dumfns(i,j,k) / dumqs(i,j,k), 1.0/ds0);
3254  dlams(i,j,k) = amrex::max(dlams(i,j,k), m_lammins);
3255  dlams(i,j,k) = amrex::min(dlams(i,j,k), m_lammaxs);
3256  }
3257 
3258  // GRAUPEL
3259  if (dumg(i,j,k) >= m_qsmall) {
3260  dlamg(i,j,k) = std::pow(m_cons2 * dumfng(i,j,k) / dumg(i,j,k), 1.0/dg0);
3261  dlamg(i,j,k) = amrex::max(dlamg(i,j,k), m_lamming);
3262  dlamg(i,j,k) = amrex::min(dlamg(i,j,k), m_lammaxg);
3263  }
3264 
3265  // Calculate number-weighted and mass-weighted terminal fall speeds
3266  // CLOUD WATER
3267  if (dumc(i,j,k) >= m_qsmall) {
3268  unc(i,j,k) = acn(i,j,k) * gamma_function(1. + m_bc + pgam(i,j,k)) /
3269  (std::pow(dlamc(i,j,k), m_bc) * gamma_function(pgam(i,j,k) + 1.));
3270  umc(i,j,k) = acn(i,j,k) * gamma_function(4. + m_bc + pgam(i,j,k)) /
3271  (std::pow(dlamc(i,j,k), m_bc) * gamma_function(pgam(i,j,k) + 4.));
3272  } else {
3273  umc(i,j,k) = 0.;
3274  unc(i,j,k) = 0.;
3275  }
3276 
3277  // CLOUD ICE
3278  if (dumi(i,j,k) >= m_qsmall) {
3279  uni(i,j,k) = ain(i,j,k) * m_cons27 / std::pow(dlami(i,j,k), m_bi);
3280  umi(i,j,k) = ain(i,j,k) * m_cons28 / std::pow(dlami(i,j,k), m_bi);
3281  } else {
3282  umi(i,j,k) = 0.;
3283  uni(i,j,k) = 0.;
3284  }
3285 
3286  // RAIN
3287  if (dumr(i,j,k) >= m_qsmall) {
3288  unr(i,j,k) = arn(i,j,k) * m_cons6 / std::pow(dlamr(i,j,k), m_br);
3289  umr(i,j,k) = arn(i,j,k) * m_cons4 / std::pow(dlamr(i,j,k), m_br);
3290  } else {
3291  umr(i,j,k) = 0.;
3292  unr(i,j,k) = 0.;
3293  }
3294 
3295  // SNOW
3296  if (dumqs(i,j,k) >= m_qsmall) {
3297  ums(i,j,k) = asn(i,j,k) * m_cons3 / std::pow(dlams(i,j,k), m_bs);
3298  uns(i,j,k) = asn(i,j,k) * m_cons5 / std::pow(dlams(i,j,k), m_bs);
3299  } else {
3300  ums(i,j,k) = 0.;
3301  uns(i,j,k) = 0.;
3302  }
3303 
3304  // GRAUPEL
3305  if (dumg(i,j,k) >= m_qsmall) {
3306  umg(i,j,k) = agn(i,j,k) * m_cons7 / std::pow(dlamg(i,j,k), m_bg);
3307  ung(i,j,k) = agn(i,j,k) * m_cons8 / std::pow(dlamg(i,j,k), m_bg);
3308  } else {
3309  umg(i,j,k) = 0.;
3310  ung(i,j,k) = 0.;
3311  }
3312 
3313  // SET REALISTIC LIMITS ON FALLSPEED
3314  // Bug fix, 10/08/09
3315  dum = std::pow(m_rhosu / rho(i,j,k), 0.54);
3316  ums(i,j,k) = std::min(ums(i,j,k), 1.2 * dum);
3317  uns(i,j,k) = std::min(uns(i,j,k), 1.2 * dum);
3318 
3319  // Fix 053011
3320  // Fix for correction by AA 4/6/11
3321  umi(i,j,k) = std::min(umi(i,j,k), 1.2 * std::pow(m_rhosu / rho(i,j,k), 0.35));
3322  uni(i,j,k) = std::min(uni(i,j,k), 1.2 * std::pow(m_rhosu / rho(i,j,k), 0.35));
3323  umr(i,j,k) = std::min(umr(i,j,k), 9.1 * dum);
3324  unr(i,j,k) = std::min(unr(i,j,k), 9.1 * dum);
3325  umg(i,j,k) = std::min(umg(i,j,k), 20. * dum);
3326  ung(i,j,k) = std::min(ung(i,j,k), 20. * dum);
3327 
3328  // Set fall speed values
3329  fr(i,j,k) = umr(i,j,k); // RAIN FALL SPEED
3330  fi(i,j,k) = umi(i,j,k); // CLOUD ICE FALL SPEED
3331  fni(i,j,k) = uni(i,j,k); // CLOUD ICE NUMBER FALL SPEED
3332  fs(i,j,k) = ums(i,j,k); // SNOW FALL SPEED
3333  fns(i,j,k) = uns(i,j,k); // SNOW NUMBER FALL SPEED
3334  fnr(i,j,k) = unr(i,j,k); // RAIN NUMBER FALL SPEED
3335  fc(i,j,k) = umc(i,j,k); // CLOUD WATER FALL SPEED
3336  fnc(i,j,k) = unc(i,j,k); // CLOUD NUMBER FALL SPEED
3337  fg(i,j,k) = umg(i,j,k); // GRAUPEL FALL SPEED
3338  fng(i,j,k) = ung(i,j,k); // GRAUPEL NUMBER FALL SPEED
3339 
3340  // V3.3 MODIFY FALLSPEED BELOW LEVEL OF PRECIP
3341  if (fr(i,j,k) < 1.e-10) {
3342  fr(i,j,k) = fr(i,j,k+1);
3343  }
3344  if (fi(i,j,k) < 1.e-10) {
3345  fi(i,j,k) = fi(i,j,k+1);
3346  }
3347  if (fni(i,j,k) < 1.e-10) {
3348  fni(i,j,k) = fni(i,j,k+1);
3349  }
3350  if (fs(i,j,k) < 1.e-10) {
3351  fs(i,j,k) = fs(i,j,k+1);
3352  }
3353  if (fns(i,j,k) < 1.e-10) {
3354  fns(i,j,k) = fns(i,j,k+1);
3355  }
3356  if (fnr(i,j,k) < 1.e-10) {
3357  fnr(i,j,k) = fnr(i,j,k+1);
3358  }
3359  if (fc(i,j,k) < 1.e-10) {
3360  fc(i,j,k) = fc(i,j,k+1);
3361  }
3362  if (fnc(i,j,k) < 1.e-10) {
3363  fnc(i,j,k) = fnc(i,j,k+1);
3364  }
3365  if (fg(i,j,k) < 1.e-10) {
3366  fg(i,j,k) = fg(i,j,k+1);
3367  }
3368  if (fng(i,j,k) < 1.e-10) {
3369  fng(i,j,k) = fng(i,j,k+1);
3370  }
3371 
3372  // CALCULATE NUMBER OF SPLIT TIME STEPS
3373  // Find maximum fall speed at this point
3374  rgvm(i,j,k) = std::max({fr(i,j,k), fi(i,j,k), fs(i,j,k), fc(i,j,k),
3375  fni(i,j,k), fnr(i,j,k), fns(i,j,k), fnc(i,j,k),
3376  fg(i,j,k), fng(i,j,k)});
3377 
3378  // Calculate number of steps (dt and nstep would need to be defined elsewhere)
3379  nstep = std::max(static_cast<int>(rgvm(i,j,k) * dt / dzq(i,j,k) + 1.), nstep);
3380  // MULTIPLY VARIABLES BY RHO
3381  dumr(i,j,k) = dumr(i,j,k) * rho(i,j,k); // Rain water content * density
3382  dumi(i,j,k) = dumi(i,j,k) * rho(i,j,k); // Cloud ice content * density
3383  dumfni(i,j,k) = dumfni(i,j,k) * rho(i,j,k); // Cloud ice number * density
3384  dumqs(i,j,k) = dumqs(i,j,k) * rho(i,j,k); // Snow content * density
3385  dumfns(i,j,k) = dumfns(i,j,k) * rho(i,j,k); // Snow number * density
3386  dumfnr(i,j,k) = dumfnr(i,j,k) * rho(i,j,k); // Rain number * density
3387  dumc(i,j,k) = dumc(i,j,k) * rho(i,j,k); // Cloud water content * density
3388  dumfnc(i,j,k) = dumfnc(i,j,k) * rho(i,j,k); // Cloud droplet number * density
3389  dumg(i,j,k) = dumg(i,j,k) * rho(i,j,k); // Graupel content * density
3390  dumfng(i,j,k) = dumfng(i,j,k) * rho(i,j,k); // Graupel number * density
3391  }
3392  // Main time stepping loop for sedimentation
3393  for (int n = 1; n <= nstep; n++) {
3394  // Calculate initial fallout for each hydrometeor type for all levels
3395  for (int k = klo; k <= khi; k++) {
3396  faloutr(i,j,k) = fr(i,j,k) * dumr(i,j,k);
3397  falouti(i,j,k) = fi(i,j,k) * dumi(i,j,k);
3398  faloutni(i,j,k) = fni(i,j,k) * dumfni(i,j,k);
3399  falouts(i,j,k) = fs(i,j,k) * dumqs(i,j,k);
3400  faloutns(i,j,k) = fns(i,j,k) * dumfns(i,j,k);
3401  faloutnr(i,j,k) = fnr(i,j,k) * dumfnr(i,j,k);
3402  faloutc(i,j,k) = fc(i,j,k) * dumc(i,j,k);
3403  faloutnc(i,j,k) = fnc(i,j,k) * dumfnc(i,j,k);
3404  faloutg(i,j,k) = fg(i,j,k) * dumg(i,j,k);
3405  faloutng(i,j,k) = fng(i,j,k) * dumfng(i,j,k);
3406  }
3407 
3408  // Process top of model level
3409  int k = khi;
3410 
3411  // Calculate tendencies at top level
3412  faltndr(i,j,k) = faloutr(i,j,k) / dzq(i,j,k);
3413  faltndi(i,j,k) = falouti(i,j,k) / dzq(i,j,k);
3414  faltndni(i,j,k) = faloutni(i,j,k) / dzq(i,j,k);
3415  faltnds(i,j,k) = falouts(i,j,k) / dzq(i,j,k);
3416  faltndns(i,j,k) = faloutns(i,j,k) / dzq(i,j,k);
3417  faltndnr(i,j,k) = faloutnr(i,j,k) / dzq(i,j,k);
3418  faltndc(i,j,k) = faloutc(i,j,k) / dzq(i,j,k);
3419  faltndnc(i,j,k) = faloutnc(i,j,k) / dzq(i,j,k);
3420  faltndg(i,j,k) = faloutg(i,j,k) / dzq(i,j,k);
3421  faltndng(i,j,k) = faloutng(i,j,k) / dzq(i,j,k);
3422 
3423  // Add fallout terms to Eulerian tendencies (scaled by time step and density)
3424  qrsten(i,j,k) = qrsten(i,j,k) - faltndr(i,j,k) / nstep / rho(i,j,k);
3425  qisten(i,j,k) = qisten(i,j,k) - faltndi(i,j,k) / nstep / rho(i,j,k);
3426  ni3dten(i,j,k) = ni3dten(i,j,k) - faltndni(i,j,k) / nstep / rho(i,j,k);
3427  qnisten(i,j,k) = qnisten(i,j,k) - faltnds(i,j,k) / nstep / rho(i,j,k);
3428  ns3dten(i,j,k) = ns3dten(i,j,k) - faltndns(i,j,k) / nstep / rho(i,j,k);
3429  nr3dten(i,j,k) = nr3dten(i,j,k) - faltndnr(i,j,k) / nstep / rho(i,j,k);
3430  qcsten(i,j,k) = qcsten(i,j,k) - faltndc(i,j,k) / nstep / rho(i,j,k);
3431  nc3dten(i,j,k) = nc3dten(i,j,k) - faltndnc(i,j,k) / nstep / rho(i,j,k);
3432  qgsten(i,j,k) = qgsten(i,j,k) - faltndg(i,j,k) / nstep / rho(i,j,k);
3433  ng3dten(i,j,k) = ng3dten(i,j,k) - faltndng(i,j,k) / nstep / rho(i,j,k);
3434 
3435  // Update temporary working variables
3436  dumr(i,j,k) = dumr(i,j,k) - faltndr(i,j,k) * dt / nstep;
3437  dumi(i,j,k) = dumi(i,j,k) - faltndi(i,j,k) * dt / nstep;
3438  dumfni(i,j,k) = dumfni(i,j,k) - faltndni(i,j,k) * dt / nstep;
3439  dumqs(i,j,k) = dumqs(i,j,k) - faltnds(i,j,k) * dt / nstep;
3440  dumfns(i,j,k) = dumfns(i,j,k) - faltndns(i,j,k) * dt / nstep;
3441  dumfnr(i,j,k) = dumfnr(i,j,k) - faltndnr(i,j,k) * dt / nstep;
3442  dumc(i,j,k) = dumc(i,j,k) - faltndc(i,j,k) * dt / nstep;
3443  dumfnc(i,j,k) = dumfnc(i,j,k) - faltndnc(i,j,k) * dt / nstep;
3444  dumg(i,j,k) = dumg(i,j,k) - faltndg(i,j,k) * dt / nstep;
3445  dumfng(i,j,k) = dumfng(i,j,k) - faltndng(i,j,k) * dt / nstep;
3446 
3447  // Process remaining levels from top to bottom
3448  for (k = khi-1; k >= klo; k--) {
3449  // Calculate tendencies based on difference between levels
3450  faltndr(i,j,k) = (faloutr(i,j,k+1) - faloutr(i,j,k)) / dzq(i,j,k);
3451  faltndi(i,j,k) = (falouti(i,j,k+1) - falouti(i,j,k)) / dzq(i,j,k);
3452  faltndni(i,j,k) = (faloutni(i,j,k+1) - faloutni(i,j,k)) / dzq(i,j,k);
3453  faltnds(i,j,k) = (falouts(i,j,k+1) - falouts(i,j,k)) / dzq(i,j,k);
3454  faltndns(i,j,k) = (faloutns(i,j,k+1) - faloutns(i,j,k)) / dzq(i,j,k);
3455  faltndnr(i,j,k) = (faloutnr(i,j,k+1) - faloutnr(i,j,k)) / dzq(i,j,k);
3456  faltndc(i,j,k) = (faloutc(i,j,k+1) - faloutc(i,j,k)) / dzq(i,j,k);
3457  faltndnc(i,j,k) = (faloutnc(i,j,k+1) - faloutnc(i,j,k)) / dzq(i,j,k);
3458  faltndg(i,j,k) = (faloutg(i,j,k+1) - faloutg(i,j,k)) / dzq(i,j,k);
3459  faltndng(i,j,k) = (faloutng(i,j,k+1) - faloutng(i,j,k)) / dzq(i,j,k);
3460 
3461  // Add fallout terms to Eulerian tendencies (positive here, as mass flows in from above)
3462  qrsten(i,j,k) = qrsten(i,j,k) + faltndr(i,j,k) / nstep / rho(i,j,k);
3463  qisten(i,j,k) = qisten(i,j,k) + faltndi(i,j,k) / nstep / rho(i,j,k);
3464  ni3dten(i,j,k) = ni3dten(i,j,k) + faltndni(i,j,k) / nstep / rho(i,j,k);
3465  qnisten(i,j,k) = qnisten(i,j,k) + faltnds(i,j,k) / nstep / rho(i,j,k);
3466  ns3dten(i,j,k) = ns3dten(i,j,k) + faltndns(i,j,k) / nstep / rho(i,j,k);
3467  nr3dten(i,j,k) = nr3dten(i,j,k) + faltndnr(i,j,k) / nstep / rho(i,j,k);
3468  qcsten(i,j,k) = qcsten(i,j,k) + faltndc(i,j,k) / nstep / rho(i,j,k);
3469  nc3dten(i,j,k) = nc3dten(i,j,k) + faltndnc(i,j,k) / nstep / rho(i,j,k);
3470  qgsten(i,j,k) = qgsten(i,j,k) + faltndg(i,j,k) / nstep / rho(i,j,k);
3471  ng3dten(i,j,k) = ng3dten(i,j,k) + faltndng(i,j,k) / nstep / rho(i,j,k);
3472  // Update temporary working variables
3473  dumr(i,j,k) = dumr(i,j,k) + faltndr(i,j,k) * dt / nstep;
3474  dumi(i,j,k) = dumi(i,j,k) + faltndi(i,j,k) * dt / nstep;
3475  dumfni(i,j,k) = dumfni(i,j,k) + faltndni(i,j,k) * dt / nstep;
3476  dumqs(i,j,k) = dumqs(i,j,k) + faltnds(i,j,k) * dt / nstep;
3477  dumfns(i,j,k) = dumfns(i,j,k) + faltndns(i,j,k) * dt / nstep;
3478  dumfnr(i,j,k) = dumfnr(i,j,k) + faltndnr(i,j,k) * dt / nstep;
3479  dumc(i,j,k) = dumc(i,j,k) + faltndc(i,j,k) * dt / nstep;
3480  dumfnc(i,j,k) = dumfnc(i,j,k) + faltndnc(i,j,k) * dt / nstep;
3481  dumg(i,j,k) = dumg(i,j,k) + faltndg(i,j,k) * dt / nstep;
3482  dumfng(i,j,k) = dumfng(i,j,k) + faltndng(i,j,k) * dt / nstep;
3483  }
3484  // Get precipitation and snowfall accumulation during the time step
3485  // Factor of 1000 converts from m to mm, but division by density
3486  // of liquid water cancels this factor of 1000
3487  int kts=klo;
3488  precrt(i,j,klo) += (faloutr(i,j,kts) + faloutc(i,j,kts) + falouts(i,j,kts) +
3489  falouti(i,j,kts) + faloutg(i,j,kts)) * dt / nstep;
3490  snowrt(i,j,klo) += (falouts(i,j,kts) + falouti(i,j,kts) + faloutg(i,j,kts)) * dt / nstep;
3491 
3492  // Added 7/13/13
3493  snowprt(i,j,klo) += (falouti(i,j,kts) + falouts(i,j,kts)) * dt / nstep;
3494  grplprt(i,j,klo) += faloutg(i,j,kts) * dt / nstep;
3495  }
3496  for(int k=klo; k<=khi; k++) {
3497  amrex::Real evs; // EVS: Saturation vapor pressure
3498  amrex::Real eis; // EIS: Ice saturation vapor pressure
3499  amrex::Real qvs; // QVS: Saturation mixing ratio
3500  amrex::Real qvi; // QVI: Ice saturation mixing ratio
3501  amrex::Real qvqvs; // QVQVS: Saturation ratio
3502  amrex::Real qvqvsi; // QVQVSI: Ice saturation ratio
3503 
3504  // ADD ON SEDIMENTATION TENDENCIES FOR MIXING RATIO TO REST OF TENDENCIES
3505  qr3dten(i,j,k) = qr3dten(i,j,k) + qrsten(i,j,k);
3506  qi3dten(i,j,k) = qi3dten(i,j,k) + qisten(i,j,k);
3507  qc3dten(i,j,k) = qc3dten(i,j,k) + qcsten(i,j,k);
3508  qg3dten(i,j,k) = qg3dten(i,j,k) + qgsten(i,j,k);
3509  qni3dten(i,j,k) = qni3dten(i,j,k) + qnisten(i,j,k);
3510  // PUT ALL CLOUD ICE IN SNOW CATEGORY IF MEAN DIAMETER EXCEEDS 2 * dcs
3511  // bug fix
3512  if (qi3d(i,j,k) >= m_qsmall && t3d(i,j,k) < 273.15 && lami(i,j,k) >= 1.e-10) {
3513  if (1.0/lami(i,j,k) >= 2.0*m_dcs) {
3514  qni3dten(i,j,k) = qni3dten(i,j,k) + qi3d(i,j,k)/dt + qi3dten(i,j,k);
3515  ns3dten(i,j,k) = ns3dten(i,j,k) + ni3d(i,j,k)/dt + ni3dten(i,j,k);
3516  qi3dten(i,j,k) = -qi3d(i,j,k)/dt;
3517  ni3dten(i,j,k) = -ni3d(i,j,k)/dt;
3518  }
3519  }
3520 
3521  // Add tendencies to ensure consistency between mixing ratio and number concentration
3522  qc3d(i,j,k) = qc3d(i,j,k) + qc3dten(i,j,k)*dt;
3523  qi3d(i,j,k) = qi3d(i,j,k) + qi3dten(i,j,k)*dt;
3524  qni3d(i,j,k) = qni3d(i,j,k) + qni3dten(i,j,k)*dt;
3525  qr3d(i,j,k) = qr3d(i,j,k) + qr3dten(i,j,k)*dt;
3526  nc3d(i,j,k) = nc3d(i,j,k) + nc3dten(i,j,k)*dt;
3527  ni3d(i,j,k) = ni3d(i,j,k) + ni3dten(i,j,k)*dt;
3528  ns3d(i,j,k) = ns3d(i,j,k) + ns3dten(i,j,k)*dt;
3529  nr3d(i,j,k) = nr3d(i,j,k) + nr3dten(i,j,k)*dt;
3530  if (m_igraup == 0) {
3531  qg3d(i,j,k) = qg3d(i,j,k) + qg3dten(i,j,k)*dt;
3532  ng3d(i,j,k) = ng3d(i,j,k) + ng3dten(i,j,k)*dt;
3533  }
3534 
3535  // ADD TEMPERATURE AND WATER VAPOR TENDENCIES FROM MICROPHYSICS
3536  t3d(i,j,k) = t3d(i,j,k) + t3dten(i,j,k)*dt;
3537  qv3d(i,j,k) = qv3d(i,j,k) + qv3dten(i,j,k)*dt;
3538  // SATURATION VAPOR PRESSURE AND MIXING RATIO
3539  // hm, add fix for low pressure, 5/12/10
3540  // Assuming POLYSVP is defined elsewhere
3541  evs = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(t3d(i,j,k), 0)); // PA
3542  eis = std::min(0.99 * pres(i,j,k), calc_saturation_vapor_pressure(t3d(i,j,k), 1)); // PA
3543 
3544  // MAKE SURE ICE SATURATION DOESN'T EXCEED WATER SAT. NEAR FREEZING
3545  if (eis > evs) {
3546  eis = evs; // temporary update: adjust ice saturation pressure
3547  }
3548 
3549  // SATURATION MIXING RATIOS
3550  qvs = m_ep_2 * evs / (pres(i,j,k) - evs); // budget equation: calculate water saturation mixing ratio
3551  qvi = m_ep_2 * eis / (pres(i,j,k) - eis); // budget equation: calculate ice saturation mixing ratio
3552 
3553  // SATURATION RATIOS
3554  qvqvs = qv3d(i,j,k) / qvs; // budget equation: calculate water saturation ratio
3555  qvqvsi = qv3d(i,j,k) / qvi; // budget equation: calculate ice saturation ratio
3556  // AT SUBSATURATION, REMOVE SMALL AMOUNTS OF CLOUD/PRECIP WATER
3557  if (qvqvs < 0.9) {
3558  if (qr3d(i,j,k) < 1.0e-8) {
3559  qv3d(i,j,k) += qr3d(i,j,k);
3560  t3d(i,j,k) -= qr3d(i,j,k) * xxlv(i,j,k) / cpm(i,j,k);
3561  qr3d(i,j,k) = 0.0;
3562  }
3563  if (qc3d(i,j,k) < 1.0e-8) {
3564  qv3d(i,j,k) += qc3d(i,j,k);
3565  t3d(i,j,k) -= qc3d(i,j,k) * xxlv(i,j,k) / cpm(i,j,k);
3566  qc3d(i,j,k) = 0.0;
3567  }
3568  }
3569  if (qvqvsi < 0.9) {
3570  if (qi3d(i,j,k) < 1.0e-8) {
3571  qv3d(i,j,k) += qi3d(i,j,k);
3572  t3d(i,j,k) -= qi3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k);
3573  qi3d(i,j,k) = 0.0;
3574  }
3575  if (qni3d(i,j,k) < 1.0e-8) {
3576  qv3d(i,j,k) += qni3d(i,j,k);
3577  t3d(i,j,k) -= qni3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k);
3578  qni3d(i,j,k) = 0.0;
3579  }
3580  if (qg3d(i,j,k) < 1.0e-8) {
3581  qv3d(i,j,k) += qg3d(i,j,k);
3582  t3d(i,j,k) -= qg3d(i,j,k) * xxls(i,j,k) / cpm(i,j,k);
3583  qg3d(i,j,k) = 0.0;
3584  }
3585  }
3586  // IF MIXING RATIO < QSMALL SET MIXING RATIO AND NUMBER CONC TO ZERO
3587  if (qc3d(i,j,k) < m_qsmall) {
3588  qc3d(i,j,k) = 0.0;
3589  nc3d(i,j,k) = 0.0;
3590  effc(i,j,k) = 0.0;
3591  }
3592  if (qr3d(i,j,k) < m_qsmall) {
3593  qr3d(i,j,k) = 0.0;
3594  nr3d(i,j,k) = 0.0;
3595  effr(i,j,k) = 0.0;
3596  }
3597  if (qi3d(i,j,k) < m_qsmall) {
3598  qi3d(i,j,k) = 0.0;
3599  ni3d(i,j,k) = 0.0;
3600  effi(i,j,k) = 0.0;
3601  }
3602  if (qni3d(i,j,k) < m_qsmall) {
3603  qni3d(i,j,k) = 0.0;
3604  ns3d(i,j,k) = 0.0;
3605  effs(i,j,k) = 0.0;
3606  }
3607  if (qg3d(i,j,k) < m_qsmall) {
3608  qg3d(i,j,k) = 0.0;
3609  ng3d(i,j,k) = 0.0;
3610  effg(i,j,k) = 0.0;
3611  }
3612  /*
3613  // Skip calculations if there is no cloud/precipitation water
3614  if ((qc3d(i,j,k) < m_qsmall && // CLOUD WATER MIXING RATIO (KG/KG)
3615  qi3d(i,j,k) < m_qsmall && // CLOUD ICE MIXING RATIO (KG/KG)
3616  qni3d(i,j,k) < m_qsmall && // SNOW MIXING RATIO (KG/KG)
3617  qr3d(i,j,k) < m_qsmall && // RAIN MIXING RATIO (KG/KG)
3618  qg3d(i,j,k) < m_qsmall)) { // GRAUPEL MIX RATIO (KG/KG)
3619  goto label_500;
3620  } else {*/
3621  if (!(qc3d(i,j,k) < m_qsmall && // CLOUD WATER MIXING RATIO (KG/KG)
3622  qi3d(i,j,k) < m_qsmall && // CLOUD ICE MIXING RATIO (KG/KG)
3623  qni3d(i,j,k) < m_qsmall && // SNOW MIXING RATIO (KG/KG)
3624  qr3d(i,j,k) < m_qsmall && // RAIN MIXING RATIO (KG/KG)
3625  qg3d(i,j,k) < m_qsmall)) { // GRAUPEL MIX RATIO (KG/KG)
3626  // CALCULATE INSTANTANEOUS PROCESSES
3627 
3628  // ADD MELTING OF CLOUD ICE TO FORM RAIN
3629  if (qi3d(i,j,k) >= m_qsmall && t3d(i,j,k) >= 273.15) {
3630  qr3d(i,j,k) = qr3d(i,j,k) + qi3d(i,j,k);
3631  t3d(i,j,k) = t3d(i,j,k) - qi3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k);
3632  qi3d(i,j,k) = 0.0;
3633  nr3d(i,j,k) = nr3d(i,j,k) + ni3d(i,j,k);
3634  ni3d(i,j,k) = 0.0;
3635  }
3636  // ****SENSITIVITY - NO ICE
3637  if ((m_iliq != 1)) {
3638 
3639  // HOMOGENEOUS FREEZING OF CLOUD WATER
3640  if (t3d(i,j,k) <= 233.15 && qc3d(i,j,k) >= m_qsmall) {
3641  qi3d(i,j,k) = qi3d(i,j,k) + qc3d(i,j,k);
3642  t3d(i,j,k) = t3d(i,j,k) + qc3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k);
3643  qc3d(i,j,k) = 0.0;
3644  ni3d(i,j,k) = ni3d(i,j,k) + nc3d(i,j,k);
3645  nc3d(i,j,k) = 0.0;
3646  }
3647  // HOMOGENEOUS FREEZING OF RAIN
3648  if (m_igraup == 0) {
3649  if (t3d(i,j,k) <= 233.15 && qr3d(i,j,k) >= m_qsmall) {
3650  qg3d(i,j,k) = qg3d(i,j,k) + qr3d(i,j,k);
3651  t3d(i,j,k) = t3d(i,j,k) + qr3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k);
3652  qr3d(i,j,k) = 0.0;
3653  ng3d(i,j,k) = ng3d(i,j,k) + nr3d(i,j,k);
3654  nr3d(i,j,k) = 0.0;
3655  }
3656  } else if (m_igraup == 1) {
3657  if (t3d(i,j,k) <= 233.15 && qr3d(i,j,k) >= m_qsmall) {
3658  qni3d(i,j,k) = qni3d(i,j,k) + qr3d(i,j,k);
3659  t3d(i,j,k) = t3d(i,j,k) + qr3d(i,j,k) * xlf(i,j,k) / cpm(i,j,k);
3660  qr3d(i,j,k) = 0.0;
3661  ns3d(i,j,k) = ns3d(i,j,k) + nr3d(i,j,k);
3662  nr3d(i,j,k) = 0.0;
3663  }
3664  }
3665 
3666  }/* else {
3667  Real dontdoanything=m_iliq;//printf("m_iliq: %d\n",m_iliq);// goto label_778;
3668  }*/
3669 
3670 // label_778:
3671  // MAKE SURE NUMBER CONCENTRATIONS AREN'T NEGATIVE
3672  ni3d(i,j,k) = std::max(0.0, ni3d(i,j,k));
3673  ns3d(i,j,k) = std::max(0.0, ns3d(i,j,k));
3674  nc3d(i,j,k) = std::max(0.0, nc3d(i,j,k));
3675  nr3d(i,j,k) = std::max(0.0, nr3d(i,j,k));
3676  ng3d(i,j,k) = std::max(0.0, ng3d(i,j,k));
3677 
3678  // CLOUD ICE
3679  if (qi3d(i,j,k) >= m_qsmall) {
3680  lami(i,j,k) = std::pow(m_cons12 * ni3d(i,j,k) / qi3d(i,j,k), 1.0/m_di);
3681  // CHECK FOR SLOPE
3682  // ADJUST VARS
3683  if (lami(i,j,k) < m_lammini) {
3684  lami(i,j,k) = m_lammini;
3685  n0i(i,j,k) = std::pow(lami(i,j,k), 4) * qi3d(i,j,k) / m_cons12;
3686  ni3d(i,j,k) = n0i(i,j,k) / lami(i,j,k);
3687  } else if (lami(i,j,k) > m_lammaxi) {
3688  lami(i,j,k) = m_lammaxi;
3689  n0i(i,j,k) = std::pow(lami(i,j,k), 4) * qi3d(i,j,k) / m_cons12;
3690  ni3d(i,j,k) = n0i(i,j,k) / lami(i,j,k);
3691  }
3692  }
3693 
3694  // RAIN
3695  if (qr3d(i,j,k) >= m_qsmall) {
3696  lamr(i,j,k) = std::pow(m_pi * m_rhow * nr3d(i,j,k) / qr3d(i,j,k), 1.0/3.0);
3697 
3698  // CHECK FOR SLOPE
3699  // ADJUST VARS
3700  if (lamr(i,j,k) < m_lamminr) {
3701  lamr(i,j,k) = m_lamminr;
3702  n0r(i,j,k) = std::pow(lamr(i,j,k), 4) * qr3d(i,j,k) / (m_pi * m_rhow);
3703  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k);
3704  } else if (lamr(i,j,k) > m_lammaxr) {
3705  lamr(i,j,k) = m_lammaxr;
3706  n0r(i,j,k) = std::pow(lamr(i,j,k), 4) * qr3d(i,j,k) / (m_pi * m_rhow);
3707  nr3d(i,j,k) = n0r(i,j,k) / lamr(i,j,k);
3708  }
3709  }
3710 
3711  // CLOUD DROPLETS
3712  // MARTIN ET AL. (1994) FORMULA FOR PGAM
3713  if (qc3d(i,j,k) >= m_qsmall) {
3714  amrex::Real dum = pres(i,j,k) / (287.15 * t3d(i,j,k));
3715  pgam(i,j,k) = 0.0005714 * (nc3d(i,j,k) / 1.0e6 * dum) + 0.2714;
3716  pgam(i,j,k) = 1.0/(std::pow(pgam(i,j,k), 2)) - 1.0;
3717  pgam(i,j,k) = std::max(pgam(i,j,k), 2.0);
3718  pgam(i,j,k) = std::min(pgam(i,j,k), 10.0);
3719 
3720  // CALCULATE LAMC
3721  lamc(i,j,k) = std::pow(m_cons26 * nc3d(i,j,k) * gamma_function(pgam(i,j,k) + 4.0) /
3722  (qc3d(i,j,k) * gamma_function(pgam(i,j,k) + 1.0)), 1.0/3.0);
3723 
3724  // LAMMIN, 60 MICRON DIAMETER
3725  // LAMMAX, 1 MICRON
3726  amrex::Real lammin = (pgam(i,j,k) + 1.0) / 60.0e-6;
3727  amrex::Real lammax = (pgam(i,j,k) + 1.0) / 1.0e-6;
3728 
3729  if (lamc(i,j,k) < lammin) {
3730  lamc(i,j,k) = lammin;
3731  nc3d(i,j,k) = std::exp(3.0 * std::log(lamc(i,j,k)) + std::log(qc3d(i,j,k)) +
3732  std::log(gamma_function(pgam(i,j,k) + 1.0)) - std::log(gamma_function(pgam(i,j,k) + 4.0))) / m_cons26;
3733  } else if (lamc(i,j,k) > lammax) {
3734  lamc(i,j,k) = lammax;
3735  nc3d(i,j,k) = std::exp(3.0 * std::log(lamc(i,j,k)) + std::log(qc3d(i,j,k)) +
3736  std::log(gamma_function(pgam(i,j,k) + 1.0)) - std::log(gamma_function(pgam(i,j,k) + 4.0))) / m_cons26;
3737  }
3738  }
3739 
3740  // SNOW
3741  if (qni3d(i,j,k) >= m_qsmall) {
3742  lams(i,j,k) = std::pow(m_cons1 * ns3d(i,j,k) / qni3d(i,j,k), 1.0/m_ds);
3743 
3744  // CHECK FOR SLOPE
3745  // ADJUST VARS
3746  if (lams(i,j,k) < m_lammins) {
3747  lams(i,j,k) = m_lammins;
3748  n0s(i,j,k) = std::pow(lams(i,j,k), 4) * qni3d(i,j,k) / m_cons1;
3749  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k);
3750  } else if (lams(i,j,k) > m_lammaxs) {
3751  lams(i,j,k) = m_lammaxs;
3752  n0s(i,j,k) = std::pow(lams(i,j,k), 4) * qni3d(i,j,k) / m_cons1;
3753  ns3d(i,j,k) = n0s(i,j,k) / lams(i,j,k);
3754  }
3755  }
3756 
3757  // GRAUPEL
3758  if (qg3d(i,j,k) >= m_qsmall) {
3759  lamg(i,j,k) = std::pow(m_cons2 * ng3d(i,j,k) / qg3d(i,j,k), 1.0/m_dg);
3760 
3761  // CHECK FOR SLOPE
3762  // ADJUST VARS
3763  if (lamg(i,j,k) < m_lamming) {
3764  lamg(i,j,k) = m_lamming;
3765  n0g(i,j,k) = std::pow(lamg(i,j,k), 4) * qg3d(i,j,k) / m_cons2;
3766  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k);
3767  } else if (lamg(i,j,k) > m_lammaxg) {
3768  lamg(i,j,k) = m_lammaxg;
3769  n0g(i,j,k) = std::pow(lamg(i,j,k), 4) * qg3d(i,j,k) / m_cons2;
3770  ng3d(i,j,k) = n0g(i,j,k) / lamg(i,j,k);
3771  }
3772  }
3773  }
3774 
3775 // label_500:
3776  // CALCULATE EFFECTIVE RADIUS
3777  if (qi3d(i,j,k) >= m_qsmall) {
3778  effi(i,j,k) = 3.0 / lami(i,j,k) / 2.0 * 1.0e6;
3779  } else {
3780  effi(i,j,k) = 25.0;
3781  }
3782 
3783  if (qni3d(i,j,k) >= m_qsmall) {
3784  effs(i,j,k) = 3.0 / lams(i,j,k) / 2.0 * 1.0e6;
3785  } else {
3786  effs(i,j,k) = 25.0;
3787  }
3788 
3789  if (qr3d(i,j,k) >= m_qsmall) {
3790  effr(i,j,k) = 3.0 / lamr(i,j,k) / 2.0 * 1.0e6;
3791  } else {
3792  effr(i,j,k) = 25.0;
3793  }
3794 
3795  if (qc3d(i,j,k) >= m_qsmall) {
3796  effc(i,j,k) = gamma_function(pgam(i,j,k) + 4.0) / gamma_function(pgam(i,j,k) + 3.0) / lamc(i,j,k) / 2.0 * 1.0e6;
3797  } else {
3798  effc(i,j,k) = 25.0;
3799  }
3800 
3801  if (qg3d(i,j,k) >= m_qsmall) {
3802  effg(i,j,k) = 3.0 / lamg(i,j,k) / 2.0 * 1.0e6;
3803  } else {
3804  effg(i,j,k) = 25.0;
3805  }
3806 
3807  // HM ADD 1/10/06, ADD UPPER BOUND ON ICE NUMBER, THIS IS NEEDED
3808  // TO PREVENT VERY LARGE ICE NUMBER DUE TO HOMOGENEOUS FREEZING
3809  // OF DROPLETS, ESPECIALLY WHEN INUM = 1, SET MAX AT 10 CM-3
3810  // HM, 12/28/12, LOWER MAXIMUM ICE CONCENTRATION TO ADDRESS PROBLEM
3811  // OF EXCESSIVE AND PERSISTENT ANVIL
3812  // NOTE: THIS MAY CHANGE/REDUCE SENSITIVITY TO AEROSOL/CCN CONCENTRATION
3813  ni3d(i,j,k) = std::min(ni3d(i,j,k), 0.3e6 / rho(i,j,k));
3814 
3815  // ADD BOUND ON DROPLET NUMBER - CANNOT EXCEED AEROSOL CONCENTRATION
3816  if (iinum == 0 && m_iact == 2) {
3817  nc3d(i,j,k) = std::min(nc3d(i,j,k), (m_nanew1 + m_nanew2) / rho(i,j,k));
3818  }
3819 
3820  // SWITCH FOR CONSTANT DROPLET NUMBER
3821  if (iinum == 1) {
3822  // CHANGE NDCNST FROM CM-3 TO KG-1
3823  nc3d(i,j,k) = m_ndcnst * 1.0e6 / rho(i,j,k);
3824  }
3825  }
3826 
3827  }/* else {
3828  goto label_400;
3829  }
3830  label_400:*/
3831  //End of _micro
3832  if(use_morr_cpp_answer) {
3833  for(int k=klo; k<=khi; k++) {
3834 
3835  // Transfer 1D variables back to 3D arrays
3836  qcl_arr(i,j,k) = qc3d(i,j,k);
3837  qci_arr(i,j,k) = qi3d(i,j,k);
3838  qps_arr(i,j,k) = qni3d(i,j,k);
3839  qpr_arr(i,j,k) = qr3d(i,j,k);
3840  ni_arr(i,j,k) = ni3d(i,j,k);
3841  ns_arr(i,j,k) = ns3d(i,j,k);
3842  nr_arr(i,j,k) = nr3d(i,j,k);
3843  qpg_arr(i,j,k) = qg3d(i,j,k);
3844  ng_arr(i,j,k) = ng3d(i,j,k);
3845 
3846  // Temperature and potential temperature conversion
3847  theta_arr(i,j,k) = t3d(i,j,k) / pii_arr(i,j,k); // Convert temp back to potential temp
3848  qv_arr(i,j,k) = qv3d(i,j,k);
3849 
3850  //Deleted wrf-check, effc, and precr type data as not used by ERF
3851  /*
3852  // NEED gpu-compatabile summation for rain_accum, check SAM or Kessler for better example
3853  rain_accum_arr(i,j,k) = rain_accum_arr(i,j,k) + precrt(i,j,k);
3854  snow_accum_arr(i,j,k) = snow_accum_arr(i,j,k) + snowprt(i,j,k);
3855  graup_accum_arr(i,j,k) = graup_accum_arr(i,j,k) + grplprt(i,j,k);*/
3856  rainncv_arr(i,j,0) = precrt(i,j,klo);
3857  snowncv_arr(i,j,0) = snowprt(i,j,klo);
3858  graupelncv_arr(i,j,0) = grplprt(i,j,klo);
3859  sr_arr(i,j,0) = snowrt(i,j,klo) / (precrt(i,j,klo) + 1.e-12);
3860  }
3861  // Update precipitation accumulation variables
3862  // These are outside the k-loop in the original code
3863  rain_accum_arr(i,j,klo) = rain_accum_arr(i,j,klo) + precrt(i,j,klo);
3864  snow_accum_arr(i,j,klo) = snow_accum_arr(i,j,klo) + snowprt(i,j,klo);
3865  graup_accum_arr(i,j,klo) = graup_accum_arr(i,j,klo) + grplprt(i,j,klo);
3866  }
3867  });
3868  // amrex::Print()<<FArrayBox(qv_arr)<<std::endl;
3869  }
3870 
3871  amrex::Print()<<"fortran should run "<<run_morr_fort<<std::endl;
3872 
3873  if(run_morr_fort) {
3874 #ifdef ERF_USE_MORR_FORT
3875  mp_morr_two_moment_c
3876  (
3877  1, // ITIMESTEP - Use 1 for simplicity
3878 
3879  // 3D arrays in Fortran expected order (assume column-major for Fortran)
3880  theta_arr.dataPtr(), // TH
3881  qv_arr.dataPtr(), // QV
3882  qcl_arr.dataPtr(), // QC
3883  qpr_arr.dataPtr(), // QR
3884  qci_arr.dataPtr(), // QI
3885  qps_arr.dataPtr(), // QS
3886  qpg_arr.dataPtr(), // QG
3887  ni_arr.dataPtr(), // NI
3888  ns_arr.dataPtr(), // NS
3889  nr_arr.dataPtr(), // NR
3890  ng_arr.dataPtr(), // NG
3891 
3892  rho_arr.dataPtr(), // RHO
3893  pii_arr.dataPtr(), // PII (Exner function)
3894  pres_arr.dataPtr(), // P (in hPa, convert if needed)
3895  dt, // DT_IN
3896  dz_arr.dataPtr(), // DZ
3897  w_arr.dataPtr(), // W (vertical velocity)
3898 
3899  // 2D arrays for precipitation accounting
3900  rain_accum_arr.dataPtr(), // RAINNC
3901  rainncv_arr.dataPtr(), // RAINNCV
3902  sr_arr.dataPtr(), // SR
3903  snow_accum_arr.dataPtr(), // SNOWNC
3904  snowncv_arr.dataPtr(), // SNOWNCV
3905  graup_accum_arr.dataPtr(),// GRAUPELNC
3906  graupelncv_arr.dataPtr(), // GRAUPELNCV
3907 
3908  // Radar reflectivity
3909  dummy_reflectivity_ptr, // refl_10cm
3910  true, // diagflag
3911  0, // do_radar_ref
3912 
3913  // Cumulus tendencies
3914  qrcuten_arr.dataPtr(), // qrcuten
3915  qscuten_arr.dataPtr(), // qscuten
3916  qicuten_arr.dataPtr(), // qicuten
3917 
3918  // WRF-Chem flags
3919  flag_qndrop, // F_QNDROP
3920  nullptr, // qndrop (not used here)
3921  ht_arr.dataPtr(), // HT (terrain height - not used)
3922 
3923  // Domain dimensions
3924  ilo, ihi, jlo, jhi, klo, khi, // IDS,IDE,JDS,JDE,KDS,KDE
3925  ilom, ihim, jlom, jhim, klom, khim, // IMS,IME,JMS,JME,KMS,KME
3926  ilo, ihi, jlo, jhi, klo, khi, // ITS,ITE,JTS,JTE,KTS,KTE
3927 
3928  // Optional WRF-Chem outputs
3929  false, // wetscav_on
3930  rainprod_arr.dataPtr(), // rainprod
3931  evapprod_arr.dataPtr(), // evapprod
3932  qlsink_arr.dataPtr(), // QLSINK
3933  precr_arr.dataPtr(), // PRECR
3934  preci_arr.dataPtr(), // PRECI
3935  precs_arr.dataPtr(), // PRECS
3936  precg_arr.dataPtr() // PRECG
3937  );
3938 #else
3939  amrex::Abort("Trying to run fortran without compiling with USE_MORR_FORT=TRUE");
3940 #endif
3941  }
3942  // amrex::Print()<<FArrayBox(qv_arr)<<std::endl;
3943  // After the call, all fields are updated
3944  // We don't need to copy results back since we passed direct pointers
3945  // to our class member arrays
3946  }
3947  }
calc_saturation_vapor_pressure
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_saturation_vapor_pressure(const amrex::Real T, const int type)
Definition: ERF_AdvanceMorrison.cpp:317
gamma_function
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real gamma_function(Real x)
Definition: ERF_AdvanceMorrison.cpp:304
pp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:219
mp_morr_two_moment_c
void mp_morr_two_moment_c(int itimestep, double *th, double *qv, double *qc, double *qr, double *qi, double *qs, double *qg, double *ni, double *ns, double *nr, double *ng, double *rho, double *pii, double *p, double dt_in, double *dz, double *w, double *rainnc, double *rainncv, double *sr, double *snownc, double *snowncv, double *graupelnc, double *graupelncv, double *refl_10cm, bool diagflag, int do_radar_ref, double *qrcuten, double *qscuten, double *qicuten, bool f_qndrop, double *qndrop, double *ht, int ids, int ide, int jds, int jde, int kds, int kde, int ims, int ime, int jms, int jme, int kms, int kme, int its, int ite, int jts, int jte, int kts, int kte, bool wetscav_on, double *rainprod, double *evapprod, double *qlsink, double *precr, double *preci, double *precs, double *precg)
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
Morrison::m_geom
amrex::Geometry m_geom
Definition: ERF_Morrison.H:278
Morrison::m_axis
int m_axis
Definition: ERF_Morrison.H:287
Morrison::m_rdOcp
amrex::Real m_rdOcp
Definition: ERF_Morrison.H:294
Morrison::mic_fab_vars
amrex::Array< FabPtr, MicVar_Morr::NumVars > mic_fab_vars
Definition: ERF_Morrison.H:301
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:25
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22
MicVar_Morr::qv
@ qv
Definition: ERF_Morrison.H:34
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:48
MicVar_Morr::nc
@ nc
Definition: ERF_Morrison.H:44
MicVar_Morr::qpg
@ qpg
Definition: ERF_Morrison.H:41
MicVar_Morr::pres
@ pres
Definition: ERF_Morrison.H:30
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:45
MicVar_Morr::qcl
@ qcl
Definition: ERF_Morrison.H:35
MicVar_Morr::tabs
@ tabs
Definition: ERF_Morrison.H:29
MicVar_Morr::theta
@ theta
Definition: ERF_Morrison.H:28
MicVar_Morr::ni
@ ni
Definition: ERF_Morrison.H:46
MicVar_Morr::ns
@ ns
Definition: ERF_Morrison.H:47
MicVar_Morr::omega
@ omega
Definition: ERF_Morrison.H:53
MicVar_Morr::qps
@ qps
Definition: ERF_Morrison.H:40
MicVar_Morr::graup_accum
@ graup_accum
Definition: ERF_Morrison.H:52
MicVar_Morr::rho
@ rho
Definition: ERF_Morrison.H:27
MicVar_Morr::qpr
@ qpr
Definition: ERF_Morrison.H:39
MicVar_Morr::qci
@ qci
Definition: ERF_Morrison.H:36
MicVar_Morr::rain_accum
@ rain_accum
Definition: ERF_Morrison.H:50
MicVar_Morr::snow_accum
@ snow_accum
Definition: ERF_Morrison.H:51
module_model_constants::p0
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
module_model_constants::xlf
real(c_double), parameter xlf
Definition: ERF_module_model_constants.F90:58
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:787
Here is the call graph for this function:

◆ Compute_Coefficients()

void Morrison::Compute_Coefficients ( )
153 {
154  auto dz = m_geom.CellSize(2);
155  auto lowz = m_geom.ProbLo(2);
156 
157  auto accrrc_t = accrrc.table();
158  auto accrsi_t = accrsi.table();
159  auto accrsc_t = accrsc.table();
160  auto coefice_t = coefice.table();
161  auto evaps1_t = evaps1.table();
162  auto evaps2_t = evaps2.table();
163  auto accrgi_t = accrgi.table();
164  auto accrgc_t = accrgc.table();
165  auto evapg1_t = evapg1.table();
166  auto evapg2_t = evapg2.table();
167  auto evapr1_t = evapr1.table();
168  auto evapr2_t = evapr2.table();
169 
170  auto rho1d_t = rho1d.table();
171  auto pres1d_t = pres1d.table();
172  auto tabs1d_t = tabs1d.table();
173 
174  auto gamaz_t = gamaz.table();
175  auto zmid_t = zmid.table();
176 
177  Real gam3 = erf_gammafff(3.0 );
178  Real gamr1 = erf_gammafff(3.0+b_rain );
179  Real gamr2 = erf_gammafff((5.0+b_rain)/2.0);
180  Real gams1 = erf_gammafff(3.0+b_snow );
181  Real gams2 = erf_gammafff((5.0+b_snow)/2.0);
182  Real gamg1 = erf_gammafff(3.0+b_grau );
183  Real gamg2 = erf_gammafff((5.0+b_grau)/2.0);
184 
185  // calculate the plane average variables
186  PlaneAverage rho_ave(mic_fab_vars[MicVar_Morr::rho].get(), m_geom, m_axis);
187  PlaneAverage theta_ave(mic_fab_vars[MicVar_Morr::theta].get(), m_geom, m_axis);
188  PlaneAverage qv_ave(mic_fab_vars[MicVar_Morr::qv].get(), m_geom, m_axis);
189  rho_ave.compute_averages(ZDir(), rho_ave.field());
190  theta_ave.compute_averages(ZDir(), theta_ave.field());
191  qv_ave.compute_averages(ZDir(), qv_ave.field());
192 
193  // get host variable rho, and rhotheta
194  int ncell = rho_ave.ncell_line();
195 
196  Gpu::HostVector<Real> rho_h(ncell), theta_h(ncell), qv_h(ncell);
197  rho_ave.line_average(0, rho_h);
198  theta_ave.line_average(0, theta_h);
199  qv_ave.line_average(0, qv_h);
200 
201  // copy data to device
202  Gpu::DeviceVector<Real> rho_d(ncell), theta_d(ncell), qv_d(ncell);
203  Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
204  Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
205  Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
206  Gpu::streamSynchronize();
207 
208  Real* rho_dptr = rho_d.data();
209  Real* theta_dptr = theta_d.data();
210  Real* qv_dptr = qv_d.data();
211 
212  Real gOcp = m_gOcp;
213 
214  ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
215  {
216  Real RhoTheta = rho_dptr[k]*theta_dptr[k];
217  Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]);
218  rho1d_t(k) = rho_dptr[k];
219  pres1d_t(k) = pressure*0.01;
220  // NOTE: Limit the temperature to the melting point of ice to avoid a divide by
221  // 0 condition when computing the cold evaporation coefficients. This should
222  // not affect results since evporation requires snow/graupel to be present
223  // and thus T<273.16
224  tabs1d_t(k) = std::min(getTgivenRandRTh(rho_dptr[k], RhoTheta, qv_dptr[k]),273.16);
225  zmid_t(k) = lowz + (k+0.5)*dz;
226  gamaz_t(k) = gOcp*zmid_t(k);
227  });
228 
229  if(round(gam3) != 2) {
230  std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl;
231  std::exit(-1);
232  }
233 
234  // Populate all the coefficients
235  ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
236  {
237  Real Prefactor;
238  Real pratio = sqrt(1.29 / rho1d_t(k));
239  //Real rrr1 = 393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5);
240  //Real rrr2 = std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k));
241  Real estw = 100.0*erf_esatw(tabs1d_t(k));
242  Real esti = 100.0*erf_esati(tabs1d_t(k));
243 
244  // accretion by snow:
245  Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0));
246  Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15));
247  accrsi_t(k) = coef1 * coef2 * esicoef;
248  accrsc_t(k) = coef1 * esccoef;
249  coefice_t(k) = coef2;
250 
251  // evaporation of snow:
252  coef1 = (lsub/(tabs1d_t(k)*R_v)-1.0)*lsub/(therco*tabs1d_t(k));
253  coef2 = R_v * R_d / (diffelq * esti);
254  Prefactor = 2.0 * PI * nzeros / (rho1d_t(k) * (coef1 + coef2));
255  Prefactor *= (2.0/PI); // Shape factor snow
256  evaps1_t(k) = Prefactor * 0.65 * sqrt(rho1d_t(k) / (PI * rhos * nzeros));
257  evaps2_t(k) = Prefactor * 0.44 * sqrt(a_snow * rho1d_t(k) / muelq) * gams2
258  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhos * nzeros) , ((5.0+b_snow)/8.0));
259 
260  // accretion by graupel:
261  coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0));
262  coef2 = exp(0.025*(tabs1d_t(k) - 273.15));
263  accrgi_t(k) = coef1 * coef2 * egicoef;
264  accrgc_t(k) = coef1 * egccoef;
265 
266  // evaporation of graupel:
267  coef1 = (lsub/(tabs1d_t(k)*R_v)-1.0)*lsub/(therco*tabs1d_t(k));
268  coef2 = R_v * R_d / (diffelq * esti);
269  Prefactor = 2.0 * PI * nzerog / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for graupel is 1
270  evapg1_t(k) = Prefactor * 0.78 * sqrt(rho1d_t(k) / (PI * rhog * nzerog));
271  evapg2_t(k) = Prefactor * 0.31 * sqrt(a_grau * rho1d_t(k) / muelq) * gamg2
272  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhog * nzerog) , ((5.0+b_grau)/8.0));
273 
274  // accretion by rain:
275  accrrc_t(k) = 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3.0+b_rain)/4.))* erccoef;
276 
277  // evaporation of rain:
278  coef1 = (lcond/(tabs1d_t(k)*R_v)-1.0)*lcond/(therco*tabs1d_t(k));
279  coef2 = R_v * R_d / (diffelq * estw);
280  Prefactor = 2.0 * PI * nzeror / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for rain is 1
281  evapr1_t(k) = Prefactor * 0.78 * sqrt(rho1d_t(k) / (PI * rhor * nzeror));
282  evapr2_t(k) = Prefactor * 0.31 * sqrt(a_rain * rho1d_t(k) / muelq) * gamr2
283  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhor * nzeror) , ((5.0+b_rain)/8.0));
284  });
285 }
rhog
constexpr amrex::Real rhog
Definition: ERF_Constants.H:30
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:11
muelq
constexpr amrex::Real muelq
Definition: ERF_Constants.H:75
nzerog
constexpr amrex::Real nzerog
Definition: ERF_Constants.H:59
a_grau
constexpr amrex::Real a_grau
Definition: ERF_Constants.H:42
lsub
constexpr amrex::Real lsub
Definition: ERF_Constants.H:68
esicoef
constexpr amrex::Real esicoef
Definition: ERF_Constants.H:53
diffelq
constexpr amrex::Real diffelq
Definition: ERF_Constants.H:73
therco
constexpr amrex::Real therco
Definition: ERF_Constants.H:74
b_grau
constexpr amrex::Real b_grau
Definition: ERF_Constants.H:43
egccoef
constexpr amrex::Real egccoef
Definition: ERF_Constants.H:54
egicoef
constexpr amrex::Real egicoef
Definition: ERF_Constants.H:55
b_rain
constexpr amrex::Real b_rain
Definition: ERF_Constants.H:39
lcond
constexpr amrex::Real lcond
Definition: ERF_Constants.H:66
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
nzeror
constexpr amrex::Real nzeror
Definition: ERF_Constants.H:57
rhos
constexpr amrex::Real rhos
Definition: ERF_Constants.H:29
esccoef
constexpr amrex::Real esccoef
Definition: ERF_Constants.H:52
rhor
constexpr amrex::Real rhor
Definition: ERF_Constants.H:28
a_rain
constexpr amrex::Real a_rain
Definition: ERF_Constants.H:38
nzeros
constexpr amrex::Real nzeros
Definition: ERF_Constants.H:58
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
a_snow
constexpr amrex::Real a_snow
Definition: ERF_Constants.H:40
b_snow
constexpr amrex::Real b_snow
Definition: ERF_Constants.H:41
erccoef
constexpr amrex::Real erccoef
Definition: ERF_Constants.H:51
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:84
getTgivenRandRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:46
erf_esatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw(amrex::Real t)
Definition: ERF_MicrophysicsUtils.H:64
erf_esati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esati(amrex::Real t)
Definition: ERF_MicrophysicsUtils.H:24
erf_gammafff
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_gammafff(amrex::Real x)
Definition: ERF_MicrophysicsUtils.H:15
Morrison::evapg2
amrex::TableData< amrex::Real, 1 > evapg2
Definition: ERF_Morrison.H:313
Morrison::pres1d
amrex::TableData< amrex::Real, 1 > pres1d
Definition: ERF_Morrison.H:319
Morrison::zmid
amrex::TableData< amrex::Real, 1 > zmid
Definition: ERF_Morrison.H:326
Morrison::m_gOcp
amrex::Real m_gOcp
Definition: ERF_Morrison.H:293
Morrison::accrgi
amrex::TableData< amrex::Real, 1 > accrgi
Definition: ERF_Morrison.H:310
Morrison::accrgc
amrex::TableData< amrex::Real, 1 > accrgc
Definition: ERF_Morrison.H:311
Morrison::tabs1d
amrex::TableData< amrex::Real, 1 > tabs1d
Definition: ERF_Morrison.H:320
Morrison::nlev
int nlev
Definition: ERF_Morrison.H:284
Morrison::gamaz
amrex::TableData< amrex::Real, 1 > gamaz
Definition: ERF_Morrison.H:325
Morrison::coefice
amrex::TableData< amrex::Real, 1 > coefice
Definition: ERF_Morrison.H:307
Morrison::accrsi
amrex::TableData< amrex::Real, 1 > accrsi
Definition: ERF_Morrison.H:305
Morrison::accrrc
amrex::TableData< amrex::Real, 1 > accrrc
Definition: ERF_Morrison.H:304
Morrison::evapg1
amrex::TableData< amrex::Real, 1 > evapg1
Definition: ERF_Morrison.H:312
Morrison::evapr1
amrex::TableData< amrex::Real, 1 > evapr1
Definition: ERF_Morrison.H:314
Morrison::evaps1
amrex::TableData< amrex::Real, 1 > evaps1
Definition: ERF_Morrison.H:308
Morrison::evaps2
amrex::TableData< amrex::Real, 1 > evaps2
Definition: ERF_Morrison.H:309
Morrison::rho1d
amrex::TableData< amrex::Real, 1 > rho1d
Definition: ERF_Morrison.H:318
Morrison::evapr2
amrex::TableData< amrex::Real, 1 > evapr2
Definition: ERF_Morrison.H:315
Morrison::accrsc
amrex::TableData< amrex::Real, 1 > accrsc
Definition: ERF_Morrison.H:306
PlaneAverage
Definition: ERF_PlaneAverage.H:14

Referenced by Update_Micro_Vars().

Here is the call graph for this function:
Here is the caller 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  // get potential total density, temperature, qt, qp
37  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
38  {
39  states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
40 
41  states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*std::max(0.0,qv_arr(i,j,k));
42  states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*std::max(0.0,qc_arr(i,j,k));
43  states_arr(i,j,k,RhoQ3_comp) = rho_arr(i,j,k)*std::max(0.0,qi_arr(i,j,k));
44 
45  states_arr(i,j,k,RhoQ4_comp) = rho_arr(i,j,k)*std::max(0.0,qpr_arr(i,j,k));
46  states_arr(i,j,k,RhoQ5_comp) = rho_arr(i,j,k)*std::max(0.0,qps_arr(i,j,k));
47  states_arr(i,j,k,RhoQ6_comp) = rho_arr(i,j,k)*std::max(0.0,qpg_arr(i,j,k));
48  });
49  }
50 
51  // Fill interior ghost cells and periodic boundaries
52  cons.FillBoundary(m_geom.periodicity());
53 }
RhoQ4_comp
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
RhoQ3_comp
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:44
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
RhoQ6_comp
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:47
RhoQ5_comp
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:46
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140

Referenced by Update_State_Vars().

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.

99 {
100  // Get the temperature, density, theta, qt and qp from input
101  for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
102  const auto& box3d = mfi.growntilebox();
103 
104  auto states_array = cons_in.array(mfi);
105 
106  // Non-precipitating
107  auto qv_array = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
108  auto qc_array = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
109  auto qi_array = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
110  auto qn_array = mic_fab_vars[MicVar_Morr::qn]->array(mfi);
111  auto qt_array = mic_fab_vars[MicVar_Morr::qt]->array(mfi);
112 
113  // Precipitating
114  auto qpr_array = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
115  auto qps_array = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
116  auto qpg_array = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
117  auto qp_array = mic_fab_vars[MicVar_Morr::qp]->array(mfi);
118 
119  auto rho_array = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
120  auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
121  auto tabs_array = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
122  auto pres_array = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
123 
124  // Get pressure, theta, temperature, density, and qt, qp
125  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
126  {
127  rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
128  theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
129 
130  qv_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp));
131  qc_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp));
132  qi_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp));
133  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
134  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
135 
136  qpr_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ4_comp)/states_array(i,j,k,Rho_comp));
137  qps_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ5_comp)/states_array(i,j,k,Rho_comp));
138  qpg_array(i,j,k) = std::max(0.0,states_array(i,j,k,RhoQ6_comp)/states_array(i,j,k,Rho_comp));
139  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
140 
141  tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
142  states_array(i,j,k,RhoTheta_comp),
143  qv_array(i,j,k));
144 
145  // NOTE: the Morrison Fortran version uses Pa not hPa so we don't divideby 100!
146  pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); // * 0.01;
147  });
148  }
149 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
MicVar_Morr::qp
@ qp
Definition: ERF_Morrison.H:38
MicVar_Morr::qn
@ qn
Definition: ERF_Morrison.H:33
MicVar_Morr::qt
@ 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_fac_cond = lcond / sc.c_p;
74  m_fac_fus = lfus / sc.c_p;
75  m_fac_sub = lsub / sc.c_p;
76  m_gOcp = CONST_GRAV / sc.c_p;
77  m_axis = sc.ave_plane;
78  m_rdOcp = sc.rdOcp;
79  }
lfus
constexpr amrex::Real lfus
Definition: ERF_Constants.H:67
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
Morrison::m_fac_cond
amrex::Real m_fac_cond
Definition: ERF_Morrison.H:290
Morrison::m_fac_sub
amrex::Real m_fac_sub
Definition: ERF_Morrison.H:292
Morrison::m_fac_fus
amrex::Real m_fac_fus
Definition: ERF_Morrison.H:291
SolverChoice::rdOcp
amrex::Real rdOcp
Definition: ERF_DataStruct.H:742
SolverChoice::c_p
amrex::Real c_p
Definition: ERF_DataStruct.H:741
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:802

◆ 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.

28 {
29  [[maybe_unused]] amrex::Real dt = dt_advance;
30  m_geom = geom;
31  m_gtoe = grids;
32 
33  m_z_phys_nd = z_phys_nd.get();
34  m_detJ_cc = detJ_cc.get();
35 
36  MicVarMap.resize(m_qmoist_size);
37  MicVarMap = {MicVar_Morr::nc, MicVar_Morr::nr, MicVar_Morr::ni, MicVar_Morr::ns, MicVar_Morr::ng,
38  MicVar_Morr::rain_accum, MicVar_Morr::snow_accum, MicVar_Morr::graup_accum};
39 
40  // initialize microphysics variables
41  for (auto ivar = 0; ivar < MicVar_Morr::NumVars; ++ivar) {
42  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
43  1, cons_in.nGrowVect());
44  mic_fab_vars[ivar]->setVal(0.);
45  }
46 
47  // Set class data members
48  for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
49  const auto& box3d = mfi.tilebox();
50 
51  const auto& lo = lbound(box3d);
52  const auto& hi = ubound(box3d);
53 
54  nlev = box3d.length(2);
55  zlo = lo.z;
56  zhi = hi.z;
57 
58  // parameters
59  accrrc.resize({zlo}, {zhi});
60  accrsi.resize({zlo}, {zhi});
61  accrsc.resize({zlo}, {zhi});
62  coefice.resize({zlo}, {zhi});
63  evaps1.resize({zlo}, {zhi});
64  evaps2.resize({zlo}, {zhi});
65  accrgi.resize({zlo}, {zhi});
66  accrgc.resize({zlo}, {zhi});
67  evapg1.resize({zlo}, {zhi});
68  evapg2.resize({zlo}, {zhi});
69  evapr1.resize({zlo}, {zhi});
70  evapr2.resize({zlo}, {zhi});
71 
72  // data (input)
73  rho1d.resize({zlo}, {zhi});
74  pres1d.resize({zlo}, {zhi});
75  tabs1d.resize({zlo}, {zhi});
76  gamaz.resize({zlo}, {zhi});
77  zmid.resize({zlo}, {zhi});
78  }
79 
80 #ifdef ERF_USE_MORR_FORT
81  int morr_rimed_ice = 0; // This is used to set something called "ihail"
82  amrex::ParmParse pp("erf");
83  MoistureType moisture_type;
84  pp.query_enum_case_insensitive("moisture_model",moisture_type);
85  int morr_noice = (moisture_type == MoistureType::Morrison_NoIce);
86  Print()<<"Setting No Ice flag in fortran to "<<morr_noice<<std::endl;
87  morr_two_moment_init_c(morr_rimed_ice, morr_noice);
88 #endif
89 }
morr_two_moment_init_c
void morr_two_moment_init_c(int morr_rimed_ice, int morr_noice)
Morrison::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_Morrison.H:275
Morrison::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_Morrison.H:297
Morrison::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_Morrison.H:298
Morrison::m_qmoist_size
int m_qmoist_size
Definition: ERF_Morrison.H:266
Morrison::zlo
int zlo
Definition: ERF_Morrison.H:284
Morrison::m_gtoe
amrex::BoxArray m_gtoe
Definition: ERF_Morrison.H:281
Morrison::zhi
int zhi
Definition: ERF_Morrison.H:284
MicVar_Morr::NumVars
@ NumVars
Definition: ERF_Morrison.H:54
Here is the call graph for this function:

◆ NewtonIterSat()

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Morrison::NewtonIterSat ( int &  i,
int &  j,
int &  k,
const amrex::Real &  fac_cond,
const amrex::Real &  fac_fus,
const amrex::Real &  fac_sub,
const amrex::Real &  an,
const amrex::Real &  bn,
const amrex::Array4< amrex::Real > &  tabs_array,
const amrex::Array4< amrex::Real > &  pres_array,
const amrex::Array4< amrex::Real > &  qv_array,
const amrex::Array4< amrex::Real > &  qc_array,
const amrex::Array4< amrex::Real > &  qi_array,
const amrex::Array4< amrex::Real > &  qn_array,
const amrex::Array4< amrex::Real > &  qt_array 
)
inlinestatic
148  {
149  // Solution tolerance
150  amrex::Real tol = 1.0e-4;
151 
152  // Saturation moisture fractions
153  amrex::Real omn, domn;
154  amrex::Real qsat, dqsat;
155  amrex::Real qsatw, dqsatw;
156  amrex::Real qsati, dqsati;
157 
158  // Newton iteration vars
159  int niter;
160  amrex::Real fff, dfff, dtabs;
161  amrex::Real lstar, dlstar;
162  amrex::Real lstarw, lstari;
163  amrex::Real delta_qv, delta_qc, delta_qi;
164 
165  // Initial guess for temperature & pressure
166  amrex::Real tabs = tabs_array(i,j,k);
167  amrex::Real pres = pres_array(i,j,k);
168 
169  niter = 0;
170  dtabs = 1;
171  //==================================================
172  // Newton iteration to qv=qsat (cloud phase only)
173  //==================================================
174  do {
175  // Latent heats and their derivatives wrt to T
176  lstarw = fac_cond;
177  lstari = fac_fus;
178  domn = 0.0;
179 
180  // Saturation moisture fractions
181  erf_qsatw(tabs, pres, qsatw);
182  erf_qsati(tabs, pres, qsati);
183  erf_dtqsatw(tabs, pres, dqsatw);
184  erf_dtqsati(tabs, pres, dqsati);
185 
186  // Cloud ice not permitted (condensation & fusion)
187  if(tabs >= tbgmax) {
188  omn = 1.0;
189  }
190  // Cloud water not permitted (sublimation & fusion)
191  else if(tabs <= tbgmin) {
192  omn = 0.0;
193  lstarw = fac_sub;
194  }
195  // Mixed cloud phase (condensation & fusion)
196  else {
197  omn = an*tabs-bn;
198  domn = an;
199  }
200 
201  // Linear combination of each component
202  qsat = omn * qsatw + (1.0-omn) * qsati;
203  dqsat = omn * dqsatw + (1.0-omn) * dqsati
204  + domn * qsatw - domn * qsati;
205  lstar = omn * lstarw + (1.0-omn) * lstari;
206  dlstar = domn * lstarw - domn * lstari;
207 
208  // Function for root finding:
209  // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
210  fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat);
211 
212  // Derivative of function (T_new iterated on)
213  dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;
214 
215  // Update the temperature
216  dtabs = -fff/dfff;
217  tabs += dtabs;
218 
219  // Update iteration
220  niter = niter+1;
221  } while(std::abs(dtabs) > tol && niter < 20);
222 
223  // Update qsat from last iteration (dq = dq/dt * dt)
224  qsat += dqsat*dtabs;
225 
226  // Changes in each component
227  delta_qv = qv_array(i,j,k) - qsat;
228  delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn);
229  delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn));
230 
231  // Partition the change in non-precipitating q
232  qv_array(i,j,k) = qsat;
233  qc_array(i,j,k) += delta_qc;
234  qi_array(i,j,k) += delta_qi;
235  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
236  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
237 
238  // Return to temperature
239  return tabs;
240  }
tbgmax
constexpr amrex::Real tbgmax
Definition: ERF_Constants.H:32
tbgmin
constexpr amrex::Real tbgmin
Definition: ERF_Constants.H:31
erf_dtqsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:170
erf_dtqsati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsati(amrex::Real t, amrex::Real p, amrex::Real &dtqsati)
Definition: ERF_MicrophysicsUtils.H:165
erf_qsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:158
erf_qsati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsati(amrex::Real t, amrex::Real p, amrex::Real &qsati)
Definition: ERF_MicrophysicsUtils.H:151
Here is the call graph for this function:

◆ Qmoist_Ptr()

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

Reimplemented from NullMoist.

118  {
119  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
120  return mic_fab_vars[MicVarMap[varIdx]].get();
121  }

◆ 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.

246  {
247  a_idx.clear();
248  a_names.clear();
249 
250  // The ordering here needs to match that in
251  // MicVarMap = {MicVar_Morr::nc, MicVar_Morr::nr, MicVar_Morr::ni, MicVar_Morr::ns, MicVar_Morr::ng,
252  // MicVar_Morr::rain_accum, MicVar_Morr::snow_accum, MicVar_Morr::graup_accum};
253  //
254  a_idx.push_back(0); a_names.push_back("Nc");
255  a_idx.push_back(1); a_names.push_back("Nr");
256  a_idx.push_back(2); a_names.push_back("Ni");
257  a_idx.push_back(3); a_names.push_back("Ns");
258  a_idx.push_back(4); a_names.push_back("Ng");
259  a_idx.push_back(5); a_names.push_back("RainAccum");
260  a_idx.push_back(6); a_names.push_back("SnowAccum");
261  a_idx.push_back(7); a_names.push_back("GraupAccum");
262  }

◆ Qmoist_Size()

int Morrison::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

127 { return Morrison::m_qmoist_size; }

◆ Qstate_Moist_Size()

int Morrison::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

130 { return Morrison::n_qstate_moist_size; }
Morrison::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_Morrison.H:269

◆ Update_Micro_Vars()

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

Reimplemented from NullMoist.

100  {
101  this->Copy_State_to_Micro(cons_in);
102  this->Compute_Coefficients();
103  }
Morrison::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitMorrison.cpp:98
Morrison::Compute_Coefficients
void Compute_Coefficients()
Definition: ERF_InitMorrison.cpp:152
Here is the call graph for this function:

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

107  {
108  this->Copy_Micro_to_State(cons_in);
109  }
Morrison::Copy_Micro_to_State
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

◆ accrgc

amrex::TableData<amrex::Real, 1> Morrison::accrgc
private

◆ accrgi

amrex::TableData<amrex::Real, 1> Morrison::accrgi
private

◆ accrrc

amrex::TableData<amrex::Real, 1> Morrison::accrrc
private

◆ accrsc

amrex::TableData<amrex::Real, 1> Morrison::accrsc
private

◆ accrsi

amrex::TableData<amrex::Real, 1> Morrison::accrsi
private

◆ CFL_MAX

constexpr amrex::Real Morrison::CFL_MAX = 0.5
staticconstexprprivate

◆ coefice

amrex::TableData<amrex::Real, 1> Morrison::coefice
private

◆ evapg1

amrex::TableData<amrex::Real, 1> Morrison::evapg1
private

◆ evapg2

amrex::TableData<amrex::Real, 1> Morrison::evapg2
private

◆ evapr1

amrex::TableData<amrex::Real, 1> Morrison::evapr1
private

◆ evapr2

amrex::TableData<amrex::Real, 1> Morrison::evapr2
private

◆ evaps1

amrex::TableData<amrex::Real, 1> Morrison::evaps1
private

◆ evaps2

amrex::TableData<amrex::Real, 1> Morrison::evaps2
private

◆ gamaz

amrex::TableData<amrex::Real, 1> Morrison::gamaz
private

◆ m_axis

int Morrison::m_axis
private

Referenced by Define().

◆ m_detJ_cc

amrex::MultiFab* Morrison::m_detJ_cc
private

◆ m_fac_cond

amrex::Real Morrison::m_fac_cond
private

Referenced by Define().

◆ m_fac_fus

amrex::Real Morrison::m_fac_fus
private

Referenced by Define().

◆ m_fac_sub

amrex::Real Morrison::m_fac_sub
private

Referenced by Define().

◆ m_geom

amrex::Geometry Morrison::m_geom
private

◆ m_gOcp

amrex::Real Morrison::m_gOcp
private

Referenced by Define().

◆ m_gtoe

amrex::BoxArray Morrison::m_gtoe
private

◆ m_qmoist_size

int Morrison::m_qmoist_size = 8
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_size

int Morrison::n_qstate_moist_size = 6
private

Referenced by Qstate_Moist_Size().

◆ nlev

int Morrison::nlev
private

◆ pres1d

amrex::TableData<amrex::Real, 1> Morrison::pres1d
private

◆ qn1d

amrex::TableData<amrex::Real, 1> Morrison::qn1d
private

◆ qt1d

amrex::TableData<amrex::Real, 1> Morrison::qt1d
private

◆ qv1d

amrex::TableData<amrex::Real, 1> Morrison::qv1d
private

◆ rho1d

amrex::TableData<amrex::Real, 1> Morrison::rho1d
private

◆ tabs1d

amrex::TableData<amrex::Real, 1> Morrison::tabs1d
private

◆ zhi

int Morrison::zhi
private

◆ zlo

int Morrison::zlo
private

◆ zmid

amrex::TableData<amrex::Real, 1> Morrison::zmid
private
The documentation for this class was generated from the following files: