ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Static Public Attributes | Private Types | Private Member Functions | Private Attributes | List of all members
WDM6 Class Reference

#include <ERF_WDM6.H>

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

Public Member Functions

 WDM6 ()
 
virtual ~WDM6 ()=default
 
void Define (SolverChoice &sc) override
 
void Init (const amrex::MultiFab &cons_in, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &z_phys_nd, std::unique_ptr< amrex::MultiFab > &detJ_cc) override
 
void Set_dzmin (const amrex::Real dz_min) override
 
void Set_Lmask (amrex::iMultiFab *lmask) override
 Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this. More...
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in, const amrex::MultiFab *base_state)
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *base_state) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &solverChoice) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
int Qstate_Moist_NumConc_Size () override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
SurfacePrecipAccumulationSources Get_Surface_Precip_Accumulation_Ptrs (const int &) 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
 
virtual void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const
 
virtual void SetCurrentLevel (const int &lev)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_RealWidth (const int)
 

Static Public Attributes

static constexpr amrex::Real dtcldcr = wdm6_literal(120.0)
 
static constexpr amrex::Real n0r = wdm6_literal(8.0e6)
 
static constexpr amrex::Real avtr = wdm6_literal(841.9)
 
static constexpr amrex::Real bvtr = wdm6_literal(0.8)
 
static constexpr amrex::Real r0 = wdm6_literal(0.8e-5)
 
static constexpr amrex::Real peaut = wdm6_literal(0.55)
 
static constexpr amrex::Real xncr = wdm6_literal(3.0e8)
 
static constexpr amrex::Real xncr0 = wdm6_literal(5.0e7)
 
static constexpr amrex::Real xncr1 = wdm6_literal(5.0e8)
 
static constexpr amrex::Real xmyu = wdm6_literal(1.718e-5)
 
static constexpr amrex::Real avts = wdm6_literal(11.72)
 
static constexpr amrex::Real bvts = wdm6_literal(0.41)
 
static constexpr amrex::Real lamdarmax = wdm6_literal(5.0e4)
 
static constexpr amrex::Real lamdarmin = wdm6_literal(2.0e3)
 
static constexpr amrex::Real lamdasmax = wdm6_literal(1.0e5)
 
static constexpr amrex::Real lamdacmax = wdm6_literal(5.0e5)
 
static constexpr amrex::Real lamdacmin = wdm6_literal(2.0e4)
 
static constexpr amrex::Real dicon = wdm6_literal(11.9)
 
static constexpr amrex::Real dimax = wdm6_literal(500.0e-6)
 
static constexpr amrex::Real pfrz1 = wdm6_literal(100.0)
 
static constexpr amrex::Real pfrz2 = wdm6_literal(0.66)
 
static constexpr amrex::Real qcrmin = wdm6_literal(1.0e-9)
 
static constexpr amrex::Real ncmin = wdm6_literal(1.0e1)
 
static constexpr amrex::Real nrmin = wdm6_literal(1.0e-2)
 
static constexpr amrex::Real eacrc = wdm6_literal(1.0)
 
static constexpr amrex::Real dens_snow = wdm6_literal(100.0)
 
static constexpr amrex::Real qs0 = wdm6_literal(6.0e-4)
 
static constexpr amrex::Real n0smax = wdm6_literal(1.0e11)
 
static constexpr amrex::Real n0s = wdm6_literal(2.0e6)
 
static constexpr amrex::Real alpha_wdm6 = wdm6_literal(0.12)
 
static constexpr amrex::Real satmax = wdm6_literal(1.0048)
 
static constexpr amrex::Real actk = wdm6_literal(0.6)
 
static constexpr amrex::Real actr = wdm6_literal(1.5)
 
static constexpr amrex::Real ncrk1 = wdm6_literal(3.03e3)
 
static constexpr amrex::Real ncrk2 = wdm6_literal(2.59e15)
 
static constexpr amrex::Real di100 = wdm6_literal(1.0e-4)
 
static constexpr amrex::Real di600 = wdm6_literal(6.0e-4)
 
static constexpr amrex::Real di2000 = wdm6_literal(2000.0e-6)
 
static constexpr amrex::Real di82 = wdm6_literal(82.0e-6)
 
static constexpr amrex::Real di15 = wdm6_literal(15.0e-6)
 

Private Types

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

Private Member Functions

void initialize_coeffs ()
 

Private Attributes

int m_qmoist_size = 3
 
int n_qstate_moist_size = 9
 
int n_qstate_moist_numconc_size = 3
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::Real dt {0.0}
 
amrex::Real m_dzmin {0.0}
 
amrex::Real m_ccn0 {100.0e6}
 
int nlev {0}
 
int zlo {0}
 
int zhi {0}
 
int m_axis {2}
 
bool m_do_cond {true}
 
amrex::Real m_rdOcp {RdoCp}
 
MoistureType m_moisture_type {MoistureType::None}
 
amrex::MultiFab * m_z_phys_nd {nullptr}
 
amrex::MultiFab * m_detJ_cc {nullptr}
 
amrex::iMultiFab * m_lmask {nullptr}
 
amrex::Array< FabPtr, MicVar_WDM6::NumVars > mic_fab_vars
 
bool m_hail_opt {false}
 
amrex::Real m_n0g {0}
 
amrex::Real m_deng {0}
 
amrex::Real m_avtg {0}
 
amrex::Real m_bvtg {0}
 
amrex::Real m_lamdagmax {0}
 
amrex::Real m_pi_wdm6 {0}
 
amrex::Real m_xlv1 {0}
 
amrex::Real m_qc0 {0}
 
amrex::Real m_qc1 {0}
 
amrex::Real m_qck1 {0}
 
amrex::Real m_pidnc {0}
 
amrex::Real m_bvtr1 {0}
 
amrex::Real m_bvtr2 {0}
 
amrex::Real m_bvtr3 {0}
 
amrex::Real m_bvtr4 {0}
 
amrex::Real m_bvtr5 {0}
 
amrex::Real m_bvtr6 {0}
 
amrex::Real m_bvtr7 {0}
 
amrex::Real m_bvtr2o5 {0}
 
amrex::Real m_bvtr3o5 {0}
 
amrex::Real m_g1pbr {0}
 
amrex::Real m_g2pbr {0}
 
amrex::Real m_g3pbr {0}
 
amrex::Real m_g4pbr {0}
 
amrex::Real m_g5pbr {0}
 
amrex::Real m_g6pbr {0}
 
amrex::Real m_g7pbr {0}
 
amrex::Real m_g5pbro2 {0}
 
amrex::Real m_g7pbro2 {0}
 
amrex::Real m_pvtr {0}
 
amrex::Real m_pvtrn {0}
 
amrex::Real m_eacrr {0}
 
amrex::Real m_pacrr {0}
 
amrex::Real m_precr1 {0}
 
amrex::Real m_precr2 {0}
 
amrex::Real m_xmmax {0}
 
amrex::Real m_roqimax {0}
 
amrex::Real m_pidn0r {0}
 
amrex::Real m_pidnr {0}
 
amrex::Real m_bvts1 {0}
 
amrex::Real m_bvts2 {0}
 
amrex::Real m_bvts3 {0}
 
amrex::Real m_bvts4 {0}
 
amrex::Real m_g1pbs {0}
 
amrex::Real m_g3pbs {0}
 
amrex::Real m_g4pbs {0}
 
amrex::Real m_g5pbso2 {0}
 
amrex::Real m_pvts {0}
 
amrex::Real m_pacrs {0}
 
amrex::Real m_precs1 {0}
 
amrex::Real m_precs2 {0}
 
amrex::Real m_pidn0s {0}
 
amrex::Real m_pacrc {0}
 
amrex::Real m_bvtg1 {0}
 
amrex::Real m_bvtg2 {0}
 
amrex::Real m_bvtg3 {0}
 
amrex::Real m_bvtg4 {0}
 
amrex::Real m_g1pbg {0}
 
amrex::Real m_g3pbg {0}
 
amrex::Real m_g4pbg {0}
 
amrex::Real m_g5pbgo2 {0}
 
amrex::Real m_pvtg {0}
 
amrex::Real m_pacrg {0}
 
amrex::Real m_precg1 {0}
 
amrex::Real m_precg2 {0}
 
amrex::Real m_pidn0g {0}
 
amrex::Real m_rslopecmax {0}
 
amrex::Real m_rslopec2max {0}
 
amrex::Real m_rslopec3max {0}
 
amrex::Real m_rslopermax {0}
 
amrex::Real m_rslopesmax {0}
 
amrex::Real m_rslopegmax {0}
 
amrex::Real m_rsloperbmax {0}
 
amrex::Real m_rslopesbmax {0}
 
amrex::Real m_rslopegbmax {0}
 
amrex::Real m_rsloper2max {0}
 
amrex::Real m_rslopes2max {0}
 
amrex::Real m_rslopeg2max {0}
 
amrex::Real m_rsloper3max {0}
 
amrex::Real m_rslopes3max {0}
 
amrex::Real m_rslopeg3max {0}
 

Additional Inherited Members

- Protected Member Functions inherited from NullMoist
void set_anelastic_reference_pressure_mode (const SolverChoice &sc)
 
void assert_base_state_available (const amrex::MultiFab *base_state) const
 
- Protected Attributes inherited from NullMoist
int m_level {0}
 
bool m_use_anelastic_reference_pressure {false}
 

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ WDM6()

WDM6::WDM6 ( )
inline
112 {}

◆ ~WDM6()

virtual WDM6::~WDM6 ( )
virtualdefault

Member Function Documentation

◆ Advance()

void WDM6::Advance ( const amrex::Real &  dt_advance,
const SolverChoice &  solverChoice 
)
overridevirtual

Reimplemented from NullMoist.

