Wsm6 microphysics initial Fortran setup

CMake/BuildERFExe.cmake

@@ -233,6 +233,17 @@ function(build_erf_lib erf_lib_name)
   target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MORR_FORT)
   endif()
 
+  if(ERF_ENABLE_WSM6_FORT)
+    target_sources(${erf_lib_name}
+       PRIVATE
+         ${SRC_DIR}/Microphysics/WSM6/ERF_module_libmassv.F90
+         ${SRC_DIR}/Microphysics/WSM6/ERF_mp_radar.F90
+         ${SRC_DIR}/Microphysics/WSM6/ERF_module_mp_wsm6.F90
+         ${SRC_DIR}/Microphysics/WSM6/ERF_module_mp_wsm6_isohelper.F90
+         )
+    target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_WSM6_FORT)
+  endif()
+
   if(ERF_ENABLE_WINDFARM)
     target_sources(${erf_lib_name} PRIVATE
       ${SRC_DIR}/Initialization/ERF_InitWindFarm.cpp
@@ -333,6 +344,9 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Microphysics/Morrison/ERF_AdvanceMorrison.cpp
        ${SRC_DIR}/Microphysics/Morrison/ERF_UpdateMorrison.cpp
        ${SRC_DIR}/Microphysics/Morrison/ERF_Morrison_Plot.cpp
+       ${SRC_DIR}/Microphysics/WSM6/ERF_InitWSM6.cpp
+       ${SRC_DIR}/Microphysics/WSM6/ERF_AdvanceWSM6.cpp
+       ${SRC_DIR}/Microphysics/WSM6/ERF_UpdateWSM6.cpp
        ${SRC_DIR}/Microphysics/SAM/ERF_InitSAM.cpp
        ${SRC_DIR}/Microphysics/SAM/ERF_CloudSAM.cpp
        ${SRC_DIR}/Microphysics/SAM/ERF_IceFall.cpp
@@ -436,7 +450,7 @@ function(build_erf_lib erf_lib_name)
 
   if(ERF_ENABLE_MPI)
     target_link_libraries(${erf_lib_name} PUBLIC $<$<BOOL:${MPI_CXX_FOUND}>:MPI::MPI_CXX>)
-    if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_NOAHMP)
+    if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_WSM6_FORT OR ERF_ENABLE_NOAHMP)
       target_link_libraries(${erf_lib_name} PUBLIC $<$<BOOL:${MPI_CXX_FOUND}>:MPI::MPI_Fortran>)
     endif()
   endif()
@@ -466,6 +480,7 @@ function(build_erf_lib erf_lib_name)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/SAM>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/Kessler>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/Morrison>)
+  target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/WSM6>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/SatAdj>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/Microphysics/SuperDropletsMoist>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/WindFarmParametrization>)

CMake/CrayCompilerDetection.cmake

@@ -96,7 +96,7 @@ else()
 endif()
 
 # Find Cray Fortran compiler wrapper (if needed)
-if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_NOAHMP)
+if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_WSM6_FORT OR ERF_ENABLE_NOAHMP)
     find_program(ERF_CRAY_FC ftn)
     if(ERF_CRAY_FC)
         set(CMAKE_Fortran_COMPILER "${ERF_CRAY_FC}" CACHE FILEPATH "Fortran compiler")
@@ -351,4 +351,4 @@ if(DEFINED ENV{MPICH_GPU_SUPPORT_ENABLED} AND "$ENV{MPICH_GPU_SUPPORT_ENABLED}"
     message(STATUS "  CMAKE_CXX_STANDARD_LIBRARIES: ${CMAKE_CXX_STANDARD_LIBRARIES}")
 else()
     message(STATUS "  GPU-aware MPI not enabled")
-endif()
+endif()

CMake/SetAmrexOptions.cmake

@@ -30,7 +30,7 @@ if(ERF_ENABLE_TESTS OR ERF_ENABLE_FCOMPARE)
 else()
   set(AMReX_PLOTFILE_TOOLS OFF)
 endif()
-if(ERF_ENABLE_MORR_FORT)
+if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_WSM6_FORT)
    set(AMReX_FORTRAN ON)
 endif()
 set(AMReX_EB ON)

CMakeLists.txt

@@ -9,7 +9,7 @@ list(APPEND CMAKE_MESSAGE_CONTEXT "ERF")
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/CMake")
 include(CrayCompilerDetection)
 
-if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_NOAHMP)
+if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_WSM6_FORT OR ERF_ENABLE_NOAHMP)
   project(ERF CXX C Fortran)
 else()
   project(ERF CXX C)
@@ -149,6 +149,7 @@ option(ERF_ENABLE_NOAHMP "Enable Noah-MP" OFF)
 
 #Other parameterizations
 option(ERF_ENABLE_WINDFARM "Enable wind farm parameterizations" OFF)
+option(ERF_ENABLE_WSM6_FORT "Enable WSM6 Fortran bridge calls" OFF)
 
 #Option to build tools
 option(ERF_ENABLE_TOOLS "Enable building of additional tools" OFF)
@@ -448,7 +449,7 @@ if(ERF_ENABLE_MPI)
     # Normal path: Use find_package with QUIET
     message(VERBOSE "Using find_package(MPI) for detection")
     set(_mpi_comps C CXX)
-    if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_NOAHMP)
+    if(ERF_ENABLE_MORR_FORT OR ERF_ENABLE_WSM6_FORT OR ERF_ENABLE_NOAHMP)
       list(APPEND _mpi_comps Fortran)
     endif()
     message(DEBUG "MPI components: ${_mpi_comps}")

Exec/Make.ERF

@@ -2,11 +2,16 @@ AMREX_HOME ?= $(ERF_HOME)/Submodules/AMReX
 EKAT_HOME  ?= $(ERF_HOME)/Submodules/ekat/build/install
 EAMXX_HOME ?= $(ERF_HOME)/external/E3SM/components/eamxx/src
 
+# Default is no Fortran build unless an explicit Fortran bridge/model is enabled.
+# Keep this TRUE when USE_MORR_FORT and USE_WSM6_FORT are both unset.
+BL_NO_FORT = TRUE
 ifeq ($(USE_MORR_FORT), TRUE)
 DEFINES += -DERF_USE_MORR_FORT
 BL_NO_FORT = FALSE
-else
-BL_NO_FORT = TRUE
+endif
+ifeq ($(USE_WSM6_FORT), TRUE)
+DEFINES += -DERF_USE_WSM6_FORT
+BL_NO_FORT = FALSE
 endif
 
 ifeq ($(USE_CUDA), TRUE)
@@ -143,6 +148,11 @@ include $(ERF_MOISTURE_MORRISON_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_MOISTURE_MORRISON_DIR)
 INCLUDE_LOCATIONS += $(ERF_MOISTURE_MORRISON_DIR)
 
+ERF_MOISTURE_WSM6_DIR = $(ERF_SOURCE_DIR)/Microphysics/WSM6
+include $(ERF_MOISTURE_WSM6_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_MOISTURE_WSM6_DIR)
+INCLUDE_LOCATIONS += $(ERF_MOISTURE_WSM6_DIR)
+
 ERF_MOISTURE_SATADJ_DIR = $(ERF_SOURCE_DIR)/Microphysics/SatAdj
 include $(ERF_MOISTURE_SATADJ_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_MOISTURE_SATADJ_DIR)

Source/DataStructs/ERF_DataStruct.H

@@ -70,7 +70,7 @@ AMREX_ENUM(MoistureModelType,
 );
 
 AMREX_ENUM(MoistureType,
-    SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, Morrison, Morrison_NoIce, SuperDroplets, None
+    SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, Morrison, Morrison_NoIce, WSM6, SuperDroplets, None
 );
 
 AMREX_ENUM(WindFarmType,
@@ -158,6 +158,7 @@ struct SolverChoice {
         }
         pp.query_enum_case_insensitive("moisture_model",moisture_type);
         if ( (moisture_type == MoistureType::Morrison) ||
+             (moisture_type == MoistureType::WSM6) ||
              (moisture_type == MoistureType::SAM) ) {
             moisture_indices = MoistureComponentIndices(
                 RhoQ1_comp,  // water vapor
@@ -211,6 +212,7 @@ struct SolverChoice {
                    moisture_type == MoistureType::SAM_NoPrecip_NoIce ||
                    moisture_type == MoistureType::Morrison ||
                    moisture_type == MoistureType::Morrison_NoIce ||
+                   moisture_type == MoistureType::WSM6 ||
                    moisture_type == MoistureType::SatAdj) )
             {
                 for (int i = 0; i <= max_level; ++i) {
@@ -906,6 +908,8 @@ struct SolverChoice {
             amrex::Print() << "Moisture Model: Morrison" << std::endl;
         } else if (moisture_type == MoistureType::Morrison_NoIce) {
             amrex::Print() << "Moisture Model: Morrison_NoIce" << std::endl;
+        } else if (moisture_type == MoistureType::WSM6) {
+            amrex::Print() << "Moisture Model: WSM6" << std::endl;
         } else if (moisture_type == MoistureType::Kessler) {
             amrex::Print() << "Moisture Model: Kessler" << std::endl;
         } else if (moisture_type == MoistureType::Kessler_NoRain) {

Source/ERF.cpp

@@ -1873,6 +1873,8 @@ ERF::InitData_post ()
 void
 ERF::Interp2DArrays (int lev, const BoxArray& my_ba2d, const DistributionMapping& my_dm)
 {
+    if (lev == 0) { return; }
+
     if (lon_m[lev-1] && !lon_m[lev]) {
         auto ngv = lon_m[lev-1]->nGrowVect(); ngv[2] = 0;
         lon_m[lev] = std::make_unique<MultiFab>(my_ba2d,my_dm,1,ngv);

Source/IO/ERF_Plotfile.cpp

@@ -1197,7 +1197,8 @@ ERF::Write3DPlotFile (int which, PlotFileType plotfile_type, Vector<std::string>
                 }
             }
             else if ( (solverChoice.moisture_type == MoistureType::SAM) ||
-                      (solverChoice.moisture_type == MoistureType::Morrison) )
+                      (solverChoice.moisture_type == MoistureType::Morrison) ||
+                      (solverChoice.moisture_type == MoistureType::WSM6) )
             {
                 int offset = (solverChoice.moisture_type == MoistureType::Morrison) ? 5 : 0;
                 if (containerHasElement(plot_var_names, "rain_accum"))

Source/Microphysics/ERF_EulerianMicrophysics.H

@@ -9,6 +9,7 @@
 #include "ERF_SAM.H"
 #include "ERF_Kessler.H"
 #include "ERF_Morrison.H"
+#include "ERF_WSM6.H"
 #include "ERF_SatAdj.H"
 #include "ERF_Microphysics.H"
 
@@ -39,6 +40,8 @@ public:
         } else if (a_model_type == MoistureType::Morrison ||
                    a_model_type == MoistureType::Morrison_NoIce) {
             SetModel<Morrison>();
+        } else if (a_model_type == MoistureType::WSM6) {
+            SetModel<WSM6>();
         } else if (a_model_type == MoistureType::SatAdj) {
             SetModel<SatAdj>();
         } else if (a_model_type == MoistureType::None) {

Source/Microphysics/ERF_Microphysics.H

@@ -103,6 +103,7 @@ public:
              || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce)
              || (a_moisture_type == MoistureType::Morrison)
              || (a_moisture_type == MoistureType::Morrison_NoIce)
+             || (a_moisture_type == MoistureType::WSM6)
              || (a_moisture_type == MoistureType::Kessler)
              || (a_moisture_type == MoistureType::Kessler_NoRain)
              || (a_moisture_type == MoistureType::SatAdj)

Source/Microphysics/WSM6/ERF_AdvanceWSM6.cpp

@@ -0,0 +1,139 @@
+#include "ERF_WSM6.H"
+#include "ERF_WSM6_Fortran_Interface.H"
+#include <cmath>
+
+using namespace amrex;
+
+void
+WSM6::Advance(const Real& dt_advance,
+              const SolverChoice&)
+{
+    dt = dt_advance;
+
+#ifdef ERF_USE_WSM6_FORT
+    static bool wsm6_inited = false;
+
+    // Minimal phase-1 initialization for single-moment WSM6.
+    if (!wsm6_inited) {
+        constexpr double den0 = 1.28;                 // Standard dry-air density (kg/m^3)
+        constexpr double denr = static_cast<double>(rhoh2o);
+        constexpr double dens = static_cast<double>(rhos);
+        constexpr double cl = static_cast<double>(Cp_l);
+        constexpr double cpv = static_cast<double>(Cp_v);
+        constexpr int hail_opt = 0;                   // Graupel mode
+        mp_wsm6_init_c(den0, denr, dens, cl, cpv, hail_opt);
+        wsm6_inited = true;
+    }
+
+    constexpr double g = static_cast<double>(CONST_GRAV);
+    constexpr double cpd = static_cast<double>(Cp_d);
+    constexpr double cpv = static_cast<double>(Cp_v);
+    constexpr double rd = static_cast<double>(R_d);
+    constexpr double rv = static_cast<double>(R_v);
+    constexpr double t0c = 273.15;
+    constexpr double ep1 = static_cast<double>(R_v / R_d - one);
+    constexpr double ep2 = static_cast<double>(R_d / R_v);
+    constexpr double qmin = 1.0e-12;
+    constexpr double xls = static_cast<double>(lsub);
+    constexpr double xlv0 = static_cast<double>(lat_vap);
+    constexpr double xlf0 = static_cast<double>(lat_ice);
+    constexpr double den0 = 1.28;
+    constexpr double denr = static_cast<double>(rhoh2o);
+    constexpr double cliq = static_cast<double>(Cp_l);
+    constexpr double cice = 2106.0;
+    constexpr double psat = 610.78;
+
+    for (MFIter mfi(*mic_fab_vars[MicVar_WSM6::qv], TileNoZ()); mfi.isValid(); ++mfi) {
+        const Box box = mfi.tilebox();
+        const Box fab_box = mfi.fabbox();
+
+        auto const& t_arr = mic_fab_vars[MicVar_WSM6::tabs]->array(mfi);
+        auto const& qv_arr = mic_fab_vars[MicVar_WSM6::qv]->array(mfi);
+        auto const& qc_arr = mic_fab_vars[MicVar_WSM6::qc]->array(mfi);
+        auto const& qi_arr = mic_fab_vars[MicVar_WSM6::qi]->array(mfi);
+        auto const& qr_arr = mic_fab_vars[MicVar_WSM6::qr]->array(mfi);
+        auto const& qs_arr = mic_fab_vars[MicVar_WSM6::qs]->array(mfi);
+        auto const& qg_arr = mic_fab_vars[MicVar_WSM6::qg]->array(mfi);
+        auto const& den_arr = mic_fab_vars[MicVar_WSM6::rho]->array(mfi);
+        auto const& p_arr = mic_fab_vars[MicVar_WSM6::pres]->array(mfi);
+        auto rain_arr = mic_fab_vars[MicVar_WSM6::rain_accum]->array(mfi);
+        auto snow_arr = mic_fab_vars[MicVar_WSM6::snow_accum]->array(mfi);
+        auto graup_arr = mic_fab_vars[MicVar_WSM6::graup_accum]->array(mfi);
+
+        const int ilo = box.smallEnd(0);
+        const int ihi = box.bigEnd(0);
+        const int jlo = box.smallEnd(1);
+        const int jhi = box.bigEnd(1);
+        const int klo = box.smallEnd(2);
+        const int khi = box.bigEnd(2);
+
+        const int imlo = fab_box.smallEnd(0);
+        const int imhi = fab_box.bigEnd(0);
+        const int jmlo = fab_box.smallEnd(1);
+        const int jmhi = fab_box.bigEnd(1);
+        const int kmlo = fab_box.smallEnd(2);
+        const int kmhi = fab_box.bigEnd(2);
+
+        const Real dz_val = m_geom.CellSize(m_axis);
+        FArrayBox delz_fab(fab_box, 1);
+        auto const& delz_arr = delz_fab.array();
+        delz_fab.setVal(dz_val);
+
+        const Array4<const Real> z_arr = (m_z_phys_nd) ? m_z_phys_nd->const_array(mfi) : Array4<const Real> {};
+        ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+            delz_arr(i,j,k) = (z_arr) ? Real(0.25) * ( (z_arr(i  ,j  ,k+1) - z_arr(i  ,j  ,k))
+                                                     + (z_arr(i+1,j  ,k+1) - z_arr(i+1,j  ,k))
+                                                     + (z_arr(i  ,j+1,k+1) - z_arr(i  ,j+1,k))
+                                                     + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz_val;
+        });
+
+        Box box2d(fab_box);
+        box2d.makeSlab(2, 0);
+        FArrayBox rainncv_fab(box2d, 1);
+        FArrayBox sr_fab(box2d, 1);
+        FArrayBox snowncv_fab(box2d, 1);
+        FArrayBox graupelncv_fab(box2d, 1);
+        FArrayBox rainacc_fab(box2d, 1);
+        FArrayBox snowacc_fab(box2d, 1);
+        FArrayBox graupacc_fab(box2d, 1);
+
+        auto const& rainncv_arr = rainncv_fab.array();
+        auto const& sr_arr = sr_fab.array();
+        auto const& snowncv_arr = snowncv_fab.array();
+        auto const& graupelncv_arr = graupelncv_fab.array();
+        auto const& rainacc_arr = rainacc_fab.array();
+        auto const& snowacc_arr = snowacc_fab.array();
+        auto const& graupacc_arr = graupacc_fab.array();
+        ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
+            rainacc_arr(i,j,0) = rain_arr(i,j,klo);
+            snowacc_arr(i,j,0) = snow_arr(i,j,klo);
+            graupacc_arr(i,j,0) = graup_arr(i,j,klo);
+            rainncv_arr(i,j,0) = Real(0.0);
+            sr_arr(i,j,0) = Real(0.0);
+            snowncv_arr(i,j,0) = Real(0.0);
+            graupelncv_arr(i,j,0) = Real(0.0);
+        });
+
+        mp_wsm6_run_c(
+            t_arr.dataPtr(),
+            qv_arr.dataPtr(), qc_arr.dataPtr(), qi_arr.dataPtr(),
+            qr_arr.dataPtr(), qs_arr.dataPtr(), qg_arr.dataPtr(),
+            den_arr.dataPtr(), p_arr.dataPtr(), delz_arr.dataPtr(),
+            static_cast<double>(dt), g, cpd, cpv, rd, rv, t0c, ep1, ep2, qmin,
+            xls, xlv0, xlf0, den0, denr, cliq, cice, psat,
+            rainacc_arr.dataPtr(), rainncv_arr.dataPtr(), sr_arr.dataPtr(),
+            snowacc_arr.dataPtr(), snowncv_arr.dataPtr(),
+            graupacc_arr.dataPtr(), graupelncv_arr.dataPtr(),
+            imlo, imhi, jmlo, jmhi, kmlo, kmhi,
+            ilo, ihi, jlo, jhi, klo, khi);
+
+        ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
+            rain_arr(i,j,klo) = rainacc_arr(i,j,0);
+            snow_arr(i,j,klo) = snowacc_arr(i,j,0);
+            graup_arr(i,j,klo) = graupacc_arr(i,j,0);
+        });
+    }
+#else
+    amrex::Abort("WSM6 Fortran bridge requested but ERF was not built with ERF_USE_WSM6_FORT");
+#endif
+}