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_State_to_Micro (const amrex::MultiFab &cons_in, const amrex::MultiFab *base_state)
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab *base_state) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in, const amrex::MultiFab &) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &solverChoice) override
 
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
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 
virtual void GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const
 
virtual void SetCurrentLevel (const int &lev)
 
virtual void InitLevel (const int, const amrex::MultiFab &)
 
virtual int getDiagnosticsInterval () const
 
virtual void Set_Lmask (amrex::iMultiFab *)
 Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this. More...
 
virtual void Set_RealWidth (const int)
 

Static Public Attributes

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

Private Types

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

Private Member Functions

void initialize_coeffs ()
 

Private Attributes

int m_qmoist_size = 3
 
int n_qstate_moist_size = 6
 
int n_qstate_moist_numconc_size = 0
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::Real dt {0.0}
 
amrex::Real m_dzmin {0.0}
 
int nlev {0}
 
int zlo {0}
 
int zhi {0}
 
int m_axis {2}
 
bool m_do_cond {true}
 
amrex::Real m_rdOcp {RdoCp}
 
MoistureType m_moisture_type {MoistureType::None}
 
amrex::MultiFab * m_z_phys_nd {nullptr}
 
amrex::MultiFab * m_detJ_cc {nullptr}
 
amrex::Array< FabPtr, MicVar_WSM6::NumVars > mic_fab_vars
 
bool m_hail_opt {false}
 
amrex::Real m_n0g {0}
 
amrex::Real m_deng {0}
 
amrex::Real m_avtg {0}
 
amrex::Real m_bvtg {0}
 
amrex::Real m_lamdagmax {0}
 
amrex::Real m_pi_wsm6 {0}
 
amrex::Real m_xlv1 {0}
 
amrex::Real m_qc0 {0}
 
amrex::Real m_qck1 {0}
 
amrex::Real m_pidnc {0}
 
amrex::Real m_bvtr1 {0}
 
amrex::Real m_bvtr2 {0}
 
amrex::Real m_bvtr3 {0}
 
amrex::Real m_bvtr4 {0}
 
amrex::Real m_bvtr6 {0}
 
amrex::Real m_g1pbr {0}
 
amrex::Real m_g3pbr {0}
 
amrex::Real m_g4pbr {0}
 
amrex::Real m_g5pbro2 {0}
 
amrex::Real m_g6pbr {0}
 
amrex::Real m_pvtr {0}
 
amrex::Real m_eacrr {0}
 
amrex::Real m_pacrr {0}
 
amrex::Real m_precr1 {0}
 
amrex::Real m_precr2 {0}
 
amrex::Real m_roqimax {0}
 
amrex::Real m_bvts1 {0}
 
amrex::Real m_bvts2 {0}
 
amrex::Real m_bvts3 {0}
 
amrex::Real m_bvts4 {0}
 
amrex::Real m_g1pbs {0}
 
amrex::Real m_g3pbs {0}
 
amrex::Real m_g4pbs {0}
 
amrex::Real m_g5pbso2 {0}
 
amrex::Real m_pvts {0}
 
amrex::Real m_pacrs {0}
 
amrex::Real m_precs1 {0}
 
amrex::Real m_precs2 {0}
 
amrex::Real m_pidn0r {0}
 
amrex::Real m_pidn0s {0}
 
amrex::Real m_pacrc {0}
 
amrex::Real m_bvtg1 {0}
 
amrex::Real m_bvtg2 {0}
 
amrex::Real m_bvtg3 {0}
 
amrex::Real m_bvtg4 {0}
 
amrex::Real m_g1pbg {0}
 
amrex::Real m_g3pbg {0}
 
amrex::Real m_g4pbg {0}
 
amrex::Real m_g5pbgo2 {0}
 
amrex::Real m_pvtg {0}
 
amrex::Real m_pacrg {0}
 
amrex::Real m_precg1 {0}
 
amrex::Real m_precg2 {0}
 
amrex::Real m_pidn0g {0}
 
amrex::Real m_rslopermax {0}
 
amrex::Real m_rslopesmax {0}
 
amrex::Real m_rslopegmax {0}
 
amrex::Real m_rsloperbmax {0}
 
amrex::Real m_rslopesbmax {0}
 
amrex::Real m_rslopegbmax {0}
 
amrex::Real m_rsloper2max {0}
 
amrex::Real m_rslopes2max {0}
 
amrex::Real m_rslopeg2max {0}
 
amrex::Real m_rsloper3max {0}
 
amrex::Real m_rslopes3max {0}
 
amrex::Real m_rslopeg3max {0}
 

Additional Inherited Members

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

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ WSM6()

WSM6::WSM6 ( )
inline
125 {}

◆ ~WSM6()

virtual WSM6::~WSM6 ( )
virtualdefault

Member Function Documentation

◆ Advance()

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

Reimplemented from NullMoist.

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

Referenced by Update_State_Vars().

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

◆ Copy_State_to_Micro() [1/2]

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

Reimplemented from NullMoist.

Referenced by Update_Micro_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro() [2/2]

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

◆ Define()

void WSM6::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

129  {
130  m_moisture_type = sc.moisture_type;
131  m_rdOcp = sc.rdOcp;
132  set_anelastic_reference_pressure_mode(sc);
133  m_axis = sc.ave_plane;
134  m_do_cond = (!sc.uses_shoc_family());
135  }
NullMoist::set_anelastic_reference_pressure_mode
void set_anelastic_reference_pressure_mode(const SolverChoice &sc)
Definition: ERF_NullMoist.H:155
WSM6::m_axis
int m_axis
Definition: ERF_WSM6.H:232
WSM6::m_moisture_type
MoistureType m_moisture_type
Definition: ERF_WSM6.H:235
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2237
SolverChoice::ave_plane
int ave_plane
Averaging plane index used by diagnostics.
Definition: ERF_DataStruct.H:2260
SolverChoice::rdOcp
amrex::Real rdOcp
Ratio of dry-air gas constant to c_p.
Definition: ERF_DataStruct.H:2062
SolverChoice::uses_shoc_family
bool uses_shoc_family() const noexcept
Query whether any SHOC-family PBL scheme is active.
Definition: ERF_DataStruct.H:2177
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.

213  {
214  SurfacePrecipAccumulationSources sources;
215  sources.total = {mic_fab_vars[MicVar_WSM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
216  sources.snow = {mic_fab_vars[MicVar_WSM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
217  sources.graupel = {mic_fab_vars[MicVar_WSM6::graup_accum].get(), rhoh2o / amrex::Real(1000.0)};
218  return sources;
219  }
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.

17 {
18  dt = dt_advance;
19  m_geom = geom;
20 
21  m_z_phys_nd = z_phys_nd.get();
22  m_detJ_cc = detJ_cc.get();
23 
24  MicVarMap.resize(m_qmoist_size);
25  MicVarMap = {MicVar_WSM6::rain_accum, MicVar_WSM6::snow_accum, MicVar_WSM6::graup_accum};
26 
27 #if defined(ERF_USE_WSM6_FORT) && defined(AMREX_USE_GPU)
28  // The host-only Fortran bridge receives dataPtr() from mic_fab_vars.
29  Arena* Arena_Used = The_Managed_Arena();
30 #else
31  Arena* Arena_Used = The_Arena();
32 #endif
33 
34  for (int ivar = 0; ivar < MicVar_WSM6::NumVars; ++ivar) {
35  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
36  1, cons_in.nGrowVect(),
37  MFInfo().SetArena(Arena_Used));
38  mic_fab_vars[ivar]->setVal(0.0);
39  }
40 
41  nlev = m_geom.Domain().length(2);
42  zlo = m_geom.Domain().smallEnd(2);
43  zhi = m_geom.Domain().bigEnd(2);
44 
45  initialize_coeffs();
46 }
WSM6::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_WSM6.H:238
WSM6::zlo
int zlo
Definition: ERF_WSM6.H:231
WSM6::nlev
int nlev
Definition: ERF_WSM6.H:231
WSM6::zhi
int zhi
Definition: ERF_WSM6.H:231
WSM6::m_qmoist_size
int m_qmoist_size
Definition: ERF_WSM6.H:222
WSM6::initialize_coeffs
void initialize_coeffs()
Definition: ERF_InitWSM6.cpp:97
WSM6::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_WSM6.H:226
MicVar_WSM6::NumVars
@ NumVars
Definition: ERF_WSM6.H:35

◆ initialize_coeffs()

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

194  {
195  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
196  return mic_fab_vars[MicVarMap[varIdx]].get();
197  }
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.

206  {
207  a_idx = {0, 1, 2};
208  a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
209  }

◆ Qmoist_Size()

int WSM6::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

199 { return m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int WSM6::Qstate_Moist_NumConc_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

201 { return n_qstate_moist_numconc_size; }
WSM6::n_qstate_moist_numconc_size
int n_qstate_moist_numconc_size
Definition: ERF_WSM6.H:224

◆ Qstate_Moist_Size()

int WSM6::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

200 { return n_qstate_moist_size; }
WSM6::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_WSM6.H:223

◆ Set_dzmin()

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

Reimplemented from NullMoist.

169 { m_dzmin = dz_min; }
WSM6::m_dzmin
amrex::Real m_dzmin
Definition: ERF_WSM6.H:230

◆ Update_Micro_Vars() [1/2]

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

Reimplemented from NullMoist.

180  {
181  Copy_State_to_Micro(cons_in);
182  }
WSM6::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Here is the call graph for this function:

◆ Update_Micro_Vars() [2/2]

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

Reimplemented from NullMoist.

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

186  {
187  Copy_Micro_to_State(cons_in);
188  }
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_rdOcp

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

Referenced by Define().

◆ m_roqimax

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

◆ m_rslopeg2max

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

◆ m_rslopeg3max

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

◆ m_rslopegbmax

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

◆ m_rslopegmax

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

◆ m_rsloper2max

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

◆ m_rsloper3max

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

◆ m_rsloperbmax

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

◆ m_rslopermax

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

◆ m_rslopes2max

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

◆ m_rslopes3max

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

◆ m_rslopesbmax

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

◆ m_rslopesmax

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

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