MAM4xx Dry Deposition process interface for EAMxx

cime_config/tests.py

@@ -725,6 +725,7 @@
             "SMS_D_Ln5.ne4pg2_oQU480.F2010-SCREAMv1-MPASSI.scream-mam4xx-optics",
             "SMS_D_Ln5.ne4pg2_oQU480.F2010-SCREAMv1-MPASSI.scream-mam4xx-aci",
             "SMS_D_Ln5.ne4pg2_oQU480.F2010-SCREAMv1-MPASSI.scream-mam4xx-wetscav",
+            "SMS_D_Ln5.ne4pg2_oQU480.F2010-SCREAMv1-MPASSI.scream-mam4xx-drydep",
         )
     },
 

components/eamxx/cime_config/namelist_defaults_scream.xml

@@ -241,6 +241,10 @@ be lost if SCREAM_HACK_XML is not enabled.
       <top_level_mam4xx type="integer" doc="level corresponding to the top of troposphere clouds" nlev="128" >0</top_level_mam4xx>
     </mam4_aci>
 
+    <!-- MAM4xx-Dry Deposition -->
+    <mam4_drydep inherit="atm_proc_base">
+    </mam4_drydep>
+
     <!-- CLD fraction -->
     <cldFraction inherit="atm_proc_base"/>
 
@@ -493,7 +497,9 @@ be lost if SCREAM_HACK_XML is not enabled.
     <surf_lw_flux_up   >0.0</surf_lw_flux_up>
     <surf_mom_flux     type="array(real)">0.0,0.0</surf_mom_flux>
     <!-- Note: MAM4xx-ACI needs dry aerosol diameter (dgnum), currently set to aitken model nominal diameter, remove it once all of MAM4xx is integrated -->
-    <dgnum type="real" doc="Dry aerosol particles diameter in meters">2.6e-08</dgnum>
+    <dgnum type="real" doc="Dry aerosol particles diameter [m]">2.6e-08</dgnum>
+    <dgnumwet type="real" doc="Wet aerosol particles diameter [m]">0.41417721820867320E-007</dgnumwet>
+    <wetdens type="real" doc="Wet density of interstitial aerosol [kg/m3]">0.15100083211582764E+004</wetdens>
     <bc_c1 type="real">0.0</bc_c1>
     <bc_c3 type="real">0.0</bc_c3>
     <bc_c4 type="real">0.0</bc_c4>

components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/drydep/shell_commands

@@ -0,0 +1,17 @@
+
+#!/bin/sh
+#------------------------------------------------------
+# MAM4xx adds additionaltracers to the simulation
+# Increase number of tracers for MAM4xx simulations
+#------------------------------------------------------
+
+$CIMEROOT/../components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/update_eamxx_num_tracers.sh -b
+
+#------------------------------------------------------
+#Update IC file and add drydep process
+#------------------------------------------------------
+$CIMEROOT/../components/eamxx/scripts/atmchange initial_conditions::Filename='$DIN_LOC_ROOT/atm/scream/init/screami_mam4xx_ne4np4L72_c20240208.nc' -b
+$CIMEROOT/../components/eamxx/scripts/atmchange physics::atm_procs_list="mac_aero_mic,rrtmgp,mam4_drydep" -b
+
+
+

components/eamxx/src/control/atmosphere_surface_coupling_importer.cpp

@@ -52,6 +52,11 @@ void SurfaceCouplingImporter::set_grids(const std::shared_ptr<const GridsManager
   add_field<Computed>("snow_depth_land",  scalar2d_layout, m,       grid_name);
   add_field<Computed>("ocnfrac",          scalar2d_layout, nondim,  grid_name);
   add_field<Computed>("landfrac",         scalar2d_layout, nondim,  grid_name);
+  add_field<Computed>("icefrac",          scalar2d_layout, nondim,  grid_name);
+  // Friction velocity [m/s]
+  add_field<Computed>("fv",               scalar2d_layout, m/s,     grid_name);
+  // Aerodynamical resistance
+  add_field<Computed>("ram1",             scalar2d_layout, s/m,     grid_name);
 }
 // =========================================================================================
   void SurfaceCouplingImporter::setup_surface_coupling_data(const SCDataManager &sc_data_manager)

components/eamxx/src/mct_coupling/scream_cpl_indices.F90

@@ -6,7 +6,7 @@ module scream_cpl_indices
   private
 
   ! Focus only on the ones that scream imports/exports (subsets of x2a and a2x)
-  integer, parameter, public :: num_scream_imports = 16
+  integer, parameter, public :: num_scream_imports = 19
   integer, parameter, public :: num_scream_exports = 17
   integer, public :: num_cpl_imports, num_cpl_exports, import_field_size, export_field_size
 
@@ -89,6 +89,9 @@ subroutine scream_set_cpl_indices (x2a, a2x)
     import_field_names(14) = 'surf_evap'
     import_field_names(15) = 'ocnfrac'
     import_field_names(16) = 'landfrac'
+    import_field_names(17) = 'icefrac'
+    import_field_names(18) = 'fv'
+    import_field_names(19) = 'ram1'
 
     ! CPL indices
     import_cpl_indices(1)  = mct_avect_indexra(x2a,'Sx_avsdr')
@@ -106,7 +109,10 @@ subroutine scream_set_cpl_indices (x2a, a2x)
     import_cpl_indices(13) = mct_avect_indexra(x2a,'Faxx_sen')
     import_cpl_indices(14) = mct_avect_indexra(x2a,'Faxx_evap')
     import_cpl_indices(15) = mct_avect_indexra(x2a,'Sf_ofrac')
-    import_cpl_indices(16) = mct_avect_indexra(x2a,'Sf_lfrac') 
+    import_cpl_indices(16) = mct_avect_indexra(x2a,'Sf_lfrac')
+    import_cpl_indices(17) = mct_avect_indexra(x2a,'Sf_ifrac') 
+    import_cpl_indices(19) = mct_avect_indexra(x2a,'Sl_fv')
+    import_cpl_indices(18) = mct_avect_indexra(x2a,'Sl_ram1')
 
     ! Vector components
     import_vector_components(11) = 0

components/eamxx/src/physics/mam/CMakeLists.txt

