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_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &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
 
virtual void Update_Micro_Vars (amrex::MultiFab &)
 
virtual void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *)
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *)
 
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 &)
 
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}
 
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}
 

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ WDM6()

WDM6::WDM6 ( )
inline
61 {}

◆ ~WDM6()

virtual WDM6::~WDM6 ( )
virtualdefault

Member Function Documentation

◆ Advance()

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

Reimplemented from NullMoist.

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

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro()

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

Reimplemented from NullMoist.

84 {
85  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
86  // Match Morrison behavior: refresh microphysics ghost zones from state.
87  const auto& box3d = mfi.growntilebox();
88  auto states = cons_in.array(mfi);
89 
90  auto rho = mic_fab_vars[MicVar_WDM6::rho]->array(mfi);
91  auto theta = mic_fab_vars[MicVar_WDM6::theta]->array(mfi);
92  auto tabs = mic_fab_vars[MicVar_WDM6::tabs]->array(mfi);
93  auto pres = mic_fab_vars[MicVar_WDM6::pres]->array(mfi);
94 
95  auto qv = mic_fab_vars[MicVar_WDM6::qv]->array(mfi);
96  auto qc = mic_fab_vars[MicVar_WDM6::qc]->array(mfi);
97  auto qi = mic_fab_vars[MicVar_WDM6::qi]->array(mfi);
98  auto qr = mic_fab_vars[MicVar_WDM6::qr]->array(mfi);
99  auto qs = mic_fab_vars[MicVar_WDM6::qs]->array(mfi);
100  auto qg = mic_fab_vars[MicVar_WDM6::qg]->array(mfi);
101 
102  auto nn = mic_fab_vars[MicVar_WDM6::nn]->array(mfi);
103  auto nc = mic_fab_vars[MicVar_WDM6::nc]->array(mfi);
104  auto nr = mic_fab_vars[MicVar_WDM6::nr]->array(mfi);
105 
106  const Real ccn0_local = m_ccn0; // CCN concentration in #/m³
107 
108  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
109  rho(i,j,k) = states(i,j,k,Rho_comp);
110  theta(i,j,k) = states(i,j,k,RhoTheta_comp) / states(i,j,k,Rho_comp);
111 
112  qv(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ1_comp) / states(i,j,k,Rho_comp));
113  qc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ2_comp) / states(i,j,k,Rho_comp));
114  qi(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ3_comp) / states(i,j,k,Rho_comp));
115  qr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ4_comp) / states(i,j,k,Rho_comp));
116  qs(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ5_comp) / states(i,j,k,Rho_comp));
117  qg(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ6_comp) / states(i,j,k,Rho_comp));
118 
119  // Number concentrations
120  nc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ7_comp) / states(i,j,k,Rho_comp));
121  nr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ9_comp) / states(i,j,k,Rho_comp));
122 
123  // nn: seed the CCN reservoir where state carries nothing yet,
124  // otherwise read it like every other state variable.
125  //
126  // This mirrors the Fortran, which does its seeding in the per-step
127  // driver wdm6 (:219) rather than in wdm6init, so this is the
128  // faithful location. The Fortran gates on itimestep == 1; RhoQ8 == 0
129  // is equivalent for every reachable case and needs no step plumbed
130  // through the Init interface. Cold start and restart from a step-0
131  // checkpoint both have RhoQ8 == 0 and seed, matching itimestep == 1;
132  // every later restart has RhoQ8 > 0 and restores. A genuine zero
133  // cannot occur after initialization because Advance bounds nn well
134  // above zero (:602). Note the Fortran block is in fact unreachable
135  // from ERF: the bridge enters wdm62D directly via mp_wdm6_run and
136  // skips the wdm6 wrapper, so BOTH legs take their nn from here.
137  const Real nn_state = states(i,j,k,RhoQ8_comp) / states(i,j,k,Rho_comp);
138  nn(i,j,k) = (nn_state > Real(0.0)) ? nn_state : ccn0_local;
139 
140  // WDM6: DO NOT enforce nc/nr minimums here!
141  // WRF starts with nc=0, nr=0 and lets CCN activation build nc naturally during
142  // the first microphysics call. Enforcing nc=10, nr=0.01 here prevents proper
143  // activation from nn. Minimums are enforced in Advance() right before physics,
144  // not during state copying.
145 
146  tabs(i,j,k) = getTgivenRandRTh(states(i,j,k,Rho_comp),
147  states(i,j,k,RhoTheta_comp),
148  qv(i,j,k));
149  pres(i,j,k) = getPgivenRTh(states(i,j,k,RhoTheta_comp), qv(i,j,k));
150  });
151  }
152 
153  // After first Copy_State_to_Micro, nn has been preserved from Init().
154  // DON'T clear the flag yet - wait until after Copy_Micro_to_State writes nn to state!
155 }
getTgivenRandRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:46
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
MicVar_Kess::tabs
@ tabs
Definition: ERF_Kessler.H:26
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:27

Referenced by Update_Micro_Vars().

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

◆ Define()

void WDM6::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

65  {
66  m_moisture_type = sc.moisture_type;
67  m_axis = sc.ave_plane;
68  m_do_cond = (!sc.uses_shoc_family());
69  }
WDM6::m_axis
int m_axis
Definition: ERF_WDM6.H:227
WDM6::m_moisture_type
MoistureType m_moisture_type
Definition: ERF_WDM6.H:229
WDM6::m_do_cond
bool m_do_cond
Definition: ERF_WDM6.H:228
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2124
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2141
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:2063
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.

207  {
208  SurfacePrecipAccumulationSources sources;
209  sources.total = {mic_fab_vars[MicVar_WDM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
210  sources.snow = {mic_fab_vars[MicVar_WDM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
211  sources.graupel = {mic_fab_vars[MicVar_WDM6::graup_accum].get(), rhoh2o / amrex::Real(1000.0)};
212  return sources;
213  }
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.

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

◆ initialize_coeffs()

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

188  {
189  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
190  return mic_fab_vars[MicVarMap[varIdx]].get();
191  }
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.

200  {
201  a_idx = {0, 1, 2};
202  a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
203  }

◆ Qmoist_Size()

int WDM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

193 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WDM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

195 { return n_qstate_moist_numconc_size; }
WDM6::n_qstate_moist_numconc_size
int n_qstate_moist_numconc_size
Definition: ERF_WDM6.H:218

◆ Qstate_Moist_Size()

int WDM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

194 { return n_qstate_moist_size; }
WDM6::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_WDM6.H:217

◆ Set_dzmin()

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

Reimplemented from NullMoist.

161 { m_dzmin = dz_min; }
WDM6::m_dzmin
amrex::Real m_dzmin
Definition: ERF_WDM6.H:224

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

166 { m_lmask = lmask; }

◆ Update_Micro_Vars() [1/3]

virtual void NullMoist::Update_Micro_Vars
inline
36 { }

◆ Update_Micro_Vars() [2/3]

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

Reimplemented from NullMoist.

174  {
175  Copy_State_to_Micro(cons_in);
176  }
WDM6::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitWDM6.cpp:83
Here is the call graph for this function:

◆ Update_Micro_Vars() [3/3]

virtual void NullMoist::Update_Micro_Vars
inline
43  {
44  Update_Micro_Vars(cons_in);
45  }
WDM6::Update_Micro_Vars
virtual void Update_Micro_Vars(amrex::MultiFab &)
Definition: ERF_NullMoist.H:36

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

180  {
181  Copy_Micro_to_State(cons_in);
182  }
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_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: