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_Lmask (amrex::iMultiFab *)
 Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this. More...
 
virtual void Set_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.queryAdd("microphysics_debug", microphysics_debug);
856  pp.queryAdd("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.queryAdd("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  // Recompute qsum from the current qs/qg (they have changed since
1864  // qsum_arr was set in G5a) -- mirrors ERF_module_mp_wsm6.F90 l.1008
1865  const Real qsum = amrex::max(qs_arr(i,j,k) + qg_arr(i,j,k), Real(1.0e-15));
1866  qsum_arr(i,j,k) = qsum;
1867  const Real vt2ave = (qsum > Real(1.0e-15))
1868  ? (vt2s * qs_arr(i,j,k) + vt2g * qg_arr(i,j,k)) / qsum
1869  : Real(0.0);
1870 
1871  if (supcol > Real(0.0) && qi_arr(i,j,k) > Real(qmin)) {
1872  if (qr_arr(i,j,k) > Real(qcrmin)) {
1873  const Real acrfac =
1874  Real(2.0) * rslope3_r_arr(i,j,k)
1875  + Real(2.0) * diameter * rslope2_r_arr(i,j,k)
1876  + diameter * diameter * rslope_r_arr(i,j,k);
1877 
1878  // WSM6-CPP TAG: PRACI
1879  // legacy_group: G13b
1880  // process: Accretion of cloud ice by rain
1881  // compare_vars: praci, qi, qr, den
1882  praci_arr(i,j,k) = Real(pi) * qi_arr(i,j,k) * Real(n0r)
1883  * std::abs(vt2r - vt2i) * acrfac / Real(4.0);
1884  praci_arr(i,j,k) *= std::pow(
1885  amrex::min(
1886  amrex::max(Real(0.0), qr_arr(i,j,k) / qi_arr(i,j,k)),
1887  Real(1.0)),
1888  Real(2.0));
1889  praci_arr(i,j,k) = amrex::min(
1890  praci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
1891 
1892  piacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(avtr)
1893  * Real(n0r) * Real(denr) * xni * denfac_arr(i,j,k)
1894  * Real(g6pbr) * rslope3_r_arr(i,j,k)
1895  * rslope3_r_arr(i,j,k) * rslopeb_r_arr(i,j,k)
1896  / Real(24.0) / den_arr(i,j,k);
1897  piacr_arr(i,j,k) *= std::pow(
1898  amrex::min(
1899  amrex::max(Real(0.0), qi_arr(i,j,k) / qr_arr(i,j,k)),
1900  Real(1.0)),
1901  Real(2.0));
1902  piacr_arr(i,j,k) = amrex::min(
1903  piacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
1904  }
1905 
1906  if (qs_arr(i,j,k) > Real(qcrmin)) {
1907  const Real acrfac =
1908  Real(2.0) * rslope3_s_arr(i,j,k)
1909  + Real(2.0) * diameter * rslope2_s_arr(i,j,k)
1910  + diameter * diameter * rslope_s_arr(i,j,k);
1911  // WSM6-CPP TAG: PSACI
1912  // legacy_group: G13e
1913  // process: Accretion of cloud ice by snow
1914  // compare_vars: psaci, qi, qs, den
1915  psaci_arr(i,j,k) = Real(pi) * qi_arr(i,j,k) * eacrs
1916  * Real(n0s) * n0sfac_arr(i,j,k)
1917  * std::abs(vt2ave - vt2i) * acrfac / Real(4.0);
1918  psaci_arr(i,j,k) = amrex::min(
1919  psaci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
1920  }
1921 
1922  if (qg_arr(i,j,k) > Real(qcrmin)) {
1923  const Real egi = std::exp(Real(0.07) * (-supcol));
1924  const Real acrfac =
1925  Real(2.0) * rslope3_g_arr(i,j,k)
1926  + Real(2.0) * diameter * rslope2_g_arr(i,j,k)
1927  + diameter * diameter * rslope_g_arr(i,j,k);
1928  pgaci_arr(i,j,k) = Real(pi) * egi * qi_arr(i,j,k)
1929  * n0g_loc * std::abs(vt2ave - vt2i) * acrfac
1930  / Real(4.0);
1931  pgaci_arr(i,j,k) = amrex::min(
1932  pgaci_arr(i,j,k), qi_arr(i,j,k) / dtcld);
1933  }
1934  }
1935 
1936  if (qs_arr(i,j,k) > Real(qcrmin) && qc_arr(i,j,k) > Real(qmin)) {
1937  psacw_arr(i,j,k) = amrex::min(
1938  Real(pacrc) * n0sfac_arr(i,j,k) * rslope3_s_arr(i,j,k)
1939  * rslopeb_s_arr(i,j,k)
1940  * std::pow(
1941  amrex::min(
1942  amrex::max(Real(0.0), qs_arr(i,j,k) / qc_arr(i,j,k)),
1943  Real(1.0)),
1944  Real(2.0))
1945  * qc_arr(i,j,k) * denfac_arr(i,j,k),
1946  qc_arr(i,j,k) / dtcld);
1947  }
1948 
1949  if (qg_arr(i,j,k) > Real(qcrmin) && qc_arr(i,j,k) > Real(qmin)) {
1950  pgacw_arr(i,j,k) = amrex::min(
1951  Real(pacrg) * rslope3_g_arr(i,j,k) * rslopeb_g_arr(i,j,k)
1952  * std::pow(
1953  amrex::min(
1954  amrex::max(Real(0.0), qg_arr(i,j,k) / qc_arr(i,j,k)),
1955  Real(1.0)),
1956  Real(2.0))
1957  * qc_arr(i,j,k) * denfac_arr(i,j,k),
1958  qc_arr(i,j,k) / dtcld);
1959  }
1960 
1961  if (qsum > Real(1.0e-15)) {
1962  paacw_arr(i,j,k) = (qs_arr(i,j,k) * psacw_arr(i,j,k)
1963  + qg_arr(i,j,k) * pgacw_arr(i,j,k))
1964  / qsum;
1965  }
1966 
1967  if (qs_arr(i,j,k) > Real(qcrmin) && qr_arr(i,j,k) > Real(qcrmin)) {
1968  if (supcol > Real(0.0)) {
1969  const Real acrfac =
1970  Real(5.0) * rslope3_s_arr(i,j,k) * rslope3_s_arr(i,j,k)
1971  * rslope_r_arr(i,j,k)
1972  + Real(2.0) * rslope3_s_arr(i,j,k) * rslope2_s_arr(i,j,k)
1973  * rslope2_r_arr(i,j,k)
1974  + Real(0.5) * rslope2_s_arr(i,j,k) * rslope2_s_arr(i,j,k)
1975  * rslope3_r_arr(i,j,k);
1976  // WSM6-CPP TAG: PRACS
1977  // legacy_group: G13f
1978  // process: Accretion of snow by rain / rain-snow interaction
1979  // compare_vars: pracs, qr, qs, den
1980  pracs_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
1981  * Real(n0s) * n0sfac_arr(i,j,k)
1982  * std::abs(vt2r - vt2ave) * (Real(dens_snow) / den_arr(i,j,k))
1983  * acrfac;
1984  pracs_arr(i,j,k) *= std::pow(
1985  amrex::min(
1986  amrex::max(Real(0.0), qr_arr(i,j,k) / qs_arr(i,j,k)),
1987  Real(1.0)),
1988  Real(2.0));
1989  pracs_arr(i,j,k) = amrex::min(
1990  pracs_arr(i,j,k), qs_arr(i,j,k) / dtcld);
1991  }
1992 
1993  {
1994  const Real acrfac =
1995  Real(5.0) * rslope3_r_arr(i,j,k) * rslope3_r_arr(i,j,k)
1996  * rslope_s_arr(i,j,k)
1997  + Real(2.0) * rslope3_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
1998  * rslope2_s_arr(i,j,k)
1999  + Real(0.5) * rslope2_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2000  * rslope3_s_arr(i,j,k);
2001  psacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
2002  * Real(n0s) * n0sfac_arr(i,j,k)
2003  * std::abs(vt2ave - vt2r) * (Real(denr) / den_arr(i,j,k))
2004  * acrfac;
2005  psacr_arr(i,j,k) *= std::pow(
2006  amrex::min(
2007  amrex::max(Real(0.0), qs_arr(i,j,k) / qr_arr(i,j,k)),
2008  Real(1.0)),
2009  Real(2.0));
2010  psacr_arr(i,j,k) = amrex::min(
2011  psacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2012  }
2013  }
2014 
2015  // pgacr: accretion of rain by graupel [HL A12] [LFO 42]
2016  // (T<T0: R->G) (T>=T0: enhance melting of graupel)
2017  // This depends only on graupel and rain, so it must not be nested inside
2018  // the (qs,qr) test above -- otherwise a rain/graupel column with
2019  // negligible snow would never accrete rain onto graupel
2020  if (qg_arr(i,j,k) > Real(qcrmin) && qr_arr(i,j,k) > Real(qcrmin)) {
2021  const Real acrfac =
2022  Real(5.0) * rslope3_r_arr(i,j,k) * rslope3_r_arr(i,j,k)
2023  * rslope_g_arr(i,j,k)
2024  + Real(2.0) * rslope3_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2025  * rslope2_g_arr(i,j,k)
2026  + Real(0.5) * rslope2_r_arr(i,j,k) * rslope2_r_arr(i,j,k)
2027  * rslope3_g_arr(i,j,k);
2028  pgacr_arr(i,j,k) = Real(pi) * Real(pi) * Real(n0r)
2029  * n0g_loc * std::abs(vt2ave - vt2r)
2030  * (Real(denr) / den_arr(i,j,k)) * acrfac;
2031  pgacr_arr(i,j,k) *= std::pow(
2032  amrex::min(
2033  amrex::max(Real(0.0), qg_arr(i,j,k) / qr_arr(i,j,k)),
2034  Real(1.0)),
2035  Real(2.0));
2036  pgacr_arr(i,j,k) = amrex::min(
2037  pgacr_arr(i,j,k), qr_arr(i,j,k) / dtcld);
2038  }
2039 
2040  if (qg_arr(i,j,k) > Real(qcrmin) && qs_arr(i,j,k) > Real(qcrmin)) {
2041  pgacs_arr(i,j,k) = Real(0.0);
2042  }
2043 
2044  if (supcol <= Real(0.0)) {
2045  const Real xlf = Real(xlf0);
2046  if (qs_arr(i,j,k) > Real(0.0)) {
2047  // WSM6-CPP TAG: PSEML
2048  // legacy_group: G13g
2049  // process: Snow evaporation/sublimation
2050  // compare_vars: pseml, qs, qv, qsat, den
2051  pseml_arr(i,j,k) = amrex::min(
2052  amrex::max(
2053  Real(cliq) * supcol
2054  * (paacw_arr(i,j,k) + psacr_arr(i,j,k)) / xlf,
2055  -qs_arr(i,j,k) / dtcld),
2056  Real(0.0));
2057  }
2058  if (qg_arr(i,j,k) > Real(0.0)) {
2059  pgeml_arr(i,j,k) = amrex::min(
2060  amrex::max(
2061  Real(cliq) * supcol
2062  * (paacw_arr(i,j,k) + pgacr_arr(i,j,k)) / xlf,
2063  -qg_arr(i,j,k) / dtcld),
2064  Real(0.0));
2065  }
2066  }
2067 
2068  if (supcol > Real(0.0)) {
2069  if (qi_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2070  // WSM6-CPP TAG: PIDEP
2071  // legacy_group: G13h
2072  // process: Ice deposition/sublimation
2073  // compare_vars: pidep, qi, qv, qsat, den, t
2074  pidep_arr(i,j,k) = Real(4.0) * diameter * xni
2075  * (rhi_arr(i,j,k) - Real(1.0))
2076  / workdiffi_arr(i,j,k);
2077  Real supice = satdt - prevp_arr(i,j,k);
2078  if (pidep_arr(i,j,k) < Real(0.0)) {
2079  pidep_arr(i,j,k) = amrex::max(
2080  amrex::max(pidep_arr(i,j,k), satdt / Real(2.0)),
2081  supice);
2082  pidep_arr(i,j,k) = amrex::max(
2083  pidep_arr(i,j,k), -qi_arr(i,j,k) / dtcld);
2084  } else {
2085  pidep_arr(i,j,k) = amrex::min(
2086  amrex::min(pidep_arr(i,j,k), satdt / Real(2.0)),
2087  supice);
2088  }
2089  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k))
2090  >= std::abs(satdt)) {
2091  ifsat = 1;
2092  }
2093  }
2094 
2095  if (qs_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2096  const Real coeres = rslope2_s_arr(i,j,k)
2097  * std::sqrt(rslope_s_arr(i,j,k)
2098  * rslopeb_s_arr(i,j,k));
2099  psdep_arr(i,j,k) = (rhi_arr(i,j,k) - Real(1.0))
2100  * n0sfac_arr(i,j,k)
2101  * (Real(precs1) * rslope2_s_arr(i,j,k)
2102  + Real(precs2) * work2_arr(i,j,k) * coeres)
2103  / workdiffi_arr(i,j,k);
2104  Real supice = satdt - prevp_arr(i,j,k) - pidep_arr(i,j,k);
2105  if (psdep_arr(i,j,k) < Real(0.0)) {
2106  psdep_arr(i,j,k) = amrex::max(
2107  psdep_arr(i,j,k), -qs_arr(i,j,k) / dtcld);
2108  psdep_arr(i,j,k) = amrex::max(
2109  amrex::max(psdep_arr(i,j,k), satdt / Real(2.0)),
2110  supice);
2111  } else {
2112  psdep_arr(i,j,k) = amrex::min(
2113  amrex::min(psdep_arr(i,j,k), satdt / Real(2.0)),
2114  supice);
2115  }
2116  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2117  + psdep_arr(i,j,k))
2118  >= std::abs(satdt)) {
2119  ifsat = 1;
2120  }
2121  }
2122 
2123  if (qg_arr(i,j,k) > Real(0.0) && ifsat != 1) {
2124  const Real coeres = rslope2_g_arr(i,j,k)
2125  * std::sqrt(rslope_g_arr(i,j,k)
2126  * rslopeb_g_arr(i,j,k));
2127  pgdep_arr(i,j,k) = (rhi_arr(i,j,k) - Real(1.0))
2128  * (Real(precg1) * rslope2_g_arr(i,j,k)
2129  + Real(precg2) * work2_arr(i,j,k) * coeres)
2130  / workdiffi_arr(i,j,k);
2131  Real supice = satdt - prevp_arr(i,j,k)
2132  - pidep_arr(i,j,k) - psdep_arr(i,j,k);
2133  if (pgdep_arr(i,j,k) < Real(0.0)) {
2134  pgdep_arr(i,j,k) = amrex::max(
2135  pgdep_arr(i,j,k), -qg_arr(i,j,k) / dtcld);
2136  pgdep_arr(i,j,k) = amrex::max(
2137  amrex::max(pgdep_arr(i,j,k), satdt / Real(2.0)),
2138  supice);
2139  } else {
2140  pgdep_arr(i,j,k) = amrex::min(
2141  amrex::min(pgdep_arr(i,j,k), satdt / Real(2.0)),
2142  supice);
2143  }
2144  if (std::abs(prevp_arr(i,j,k) + pidep_arr(i,j,k)
2145  + psdep_arr(i,j,k) + pgdep_arr(i,j,k))
2146  >= std::abs(satdt)) {
2147  ifsat = 1;
2148  }
2149  }
2150 
2151  if (supsat > Real(0.0) && ifsat != 1) {
2152  const Real supice = satdt - prevp_arr(i,j,k)
2153  - pidep_arr(i,j,k)
2154  - psdep_arr(i,j,k)
2155  - pgdep_arr(i,j,k);
2156  const Real xni0 = Real(1.0e3) * std::exp(Real(0.1) * supcol);
2157  const Real roqi0 = Real(4.92e-11) * std::pow(xni0, Real(1.33));
2158  pigen_arr(i,j,k) = amrex::max(
2159  Real(0.0),
2160  (roqi0 / den_arr(i,j,k)
2161  - amrex::max(qi_arr(i,j,k), Real(0.0))) / dtcld);
2162  pigen_arr(i,j,k) = amrex::min(
2163  amrex::min(pigen_arr(i,j,k), satdt), supice);
2164  }
2165 
2166  if (qi_arr(i,j,k) > Real(0.0)) {
2167  const Real qimax = Real(roqimax) / den_arr(i,j,k);
2168  // WSM6-CPP TAG: PSAUT
2169  // legacy_group: G13i
2170  // process: Autoconversion to snow
2171  // compare_vars: psaut, qi, qs, den
2172  psaut_arr(i,j,k) = amrex::max(
2173  Real(0.0), (qi_arr(i,j,k) - qimax) / dtcld);
2174  }
2175 
2176  if (qs_arr(i,j,k) > Real(0.0)) {
2177  const Real alpha2 = Real(1.0e-3)
2178  * std::exp(Real(0.09) * (-supcol));
2179  pgaut_arr(i,j,k) = amrex::min(
2180  amrex::max(
2181  Real(0.0), alpha2 * (qs_arr(i,j,k) - Real(qs0))),
2182  qs_arr(i,j,k) / dtcld);
2183  }
2184  }
2185 
2186  if (supcol < Real(0.0)) {
2187  if (qs_arr(i,j,k) > Real(0.0)
2188  && rhw_arr(i,j,k) < Real(1.0)) {
2189  const Real coeres = rslope2_s_arr(i,j,k)
2190  * std::sqrt(rslope_s_arr(i,j,k)
2191  * rslopeb_s_arr(i,j,k));
2192  // WSM6-CPP TAG: PSEVP
2193  // legacy_group: G13j
2194  // process: Graupel evaporation/sublimation
2195  // compare_vars: psevp, qg, qv, qsat, den
2196  psevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2197  * n0sfac_arr(i,j,k)
2198  * (Real(precs1) * rslope2_s_arr(i,j,k)
2199  + Real(precs2) * work2_arr(i,j,k) * coeres)
2200  / workdiffw_arr(i,j,k);
2201  psevp_arr(i,j,k) = amrex::min(
2202  amrex::max(psevp_arr(i,j,k),
2203  -qs_arr(i,j,k) / dtcld),
2204  Real(0.0));
2205  }
2206 
2207  if (qg_arr(i,j,k) > Real(0.0)
2208  && rhw_arr(i,j,k) < Real(1.0)) {
2209  const Real coeres = rslope2_g_arr(i,j,k)
2210  * std::sqrt(rslope_g_arr(i,j,k)
2211  * rslopeb_g_arr(i,j,k));
2212  pgevp_arr(i,j,k) = (rhw_arr(i,j,k) - Real(1.0))
2213  * (Real(precg1) * rslope2_g_arr(i,j,k)
2214  + Real(precg2) * work2_arr(i,j,k) * coeres)
2215  / workdiffw_arr(i,j,k);
2216  pgevp_arr(i,j,k) = amrex::min(
2217  amrex::max(pgevp_arr(i,j,k),
2218  -qg_arr(i,j,k) / dtcld),
2219  Real(0.0));
2220  }
2221  }
2222  });
2223  // G14: mass conservation check and state update [lines 1200-1388]
2224  // WSM6-CPP TAG: UPDATE
2225  // legacy_group: G14
2226  // process: Mass conservation and state update
2227  // compare_vars: t, q, qci, qrs, qv
2228  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2229  const Real qmin_l = Real(qmin);
2230  const Real qcrmin_l= Real(qcrmin);
2231  const Real t0c_l = Real(t0c);
2232 
2233  const Real delta2 =
2234  (qr_arr(i,j,k) < Real(1.0e-4) && qs_arr(i,j,k) < Real(1.0e-4))
2235  ? Real(1.0) : Real(0.0);
2236  const Real delta3 =
2237  (qr_arr(i,j,k) < Real(1.0e-4)) ? Real(1.0) : Real(0.0);
2238 
2239  if (t_arr(i,j,k) <= t0c_l) {
2240  Real value, source, factor, xlf, xlwork2;
2241 
2242  value = amrex::max(qmin_l, qc_arr(i,j,k));
2243  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2244  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2245 // + paacw_arr(i,j,k)) * dtcld;
2246  if (source > value) {
2247  factor = value / source;
2248  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2249  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2250  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2251  }
2252 
2253  value = amrex::max(qmin_l, qi_arr(i,j,k));
2254  source = (psaut_arr(i,j,k) - pigen_arr(i,j,k)
2255  - pidep_arr(i,j,k) + praci_arr(i,j,k)
2256  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)) * dtcld;
2257  if (source > value) {
2258  factor = value / source;
2259  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2260  pigen_arr(i,j,k) = pigen_arr(i,j,k) * factor;
2261  pidep_arr(i,j,k) = pidep_arr(i,j,k) * factor;
2262  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2263  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2264  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2265  }
2266 
2267  value = amrex::max(qmin_l, qr_arr(i,j,k));
2268  source = (-praut_arr(i,j,k) - prevp_arr(i,j,k)
2269  - pracw_arr(i,j,k) + piacr_arr(i,j,k)
2270  + psacr_arr(i,j,k) + pgacr_arr(i,j,k)) * dtcld;
2271  if (source > value) {
2272  factor = value / source;
2273  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2274  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2275  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2276  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2277  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2278  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2279  }
2280 
2281  value = amrex::max(qmin_l, qs_arr(i,j,k));
2282  source = -(psdep_arr(i,j,k) + psaut_arr(i,j,k)
2283  - pgaut_arr(i,j,k) + paacw_arr(i,j,k)
2284  + piacr_arr(i,j,k) * delta3
2285  + praci_arr(i,j,k) * delta3
2286  - pracs_arr(i,j,k) * (Real(1.0) - delta2)
2287  + psacr_arr(i,j,k) * delta2
2288  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)) * dtcld;
2289  if (source > value) {
2290  factor = value / source;
2291  psdep_arr(i,j,k) = psdep_arr(i,j,k) * factor;
2292  psaut_arr(i,j,k) = psaut_arr(i,j,k) * factor;
2293  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2294  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2295  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2296  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2297  psaci_arr(i,j,k) = psaci_arr(i,j,k) * factor;
2298  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2299  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2300  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2301  }
2302 
2303  value = amrex::max(qmin_l, qg_arr(i,j,k));
2304  source = -(pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2305  + piacr_arr(i,j,k) * (Real(1.0) - delta3)
2306  + praci_arr(i,j,k) * (Real(1.0) - delta3)
2307  + psacr_arr(i,j,k) * (Real(1.0) - delta2)
2308  + pracs_arr(i,j,k) * (Real(1.0) - delta2)
2309  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2310  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k)) * dtcld;
2311  if (source > value) {
2312  factor = value / source;
2313  pgdep_arr(i,j,k) = pgdep_arr(i,j,k) * factor;
2314  pgaut_arr(i,j,k) = pgaut_arr(i,j,k) * factor;
2315  piacr_arr(i,j,k) = piacr_arr(i,j,k) * factor;
2316  praci_arr(i,j,k) = praci_arr(i,j,k) * factor;
2317  psacr_arr(i,j,k) = psacr_arr(i,j,k) * factor;
2318  pracs_arr(i,j,k) = pracs_arr(i,j,k) * factor;
2319  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2320  pgaci_arr(i,j,k) = pgaci_arr(i,j,k) * factor;
2321  pgacr_arr(i,j,k) = pgacr_arr(i,j,k) * factor;
2322  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2323  }
2324 
2325  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psdep_arr(i,j,k)
2326  + pgdep_arr(i,j,k) + pigen_arr(i,j,k)
2327  + pidep_arr(i,j,k));
2328  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
2329 // + paacw_arr(i,j,k))
2330  qc_arr(i,j,k) = amrex::max(
2331  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
2332  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
2333  * dtcld,
2334  Real(0.0));
2335  qr_arr(i,j,k) = amrex::max(
2336  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
2337  + prevp_arr(i,j,k) - piacr_arr(i,j,k)
2338  - pgacr_arr(i,j,k) - psacr_arr(i,j,k))
2339  * dtcld,
2340  Real(0.0));
2341  qi_arr(i,j,k) = amrex::max(
2342  qi_arr(i,j,k) - (psaut_arr(i,j,k) + praci_arr(i,j,k)
2343  + psaci_arr(i,j,k) + pgaci_arr(i,j,k)
2344  - pigen_arr(i,j,k) - pidep_arr(i,j,k))
2345  * dtcld,
2346  Real(0.0));
2347  qs_arr(i,j,k) = amrex::max(
2348  qs_arr(i,j,k) + (psdep_arr(i,j,k) + psaut_arr(i,j,k)
2349  + paacw_arr(i,j,k) - pgaut_arr(i,j,k)
2350  + piacr_arr(i,j,k) * delta3
2351  + praci_arr(i,j,k) * delta3
2352  + psaci_arr(i,j,k) - pgacs_arr(i,j,k)
2353  - pracs_arr(i,j,k) * (Real(1.0) - delta2)
2354  + psacr_arr(i,j,k) * delta2) * dtcld,
2355  Real(0.0));
2356  qg_arr(i,j,k) = amrex::max(
2357  qg_arr(i,j,k) + (pgdep_arr(i,j,k) + pgaut_arr(i,j,k)
2358  + piacr_arr(i,j,k) * (Real(1.0) - delta3)
2359  + praci_arr(i,j,k) * (Real(1.0) - delta3)
2360  + psacr_arr(i,j,k) * (Real(1.0) - delta2)
2361  + pracs_arr(i,j,k) * (Real(1.0) - delta2)
2362  + pgaci_arr(i,j,k) + paacw_arr(i,j,k)
2363  + pgacr_arr(i,j,k) + pgacs_arr(i,j,k))
2364  * dtcld,
2365  Real(0.0));
2366  xlf = Real(xls) - xl_arr(i,j,k);
2367  xlwork2 = -Real(xls) * (psdep_arr(i,j,k) + pgdep_arr(i,j,k)
2368  + pidep_arr(i,j,k) + pigen_arr(i,j,k))
2369  - xl_arr(i,j,k) * prevp_arr(i,j,k)
2370  - xlf * (piacr_arr(i,j,k) + paacw_arr(i,j,k)
2371  + paacw_arr(i,j,k) + pgacr_arr(i,j,k)
2372  + psacr_arr(i,j,k));
2373  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
2374  } else {
2375  Real value, source, factor, xlf, xlwork2;
2376 
2377  value = amrex::max(qmin_l, qc_arr(i,j,k));
2378  source = (praut_arr(i,j,k) + pracw_arr(i,j,k)
2379  + paacw_arr(i,j,k) + paacw_arr(i,j,k)) * dtcld;
2380 // + paacw_arr(i,j,k)) * dtcld;
2381  if (source > value) {
2382  factor = value / source;
2383  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2384  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2385  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2386  }
2387 
2388  value = amrex::max(qmin_l, qr_arr(i,j,k));
2389  source = (-paacw_arr(i,j,k) - praut_arr(i,j,k)
2390  + pseml_arr(i,j,k) + pgeml_arr(i,j,k)
2391  - pracw_arr(i,j,k) - paacw_arr(i,j,k)
2392  - prevp_arr(i,j,k)) * dtcld;
2393  if (source > value) {
2394  factor = value / source;
2395  praut_arr(i,j,k) = praut_arr(i,j,k) * factor;
2396  prevp_arr(i,j,k) = prevp_arr(i,j,k) * factor;
2397  pracw_arr(i,j,k) = pracw_arr(i,j,k) * factor;
2398  paacw_arr(i,j,k) = paacw_arr(i,j,k) * factor;
2399  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2400  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2401  }
2402 
2403  value = amrex::max(qcrmin_l, qs_arr(i,j,k));
2404  source = (pgacs_arr(i,j,k) - pseml_arr(i,j,k)
2405  - psevp_arr(i,j,k)) * dtcld;
2406  if (source > value) {
2407  factor = value / source;
2408  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2409  psevp_arr(i,j,k) = psevp_arr(i,j,k) * factor;
2410  pseml_arr(i,j,k) = pseml_arr(i,j,k) * factor;
2411  }
2412 
2413  value = amrex::max(qcrmin_l, qg_arr(i,j,k));
2414  source = -(pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
2415  + pgeml_arr(i,j,k)) * dtcld;
2416  if (source > value) {
2417  factor = value / source;
2418  pgacs_arr(i,j,k) = pgacs_arr(i,j,k) * factor;
2419  pgevp_arr(i,j,k) = pgevp_arr(i,j,k) * factor;
2420  pgeml_arr(i,j,k) = pgeml_arr(i,j,k) * factor;
2421  }
2422 
2423  work2_arr(i,j,k) = -(prevp_arr(i,j,k) + psevp_arr(i,j,k)
2424  + pgevp_arr(i,j,k));
2425  qv_arr(i,j,k) = qv_arr(i,j,k) + work2_arr(i,j,k) * dtcld;
2426 // + paacw_arr(i,j,k))
2427  qc_arr(i,j,k) = amrex::max(
2428  qc_arr(i,j,k) - (praut_arr(i,j,k) + pracw_arr(i,j,k)
2429  + paacw_arr(i,j,k) + paacw_arr(i,j,k))
2430  * dtcld,
2431  Real(0.0));
2432  qr_arr(i,j,k) = amrex::max(
2433  qr_arr(i,j,k) + (praut_arr(i,j,k) + pracw_arr(i,j,k)
2434  + prevp_arr(i,j,k) + paacw_arr(i,j,k)
2435  + paacw_arr(i,j,k) - pseml_arr(i,j,k)
2436  - pgeml_arr(i,j,k)) * dtcld,
2437  Real(0.0));
2438  qs_arr(i,j,k) = amrex::max(
2439  qs_arr(i,j,k) + (psevp_arr(i,j,k) - pgacs_arr(i,j,k)
2440  + pseml_arr(i,j,k)) * dtcld,
2441  Real(0.0));
2442  qg_arr(i,j,k) = amrex::max(
2443  qg_arr(i,j,k) + (pgacs_arr(i,j,k) + pgevp_arr(i,j,k)
2444  + pgeml_arr(i,j,k)) * dtcld,
2445  Real(0.0));
2446  xlf = Real(xls) - xl_arr(i,j,k);
2447  xlwork2 = -xl_arr(i,j,k) * (prevp_arr(i,j,k)
2448  + psevp_arr(i,j,k)
2449  + pgevp_arr(i,j,k))
2450  - xlf * (pseml_arr(i,j,k) + pgeml_arr(i,j,k));
2451  t_arr(i,j,k) = t_arr(i,j,k) - xlwork2 / cpm_arr(i,j,k) * dtcld;
2452  }
2453  });
2454  // G15: second qsat computation [lines 1390-1420]
2455  // WSM6-CPP TAG: QSAT2
2456  // legacy_group: G15
2457  // process: Second saturation mixing ratio computation
2458  // compare_vars: qs, qvs, den, denfac, t, p
2459  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2460  const Real ttp = Real(t0c) + Real(0.01);
2461  const Real dldt = Real(cpv) - Real(cliq);
2462  const Real xa = -dldt / Real(rv);
2463  const Real xb = xa + Real(xlv0) / (Real(rv) * ttp);
2464 
2465  Real tr = ttp / t_arr(i,j,k);
2466  Real qsw = Real(psat) * std::exp(std::log(tr) * xa)
2467  * std::exp(xb * (Real(1.0) - tr));
2468  qsw = amrex::min(qsw, Real(0.99) * p_arr(i,j,k));
2469  qsatw_arr(i,j,k) = Real(ep2) * qsw / (p_arr(i,j,k) - qsw);
2470  qsatw_arr(i,j,k) = amrex::max(qsatw_arr(i,j,k), Real(qmin));
2471  });
2472  // G16: pcond condensational/evaporational update [lines 1427-1437]
2473  // WSM6-CPP TAG: PCOND
2474  // legacy_group: G16
2475  // process: Condensation/evaporation update
2476  // compare_vars: pcond, t, qv, qc, qsat
2477  if (m_do_cond) {
2478  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2479  const Real workcond = wsm6_conden(
2480  t_arr(i,j,k), qv_arr(i,j,k), qsatw_arr(i,j,k),
2481  xl_arr(i,j,k), cpm_arr(i,j,k), Real(qmin), Real(rv));
2482  const Real work2loc = qc_arr(i,j,k) + workcond;
2483  static_cast<void>(work2loc);
2484  pcond_arr(i,j,k) = amrex::min(
2485  amrex::max(workcond / dtcld, Real(0.0)),
2486  amrex::max(qv_arr(i,j,k), Real(0.0)) / dtcld);
2487  if (qc_arr(i,j,k) > Real(0.0) && workcond < Real(0.0)) {
2488  pcond_arr(i,j,k) = amrex::max(workcond, -qc_arr(i,j,k)) / dtcld;
2489  }
2490  qv_arr(i,j,k) = qv_arr(i,j,k) - pcond_arr(i,j,k) * dtcld;
2491  qc_arr(i,j,k) = amrex::max(
2492  qc_arr(i,j,k) + pcond_arr(i,j,k) * dtcld,
2493  Real(0.0));
2494  t_arr(i,j,k) = t_arr(i,j,k)
2495  + pcond_arr(i,j,k) * xl_arr(i,j,k)
2496  / cpm_arr(i,j,k) * dtcld;
2497  });
2498  }
2499  // G17: padding for small values [lines 1444-1449]
2500  // WSM6-CPP TAG: CLIP
2501  // legacy_group: G17
2502  // process: Padding/clipping for small values
2503  // compare_vars: qv, qc, qr, qi, qs, qg
2504  ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
2505  if (qc_arr(i,j,k) <= Real(qmin)) qc_arr(i,j,k) = Real(0.0);
2506  if (qi_arr(i,j,k) <= Real(qmin)) qi_arr(i,j,k) = Real(0.0);
2507  });
2508 
2509  }
2510 #ifdef ERF_USE_WSM6_FORT
2511  }
2512 #endif
2513  ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
2514  rain_arr(i,j,klo) = rainacc_arr(i,j,0);
2515  snow_arr(i,j,klo) = snowacc_arr(i,j,0);
2516  graup_arr(i,j,klo) = graupacc_arr(i,j,0);
2517  });
2518  }
2519 }
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
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_wdm6::pacrc
real(kind=kind_phys), save pacrc
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precs1
real(kind=kind_phys), save precs1
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::pacrr
real(kind=kind_phys), save pacrr
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precr2
real(kind=kind_phys), save precr2
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::pvtr
real(kind=kind_phys), save pvtr
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precr1
real(kind=kind_phys), save precr1
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::pvtg
real(kind=kind_phys), save pvtg
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::g6pbr
real(kind=kind_phys), save g6pbr
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::dens
real(kind=kind_phys), parameter, private dens
Definition: ERF_module_mp_wdm6.F90:61
mp_wdm6::roqimax
real(kind=kind_phys), save roqimax
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::qck1
real(kind=kind_phys), save qck1
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::qc0
real(kind=kind_phys), save qc0
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::pvts
real(kind=kind_phys), save pvts
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precs2
real(kind=kind_phys), save precs2
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::pacrg
real(kind=kind_phys), save pacrg
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precg1
real(kind=kind_phys), save precg1
Definition: ERF_module_mp_wdm6.F90:100
mp_wdm6::precg2
real(kind=kind_phys), save precg2
Definition: ERF_module_mp_wdm6.F90:100
Here is the call graph for this function:

◆ Copy_Micro_to_State()

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

Reimplemented from NullMoist.

8 {
9  for (MFIter mfi(cons, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
10  const auto& box3d = mfi.tilebox();
11  auto states = cons.array(mfi);
12 
13  auto rho = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
14  auto theta = mic_fab_vars[MicVar_WSM6::theta]->array(mfi);
15  auto tabs = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
16 
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:48
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:46
RhoQ3_comp
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:47
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
RhoQ6_comp
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:50
RhoQ5_comp
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:49
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:20
MicVar_Kess::tabs
@ tabs
Definition: ERF_Kessler.H:26
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:24
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:30
MicVar_SatAdj::qc
@ qc
Definition: ERF_SatAdj.H:41
MicVar_WDM6::qi
@ qi
Definition: ERF_WDM6.H:27
MicVar_WDM6::qg
@ qg
Definition: ERF_WDM6.H:30
MicVar_WDM6::qs
@ qs
Definition: ERF_WDM6.H:29
MicVar_WSM6::theta
@ theta
Definition: ERF_WSM6.H:22
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
WDM6SedCellScratch::qr
@ qr
Definition: ERF_AdvanceWDM6.cpp:269

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro()

void 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:39
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:2124
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2141
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:2063
Here is the call graph for this function:

◆ Get_Surface_Precip_Accumulation_Ptrs()

SurfacePrecipAccumulationSources 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:34
SurfacePrecipAccumulationSources::snow
SurfacePrecipAccumulationSource snow
Definition: ERF_SurfacePrecipitation.H:37
SurfacePrecipAccumulationSources::total
SurfacePrecipAccumulationSource total
Definition: ERF_SurfacePrecipitation.H:35
SurfacePrecipAccumulationSources::graupel
SurfacePrecipAccumulationSource graupel
Definition: ERF_SurfacePrecipitation.H:38

◆ Init()

void 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_wdm6::rgmma
real(kind=kind_phys) function rgmma(x)
Definition: ERF_module_mp_wdm6.F90:2174
Here is the call graph for this function:

◆ Qmoist_Ptr()

amrex::MultiFab* 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
36 { }

◆ 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
43  {
44  Update_Micro_Vars(cons_in);
45  }
WSM6::Update_Micro_Vars
virtual void Update_Micro_Vars(amrex::MultiFab &)
Definition: ERF_NullMoist.H:36

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