Major MMF2/PAM updates

cime_config/tests.py

@@ -370,7 +370,7 @@
             "ERP_Ln9.ne4pg2_oQU480.WCYCL20TRNS-MMF1.allactive-mmf_fixed_subcycle",
             "ERS_Ln9.ne4pg2_ne4pg2.FRCE-MMF1.eam-cosp_nhtfrq9",
             "SMS_Ln5.ne4_ne4.FSCM-ARM97-MMF1",
-            # "SMS_Ln3.ne4pg2_ne4pg2.F2010-MMF2", - temporarily disabled due to ongoing dev issues
+            "SMS_Ln3.ne4pg2_oQU480.F2010-MMF2",
             )
         },
 

components/eam/cime_config/config_component.xml

@@ -71,14 +71,15 @@
 
       <!-- Multiscale modeling framework (MMF) -->
       <value compset="_EAM%.*MMF1"     >-crm samxx -crm_dt 10 </value>
-      <value compset="_EAM%.*MMF2_"    >-crm pam -pam_dycor spam -crm_dt 8 </value>
-      <value compset="_EAM%.*MMF2-AWFL">-crm pam -pam_dycor awfl -crm_dt 8 </value>
+      <value compset="_EAM%.*MMF2_"    >-crm pam -pam_dycor spam -crm_dt 10 </value>
+      <value compset="_EAM%.*MMF2-AWFL">-crm pam -pam_dycor awfl -crm_dt 10 </value>
       <value compset="_EAM%.*MMF"      >-use_MMF -nlev 60 -crm_nz 50</value>
-      <value compset="_EAM%.*MMF"      >-crm_dx 2000 -crm_nx 64 -crm_ny 1 </value>
+      <value compset="_EAM%.*MMF1"     >-crm_dx 2000 -crm_nx 64 -crm_ny 1 -crm_nx_rad 4 -crm_ny_rad 1</value>
+      <value compset="_EAM%.*MMF2"     >-crm_dx 3000 -crm_nx 45 -crm_ny 1 -crm_nx_rad 5 -crm_ny_rad 1</value>
       <value compset="_EAM%.*MMF1"     >-MMF_microphysics_scheme sam1mom -chem none</value>
       <value compset="_EAM%.*MMF2"     >-MMF_microphysics_scheme p3 -chem none</value>
-      <value compset="_EAM%.*MMF"      >-crm_nx_rad 4 -crm_ny_rad 1 -rad rrtmgp -rrtmgpxx</value>
-      <value compset="_EAM%.*MMF1"     >-use_MMF_VT</value>   <!-- CRM variance transport -->
+      <value compset="_EAM%.*MMF"      >-rad rrtmgp -rrtmgpxx</value>
+      <value compset="_EAM%.*MMF"      >-use_MMF_VT</value>   <!-- CRM variance transport -->
       <value compset="_EAM%.*MMF1"     >-use_MMF_ESMT</value> <!-- explicit scalar momentum transport -->
       <value compset="_EAM%AQP-MMF.*"  >-aquaplanet</value>
       <value compset="_EAM%RCE-MMF.*"  >-aquaplanet -rce</value>

components/eam/src/physics/crm/crm_physics.F90

@@ -1379,8 +1379,8 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
       call pam_mirror_array_readwrite( 'output_ni_mean',         crm_output%ni_mean,         '' )
       call pam_mirror_array_readwrite( 'output_qm_mean',         crm_output%qm_mean,         '' )
       call pam_mirror_array_readwrite( 'output_bm_mean',         crm_output%bm_mean,         '' )
-      call pam_mirror_array_readwrite( 'output_rho_d_mean',         crm_output%rho_d_mean,         '' )
-      call pam_mirror_array_readwrite( 'output_rho_v_mean',         crm_output%rho_v_mean,         '' )
+      call pam_mirror_array_readwrite( 'output_rho_d_mean',      crm_output%rho_d_mean,         '' )
+      call pam_mirror_array_readwrite( 'output_rho_v_mean',      crm_output%rho_v_mean,         '' )
 
       call pam_mirror_array_readwrite( 'output_qt_ls',       crm_output%qt_ls,       '' )
       call pam_mirror_array_readwrite( 'output_t_ls',        crm_output%t_ls,        '' )

components/eam/src/physics/crm/pam/CMakeLists.txt

@@ -11,6 +11,9 @@ set(PAM_DRIVER_SRC
 add_library(pam_driver 
            ${PAM_DRIVER_SRC})
 
+set(SCREAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../../../../..)
+add_library(eamxx_physics INTERFACE ${SCREAM_HOME}/components/eamxx/src/physics/)
+
 if (${USE_CUDA})
   # PAM will be CUDA-linked with device symbols resolved at library creation
   set_target_properties(pam_driver
@@ -37,7 +40,6 @@ set(PAM_MICRO   p3)
 set(PAM_SGS     shoc)
 set(PAM_RAD     forced)
 set(PAM_SCREAM_USE_CXX True)
-set(SCREAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../../../../..)
 
 # removing this once the cime cmake_macros aren't using obsolete
 # Kokkos settings like KOKKOS_OPTIONS.
@@ -124,7 +126,7 @@ target_compile_options(pam_driver PUBLIC
 )
 
 if (PAM_SCREAM_USE_CXX)
-  target_link_libraries(pam_driver pam_core physics dynamics pam_scream_cxx_interfaces ekat p3 shoc physics_share scream_share)
+  target_link_libraries(pam_driver pam_core physics dynamics pam_scream_cxx_interfaces ekat p3 shoc physics_share scream_share eamxx_physics)
 else()
   target_link_libraries(pam_driver pam_core physics dynamics )
 endif()

components/eam/src/physics/crm/pam/pam_accelerate.h

@@ -17,16 +17,22 @@ inline void pam_accelerate_init( pam::PamCoupler &coupler ) {
   using yakl::c::parallel_for;
   using yakl::c::SimpleBounds;
   using yakl::atomicAdd;
+  //------------------------------------------------------------------------------------------------
+  // The ability to accelerate dry density is implemented, but tests suggested 
+  // that this was contributing to instabilities in PAM. Combined with the 
+  // uncertainty of how to handle dry density in the anelastic case, we will
+  // disable dry density mean-state acceleration by default.
+  coupler.set_option<bool>("crm_accel_rd",false);
+  //------------------------------------------------------------------------------------------------
   auto &dm_device   = coupler.get_data_manager_device_readwrite();
   auto nens         = coupler.get_option<int>("ncrms");
   auto nz           = coupler.get_option<int>("crm_nz");
   auto ny           = coupler.get_option<int>("crm_ny");
   auto nx           = coupler.get_option<int>("crm_nx");
-  auto crm_accel_uv = coupler.get_option<bool>("crm_accel_uv");
   //------------------------------------------------------------------------------------------------
   dm_device.register_and_allocate<real>("accel_save_t", "saved temperature for MSA", {nz,nens}, {"z","nens"} );
-  dm_device.register_and_allocate<real>("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} );
   dm_device.register_and_allocate<real>("accel_save_q", "saved total water for MSA", {nz,nens}, {"z","nens"} );
+  dm_device.register_and_allocate<real>("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} );
   dm_device.register_and_allocate<real>("accel_save_u", "saved uvel for MSA",        {nz,nens}, {"z","nens"} );
   dm_device.register_and_allocate<real>("accel_save_v", "saved vvel for MSA",        {nz,nens}, {"z","nens"} );
   //------------------------------------------------------------------------------------------------
@@ -42,40 +48,37 @@ inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) {
   auto nz         = coupler.get_option<int>("crm_nz");
   auto ny         = coupler.get_option<int>("crm_ny");
   auto nx         = coupler.get_option<int>("crm_nx");
+  auto crm_accel_rd = coupler.get_option<bool>("crm_accel_rd");
   auto crm_accel_uv = coupler.get_option<bool>("crm_accel_uv");
   //------------------------------------------------------------------------------------------------
   auto temp = dm_device.get<real,4>("temp"       );
-  auto rhod = dm_device.get<real,4>("density_dry");
   auto rhov = dm_device.get<real,4>("water_vapor");
   auto rhol = dm_device.get<real,4>("cloud_water");
   auto rhoi = dm_device.get<real,4>("ice"        );
+  auto rhod = dm_device.get<real,4>("density_dry");
   auto uvel = dm_device.get<real,4>("uvel"       );
   auto vvel = dm_device.get<real,4>("vvel"       );
   auto accel_save_t = dm_device.get<real,2>("accel_save_t");
-  auto accel_save_r = dm_device.get<real,2>("accel_save_r");
   auto accel_save_q = dm_device.get<real,2>("accel_save_q");
+  auto accel_save_r = dm_device.get<real,2>("accel_save_r");
   auto accel_save_u = dm_device.get<real,2>("accel_save_u");
   auto accel_save_v = dm_device.get<real,2>("accel_save_v");
   //------------------------------------------------------------------------------------------------
   // compute horizontal means needed later for mean-state acceleration
   parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
     accel_save_t(k,n) = 0.0;
-    accel_save_r(k,n) = 0.0;
     accel_save_q(k,n) = 0.0;
-    if (crm_accel_uv) {
-      accel_save_u(k,n) = 0.0;
-      accel_save_v(k,n) = 0.0;
-    }
+    if (crm_accel_rd) { accel_save_r(k,n) = 0.0; }
+    if (crm_accel_uv) { accel_save_u(k,n) = 0.0; }
+    if (crm_accel_uv) { accel_save_v(k,n) = 0.0; }
   });
   real r_nx_ny  = 1._fp/(nx*ny);  // precompute reciprocal to avoid costly divisions
   parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
     yakl::atomicAdd( accel_save_t(k,n), temp(k,j,i,n) * r_nx_ny );
-    yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny );
     yakl::atomicAdd( accel_save_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny );
-    if (crm_accel_uv) {
-      yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny );
-      yakl::atomicAdd( accel_save_v(k,n), vvel(k,j,i,n) * r_nx_ny );
-    }
+    if (crm_accel_rd) { yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny ); }
+    if (crm_accel_uv) { yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny ); }
+    if (crm_accel_uv) { yakl::atomicAdd( accel_save_v(k,n), vvel(k,j,i,n) * r_nx_ny ); }
   });
   //------------------------------------------------------------------------------------------------
 }
@@ -86,28 +89,28 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) {
   using yakl::c::SimpleBounds;
   using yakl::atomicAdd;
   using yakl::ScalarLiveOut;
-  auto &dm_device = coupler.get_data_manager_device_readwrite();
-  auto nens       = coupler.get_option<int>("ncrms");
-  auto nz         = coupler.get_option<int>("crm_nz");
-  auto ny         = coupler.get_option<int>("crm_ny");
-  auto nx         = coupler.get_option<int>("crm_nx");
+  auto &dm_device       = coupler.get_data_manager_device_readwrite();
+  auto nens             = coupler.get_option<int>("ncrms");
+  auto nz               = coupler.get_option<int>("crm_nz");
+  auto ny               = coupler.get_option<int>("crm_ny");
+  auto nx               = coupler.get_option<int>("crm_nx");
+  auto crm_accel_rd     = coupler.get_option<bool>("crm_accel_rd");
+  auto crm_accel_uv     = coupler.get_option<bool>("crm_accel_uv");
+  real crm_accel_factor = coupler.get_option<real>("crm_accel_factor");
   //------------------------------------------------------------------------------------------------
   auto temp = dm_device.get<real,4>("temp"       );
-  auto rhod = dm_device.get<real,4>("density_dry");
   auto rhov = dm_device.get<real,4>("water_vapor");
   auto rhol = dm_device.get<real,4>("cloud_water");
   auto rhoi = dm_device.get<real,4>("ice"        );
+  auto rhod = dm_device.get<real,4>("density_dry");
   auto uvel = dm_device.get<real,4>("uvel"       );
   auto vvel = dm_device.get<real,4>("vvel"       );
   auto accel_save_t = dm_device.get<real,2>("accel_save_t");
-  auto accel_save_r = dm_device.get<real,2>("accel_save_r");
   auto accel_save_q = dm_device.get<real,2>("accel_save_q");
+  auto accel_save_r = dm_device.get<real,2>("accel_save_r");
   auto accel_save_u = dm_device.get<real,2>("accel_save_u");
   auto accel_save_v = dm_device.get<real,2>("accel_save_v");
   //------------------------------------------------------------------------------------------------
-  bool crm_accel_uv     = coupler.get_option<bool>("crm_accel_uv");
-  real crm_accel_factor = coupler.get_option<real>("crm_accel_factor");
-  //------------------------------------------------------------------------------------------------
   real2d hmean_t  ("hmean_t",   nz,nens);
   real2d hmean_r  ("hmean_r",   nz,nens);
   real2d hmean_q  ("hmean_q",   nz,nens);
@@ -127,35 +130,29 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) {
   // Compute the horizontal mean for each variable
   parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
     hmean_t(k,n) = 0.0;
-    hmean_r(k,n) = 0.0;
     hmean_q(k,n) = 0.0;
-    if (crm_accel_uv) {
-      hmean_u(k,n) = 0.0;
-      hmean_v(k,n) = 0.0;
-    }
+    if (crm_accel_rd) { hmean_r(k,n) = 0.0; }
+    if (crm_accel_uv) { hmean_u(k,n) = 0.0; }
+    if (crm_accel_uv) { hmean_v(k,n) = 0.0; }
   });
   real r_nx_ny  = 1._fp/(nx*ny);  // precompute reciprocal to avoid costly divisions
   parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
     yakl::atomicAdd( hmean_t(k,n), temp(k,j,i,n) * r_nx_ny );
-    yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny );
     yakl::atomicAdd( hmean_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny );
-    if (crm_accel_uv) {
-      yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny );
-      yakl::atomicAdd( hmean_v(k,n), vvel(k,j,i,n) * r_nx_ny );
-    }
+    if (crm_accel_rd) { yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny ); }
+    if (crm_accel_uv) { yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny ); }
+    if (crm_accel_uv) { yakl::atomicAdd( hmean_v(k,n), vvel(k,j,i,n) * r_nx_ny ); }
   });
   //------------------------------------------------------------------------------------------------
   // Compute the accelerated tendencies
   // NOTE - these are tendencies multiplied by the time step (i.e. d/dt * crm_dt)
   ScalarLiveOut<bool> ceaseflag_liveout(false);
   parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
     ttend_acc(k,n) = hmean_t(k,n) - accel_save_t(k,n);
-    rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n);
     qtend_acc(k,n) = hmean_q(k,n) - accel_save_q(k,n);
-    if (crm_accel_uv) {
-      utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n);
-      vtend_acc(k,n) = hmean_v(k,n) - accel_save_v(k,n);
-    }
+    if (crm_accel_rd) { rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n); }
+    if (crm_accel_uv) { utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n); }
+    if (crm_accel_uv) { vtend_acc(k,n) = hmean_v(k,n) - accel_save_v(k,n); }
     if (abs(ttend_acc(k,n)) > dtemp_max) {
       ceaseflag_liveout = true;
     }
@@ -191,14 +188,11 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) {
   //------------------------------------------------------------------------------------------------
   // Apply the accelerated tendencies
   parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
-
     temp(k,j,i,n) = temp(k,j,i,n) + crm_accel_factor * ttend_acc(k,n);
-    rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n);
     rhov(k,j,i,n) = rhov(k,j,i,n) + crm_accel_factor * qtend_acc(k,n);
-    if (crm_accel_uv) {
-      uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n); 
-      vvel(k,j,i,n) = vvel(k,j,i,n) + crm_accel_factor * vtend_acc(k,n); 
-    }
+    if (crm_accel_rd) { rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n); }
+    if (crm_accel_uv) { uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n); }
+    if (crm_accel_uv) { vvel(k,j,i,n) = vvel(k,j,i,n) + crm_accel_factor * vtend_acc(k,n); }
     // apply limiter on temperature
     temp(k,j,i,n) = std::max( temp_min, temp(k,j,i,n) );
   });
@@ -231,4 +225,3 @@ inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) {
   });
   //------------------------------------------------------------------------------------------------
 }
-