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

#include <ERF_WSM6.H>

Inheritance diagram for WSM6:
Collaboration diagram for WSM6:

Public Member Functions

 WSM6 ()
 
virtual ~WSM6 ()=default
 
void Define (SolverChoice &sc) override
 
void Init (const amrex::MultiFab &cons_in, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &z_phys_nd, std::unique_ptr< amrex::MultiFab > &detJ_cc) override
 
void Set_dzmin (const amrex::Real dz_min) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in, const amrex::MultiFab *base_state)
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *base_state) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &solverChoice) override
 
const SedimentationFlux * Sed_Flux () const override
 The sedimentation face fluxes of the most recent Advance, or null for a scheme that does not sediment (see ERF_SedimentationFlux.H) More...
 
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
 
void AdvanceColumns (const amrex::Real &dt_advance, const SolverChoice &sc)
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 
virtual void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const
 
virtual void SetCurrentLevel (const int &lev)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_Lmask (amrex::iMultiFab *)
 Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this. More...
 
virtual void Set_Sed_Inflow (const amrex::MultiFab *inflow)
 The inflow through the top of a fine patch for the next Advance, or null (see ERF_SedimentationFlux.H); the caller keeps it alive for that Advance. More...
 
virtual void Set_RealWidth (const int)
 

Static Public Attributes

static constexpr amrex::Real dtcldcr = amrex::Real(120.0)
 
static constexpr amrex::Real n0r = amrex::Real(8.0e6)
 
static constexpr amrex::Real avtr = amrex::Real(841.9)
 
static constexpr amrex::Real bvtr = amrex::Real(0.8)
 
static constexpr amrex::Real r0 = amrex::Real(0.8e-5)
 
static constexpr amrex::Real peaut = amrex::Real(0.55)
 
static constexpr amrex::Real xncr = amrex::Real(3.0e8)
 
static constexpr amrex::Real xmyu = amrex::Real(1.718e-5)
 
static constexpr amrex::Real avts = amrex::Real(11.72)
 
static constexpr amrex::Real bvts = amrex::Real(0.41)
 
static constexpr amrex::Real lamdarmax = amrex::Real(8.0e4)
 
static constexpr amrex::Real lamdasmax = amrex::Real(1.0e5)
 
static constexpr amrex::Real dicon = amrex::Real(11.9)
 
static constexpr amrex::Real dimax = amrex::Real(500.0e-6)
 
static constexpr amrex::Real pfrz1 = amrex::Real(100.0)
 
static constexpr amrex::Real pfrz2 = amrex::Real(0.66)
 
static constexpr amrex::Real qcrmin = amrex::Real(1.0e-9)
 
static constexpr amrex::Real eacrc = amrex::Real(1.0)
 
static constexpr amrex::Real dens_snow = amrex::Real(100.0)
 
static constexpr amrex::Real qs0 = amrex::Real(6.0e-4)
 
static constexpr amrex::Real n0smax = amrex::Real(1.0e11)
 
static constexpr amrex::Real n0s = amrex::Real(2.0e6)
 
static constexpr amrex::Real alpha_wsm6 = amrex::Real(0.12)
 

Private Types

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

Private Member Functions

void initialize_coeffs ()
 

Private Attributes

int m_qmoist_size = 4
 
int n_qstate_moist_size = 6
 
int n_qstate_moist_numconc_size = 0
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::Real dt {0.0}
 
amrex::Real m_dzmin {0.0}
 
int nlev {0}
 
int zlo {0}
 
int zhi {0}
 
int m_axis {2}
 
bool m_do_cond {true}
 
amrex::Real m_rdOcp {RdoCp}
 
MoistureType m_moisture_type {MoistureType::None}
 
amrex::MultiFab * m_z_phys_nd {nullptr}
 
amrex::MultiFab * m_detJ_cc {nullptr}
 
amrex::Array< FabPtr, MicVar_WSM6::NumVars > mic_fab_vars
 
MicrophysicsColumns m_columns
 
SedimentationFlux m_sed
 
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_wsm6 {0}
 
amrex::Real m_xlv1 {0}
 
amrex::Real m_qc0 {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_bvtr6 {0}
 
amrex::Real m_g1pbr {0}
 
amrex::Real m_g3pbr {0}
 
amrex::Real m_g4pbr {0}
 
amrex::Real m_g5pbro2 {0}
 
amrex::Real m_g6pbr {0}
 
amrex::Real m_pvtr {0}
 
amrex::Real m_eacrr {0}
 
amrex::Real m_pacrr {0}
 
amrex::Real m_precr1 {0}
 
amrex::Real m_precr2 {0}
 
amrex::Real m_roqimax {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_pidn0r {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_rslopermax {0}
 
amrex::Real m_rslopesmax {0}
 
amrex::Real m_rslopegmax {0}
 
amrex::Real m_rsloperbmax {0}
 
amrex::Real m_rslopesbmax {0}
 
amrex::Real m_rslopegbmax {0}
 
amrex::Real m_rsloper2max {0}
 
amrex::Real m_rslopes2max {0}
 
amrex::Real m_rslopeg2max {0}
 
amrex::Real m_rsloper3max {0}
 
amrex::Real m_rslopes3max {0}
 
amrex::Real m_rslopeg3max {0}
 

Additional Inherited Members

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

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ WSM6()

WSM6::WSM6 ( )
inline
132 {}

◆ ~WSM6()

virtual WSM6::~WSM6 ( )
virtualdefault

Member Function Documentation

◆ Advance()

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

Reimplemented from NullMoist.

862 {
863  m_sed.reset(dt_advance);
864 
865  if (!m_columns.needs_gather()) {
866  AdvanceColumns(dt_advance, sc);
867  return;
868  }
869 
870  // Some column is split between boxes stacked in z. Sedimentation treats each box's
871  // k range as the whole column, so run on a copy of the fields in which every column
872  // lies in one box, then copy the result back.
873  amrex::Array<FabPtr, MicVar_WSM6::NumVars> saved;
874  m_columns.gather(mic_fab_vars, saved);
875  MultiFab* z_phys_nd_saved = m_z_phys_nd;
876  auto z_phys_nd_col = m_columns.gather_nodal(m_z_phys_nd);
877  m_z_phys_nd = z_phys_nd_col.get();
878  auto sed_saved = m_sed.down;
879  m_sed.down = m_columns.gather_one(*sed_saved);
880  const MultiFab* inflow_saved = m_sed_inflow;
881  auto inflow_col = (m_sed_inflow) ? m_columns.gather_one(*m_sed_inflow) : nullptr;
882  m_sed_inflow = inflow_col.get();
883 
884  AdvanceColumns(dt_advance, sc);
885 
886  m_sed_inflow = inflow_saved;
887  sed_saved->ParallelCopy(*m_sed.down);
888  m_sed.down = sed_saved;
889  m_z_phys_nd = z_phys_nd_saved;
891 }
void gather(std::array< FabPtr, N > &vars, std::array< FabPtr, N > &saved) const
Swap each non-null field for a copy on the column layout; the originals go to saved.
Definition: ERF_MicrophysicsColumns.H:48
FabPtr gather_one(const amrex::MultiFab &src) const
Copy of a field (cell-centered or z-face; with its ghost cells) on the column layout.
Definition: ERF_MicrophysicsColumns.cpp:100
void scatter(std::array< FabPtr, N > &vars, std::array< FabPtr, N > &saved) const
Copy the column-layout fields back into the originals and restore them.
Definition: ERF_MicrophysicsColumns.H:58
std::unique_ptr< amrex::MultiFab > gather_nodal(const amrex::MultiFab *src) const
Copy of a nodal field (e.g. z_phys_nd) on the column layout, or null for null input.
Definition: ERF_MicrophysicsColumns.cpp:116
bool needs_gather() const
True when some column is split between boxes, so the fields must be gathered.
Definition: ERF_MicrophysicsColumns.H:41
const amrex::MultiFab * m_sed_inflow
Definition: ERF_NullMoist.H:168
amrex::MultiFab * m_z_phys_nd
Definition: ERF_WSM6.H:252
MicrophysicsColumns m_columns
Definition: ERF_WSM6.H:258
void AdvanceColumns(const amrex::Real &dt_advance, const SolverChoice &sc)
Definition: ERF_AdvanceWSM6.cpp:894
SedimentationFlux m_sed
Definition: ERF_WSM6.H:261
amrex::Array< FabPtr, MicVar_WSM6::NumVars > mic_fab_vars
Definition: ERF_WSM6.H:255
std::shared_ptr< amrex::MultiFab > down
Definition: ERF_SedimentationFlux.H:39
void reset(amrex::Real a_dt)
Zero the fluxes at the start of an Advance of length a_dt.
Definition: ERF_SedimentationFlux.H:57

◆ AdvanceColumns()

void WSM6::AdvanceColumns ( const amrex::Real &  dt_advance,
const SolverChoice &  sc 
)
896 {
897  dt = dt_advance;
898 
899  int microphysics_debug = 0;
900  std::string micro_diag_mode = "canonical";
901  std::vector<std::string> micro_diag_tags = {"standing"};
902  std::vector<std::string> micro_diag_expr = {"standing"};
903  std::vector<std::string> micro_diag_store = {"standing"};
904  std::vector<int> micro_diag_target_column;
905  {
906  amrex::ParmParse pp("erf");
907  pp.queryAdd("microphysics_debug", microphysics_debug);
908  pp.queryAdd("micro_diag_mode", micro_diag_mode);
909  pp.queryarr("micro_diag_tags", micro_diag_tags);
910  pp.queryarr("micro_diag_expr", micro_diag_expr);
911  pp.queryarr("micro_diag_store", micro_diag_store);
912  pp.queryarr("micro_diag_target_column", micro_diag_target_column);
913  }
914  microphysics_debug = std::max(0, std::min(2, microphysics_debug));
915 #ifdef ERF_USE_WSM6_FORT
916  const std::string micro_diag_mode_lower = [&]{ std::string m = micro_diag_mode; std::transform(m.begin(), m.end(), m.begin(), [](unsigned char c){ return static_cast<char>(std::tolower(c)); }); return m; }();
917  const bool micro_diag_forensic = (micro_diag_mode_lower == "forensic" || micro_diag_mode_lower == "both");
918  const int microphysics_debug_bridge = micro_diag_forensic
919  ? microphysics_debug
920  : std::min(microphysics_debug, 1);
921  bool use_wsm6_cpp_answer = false;
922  { amrex::ParmParse pp("erf");
923  pp.queryAdd("use_wsm6_cpp_answer", use_wsm6_cpp_answer); }
924  bool run_wsm6_fort = !use_wsm6_cpp_answer;
925 
926  // The Fortran bridge takes each box's k range as the whole column (the box top as the
927  // model top, its bottom as the ground), so it needs boxes that do span it.
929  "WSM6: the Fortran bridge (erf.use_wsm6_cpp_answer = false) requires every box to span "
930  "the domain in z (no grids split in z, and no refined levels short of the domain top or bottom)");
931 
932  static bool wsm6_inited = false;
933 
934  // Minimal phase-1 initialization for single-moment WSM6.
935  if (!wsm6_inited) {
936  constexpr double den0 = 1.28; // Standard dry-air density (kg/m^3)
937  constexpr double denr = static_cast<double>(rhoh2o);
938  constexpr double dens = static_cast<double>(rhos);
939  constexpr double cl = static_cast<double>(Cp_l);
940  constexpr double cpv = static_cast<double>(Cp_v);
941  constexpr int hail_opt = 0; // Graupel mode
942  mp_wsm6_init_c(den0, denr, dens, cl, cpv, hail_opt);
943  wsm6_inited = true;
944  }
945 #endif
946 
947  constexpr double g = static_cast<double>(CONST_GRAV);
948  constexpr double cpd = static_cast<double>(Cp_d);
949  constexpr double cpv = static_cast<double>(Cp_v);
950  constexpr double rd = static_cast<double>(R_d);
951  constexpr double rv = static_cast<double>(R_v);
952  constexpr double t0c = 273.15;
953  constexpr double ep1 = static_cast<double>(R_v / R_d - one);
954  amrex::ignore_unused(g, rd, ep1);
955  constexpr double ep2 = static_cast<double>(R_d / R_v);
956  constexpr double qmin = 1.0e-12;
957  constexpr double xls = static_cast<double>(lsub);
958  constexpr double xlv0 = static_cast<double>(lat_vap);
959  constexpr double xlf0 = static_cast<double>(lat_ice);
960  constexpr double den0 = 1.28;
961  constexpr double denr = static_cast<double>(rhoh2o);
962  constexpr double cliq = static_cast<double>(Cp_l);
963  constexpr double cice = 2106.0;
964  constexpr double psat = 610.78;
965  for (MFIter mfi(*mic_fab_vars[MicVar_WSM6::qv], TileNoZ()); mfi.isValid(); ++mfi) {
966  const Box box = mfi.tilebox();
967  const Box fab_box = mfi.fabbox();
968 
969  auto const& t_arr = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
970  auto const& qv_arr = mic_fab_vars[MicVar_WSM6::qv]->array(mfi);
971  auto const& qc_arr = mic_fab_vars[MicVar_WSM6::qc]->array(mfi);
972  auto const& qi_arr = mic_fab_vars[MicVar_WSM6::qi]->array(mfi);
973  auto const& qr_arr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi);
974  auto const& qs_arr = mic_fab_vars[MicVar_WSM6::qs]->array(mfi);
975  auto const& qg_arr = mic_fab_vars[MicVar_WSM6::qg]->array(mfi);
976  auto const& den_arr = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
977  auto const& p_arr = mic_fab_vars[MicVar_WSM6::pres]->array(mfi);
978  auto rain_arr = mic_fab_vars[MicVar_WSM6::rain_accum]->array(mfi);
979  auto snow_arr = mic_fab_vars[MicVar_WSM6::snow_accum]->array(mfi);
980  auto graupel_arr = mic_fab_vars[MicVar_WSM6::graupel_accum]->array(mfi);
981  auto precip_arr = mic_fab_vars[MicVar_WSM6::precip_accum]->array(mfi);
982  auto const& sed_arr = m_sed.down->array(mfi);
983 
984  const int ilo = box.smallEnd(0);
985  const int ihi = box.bigEnd(0);
986  const int jlo = box.smallEnd(1);
987  const int jhi = box.bigEnd(1);
988  const int klo = box.smallEnd(2);
989  const int khi = box.bigEnd(2);
990  const bool has_target_override = (micro_diag_target_column.size() == 2);
991  const int diag_i = has_target_override ? micro_diag_target_column[0] : ilo;
992  const int diag_j = has_target_override ? micro_diag_target_column[1] : jlo;
993 
994  const int imlo = fab_box.smallEnd(0);
995  const int imhi = fab_box.bigEnd(0);
996  const int jmlo = fab_box.smallEnd(1);
997  const int jmhi = fab_box.bigEnd(1);
998  const int kmlo = fab_box.smallEnd(2);
999  const int kmhi = fab_box.bigEnd(2);
1000  amrex::ignore_unused(ihi, jhi, diag_i, diag_j, imlo, imhi, jmlo, jmhi, kmlo, kmhi);
1001 
1002 #if defined(ERF_USE_WSM6_FORT) && defined(AMREX_USE_GPU)
1003  Arena* Arena_Used = (run_wsm6_fort) ? The_Pinned_Arena() : The_Async_Arena();
1004 #else
1005  Arena* Arena_Used = The_Async_Arena();
1006 #endif
1007 
1008  const Real dz_val = m_geom.CellSize(2);
1009  FArrayBox delz_fab(fab_box, 1, Arena_Used);
1010  auto const& delz_arr = delz_fab.array();
1011  ParallelFor(fab_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1012  delz_arr(i,j,k) = dz_val;
1013  });
1014 
1015  const Array4<const Real> z_arr = (m_z_phys_nd) ? m_z_phys_nd->const_array(mfi) : Array4<const Real> {};
1016  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1017  delz_arr(i,j,k) = (z_arr) ? Real(0.25) * ( (z_arr(i ,j ,k+1) - z_arr(i ,j ,k))
1018  + (z_arr(i+1,j ,k+1) - z_arr(i+1,j ,k))
1019  + (z_arr(i ,j+1,k+1) - z_arr(i ,j+1,k))
1020  + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz_val;
1021  });
1022 
1023  Box box2d(box);
1024  box2d.makeSlab(2, 0);
1025  Box fab_box2d(fab_box);
1026  fab_box2d.makeSlab(2, 0);
1027 
1028  // Fortran bridge uses ims:ime, jms:jme storage bounds; these buffers must
1029  // therefore be allocated on fab_box extents even if C++ kernels only
1030  // update the valid tile slab (box2d).
1031  FArrayBox rainncv_fab(fab_box2d, 1, Arena_Used);
1032  FArrayBox sr_fab(fab_box2d, 1, Arena_Used);
1033  FArrayBox snowncv_fab(fab_box2d, 1, Arena_Used);
1034  FArrayBox graupelncv_fab(fab_box2d, 1, Arena_Used);
1035  FArrayBox rainacc_fab(fab_box2d, 1, Arena_Used);
1036  FArrayBox snowacc_fab(fab_box2d, 1, Arena_Used);
1037  FArrayBox graupacc_fab(fab_box2d, 1, Arena_Used);
1038  FArrayBox rainonly_fab(fab_box2d, 1, Arena_Used);
1039  FArrayBox snowonly_fab(fab_box2d, 1, Arena_Used);
1040 
1041  auto const& rainncv_arr = rainncv_fab.array();
1042  auto const& sr_arr = sr_fab.array();
1043  auto const& snowncv_arr = snowncv_fab.array();
1044  auto const& graupelncv_arr = graupelncv_fab.array();
1045  auto const& rainacc_arr = rainacc_fab.array();
1046  auto const& snowacc_arr = snowacc_fab.array();
1047  auto const& graupacc_arr = graupacc_fab.array();
1048  auto const& rainonly_arr = rainonly_fab.array();
1049  auto const& snowonly_arr = snowonly_fab.array();
1050  ParallelFor(fab_box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1051  rainncv_arr(i,j,k) = Real(0.0);
1052  sr_arr(i,j,k) = Real(0.0);
1053  snowncv_arr(i,j,k) = Real(0.0);
1054  graupelncv_arr(i,j,k) = Real(0.0);
1055  rainacc_arr(i,j,k) = Real(0.0);
1056  snowacc_arr(i,j,k) = Real(0.0);
1057  graupacc_arr(i,j,k) = Real(0.0);
1058  rainonly_arr(i,j,k) = Real(0.0);
1059  snowonly_arr(i,j,k) = Real(0.0);
1060  });
1061  // The WRF-style accumulators start this step at zero, so afterwards they hold
1062  // what fell during it: rainacc the total, snowacc the snow and cloud ice,
1063  // graupacc the graupel, and rainonly and snowonly the rain and snow alone.
1064  // They are added to the running accumulations at the end.
1065 
1066 #ifdef ERF_USE_WSM6_FORT
1067  if (run_wsm6_fort) {
1068  // Host-only Fortran reads mic_fab_vars via dataPtr(), so GPU writes
1069  // must complete before crossing the language boundary.
1070  Gpu::streamSynchronize();
1071  mp_wsm6_run_c(
1072  t_arr.dataPtr(),
1073  qv_arr.dataPtr(), qc_arr.dataPtr(), qi_arr.dataPtr(),
1074  qr_arr.dataPtr(), qs_arr.dataPtr(), qg_arr.dataPtr(),
1075  den_arr.dataPtr(), p_arr.dataPtr(), delz_arr.dataPtr(),
1076  static_cast<double>(dt), g, cpd, cpv, rd, rv, t0c, ep1, ep2, qmin,
1077  xls, xlv0, xlf0, den0, denr, cliq, cice, psat,
1078  rainacc_arr.dataPtr(), rainncv_arr.dataPtr(), sr_arr.dataPtr(),
1079  snowacc_arr.dataPtr(), snowncv_arr.dataPtr(),
1080  graupacc_arr.dataPtr(), graupelncv_arr.dataPtr(),
1081  rainonly_arr.dataPtr(), snowonly_arr.dataPtr(),
1082  imlo, imhi, jmlo, jmhi, kmlo, kmhi,
1083  ilo, ihi, jlo, jhi, klo, khi, microphysics_debug_bridge, diag_i, diag_j);
1084  } else {
1085 #endif
1086  // --- Phase 4 native C++ kernel ---
1087 
1088  // box2d for 1D per-column arrays (already defined above)
1089  // delqrs1/2/3, delqi: surface precipitation flux accumulators
1090  FArrayBox delqrs1_fab(box2d,1, Arena_Used);
1091  FArrayBox delqrs2_fab(box2d,1, Arena_Used);
1092  FArrayBox delqrs3_fab(box2d,1, Arena_Used);
1093  FArrayBox delqi_fab(box2d,1, Arena_Used);
1094  FArrayBox tstepsnow_fab(box2d,1, Arena_Used);
1095  FArrayBox tstepgraup_fab(box2d,1, Arena_Used);
1096  auto const& delqrs1_arr = delqrs1_fab.array();
1097  auto const& delqrs2_arr = delqrs2_fab.array();
1098  auto const& delqrs3_arr = delqrs3_fab.array();
1099  auto const& delqi_arr = delqi_fab.array();
1100  auto const& tstepsnow_arr = tstepsnow_fab.array();
1101  auto const& tstepgraup_arr = tstepgraup_fab.array();
1102 
1103  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1104  delqrs1_arr(i,j,k) = Real(0.0);
1105  delqrs2_arr(i,j,k) = Real(0.0);
1106  delqrs3_arr(i,j,k) = Real(0.0);
1107  delqi_arr(i,j,k) = Real(0.0);
1108  tstepsnow_arr(i,j,k) = Real(0.0);
1109  tstepgraup_arr(i,j,k) = Real(0.0);
1110  });
1111  // 3D working FABs
1112  FArrayBox denfac_fab(fab_box,1, Arena_Used); FArrayBox xni_fab(fab_box,1, Arena_Used);
1113  FArrayBox cpm_fab(fab_box,1, Arena_Used); FArrayBox xl_fab(fab_box,1, Arena_Used);
1114  FArrayBox qsatw_fab(fab_box,1, Arena_Used); FArrayBox qsati_fab(fab_box,1, Arena_Used);
1115  FArrayBox rhw_fab(fab_box,1, Arena_Used); FArrayBox rhi_fab(fab_box,1, Arena_Used);
1116  FArrayBox den_tmp_fab(fab_box,1, Arena_Used); FArrayBox delz_tmp_fab(fab_box,1, Arena_Used);
1117  FArrayBox n0sfac_fab(fab_box,1, Arena_Used);
1118  FArrayBox qrs_tmp_r_fab(fab_box,1, Arena_Used); FArrayBox qrs_tmp_s_fab(fab_box,1, Arena_Used);
1119  FArrayBox qrs_tmp_g_fab(fab_box,1, Arena_Used);
1120  FArrayBox rslope_r_fab(fab_box,1, Arena_Used); FArrayBox rslope_s_fab(fab_box,1, Arena_Used);
1121  FArrayBox rslope_g_fab(fab_box,1, Arena_Used);
1122  FArrayBox rslope2_r_fab(fab_box,1, Arena_Used); FArrayBox rslope2_s_fab(fab_box,1, Arena_Used);
1123  FArrayBox rslope2_g_fab(fab_box,1, Arena_Used);
1124  FArrayBox rslope3_r_fab(fab_box,1, Arena_Used); FArrayBox rslope3_s_fab(fab_box,1, Arena_Used);
1125  FArrayBox rslope3_g_fab(fab_box,1, Arena_Used);
1126  FArrayBox rslopeb_r_fab(fab_box,1, Arena_Used); FArrayBox rslopeb_s_fab(fab_box,1, Arena_Used);
1127  FArrayBox rslopeb_g_fab(fab_box,1, Arena_Used);
1128  FArrayBox work1_r_fab(fab_box,1, Arena_Used); FArrayBox work1_s_fab(fab_box,1, Arena_Used);
1129  FArrayBox work1_g_fab(fab_box,1, Arena_Used);
1130  FArrayBox work2_fab(fab_box,1, Arena_Used); FArrayBox workdiffw_fab(fab_box,1, Arena_Used);
1131  FArrayBox workdiffi_fab(fab_box,1, Arena_Used);
1132  FArrayBox workr_fab(fab_box,1, Arena_Used); FArrayBox worka_fab(fab_box,1, Arena_Used);
1133  FArrayBox work1c_fab(fab_box,1, Arena_Used);
1134  FArrayBox denqrs1_fab(fab_box,1, Arena_Used); FArrayBox denqrs2_fab(fab_box,1, Arena_Used);
1135  FArrayBox denqrs3_fab(fab_box,1, Arena_Used); FArrayBox denqci_fab(fab_box,1, Arena_Used);
1136  FArrayBox fall_r_fab(fab_box,1, Arena_Used); FArrayBox fall_s_fab(fab_box,1, Arena_Used);
1137  FArrayBox fall_g_fab(fab_box,1, Arena_Used); FArrayBox fallc_fab(fab_box,1, Arena_Used);
1138  FArrayBox qsum_fab(fab_box,1, Arena_Used);
1139  FArrayBox nislfv_r_diag_fab(fab_box,6, Arena_Used);
1140  FArrayBox nislfv_sg_diag_fab(fab_box,6, Arena_Used);
1141  FArrayBox sed_cell_scratch_fab(fab_box, WSM6SedCellScratch::NumComps, Arena_Used);
1142  Box sed_node_box = amrex::surroundingNodes(fab_box, 2);
1143  FArrayBox sed_node_scratch_fab(sed_node_box, WSM6SedNodeScratch::NumComps, Arena_Used);
1144  // process rates
1145  FArrayBox praut_fab(fab_box,1, Arena_Used); FArrayBox pracw_fab(fab_box,1, Arena_Used);
1146  FArrayBox prevp_fab(fab_box,1, Arena_Used); FArrayBox psdep_fab(fab_box,1, Arena_Used);
1147  FArrayBox pgdep_fab(fab_box,1, Arena_Used); FArrayBox psaut_fab(fab_box,1, Arena_Used);
1148  FArrayBox pgaut_fab(fab_box,1, Arena_Used); FArrayBox praci_fab(fab_box,1, Arena_Used);
1149  FArrayBox piacr_fab(fab_box,1, Arena_Used); FArrayBox psaci_fab(fab_box,1, Arena_Used);
1150  FArrayBox psacw_fab(fab_box,1, Arena_Used); FArrayBox pgacw_fab(fab_box,1, Arena_Used);
1151  FArrayBox pgaci_fab(fab_box,1, Arena_Used); FArrayBox paacw_fab(fab_box,1, Arena_Used);
1152  FArrayBox pracs_fab(fab_box,1, Arena_Used); FArrayBox psacr_fab(fab_box,1, Arena_Used);
1153  FArrayBox pgacr_fab(fab_box,1, Arena_Used); FArrayBox pgacs_fab(fab_box,1, Arena_Used);
1154  FArrayBox pigen_fab(fab_box,1, Arena_Used); FArrayBox pidep_fab(fab_box,1, Arena_Used);
1155  FArrayBox pcond_fab(fab_box,1, Arena_Used); FArrayBox psmlt_fab(fab_box,1, Arena_Used);
1156  FArrayBox pgmlt_fab(fab_box,1, Arena_Used); FArrayBox pseml_fab(fab_box,1, Arena_Used);
1157  FArrayBox pgeml_fab(fab_box,1, Arena_Used); FArrayBox psevp_fab(fab_box,1, Arena_Used);
1158  FArrayBox pgevp_fab(fab_box,1, Arena_Used);
1159  FArrayBox pimlt_fab(fab_box,1, Arena_Used); FArrayBox pihmf_fab(fab_box,1, Arena_Used);
1160  FArrayBox pihtf_fab(fab_box,1, Arena_Used); FArrayBox pgfrz_fab(fab_box,1, Arena_Used);
1161 
1162  auto const& denfac_arr = denfac_fab.array();
1163  auto const& xni_arr = xni_fab.array();
1164  auto const& cpm_arr = cpm_fab.array();
1165  auto const& xl_arr = xl_fab.array();
1166  auto const& qsatw_arr = qsatw_fab.array();
1167  auto const& qsati_arr = qsati_fab.array();
1168  auto const& rhw_arr = rhw_fab.array();
1169  auto const& rhi_arr = rhi_fab.array();
1170  auto const& den_tmp_arr = den_tmp_fab.array();
1171  auto const& delz_tmp_arr = delz_tmp_fab.array();
1172  auto const& n0sfac_arr = n0sfac_fab.array();
1173  auto const& qrs_tmp_r_arr = qrs_tmp_r_fab.array();
1174  auto const& qrs_tmp_s_arr = qrs_tmp_s_fab.array();
1175  auto const& qrs_tmp_g_arr = qrs_tmp_g_fab.array();
1176  auto const& rslope_r_arr = rslope_r_fab.array();
1177  auto const& rslope_s_arr = rslope_s_fab.array();
1178  auto const& rslope_g_arr = rslope_g_fab.array();
1179  auto const& rslope2_r_arr = rslope2_r_fab.array();
1180  auto const& rslope2_s_arr = rslope2_s_fab.array();
1181  auto const& rslope2_g_arr = rslope2_g_fab.array();
1182  auto const& rslope3_r_arr = rslope3_r_fab.array();
1183  auto const& rslope3_s_arr = rslope3_s_fab.array();
1184  auto const& rslope3_g_arr = rslope3_g_fab.array();
1185  auto const& rslopeb_r_arr = rslopeb_r_fab.array();
1186  auto const& rslopeb_s_arr = rslopeb_s_fab.array();
1187  auto const& rslopeb_g_arr = rslopeb_g_fab.array();
1188  auto const& work1_r_arr = work1_r_fab.array();
1189  auto const& work1_s_arr = work1_s_fab.array();
1190  auto const& work1_g_arr = work1_g_fab.array();
1191  auto const& work2_arr = work2_fab.array();
1192  auto const& workdiffw_arr = workdiffw_fab.array();
1193  auto const& workdiffi_arr = workdiffi_fab.array();
1194  auto const& workr_arr = workr_fab.array();
1195  auto const& worka_arr = worka_fab.array();
1196  auto const& work1c_arr = work1c_fab.array();
1197  auto const& denqrs1_arr = denqrs1_fab.array();
1198  auto const& denqrs2_arr = denqrs2_fab.array();
1199  auto const& denqrs3_arr = denqrs3_fab.array();
1200  auto const& denqci_arr = denqci_fab.array();
1201  auto const& fall_r_arr = fall_r_fab.array();
1202  auto const& fall_s_arr = fall_s_fab.array();
1203  auto const& fall_g_arr = fall_g_fab.array();
1204  auto const& fallc_arr = fallc_fab.array();
1205  auto const& qsum_arr = qsum_fab.array();
1206  auto const& nislfv_r_diag_arr = nislfv_r_diag_fab.array();
1207  auto const& nislfv_sg_diag_arr = nislfv_sg_diag_fab.array();
1208  auto const& sed_cell_scratch_arr = sed_cell_scratch_fab.array();
1209  auto const& sed_node_scratch_arr = sed_node_scratch_fab.array();
1210 
1211  ParallelFor(fab_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1212  work1c_arr(i,j,k) = Real(0.0);
1213  });
1214  ParallelFor(fab_box, nislfv_r_diag_fab.nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
1215  nislfv_r_diag_arr(i,j,k,n) = Real(0.0);
1216  nislfv_sg_diag_arr(i,j,k,n) = Real(0.0);
1217  });
1218  auto const& praut_arr = praut_fab.array();
1219  auto const& pracw_arr = pracw_fab.array();
1220  auto const& prevp_arr = prevp_fab.array();
1221  auto const& psdep_arr = psdep_fab.array();
1222  auto const& pgdep_arr = pgdep_fab.array();
1223  auto const& psaut_arr = psaut_fab.array();
1224  auto const& pgaut_arr = pgaut_fab.array();
1225  auto const& praci_arr = praci_fab.array();
1226  auto const& piacr_arr = piacr_fab.array();
1227  auto const& psaci_arr = psaci_fab.array();
1228  auto const& psacw_arr = psacw_fab.array();
1229  auto const& pgacw_arr = pgacw_fab.array();
1230  auto const& pgaci_arr = pgaci_fab.array();
1231  auto const& paacw_arr = paacw_fab.array();
1232  auto const& pracs_arr = pracs_fab.array();
1233  auto const& psacr_arr = psacr_fab.array();
1234  auto const& pgacr_arr = pgacr_fab.array();
1235  auto const& pgacs_arr = pgacs_fab.array();
1236  auto const& pigen_arr = pigen_fab.array();
1237  auto const& pidep_arr = pidep_fab.array();
1238  auto const& pcond_arr = pcond_fab.array();
1239  auto const& psmlt_arr = psmlt_fab.array();
1240  auto const& pgmlt_arr = pgmlt_fab.array();
1241  auto const& pseml_arr = pseml_fab.array();
1242  auto const& pgeml_arr = pgeml_fab.array();
1243  auto const& psevp_arr = psevp_fab.array();
1244  auto const& pgevp_arr = pgevp_fab.array();
1245  auto const& pimlt_arr = pimlt_fab.array();
1246  auto const& pihmf_arr = pihmf_fab.array();
1247  auto const& pihtf_arr = pihtf_fab.array();
1248  auto const& pgfrz_arr = pgfrz_fab.array();
1249 
1250  // Groups A-E: pre-loop setup
1251  // Clamp negative values (Group A)
1252  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1253  qc_arr(i,j,k) = amrex::max(qc_arr(i,j,k), Real(0.0));
1254  qr_arr(i,j,k) = amrex::max(qr_arr(i,j,k), Real(0.0));
1255  qi_arr(i,j,k) = amrex::max(qi_arr(i,j,k), Real(0.0));
1256  qs_arr(i,j,k) = amrex::max(qs_arr(i,j,k), Real(0.0));
1257  qg_arr(i,j,k) = amrex::max(qg_arr(i,j,k), Real(0.0));
1258  den_tmp_arr(i,j,k) = den_arr(i,j,k);
1259  delz_tmp_arr(i,j,k) = delz_arr(i,j,k);
1260  });
1261 
1262  // Group B: cpm, xl — computed once from initial state [lines 455-460]
1263  const Real xlv1_loc = m_xlv1;
1264  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1265  cpm_arr(i,j,k) = wsm6_cpmcal(qv_arr(i,j,k), Real(qmin), Real(cpd), Real(cpv));
1266  xl_arr(i,j,k) = wsm6_xlcal(t_arr(i,j,k), Real(xlv0), xlv1_loc, Real(t0c));
1267  });
1268 
1269  // Outer minor timestep loop (Rule 29)
1270  const int wsm6_loops = std::max(
1271  static_cast<int>(std::round(dt / dtcldcr)), 1);
1272  const Real dtcld = dt / static_cast<Real>(wsm6_loops);
1273  const Real qc0 = m_qc0;
1274  const Real qck1 = m_qck1;
1275  const Real pvtr = m_pvtr;
1276  const Real pacrr = m_pacrr;
1277  const Real precr1 = m_precr1;
1278  const Real precr2 = m_precr2;
1279  const Real roqimax= m_roqimax;
1280  const Real pvts = m_pvts;
1281  const Real pacrc = m_pacrc;
1282  const Real precs1 = m_precs1;
1283  const Real precs2 = m_precs2;
1284  const Real g6pbr = m_g6pbr;
1285  const Real pvtg = m_pvtg;
1286  const Real pacrg = m_pacrg;
1287  const Real precg1 = m_precg1;
1288  const Real precg2 = m_precg2;
1289 
1290  // At the top of a fine patch, rain, snow, ... come in from the coarse level: add
1291  // each minor step's share of it to the top cell before that step, so it falls with
1292  // what is already there, and record it on the top face.
1293  const auto inflow_arr = (m_sed_inflow) ? m_sed_inflow->const_array(mfi)
1294  : Array4<const Real>{};
1295  const Real inflow_share = Real(1) / static_cast<Real>(wsm6_loops);
1296 
1297  for (int loop = 0; loop < wsm6_loops; ++loop) {
1298  if (inflow_arr) {
1299  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
1300  if (inflow_arr(i,j,khi+1,WSM6Sed::NumComps) < Real(0.5)) { return; }
1301  const Real fac = inflow_share / (den_arr(i,j,khi) * delz_arr(i,j,khi));
1302  qi_arr(i,j,khi) += inflow_arr(i,j,khi+1,WSM6Sed::qi) * fac;
1303  qr_arr(i,j,khi) += inflow_arr(i,j,khi+1,WSM6Sed::qr) * fac;
1304  qs_arr(i,j,khi) += inflow_arr(i,j,khi+1,WSM6Sed::qs) * fac;
1305  qg_arr(i,j,khi) += inflow_arr(i,j,khi+1,WSM6Sed::qg) * fac;
1306  for (int n = 0; n < WSM6Sed::NumComps; ++n) {
1307  sed_arr(i,j,khi+1,n) += inflow_arr(i,j,khi+1,n) * inflow_share;
1308  }
1309  });
1310  }
1311  // G1b: denfac = sqrt(den0/den) [lines 503-515]
1312  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1313  const Real invden = Real(1.0) / den_arr(i,j,k);
1314  denfac_arr(i,j,k) = std::sqrt(invden * Real(den0));
1315  });
1316  // G1c: qsatw, qsati, rhw, rhi [lines 517-549]
1317  {
1318  const Real ttp = Real(t0c) + Real(0.01);
1319  const Real dldt = Real(cpv) - Real(cliq);
1320  const Real xa = -dldt / Real(rv);
1321  const Real xb = xa + Real(xlv0) / (Real(rv)*ttp);
1322  const Real dldti= Real(cpv) - Real(cice);
1323  const Real xai = -dldti / Real(rv);
1324  const Real xbi = xai + Real(xls) / (Real(rv)*ttp);
1325  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1326  const Real tr = ttp / t_arr(i,j,k);
1327  Real qsw = Real(psat)*std::exp(std::log(tr)*xa)*std::exp(xb*(Real(1.0)-tr));
1328  qsw = amrex::min(qsw, Real(0.99)*p_arr(i,j,k));
1329  qsw = Real(ep2)*qsw / (p_arr(i,j,k) - qsw);
1330  qsw = amrex::max(qsw, Real(qmin));
1331  qsatw_arr(i,j,k) = qsw;
1332  rhw_arr(i,j,k) = amrex::max(qv_arr(i,j,k)/qsw, Real(qmin));
1333  Real qsi = (t_arr(i,j,k) < ttp)
1334  ? Real(psat)*std::exp(std::log(tr)*xai)*std::exp(xbi*(Real(1.0)-tr))
1335  : Real(psat)*std::exp(std::log(tr)*xa )*std::exp(xb *(Real(1.0)-tr));
1336  qsi = amrex::min(qsi, Real(0.99)*p_arr(i,j,k));
1337  qsi = Real(ep2)*qsi / (p_arr(i,j,k) - qsi);
1338  qsi = amrex::max(qsi, Real(qmin));
1339  qsati_arr(i,j,k) = qsi;
1340  rhi_arr(i,j,k) = amrex::max(qv_arr(i,j,k)/qsi, Real(qmin));
1341  });
1342  }
1343  // G2: zero all process rates each sub-step [lines 555-594]
1344  // WSM6-CPP TAG: RATES_ZERO
1345  // legacy_group: G2
1346  // process: Zero all process rates each sub-step
1347  // compare_vars: prevp, psdep, pgdep, praut, psaut, pgaut, pracw, praci, psaci, pracs, pidep, pcond, psmlt, pgmlt, pseml, psevp, fall_r, fall_s, fall_g, fallc
1348  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1349  prevp_arr(i,j,k) = Real(0.0); psdep_arr(i,j,k) = Real(0.0);
1350  pgdep_arr(i,j,k) = Real(0.0); praut_arr(i,j,k) = Real(0.0);
1351  psaut_arr(i,j,k) = Real(0.0); pgaut_arr(i,j,k) = Real(0.0);
1352  pracw_arr(i,j,k) = Real(0.0); praci_arr(i,j,k) = Real(0.0);
1353  piacr_arr(i,j,k) = Real(0.0); psaci_arr(i,j,k) = Real(0.0);
1354  psacw_arr(i,j,k) = Real(0.0); pracs_arr(i,j,k) = Real(0.0);
1355  psacr_arr(i,j,k) = Real(0.0); pgacw_arr(i,j,k) = Real(0.0);
1356  paacw_arr(i,j,k) = Real(0.0); pgaci_arr(i,j,k) = Real(0.0);
1357  pgacr_arr(i,j,k) = Real(0.0); pgacs_arr(i,j,k) = Real(0.0);
1358  pigen_arr(i,j,k) = Real(0.0); pidep_arr(i,j,k) = Real(0.0);
1359  pcond_arr(i,j,k) = Real(0.0); psmlt_arr(i,j,k) = Real(0.0);
1360  pgmlt_arr(i,j,k) = Real(0.0); pseml_arr(i,j,k) = Real(0.0);
1361  pgeml_arr(i,j,k) = Real(0.0); psevp_arr(i,j,k) = Real(0.0);
1362  pgevp_arr(i,j,k) = Real(0.0);
1363  fall_r_arr(i,j,k) = Real(0.0); fall_s_arr(i,j,k) = Real(0.0);
1364  fall_g_arr(i,j,k) = Real(0.0); fallc_arr(i,j,k) = Real(0.0);
1365  });
1366 
1367  // G3: xni ice crystal number concentration [lines 598-604]
1368  // WSM6-CPP TAG: XNI
1369  // legacy_group: G3
1370  // process: Ice crystal number concentration
1371  // compare_vars: xni, qi, den
1372  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1373  const Real tmp = den_arr(i,j,k)*amrex::max(qi_arr(i,j,k), Real(qmin));
1374  xni_arr(i,j,k) = amrex::min(
1375  amrex::max(Real(5.38e7)*std::sqrt(std::sqrt(tmp*tmp*tmp)), Real(1.e3)),
1376  Real(1.e6));
1377  });
1378 
1379  // G4: pack qrs_tmp, first slope_wsm6 [lines 610-618]
1380  // WSM6-CPP TAG: SLOPE1
1381  // legacy_group: G4
1382  // process: First slope calculation
1383  // compare_vars: rslope, rslope2, rslope3, rslopeb, falk, fall, work1
1384  const Real pidn0r_loc = m_pidn0r;
1385  const Real rslopermax_loc = m_rslopermax;
1386  const Real rsloperbmax_loc = m_rsloperbmax;
1387  const Real rsloper2max_loc = m_rsloper2max;
1388  const Real rsloper3max_loc = m_rsloper3max;
1389  const Real pidn0s_loc = m_pidn0s;
1390  const Real rslopesmax_loc = m_rslopesmax;
1391  const Real rslopesbmax_loc = m_rslopesbmax;
1392  const Real rslopes2max_loc = m_rslopes2max;
1393  const Real rslopes3max_loc = m_rslopes3max;
1394  const Real pidn0g_loc = m_pidn0g;
1395  const Real rslopegmax_loc = m_rslopegmax;
1396  const Real rslopegbmax_loc = m_rslopegbmax;
1397  const Real rslopeg2max_loc = m_rslopeg2max;
1398  const Real rslopeg3max_loc = m_rslopeg3max;
1399  const Real bvtg_loc = m_bvtg;
1400  const Real pi_wsm6_loc = m_pi_wsm6;
1401  const Real n0g_loc = m_n0g;
1402 
1403  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1404  qrs_tmp_r_arr(i,j,k) = qr_arr(i,j,k);
1405  qrs_tmp_s_arr(i,j,k) = qs_arr(i,j,k);
1406  qrs_tmp_g_arr(i,j,k) = qg_arr(i,j,k);
1407  Real dummy_n0sfac;
1409  qrs_tmp_r_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1410  pidn0r_loc, Real(qcrmin), rslopermax_loc, rsloperbmax_loc,
1411  rsloper2max_loc, rsloper3max_loc, Real(bvtr), Real(pvtr),
1412  rslope_r_arr(i,j,k), rslopeb_r_arr(i,j,k),
1413  rslope2_r_arr(i,j,k), rslope3_r_arr(i,j,k),
1414  work1_r_arr(i,j,k));
1416  qrs_tmp_s_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1417  t_arr(i,j,k), pidn0s_loc, Real(alpha_wsm6),
1418  Real(n0smax), Real(n0s), Real(t0c), Real(qcrmin),
1419  rslopesmax_loc, rslopesbmax_loc,
1420  rslopes2max_loc, rslopes3max_loc,
1421  Real(bvts), Real(pvts),
1422  rslope_s_arr(i,j,k), rslopeb_s_arr(i,j,k),
1423  rslope2_s_arr(i,j,k), rslope3_s_arr(i,j,k),
1424  work1_s_arr(i,j,k), dummy_n0sfac);
1426  qrs_tmp_g_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1427  pidn0g_loc, Real(qcrmin),
1428  rslopegmax_loc, rslopegbmax_loc,
1429  rslopeg2max_loc, rslopeg3max_loc,
1430  bvtg_loc, Real(pvtg),
1431  rslope_g_arr(i,j,k), rslopeb_g_arr(i,j,k),
1432  rslope2_g_arr(i,j,k), rslope3_g_arr(i,j,k),
1433  work1_g_arr(i,j,k));
1434  n0sfac_arr(i,j,k) = dummy_n0sfac;
1435  });
1436  // G5a-G5e: sedimentation setup, nislfv calls, and flux updates
1437  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
1438  const int km_local = khi - klo + 1;
1439 
1440  constexpr Real qsum_min = Real(1.0e-15);
1441  auto den_col = [&](int k) -> amrex::Real& {
1442  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::den);
1443  };
1444  auto denfac_col = [&](int k) -> amrex::Real& {
1445  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denfac);
1446  };
1447  auto t_col = [&](int k) -> amrex::Real& {
1448  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::tk);
1449  };
1450  auto dz_col = [&](int k) -> amrex::Real& {
1451  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::dz);
1452  };
1453  auto workr_col = [&](int k) -> amrex::Real& {
1454  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::workr_col);
1455  };
1456  auto worka_col = [&](int k) -> amrex::Real& {
1457  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::worka_col);
1458  };
1459  auto denqrs1_col = [&](int k) -> amrex::Real& {
1460  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denqrs1_col);
1461  };
1462  auto denqrs2_col = [&](int k) -> amrex::Real& {
1463  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denqrs2_col);
1464  };
1465  auto denqrs3_col = [&](int k) -> amrex::Real& {
1466  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denqrs3_col);
1467  };
1468  auto qsum_col = [&](int k) -> amrex::Real& {
1469  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::qsum_col);
1470  };
1471  Real delqrs1_col = Real(0.0);
1472  Real delqrs2_col = Real(0.0);
1473  Real delqrs3_col = Real(0.0);
1474 
1475  // G5a: pack sedimentation work arrays
1476  for (int k = klo; k <= khi; ++k) {
1477  const int kk = k - klo;
1478  den_col(kk) = den_arr(i,j,k);
1479  denfac_col(kk) = denfac_arr(i,j,k);
1480  t_col(kk) = t_arr(i,j,k);
1481  dz_col(kk) = delz_tmp_arr(i,j,k);
1482  workr_col(kk) = work1_r_arr(i,j,k);
1483  qsum_col(kk) = amrex::max(qs_arr(i,j,k) + qg_arr(i,j,k), qsum_min);
1484  if (qsum_col(kk) > qsum_min) {
1485  worka_col(kk) = (work1_s_arr(i,j,k) * qs_arr(i,j,k)
1486  + work1_g_arr(i,j,k) * qg_arr(i,j,k))
1487  / qsum_col(kk);
1488  } else {
1489  worka_col(kk) = Real(0.0);
1490  }
1491  denqrs1_col(kk) = den_col(kk) * qr_arr(i,j,k);
1492  denqrs2_col(kk) = den_col(kk) * qs_arr(i,j,k);
1493  denqrs3_col(kk) = den_col(kk) * qg_arr(i,j,k);
1494  if (qr_arr(i,j,k) <= Real(0.0)) {
1495  workr_col(kk) = Real(0.0);
1496  }
1497  }
1498 
1499  // The face fluxes are what the cells above each face lose: add the
1500  // top-down sums of rho*q*dz now, and subtract them after G5d.
1501  sedimentation_add_column_sums(sed_arr, WSM6Sed::qr, i, j, klo, khi, Real(1.0),
1502  [&](int k) { return den_arr(i,j,k) * qr_arr(i,j,k) * delz_arr(i,j,k); });
1503  sedimentation_add_column_sums(sed_arr, WSM6Sed::qs, i, j, klo, khi, Real(1.0),
1504  [&](int k) { return den_arr(i,j,k) * qs_arr(i,j,k) * delz_arr(i,j,k); });
1505  sedimentation_add_column_sums(sed_arr, WSM6Sed::qg, i, j, klo, khi, Real(1.0),
1506  [&](int k) { return den_arr(i,j,k) * qg_arr(i,j,k) * delz_arr(i,j,k); });
1507 
1508  // G5b: rain sedimentation
1510  km_local,
1513  &delqrs1_col, dtcld, 1,
1514  sed_cell_scratch_arr, sed_node_scratch_arr, i, j, klo);
1515  // Strict Rule 30 snapshot: immediately after G5b
1516  for (int k = klo; k <= khi; ++k) {
1517  const int kk = k - klo;
1518  nislfv_r_diag_arr(i,j,k,0) = amrex::max(denqrs1_col(kk) / den_col(kk), Real(0.0));
1519  nislfv_r_diag_arr(i,j,k,1) = denqrs1_col(kk) * workr_col(kk) / delz_arr(i,j,k);
1520  nislfv_r_diag_arr(i,j,k,2) = workr_col(kk);
1521  nislfv_r_diag_arr(i,j,k,3) = denqrs1_col(kk);
1522  nislfv_r_diag_arr(i,j,k,4) = den_col(kk);
1523  nislfv_r_diag_arr(i,j,k,5) = denfac_col(kk);
1524  }
1525 
1526  // G5c: snow + graupel sedimentation
1528  km_local,
1532  &delqrs2_col, &delqrs3_col, dtcld, 1,
1533  sed_cell_scratch_arr, sed_node_scratch_arr, i, j, klo);
1534  // Strict Rule 30 snapshot: immediately after G5c
1535  for (int k = klo; k <= khi; ++k) {
1536  const int kk = k - klo;
1537  nislfv_sg_diag_arr(i,j,k,0) = amrex::max(denqrs2_col(kk) / den_col(kk), Real(0.0));
1538  nislfv_sg_diag_arr(i,j,k,1) = amrex::max(denqrs3_col(kk) / den_col(kk), Real(0.0));
1539  nislfv_sg_diag_arr(i,j,k,2) = denqrs2_col(kk) * worka_col(kk) / delz_arr(i,j,k);
1540  nislfv_sg_diag_arr(i,j,k,3) = denqrs3_col(kk) * worka_col(kk) / delz_arr(i,j,k);
1541  nislfv_sg_diag_arr(i,j,k,4) = denqrs2_col(kk);
1542  nislfv_sg_diag_arr(i,j,k,5) = denqrs3_col(kk);
1543  }
1544 
1545  // G5d: update species and fall speeds
1546  for (int k = klo; k <= khi; ++k) {
1547  const int kk = k - klo;
1548  qsum_arr(i,j,k) = qsum_col(kk);
1549  workr_arr(i,j,k) = workr_col(kk);
1550  worka_arr(i,j,k) = worka_col(kk);
1551  denqrs1_arr(i,j,k) = denqrs1_col(kk);
1552  denqrs2_arr(i,j,k) = denqrs2_col(kk);
1553  denqrs3_arr(i,j,k) = denqrs3_col(kk);
1554  qr_arr(i,j,k) = amrex::max(denqrs1_col(kk) / den_col(kk), Real(0.0));
1555  qs_arr(i,j,k) = amrex::max(denqrs2_col(kk) / den_col(kk), Real(0.0));
1556  qg_arr(i,j,k) = amrex::max(denqrs3_col(kk) / den_col(kk), Real(0.0));
1557  fall_r_arr(i,j,k) = denqrs1_col(kk) * workr_col(kk) / delz_arr(i,j,k);
1558  fall_s_arr(i,j,k) = denqrs2_col(kk) * worka_col(kk) / delz_arr(i,j,k);
1559  fall_g_arr(i,j,k) = denqrs3_col(kk) * worka_col(kk) / delz_arr(i,j,k);
1560  }
1561 
1562  sedimentation_add_column_sums(sed_arr, WSM6Sed::qr, i, j, klo, khi, Real(-1.0),
1563  [&](int k) { return den_arr(i,j,k) * qr_arr(i,j,k) * delz_arr(i,j,k); });
1564  sedimentation_add_column_sums(sed_arr, WSM6Sed::qs, i, j, klo, khi, Real(-1.0),
1565  [&](int k) { return den_arr(i,j,k) * qs_arr(i,j,k) * delz_arr(i,j,k); });
1566  sedimentation_add_column_sums(sed_arr, WSM6Sed::qg, i, j, klo, khi, Real(-1.0),
1567  [&](int k) { return den_arr(i,j,k) * qg_arr(i,j,k) * delz_arr(i,j,k); });
1568 
1569  // G5e: slab fall fluxes at the lower boundary
1570  delqrs1_arr(i,j,0) = delqrs1_col / delz_arr(i,j,klo) / dtcld;
1571  delqrs2_arr(i,j,0) = delqrs2_col / delz_arr(i,j,klo) / dtcld;
1572  delqrs3_arr(i,j,0) = delqrs3_col / delz_arr(i,j,klo) / dtcld;
1573  fall_r_arr(i,j,klo) = delqrs1_arr(i,j,0);
1574  fall_s_arr(i,j,klo) = delqrs2_arr(i,j,0);
1575  fall_g_arr(i,j,klo) = delqrs3_arr(i,j,0);
1576  });
1577 
1578  // G6: repack qrs_tmp, second slope_wsm6 [lines 655-663]
1579  // slope params updated after sedimentation moved mass
1580  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1581  qrs_tmp_r_arr(i,j,k) = qr_arr(i,j,k);
1582  qrs_tmp_s_arr(i,j,k) = qs_arr(i,j,k);
1583  qrs_tmp_g_arr(i,j,k) = qg_arr(i,j,k);
1584  Real dummy_n0sfac;
1586  qrs_tmp_r_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1587  pidn0r_loc, Real(qcrmin), rslopermax_loc, rsloperbmax_loc,
1588  rsloper2max_loc, rsloper3max_loc, Real(bvtr), Real(pvtr),
1589  rslope_r_arr(i,j,k), rslopeb_r_arr(i,j,k),
1590  rslope2_r_arr(i,j,k), rslope3_r_arr(i,j,k),
1591  work1_r_arr(i,j,k));
1593  qrs_tmp_s_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1594  t_arr(i,j,k), pidn0s_loc, Real(alpha_wsm6),
1595  Real(n0smax), Real(n0s), Real(t0c), Real(qcrmin),
1596  rslopesmax_loc, rslopesbmax_loc,
1597  rslopes2max_loc, rslopes3max_loc,
1598  Real(bvts), Real(pvts),
1599  rslope_s_arr(i,j,k), rslopeb_s_arr(i,j,k),
1600  rslope2_s_arr(i,j,k), rslope3_s_arr(i,j,k),
1601  work1_s_arr(i,j,k), dummy_n0sfac);
1603  qrs_tmp_g_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1604  pidn0g_loc, Real(qcrmin),
1605  rslopegmax_loc, rslopegbmax_loc,
1606  rslopeg2max_loc, rslopeg3max_loc,
1607  bvtg_loc, Real(pvtg),
1608  rslope_g_arr(i,j,k), rslopeb_g_arr(i,j,k),
1609  rslope2_g_arr(i,j,k), rslope3_g_arr(i,j,k),
1610  work1_g_arr(i,j,k));
1611  });
1612 
1613  // G7: melting (T>T0 only) [lines 665-704]
1614  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1615  const Real supcol = Real(t0c) - t_arr(i,j,k);
1616  n0sfac_arr(i,j,k) = amrex::max(
1617  amrex::min(std::exp(Real(alpha_wsm6) * supcol),
1618  Real(n0smax) / Real(n0s)),
1619  Real(1.0));
1620 
1621  if (t_arr(i,j,k) > Real(t0c)) {
1622  const Real xlf = Real(xlf0);
1623  work2_arr(i,j,k) = wsm6_venfac(
1624  p_arr(i,j,k), t_arr(i,j,k), den_arr(i,j,k), Real(den0));
1625 
1626  if (qs_arr(i,j,k) > Real(0.0)) {
1627  const Real coeres =
1628  rslope2_s_arr(i,j,k) *
1629  std::sqrt(rslope_s_arr(i,j,k) * rslopeb_s_arr(i,j,k));
1630  psmlt_arr(i,j,k) =
1631  wsm6_xka(t_arr(i,j,k), den_arr(i,j,k)) / xlf *
1632  (Real(t0c) - t_arr(i,j,k)) * pi_wsm6_loc * Real(0.5) *
1633  n0sfac_arr(i,j,k) *
1634  (Real(precs1) * rslope2_s_arr(i,j,k) +
1635  Real(precs2) * work2_arr(i,j,k) * coeres) /
1636  den_arr(i,j,k);
1637  psmlt_arr(i,j,k) = amrex::min(
1638  amrex::max(psmlt_arr(i,j,k) * dtcld,
1639  -qs_arr(i,j,k)),
1640  Real(0.0));
1641  qs_arr(i,j,k) = qs_arr(i,j,k) + psmlt_arr(i,j,k);
1642  qr_arr(i,j,k) = qr_arr(i,j,k) - psmlt_arr(i,j,k);
1643  t_arr(i,j,k) = t_arr(i,j,k) + xlf / cpm_arr(i,j,k) * psmlt_arr(i,j,k);
1644  }
1645 
1646  if (qg_arr(i,j,k) > Real(0.0)) {
1647  const Real coeres =
1648  rslope2_g_arr(i,j,k) *
1649  std::sqrt(rslope_g_arr(i,j,k) * rslopeb_g_arr(i,j,k));
1650  pgmlt_arr(i,j,k) =
1651  wsm6_xka(t_arr(i,j,k), den_arr(i,j,k)) / xlf *
1652  (Real(t0c) - t_arr(i,j,k)) *
1653  (Real(precg1) * rslope2_g_arr(i,j,k) +
1654  Real(precg2) * work2_arr(i,j,k) * coeres) /
1655  den_arr(i,j,k);
1656  pgmlt_arr(i,j,k) = amrex::min(
1657  amrex::max(pgmlt_arr(i,j,k) * dtcld,
1658  -qg_arr(i,j,k)),
1659  Real(0.0));
1660  qg_arr(i,j,k) = qg_arr(i,j,k) + pgmlt_arr(i,j,k);
1661  qr_arr(i,j,k) = qr_arr(i,j,k) - pgmlt_arr(i,j,k);
1662  t_arr(i,j,k) = t_arr(i,j,k) + xlf / cpm_arr(i,j,k) * pgmlt_arr(i,j,k);
1663  }
1664  }
1665  });
1666 
1667  // G8: cloud ice sedimentation/fallout [lines 708-735]
1668  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1669  if (qi_arr(i,j,k) <= Real(0.0)) {
1670  work1c_arr(i,j,k) = Real(0.0);
1671  } else {
1672  const Real tmp = den_arr(i,j,k) *
1673  amrex::max(qi_arr(i,j,k), Real(qmin));
1674  const Real xni = amrex::min(
1675  amrex::max(Real(5.38e7) *
1676  std::sqrt(std::sqrt(tmp*tmp*tmp)),
1677  Real(1.0e3)),
1678  Real(1.0e6));
1679  const Real xmi = den_arr(i,j,k) * qi_arr(i,j,k) / xni;
1680  const Real diameter = amrex::max(
1681  amrex::min(Real(dicon) * std::sqrt(xmi), Real(dimax)),
1682  Real(1.0e-25));
1683  work1c_arr(i,j,k) = Real(1.49e4) *
1684  std::exp(std::log(diameter) * Real(1.31));
1685  }
1686  denqci_arr(i,j,k) = den_arr(i,j,k) * qi_arr(i,j,k);
1687  });
1688 
1689  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
1690  const int km_local = khi - klo + 1;
1691  auto den_col = [&](int k) -> amrex::Real& {
1692  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::den);
1693  };
1694  auto denfac_col = [&](int k) -> amrex::Real& {
1695  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denfac);
1696  };
1697  auto t_col = [&](int k) -> amrex::Real& {
1698  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::tk);
1699  };
1700  auto dz_col = [&](int k) -> amrex::Real& {
1701  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::dz);
1702  };
1703  auto work1c_col = [&](int k) -> amrex::Real& {
1704  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::work1c_col);
1705  };
1706  auto denqci_col = [&](int k) -> amrex::Real& {
1707  return sed_cell_scratch_arr(i, j, klo + k, WSM6SedCellScratch::denqci_col);
1708  };
1709  Real delqi_col = Real(0.0);
1710 
1711  for (int k = klo; k <= khi; ++k) {
1712  const int kk = k - klo;
1713  den_col(kk) = den_arr(i,j,k);
1714  denfac_col(kk) = denfac_arr(i,j,k);
1715  t_col(kk) = t_arr(i,j,k);
1716  dz_col(kk) = delz_tmp_arr(i,j,k);
1717  work1c_col(kk) = work1c_arr(i,j,k);
1718  denqci_col(kk) = denqci_arr(i,j,k);
1719  }
1720 
1721  // Face fluxes, as for G5: the top-down sums before, minus those after
1722  sedimentation_add_column_sums(sed_arr, WSM6Sed::qi, i, j, klo, khi, Real(1.0),
1723  [&](int k) { return den_arr(i,j,k) * qi_arr(i,j,k) * delz_arr(i,j,k); });
1724 
1726  km_local,
1729  &delqi_col, dtcld, 0,
1730  sed_cell_scratch_arr, sed_node_scratch_arr, i, j, klo);
1731 
1732  for (int k = klo; k <= khi; ++k) {
1733  const int kk = k - klo;
1734  work1c_arr(i,j,k) = work1c_col(kk);
1735  denqci_arr(i,j,k) = denqci_col(kk);
1736  qi_arr(i,j,k) = amrex::max(
1737  denqci_col(kk) / den_col(kk), Real(0.0));
1738  }
1739  sedimentation_add_column_sums(sed_arr, WSM6Sed::qi, i, j, klo, khi, Real(-1.0),
1740  [&](int k) { return den_arr(i,j,k) * qi_arr(i,j,k) * delz_arr(i,j,k); });
1741 
1742  delqi_arr(i,j,0) = delqi_col / delz_arr(i,j,klo) / dtcld;
1743  fallc_arr(i,j,klo) = delqi_arr(i,j,0);
1744  });
1745 
1746  // G9: surface precipitation accumulation [lines 741-770]
1747  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
1748  const Real fallsum =
1749  fall_r_arr(i,j,klo) + fall_s_arr(i,j,klo) +
1750  fall_g_arr(i,j,klo) + fallc_arr(i,j,klo);
1751  const Real fallsum_qsi = fall_s_arr(i,j,klo) + fallc_arr(i,j,klo);
1752  const Real fallsum_qg = fall_g_arr(i,j,klo);
1753  const Real precip = delz_arr(i,j,klo) / denr * dtcld * Real(1000.0);
1754 
1755  if (fallsum > Real(0.0)) {
1756  rainncv_arr(i,j,0) += fallsum * precip;
1757  rainacc_arr(i,j,0) += fallsum * precip;
1758  }
1759  if (fallsum_qsi > Real(0.0)) {
1760  tstepsnow_arr(i,j,0) += fallsum_qsi * precip;
1761  snowncv_arr(i,j,0) += fallsum_qsi * precip;
1762  snowacc_arr(i,j,0) += fallsum_qsi * precip;
1763  }
1764  if (fallsum_qg > Real(0.0)) {
1765  tstepgraup_arr(i,j,0) += fallsum_qg * precip;
1766  graupelncv_arr(i,j,0) += fallsum_qg * precip;
1767  graupacc_arr(i,j,0) += fallsum_qg * precip;
1768  }
1769  rainonly_arr(i,j,0) += fall_r_arr(i,j,klo) * precip;
1770  snowonly_arr(i,j,0) += fall_s_arr(i,j,klo) * precip;
1771  if (fallsum > Real(0.0)) {
1772  sr_arr(i,j,0) =
1773  (snowncv_arr(i,j,0) + graupelncv_arr(i,j,0)) /
1774  (rainncv_arr(i,j,0) + Real(1.0e-12));
1775  }
1776  });
1777 
1778  constexpr Real pi = Real(3.141592653589793238462643383279502884);
1779 
1780  // G10: instantaneous phase changes [lines 774-830]
1781  // WSM6-CPP TAG: PHASE
1782  // legacy_group: G10
1783  // process: Instantaneous phase changes
1784  // compare_vars: t, q, qci, qrs
1785  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1786  const Real supcol = Real(t0c) - t_arr(i,j,k);
1787  const Real xlf = (supcol < Real(0.0))
1788  ? Real(xlf0)
1789  : Real(xls) - xl_arr(i,j,k);
1790  pimlt_arr(i,j,k) = Real(0.0);
1791  pihmf_arr(i,j,k) = Real(0.0);
1792  pihtf_arr(i,j,k) = Real(0.0);
1793  pgfrz_arr(i,j,k) = Real(0.0);
1794 
1795  if (supcol < Real(0.0) && qi_arr(i,j,k) > Real(0.0)) {
1796  pimlt_arr(i,j,k) = qi_arr(i,j,k);
1797  qc_arr(i,j,k) = qc_arr(i,j,k) + qi_arr(i,j,k);
1798  t_arr(i,j,k) = t_arr(i,j,k) - xlf / cpm_arr(i,j,k) * qi_arr(i,j,k);
1799  qi_arr(i,j,k) = Real(0.0);
1800  }
1801 
1802  if (supcol > Real(40.0) && qc_arr(i,j,k) > Real(0.0)) {
1803  pihmf_arr(i,j,k) = qc_arr(i,j,k);
1804  qi_arr(i,j,k) = qi_arr(i,j,k) + qc_arr(i,j,k);
1805  t_arr(i,j,k) = t_arr(i,j,k) + xlf / cpm_arr(i,j,k) * qc_arr(i,j,k);
1806  qc_arr(i,j,k) = Real(0.0);
1807  }
1808 
1809  if (supcol > Real(0.0) && qc_arr(i,j,k) > Real(qmin)) {
1810  const Real supcolt = amrex::min(supcol, Real(50.0));
1811  const Real pfrzdtc = amrex::min(
1812  Real(pfrz1) *
1813  (std::exp(Real(pfrz2) * supcolt) - Real(1.0)) *
1814  den_arr(i,j,k) / Real(denr) / Real(WSM6::xncr) *
1815  qc_arr(i,j,k) * qc_arr(i,j,k) * dtcld,
1816  qc_arr(i,j,k));
1817  pihtf_arr(i,j,k) = pfrzdtc;
1818  qi_arr(i,j,k) = qi_arr(i,j,k) + pfrzdtc;
1819  t_arr(i,j,k) = t_arr(i,j,k) + xlf / cpm_arr(i,j,k) * pfrzdtc;
1820  qc_arr(i,j,k) = qc_arr(i,j,k) - pfrzdtc;
1821  }
1822 
1823  if (supcol > Real(0.0) && qr_arr(i,j,k) > Real(0.0)) {
1824  Real temp = rslope3_r_arr(i,j,k);
1825  temp = temp * temp * rslope_r_arr(i,j,k);
1826  const Real supcolt = amrex::min(supcol, Real(50.0));
1827  const Real pfrzdtr = amrex::min(
1828  Real(20.0) * pi * pi * Real(pfrz1) * Real(WSM6::n0r) *
1829  Real(denr) / den_arr(i,j,k) *
1830  (std::exp(Real(pfrz2) * supcolt) - Real(1.0)) *
1831  temp * dtcld,
1832  qr_arr(i,j,k));
1833  pgfrz_arr(i,j,k) = pfrzdtr;
1834  qg_arr(i,j,k) = qg_arr(i,j,k) + pfrzdtr;
1835  t_arr(i,j,k) = t_arr(i,j,k) + xlf / cpm_arr(i,j,k) * pfrzdtr;
1836  qr_arr(i,j,k) = qr_arr(i,j,k) - pfrzdtr;
1837  }
1838  });
1839  // G11: third slope_wsm6 call [lines 836-844]
1840  // WSM6-CPP TAG: SLOPE3
1841  // legacy_group: G11
1842  // process: Third slope calculation
1843  // compare_vars: rslope, rslope2, rslope3, rslopeb, falk, fall, work1
1844  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1845  qrs_tmp_r_arr(i,j,k) = qr_arr(i,j,k);
1846  qrs_tmp_s_arr(i,j,k) = qs_arr(i,j,k);
1847  qrs_tmp_g_arr(i,j,k) = qg_arr(i,j,k);
1848  Real dummy_n0sfac;
1850  qrs_tmp_r_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1851  pidn0r_loc, Real(qcrmin), rslopermax_loc, rsloperbmax_loc,
1852  rsloper2max_loc, rsloper3max_loc, Real(bvtr), Real(pvtr),
1853  rslope_r_arr(i,j,k), rslopeb_r_arr(i,j,k),
1854  rslope2_r_arr(i,j,k), rslope3_r_arr(i,j,k),
1855  work1_r_arr(i,j,k));
1857  qrs_tmp_s_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1858  t_arr(i,j,k), pidn0s_loc, Real(alpha_wsm6),
1859  Real(n0smax), Real(n0s), Real(t0c), Real(qcrmin),
1860  rslopesmax_loc, rslopesbmax_loc,
1861  rslopes2max_loc, rslopes3max_loc,
1862  Real(bvts), Real(pvts),
1863  rslope_s_arr(i,j,k), rslopeb_s_arr(i,j,k),
1864  rslope2_s_arr(i,j,k), rslope3_s_arr(i,j,k),
1865  work1_s_arr(i,j,k), dummy_n0sfac);
1867  qrs_tmp_g_arr(i,j,k), den_arr(i,j,k), denfac_arr(i,j,k),
1868  pidn0g_loc, Real(qcrmin),
1869  rslopegmax_loc, rslopegbmax_loc,
1870  rslopeg2max_loc, rslopeg3max_loc,
1871  bvtg_loc, Real(pvtg),
1872  rslope_g_arr(i,j,k), rslopeb_g_arr(i,j,k),
1873  rslope2_g_arr(i,j,k), rslope3_g_arr(i,j,k),
1874  work1_g_arr(i,j,k));
1875  n0sfac_arr(i,j,k) = dummy_n0sfac;
1876  });
1877  // G12: workdiffw, workdiffi, work2 [lines 851-857]
1878  // WSM6-CPP TAG: DIFF_PREP
1879  // legacy_group: G12
1880  // process: Prepare diffusion/work terms
1881  // compare_vars: workdiffw, workdiffi, work2
1882  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1883  workdiffw_arr(i,j,k) = wsm6_diffac(
1884  xl_arr(i,j,k), p_arr(i,j,k), t_arr(i,j,k),
1885  den_arr(i,j,k), qsatw_arr(i,j,k), Real(rv));
1886  workdiffi_arr(i,j,k) = wsm6_diffac(
1887  Real(xls), p_arr(i,j,k), t_arr(i,j,k),
1888  den_arr(i,j,k), qsati_arr(i,j,k), Real(rv));
1889  work2_arr(i,j,k) = wsm6_venfac(
1890  p_arr(i,j,k), t_arr(i,j,k), den_arr(i,j,k), Real(den0));
1891  });
1892 
1893  // G13a: warm rain — praut, pracw, prevp [lines 867-903]
1894  {
1895  const Real dtcld_l = dtcld;
1896  const Real qmin_l = Real(qmin);
1897  const Real qcrmin_l= Real(qcrmin);
1898  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1899  Real supsat = amrex::max(qv_arr(i,j,k), qmin_l)
1900  - qsatw_arr(i,j,k);
1901  Real satdt = supsat / dtcld_l;
1902  // praut: C->R autoconversion
1903  if (qc_arr(i,j,k) > Real(qc0)) {
1904  praut_arr(i,j,k) = amrex::min(
1905  Real(qck1)*std::pow(qc_arr(i,j,k), Real(7.0/3.0)),
1906  qc_arr(i,j,k)/dtcld_l);
1907  }
1908  // pracw: C->R accretion by rain
1909  if (qr_arr(i,j,k) > qcrmin_l && qc_arr(i,j,k) > qmin_l) {
1910  pracw_arr(i,j,k) = amrex::min(
1911  Real(pacrr)*rslope3_r_arr(i,j,k)
1912  *rslopeb_r_arr(i,j,k)
1913  *qc_arr(i,j,k)*denfac_arr(i,j,k),
1914  qc_arr(i,j,k)/dtcld_l);
1915  }
1916  // prevp: R evaporation/condensation
1917  if (qr_arr(i,j,k) > Real(0.0)) {
1918  Real coeres = rslope2_r_arr(i,j,k)*std::sqrt(
1919  rslope_r_arr(i,j,k)*rslopeb_r_arr(i,j,k));
1920  Real rate = (rhw_arr(i,j,k)-Real(1.0))
1921  *(Real(precr1)*rslope2_r_arr(i,j,k)
1922  + Real(precr2)*work2_arr(i,j,k)*coeres)
1923  / workdiffw_arr(i,j,k);
1924  if (rate < Real(0.0)) {
1925  rate = amrex::max(rate, -qr_arr(i,j,k)/dtcld_l);
1926  rate = amrex::max(rate, satdt/Real(2.0));
1927  } else {
1928  rate = amrex::min(rate, satdt/Real(2.0));
1929  }
1930  prevp_arr(i,j,k) = rate;
1931  }
1932  });
1933  }
1934  // G13b: cold-rain, mixed-phase, and ice deposition/nucleation
1935  // [lines 904-1062, 1075-1192]
1936  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
1937  const Real supcol = Real(t0c) - t_arr(i,j,k);
1938  n0sfac_arr(i,j,k) = amrex::max(
1939  amrex::min(std::exp(Real(alpha_wsm6) * supcol),
1940  Real(n0smax) / Real(n0s)),
1941  Real(1.0));
1942 
1943  const Real supsat = amrex::max(qv_arr(i,j,k), Real(qmin))
1944  - qsati_arr(i,j,k);
1945  const Real satdt = supsat / dtcld;
1946  int ifsat = 0;
1947 
1948  const Real tmp = den_arr(i,j,k)
1949  * amrex::max(qi_arr(i,j,k), Real(qmin));
1950  const Real temp = std::sqrt(std::sqrt(tmp * tmp * tmp));
1951  xni_arr(i,j,k) = amrex::min(
1952  amrex::max(Real(5.38e7) * temp, Real(1.0e3)),
1953  Real(1.0e6));
1954 
1955  const Real xni = xni_arr(i,j,k);
1956  const Real eacrs = std::exp(Real(0.07) * (-supcol));
1957 
1958  const Real xmi = den_arr(i,j,k) * qi_arr(i,j,k) / xni;
1959  const Real diameter = amrex::min(
1960  Real(dicon) * std::sqrt(xmi), Real(dimax));
1961  const Real vt2i = Real(1.49e4) * std::pow(diameter, Real(1.31));
1962  const Real vt2r = Real(pvtr) * rslopeb_r_arr(i,j,k)
1963  * denfac_arr(i,j,k);
1964  const Real vt2s = Real(pvts) * rslopeb_s_arr(i,j,k)
1965  * denfac_arr(i,j,k);
1966  const Real vt2g = Real(pvtg) * rslopeb_g_arr(i,j,k)
1967  * denfac_arr(i,j,k);
1968 
1969  // Recompute qsum from the current qs/qg (they have changed since
1970  // qsum_arr was set in G5a) -- mirrors ERF_module_mp_wsm6.F90 l.1008
1971  const Real qsum = amrex::max(qs_arr(i,j,k) + qg_arr(i,j,k), Real(1.0e-15));
1972  qsum_arr(i,j,k) = qsum;
1973  const Real vt2ave = (qsum > Real(1.0e-15))
1974  ? (vt2s * qs_arr(i,j,k) + vt2g * qg_arr(i,j,k)) / qsum
1975  : Real(0.0);
1976 
1977  if (supcol > Real(0.0) && qi_arr(i,j,k) > Real(qmin)) {
1978  if (qr_arr(i,j,k) > Real(qcrmin)) {
1979  const Real acrfac =
1980  Real(2.0) * rslope3_r_arr(i,j,k)
1981  + Real(2.0) * diameter * rslope2_r_arr(i,j,k)
1982  + diameter * diameter * rslope_r_arr(i,j,k);
1983 
1984  // WSM6-CPP TAG: PRACI
1985  // legacy_group: G13b
1986  // process: Accretion of cloud ice by rain
1987  // compare_vars: praci, qi, qr, den
1988  praci_arr(i,j,k) = Real(pi) * qi_arr(i,j,k) * Real(n0r)
1989  * std::abs(vt2r - vt2i) * acrfac / Real(4.0);
1990  praci_arr(i,j,k) *= std::pow(
1991  amrex::min(
1992  amrex::max(Real(0.0), qr_arr(i,j,k) / qi_arr(i,j,k)),
1993  Real(1.0)),
1994  Real(2.0));
1995  praci_arr(i,j,k) = amrex::min(
1996  praci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
1997 
1998  piacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(avtr)
1999  * Real(n0r) * Real(denr) * xni * denfac_arr(i,j,k)
2000  * Real(g6pbr) * rslope3_r_arr(i,j,k)
2001  * rslope3_r_arr(i,j,k) * rslopeb_r_arr(i,j,k)
2002  / Real(24.0) / den_arr(i,j,k);
2003  piacr_arr(i,j,k) *= std::pow(
2004  amrex::min(
2005  amrex::max(Real(0.0), qi_arr(i,j,k) / qr_arr(i,j,k)),
2006  Real(1.0)),
2007  Real(2.0));
2008  piacr_arr(i,j,k) = amrex::min(
2009  piacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2010  }
2011 
2012  if (qs_arr(i,j,k) > Real(qcrmin)) {
2013  const Real acrfac =
2014  Real(2.0) * rslope3_s_arr(i,j,k)
2015  + Real(2.0) * diameter * rslope2_s_arr(i,j,k)
2016  + diameter * diameter * rslope_s_arr(i,j,k);
2017  // WSM6-CPP TAG: PSACI
2018  // legacy_group: G13e
2019  // process: Accretion of cloud ice by snow
2020  // compare_vars: psaci, qi, qs, den
2021  psaci_arr(i,j,k) = Real(pi) * qi_arr(i,j,k) * eacrs
2022  * Real(n0s) * n0sfac_arr(i,j,k)
2023  * std::abs(vt2ave - vt2i) * acrfac / Real(4.0);
2024  psaci_arr(i,j,k) = amrex::min(
2025  psaci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
2026  }
2027 
2028  if (qg_arr(i,j,k) > Real(qcrmin)) {
2029  const Real egi = std::exp(Real(0.07) * (-supcol));
2030  const Real acrfac =
2031  Real(2.0) * rslope3_g_arr(i,j,k)
2032  + Real(2.0) * diameter * rslope2_g_arr(i,j,k)
2033  + diameter * diameter * rslope_g_arr(i,j,k);
2034  pgaci_arr(i,j,k) = Real(pi) * egi * qi_arr(i,j,k)
2035  * n0g_loc * std::abs(vt2ave - vt2i) * acrfac
2036  / Real(4.0);
2037  pgaci_arr(i,j,k) = amrex::min(
2038  pgaci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
2039  }
2040  }
2041 
2042  if (qs_arr(i,j,k) > Real(qcrmin) && qc_arr(i,j,k) > Real(qmin)) {
2043  psacw_arr(i,j,k) = amrex::min(
2044  Real(pacrc) * n0sfac_arr(i,j,k) * rslope3_s_arr(i,j,k)
2045  * rslopeb_s_arr(i,j,k)
2046  * std::pow(
2047  amrex::min(
2048  amrex::max(Real(0.0), qs_arr(i,j,k) / qc_arr(i,j,k)),
2049  Real(1.0)),
2050  Real(2.0))
2051  * qc_arr(i,j,k) * denfac_arr(i,j,k),
2052  qc_arr(i,j,k) / dtcld);
2053  }
2054 
2055  if (qg_arr(i,j,k) > Real(qcrmin) && qc_arr(i,j,k) > Real(qmin)) {
2056  pgacw_arr(i,j,k) = amrex::min(
2057  Real(pacrg) * rslope3_g_arr(i,j,k) * rslopeb_g_arr(i,j,k)
2058  * std::pow(
2059  amrex::min(
2060  amrex::max(Real(0.0), qg_arr(i,j,k) / qc_arr(i,j,k)),
2061  Real(1.0)),
2062  Real(2.0))
2063  * qc_arr(i,j,k) * denfac_arr(i,j,k),
2064  qc_arr(i,j,k) / dtcld);
2065  }
2066 
2067  if (qsum > Real(1.0e-15)) {
2068  paacw_arr(i,j,k) = (qs_arr(i,j,k) * psacw_arr(i,j,k)
2069  + qg_arr(i,j,k) * pgacw_arr(i,j,k))
2070  / qsum;
2071  }
2072 
2073  if (qs_arr(i,j,k) > Real(qcrmin) && qr_arr(i,j,k) > Real(qcrmin)) {
2074  if (supcol > Real(0.0)) {
2075  const Real acrfac =
2076  Real(5.0) * rslope3_s_arr(i,j,k) * rslope3_s_arr(i,j,k)
2077  * rslope_r_arr(i,j,k)
2078  + Real(2.0) * rslope3_s_arr(i,j,k) * rslope2_s_arr(i,j,k)
2079  * rslope2_r_arr(i,j,k)
2080  + Real(0.5) * rslope2_s_arr(i,j,k) * rslope2_s_arr(i,j,k)
2081  * rslope3_r_arr(i,j,k);
2082  // WSM6-CPP TAG: PRACS
2083  // legacy_group: G13f
2084  // process: Accretion of snow by rain / rain-snow interaction
2085  // compare_vars: pracs, qr, qs, den
2086  pracs_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
2087  * Real(n0s) * n0sfac_arr(i,j,k)
2088  * std::abs(vt2r - vt2ave) * (Real(dens_snow) / den_arr(i,j,k))
2089  * acrfac;
2090  pracs_arr(i,j,k) *= std::pow(
2091  amrex::min(
2092  amrex::max(Real(0.0), qr_arr(i,j,k) / qs_arr(i,j,k)),
2093  Real(1.0)),
2094  Real(2.0));
2095  pracs_arr(i,j,k) = amrex::min(
2096  pracs_arr(i,j,k), qs_arr(i,j,k) / dtcld);
2097  }
2098 
2099  {
2100  const Real acrfac =
2101  Real(5.0) * rslope3_r_arr(i,j,k) * rslope3_r_arr(i,j,k)
2102  * rslope_s_arr(i,j,k)
2103  + Real(2.0) * rslope3_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2104  * rslope2_s_arr(i,j,k)
2105  + Real(0.5) * rslope2_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2106  * rslope3_s_arr(i,j,k);
2107  psacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
2108  * Real(n0s) * n0sfac_arr(i,j,k)
2109  * std::abs(vt2ave - vt2r) * (Real(denr) / den_arr(i,j,k))
2110  * acrfac;
2111  psacr_arr(i,j,k) *= std::pow(
2112  amrex::min(
2113  amrex::max(Real(0.0), qs_arr(i,j,k) / qr_arr(i,j,k)),
2114  Real(1.0)),
2115  Real(2.0));
2116  psacr_arr(i,j,k) = amrex::min(
2117  psacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2118  }
2119  }
2120 
2121  // pgacr: accretion of rain by graupel [HL A12] [LFO 42]
2122  // (T<T0: R->G) (T>=T0: enhance melting of graupel)
2123  // This depends only on graupel and rain, so it must not be nested inside
2124  // the (qs,qr) test above -- otherwise a rain/graupel column with
2125  // negligible snow would never accrete rain onto graupel
2126  if (qg_arr(i,j,k) > Real(qcrmin) && qr_arr(i,j,k) > Real(qcrmin)) {
2127  const Real acrfac =
2128  Real(5.0) * rslope3_r_arr(i,j,k) * rslope3_r_arr(i,j,k)
2129  * rslope_g_arr(i,j,k)
2130  + Real(2.0) * rslope3_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2131  * rslope2_g_arr(i,j,k)
2132  + Real(0.5) * rslope2_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2133  * rslope3_g_arr(i,j,k);
2134  pgacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
2135  * n0g_loc * std::abs(vt2ave - vt2r)
2136  * (Real(denr) / den_arr(i,j,k)) * acrfac;
2137  pgacr_arr(i,j,k) *= std::pow(
2138  amrex::min(
2139  amrex::max(Real(0.0), qg_arr(i,j,k) / qr_arr(i,j,k)),
2140  Real(1.0)),
2141  Real(2.0));
2142  pgacr_arr(i,j,k) = amrex::min(
2143  pgacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2144  }
2145 
2146  if (qg_arr(i,j,k) > Real(qcrmin) && qs_arr(i,j,k) > Real(qcrmin)) {
2147  pgacs_arr(i,j,k) = Real(0.0);
2148  }
2149 
2150  if (supcol <= Real(0.0)) {
2151  const Real xlf = Real(xlf0);
2152  if (qs_arr(i,j,k) > Real(0.0)) {
2153  // WSM6-CPP TAG: PSEML
2154  // legacy_group: G13g
2155  // process: Snow evaporation/sublimation
2156  // compare_vars: pseml, qs, qv, qsat, den
2157  pseml_arr(i,j,k) = amrex::min(
2158  amrex::max(
2159  Real(cliq) * supcol
2160  * (paacw_arr(i,j,k) + psacr_arr(i,j,k)) / xlf,
2161  -qs_arr(i,j,k) / dtcld),
2162  Real(0.0));
2163  }
2164  if (qg_arr(i,j,k) > Real(0.0)) {
2165  pgeml_arr(i,j,k) = amrex::min(
2166  amrex::max(
2167  Real(cliq) * supcol
2168  * (paacw_arr(i,j,k) + pgacr_arr(i,j,k)) / xlf,
2169  -qg_arr(i,j,k) / dtcld),
2170  Real(0.0));
2171  }
2172  }
2173 
2174  if (supcol > Real(0.0)) {
2175  if (qi_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2176  // WSM6-CPP TAG: PIDEP
2177  // legacy_group: G13h
2178  // process: Ice deposition/sublimation
2179  // compare_vars: pidep, qi, qv, qsat, den, t
2180  pidep_arr(i,j,k) = Real(4.0) * diameter * xni
2181  * (rhi_arr(i,j,k) - Real(1.0))
2182  / workdiffi_arr(i,j,k);
2183  Real supice = satdt - prevp_arr(i,j,k);
2184  if (pidep_arr(i,j,k) < Real(0.0)) {
2185  pidep_arr(i,j,k) = amrex::max(
2186  amrex::max(pidep_arr(i,j,k), satdt / Real(2.0)),
2187  supice);
2188  pidep_arr(i,j,k) = amrex::max(
2189  pidep_arr(i,j,k), -qi_arr(i,j,k) / dtcld);
2190  } else {
2191  pidep_arr(i,j,k) = amrex::min(
2192  amrex::min(pidep_arr(i,j,k), satdt / Real(2.0)),
2193  supice);
2194  }
2195  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k))
2196  >= std::abs(satdt)) {
2197  ifsat = 1;
2198  }
2199  }
2200 
2201  if (qs_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2202  const Real coeres = rslope2_s_arr(i,j,k)
2203  * std::sqrt(rslope_s_arr(i,j,k)
2204  * rslopeb_s_arr(i,j,k));
2205  psdep_arr(i,j,k) = (rhi_arr(i,j,k) - Real(1.0))
2206  * n0sfac_arr(i,j,k)
2207  * (Real(precs1) * rslope2_s_arr(i,j,k)
2208  + Real(precs2) * work2_arr(i,j,k) * coeres)
2209  / workdiffi_arr(i,j,k);
2210  Real supice = satdt - prevp_arr(i,j,k) - pidep_arr(i,j,k);
2211  if (psdep_arr(i,j,k) < Real(0.0)) {
2212  psdep_arr(i,j,k) = amrex::max(
2213  psdep_arr(i,j,k), -qs_arr(i,j,k) / dtcld);
2214  psdep_arr(i,j,k) = amrex::max(
2215  amrex::max(psdep_arr(i,j,k), satdt / Real(2.0)),
2216  supice);
2217  } else {
2218  psdep_arr(i,j,k) = amrex::min(
2219  amrex::min(psdep_arr(i,j,k), satdt / Real(2.0)),
2220  supice);
2221  }
2222  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2223  + psdep_arr(i,j,k))
2224  >= std::abs(satdt)) {
2225  ifsat = 1;
2226  }
2227  }
2228 
2229  if (qg_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2230  const Real coeres = rslope2_g_arr(i,j,k)
2231  * std::sqrt(rslope_g_arr(i,j,k)
2232  * rslopeb_g_arr(i,j,k));
2233  pgdep_arr(i,j,k) = (rhi_arr(i,j,k) - Real(1.0))
2234  * (Real(precg1) * rslope2_g_arr(i,j,k)
2235  + Real(precg2) * work2_arr(i,j,k) * coeres)
2236  / workdiffi_arr(i,j,k);
2237  Real supice = satdt - prevp_arr(i,j,k)
2238  - pidep_arr(i,j,k) - psdep_arr(i,j,k);
2239  if (pgdep_arr(i,j,k) < Real(0.0)) {
2240  pgdep_arr(i,j,k) = amrex::max(
2241  pgdep_arr(i,j,k), -qg_arr(i,j,k) / dtcld);
2242  pgdep_arr(i,j,k) = amrex::max(
2243  amrex::max(pgdep_arr(i,j,k), satdt / Real(2.0)),
2244  supice);
2245  } else {
2246  pgdep_arr(i,j,k) = amrex::min(
2247  amrex::min(pgdep_arr(i,j,k), satdt / Real(2.0)),
2248  supice);
2249  }
2250  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2251  + psdep_arr(i,j,k) + pgdep_arr(i,j,k))
2252  >= std::abs(satdt)) {
2253  ifsat = 1;
2254  }
2255  }
2256 
2257  if (supsat > Real(0.0) && ifsat != 1) {
2258  const Real supice = satdt - prevp_arr(i,j,k)
2259  - pidep_arr(i,j,k)
2260  - psdep_arr(i,j,k)
2261  - pgdep_arr(i,j,k);
2262  const Real xni0 = Real(1.0e3) * std::exp(Real(0.1) * supcol);
2263  const Real roqi0 = Real(4.92e-11) * std::pow(xni0, Real(1.33));
2264  pigen_arr(i,j,k) = amrex::max(
2265  Real(0.0),
2266  (roqi0 / den_arr(i,j,k)
2267  - amrex::max(qi_arr(i,j,k), Real(0.0))) / dtcld);
2268  pigen_arr(i,j,k) = amrex::min(
2269  amrex::min(pigen_arr(i,j,k), satdt), supice);
2270  }
2271 
2272  if (qi_arr(i,j,k) > Real(0.0)) {
2273  const Real qimax = Real(roqimax) / den_arr(i,j,k);
2274  // WSM6-CPP TAG: PSAUT
2275  // legacy_group: G13i
2276  // process: Autoconversion to snow
2277  // compare_vars: psaut, qi, qs, den
2278  psaut_arr(i,j,k) = amrex::max(
2279  Real(0.0), (qi_arr(i,j,k) - qimax) / dtcld);
2280  }
2281 
2282  if (qs_arr(i,j,k) > Real(0.0)) {
2283  const Real alpha2 = Real(1.0e-3)
2284  * std::exp(Real(0.09) * (-supcol));
2285  pgaut_arr(i,j,k) = amrex::min(
2286  amrex::max(
2287  Real(0.0), alpha2 * (qs_arr(i,j,k) - Real(qs0))),
2288  qs_arr(i,j,k) / dtcld);
2289  }
2290  }
2291 
2292  if (supcol < Real(0.0)) {
2293  if (qs_arr(i,j,k) > Real(0.0)
2294  && rhw_arr(i,j,k) < Real(1.0)) {
2295  const Real coeres = rslope2_s_arr(i,j,k)
2296  * std::sqrt(rslope_s_arr(i,j,k)
2297  * rslopeb_s_arr(i,j,k));
2298  // WSM6-CPP TAG: PSEVP
2299  // legacy_group: G13j
2300  // process: Graupel evaporation/sublimation
2301  // compare_vars: psevp, qg, qv, qsat, den
2302  psevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2303  * n0sfac_arr(i,j,k)
2304  * (Real(precs1) * rslope2_s_arr(i,j,k)
2305  + Real(precs2) * work2_arr(i,j,k) * coeres)
2306  / workdiffw_arr(i,j,k);
2307  psevp_arr(i,j,k) = amrex::min(
2308  amrex::max(psevp_arr(i,j,k),
2309  -qs_arr(i,j,k) / dtcld),
2310  Real(0.0));
2311  }
2312 
2313  if (qg_arr(i,j,k) > Real(0.0)
2314  && rhw_arr(i,j,k) < Real(1.0)) {
2315  const Real coeres = rslope2_g_arr(i,j,k)
2316  * std::sqrt(rslope_g_arr(i,j,k)
2317  * rslopeb_g_arr(i,j,k));
2318  pgevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2319  * (Real(precg1) * rslope2_g_arr(i,j,k)
2320  + Real(precg2) * work2_arr(i,j,k) * coeres)
2321  / workdiffw_arr(i,j,k);
2322  pgevp_arr(i,j,k) = amrex::min(
2323  amrex::max(pgevp_arr(i,j,k),
2324  -qg_arr(i,j,k) / dtcld),
2325  Real(0.0));
2326  }
2327  }
2328  });
2329  // G14: mass conservation check and state update [lines 1200-1388]
2330  // WSM6-CPP TAG: UPDATE
2331  // legacy_group: G14
2332  // process: Mass conservation and state update
2333  // compare_vars: t, q, qci, qrs, qv
2334  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2335  const Real qmin_l = Real(qmin);
2336  const Real qcrmin_l= Real(qcrmin);
2337  const Real t0c_l = Real(t0c);
2338 
2339  const Real delta2 =
2340  (qr_arr(i,j,k) < Real(1.0e-4) && qs_arr(i,j,k) < Real(1.0e-4))
2341  ? Real(1.0) : Real(0.0);
2342  const Real delta3 =
2343  (qr_arr(i,j,k) < Real(1.0e-4)) ? Real(1.0) : Real(0.0);
2344 
2345  if (t_arr(i,j,k) <= t0c_l) {
2346  Real value, source, factor, xlf, xlwork2;
2347 
2348  value = amrex::max(qmin_l, qc_arr(i,j,k));
2349  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2350  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2351 // + paacw_arr(i,j,k)) * dtcld;
2352  if (source > value) {
2353  factor = value / source;
2354  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2355  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2356  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2357  }
2358 
2359  value = amrex::max(qmin_l, qi_arr(i,j,k));
2360  source = (psaut_arr(i,j,k) - pigen_arr(i,j,k)
2361  - pidep_arr(i,j,k) + praci_arr(i,j,k)
2362  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)) * dtcld;
2363  if (source > value) {
2364  factor = value / source;
2365  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2366  pigen_arr(i,j,k) = pigen_arr(i,j,k) * factor;
2367  pidep_arr(i,j,k) = pidep_arr(i,j,k) * factor;
2368  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2369  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2370  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2371  }
2372 
2373  value = amrex::max(qmin_l, qr_arr(i,j,k));
2374  source = (-praut_arr(i,j,k) - prevp_arr(i,j,k)
2375  - pracw_arr(i,j,k) + piacr_arr(i,j,k)
2376  + psacr_arr(i,j,k) + pgacr_arr(i,j,k)) * dtcld;
2377  if (source > value) {
2378  factor = value / source;
2379  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2380  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2381  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2382  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2383  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2384  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2385  }
2386 
2387  value = amrex::max(qmin_l, qs_arr(i,j,k));
2388  source = -(psdep_arr(i,j,k) + psaut_arr(i,j,k)
2389  - pgaut_arr(i,j,k) + paacw_arr(i,j,k)
2390  + piacr_arr(i,j,k) * delta3
2391  + praci_arr(i,j,k) * delta3
2392  - pracs_arr(i,j,k) * (Real(1.0) - delta2)
2393  + psacr_arr(i,j,k) * delta2
2394  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)) * dtcld;
2395  if (source > value) {
2396  factor = value / source;
2397  psdep_arr(i,j,k) = psdep_arr(i,j,k) * factor;
2398  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2399  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2400  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2401  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2402  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2403  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2404  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2405  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2406  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2407  }
2408 
2409  value = amrex::max(qmin_l, qg_arr(i,j,k));
2410  source = -(pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2411  + piacr_arr(i,j,k) * (Real(1.0) - delta3)
2412  + praci_arr(i,j,k) * (Real(1.0) - delta3)
2413  + psacr_arr(i,j,k) * (Real(1.0) - delta2)
2414  + pracs_arr(i,j,k) * (Real(1.0) - delta2)
2415  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2416  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k)) * dtcld;
2417  if (source > value) {
2418  factor = value / source;
2419  pgdep_arr(i,j,k) = pgdep_arr(i,j,k) * factor;
2420  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2421  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2422  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2423  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2424  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2425  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2426  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2427  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2428  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2429  }
2430 
2431  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psdep_arr(i,j,k)
2432  + pgdep_arr(i,j,k) + pigen_arr(i,j,k)
2433  + pidep_arr(i,j,k));
2434  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
2435 // + paacw_arr(i,j,k))
2436  qc_arr(i,j,k) = amrex::max(
2437  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
2438  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
2439  * dtcld,
2440  Real(0.0));
2441  qr_arr(i,j,k) = amrex::max(
2442  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
2443  + prevp_arr(i,j,k) - piacr_arr(i,j,k)
2444  - pgacr_arr(i,j,k) - psacr_arr(i,j,k))
2445  * dtcld,
2446  Real(0.0));
2447  qi_arr(i,j,k) = amrex::max(
2448  qi_arr(i,j,k) - (psaut_arr(i,j,k) + praci_arr(i,j,k)
2449  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)
2450  - pigen_arr(i,j,k) - pidep_arr(i,j,k))
2451  * dtcld,
2452  Real(0.0));
2453  qs_arr(i,j,k) = amrex::max(
2454  qs_arr(i,j,k) + (psdep_arr(i,j,k) + psaut_arr(i,j,k)
2455  + paacw_arr(i,j,k) - pgaut_arr(i,j,k)
2456  + piacr_arr(i,j,k) * delta3
2457  + praci_arr(i,j,k) * delta3
2458  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)
2459  - pracs_arr(i,j,k) * (Real(1.0) - delta2)
2460  + psacr_arr(i,j,k) * delta2) * dtcld,
2461  Real(0.0));
2462  qg_arr(i,j,k) = amrex::max(
2463  qg_arr(i,j,k) + (pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2464  + piacr_arr(i,j,k) * (Real(1.0) - delta3)
2465  + praci_arr(i,j,k) * (Real(1.0) - delta3)
2466  + psacr_arr(i,j,k) * (Real(1.0) - delta2)
2467  + pracs_arr(i,j,k) * (Real(1.0) - delta2)
2468  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2469  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k))
2470  * dtcld,
2471  Real(0.0));
2472  xlf = Real(xls) - xl_arr(i,j,k);
2473  xlwork2 = -Real(xls) * (psdep_arr(i,j,k) + pgdep_arr(i,j,k)
2474  + pidep_arr(i,j,k) + pigen_arr(i,j,k))
2475  - xl_arr(i,j,k) * prevp_arr(i,j,k)
2476  - xlf * (piacr_arr(i,j,k) + paacw_arr(i,j,k)
2477  + paacw_arr(i,j,k) + pgacr_arr(i,j,k)
2478  + psacr_arr(i,j,k));
2479  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
2480  } else {
2481  Real value, source, factor, xlf, xlwork2;
2482 
2483  value = amrex::max(qmin_l, qc_arr(i,j,k));
2484  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2485  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2486 // + paacw_arr(i,j,k)) * dtcld;
2487  if (source > value) {
2488  factor = value / source;
2489  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2490  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2491  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2492  }
2493 
2494  value = amrex::max(qmin_l, qr_arr(i,j,k));
2495  source = (-paacw_arr(i,j,k) - praut_arr(i,j,k)
2496  + pseml_arr(i,j,k) + pgeml_arr(i,j,k)
2497  - pracw_arr(i,j,k) - paacw_arr(i,j,k)
2498  - prevp_arr(i,j,k)) * dtcld;
2499  if (source > value) {
2500  factor = value / source;
2501  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2502  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2503  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2504  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2505  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2506  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2507  }
2508 
2509  value = amrex::max(qcrmin_l, qs_arr(i,j,k));
2510  source = (pgacs_arr(i,j,k) - pseml_arr(i,j,k)
2511  - psevp_arr(i,j,k)) * dtcld;
2512  if (source > value) {
2513  factor = value / source;
2514  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2515  psevp_arr(i,j,k) = psevp_arr(i,j,k) * factor;
2516  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2517  }
2518 
2519  value = amrex::max(qcrmin_l, qg_arr(i,j,k));
2520  source = -(pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
2521  + pgeml_arr(i,j,k)) * dtcld;
2522  if (source > value) {
2523  factor = value / source;
2524  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2525  pgevp_arr(i,j,k) = pgevp_arr(i,j,k) * factor;
2526  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2527  }
2528 
2529  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psevp_arr(i,j,k)
2530  + pgevp_arr(i,j,k));
2531  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
2532 // + paacw_arr(i,j,k))
2533  qc_arr(i,j,k) = amrex::max(
2534  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
2535  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
2536  * dtcld,
2537  Real(0.0));
2538  qr_arr(i,j,k) = amrex::max(
2539  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
2540  + prevp_arr(i,j,k) + paacw_arr(i,j,k)
2541  + paacw_arr(i,j,k) - pseml_arr(i,j,k)
2542  - pgeml_arr(i,j,k)) * dtcld,
2543  Real(0.0));
2544  qs_arr(i,j,k) = amrex::max(
2545  qs_arr(i,j,k) + (psevp_arr(i,j,k) - pgacs_arr(i,j,k)
2546  + pseml_arr(i,j,k)) * dtcld,
2547  Real(0.0));
2548  qg_arr(i,j,k) = amrex::max(
2549  qg_arr(i,j,k) + (pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
2550  + pgeml_arr(i,j,k)) * dtcld,
2551  Real(0.0));
2552  xlf = Real(xls) - xl_arr(i,j,k);
2553  xlwork2 = -xl_arr(i,j,k) * (prevp_arr(i,j,k)
2554  + psevp_arr(i,j,k)
2555  + pgevp_arr(i,j,k))
2556  - xlf * (pseml_arr(i,j,k) + pgeml_arr(i,j,k));
2557  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
2558  }
2559  });
2560  // G15: second qsat computation [lines 1390-1420]
2561  // WSM6-CPP TAG: QSAT2
2562  // legacy_group: G15
2563  // process: Second saturation mixing ratio computation
2564  // compare_vars: qs, qvs, den, denfac, t, p
2565  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2566  const Real ttp = Real(t0c) + Real(0.01);
2567  const Real dldt = Real(cpv) - Real(cliq);
2568  const Real xa = -dldt / Real(rv);
2569  const Real xb = xa + Real(xlv0) / (Real(rv) * ttp);
2570 
2571  Real tr = ttp / t_arr(i,j,k);
2572  Real qsw = Real(psat) * std::exp(std::log(tr) * xa)
2573  * std::exp(xb * (Real(1.0) - tr));
2574  qsw = amrex::min(qsw, Real(0.99) * p_arr(i,j,k));
2575  qsatw_arr(i,j,k) = Real(ep2) * qsw / (p_arr(i,j,k) - qsw);
2576  qsatw_arr(i,j,k) = amrex::max(qsatw_arr(i,j,k), Real(qmin));
2577  });
2578  // G16: pcond condensational/evaporational update [lines 1427-1437]
2579  // WSM6-CPP TAG: PCOND
2580  // legacy_group: G16
2581  // process: Condensation/evaporation update
2582  // compare_vars: pcond, t, qv, qc, qsat
2583  if (m_do_cond) {
2584  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2585  const Real workcond = wsm6_conden(
2586  t_arr(i,j,k), qv_arr(i,j,k), qsatw_arr(i,j,k),
2587  xl_arr(i,j,k), cpm_arr(i,j,k), Real(qmin), Real(rv));
2588  const Real work2loc = qc_arr(i,j,k) + workcond;
2589  static_cast<void>(work2loc);
2590  pcond_arr(i,j,k) = amrex::min(
2591  amrex::max(workcond / dtcld, Real(0.0)),
2592  amrex::max(qv_arr(i,j,k), Real(0.0)) / dtcld);
2593  if (qc_arr(i,j,k) > Real(0.0) && workcond < Real(0.0)) {
2594  pcond_arr(i,j,k) = amrex::max(workcond, -qc_arr(i,j,k)) / dtcld;
2595  }
2596  qv_arr(i,j,k) = qv_arr(i,j,k) - pcond_arr(i,j,k) * dtcld;
2597  qc_arr(i,j,k) = amrex::max(
2598  qc_arr(i,j,k) + pcond_arr(i,j,k) * dtcld,
2599  Real(0.0));
2600  t_arr(i,j,k) = t_arr(i,j,k)
2601  + pcond_arr(i,j,k) * xl_arr(i,j,k)
2602  / cpm_arr(i,j,k) * dtcld;
2603  });
2604  }
2605 
2606  // G17: padding for small values [lines 1444-1449]
2607  // WSM6-CPP TAG: CLIP
2608  // legacy_group: G17
2609  // process: Padding/clipping for small values
2610  // compare_vars: qv, qc, qr, qi, qs, qg
2611  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2612  if (qc_arr(i,j,k) <= Real(qmin)) qc_arr(i,j,k) = Real(0.0);
2613  if (qi_arr(i,j,k) <= Real(qmin)) qi_arr(i,j,k) = Real(0.0);
2614  });
2615 
2616  }
2617 #ifdef ERF_USE_WSM6_FORT
2618  }
2619 #endif
2620  // Only a column standing on the domain bottom has a ground to accumulate precipitation on
2621  if (klo == m_geom.Domain().smallEnd(2)) {
2622  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
2623  rain_arr(i,j,klo) += rainonly_arr(i,j,0);
2624  snow_arr(i,j,klo) += snowonly_arr(i,j,0);
2625  graupel_arr(i,j,klo) += graupacc_arr(i,j,0);
2626  precip_arr(i,j,klo) += rainacc_arr(i,j,0);
2627  });
2628  }
2629  }
2630 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_xka(Real x, Real y)
Definition: ERF_AdvanceWSM6.cpp:46
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_conden(Real a, Real b, Real c, Real d, Real e, Real qmin_arg, Real rv_arg)
Definition: ERF_AdvanceWSM6.cpp:66
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_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_AdvanceWSM6.cpp:172
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_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_AdvanceWSM6.cpp:202
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_slope_rain_cell(Real qr, Real den, Real denfac, Real pidn0r_arg, Real qcrmin_arg, Real rslopermax_arg, Real rsloperbmax_arg, Real rsloper2max_arg, Real rsloper3max_arg, Real bvtr_arg, Real pvtr_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt)
Definition: ERF_AdvanceWSM6.cpp:147
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_nislfv_rain_plm_scratch(int km, int ww_comp, int rq_comp, Real *precip, Real dt, int iter, Array4< Real > const &sed_cell, Array4< Real > const &sed_node, int i_s, int j_s, int klo_s)
Definition: ERF_AdvanceWSM6.cpp:226
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_xlcal(Real x, Real xlv0_arg, Real xlv1_arg, Real t0c_arg)
Definition: ERF_AdvanceWSM6.cpp:31
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_venfac(Real a, Real b, Real c, Real den0_arg)
Definition: ERF_AdvanceWSM6.cpp:58
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_cpmcal(Real x, Real qmin_arg, Real cpd_arg, Real cpv_arg)
Definition: ERF_AdvanceWSM6.cpp:25
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_nislfv_rain_plm6_scratch(int km, int ww_comp, int rq_comp, int rq2_comp, Real *precip1, Real *precip2, Real dt, int iter, Array4< Real > const &sed_cell, Array4< Real > const &sed_node, int i_s, int j_s, int klo_s)
Definition: ERF_AdvanceWSM6.cpp:508
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_diffac(Real a, Real b, Real c, Real d, Real e, Real rv_arg)
Definition: ERF_AdvanceWSM6.cpp:51
constexpr amrex::Real R_v
Definition: ERF_Constants.H:50
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:51
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:72
constexpr amrex::Real Cp_l
Definition: ERF_Constants.H:53
constexpr amrex::Real Cp_v
Definition: ERF_Constants.H:52
constexpr amrex::Real R_d
Definition: ERF_Constants.H:49
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(axis >=0 &&axis< AMREX_SPACEDIM, "prob.axis must be 0, 1 or 2")
ParmParse pp("prob")
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
const int klo
Definition: ERF_InitCustomPert_KE.H:9
auto qv_arr
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:222
constexpr amrex::Real lat_vap
Definition: ERF_MicrophysicsConstants.H:113
constexpr amrex::Real lsub
Definition: ERF_MicrophysicsConstants.H:111
constexpr amrex::Real rhoh2o
Definition: ERF_MicrophysicsConstants.H:42
constexpr amrex::Real lat_ice
Definition: ERF_MicrophysicsConstants.H:114
constexpr amrex::Real rhos
Definition: ERF_MicrophysicsConstants.H:40
Arena * Arena_Used
Definition: ERF_Morrison_Advance_F.H:23
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);})
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void sedimentation_add_column_sums(const amrex::Array4< amrex::Real > &sed, int n, int i, int j, int klo, int khi, amrex::Real sign, F &&w) noexcept
Definition: ERF_SedimentationFlux.H:74
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
void mp_wsm6_init_c(double den0, double denr, double dens, double cl, double cpv, int hail_opt)
void mp_wsm6_run_c(double *t, double *qv, double *qc, double *qi, double *qr, double *qs, double *qg, 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 *rain, double *rainncv, double *sr, double *snow, double *snowncv, double *graupel, double *graupelncv, double *rainonly, double *snowonly, 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)
bool boxes_span_domain() const
True when every box spans the domain in z (required by the Fortran bridges)
Definition: ERF_MicrophysicsColumns.H:44
amrex::Real m_roqimax
Definition: ERF_WSM6.H:276
amrex::Real m_rslopesbmax
Definition: ERF_WSM6.H:285
static constexpr amrex::Real avtr
Definition: ERF_WSM6.H:147
static constexpr amrex::Real dicon
Definition: ERF_WSM6.H:157
amrex::Real m_rslopes2max
Definition: ERF_WSM6.H:286
amrex::Real m_rsloper2max
Definition: ERF_WSM6.H:286
static constexpr amrex::Real qcrmin
Definition: ERF_WSM6.H:161
amrex::Real m_rslopermax
Definition: ERF_WSM6.H:284
amrex::Real m_pacrc
Definition: ERF_WSM6.H:280
static constexpr amrex::Real pfrz1
Definition: ERF_WSM6.H:159
amrex::Real m_rsloperbmax
Definition: ERF_WSM6.H:285
amrex::Real m_pi_wsm6
Definition: ERF_WSM6.H:271
static constexpr amrex::Real n0smax
Definition: ERF_WSM6.H:165
bool m_do_cond
Definition: ERF_WSM6.H:248
amrex::Real m_pidn0s
Definition: ERF_WSM6.H:280
static constexpr amrex::Real dimax
Definition: ERF_WSM6.H:158
static constexpr amrex::Real pfrz2
Definition: ERF_WSM6.H:160
static constexpr amrex::Real alpha_wsm6
Definition: ERF_WSM6.H:167
amrex::Real m_pvtg
Definition: ERF_WSM6.H:283
amrex::Real m_qck1
Definition: ERF_WSM6.H:272
amrex::Real m_precg2
Definition: ERF_WSM6.H:283
amrex::Real m_precg1
Definition: ERF_WSM6.H:283
amrex::Real m_rsloper3max
Definition: ERF_WSM6.H:287
static constexpr amrex::Real n0s
Definition: ERF_WSM6.H:166
amrex::Real m_precr2
Definition: ERF_WSM6.H:276
static constexpr amrex::Real n0r
Definition: ERF_WSM6.H:146
amrex::Real m_n0g
Definition: ERF_WSM6.H:266
amrex::Real m_pidn0g
Definition: ERF_WSM6.H:283
static constexpr amrex::Real dens_snow
Definition: ERF_WSM6.H:163
static constexpr amrex::Real dtcldcr
Definition: ERF_WSM6.H:145
amrex::Real m_rslopegmax
Definition: ERF_WSM6.H:284
static constexpr amrex::Real bvts
Definition: ERF_WSM6.H:154
amrex::Real m_g6pbr
Definition: ERF_WSM6.H:274
amrex::Real m_pacrg
Definition: ERF_WSM6.H:283
amrex::Real m_precs2
Definition: ERF_WSM6.H:279
amrex::Real m_rslopeg3max
Definition: ERF_WSM6.H:287
amrex::Real m_pvts
Definition: ERF_WSM6.H:279
amrex::Real m_rslopegbmax
Definition: ERF_WSM6.H:285
static constexpr amrex::Real qs0
Definition: ERF_WSM6.H:164
amrex::Geometry m_geom
Definition: ERF_WSM6.H:243
amrex::Real m_precs1
Definition: ERF_WSM6.H:279
amrex::Real dt
Definition: ERF_WSM6.H:244
static constexpr amrex::Real xncr
Definition: ERF_WSM6.H:151
static constexpr amrex::Real bvtr
Definition: ERF_WSM6.H:148
amrex::Real m_rslopesmax
Definition: ERF_WSM6.H:284
amrex::Real m_qc0
Definition: ERF_WSM6.H:272
amrex::Real m_bvtg
Definition: ERF_WSM6.H:266
amrex::Real m_pvtr
Definition: ERF_WSM6.H:275
amrex::Real m_pacrr
Definition: ERF_WSM6.H:275
amrex::Real m_precr1
Definition: ERF_WSM6.H:276
amrex::Real m_pidn0r
Definition: ERF_WSM6.H:280
amrex::Real m_rslopes3max
Definition: ERF_WSM6.H:287
amrex::Real m_xlv1
Definition: ERF_WSM6.H:271
amrex::Real m_rslopeg2max
Definition: ERF_WSM6.H:286
@ xlf
Definition: ERF_AdvanceMorrison.cpp:159
@ precip_accum
Definition: ERF_WSM6.H:36
@ qr
Definition: ERF_WSM6.H:30
@ qi
Definition: ERF_WSM6.H:29
@ graupel_accum
Definition: ERF_WSM6.H:35
@ qs
Definition: ERF_WSM6.H:31
@ qv
Definition: ERF_WSM6.H:27
@ qc
Definition: ERF_WSM6.H:28
@ rain_accum
Definition: ERF_WSM6.H:33
@ snow_accum
Definition: ERF_WSM6.H:34
@ rho
Definition: ERF_WSM6.H:23
@ qg
Definition: ERF_WSM6.H:32
@ pres
Definition: ERF_WSM6.H:26
@ tabs
Definition: ERF_WSM6.H:25
@ precip
Definition: ERF_SLM.H:133
@ qsum
Definition: ERF_WSM6.H:344
@ xni
Definition: ERF_WSM6.H:349
@ NumComps
Definition: ERF_AdvanceWSM6.cpp:128
@ workr_col
Definition: ERF_AdvanceWSM6.cpp:120
@ tmp
Definition: ERF_AdvanceWSM6.cpp:116
@ denqrs2_col
Definition: ERF_AdvanceWSM6.cpp:123
@ qsum_col
Definition: ERF_AdvanceWSM6.cpp:125
@ denqci_col
Definition: ERF_AdvanceWSM6.cpp:127
@ denqrs1_col
Definition: ERF_AdvanceWSM6.cpp:122
@ worka_col
Definition: ERF_AdvanceWSM6.cpp:121
@ denqrs3_col
Definition: ERF_AdvanceWSM6.cpp:124
@ den
Definition: ERF_AdvanceWSM6.cpp:111
@ dz
Definition: ERF_AdvanceWSM6.cpp:106
@ tk
Definition: ERF_AdvanceWSM6.cpp:113
@ denfac
Definition: ERF_AdvanceWSM6.cpp:112
@ work1c_col
Definition: ERF_AdvanceWSM6.cpp:126
@ NumComps
Definition: ERF_AdvanceWSM6.cpp:142
@ qs
Definition: ERF_WSM6.H:125
@ qr
Definition: ERF_WSM6.H:125
@ qg
Definition: ERF_WSM6.H:125
@ NumComps
Definition: ERF_WSM6.H:125
@ qi
Definition: ERF_WSM6.H:125
real(c_double), parameter cice
Definition: ERF_module_model_constants.F90:30
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
real(c_double), parameter cliq
Definition: ERF_module_model_constants.F90:29
real(c_double), parameter xlv0
Definition: ERF_module_model_constants.F90:51
real(c_double), parameter cpv
Definition: ERF_module_model_constants.F90:26
real(c_double), parameter xls
Definition: ERF_module_model_constants.F90:56
real(c_double), parameter psat
Definition: ERF_module_model_constants.F90:31
real(c_double), parameter, private pi
Definition: ERF_module_mp_morr_two_moment.F90:100
real(kind=kind_phys), save pacrc
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precs1
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save pacrr
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precr2
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save pvtr
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precr1
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save pvtg
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save g6pbr
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), parameter, private dens
Definition: ERF_module_mp_wdm6.F90:61
real(kind=kind_phys), save roqimax
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save qck1
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save qc0
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save pvts
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precs2
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save pacrg
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precg1
Definition: ERF_module_mp_wdm6.F90:100
real(kind=kind_phys), save precg2
Definition: ERF_module_mp_wdm6.F90:100
Here is the call graph for this function:

◆ Copy_Micro_to_State()

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

Reimplemented from NullMoist.

8 {
9  for (MFIter mfi(cons, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
10  const auto& box3d = mfi.tilebox();
11  auto states = cons.array(mfi);
12 
13  auto rho = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
14  auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi);
15  auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
16  auto pres = mic_fab_vars[MicVar_WSM6::pres]->const_array(mfi);
17 
18  auto qv = mic_fab_vars[MicVar_WSM6::qv]->array(mfi);
19  auto qc = mic_fab_vars[MicVar_WSM6::qc]->array(mfi);
20  auto qi = mic_fab_vars[MicVar_WSM6::qi]->array(mfi);
21  auto qr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi);
22  auto qs = mic_fab_vars[MicVar_WSM6::qs]->array(mfi);
23  auto qg = mic_fab_vars[MicVar_WSM6::qg]->array(mfi);
24  const Real rdOcp = m_rdOcp;
25  const bool use_anelastic_reference_pressure =
27 
28  ParallelFor(box3d, [=]
29  AMREX_GPU_DEVICE(int i, int j, int k) {
31  states, theta, rho, tabs, pres, qv, qc, qi, qr, qs, qg,
32  use_anelastic_reference_pressure, rdOcp, i, j, k);
33  });
34  }
35 
36  cons.FillBoundary(m_geom.periodicity());
37 }
const Real rdOcp
Definition: ERF_InitCustomPert_ABL.H:72
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_copy_micro_to_state_cell(const amrex::Array4< amrex::Real > &states, const amrex::Array4< amrex::Real > &theta, const amrex::Array4< const amrex::Real > &rho, const amrex::Array4< const amrex::Real > &tabs, const amrex::Array4< const amrex::Real > &pres, const amrex::Array4< const amrex::Real > &qv, const amrex::Array4< const amrex::Real > &qc, const amrex::Array4< const amrex::Real > &qi, const amrex::Array4< const amrex::Real > &qr, const amrex::Array4< const amrex::Real > &qs, const amrex::Array4< const amrex::Real > &qg, const bool use_anelastic_reference_pressure, const amrex::Real rdOcp, const int i, const int j, const int k) noexcept
Definition: ERF_WSM6.H:94
bool m_use_anelastic_reference_pressure
Definition: ERF_NullMoist.H:191
amrex::Real m_rdOcp
Definition: ERF_WSM6.H:249
@ tabs
Definition: ERF_Kessler.H:27
@ pres
Definition: ERF_Kessler.H:28
@ rho
Definition: ERF_Kessler.H:25
@ theta
Definition: ERF_Kessler.H:26
@ qv
Definition: ERF_Kessler.H:31
@ theta
Definition: ERF_WSM6.H:24
@ qg
Definition: ERF_Morrison.H:114
@ qr
Definition: ERF_Morrison.H:112
@ qi
Definition: ERF_Morrison.H:111
@ qs
Definition: ERF_Morrison.H:113
@ qc
Definition: ERF_Morrison.H:110
@ cons
Definition: ERF_IndexDefines.H:254

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro() [1/2]

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

Reimplemented from NullMoist.

Referenced by Update_Micro_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro() [2/2]

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

◆ Define()

void WSM6::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

136  {
138  m_rdOcp = sc.rdOcp;
140  m_axis = sc.ave_plane;
141  m_do_cond = (!sc.uses_shoc_family());
142  }
void set_anelastic_reference_pressure_mode(const SolverChoice &sc)
Definition: ERF_NullMoist.H:171
int m_axis
Definition: ERF_WSM6.H:247
MoistureType m_moisture_type
Definition: ERF_WSM6.H:250
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2942
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2990
amrex::Real rdOcp
Ratio of dry-air gas constant to c_p.
Definition: ERF_DataStruct.H:2738
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:2864
Here is the call graph for this function:

◆ Get_Surface_Precip_Accumulation_Ptrs()

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

Reimplemented from NullMoist.

222  {
224  sources.total = {mic_fab_vars[MicVar_WSM6::precip_accum].get(), rhoh2o / amrex::Real(1000.0)};
225  sources.rain = {mic_fab_vars[MicVar_WSM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
226  sources.snow = {mic_fab_vars[MicVar_WSM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
227  sources.graupel = {mic_fab_vars[MicVar_WSM6::graupel_accum].get(), rhoh2o / amrex::Real(1000.0)};
228  return sources;
229  }
Definition: ERF_SurfacePrecipitation.H:34
SurfacePrecipAccumulationSource rain
Definition: ERF_SurfacePrecipitation.H:36
SurfacePrecipAccumulationSource snow
Definition: ERF_SurfacePrecipitation.H:37
SurfacePrecipAccumulationSource total
Definition: ERF_SurfacePrecipitation.H:35
SurfacePrecipAccumulationSource graupel
Definition: ERF_SurfacePrecipitation.H:38

◆ Init()

void WSM6::Init ( const amrex::MultiFab &  cons_in,
const amrex::BoxArray &  grids,
const amrex::Geometry &  geom,
const amrex::Real &  dt_advance,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
std::unique_ptr< amrex::MultiFab > &  detJ_cc 
)
overridevirtual

Reimplemented from NullMoist.

17 {
18  dt = dt_advance;
19  m_geom = geom;
20 
21  m_z_phys_nd = z_phys_nd.get();
22  m_detJ_cc = detJ_cc.get();
23 
24  m_columns.define(cons_in.boxArray(), geom.Domain());
25 
26  MicVarMap.resize(m_qmoist_size);
29 
30 #if defined(ERF_USE_WSM6_FORT) && defined(AMREX_USE_GPU)
31  // The host-only Fortran bridge receives dataPtr() from mic_fab_vars.
32  Arena* Arena_Used = The_Managed_Arena();
33 #else
34  Arena* Arena_Used = The_Arena();
35 #endif
36 
37  for (int ivar = 0; ivar < MicVar_WSM6::NumVars; ++ivar) {
38  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
39  1, cons_in.nGrowVect(),
40  MFInfo().SetArena(Arena_Used));
41  mic_fab_vars[ivar]->setVal(0.0);
42  }
43 
44  m_sed.define(cons_in.boxArray(), cons_in.DistributionMap(),
45  {RhoQ3_comp, RhoQ4_comp, RhoQ5_comp, RhoQ6_comp});
46 
47  nlev = m_geom.Domain().length(2);
48  zlo = m_geom.Domain().smallEnd(2);
49  zhi = m_geom.Domain().bigEnd(2);
50 
52 }
void define(const amrex::BoxArray &ba, const amrex::Box &domain)
Build from a level's grids; call from the scheme's Init (which also runs on regrid)
Definition: ERF_MicrophysicsColumns.cpp:9
amrex::MultiFab * m_detJ_cc
Definition: ERF_WSM6.H:253
int zlo
Definition: ERF_WSM6.H:246
int nlev
Definition: ERF_WSM6.H:246
int zhi
Definition: ERF_WSM6.H:246
int m_qmoist_size
Definition: ERF_WSM6.H:237
void initialize_coeffs()
Definition: ERF_InitWSM6.cpp:103
amrex::Vector< int > MicVarMap
Definition: ERF_WSM6.H:241
@ NumVars
Definition: ERF_WSM6.H:37
void define(const amrex::BoxArray &cc_ba, const amrex::DistributionMapping &dm, const amrex::Vector< int > &a_cons_comp)
Definition: ERF_SedimentationFlux.H:43

◆ initialize_coeffs()

void WSM6::initialize_coeffs ( )
private
104 {
105  using amrex::Real;
106 
107  // Exact port of Fortran rgmma() — ERF_module_mp_wsm6.F90 lines 1472-1495
108  // Weierstrass infinite product form for Gamma(x).
109  // Special case x==1 returns 0.0 matching Fortran exactly.
110  // Never triggered in practice (bvtr1=1.8, bvts1=1.41, etc.)
111  // but preserved for bit-compatible validation against Fortran reference.
112  // CPU-only: called from initialize_coeffs(), not from GPU kernels.
113  auto rgmma = [](Real x) -> Real {
114  if (x == Real(1.0)) return Real(0.0);
115  constexpr Real euler = Real(0.577215664901532);
116  Real result = x * std::exp(euler * x);
117  for (int i = 1; i <= 10000; ++i) {
118  Real y = Real(i);
119  result = result * (Real(1.0) + x/y) * std::exp(-x/y);
120  }
121  return Real(1.0) / result;
122  };
123 
124  // Physical constants matching mp_wsm6_init argument list
125  // den0: reference air density at 850mb (kg/m^3)
126  // denr: liquid water density — rhoh2o from ERF_Constants.H
127  // dens_arg: snow density — dens_snow constexpr = 100.0
128  // cl: specific heat liquid water — Cp_l from ERF_Constants.H
129  // cpv: specific heat water vapor — Cp_v from ERF_Constants.H
130  const Real den0 = Real(1.28);
131  const Real denr = Real(rhoh2o);
132  const Real dens_arg = dens_snow;
133  const Real cl = Real(Cp_l);
134  const Real cpv_loc = Real(Cp_v);
135 
136  // hail_opt branch: 5 regime-dependent coefficients
137  if (m_hail_opt) {
138  m_n0g = Real(4.0e4);
139  m_deng = Real(700.0);
140  m_avtg = Real(285.0);
141  m_bvtg = Real(0.8);
142  m_lamdagmax = Real(2.0e4);
143  } else {
144  m_n0g = Real(4.0e6);
145  m_deng = Real(500.0);
146  m_avtg = Real(330.0);
147  m_bvtg = Real(0.8);
148  m_lamdagmax = Real(6.0e4);
149  }
150 
151  m_pi_wsm6 = Real(4.0) * std::atan(Real(1.0));
152  m_xlv1 = cl - cpv_loc;
153 
154  m_qc0 = Real(4.0)/Real(3.0) * m_pi_wsm6 * denr
155  * std::pow(r0, Real(3.0)) * xncr / den0;
156  m_qck1 = Real(0.104) * Real(9.8) * peaut
157  / std::pow(xncr * denr, Real(1.0)/Real(3.0))
158  / xmyu * std::pow(den0, Real(4.0)/Real(3.0));
159  m_pidnc = m_pi_wsm6 * denr / Real(6.0);
160 
161  // Rain coefficients
162  m_bvtr1 = Real(1.0) + bvtr;
163  m_bvtr2 = Real(2.5) + Real(0.5) * bvtr;
164  m_bvtr3 = Real(3.0) + bvtr;
165  m_bvtr4 = Real(4.0) + bvtr;
166  m_bvtr6 = Real(6.0) + bvtr;
167  m_g1pbr = rgmma(m_bvtr1);
168  m_g3pbr = rgmma(m_bvtr3);
169  m_g4pbr = rgmma(m_bvtr4);
170  m_g6pbr = rgmma(m_bvtr6);
172  m_pvtr = avtr * m_g4pbr / Real(6.0);
173  m_eacrr = Real(1.0);
174  m_pacrr = m_pi_wsm6 * n0r * avtr * m_g3pbr * Real(0.25) * m_eacrr;
175  m_precr1 = Real(2.0) * m_pi_wsm6 * n0r * Real(0.78);
176  m_precr2 = Real(2.0) * m_pi_wsm6 * n0r * Real(0.31)
177  * std::pow(avtr, Real(0.5)) * m_g5pbro2;
178  m_roqimax = Real(2.08e22) * std::pow(dimax, Real(8.0));
179 
180  // Snow coefficients
181  m_bvts1 = Real(1.0) + bvts;
182  m_bvts2 = Real(2.5) + Real(0.5) * bvts;
183  m_bvts3 = Real(3.0) + bvts;
184  m_bvts4 = Real(4.0) + bvts;
185  m_g1pbs = rgmma(m_bvts1);
186  m_g3pbs = rgmma(m_bvts3);
187  m_g4pbs = rgmma(m_bvts4);
189  m_pvts = avts * m_g4pbs / Real(6.0);
190  m_pacrs = m_pi_wsm6 * n0s * avts * m_g3pbs * Real(0.25);
191  m_precs1 = Real(4.0) * n0s * Real(0.65);
192  m_precs2 = Real(4.0) * n0s * Real(0.44)
193  * std::pow(avts, Real(0.5)) * m_g5pbso2;
194  m_pidn0r = m_pi_wsm6 * denr * n0r;
195  m_pidn0s = m_pi_wsm6 * dens_arg * n0s;
196  m_pacrc = m_pi_wsm6 * n0s * avts * m_g3pbs * Real(0.25) * eacrc;
197 
198  // Graupel/hail coefficients
199  m_bvtg1 = Real(1.0) + m_bvtg;
200  m_bvtg2 = Real(2.5) + Real(0.5) * m_bvtg;
201  m_bvtg3 = Real(3.0) + m_bvtg;
202  m_bvtg4 = Real(4.0) + m_bvtg;
203  m_g1pbg = rgmma(m_bvtg1);
204  m_g3pbg = rgmma(m_bvtg3);
205  m_g4pbg = rgmma(m_bvtg4);
206  m_pacrg = m_pi_wsm6 * m_n0g * m_avtg * m_g3pbg * Real(0.25);
208  m_pvtg = m_avtg * m_g4pbg / Real(6.0);
209  m_precg1 = Real(2.0) * m_pi_wsm6 * m_n0g * Real(0.78);
210  m_precg2 = Real(2.0) * m_pi_wsm6 * m_n0g * Real(0.31)
211  * std::pow(m_avtg, Real(0.5)) * m_g5pbgo2;
213 
214  // Slope parameter limits
215  m_rslopermax = Real(1.0) / lamdarmax;
216  m_rslopesmax = Real(1.0) / lamdasmax;
217  m_rslopegmax = Real(1.0) / m_lamdagmax;
218  m_rsloperbmax = std::pow(m_rslopermax, bvtr);
219  m_rslopesbmax = std::pow(m_rslopesmax, bvts);
220  m_rslopegbmax = std::pow(m_rslopegmax, m_bvtg);
227 }
amrex::Real m_g5pbgo2
Definition: ERF_WSM6.H:282
amrex::Real m_bvts3
Definition: ERF_WSM6.H:277
bool m_hail_opt
Definition: ERF_WSM6.H:265
amrex::Real m_bvtr1
Definition: ERF_WSM6.H:273
amrex::Real m_g4pbs
Definition: ERF_WSM6.H:278
amrex::Real m_g1pbg
Definition: ERF_WSM6.H:282
amrex::Real m_g4pbg
Definition: ERF_WSM6.H:282
amrex::Real m_g5pbro2
Definition: ERF_WSM6.H:274
amrex::Real m_g3pbr
Definition: ERF_WSM6.H:274
amrex::Real m_bvtg1
Definition: ERF_WSM6.H:281
amrex::Real m_bvts4
Definition: ERF_WSM6.H:277
amrex::Real m_deng
Definition: ERF_WSM6.H:266
amrex::Real m_g3pbg
Definition: ERF_WSM6.H:282
static constexpr amrex::Real lamdasmax
Definition: ERF_WSM6.H:156
amrex::Real m_bvtr6
Definition: ERF_WSM6.H:273
amrex::Real m_bvtr2
Definition: ERF_WSM6.H:273
amrex::Real m_bvts1
Definition: ERF_WSM6.H:277
amrex::Real m_pacrs
Definition: ERF_WSM6.H:279
amrex::Real m_eacrr
Definition: ERF_WSM6.H:275
amrex::Real m_g1pbs
Definition: ERF_WSM6.H:278
amrex::Real m_lamdagmax
Definition: ERF_WSM6.H:266
amrex::Real m_bvtg2
Definition: ERF_WSM6.H:281
static constexpr amrex::Real r0
Definition: ERF_WSM6.H:149
amrex::Real m_bvtg4
Definition: ERF_WSM6.H:281
amrex::Real m_bvtr4
Definition: ERF_WSM6.H:273
amrex::Real m_bvts2
Definition: ERF_WSM6.H:277
amrex::Real m_bvtg3
Definition: ERF_WSM6.H:281
amrex::Real m_pidnc
Definition: ERF_WSM6.H:272
amrex::Real m_g3pbs
Definition: ERF_WSM6.H:278
static constexpr amrex::Real lamdarmax
Definition: ERF_WSM6.H:155
amrex::Real m_avtg
Definition: ERF_WSM6.H:266
amrex::Real m_g4pbr
Definition: ERF_WSM6.H:274
amrex::Real m_g1pbr
Definition: ERF_WSM6.H:274
static constexpr amrex::Real peaut
Definition: ERF_WSM6.H:150
static constexpr amrex::Real avts
Definition: ERF_WSM6.H:153
static constexpr amrex::Real eacrc
Definition: ERF_WSM6.H:162
amrex::Real m_g5pbso2
Definition: ERF_WSM6.H:278
static constexpr amrex::Real xmyu
Definition: ERF_WSM6.H:152
amrex::Real m_bvtr3
Definition: ERF_WSM6.H:273
real(kind=kind_phys) function rgmma(x)
Definition: ERF_module_mp_wdm6.F90:2185
Here is the call graph for this function:

◆ Qmoist_Ptr()

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

Reimplemented from NullMoist.

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

◆ Qmoist_Restart_Vars()

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

Reimplemented from NullMoist.

215  {
216  a_idx = {0, 1, 2, 3};
217  a_names = {"RainAccum", "SnowAccum", "GraupAccum", "PrecipAccum"};
218  }

◆ Qmoist_Size()

int WSM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

208 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WSM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

210 { return n_qstate_moist_numconc_size; }
int n_qstate_moist_numconc_size
Definition: ERF_WSM6.H:239

◆ Qstate_Moist_Size()

int WSM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

209 { return n_qstate_moist_size; }
int n_qstate_moist_size
Definition: ERF_WSM6.H:238

◆ Sed_Flux()

const SedimentationFlux* WSM6::Sed_Flux ( ) const
inlineoverridevirtual

The sedimentation face fluxes of the most recent Advance, or null for a scheme that does not sediment (see ERF_SedimentationFlux.H)

Reimplemented from NullMoist.

200 { return &m_sed; }

◆ Set_dzmin()

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

Reimplemented from NullMoist.

176 { m_dzmin = dz_min; }
amrex::Real m_dzmin
Definition: ERF_WSM6.H:245

◆ Update_Micro_Vars() [1/2]

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

Reimplemented from NullMoist.

187  {
188  Copy_State_to_Micro(cons_in);
189  }
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Here is the call graph for this function:

◆ Update_Micro_Vars() [2/2]

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

Reimplemented from NullMoist.

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

193  {
194  Copy_Micro_to_State(cons_in);
195  }
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateWSM6.cpp:7
Here is the call graph for this function:

Member Data Documentation

◆ alpha_wsm6

constexpr amrex::Real WSM6::alpha_wsm6 = amrex::Real(0.12)
staticconstexpr

◆ avtr

constexpr amrex::Real WSM6::avtr = amrex::Real(841.9)
staticconstexpr

◆ avts

constexpr amrex::Real WSM6::avts = amrex::Real(11.72)
staticconstexpr

◆ bvtr

constexpr amrex::Real WSM6::bvtr = amrex::Real(0.8)
staticconstexpr

◆ bvts

constexpr amrex::Real WSM6::bvts = amrex::Real(0.41)
staticconstexpr

◆ dens_snow

constexpr amrex::Real WSM6::dens_snow = amrex::Real(100.0)
staticconstexpr

◆ dicon

constexpr amrex::Real WSM6::dicon = amrex::Real(11.9)
staticconstexpr

◆ dimax

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

◆ dt

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

◆ dtcldcr

constexpr amrex::Real WSM6::dtcldcr = amrex::Real(120.0)
staticconstexpr

◆ eacrc

constexpr amrex::Real WSM6::eacrc = amrex::Real(1.0)
staticconstexpr

◆ lamdarmax

constexpr amrex::Real WSM6::lamdarmax = amrex::Real(8.0e4)
staticconstexpr

◆ lamdasmax

constexpr amrex::Real WSM6::lamdasmax = amrex::Real(1.0e5)
staticconstexpr

◆ m_avtg

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

◆ m_axis

int WSM6::m_axis {2}
private

Referenced by Define().

◆ m_bvtg

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

◆ m_bvtg1

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

◆ m_bvtg2

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

◆ m_bvtg3

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

◆ m_bvtg4

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

◆ m_bvtr1

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

◆ m_bvtr2

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

◆ m_bvtr3

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

◆ m_bvtr4

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

◆ m_bvtr6

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

◆ m_bvts1

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

◆ m_bvts2

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

◆ m_bvts3

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

◆ m_bvts4

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

◆ m_columns

MicrophysicsColumns WSM6::m_columns
private

◆ m_deng

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

◆ m_detJ_cc

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

◆ m_do_cond

bool WSM6::m_do_cond {true}
private

Referenced by Define().

◆ m_dzmin

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

Referenced by Set_dzmin().

◆ m_eacrr

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

◆ m_g1pbg

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

◆ m_g1pbr

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

◆ m_g1pbs

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

◆ m_g3pbg

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

◆ m_g3pbr

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

◆ m_g3pbs

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

◆ m_g4pbg

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

◆ m_g4pbr

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

◆ m_g4pbs

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

◆ m_g5pbgo2

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

◆ m_g5pbro2

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

◆ m_g5pbso2

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

◆ m_g6pbr

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

◆ m_geom

amrex::Geometry WSM6::m_geom
private

◆ m_hail_opt

bool WSM6::m_hail_opt {false}
private

◆ m_lamdagmax

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

◆ m_moisture_type

MoistureType WSM6::m_moisture_type {MoistureType::None}
private

Referenced by Define().

◆ m_n0g

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

◆ m_pacrc

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

◆ m_pacrg

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

◆ m_pacrr

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

◆ m_pacrs

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

◆ m_pi_wsm6

amrex::Real WSM6::m_pi_wsm6 {0}
private

◆ m_pidn0g

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

◆ m_pidn0r

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

◆ m_pidn0s

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

◆ m_pidnc

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

◆ m_precg1

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

◆ m_precg2

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

◆ m_precr1

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

◆ m_precr2

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

◆ m_precs1

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

◆ m_precs2

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

◆ m_pvtg

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

◆ m_pvtr

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

◆ m_pvts

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

◆ m_qc0

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

◆ m_qck1

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

◆ m_qmoist_size

int WSM6::m_qmoist_size = 4
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

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

Referenced by Define().

◆ m_roqimax

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

◆ m_rslopeg2max

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

◆ m_rslopeg3max

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

◆ m_rslopegbmax

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

◆ m_rslopegmax

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

◆ m_rsloper2max

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

◆ m_rsloper3max

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

◆ m_rsloperbmax

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

◆ m_rslopermax

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

◆ m_rslopes2max

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

◆ m_rslopes3max

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

◆ m_rslopesbmax

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

◆ m_rslopesmax

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

◆ m_sed

SedimentationFlux WSM6::m_sed
private

Referenced by Sed_Flux().

◆ m_xlv1

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

◆ m_z_phys_nd

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

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar_WSM6::NumVars> WSM6::mic_fab_vars
private

◆ MicVarMap

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

Referenced by Qmoist_Ptr().

◆ n0r

constexpr amrex::Real WSM6::n0r = amrex::Real(8.0e6)
staticconstexpr

◆ n0s

constexpr amrex::Real WSM6::n0s = amrex::Real(2.0e6)
staticconstexpr

◆ n0smax

constexpr amrex::Real WSM6::n0smax = amrex::Real(1.0e11)
staticconstexpr

◆ n_qstate_moist_numconc_size

int WSM6::n_qstate_moist_numconc_size = 0
private

◆ n_qstate_moist_size

int WSM6::n_qstate_moist_size = 6
private

Referenced by Qstate_Moist_Size().

◆ nlev

int WSM6::nlev {0}
private

◆ peaut

constexpr amrex::Real WSM6::peaut = amrex::Real(0.55)
staticconstexpr

◆ pfrz1

constexpr amrex::Real WSM6::pfrz1 = amrex::Real(100.0)
staticconstexpr

◆ pfrz2

constexpr amrex::Real WSM6::pfrz2 = amrex::Real(0.66)
staticconstexpr

◆ qcrmin

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

◆ qs0

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

◆ r0

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

◆ xmyu

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

◆ xncr

constexpr amrex::Real WSM6::xncr = amrex::Real(3.0e8)
staticconstexpr

Referenced by AdvanceColumns().

◆ zhi

int WSM6::zhi {0}
private

◆ zlo

int WSM6::zlo {0}
private

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