diff --git a/components/eamxx/cime_config/namelist_defaults_scream.xml b/components/eamxx/cime_config/namelist_defaults_scream.xml index f09d10b0b32..dadb6dbba35 100644 --- a/components/eamxx/cime_config/namelist_defaults_scream.xml +++ b/components/eamxx/cime_config/namelist_defaults_scream.xml @@ -242,8 +242,7 @@ be lost if SCREAM_HACK_XML is not enabled. - - + @@ -266,8 +265,41 @@ be lost if SCREAM_HACK_XML is not enabled. - - + + + + + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/DMSflux.2010.ne4pg2_conserv.POPmonthlyClimFromACES4BGC_c20240814.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_so2_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_bc_a4_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_num_a1_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_num_a2_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_num_a4_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_pom_a4_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_so4_a1_surf_ne4pg2_2010_clim_c20240815.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/surface/cmip6_mam4_so4_a2_surf_ne4pg2_2010_clim_c20240815.nc + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/DMSflux.2010.ne30pg2_conserv.POPmonthlyClimFromACES4BGC_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_so2_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_bc_a4_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_num_a1_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_num_a2_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_num_a4_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_pom_a4_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_so4_a1_surf_ne30pg2_2010_clim_c20240816.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/surface/cmip6_mam4_so4_a2_surf_ne30pg2_2010_clim_c20240816.nc + + + + ${DIN_LOC_ROOT}/atm/scream/maps/map_ne30pg2_to_ne120pg2_20231201.nc + ${DIN_LOC_ROOT}/atm/scream/maps/map_ne30pg2_to_ne256pg2_20231201.nc + ${DIN_LOC_ROOT}/atm/scream/maps/map_ne30pg2_to_ne512pg2_20231201.nc + ${DIN_LOC_ROOT}/atm/scream/maps/map_ne30pg2_to_ne1024pg2_20231201.nc + + + + diff --git a/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/srf_online_emiss_constituent_fluxes/shell_commands b/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/srf_online_emiss_constituent_fluxes/shell_commands new file mode 100644 index 00000000000..f47a1729ed7 --- /dev/null +++ b/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/srf_online_emiss_constituent_fluxes/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 the processes +#------------------------------------------------------ +$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="mam4_constituent_fluxes,mac_aero_mic,rrtmgp,mam4_srf_online_emiss" -b + + + diff --git a/components/eamxx/src/physics/mam/CMakeLists.txt b/components/eamxx/src/physics/mam/CMakeLists.txt index c5feb52e12a..9874f79f9eb 100644 --- a/components/eamxx/src/physics/mam/CMakeLists.txt +++ b/components/eamxx/src/physics/mam/CMakeLists.txt @@ -45,7 +45,9 @@ add_library(mam 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) + eamxx_mam_wetscav_process_interface.cpp + eamxx_mam_srf_and_online_emissions_process_interface.cpp + eamxx_mam_constituent_fluxes_interface.cpp) target_compile_definitions(mam PUBLIC EAMXX_HAS_MAM) add_dependencies(mam mam4xx) target_include_directories(mam PUBLIC diff --git a/components/eamxx/src/physics/mam/eamxx_mam_aci_functions.hpp b/components/eamxx/src/physics/mam/eamxx_mam_aci_functions.hpp index 03b04ac739a..8bb52c918d2 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_aci_functions.hpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_aci_functions.hpp @@ -198,23 +198,6 @@ void store_liquid_cloud_fraction( }); } -void compute_recipical_pseudo_density(haero::ThreadTeamPolicy team_policy, - MAMAci::const_view_2d pdel, - const int nlev, - // output - MAMAci::view_2d rpdel) { - Kokkos::parallel_for( - team_policy, KOKKOS_LAMBDA(const haero::ThreadTeam &team) { - const int icol = team.league_rank(); - Kokkos::parallel_for( - Kokkos::TeamVectorRange(team, 0, nlev), [&](int kk) { - EKAT_KERNEL_ASSERT_MSG(0 < pdel(icol, kk), - "Error: pdel should be > 0.\n"); - rpdel(icol, kk) = 1 / pdel(icol, kk); - }); - }); -} - void call_function_dropmixnuc( haero::ThreadTeamPolicy team_policy, const Real dt, mam_coupling::DryAtmosphere &dry_atmosphere, const MAMAci::view_2d rpdel, @@ -397,7 +380,7 @@ void call_function_dropmixnuc( progs_at_col, haero_atm, state_q_at_lev_col, klev); // get the start index for aerosols species in the state_q array - int istart = mam4::aero_model::pcnst - mam4::ndrop::ncnst_tot; + int istart = mam4::utils::aero_start_ind(); // create colum views of state_q for(int icnst = istart; icnst < mam4::aero_model::pcnst; diff --git a/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.cpp b/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.cpp index e919aa066b4..ef2c0f2b288 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.cpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.cpp @@ -253,7 +253,7 @@ void MAMAci::set_grids( frz_unit, grid_name); // heterogeneous freezing by deposition nucleation [cm^-3 s^-1] - add_field("hetfrz_depostion_nucleation_tend", scalar3d_layout_mid, + add_field("hetfrz_deposition_nucleation_tend", scalar3d_layout_mid, frz_unit, grid_name); } // function set_grids ends @@ -390,8 +390,8 @@ void MAMAci::initialize_impl(const RunType run_type) { get_field_out("hetfrz_immersion_nucleation_tend").get_view(); hetfrz_contact_nucleation_tend_ = get_field_out("hetfrz_contact_nucleation_tend").get_view(); - hetfrz_depostion_nucleation_tend_ = - get_field_out("hetfrz_depostion_nucleation_tend").get_view(); + hetfrz_deposition_nucleation_tend_ = + get_field_out("hetfrz_deposition_nucleation_tend").get_view(); //--------------------------------------------------------------------------------- // Allocate memory for the class members @@ -597,7 +597,7 @@ void MAMAci::run_impl(const double dt) { // output cloud_frac_, cloud_frac_prev_); - compute_recipical_pseudo_density(team_policy, dry_atm_.p_del, nlev_, + mam_coupling::compute_recipical_pseudo_density(team_policy, dry_atm_.p_del, nlev_, // output rpdel_); @@ -642,7 +642,7 @@ void MAMAci::run_impl(const double dt) { team_policy, hetfrz_, dry_atm_, dry_aero_, factnum_, dt, nlev_, // ## output to be used by the other processes ## hetfrz_immersion_nucleation_tend_, hetfrz_contact_nucleation_tend_, - hetfrz_depostion_nucleation_tend_, + hetfrz_deposition_nucleation_tend_, // work arrays diagnostic_scratch_); diff --git a/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.hpp b/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.hpp index fc930ea11d7..63f87dd9f24 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.hpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_aci_process_interface.hpp @@ -123,7 +123,7 @@ class MAMAci final : public scream::AtmosphereProcess { // added correctly to the cloud-micorphysics scheme. view_2d hetfrz_immersion_nucleation_tend_; view_2d hetfrz_contact_nucleation_tend_; - view_2d hetfrz_depostion_nucleation_tend_; + view_2d hetfrz_deposition_nucleation_tend_; view_2d diagnostic_scratch_[hetro_scratch_]; diff --git a/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_functions.hpp b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_functions.hpp new file mode 100644 index 00000000000..602fc9ebf89 --- /dev/null +++ b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_functions.hpp @@ -0,0 +1,84 @@ +#ifndef EAMXX_MAM_CONSTITUTE_FLUXES_FUNCTIONS_HPP +#define EAMXX_MAM_CONSTITUTE_FLUXES_FUNCTIONS_HPP + +#include +#include + +namespace scream { + +namespace { + +void update_gas_aerosols_using_constituents( + const int ncol, const int nlev, const double dt, + const mam_coupling::DryAtmosphere &dry_atm, + const MAMConstituentFluxes::const_view_2d &constituent_fluxes, + // output + const mam_coupling::AerosolState &wet_aero) { + using C = physics::Constants; + static constexpr auto gravit = C::gravit; // Gravity [m/s2] + static constexpr int pcnst = mam4::aero_model::pcnst; + + // Declare local variables + const int surface_lev = nlev - 1; + + // get the start index for gas species in the state_q array + int istart = mam4::utils::gasses_start_ind(); + + // number of constituents to update (currently updating only MAM4xx + // constituents) + const int nconstituents = pcnst - istart; + + // Create a policy to loop over columns annd number of constituents + // to update + // FIXME: TODO:We don't need a team for "nconstituents", so we can make the + // kookos_for simple by using just ncols + const auto policy = ekat::ExeSpaceUtils:: + get_default_team_policy(ncol, nconstituents); + + // Loop through all columns to update tracer mixing rations + Kokkos::parallel_for( + policy, KOKKOS_LAMBDA(const haero::ThreadTeam &team) { + const int icol = team.league_rank(); + + //---------------------------------------------------------------------- + // To form EAM like state%q array, we need prognostics (gas and aerosol + // mmrs) atmosphere (qv, qc, nc, ni, etc.) + //---------------------------------------------------------------------- + + // Get prognostics + mam4::Prognostics progs_at_col = + mam_coupling::aerosols_for_column(wet_aero, // output + icol); // input + // Get atmospheric quantities + const haero::Atmosphere haero_atm = + atmosphere_for_column(dry_atm, // output + icol); // input + + // Form state%q like array at surface level + Real state_q_at_surf_lev[pcnst] = {}; + mam4::utils::extract_stateq_from_prognostics( + progs_at_col, haero_atm, // input + state_q_at_surf_lev, // output + surface_lev); // input + + // Compute the units conversion factor (kg/m2/s to kg/kg) + EKAT_KERNEL_ASSERT_MSG(dry_atm.p_del(icol, surface_lev) != 0, + "Error! dry_atm.pdel must be non-zero!\n"); + const Real rpdel = 1.0 / dry_atm.p_del(icol, surface_lev); + const Real unit_factor = dt * gravit * rpdel; + + // Update state vector with constituent fluxes + for(int icnst = istart; icnst < pcnst; ++icnst) { + state_q_at_surf_lev[icnst] += + constituent_fluxes(icol, icnst) * unit_factor; + } + mam4::utils::inject_stateq_to_prognostics(state_q_at_surf_lev, // input + progs_at_col, // output + surface_lev); // input + }); // icol loop +} + +} // namespace +} // namespace scream + +#endif // ifdef EAMXX_MAM_CONSTITUTE_FLUXES_FUNCTIONS_HPP diff --git a/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.cpp b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.cpp new file mode 100644 index 00000000000..c168c1ef0e3 --- /dev/null +++ b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.cpp @@ -0,0 +1,306 @@ +#include +#include +namespace scream { + +// ================================================================ +// Constructor +// ================================================================ +MAMConstituentFluxes::MAMConstituentFluxes(const ekat::Comm &comm, + const ekat::ParameterList ¶ms) + : AtmosphereProcess(comm, params) { + /* Anything that can be initialized without grid information can be + * initialized here. Like universal constants, mam wetscav options. + */ +} + +// ================================================================ +// SET_GRIDS +// ================================================================ +void MAMConstituentFluxes::set_grids( + const std::shared_ptr grids_manager) { + grid_ = grids_manager->get_grid("Physics"); + const auto &grid_name = grid_->name(); + + ncol_ = grid_->get_num_local_dofs(); // Number of columns on this rank + nlev_ = grid_->get_num_vertical_levels(); // Number of levels per column + + using namespace ekat::units; + auto q_unit = kg / kg; // units of mass mixing ratios of tracers + auto n_unit = 1 / kg; // units of number mixing ratios of tracers + auto nondim = ekat::units::Units::nondimensional(); + + FieldLayout scalar2d = grid_->get_2d_scalar_layout(); + FieldLayout scalar3d_mid = grid_->get_3d_scalar_layout(true); + FieldLayout scalar3d_int = grid_->get_3d_scalar_layout(false); + + static constexpr int pcnst = mam4::aero_model::pcnst; + + const FieldLayout scalar2d_pcnct = + grid_->get_2d_vector_layout(pcnst, "num_phys_constituents"); + + // -------------------------------------------------------------------------- + // These variables are "Required" or pure inputs for the process + // -------------------------------------------------------------------------- + // ----------- Atmospheric quantities ------------- + // Specific humidity [kg/kg](Require only for building DS) + add_field("qv", scalar3d_mid, q_unit, grid_name, "tracers"); + + // Cloud liquid mass mixing ratio [kg/kg](Require only for building DS) + add_field("qc", scalar3d_mid, q_unit, grid_name, "tracers"); + + // Cloud ice mass mixing ratio [kg/kg](Require only for building DS) + add_field("qi", scalar3d_mid, q_unit, grid_name, "tracers"); + + // Cloud liquid number mixing ratio [1/kg](Require only for building DS) + add_field("nc", scalar3d_mid, n_unit, grid_name, "tracers"); + + // Cloud ice number mixing ratio [1/kg](Require only for building DS) + add_field("ni", scalar3d_mid, n_unit, grid_name, "tracers"); + + // Temperature[K] at midpoints + add_field("T_mid", scalar3d_mid, K, grid_name); + + // Vertical pressure velocity [Pa/s] at midpoints (Require only for building + // DS) + add_field("omega", scalar3d_mid, Pa / s, grid_name); + + // Total pressure [Pa] at midpoints + add_field("p_mid", scalar3d_mid, Pa, grid_name); + + // Total pressure [Pa] at interfaces + add_field("p_int", scalar3d_int, Pa, grid_name); + + // Layer thickness(pdel) [Pa] at midpoints + add_field("pseudo_density", scalar3d_mid, Pa, grid_name); + + // Planetary boundary layer height [m] (Require only for building DS) + add_field("pbl_height", scalar2d, m, grid_name); + + // cloud fraction [nondimensional] computed by eamxx_cld_fraction_process + add_field("cldfrac_tot", scalar3d_mid, nondim, grid_name); + + static constexpr Units m2(m * m, "m2"); + static constexpr Units s2(s * s, "s2"); + + // Surface geopotential [m2/s2] (Require only for building DS) + add_field("phis", scalar2d, m2 / s2, grid_name); + + // Constituent fluxes at the surface (gasses and aerosols) + //[units: kg/m2/s (mass) or #/m2/s (number)] + add_field("constituent_fluxes", scalar2d_pcnct, kg / m2 / s, + grid_name); + + // --------------------------------------------------------------------- + // These variables are "Updated" or inputs/outputs for the process + // --------------------------------------------------------------------- + // NOTE: Cloud borne aerosols are not updated in this process but are included + // to create data structures. + + // interstitial and cloudborne aerosol tracers of interest: mass (q) and + // number (n) mixing ratios + for(int mode = 0; mode < mam_coupling::num_aero_modes(); ++mode) { + // interstitial aerosol tracers of interest: number (n) mixing ratios + const std::string int_nmr_field_name = + mam_coupling::int_aero_nmr_field_name(mode); + add_field(int_nmr_field_name, scalar3d_mid, n_unit, grid_name, + "tracers"); + + // cloudborne aerosol tracers of interest: number (n) mixing ratios + // NOTE: DO NOT add cld borne aerosols to the "tracer" group as these are + // NOT advected + const std::string cld_nmr_field_name = + mam_coupling::cld_aero_nmr_field_name(mode); + add_field(cld_nmr_field_name, scalar3d_mid, n_unit, grid_name); + + for(int a = 0; a < mam_coupling::num_aero_species(); ++a) { + // (interstitial) aerosol tracers of interest: mass (q) mixing ratios + const std::string int_mmr_field_name = + mam_coupling::int_aero_mmr_field_name(mode, a); + if(not int_mmr_field_name.empty()) { + add_field(int_mmr_field_name, scalar3d_mid, q_unit, grid_name, + "tracers"); + } + // (cloudborne) aerosol tracers of interest: mass (q) mixing ratios + // NOTE: DO NOT add cld borne aerosols to the "tracer" group as these are + // NOT advected + const std::string cld_mmr_field_name = + mam_coupling::cld_aero_mmr_field_name(mode, a); + if(not cld_mmr_field_name.empty()) { + add_field(cld_mmr_field_name, scalar3d_mid, q_unit, grid_name); + } + } // end for loop num species + } // end for loop for num modes + + for(int g = 0; g < mam_coupling::num_aero_gases(); ++g) { + const std::string gas_mmr_field_name = mam_coupling::gas_mmr_field_name(g); + add_field(gas_mmr_field_name, scalar3d_mid, q_unit, grid_name, + "tracers"); + } // end for loop num gases + +} // set_grid + +// ================================================================ +// REQUEST_BUFFER_SIZE_IN_BYTES +// ================================================================ +// ON HOST, returns the number of bytes of device memory needed by the above +// Buffer type given the number of columns and vertical levels +size_t MAMConstituentFluxes::requested_buffer_size_in_bytes() const { + return mam_coupling::buffer_size(ncol_, nlev_); +} + +// ================================================================ +// INIT_BUFFERS +// ================================================================ +// ON HOST, initializes the Buffer type with sufficient memory to store +// intermediate (dry) quantities on the given number of columns with the given +// number of vertical levels. Returns the number of bytes allocated. +void MAMConstituentFluxes::init_buffers( + const ATMBufferManager &buffer_manager) { + EKAT_REQUIRE_MSG( + buffer_manager.allocated_bytes() >= requested_buffer_size_in_bytes(), + "Error! Insufficient buffer size.\n"); + + size_t used_mem = + mam_coupling::init_buffer(buffer_manager, ncol_, nlev_, buffer_); + EKAT_REQUIRE_MSG( + used_mem == requested_buffer_size_in_bytes(), + "Error! Used memory != requested memory for MAMConstituentFluxes."); +} + +// ================================================================ +// INITIALIZE_IMPL +// ================================================================ +void MAMConstituentFluxes::initialize_impl(const RunType run_type) { + // --------------------------------------------------------------- + // Input fields read in from IC file, namelist or other processes + // --------------------------------------------------------------- + + // Populate the wet atmosphere state with views from fields + // FIMXE: specifically look which among these are actually used by the process + wet_atm_.qv = get_field_in("qv").get_view(); + + // Following wet_atm vars are required only for building DS + wet_atm_.qc = get_field_in("qc").get_view(); + wet_atm_.nc = get_field_in("nc").get_view(); + wet_atm_.qi = get_field_in("qi").get_view(); + wet_atm_.ni = get_field_in("ni").get_view(); + + // Populate the dry atmosphere state with views from fields + dry_atm_.T_mid = get_field_in("T_mid").get_view(); + dry_atm_.p_mid = get_field_in("p_mid").get_view(); + dry_atm_.p_del = get_field_in("pseudo_density").get_view(); + dry_atm_.p_int = get_field_in("p_int").get_view(); + + // Following dry_atm vars are required only for building DS + dry_atm_.cldfrac = get_field_in("cldfrac_tot").get_view(); + dry_atm_.pblh = get_field_in("pbl_height").get_view(); + dry_atm_.omega = get_field_in("omega").get_view(); + dry_atm_.p_del = get_field_in("pseudo_density").get_view(); + + // store fields converted to dry mmr from wet mmr in dry_atm_ + dry_atm_.z_mid = buffer_.z_mid; + dry_atm_.z_iface = buffer_.z_iface; + dry_atm_.dz = buffer_.dz; + dry_atm_.qv = buffer_.qv_dry; + dry_atm_.qc = buffer_.qc_dry; + dry_atm_.nc = buffer_.nc_dry; + dry_atm_.qi = buffer_.qi_dry; + dry_atm_.ni = buffer_.ni_dry; + dry_atm_.w_updraft = buffer_.w_updraft; + dry_atm_.z_surf = 0.0; // FIXME: for now + + // Constituent fluxes at the surface (gasses and aerosols) [kg/m2/s] + constituent_fluxes_ = + get_field_in("constituent_fluxes").get_view(); + + // interstitial and cloudborne aerosol tracers of interest: mass (q) and + // number (n) mixing ratios + for(int m = 0; m < mam_coupling::num_aero_modes(); ++m) { + // interstitial aerosol tracers of interest: number (n) mixing ratios + const std::string int_nmr_field_name = + mam_coupling::int_aero_nmr_field_name(m); + wet_aero_.int_aero_nmr[m] = + get_field_out(int_nmr_field_name).get_view(); + + // cloudborne aerosol tracers of interest: number (n) mixing ratios + const std::string cld_nmr_field_name = + mam_coupling::cld_aero_nmr_field_name(m); + wet_aero_.cld_aero_nmr[m] = + get_field_out(cld_nmr_field_name).get_view(); + + for(int a = 0; a < mam_coupling::num_aero_species(); ++a) { + // (interstitial) aerosol tracers of interest: mass (q) mixing ratios + const std::string int_mmr_field_name = + mam_coupling::int_aero_mmr_field_name(m, a); + + if(not int_mmr_field_name.empty()) { + wet_aero_.int_aero_mmr[m][a] = + get_field_out(int_mmr_field_name).get_view(); + } + + // (cloudborne) aerosol tracers of interest: mass (q) mixing ratios + const std::string cld_mmr_field_name = + mam_coupling::cld_aero_mmr_field_name(m, a); + if(not cld_mmr_field_name.empty()) { + wet_aero_.cld_aero_mmr[m][a] = + get_field_out(cld_mmr_field_name).get_view(); + } + } + } + for(int g = 0; g < mam_coupling::num_aero_gases(); ++g) { + const std::string gas_mmr_field_name = mam_coupling::gas_mmr_field_name(g); + wet_aero_.gas_mmr[g] = + get_field_out(gas_mmr_field_name).get_view(); + } + +} // end initialize_impl() + +// ================================================================ +// RUN_IMPL +// ================================================================ +void MAMConstituentFluxes::run_impl(const double dt) { + // ------------------------------------------------------------------- + // (LONG) NOTE: The following code is an adaptation of cflx.F90 code in + // E3SM. In EAMxx, all constituents are considered "wet" (or have wet + // mixing ratios), we are *not* doing any wet to dry conversions in the + // for this process. We are simply updating the MAM4xx tracers using the + // "constituent fluxes". + // We are converting wet atm to dry atm. Since we do not use or update + // any of the water constituents (qc, qv, qi etc.), we should be okay + // to do this conversion. We need to do this conversion as our function + // are built following HAERO data structures. + // ------------------------------------------------------------------- + + // Compute vertical layer heights and updraft velocity. We need these to fully + // populate dry_atm_, so that we can form a HAERO atmosphere object. HAERO + // atmosphere object is used to for state%q like array. + auto lambda = + KOKKOS_LAMBDA(const Kokkos::TeamPolicy::member_type &team) { + const int icol = team.league_rank(); // column index + compute_dry_mixing_ratios(team, wet_atm_, // in + dry_atm_, // out + icol); // in + team.team_barrier(); + // vertical heights has to be computed after computing dry mixing ratios + // for atmosphere + compute_vertical_layer_heights(team, // in + dry_atm_, // out + icol); // in + compute_updraft_velocities(team, wet_atm_, // in + dry_atm_, // out + icol); // in + }; + // policy + const auto scan_policy = ekat::ExeSpaceUtils< + KT::ExeSpace>::get_thread_range_parallel_scan_team_policy(ncol_, nlev_); + + Kokkos::parallel_for("mam_cfi_compute_updraft", scan_policy, lambda); + + update_gas_aerosols_using_constituents(ncol_, nlev_, dt, dry_atm_, + constituent_fluxes_, + // output + wet_aero_); +} // run_impl ends + +// ============================================================================= +} // namespace scream diff --git a/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.hpp b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.hpp new file mode 100644 index 00000000000..7d54761a297 --- /dev/null +++ b/components/eamxx/src/physics/mam/eamxx_mam_constituent_fluxes_interface.hpp @@ -0,0 +1,76 @@ +#ifndef EAMXX_MAM_CONSTITUENT_FLUXES_HPP +#define EAMXX_MAM_CONSTITUENT_FLUXES_HPP + +// For declaring contituent fluxes class derived from atm process +// class +#include + +// For MAM4 aerosol configuration +#include +#include + +namespace scream { + +// The process responsible for applying MAM4 constituent fluxes. The +// AD stores exactly ONE instance of this class in its list of subcomponents. +class MAMConstituentFluxes final : public scream::AtmosphereProcess { + public: + using KT = ekat::KokkosTypes; + using const_view_2d = Field::view_dev_t; + + private: + // number of horizontal columns + int ncol_, nlev_; + + // Wet and dry states of atmosphere + mam_coupling::WetAtmosphere wet_atm_; + mam_coupling::DryAtmosphere dry_atm_; + + // aerosol state variables + mam_coupling::AerosolState wet_aero_; + + // buffer for sotring temporary variables + mam_coupling::Buffer buffer_; + + // physics grid for column information + std::shared_ptr grid_; + + const_view_2d constituent_fluxes_; + + public: + // Constructor + MAMConstituentFluxes(const ekat::Comm &comm, + const ekat::ParameterList ¶ms); + + // -------------------------------------------------------------------------- + // AtmosphereProcess overrides (see share/atm_process/atmosphere_process.hpp) + // -------------------------------------------------------------------------- + + // The type of subcomponent + AtmosphereProcessType type() const { return AtmosphereProcessType::Physics; } + + // The name of the subcomponent + std::string name() const { return "mam_constituent_fluxes"; } + + // grid + void set_grids( + const std::shared_ptr grids_manager) override; + + // management of common atm process memory + size_t requested_buffer_size_in_bytes() const override; + void init_buffers(const ATMBufferManager &buffer_manager) override; + + // Initialize variables + void initialize_impl(const RunType run_type) override; + + // Run the process by one time step + void run_impl(const double dt) override; + + // Finalize + void finalize_impl(){/*Do nothing*/}; + +}; // MAMConstituentFluxes + +} // namespace scream + +#endif // EAMXX_MAM_CONSTITUENT_FLUXES_HPP diff --git a/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.cpp b/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.cpp index 8b93936148c..f8a34bd90d5 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.cpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.cpp @@ -450,7 +450,7 @@ void MAMMicrophysics::run_impl(const double dt) { impl::compute_water_content(progs, k, qv, temp, pmid, dgncur_a, dgncur_awet, wetdens, qaerwat); // do aerosol microphysics (gas-aerosol exchange, nucleation, coagulation) - impl::modal_aero_amicphys_intr(config.amicphys, some_step, dt, t, pmid, pdel, + impl::modal_aero_amicphys_intr(config.amicphys, step_, dt, t, pmid, pdel, zm, pblh, qv, cldfrac, vmr, vmrcw, vmr_pregaschem, vmr_precldchem, vmrcw_precldchem, vmr_tendbb, vmrcw_tendbb, dgncur_a, dgncur_awet, diff --git a/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.cpp b/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.cpp new file mode 100644 index 00000000000..850a82d0896 --- /dev/null +++ b/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.cpp @@ -0,0 +1,275 @@ +//#include +#include + +namespace scream { + +// ================================================================ +// Constructor +// ================================================================ +MAMSrfOnlineEmiss::MAMSrfOnlineEmiss(const ekat::Comm &comm, + const ekat::ParameterList ¶ms) + : AtmosphereProcess(comm, params) { + /* Anything that can be initialized without grid information can be + * initialized here. Like universal constants. + */ +} + +// ================================================================ +// SET_GRIDS +// ================================================================ +void MAMSrfOnlineEmiss::set_grids( + const std::shared_ptr grids_manager) { + grid_ = grids_manager->get_grid("Physics"); + const auto &grid_name = grid_->name(); + + ncol_ = grid_->get_num_local_dofs(); // Number of columns on this rank + nlev_ = grid_->get_num_vertical_levels(); // Number of levels per column + + using namespace ekat::units; + + static constexpr int pcnst = mam4::aero_model::pcnst; + const FieldLayout scalar2d_pcnct = + grid_->get_2d_vector_layout(pcnst, "num_phys_constituents"); + + // ------------------------------------------------------------- + // These variables are "Computed" or outputs for the process + // ------------------------------------------------------------- + static constexpr Units m2(m * m, "m2"); + // Constituent fluxes of species in [kg/m2/s] + add_field("constituent_fluxes", scalar2d_pcnct, kg / m2 / s, + grid_name); + + // Surface emissions remapping file + auto srf_map_file = m_params.get("srf_remap_file", ""); + + // FIXME: We can extract the following info about each species + // in a separate hpp file + //-------------------------------------------------------------------- + // Init dms srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ dms; + // File name, name and sectors + dms.data_file = m_params.get("srf_emis_specifier_for_DMS"); + dms.species_name = "dms"; + dms.sectors = {"DMS"}; + srf_emiss_species_.push_back(dms); // add to the vector + + //-------------------------------------------------------------------- + // Init so2 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ so2; + // File name, name and sectors + so2.data_file = m_params.get("srf_emis_specifier_for_SO2"); + so2.species_name = "so2"; + so2.sectors = {"AGR", "RCO", "SHP", "SLV", "TRA", "WST"}; + srf_emiss_species_.push_back(so2); // add to the vector + //-------------------------------------------------------------------- + // Init bc_a4 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ bc_a4; + // File name, name and sectors + bc_a4.data_file = m_params.get("srf_emis_specifier_for_bc_a4"); + bc_a4.species_name = "bc_a4"; + bc_a4.sectors = {"AGR", "ENE", "IND", "RCO", "SHP", "SLV", "TRA", "WST"}; + srf_emiss_species_.push_back(bc_a4); // add to the vector + + //-------------------------------------------------------------------- + // Init num_a1 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ num_a1; + // File name, name and sectors + num_a1.data_file = m_params.get("srf_emis_specifier_for_num_a1"); + num_a1.species_name = "num_a1"; + num_a1.sectors = {"num_a1_SO4_AGR", "num_a1_SO4_SHP", "num_a1_SO4_SLV", + "num_a1_SO4_WST"}; + srf_emiss_species_.push_back(num_a1); // add to the vector + + //-------------------------------------------------------------------- + // Init num_a2 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ num_a2; + // File name, name and sectors + num_a2.data_file = m_params.get("srf_emis_specifier_for_num_a2"); + num_a2.species_name = "num_a2"; + num_a2.sectors = {"num_a2_SO4_RCO", "num_a2_SO4_TRA"}; + srf_emiss_species_.push_back(num_a2); // add to the vector + + //-------------------------------------------------------------------- + // Init num_a4 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ num_a4; + // File name, name and sectors + num_a4.data_file = m_params.get("srf_emis_specifier_for_num_a4"); + num_a4.species_name = "num_a4"; + num_a4.sectors = { + "num_a1_BC_AGR", "num_a1_BC_ENE", "num_a1_BC_IND", "num_a1_BC_RCO", + "num_a1_BC_SHP", "num_a1_BC_SLV", "num_a1_BC_TRA", "num_a1_BC_WST", + "num_a1_POM_AGR", "num_a1_POM_ENE", "num_a1_POM_IND", "num_a1_POM_RCO", + "num_a1_POM_SHP", "num_a1_POM_SLV", "num_a1_POM_TRA", "num_a1_POM_WST"}; + srf_emiss_species_.push_back(num_a4); // add to the vector + + //-------------------------------------------------------------------- + // Init pom_a4 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ pom_a4; + // File name, name and sectors + pom_a4.data_file = m_params.get("srf_emis_specifier_for_pom_a4"); + pom_a4.species_name = "pom_a4"; + pom_a4.sectors = {"AGR", "ENE", "IND", "RCO", "SHP", "SLV", "TRA", "WST"}; + srf_emiss_species_.push_back(pom_a4); // add to the vector + + //-------------------------------------------------------------------- + // Init so4_a1 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ so4_a1; + // File name, name and sectors + so4_a1.data_file = m_params.get("srf_emis_specifier_for_so4_a1"); + so4_a1.species_name = "so4_a1"; + so4_a1.sectors = {"AGR", "SHP", "SLV", "WST"}; + srf_emiss_species_.push_back(so4_a1); + + //-------------------------------------------------------------------- + // Init so4_a2 srf emiss data structures + //-------------------------------------------------------------------- + srf_emiss_ so4_a2; + // File name, name and sectors + so4_a2.data_file = m_params.get("srf_emis_specifier_for_so4_a2"); + so4_a2.species_name = "so4_a2"; + so4_a2.sectors = {"RCO", "TRA"}; + srf_emiss_species_.push_back(so4_a2); + + //-------------------------------------------------------------------- + // Init data structures to read and interpolate + //-------------------------------------------------------------------- + for(srf_emiss_ &ispec_srf : srf_emiss_species_) { + srfEmissFunc::init_srf_emiss_objects( + ncol_, grid_, ispec_srf.data_file, ispec_srf.sectors, srf_map_file, + // output + ispec_srf.horizInterp_, ispec_srf.data_start_, ispec_srf.data_end_, + ispec_srf.data_out_, ispec_srf.dataReader_); + } + +} // set_grid ends + +// ================================================================ +// REQUEST_BUFFER_SIZE_IN_BYTES +// ================================================================ +// ON HOST, returns the number of bytes of device memory needed by the above +// Buffer type given the number of columns and vertical levels +size_t MAMSrfOnlineEmiss::requested_buffer_size_in_bytes() const { + return mam_coupling::buffer_size(ncol_, nlev_); +} + +// ================================================================ +// INIT_BUFFERS +// ================================================================ +// ON HOST, initializes the Buffer type with sufficient memory to store +// intermediate (dry) quantities on the given number of columns with the given +// number of vertical levels. Returns the number of bytes allocated. +void MAMSrfOnlineEmiss::init_buffers(const ATMBufferManager &buffer_manager) { + EKAT_REQUIRE_MSG( + buffer_manager.allocated_bytes() >= requested_buffer_size_in_bytes(), + "Error! Insufficient buffer size.\n"); + + size_t used_mem = + mam_coupling::init_buffer(buffer_manager, ncol_, nlev_, buffer_); + EKAT_REQUIRE_MSG( + used_mem == requested_buffer_size_in_bytes(), + "Error! Used memory != requested memory for MAMSrfOnlineEmiss."); +} + +// ================================================================ +// INITIALIZE_IMPL +// ================================================================ +void MAMSrfOnlineEmiss::initialize_impl(const RunType run_type) { + // --------------------------------------------------------------- + // Output fields + // --------------------------------------------------------------- + // Constituent fluxes of species in [kg/m2/s] + constituent_fluxes_ = get_field_out("constituent_fluxes").get_view(); + + // --------------------------------------------------------------- + // Allocate memory for local and work arrays + // --------------------------------------------------------------- + + // Work array to store fluxes after unit conversions to kg/m2/s + fluxes_in_mks_units_ = view_1d("fluxes_in_mks_units", ncol_); + + // Current month ( 0-based) + const int curr_month = timestamp().get_month() - 1; + + // Load the first month into data_end. + + // Note: At the first time step, the data will be moved into data_beg, + // and data_end will be reloaded from file with the new month. + + //-------------------------------------------------------------------- + // Update surface emissions from file + //-------------------------------------------------------------------- + for(srf_emiss_ &ispec_srf : srf_emiss_species_) { + srfEmissFunc::update_srfEmiss_data_from_file( + ispec_srf.dataReader_, timestamp(), curr_month, *ispec_srf.horizInterp_, + ispec_srf.data_end_); // output + } + + //----------------------------------------------------------------- + // Setup preprocessing and post processing + //----------------------------------------------------------------- + preprocess_.initialize(constituent_fluxes_); + +} // end initialize_impl() + +// ================================================================ +// RUN_IMPL +// ================================================================ +void MAMSrfOnlineEmiss::run_impl(const double dt) { + // Zero-out output + Kokkos::deep_copy(preprocess_.constituent_fluxes_pre_, 0); + + // Gather time and state information for interpolation + auto ts = timestamp() + dt; + + //-------------------------------------------------------------------- + // Interpolate srf emiss data + //-------------------------------------------------------------------- + + for(srf_emiss_ &ispec_srf : srf_emiss_species_) { + // Update TimeState, note the addition of dt + ispec_srf.timeState_.t_now = ts.frac_of_year_in_days(); + + // Update time state and if the month has changed, update the data. + srfEmissFunc::update_srfEmiss_timestate( + ispec_srf.dataReader_, ts, *ispec_srf.horizInterp_, + // output + ispec_srf.timeState_, ispec_srf.data_start_, ispec_srf.data_end_); + + // Call the main srfEmiss routine to get interpolated aerosol forcings. + srfEmissFunc::srfEmiss_main(ispec_srf.timeState_, ispec_srf.data_start_, + ispec_srf.data_end_, ispec_srf.data_out_); + + //-------------------------------------------------------------------- + // Modify units to MKS units (from molecules/cm2/s to kg/m2/s) + //-------------------------------------------------------------------- + // Get species index in array with pcnst dimension (e.g., state_q or + // constituent_fluxes_) + const int species_index = spcIndex_in_pcnst_.at(ispec_srf.species_name); + + // modify units from molecules/cm2/s to kg/m2/s + auto fluxes_in_mks_units = this->fluxes_in_mks_units_; + auto constituent_fluxes = this->constituent_fluxes_; + const Real mfactor = + amufac * mam4::gas_chemistry::adv_mass[species_index - offset_]; + // Parallel loop over all the columns to update units + Kokkos::parallel_for( + "fluxes", ncol_, KOKKOS_LAMBDA(int icol) { + fluxes_in_mks_units(icol) = + ispec_srf.data_out_.emiss_sectors(0, icol) * mfactor; + constituent_fluxes(icol, species_index) = fluxes_in_mks_units(icol); + }); + + } // for loop for species + Kokkos::fence(); +} // run_imple ends + +// ============================================================================= +} // namespace scream diff --git a/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.hpp b/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.hpp new file mode 100644 index 00000000000..031fb62d8b7 --- /dev/null +++ b/components/eamxx/src/physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.hpp @@ -0,0 +1,133 @@ +#ifndef EAMXX_MAM_SRF_ONLINE_EMISS_HPP +#define EAMXX_MAM_SRF_ONLINE_EMISS_HPP + +#include "share/grid/remap/abstract_remapper.hpp" +#include "share/io/scorpio_input.hpp" + +// For MAM4 aerosol configuration +#include +#include + +// For declaring surface and online emission class derived from atm process +// class +#include + +// #include +#include + +namespace scream { + +// The process responsible for handling MAM4 surface and online emissions. The +// AD stores exactly ONE instance of this class in its list of subcomponents. +class MAMSrfOnlineEmiss final : public scream::AtmosphereProcess { + using KT = ekat::KokkosTypes; + using view_1d = typename KT::template view_1d; + using view_2d = typename KT::template view_2d; + + // number of horizontal columns and vertical levels + int ncol_, nlev_; + + // buffer for sotring temporary variables + mam_coupling::Buffer buffer_; + + // physics grid for column information + std::shared_ptr grid_; + + // Constituent fluxes of species in [kg/m2/s] + view_2d constituent_fluxes_; + + // Work array to store fluxes after unit conversions to kg/m2/s + view_1d fluxes_in_mks_units_; + + // Unified atomic mass unit used for unit conversion (BAD constant) + static constexpr Real amufac = 1.65979e-23; // 1.e4* kg / amu + + public: + using srfEmissFunc = mam_coupling::srfEmissFunctions; + + // Constructor + MAMSrfOnlineEmiss(const ekat::Comm &comm, const ekat::ParameterList ¶ms); + + // -------------------------------------------------------------------------- + // AtmosphereProcess overrides (see share/atm_process/atmosphere_process.hpp) + // -------------------------------------------------------------------------- + + // The type of subcomponent + AtmosphereProcessType type() const { return AtmosphereProcessType::Physics; } + + // The name of the subcomponent + std::string name() const { return "mam_srf_online_emissions"; } + + // grid + void set_grids( + const std::shared_ptr grids_manager) override; + + // management of common atm process memory + size_t requested_buffer_size_in_bytes() const override; + void init_buffers(const ATMBufferManager &buffer_manager) override; + + // Initialize variables + void initialize_impl(const RunType run_type) override; + + // Run the process by one time step + void run_impl(const double dt) override; + + // Finalize + void finalize_impl(){/*Do nothing*/}; + // Atmosphere processes often have a pre-processing step that constructs + // required variables from the set of fields stored in the field manager. + // This functor implements this step, which is called during run_impl. + struct Preprocess { + Preprocess() = default; + // on host: initializes preprocess functor with necessary state data + void initialize(const view_2d &constituent_fluxes) { + constituent_fluxes_pre_ = constituent_fluxes; + } + // local variables for preprocess struct + view_2d constituent_fluxes_pre_; + }; // MAMSrfOnlineEmiss::Preprocess + + private: + // preprocessing scratch pad + Preprocess preprocess_; + + // Species index (zero-based) in tracer array with "pcnst" dimension + // FIXME: Remove the hardwired indices and use a function + // to find them from an array. + const std::map spcIndex_in_pcnst_ = { + {"so2", 12}, {"dms", 13}, {"so4_a1", 15}, + {"num_a1", 22}, {"so4_a2", 23}, {"num_a2", 27}, + {"pom_a4", 36}, {"bc_a4", 37}, {"num_a4", 39}}; + + // A struct carrying all the fields needed to read + // surface emissions of a species + struct srf_emiss_ { + // species name + std::string species_name; + + // Data file name + std::string data_file; + + // Sector names in file + std::vector sectors; + + // Data structure for reading interpolation + std::shared_ptr horizInterp_; + std::shared_ptr dataReader_; + srfEmissFunc::srfEmissTimeState timeState_; + srfEmissFunc::srfEmissInput data_start_, data_end_; + srfEmissFunc::srfEmissOutput data_out_; + }; + + // A vector for carrying emissions for all the species + std::vector srf_emiss_species_; + + // offset for converting pcnst index to gas_pcnst index + static constexpr int offset_ = + mam4::aero_model::pcnst - mam4::gas_chemistry::gas_pcnst; + +}; // MAMSrfOnlineEmiss + +} // namespace scream + +#endif // EAMXX_MAM_SRF_ONLINE_EMISS_HPP diff --git a/components/eamxx/src/physics/mam/mam_coupling.hpp b/components/eamxx/src/physics/mam/mam_coupling.hpp index d490ab76155..a3a5f746aa3 100644 --- a/components/eamxx/src/physics/mam/mam_coupling.hpp +++ b/components/eamxx/src/physics/mam/mam_coupling.hpp @@ -776,6 +776,24 @@ void compute_wet_mixing_ratios(const Team& team, }); } +// Computes the reciprocal of pseudo density for a column +inline +void compute_recipical_pseudo_density(haero::ThreadTeamPolicy team_policy, + const_view_2d pdel, + const int nlev, + // output + view_2d rpdel) { + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const haero::ThreadTeam &team) { + const int icol = team.league_rank(); + Kokkos::parallel_for( + Kokkos::TeamVectorRange(team, 0, nlev), [&](int kk) { + EKAT_KERNEL_ASSERT_MSG(0 < pdel(icol, kk), + "Error: pdel should be > 0.\n"); + rpdel(icol, kk) = 1 / pdel(icol, kk); + }); + }); +} // Scream (or EAMxx) can sometimes extend views beyond model levels (nlev) as it uses // "packs". Following function copies a 2d view till model levels inline diff --git a/components/eamxx/src/physics/mam/srf_emission.hpp b/components/eamxx/src/physics/mam/srf_emission.hpp new file mode 100644 index 00000000000..29aaca421ea --- /dev/null +++ b/components/eamxx/src/physics/mam/srf_emission.hpp @@ -0,0 +1,138 @@ +#ifndef SRF_EMISSION_HPP +#define SRF_EMISSION_HPP + +#include "share/util/scream_timing.hpp" + +namespace scream::mam_coupling { +namespace { + +template +struct srfEmissFunctions { + using Device = DeviceType; + + using KT = KokkosTypes; + using MemberType = typename KT::MemberType; + + struct srfEmissTimeState { + srfEmissTimeState() = default; + // Whether the timestate has been initialized. + // The current month + int current_month = -1; + // Julian Date for the beginning of the month, as defined in + // /src/share/util/scream_time_stamp.hpp + // See this file for definition of Julian Date. + Real t_beg_month; + // Current simulation Julian Date + Real t_now; + // Number of days in the current month, cast as a Real + Real days_this_month; + }; // srfEmissTimeState + + struct srfEmissData { + srfEmissData() = default; + srfEmissData(const int ncol_, const int nsectors_) + : ncols(ncol_), nsectors(nsectors_) { + init(ncols, nsectors, true); + } + + void init(const int ncol, const int nsector, const bool allocate) { + ncols = ncol; + nsectors = nsector; + if(allocate) emiss_sectors = view_2d("AllSectors", nsectors, ncols); + } // srfEmissData init + + // Basic spatial dimensions of the data + int ncols, nsectors; + view_2d emiss_sectors; + }; // srfEmissData + + struct srfEmissInput { + srfEmissInput() = default; + srfEmissInput(const int ncols_, const int nsectors_) { + init(ncols_, nsectors_); + } + + void init(const int ncols_, const int nsectors_) { + data.init(ncols_, nsectors_, true); + } + srfEmissData data; // All srfEmiss fields + }; // srfEmissInput + + // The output is really just srfEmissData, but for clarity it might + // help to see a srfEmissOutput along a srfEmissInput in functions signatures + using srfEmissOutput = srfEmissData; + + /* ------------------------------------------------------------------------------------------- + */ + // Surface emissions routines + template + static std::shared_ptr create_horiz_remapper( + const std::shared_ptr &model_grid, + const std::string &srfEmiss_data_file, + const std::array &field_names, + const std::string &map_file); + + static std::shared_ptr create_horiz_remapper( + const std::shared_ptr &model_grid, + const std::string &srfEmiss_data_file, + const std::vector &field_names, const std::string &map_file); + + static std::shared_ptr create_srfEmiss_data_reader( + const std::shared_ptr &horiz_remapper, + const std::string &srfEmiss_data_file); + + static void srfEmiss_main(const srfEmissTimeState &time_state, + const srfEmissInput &data_beg, + const srfEmissInput &data_end, + const srfEmissOutput &data_out); + + static void update_srfEmiss_data_from_file( + std::shared_ptr &scorpio_reader, + const util::TimeStamp &ts, + const int time_index, // zero-based + AbstractRemapper &srfEmiss_horiz_interp, srfEmissInput &srfEmiss_input); + static void update_srfEmiss_timestate( + std::shared_ptr &scorpio_reader, + const util::TimeStamp &ts, AbstractRemapper &srfEmiss_horiz_interp, + srfEmissTimeState &time_state, srfEmissInput &srfEmiss_beg, + srfEmissInput &srfEmiss_end); + + // The following three are called during srfEmiss_main + static void perform_time_interpolation(const srfEmissTimeState &time_state, + const srfEmissInput &data_beg, + const srfEmissInput &data_end, + const srfEmissOutput &data_out); + + // Performs convex interpolation of x0 and x1 at point t + template + KOKKOS_INLINE_FUNCTION static ScalarX linear_interp(const ScalarX &x0, + const ScalarX &x1, + const ScalarT &t); + template + static void init_srf_emiss_objects( + const int ncol, const std::shared_ptr &grid, + const std::string &data_file, + const std::array §ors, + const std::string &srf_map_file, + // output + std::shared_ptr &SrfEmissHorizInterp, + srfEmissInput &SrfEmissData_start, srfEmissInput &SrfEmissData_end, + srfEmissOutput &SrfEmissData_out, + std::shared_ptr &SrfEmissDataReader); + + static void init_srf_emiss_objects( + const int ncol, const std::shared_ptr &grid, + const std::string &data_file, const std::vector §ors, + const std::string &srf_map_file, + // output + std::shared_ptr &SrfEmissHorizInterp, + srfEmissInput &SrfEmissData_start, srfEmissInput &SrfEmissData_end, + srfEmissOutput &SrfEmissData_out, + std::shared_ptr &SrfEmissDataReader); + +}; // struct srfEmissFunctions +} // namespace +} // namespace scream::mam_coupling +#endif // SRF_EMISSION_HPP + +#include "srf_emission_impl.hpp" diff --git a/components/eamxx/src/physics/mam/srf_emission_impl.hpp b/components/eamxx/src/physics/mam/srf_emission_impl.hpp new file mode 100644 index 00000000000..b8ac7de768d --- /dev/null +++ b/components/eamxx/src/physics/mam/srf_emission_impl.hpp @@ -0,0 +1,292 @@ +#ifndef SRF_EMISSION_IMPL_HPP +#define SRF_EMISSION_IMPL_HPP + +#include "share/grid/remap/identity_remapper.hpp" +#include "share/grid/remap/refining_remapper_p2p.hpp" +#include "share/io/scream_scorpio_interface.hpp" + +namespace scream::mam_coupling { +namespace { + +template +std::shared_ptr +srfEmissFunctions::create_horiz_remapper( + const std::shared_ptr &model_grid, + const std::string &data_file, const std::vector §or_names, + const std::string &map_file) { + using namespace ShortFieldTagsNames; + + scorpio::register_file(data_file, scorpio::Read); + const int ncols_data = scorpio::get_dimlen(data_file, "ncol"); + scorpio::release_file(data_file); + + // We could use model_grid directly if using same num levels, + // but since shallow clones are cheap, we may as well do it (less lines of + // code) + auto horiz_interp_tgt_grid = + model_grid->clone("srf_emiss_horiz_interp_tgt_grid", true); + + const int ncols_model = model_grid->get_num_global_dofs(); + std::shared_ptr remapper; + // if the file's grid is same as model's native grid, we identity remapper + // (i.e., no interpolation) + if(ncols_data == ncols_model) { + remapper = std::make_shared( + horiz_interp_tgt_grid, IdentityRemapper::SrcAliasTgt); + } else { + EKAT_REQUIRE_MSG(ncols_data <= ncols_model, + "Error! We do not allow to coarsen srfEmiss data to fit " + "the model. We only allow\n" + "srfEmiss data to be at the same or coarser resolution as " + "the model.\n"); + // We must have a valid map file + EKAT_REQUIRE_MSG( + map_file != "", + "ERROR: srfEmiss data is on a different grid than the model one,\n" + "but srfEmiss_remap_file is missing from srfEmiss parameter " + "list."); + + remapper = + std::make_shared(horiz_interp_tgt_grid, map_file); + } + + remapper->registration_begins(); + + const auto tgt_grid = remapper->get_tgt_grid(); + + const auto layout_2d = tgt_grid->get_2d_scalar_layout(); + const auto nondim = ekat::units::Units::nondimensional(); + + std::vector field_emiss_sectors; + + const int sector_size = sector_names.size(); + for(int icomp = 0; icomp < sector_size; ++icomp) { + auto comp_name = sector_names[icomp]; + // set and allocate fields + Field f(FieldIdentifier(comp_name, layout_2d, nondim, tgt_grid->name())); + f.allocate_view(); + field_emiss_sectors.push_back(f); + remapper->register_field_from_tgt(f); + } + + remapper->registration_ends(); + + return remapper; +} // create_horiz_remapper + +template +std::shared_ptr +srfEmissFunctions::create_srfEmiss_data_reader( + const std::shared_ptr &horiz_remapper, + const std::string &srfEmiss_data_file) { + std::vector field_emiss_sectors; + for(int i = 0; i < horiz_remapper->get_num_fields(); ++i) { + field_emiss_sectors.push_back(horiz_remapper->get_src_field(i)); + } + const auto io_grid = horiz_remapper->get_src_grid(); + return std::make_shared(srfEmiss_data_file, io_grid, + field_emiss_sectors, true); +} // create_srfEmiss_data_reader + +template +template +ScalarX srfEmissFunctions::linear_interp(const ScalarX &x0, + const ScalarX &x1, + const ScalarT &t) { + return (1 - t) * x0 + t * x1; +} // linear_interp + +template +void srfEmissFunctions::perform_time_interpolation( + const srfEmissTimeState &time_state, const srfEmissInput &data_beg, + const srfEmissInput &data_end, const srfEmissOutput &data_out) { + // NOTE: we *assume* data_beg and data_end have the *same* hybrid v coords. + // IF this ever ceases to be the case, you can interp those too. + + using ExeSpace = typename KT::ExeSpace; + using ESU = ekat::ExeSpaceUtils; + + // Gather time stamp info + auto &t_now = time_state.t_now; + auto &t_beg = time_state.t_beg_month; + auto &delta_t = time_state.days_this_month; + + // At this stage, begin/end must have the same dimensions + EKAT_REQUIRE(data_end.data.ncols == data_beg.data.ncols); + + auto delta_t_fraction = (t_now - t_beg) / delta_t; + + EKAT_REQUIRE_MSG(delta_t_fraction >= 0 && delta_t_fraction <= 1, + "Error! Convex interpolation with coefficient out of " + "[0,1].\n t_now : " + + std::to_string(t_now) + + "\n" + " t_beg : " + + std::to_string(t_beg) + + "\n delta_t: " + std::to_string(delta_t) + "\n"); + + const int nsectors = data_beg.data.nsectors; + const int ncols = data_beg.data.ncols; + using ExeSpace = typename KT::ExeSpace; + using ESU = ekat::ExeSpaceUtils; + const auto policy = ESU::get_default_team_policy(ncols, nsectors); + + Kokkos::parallel_for( + policy, KOKKOS_LAMBDA(const MemberType &team) { + const int icol = team.league_rank(); // column index + Real accum = 0; + // Parallel reduction over sectors + // FIXME: Do we need to use Kokkos::Single for each team here??? + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange(team, nsectors), + [&](const int i, Real &update) { + const auto beg = data_beg.data.emiss_sectors(i, icol); + const auto end = data_end.data.emiss_sectors(i, icol); + update += linear_interp(beg, end, delta_t_fraction); + }, + accum); + // Assign the accumulated value to the output + data_out.emiss_sectors(0, icol) = accum; + }); + Kokkos::fence(); +} // perform_time_interpolation + +template +void srfEmissFunctions::srfEmiss_main(const srfEmissTimeState &time_state, + const srfEmissInput &data_beg, + const srfEmissInput &data_end, + const srfEmissOutput &data_out) { + // Beg/End/Tmp month must have all sizes matching + + EKAT_REQUIRE_MSG( + data_end.data.ncols == data_beg.data.ncols, + "Error! srfEmissInput data structs must have the same number of " + "columns/levels.\n"); + + // Horiz interpolation can be expensive, and does not depend on the particular + // time of the month, so it can be done ONCE per month, *outside* + // srfEmiss_main (when updating the beg/end states, reading them from file). + EKAT_REQUIRE_MSG( + data_end.data.ncols == data_out.ncols, + "Error! Horizontal interpolation is performed *before* " + "calling srfEmiss_main,\n" + " srfEmissInput and srfEmissOutput data structs must have the " + "same number columns.\n"); + + // Step 1. Perform time interpolation + perform_time_interpolation(time_state, data_beg, data_end, data_out); +} // srfEmiss_main + +template +void srfEmissFunctions::update_srfEmiss_data_from_file( + std::shared_ptr &scorpio_reader, const util::TimeStamp &ts, + const int time_index, // zero-based + AbstractRemapper &srfEmiss_horiz_interp, srfEmissInput &srfEmiss_input) { + using namespace ShortFieldTagsNames; + using ESU = ekat::ExeSpaceUtils; + using Member = typename KokkosTypes::MemberType; + + start_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file"); + + // 1. Read from file + start_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::read_data"); + scorpio_reader->read_variables(time_index); + stop_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::read_data"); + + // 2. Run the horiz remapper (it is a do-nothing op if srfEmiss data is on + // same grid as model) + start_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::horiz_remap"); + srfEmiss_horiz_interp.remap(/*forward = */ true); + stop_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::horiz_remap"); + + // 3. Copy from the tgt field of the remapper into the srfEmiss_data, padding + // data if necessary + start_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::copy_and_pad"); + // Recall, the fields are registered in the order: ps, ccn3, g_sw, ssa_sw, + // tau_sw, tau_lw + + const auto &layout = srfEmiss_horiz_interp.get_tgt_field(0) + .get_header() + .get_identifier() + .get_layout(); + + const int ncols = layout.dim(COL); + + // Read fields from the file + for(int i = 0; i < srfEmiss_horiz_interp.get_num_fields(); ++i) { + auto sector = + srfEmiss_horiz_interp.get_tgt_field(i).get_view(); + const auto emiss = + Kokkos::subview(srfEmiss_input.data.emiss_sectors, i, Kokkos::ALL()); + Kokkos::deep_copy(emiss, sector); + } + + Kokkos::fence(); + stop_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file::copy_and_pad"); + + stop_timer("EAMxx::srfEmiss::update_srfEmiss_data_from_file"); + +} // END update_srfEmiss_data_from_file + +template +void srfEmissFunctions::update_srfEmiss_timestate( + std::shared_ptr &scorpio_reader, const util::TimeStamp &ts, + AbstractRemapper &srfEmiss_horiz_interp, srfEmissTimeState &time_state, + srfEmissInput &srfEmiss_beg, srfEmissInput &srfEmiss_end) { + // Now we check if we have to update the data that changes monthly + // NOTE: This means that srfEmiss assumes monthly data to update. Not + // any other frequency. + const auto month = ts.get_month() - 1; // Make it 0-based + if(month != time_state.current_month) { + // Update the srfEmiss time state information + time_state.current_month = month; + time_state.t_beg_month = + util::TimeStamp({ts.get_year(), month + 1, 1}, {0, 0, 0}) + .frac_of_year_in_days(); + time_state.days_this_month = util::days_in_month(ts.get_year(), month + 1); + + // Copy srfEmiss_end'data into srfEmiss_beg'data, and read in the new + // srfEmiss_end + std::swap(srfEmiss_beg, srfEmiss_end); + + // Update the srfEmiss forcing data for this month and next month + // Start by copying next months data to this months data structure. + // NOTE: If the timestep is bigger than monthly this could cause the wrong + // values + // to be assigned. A timestep greater than a month is very unlikely + // so we will proceed. + int next_month = (time_state.current_month + 1) % 12; + update_srfEmiss_data_from_file(scorpio_reader, ts, next_month, + srfEmiss_horiz_interp, srfEmiss_end); + } + +} // END updata_srfEmiss_timestate + +template +void srfEmissFunctions::init_srf_emiss_objects( + const int ncol, const std::shared_ptr &grid, + const std::string &data_file, const std::vector §ors, + const std::string &srf_map_file, + // output + std::shared_ptr &SrfEmissHorizInterp, + srfEmissInput &SrfEmissData_start, srfEmissInput &SrfEmissData_end, + srfEmissOutput &SrfEmissData_out, + std::shared_ptr &SrfEmissDataReader) { + // Init horizontal remap + SrfEmissHorizInterp = + create_horiz_remapper(grid, data_file, sectors, srf_map_file); + + // Initialize the size of start/end/out data structures + SrfEmissData_start = srfEmissInput(ncol, sectors.size()); + SrfEmissData_end = srfEmissInput(ncol, sectors.size()); + SrfEmissData_out.init(ncol, 1, true); + + // Create reader (an AtmosphereInput object) + SrfEmissDataReader = + create_srfEmiss_data_reader(SrfEmissHorizInterp, data_file); +} // init_srf_emiss_objects + +} // namespace +} // namespace scream::mam_coupling + +#endif // SRF_EMISSION_IMPL_HPP diff --git a/components/eamxx/src/physics/register_physics.hpp b/components/eamxx/src/physics/register_physics.hpp index aea9f2f543a..d837e96311f 100644 --- a/components/eamxx/src/physics/register_physics.hpp +++ b/components/eamxx/src/physics/register_physics.hpp @@ -29,6 +29,8 @@ #include "physics/mam/eamxx_mam_dry_deposition_process_interface.hpp" #include "physics/mam/eamxx_mam_aci_process_interface.hpp" #include "physics/mam/eamxx_mam_wetscav_process_interface.hpp" +#include "physics/mam/eamxx_mam_srf_and_online_emissions_process_interface.hpp" +#include "physics/mam/eamxx_mam_constituent_fluxes_interface.hpp" #endif #ifdef EAMXX_HAS_COSP #include "physics/cosp/eamxx_cosp.hpp" @@ -68,6 +70,8 @@ inline void register_physics () { proc_factory.register_product("mam4_drydep",&create_atmosphere_process); proc_factory.register_product("mam4_aci",&create_atmosphere_process); proc_factory.register_product("mam4_wetscav",&create_atmosphere_process); + proc_factory.register_product("mam4_srf_online_emiss",&create_atmosphere_process); + proc_factory.register_product("mam4_constituent_fluxes",&create_atmosphere_process); #endif #ifdef EAMXX_HAS_COSP proc_factory.register_product("Cosp",&create_atmosphere_process); diff --git a/components/eamxx/tests/multi-process/physics_only/CMakeLists.txt b/components/eamxx/tests/multi-process/physics_only/CMakeLists.txt index d6ed13dfb0a..bb23911ff1c 100644 --- a/components/eamxx/tests/multi-process/physics_only/CMakeLists.txt +++ b/components/eamxx/tests/multi-process/physics_only/CMakeLists.txt @@ -11,6 +11,7 @@ if (SCREAM_DOUBLE_PRECISION) add_subdirectory(mam/shoc_cldfrac_mam4_aci_p3_mam4_optics_rrtmgp) add_subdirectory(mam/p3_mam4_wetscav) add_subdirectory(mam/shoc_cldfrac_p3_wetscav) + add_subdirectory(mam/mam4_srf_online_emiss_mam4_constituent_fluxes) endif() endif() diff --git a/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/CMakeLists.txt b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/CMakeLists.txt new file mode 100644 index 00000000000..886936ca084 --- /dev/null +++ b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/CMakeLists.txt @@ -0,0 +1,53 @@ +INCLUDE (ScreamUtils) + +set (TEST_BASE_NAME mam4_srf_online_emiss_mam4_constituent_fluxes) +set (FIXTURES_BASE_NAME ${TEST_BASE_NAME}_generate_output_nc_files) + +# Create the test +CreateADUnitTest(${TEST_BASE_NAME} + LIBS mam + LABELS physics mam4_srf_online_emiss mam4_constituent_fluxes + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + FIXTURES_SETUP_INDIVIDUAL ${FIXTURES_BASE_NAME} +) + +# Set AD configurable options +set (ATM_TIME_STEP 1800) +SetVarDependingOnTestSize(NUM_STEPS 2 5 48) # 1h 4h 24h +set (RUN_T0 2021-10-12-45000) + +# Ensure test input files are present in the data dir +GetInputFile(scream/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev}) +GetInputFile(cam/topo/${EAMxx_tests_TOPO_FILE}) + +# Ensure test input files are present in the data dir +set (TEST_INPUT_FILES + scream/mam4xx/emissions/ne2np4/surface/DMSflux.2010.ne2np4_conserv.POPmonthlyClimFromACES4BGC_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so2_surf_ne2np4_2010_clim_c20240723.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_bc_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a1_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a2_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_pom_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a1_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a2_surf_ne2np4_2010_clim_c20240726.nc +) +foreach (file IN ITEMS ${TEST_INPUT_FILES}) + GetInputFile(${file}) +endforeach() + +## Copy (and configure) yaml files needed by tests +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input.yaml) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/output.yaml + ${CMAKE_CURRENT_BINARY_DIR}/output.yaml) + +# Compare output files produced by npX tests, to ensure they are bfb +include (CompareNCFiles) +CompareNCFilesFamilyMpi ( + TEST_BASE_NAME ${TEST_BASE_NAME} + FILE_META_NAME ${TEST_BASE_NAME}_output.INSTANT.nsteps_x${NUM_STEPS}.npMPIRANKS.${RUN_T0}.nc + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + LABELS physics mam4_srf_online_emiss mam4_constituent_fluxes + META_FIXTURES_REQUIRED ${FIXTURES_BASE_NAME}_npMPIRANKS_omp1 +) diff --git a/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/input.yaml b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/input.yaml new file mode 100644 index 00000000000..55c62d69aa2 --- /dev/null +++ b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/input.yaml @@ -0,0 +1,47 @@ +%YAML 1.1 +--- +driver_options: + atmosphere_dag_verbosity_level: 5 + +time_stepping: + time_step: ${ATM_TIME_STEP} + run_t0: ${RUN_T0} # YYYY-MM-DD-XXXXX + number_of_steps: ${NUM_STEPS} + +atmosphere_processes: + atm_procs_list: [mam4_srf_online_emiss, mam4_constituent_fluxes] + schedule_type: Sequential + mam4_srf_online_emiss: + # MAM4xx-Surface-Emissions + srf_remap_file: "" + srf_emis_specifier_for_DMS: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/DMSflux.2010.ne2np4_conserv.POPmonthlyClimFromACES4BGC_c20240726.nc + srf_emis_specifier_for_SO2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so2_surf_ne2np4_2010_clim_c20240723.nc + srf_emis_specifier_for_bc_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_bc_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a1: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a1_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a2_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_pom_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_pom_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_so4_a1: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a1_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_so4_a2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a2_surf_ne2np4_2010_clim_c20240726.nc + +grids_manager: + Type: Mesh Free + geo_data_source: IC_FILE + grids_names: [Physics GLL] + Physics GLL: + aliases: [Physics] + type: point_grid + number_of_global_columns: 218 + number_of_vertical_levels: 72 + +initial_conditions: + # The name of the file containing the initial conditions for this test. + Filename: ${SCREAM_DATA_DIR}/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev} + topography_filename: ${TOPO_DATA_DIR}/${EAMxx_tests_TOPO_FILE} + + pbl_height: 0.0 + +# The parameters for I/O control +Scorpio: + output_yaml_files: ["output.yaml"] +... diff --git a/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/output.yaml b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/output.yaml new file mode 100644 index 00000000000..c9cad5e2761 --- /dev/null +++ b/components/eamxx/tests/multi-process/physics_only/mam/mam4_srf_online_emiss_mam4_constituent_fluxes/output.yaml @@ -0,0 +1,41 @@ +%YAML 1.1 +--- +filename_prefix: mam4_srf_online_emiss_mam4_constituent_fluxes_output +Averaging Type: Instant +Field Names: + - bc_a1 + - bc_a3 + - bc_a4 + - dst_a1 + - dst_a3 + - so4_a1 + - so4_a2 + - so4_a3 + - pom_a1 + - pom_a3 + - pom_a4 + - soa_a1 + - soa_a2 + - soa_a3 + - nacl_a1 + - nacl_a2 + - nacl_a3 + - mom_a1 + - mom_a2 + - mom_a3 + - mom_a4 + - num_a1 + - num_a2 + - num_a3 + - num_a4 + - constituent_fluxes + - O3 + - H2O2 + - H2SO4 + - SO2 + - DMS + - SOAG +output_control: + Frequency: ${NUM_STEPS} + frequency_units: nsteps +... diff --git a/components/eamxx/tests/single-process/CMakeLists.txt b/components/eamxx/tests/single-process/CMakeLists.txt index 95e4d81ccbb..40565d94bf8 100644 --- a/components/eamxx/tests/single-process/CMakeLists.txt +++ b/components/eamxx/tests/single-process/CMakeLists.txt @@ -23,6 +23,8 @@ if (SCREAM_ENABLE_MAM) add_subdirectory(mam/aci) add_subdirectory(mam/drydep) add_subdirectory(mam/wet_scav) + add_subdirectory(mam/emissions) + add_subdirectory(mam/constituent_fluxes) endif() if (SCREAM_TEST_LEVEL GREATER_EQUAL SCREAM_TEST_LEVEL_EXPERIMENTAL) add_subdirectory(zm) diff --git a/components/eamxx/tests/single-process/mam/constituent_fluxes/CMakeLists.txt b/components/eamxx/tests/single-process/mam/constituent_fluxes/CMakeLists.txt new file mode 100644 index 00000000000..2615f152b26 --- /dev/null +++ b/components/eamxx/tests/single-process/mam/constituent_fluxes/CMakeLists.txt @@ -0,0 +1,44 @@ +include (ScreamUtils) + +set (TEST_BASE_NAME mam4_constituent_fluxes_standalone) +set (FIXTURES_BASE_NAME ${TEST_BASE_NAME}_generate_output_nc_files) + +# Create the test +CreateADUnitTest(${TEST_BASE_NAME} + LABELS mam4_constituent_fluxes physics + LIBS mam + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + FIXTURES_SETUP_INDIVIDUAL ${FIXTURES_BASE_NAME} +) + +# Set AD configurable options +set (ATM_TIME_STEP 1800) +SetVarDependingOnTestSize(NUM_STEPS 2 5 48) # 1h 2.5h 24h +set (RUN_T0 2021-10-12-45000) + +## Copy (and configure) yaml files needed by tests +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input.yaml) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/output.yaml + ${CMAKE_CURRENT_BINARY_DIR}/output.yaml) + +# Ensure test input files are present in the data dir +GetInputFile(scream/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev}) + +# Compare output files produced by npX tests, to ensure they are bfb +include (CompareNCFiles) + +CompareNCFilesFamilyMpi ( + TEST_BASE_NAME ${TEST_BASE_NAME} + FILE_META_NAME ${TEST_BASE_NAME}_output.INSTANT.nsteps_x1.npMPIRANKS.${RUN_T0}.nc + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + LABELS mam4_constituent_fluxes physics + META_FIXTURES_REQUIRED ${FIXTURES_BASE_NAME}_npMPIRANKS_omp1 +) + +if (SCREAM_ENABLE_BASELINE_TESTS) + # Compare one of the output files with the baselines. + # Note: one is enough, since we already check that np1 is BFB with npX + set (OUT_FILE ${TEST_BASE_NAME}_output.INSTANT.nsteps_x1.np${TEST_RANK_END}.${RUN_T0}.nc) + CreateBaselineTest(${TEST_BASE_NAME} ${TEST_RANK_END} ${OUT_FILE} ${FIXTURES_BASE_NAME}) +endif() diff --git a/components/eamxx/tests/single-process/mam/constituent_fluxes/input.yaml b/components/eamxx/tests/single-process/mam/constituent_fluxes/input.yaml new file mode 100644 index 00000000000..d012274d146 --- /dev/null +++ b/components/eamxx/tests/single-process/mam/constituent_fluxes/input.yaml @@ -0,0 +1,37 @@ +%YAML 1.1 +--- +driver_options: + atmosphere_dag_verbosity_level: 5 + +time_stepping: + time_step: ${ATM_TIME_STEP} + run_t0: ${RUN_T0} # YYYY-MM-DD-XXXXX + number_of_steps: ${NUM_STEPS} + +atmosphere_processes: + atm_procs_list: [mam4_constituent_fluxes] +grids_manager: + Type: Mesh Free + geo_data_source: IC_FILE + grids_names: [Physics GLL] + Physics GLL: + type: point_grid + aliases: [Physics] + number_of_global_columns: 218 + number_of_vertical_levels: 72 + +initial_conditions: + # The name of the file containing the initial conditions for this test. + Filename: ${SCREAM_DATA_DIR}/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev} + topography_filename: ${TOPO_DATA_DIR}/${EAMxx_tests_TOPO_FILE} + phis : 1.0 + #These should come from the input file + + #we should get the following variables from other processes + pbl_height : 1.0 + constituent_fluxes: 1.e-5 + +# The parameters for I/O control +Scorpio: + output_yaml_files: ["output.yaml"] +... diff --git a/components/eamxx/tests/single-process/mam/constituent_fluxes/output.yaml b/components/eamxx/tests/single-process/mam/constituent_fluxes/output.yaml new file mode 100644 index 00000000000..b237056b5cc --- /dev/null +++ b/components/eamxx/tests/single-process/mam/constituent_fluxes/output.yaml @@ -0,0 +1,43 @@ +%YAML 1.1 +--- +filename_prefix: mam4_constituent_fluxes_standalone_output +Averaging Type: Instant +Fields: + Physics: + Field Names: + - O3 + - H2O2 + - H2SO4 + - SO2 + - DMS + - SOAG + - bc_a1 + - bc_a3 + - bc_a4 + - dst_a1 + - dst_a3 + - so4_a1 + - so4_a2 + - so4_a3 + - pom_a1 + - pom_a3 + - pom_a4 + - soa_a1 + - soa_a2 + - soa_a3 + - nacl_a1 + - nacl_a2 + - nacl_a3 + - mom_a1 + - mom_a2 + - mom_a3 + - mom_a4 + - num_a1 + - num_a2 + - num_a3 + - num_a4 + +output_control: + Frequency: 1 + frequency_units: nsteps +... diff --git a/components/eamxx/tests/single-process/mam/emissions/CMakeLists.txt b/components/eamxx/tests/single-process/mam/emissions/CMakeLists.txt new file mode 100644 index 00000000000..c8e690b929f --- /dev/null +++ b/components/eamxx/tests/single-process/mam/emissions/CMakeLists.txt @@ -0,0 +1,60 @@ +include (ScreamUtils) + +set (TEST_BASE_NAME mam4_srf_online_emiss_standalone) +set (FIXTURES_BASE_NAME ${TEST_BASE_NAME}_generate_output_nc_files) + +# Create the test +CreateADUnitTest(${TEST_BASE_NAME} + LABELS mam4_srf_online_emiss physics + LIBS mam + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + FIXTURES_SETUP_INDIVIDUAL ${FIXTURES_BASE_NAME} +) + +# Set AD configurable options +set (ATM_TIME_STEP 1800) +SetVarDependingOnTestSize(NUM_STEPS 2 5 48) # 1h 2.5h 24h +set (RUN_T0 2021-10-12-45000) + +## Copy (and configure) yaml files needed by tests +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input.yaml) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/output.yaml + ${CMAKE_CURRENT_BINARY_DIR}/output.yaml) + +# Ensure test input files are present in the data dir +GetInputFile(scream/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev}) + +# Ensure test input files are present in the data dir +set (TEST_INPUT_FILES + scream/mam4xx/emissions/ne2np4/surface/DMSflux.2010.ne2np4_conserv.POPmonthlyClimFromACES4BGC_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so2_surf_ne2np4_2010_clim_c20240723.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_bc_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a1_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a2_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_pom_a4_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a1_surf_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a2_surf_ne2np4_2010_clim_c20240726.nc +) +foreach (file IN ITEMS ${TEST_INPUT_FILES}) + GetInputFile(${file}) +endforeach() + +# Compare output files produced by npX tests, to ensure they are bfb +include (CompareNCFiles) + +CompareNCFilesFamilyMpi ( + TEST_BASE_NAME ${TEST_BASE_NAME} + FILE_META_NAME ${TEST_BASE_NAME}_output.INSTANT.nsteps_x1.npMPIRANKS.${RUN_T0}.nc + MPI_RANKS ${TEST_RANK_START} ${TEST_RANK_END} + LABELS mam4_srf_online_emiss physics + META_FIXTURES_REQUIRED ${FIXTURES_BASE_NAME}_npMPIRANKS_omp1 +) + +if (SCREAM_ENABLE_BASELINE_TESTS) + # Compare one of the output files with the baselines. + # Note: one is enough, since we already check that np1 is BFB with npX + set (OUT_FILE ${TEST_BASE_NAME}_output.INSTANT.nsteps_x1.np${TEST_RANK_END}.${RUN_T0}.nc) + CreateBaselineTest(${TEST_BASE_NAME} ${TEST_RANK_END} ${OUT_FILE} ${FIXTURES_BASE_NAME}) +endif() diff --git a/components/eamxx/tests/single-process/mam/emissions/input.yaml b/components/eamxx/tests/single-process/mam/emissions/input.yaml new file mode 100644 index 00000000000..e7af45178aa --- /dev/null +++ b/components/eamxx/tests/single-process/mam/emissions/input.yaml @@ -0,0 +1,49 @@ +%YAML 1.1 +--- +driver_options: + atmosphere_dag_verbosity_level: 5 + +time_stepping: + time_step: ${ATM_TIME_STEP} + run_t0: ${RUN_T0} # YYYY-MM-DD-XXXXX + number_of_steps: ${NUM_STEPS} + +atmosphere_processes: + atm_procs_list: [mam4_srf_online_emiss] + mam4_srf_online_emiss: + # MAM4xx-Surface-Emissions + srf_remap_file: "" + srf_emis_specifier_for_DMS: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/DMSflux.2010.ne2np4_conserv.POPmonthlyClimFromACES4BGC_c20240726.nc + srf_emis_specifier_for_SO2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so2_surf_ne2np4_2010_clim_c20240723.nc + srf_emis_specifier_for_bc_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_bc_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a1: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a1_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a2_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_num_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_num_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_pom_a4: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_pom_a4_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_so4_a1: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a1_surf_ne2np4_2010_clim_c20240726.nc + srf_emis_specifier_for_so4_a2: ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/surface/cmip6_mam4_so4_a2_surf_ne2np4_2010_clim_c20240726.nc + +grids_manager: + Type: Mesh Free + geo_data_source: IC_FILE + grids_names: [Physics GLL] + Physics GLL: + type: point_grid + aliases: [Physics] + number_of_global_columns: 218 + number_of_vertical_levels: 72 + +initial_conditions: + # The name of the file containing the initial conditions for this test. + Filename: ${SCREAM_DATA_DIR}/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev} + topography_filename: ${TOPO_DATA_DIR}/${EAMxx_tests_TOPO_FILE} + phis : 1.0 + #These should come from the input file + + #we should get the following variables from other processes + pbl_height : 1.0 + +# The parameters for I/O control +Scorpio: + output_yaml_files: ["output.yaml"] +... diff --git a/components/eamxx/tests/single-process/mam/emissions/output.yaml b/components/eamxx/tests/single-process/mam/emissions/output.yaml new file mode 100644 index 00000000000..433f7d58fa4 --- /dev/null +++ b/components/eamxx/tests/single-process/mam/emissions/output.yaml @@ -0,0 +1,12 @@ +%YAML 1.1 +--- +filename_prefix: mam4_srf_online_emiss_standalone_output +Averaging Type: Instant +Fields: + Physics: + Field Names: + - constituent_fluxes +output_control: + Frequency: 1 + frequency_units: nsteps +... diff --git a/components/eamxx/tests/single-process/mam/optics/CMakeLists.txt b/components/eamxx/tests/single-process/mam/optics/CMakeLists.txt index a3e6f64ec7d..0a3a0b877c5 100644 --- a/components/eamxx/tests/single-process/mam/optics/CMakeLists.txt +++ b/components/eamxx/tests/single-process/mam/optics/CMakeLists.txt @@ -40,7 +40,6 @@ set (TEST_INPUT_FILES scream/mam4xx/physprops/ocpho_rrtmg_c20240206.nc scream/mam4xx/physprops/bcpho_rrtmg_c20240206.nc scream/mam4xx/physprops/poly_rrtmg_c20240206.nc - scream/init/scream_unit_tests_aerosol_optics_ne2np4L72_20220822.nc ) foreach (file IN ITEMS ${TEST_INPUT_FILES}) @@ -63,4 +62,4 @@ if (SCREAM_ENABLE_BASELINE_TESTS) # Note: one is enough, since we already check that np1 is BFB with npX set (OUT_FILE ${TEST_BASE_NAME}_output.INSTANT.nsteps_x2.np${TEST_RANK_END}.${RUN_T0}.nc) CreateBaselineTest(${TEST_BASE_NAME} ${TEST_RANK_END} ${OUT_FILE} ${FIXTURES_BASE_NAME}) -endif() \ No newline at end of file +endif() diff --git a/externals/mam4xx b/externals/mam4xx index 5495c27a6df..2eee1988b1e 160000 --- a/externals/mam4xx +++ b/externals/mam4xx @@ -1 +1 @@ -Subproject commit 5495c27a6df7ae4e958d67077abc26d4c7e5765d +Subproject commit 2eee1988b1e6488c904e0a28564bdcadfa190d34