654 {
655  // ---------------------------------------------------------------
656  // Dual-mode implementation following WSM6 pattern:
657  // - With ERF_USE_WDM6_FORT: Call Fortran bridge (CPU-only)
658  // - Without: Use C++ GPU kernels (not yet implemented)
659  // ---------------------------------------------------------------
660 
661 #ifdef ERF_USE_WDM6_FORT
662  static int call_count = 0;
663  call_count++;
664  [[maybe_unused]] const bool first_call = (call_count == 1);
665 
666  // Fortran bridge mode - initialize once
667  static bool wdm6_inited = false;
668  if (!wdm6_inited) {
669  constexpr double den0 = 1.28; // Standard dry-air density (kg/m^3)
670  constexpr double denr = static_cast<double>(rhoh2o);
671  constexpr double dens = static_cast<double>(rhos);
672  constexpr double cl = static_cast<double>(Cp_l);
673  constexpr double cpv = static_cast<double>(Cp_v);
674  const double ccn0 = static_cast<double>(m_ccn0);
675  // Honour wdm6.hail_opt on this leg too. m_hail_opt is resolved in Init
676  // from the same input the native coefficients branch on, so the two
677  // legs cannot disagree about the regime. The bool maps back to {0,1}
678  // losslessly with respect to the Fortran, whose only test is
679  // `hail_opt .eq. 1` (ERF_module_mp_wdm6.F90:3259).
680  const int hail_opt = m_hail_opt ? 1 : 0;
681  mp_wdm6_init_c(den0, denr, dens, cl, cpv, ccn0, hail_opt);
682  wdm6_inited = true;
683  if (first_call) {
684  amrex::Print() << "WDM6 Fortran bridge initialized\n";
685  }
686  }
687 #endif
688 
689  int microphysics_debug = 0;
690  std::vector<int> micro_diag_target_column;
691  {
692  amrex::ParmParse pp("erf");
693  pp.queryAdd("microphysics_debug", microphysics_debug);
694  pp.queryarr("micro_diag_target_column", micro_diag_target_column);
695  }
696  microphysics_debug = std::max(0, std::min(2, microphysics_debug));
697 #ifdef ERF_USE_WDM6_FORT
698  bool use_wdm6_cpp_answer = false;
699  {
700  amrex::ParmParse pp("erf");
701  pp.queryAdd("use_wdm6_cpp_answer", use_wdm6_cpp_answer);
702  }
703  const bool run_wdm6_fort = !use_wdm6_cpp_answer;
704 #endif
705 
706  // Physical constants
707  [[maybe_unused]] constexpr double g = static_cast<double>(CONST_GRAV);
708  constexpr double cpd = static_cast<double>(Cp_d);
709  constexpr double cpv = static_cast<double>(Cp_v);
710  [[maybe_unused]] constexpr double rd = static_cast<double>(R_d);
711  constexpr double rv = static_cast<double>(R_v);
712  constexpr double t0c = 273.15;
713  [[maybe_unused]] constexpr double ep1 = static_cast<double>(R_v / R_d - one);
714  constexpr double ep2 = static_cast<double>(R_d / R_v);
715  constexpr double qmin = 1.0e-12;
716  constexpr double xls = static_cast<double>(lsub);
717  constexpr double xlv0 = static_cast<double>(lat_vap);
718  constexpr double xlf0 = static_cast<double>(lat_ice);
719  constexpr double den0 = 1.28;
720  constexpr double denr = static_cast<double>(rhoh2o);
721  constexpr double dens = static_cast<double>(rhos);
722  constexpr double cliq = static_cast<double>(Cp_l);
723  constexpr double cice = 2106.0;
724  constexpr double psat = 610.78;
725  // Only the Fortran bridge consumes these; the native C++ path below derives
726  // its own. Mirrors ERF_AdvanceWSM6.cpp.
727  amrex::ignore_unused(g, rd, ep1);
728  [[maybe_unused]] const double ccn0 = static_cast<double>(m_ccn0);
729 
730  for (MFIter mfi(*mic_fab_vars[MicVar_WDM6::qv], TileNoZ()); mfi.isValid(); ++mfi) {
731  const Box box = mfi.tilebox();
732  const Box fab_box = mfi.fabbox();
733 
734  // Get array pointers for all WDM6 variables
735  auto const& t_arr = mic_fab_vars[MicVar_WDM6::tabs]->array(mfi);
736  auto const& qv_arr = mic_fab_vars[MicVar_WDM6::qv]->array(mfi);
737  auto const& qc_arr = mic_fab_vars[MicVar_WDM6::qc]->array(mfi);
738  auto const& qi_arr = mic_fab_vars[MicVar_WDM6::qi]->array(mfi);
739  auto const& qr_arr = mic_fab_vars[MicVar_WDM6::qr]->array(mfi);
740  auto const& qs_arr = mic_fab_vars[MicVar_WDM6::qs]->array(mfi);
741  auto const& qg_arr = mic_fab_vars[MicVar_WDM6::qg]->array(mfi);
742  auto const& nn_arr = mic_fab_vars[MicVar_WDM6::nn]->array(mfi); // Aerosol number
743  auto const& nc_arr = mic_fab_vars[MicVar_WDM6::nc]->array(mfi); // Cloud droplet number
744  auto const& nr_arr = mic_fab_vars[MicVar_WDM6::nr]->array(mfi); // Rain drop number
745  auto const& den_arr = mic_fab_vars[MicVar_WDM6::rho]->array(mfi);
746  auto const& p_arr = mic_fab_vars[MicVar_WDM6::pres]->array(mfi);
747  auto rain_arr = mic_fab_vars[MicVar_WDM6::rain_accum]->array(mfi);
748  auto snow_arr = mic_fab_vars[MicVar_WDM6::snow_accum]->array(mfi);
749  auto graup_arr = mic_fab_vars[MicVar_WDM6::graup_accum]->array(mfi);
750 
751  const int ilo = box.smallEnd(0);
752  const int ihi = box.bigEnd(0);
753  const int jlo = box.smallEnd(1);
754  const int jhi = box.bigEnd(1);
755  const int klo = box.smallEnd(2);
756  const int khi = box.bigEnd(2);
757 
758  [[maybe_unused]] const int imlo = fab_box.smallEnd(0);
759  [[maybe_unused]] const int imhi = fab_box.bigEnd(0);
760  [[maybe_unused]] const int jmlo = fab_box.smallEnd(1);
761  [[maybe_unused]] const int jmhi = fab_box.bigEnd(1);
762  [[maybe_unused]] const int kmlo = fab_box.smallEnd(2);
763  [[maybe_unused]] const int kmhi = fab_box.bigEnd(2);
764  const bool has_target_override = (micro_diag_target_column.size() == 2);
765  const int diag_i = has_target_override ? micro_diag_target_column[0] : ilo;
766  const int diag_j = has_target_override ? micro_diag_target_column[1] : jlo;
767 
768  auto const& w1_theta = mic_fab_vars[MicVar_WDM6::theta]->array(mfi);
769 
770 #if defined(ERF_USE_WDM6_FORT) && defined(AMREX_USE_GPU)
771  Arena* Arena_Used = run_wdm6_fort ? The_Pinned_Arena() : The_Async_Arena();
772 #else
773  Arena* Arena_Used = The_Async_Arena();
774 #endif
775 
776 #ifdef ERF_USE_WDM6_FORT
777  if (run_wdm6_fort) {
778  // Fortran bridge path
779  // Create delz array (cell thickness). Uniform dz over the storage box,
780  // then the physical thickness over the valid tile where a terrain /
781  // stretched-mesh nodal height field exists. Mirrors the WSM6 delz_arr
782  // fill (ERF_AdvanceWSM6.cpp:948-961) including the four-corner average.
783  // The ghost entries keep dz_val because the isohelper repacks delz only
784  // over its:ite / kts:kte (ERF_module_mp_wdm6_isohelper.F90:140), so no
785  // ghost value ever reaches the Fortran.
786  const Real dz_val = m_geom.CellSize(2);
787  FArrayBox delz_fab(fab_box, 1, Arena_Used);
788  auto const& delz_arr = delz_fab.array();
789  ParallelFor(fab_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
790  delz_arr(i,j,k) = dz_val;
791  });
792 
793  const Array4<const Real> z_arr = (m_z_phys_nd) ? m_z_phys_nd->const_array(mfi) : Array4<const Real> {};
794  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
795  delz_arr(i,j,k) = (z_arr) ? Real(0.25) * ( (z_arr(i ,j ,k+1) - z_arr(i ,j ,k))
796  + (z_arr(i+1,j ,k+1) - z_arr(i+1,j ,k))
797  + (z_arr(i ,j+1,k+1) - z_arr(i ,j+1,k))
798  + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz_val;
799  });
800 
801  // Surface slabs are anchored at klo, NOT at k = 0. The 2D buffers below
802  // are addressed with the ParallelFor's own k, while the 3D accumulation
803  // fields (rain_arr and friends, allocated on cons_in.boxArray()) are
804  // addressed at klo. Collapsing the slab to 0 makes those two agree only
805  // when klo happens to be 0, which is every case ERF can currently run
806  // because ERF_InitCustomPert_Bubble.H asserts each box spans the whole
807  // column. Anchoring at klo makes the correspondence hold by
808  // construction, and is bit-identical wherever klo == 0.
809  Box box2d(box);
810  box2d.makeSlab(2, klo);
811  Box fab_box2d(fab_box);
812  fab_box2d.makeSlab(2, klo);
813 
814  // Create landmask array in the WRF xland encoding: 1 = land, 2 = water.
815  // xland is handed to wdm62D's slmsk dummy unconverted (see the comment
816  // in mp_wdm6_run), and its only consumer is the maritime/continental
817  // autoconversion threshold at ERF_module_mp_wdm6.F90:628-633, which
818  // tests `slmsk .eq. 2`.
819  //
820  // ERF's lmask uses a DIFFERENT encoding -- 0 = water, 1 = land,
821  // 2 = building (ERF_MakeNewArrays.cpp:488 and :804) -- so the two
822  // cannot be passed through unmapped. Only lmask == 0 is water; a
823  // building cell is land as far as CCN is concerned.
824  FArrayBox xland_fab(fab_box2d, 1, Arena_Used);
825  auto const& xland_arr = xland_fab.array();
826  if (m_lmask != nullptr) {
827  auto const& lmask_arr = m_lmask->const_array(mfi);
828  ParallelFor(fab_box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
829  xland_arr(i,j,k) = (lmask_arr(i,j,0) == 0) ? Real(2.0) : Real(1.0);
830  });
831  } else {
832  // No mask wired (e.g. a restart path that has not set it): land.
833  ParallelFor(fab_box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
834  xland_arr(i,j,k) = Real(1.0);
835  });
836  }
837 
838  // Create 2D accumulation arrays
839  // Fortran bridge uses ims:ime, jms:jme storage bounds; these buffers must
840  // therefore be allocated on fab_box extents even if C++ kernels only
841  // update the valid tile slab (box2d).
842  FArrayBox rainacc_fab(fab_box2d, 1, Arena_Used);
843  FArrayBox rainncv_fab(fab_box2d, 1, Arena_Used);
844  FArrayBox sr_fab(fab_box2d, 1, Arena_Used);
845  FArrayBox snowacc_fab(fab_box2d, 1, Arena_Used);
846  FArrayBox snowncv_fab(fab_box2d, 1, Arena_Used);
847  FArrayBox graupacc_fab(fab_box2d, 1, Arena_Used);
848  FArrayBox graupelncv_fab(fab_box2d, 1, Arena_Used);
849 
850  auto const& rainacc_arr = rainacc_fab.array();
851  auto const& rainncv_arr = rainncv_fab.array();
852  auto const& sr_arr = sr_fab.array();
853  auto const& snowacc_arr = snowacc_fab.array();
854  auto const& snowncv_arr = snowncv_fab.array();
855  auto const& graupacc_arr = graupacc_fab.array();
856  auto const& graupelncv_arr = graupelncv_fab.array();
857 
858  // Initialize 2D arrays to zero over the full storage slab so the Fortran
859  // bridge never reads ghost entries that were never written.
860  ParallelFor(fab_box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
861  rainacc_arr(i,j,k) = Real(0.0);
862  rainncv_arr(i,j,k) = Real(0.0);
863  sr_arr(i,j,k) = Real(0.0);
864  snowacc_arr(i,j,k) = Real(0.0);
865  snowncv_arr(i,j,k) = Real(0.0);
866  graupacc_arr(i,j,k) = Real(0.0);
867  graupelncv_arr(i,j,k) = Real(0.0);
868  });
869 
870  // (Tile-based diagnostics removed - using global diagnostics instead)
871 
872  // Host-only Fortran reads mic_fab_vars and the buffers filled above via
873  // dataPtr(), so the GPU writes must complete before crossing the
874  // language boundary. Without this the Fortran can read delz, xland and
875  // the zeroed accumulators while their kernels are still in flight.
876  // Mirrors ERF_AdvanceWSM6.cpp.
877  Gpu::streamSynchronize();
878 
879  // Call Fortran WDM6
880  mp_wdm6_run_c(
881  t_arr.dataPtr(),
882  qv_arr.dataPtr(), qc_arr.dataPtr(), qi_arr.dataPtr(),
883  qr_arr.dataPtr(), qs_arr.dataPtr(), qg_arr.dataPtr(),
884  nn_arr.dataPtr(), nc_arr.dataPtr(), nr_arr.dataPtr(), // WDM6 number concentrations
885  den_arr.dataPtr(), p_arr.dataPtr(), delz_arr.dataPtr(),
886  static_cast<double>(dt_advance), g, cpd, cpv, rd, rv, t0c, ep1, ep2, qmin,
887  xls, xlv0, xlf0, den0, denr, cliq, cice, psat,
888  ccn0, xland_arr.dataPtr(), // WDM6 parameters
889  rainacc_arr.dataPtr(), rainncv_arr.dataPtr(), sr_arr.dataPtr(),
890  snowacc_arr.dataPtr(), snowncv_arr.dataPtr(),
891  graupacc_arr.dataPtr(), graupelncv_arr.dataPtr(),
892  imlo, imhi, jmlo, jmhi, kmlo, kmhi,
893  ilo, ihi, jlo, jhi, klo, khi,
894  microphysics_debug, diag_i, diag_j);
895 
896  // CRITICAL: Convert updated temperature back to potential temperature
897  // The Fortran WDM6 modifies t_arr (absolute temperature) due to latent heating/cooling
898  // from condensation, evaporation, freezing, melting, etc.
899  // ERF stores theta (potential temperature), so we must convert back: theta = T / exner
900  // This matches WRF's conversion: th(i,k,j) = t(i,k) / pii(i,k,j)
901  auto const& theta_arr = mic_fab_vars[MicVar_WDM6::theta]->array(mfi);
902  const Real configured_rdOcp = m_rdOcp;
903  const bool use_anelastic_reference_pressure =
904  m_use_anelastic_reference_pressure;
905  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
906  // Use the same held pressure that was passed to WDM6. In
907  // anelastic mode this is the BaseState p0; in compressible mode
908  // it is the pressure diagnosed from the incoming state.
909  theta_arr(i,j,k) = wdm6_theta_from_temperature_and_pressure(
910  t_arr(i,j,k), p_arr(i,j,k), configured_rdOcp,
911  use_anelastic_reference_pressure);
912  });
913 
914  // (Tile-based precipitation diagnostics removed - using global diagnostics instead)
915 
916  // Accumulate precipitation
917  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
918  rain_arr(i,j,k) += rainacc_arr(i,j,k);
919  snow_arr(i,j,k) += snowacc_arr(i,j,k);
920  graup_arr(i,j,k) += graupacc_arr(i,j,k);
921  });
922  } else {
923 #endif
924  // ===================================================================
925  // WDM6 C++ GPU kernel path (adapted from WSM6 with double-moment)
926  // ===================================================================
927 
928  // Working FABs (similar to WSM6 but with WDM6-specific additions)
929  // 3D working arrays
930  const Real dz_val = m_geom.CellSize(2);
931  // Nodal physical heights, when a terrain / stretched mesh supplies them.
932  // Consumed in the delz_arr fill below; kept identical to the bridge
933  // leg's four-corner average so the two legs see the same thickness.
934  const Array4<const Real> z_arr = (m_z_phys_nd) ? m_z_phys_nd->const_array(mfi) : Array4<const Real> {};
935  FArrayBox delz_fab(fab_box,1, Arena_Used);
936  FArrayBox denfac_fab(fab_box,1, Arena_Used);
937  FArrayBox xni_fab(fab_box,1, Arena_Used);
938  FArrayBox rslopec_fab(fab_box,1, Arena_Used);
939  FArrayBox rslopec2_fab(fab_box,1, Arena_Used);
940  FArrayBox rslopec3_fab(fab_box,1, Arena_Used);
941  FArrayBox rslope_fab(fab_box,3, Arena_Used);
942  FArrayBox rslopeb_fab(fab_box,3, Arena_Used);
943  FArrayBox rslope2_fab(fab_box,3, Arena_Used);
944  FArrayBox rslope3_fab(fab_box,3, Arena_Used);
945  FArrayBox work1_fab(fab_box,3, Arena_Used);
946  FArrayBox workn_fab(fab_box,1, Arena_Used);
947  FArrayBox work2_fab(fab_box,1, Arena_Used); // Ventilation factor for diffusion (G11+)
948  // Anchored at klo, not 0. mstep/numdt/sr/fallc/delqi all live on this
949  // slab and are written with the ParallelFor's own k, while the G9
950  // accumulation block reads them back at klo alongside the 3D surface
951  // fields (rain_arr, delz_arr, work1_arr). With the slab at 0 those two
952  // conventions only coincide when klo == 0; sr_arr(i,j,klo) and
953  // fallc_arr(i,j,klo) were otherwise out of bounds on this FAB.
954  Box box2d(IntVect(ilo,jlo,klo), IntVect(ihi,jhi,klo));
955  IArrayBox mstep_fab(box2d,1, Arena_Used);
956  IArrayBox numdt_fab(box2d,1, Arena_Used);
957  FArrayBox sr_fab(box2d, 1, Arena_Used); // Snow ratio for G9 output
958  FArrayBox cpm_fab(fab_box,1, Arena_Used);
959  FArrayBox xl_fab(fab_box,1, Arena_Used);
960  FArrayBox qsatw_fab(fab_box,1, Arena_Used);
961  FArrayBox qsati_fab(fab_box,1, Arena_Used);
962  FArrayBox rhw_fab(fab_box,1, Arena_Used);
963  FArrayBox rhi_fab(fab_box,1, Arena_Used);
964  FArrayBox qcr_fab(fab_box,1, Arena_Used);
965  FArrayBox sed_cell_scratch_fab(fab_box, WDM6SedCellScratch::NumComps, Arena_Used);
966  Box sed_node_box = amrex::surroundingNodes(fab_box, 2);
967  FArrayBox sed_node_scratch_fab(sed_node_box, WDM6SedNodeScratch::NumComps, Arena_Used);
968 
969  // Process rate arrays
970  FArrayBox praut_fab(fab_box,1, Arena_Used);
971  FArrayBox pracw_fab(fab_box,1, Arena_Used);
972  FArrayBox prevp_fab(fab_box,1, Arena_Used);
973  FArrayBox pidep_fab(fab_box,1, Arena_Used);
974  FArrayBox psdep_fab(fab_box,1, Arena_Used);
975  FArrayBox pgdep_fab(fab_box,1, Arena_Used);
976  FArrayBox pigen_fab(fab_box,1, Arena_Used);
977  FArrayBox psaut_fab(fab_box,1, Arena_Used);
978  FArrayBox pgaut_fab(fab_box,1, Arena_Used);
979  FArrayBox pcact_fab(fab_box,1, Arena_Used);
980  FArrayBox pcond_fab(fab_box,1, Arena_Used);
981  FArrayBox praci_fab(fab_box,1, Arena_Used);
982  FArrayBox piacr_fab(fab_box,1, Arena_Used);
983  FArrayBox niacr_fab(fab_box,1, Arena_Used);
984  FArrayBox psaci_fab(fab_box,1, Arena_Used);
985  FArrayBox pgaci_fab(fab_box,1, Arena_Used);
986  FArrayBox psacw_fab(fab_box,1, Arena_Used);
987  FArrayBox nsacw_fab(fab_box,1, Arena_Used);
988  FArrayBox pgacw_fab(fab_box,1, Arena_Used);
989  FArrayBox ngacw_fab(fab_box,1, Arena_Used);
990  FArrayBox paacw_fab(fab_box,1, Arena_Used);
991  FArrayBox naacw_fab(fab_box,1, Arena_Used);
992  FArrayBox pracs_fab(fab_box,1, Arena_Used);
993  FArrayBox psacr_fab(fab_box,1, Arena_Used);
994  FArrayBox nsacr_fab(fab_box,1, Arena_Used);
995  FArrayBox pgacr_fab(fab_box,1, Arena_Used);
996  FArrayBox ngacr_fab(fab_box,1, Arena_Used);
997  FArrayBox pgacs_fab(fab_box,1, Arena_Used);
998  FArrayBox pseml_fab(fab_box,1, Arena_Used);
999  FArrayBox nseml_fab(fab_box,1, Arena_Used);
1000  FArrayBox pgeml_fab(fab_box,1, Arena_Used);
1001  FArrayBox ngeml_fab(fab_box,1, Arena_Used);
1002  FArrayBox psevp_fab(fab_box,1, Arena_Used);
1003  FArrayBox pgevp_fab(fab_box,1, Arena_Used);
1004 
1005  // Number concentration process rates (WDM6-specific)
1006  FArrayBox ncauto_fab(fab_box,1, Arena_Used); // nc lost to autoconversion
1007  FArrayBox ncaccr_fab(fab_box,1, Arena_Used); // nc lost to accretion by rain
1008  FArrayBox nrauto_fab(fab_box,1, Arena_Used); // nr gained from autoconversion
1009  FArrayBox nraccr_fab(fab_box,1, Arena_Used); // nr gained from accretion
1010  FArrayBox nrevp_fab(fab_box,1, Arena_Used); // nr lost to evaporation
1011  FArrayBox ncact_fab(fab_box,1, Arena_Used); // nn -> nc activation rate
1012  FArrayBox act_ratio_fab(fab_box,1, Arena_Used);
1013  FArrayBox act_fraction_fab(fab_box,1, Arena_Used);
1014  FArrayBox act_raw_fab(fab_box,1, Arena_Used);
1015  FArrayBox act_cap_fab(fab_box,1, Arena_Used);
1016  FArrayBox nccol_fab(fab_box,1, Arena_Used); // cloud self-collection number sink
1017  FArrayBox nrcol_fab(fab_box,1, Arena_Used); // rain self-collection number sink
1018 
1019  // G6 temporary arrays for slope_wdm6 call
1020  FArrayBox qrs_tmp_fab(fab_box,3, Arena_Used); // Temporary qr, qs, qg (components 0,1,2)
1021  FArrayBox ncr_tmp_fab(fab_box,1, Arena_Used); // Temporary nr
1022 
1023  // G11 arrays (particle diameter work array)
1024  FArrayBox avedia_fab(fab_box,2, Arena_Used); // avedia(:,:,1:2) for rain and cloud slopes
1025 
1026  // G8 ice sedimentation arrays
1027  FArrayBox work1c_fab(fab_box,1, Arena_Used); // Ice crystal fall speed
1028  FArrayBox fallc_fab(box2d,1, Arena_Used); // Ice fallout (2D surface)
1029  FArrayBox delqi_fab(box2d,1, Arena_Used); // Ice precipitation at bottom (2D)
1030 
1031  auto const& delz_arr = delz_fab.array();
1032  auto const& denfac_arr = denfac_fab.array();
1033  auto const& xni_arr = xni_fab.array();
1034  auto const& rslopec_arr = rslopec_fab.array();
1035  auto const& rslopec2_arr = rslopec2_fab.array();
1036  auto const& rslopec3_arr = rslopec3_fab.array();
1037  auto const& rslope_arr = rslope_fab.array();
1038  auto const& rslopeb_arr = rslopeb_fab.array();
1039  auto const& rslope2_arr = rslope2_fab.array();
1040  auto const& rslope3_arr = rslope3_fab.array();
1041  auto const& work1_arr = work1_fab.array();
1042  auto const& workn_arr = workn_fab.array();
1043  auto const& work2_arr = work2_fab.array(); // Ventilation factor for diffusion
1044  auto const& mstep_arr = mstep_fab.array();
1045  auto const& numdt_arr = numdt_fab.array();
1046  auto const& sr_arr = sr_fab.array(); // Snow ratio for G9
1047  auto const& cpm_arr = cpm_fab.array();
1048  auto const& xl_arr = xl_fab.array();
1049  auto const& qsatw_arr = qsatw_fab.array();
1050  auto const& qsati_arr = qsati_fab.array();
1051  auto const& rhw_arr = rhw_fab.array();
1052  auto const& rhi_arr = rhi_fab.array();
1053  auto const& qcr_arr = qcr_fab.array();
1054  auto const& sed_cell_scratch_arr = sed_cell_scratch_fab.array();
1055  auto const& sed_node_scratch_arr = sed_node_scratch_fab.array();
1056  auto const& praut_arr = praut_fab.array();
1057  auto const& pracw_arr = pracw_fab.array();
1058  auto const& prevp_arr = prevp_fab.array();
1059  auto const& pidep_arr = pidep_fab.array();
1060  auto const& psdep_arr = psdep_fab.array();
1061  auto const& pgdep_arr = pgdep_fab.array();
1062  auto const& pigen_arr = pigen_fab.array();
1063  auto const& psaut_arr = psaut_fab.array();
1064  auto const& pgaut_arr = pgaut_fab.array();
1065  auto const& pcact_arr = pcact_fab.array();
1066  auto const& pcond_arr = pcond_fab.array();
1067  auto const& praci_arr = praci_fab.array();
1068  auto const& piacr_arr = piacr_fab.array();
1069  auto const& niacr_arr = niacr_fab.array();
1070  auto const& psaci_arr = psaci_fab.array();
1071  auto const& pgaci_arr = pgaci_fab.array();
1072  auto const& psacw_arr = psacw_fab.array();
1073  auto const& nsacw_arr = nsacw_fab.array();
1074  auto const& pgacw_arr = pgacw_fab.array();
1075  auto const& ngacw_arr = ngacw_fab.array();
1076  auto const& paacw_arr = paacw_fab.array();
1077  auto const& naacw_arr = naacw_fab.array();
1078  auto const& pracs_arr = pracs_fab.array();
1079  auto const& psacr_arr = psacr_fab.array();
1080  auto const& nsacr_arr = nsacr_fab.array();
1081  auto const& pgacr_arr = pgacr_fab.array();
1082  auto const& ngacr_arr = ngacr_fab.array();
1083  auto const& pgacs_arr = pgacs_fab.array();
1084  auto const& pseml_arr = pseml_fab.array();
1085  auto const& nseml_arr = nseml_fab.array();
1086  auto const& pgeml_arr = pgeml_fab.array();
1087  auto const& ngeml_arr = ngeml_fab.array();
1088  auto const& psevp_arr = psevp_fab.array();
1089  auto const& pgevp_arr = pgevp_fab.array();
1090  auto const& ncauto_arr = ncauto_fab.array();
1091  auto const& ncaccr_arr = ncaccr_fab.array();
1092  auto const& nrauto_arr = nrauto_fab.array();
1093  auto const& nraccr_arr = nraccr_fab.array();
1094  auto const& nrevp_arr = nrevp_fab.array();
1095  auto const& ncact_arr = ncact_fab.array();
1096  auto const& act_ratio_arr = act_ratio_fab.array();
1097  auto const& act_fraction_arr = act_fraction_fab.array();
1098  auto const& act_raw_arr = act_raw_fab.array();
1099  auto const& act_cap_arr = act_cap_fab.array();
1100  auto const& nccol_arr = nccol_fab.array();
1101  auto const& nrcol_arr = nrcol_fab.array();
1102  auto const& qrs_tmp_arr = qrs_tmp_fab.array(); // G6: temporary qr, qs, qg
1103  auto const& ncr_tmp_arr = ncr_tmp_fab.array(); // G6: temporary nr
1104  auto const& avedia_arr = avedia_fab.array(); // G11: particle diameter work array
1105 
1106  // G8 array references
1107  auto const& work1c_arr = work1c_fab.array(); // Ice crystal fall speed
1108  auto const& fallc_arr = fallc_fab.array(); // Ice fallout (2D)
1109  auto const& delqi_arr = delqi_fab.array(); // Ice precip at bottom (2D)
1110 
1111  // Clamp negative values and enforce minimums
1112  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1113  delz_arr(i,j,k) = (z_arr) ? Real(0.25) * ( (z_arr(i ,j ,k+1) - z_arr(i ,j ,k))
1114  + (z_arr(i+1,j ,k+1) - z_arr(i+1,j ,k))
1115  + (z_arr(i ,j+1,k+1) - z_arr(i ,j+1,k))
1116  + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz_val;
1117  qc_arr(i,j,k) = amrex::max(qc_arr(i,j,k), Real(0.0));
1118  qr_arr(i,j,k) = amrex::max(qr_arr(i,j,k), Real(0.0));
1119  qi_arr(i,j,k) = amrex::max(qi_arr(i,j,k), Real(0.0));
1120  qs_arr(i,j,k) = amrex::max(qs_arr(i,j,k), Real(0.0));
1121  qg_arr(i,j,k) = amrex::max(qg_arr(i,j,k), Real(0.0));
1122 
1123  // Match Fortran pre-G3 behavior: nc is non-negative here, but not floored to ncmin yet.
1124  nc_arr(i,j,k) = amrex::max(nc_arr(i,j,k), Real(0.0));
1125  nr_arr(i,j,k) = amrex::max(nr_arr(i,j,k), Real(0.0));
1126  // CCN is clamped into [1.e8, 2.e10], not merely floored at zero:
1127  // ncr(i,k,1) = min(max(ncr(i,k,1),1.e8),2.e10)
1128  // at ERF_module_mp_wdm6.F90:604. This re-seeds the aerosol field to
1129  // the background concentration at the top of every macro timestep,
1130  // so activation depletion never carries across steps. A plain
1131  // max(.,0) is a no-op on step 1, where Init() has just set nn to
1132  // ccn0 = 1.e8 everywhere, and only diverges from step 2 onward once
1133  // activation has drawn nn below the floor. Both literals are exact
1134  // in float32, so the precision contract is satisfied either way.
1135  nn_arr(i,j,k) = amrex::min(amrex::max(nn_arr(i,j,k), wdm6_literal(1.0e8)),
1136  wdm6_literal(2.0e10));
1137  });
1138 
1139  // Compute cpm and xl once from initial state
1140  const Real xlv1_loc = m_xlv1;
1141  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1142  cpm_arr(i,j,k) = wdm6_cpmcal(qv_arr(i,j,k), Real(qmin), Real(cpd), Real(cpv));
1143  xl_arr(i,j,k) = wdm6_xlcal(t_arr(i,j,k), Real(xlv0), xlv1_loc, Real(t0c));
1144  });
1145 
1146  // Minor timestep loop (match WSM6 structure)
1147  const int wdm6_loops = std::max(
1148  static_cast<int>(std::round(dt_advance / Real(120.0))), 1); // dtcldcr = 120s
1149  const Real dtcld = dt_advance / static_cast<Real>(wdm6_loops);
1150 
1151  // Extract WDM6 coefficients
1152  const Real qc0_loc = m_qc0;
1153  const Real qc1_loc = m_qc1;
1154  const Real qck1_loc = m_qck1;
1155  const Real pidnc_loc = m_pidnc;
1156  const Real pidnr_loc = m_pidnr;
1157  const Real pvtr_loc = m_pvtr;
1158  const Real pvtrn_loc = m_pvtrn;
1159  const Real pvts_loc = m_pvts;
1160  const Real pvtg_loc = m_pvtg;
1161  const Real slope_bvtg_loc = m_bvtg;
1162  const Real pidn0s_loc = m_pidn0s;
1163  const Real pidn0g_loc = m_pidn0g;
1164  const Real rslopermax_loc = m_rslopermax;
1165  const Real rsloperbmax_loc = m_rsloperbmax;
1166  const Real rsloper2max_loc = m_rsloper2max;
1167  const Real rsloper3max_loc = m_rsloper3max;
1168  const Real rslopesmax_loc = m_rslopesmax;
1169  const Real rslopesbmax_loc = m_rslopesbmax;
1170  const Real rslopes2max_loc = m_rslopes2max;
1171  const Real rslopes3max_loc = m_rslopes3max;
1172  const Real rslopegmax_loc = m_rslopegmax;
1173  const Real rslopegbmax_loc = m_rslopegbmax;
1174  const Real rslopeg2max_loc = m_rslopeg2max;
1175  const Real rslopeg3max_loc = m_rslopeg3max;
1176  const Real rslopecmax_loc = m_rslopecmax;
1177  const Real rslopec2max_loc = m_rslopec2max;
1178  const Real rslopec3max_loc = m_rslopec3max;
1179  const Real precs1_loc = m_precs1;
1180  const Real precs2_loc = m_precs2;
1181  const Real precg1_loc = m_precg1;
1182  const Real precg2_loc = m_precg2;
1183  const Real precr1_loc = m_precr1;
1184  const Real precr2_loc = m_precr2;
1185  const Real n0g_loc = m_n0g;
1186  const Real pi_wdm6_loc = m_pi_wdm6;
1187  const Real pacrc_loc = m_pacrc;
1188  const Real pacrg_loc = m_pacrg;
1189  const Real g4pbr_loc = m_g4pbr;
1190  const Real g7pbr_loc = m_g7pbr;
1191  constexpr Real pfrz1_loc = pfrz1;
1192  constexpr Real pfrz2_loc = pfrz2;
1193  const int diag_k = klo;
1194 
1195  // Maritime/continental autoconversion threshold. The Fortran selects on
1196  // `slmsk .eq. 2` (ERF_module_mp_wdm6.F90:628-633) where slmsk is xland
1197  // in the WRF encoding, so 2 means WATER and picks qc0. ERF's lmask
1198  // encodes water as 0, not 2 (ERF_MakeNewArrays.cpp:488, :804), so the
1199  // test must be against 0 here. Testing == 2 would select maritime for
1200  // building cells and would disagree with the bridge leg, which maps
1201  // lmask == 0 to xland = 2.0.
1202  if (m_lmask != nullptr) {
1203  auto const& lmask_arr = m_lmask->const_array(mfi);
1204  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1205  qcr_arr(i,j,k) = (lmask_arr(i,j,0) == 0) ? qc0_loc : qc1_loc;
1206  });
1207  } else {
1208  // Match the bridge leg's no-mask posture: default to land.
1209  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1210  qcr_arr(i,j,k) = qc1_loc;
1211  });
1212  }
1213 
1214  for (int loop = 0; loop < wdm6_loops; ++loop) {
1215  // ============================================================
1216  // Step 1: Density factor (G1b / DENFAC)
1217  // ============================================================
1218 
1219  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1220  denfac_arr(i,j,k) = std::sqrt(Real(den0) / den_arr(i,j,k));
1221  });
1222 
1223 
1224  // ============================================================
1225  // Step 2: Saturation calculations (G1c / QSAT)
1226  // ============================================================
1227  {
1228  const Real ttp = Real(t0c) + wdm6_literal(0.01);
1229  const Real dldt = Real(cpv) - Real(cliq);
1230  const Real xa = -dldt / Real(rv);
1231  const Real xb = xa + Real(xlv0) / (Real(rv) * ttp);
1232  const Real dldti = Real(cpv) - Real(cice);
1233  const Real xai = -dldti / Real(rv);
1234  const Real xbi = xai + Real(xls) / (Real(rv) * ttp);
1235 
1236 
1237  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1238  const Real tr = ttp / t_arr(i,j,k);
1239 
1240  // Saturation over water
1241  Real qsw = Real(psat) * std::exp(std::log(tr) * xa) * std::exp(xb * (Real(1.0) - tr));
1242  qsw = amrex::min(qsw, wdm6_literal(0.99) * p_arr(i,j,k));
1243  qsw = Real(ep2) * qsw / (p_arr(i,j,k) - qsw);
1244  qsw = amrex::max(qsw, Real(qmin));
1245  qsatw_arr(i,j,k) = qsw;
1246  rhw_arr(i,j,k) = amrex::max(qv_arr(i,j,k) / qsw, Real(qmin));
1247 
1248  // Saturation over ice
1249  Real qsi = (t_arr(i,j,k) < ttp)
1250  ? Real(psat) * std::exp(std::log(tr) * xai) * std::exp(xbi * (Real(1.0) - tr))
1251  : Real(psat) * std::exp(std::log(tr) * xa) * std::exp(xb * (Real(1.0) - tr));
1252  qsi = amrex::min(qsi, wdm6_literal(0.99) * p_arr(i,j,k));
1253  qsi = Real(ep2) * qsi / (p_arr(i,j,k) - qsi);
1254  qsi = amrex::max(qsi, Real(qmin));
1255  qsati_arr(i,j,k) = qsi;
1256  rhi_arr(i,j,k) = amrex::max(qv_arr(i,j,k) / qsi, Real(qmin));
1257  });
1258  }
1259 
1260 
1261  // ============================================================
1262  // Step 3: Zero process rates (G2 / RATES_ZERO)
1263  // ============================================================
1264 
1265  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1266  praut_arr(i,j,k) = Real(0.0);
1267  pracw_arr(i,j,k) = Real(0.0);
1268  prevp_arr(i,j,k) = Real(0.0);
1269  pidep_arr(i,j,k) = Real(0.0);
1270  psdep_arr(i,j,k) = Real(0.0);
1271  pgdep_arr(i,j,k) = Real(0.0);
1272  pigen_arr(i,j,k) = Real(0.0);
1273  psaut_arr(i,j,k) = Real(0.0);
1274  pgaut_arr(i,j,k) = Real(0.0);
1275  pcact_arr(i,j,k) = Real(0.0);
1276  pcond_arr(i,j,k) = Real(0.0);
1277  praci_arr(i,j,k) = Real(0.0);
1278  piacr_arr(i,j,k) = Real(0.0);
1279  niacr_arr(i,j,k) = Real(0.0);
1280  psaci_arr(i,j,k) = Real(0.0);
1281  pgaci_arr(i,j,k) = Real(0.0);
1282  psacw_arr(i,j,k) = Real(0.0);
1283  nsacw_arr(i,j,k) = Real(0.0);
1284  pgacw_arr(i,j,k) = Real(0.0);
1285  ngacw_arr(i,j,k) = Real(0.0);
1286  paacw_arr(i,j,k) = Real(0.0);
1287  naacw_arr(i,j,k) = Real(0.0);
1288  pracs_arr(i,j,k) = Real(0.0);
1289  psacr_arr(i,j,k) = Real(0.0);
1290  nsacr_arr(i,j,k) = Real(0.0);
1291  pgacr_arr(i,j,k) = Real(0.0);
1292  ngacr_arr(i,j,k) = Real(0.0);
1293  pgacs_arr(i,j,k) = Real(0.0);
1294  pseml_arr(i,j,k) = Real(0.0);
1295  nseml_arr(i,j,k) = Real(0.0);
1296  pgeml_arr(i,j,k) = Real(0.0);
1297  ngeml_arr(i,j,k) = Real(0.0);
1298  psevp_arr(i,j,k) = Real(0.0);
1299  pgevp_arr(i,j,k) = Real(0.0);
1300  ncauto_arr(i,j,k) = Real(0.0);
1301  ncaccr_arr(i,j,k) = Real(0.0);
1302  nrauto_arr(i,j,k) = Real(0.0);
1303  nraccr_arr(i,j,k) = Real(0.0);
1304  nrevp_arr(i,j,k) = Real(0.0);
1305  ncact_arr(i,j,k) = Real(0.0);
1306  act_ratio_arr(i,j,k) = Real(0.0);
1307  act_fraction_arr(i,j,k) = Real(0.0);
1308  act_raw_arr(i,j,k) = Real(0.0);
1309  act_cap_arr(i,j,k) = Real(0.0);
1310  nccol_arr(i,j,k) = Real(0.0);
1311  nrcol_arr(i,j,k) = Real(0.0);
1312  });
1313 
1314 
1315  // ============================================================
1316  // Step 3b: CLOUD_SETUP (G3)
1317  // Exact port of the frozen Fortran block:
1318  // - cloud droplet slope parameter rslopec{,2,3}
1319  // - ice number concentration xni
1320  // ============================================================
1321  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1322  if (qc_arr(i,j,k) <= Real(qmin) || nc_arr(i,j,k) <= Real(1.e1)) {
1323  rslopec_arr(i,j,k) = rslopecmax_loc;
1324  rslopec2_arr(i,j,k) = rslopec2max_loc;
1325  rslopec3_arr(i,j,k) = rslopec3max_loc;
1326  } else {
1327  rslopec_arr(i,j,k) = wdm6_rslopec_exact(qc_arr(i,j,k), den_arr(i,j,k),
1328  nc_arr(i,j,k), pidnc_loc);
1329  rslopec2_arr(i,j,k) = rslopec_arr(i,j,k) * rslopec_arr(i,j,k);
1330  rslopec3_arr(i,j,k) = rslopec2_arr(i,j,k) * rslopec_arr(i,j,k);
1331  }
1332  xni_arr(i,j,k) = wdm6_xni_exact(qi_arr(i,j,k), den_arr(i,j,k), Real(qmin));
1333  });
1334 
1335 
1336  // ============================================================
1337  // Step 3c: SLOPE1 (G4)
1338  // Exact first packed slope_wdm6 surface for rain/snow/graupel.
1339  // ============================================================
1340  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1341  Real rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3, rain_vt, rain_vtn;
1342  Real snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt, snow_n0sfac;
1343  Real graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt;
1344 
1345  wdm6_slope_rain_cell(qr_arr(i,j,k), nr_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1346  Real(qcrmin), Real(nrmin),
1347  rslopermax_loc, rsloperbmax_loc, rsloper2max_loc, rsloper3max_loc,
1348  Real(bvtr), pvtr_loc, pvtrn_loc, pidnr_loc,
1349  rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3,
1350  rain_vt, rain_vtn);
1351  wdm6_slope_snow_cell(qs_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k), t_arr(i,j,k),
1352  pidn0s_loc, Real(alpha_wdm6), Real(n0smax), Real(n0s),
1353  wdm6_slope_t0c, Real(qcrmin),
1354  rslopesmax_loc, rslopesbmax_loc, rslopes2max_loc, rslopes3max_loc,
1355  Real(bvts), pvts_loc,
1356  snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt,
1357  snow_n0sfac);
1358  wdm6_slope_graup_cell(qg_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1359  pidn0g_loc, Real(qcrmin),
1360  rslopegmax_loc, rslopegbmax_loc, rslopeg2max_loc, rslopeg3max_loc,
1361  slope_bvtg_loc, pvtg_loc,
1362  graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt);
1363 
1364  rslope_arr(i,j,k,0) = rain_rslope;
1365  rslope_arr(i,j,k,1) = snow_rslope;
1366  rslope_arr(i,j,k,2) = graup_rslope;
1367  rslopeb_arr(i,j,k,0) = rain_rslopeb;
1368  rslopeb_arr(i,j,k,1) = snow_rslopeb;
1369  rslopeb_arr(i,j,k,2) = graup_rslopeb;
1370  rslope2_arr(i,j,k,0) = rain_rslope2;
1371  rslope2_arr(i,j,k,1) = snow_rslope2;
1372  rslope2_arr(i,j,k,2) = graup_rslope2;
1373  rslope3_arr(i,j,k,0) = rain_rslope3;
1374  rslope3_arr(i,j,k,1) = snow_rslope3;
1375  rslope3_arr(i,j,k,2) = graup_rslope3;
1376  work1_arr(i,j,k,0) = rain_vt;
1377  work1_arr(i,j,k,1) = snow_vt;
1378  work1_arr(i,j,k,2) = graup_vt;
1379  workn_arr(i,j,k) = rain_vtn;
1380  });
1381 
1382  // ============================================================
1383  // Step 3d: Rain sedimentation setup (G5a)
1384  // Determine rain/number substep count from normalized work1/workn.
1385  // ============================================================
1386  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1387  mstep_arr(i,j,k) = 1;
1388  numdt_arr(i,j,k) = 1;
1389  sr_arr(i,j,k) = Real(0.0); // Initialize snow ratio to zero
1390  });
1391  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1392  int mstep_loc = 1;
1393  int numdt_loc = 1;
1394  for (int kk = khi; kk >= klo; --kk) {
1395  work1_arr(i,j,kk,0) = work1_arr(i,j,kk,0) / delz_arr(i,j,kk);
1396  workn_arr(i,j,kk) = workn_arr(i,j,kk) / delz_arr(i,j,kk);
1397  numdt_loc = amrex::max(static_cast<int>(amrex::max(work1_arr(i,j,kk,0),
1398  workn_arr(i,j,kk)) * dtcld + Real(0.5)),
1399  1);
1400  if (numdt_loc >= mstep_loc) {
1401  mstep_loc = numdt_loc;
1402  }
1403  }
1404  mstep_arr(i,j,k) = mstep_loc;
1405  numdt_arr(i,j,k) = numdt_loc;
1406  });
1407  ReduceOps<ReduceOpMax> reduce_op;
1408  ReduceData<int> reduce_data(reduce_op);
1409  reduce_op.eval(box2d, reduce_data,
1410  [=] AMREX_GPU_DEVICE (int i, int j, int k) -> GpuTuple<int> {
1411  return {mstep_arr(i,j,k)};
1412  });
1413  int mstepmax = amrex::get<0>(reduce_data.value());
1414  amrex::ignore_unused(mstepmax);
1415 
1416  // ============================================================
1417  // Step 3e: Rain sedimentation substeps (G5b)
1418  // Match the bounded top-down rain/nr transport loop and
1419  // slope_rain refresh that immediately follows G5a in Fortran.
1420  // ============================================================
1421  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1422  amrex::ignore_unused(k);
1423  const int col_mstep = mstep_arr(i,j,klo);
1424  for (int n = 1; n <= mstepmax; ++n) {
1425  if (n > col_mstep) {
1426  continue;
1427  }
1428 
1429  const int kk_top = khi;
1430  const Real top_flux_qr = den_arr(i,j,kk_top) * qr_arr(i,j,kk_top)
1431  * work1_arr(i,j,kk_top,0) / static_cast<Real>(col_mstep);
1432  const Real top_flux_nr = nr_arr(i,j,kk_top)
1433  * workn_arr(i,j,kk_top) / static_cast<Real>(col_mstep);
1434 
1435  qr_arr(i,j,kk_top) = amrex::max(
1436  qr_arr(i,j,kk_top) - top_flux_qr * dtcld / den_arr(i,j,kk_top),
1437  Real(0.0));
1438  nr_arr(i,j,kk_top) = amrex::max(
1439  nr_arr(i,j,kk_top) - top_flux_nr * dtcld,
1440  Real(0.0));
1441 
1442  Real flux_qr_above = top_flux_qr;
1443  Real flux_nr_above = top_flux_nr;
1444  for (int kk = khi - 1; kk >= klo; --kk) {
1445  const Real flux_qr = den_arr(i,j,kk) * qr_arr(i,j,kk)
1446  * work1_arr(i,j,kk,0) / static_cast<Real>(col_mstep);
1447  const Real flux_nr = nr_arr(i,j,kk)
1448  * workn_arr(i,j,kk) / static_cast<Real>(col_mstep);
1449 
1450  const Real dqr_self = amrex::min(
1451  flux_qr * dtcld / den_arr(i,j,kk),
1452  qr_arr(i,j,kk));
1453  const Real dqr_from_above = amrex::min(
1454  flux_qr_above * delz_arr(i,j,kk+1) / delz_arr(i,j,kk)
1455  * dtcld / den_arr(i,j,kk),
1456  qr_arr(i,j,kk+1));
1457  const Real dnr_self = amrex::min(
1458  flux_nr * dtcld,
1459  nr_arr(i,j,kk));
1460  const Real dnr_from_above = amrex::min(
1461  flux_nr_above * delz_arr(i,j,kk+1) / delz_arr(i,j,kk)
1462  * dtcld,
1463  nr_arr(i,j,kk+1));
1464 
1465  qr_arr(i,j,kk) = amrex::max(
1466  qr_arr(i,j,kk) - dqr_self + dqr_from_above,
1467  Real(0.0));
1468  nr_arr(i,j,kk) = amrex::max(
1469  nr_arr(i,j,kk) - dnr_self + dnr_from_above,
1470  Real(0.0));
1471 
1472  flux_qr_above = flux_qr;
1473  flux_nr_above = flux_nr;
1474  }
1475 
1476  for (int kk = klo; kk <= khi; ++kk) {
1477  Real rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3;
1478  Real rain_vt, rain_vtn;
1479  wdm6_slope_rain_cell(
1480  qr_arr(i,j,kk), nr_arr(i,j,kk),
1481  den_arr(i,j,kk), denfac_arr(i,j,kk),
1482  Real(qcrmin), Real(nrmin),
1483  rslopermax_loc, rsloperbmax_loc,
1484  rsloper2max_loc, rsloper3max_loc,
1485  Real(bvtr), pvtr_loc, pvtrn_loc, pidnr_loc,
1486  rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3,
1487  rain_vt, rain_vtn);
1488  rslope_arr(i,j,kk,0) = rain_rslope;
1489  rslopeb_arr(i,j,kk,0) = rain_rslopeb;
1490  rslope2_arr(i,j,kk,0) = rain_rslope2;
1491  rslope3_arr(i,j,kk,0) = rain_rslope3;
1492  work1_arr(i,j,kk,0) = rain_vt / delz_arr(i,j,kk);
1493  workn_arr(i,j,kk) = rain_vtn / delz_arr(i,j,kk);
1494  }
1495  }
1496  });
1497 
1498  // ============================================================
1499  // Step 3f: Snow/graupel sedimentation (G5c)
1500  // Exact bounded NISLFV_SG slice plus lower-boundary slab fall.
1501  // ============================================================
1502 
1503  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1504  amrex::ignore_unused(k);
1505  const int km = khi - klo + 1;
1506  auto dz_col = [&](int kk) -> amrex::Real& {
1507  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::dz);
1508  };
1509  auto den_col = [&](int kk) -> amrex::Real& {
1510  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::den);
1511  };
1512  auto denfac_col = [&](int kk) -> amrex::Real& {
1513  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::denfac);
1514  };
1515  auto tk_col = [&](int kk) -> amrex::Real& {
1516  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::tk);
1517  };
1518  auto work_col = [&](int kk) -> amrex::Real& {
1519  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::work_col);
1520  };
1521  auto rq_col = [&](int kk) -> amrex::Real& {
1522  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::rq_col);
1523  };
1524  auto rq2_col = [&](int kk) -> amrex::Real& {
1525  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::rq2_col);
1526  };
1527  for (int kk = 0; kk < km; ++kk) {
1528  const int k3 = klo + kk;
1529  dz_col(kk) = delz_arr(i,j,k3);
1530  den_col(kk) = den_arr(i,j,k3);
1531  denfac_col(kk) = denfac_arr(i,j,k3);
1532  tk_col(kk) = t_arr(i,j,k3);
1533  rq_col(kk) = den_col(kk) * qs_arr(i,j,k3);
1534  rq2_col(kk) = den_col(kk) * qg_arr(i,j,k3);
1535  const Real qsum = amrex::max(qs_arr(i,j,k3) + qg_arr(i,j,k3), wdm6_literal(1.0e-15));
1536  work_col(kk) = (qsum > wdm6_literal(1.0e-15))
1537  ? (work1_arr(i,j,k3,1) * qs_arr(i,j,k3) + work1_arr(i,j,k3,2) * qg_arr(i,j,k3)) / qsum
1538  : Real(0.0);
1539  }
1540 
1541  Real delqrs2 = Real(0.0);
1542  Real delqrs3 = Real(0.0);
1543  wdm6_nislfv_rain_plm6_column(
1544  km, &delqrs2, &delqrs3, dtcld, 1,
1545  pidn0s_loc, pidn0g_loc, Real(qcrmin), Real(alpha_wdm6), Real(n0smax), Real(n0s), wdm6_slope_t0c,
1546  rslopesmax_loc, rslopesbmax_loc, rslopes2max_loc, rslopes3max_loc, Real(bvts), pvts_loc,
1547  rslopegmax_loc, rslopegbmax_loc, rslopeg2max_loc, rslopeg3max_loc, slope_bvtg_loc, pvtg_loc,
1548  sed_cell_scratch_arr, sed_node_scratch_arr, i, j, klo);
1549 
1550  for (int kk = 0; kk < km; ++kk) {
1551  const int k3 = klo + kk;
1552  qs_arr(i,j,k3) = amrex::max(rq_col(kk) / den_col(kk), Real(0.0));
1553  qg_arr(i,j,k3) = amrex::max(rq2_col(kk) / den_col(kk), Real(0.0));
1554  }
1555 
1556  work1_arr(i,j,klo,1) = delqrs2 / dz_col(0) / dtcld;
1557  work1_arr(i,j,klo,2) = delqrs3 / dz_col(0) / dtcld;
1558  });
1559 
1560 
1561  // ============================================================
1562  // Step 3g: Second slope_wdm6 call after sedimentation (G6)
1563  // Recompute slope parameters for all species after G5c sed.
1564  // ============================================================
1565 
1566  // Copy current state to temporary arrays (qrs_tmp, ncr_tmp)
1567  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1568  qrs_tmp_arr(i,j,k,0) = qr_arr(i,j,k); // Fortran qrs(:,:,1)
1569  qrs_tmp_arr(i,j,k,1) = qs_arr(i,j,k); // Fortran qrs(:,:,2)
1570  qrs_tmp_arr(i,j,k,2) = qg_arr(i,j,k); // Fortran qrs(:,:,3)
1571  ncr_tmp_arr(i,j,k) = nr_arr(i,j,k); // Fortran ncr(:,:,3)
1572  });
1573 
1574  // Call slope_wdm6: compute slope parameters for rain/snow/graupel
1575  // This parallels the Fortran call: slope_wdm6(qrs_tmp,ncr_tmp,...)
1576  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1577  Real rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3, rain_vt, rain_vtn;
1578  Real snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt, snow_n0sfac;
1579  Real graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt;
1580 
1581  // Compute slope parameters using temporary (post-sedimentation) values
1582  wdm6_slope_rain_cell(qrs_tmp_arr(i,j,k,0), ncr_tmp_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1583  Real(qcrmin), Real(nrmin),
1584  rslopermax_loc, rsloperbmax_loc, rsloper2max_loc, rsloper3max_loc,
1585  Real(bvtr), pvtr_loc, pvtrn_loc, pidnr_loc,
1586  rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3,
1587  rain_vt, rain_vtn);
1588  wdm6_slope_snow_cell(qrs_tmp_arr(i,j,k,1), den_arr(i,j,k), denfac_arr(i,j,k), t_arr(i,j,k),
1589  pidn0s_loc, Real(alpha_wdm6), Real(n0smax), Real(n0s),
1590  wdm6_slope_t0c, Real(qcrmin),
1591  rslopesmax_loc, rslopesbmax_loc, rslopes2max_loc, rslopes3max_loc,
1592  Real(bvts), pvts_loc,
1593  snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt,
1594  snow_n0sfac);
1595  wdm6_slope_graup_cell(qrs_tmp_arr(i,j,k,2), den_arr(i,j,k), denfac_arr(i,j,k),
1596  pidn0g_loc, Real(qcrmin),
1597  rslopegmax_loc, rslopegbmax_loc, rslopeg2max_loc, rslopeg3max_loc,
1598  slope_bvtg_loc, pvtg_loc,
1599  graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt);
1600 
1601  // Store results in output slope arrays
1602  rslope_arr(i,j,k,0) = rain_rslope;
1603  rslope_arr(i,j,k,1) = snow_rslope;
1604  rslope_arr(i,j,k,2) = graup_rslope;
1605  rslopeb_arr(i,j,k,0) = rain_rslopeb;
1606  rslopeb_arr(i,j,k,1) = snow_rslopeb;
1607  rslopeb_arr(i,j,k,2) = graup_rslopeb;
1608  rslope2_arr(i,j,k,0) = rain_rslope2;
1609  rslope2_arr(i,j,k,1) = snow_rslope2;
1610  rslope2_arr(i,j,k,2) = graup_rslope2;
1611  rslope3_arr(i,j,k,0) = rain_rslope3;
1612  rslope3_arr(i,j,k,1) = snow_rslope3;
1613  rslope3_arr(i,j,k,2) = graup_rslope3;
1614  work1_arr(i,j,k,0) = rain_vt;
1615  work1_arr(i,j,k,1) = snow_vt;
1616  work1_arr(i,j,k,2) = graup_vt;
1617  workn_arr(i,j,k) = rain_vtn;
1618  });
1619 
1620 
1621  // ============================================================
1622  // WDM6-CPP_PRE_G7
1623  // ============================================================
1624 
1625  // ============================================================
1626  // Step 3h: G7 warm-phase snow/graupel melting
1627  // ============================================================
1628  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1629  // Warm-phase snow and graupel melting (t > t0c)
1630  if (t_arr(i,j,k) > t0c) {
1631  const Real supcol = t0c - t_arr(i,j,k);
1632  const Real n0sfac = amrex::max(
1633  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
1634  Real(n0smax) / Real(n0s)),
1635  Real(1.0));
1636 
1637  const Real xlf = xlf0;
1638  const Real work2 = wdm6_venfac(p_arr(i,j,k), t_arr(i,j,k),
1639  den_arr(i,j,k), Real(den0));
1640 
1641  // --- Snow melting: psmlt ---
1642  if (qs_arr(i,j,k) > Real(0.0)) {
1643  const Real coeres_s = rslope2_arr(i,j,k,1) *
1644  std::sqrt(rslope_arr(i,j,k,1) * rslopeb_arr(i,j,k,1));
1645 
1646  Real psmlt = wdm6_xka(t_arr(i,j,k), den_arr(i,j,k)) / xlf *
1647  (t0c - t_arr(i,j,k)) * pi_wdm6_loc * Real(0.5) * n0sfac *
1648  (precs1_loc * rslope2_arr(i,j,k,1) +
1649  precs2_loc * work2 * coeres_s) / den_arr(i,j,k);
1650 
1651  psmlt = amrex::min(amrex::max(psmlt * dtcld, -qs_arr(i,j,k)),
1652  Real(0.0));
1653 
1654  if (qs_arr(i,j,k) > Real(qcrmin)) {
1655  const Real sfac = rslope_arr(i,j,k,1) * Real(n0s) * n0sfac /
1656  qs_arr(i,j,k);
1657  nr_arr(i,j,k) = amrex::max(nr_arr(i,j,k) - sfac * psmlt, Real(0.0));
1658  }
1659 
1660  qs_arr(i,j,k) += psmlt;
1661  qr_arr(i,j,k) -= psmlt;
1662  t_arr(i,j,k) += xlf / cpm_arr(i,j,k) * psmlt;
1663  }
1664 
1665  // --- Graupel melting: pgmlt ---
1666  if (qg_arr(i,j,k) > Real(0.0)) {
1667  const Real coeres_g = rslope2_arr(i,j,k,2) *
1668  std::sqrt(rslope_arr(i,j,k,2) * rslopeb_arr(i,j,k,2));
1669 
1670  Real pgmlt = wdm6_xka(t_arr(i,j,k), den_arr(i,j,k)) / xlf *
1671  (t0c - t_arr(i,j,k)) *
1672  (precg1_loc * rslope2_arr(i,j,k,2) +
1673  precg2_loc * work2 * coeres_g) / den_arr(i,j,k);
1674 
1675  pgmlt = amrex::min(amrex::max(pgmlt * dtcld, -qg_arr(i,j,k)),
1676  Real(0.0));
1677 
1678  if (qg_arr(i,j,k) > Real(qcrmin)) {
1679  const Real gfac = rslope_arr(i,j,k,2) * n0g_loc /
1680  qg_arr(i,j,k);
1681  nr_arr(i,j,k) = amrex::max(nr_arr(i,j,k) - gfac * pgmlt, Real(0.0));
1682  }
1683 
1684  qg_arr(i,j,k) += pgmlt;
1685  qr_arr(i,j,k) -= pgmlt;
1686  t_arr(i,j,k) += xlf / cpm_arr(i,j,k) * pgmlt;
1687  }
1688  }
1689  });
1690 
1691  // ============================================================
1692  // WDM6-CPP_POST_G7
1693  // ============================================================
1694 
1695  // ============================================================
1696  // Step 3i: G8 ice fall speed + ice sedimentation
1697  // VICE: Ice crystal terminal velocity + nislfv_rain_plmr
1698  // ============================================================
1699 
1700  // Compute ice crystal fall speed (work1c). Fortran:
1701  // xmi = den*qci(:,:,2)/xni
1702  // diameter = max(min(dicon*sqrt(xmi),dimax),1.e-25)
1703  // work1c = 1.49e4*exp(log(diameter)*(1.31))
1704  // dicon and dimax are the CLASS constants, already routed through
1705  // wdm6_literal. This block previously declared its own locals that
1706  // shadowed them, and both were wrong: dicon was 11.45e-9 against the
1707  // Fortran's 11.9, nine orders of magnitude small, which collapsed the
1708  // ice fall speed to ~5e-13 m/s so the native path effectively did not
1709  // sediment ice at all; and dimax was a true double, reintroducing at
1710  // this site exactly the literal-precision defect already closed in
1711  // G13e. Hoisted out of the lambda in the established _loc style so no
1712  // capture of `this` is needed on device.
1713  constexpr Real dicon_loc = dicon;
1714  constexpr Real dimax_loc = dimax;
1715  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1716  Real work1c = Real(0.0);
1717  if (qi_arr(i,j,k) > Real(0.0)) {
1718  const Real xni_safe = amrex::max(xni_arr(i,j,k), Real(1.0e-30));
1719  const Real xmi = den_arr(i,j,k) * qi_arr(i,j,k) / xni_safe;
1720  const Real diameter = amrex::max(amrex::min(dicon_loc * std::sqrt(xmi), dimax_loc),
1721  wdm6_literal(1.0e-25));
1722  work1c = wdm6_literal(1.49e4) * std::exp(std::log(diameter) * wdm6_literal(1.31));
1723  }
1724  work1c_arr(i,j,k) = work1c;
1725  });
1726 
1727  // Ice sedimentation via simplified PLM6 scheme
1728  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1729  amrex::ignore_unused(k);
1730  const int km = khi - klo + 1;
1731  auto dz_col = [&](int kk) -> amrex::Real& {
1732  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::dz);
1733  };
1734  auto den_col = [&](int kk) -> amrex::Real& {
1735  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::den);
1736  };
1737  auto denfac_col = [&](int kk) -> amrex::Real& {
1738  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::denfac);
1739  };
1740  auto tk_col = [&](int kk) -> amrex::Real& {
1741  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::tk);
1742  };
1743  auto work_col = [&](int kk) -> amrex::Real& {
1744  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::work_col);
1745  };
1746  auto rq_col = [&](int kk) -> amrex::Real& {
1747  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::rq_col);
1748  };
1749  auto rq2_col = [&](int kk) -> amrex::Real& {
1750  return sed_cell_scratch_arr(i, j, klo + kk, WDM6SedCellScratch::rq2_col);
1751  };
1752 
1753  // Fill column arrays
1754  for (int kk = 0; kk < km; ++kk) {
1755  const int k3 = klo + kk;
1756  dz_col(kk) = delz_arr(i,j,k3);
1757  den_col(kk) = den_arr(i,j,k3);
1758  denfac_col(kk) = denfac_arr(i,j,k3);
1759  tk_col(kk) = t_arr(i,j,k3);
1760  work_col(kk) = work1c_arr(i,j,k3);
1761  rq_col(kk) = den_col(kk) * qi_arr(i,j,k3);
1762  rq2_col(kk) = den_col(kk) * qi_arr(i,j,k3);
1763  }
1764 
1765  Real delqi_col = Real(0.0);
1766  Real delqi2_col = Real(0.0);
1767 
1768  // Ice sedimentation. ITER MUST BE 0 HERE. The Fortran calls
1769  // nislfv_rain_plmr(...,work1c,denqci,denqci,delqi,dtcld,1,0,0)
1770  // whose trailing arguments are dt, id, iter, rid -- so id=1 and
1771  // iter=0. This helper's signature has no id parameter and reads
1772  // (..., dt, iter, ...), so passing the Fortran's id of 1 in that
1773  // slot silently set iter=1 and ran the velocity-iteration block
1774  // that G8 is supposed to skip, recomputing the ICE fall speed
1775  // from the snow and graupel slope routines. Contrast G5c, which
1776  // genuinely passes iter=1 (Fortran :1122 ends dtcld,1,1).
1777  // This was verified against the Fortran bridge by tracing every
1778  // input with prints archived on a branch: qi, xni, den, xmi,
1779  // dicon, diameter and work1c are all bitwise equal at every k in
1780  // the 85..100 window, so the kernel invocation was the only
1781  // remaining candidate.
1782  wdm6_nislfv_rain_plm6_column(
1783  km, &delqi_col, &delqi2_col, dtcld, 0,
1784  pidn0s_loc, pidn0g_loc, Real(qcrmin), Real(alpha_wdm6), Real(n0smax), Real(n0s), wdm6_slope_t0c,
1785  rslopesmax_loc, rslopesbmax_loc, rslopes2max_loc, rslopes3max_loc, Real(bvts), pvts_loc,
1786  rslopegmax_loc, rslopegbmax_loc, rslopeg2max_loc, rslopeg3max_loc, slope_bvtg_loc, pvtg_loc,
1787  sed_cell_scratch_arr, sed_node_scratch_arr, i, j, klo);
1788 
1789  // Update ice concentrations from sedimented state
1790  for (int kk = 0; kk < km; ++kk) {
1791  const int k3 = klo + kk;
1792  qi_arr(i,j,k3) = amrex::max(rq_col(kk) / den_col(kk), Real(0.0));
1793  }
1794 
1795  // Ice fallout at surface. Indexed at klo to match the G9 block
1796  // below, which reads fallc_arr(i,j,klo); this lambda discards
1797  // its own k (see ignore_unused above), so spelling the index
1798  // out keeps the two sites textually identical.
1799  fallc_arr(i,j,klo) = delqi_col / dz_col(0) / dtcld;
1800  delqi_arr(i,j,klo) = delqi_col;
1801  });
1802 
1803 
1804  // ============================================================
1805  // Step 3h: G9 Surface precipitation accumulation and sr update
1806  // ============================================================
1807 
1808  // Surface precipitation accumulation from sedimentation fallout
1809  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
1810  {
1811  const Real fall_r = work1_arr(i,j,klo,0); // rain fallout
1812  const Real fall_s = work1_arr(i,j,klo,1); // snow fallout
1813  const Real fall_g = work1_arr(i,j,klo,2); // graupel fallout
1814  const Real fall_c = fallc_arr(i,j,klo); // ice fallout
1815 
1816  const Real fallsum = fall_r + fall_s + fall_g + fall_c;
1817  const Real fallsum_qsi = fall_s + fall_c;
1818  const Real fallsum_qg = fall_g;
1819 
1820  const Real conv = delz_arr(i,j,klo) / Real(rhoh2o) * dtcld * Real(1000.0);
1821 
1822  if (fallsum > Real(0.0)) {
1823  rain_arr(i,j,klo) += fallsum * conv;
1824  }
1825 
1826  if (fallsum_qsi > Real(0.0)) {
1827  snow_arr(i,j,klo) += fallsum_qsi * conv;
1828  }
1829 
1830  if (fallsum_qg > Real(0.0)) {
1831  graup_arr(i,j,klo) += fallsum_qg * conv;
1832  }
1833 
1834  if (fallsum > Real(0.0)) {
1835  sr_arr(i,j,klo) = (snow_arr(i,j,klo) + graup_arr(i,j,klo))
1836  / (rain_arr(i,j,klo) + wdm6_literal(1.0e-12));
1837  }
1838  });
1839 
1840 
1841  // ============================================================
1842  // Step 3j: G10a Cloud-Ice Melt to Cloud Water
1843  // Exact port of bounded Fortran block (lines 1189-1197)
1844  // ============================================================
1845 
1846  // ================================================================
1847  // G10a-G10d FUSED, matching the Fortran's single loop nest.
1848  //
1849  // ERF_module_mp_wdm6.F90:1420-1421 opens one do k / do i nest and
1850  // computes supcol (:1422) and xlf (:1423-1424) ONCE per cell; all
1851  // four phase blocks then reuse them while t is progressively
1852  // updated by latent heat (:1432 cools, :1459 warms, :1495 warms).
1853  // This path previously ran the four blocks as four separate
1854  // ParallelFor kernels, each recomputing supcol and xlf from the
1855  // already-updated t_arr. That made pgfrz consume a fresher supcol
1856  // than the Fortran does: measured 1.048787e-06 on supcolt at
1857  // (107,3,46), amplified to 4.908569e-06 in expterm by the pfrz2
1858  // factor and inherited by pfrzdtr.
1859  //
1860  // xlf is the more dangerous of the two because it is set by a
1861  // BRANCH on supcol. Melting cools and freezing warms, so t crossing
1862  // t0c mid-section is exactly what this code does, and a per-kernel
1863  // re-evaluation would flip the branch where the Fortran does not.
1864  // Note G10b did not even carry the branch, only xlf = xls - xl.
1865  //
1866  // The group-boundary tags stay at their group boundaries per
1867  // Rule 30; they are emitted per cell from inside the kernel now
1868  // rather than from a column sweep afterwards. For these purely
1869  // per-cell blocks that is the same state, and it is strictly closer
1870  // to the "immediate process-group boundary" the rule asks for.
1871  // ================================================================
1872  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1873  const Real supcol = t0c - t_arr(i,j,k);
1874  Real xlf = xls - xl_arr(i,j,k);
1875  if (supcol < Real(0.0)) xlf = xlf0;
1876 
1877  // Cloud-ice melt to cloud water (T > t0c)
1878  if (supcol < Real(0.0) && qi_arr(i,j,k) > Real(0.0)) {
1879  const Real qim = qi_arr(i,j,k); // preserve old ice amount
1880 
1881  qc_arr(i,j,k) += qim; // qci(:,:,1) += qci(:,:,2)
1882  if (qim > Real(qmin)) {
1883  nc_arr(i,j,k) += xni_arr(i,j,k); // ncr(:,:,2) += xni
1884  }
1885  t_arr(i,j,k) -= xlf / cpm_arr(i,j,k) * qim; // latent heat release
1886  qi_arr(i,j,k) = Real(0.0); // zero out ice
1887  }
1888 
1889 
1890 
1891  // ============================================================
1892  // Step 3k: G10b Cloud-Water Homogeneous Freezing
1893  // Exact port of bounded Fortran block (lines 1215-1224)
1894  // ============================================================
1895 
1896 
1897  // Homogeneous freezing of cloud water when supcol > 40K (T < -40C)
1898  if (supcol > Real(40.0) && qc_arr(i,j,k) > Real(0.0)) {
1899  const Real qc_old = qc_arr(i,j,k);
1900 
1901  qi_arr(i,j,k) += qc_old;
1902 
1903  if (nc_arr(i,j,k) > Real(0.0)) {
1904  nc_arr(i,j,k) = Real(0.0);
1905  }
1906 
1907  t_arr(i,j,k) += xlf / cpm_arr(i,j,k) * qc_old;
1908  qc_arr(i,j,k) = Real(0.0);
1909  }
1910 
1911 
1912  // ============================================================
1913  // Step 3l: G10c Cloud-Water Heterogeneous Freezing
1914  // Exact port of bounded Fortran block (pihtf, lines 1241-1258)
1915  // ============================================================
1916 
1917  // Latent heat of FUSION, per the Fortran's per-k prologue:
1918  // xlf = xls-xl(i,k) ; if(supcol.lt.0.) xlf = xlf0
1919  // shared by G10a-G10d. Freezing releases xlf, not xl.
1920 
1921  // Heterogeneous freezing (Biggs, contact/immersion): 0 > T > -40C
1922  // Trigger: supcol > 0 (T < T0c) AND qc > qmin
1923  if (supcol > Real(0.0) && qc_arr(i,j,k) > Real(qmin)) {
1924  const Real supcolt = amrex::min(supcol, Real(70.0));
1925  const Real expterm = std::exp(pfrz2_loc * supcolt) - Real(1.0);
1926 
1927  // Use the stored rslopec3, as the Fortran does: in the
1928  // rslopecmax branch it is the precomputed rslopec3max
1929  // constant, which need not equal rslopecmax cubed.
1930  const Real rs3 = rslopec3_arr(i,j,k);
1931 
1932  // pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den
1933  // *ncr(:,:,2)*rslopec3*rslopec3/18.*dtcld, qci(:,:,1))
1934  // rslopec3 appears TWICE. Using it once left the rate about
1935  // 1/rslopec3 too large -- roughly 8e12 for a typical cloud
1936  // droplet slope -- so pfrzdtc always saturated at its qc cap
1937  // and every cell below freezing lost all its cloud water to
1938  // ice in a single substep. Factor order follows the Fortran
1939  // left to right so the roundings agree.
1940  Real pfrzdtc = pi_wdm6_loc * pi_wdm6_loc * pfrz1_loc * expterm
1941  * denr / den_arr(i,j,k) * nc_arr(i,j,k) * rs3 * rs3
1942  / Real(18.0) * dtcld;
1943  pfrzdtc = amrex::min(pfrzdtc, qc_arr(i,j,k));
1944 
1945  // nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(:,:,2)
1946  // *rslopec3/6.*dtcld, ncr(:,:,2)) -- one factor here.
1947  Real nfrzdtc = pi_wdm6_loc * pfrz1_loc * expterm
1948  * nc_arr(i,j,k) * rs3
1949  / Real(6.0) * dtcld;
1950  nfrzdtc = amrex::min(nfrzdtc, nc_arr(i,j,k));
1951 
1952  // Apply number loss (only if nc > ncmin)
1953  if (nc_arr(i,j,k) > Real(ncmin)) {
1954  nc_arr(i,j,k) -= nfrzdtc;
1955  }
1956 
1957  // Apply mass and temperature updates
1958  qi_arr(i,j,k) += pfrzdtc;
1959  t_arr(i,j,k) += xlf / cpm_arr(i,j,k) * pfrzdtc;
1960  qc_arr(i,j,k) -= pfrzdtc;
1961  }
1962 
1963 
1964  // ============================================================
1965  // Step 3m: G10d Rain-to-Graupel Freezing
1966  // Exact port of bounded Fortran block (pgfrz, rain freezing)
1967  // ============================================================
1968 
1969  // Latent heat of fusion, as in G10a-G10c above.
1970 
1971  // Rain freezing to graupel: trigger when T < t0c and qr > 0
1972  if (supcol > Real(0.0) && qr_arr(i,j,k) > Real(0.0)) {
1973  const Real supcolt = amrex::min(supcol, Real(70.0));
1974  const Real expterm = std::exp(pfrz2_loc * supcolt) - Real(1.0);
1975 
1976  // Rain slope cubed (from G4/G6 computation)
1977  const Real rs3 = rslope3_arr(i,j,k,0);
1978 
1979  // pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(:,:,3)*denr/den
1980  // *(exp(pfrz2*supcolt)-1.)*rslope3*rslope3*dtcld, qrs(:,:,1))
1981  // Factor order follows the Fortran left to right.
1982  Real pfrzdtr = Real(140.0) * (pi_wdm6_loc * pi_wdm6_loc)
1983  * pfrz1_loc * nr_arr(i,j,k)
1984  * denr / den_arr(i,j,k)
1985  * expterm * rs3 * rs3 * dtcld;
1986  pfrzdtr = amrex::min(pfrzdtr, qr_arr(i,j,k));
1987 
1988 
1989  // Number freezing rate (conditional on nr > nrmin)
1990  if (nr_arr(i,j,k) > Real(nrmin)) {
1991  Real nfrzdtr = Real(4.0) * pi_wdm6_loc * pfrz1_loc
1992  * nr_arr(i,j,k) * expterm * rs3 * dtcld;
1993  nfrzdtr = amrex::min(nfrzdtr, nr_arr(i,j,k));
1994  nr_arr(i,j,k) -= nfrzdtr;
1995  }
1996 
1997  // Apply mass and temperature updates
1998  qg_arr(i,j,k) += pfrzdtr;
1999  t_arr(i,j,k) += xlf / cpm_arr(i,j,k) * pfrzdtr;
2000  qr_arr(i,j,k) -= pfrzdtr;
2001  }
2002  });
2003 
2004 
2005  // ============================================================
2006  // Step 3k: G10e Phase cleanup — clamp number concentrations non-negative
2007  // ============================================================
2008 
2009  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2010  nc_arr(i,j,k) = amrex::max(nc_arr(i,j,k), Real(0.0));
2011  nr_arr(i,j,k) = amrex::max(nr_arr(i,j,k), Real(0.0));
2012  });
2013 
2014 
2015  // ============================================================
2016  // G11: SLOPE3 — Third slope_wdm6 call + avedia/rslopec recompute
2017  // ============================================================
2018 
2019 
2020  // G11(a): Repack fields for slope recomputation (Fortran qrs_tmp/ncr_tmp)
2021  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2022  qrs_tmp_arr(i,j,k,0) = qr_arr(i,j,k); // qrs(:,:,1)
2023  qrs_tmp_arr(i,j,k,1) = qs_arr(i,j,k); // qrs(:,:,2)
2024  qrs_tmp_arr(i,j,k,2) = qg_arr(i,j,k); // qrs(:,:,3)
2025  ncr_tmp_arr(i,j,k) = nr_arr(i,j,k); // ncr(:,:,3)
2026  });
2027 
2028  // G11(b): slope_wdm6-equivalent recompute for rain/snow/graupel
2029  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2030  Real rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3, rain_vt, rain_vtn;
2031  Real snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt, snow_n0sfac;
2032  Real graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt;
2033 
2034  // Compute slope parameters using qrs_tmp/ncr_tmp values
2035  wdm6_slope_rain_cell(qrs_tmp_arr(i,j,k,0), ncr_tmp_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
2036  Real(qcrmin), Real(nrmin),
2037  rslopermax_loc, rsloperbmax_loc, rsloper2max_loc, rsloper3max_loc,
2038  Real(bvtr), pvtr_loc, pvtrn_loc, pidnr_loc,
2039  rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3,
2040  rain_vt, rain_vtn);
2041  wdm6_slope_snow_cell(qrs_tmp_arr(i,j,k,1), den_arr(i,j,k), denfac_arr(i,j,k), t_arr(i,j,k),
2042  pidn0s_loc, Real(alpha_wdm6), Real(n0smax), Real(n0s),
2043  wdm6_slope_t0c, Real(qcrmin),
2044  rslopesmax_loc, rslopesbmax_loc, rslopes2max_loc, rslopes3max_loc,
2045  Real(bvts), pvts_loc,
2046  snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt,
2047  snow_n0sfac);
2048  wdm6_slope_graup_cell(qrs_tmp_arr(i,j,k,2), den_arr(i,j,k), denfac_arr(i,j,k),
2049  pidn0g_loc, Real(qcrmin),
2050  rslopegmax_loc, rslopegbmax_loc, rslopeg2max_loc, rslopeg3max_loc,
2051  slope_bvtg_loc, pvtg_loc,
2052  graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt);
2053 
2054  // Store results in output slope arrays
2055  rslope_arr(i,j,k,0) = rain_rslope;
2056  rslope_arr(i,j,k,1) = snow_rslope;
2057  rslope_arr(i,j,k,2) = graup_rslope;
2058  rslopeb_arr(i,j,k,0) = rain_rslopeb;
2059  rslopeb_arr(i,j,k,1) = snow_rslopeb;
2060  rslopeb_arr(i,j,k,2) = graup_rslopeb;
2061  rslope2_arr(i,j,k,0) = rain_rslope2;
2062  rslope2_arr(i,j,k,1) = snow_rslope2;
2063  rslope2_arr(i,j,k,2) = graup_rslope2;
2064  rslope3_arr(i,j,k,0) = rain_rslope3;
2065  rslope3_arr(i,j,k,1) = snow_rslope3;
2066  rslope3_arr(i,j,k,2) = graup_rslope3;
2067  work1_arr(i,j,k,0) = rain_vt;
2068  work1_arr(i,j,k,1) = snow_vt;
2069  work1_arr(i,j,k,2) = graup_vt;
2070  workn_arr(i,j,k) = rain_vtn;
2071  });
2072 
2073  // G11(c): avedia + rslopec/2/3 recompute (lamdac fallback branch)
2074  const Real cbrt24 = wdm6_default_real_pow(24.0, 0.3333333);
2075  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2076  // avedia component from rain slope: avedia(:,:,2) = rslope(:,:,1) * (24.)^(1/3)
2077  avedia_arr(i,j,k,1) = rslope_arr(i,j,k,0) * cbrt24;
2078 
2079  // Cloud water slope parameters: rslopec = 1/lamdac(qci(:,:,1), den, ncr(:,:,2))
2080  // Fortran gate: if (qci(:,:,1) <= qmin or ncr(:,:,2) <= ncmin) use maxima, else use lamdac
2081  const Real qci_for_lamdac = qc_arr(i,j,k);
2082  const Real nc_for_lamdac = nc_arr(i,j,k);
2083 
2084  if (qci_for_lamdac <= Real(qmin) || nc_for_lamdac <= Real(ncmin)) {
2085  rslopec_arr(i,j,k) = rslopecmax_loc;
2086  rslopec2_arr(i,j,k) = rslopec2max_loc;
2087  rslopec3_arr(i,j,k) = rslopec3max_loc;
2088  } else {
2089  // Same expression as G3; share the one exact port. The
2090  // previous wdm6_lamdac() here differed from the Fortran on
2091  // three counts: it cancelled den out of the argument, it
2092  // used std::pow with a true 1/3 rather than exp/log with the
2093  // unsuffixed .33333333, and it carried a qc/nc guard of its
2094  // own that fired inside this block's own gate.
2095  const Real rslc = wdm6_rslopec_exact(qci_for_lamdac, den_arr(i,j,k),
2096  nc_for_lamdac, pidnc_loc);
2097  rslopec_arr(i,j,k) = rslc;
2098  rslopec2_arr(i,j,k) = rslc * rslc;
2099  rslopec3_arr(i,j,k) = rslopec2_arr(i,j,k) * rslc;
2100  }
2101 
2102  // avedia component from cloud slope: avedia(:,:,1) = rslopec(:,:)
2103  avedia_arr(i,j,k,0) = rslopec_arr(i,j,k);
2104  });
2105 
2106  // G11(d): work arrays recompute via diffac/venfac
2107  // work1(:,:,1) = diffac(xl, p, t, den, qs(:,:,1), rv)
2108  // work1(:,:,2) = diffac(xls, p, t, den, qs(:,:,2), rv)
2109  // work2(:,:) = venfac(p, t, den, den0) [stored for G13a warm-rain rates]
2110  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2111  // Recompute work1 diffusion coefficients using current state
2112  if (i == diag_i && j == diag_j && k == diag_k) {
2113  }
2114  work1_arr(i,j,k,0) = wdm6_diffac(xl_arr(i,j,k), p_arr(i,j,k), t_arr(i,j,k),
2115  den_arr(i,j,k), qsatw_arr(i,j,k), Real(rv)); // qs(:,:,1)
2116  if (i == diag_i && j == diag_j && k == diag_k) {
2117  }
2118  if (i == diag_i && j == diag_j && k == diag_k) {
2119  }
2120  work1_arr(i,j,k,1) = wdm6_diffac(Real(xls), p_arr(i,j,k), t_arr(i,j,k),
2121  den_arr(i,j,k), qsati_arr(i,j,k), Real(rv)); // qs(:,:,2)
2122  if (i == diag_i && j == diag_j && k == diag_k) {
2123  }
2124  // Compute work2 ventilation factor for diffusion (used in G13a warm-rain rates)
2125  work2_arr(i,j,k) = wdm6_venfac(p_arr(i,j,k), t_arr(i,j,k), den_arr(i,j,k), Real(den0));
2126  });
2127 
2128 
2129 
2130 
2131  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2132  const Real supsat = amrex::max(qv_arr(i,j,k), Real(qmin)) - qsatw_arr(i,j,k);
2133  const Real satdt = supsat / dtcld;
2134  const Real lencon = wdm6_literal(2.7e-2) * den_arr(i,j,k) * qc_arr(i,j,k)
2135  * (wdm6_literal(1.0e20/16.0) * rslopec2_arr(i,j,k) * rslopec2_arr(i,j,k)
2136  - wdm6_literal(0.4));
2137  const Real lenconcr = amrex::max(wdm6_literal(1.2) * lencon, Real(qcrmin));
2138 
2139  if (qc_arr(i,j,k) > qcr_arr(i,j,k) && nc_arr(i,j,k) > Real(ncmin)) {
2140  // Fortran :1889
2141  // praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.)
2142  // Both exponents are unsuffixed, so they are float32 values
2143  // widened, not true thirds: 7./3. is 2.3333332538604736
2144  // against 2.3333333333333335, and 1./3. is 0.3333333432674408
2145  // against 0.33333333333333331. Exponents amplify: the error
2146  // enters as ln(base)*delta, and with qc about 1.6e-02 and nc
2147  // about 1.8e+08 those logs are -4.2 and +19.0, which is why a
2148  // 3e-08 exponent difference produced 1.86e-07 on praut.
2149  praut_arr(i,j,k) = qck1_loc * std::pow(qc_arr(i,j,k), wdm6_literal(7.0/3.0))
2150  * std::pow(nc_arr(i,j,k), wdm6_literal(-1.0/3.0));
2151  praut_arr(i,j,k) = amrex::min(praut_arr(i,j,k), qc_arr(i,j,k) / dtcld);
2152 
2153  nrauto_arr(i,j,k) = Real(3.5e9) * den_arr(i,j,k) * praut_arr(i,j,k);
2154  if (qr_arr(i,j,k) > lenconcr) {
2155  nrauto_arr(i,j,k) = nr_arr(i,j,k) / qr_arr(i,j,k) * praut_arr(i,j,k);
2156  }
2157  nrauto_arr(i,j,k) = amrex::min(nrauto_arr(i,j,k), nc_arr(i,j,k) / dtcld);
2158  }
2159 
2160  if (qr_arr(i,j,k) >= lenconcr) {
2161  if (avedia_arr(i,j,k,1) >= Real(di100)) {
2162  nraccr_arr(i,j,k) = amrex::min(
2163  Real(ncrk1) * nc_arr(i,j,k) * nr_arr(i,j,k)
2164  * (rslopec3_arr(i,j,k) + Real(24.0) * rslope3_arr(i,j,k,0)),
2165  nc_arr(i,j,k) / dtcld);
2166  pracw_arr(i,j,k) = amrex::min(
2167  pi_wdm6_loc / Real(6.0) * (Real(denr) / den_arr(i,j,k))
2168  * Real(ncrk1) * nc_arr(i,j,k) * nr_arr(i,j,k)
2169  * rslopec3_arr(i,j,k)
2170  * (Real(2.0) * rslopec3_arr(i,j,k) + Real(24.0) * rslope3_arr(i,j,k,0)),
2171  qc_arr(i,j,k) / dtcld);
2172  } else {
2173  nraccr_arr(i,j,k) = amrex::min(
2174  Real(ncrk2) * nc_arr(i,j,k) * nr_arr(i,j,k)
2175  * (Real(2.0) * rslopec3_arr(i,j,k) * rslopec3_arr(i,j,k)
2176  + Real(5040.0) * rslope3_arr(i,j,k,0) * rslope3_arr(i,j,k,0)),
2177  nc_arr(i,j,k) / dtcld);
2178  pracw_arr(i,j,k) = amrex::min(
2179  pi_wdm6_loc / Real(6.0) * (Real(denr) / den_arr(i,j,k))
2180  * Real(ncrk2) * nc_arr(i,j,k) * nr_arr(i,j,k)
2181  * rslopec3_arr(i,j,k)
2182  * (Real(6.0) * rslopec3_arr(i,j,k) * rslopec3_arr(i,j,k)
2183  + Real(5040.0) * rslope3_arr(i,j,k,0) * rslope3_arr(i,j,k,0)),
2184  qc_arr(i,j,k) / dtcld);
2185  }
2186  }
2187 
2188  if (avedia_arr(i,j,k,0) >= Real(di100)) {
2189  nccol_arr(i,j,k) = Real(ncrk1) * nc_arr(i,j,k) * nc_arr(i,j,k) * rslopec3_arr(i,j,k);
2190  } else {
2191  nccol_arr(i,j,k) = Real(2.0) * Real(ncrk2) * nc_arr(i,j,k) * nc_arr(i,j,k)
2192  * rslopec3_arr(i,j,k) * rslopec3_arr(i,j,k);
2193  }
2194 
2195  if (qr_arr(i,j,k) >= lenconcr) {
2196  if (avedia_arr(i,j,k,1) < Real(di100)) {
2197  nrcol_arr(i,j,k) = Real(5040.0) * Real(ncrk2) * nr_arr(i,j,k) * nr_arr(i,j,k)
2198  * rslope3_arr(i,j,k,0) * rslope3_arr(i,j,k,0);
2199  } else if (avedia_arr(i,j,k,1) < Real(di600)) {
2200  nrcol_arr(i,j,k) = Real(24.0) * Real(ncrk1) * nr_arr(i,j,k) * nr_arr(i,j,k)
2201  * rslope3_arr(i,j,k,0);
2202  } else if (avedia_arr(i,j,k,1) < Real(di2000)) {
2203  const Real coecol = -Real(2.5e3) * (avedia_arr(i,j,k,1) - Real(di600));
2204  nrcol_arr(i,j,k) = Real(24.0) * std::exp(coecol) * Real(ncrk1)
2205  * nr_arr(i,j,k) * nr_arr(i,j,k) * rslope3_arr(i,j,k,0);
2206  } else {
2207  nrcol_arr(i,j,k) = Real(0.0);
2208  }
2209  }
2210 
2211  if (qr_arr(i,j,k) > Real(0.0)) {
2212  const Real coeres = rslope_arr(i,j,k,0)
2213  * std::sqrt(rslope_arr(i,j,k,0) * rslopeb_arr(i,j,k,0));
2214  prevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0)) * nr_arr(i,j,k)
2215  * (precr1_loc * rslope_arr(i,j,k,0) + precr2_loc * work2_arr(i,j,k) * coeres)
2216  / work1_arr(i,j,k,0);
2217  if (prevp_arr(i,j,k) < Real(0.0)) {
2218  prevp_arr(i,j,k) = amrex::max(prevp_arr(i,j,k), -qr_arr(i,j,k) / dtcld);
2219  prevp_arr(i,j,k) = amrex::max(prevp_arr(i,j,k), satdt / Real(2.0));
2220 
2221  if (prevp_arr(i,j,k) == -qr_arr(i,j,k) / dtcld) {
2222  nn_arr(i,j,k) = nn_arr(i,j,k) + nr_arr(i,j,k);
2223  nr_arr(i,j,k) = Real(0.0);
2224  }
2225  } else if (prevp_arr(i,j,k) == Real(0.0)) {
2226  // A zero kinetic rate must not become evaporation
2227  // through a negative saturation-rate limiter.
2228  prevp_arr(i,j,k) = Real(0.0);
2229  } else {
2230  prevp_arr(i,j,k) = amrex::min(prevp_arr(i,j,k), satdt / Real(2.0));
2231  }
2232  }
2233  });
2234 
2235 
2236 
2237  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2238  const Real supcol = Real(t0c) - t_arr(i,j,k);
2239  const Real n0sfac = amrex::max(
2240  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
2241  Real(n0smax) / Real(n0s)),
2242  Real(1.0));
2243  const Real supsat = amrex::max(qv_arr(i,j,k), Real(qmin)) - qsati_arr(i,j,k);
2244  amrex::ignore_unused(supsat);
2245  const Real satdt = supsat / dtcld;
2246  amrex::ignore_unused(satdt);
2247 
2248  const Real qi_val = qi_arr(i,j,k);
2249  if (!(supcol > Real(0.0) && qi_val > Real(qmin))) {
2250  return;
2251  }
2252 
2253  Real temp = den_arr(i,j,k) * amrex::max(qi_val, Real(qmin));
2254  temp = std::sqrt(std::sqrt(temp * temp * temp));
2255  xni_arr(i,j,k) = amrex::min(amrex::max(Real(5.38e7) * temp, Real(1.e3)), Real(1.e6));
2256 
2257  // Fortran :2021 eacrs = exp(0.07*(-supcol)); 0.07 is unsuffixed
2258  // so it carries the float assumption. Verified by tracing psaci's
2259  // inputs: eacrs diverged 1.03e-08 to 1.14e-08 relative across
2260  // k=89..95, matching delta(0.07)*(-supcol).
2261  const Real eacrs = std::exp(wdm6_literal(0.07) * (-supcol));
2262  const Real xni_safe = amrex::max(xni_arr(i,j,k), Real(1.0e-30));
2263  const Real xmi = den_arr(i,j,k) * qi_val / xni_safe;
2264  const Real diameter = amrex::min(Real(dicon) * std::sqrt(xmi), Real(dimax));
2265  // Fortran :2025 vt2i = 1.49e4*diameter**1.31. 1.49e4 is exactly
2266  // representable in float32 so it needs no routing, but the 1.31
2267  // EXPONENT does, and exponents amplify: the error enters as
2268  // ln(diameter)*delta, and with diameter at the dimax cap of
2269  // 5e-4 that is ln(5e-4) = -7.60 times -5.722e-08, i.e. 4.349e-07.
2270  // Measured: exactly 4.349270e-07 on vt2i across k=91..95.
2271  // It then amplifies again through the |vt2ave - vt2i|
2272  // cancellation in psaci: at k=89 vt2diff is 0.139 against vt2i
2273  // 0.681, a 4.90x factor giving the observed 2.14e-06.
2274  const Real vt2i = Real(1.49e4) * std::pow(diameter, wdm6_literal(1.31));
2275  const Real vt2r = pvtr_loc * rslopeb_arr(i,j,k,0) * denfac_arr(i,j,k);
2276  const Real vt2s = pvts_loc * rslopeb_arr(i,j,k,1) * denfac_arr(i,j,k);
2277  const Real vt2g = pvtg_loc * rslopeb_arr(i,j,k,2) * denfac_arr(i,j,k);
2278  // The 1.e-15 qsum floor is unsuffixed in the Fortran (:1110, :2072,
2279  // :2223, :3979), so it is float32(1e-15) = 1.0000000036274937e-15,
2280  // which is LARGER than the exact double by 3.63e-09 relative.
2281  //
2282  // It was previously judged inert because the floor and the gate use
2283  // the same literal, so a floored qsum can never pass `> literal`.
2284  // That reasoning holds WITHIN a leg but misses the cross-leg case:
2285  // the THRESHOLD itself differed between legs, so a cell with
2286  // qs+qg in (1e-15, 1.0000000036274937e-15) took the branch on the
2287  // port and skipped it on the oracle. A narrow window, but a branch
2288  // divergence rather than a rounding difference, so it is routed
2289  // rather than argued about.
2290  const Real qsum = amrex::max(qs_arr(i,j,k) + qg_arr(i,j,k), wdm6_literal(1.0e-15));
2291  const Real vt2ave = (qsum > wdm6_literal(1.0e-15))
2292  ? (vt2s * qs_arr(i,j,k) + vt2g * qg_arr(i,j,k)) / qsum
2293  : Real(0.0);
2294 
2295  if (qr_arr(i,j,k) > Real(qcrmin)) {
2296  const Real acrfac = Real(6.0) * rslope2_arr(i,j,k,0)
2297  + Real(4.0) * diameter * rslope_arr(i,j,k,0)
2298  + diameter * diameter;
2299  praci_arr(i,j,k) = pi_wdm6_loc * qi_val * nr_arr(i,j,k)
2300  * std::abs(vt2r - vt2i) * acrfac / Real(4.0);
2301  praci_arr(i,j,k) *= std::pow(
2302  amrex::min(amrex::max(Real(0.0), qr_arr(i,j,k) / qi_val), Real(1.0)),
2303  Real(2.0));
2304  praci_arr(i,j,k) = amrex::min(praci_arr(i,j,k), qi_val / dtcld);
2305 
2306  piacr_arr(i,j,k) = pi_wdm6_loc * pi_wdm6_loc * Real(WDM6::avtr)
2307  * nr_arr(i,j,k) * Real(denr) * xni_arr(i,j,k) * denfac_arr(i,j,k)
2308  * g7pbr_loc * rslope3_arr(i,j,k,0) * rslope2_arr(i,j,k,0)
2309  * rslopeb_arr(i,j,k,0) / (Real(24.0) * den_arr(i,j,k));
2310  piacr_arr(i,j,k) *= std::pow(
2311  amrex::min(amrex::max(Real(0.0), qi_val / qr_arr(i,j,k)), Real(1.0)),
2312  Real(2.0));
2313  piacr_arr(i,j,k) = amrex::min(piacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2314  }
2315 
2316  if (nr_arr(i,j,k) > Real(nrmin)) {
2317  niacr_arr(i,j,k) = pi_wdm6_loc * Real(WDM6::avtr) * nr_arr(i,j,k)
2318  * xni_arr(i,j,k) * denfac_arr(i,j,k) * g4pbr_loc
2319  * rslope2_arr(i,j,k,0) * rslopeb_arr(i,j,k,0) / Real(4.0);
2320  niacr_arr(i,j,k) *= std::pow(
2321  amrex::min(amrex::max(Real(0.0), qi_val / qr_arr(i,j,k)), Real(1.0)),
2322  Real(2.0));
2323  niacr_arr(i,j,k) = amrex::min(niacr_arr(i,j,k), nr_arr(i,j,k) / dtcld);
2324  }
2325 
2326  if (qs_arr(i,j,k) > Real(qcrmin)) {
2327  const Real acrfac = Real(2.0) * rslope3_arr(i,j,k,1)
2328  + Real(2.0) * diameter * rslope2_arr(i,j,k,1)
2329  + diameter * diameter * rslope_arr(i,j,k,1);
2330  psaci_arr(i,j,k) = pi_wdm6_loc * qi_val * eacrs * Real(n0s) * n0sfac
2331  * std::abs(vt2ave - vt2i) * acrfac / Real(4.0);
2332  psaci_arr(i,j,k) = amrex::min(psaci_arr(i,j,k), qi_val / dtcld);
2333 
2334  }
2335 
2336  if (qg_arr(i,j,k) > Real(qcrmin)) {
2337  // Fortran :2114, the same unsuffixed 0.07 as eacrs above.
2338  // pgaci reads bitwise-equal in the current column only
2339  // because the branch is inactive there; routed for the same
2340  // reason, not on separate evidence.
2341  const Real egi = std::exp(wdm6_literal(0.07) * (-supcol));
2342  const Real acrfac = Real(2.0) * rslope3_arr(i,j,k,2)
2343  + Real(2.0) * diameter * rslope2_arr(i,j,k,2)
2344  + diameter * diameter * rslope_arr(i,j,k,2);
2345  pgaci_arr(i,j,k) = pi_wdm6_loc * egi * qi_val * n0g_loc
2346  * std::abs(vt2ave - vt2i) * acrfac / Real(4.0);
2347  pgaci_arr(i,j,k) = amrex::min(pgaci_arr(i,j,k), qi_val / dtcld);
2348  }
2349  });
2350 
2351 
2352 
2353  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2354  const Real supcol = Real(t0c) - t_arr(i,j,k);
2355  const Real n0sfac = amrex::max(
2356  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
2357  Real(n0smax) / Real(n0s)),
2358  Real(1.0));
2359  const Real qs_val = qs_arr(i,j,k);
2360  const Real qg_val = qg_arr(i,j,k);
2361  const Real qc_val = qc_arr(i,j,k);
2362  const Real nc_val = nc_arr(i,j,k);
2363 
2364  // The Fortran (:2143, :2152) writes these ratios UNGUARDED:
2365  // min(max(0.0, qrs(i,k,2)/qci(i,k,1)), 1.)
2366  // so when qc == 0 and qs > 0 the division yields IEEE +Inf,
2367  // max(0,Inf) is Inf, and the clamp returns 1.0. Returning 0.0
2368  // for that case instead -- as a defensive divide-by-zero guard
2369  // naturally does -- lands on the OPPOSITE end of the clamp and
2370  // silently changes the physics.
2371  //
2372  // It hides in psacw, which is gated on qc > qmin and so never
2373  // sees qc == 0, but nsacw is gated on nc > ncmin, which passes
2374  // with qc == 0, and there the ratio is load-bearing. Measured
2375  // at (199,3,88) from a bitwise-identical PRE_G13C with qc = 0,
2376  // qs = 1.935e-02 and nc = 1.048e+06: the bridge produced
2377  // nsacw = 9.6929949510183156e+04 and the native leg exactly 0.
2378  // That single term is the entire step-6 nn discontinuity.
2379  //
2380  // The division is still avoided rather than allowed to produce
2381  // Inf, so builds with FP-exception trapping stay clean. The
2382  // qc == 0 and qs == 0 case cannot be reached by any consumer --
2383  // every use is gated on qs > qcrmin, qg > qcrmin or qc > qmin --
2384  // and 0.0 is returned there only for definiteness.
2385  const Real ratio_s = (qc_val > Real(0.0))
2386  ? amrex::min(amrex::max(Real(0.0), qs_val / qc_val), Real(1.0))
2387  : (qs_val > Real(0.0) ? Real(1.0) : Real(0.0));
2388  // Same treatment as ratio_s above, and for the same reason. This
2389  // was MISSED when ratio_s was fixed: that edit replaced only the
2390  // ratio_s else-branch while its comment claimed both, so ngacw
2391  // kept returning 0 where the Fortran returns 1.0. Surfaced at
2392  // step 8 on ngacw, at relative error exactly 1.0: bridge
2393  // 8.85033631487811001e-02 against native 0, which was the whole
2394  // remaining nn residual at 10 steps.
2395  const Real ratio_g = (qc_val > Real(0.0))
2396  ? amrex::min(amrex::max(Real(0.0), qg_val / qc_val), Real(1.0))
2397  : (qg_val > Real(0.0) ? Real(1.0) : Real(0.0));
2398 
2399  if (qs_val > Real(qcrmin) && qc_val > Real(qmin)) {
2400  psacw_arr(i,j,k) = amrex::min(
2401  pacrc_loc * n0sfac * rslope3_arr(i,j,k,1) * rslopeb_arr(i,j,k,1)
2402  * ratio_s * ratio_s * qc_val * denfac_arr(i,j,k),
2403  qc_val / dtcld);
2404  }
2405  if (qs_val > Real(qcrmin) && nc_val > Real(ncmin)) {
2406  nsacw_arr(i,j,k) = amrex::min(
2407  pacrc_loc * n0sfac * rslope3_arr(i,j,k,1) * rslopeb_arr(i,j,k,1)
2408  * ratio_s * ratio_s * nc_val * denfac_arr(i,j,k),
2409  nc_val / dtcld);
2410  }
2411  if (qg_val > Real(qcrmin) && qc_val > Real(qmin)) {
2412  pgacw_arr(i,j,k) = amrex::min(
2413  pacrg_loc * rslope3_arr(i,j,k,2) * rslopeb_arr(i,j,k,2)
2414  * qc_val * ratio_g * ratio_g * denfac_arr(i,j,k),
2415  qc_val / dtcld);
2416  }
2417  if (qg_val > Real(qcrmin) && nc_val > Real(ncmin)) {
2418  ngacw_arr(i,j,k) = amrex::min(
2419  pacrg_loc * rslope3_arr(i,j,k,2) * rslopeb_arr(i,j,k,2)
2420  * nc_val * ratio_g * ratio_g * denfac_arr(i,j,k),
2421  nc_val / dtcld);
2422  }
2423 
2424  const Real qsum = amrex::max(qs_val + qg_val, wdm6_literal(1.0e-15));
2425  if (qsum > wdm6_literal(1.0e-15)) {
2426  paacw_arr(i,j,k) = (qs_val * psacw_arr(i,j,k) + qg_val * pgacw_arr(i,j,k)) / qsum;
2427  naacw_arr(i,j,k) = (qs_val * nsacw_arr(i,j,k) + qg_val * ngacw_arr(i,j,k)) / qsum;
2428  }
2429  });
2430 
2431 
2432 
2433  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2434  const Real supcol = Real(t0c) - t_arr(i,j,k);
2435  const Real n0sfac = amrex::max(
2436  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
2437  Real(n0smax) / Real(n0s)),
2438  Real(1.0));
2439  const Real qr_val = qr_arr(i,j,k);
2440  const Real qs_val = qs_arr(i,j,k);
2441  const Real qg_val = qg_arr(i,j,k);
2442  const Real nr_val = nr_arr(i,j,k);
2443  const Real vt2r = pvtr_loc * rslopeb_arr(i,j,k,0) * denfac_arr(i,j,k);
2444  const Real vt2s = pvts_loc * rslopeb_arr(i,j,k,1) * denfac_arr(i,j,k);
2445  const Real vt2g = pvtg_loc * rslopeb_arr(i,j,k,2) * denfac_arr(i,j,k);
2446  const Real qsum = qs_val + qg_val;
2447  const Real vt2ave = (qsum > Real(0.0))
2448  ? (vt2s * qs_val + vt2g * qg_val) / qsum
2449  : Real(0.0);
2450 
2451  if (qs_val > Real(qcrmin) && qr_val > Real(qcrmin)) {
2452  if (supcol > Real(0.0)) {
2453  const Real acrfac =
2454  Real(5.0) * rslope3_arr(i,j,k,1) * rslope3_arr(i,j,k,1)
2455  + Real(4.0) * rslope3_arr(i,j,k,1) * rslope2_arr(i,j,k,1)
2456  * rslope_arr(i,j,k,0)
2457  + Real(1.5) * rslope2_arr(i,j,k,1) * rslope2_arr(i,j,k,1)
2458  * rslope2_arr(i,j,k,0);
2459  pracs_arr(i,j,k) = pi_wdm6_loc * pi_wdm6_loc * nr_val * Real(n0s)
2460  * n0sfac * std::abs(vt2r - vt2ave)
2461  * (Real(dens) / den_arr(i,j,k)) * acrfac;
2462  const Real ratio = amrex::min(
2463  amrex::max(Real(0.0), qr_val / qs_val), Real(1.0));
2464  pracs_arr(i,j,k) *= ratio * ratio;
2465  pracs_arr(i,j,k) = amrex::min(pracs_arr(i,j,k), qs_val / dtcld);
2466  }
2467 
2468  const Real acrfac =
2469  Real(30.0) * rslope3_arr(i,j,k,0) * rslope2_arr(i,j,k,0)
2470  * rslope_arr(i,j,k,1)
2471  + Real(10.0) * rslope2_arr(i,j,k,0) * rslope2_arr(i,j,k,0)
2472  * rslope2_arr(i,j,k,1)
2473  + Real(2.0) * rslope3_arr(i,j,k,0) * rslope3_arr(i,j,k,1);
2474  psacr_arr(i,j,k) = pi_wdm6_loc * pi_wdm6_loc * nr_val * Real(n0s)
2475  * n0sfac * std::abs(vt2ave - vt2r)
2476  * (Real(denr) / den_arr(i,j,k)) * acrfac;
2477  const Real ratio = amrex::min(
2478  amrex::max(Real(0.0), qs_val / qr_val), Real(1.0));
2479  psacr_arr(i,j,k) *= ratio * ratio;
2480  psacr_arr(i,j,k) = amrex::min(psacr_arr(i,j,k), qr_val / dtcld);
2481  }
2482 
2483  if (qs_val > Real(qcrmin) && nr_val > Real(nrmin)) {
2484  const Real acrfac =
2485  Real(1.5) * rslope2_arr(i,j,k,0) * rslope_arr(i,j,k,1)
2486  + rslope_arr(i,j,k,0) * rslope2_arr(i,j,k,1)
2487  + Real(0.5) * rslope3_arr(i,j,k,1);
2488  nsacr_arr(i,j,k) = pi_wdm6_loc * nr_val * Real(n0s) * n0sfac
2489  * std::abs(vt2ave - vt2r) * acrfac;
2490  const Real ratio = amrex::min(
2491  amrex::max(Real(0.0), qs_val / qr_val), Real(1.0));
2492  nsacr_arr(i,j,k) *= ratio * ratio;
2493  nsacr_arr(i,j,k) = amrex::min(nsacr_arr(i,j,k), nr_val / dtcld);
2494  }
2495 
2496  if (qg_val > Real(qcrmin) && qr_val > Real(qcrmin)) {
2497  const Real acrfac =
2498  Real(30.0) * rslope3_arr(i,j,k,0) * rslope2_arr(i,j,k,0)
2499  * rslope_arr(i,j,k,2)
2500  + Real(10.0) * rslope2_arr(i,j,k,0) * rslope2_arr(i,j,k,0)
2501  * rslope2_arr(i,j,k,2)
2502  + Real(2.0) * rslope3_arr(i,j,k,0) * rslope3_arr(i,j,k,2);
2503  pgacr_arr(i,j,k) = pi_wdm6_loc * pi_wdm6_loc * nr_val * n0g_loc
2504  * std::abs(vt2ave - vt2r) * (Real(denr) / den_arr(i,j,k))
2505  * acrfac;
2506  const Real ratio = amrex::min(
2507  amrex::max(Real(0.0), qg_val / qr_val), Real(1.0));
2508  pgacr_arr(i,j,k) *= ratio * ratio;
2509  pgacr_arr(i,j,k) = amrex::min(pgacr_arr(i,j,k), qr_val / dtcld);
2510  }
2511 
2512  if (qg_val > Real(qcrmin) && nr_val > Real(nrmin)) {
2513  const Real acrfac =
2514  Real(1.5) * rslope2_arr(i,j,k,0) * rslope_arr(i,j,k,2)
2515  + rslope_arr(i,j,k,0) * rslope2_arr(i,j,k,2)
2516  + Real(0.5) * rslope3_arr(i,j,k,2);
2517  ngacr_arr(i,j,k) = pi_wdm6_loc * nr_val * n0g_loc
2518  * std::abs(vt2ave - vt2r) * acrfac;
2519  const Real ratio = amrex::min(
2520  amrex::max(Real(0.0), qg_val / qr_val), Real(1.0));
2521  ngacr_arr(i,j,k) *= ratio * ratio;
2522  ngacr_arr(i,j,k) = amrex::min(ngacr_arr(i,j,k), nr_val / dtcld);
2523  }
2524 
2525  if (qg_val > Real(qcrmin) && qs_val > Real(qcrmin)) {
2526  pgacs_arr(i,j,k) = Real(0.0);
2527  }
2528 
2529  if (supcol <= Real(0.0)) {
2530  const Real xlf = Real(xlf0);
2531  if (qs_val > Real(0.0)) {
2532  pseml_arr(i,j,k) = amrex::min(
2533  amrex::max(
2534  Real(cliq) * supcol
2535  * (paacw_arr(i,j,k) + psacr_arr(i,j,k)) / xlf,
2536  -qs_val / dtcld),
2537  Real(0.0));
2538  }
2539  if (qs_val > Real(qcrmin)) {
2540  const Real sfac = rslope_arr(i,j,k,1) * Real(n0s) * n0sfac / qs_val;
2541  nseml_arr(i,j,k) = -sfac * pseml_arr(i,j,k);
2542  }
2543 
2544  if (qg_val > Real(0.0)) {
2545  pgeml_arr(i,j,k) = amrex::min(
2546  amrex::max(
2547  Real(cliq) * supcol
2548  * (paacw_arr(i,j,k) + pgacr_arr(i,j,k)) / xlf,
2549  -qg_val / dtcld),
2550  Real(0.0));
2551  }
2552  if (qg_val > Real(qcrmin)) {
2553  const Real gfac = rslope_arr(i,j,k,2) * n0g_loc / qg_val;
2554  ngeml_arr(i,j,k) = -gfac * pgeml_arr(i,j,k);
2555  }
2556  }
2557  });
2558 
2559 
2560 
2561  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2562  const Real supcol = Real(t0c) - t_arr(i,j,k);
2563  if (supcol <= Real(0.0)) {
2564  return;
2565  }
2566 
2567  const Real n0sfac = amrex::max(
2568  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
2569  Real(n0smax) / Real(n0s)),
2570  Real(1.0));
2571  const Real supsat =
2572  amrex::max(qv_arr(i,j,k), Real(qmin)) - qsati_arr(i,j,k);
2573  const Real satdt = supsat / dtcld;
2574  int ifsat = 0;
2575 
2576  const Real qi_val = qi_arr(i,j,k);
2577  const Real qs_val = qs_arr(i,j,k);
2578  const Real qg_val = qg_arr(i,j,k);
2579  const Real xni = xni_arr(i,j,k);
2580  const Real xni_safe = amrex::max(xni, Real(1.0e-30));
2581  const Real xmi = den_arr(i,j,k) * qi_val / xni_safe;
2582  const Real diameter = amrex::min(
2583  Real(dicon) * std::sqrt(xmi), Real(dimax));
2584  const Real rhi = rhi_arr(i,j,k);
2585  const Real work1i = work1_arr(i,j,k,1);
2586 
2587  if (qi_val > Real(0.0) && ifsat != 1) {
2588  pidep_arr(i,j,k) = Real(4.0) * diameter * xni
2589  * (rhi - Real(1.0)) / work1i;
2590  Real supice = satdt - prevp_arr(i,j,k);
2591  if (pidep_arr(i,j,k) < Real(0.0)) {
2592  pidep_arr(i,j,k) = amrex::max(
2593  amrex::max(pidep_arr(i,j,k), satdt / Real(2.0)),
2594  supice);
2595  pidep_arr(i,j,k) = amrex::max(
2596  pidep_arr(i,j,k), -qi_val / dtcld);
2597  } else {
2598  pidep_arr(i,j,k) = amrex::min(
2599  amrex::min(pidep_arr(i,j,k), satdt / Real(2.0)),
2600  supice);
2601  }
2602  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k))
2603  >= std::abs(satdt)) {
2604  ifsat = 1;
2605  }
2606  }
2607 
2608  if (qs_val > Real(0.0) && ifsat != 1) {
2609  const Real coeres_s = rslope2_arr(i,j,k,1)
2610  * std::sqrt(rslope_arr(i,j,k,1) * rslopeb_arr(i,j,k,1));
2611  psdep_arr(i,j,k) = (rhi - Real(1.0)) * n0sfac
2612  * (precs1_loc * rslope2_arr(i,j,k,1)
2613  + precs2_loc * work2_arr(i,j,k) * coeres_s)
2614  / work1i;
2615  Real supice = satdt - prevp_arr(i,j,k) - pidep_arr(i,j,k);
2616  if (psdep_arr(i,j,k) < Real(0.0)) {
2617  psdep_arr(i,j,k) = amrex::max(
2618  psdep_arr(i,j,k), -qs_val / dtcld);
2619  psdep_arr(i,j,k) = amrex::max(
2620  amrex::max(psdep_arr(i,j,k), satdt / Real(2.0)),
2621  supice);
2622  } else {
2623  psdep_arr(i,j,k) = amrex::min(
2624  amrex::min(psdep_arr(i,j,k), satdt / Real(2.0)),
2625  supice);
2626  }
2627  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2628  + psdep_arr(i,j,k)) >= std::abs(satdt)) {
2629  ifsat = 1;
2630  }
2631  }
2632 
2633  if (qg_val > Real(0.0) && ifsat != 1) {
2634  const Real coeres_g = rslope2_arr(i,j,k,2)
2635  * std::sqrt(rslope_arr(i,j,k,2) * rslopeb_arr(i,j,k,2));
2636  pgdep_arr(i,j,k) = (rhi - Real(1.0))
2637  * (precg1_loc * rslope2_arr(i,j,k,2)
2638  + precg2_loc * work2_arr(i,j,k) * coeres_g)
2639  / work1i;
2640  Real supice = satdt - prevp_arr(i,j,k) - pidep_arr(i,j,k)
2641  - psdep_arr(i,j,k);
2642  if (pgdep_arr(i,j,k) < Real(0.0)) {
2643  pgdep_arr(i,j,k) = amrex::max(
2644  pgdep_arr(i,j,k), -qg_val / dtcld);
2645  pgdep_arr(i,j,k) = amrex::max(
2646  amrex::max(pgdep_arr(i,j,k), satdt / Real(2.0)),
2647  supice);
2648  } else {
2649  pgdep_arr(i,j,k) = amrex::min(
2650  amrex::min(pgdep_arr(i,j,k), satdt / Real(2.0)),
2651  supice);
2652  }
2653  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2654  + psdep_arr(i,j,k) + pgdep_arr(i,j,k))
2655  >= std::abs(satdt)) {
2656  ifsat = 1;
2657  }
2658  }
2659 
2660  if (supsat > Real(0.0) && ifsat != 1) {
2661  const Real supice = satdt - prevp_arr(i,j,k) - pidep_arr(i,j,k)
2662  - psdep_arr(i,j,k) - pgdep_arr(i,j,k);
2663  const Real xni0 = wdm6_literal(1.0e3) * std::exp(wdm6_literal(0.1) * supcol);
2664  const Real roqi0 = wdm6_literal(4.92e-11) * std::pow(xni0, wdm6_literal(1.33));
2665  const Real pigen_raw = amrex::max(
2666  Real(0.0),
2667  (roqi0 / den_arr(i,j,k) - amrex::max(qi_val, Real(0.0))) / dtcld);
2668  pigen_arr(i,j,k) = amrex::min(
2669  amrex::min(pigen_raw, satdt), supice);
2670  }
2671  });
2672 
2673 
2674 
2675  const Real roqimax_loc = m_roqimax;
2676  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2677  const Real supcol = Real(t0c) - t_arr(i,j,k);
2678 
2679  if (qi_arr(i,j,k) > Real(0.0)) {
2680  const Real qimax = roqimax_loc / den_arr(i,j,k);
2681  psaut_arr(i,j,k) = amrex::max(
2682  Real(0.0), (qi_arr(i,j,k) - qimax) / dtcld);
2683  }
2684 
2685  if (qs_arr(i,j,k) > Real(0.0)) {
2686  // Fortran :2510 alpha2 = 1.e-3*exp(0.09*(-supcol)). Both
2687  // literals are unsuffixed and they push OPPOSITE ways, so
2688  // neither can be judged alone: 1.e-3 contributes +4.749745e-08
2689  // and 0.09 contributes -3.5763e-09*supcol, and the residual
2690  // is their signed sum. Verified with zero free parameters
2691  // against measured supcol at (199,3), k=91..98: predicted
2692  // and measured relative error agree to 7 digits at every
2693  // cell, ratio 1.0000. At k=91 that value is 7.652614e-08,
2694  // which is exactly the step-2 rhoQ6/qgraup plotfile
2695  // residual, so this literal pair is that residual.
2696  const Real alpha2 = wdm6_literal(1.e-3)
2697  * std::exp(wdm6_literal(0.09) * (-supcol));
2698  pgaut_arr(i,j,k) = amrex::min(
2699  amrex::max(
2700  Real(0.0), alpha2 * (qs_arr(i,j,k) - Real(qs0))),
2701  qs_arr(i,j,k) / dtcld);
2702  }
2703  });
2704 
2705 
2706 
2707  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2708  const Real supcol = Real(t0c) - t_arr(i,j,k);
2709  if (supcol < Real(0.0)) {
2710  const Real n0sfac = amrex::max(
2711  amrex::min(std::exp(Real(alpha_wdm6) * supcol),
2712  Real(n0smax) / Real(n0s)),
2713  Real(1.0));
2714 
2715  if (qs_arr(i,j,k) > Real(0.0)
2716  && rhw_arr(i,j,k) < Real(1.0)) {
2717  const Real coeres = rslope2_arr(i,j,k,1)
2718  * std::sqrt(rslope_arr(i,j,k,1)
2719  * rslopeb_arr(i,j,k,1));
2720  psevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2721  * n0sfac
2722  * (precs1_loc * rslope2_arr(i,j,k,1)
2723  + precs2_loc * work2_arr(i,j,k) * coeres)
2724  / work1_arr(i,j,k,0);
2725  psevp_arr(i,j,k) = amrex::min(
2726  amrex::max(psevp_arr(i,j,k),
2727  -qs_arr(i,j,k) / dtcld),
2728  Real(0.0));
2729  }
2730 
2731  if (qg_arr(i,j,k) > Real(0.0)
2732  && rhw_arr(i,j,k) < Real(1.0)) {
2733  const Real coeres = rslope2_arr(i,j,k,2)
2734  * std::sqrt(rslope_arr(i,j,k,2)
2735  * rslopeb_arr(i,j,k,2));
2736  pgevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2737  * (precg1_loc * rslope2_arr(i,j,k,2)
2738  + precg2_loc * work2_arr(i,j,k) * coeres)
2739  / work1_arr(i,j,k,0);
2740  pgevp_arr(i,j,k) = amrex::min(
2741  amrex::max(pgevp_arr(i,j,k),
2742  -qg_arr(i,j,k) / dtcld),
2743  Real(0.0));
2744  }
2745  }
2746  });
2747 
2748 
2749 
2750  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2751  const Real qmin_l = Real(qmin);
2752  const Real qcrmin_l = Real(qcrmin);
2753  const Real ncmin_l = Real(ncmin);
2754  const Real nrmin_l = Real(nrmin);
2755  const Real t0c_l = Real(t0c);
2756  const Real one_l = Real(1.0);
2757  const Real zero_l = Real(0.0);
2758  const Real delta2 =
2759  (qr_arr(i,j,k) < wdm6_literal(1.0e-4) && qs_arr(i,j,k) < wdm6_literal(1.0e-4))
2760  ? one_l : zero_l;
2761  const Real delta3 =
2762  (qr_arr(i,j,k) < wdm6_literal(1.0e-4)) ? one_l : zero_l;
2763 
2764  if (t_arr(i,j,k) <= t0c_l) {
2765  Real value, source, factor, xlf, xlwork2;
2766 
2767  value = amrex::max(qmin_l, qc_arr(i,j,k));
2768  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2769  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2770  if (source > value) {
2771  factor = value / source;
2772  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2773  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2774  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2775  }
2776 
2777  value = amrex::max(qmin_l, qi_arr(i,j,k));
2778  source = (psaut_arr(i,j,k) - pigen_arr(i,j,k)
2779  - pidep_arr(i,j,k) + praci_arr(i,j,k)
2780  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)) * dtcld;
2781  if (source > value) {
2782  factor = value / source;
2783  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2784  pigen_arr(i,j,k) = pigen_arr(i,j,k) * factor;
2785  pidep_arr(i,j,k) = pidep_arr(i,j,k) * factor;
2786  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2787  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2788  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2789  }
2790 
2791  value = amrex::max(qmin_l, qr_arr(i,j,k));
2792  source = (-praut_arr(i,j,k) - prevp_arr(i,j,k)
2793  - pracw_arr(i,j,k) + piacr_arr(i,j,k)
2794  + psacr_arr(i,j,k) + pgacr_arr(i,j,k)) * dtcld;
2795  if (source > value) {
2796  factor = value / source;
2797  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2798  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2799  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2800  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2801  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2802  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2803  }
2804 
2805  value = amrex::max(qmin_l, qs_arr(i,j,k));
2806  source = -(psdep_arr(i,j,k) + psaut_arr(i,j,k)
2807  - pgaut_arr(i,j,k) + paacw_arr(i,j,k)
2808  + piacr_arr(i,j,k) * delta3
2809  + praci_arr(i,j,k) * delta3
2810  - pracs_arr(i,j,k) * (one_l - delta2)
2811  + psacr_arr(i,j,k) * delta2
2812  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)) * dtcld;
2813  if (source > value) {
2814  factor = value / source;
2815  psdep_arr(i,j,k) = psdep_arr(i,j,k) * factor;
2816  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2817  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2818  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2819  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2820  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2821  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2822  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2823  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2824  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2825  }
2826 
2827  value = amrex::max(qmin_l, qg_arr(i,j,k));
2828  source = -(pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2829  + piacr_arr(i,j,k) * (one_l - delta3)
2830  + praci_arr(i,j,k) * (one_l - delta3)
2831  + psacr_arr(i,j,k) * (one_l - delta2)
2832  + pracs_arr(i,j,k) * (one_l - delta2)
2833  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2834  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k)) * dtcld;
2835  if (source > value) {
2836  factor = value / source;
2837  pgdep_arr(i,j,k) = pgdep_arr(i,j,k) * factor;
2838  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2839  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2840  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2841  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2842  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2843  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2844  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2845  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2846  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2847  }
2848 
2849  value = amrex::max(ncmin_l, nc_arr(i,j,k));
2850  source = (nrauto_arr(i,j,k) + nccol_arr(i,j,k)
2851  + nraccr_arr(i,j,k) + naacw_arr(i,j,k)
2852  + naacw_arr(i,j,k)) * dtcld;
2853  if (source > value) {
2854  factor = value / source;
2855  nrauto_arr(i,j,k) = nrauto_arr(i,j,k) * factor;
2856  nccol_arr(i,j,k) = nccol_arr(i,j,k) * factor;
2857  nraccr_arr(i,j,k) = nraccr_arr(i,j,k) * factor;
2858  naacw_arr(i,j,k) = naacw_arr(i,j,k) * factor;
2859  }
2860 
2861  value = amrex::max(nrmin_l, nr_arr(i,j,k));
2862  source = (-nrauto_arr(i,j,k) + nrcol_arr(i,j,k)
2863  + niacr_arr(i,j,k) + nsacr_arr(i,j,k)
2864  + ngacr_arr(i,j,k)) * dtcld;
2865  if (source > value) {
2866  factor = value / source;
2867  nrauto_arr(i,j,k) = nrauto_arr(i,j,k) * factor;
2868  nrcol_arr(i,j,k) = nrcol_arr(i,j,k) * factor;
2869  niacr_arr(i,j,k) = niacr_arr(i,j,k) * factor;
2870  nsacr_arr(i,j,k) = nsacr_arr(i,j,k) * factor;
2871  ngacr_arr(i,j,k) = ngacr_arr(i,j,k) * factor;
2872  }
2873 
2874  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psdep_arr(i,j,k)
2875  + pgdep_arr(i,j,k) + pigen_arr(i,j,k)
2876  + pidep_arr(i,j,k));
2877  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
2878  qc_arr(i,j,k) = amrex::max(
2879  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
2880  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
2881  * dtcld,
2882  zero_l);
2883  qr_arr(i,j,k) = amrex::max(
2884  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
2885  + prevp_arr(i,j,k) - piacr_arr(i,j,k)
2886  - pgacr_arr(i,j,k) - psacr_arr(i,j,k))
2887  * dtcld,
2888  zero_l);
2889  qi_arr(i,j,k) = amrex::max(
2890  qi_arr(i,j,k) - (psaut_arr(i,j,k) + praci_arr(i,j,k)
2891  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)
2892  - pigen_arr(i,j,k) - pidep_arr(i,j,k))
2893  * dtcld,
2894  zero_l);
2895  qs_arr(i,j,k) = amrex::max(
2896  qs_arr(i,j,k) + (psdep_arr(i,j,k) + psaut_arr(i,j,k)
2897  + paacw_arr(i,j,k) - pgaut_arr(i,j,k)
2898  + piacr_arr(i,j,k) * delta3
2899  + praci_arr(i,j,k) * delta3
2900  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)
2901  - pracs_arr(i,j,k) * (one_l - delta2)
2902  + psacr_arr(i,j,k) * delta2) * dtcld,
2903  zero_l);
2904  qg_arr(i,j,k) = amrex::max(
2905  qg_arr(i,j,k) + (pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2906  + piacr_arr(i,j,k) * (one_l - delta3)
2907  + praci_arr(i,j,k) * (one_l - delta3)
2908  + psacr_arr(i,j,k) * (one_l - delta2)
2909  + pracs_arr(i,j,k) * (one_l - delta2)
2910  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2911  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k))
2912  * dtcld,
2913  zero_l);
2914  nc_arr(i,j,k) = amrex::max(
2915  nc_arr(i,j,k) + (-nrauto_arr(i,j,k) - nccol_arr(i,j,k)
2916  - nraccr_arr(i,j,k) - naacw_arr(i,j,k)
2917  - naacw_arr(i,j,k)) * dtcld,
2918  zero_l);
2919  nr_arr(i,j,k) = amrex::max(
2920  nr_arr(i,j,k) + (nrauto_arr(i,j,k) - nrcol_arr(i,j,k)
2921  - niacr_arr(i,j,k) - nsacr_arr(i,j,k)
2922  - ngacr_arr(i,j,k)) * dtcld,
2923  zero_l);
2924  xlf = Real(xls) - xl_arr(i,j,k);
2925  xlwork2 = -Real(xls) * (psdep_arr(i,j,k) + pgdep_arr(i,j,k)
2926  + pidep_arr(i,j,k) + pigen_arr(i,j,k))
2927  - xl_arr(i,j,k) * prevp_arr(i,j,k)
2928  - xlf * (piacr_arr(i,j,k) + paacw_arr(i,j,k)
2929  + paacw_arr(i,j,k) + pgacr_arr(i,j,k)
2930  + psacr_arr(i,j,k));
2931  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
2932  } else {
2933  Real value, source, factor, xlf, xlwork2;
2934 
2935  value = amrex::max(qmin_l, qc_arr(i,j,k));
2936  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2937  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2938  if (source > value) {
2939  factor = value / source;
2940  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2941  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2942  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2943  }
2944 
2945  value = amrex::max(qmin_l, qr_arr(i,j,k));
2946  source = (-paacw_arr(i,j,k) - praut_arr(i,j,k)
2947  + pseml_arr(i,j,k) + pgeml_arr(i,j,k)
2948  - pracw_arr(i,j,k) - paacw_arr(i,j,k)
2949  - prevp_arr(i,j,k)) * dtcld;
2950  if (source > value) {
2951  factor = value / source;
2952  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2953  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2954  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2955  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2956  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2957  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2958  }
2959 
2960  value = amrex::max(qcrmin_l, qs_arr(i,j,k));
2961  source = (pgacs_arr(i,j,k) - pseml_arr(i,j,k)
2962  - psevp_arr(i,j,k)) * dtcld;
2963  if (source > value) {
2964  factor = value / source;
2965  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2966  psevp_arr(i,j,k) = psevp_arr(i,j,k) * factor;
2967  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2968  }
2969 
2970  value = amrex::max(qcrmin_l, qg_arr(i,j,k));
2971  source = -(pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
2972  + pgeml_arr(i,j,k)) * dtcld;
2973  if (source > value) {
2974  factor = value / source;
2975  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2976  pgevp_arr(i,j,k) = pgevp_arr(i,j,k) * factor;
2977  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2978  }
2979 
2980  value = amrex::max(ncmin_l, nc_arr(i,j,k));
2981  source = (nrauto_arr(i,j,k) + nccol_arr(i,j,k)
2982  + nraccr_arr(i,j,k) + naacw_arr(i,j,k)
2983  + naacw_arr(i,j,k)) * dtcld;
2984  if (source > value) {
2985  factor = value / source;
2986  nrauto_arr(i,j,k) = nrauto_arr(i,j,k) * factor;
2987  nccol_arr(i,j,k) = nccol_arr(i,j,k) * factor;
2988  nraccr_arr(i,j,k) = nraccr_arr(i,j,k) * factor;
2989  naacw_arr(i,j,k) = naacw_arr(i,j,k) * factor;
2990  }
2991 
2992  value = amrex::max(nrmin_l, nr_arr(i,j,k));
2993  source = (-nrauto_arr(i,j,k) + nrcol_arr(i,j,k)
2994  - nseml_arr(i,j,k) - ngeml_arr(i,j,k)) * dtcld;
2995  if (source > value) {
2996  factor = value / source;
2997  nrauto_arr(i,j,k) = nrauto_arr(i,j,k) * factor;
2998  nrcol_arr(i,j,k) = nrcol_arr(i,j,k) * factor;
2999  nseml_arr(i,j,k) = nseml_arr(i,j,k) * factor;
3000  ngeml_arr(i,j,k) = ngeml_arr(i,j,k) * factor;
3001  }
3002 
3003  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psevp_arr(i,j,k)
3004  + pgevp_arr(i,j,k));
3005  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
3006  qc_arr(i,j,k) = amrex::max(
3007  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
3008  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
3009  * dtcld,
3010  zero_l);
3011  qr_arr(i,j,k) = amrex::max(
3012  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
3013  + prevp_arr(i,j,k) + paacw_arr(i,j,k)
3014  + paacw_arr(i,j,k) - pseml_arr(i,j,k)
3015  - pgeml_arr(i,j,k)) * dtcld,
3016  zero_l);
3017  qs_arr(i,j,k) = amrex::max(
3018  qs_arr(i,j,k) + (psevp_arr(i,j,k) - pgacs_arr(i,j,k)
3019  + pseml_arr(i,j,k)) * dtcld,
3020  zero_l);
3021  qg_arr(i,j,k) = amrex::max(
3022  qg_arr(i,j,k) + (pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
3023  + pgeml_arr(i,j,k)) * dtcld,
3024  zero_l);
3025  nc_arr(i,j,k) = amrex::max(
3026  nc_arr(i,j,k) + (-nrauto_arr(i,j,k) - nccol_arr(i,j,k)
3027  - nraccr_arr(i,j,k) - naacw_arr(i,j,k)
3028  - naacw_arr(i,j,k)) * dtcld,
3029  zero_l);
3030  nr_arr(i,j,k) = amrex::max(
3031  nr_arr(i,j,k) + (nrauto_arr(i,j,k) - nrcol_arr(i,j,k)
3032  + nseml_arr(i,j,k) + ngeml_arr(i,j,k))
3033  * dtcld,
3034  zero_l);
3035  xlf = Real(xls) - xl_arr(i,j,k);
3036  xlwork2 = -xl_arr(i,j,k) * (prevp_arr(i,j,k)
3037  + psevp_arr(i,j,k)
3038  + pgevp_arr(i,j,k))
3039  - xlf * (pseml_arr(i,j,k) + pgeml_arr(i,j,k));
3040  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
3041  }
3042  });
3043 
3044 
3045 
3046  const Real hsub = Real(xls);
3047  const Real hvap = Real(xlv0);
3048  const Real cvap = Real(cpv);
3049  const Real ttp = Real(t0c) + wdm6_literal(0.01);
3050  const Real dldt = cvap - Real(cliq);
3051  const Real xa = -dldt / Real(rv);
3052  const Real xb = xa + hvap / (Real(rv) * ttp);
3053  const Real dldti = cvap - Real(cice);
3054  const Real xai = -dldti / Real(rv);
3055  const Real xbi = xai + hsub / (Real(rv) * ttp);
3056 
3057  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3058  const Real tr = ttp / t_arr(i,j,k);
3059  Real qsw = Real(psat) * std::exp(std::log(tr) * xa)
3060  * std::exp(xb * (Real(1.) - tr));
3061  qsw = amrex::min(qsw, wdm6_literal(0.99) * p_arr(i,j,k));
3062  qsw = Real(ep2) * qsw / (p_arr(i,j,k) - qsw);
3063  qsw = amrex::max(qsw, Real(qmin));
3064  qsatw_arr(i,j,k) = qsw;
3065 
3066  Real qsi;
3067  if (t_arr(i,j,k) < ttp) {
3068  qsi = Real(psat) * std::exp(std::log(tr) * xai)
3069  * std::exp(xbi * (Real(1.) - tr));
3070  } else {
3071  qsi = Real(psat) * std::exp(std::log(tr) * xa)
3072  * std::exp(xb * (Real(1.) - tr));
3073  }
3074  qsi = amrex::min(qsi, wdm6_literal(0.99) * p_arr(i,j,k));
3075  qsi = Real(ep2) * qsi / (p_arr(i,j,k) - qsi);
3076  qsi = amrex::max(qsi, Real(qmin));
3077  qsati_arr(i,j,k) = qsi;
3078 
3079  rhw_arr(i,j,k) = amrex::max(qv_arr(i,j,k) / qsw, Real(qmin));
3080  });
3081 
3082 
3083  const Real g16a_cbrt24 = wdm6_default_real_pow(24.0, 0.3333333);
3084 
3085 
3086  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3087  qrs_tmp_arr(i,j,k,0) = qr_arr(i,j,k);
3088  qrs_tmp_arr(i,j,k,1) = qs_arr(i,j,k);
3089  qrs_tmp_arr(i,j,k,2) = qg_arr(i,j,k);
3090  ncr_tmp_arr(i,j,k) = nr_arr(i,j,k);
3091  });
3092 
3093  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3094  Real rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3, rain_vt, rain_vtn;
3095  Real snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3, snow_vt, snow_n0sfac;
3096  Real graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3, graup_vt;
3097 
3098  wdm6_slope_rain_cell(qrs_tmp_arr(i,j,k,0), ncr_tmp_arr(i,j,k),
3099  den_arr(i,j,k), denfac_arr(i,j,k),
3100  Real(qcrmin), Real(nrmin),
3101  rslopermax_loc, rsloperbmax_loc,
3102  rsloper2max_loc, rsloper3max_loc,
3103  Real(bvtr), pvtr_loc, pvtrn_loc, pidnr_loc,
3104  rain_rslope, rain_rslopeb, rain_rslope2, rain_rslope3,
3105  rain_vt, rain_vtn);
3106  wdm6_slope_snow_cell(qrs_tmp_arr(i,j,k,1), den_arr(i,j,k), denfac_arr(i,j,k),
3107  t_arr(i,j,k), pidn0s_loc, Real(alpha_wdm6),
3108  Real(n0smax), Real(n0s), wdm6_slope_t0c, Real(qcrmin),
3109  rslopesmax_loc, rslopesbmax_loc,
3110  rslopes2max_loc, rslopes3max_loc,
3111  Real(bvts), pvts_loc,
3112  snow_rslope, snow_rslopeb, snow_rslope2, snow_rslope3,
3113  snow_vt, snow_n0sfac);
3114  wdm6_slope_graup_cell(qrs_tmp_arr(i,j,k,2), den_arr(i,j,k), denfac_arr(i,j,k),
3115  pidn0g_loc, Real(qcrmin),
3116  rslopegmax_loc, rslopegbmax_loc,
3117  rslopeg2max_loc, rslopeg3max_loc,
3118  slope_bvtg_loc, pvtg_loc,
3119  graup_rslope, graup_rslopeb, graup_rslope2, graup_rslope3,
3120  graup_vt);
3121 
3122  rslope_arr(i,j,k,0) = rain_rslope;
3123  rslope_arr(i,j,k,1) = snow_rslope;
3124  rslope_arr(i,j,k,2) = graup_rslope;
3125  rslopeb_arr(i,j,k,0) = rain_rslopeb;
3126  rslopeb_arr(i,j,k,1) = snow_rslopeb;
3127  rslopeb_arr(i,j,k,2) = graup_rslopeb;
3128  rslope2_arr(i,j,k,0) = rain_rslope2;
3129  rslope2_arr(i,j,k,1) = snow_rslope2;
3130  rslope2_arr(i,j,k,2) = graup_rslope2;
3131  rslope3_arr(i,j,k,0) = rain_rslope3;
3132  rslope3_arr(i,j,k,1) = snow_rslope3;
3133  rslope3_arr(i,j,k,2) = graup_rslope3;
3134  work1_arr(i,j,k,0) = rain_vt;
3135  work1_arr(i,j,k,1) = snow_vt;
3136  work1_arr(i,j,k,2) = graup_vt;
3137  workn_arr(i,j,k) = rain_vtn;
3138 
3139  avedia_arr(i,j,k,1) = rslope_arr(i,j,k,0) * g16a_cbrt24;
3140  if (avedia_arr(i,j,k,1) <= Real(di82)) {
3141  nc_arr(i,j,k) += nr_arr(i,j,k);
3142  nr_arr(i,j,k) = Real(0.0);
3143  qc_arr(i,j,k) += qr_arr(i,j,k);
3144  qr_arr(i,j,k) = Real(0.0);
3145  }
3146  });
3147 
3148 
3149 
3150  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3151  if (rhw_arr(i,j,k) > Real(1.0)) {
3152  const Real ratio = rhw_arr(i,j,k) / Real(satmax);
3153  const Real fraction = amrex::min(Real(1.0),
3154  std::exp(std::log(ratio) * Real(actk)));
3155  Real ncact_raw = (nn_arr(i,j,k) + nc_arr(i,j,k)) * fraction - nc_arr(i,j,k);
3156  Real ncact = amrex::max(Real(0.0), ncact_raw);
3157  ncact /= dtcld;
3158  const Real ncact_cap = amrex::max(nn_arr(i,j,k), Real(0.0)) / dtcld;
3159  ncact = amrex::min(ncact, ncact_cap);
3160  const Real actr_um = Real(actr) * wdm6_literal(1.0e-6);
3161  const Real pcact = amrex::min(
3162  Real(4.0) * pi_wdm6_loc * Real(denr)
3163  * actr_um * actr_um * actr_um * ncact
3164  / (Real(3.0) * den_arr(i,j,k)),
3165  amrex::max(qv_arr(i,j,k), Real(0.0)) / dtcld);
3166 
3167  ncact_arr(i,j,k) = ncact;
3168  act_ratio_arr(i,j,k) = ratio;
3169  act_fraction_arr(i,j,k) = fraction;
3170  act_raw_arr(i,j,k) = ncact_raw;
3171  act_cap_arr(i,j,k) = ncact_cap;
3172  pcact_arr(i,j,k) = pcact;
3173  qv_arr(i,j,k) = amrex::max(qv_arr(i,j,k) - pcact * dtcld, Real(0.0));
3174  qc_arr(i,j,k) = amrex::max(qc_arr(i,j,k) + pcact * dtcld, Real(0.0));
3175  nn_arr(i,j,k) = amrex::max(nn_arr(i,j,k) - ncact * dtcld, Real(0.0));
3176  nc_arr(i,j,k) = amrex::max(nc_arr(i,j,k) + ncact * dtcld, Real(0.0));
3177  t_arr(i,j,k) += pcact * xl_arr(i,j,k) / cpm_arr(i,j,k) * dtcld;
3178  }
3179 
3180  const Real tr = ttp / t_arr(i,j,k);
3181  Real qsw = Real(psat) * std::exp(std::log(tr) * xa)
3182  * std::exp(xb * (Real(1.0) - tr));
3183  qsw = amrex::min(qsw, wdm6_literal(0.99) * p_arr(i,j,k));
3184  qsw = Real(ep2) * qsw / (p_arr(i,j,k) - qsw);
3185  qsw = amrex::max(qsw, Real(qmin));
3186  qsatw_arr(i,j,k) = qsw;
3187 
3188  work1_arr(i,j,k,0) = wdm6_conden(
3189  t_arr(i,j,k), qv_arr(i,j,k), qsw, xl_arr(i,j,k), cpm_arr(i,j,k),
3190  Real(qmin), Real(rv));
3191  work2_arr(i,j,k) = qc_arr(i,j,k) + work1_arr(i,j,k,0);
3192 
3193  Real pcond = amrex::min(
3194  amrex::max(work1_arr(i,j,k,0) / dtcld, Real(0.0)),
3195  amrex::max(qv_arr(i,j,k), Real(0.0)) / dtcld);
3196  if (qc_arr(i,j,k) > Real(0.0) && work1_arr(i,j,k,0) < Real(0.0)) {
3197  pcond = amrex::max(work1_arr(i,j,k,0), -qc_arr(i,j,k)) / dtcld;
3198  }
3199  pcond_arr(i,j,k) = pcond;
3200 
3201  if (pcond == -qc_arr(i,j,k) / dtcld) {
3202  nn_arr(i,j,k) += nc_arr(i,j,k);
3203  nc_arr(i,j,k) = Real(0.0);
3204  }
3205 
3206  qv_arr(i,j,k) = amrex::max(qv_arr(i,j,k) - pcond * dtcld, Real(0.0));
3207  qc_arr(i,j,k) = amrex::max(qc_arr(i,j,k) + pcond * dtcld, Real(0.0));
3208  t_arr(i,j,k) += pcond * xl_arr(i,j,k) / cpm_arr(i,j,k) * dtcld;
3209  });
3210 
3211 
3212  const Real g17_pidnc = pi_wdm6_loc * Real(denr) / Real(6.0);
3213  const Real g17_pidnr = Real(4.0) * pi_wdm6_loc * Real(denr);
3214 
3215 
3216  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3217  if (qc_arr(i,j,k) <= Real(qmin)) qc_arr(i,j,k) = Real(0.0);
3218  if (qi_arr(i,j,k) <= Real(qmin)) qi_arr(i,j,k) = Real(0.0);
3219 
3220  if (qr_arr(i,j,k) >= Real(qcrmin) && nr_arr(i,j,k) >= Real(nrmin)) {
3221  Real lamdr = std::exp(std::log(
3222  (g17_pidnr * nr_arr(i,j,k)) / (den_arr(i,j,k) * qr_arr(i,j,k))
3223  ) * wdm6_literal(0.33333333));
3224  if (lamdr <= Real(lamdarmin)) {
3225  lamdr = Real(lamdarmin);
3226  nr_arr(i,j,k) = den_arr(i,j,k) * qr_arr(i,j,k)
3227  * std::pow(lamdr, Real(3.0)) / g17_pidnr;
3228  } else if (lamdr >= Real(lamdarmax)) {
3229  lamdr = Real(lamdarmax);
3230  nr_arr(i,j,k) = den_arr(i,j,k) * qr_arr(i,j,k)
3231  * std::pow(lamdr, Real(3.0)) / g17_pidnr;
3232  }
3233  }
3234 
3235  if (qc_arr(i,j,k) >= Real(qmin) && nc_arr(i,j,k) >= Real(ncmin)) {
3236  Real lamdc = std::exp(std::log(
3237  (g17_pidnc * nc_arr(i,j,k)) / (den_arr(i,j,k) * qc_arr(i,j,k))
3238  ) * wdm6_literal(0.33333333));
3239  if (lamdc <= Real(lamdacmin)) {
3240  lamdc = Real(lamdacmin);
3241  nc_arr(i,j,k) = den_arr(i,j,k) * qc_arr(i,j,k)
3242  * std::pow(lamdc, Real(3.0)) / g17_pidnc;
3243  } else if (lamdc >= Real(lamdacmax)) {
3244  lamdc = Real(lamdacmax);
3245  nc_arr(i,j,k) = den_arr(i,j,k) * qc_arr(i,j,k)
3246  * std::pow(lamdc, Real(3.0)) / g17_pidnc;
3247  }
3248  }
3249  });
3250 
3251 
3252 
3253 
3254  } // End minor timestep loop
3255 
3256 
3257  // Convert updated temperature back to potential temperature, mirroring
3258  // the bridge leg above. Without this the native path left
3259  // mic_fab_vars[theta] at its pre-microphysics value while the bridge
3260  // leg advanced it, and Copy_Micro_to_State then formed
3261  // RhoTheta = rho * theta from the stale field, dropping microphysical
3262  // latent heating from the thermodynamic state. Confirmed by comparing
3263  // theta before and after the bridge call: native PRE theta equalled
3264  // POST theta at 100/100 levels, and the resulting
3265  // divergence 6.625510053 at k=90 matched the step-1 theta error exactly.
3266  {
3267  const Real configured_rdOcp = m_rdOcp;
3268  const bool use_anelastic_reference_pressure =
3269  m_use_anelastic_reference_pressure;
3270  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
3271  w1_theta(i,j,k) = wdm6_theta_from_temperature_and_pressure(
3272  t_arr(i,j,k), p_arr(i,j,k), configured_rdOcp,
3273  use_anelastic_reference_pressure);
3274  });
3275  }
3276 
3277 #ifdef ERF_USE_WDM6_FORT
3278  }
3279 #endif
3280 
3281  } // MFIter loop
3282 }
wdm6_diffac
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_diffac(Real a, Real b, Real c, Real d, Real e, Real rv_arg)
Definition: ERF_AdvanceWDM6.cpp:63
wdm6_xka
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_xka(Real x, Real y)
Definition: ERF_AdvanceWDM6.cpp:58
wdm6_venfac
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_venfac(Real a, Real b, Real c, Real den0_arg)
Definition: ERF_AdvanceWDM6.cpp:70
wdm6_nislfv_rain_plm6_column
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_nislfv_rain_plm6_column(int km, Real *precip1, Real *precip2, Real dt, int iter, Real pidn0s, Real pidn0g, Real qcrmin, Real alpha, Real n0smax, Real n0s, Real t0c, Real rslopesmax, Real rslopesbmax, Real rslopes2max, Real rslopes3max, Real bvts, Real pvts, Real rslopegmax, Real rslopegbmax, Real rslopeg2max, Real rslopeg3max, Real bvtg, Real pvtg, Array4< Real > const &sed_cell, Array4< Real > const &sed_node, int i_s, int j_s, int klo_s)
Definition: ERF_AdvanceWDM6.cpp:280
wdm6_xlcal
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_xlcal(Real x, Real xlv0_arg, Real xlv1_arg, Real t0c_arg)
Definition: ERF_AdvanceWDM6.cpp:34
wdm6_slope_t0c
constexpr amrex::Real wdm6_slope_t0c
Definition: ERF_AdvanceWDM6.cpp:198
wdm6_slope_rain_cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_slope_rain_cell(Real qr, Real nr, Real den, Real denfac, Real qcrmin_arg, Real nrmin_arg, Real rslopermax_arg, Real rsloperbmax_arg, Real rsloper2max_arg, Real rsloper3max_arg, Real bvtr_arg, Real pvtr_arg, Real pvtrn_arg, Real pidnr_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt, Real &vtn)
Definition: ERF_AdvanceWDM6.cpp:140
wdm6_slope_snow_cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_slope_snow_cell(Real qs, Real den, Real denfac, Real t, Real pidn0s_arg, Real alpha_arg, Real n0smax_arg, Real n0s_arg, Real t0c_arg, Real qcrmin_arg, Real rslopesmax_arg, Real rslopesbmax_arg, Real rslopes2max_arg, Real rslopes3max_arg, Real bvts_arg, Real pvts_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt, Real &n0sfac)
Definition: ERF_AdvanceWDM6.cpp:217
wdm6_cpmcal
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_cpmcal(Real x, Real qmin_arg, Real cpd_arg, Real cpv_arg)
Definition: ERF_AdvanceWDM6.cpp:28
wdm6_xni_exact
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_xni_exact(Real qi, Real den, Real qmin_arg)
Definition: ERF_AdvanceWDM6.cpp:133
wdm6_default_real_pow
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_default_real_pow(double base, double exponent)
Definition: ERF_AdvanceWDM6.cpp:39
wdm6_slope_graup_cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_slope_graup_cell(Real qg, Real den, Real denfac, Real pidn0g_arg, Real qcrmin_arg, Real rslopegmax_arg, Real rslopegbmax_arg, Real rslopeg2max_arg, Real rslopeg3max_arg, Real bvtg_arg, Real pvtg_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt)
Definition: ERF_AdvanceWDM6.cpp:247
wdm6_conden
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_conden(Real a, Real b, Real c, Real d, Real e, Real qmin_arg, Real rv_arg)
Definition: ERF_AdvanceWDM6.cpp:79
wdm6_rslopec_exact
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_rslopec_exact(Real qc, Real den, Real nc, Real pidnc_arg)
Definition: ERF_AdvanceWDM6.cpp:122
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:35
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:36
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:56
Cp_l
constexpr amrex::Real Cp_l
Definition: ERF_Constants.H:38
Cp_v
constexpr amrex::Real Cp_v
Definition: ERF_Constants.H:37
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:34
klo
const int klo
Definition: ERF_InitCustomPert_ABL.H:75
pp
ParmParse pp("prob")
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
qv_arr
auto qv_arr
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:210
lat_vap
constexpr amrex::Real lat_vap
Definition: ERF_MicrophysicsConstants.H:113
lsub
constexpr amrex::Real lsub
Definition: ERF_MicrophysicsConstants.H:111
rhoh2o
constexpr amrex::Real rhoh2o
Definition: ERF_MicrophysicsConstants.H:42
lat_ice
constexpr amrex::Real lat_ice
Definition: ERF_MicrophysicsConstants.H:114
rhos
constexpr amrex::Real rhos
Definition: ERF_MicrophysicsConstants.H:40
Arena_Used
Arena * Arena_Used
Definition: ERF_Morrison_Advance_F.H:23
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
wdm6_literal
constexpr amrex::Real wdm6_literal(double d)
Definition: ERF_WDM6.H:96
wdm6_theta_from_temperature_and_pressure
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real wdm6_theta_from_temperature_and_pressure(const amrex::Real temperature, const amrex::Real pressure, const amrex::Real configured_rdOcp, const bool use_anelastic_reference_pressure) noexcept
Definition: ERF_WDM6.H:81
mp_wdm6_run_c
void mp_wdm6_run_c(double *t, double *qv, double *qc, double *qi, double *qr, double *qs, double *qg, double *nn, double *nc, double *nr, double *den, double *p, double *delz, double delt, double g, double cpd, double cpv, double rd, double rv, double t0c, double ep1, double ep2, double qmin, double xls, double xlv0, double xlf0, double den0, double denr, double cliq, double cice, double psat, double ccn0, double *xland, double *rain, double *rainncv, double *sr, double *snow, double *snowncv, double *graupel, double *graupelncv, int ims, int ime, int jms, int jme, int kms, int kme, int its, int ite, int jts, int jte, int kts, int kte, int microphysics_debug, int diag_i_dbg, int diag_j_dbg)
mp_wdm6_init_c
void mp_wdm6_init_c(double den0, double denr, double dens, double cl, double cpv, double ccn0, int hail_opt)
NullMoist::m_use_anelastic_reference_pressure
bool m_use_anelastic_reference_pressure
Definition: ERF_NullMoist.H:175
WDM6::m_precs2
amrex::Real m_precs2
Definition: ERF_WDM6.H:315
WDM6::ncmin
static constexpr amrex::Real ncmin
Definition: ERF_WDM6.H:188
WDM6::qcrmin
static constexpr amrex::Real qcrmin
Definition: ERF_WDM6.H:187
WDM6::m_pidn0g
amrex::Real m_pidn0g
Definition: ERF_WDM6.H:319
WDM6::m_precr2
amrex::Real m_precr2
Definition: ERF_WDM6.H:311
WDM6::m_hail_opt
bool m_hail_opt
Definition: ERF_WDM6.H:298
WDM6::m_lmask
amrex::iMultiFab * m_lmask
Definition: ERF_WDM6.H:290
WDM6::nrmin
static constexpr amrex::Real nrmin
Definition: ERF_WDM6.H:189
WDM6::m_pvtrn
amrex::Real m_pvtrn
Definition: ERF_WDM6.H:310
WDM6::m_rslopecmax
amrex::Real m_rslopecmax
Definition: ERF_WDM6.H:320
WDM6::lamdarmax
static constexpr amrex::Real lamdarmax
Definition: ERF_WDM6.H:178
WDM6::m_pidn0s
amrex::Real m_pidn0s
Definition: ERF_WDM6.H:316
WDM6::m_qc1
amrex::Real m_qc1
Definition: ERF_WDM6.H:305
WDM6::di100
static constexpr amrex::Real di100
Definition: ERF_WDM6.H:201
WDM6::m_pvts
amrex::Real m_pvts
Definition: ERF_WDM6.H:315
WDM6::pfrz1
static constexpr amrex::Real pfrz1
Definition: ERF_WDM6.H:185
WDM6::m_rslopec3max
amrex::Real m_rslopec3max
Definition: ERF_WDM6.H:320
WDM6::dicon
static constexpr amrex::Real dicon
Definition: ERF_WDM6.H:183
WDM6::m_pvtr
amrex::Real m_pvtr
Definition: ERF_WDM6.H:310
WDM6::m_rslopegmax
amrex::Real m_rslopegmax
Definition: ERF_WDM6.H:321
WDM6::bvtr
static constexpr amrex::Real bvtr
Definition: ERF_WDM6.H:169
WDM6::m_rslopesmax
amrex::Real m_rslopesmax
Definition: ERF_WDM6.H:321
WDM6::m_pvtg
amrex::Real m_pvtg
Definition: ERF_WDM6.H:319
WDM6::mic_fab_vars
amrex::Array< FabPtr, MicVar_WDM6::NumVars > mic_fab_vars
Definition: ERF_WDM6.H:292
WDM6::lamdacmax
static constexpr amrex::Real lamdacmax
Definition: ERF_WDM6.H:181
WDM6::alpha_wdm6
static constexpr amrex::Real alpha_wdm6
Definition: ERF_WDM6.H:195
WDM6::m_pi_wdm6
amrex::Real m_pi_wdm6
Definition: ERF_WDM6.H:304
WDM6::m_rslopes2max
amrex::Real m_rslopes2max
Definition: ERF_WDM6.H:323
WDM6::m_rslopegbmax
amrex::Real m_rslopegbmax
Definition: ERF_WDM6.H:322
WDM6::m_rslopes3max
amrex::Real m_rslopes3max
Definition: ERF_WDM6.H:324
WDM6::m_qck1
amrex::Real m_qck1
Definition: ERF_WDM6.H:305
WDM6::m_rdOcp
amrex::Real m_rdOcp
Definition: ERF_WDM6.H:285
WDM6::qs0
static constexpr amrex::Real qs0
Definition: ERF_WDM6.H:192
WDM6::m_rslopeg2max
amrex::Real m_rslopeg2max
Definition: ERF_WDM6.H:323
WDM6::n0s
static constexpr amrex::Real n0s
Definition: ERF_WDM6.H:194
WDM6::m_rslopesbmax
amrex::Real m_rslopesbmax
Definition: ERF_WDM6.H:322
WDM6::m_geom
amrex::Geometry m_geom
Definition: ERF_WDM6.H:278
WDM6::di2000
static constexpr amrex::Real di2000
Definition: ERF_WDM6.H:203
WDM6::m_precg1
amrex::Real m_precg1
Definition: ERF_WDM6.H:319
WDM6::m_rsloperbmax
amrex::Real m_rsloperbmax
Definition: ERF_WDM6.H:322
WDM6::m_precr1
amrex::Real m_precr1
Definition: ERF_WDM6.H:311
WDM6::di82
static constexpr amrex::Real di82
Definition: ERF_WDM6.H:204
WDM6::m_g7pbr
amrex::Real m_g7pbr
Definition: ERF_WDM6.H:308
WDM6::actr
static constexpr amrex::Real actr
Definition: ERF_WDM6.H:198
WDM6::lamdarmin
static constexpr amrex::Real lamdarmin
Definition: ERF_WDM6.H:179
WDM6::lamdacmin
static constexpr amrex::Real lamdacmin
Definition: ERF_WDM6.H:182
WDM6::m_ccn0
amrex::Real m_ccn0
Definition: ERF_WDM6.H:281
WDM6::m_pidnc
amrex::Real m_pidnc
Definition: ERF_WDM6.H:305
WDM6::m_qc0
amrex::Real m_qc0
Definition: ERF_WDM6.H:305
WDM6::m_bvtg
amrex::Real m_bvtg
Definition: ERF_WDM6.H:299
WDM6::m_rsloper3max
amrex::Real m_rsloper3max
Definition: ERF_WDM6.H:324
WDM6::m_rslopeg3max
amrex::Real m_rslopeg3max
Definition: ERF_WDM6.H:324
WDM6::m_pidnr
amrex::Real m_pidnr
Definition: ERF_WDM6.H:312
WDM6::m_g4pbr
amrex::Real m_g4pbr
Definition: ERF_WDM6.H:308
WDM6::bvts
static constexpr amrex::Real bvts
Definition: ERF_WDM6.H:177
WDM6::m_n0g
amrex::Real m_n0g
Definition: ERF_WDM6.H:299
WDM6::di600
static constexpr amrex::Real di600
Definition: ERF_WDM6.H:202
WDM6::m_precg2
amrex::Real m_precg2
Definition: ERF_WDM6.H:319
WDM6::m_rslopermax
amrex::Real m_rslopermax
Definition: ERF_WDM6.H:321
WDM6::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_WDM6.H:288
WDM6::ncrk2
static constexpr amrex::Real ncrk2
Definition: ERF_WDM6.H:200
WDM6::m_xlv1
amrex::Real m_xlv1
Definition: ERF_WDM6.H:304
WDM6::m_pacrg
amrex::Real m_pacrg
Definition: ERF_WDM6.H:319
WDM6::ncrk1
static constexpr amrex::Real ncrk1
Definition: ERF_WDM6.H:199
WDM6::avtr
static constexpr amrex::Real avtr
Definition: ERF_WDM6.H:168
WDM6::m_pacrc
amrex::Real m_pacrc
Definition: ERF_WDM6.H:316
WDM6::m_roqimax
amrex::Real m_roqimax
Definition: ERF_WDM6.H:311
WDM6::dimax
static constexpr amrex::Real dimax
Definition: ERF_WDM6.H:184
WDM6::actk
static constexpr amrex::Real actk
Definition: ERF_WDM6.H:197
WDM6::pfrz2
static constexpr amrex::Real pfrz2
Definition: ERF_WDM6.H:186
WDM6::m_precs1
amrex::Real m_precs1
Definition: ERF_WDM6.H:315
WDM6::satmax
static constexpr amrex::Real satmax
Definition: ERF_WDM6.H:196
WDM6::n0smax
static constexpr amrex::Real n0smax
Definition: ERF_WDM6.H:193
WDM6::m_rsloper2max
amrex::Real m_rsloper2max
Definition: ERF_WDM6.H:323
WDM6::m_rslopec2max
amrex::Real m_rslopec2max
Definition: ERF_WDM6.H:320
MORRInd::xlf
@ xlf
Definition: ERF_AdvanceMorrison.cpp:157
MicVar_WDM6::qr
@ qr
Definition: ERF_WDM6.H:29
MicVar_WDM6::qv
@ qv
Definition: ERF_WDM6.H:26
MicVar_WDM6::qc
@ qc
Definition: ERF_WDM6.H:27
MicVar_WDM6::qi
@ qi
Definition: ERF_WDM6.H:28
MicVar_WDM6::graup_accum
@ graup_accum
Definition: ERF_WDM6.H:37
MicVar_WDM6::rain_accum
@ rain_accum
Definition: ERF_WDM6.H:35
MicVar_WDM6::pres
@ pres
Definition: ERF_WDM6.H:25
MicVar_WDM6::nr
@ nr
Definition: ERF_WDM6.H:34
MicVar_WDM6::qg
@ qg
Definition: ERF_WDM6.H:31
MicVar_WDM6::theta
@ theta
Definition: ERF_WDM6.H:23
MicVar_WDM6::qs
@ qs
Definition: ERF_WDM6.H:30
MicVar_WDM6::nc
@ nc
Definition: ERF_WDM6.H:33
MicVar_WDM6::nn
@ nn
Definition: ERF_WDM6.H:32
MicVar_WDM6::rho
@ rho
Definition: ERF_WDM6.H:22
MicVar_WDM6::tabs
@ tabs
Definition: ERF_WDM6.H:24
MicVar_WDM6::snow_accum
@ snow_accum
Definition: ERF_WDM6.H:36
WDM6SedCellScratch::tk
@ tk
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::work_col
@ work_col
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::den
@ den
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::denfac
@ denfac
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::rq2_col
@ rq2_col
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::rq_col
@ rq_col
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::NumComps
@ NumComps
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
WDM6SedNodeScratch::NumComps
@ NumComps
Definition: ERF_AdvanceWDM6.cpp:276
WSM6Ind::fall_s
@ fall_s
Definition: ERF_WSM6.H:347
WSM6Ind::n0sfac
@ n0sfac
Definition: ERF_WSM6.H:333
WSM6Ind::work2
@ work2
Definition: ERF_WSM6.H:326
WSM6Ind::pcond
@ pcond
Definition: ERF_WSM6.H:296
WSM6Ind::qsum
@ qsum
Definition: ERF_WSM6.H:323
WSM6Ind::psmlt
@ psmlt
Definition: ERF_WSM6.H:317
WSM6Ind::fall_g
@ fall_g
Definition: ERF_WSM6.H:347
WSM6Ind::work1c
@ work1c
Definition: ERF_WSM6.H:286
WSM6Ind::pgmlt
@ pgmlt
Definition: ERF_WSM6.H:318
WSM6Ind::rhi
@ rhi
Definition: ERF_WSM6.H:339
WSM6Ind::fall_r
@ fall_r
Definition: ERF_WSM6.H:347
WSM6Ind::xni
@ xni
Definition: ERF_WSM6.H:328
module_model_constants::cice
real(c_double), parameter cice
Definition: ERF_module_model_constants.F90:30
module_model_constants::g
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
module_model_constants::cliq
real(c_double), parameter cliq
Definition: ERF_module_model_constants.F90:29
module_model_constants::xlv0
real(c_double), parameter xlv0
Definition: ERF_module_model_constants.F90:51
module_model_constants::cpv
real(c_double), parameter cpv
Definition: ERF_module_model_constants.F90:26
module_model_constants::xls
real(c_double), parameter xls
Definition: ERF_module_model_constants.F90:56
module_model_constants::psat
real(c_double), parameter psat
Definition: ERF_module_model_constants.F90:31
mp_wdm6::dens
real(kind=kind_phys), parameter, private dens
Definition: ERF_module_mp_wdm6.F90:61
Here is the call graph for this function:

◆ Copy_Micro_to_State()

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

Reimplemented from NullMoist.

10 {
11  // Conservative update of all fields
12  for (MFIter mfi(cons, TilingIfNotGPU()); mfi.isValid(); ++mfi)
13  {
14  const auto& box3d = mfi.tilebox();
15  auto states = cons.array(mfi);
16 
17  const auto& rho = mic_fab_vars[MicVar_WDM6::rho]->array(mfi);
18  const auto& theta = mic_fab_vars[MicVar_WDM6::theta]->array(mfi);
19  const auto& qv = mic_fab_vars[MicVar_WDM6::qv]->array(mfi);
20  const auto& qc = mic_fab_vars[MicVar_WDM6::qc]->array(mfi);
21  const auto& qi = mic_fab_vars[MicVar_WDM6::qi]->array(mfi);
22  const auto& qr = mic_fab_vars[MicVar_WDM6::qr]->array(mfi);
23  const auto& qs = mic_fab_vars[MicVar_WDM6::qs]->array(mfi);
24  const auto& qg = mic_fab_vars[MicVar_WDM6::qg]->array(mfi);
25  const auto& nn = mic_fab_vars[MicVar_WDM6::nn]->array(mfi);
26  const auto& nc = mic_fab_vars[MicVar_WDM6::nc]->array(mfi);
27  const auto& nr = mic_fab_vars[MicVar_WDM6::nr]->array(mfi);
28 
29  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
30  {
31  states(i,j,k,RhoTheta_comp) = rho(i,j,k) * theta(i,j,k);
32 
33  states(i,j,k,RhoQ1_comp) = rho(i,j,k) * amrex::max(Real(0), qv(i,j,k));
34  states(i,j,k,RhoQ2_comp) = rho(i,j,k) * amrex::max(Real(0), qc(i,j,k));
35  states(i,j,k,RhoQ3_comp) = rho(i,j,k) * amrex::max(Real(0), qi(i,j,k));
36  states(i,j,k,RhoQ4_comp) = rho(i,j,k) * amrex::max(Real(0), qr(i,j,k));
37  states(i,j,k,RhoQ5_comp) = rho(i,j,k) * amrex::max(Real(0), qs(i,j,k));
38  states(i,j,k,RhoQ6_comp) = rho(i,j,k) * amrex::max(Real(0), qg(i,j,k));
39 
40  // Number concentrations
41  states(i,j,k,RhoQ7_comp) = rho(i,j,k) * amrex::max(Real(0), nc(i,j,k));
42  states(i,j,k,RhoQ8_comp) = rho(i,j,k) * amrex::max(Real(0), nn(i,j,k));
43  states(i,j,k,RhoQ9_comp) = rho(i,j,k) * amrex::max(Real(0), nr(i,j,k));
44  });
45  }
46 
47  // No nn bookkeeping needed here. Copy_State_to_Micro decides per cell from
48  // the state itself, so once RhoQ8 has been written by the copyback above it
49  // is read back like any other variable.
50 
51  cons.FillBoundary(m_geom.periodicity());
52 }
RhoQ4_comp
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:48
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:46
RhoQ3_comp
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:47
RhoQ9_comp
#define RhoQ9_comp
Definition: ERF_IndexDefines.H:53
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
RhoQ6_comp
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:50
RhoQ8_comp
#define RhoQ8_comp
Definition: ERF_IndexDefines.H:52
RhoQ5_comp
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:49
RhoQ7_comp
#define RhoQ7_comp
Definition: ERF_IndexDefines.H:51
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:19
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:31
MicVar_Morr::nc
@ nc
Definition: ERF_Morrison.H:46
MicVar_Morr::nr
@ nr
Definition: ERF_Morrison.H:47
MicVar_SatAdj::qc
@ qc
Definition: ERF_SatAdj.H:42
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
WDM6SedCellScratch::qr
@ qr
Definition: ERF_AdvanceWDM6.cpp:271

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro() [1/2]

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

Reimplemented from NullMoist.

Referenced by Update_Micro_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro() [2/2]

void WDM6::Copy_State_to_Micro ( const amrex::MultiFab &  cons_in,
const amrex::MultiFab *  base_state 
)

◆ Define()

void WDM6::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

116  {
117  m_moisture_type = sc.moisture_type;
118  m_rdOcp = sc.rdOcp;
119  set_anelastic_reference_pressure_mode(sc);
120  m_axis = sc.ave_plane;
121  m_do_cond = (!sc.uses_shoc_family());
122  }
NullMoist::set_anelastic_reference_pressure_mode
void set_anelastic_reference_pressure_mode(const SolverChoice &sc)
Definition: ERF_NullMoist.H:155
WDM6::m_axis
int m_axis
Definition: ERF_WDM6.H:283
WDM6::m_moisture_type
MoistureType m_moisture_type
Definition: ERF_WDM6.H:286
WDM6::m_do_cond
bool m_do_cond
Definition: ERF_WDM6.H:284
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2237
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2260
SolverChoice::rdOcp
amrex::Real rdOcp
Ratio of dry-air gas constant to c_p.
Definition: ERF_DataStruct.H:2062
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:2177
Here is the call graph for this function:

◆ Get_Surface_Precip_Accumulation_Ptrs()

SurfacePrecipAccumulationSources WDM6::Get_Surface_Precip_Accumulation_Ptrs ( const int &  ) const
inlineoverridevirtual

Reimplemented from NullMoist.

263  {
264  SurfacePrecipAccumulationSources sources;
265  sources.total = {mic_fab_vars[MicVar_WDM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
266  sources.snow = {mic_fab_vars[MicVar_WDM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
267  sources.graupel = {mic_fab_vars[MicVar_WDM6::graup_accum].get(), rhoh2o / amrex::Real(1000.0)};
268  return sources;
269  }
SurfacePrecipAccumulationSources
Definition: ERF_SurfacePrecipitation.H:34
SurfacePrecipAccumulationSources::snow
SurfacePrecipAccumulationSource snow
Definition: ERF_SurfacePrecipitation.H:37
SurfacePrecipAccumulationSources::total
SurfacePrecipAccumulationSource total
Definition: ERF_SurfacePrecipitation.H:35
SurfacePrecipAccumulationSources::graupel
SurfacePrecipAccumulationSource graupel
Definition: ERF_SurfacePrecipitation.H:38

◆ Init()

void WDM6::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

Reimplemented from NullMoist.

17 {
18  dt = dt_advance;
19  m_geom = geom;
20 
21  m_z_phys_nd = z_phys_nd.get();
22  m_detJ_cc = detJ_cc.get();
23 
24  // Read CCN concentration from input
25  amrex::ParmParse pp("wdm6");
26  pp.queryAdd("ccn0", m_ccn0); // default 100.0e6 m^-3
27 
28  // Graupel/hail regime selector. The Fortran gates on `hail_opt .eq. 1`
29  // (ERF_module_mp_wdm6.F90:3259), so every value other than 1 selects the
30  // graupel branch; the bool below reproduces that test exactly rather than
31  // inventing a stricter 0/1 validation the Fortran does not perform.
32  // Must be resolved before initialize_coeffs(), which branches on it.
33  int hail_opt_in = 0;
34  pp.queryAdd("hail_opt", hail_opt_in); // default 0 = graupel
35  m_hail_opt = (hail_opt_in == 1);
36 
37  amrex::Print() << "WDM6 Initialization: CCN0 = " << m_ccn0 << " #/m^3"
38  << ", hail_opt = " << (m_hail_opt ? 1 : 0)
39  << (m_hail_opt ? " (hail)" : " (graupel)") << "\n";
40 
41  MicVarMap.resize(m_qmoist_size);
42  MicVarMap = {MicVar_WDM6::rain_accum, MicVar_WDM6::snow_accum, MicVar_WDM6::graup_accum};
43 
44  // Select appropriate Arena based on execution mode
45  // - Fortran bridge + GPU: Use managed memory for CPU (Fortran) ↔ GPU data transfer
46  // - Fortran bridge CPU-only: Use pinned memory for CPU efficiency
47  // - C++ GPU kernels: Use standard device memory
48 #if defined(ERF_USE_WDM6_FORT) && defined(AMREX_USE_GPU)
49  Arena* Arena_Used = The_Managed_Arena();
50 #elif defined(ERF_USE_WDM6_FORT)
51  Arena* Arena_Used = The_Pinned_Arena();
52 #else
53  Arena* Arena_Used = The_Arena();
54 #endif
55 
56  for (int ivar = 0; ivar < MicVar_WDM6::NumVars; ++ivar) {
57  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
58  1, cons_in.nGrowVect(),
59  MFInfo().SetArena(Arena_Used));
60  mic_fab_vars[ivar]->setVal(0.0);
61  }
62 
63  // Initialize nn to ccn0 / rho (matches WRF itimestep==1 behavior)
64  // This must be done here so nn is available before the first microphysics call
65  // NOTE: Can't write to cons_in (it's const), so only initialize mic_fab_vars.
66  // Copy_Micro_to_State will write it back to state after the first microphysics call.
67  // IMPORTANT: Use growntilebox to include ghost zones, since Copy_State_to_Micro will
68  // skip reading nn and expect it to be initialized everywhere.
69  // nn is NOT seeded here. The Fortran seeds it inside the per-step driver
70  // wdm6 (:219), immediately before the kernel, not in wdm6init, so the
71  // faithful ERF home is Copy_State_to_Micro, which is the per-step pack.
72  // Seeding here instead was also wrong on restart: Init runs again on a
73  // restart, so it discarded the checkpoint's evolved RhoQ8 and the first
74  // Copy_State_to_Micro kept the uniform background. See the state-based
75  // seeding in Copy_State_to_Micro below for the condition and its rationale.
76 
77  nlev = m_geom.Domain().length(2);
78  zlo = m_geom.Domain().smallEnd(2);
79  zhi = m_geom.Domain().bigEnd(2);
80 
81  initialize_coeffs();
82 }
WDM6::nlev
int nlev
Definition: ERF_WDM6.H:282
WDM6::dt
amrex::Real dt
Definition: ERF_WDM6.H:279
WDM6::m_qmoist_size
int m_qmoist_size
Definition: ERF_WDM6.H:272
WDM6::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_WDM6.H:276
WDM6::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_WDM6.H:289
WDM6::zlo
int zlo
Definition: ERF_WDM6.H:282
WDM6::zhi
int zhi
Definition: ERF_WDM6.H:282
WDM6::initialize_coeffs
void initialize_coeffs()
Definition: ERF_InitWDM6.cpp:170
MicVar_WDM6::NumVars
@ NumVars
Definition: ERF_WDM6.H:38
Here is the call graph for this function:

◆ initialize_coeffs()

void WDM6::initialize_coeffs ( )
private
171 {
172  using amrex::Real;
173 
174  // Exact port of Fortran rgmma() function from WSM6
175  auto rgmma = [](Real x) -> Real {
176  if (x == Real(1.0)) return Real(0.0);
177  constexpr Real euler = Real(0.577215664901532);
178  Real result = x * std::exp(euler * x);
179  for (int i = 1; i <= 10000; ++i) {
180  Real y = Real(i);
181  result = result * (Real(1.0) + x/y) * std::exp(-x/y);
182  }
183  return Real(1.0) / result;
184  };
185 
186  // Physical constants
187  const Real den0 = Real(1.28);
188  const Real denr = Real(rhoh2o);
189  const Real dens_arg = dens_snow;
190  const Real cl = Real(Cp_l);
191  const Real cpv_loc = Real(Cp_v);
192 
193  // hail_opt branch
194  if (m_hail_opt) {
195  m_n0g = Real(4.0e4);
196  m_deng = Real(700.0);
197  m_avtg = Real(285.0);
198  m_bvtg = wdm6_literal(0.8);
199  m_lamdagmax = Real(2.0e4);
200  } else {
201  m_n0g = Real(4.0e6);
202  m_deng = Real(500.0);
203  m_avtg = Real(330.0);
204  m_bvtg = wdm6_literal(0.8);
205  m_lamdagmax = Real(6.0e4);
206  }
207 
208  // LITERAL PRECISION CONTRACT, second hazard class. The Fortran writes
209  // pi = 4.*atan(1.)
210  // at ERF_module_mp_wdm6.F90:3111. Neither literal carries a kind suffix, so
211  // `1.` is a default REAL and the generic ATAN resolves to its SINGLE
212  // precision specific; the whole expression is evaluated in float32 and only
213  // then widened to kind_phys. Measured with gfortran:
214  // 4.*atan(1.) = 3.14159274101257324219
215  // 4._kp*atan(1._kp) = 3.14159265358979311600
216  // relative difference = 2.7828e-08
217  // This is distinct from the unsuffixed-parameter hazard the contract already
218  // covers: nothing here is a rounded constant, it is an intrinsic whose result
219  // kind is chosen by its argument kind. pi is load-bearing -- it feeds pidn0s,
220  // pidnc, pidnr, pidn0r, pidn0s, pidn0g, pacrs, pacrc and the pi*pi in G10c's
221  // pfrzdtc -- so the error propagates broadly.
222  // Routed through wdm6_literal like every other constant: it is constexpr but
223  // not restricted to constant expressions, so it applies the same rounding to
224  // a runtime value. Do not re-expand the toggle inline here; the contract is
225  // that ERF_WDM6_F32_LITERALS is honoured in exactly one place.
226  m_pi_wdm6 = wdm6_literal(Real(4.0) * std::atan(Real(1.0)));
227  m_xlv1 = cl - cpv_loc;
228 
229  // Cloud droplet parameters. The Fortran, ERF_module_mp_wdm6.F90:3114-3116:
230  // qc0 = 4./3.*pi*denr*r0**3.*xncr0/den0
231  // qc1 = 4./3.*pi*denr*r0**3.*xncr1/den0
232  // qck1 = .104*9.8*peaut/(denr)**(1./3.)/xmyu*den0**(4./3.)
233  // qc0 takes xncr0 and qc1 takes xncr1, unconditionally: the maritime versus
234  // continental choice is made downstream by the slmsk branch that sets
235  // qcr = qc0 or qc1 (:628-632), not here. The previous code used xncr for qc0
236  // and a ccn0-dependent xncr_use for qc1, so with the default ccn0 = 1.0e8 it
237  // selected xncr0 for BOTH.
238  //
239  // qck1 is the consequential one. WDM6 divides by denr**(1./3.) alone; the
240  // previous code divided by (xncr*denr)**(1./3.), which is the WSM6 form at
241  // ERF_module_mp_wsm6.F90:105. WSM6 is single-moment and folds its fixed
242  // droplet number xncr into qck1; WDM6 is double-moment and already carries
243  // number explicitly as ncr(:,:,2)**(-1./3.) in praut, so the WSM6 form
244  // double-counts it. The error is exactly xncr**(1./3.) = (3.0e8)**(1/3) =
245  // 669.4, and praut and nraut were both measured low by that factor at
246  // (199,3,43) from bitwise-equal inputs.
247  m_qc0 = wdm6_literal(4.0/3.0) * m_pi_wdm6 * denr
248  * std::pow(r0, Real(3.0)) * xncr0 / den0;
249  m_qc1 = wdm6_literal(4.0/3.0) * m_pi_wdm6 * denr
250  * std::pow(r0, Real(3.0)) * xncr1 / den0;
251  // A FOURTH HAZARD CLASS under the literal precision contract: the float
252  // assumptions apply to the whole literal PRODUCT, not merely to each
253  // literal. In the Fortran, .104 and 9.8 are both default REAL, so
254  // `.104*9.8` is itself a default-REAL expression and the PRODUCT is
255  // evaluated in single precision before anything widens it:
256  // float32(float32(.104)*float32(9.8)) = 1.01920008659362792969
257  // Applying the float assumption to each operand separately and then
258  // multiplying performs the product in double:
259  // float32(.104)*float32(9.8) = 1.01920004086494486728
260  // a 4.486723e-08 relative difference. It lands on qck1 as a constant
261  // multiplicative factor, so praut carries exactly that relative error at
262  // every k, and nraut = 3.5e9*den*praut inherits it unchanged.
263  //
264  // Proven by reconstruction at clean SHA ac1e445fa from bitwise-equal
265  // PRE_G13A inputs, no Tier 2 instrumentation: of 61 live G13a cells across
266  // columns (108,3) and (199,3) over two steps, 59 reproduce BOTH legs'
267  // printed praut bitwise to 20 digits from the two qck1 variants; the other
268  // 2 are cells where the min(praut, qc/dtcld) clamp binds and the legs
269  // already agree bitwise.
270  //
271  // The nesting is load-bearing. wdm6_literal(0.104 * 9.8) is NOT equivalent
272  // and is worse at 1.169636e-07: it rounds the double product, whereas the
273  // Fortran rounds the product of already-rounded operands. Round the
274  // operands first, then round their product. Under Option B the helper is
275  // the identity, so this collapses to 0.104*9.8 in double, which is what
276  // -fdefault-real-8 gives on the Fortran side.
277  m_qck1 = wdm6_literal(wdm6_literal(0.104) * wdm6_literal(9.8)) * peaut
278  / std::pow(denr, wdm6_literal(1.0/3.0))
279  / xmyu * std::pow(den0, wdm6_literal(4.0/3.0));
280  m_pidnc = m_pi_wdm6 * denr / Real(6.0);
281 
282  // Rain coefficients
283  m_bvtr1 = Real(1.0) + bvtr;
284  m_bvtr2 = Real(2.0) + bvtr;
285  m_bvtr3 = Real(3.0) + bvtr;
286  m_bvtr4 = Real(4.0) + bvtr;
287  m_bvtr5 = Real(5.0) + bvtr;
288  m_bvtr6 = Real(6.0) + bvtr;
289  m_bvtr7 = Real(7.0) + bvtr;
290  m_bvtr2o5 = Real(2.5) + Real(0.5) * bvtr;
291  m_bvtr3o5 = Real(3.5) + Real(0.5) * bvtr;
292 
293  m_g1pbr = rgmma(m_bvtr1);
294  m_g2pbr = rgmma(m_bvtr2);
295  m_g3pbr = rgmma(m_bvtr3);
296  m_g4pbr = rgmma(m_bvtr4);
297  m_g5pbr = rgmma(m_bvtr5);
298  m_g6pbr = rgmma(m_bvtr6);
299  m_g7pbr = rgmma(m_bvtr7);
300  m_g5pbro2 = rgmma(m_bvtr2o5);
301  m_g7pbro2 = rgmma(m_bvtr3o5);
302 
303  m_pvtr = avtr * m_g5pbr / Real(24.0);
304  m_pvtrn = avtr * m_g2pbr;
305  m_eacrr = Real(1.0);
306  m_pacrr = m_pi_wdm6 * n0r * avtr * m_g3pbr * Real(0.25) * m_eacrr;
307  // Fortran WDM6 :3200-3201
308  // precr1 = 2.*pi*1.56
309  // precr2 = 2.*pi*.31*avtr**.5*g7pbro2
310  // NEITHER carries n0r. WSM6 :121-122 does -- precr1 = 2.*pi*n0r*.78 and
311  // precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2 -- for the same reason qck1
312  // differs between the schemes: WSM6 is single-moment and folds the fixed
313  // intercept n0r into its coefficients, while WDM6 is double-moment and
314  // multiplies by ncr(:,:,3) explicitly at the use site in prevp. precr1 had
315  // already been ported from the WDM6 form; precr2 kept WSM6's n0r, making it
316  // 8.0e6 times too large. That drove the raw prevp far enough negative to
317  // hit its satdt/2 clamp on the native leg only: at (103,3,34) on step 3 the
318  // bridge held prevp at -8.3286e-09 while the native leg sat exactly on the
319  // clamp at -2.5824185419052320E-06, 310x larger.
320  m_precr1 = wdm6_literal(2.0) * m_pi_wdm6 * wdm6_literal(1.56);
321  m_precr2 = wdm6_literal(2.0) * m_pi_wdm6 * wdm6_literal(0.31)
322  * std::pow(avtr, wdm6_literal(0.5)) * m_g7pbro2;
323  // Fortran: roqimax = 2.08e22*dimax**8
324  // Two mismatches, both confirmed against the Tier 2 psaut trace:
325  // - 2.08e22 is unsuffixed there, so it is float32-rounded before being
326  // promoted. Route it through wdm6_literal so it follows the same
327  // switch as the ERF_WDM6.H constants.
328  // - dimax**8 has an INTEGER exponent, which gfortran expands to repeated
329  // multiplication, not a libm pow call. Written as explicit squaring
330  // here so the two sides agree structurally rather than by luck. For
331  // this value the two routes happen to agree bitwise, but that is not
332  // guaranteed in general and should not be relied on.
333  {
334  const Real d2 = dimax * dimax;
335  const Real d4 = d2 * d2;
336  m_roqimax = wdm6_literal(2.08e22) * (d4 * d4);
337  }
338  // Fortran: xmmax = (dimax/dicon)**2, again an integer exponent.
339  {
340  const Real r = dimax / dicon;
341  m_xmmax = r * r;
342  }
343 
344  m_pidn0r = m_pi_wdm6 * denr * n0r;
345  m_pidnr = Real(4.0) * m_pi_wdm6 * denr;
346 
347  // Snow coefficients
348  m_bvts1 = Real(1.0) + bvts;
349  m_bvts2 = Real(2.5) + Real(0.5) * bvts;
350  m_bvts3 = Real(3.0) + bvts;
351  m_bvts4 = Real(4.0) + bvts;
352  m_g1pbs = rgmma(m_bvts1);
353  m_g3pbs = rgmma(m_bvts3);
354  m_g4pbs = rgmma(m_bvts4);
355  m_g5pbso2 = rgmma(m_bvts2);
356  m_pvts = avts * m_g4pbs / Real(6.0);
357  m_pacrs = m_pi_wdm6 * n0s * avts * m_g3pbs * Real(0.25);
358  // Fortran :3248-3249
359  // precs1 = 4.*n0s*.65
360  // precs2 = 4.*n0s*.44*avts**.5*g5pbso2
361  // .65 and .44 are unsuffixed so they carry the float assumption; 4. and
362  // avts**.5 are binary-exact and need no routing. Proven at Tier 2 on the
363  // G13E psdep decomposition at (199,3), k=85..98: coeres, rh-1 and work1
364  // are all bitwise, precs1 diverged exactly 3.667978e-08 and precs2 exactly
365  // 5.418604e-09, matching float32(.65) and float32(.44), and psdep's
366  // 2.523449e-08 is the term-weighted mix of the two.
367  m_precs1 = Real(4.0) * n0s * wdm6_literal(0.65);
368  m_precs2 = Real(4.0) * n0s * wdm6_literal(0.44)
369  * std::pow(avts, Real(0.5)) * m_g5pbso2;
370  m_pidn0s = m_pi_wdm6 * dens_arg * n0s;
371  m_pacrc = m_pi_wdm6 * n0s * avts * m_g3pbs * Real(0.25) * eacrc;
372 
373  // Graupel/hail coefficients
374  m_bvtg1 = Real(1.0) + m_bvtg;
375  m_bvtg2 = Real(2.5) + Real(0.5) * m_bvtg;
376  m_bvtg3 = Real(3.0) + m_bvtg;
377  m_bvtg4 = Real(4.0) + m_bvtg;
378  m_g1pbg = rgmma(m_bvtg1);
379  m_g3pbg = rgmma(m_bvtg3);
380  m_g4pbg = rgmma(m_bvtg4);
381  m_pacrg = m_pi_wdm6 * m_n0g * m_avtg * m_g3pbg * Real(0.25);
382  m_g5pbgo2 = rgmma(m_bvtg2);
383  m_pvtg = m_avtg * m_g4pbg / Real(6.0);
384  // Fortran :3264-3265, the graupel analogues of precs1/precs2 with the same
385  // unsuffixed .78 and .31. Routed on the same evidence class as precs1 and
386  // precs2 rather than on separate measurements: pgdep reads bitwise in the
387  // current column only because its branch is inactive there.
388  m_precg1 = Real(2.0) * m_pi_wdm6 * m_n0g * wdm6_literal(0.78);
389  m_precg2 = Real(2.0) * m_pi_wdm6 * m_n0g * wdm6_literal(0.31)
390  * std::pow(m_avtg, Real(0.5)) * m_g5pbgo2;
391  m_pidn0g = m_pi_wdm6 * m_deng * m_n0g;
392 
393  // Slope parameter limits
394  m_rslopecmax = Real(1.0) / lamdacmax;
395  m_rslopec2max = m_rslopecmax * m_rslopecmax;
396  m_rslopec3max = m_rslopec2max * m_rslopecmax;
397 
398  m_rslopermax = Real(1.0) / lamdarmax;
399  m_rslopesmax = Real(1.0) / lamdasmax;
400  m_rslopegmax = Real(1.0) / m_lamdagmax;
401  m_rsloperbmax = std::pow(m_rslopermax, bvtr);
402  m_rslopesbmax = std::pow(m_rslopesmax, bvts);
403  m_rslopegbmax = std::pow(m_rslopegmax, m_bvtg);
404  m_rsloper2max = m_rslopermax * m_rslopermax;
405  m_rslopes2max = m_rslopesmax * m_rslopesmax;
406  m_rslopeg2max = m_rslopegmax * m_rslopegmax;
407  m_rsloper3max = m_rsloper2max * m_rslopermax;
408  m_rslopes3max = m_rslopes2max * m_rslopesmax;
409  m_rslopeg3max = m_rslopeg2max * m_rslopegmax;
410 }
Coord::y
@ y
Coord::x
@ x
WDM6::m_g5pbro2
amrex::Real m_g5pbro2
Definition: ERF_WDM6.H:309
WDM6::m_bvtg1
amrex::Real m_bvtg1
Definition: ERF_WDM6.H:317
WDM6::n0r
static constexpr amrex::Real n0r
Definition: ERF_WDM6.H:167
WDM6::m_bvtr3o5
amrex::Real m_bvtr3o5
Definition: ERF_WDM6.H:307
WDM6::m_g5pbr
amrex::Real m_g5pbr
Definition: ERF_WDM6.H:308
WDM6::m_g3pbs
amrex::Real m_g3pbs
Definition: ERF_WDM6.H:314
WDM6::m_pacrs
amrex::Real m_pacrs
Definition: ERF_WDM6.H:315
WDM6::eacrc
static constexpr amrex::Real eacrc
Definition: ERF_WDM6.H:190
WDM6::m_bvtg4
amrex::Real m_bvtg4
Definition: ERF_WDM6.H:317
WDM6::m_bvtr7
amrex::Real m_bvtr7
Definition: ERF_WDM6.H:306
WDM6::m_bvtr3
amrex::Real m_bvtr3
Definition: ERF_WDM6.H:306
WDM6::m_pidn0r
amrex::Real m_pidn0r
Definition: ERF_WDM6.H:312
WDM6::avts
static constexpr amrex::Real avts
Definition: ERF_WDM6.H:176
WDM6::xmyu
static constexpr amrex::Real xmyu
Definition: ERF_WDM6.H:175
WDM6::m_bvtr1
amrex::Real m_bvtr1
Definition: ERF_WDM6.H:306
WDM6::m_g1pbs
amrex::Real m_g1pbs
Definition: ERF_WDM6.H:314
WDM6::m_bvtg3
amrex::Real m_bvtg3
Definition: ERF_WDM6.H:317
WDM6::m_bvtr2o5
amrex::Real m_bvtr2o5
Definition: ERF_WDM6.H:307
WDM6::m_avtg
amrex::Real m_avtg
Definition: ERF_WDM6.H:299
WDM6::m_g3pbg
amrex::Real m_g3pbg
Definition: ERF_WDM6.H:318
WDM6::m_bvts1
amrex::Real m_bvts1
Definition: ERF_WDM6.H:313
WDM6::dens_snow
static constexpr amrex::Real dens_snow
Definition: ERF_WDM6.H:191
WDM6::m_g1pbg
amrex::Real m_g1pbg
Definition: ERF_WDM6.H:318
WDM6::m_pacrr
amrex::Real m_pacrr
Definition: ERF_WDM6.H:310
WDM6::r0
static constexpr amrex::Real r0
Definition: ERF_WDM6.H:170
WDM6::m_g5pbgo2
amrex::Real m_g5pbgo2
Definition: ERF_WDM6.H:318
WDM6::m_g4pbs
amrex::Real m_g4pbs
Definition: ERF_WDM6.H:314
WDM6::m_g3pbr
amrex::Real m_g3pbr
Definition: ERF_WDM6.H:308
WDM6::lamdasmax
static constexpr amrex::Real lamdasmax
Definition: ERF_WDM6.H:180
WDM6::m_bvtr2
amrex::Real m_bvtr2
Definition: ERF_WDM6.H:306
WDM6::xncr0
static constexpr amrex::Real xncr0
Definition: ERF_WDM6.H:173
WDM6::peaut
static constexpr amrex::Real peaut
Definition: ERF_WDM6.H:171
WDM6::m_bvtr6
amrex::Real m_bvtr6
Definition: ERF_WDM6.H:306
WDM6::m_g6pbr
amrex::Real m_g6pbr
Definition: ERF_WDM6.H:308
WDM6::xncr1
static constexpr amrex::Real xncr1
Definition: ERF_WDM6.H:174
WDM6::m_bvts4
amrex::Real m_bvts4
Definition: ERF_WDM6.H:313
WDM6::m_lamdagmax
amrex::Real m_lamdagmax
Definition: ERF_WDM6.H:299
WDM6::m_g1pbr
amrex::Real m_g1pbr
Definition: ERF_WDM6.H:308
WDM6::m_g2pbr
amrex::Real m_g2pbr
Definition: ERF_WDM6.H:308
WDM6::m_bvts3
amrex::Real m_bvts3
Definition: ERF_WDM6.H:313
WDM6::m_bvtr5
amrex::Real m_bvtr5
Definition: ERF_WDM6.H:306
WDM6::m_xmmax
amrex::Real m_xmmax
Definition: ERF_WDM6.H:311
WDM6::m_bvtr4
amrex::Real m_bvtr4
Definition: ERF_WDM6.H:306
WDM6::m_bvtg2
amrex::Real m_bvtg2
Definition: ERF_WDM6.H:317
WDM6::m_g5pbso2
amrex::Real m_g5pbso2
Definition: ERF_WDM6.H:314
WDM6::m_g4pbg
amrex::Real m_g4pbg
Definition: ERF_WDM6.H:318
WDM6::m_deng
amrex::Real m_deng
Definition: ERF_WDM6.H:299
WDM6::m_eacrr
amrex::Real m_eacrr
Definition: ERF_WDM6.H:310
WDM6::m_g7pbro2
amrex::Real m_g7pbro2
Definition: ERF_WDM6.H:309
WDM6::m_bvts2
amrex::Real m_bvts2
Definition: ERF_WDM6.H:313
mp_wdm6::rgmma
real(kind=kind_phys) function rgmma(x)
Definition: ERF_module_mp_wdm6.F90:2174
Here is the call graph for this function:

◆ Qmoist_Ptr()

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

Reimplemented from NullMoist.

244  {
245  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
246  return mic_fab_vars[MicVarMap[varIdx]].get();
247  }
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
Here is the call graph for this function:

◆ Qmoist_Restart_Vars()

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

Reimplemented from NullMoist.

256  {
257  a_idx = {0, 1, 2};
258  a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
259  }

◆ Qmoist_Size()

int WDM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

249 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WDM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

251 { return n_qstate_moist_numconc_size; }
WDM6::n_qstate_moist_numconc_size
int n_qstate_moist_numconc_size
Definition: ERF_WDM6.H:274

◆ Qstate_Moist_Size()

int WDM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

250 { return n_qstate_moist_size; }
WDM6::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_WDM6.H:273

◆ Set_dzmin()

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

Reimplemented from NullMoist.

214 { m_dzmin = dz_min; }
WDM6::m_dzmin
amrex::Real m_dzmin
Definition: ERF_WDM6.H:280

◆ Set_Lmask()

void WDM6::Set_Lmask ( amrex::iMultiFab *  )
inlineoverridevirtual

Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this.

Reimplemented from NullMoist.

219 { m_lmask = lmask; }

◆ Update_Micro_Vars() [1/2]

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

Reimplemented from NullMoist.

230  {
231  Copy_State_to_Micro(cons_in);
232  }
WDM6::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Here is the call graph for this function:

◆ Update_Micro_Vars() [2/2]

void WDM6::Update_Micro_Vars ( amrex::MultiFab &  cons_in,
const amrex::MultiFab *  base_state 
)
overridevirtual

Reimplemented from NullMoist.

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

236  {
237  Copy_Micro_to_State(cons_in);
238  }
WDM6::Copy_Micro_to_State
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateWDM6.cpp:9
Here is the call graph for this function:

Member Data Documentation

◆ actk

constexpr amrex::Real WDM6::actk = wdm6_literal(0.6)
staticconstexpr

◆ actr

constexpr amrex::Real WDM6::actr = wdm6_literal(1.5)
staticconstexpr

◆ alpha_wdm6

constexpr amrex::Real WDM6::alpha_wdm6 = wdm6_literal(0.12)
staticconstexpr

◆ avtr

constexpr amrex::Real WDM6::avtr = wdm6_literal(841.9)
staticconstexpr

Referenced by Advance().

◆ avts

constexpr amrex::Real WDM6::avts = wdm6_literal(11.72)
staticconstexpr

◆ bvtr

constexpr amrex::Real WDM6::bvtr = wdm6_literal(0.8)
staticconstexpr

◆ bvts

constexpr amrex::Real WDM6::bvts = wdm6_literal(0.41)
staticconstexpr

◆ dens_snow

constexpr amrex::Real WDM6::dens_snow = wdm6_literal(100.0)
staticconstexpr

◆ di100

constexpr amrex::Real WDM6::di100 = wdm6_literal(1.0e-4)
staticconstexpr

◆ di15

constexpr amrex::Real WDM6::di15 = wdm6_literal(15.0e-6)
staticconstexpr

◆ di2000

constexpr amrex::Real WDM6::di2000 = wdm6_literal(2000.0e-6)
staticconstexpr

◆ di600

constexpr amrex::Real WDM6::di600 = wdm6_literal(6.0e-4)
staticconstexpr

◆ di82

constexpr amrex::Real WDM6::di82 = wdm6_literal(82.0e-6)
staticconstexpr

◆ dicon

constexpr amrex::Real WDM6::dicon = wdm6_literal(11.9)
staticconstexpr

◆ dimax

constexpr amrex::Real WDM6::dimax = wdm6_literal(500.0e-6)
staticconstexpr

◆ dt

amrex::Real WDM6::dt {0.0}
private

◆ dtcldcr

constexpr amrex::Real WDM6::dtcldcr = wdm6_literal(120.0)
staticconstexpr

◆ eacrc

constexpr amrex::Real WDM6::eacrc = wdm6_literal(1.0)
staticconstexpr

◆ lamdacmax

constexpr amrex::Real WDM6::lamdacmax = wdm6_literal(5.0e5)
staticconstexpr

◆ lamdacmin

constexpr amrex::Real WDM6::lamdacmin = wdm6_literal(2.0e4)
staticconstexpr

◆ lamdarmax

constexpr amrex::Real WDM6::lamdarmax = wdm6_literal(5.0e4)
staticconstexpr

◆ lamdarmin

constexpr amrex::Real WDM6::lamdarmin = wdm6_literal(2.0e3)
staticconstexpr

◆ lamdasmax

constexpr amrex::Real WDM6::lamdasmax = wdm6_literal(1.0e5)
staticconstexpr

◆ m_avtg

amrex::Real WDM6::m_avtg {0}
private

◆ m_axis

int WDM6::m_axis {2}
private

Referenced by Define().

◆ m_bvtg

amrex::Real WDM6::m_bvtg {0}
private

◆ m_bvtg1

amrex::Real WDM6::m_bvtg1 {0}
private

◆ m_bvtg2

amrex::Real WDM6::m_bvtg2 {0}
private

◆ m_bvtg3

amrex::Real WDM6::m_bvtg3 {0}
private

◆ m_bvtg4

amrex::Real WDM6::m_bvtg4 {0}
private

◆ m_bvtr1

amrex::Real WDM6::m_bvtr1 {0}
private

◆ m_bvtr2

amrex::Real WDM6::m_bvtr2 {0}
private

◆ m_bvtr2o5

amrex::Real WDM6::m_bvtr2o5 {0}
private

◆ m_bvtr3

amrex::Real WDM6::m_bvtr3 {0}
private

◆ m_bvtr3o5

amrex::Real WDM6::m_bvtr3o5 {0}
private

◆ m_bvtr4

amrex::Real WDM6::m_bvtr4 {0}
private

◆ m_bvtr5

amrex::Real WDM6::m_bvtr5 {0}
private

◆ m_bvtr6

amrex::Real WDM6::m_bvtr6 {0}
private

◆ m_bvtr7

amrex::Real WDM6::m_bvtr7 {0}
private

◆ m_bvts1

amrex::Real WDM6::m_bvts1 {0}
private

◆ m_bvts2

amrex::Real WDM6::m_bvts2 {0}
private

◆ m_bvts3

amrex::Real WDM6::m_bvts3 {0}
private

◆ m_bvts4

amrex::Real WDM6::m_bvts4 {0}
private

◆ m_ccn0

amrex::Real WDM6::m_ccn0 {100.0e6}
private

◆ m_deng

amrex::Real WDM6::m_deng {0}
private

◆ m_detJ_cc

amrex::MultiFab* WDM6::m_detJ_cc {nullptr}
private

◆ m_do_cond

bool WDM6::m_do_cond {true}
private

Referenced by Define().

◆ m_dzmin

amrex::Real WDM6::m_dzmin {0.0}
private

Referenced by Set_dzmin().

◆ m_eacrr

amrex::Real WDM6::m_eacrr {0}
private

◆ m_g1pbg

amrex::Real WDM6::m_g1pbg {0}
private

◆ m_g1pbr

amrex::Real WDM6::m_g1pbr {0}
private

◆ m_g1pbs

amrex::Real WDM6::m_g1pbs {0}
private

◆ m_g2pbr

amrex::Real WDM6::m_g2pbr {0}
private

◆ m_g3pbg

amrex::Real WDM6::m_g3pbg {0}
private

◆ m_g3pbr

amrex::Real WDM6::m_g3pbr {0}
private

◆ m_g3pbs

amrex::Real WDM6::m_g3pbs {0}
private

◆ m_g4pbg

amrex::Real WDM6::m_g4pbg {0}
private

◆ m_g4pbr

amrex::Real WDM6::m_g4pbr {0}
private

◆ m_g4pbs

amrex::Real WDM6::m_g4pbs {0}
private

◆ m_g5pbgo2

amrex::Real WDM6::m_g5pbgo2 {0}
private

◆ m_g5pbr

amrex::Real WDM6::m_g5pbr {0}
private

◆ m_g5pbro2

amrex::Real WDM6::m_g5pbro2 {0}
private

◆ m_g5pbso2

amrex::Real WDM6::m_g5pbso2 {0}
private

◆ m_g6pbr

amrex::Real WDM6::m_g6pbr {0}
private

◆ m_g7pbr

amrex::Real WDM6::m_g7pbr {0}
private

◆ m_g7pbro2

amrex::Real WDM6::m_g7pbro2 {0}
private

◆ m_geom

amrex::Geometry WDM6::m_geom
private

◆ m_hail_opt

bool WDM6::m_hail_opt {false}
private

◆ m_lamdagmax

amrex::Real WDM6::m_lamdagmax {0}
private

◆ m_lmask

amrex::iMultiFab* WDM6::m_lmask {nullptr}
private

Referenced by Set_Lmask().

◆ m_moisture_type

MoistureType WDM6::m_moisture_type {MoistureType::None}
private

Referenced by Define().

◆ m_n0g

amrex::Real WDM6::m_n0g {0}
private

◆ m_pacrc

amrex::Real WDM6::m_pacrc {0}
private

◆ m_pacrg

amrex::Real WDM6::m_pacrg {0}
private

◆ m_pacrr

amrex::Real WDM6::m_pacrr {0}
private

◆ m_pacrs

amrex::Real WDM6::m_pacrs {0}
private

◆ m_pi_wdm6

amrex::Real WDM6::m_pi_wdm6 {0}
private

◆ m_pidn0g

amrex::Real WDM6::m_pidn0g {0}
private

◆ m_pidn0r

amrex::Real WDM6::m_pidn0r {0}
private

◆ m_pidn0s

amrex::Real WDM6::m_pidn0s {0}
private

◆ m_pidnc

amrex::Real WDM6::m_pidnc {0}
private

◆ m_pidnr

amrex::Real WDM6::m_pidnr {0}
private

◆ m_precg1

amrex::Real WDM6::m_precg1 {0}
private

◆ m_precg2

amrex::Real WDM6::m_precg2 {0}
private

◆ m_precr1

amrex::Real WDM6::m_precr1 {0}
private

◆ m_precr2

amrex::Real WDM6::m_precr2 {0}
private

◆ m_precs1

amrex::Real WDM6::m_precs1 {0}
private

◆ m_precs2

amrex::Real WDM6::m_precs2 {0}
private

◆ m_pvtg

amrex::Real WDM6::m_pvtg {0}
private

◆ m_pvtr

amrex::Real WDM6::m_pvtr {0}
private

◆ m_pvtrn

amrex::Real WDM6::m_pvtrn {0}
private

◆ m_pvts

amrex::Real WDM6::m_pvts {0}
private

◆ m_qc0

amrex::Real WDM6::m_qc0 {0}
private

◆ m_qc1

amrex::Real WDM6::m_qc1 {0}
private

◆ m_qck1

amrex::Real WDM6::m_qck1 {0}
private

◆ m_qmoist_size

int WDM6::m_qmoist_size = 3
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

amrex::Real WDM6::m_rdOcp {RdoCp}
private

Referenced by Define().

◆ m_roqimax

amrex::Real WDM6::m_roqimax {0}
private

◆ m_rslopec2max

amrex::Real WDM6::m_rslopec2max {0}
private

◆ m_rslopec3max

amrex::Real WDM6::m_rslopec3max {0}
private

◆ m_rslopecmax

amrex::Real WDM6::m_rslopecmax {0}
private

◆ m_rslopeg2max

amrex::Real WDM6::m_rslopeg2max {0}
private

◆ m_rslopeg3max

amrex::Real WDM6::m_rslopeg3max {0}
private

◆ m_rslopegbmax

amrex::Real WDM6::m_rslopegbmax {0}
private

◆ m_rslopegmax

amrex::Real WDM6::m_rslopegmax {0}
private

◆ m_rsloper2max

amrex::Real WDM6::m_rsloper2max {0}
private

◆ m_rsloper3max

amrex::Real WDM6::m_rsloper3max {0}
private

◆ m_rsloperbmax

amrex::Real WDM6::m_rsloperbmax {0}
private

◆ m_rslopermax

amrex::Real WDM6::m_rslopermax {0}
private

◆ m_rslopes2max

amrex::Real WDM6::m_rslopes2max {0}
private

◆ m_rslopes3max

amrex::Real WDM6::m_rslopes3max {0}
private

◆ m_rslopesbmax

amrex::Real WDM6::m_rslopesbmax {0}
private

◆ m_rslopesmax

amrex::Real WDM6::m_rslopesmax {0}
private

◆ m_xlv1

amrex::Real WDM6::m_xlv1 {0}
private

◆ m_xmmax

amrex::Real WDM6::m_xmmax {0}
private

◆ m_z_phys_nd

amrex::MultiFab* WDM6::m_z_phys_nd {nullptr}
private

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar_WDM6::NumVars> WDM6::mic_fab_vars
private

Referenced by Get_Surface_Precip_Accumulation_Ptrs(), and Qmoist_Ptr().

◆ MicVarMap

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

Referenced by Qmoist_Ptr().

◆ n0r

constexpr amrex::Real WDM6::n0r = wdm6_literal(8.0e6)
staticconstexpr

◆ n0s

constexpr amrex::Real WDM6::n0s = wdm6_literal(2.0e6)
staticconstexpr

◆ n0smax

constexpr amrex::Real WDM6::n0smax = wdm6_literal(1.0e11)
staticconstexpr

◆ n_qstate_moist_numconc_size

int WDM6::n_qstate_moist_numconc_size = 3
private

Referenced by Qstate_Moist_NumConc_Size().

◆ n_qstate_moist_size

int WDM6::n_qstate_moist_size = 9
private

Referenced by Qstate_Moist_Size().

◆ ncmin

constexpr amrex::Real WDM6::ncmin = wdm6_literal(1.0e1)
staticconstexpr

◆ ncrk1

constexpr amrex::Real WDM6::ncrk1 = wdm6_literal(3.03e3)
staticconstexpr

◆ ncrk2

constexpr amrex::Real WDM6::ncrk2 = wdm6_literal(2.59e15)
staticconstexpr

◆ nlev

int WDM6::nlev {0}
private

◆ nrmin

constexpr amrex::Real WDM6::nrmin = wdm6_literal(1.0e-2)
staticconstexpr

◆ peaut

constexpr amrex::Real WDM6::peaut = wdm6_literal(0.55)
staticconstexpr

◆ pfrz1

constexpr amrex::Real WDM6::pfrz1 = wdm6_literal(100.0)
staticconstexpr

◆ pfrz2

constexpr amrex::Real WDM6::pfrz2 = wdm6_literal(0.66)
staticconstexpr

◆ qcrmin

constexpr amrex::Real WDM6::qcrmin = wdm6_literal(1.0e-9)
staticconstexpr

◆ qs0

constexpr amrex::Real WDM6::qs0 = wdm6_literal(6.0e-4)
staticconstexpr

◆ r0

constexpr amrex::Real WDM6::r0 = wdm6_literal(0.8e-5)
staticconstexpr

◆ satmax

constexpr amrex::Real WDM6::satmax = wdm6_literal(1.0048)
staticconstexpr

◆ xmyu

constexpr amrex::Real WDM6::xmyu = wdm6_literal(1.718e-5)
staticconstexpr

◆ xncr

constexpr amrex::Real WDM6::xncr = wdm6_literal(3.0e8)
staticconstexpr

◆ xncr0

constexpr amrex::Real WDM6::xncr0 = wdm6_literal(5.0e7)
staticconstexpr

◆ xncr1

constexpr amrex::Real WDM6::xncr1 = wdm6_literal(5.0e8)
staticconstexpr

◆ zhi

int WDM6::zhi {0}
private

◆ zlo

int WDM6::zlo {0}
private
The documentation for this class was generated from the following files: