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
 
SurfacePrecipAccumulationSources Get_Surface_Precip_Accumulation_Ptrs (const int &) const override
 
virtual void Update_Micro_Vars (amrex::MultiFab &)
 
virtual void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *)
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *)
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 
virtual void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const
 
virtual void SetCurrentLevel (const int &)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_RealWidth (const int)
 

Static Public Attributes

static constexpr amrex::Real dtcldcr = 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}
 
MoistureType m_moisture_type {MoistureType::None}
 
amrex::MultiFab * m_z_phys_nd {nullptr}
 
amrex::MultiFab * m_detJ_cc {nullptr}
 
amrex::Array< FabPtr, MicVar_WSM6::NumVars > mic_fab_vars
 
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
42 {}

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

48 {
49  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
50  // Match Morrison behavior: refresh microphysics ghost zones from state.
51  // WSM6 Fortran reads the full (ims:ime, jms:jme, kms:kme) slab.
52  const auto& box3d = mfi.growntilebox();
53  auto states = cons_in.array(mfi);
54 
55  auto rho = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
56  auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi);
57  auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
58  auto pres = mic_fab_vars[MicVar_WSM6::pres]->array(mfi);
59 
60  auto qv = mic_fab_vars[MicVar_WSM6::qv]->array(mfi);
61  auto qc = mic_fab_vars[MicVar_WSM6::qc]->array(mfi);
62  auto qi = mic_fab_vars[MicVar_WSM6::qi]->array(mfi);
63  auto qr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi);
64  auto qs = mic_fab_vars[MicVar_WSM6::qs]->array(mfi);
65  auto qg = mic_fab_vars[MicVar_WSM6::qg]->array(mfi);
66 
67  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
68  rho(i,j,k) = states(i,j,k,Rho_comp);
69  theta(i,j,k) = states(i,j,k,RhoTheta_comp) / states(i,j,k,Rho_comp);
70 
71  qv(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ1_comp) / states(i,j,k,Rho_comp));
72  qc(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ2_comp) / states(i,j,k,Rho_comp));
73  qi(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ3_comp) / states(i,j,k,Rho_comp));
74  qr(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ4_comp) / states(i,j,k,Rho_comp));
75  qs(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ5_comp) / states(i,j,k,Rho_comp));
76  qg(i,j,k) = amrex::max(Real(0.0), states(i,j,k,RhoQ6_comp) / states(i,j,k,Rho_comp));
77 
78  tabs(i,j,k) = getTgivenRandRTh(states(i,j,k,Rho_comp),
79  states(i,j,k,RhoTheta_comp),
80  qv(i,j,k));
81  pres(i,j,k) = getPgivenRTh(states(i,j,k,RhoTheta_comp), qv(i,j,k));
82  });
83  }
84 }
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:27

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.

46  {
47  m_moisture_type = sc.moisture_type;
48  m_axis = sc.ave_plane;
49  m_do_cond = (!sc.uses_shoc_family());
50  }
WSM6::m_axis
int m_axis
Definition: ERF_WSM6.H:144
WSM6::m_moisture_type
MoistureType m_moisture_type
Definition: ERF_WSM6.H:146
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:1604
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:1618
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:1566
Here is the call graph for this function:

◆ Get_Surface_Precip_Accumulation_Ptrs()

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

Reimplemented from NullMoist.

125  {
126  SurfacePrecipAccumulationSources sources;
127  sources.total = {mic_fab_vars[MicVar_WSM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
128  sources.snow = {mic_fab_vars[MicVar_WSM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
129  sources.graupel = {mic_fab_vars[MicVar_WSM6::graup_accum].get(), rhoh2o / amrex::Real(1000.0)};
130  return sources;
131  }
SurfacePrecipAccumulationSources
Definition: ERF_SurfacePrecipitation.H:21
SurfacePrecipAccumulationSources::snow
SurfacePrecipAccumulationSource snow
Definition: ERF_SurfacePrecipitation.H:24
SurfacePrecipAccumulationSources::total
SurfacePrecipAccumulationSource total
Definition: ERF_SurfacePrecipitation.H:22
SurfacePrecipAccumulationSources::graupel
SurfacePrecipAccumulationSource graupel
Definition: ERF_SurfacePrecipitation.H:25

◆ 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 #if defined(ERF_USE_WSM6_FORT) && defined(AMREX_USE_GPU)
26  // The host-only Fortran bridge receives dataPtr() from mic_fab_vars.
27  Arena* Arena_Used = The_Managed_Arena();
28 #else
29  Arena* Arena_Used = The_Arena();
30 #endif
31 
32  for (int ivar = 0; ivar < MicVar_WSM6::NumVars; ++ivar) {
33  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
34  1, cons_in.nGrowVect(),
35  MFInfo().SetArena(Arena_Used));
36  mic_fab_vars[ivar]->setVal(0.0);
37  }
38 
39  nlev = m_geom.Domain().length(2);
40  zlo = m_geom.Domain().smallEnd(2);
41  zhi = m_geom.Domain().bigEnd(2);
42 
43  initialize_coeffs();
44 }
WSM6::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_WSM6.H:149
WSM6::zlo
int zlo
Definition: ERF_WSM6.H:143
WSM6::nlev
int nlev
Definition: ERF_WSM6.H:143
WSM6::zhi
int zhi
Definition: ERF_WSM6.H:143
WSM6::m_qmoist_size
int m_qmoist_size
Definition: ERF_WSM6.H:134
WSM6::initialize_coeffs
void initialize_coeffs()
Definition: ERF_InitWSM6.cpp:87
WSM6::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_WSM6.H:138
MicVar_WSM6::NumVars
@ NumVars
Definition: ERF_WSM6.H:34

◆ initialize_coeffs()

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

106  {
107  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
108  return mic_fab_vars[MicVarMap[varIdx]].get();
109  }
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.

118  {
119  a_idx = {0, 1, 2};
120  a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
121  }

◆ Qmoist_Size()

int WSM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

111 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WSM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

113 { return n_qstate_moist_numconc_size; }
WSM6::n_qstate_moist_numconc_size
int n_qstate_moist_numconc_size
Definition: ERF_WSM6.H:136

◆ Qstate_Moist_Size()

int WSM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

112 { return n_qstate_moist_size; }
WSM6::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_WSM6.H:135

◆ Set_dzmin()

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

Reimplemented from NullMoist.

84 { m_dzmin = dz_min; }
WSM6::m_dzmin
amrex::Real m_dzmin
Definition: ERF_WSM6.H:142

◆ Update_Micro_Vars() [1/3]

virtual void NullMoist::Update_Micro_Vars
inline
35 { }

◆ Update_Micro_Vars() [2/3]

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

Reimplemented from NullMoist.

92  {
93  Copy_State_to_Micro(cons_in);
94  }
WSM6::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitWSM6.cpp:47
Here is the call graph for this function:

◆ Update_Micro_Vars() [3/3]

virtual void NullMoist::Update_Micro_Vars
inline
42  {
43  Update_Micro_Vars(cons_in);
44  }
WSM6::Update_Micro_Vars
virtual void Update_Micro_Vars(amrex::MultiFab &)
Definition: ERF_NullMoist.H:35

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

98  {
99  Copy_Micro_to_State(cons_in);
100  }
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_moisture_type

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

Referenced by Define().

◆ m_n0g

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

◆ m_pacrc

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

◆ m_pacrg

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

◆ m_pacrr

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

◆ m_pacrs

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

◆ m_pi_wsm6

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

◆ m_pidn0g

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

◆ m_pidn0r

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

◆ m_pidn0s

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

◆ m_pidnc

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

◆ m_precg1

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

◆ m_precg2

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

◆ m_precr1

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

◆ m_precr2

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

◆ m_precs1

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

◆ m_precs2

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

◆ m_pvtg

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

◆ m_pvtr

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

◆ m_pvts

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

◆ m_qc0

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

◆ m_qck1

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

◆ m_qmoist_size

int WSM6::m_qmoist_size = 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 Get_Surface_Precip_Accumulation_Ptrs(), and 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: