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

#include <ERF_WSM6.H>

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

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_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &solverChoice) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
int Qstate_Moist_NumConc_Size () override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
- 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 &)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_RealWidth (const int)
 

Static Public Attributes

static constexpr amrex::Real dtcldcr = 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 = 3
 
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::MultiFab * m_z_phys_nd {nullptr}
 
amrex::MultiFab * m_detJ_cc {nullptr}
 
amrex::Array< FabPtr, MicVar_WSM6::NumVars > mic_fab_vars
 
bool m_hail_opt {false}
 
amrex::Real m_n0g {0}
 
amrex::Real m_deng {0}
 
amrex::Real m_avtg {0}
 
amrex::Real m_bvtg {0}
 
amrex::Real m_lamdagmax {0}
 
amrex::Real m_pi_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}
 

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ WSM6()

WSM6::WSM6 ( )
inline
41 {}

◆ ~WSM6()

virtual WSM6::~WSM6 ( )
virtualdefault

Member Function Documentation

◆ Advance()

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

Reimplemented from NullMoist.

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

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro()

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

Reimplemented from NullMoist.

40 {
41  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
42  // Match Morrison behavior: refresh microphysics ghost zones from state.
43  // WSM6 Fortran reads the full (ims:ime, jms:jme, kms:kme) slab.
44  const auto& box3d = mfi.growntilebox();
45  auto states = cons_in.array(mfi);
46 
47  auto rho = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
48  auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi);
49  auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
50  auto pres = mic_fab_vars[MicVar_WSM6::pres]->array(mfi);
51 
52  auto qv = mic_fab_vars[MicVar_WSM6::qv]->array(mfi);
53  auto qc = mic_fab_vars[MicVar_WSM6::qc]->array(mfi);
54  auto qi = mic_fab_vars[MicVar_WSM6::qi]->array(mfi);
55  auto qr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi);
56  auto qs = mic_fab_vars[MicVar_WSM6::qs]->array(mfi);
57  auto qg = mic_fab_vars[MicVar_WSM6::qg]->array(mfi);
58 
59  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
60  rho(i,j,k) = states(i,j,k,Rho_comp);
61  theta(i,j,k) = states(i,j,k,RhoTheta_comp) / states(i,j,k,Rho_comp);
62 
63  qv(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ1_comp) / states(i,j,k,Rho_comp));
64  qc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ2_comp) / states(i,j,k,Rho_comp));
65  qi(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ3_comp) / states(i,j,k,Rho_comp));
66  qr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ4_comp) / states(i,j,k,Rho_comp));
67  qs(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ5_comp) / states(i,j,k,Rho_comp));
68  qg(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ6_comp) / states(i,j,k,Rho_comp));
69 
70  tabs(i,j,k) = getTgivenRandRTh(states(i,j,k,Rho_comp),
71  states(i,j,k,RhoTheta_comp),
72  qv(i,j,k));
73  pres(i,j,k) = getPgivenRTh(states(i,j,k,RhoTheta_comp), qv(i,j,k));
74  });
75  }
76 }
getTgivenRandRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:46
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:26

Referenced by Update_Micro_Vars().

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

◆ Define()

void WSM6::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

45  {
46  m_axis = sc.ave_plane;
47  m_do_cond = (!sc.uses_shoc_family());
48  }
WSM6::m_axis
int m_axis
Definition: ERF_WSM6.H:130
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Definition: ERF_DataStruct.H:1359
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:1403
Here is the call graph for this function:

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

15 {
16  dt = dt_advance;
17  m_geom = geom;
18 
19  m_z_phys_nd = z_phys_nd.get();
20  m_detJ_cc = detJ_cc.get();
21 
22  MicVarMap.resize(m_qmoist_size);
23  MicVarMap = {MicVar_WSM6::rain_accum, MicVar_WSM6::snow_accum, MicVar_WSM6::graup_accum};
24 
25  for (int ivar = 0; ivar < MicVar_WSM6::NumVars; ++ivar) {
26  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
27  1, cons_in.nGrowVect());
28  mic_fab_vars[ivar]->setVal(0.0);
29  }
30 
31  nlev = m_geom.Domain().length(2);
32  zlo = m_geom.Domain().smallEnd(2);
33  zhi = m_geom.Domain().bigEnd(2);
34 
35  initialize_coeffs();
36 }
WSM6::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_WSM6.H:134
WSM6::zlo
int zlo
Definition: ERF_WSM6.H:129
WSM6::nlev
int nlev
Definition: ERF_WSM6.H:129
WSM6::zhi
int zhi
Definition: ERF_WSM6.H:129
WSM6::m_qmoist_size
int m_qmoist_size
Definition: ERF_WSM6.H:120
WSM6::initialize_coeffs
void initialize_coeffs()
Definition: ERF_InitWSM6.cpp:79
WSM6::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_WSM6.H:124
MicVar_WSM6::NumVars
@ NumVars
Definition: ERF_WSM6.H:33

◆ initialize_coeffs()

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

◆ Qmoist_Ptr()

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

Reimplemented from NullMoist.

102  {
103  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
104  return mic_fab_vars[MicVarMap[varIdx]].get();
105  }
AMREX_ALWAYS_ASSERT
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.

114  {
115  a_idx = {0, 1, 2};
116  a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
117  }

◆ Qmoist_Size()

int WSM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

107 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WSM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

109 { return n_qstate_moist_numconc_size; }
WSM6::n_qstate_moist_numconc_size
int n_qstate_moist_numconc_size
Definition: ERF_WSM6.H:122

◆ Qstate_Moist_Size()

int WSM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

108 { return n_qstate_moist_size; }
WSM6::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_WSM6.H:121

◆ Set_dzmin()

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

Reimplemented from NullMoist.

82 { m_dzmin = dz_min; }
WSM6::m_dzmin
amrex::Real m_dzmin
Definition: ERF_WSM6.H:128

◆ Update_Micro_Vars()

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

Reimplemented from NullMoist.

88  {
89  Copy_State_to_Micro(cons_in);
90  }
WSM6::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitWSM6.cpp:39
Here is the call graph for this function:

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

94  {
95  Copy_Micro_to_State(cons_in);
96  }
WSM6::Copy_Micro_to_State
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

Referenced by wsm6_nislfv_rain_plm_scratch().

◆ avts

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

Referenced by wsm6_nislfv_rain_plm6_scratch().

◆ bvtr

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

Referenced by wsm6_nislfv_rain_plm_scratch().

◆ bvts

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

Referenced by wsm6_nislfv_rain_plm6_scratch().

◆ dens_snow

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

Referenced by wsm6_nislfv_rain_plm6_scratch().

◆ 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

Referenced by wsm6_nislfv_rain_plm_scratch().

◆ lamdasmax

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

Referenced by wsm6_nislfv_rain_plm6_scratch().

◆ 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_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_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 = 3
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

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

Referenced by Qmoist_Ptr().

◆ MicVarMap

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

Referenced by Qmoist_Ptr().

◆ n0r

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

Referenced by Advance(), and wsm6_nislfv_rain_plm_scratch().

◆ n0s

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

Referenced by wsm6_nislfv_rain_plm6_scratch().

◆ 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

Referenced by Qstate_Moist_NumConc_Size().

◆ 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

Referenced by wsm6_nislfv_rain_plm6_scratch(), and wsm6_nislfv_rain_plm_scratch().

◆ 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 Advance().

◆ zhi

int WSM6::zhi {0}
private

◆ zlo

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