@@ -43,6 +43,7 @@ add_subdirectory(${EXTERNALS_SOURCE_DIR}/mam4xx ${CMAKE_BINARY_DIR}/externals/ma
 add_library(mam
   eamxx_mam_microphysics_process_interface.cpp
   eamxx_mam_optics_process_interface.cpp
+  eamxx_mam_dry_deposition_process_interface.cpp
   eamxx_mam_aci_process_interface.cpp
   eamxx_mam_wetscav_process_interface.cpp)
 target_compile_definitions(mam PUBLIC EAMXX_HAS_MAM)

components/eamxx/src/physics/mam/eamxx_mam_dry_deposition_functions.hpp

@@ -0,0 +1,205 @@
+#ifndef EAMXX_MAM_DRY_DEPOSITION_FUNCTIONS_HPP
+#define EAMXX_MAM_DRY_DEPOSITION_FUNCTIONS_HPP
+
+#include <ekat/kokkos/ekat_subview_utils.hpp>
+#include <mam4xx/aero_config.hpp>
+#include <mam4xx/convproc.hpp>
+#include <mam4xx/mam4.hpp>
+
+namespace scream {
+
+namespace {
+void compute_tendencies(
+    // inputs
+    const int ncol, const int nlev, const double dt,
+    const MAMDryDep::const_view_1d obklen,
+    const MAMDryDep::const_view_1d surfric,
+    const MAMDryDep::const_view_1d landfrac,
+    const MAMDryDep::const_view_1d icefrac,
+    const MAMDryDep::const_view_1d ocnfrac,
+    const MAMDryDep::const_view_1d friction_velocity,
+    const MAMDryDep::const_view_1d aerodynamical_resistance,
+    MAMDryDep::view_3d qtracers, MAMDryDep::view_2d fraction_landuse_,
+    const MAMDryDep::const_view_3d dgncur_awet_,
+    const MAMDryDep::const_view_3d wet_dens_,
+    const mam_coupling::DryAtmosphere dry_atm,
+    const mam_coupling::AerosolState dry_aero,
+
+    // input-outputs
+    MAMDryDep::view_3d qqcw_,
+
+    // outputs
+    MAMDryDep::view_3d ptend_q, MAMDryDep::view_2d aerdepdrycw,
+    MAMDryDep::view_2d aerdepdryis,
+
+    // work arrays
+    MAMDryDep::view_2d rho_, MAMDryDep::view_4d vlc_dry_,
+    MAMDryDep::view_3d vlc_trb_, MAMDryDep::view_4d vlc_grv_,
+    MAMDryDep::view_3d dqdt_tmp_) {
+  static constexpr int num_aero_modes = mam_coupling::num_aero_modes();
+  const auto policy =
+      ekat::ExeSpaceUtils<MAMDryDep::KT::ExeSpace>::get_default_team_policy(
+          ncol, nlev);
+
+  // Parallel loop over all the columns
+  Kokkos::parallel_for(
+      policy, KOKKOS_LAMBDA(const MAMDryDep::KT::MemberType &team) {
+        static constexpr int num_aero_species =
+            mam_coupling::num_aero_species();
+
+        const int icol = team.league_rank();
+        // Parallel loop over all the levels to populate qtracers array using
+        // dry_aero
+        Kokkos::parallel_for(
+            Kokkos::TeamVectorRange(team, nlev), [&](const int lev) {
+              for(int mode = 0; mode < num_aero_modes; ++mode) {
+                int icnst = mam4::ConvProc::numptrcw_amode(mode);
+                qtracers(icol, lev, icnst) =
+                    dry_aero.int_aero_nmr[mode](icol, lev);
+                for(int species = 0; species < num_aero_species; ++species) {
+                  icnst = mam4::ConvProc::lmassptrcw_amode(species, mode);
+                  if(-1 < icnst) {
+                    qtracers(icol, lev, icnst) =
+                        dry_aero.int_aero_mmr[mode][species](icol, lev);
+                  }
+                }
+              }
+            });  // parallel_for for nlevs
+        team.team_barrier();
+
+        // Create atm and progs objects
+        mam4::Atmosphere atm    = atmosphere_for_column(dry_atm, icol);
+        mam4::Prognostics progs = aerosols_for_column(dry_aero, icol);
+
+        // Extract column data (or 1d view) from 2d views of data
+        mam4::ConstColumnView dgncur_awet[num_aero_modes];
+        mam4::ConstColumnView wet_dens[num_aero_modes];
+
+        for(int i = 0; i < num_aero_modes; ++i) {
+          dgncur_awet[i] = ekat::subview(dgncur_awet_, icol, i);
+          wet_dens[i]    = ekat::subview(wet_dens_, icol, i);
+        }
+
+        mam4::ColumnView rho;
+        rho = ekat::subview(rho_, icol);
+
+        static constexpr int n_land_type = MAMDryDep::n_land_type;
+        Real fraction_landuse[n_land_type];
+        for(int i = 0; i < n_land_type; ++i) {
+          fraction_landuse[i] = fraction_landuse_(i, icol);
+        }
+
+        static constexpr int nmodes = mam4::AeroConfig::num_modes();
+        mam4::ColumnView vlc_dry[nmodes][MAMDryDep::aerosol_categories_];
+        mam4::ColumnView vlc_grv[nmodes][MAMDryDep::aerosol_categories_];
+        Real vlc_trb[nmodes][MAMDryDep::aerosol_categories_];
+
+        for(int i = 0; i < nmodes; ++i) {
+          for(int j = 0; j < MAMDryDep::aerosol_categories_; ++j) {
+            vlc_dry[i][j] = ekat::subview(vlc_dry_, i, j, icol);
+            vlc_trb[i][j] = vlc_trb_(i, j, icol);
+            vlc_grv[i][j] = ekat::subview(vlc_grv_, i, j, icol);
+          }
+        }
+        static constexpr int pcnst = mam4::aero_model::pcnst;
+        mam4::ColumnView qqcw[pcnst];
+        mam4::ColumnView dqdt_tmp[pcnst];
+        for(int i = 0; i < pcnst; ++i) {
+          qqcw[i]     = ekat::subview(qqcw_, i, icol);
+          dqdt_tmp[i] = ekat::subview(dqdt_tmp_, i, icol);
+        }
+        // Extract qqcw from Prognostics
+        Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev), [&](int kk) {
+          for(int m = 0; m < nmodes; ++m) {
+            qqcw[mam4::ConvProc::numptrcw_amode(m)][kk] = progs.n_mode_c[m][kk];
+            for(int a = 0; a < mam4::AeroConfig::num_aerosol_ids(); ++a)
+              if(-1 < mam4::ConvProc::lmassptrcw_amode(a, m))
+                qqcw[mam4::ConvProc::lmassptrcw_amode(a, m)][kk] =
+                    progs.q_aero_c[m][a][kk];
+          }
+        });  // parallel_for nlevs
+        team.team_barrier();
+        bool ptend_lq[pcnst];  // currently unused
+        mam4::aero_model_drydep(
+            // inputs
+            team, fraction_landuse, atm.temperature, atm.pressure,
+            atm.interface_pressure, atm.hydrostatic_dp,
+            ekat::subview(qtracers, icol), dgncur_awet, wet_dens, obklen[icol],
+            surfric[icol], landfrac[icol], icefrac[icol], ocnfrac[icol],
+            friction_velocity[icol], aerodynamical_resistance[icol], dt,
+            // input-outputs
+            qqcw,
+            // outputs
+            ekat::subview(ptend_q, icol), ptend_lq,
+            ekat::subview(aerdepdrycw, icol), ekat::subview(aerdepdryis, icol),
+            // work arrays
+            rho, vlc_dry, vlc_trb, vlc_grv, dqdt_tmp);
+      });  // parallel_for for ncols
+}  // Compute_tendencies ends
+
+// Update interstitial aerosols using ptend_q tendencies
+void update_interstitial_mmrs(const MAMDryDep::view_3d ptend_q, const double dt,
+                              const int ncol, const int nlev,
+                              // output
+                              const mam_coupling::AerosolState dry_aero) {
+  const auto policy =
+      ekat::ExeSpaceUtils<MAMDryDep::KT::ExeSpace>::get_default_team_policy(
+          ncol, nlev);
+  static constexpr int nmodes = mam4::AeroConfig::num_modes();
+  Kokkos::parallel_for(
+      policy, KOKKOS_LAMBDA(const MAMDryDep::KT::MemberType &team) {
+        const int icol = team.league_rank();
+        Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev), [&](int kk) {
+          for(int m = 0; m < nmodes; ++m) {
+            dry_aero.int_aero_nmr[m](icol, kk) +=
+                ptend_q(icol, kk, mam4::ConvProc::numptrcw_amode(m)) * dt;
+            for(int a = 0; a < mam4::AeroConfig::num_aerosol_ids(); ++a)
+              if(-1 < mam4::ConvProc::lmassptrcw_amode(a, m))
+                dry_aero.int_aero_mmr[m][a](icol, kk) +=
+                    ptend_q(icol, kk, mam4::ConvProc::lmassptrcw_amode(a, m)) *
+                    dt;
+          }
+        });  // parallel_for nlevs
+      });    // parallel_for icol
+}  // Update interstitial aerosols ends
+
+// Update cloud borne aerosols using qqcw
+void update_cloudborne_mmrs(const MAMDryDep::view_3d qqcw, const double dt,
+                            const int nlev_,
+                            // output
+                            const mam_coupling::AerosolState dry_aero) {
+  for(int m = 0; m < mam_coupling::num_aero_modes(); ++m) {
+    Kokkos::deep_copy(dry_aero.cld_aero_nmr[m],
+                      ekat::subview(qqcw, mam4::ConvProc::numptrcw_amode(m)));
+    for(int a = 0; a < mam_coupling::num_aero_species(); ++a) {
+      if(dry_aero.cld_aero_mmr[m][a].data()) {
+        Kokkos::deep_copy(
+            dry_aero.cld_aero_mmr[m][a],
+            ekat::subview(qqcw, mam4::ConvProc::lmassptrcw_amode(a, m)));
+      }
+    }
+  }
+}  // Update cloud borne aerosols ends
+
+// FIXME: remove the following function after implementing file read for landuse
+void populated_fraction_landuse(MAMDryDep::view_2d flu, const int ncol) {
+  Kokkos::parallel_for(
+      "populated_fraction_landuse", 1, KOKKOS_LAMBDA(int) {
+        static constexpr int n_land_type = MAMDryDep::n_land_type;
+        const Real temp[n_land_type]     = {
+                0.28044346587077795E-003, 0.26634987180780171E-001,
+                0.16803558403621365E-001, 0.18076055155371872E-001,
+                0.00000000000000000E+000, 0.00000000000000000E+000,
+                0.91803784897907303E+000, 0.17186036997038400E-002,
+                0.00000000000000000E+000, 0.00000000000000000E+000,
+                0.18448503115578840E-001};
+        for(int i = 0; i < n_land_type; ++i)
+          for(int j = 0; j < ncol; ++j) flu(i, j) = temp[i];
+      });
+  Kokkos::fence();
+}
+
+}  // namespace
+}  // namespace scream
+
+#endif