diff --git a/components/eamxx/cime_config/namelist_defaults_scream.xml b/components/eamxx/cime_config/namelist_defaults_scream.xml index 9a7975f01ba..a1435b03de1 100644 --- a/components/eamxx/cime_config/namelist_defaults_scream.xml +++ b/components/eamxx/cime_config/namelist_defaults_scream.xml @@ -268,6 +268,64 @@ be lost if SCREAM_HACK_XML is not enabled. + + + + + true + true + true + true + + 172800.0 + 3.0E-008 + 4 + 193.0 + 20100101 + ${DIN_LOC_ROOT}/atm/scream/mam4xx/linoz/ne30pg2/linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne30pg2_c20240724.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/linoz/ne4pg2/linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne4pg2_c20240724.nc + + 20150101 + ${DIN_LOC_ROOT}/atm/scream/mam4xx/invariants/ne30pg2/oxid_1.9x2.5_L26_1850-2015_ne30pg2_c20241009.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/invariants/ne4pg2/oxid_1.9x2.5_L26_1850-2015_ne4pg2_c20241009.nc + 20100101 + ${DIN_LOC_ROOT}/atm/scream/mam4xx/linoz/Linoz_Chlorine_Loading_CMIP6_0003-2017_c20171114.nc + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/photolysis/RSF_GT200nm_v3.0_c080811.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/photolysis/temp_prs_GT200nm_JPL10_c130206.nc + + 20100101 + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_so2_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_so4_a1_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_so4_a2_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_pom_a4_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_bc_a4_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_num_a1_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_num_a2_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_num_a4_elev_1x1_2010_clim_ne30pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne30pg2/elevated/cmip6_mam4_soag_elev_1x1_2010_clim_ne30pg2_c20241008.nc + + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_so2_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_so4_a1_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_so4_a2_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_pom_a4_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_bc_a4_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_num_a1_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_num_a2_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_num_a4_elev_1x1_2010_clim_ne4pg2_c20241008.nc + ${DIN_LOC_ROOT}/atm/scream/mam4xx/emissions/ne4pg2/elevated/cmip6_mam4_soag_elev_1x1_2010_clim_ne4pg2_c20241008.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 + + + ${DIN_LOC_ROOT}/atm/scream/mam4xx/physprops/mam4_mode1_rrtmg_aeronetdust_c20240206.nc @@ -290,7 +348,7 @@ be lost if SCREAM_HACK_XML is not enabled. - + ${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 @@ -300,8 +358,8 @@ be lost if SCREAM_HACK_XML is not enabled. ${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/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 @@ -311,7 +369,7 @@ be lost if SCREAM_HACK_XML is not enabled. ${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/ne4pg2/surface/DMSflux.2010.ne4pg2_conserv.POPmonthlyClimFromACES4BGC_c20240814.nc @@ -576,9 +634,9 @@ be lost if SCREAM_HACK_XML is not enabled. 0.0 0.0,0.0 - 2.6e-08 - 0.41417721820867320E-007 - 0.15100083211582764E+004 + 1.37146e-07 ,3.45899e-08 ,1.00000e-06 ,9.99601e-08 + 1.37452e-07 ,3.46684e-08 ,1.00900e-06 ,9.99601e-08 + 5.08262e-12 ,1.54035e-13 ,3.09018e-13 ,9.14710e-22 0.0 0.0 0.0 @@ -646,7 +704,7 @@ be lost if SCREAM_HACK_XML is not enabled. + --> diff --git a/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/aero_microphysics/shell_commands b/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/aero_microphysics/shell_commands new file mode 100644 index 00000000000..1d6757a5bd9 --- /dev/null +++ b/components/eamxx/cime_config/testdefs/testmods_dirs/scream/mam4xx/aero_microphysics/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_aero_microphys" -b + + + diff --git a/components/eamxx/src/physics/mam/CMakeLists.txt b/components/eamxx/src/physics/mam/CMakeLists.txt index 9874f79f9eb..fb6a48733f8 100644 --- a/components/eamxx/src/physics/mam/CMakeLists.txt +++ b/components/eamxx/src/physics/mam/CMakeLists.txt @@ -42,6 +42,7 @@ add_subdirectory(${EXTERNALS_SOURCE_DIR}/mam4xx ${CMAKE_BINARY_DIR}/externals/ma # EAMxx mam4xx-based atmospheric processes add_library(mam eamxx_mam_microphysics_process_interface.cpp + ${SCREAM_BASE_DIR}/src/physics/rrtmgp/shr_orb_mod_c2f.F90 eamxx_mam_optics_process_interface.cpp eamxx_mam_dry_deposition_process_interface.cpp eamxx_mam_aci_process_interface.cpp @@ -54,7 +55,7 @@ target_include_directories(mam PUBLIC ${EXTERNALS_SOURCE_DIR}/haero ${EXTERNALS_SOURCE_DIR}/mam4xx/src ) -target_link_libraries(mam PUBLIC physics_share scream_share mam4xx haero) +target_link_libraries(mam PUBLIC physics_share csm_share scream_share mam4xx haero) #if (NOT SCREAM_LIB_ONLY) # add_subdirectory(tests) 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 f2d58c54e44..449ba4a55d8 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 @@ -63,26 +63,22 @@ void MAMAci::set_grids( // Define the different field layouts that will be used for this process using namespace ShortFieldTagsNames; - // Layout for 3D (2d horiz X 1d vertical) variables - // mid points - FieldLayout scalar3d_layout_mid{{COL, LEV}, {ncol_, nlev_}}; + // Layout for 2D (2d horiz) variable + const FieldLayout scalar2d = grid_->get_2d_scalar_layout(); + + // Layout for 3D (2d horiz X 1d vertical) variable defined at mid-level and // interfaces - FieldLayout scalar3d_layout_int{{COL, ILEV}, {ncol_, nlev_ + 1}}; + const FieldLayout scalar3d_mid = grid_->get_3d_scalar_layout(true); + const FieldLayout scalar3d_int = grid_->get_3d_scalar_layout(false); // layout for 2D (1d horiz X 1d vertical) variable FieldLayout scalar2d_layout_col{{COL}, {ncol_}}; - auto make_layout = [](const std::vector &extents, - const std::vector &names) { - std::vector tags(extents.size(), CMP); - return FieldLayout(tags, extents, names); - }; - 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 + constexpr auto q_unit = kg / kg; // units of mass mixing ratios of tracers + constexpr auto n_unit = 1 / kg; // units of number mixing ratios of tracers - auto nondim = ekat::units::Units::nondimensional(); + constexpr auto nondim = ekat::units::Units::nondimensional(); // atmospheric quantities // specific humidity [kg/kg] @@ -101,28 +97,28 @@ void MAMAci::set_grids( add_tracer("ni", grid_, n_unit); // Temperature[K] at midpoints - add_field("T_mid", scalar3d_layout_mid, K, grid_name); + add_field("T_mid", scalar3d_mid, K, grid_name); // Vertical pressure velocity [Pa/s] at midpoints - add_field("omega", scalar3d_layout_mid, Pa / s, grid_name); + add_field("omega", scalar3d_mid, Pa / s, grid_name); // Total pressure [Pa] at midpoints - add_field("p_mid", scalar3d_layout_mid, Pa, grid_name); + add_field("p_mid", scalar3d_mid, Pa, grid_name); // Total pressure [Pa] at interfaces - add_field("p_int", scalar3d_layout_int, Pa, grid_name); + add_field("p_int", scalar3d_int, Pa, grid_name); // Layer thickness(pdel) [Pa] at midpoints - add_field("pseudo_density", scalar3d_layout_mid, Pa, grid_name); + add_field("pseudo_density", scalar3d_mid, Pa, grid_name); // planetary boundary layer height add_field("pbl_height", scalar2d_layout_col, m, grid_name); // cloud fraction [nondimensional] computed by eamxx_cld_fraction_process - add_field("cldfrac_tot", scalar3d_layout_mid, nondim, grid_name); + add_field("cldfrac_tot", scalar3d_mid, nondim, grid_name); - auto m2 = pow(m, 2); - auto s2 = pow(s, 2); + constexpr auto m2 = pow(m, 2); + constexpr auto s2 = pow(s, 2); // NOTE: w_variance im microp_aero.F90 in EAM is at "itim_old" dynamics time // step. Since, we are using SE dycore, itim_old is 1 which is equivalent to @@ -130,29 +126,30 @@ void MAMAci::set_grids( // and we might need to revisit this. // Vertical velocity variance at midpoints - add_field("w_variance", scalar3d_layout_mid, m2 / s2, grid_name); + add_field("w_variance", scalar3d_mid, m2 / s2, grid_name); // NOTE: "cldfrac_liq" is updated in SHOC. "cldfrac_liq" in C++ code is // equivalent to "alst" in the shoc_intr.F90. In the C++ code, it is used as // "shoc_cldfrac" and in the F90 code it is called "cloud_frac" // Liquid stratiform cloud fraction at midpoints - add_field("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name); + add_field("cldfrac_liq", scalar3d_mid, nondim, grid_name); // Previous value of liquid stratiform cloud fraction at midpoints - add_field("cldfrac_liq_prev", scalar3d_layout_mid, nondim, - grid_name); + add_field("cldfrac_liq_prev", scalar3d_mid, nondim, grid_name); // Eddy diffusivity for heat - add_field("eddy_diff_heat", scalar3d_layout_mid, m2 / s, grid_name); + add_field("eddy_diff_heat", scalar3d_mid, m2 / s, grid_name); - // Layout for 4D (2d horiz X 1d vertical x number of modes) variables - const int num_aero_modes = mam_coupling::num_aero_modes(); - FieldLayout scalar4d_layout_mid = - make_layout({ncol_, num_aero_modes, nlev_}, {"COL", "NMODES", "LEV"}); + // Number of modes + constexpr int nmodes = mam4::AeroConfig::num_modes(); + + // layout for 3D (ncol, nmodes, nlevs) + FieldLayout scalar3d_mid_nmodes = + grid_->get_3d_vector_layout(true, nmodes, "nmodes"); // dry diameter of aerosols [m] - add_field("dgnum", scalar4d_layout_mid, m, grid_name); + add_field("dgnum", scalar3d_mid_nmodes, m, grid_name); // ======================================================================== // Output from this whole process @@ -171,8 +168,7 @@ void MAMAci::set_grids( // NOT advected const char *cld_nmr_field_name = mam_coupling::cld_aero_nmr_field_name(mode); - add_field(cld_nmr_field_name, scalar3d_layout_mid, n_unit, - grid_name); + 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 @@ -187,8 +183,7 @@ void MAMAci::set_grids( const char *cld_mmr_field_name = mam_coupling::cld_aero_mmr_field_name(mode, a); if(strlen(cld_mmr_field_name) > 0) { - add_field(cld_mmr_field_name, scalar3d_layout_mid, q_unit, - grid_name); + add_field(cld_mmr_field_name, scalar3d_mid, q_unit, grid_name); } } // end for loop num species } // end for loop for num modes @@ -203,54 +198,52 @@ void MAMAci::set_grids( // ------------------------------------------------------------------------ // number of activated aerosol for ice nucleation[#/kg] - add_field("ni_activated", scalar3d_layout_mid, n_unit, grid_name); + add_field("ni_activated", scalar3d_mid, n_unit, grid_name); // ------------------------------------------------------------------------ // Output from droplet activation process (dropmixnuc) // ------------------------------------------------------------------------ // tendency in droplet number mixing ratio [#/kg/s] - add_field("nc_nuceat_tend", scalar3d_layout_mid, n_unit / s, - grid_name); + add_field("nc_nuceat_tend", scalar3d_mid, n_unit / s, grid_name); // FIXME: [TEMPORARY]droplet number mixing ratio source tendency [#/kg/s] - add_field("nsource", scalar3d_layout_mid, n_unit / s, grid_name); + add_field("nsource", scalar3d_mid, n_unit / s, grid_name); // FIXME: [TEMPORARY]droplet number mixing ratio tendency due to mixing // [#/kg/s] - add_field("ndropmix", scalar3d_layout_mid, n_unit / s, grid_name); + add_field("ndropmix", scalar3d_mid, n_unit / s, grid_name); // FIXME: [TEMPORARY]droplet number as seen by ACI [#/kg] - add_field("nc_inp_to_aci", scalar3d_layout_mid, n_unit / s, - grid_name); - const auto cm_tmp = m / 100; // FIXME: [TEMPORARY] remove this - const auto cm3 = cm_tmp * cm_tmp * cm_tmp; // FIXME: [TEMPORARY] remove this + add_field("nc_inp_to_aci", scalar3d_mid, n_unit / s, grid_name); + constexpr auto cm_tmp = m / 100; // FIXME: [TEMPORARY] remove this + constexpr auto cm3 = pow(cm_tmp, 3); // FIXME: [TEMPORARY] remove this // FIXME: [TEMPORARY] remove the following ccn outputs - add_field("ccn_0p02", scalar3d_layout_mid, cm3, grid_name); - add_field("ccn_0p05", scalar3d_layout_mid, cm3, grid_name); - add_field("ccn_0p1", scalar3d_layout_mid, cm3, grid_name); - add_field("ccn_0p2", scalar3d_layout_mid, cm3, grid_name); - add_field("ccn_0p5", scalar3d_layout_mid, cm3, grid_name); - add_field("ccn_1p0", scalar3d_layout_mid, cm3, grid_name); + add_field("ccn_0p02", scalar3d_mid, cm3, grid_name); + add_field("ccn_0p05", scalar3d_mid, cm3, grid_name); + add_field("ccn_0p1", scalar3d_mid, cm3, grid_name); + add_field("ccn_0p2", scalar3d_mid, cm3, grid_name); + add_field("ccn_0p5", scalar3d_mid, cm3, grid_name); + add_field("ccn_1p0", scalar3d_mid, cm3, grid_name); // ------------------------------------------------------------------------ // Output from hetrozenous freezing // ------------------------------------------------------------------------ - const auto cm = m / 100; + constexpr auto cm = m / 100; // units of number mixing ratios of tracers - auto frz_unit = 1 / (cm * cm * cm * s); + constexpr auto frz_unit = 1 / (cm * cm * cm * s); // heterogeneous freezing by immersion nucleation [cm^-3 s^-1] - add_field("hetfrz_immersion_nucleation_tend", scalar3d_layout_mid, + add_field("hetfrz_immersion_nucleation_tend", scalar3d_mid, frz_unit, grid_name); // heterogeneous freezing by contact nucleation [cm^-3 s^-1] - add_field("hetfrz_contact_nucleation_tend", scalar3d_layout_mid, - frz_unit, grid_name); + add_field("hetfrz_contact_nucleation_tend", scalar3d_mid, frz_unit, + grid_name); // heterogeneous freezing by deposition nucleation [cm^-3 s^-1] - add_field("hetfrz_deposition_nucleation_tend", scalar3d_layout_mid, + add_field("hetfrz_deposition_nucleation_tend", scalar3d_mid, frz_unit, grid_name); } // function set_grids ends @@ -299,11 +292,11 @@ void MAMAci::initialize_impl(const RunType run_type) { // store rest fo the atm fields in dry_atm_in dry_atm_.z_surf = 0; - 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_int = get_field_in("p_int").get_view(); - dry_atm_.p_del = get_field_in("pseudo_density").get_view(); - dry_atm_.omega = get_field_in("omega").get_view(); + 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_int = get_field_in("p_int").get_view(); + dry_atm_.p_del = get_field_in("pseudo_density").get_view(); + dry_atm_.omega = get_field_in("omega").get_view(); // store fields converted to dry mmr from wet mmr in dry_atm_ dry_atm_.qv = buffer_.qv_dry; @@ -561,7 +554,8 @@ void MAMAci::run_impl(const double dt) { // output w0_, rho_); - compute_tke_at_interfaces(team_policy, w_sec_mid_, dry_atm_.dz, nlev_, w_sec_int_, + compute_tke_at_interfaces(team_policy, w_sec_mid_, dry_atm_.dz, nlev_, + w_sec_int_, // output tke_); @@ -594,25 +588,26 @@ void MAMAci::run_impl(const double dt) { // output cloud_frac_, cloud_frac_prev_); - mam_coupling::compute_recipical_pseudo_density(team_policy, dry_atm_.p_del, nlev_, - // output - rpdel_); + mam_coupling::compute_recipical_pseudo_density(team_policy, dry_atm_.p_del, + nlev_, + // output + rpdel_); Kokkos::fence(); // wait for rpdel_ to be computed. // Compute activated CCN number tendency (tendnd_) and updated // cloud borne aerosols (stored in a work array) and interstitial // aerosols tendencies - call_function_dropmixnuc(team_policy, dt, dry_atm_, rpdel_, kvh_mid_, kvh_int_, wsub_, - cloud_frac_, cloud_frac_prev_, dry_aero_, nlev_, - // output - coltend_, coltend_cw_, qcld_, ndropcol_, ndropmix_, - nsource_, wtke_, ccn_, - // ## output to be used by the other processes ## - qqcw_fld_work_, ptend_q_, factnum_, tendnd_, - // work arrays - raercol_cw_, raercol_, state_q_work_, nact_, mact_, - dropmixnuc_scratch_mem_); + call_function_dropmixnuc( + team_policy, dt, dry_atm_, rpdel_, kvh_mid_, kvh_int_, wsub_, cloud_frac_, + cloud_frac_prev_, dry_aero_, nlev_, + // output + coltend_, coltend_cw_, qcld_, ndropcol_, ndropmix_, nsource_, wtke_, ccn_, + // ## output to be used by the other processes ## + qqcw_fld_work_, ptend_q_, factnum_, tendnd_, + // work arrays + raercol_cw_, raercol_, state_q_work_, nact_, mact_, + dropmixnuc_scratch_mem_); Kokkos::fence(); // wait for ptend_q_ to be computed. Kokkos::deep_copy(ccn_0p02_, 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 090fdfcb731..b0ab232c255 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 @@ -1,215 +1,403 @@ #include -#include -#include -#include -#include "scream_config.h" // for SCREAM_CIME_BUILD +// impl namespace for some driver level functions for microphysics -#include // for serial NetCDF file reads on MPI root -#include +#include "readfiles/photo_table_utils.cpp" +#include "physics/rrtmgp/shr_orb_mod_c2f.hpp" -// NOTE: see the impl/ directory for the contents of the impl namespace -#include "impl/compute_o3_column_density.cpp" -#include "impl/compute_water_content.cpp" -#include "impl/gas_phase_chemistry.cpp" +namespace scream { -namespace scream -{ +MAMMicrophysics::MAMMicrophysics(const ekat::Comm &comm, + const ekat::ParameterList ¶ms) + : AtmosphereProcess(comm, params), aero_config_() { + config_.amicphys.do_cond = m_params.get("mam4_do_cond"); + config_.amicphys.do_rename = m_params.get("mam4_do_rename"); + config_.amicphys.do_newnuc = m_params.get("mam4_do_newnuc"); + config_.amicphys.do_coag = m_params.get("mam4_do_coag"); -MAMMicrophysics::MAMMicrophysics( - const ekat::Comm& comm, - const ekat::ParameterList& params) - : AtmosphereProcess(comm, params), - aero_config_() { - configure(params); + // these parameters guide the coupling between parameterizations + // NOTE: mam4xx was ported with these parameters fixed, so it's probably not + // NOTE: safe to change these without code modifications. + + // controls treatment of h2so4 condensation in mam_gasaerexch_1subarea + // 1 = sequential calc. of gas-chem prod then condensation loss + // 2 = simultaneous calc. of gas-chem prod and condensation loss + config_.amicphys.gaexch_h2so4_uptake_optaa = 2; + + // controls treatment of h2so4 concentrationin mam_newnuc_1subarea + // 1 = use average value calculated in standard cam5.2.10 and earlier + // 2 = use average value calculated in mam_gasaerexch_1subarea + // 11 = use average of initial and final values from mam_gasaerexch_1subarea + // 12 = use final value from mam_gasaerexch_1subarea + config_.amicphys.newnuc_h2so4_conc_optaa = 2; + + // LINOZ namelist parameters + config_.linoz.o3_lbl = m_params.get("mam4_o3_lbl"); + config_.linoz.o3_tau = m_params.get("mam4_o3_tau"); + config_.linoz.o3_sfc = m_params.get("mam4_o3_sfc"); + config_.linoz.psc_T = m_params.get("mam4_psc_T"); } AtmosphereProcessType MAMMicrophysics::type() const { return AtmosphereProcessType::Physics; } -std::string MAMMicrophysics::name() const { - return "mam4_micro"; -} +// ================================================================ +// SET_GRIDS +// ================================================================ +void MAMMicrophysics::set_grids( + const std::shared_ptr grids_manager) { + // set grid for all the inputs and outputs + // use physics grid + grid_ = grids_manager->get_grid("Physics"); + const auto &grid_name = grid_->name(); -namespace { + ncol_ = grid_->get_num_local_dofs(); // number of columns on this rank + nlev_ = grid_->get_num_vertical_levels(); // number of levels per column -void set_data_file(const char *name, const char *path, char location[MAX_FILENAME_LEN]) { - EKAT_REQUIRE_MSG(strlen(SCREAM_DATA_DIR) + strlen(path) < MAX_FILENAME_LEN, - "Error! " << name << " path is too long (must be < " << MAX_FILENAME_LEN << " characters)"); - sprintf(location, "%s/%s", SCREAM_DATA_DIR, path); -} + // get column geometry and locations + col_latitudes_ = grid_->get_geometry_data("lat").get_view(); -} + // define the different field layouts that will be used for this process + using namespace ShortFieldTagsNames; -#define set_file_location(data_file, path) set_data_file(#data_file, path, data_file) + // Layout for 2D (2d horiz) variable + const FieldLayout scalar2d = grid_->get_2d_scalar_layout(); -void MAMMicrophysics::set_defaults_() { - config_.amicphys.do_cond = true; - config_.amicphys.do_rename = true; - config_.amicphys.do_newnuc = true; - config_.amicphys.do_coag = true; + // Layout for 3D (2d horiz X 1d vertical) variable defined at mid-level and + // interfaces + const FieldLayout scalar3d_mid = grid_->get_3d_scalar_layout(true); + const FieldLayout scalar3d_int = grid_->get_3d_scalar_layout(false); - config_.amicphys.nucleation = {}; - config_.amicphys.nucleation.dens_so4a_host = 1770.0; - config_.amicphys.nucleation.mw_so4a_host = 115.0; - config_.amicphys.nucleation.newnuc_method_user_choice = 2; - config_.amicphys.nucleation.pbl_nuc_wang2008_user_choice = 1; - config_.amicphys.nucleation.adjust_factor_pbl_ratenucl = 1.0; - config_.amicphys.nucleation.accom_coef_h2so4 = 1.0; - config_.amicphys.nucleation.newnuc_adjust_factor_dnaitdt = 1.0; + using namespace ekat::units; + constexpr auto q_unit = kg / kg; // units of mass mixing ratios of tracers + constexpr auto n_unit = 1 / kg; // units of number mixing ratios of tracers - // these parameters guide the coupling between parameterizations - // NOTE: mam4xx was ported with these parameters fixed, so it's probably not - // NOTE: safe to change these without code modifications. - config_.amicphys.gaexch_h2so4_uptake_optaa = 2; - config_.amicphys.newnuc_h2so4_conc_optaa = 2; + // -------------------------------------------------------------------------- + // These variables are "Required" or pure inputs for the process + // -------------------------------------------------------------------------- - //=========================================================== - // default data file locations (relative to SCREAM_DATA_DIR) - //=========================================================== + // ----------- Atmospheric quantities ------------- - // many of these paths were extracted from - // e3smv2/bld/namelist_files/namelist_defaults_eam.xml + // Specific humidity [kg/kg](Require only for building DS) + add_tracer("qv", grid_, kg/kg); // specific humidity - // photolysis - set_file_location(config_.photolysis.rsf_file, "../waccm/phot/RSF_GT200nm_v3.0_c080811.nc"); - set_file_location(config_.photolysis.xs_long_file, "../waccm/phot/temp_prs_GT200nm_JPL10_c130206.nc"); + // Cloud liquid mass mixing ratio [kg/kg](Require only for building DS) + add_tracer("qc", grid_, kg/kg); // cloud liquid wet mixing ratio - // stratospheric chemistry - set_file_location(config_.linoz.chlorine_loading_file, "../cam/chem/trop_mozart/ub/Linoz_Chlorine_Loading_CMIP6_0003-2017_c20171114.nc"); + // Cloud ice mass mixing ratio [kg/kg](Require only for building DS) + add_tracer("qi", grid_, kg/kg); // ice wet mixing ratio - // dry deposition - set_file_location(config_.drydep.srf_file, "../cam/chem/trop_mam/atmsrf_ne4pg2_200527.nc"); -} + // Cloud liquid number mixing ratio [1/kg](Require only for building DS) + add_tracer("nc", grid_, n_unit); // cloud liquid wet number mixing ratio -void MAMMicrophysics::configure(const ekat::ParameterList& params) { - set_defaults_(); - // FIXME: implement "namelist" parsing -} + // Cloud ice number mixing ratio [1/kg](Require only for building DS) + add_tracer("ni", grid_, n_unit); // ice number mixing ratio -void MAMMicrophysics::set_grids(const std::shared_ptr grids_manager) { - using namespace ekat::units; + // Temperature[K] at midpoints + add_field("T_mid", scalar3d_mid, K, grid_name); - Units nondim = Units::nondimensional(); - Units n_unit (1/kg,"#/kg"); // number mixing ratios [# / kg air] - const auto m2 = pow(m,2); - const auto s2 = pow(s,2); + // Vertical pressure velocity [Pa/s] at midpoints (Require only for building + // DS) + add_field("omega", scalar3d_mid, Pa / s, grid_name); - grid_ = grids_manager->get_grid("Physics"); - const auto& grid_name = grid_->name(); + // Total pressure [Pa] at midpoints + add_field("p_mid", scalar3d_mid, Pa, 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 + // Total pressure [Pa] at interfaces + add_field("p_int", scalar3d_int, Pa, grid_name); - // get column geometry and locations - col_areas_ = grid_->get_geometry_data("area").get_view(); - col_latitudes_ = grid_->get_geometry_data("lat").get_view(); - col_longitudes_ = grid_->get_geometry_data("lon").get_view(); + // Layer thickness(pdel) [Pa] at midpoints + add_field("pseudo_density", scalar3d_mid, Pa, grid_name); - // define the different field layouts that will be used for this process - using namespace ShortFieldTagsNames; + // Planetary boundary layer height [m] + add_field("pbl_height", scalar2d, m, grid_name); - // layout for 2D (1d horiz X 1d vertical) variable - FieldLayout scalar2d_layout_col{ {COL}, {ncol_} }; + constexpr auto m2 = pow(m, 2); + constexpr auto s2 = pow(s, 2); - // layout for 3D (2d horiz X 1d vertical) variables - FieldLayout scalar3d_layout_mid{ {COL, LEV}, {ncol_, nlev_} }; + // Surface geopotential [m2/s2] + add_field("phis", scalar2d, m2 / s2, grid_name); - // define fields needed in mam4xx + //----------- Variables from microphysics scheme ------------- + constexpr auto nondim = ekat::units::Units::nondimensional(); + // Total cloud fraction [fraction] + add_field("cldfrac_liq", scalar3d_mid, nondim, grid_name); - // atmospheric quantities - add_field("omega", scalar3d_layout_mid, Pa/s, grid_name); // vertical pressure velocity - add_field("T_mid", scalar3d_layout_mid, K, grid_name); // Temperature - add_field("p_mid", scalar3d_layout_mid, Pa, grid_name); // total pressure - add_field("pbl_height", scalar2d_layout_col, m, grid_name); // planetary boundary layer height - add_field("pseudo_density", scalar3d_layout_mid, Pa, grid_name); // p_del, hydrostatic pressure - add_field("phis", scalar2d_layout_col, m2/s2, grid_name); - add_field("cldfrac_tot", scalar3d_layout_mid, nondim, grid_name); // cloud fraction - add_tracer("qv", grid_, kg/kg); // specific humidity - add_tracer("qi", grid_, kg/kg); // ice wet mixing ratio - add_tracer("ni", grid_, n_unit); // ice number mixing ratio + //----------- Variables from other mam4xx processes ------------ + // Number of modes + constexpr int nmodes = mam4::AeroConfig::num_modes(); - // droplet activation can alter cloud liquid and number mixing ratios - add_tracer("qc", grid_, kg/kg); // cloud liquid wet mixing ratio - add_tracer("nc", grid_, n_unit); // cloud liquid wet number mixing ratio + // layout for 3D (ncol, nmodes, nlevs) + FieldLayout scalar3d_mid_nmodes = + grid_->get_3d_vector_layout(true, nmodes, "nmodes"); + + // Geometric mean dry diameter for number distribution [m] + add_field("dgnum", scalar3d_mid_nmodes, m, grid_name); + // Geometric mean wet diameter for number distribution [m] + add_field("dgnumwet", scalar3d_mid_nmodes, m, grid_name); + + constexpr auto m3 = pow(m, 3); + // Wet density of interstitial aerosol [kg/m3] + add_field("wetdens", scalar3d_mid_nmodes, kg / m3, grid_name); + + //----------- Variables from coupler (land component)--------- + // surface albedo shortwave, direct + add_field("sfc_alb_dir_vis", scalar2d, nondim, grid_name); + + // --------------------------------------------------------------------- + // These variables are "updated" or inputs/outputs for the process + // --------------------------------------------------------------------- + + // (interstitial) aerosol tracers of interest: mass (q) and number (n) mixing + // ratios + for(int m = 0; m < nmodes; ++m) { + const char *int_nmr_field_name = mam_coupling::int_aero_nmr_field_name(m); - // (interstitial) aerosol tracers of interest: mass (q) and number (n) mixing ratios - for (int m = 0; m < mam_coupling::num_aero_modes(); ++m) { - const char* int_nmr_field_name = mam_coupling::int_aero_nmr_field_name(m); add_tracer(int_nmr_field_name, grid_, n_unit); - for (int a = 0; a < mam_coupling::num_aero_species(); ++a) { - const char* int_mmr_field_name = mam_coupling::int_aero_mmr_field_name(m, a); - if (strlen(int_mmr_field_name) > 0) { + for(int a = 0; a < mam_coupling::num_aero_species(); ++a) { + const char *int_mmr_field_name = + mam_coupling::int_aero_mmr_field_name(m, a); + + if(strlen(int_mmr_field_name) > 0) { add_tracer(int_mmr_field_name, grid_, kg/kg); } - } - } + } // for loop species + } // for loop nmodes interstitial + // (cloud) aerosol tracers of interest: mass (q) and number (n) mixing ratios + for(int m = 0; m < nmodes; ++m) { + const char *cld_nmr_field_name = mam_coupling::cld_aero_nmr_field_name(m); + + add_field(cld_nmr_field_name, scalar3d_mid, n_unit, grid_name); + for(int a = 0; a < mam_coupling::num_aero_species(); ++a) { + const char *cld_mmr_field_name = + mam_coupling::cld_aero_mmr_field_name(m, a); + + if(strlen(cld_mmr_field_name) > 0) { + add_field(cld_mmr_field_name, scalar3d_mid, q_unit, grid_name); + } + } // for loop species + } // for loop nmodes cld borne // aerosol-related gases: mass mixing ratios - for (int g = 0; g < mam_coupling::num_aero_gases(); ++g) { - const char* gas_mmr_field_name = mam_coupling::gas_mmr_field_name(g); + for(int g = 0; g < mam_coupling::num_aero_gases(); ++g) { + const char *gas_mmr_field_name = mam_coupling::gas_mmr_field_name(g); add_tracer(gas_mmr_field_name, grid_, kg/kg); } - // Tracers group -- do we need this in addition to the tracers above? In any - // case, this call should be idempotent, so it can't hurt. - add_group("tracers", grid_name, 1, Bundling::Required); -} - -// this checks whether we have the tracers we expect -void MAMMicrophysics:: -set_computed_group_impl(const FieldGroup& group) { - const auto& name = group.m_info->m_group_name; - EKAT_REQUIRE_MSG(name=="tracers", - "Error! MAM4 expects a 'tracers' field group (got '" << name << "')\n"); - - EKAT_REQUIRE_MSG(group.m_info->m_bundled, - "Error! MAM4 expects bundled fields for tracers.\n"); - - // how many aerosol/gas tracers do we expect? - int num_tracers = 2 * (mam_coupling::num_aero_modes() + - mam_coupling::num_aero_tracers()) + - mam_coupling::num_aero_gases(); - EKAT_REQUIRE_MSG(group.m_info->size() >= num_tracers, - "Error! MAM4 requires at least " << num_tracers << " aerosol tracers."); -} - -size_t MAMMicrophysics::requested_buffer_size_in_bytes() const -{ + // Creating a Linoz reader and setting Linoz parameters involves reading data + // from a file and configuring the necessary parameters for the Linoz model. + { + linoz_file_name_ = m_params.get("mam4_linoz_file_name"); + const std::string linoz_map_file = + m_params.get("aero_microphys_remap_file", ""); + const std::vector var_names{ + "o3_clim", "o3col_clim", "t_clim", "PmL_clim", + "dPmL_dO3", "dPmL_dT", "dPmL_dO3col", "cariolle_pscs"}; + + // in format YYYYMMDD + const int linoz_cyclical_ymd = m_params.get("mam4_linoz_ymd"); + scream::mam_coupling::setup_tracer_data(linoz_data_, linoz_file_name_, + linoz_cyclical_ymd); + LinozHorizInterp_ = scream::mam_coupling::create_horiz_remapper( + grid_, linoz_file_name_, linoz_map_file, var_names, linoz_data_); + LinozDataReader_ = scream::mam_coupling::create_tracer_data_reader( + LinozHorizInterp_, linoz_file_name_); + + // linoz reader + const auto io_grid_linoz = LinozHorizInterp_->get_src_grid(); + const int num_cols_io_linoz = + io_grid_linoz->get_num_local_dofs(); // Number of columns on this rank + const int num_levs_io_linoz = + io_grid_linoz + ->get_num_vertical_levels(); // Number of levels per column + const int nvars = int(var_names.size()); + linoz_data_.init(num_cols_io_linoz, num_levs_io_linoz, nvars); + linoz_data_.allocate_temporal_views(); + } // LINOZ reader + + { + oxid_file_name_ = m_params.get("mam4_oxid_file_name"); + const std::string oxid_map_file = + m_params.get("aero_microphys_remap_file", ""); + // NOTE: order matches mam4xx: + const std::vector var_names{"O3", "OH", "NO3", "HO2"}; + + // in format YYYYMMDD + const int oxid_ymd = m_params.get("mam4_oxid_ymd"); + scream::mam_coupling::setup_tracer_data(tracer_data_, oxid_file_name_, + oxid_ymd); + TracerHorizInterp_ = scream::mam_coupling::create_horiz_remapper( + grid_, oxid_file_name_, oxid_map_file, var_names, tracer_data_); + TracerDataReader_ = scream::mam_coupling::create_tracer_data_reader( + TracerHorizInterp_, oxid_file_name_); + + const int nvars = int(var_names.size()); + const auto io_grid = TracerHorizInterp_->get_src_grid(); + const int num_cols_io = + io_grid->get_num_local_dofs(); // Number of columns on this rank + const int num_levs_io = + io_grid->get_num_vertical_levels(); // Number of levels per column + tracer_data_.init(num_cols_io, num_levs_io, nvars); + tracer_data_.allocate_temporal_views(); + + for(int ivar = 0; ivar < nvars; ++ivar) { + cnst_offline_[ivar] = view_2d("cnst_offline_", ncol_, nlev_); + } + } // oxid file reader + + { + const std::string extfrc_map_file = + m_params.get("aero_microphys_remap_file", ""); + // NOTE: order of forcing species is important. + // extfrc_lst(: 9) = {'SO2 ','so4_a1 ','so4_a2 + // ','pom_a4 ','bc_a4 ', 'num_a1 ','num_a2 + // ','num_a4 ','SOAG ' } + // This order corresponds to files in namelist e3smv2 + extfrc_lst_ = {"so2", "so4_a1", "so4_a2", "pom_a4", "bc_a4", + "num_a1", "num_a2", "num_a4", "soag"}; + + for(const auto &var_name : extfrc_lst_) { + std::string item_name = "mam4_" + var_name + "_verti_emiss_file_name"; + const auto file_name = m_params.get(item_name); + vert_emis_file_name_[var_name] = file_name; + } + vert_emis_var_names_["so2"] = {"BB", "ENE_ELEV", "IND_ELEV", "contvolc"}; + vert_emis_var_names_["so4_a1"] = {"BB", "ENE_ELEV", "IND_ELEV", "contvolc"}; + vert_emis_var_names_["so4_a2"] = {"contvolc"}; + vert_emis_var_names_["pom_a4"] = {"BB"}; + vert_emis_var_names_["bc_a4"] = {"BB"}; + vert_emis_var_names_["num_a1"] = { + "num_a1_SO4_ELEV_BB", "num_a1_SO4_ELEV_ENE", "num_a1_SO4_ELEV_IND", + "num_a1_SO4_ELEV_contvolc"}; + vert_emis_var_names_["num_a2"] = {"num_a2_SO4_ELEV_contvolc"}; + // num_a4 + // FIXME: why the sectors in this files are num_a1; + // I guess this should be num_a4? Is this a bug in the orginal nc files? + vert_emis_var_names_["num_a4"] = {"num_a1_BC_ELEV_BB", + "num_a1_POM_ELEV_BB"}; + vert_emis_var_names_["soag"] = {"SOAbb_src", "SOAbg_src", "SOAff_src"}; + + int verti_emiss_cyclical_ymd = m_params.get("verti_emiss_ymd"); + + for(const auto &var_name : extfrc_lst_) { + const auto file_name = vert_emis_file_name_[var_name]; + const auto var_names = vert_emis_var_names_[var_name]; + + scream::mam_coupling::TracerData data_tracer; + scream::mam_coupling::setup_tracer_data(data_tracer, file_name, + verti_emiss_cyclical_ymd); + auto hor_rem = scream::mam_coupling::create_horiz_remapper( + grid_, file_name, extfrc_map_file, var_names, data_tracer); + auto file_reader = + scream::mam_coupling::create_tracer_data_reader(hor_rem, file_name); + VertEmissionsHorizInterp_.push_back(hor_rem); + VertEmissionsDataReader_.push_back(file_reader); + vert_emis_data_.push_back(data_tracer); + } // var_name vert emissions + int i = 0; + int offset_emis_ver = 0; + for(const auto &var_name : extfrc_lst_) { + const auto file_name = vert_emis_file_name_[var_name]; + const auto var_names = vert_emis_var_names_[var_name]; + const int nvars = static_cast(var_names.size()); + + forcings_[i].nsectors = nvars; + // I am assuming the order of species in extfrc_lst_. + // Indexing in mam4xx is fortran. + forcings_[i].frc_ndx = i + 1; + const auto io_grid_emis = VertEmissionsHorizInterp_[i]->get_src_grid(); + const int num_cols_io_emis = + io_grid_emis->get_num_local_dofs(); // Number of columns on this rank + const int num_levs_io_emis = + io_grid_emis + ->get_num_vertical_levels(); // Number of levels per column + vert_emis_data_[i].init(num_cols_io_emis, num_levs_io_emis, nvars); + vert_emis_data_[i].allocate_temporal_views(); + forcings_[i].file_alt_data = vert_emis_data_[i].has_altitude_; + for(int isp = 0; isp < nvars; ++isp) { + forcings_[i].offset = offset_emis_ver; + vert_emis_output_[isp + offset_emis_ver] = + view_2d("vert_emis_output_", ncol_, nlev_); + } + offset_emis_ver += nvars; + ++i; + } // end i + EKAT_REQUIRE_MSG( + offset_emis_ver <= int(mam_coupling::MAX_NUM_VERT_EMISSION_FIELDS), + "Error! Number of fields is bigger than " + "MAX_NUM_VERT_EMISSION_FIELDS. Increase the " + "MAX_NUM_VERT_EMISSION_FIELDS in tracer_reader_utils.hpp \n"); + + } // Tracer external forcing data +} // set_grids + +// ================================================================ +// 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 MAMMicrophysics::requested_buffer_size_in_bytes() const { return mam_coupling::buffer_size(ncol_, nlev_); } -void MAMMicrophysics::init_buffers(const ATMBufferManager &buffer_manager) { - EKAT_REQUIRE_MSG(buffer_manager.allocated_bytes() >= requested_buffer_size_in_bytes(), - "Error! Insufficient buffer size.\n"); +// ================================================================ +// INIT_BUFFERS +// ================================================================ +// ON HOST, initializeŃ• 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. - 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 MAMMicrophysics."); +void MAMMicrophysics::init_buffers(const ATMBufferManager &buffer_manager) { + 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 MAMMicrophysics." + " Used memory: " + << std::to_string(used_mem) + << "." + " Requested memory: " + << std::to_string(requested_buffer_size_in_bytes()) + << ". \n"); } - +// ================================================================ +// INITIALIZE_IMPL +// ================================================================ void MAMMicrophysics::initialize_impl(const RunType run_type) { - - step_ = 0; - - // populate the wet and dry atmosphere states with views from fields and - // the buffer - wet_atm_.qv = get_field_in("qv").get_view(); - wet_atm_.qc = get_field_out("qc").get_view(); - wet_atm_.nc = get_field_out("nc").get_view(); - wet_atm_.qi = get_field_in("qi").get_view(); - wet_atm_.ni = get_field_in("ni").get_view(); - - - 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_.cldfrac = get_field_in("cldfrac_tot").get_view(); // FIXME: tot or liq? - dry_atm_.pblh = get_field_in("pbl_height").get_view(); - dry_atm_.phis = get_field_in("phis").get_view(); - dry_atm_.omega = get_field_in("omega").get_view(); + // Determine orbital year. If orbital_year is negative, use current year + // from timestamp for orbital year; if positive, use provided orbital year + // for duration of simulation. + m_orbital_year = m_params.get("orbital_year", -9999); + + // Get orbital parameters from yaml file + m_orbital_eccen = m_params.get("orbital_eccentricity", -9999); + m_orbital_obliq = m_params.get("orbital_obliquity", -9999); + m_orbital_mvelp = m_params.get("orbital_mvelp", -9999); + + // --------------------------------------------------------------- + // 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(); + 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(); + + 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_int = get_field_in("p_int").get_view(); + dry_atm_.p_del = get_field_in("pseudo_density").get_view(); + dry_atm_.cldfrac = get_field_in("cldfrac_liq").get_view(); + dry_atm_.pblh = get_field_in("pbl_height").get_view(); + dry_atm_.phis = get_field_in("phis").get_view(); + dry_atm_.omega = get_field_in("omega").get_view(); dry_atm_.z_mid = buffer_.z_mid; dry_atm_.dz = buffer_.dz; dry_atm_.z_iface = buffer_.z_iface; @@ -219,43 +407,60 @@ void MAMMicrophysics::initialize_impl(const RunType run_type) { 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 - - // perform any initialization work - if (run_type==RunType::Initial) { - } - - // set wet/dry aerosol state data (interstitial aerosols only) - for (int m = 0; m < mam_coupling::num_aero_modes(); ++m) { - const char* 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(); + dry_atm_.z_surf = 0.0; // It is always zero. + + // get surface albedo: shortwave, direct + d_sfc_alb_dir_vis_ = get_field_in("sfc_alb_dir_vis").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 char *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(); dry_aero_.int_aero_nmr[m] = buffer_.dry_int_aero_nmr[m]; - for (int a = 0; a < mam_coupling::num_aero_species(); ++a) { - const char* int_mmr_field_name = mam_coupling::int_aero_mmr_field_name(m, a); - if (strlen(int_mmr_field_name) > 0) { - wet_aero_.int_aero_mmr[m][a] = get_field_out(int_mmr_field_name).get_view(); + + // cloudborne aerosol tracers of interest: number (n) mixing ratios + const char *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(); + dry_aero_.cld_aero_nmr[m] = buffer_.dry_cld_aero_nmr[m]; + + for(int a = 0; a < mam_coupling::num_aero_species(); ++a) { + // (interstitial) aerosol tracers of interest: mass (q) mixing ratios + const char *int_mmr_field_name = + mam_coupling::int_aero_mmr_field_name(m, a); + if(strlen(int_mmr_field_name) > 0) { + wet_aero_.int_aero_mmr[m][a] = + get_field_out(int_mmr_field_name).get_view(); dry_aero_.int_aero_mmr[m][a] = buffer_.dry_int_aero_mmr[m][a]; } - } - } + + // (cloudborne) aerosol tracers of interest: mass (q) mixing ratios + const char *cld_mmr_field_name = + mam_coupling::cld_aero_mmr_field_name(m, a); + if(strlen(cld_mmr_field_name) > 0) { + wet_aero_.cld_aero_mmr[m][a] = + get_field_out(cld_mmr_field_name).get_view(); + dry_aero_.cld_aero_mmr[m][a] = buffer_.dry_cld_aero_mmr[m][a]; + } + } // for loop species + } // for loop num_aero_modes() // set wet/dry aerosol-related gas state data - for (int g = 0; g < mam_coupling::num_aero_gases(); ++g) { - const char* mmr_field_name = mam_coupling::gas_mmr_field_name(g); - wet_aero_.gas_mmr[g] = get_field_out(mmr_field_name).get_view(); + for(int g = 0; g < mam_coupling::num_aero_gases(); ++g) { + const char *mmr_field_name = mam_coupling::gas_mmr_field_name(g); + wet_aero_.gas_mmr[g] = get_field_out(mmr_field_name).get_view(); dry_aero_.gas_mmr[g] = buffer_.dry_gas_mmr[g]; } // create our photolysis rate calculation table - photo_table_ = impl::read_photo_table(get_comm(), - config_.photolysis.rsf_file, - config_.photolysis.xs_long_file); - - // FIXME: read relevant land use data from drydep surface file + const std::string rsf_file = m_params.get("mam4_rsf_file"); + const std::string xs_long_file = + m_params.get("mam4_xs_long_file"); - // set up our preprocess/postprocess functors - preprocess_.initialize(ncol_, nlev_, wet_atm_, wet_aero_, dry_atm_, dry_aero_); - postprocess_.initialize(ncol_, nlev_, wet_atm_, wet_aero_, dry_atm_, dry_aero_); + photo_table_ = impl::read_photo_table(rsf_file, xs_long_file); // set field property checks for the fields in this process /* e.g. @@ -267,249 +472,342 @@ void MAMMicrophysics::initialize_impl(const RunType run_type) { add_postcondition_check(get_field_out("tke"),m_grid,0); */ - // set up WSM for internal local variables - // (we'll probably need this later, but we'll just use ATMBufferManager for now) - //const auto default_policy = ekat::ExeSpaceUtils::get_default_team_policy(ncol_, nlev_); - //workspace_mgr_.setup(buffer_.wsm_data, nlev_+1, 13+(n_wind_slots+n_trac_slots), default_policy); -} + { + // climatology data for linear stratospheric chemistry + auto linoz_o3_clim = buffer_.scratch[0]; // ozone (climatology) [vmr] + auto linoz_o3col_clim = + buffer_.scratch[1]; // column o3 above box (climatology) [Dobson Units + // (DU)] + auto linoz_t_clim = buffer_.scratch[2]; // temperature (climatology) [K] + auto linoz_PmL_clim = + buffer_.scratch[3]; // P minus L (climatology) [vmr/s] + auto linoz_dPmL_dO3 = + buffer_.scratch[4]; // sensitivity of P minus L to O3 [1/s] + auto linoz_dPmL_dT = + buffer_.scratch[5]; // sensitivity of P minus L to T3 [K] + auto linoz_dPmL_dO3col = buffer_.scratch[6]; // sensitivity of P minus L to + // overhead O3 column [vmr/DU] + auto linoz_cariolle_pscs = + buffer_.scratch[7]; // Cariolle parameter for PSC loss of ozone [1/s] + + auto ts = timestamp(); + std::string linoz_chlorine_file = + m_params.get("mam4_linoz_chlorine_file"); + int chlorine_loading_ymd = m_params.get("mam4_chlorine_loading_ymd"); + scream::mam_coupling::create_linoz_chlorine_reader( + linoz_chlorine_file, ts, chlorine_loading_ymd, chlorine_values_, + chlorine_time_secs_); + } // LINOZ + + const int photo_table_len = get_photo_table_work_len(photo_table_); + work_photo_table_ = view_2d("work_photo_table", ncol_, photo_table_len); -void MAMMicrophysics::run_impl(const double dt) { + // here's where we store per-column photolysis rates + photo_rates_ = view_3d("photo_rates", ncol_, nlev_, mam4::mo_photo::phtcnt); + + // Load the first month into extfrc_lst_end. + // Note: At the first time step, the data will be moved into extfrc_lst_beg, + // and extfrc_lst_end will be reloaded from file with the new month. + const int curr_month = timestamp().get_month() - 1; // 0-based + + scream::mam_coupling::update_tracer_data_from_file( + LinozDataReader_, curr_month, *LinozHorizInterp_, linoz_data_); + + scream::mam_coupling::update_tracer_data_from_file( + TracerDataReader_, curr_month, *TracerHorizInterp_, tracer_data_); + + for(int i = 0; i < static_cast(extfrc_lst_.size()); ++i) { + scream::mam_coupling::update_tracer_data_from_file( + VertEmissionsDataReader_[i], curr_month, *VertEmissionsHorizInterp_[i], + vert_emis_data_[i]); + } + + invariants_ = view_3d("invarians", ncol_, nlev_, mam4::gas_chemistry::nfs); + + constexpr int extcnt = mam4::gas_chemistry::extcnt; + extfrc_ = view_3d("extfrc_", ncol_, nlev_, extcnt); + + // + acos_cosine_zenith_host_ = view_1d_host("host_acos(cosine_zenith)", ncol_); + acos_cosine_zenith_ = view_1d("device_acos(cosine_zenith)", ncol_); + + //----------------------------------------------------------------- + // Setup preprocessing and post processing + //----------------------------------------------------------------- + preprocess_.initialize(ncol_, nlev_, wet_atm_, wet_aero_, dry_atm_, + dry_aero_); + postprocess_.initialize(ncol_, nlev_, wet_atm_, wet_aero_, dry_atm_, + dry_aero_); + +} // initialize_impl - const auto scan_policy = ekat::ExeSpaceUtils::get_thread_range_parallel_scan_team_policy(ncol_, nlev_); - const auto policy = ekat::ExeSpaceUtils::get_default_team_policy(ncol_, nlev_); +// ================================================================ +// RUN_IMPL +// ================================================================ +void MAMMicrophysics::run_impl(const double dt) { + const auto scan_policy = ekat::ExeSpaceUtils< + KT::ExeSpace>::get_thread_range_parallel_scan_team_policy(ncol_, nlev_); + const auto policy = + ekat::ExeSpaceUtils::get_default_team_policy(ncol_, nlev_); // preprocess input -- needs a scan for the calculation of atm height Kokkos::parallel_for("preprocess", scan_policy, preprocess_); Kokkos::fence(); - // reset internal WSM variables - //workspace_mgr_.reset_internals(); - - // NOTE: nothing depends on simulation time (yet), so we can just use zero for now - double t = 0.0; - - // here's where we store per-column photolysis rates - using View2D = haero::DeviceType::view_2d; - View2D photo_rates("photo_rates", nlev_, mam4::mo_photo::phtcnt); + const auto wet_geometric_mean_diameter_i = + get_field_in("dgnumwet").get_view(); + const auto dry_geometric_mean_diameter_i = + get_field_in("dgnum").get_view(); + const auto wetdens = get_field_in("wetdens").get_view(); // climatology data for linear stratospheric chemistry - auto linoz_o3_clim = buffer_.scratch[0]; // ozone (climatology) [vmr] - auto linoz_o3col_clim = buffer_.scratch[1]; // column o3 above box (climatology) [Dobson Units (DU)] - auto linoz_t_clim = buffer_.scratch[2]; // temperature (climatology) [K] - auto linoz_PmL_clim = buffer_.scratch[3]; // P minus L (climatology) [vmr/s] - auto linoz_dPmL_dO3 = buffer_.scratch[4]; // sensitivity of P minus L to O3 [1/s] - auto linoz_dPmL_dT = buffer_.scratch[5]; // sensitivity of P minus L to T3 [K] - auto linoz_dPmL_dO3col = buffer_.scratch[6]; // sensitivity of P minus L to overhead O3 column [vmr/DU] - auto linoz_cariolle_psc = buffer_.scratch[7]; // Cariolle parameter for PSC loss of ozone [1/s] - + // ozone (climatology) [vmr] + auto linoz_o3_clim = buffer_.scratch[0]; + // column o3 above box (climatology) [Dobson Units (DU)] + auto linoz_o3col_clim = buffer_.scratch[1]; + auto linoz_t_clim = buffer_.scratch[2]; // temperature (climatology) [K] + auto linoz_PmL_clim = buffer_.scratch[3]; // P minus L (climatology) [vmr/s] + // sensitivity of P minus L to O3 [1/s] + auto linoz_dPmL_dO3 = buffer_.scratch[4]; + // sensitivity of P minus L to T3 [K] + auto linoz_dPmL_dT = buffer_.scratch[5]; + // sensitivity of P minus L to overhead O3 column [vmr/DU] + auto linoz_dPmL_dO3col = buffer_.scratch[6]; + // Cariolle parameter for PSC loss of ozone [1/s] + auto linoz_cariolle_pscs = buffer_.scratch[7]; + + view_2d linoz_output[8]; + linoz_output[0] = linoz_o3_clim; + linoz_output[1] = linoz_o3col_clim; + linoz_output[2] = linoz_t_clim; + linoz_output[3] = linoz_PmL_clim; + linoz_output[4] = linoz_dPmL_dO3; + linoz_output[5] = linoz_dPmL_dT; + linoz_output[6] = linoz_dPmL_dO3col; + linoz_output[7] = linoz_cariolle_pscs; // it's a bit wasteful to store this for all columns, but simpler from an // allocation perspective auto o3_col_dens = buffer_.scratch[8]; - const_view_1d &col_latitudes = col_latitudes_; - mam_coupling::DryAtmosphere &dry_atm = dry_atm_; - mam_coupling::AerosolState &dry_aero = dry_aero_; - mam4::mo_photo::PhotoTableData &photo_table = photo_table_; - const int nlev = nlev_; - const Config &config = config_; - // FIXME: read relevant linoz climatology data from file(s) based on time - - // FIXME: read relevant chlorine loading data from file based on time - - // loop over atmosphere columns and compute aerosol microphyscs - auto some_step = step_; - - Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const ThreadTeam& team) { - const int icol = team.league_rank(); // column index - - Real col_lat = col_latitudes(icol); // column latitude (degrees?) - - // fetch column-specific atmosphere state data - auto atm = mam_coupling::atmosphere_for_column(dry_atm, icol); - auto z_iface = ekat::subview(dry_atm.z_iface, icol); - Real phis = dry_atm.phis(icol); - - // set surface state data - haero::Surface sfc{}; - - // fetch column-specific subviews into aerosol prognostics - mam4::Prognostics progs = mam_coupling::interstitial_aerosols_for_column(dry_aero, icol); - - // set up diagnostics - mam4::Diagnostics diags(nlev); - - // calculate o3 column densities (first component of col_dens in Fortran code) - auto o3_col_dens_i = ekat::subview(o3_col_dens, icol); - impl::compute_o3_column_density(team, atm, progs, o3_col_dens_i); - - // set up photolysis work arrays for this column. - mam4::mo_photo::PhotoTableWorkArrays photo_work_arrays; - // FIXME: set views here - - // ... look up photolysis rates from our table - // NOTE: the table interpolation operates on an entire column of data, so we - // NOTE: must do it before dispatching to individual vertical levels - Real zenith_angle = 0.0; // FIXME: need to get this from EAMxx [radians] - Real surf_albedo = 0.0; // FIXME: surface albedo - Real esfact = 0.0; // FIXME: earth-sun distance factor - mam4::ColumnView lwc; // FIXME: liquid water cloud content: where do we get this? - mam4::mo_photo::table_photo(photo_rates, atm.pressure, atm.hydrostatic_dp, - atm.temperature, o3_col_dens_i, zenith_angle, surf_albedo, lwc, - atm.cloud_fraction, esfact, photo_table, photo_work_arrays); - - // compute external forcings at time t(n+1) [molecules/cm^3/s] - constexpr int extcnt = mam4::gas_chemistry::extcnt; - view_2d extfrc; // FIXME: where to allocate? (nlev, extcnt) - mam4::mo_setext::Forcing forcings[extcnt]; // FIXME: forcings seem to require file data - mam4::mo_setext::extfrc_set(forcings, extfrc); - - // compute aerosol microphysics on each vertical level within this column - Kokkos::parallel_for(Kokkos::TeamThreadRange(team, nlev), [&](const int k) { - - constexpr int num_modes = mam4::AeroConfig::num_modes(); - constexpr int gas_pcnst = mam_coupling::gas_pcnst(); - constexpr int nqtendbb = mam_coupling::nqtendbb(); - - // extract atm state variables (input) - Real temp = atm.temperature(k); - Real pmid = atm.pressure(k); - Real pdel = atm.hydrostatic_dp(k); - Real zm = atm.height(k); - Real zi = z_iface(k); - Real pblh = atm.planetary_boundary_layer_height; - Real qv = atm.vapor_mixing_ratio(k); - Real cldfrac = atm.cloud_fraction(k); - - // extract aerosol state variables into "working arrays" (mass mixing ratios) - // (in EAM, this is done in the gas_phase_chemdr subroutine defined within - // mozart/mo_gas_phase_chemdr.F90) - Real q[gas_pcnst] = {}; - Real qqcw[gas_pcnst] = {}; - mam_coupling::transfer_prognostics_to_work_arrays(progs, k, q, qqcw); - - // convert mass mixing ratios to volume mixing ratios (VMR), equivalent - // to tracer mixing ratios (TMR)) - Real vmr[gas_pcnst], vmrcw[gas_pcnst]; - mam_coupling::convert_work_arrays_to_vmr(q, qqcw, vmr, vmrcw); - - // aerosol/gas species tendencies (output) - Real vmr_tendbb[gas_pcnst][nqtendbb] = {}; - Real vmrcw_tendbb[gas_pcnst][nqtendbb] = {}; - - // create work array copies to retain "pre-chemistry" values - Real vmr_pregaschem[gas_pcnst] = {}; - Real vmr_precldchem[gas_pcnst] = {}; - Real vmrcw_precldchem[gas_pcnst] = {}; - for (int i = 0; i < gas_pcnst; ++i) { - vmr_pregaschem[i] = vmr[i]; - vmr_precldchem[i] = vmr[i]; - vmrcw_precldchem[i] = vmrcw[i]; - } + /* Gather time and state information for interpolation */ + const auto ts = timestamp() + dt; + + const Real chlorine_loading = scream::mam_coupling::chlorine_loading_advance( + ts, chlorine_values_, chlorine_time_secs_); + + // /* Update the TracerTimeState to reflect the current time, note the + // addition of dt */ + trace_time_state_.t_now = ts.frac_of_year_in_days(); + scream::mam_coupling::advance_tracer_data( + TracerDataReader_, // in + *TracerHorizInterp_, // out + ts, // in + trace_time_state_, tracer_data_, // out + dry_atm_.p_mid, dry_atm_.z_iface, // in + cnst_offline_); // out + Kokkos::fence(); - //--------------------- - // Gas Phase Chemistry - //--------------------- - Real photo_rates_k[mam4::mo_photo::phtcnt]; - for (int i = 0; i < mam4::mo_photo::phtcnt; ++i) { - photo_rates_k[i] = photo_rates(k, i); - } - Real extfrc_k[extcnt]; - for (int i = 0; i < extcnt; ++i) { - extfrc_k[i] = extfrc(k, i); - } - constexpr int nfs = mam4::gas_chemistry::nfs; // number of "fixed species" - // NOTE: we compute invariants here and pass them out to use later with - // NOTE: setsox - Real invariants[nfs]; - impl::gas_phase_chemistry(zm, zi, phis, temp, pmid, pdel, dt, - photo_rates_k, extfrc_k, vmr, invariants); - - //---------------------- - // Aerosol microphysics - //---------------------- - // the logic below is taken from the aero_model_gasaerexch subroutine in - // eam/src/chemistry/modal_aero/aero_model.F90 - - // aqueous chemistry ... - const int loffset = 8; // offset of first tracer in work arrays - // (taken from mam4xx setsox validation test) - const Real mbar = haero::Constants::molec_weight_dry_air; - constexpr int indexm = 0; // FIXME: index of xhnm in invariants array (??) - Real cldnum = 0.0; // FIXME: droplet number concentration: where do we get this? - setsox_single_level(loffset, dt, pmid, pdel, temp, mbar, lwc(k), - cldfrac, cldnum, invariants[indexm], config.setsox, vmrcw, vmr); - - // calculate aerosol water content using water uptake treatment - // * dry and wet diameters [m] - // * wet densities [kg/m3] - // * aerosol water mass mixing ratio [kg/kg] - Real dgncur_a[num_modes] = {}; - Real dgncur_awet[num_modes] = {}; - Real wetdens[num_modes] = {}; - Real qaerwat[num_modes] = {}; - 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, 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, - wetdens, qaerwat); - - //----------------- - // LINOZ chemistry - //----------------- - - // the following things are diagnostics, which we're not - // including in the first rev - Real do3_linoz, do3_linoz_psc, ss_o3, o3col_du_diag, o3clim_linoz_diag, - zenith_angle_degrees; - - // FIXME: Need to get chlorine loading data from file - Real chlorine_loading = 0.0; - - Real rlats = col_lat * M_PI / 180.0; // convert column latitude to radians - int o3_ndx = 0; // index of "O3" in solsym array (in EAM) - mam4::lin_strat_chem::lin_strat_chem_solve_kk(o3_col_dens_i(k), temp, - zenith_angle, pmid, dt, rlats, - linoz_o3_clim(icol, k), linoz_t_clim(icol, k), linoz_o3col_clim(icol, k), - linoz_PmL_clim(icol, k), linoz_dPmL_dO3(icol, k), linoz_dPmL_dT(icol, k), - linoz_dPmL_dO3col(icol, k), linoz_cariolle_psc(icol, k), - chlorine_loading, config.linoz.psc_T, vmr[o3_ndx], - do3_linoz, do3_linoz_psc, ss_o3, - o3col_du_diag, o3clim_linoz_diag, zenith_angle_degrees); - - // update source terms above the ozone decay threshold - if (k > nlev - config.linoz.o3_lbl - 1) { - Real do3_mass; // diagnostic, not needed - mam4::lin_strat_chem::lin_strat_sfcsink_kk(dt, pdel, vmr[o3_ndx], config.linoz.o3_sfc, - config.linoz.o3_tau, do3_mass); - } + scream::mam_coupling::advance_tracer_data( + LinozDataReader_, // in + *LinozHorizInterp_, // out + ts, // in + linoz_time_state_, linoz_data_, // out + dry_atm_.p_mid, dry_atm_.z_iface, // in + linoz_output); // out + Kokkos::fence(); - // ... check for negative values and reset to zero - for (int i = 0; i < gas_pcnst; ++i) { - if (vmr[i] < 0.0) vmr[i] = 0.0; - } + vert_emiss_time_state_.t_now = ts.frac_of_year_in_days(); + int i = 0; + for(const auto &var_name : extfrc_lst_) { + const auto file_name = vert_emis_file_name_[var_name]; + const auto var_names = vert_emis_var_names_[var_name]; + const int nsectors = int(var_names.size()); + view_2d vert_emis_output[nsectors]; + for(int isp = 0; isp < nsectors; ++isp) { + vert_emis_output[isp] = vert_emis_output_[isp + forcings_[i].offset]; + } + scream::mam_coupling::advance_tracer_data( + VertEmissionsDataReader_[i], *VertEmissionsHorizInterp_[i], ts, + vert_emiss_time_state_, vert_emis_data_[i], dry_atm_.p_mid, + dry_atm_.z_iface, vert_emis_output); + i++; + Kokkos::fence(); + } - //---------------------- - // Dry deposition (gas) - //---------------------- + const_view_1d &col_latitudes = col_latitudes_; + const_view_1d &d_sfc_alb_dir_vis = d_sfc_alb_dir_vis_; - // FIXME: C++ port in progress! - //mam4::drydep::drydep_xactive(...); + mam_coupling::DryAtmosphere &dry_atm = dry_atm_; + mam_coupling::AerosolState &dry_aero = dry_aero_; - // transfer updated prognostics from work arrays - mam_coupling::convert_work_arrays_to_mmr(vmr, vmrcw, q, qqcw); - mam_coupling::transfer_work_arrays_to_prognostics(q, qqcw, progs, k); - }); - }); + mam4::mo_photo::PhotoTableData &photo_table = photo_table_; + const Config &config = config_; + const auto &work_photo_table = work_photo_table_; + const auto &photo_rates = photo_rates_; + + const auto &invariants = invariants_; + const auto &cnst_offline = cnst_offline_; + + // Compute orbital parameters; these are used both for computing + // the solar zenith angle. + // Note: We are following the RRTMGP EAMxx interface to compute the zenith + // angle. This operation is performed on the host because the routine + // shr_orb_cosz_c2f has not been ported to C++. + auto ts2 = timestamp(); + auto orbital_year = m_orbital_year; + // Note: We need double precision because + // shr_orb_params_c2f and shr_orb_decl_c2f only support double precision. + double obliqr, lambm0, mvelpp; + double eccen = m_orbital_eccen; + double obliq = m_orbital_obliq; + double mvelp = m_orbital_mvelp; + // Use the orbital parameters to calculate the solar declination and + // eccentricity factor + double delta, eccf; + if(eccen >= 0 && obliq >= 0 && mvelp >= 0) { + // use fixed orbital parameters; to force this, we need to set + // orbital_year to SHR_ORB_UNDEF_INT, which is exposed through + // our c2f bridge as shr_orb_undef_int_c2f + orbital_year = shr_orb_undef_int_c2f; + } else if(orbital_year < 0) { + // compute orbital parameters based on current year + orbital_year = ts2.get_year(); + } + shr_orb_params_c2f(&orbital_year, // in + &eccen, &obliq, &mvelp, &obliqr, &lambm0, &mvelpp); // out + + // Want day + fraction; calday 1 == Jan 1 0Z + auto calday = ts2.frac_of_year_in_days() + 1; + shr_orb_decl_c2f(calday, eccen, mvelpp, lambm0, obliqr, // in + &delta, &eccf); // out + { + const auto col_latitudes_host = + grid_->get_geometry_data("lat").get_view(); + const auto col_longitudes_host = + grid_->get_geometry_data("lon").get_view(); + // get a host copy of lat/lon + // Determine the cosine zenith angle + // NOTE: Since we are bridging to F90 arrays this must be done on HOST and + // then deep copied to a device view. + + // Now use solar declination to calculate zenith angle for all points + for(int i = 0; i < ncol_; i++) { + Real lat = + col_latitudes_host(i) * M_PI / 180.0; // Convert lat/lon to radians + Real lon = col_longitudes_host(i) * M_PI / 180.0; + // what's the aerosol microphys frequency? + Real temp = shr_orb_cosz_c2f(calday, lat, lon, delta, dt); + acos_cosine_zenith_host_(i) = acos(temp); + } + Kokkos::deep_copy(acos_cosine_zenith_, acos_cosine_zenith_host_); + } + const auto zenith_angle = acos_cosine_zenith_; + constexpr int gas_pcnst = mam_coupling::gas_pcnst(); + + const auto& vert_emis_output = vert_emis_output_; + const auto& extfrc = extfrc_; + const auto& forcings = forcings_; + constexpr int extcnt = mam4::gas_chemistry::extcnt; + + const int offset_aerosol = mam4::utils::gasses_start_ind(); + Real adv_mass_kg_per_moles[gas_pcnst]; + // NOTE: Making copies of clsmap_4 and permute_4 to fix undefined arrays on + // the device. + int clsmap_4[gas_pcnst], permute_4[gas_pcnst]; + for(int i = 0; i < gas_pcnst; ++i) { + // NOTE: state_q is kg/kg-dry-air; adv_mass is in g/mole. + // Convert adv_mass to kg/mole as vmr_from_mmr function uses + // molec_weight_dry_air with kg/mole units + adv_mass_kg_per_moles[i] = mam4::gas_chemistry::adv_mass[i] / 1e3; + clsmap_4[i] = mam4::gas_chemistry::clsmap_4[i]; + permute_4[i] = mam4::gas_chemistry::permute_4[i]; + } + // loop over atmosphere columns and compute aerosol microphyscs + Kokkos::parallel_for( + policy, KOKKOS_LAMBDA(const ThreadTeam &team) { + const int icol = team.league_rank(); // column index + const Real col_lat = col_latitudes(icol); // column latitude (degrees?) + + // convert column latitude to radians + const Real rlats = col_lat * M_PI / 180.0; + + // fetch column-specific atmosphere state data + const auto atm = mam_coupling::atmosphere_for_column(dry_atm, icol); + const auto wet_diameter_icol = + ekat::subview(wet_geometric_mean_diameter_i, icol); + const auto dry_diameter_icol = + ekat::subview(dry_geometric_mean_diameter_i, icol); + const auto wetdens_icol = ekat::subview(wetdens, icol); + + // fetch column-specific subviews into aerosol prognostics + mam4::Prognostics progs = + mam_coupling::aerosols_for_column(dry_aero, icol); + + const auto invariants_icol = ekat::subview(invariants, icol); + mam4::mo_setext::Forcing forcings_in[extcnt]; + + for(int i = 0; i < extcnt; ++i) { + const int nsectors = forcings[i].nsectors; + const int frc_ndx = forcings[i].frc_ndx; + const auto file_alt_data = forcings[i].file_alt_data; + + forcings_in[i].nsectors = nsectors; + forcings_in[i].frc_ndx = frc_ndx; + // We may need to move this line where we read files. + forcings_in[i].file_alt_data = file_alt_data; + for(int isec = 0; isec < forcings[i].nsectors; ++isec) { + const auto field = vert_emis_output[isec + forcings[i].offset]; + forcings_in[i].fields_data[isec] = ekat::subview(field, icol); + } + } // extcnt for loop + + const auto extfrc_icol = ekat::subview(extfrc, icol); + + view_1d cnst_offline_icol[mam4::mo_setinv::num_tracer_cnst]; + for(int i = 0; i < mam4::mo_setinv::num_tracer_cnst; ++i) { + cnst_offline_icol[i] = ekat::subview(cnst_offline[i], icol); + } + + // calculate o3 column densities (first component of col_dens in Fortran + // code) + auto o3_col_dens_i = ekat::subview(o3_col_dens, icol); + const auto &work_photo_table_icol = + ekat::subview(work_photo_table, icol); + + const auto &photo_rates_icol = ekat::subview(photo_rates, icol); + + const auto linoz_o3_clim_icol = ekat::subview(linoz_o3_clim, icol); + const auto linoz_t_clim_icol = ekat::subview(linoz_t_clim, icol); + const auto linoz_o3col_clim_icol = + ekat::subview(linoz_o3col_clim, icol); + const auto linoz_PmL_clim_icol = ekat::subview(linoz_PmL_clim, icol); + const auto linoz_dPmL_dO3_icol = ekat::subview(linoz_dPmL_dO3, icol); + const auto linoz_dPmL_dT_icol = ekat::subview(linoz_dPmL_dT, icol); + const auto linoz_dPmL_dO3col_icol = + ekat::subview(linoz_dPmL_dO3col, icol); + const auto linoz_cariolle_pscs_icol = + ekat::subview(linoz_cariolle_pscs, icol); + // Note: All variables are inputs, except for progs, which is an + // input/output variable. + mam4::microphysics::perform_atmospheric_chemistry_and_microphysics( + team, dt, rlats, cnst_offline_icol, forcings_in, atm, progs, + photo_table, chlorine_loading, config.setsox, config.amicphys, + config.linoz.psc_T, zenith_angle(icol), d_sfc_alb_dir_vis(icol), + o3_col_dens_i, photo_rates_icol, extfrc_icol, invariants_icol, + work_photo_table_icol, linoz_o3_clim_icol, linoz_t_clim_icol, + linoz_o3col_clim_icol, linoz_PmL_clim_icol, linoz_dPmL_dO3_icol, + linoz_dPmL_dT_icol, linoz_dPmL_dO3col_icol, + linoz_cariolle_pscs_icol, eccf, adv_mass_kg_per_moles, clsmap_4, + permute_4, offset_aerosol, + config.linoz.o3_sfc, config.linoz.o3_tau, config.linoz.o3_lbl, + dry_diameter_icol, wet_diameter_icol, wetdens_icol); + }); // parallel_for for the column loop + Kokkos::fence(); // postprocess output Kokkos::parallel_for("postprocess", policy, postprocess_); Kokkos::fence(); -} -void MAMMicrophysics::finalize_impl() { -} +} // MAMMicrophysics::run_impl -} // namespace scream +} // namespace scream diff --git a/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.hpp b/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.hpp index 5f5a44b846b..6b1dd33dfaa 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.hpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_microphysics_process_interface.hpp @@ -4,25 +4,12 @@ #include #include #include - -#include "impl/mam4_amicphys.cpp" // mam4xx top-level microphysics function(s) - -#include -#include +#include "readfiles/tracer_reader_utils.hpp" +// For calling MAM4 processes #include - #include -#ifndef KOKKOS_ENABLE_CUDA -#define protected_except_cuda public -#define private_except_cuda public -#else -#define protected_except_cuda protected -#define private_except_cuda private -#endif - -namespace scream -{ +namespace scream { // The process responsible for handling MAM4 aerosol microphysics. The AD // stores exactly ONE instance of this class in its list of subcomponents. @@ -31,28 +18,19 @@ class MAMMicrophysics final : public scream::AtmosphereProcess { using KT = ekat::KokkosTypes; // views for single- and multi-column data - using view_1d_int = typename KT::template view_1d; using view_1d = typename KT::template view_1d; using view_2d = typename KT::template view_2d; + using view_3d = typename KT::template view_3d; using const_view_1d = typename KT::template view_1d; - using const_view_2d = typename KT::template view_2d; - - // unmanaged views (for buffer and workspace manager) - using uview_1d = Unmanaged>; - using uview_2d = Unmanaged>; - // a quantity stored in a single vertical column with a single index - using ColumnView = mam4::ColumnView; + using view_1d_host = typename KT::view_1d::HostMirror; // a thread team dispatched to a single vertical column using ThreadTeam = mam4::ThreadTeam; -public: - + public: // Constructor - MAMMicrophysics(const ekat::Comm& comm, const ekat::ParameterList& params); - -protected_except_cuda: + MAMMicrophysics(const ekat::Comm &comm, const ekat::ParameterList ¶ms); // -------------------------------------------------------------------------- // AtmosphereProcess overrides (see share/atm_process/atmosphere_process.hpp) @@ -60,59 +38,59 @@ class MAMMicrophysics final : public scream::AtmosphereProcess { // process metadata AtmosphereProcessType type() const override; - std::string name() const override; - // set aerosol microphysics configuration parameters (called by constructor) - void configure(const ekat::ParameterList& params); + // The name of the subcomponent + std::string name() const { return "mam_aero_microphysics"; } // grid - void set_grids(const std::shared_ptr grids_manager) override; + 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; - // process behavior + // Initialize variables void initialize_impl(const RunType run_type) override; - void run_impl(const double dt) override; - void finalize_impl() override; - // performs some checks on the tracers group - void set_computed_group_impl(const FieldGroup& group) override; + // Run the process by one time step + void run_impl(const double dt) override; -private_except_cuda: + // Finalize + void finalize_impl(){/*Do nothing*/}; + private: // number of horizontal columns and vertical levels int ncol_, nlev_; - // configuration data (for the moment, we plan to be able to move this to - // the device, so we can't use C++ strings) + // The orbital year, used for zenith angle calculations: + // If > 0, use constant orbital year for duration of simulation + // If < 0, use year from timestamp for orbital parameters + int m_orbital_year; + + // Orbital parameters, used for zenith angle calculations. + // If >= 0, bypass computation based on orbital year and use fixed parameters + // If < 0, compute based on orbital year, specified above + // These variables are required to be double. + double m_orbital_eccen; // Eccentricity + double m_orbital_obliq; // Obliquity + double m_orbital_mvelp; // Vernal Equinox Mean Longitude of Perihelion + struct Config { - // photolysis parameters - struct { - char rsf_file[MAX_FILENAME_LEN]; - char xs_long_file[MAX_FILENAME_LEN]; - } photolysis; // stratospheric chemistry parameters struct { - int o3_lbl; // number of layers with ozone decay from the surface - int o3_sfc; // set from namelist input linoz_sfc - int o3_tau; // set from namelist input linoz_tau - Real psc_T; // set from namelist input linoz_psc_T - char chlorine_loading_file[MAX_FILENAME_LEN]; + int o3_lbl; // number of layers with ozone decay from the surface + Real o3_sfc; // set from namelist input linoz_sfc + Real o3_tau; // set from namelist input linoz_tau + Real psc_T; // set from namelist input linoz_psc_T } linoz; // aqueous chemistry parameters mam4::mo_setsox::Config setsox; // aero microphysics configuration (see impl/mam4_amicphys.cpp) - impl::AmicPhysConfig amicphys; - - // dry deposition parameters - struct { - char srf_file[MAX_FILENAME_LEN]; - } drydep; + mam4::microphysics::AmicPhysConfig amicphys; }; Config config_; @@ -124,39 +102,41 @@ class MAMMicrophysics final : public scream::AtmosphereProcess { // on host: initializes preprocess functor with necessary state data void initialize(const int ncol, const int nlev, - const mam_coupling::WetAtmosphere& wet_atm, - const mam_coupling::AerosolState& wet_aero, - const mam_coupling::DryAtmosphere& dry_atm, - const mam_coupling::AerosolState& dry_aero) { - ncol_ = ncol; - nlev_ = nlev; - wet_atm_ = wet_atm; - wet_aero_ = wet_aero; - dry_atm_ = dry_atm; - dry_aero_ = dry_aero; + const mam_coupling::WetAtmosphere &wet_atm, + const mam_coupling::AerosolState &wet_aero, + const mam_coupling::DryAtmosphere &dry_atm, + const mam_coupling::AerosolState &dry_aero) { + ncol_pre_ = ncol; + nlev_pre_ = nlev; + wet_atm_pre_ = wet_atm; + wet_aero_pre_ = wet_aero; + dry_atm_pre_ = dry_atm; + dry_aero_pre_ = dry_aero; } KOKKOS_INLINE_FUNCTION - void operator()(const Kokkos::TeamPolicy::member_type& team) const { - const int i = team.league_rank(); // column index - - compute_vertical_layer_heights(team, dry_atm_, i); - team.team_barrier(); // allows kernels below to use layer heights - compute_updraft_velocities(team, wet_atm_, dry_atm_, i); - compute_dry_mixing_ratios(team, wet_atm_, dry_atm_, i); - compute_dry_mixing_ratios(team, wet_atm_, wet_aero_, dry_aero_, i); + void operator()( + const Kokkos::TeamPolicy::member_type &team) const { + const int i = team.league_rank(); // column index + + compute_dry_mixing_ratios(team, wet_atm_pre_, dry_atm_pre_, i); + compute_dry_mixing_ratios(team, wet_atm_pre_, wet_aero_pre_, + dry_aero_pre_, i); team.team_barrier(); - } // operator() + + compute_vertical_layer_heights(team, dry_atm_pre_, i); + compute_updraft_velocities(team, wet_atm_pre_, dry_atm_pre_, i); + } // operator() // number of horizontal columns and vertical levels - int ncol_, nlev_; + int ncol_pre_, nlev_pre_; // local atmospheric and aerosol state data - mam_coupling::WetAtmosphere wet_atm_; - mam_coupling::DryAtmosphere dry_atm_; - mam_coupling::AerosolState wet_aero_, dry_aero_; + mam_coupling::WetAtmosphere wet_atm_pre_; + mam_coupling::DryAtmosphere dry_atm_pre_; + mam_coupling::AerosolState wet_aero_pre_, dry_aero_pre_; - }; // MAMMicrophysics::Preprocess + }; // MAMMicrophysics::Preprocess // Postprocessing functor struct Postprocess { @@ -164,33 +144,35 @@ class MAMMicrophysics final : public scream::AtmosphereProcess { // on host: initializes postprocess functor with necessary state data void initialize(const int ncol, const int nlev, - const mam_coupling::WetAtmosphere& wet_atm, - const mam_coupling::AerosolState& wet_aero, - const mam_coupling::DryAtmosphere& dry_atm, - const mam_coupling::AerosolState& dry_aero) { - ncol_ = ncol; - nlev_ = nlev; - wet_atm_ = wet_atm; - wet_aero_ = wet_aero; - dry_atm_ = dry_atm; - dry_aero_ = dry_aero; + const mam_coupling::WetAtmosphere &wet_atm, + const mam_coupling::AerosolState &wet_aero, + const mam_coupling::DryAtmosphere &dry_atm, + const mam_coupling::AerosolState &dry_aero) { + ncol_post_ = ncol; + nlev_post_ = nlev; + wet_atm_post_ = wet_atm; + wet_aero_post_ = wet_aero; + dry_atm_post_ = dry_atm; + dry_aero_post_ = dry_aero; } KOKKOS_INLINE_FUNCTION - void operator()(const Kokkos::TeamPolicy::member_type& team) const { - const int i = team.league_rank(); // column index - compute_wet_mixing_ratios(team, dry_atm_, dry_aero_, wet_aero_, i); + void operator()( + const Kokkos::TeamPolicy::member_type &team) const { + const int i = team.league_rank(); // column index + compute_wet_mixing_ratios(team, dry_atm_post_, dry_aero_post_, + wet_aero_post_, i); team.team_barrier(); - } // operator() + } // operator() // number of horizontal columns and vertical levels - int ncol_, nlev_; + int ncol_post_, nlev_post_; // local atmospheric and aerosol state data - mam_coupling::WetAtmosphere wet_atm_; - mam_coupling::DryAtmosphere dry_atm_; - mam_coupling::AerosolState wet_aero_, dry_aero_; - }; // MAMMicrophysics::Postprocess + mam_coupling::WetAtmosphere wet_atm_post_; + mam_coupling::DryAtmosphere dry_atm_post_; + mam_coupling::AerosolState wet_aero_post_, dry_aero_post_; + }; // MAMMicrophysics::Postprocess // MAM4 aerosol particle size description mam4::AeroConfig aero_config_; @@ -202,29 +184,63 @@ class MAMMicrophysics final : public scream::AtmosphereProcess { // atmospheric and aerosol state variables mam_coupling::WetAtmosphere wet_atm_; mam_coupling::DryAtmosphere dry_atm_; - mam_coupling::AerosolState wet_aero_, dry_aero_; + mam_coupling::AerosolState wet_aero_, dry_aero_; // photolysis rate table (column-independent) mam4::mo_photo::PhotoTableData photo_table_; // column areas, latitudes, longitudes - const_view_1d col_areas_, col_latitudes_, col_longitudes_; + const_view_1d col_latitudes_; - // time step number - int step_; + // surface albedo: shortwave, direct + const_view_1d d_sfc_alb_dir_vis_; // workspace manager for internal local variables - //ekat::WorkspaceManager workspace_mgr_; + // ekat::WorkspaceManager workspace_mgr_; mam_coupling::Buffer buffer_; // physics grid for column information std::shared_ptr grid_; - // sets defaults for "namelist parameters" - void set_defaults_(); - -}; // MAMMicrophysics - -} // namespace scream - -#endif // EAMXX_MAM_MICROPHYSICS_HPP + mam_coupling::TracerTimeState linoz_time_state_; + view_2d work_photo_table_; + std::vector chlorine_values_; + std::vector chlorine_time_secs_; + view_3d photo_rates_; + + // invariants members + mam_coupling::TracerTimeState trace_time_state_; + std::shared_ptr TracerDataReader_; + std::shared_ptr TracerHorizInterp_; + mam_coupling::TracerData tracer_data_; + view_3d invariants_; + std::string oxid_file_name_; + view_2d cnst_offline_[4]; + + // linoz reader + std::shared_ptr LinozDataReader_; + std::shared_ptr LinozHorizInterp_; + mam_coupling::TracerData linoz_data_; + std::string linoz_file_name_; + + // Vertical emission uses 9 files, here I am using std::vector to stote + // instance of each file. + mam_coupling::TracerTimeState vert_emiss_time_state_; + std::vector> VertEmissionsDataReader_; + std::vector> VertEmissionsHorizInterp_; + std::vector extfrc_lst_; + std::vector vert_emis_data_; + std::map vert_emis_file_name_; + std::map> vert_emis_var_names_; + view_2d vert_emis_output_[mam_coupling::MAX_NUM_VERT_EMISSION_FIELDS]; + view_3d extfrc_; + mam_coupling::ForcingHelper forcings_[mam4::gas_chemistry::extcnt]; + + view_1d_host acos_cosine_zenith_host_; + view_1d acos_cosine_zenith_; + +}; // MAMMicrophysics + +} // namespace scream + +#endif // EAMXX_MAM_MICROPHYSICS_HPP diff --git a/components/eamxx/src/physics/mam/eamxx_mam_wetscav_process_interface.cpp b/components/eamxx/src/physics/mam/eamxx_mam_wetscav_process_interface.cpp index 32148d580fe..863da8021dd 100644 --- a/components/eamxx/src/physics/mam/eamxx_mam_wetscav_process_interface.cpp +++ b/components/eamxx/src/physics/mam/eamxx_mam_wetscav_process_interface.cpp @@ -28,8 +28,8 @@ void MAMWetscav::set_grids( // The units of mixing ratio Q are technically non-dimensional. // Nevertheless, for output reasons, we like to see 'kg/kg'. - auto q_unit = kg / kg; - auto n_unit = 1 / kg; // units of number mixing ratios of tracers + auto q_unit = kg / kg; + auto n_unit = 1 / kg; // units of number mixing ratios of tracers m_grid = grids_manager->get_grid("Physics"); const auto &grid_name = m_grid->name(); @@ -205,7 +205,7 @@ void MAMWetscav::set_grids( static constexpr auto m3 = m * m * m; // Aerosol dry particle diameter [m] - add_field("dgncur_a", scalar3d_mid_nmodes, m, grid_name); + add_field("dgnum", scalar3d_mid_nmodes, m, grid_name); // Wet aerosol density [kg/m3] add_field("wetdens", scalar3d_mid_nmodes, kg / m3, grid_name); @@ -475,7 +475,7 @@ void MAMWetscav::run_impl(const double dt) { const auto wet_geometric_mean_diameter_i = get_field_out("dgnumwet").get_view(); const auto dry_geometric_mean_diameter_i = - get_field_out("dgncur_a").get_view(); + get_field_out("dgnum").get_view(); const auto qaerwat = get_field_out("qaerwat").get_view(); const auto wetdens = get_field_out("wetdens").get_view(); diff --git a/components/eamxx/src/physics/mam/impl/README.md b/components/eamxx/src/physics/mam/impl/README.md deleted file mode 100644 index f05a484a82a..00000000000 --- a/components/eamxx/src/physics/mam/impl/README.md +++ /dev/null @@ -1,11 +0,0 @@ -# MAM4 Integration Code - -This folder contains C++ implementations of the higher-level MAM4 interface -routines for aerosol microphysics, cloud-aerosol interactions, and optical -properties. We've retained the overall structure of the original Fortran code -to make it easier to understand for folks who are more familiar with the -original implementation of MAM4. - -## Contents - -* `mam4_amicphys.cpp` - high-level MAM4 microphysics interface code diff --git a/components/eamxx/src/physics/mam/impl/compute_o3_column_density.cpp b/components/eamxx/src/physics/mam/impl/compute_o3_column_density.cpp deleted file mode 100644 index ceb992bc886..00000000000 --- a/components/eamxx/src/physics/mam/impl/compute_o3_column_density.cpp +++ /dev/null @@ -1,43 +0,0 @@ -namespace scream::impl { - -KOKKOS_INLINE_FUNCTION -void compute_o3_column_density(const ThreadTeam& team, const haero::Atmosphere& atm, - const mam4::Prognostics &progs, ColumnView o3_col_dens) { - constexpr int gas_pcnst = mam4::gas_chemistry::gas_pcnst; // number of gas phase species - constexpr int nfs = mam4::gas_chemistry::nfs; // number of "fixed species" - // constexpr Real mwdry = 1.0/haero::Constants::molec_weight_dry_air; - - Real o3_col_deltas[mam4::nlev+1] = {}; // o3 column density above model [1/cm^2] - // NOTE: if we need o2 column densities, set_ub_col and setcol must be changed - Kokkos::parallel_for(Kokkos::TeamThreadRange(team, atm.num_levels()), [&](const int k) { - - // Real temp = atm.temperature(k); - // Real pmid = atm.pressure(k); - Real pdel = atm.hydrostatic_dp(k); - // Real qv = atm.vapor_mixing_ratio(k); - - // ... map incoming mass mixing ratios to working array - Real q[gas_pcnst], qqcw[gas_pcnst]; - mam_coupling::transfer_prognostics_to_work_arrays(progs, k, q, qqcw); - - // ... set atmosphere mean mass to the molecular weight of dry air - // and compute water vapor vmr - // Real mbar = mwdry; - // Real h2ovmr = mam4::conversions::vmr_from_mmr(qv, mbar); - - // ... Xform from mmr to vmr - Real vmr[gas_pcnst], vmrcw[gas_pcnst]; - mam_coupling::convert_work_arrays_to_vmr(q, qqcw, vmr, vmrcw); - - // ... compute invariants for this level - Real invariants[nfs]; - // setinv(invariants, temp, h2ovmr, vmr, pmid); FIXME: not yet ported - - // compute the change in o3 density for this column above its neighbor - mam4::mo_photo::set_ub_col(o3_col_deltas[k+1], vmr, invariants, pdel); - }); - // sum the o3 column deltas to densities - mam4::mo_photo::setcol(o3_col_deltas, o3_col_dens); -} - -} // namespace scream::impl diff --git a/components/eamxx/src/physics/mam/impl/compute_water_content.cpp b/components/eamxx/src/physics/mam/impl/compute_water_content.cpp deleted file mode 100644 index ea7190afff9..00000000000 --- a/components/eamxx/src/physics/mam/impl/compute_water_content.cpp +++ /dev/null @@ -1,99 +0,0 @@ -#include - -namespace scream::impl { - -KOKKOS_INLINE_FUNCTION -void compute_water_content(const mam4::Prognostics &progs, int k, - Real qv, Real temp, Real pmid, - Real dgncur_a[mam4::AeroConfig::num_modes()], - Real dgncur_awet[mam4::AeroConfig::num_modes()], - Real wetdens[mam4::AeroConfig::num_modes()], - Real qaerwat[mam4::AeroConfig::num_modes()]) { - constexpr int num_modes = mam4::AeroConfig::num_modes(); - constexpr int num_aero_ids = mam4::AeroConfig::num_aerosol_ids(); - - // get some information about aerosol species - // FIXME: this isn't great! - constexpr int maxd_aspectype = mam4::water_uptake::maxd_aspectype; - int nspec_amode[num_modes], lspectype_amode[maxd_aspectype][num_modes]; - Real specdens_amode[maxd_aspectype], spechygro[maxd_aspectype]; - mam4::water_uptake::get_e3sm_parameters(nspec_amode, lspectype_amode, - specdens_amode, spechygro); - - // extract aerosol tracers for this level into state_q, which is needed - // for computing dry aerosol properties below - // FIXME: we should eliminate this index translation stuff - constexpr int nvars = aero_model::pcnst; - Real state_q[nvars]; // aerosol tracers for level k - for (int imode = 0; imode < num_modes; ++imode) { - int la, lc; // interstitial and cloudborne indices within state_q - - // number mixing ratios - mam4::convproc::assign_la_lc(imode, -1, la, lc); - state_q[la] = progs.n_mode_i[imode](k); - state_q[lc] = progs.n_mode_c[imode](k); - // aerosol mass mixing ratios - for (int iaero = 0; iaero < num_aero_ids; ++iaero) { - mam4::convproc::assign_la_lc(imode, iaero, la, lc); - auto mode = static_cast(imode); - auto aero = static_cast(iaero); - int ispec = mam4::aerosol_index_for_mode(mode, aero); - if (ispec != -1) { - state_q[la] = progs.q_aero_i[imode][ispec](k); - state_q[lc] = progs.q_aero_c[imode][ispec](k); - } - } - } - - // compute the dry volume for each mode, and from it the current dry - // geometric nominal particle diameter. - // FIXME: We have to do some gymnastics here to set up the calls to - // FIXME: calcsize. This could be improved. - Real inv_densities[num_modes][num_aero_ids] = {}; - for (int imode = 0; imode < num_modes; ++imode) { - const int n_spec = mam4::num_species_mode(imode); - for (int ispec = 0; ispec < n_spec; ++ispec) { - const int iaer = static_cast(mam4::mode_aero_species(imode, ispec)); - const Real density = mam4::aero_species(iaer).density; - inv_densities[imode][ispec] = 1.0 / density; - } - } - for (int imode = 0; imode < num_modes; ++imode) { - Real dryvol_i, dryvol_c; // interstitial and cloudborne dry volumes - mam4::calcsize::compute_dry_volume_k(k, imode, inv_densities, progs, - dryvol_i, dryvol_c); - - // NOTE: there's some disagreement over whether vol2num should be called - // NOTE: num2vol here, so I'm just adopting the nomenclature used by - // NOTE: the following call to calcsize) - const mam4::Mode& mode = mam4::modes(imode); - Real vol2num_min = 1.0/mam4::conversions::mean_particle_volume_from_diameter( - mode.max_diameter, mode.mean_std_dev); - Real vol2num_max = 1.0/mam4::conversions::mean_particle_volume_from_diameter( - mode.min_diameter, mode.mean_std_dev); - Real vol2num; - mam4::calcsize::update_diameter_and_vol2num(dryvol_i, - progs.n_mode_i[imode](k), vol2num_min, vol2num_max, - mode.min_diameter, mode.max_diameter, mode.mean_std_dev, - dgncur_a[imode], vol2num); - } - - // calculate dry aerosol properties - Real hygro[num_modes], naer[num_modes], dryrad[num_modes], - dryvol[num_modes], drymass[num_modes], - rhcrystal[num_modes], rhdeliques[num_modes], specdens_1[num_modes]; - mam4::water_uptake::modal_aero_water_uptake_dryaer(nspec_amode, specdens_amode, - spechygro, lspectype_amode, state_q, dgncur_a, hygro, - naer, dryrad, dryvol, drymass, rhcrystal, rhdeliques, specdens_1); - - // calculate wet aerosol properties - Real rh = mam4::conversions::relative_humidity_from_vapor_mixing_ratio(qv, temp, pmid); - Real wetrad[num_modes], wetvol[num_modes], wtrvol[num_modes]; - mam4::water_uptake::modal_aero_water_uptake_wetaer(rhcrystal, rhdeliques, dgncur_a, - dryrad, hygro, rh, naer, dryvol, wetrad, wetvol, wtrvol, dgncur_awet, - qaerwat); - mam4::water_uptake::modal_aero_water_uptake_wetdens(wetvol, wtrvol, - drymass, specdens_1, wetdens); -} - -} // namespace scream::impl diff --git a/components/eamxx/src/physics/mam/impl/gas_phase_chemistry.cpp b/components/eamxx/src/physics/mam/impl/gas_phase_chemistry.cpp deleted file mode 100644 index 5401fa5e3da..00000000000 --- a/components/eamxx/src/physics/mam/impl/gas_phase_chemistry.cpp +++ /dev/null @@ -1,381 +0,0 @@ -#include - -namespace scream::impl { - -using mam4::utils::min_max_bound; - -using HostView1D = haero::DeviceType::view_1d::HostMirror; -using HostViewInt1D = haero::DeviceType::view_1d::HostMirror; - -//------------------------------------------------------------------------- -// Reading the photolysis table -//------------------------------------------------------------------------- -// This logic is currently implemented using serial NetCDF calls for -// clarity of purpose. We should probably read the data for the photolysis -// table using SCREAM's SCORPIO interface instead, but I wanted to make -// clear what we're trying to do in terms of "elementary" operations first. - -// ON HOST (MPI root rank only), reads the dimension of a NetCDF variable from -// the file with the given ID -int nc_dimension(const char *file, int nc_id, const char *dim_name) { - int dim_id; - int result = nc_inq_dimid(nc_id, dim_name, &dim_id); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't fetch " << dim_name << - " dimension ID from NetCDF file '" << file << "'\n"); - size_t dim; - result = nc_inq_dimlen(nc_id, dim_id, &dim); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't fetch " << dim_name << - " dimension from NetCDF file '" << file << "'\n"); - return static_cast(dim); -} - -// ON HOST (MPI root rank only), reads data from the given NetCDF variable from -// the file with the given ID into the given Kokkos host View -template -void read_nc_var(const char *file, int nc_id, const char *var_name, V host_view) { - int var_id; - int result = nc_inq_varid(nc_id, var_name, &var_id); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't fetch ID for variable '" << var_name << - "' from NetCDF file '" << file << "'\n"); - result = nc_get_var(nc_id, var_id, host_view.data()); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't read data for variable '" << var_name << - "' from NetCDF file '" << file << "'\n"); -} - -// ON HOST (MPI root rank only), reads data from the NetCDF variable with the -// given ID, from the file with the given ID, into the given Kokkos host View -template -void read_nc_var(const char *file, int nc_id, int var_id, V host_view) { - int result = nc_get_var(nc_id, var_id, host_view.data()); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't read data for variable with ID " << - var_id << " from NetCDF file '" << file << "'\n"); -} - -// ON HOST (MPI root only), sets the lng_indexer and pht_alias_mult_1 host views -// according to parameters in our (hardwired) chemical mechanism -void set_lng_indexer_and_pht_alias_mult_1(const char *file, int nc_id, - HostViewInt1D lng_indexer, - HostView1D pht_alias_mult_1) { - // NOTE: it seems that the chemical mechanism we're using - // NOTE: 1. sets pht_alias_lst to a blank string [1] - // NOTE: 2. sets pht_alias_mult_1 to 1.0 [1] - // NOTE: 3. sets rxt_tag_lst to ['jh2o2', 'usr_HO2_HO2', 'usr_SO2_OH', 'usr_DMS_OH'] [2] - // NOTE: References: - // NOTE: [1] (https://github.com/eagles-project/e3sm_mam4_refactor/blob/refactor-maint-2.0/components/eam/src/chemistry/pp_linoz_mam4_resus_mom_soag/mo_sim_dat.F90#L117) - // NOTE: [2] (https://github.com/eagles-project/e3sm_mam4_refactor/blob/refactor-maint-2.0/components/eam/src/chemistry/pp_linoz_mam4_resus_mom_soag/mo_sim_dat.F90#L99) - - // populate lng_indexer (see https://github.com/eagles-project/e3sm_mam4_refactor/blob/refactor-maint-2.0/components/eam/src/chemistry/mozart/mo_jlong.F90#L180) - static const char *var_names[4] = {"jh2o2", "usr_HO2_HO2", "usr_SO2_OH", "usr_DMS_OH"}; - for (int m = 0; m < mam4::mo_photo::phtcnt; ++m) { - int var_id; - int result = nc_inq_varid(nc_id, var_names[m], &var_id); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't fetch ID for variable '" - << var_names[m] << "' from NetCDF file '" << file << "'\n"); - lng_indexer(m) = var_id; - } - - // set pht_alias_mult_1 to 1 - Kokkos::deep_copy(pht_alias_mult_1, 1.0); -} - -// ON HOST (MPI root only), populates the etfphot view using rebinned -// solar data from our solar_data_file -void populate_etfphot(HostView1D we, HostView1D etfphot) { - // FIXME: It looks like EAM is relying on a piece of infrastructure that - // FIXME: we just don't have in EAMxx (eam/src/chemistry/utils/solar_data.F90). - // FIXME: I have no idea whether EAMxx has a plan for supporting this - // FIXME: solar irradiance / photon flux data, and I'm not going to recreate - // FIXME: that capability here. So this is an unplugged hole. - // FIXME: - // FIXME: If we are going to do this the way EAM does it, the relevant logic - // FIXME: is the call to rebin() in eam/src/chemistry/mozart/mo_jlong.F90, - // FIXME: around line 104. - - // FIXME: zero the photon flux for now - Kokkos::deep_copy(etfphot, 0); -} - -// ON HOST, reads the photolysis table (used for gas phase chemistry) from the -// files with the given names -mam4::mo_photo::PhotoTableData read_photo_table(const ekat::Comm& comm, - const char *rsf_file, - const char* xs_long_file) { - // NOTE: at the time of development, SCREAM's SCORPIO interface seems intended - // NOTE: for domain-decomposed grid data. The files we're reading here are not - // NOTE: spatial data, and should be the same everywhere, so we read them - // NOTE: using serial NetCDF calls on MPI rank 0 and broadcast to other ranks. - const int mpi_root = 0; - int rsf_id, xs_long_id; // NetCDF file IDs (used only on MPI root) - int nw, nump, numsza, numcolo3, numalb, nt, np_xs; // table dimensions - if (comm.rank() == mpi_root) { // read dimension data from files and broadcast - // open files - int result = nc_open(rsf_file, NC_NOWRITE, &rsf_id); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't open rsf_file '" << rsf_file << "'\n"); - result = nc_open(xs_long_file, NC_NOWRITE, &xs_long_id); - EKAT_REQUIRE_MSG(result == 0, "Error! Couldn't open xs_long_file '" << xs_long_file << "'\n"); - - // read and broadcast dimension data - nump = nc_dimension(rsf_file, rsf_id, "numz"); - numsza = nc_dimension(rsf_file, rsf_id, "numsza"); - numalb = nc_dimension(rsf_file, rsf_id, "numalb"); - numcolo3 = nc_dimension(rsf_file, rsf_id, "numcolo3fact"); - nt = nc_dimension(xs_long_file, xs_long_id, "numtemp"); - nw = nc_dimension(xs_long_file, xs_long_id, "numwl"); - np_xs = nc_dimension(xs_long_file, xs_long_id, "numprs"); - - int dim_data[7] = {nump, numsza, numcolo3, numalb, nt, nw, np_xs}; - comm.broadcast(dim_data, 7, mpi_root); - } else { // receive broadcasted dimension data from root rank - int dim_data[7]; - comm.broadcast(dim_data, 7, mpi_root); - nump = dim_data[0]; - numsza = dim_data[1]; - numcolo3 = dim_data[2]; - numalb = dim_data[3]; - nt = dim_data[4]; - nw = dim_data[5]; - np_xs = dim_data[6]; - } - - // set up the lng_indexer and pht_alias_mult_1 views based on our - // (hardwired) chemical mechanism - HostViewInt1D lng_indexer_h("lng_indexer(host)", mam4::mo_photo::phtcnt); - HostView1D pht_alias_mult_1_h("pht_alias_mult_1(host)", 2); - if (comm.rank() == mpi_root) { - set_lng_indexer_and_pht_alias_mult_1(xs_long_file, xs_long_id, - lng_indexer_h, pht_alias_mult_1_h); - } - comm.broadcast(lng_indexer_h.data(), mam4::mo_photo::phtcnt, mpi_root); - comm.broadcast(pht_alias_mult_1_h.data(), 2, mpi_root); - - // compute the size of the foremost dimension of xsqy using lng_indexer - int numj = 0; - for (int m = 0; m < mam4::mo_photo::phtcnt; ++m) { - if (lng_indexer_h(m) > 0) { - for (int mm = 0; mm < m; ++mm) { - if (lng_indexer_h(mm) == lng_indexer_h(m)) { - break; - } - ++numj; - } - } - } - - // allocate the photolysis table - auto table = mam4::mo_photo::create_photo_table_data(nw, nt, np_xs, numj, - nump, numsza, numcolo3, - numalb); - - // allocate host views for table data - auto rsf_tab_h = Kokkos::create_mirror_view(table.rsf_tab); - auto xsqy_h = Kokkos::create_mirror_view(table.xsqy); - auto sza_h = Kokkos::create_mirror_view(table.sza); - auto alb_h = Kokkos::create_mirror_view(table.alb); - auto press_h = Kokkos::create_mirror_view(table.press); - auto colo3_h = Kokkos::create_mirror_view(table.colo3); - auto o3rat_h = Kokkos::create_mirror_view(table.o3rat); - auto etfphot_h = Kokkos::create_mirror_view(table.etfphot); - auto prs_h = Kokkos::create_mirror_view(table.prs); - - if (comm.rank() == mpi_root) { // read data from files and broadcast - // read file data into our host views - read_nc_var(rsf_file, rsf_id, "pm", press_h); - read_nc_var(rsf_file, rsf_id, "sza", sza_h); - read_nc_var(rsf_file, rsf_id, "alb", alb_h); - read_nc_var(rsf_file, rsf_id, "colo3fact", o3rat_h); - read_nc_var(rsf_file, rsf_id, "colo3", colo3_h); - read_nc_var(rsf_file, rsf_id, "RSF", rsf_tab_h); - - read_nc_var(xs_long_file, xs_long_id, "pressure", prs_h); - - // read xsqy data (using lng_indexer_h for the first index) - int ndx = 0; - for (int m = 0; m < mam4::mo_photo::phtcnt; ++m) { - if (lng_indexer_h(m) > 0) { - auto xsqy_ndx_h = ekat::subview(xsqy_h, ndx); - read_nc_var(xs_long_file, xs_long_id, lng_indexer_h(m), xsqy_ndx_h); - ++ndx; - } - } - - // populate etfphot by rebinning solar data - HostView1D wc_h("wc", nw), wlintv_h("wlintv", nw), we_h("we", nw+1); - read_nc_var(rsf_file, rsf_id, "wc", wc_h); - read_nc_var(rsf_file, rsf_id, "wlintv", wlintv_h); - for (int i = 0; i < nw; ++i) { - we_h(i) = wc_h(i) - 0.5 * wlintv_h(i); - } - we_h(nw) = wc_h(nw-1) - 0.5 * wlintv_h(nw-1); - populate_etfphot(we_h, etfphot_h); - - // close the files - nc_close(rsf_id); - nc_close(xs_long_id); - } - - // broadcast host views from MPI root to others - comm.broadcast(rsf_tab_h.data(), nw*numalb*numcolo3*numsza*nump, mpi_root); - comm.broadcast(xsqy_h.data(), numj*nw*nt*np_xs, mpi_root); - comm.broadcast(sza_h.data(), numsza, mpi_root); - comm.broadcast(alb_h.data(), numalb, mpi_root); - comm.broadcast(press_h.data(), nump, mpi_root); - comm.broadcast(o3rat_h.data(), numcolo3, mpi_root); - comm.broadcast(colo3_h.data(), nump, mpi_root); - comm.broadcast(etfphot_h.data(), nw, mpi_root); - comm.broadcast(prs_h.data(), np_xs, mpi_root); - - // copy host photolysis table into place on device - Kokkos::deep_copy(table.rsf_tab, rsf_tab_h); - Kokkos::deep_copy(table.xsqy, xsqy_h); - Kokkos::deep_copy(table.sza, sza_h); - Kokkos::deep_copy(table.alb, alb_h); - Kokkos::deep_copy(table.press, press_h); - Kokkos::deep_copy(table.colo3, colo3_h); - Kokkos::deep_copy(table.o3rat, o3rat_h); - Kokkos::deep_copy(table.etfphot, etfphot_h); - Kokkos::deep_copy(table.prs, prs_h); - Kokkos::deep_copy(table.pht_alias_mult_1, pht_alias_mult_1_h); - Kokkos::deep_copy(table.lng_indexer, lng_indexer_h); - - // compute gradients (on device) - Kokkos::parallel_for("del_p", nump-1, KOKKOS_LAMBDA(int i) { - table.del_p(i) = 1.0/::abs(table.press(i)- table.press(i+1)); - }); - Kokkos::parallel_for("del_sza", numsza-1, KOKKOS_LAMBDA(int i) { - table.del_sza(i) = 1.0/(table.sza(i+1) - table.sza(i)); - }); - Kokkos::parallel_for("del_alb", numalb-1, KOKKOS_LAMBDA(int i) { - table.del_alb(i) = 1.0/(table.alb(i+1) - table.alb(i)); - }); - Kokkos::parallel_for("del_o3rat", numcolo3-1, KOKKOS_LAMBDA(int i) { - table.del_o3rat(i) = 1.0/(table.o3rat(i+1) - table.o3rat(i)); - }); - Kokkos::parallel_for("dprs", np_xs-1, KOKKOS_LAMBDA(int i) { - table.dprs(i) = 1.0/(table.prs(i) - table.prs(i+1)); - }); - - return table; -} - -// performs gas phase chemistry calculations on a single level of a single -// atmospheric column -KOKKOS_INLINE_FUNCTION -void gas_phase_chemistry(Real zm, Real zi, Real phis, Real temp, Real pmid, Real pdel, Real dt, - const Real photo_rates[mam4::mo_photo::phtcnt], // in - const Real extfrc[mam4::gas_chemistry::extcnt], // in - Real q[mam4::gas_chemistry::gas_pcnst], // VMRs, inout - Real invariants[mam4::gas_chemistry::nfs]) { // out - // constexpr Real rga = 1.0/haero::Constants::gravity; - // constexpr Real m2km = 0.01; // converts m -> km - - // The following things are chemical mechanism dependent! See mam4xx/src/mam4xx/gas_chem_mechanism.hpp) - constexpr int gas_pcnst = mam4::gas_chemistry::gas_pcnst; // number of gas phase species - constexpr int rxntot = mam4::gas_chemistry::rxntot; // number of chemical reactions - constexpr int extcnt = mam4::gas_chemistry::extcnt; // number of species with external forcing - constexpr int indexm = 0; // FIXME: index of total atm density in invariants array - - constexpr int phtcnt = mam4::mo_photo::phtcnt; // number of photolysis reactions - - constexpr int itermax = mam4::gas_chemistry::itermax; - constexpr int clscnt4 = mam4::gas_chemistry::clscnt4; - - // NOTE: vvv these arrays were copied from mam4xx/gas_chem_mechanism.hpp vvv - constexpr int permute_4[gas_pcnst] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, - 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, - 20, 21, 22, 23, 24, 25, 26, 27, 28, 29}; - constexpr int clsmap_4[gas_pcnst] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, - 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, - 21, 22, 23, 24, 25, 26, 27, 28, 29, 30}; - - // These indices for species are fixed by the chemical mechanism - // std::string solsym[] = {"O3", "H2O2", "H2SO4", "SO2", "DMS", "SOAG", - // "so4_a1", "pom_a1", "soa_a1", "bc_a1", "dst_a1", - // "ncl_a1", "mom_a1", "num_a1", "so4_a2", "soa_a2", - // "ncl_a2", "mom_a2", "num_a2", "dst_a3", "ncl_a3", - // "so4_a3", "bc_a3", "pom_a3", "soa_a3", "mom_a3", - // "num_a3", "pom_a4", "bc_a4", "mom_a4", "num_a4"}; - constexpr int ndx_h2so4 = 2; - // std::string extfrc_list[] = {"SO2", "so4_a1", "so4_a2", "pom_a4", "bc_a4", - // "num_a1", "num_a2", "num_a3", "num_a4", "SOAG"}; - constexpr int synoz_ndx = -1; - - // fetch the zenith angle (not its cosine!) in degrees for this column. - // FIXME: For now, we fix the zenith angle. At length, we need to compute it - // FIXME: from EAMxx's current set of orbital parameters, which requires some - // FIXME: conversation with the EAMxx team. - - // xform geopotential height from m to km and pressure from Pa to mb - // Real zsurf = rga * phis; - // Real zmid = m2km * (zm + zsurf); - - // ... compute the column's invariants - // Real h2ovmr = q[0]; - // setinv(invariants, temp, h2ovmr, q, pmid); FIXME: not ported yet - - // ... set rates for "tabular" and user specified reactions - Real reaction_rates[rxntot]; - mam4::gas_chemistry::setrxt(reaction_rates, temp); - - // set reaction rates based on chemical invariants - // (indices (ndxes?) are taken from mam4 validation data and translated from - // 1-based indices to 0-based indices) - int usr_HO2_HO2_ndx = 1, usr_DMS_OH_ndx = 5, - usr_SO2_OH_ndx = 3, inv_h2o_ndx = 3; - mam4::gas_chemistry::usrrxt(reaction_rates, temp, invariants, invariants[indexm], - usr_HO2_HO2_ndx, usr_DMS_OH_ndx, - usr_SO2_OH_ndx, inv_h2o_ndx); - mam4::gas_chemistry::adjrxt(reaction_rates, invariants, invariants[indexm]); - - //=================================== - // Photolysis rates at time = t(n+1) - //=================================== - - // compute the rate of change from forcing - Real extfrc_rates[extcnt]; // [1/cm^3/s] - for (int mm = 0; mm < extcnt; ++mm) { - if (mm != synoz_ndx) { - extfrc_rates[mm] = extfrc[mm] / invariants[indexm]; - } - } - - // ... Form the washout rates - Real het_rates[gas_pcnst]; - // FIXME: not ported yet - //sethet(het_rates, pmid, zmid, phis, temp, cmfdqr, prain, nevapr, delt, - // invariants[indexm], q); - - - // save h2so4 before gas phase chem (for later new particle nucleation) - Real del_h2so4_gasprod = q[ndx_h2so4]; - - //=========================== - // Class solution algorithms - //=========================== - - // copy photolysis rates into reaction_rates (assumes photolysis rates come first) - for (int i = 0; i < phtcnt; ++i) { - reaction_rates[i] = photo_rates[i]; - } - - // ... solve for "Implicit" species - bool factor[itermax]; - for (int i = 0; i < itermax; ++i) { - factor[i] = true; - } - - // initialize error tolerances - Real epsilon[clscnt4]; - mam4::gas_chemistry::imp_slv_inti(epsilon); - - // solve chemical system implicitly - Real prod_out[clscnt4], loss_out[clscnt4]; - mam4::gas_chemistry::imp_sol(q, reaction_rates, het_rates, extfrc_rates, dt, - permute_4, clsmap_4, factor, epsilon, prod_out, loss_out); - - // save h2so4 change by gas phase chem (for later new particle nucleation) - if (ndx_h2so4 > 0) { - del_h2so4_gasprod = q[ndx_h2so4] - del_h2so4_gasprod; - } -} - -} // namespace scream::impl diff --git a/components/eamxx/src/physics/mam/impl/mam4_amicphys.cpp b/components/eamxx/src/physics/mam/impl/mam4_amicphys.cpp deleted file mode 100644 index 71cee1eeb40..00000000000 --- a/components/eamxx/src/physics/mam/impl/mam4_amicphys.cpp +++ /dev/null @@ -1,1769 +0,0 @@ -#include -#include -#include -#include -#include -#include - -namespace scream::impl { - -#define MAX_FILENAME_LEN 256 - -using namespace mam4; - -// number of constituents in gas chemistry "work arrays" -KOKKOS_INLINE_FUNCTION -constexpr int gas_pcnst() { - constexpr int gas_pcnst_ = mam4::gas_chemistry::gas_pcnst; - return gas_pcnst_; -} - -// number of aerosol/gas species tendencies -KOKKOS_INLINE_FUNCTION -constexpr int nqtendbb() { return 4; } - -// MAM4 aerosol microphysics configuration data -struct AmicPhysConfig { - // these switches activate various aerosol microphysics processes - bool do_cond; // condensation (a.k.a gas-aerosol exchange) - bool do_rename; // mode "renaming" - bool do_newnuc; // gas -> aerosol nucleation - bool do_coag; // aerosol coagulation - - // configurations for specific aerosol microphysics - mam4::GasAerExchProcess::ProcessConfig condensation; - mam4::NucleationProcess::ProcessConfig nucleation; - - // controls treatment of h2so4 condensation in mam_gasaerexch_1subarea - // 1 = sequential calc. of gas-chem prod then condensation loss - // 2 = simultaneous calc. of gas-chem prod and condensation loss - int gaexch_h2so4_uptake_optaa; - - // controls how nucleation interprets h2so4 concentrations - int newnuc_h2so4_conc_optaa; -}; - -namespace { - -KOKKOS_INLINE_FUNCTION constexpr int nqtendaa() { return 5; } -KOKKOS_INLINE_FUNCTION constexpr int nqqcwtendaa() { return 1; } -KOKKOS_INLINE_FUNCTION constexpr int nqqcwtendbb() { return 1; } -KOKKOS_INLINE_FUNCTION constexpr int iqtend_cond() { return 0; } -KOKKOS_INLINE_FUNCTION constexpr int iqtend_rnam() { return 1; } -KOKKOS_INLINE_FUNCTION constexpr int iqtend_nnuc() { return 2; } -KOKKOS_INLINE_FUNCTION constexpr int iqtend_coag() { return 3; } -KOKKOS_INLINE_FUNCTION constexpr int iqtend_cond_only() { return 4; } -KOKKOS_INLINE_FUNCTION constexpr int iqqcwtend_rnam() { return 0; } -KOKKOS_INLINE_FUNCTION constexpr int maxsubarea() { return 2; } - -// conversion factors -KOKKOS_INLINE_FUNCTION Real fcvt_gas(int gas_id) { - static const Real fcvt_gas_[AeroConfig::num_gas_ids()] = {1, 1, 1}; - return fcvt_gas_[gas_id]; -} -KOKKOS_INLINE_FUNCTION Real fcvt_aer(int aero_id) { - static const Real fcvt_aer_[AeroConfig::num_aerosol_ids()] = {1, 1, 1, 1, 1, 1, 1}; - return fcvt_aer_[aero_id]; -} - -// leave number mix-ratios unchanged (#/kmol-air) -KOKKOS_INLINE_FUNCTION Real fcvt_num() { return 1.0; } -// factor for converting aerosol water mix-ratios from (kg/kg) to (mol/mol) -KOKKOS_INLINE_FUNCTION Real fcvt_wtr() { return 1.0; } - -KOKKOS_INLINE_FUNCTION constexpr int lmapcc_val_nul() { return 0; } -KOKKOS_INLINE_FUNCTION constexpr int lmapcc_val_gas() { return 1; } -KOKKOS_INLINE_FUNCTION constexpr int lmapcc_val_aer() { return 2; } -KOKKOS_INLINE_FUNCTION constexpr int lmapcc_val_num() { return 3; } -KOKKOS_INLINE_FUNCTION int lmapcc_all(int index) { - static const int lmapcc_all_[gas_pcnst()] = { - lmapcc_val_nul(), lmapcc_val_gas(), lmapcc_val_nul(), lmapcc_val_nul(), - lmapcc_val_gas(), lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_num(), lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_aer(), lmapcc_val_num(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_aer(), lmapcc_val_num(), lmapcc_val_aer(), lmapcc_val_aer(), - lmapcc_val_aer(), lmapcc_val_num()}; - return lmapcc_all_[index]; -} - -// Where lmapcc_val_num are defined in lmapcc_all -KOKKOS_INLINE_FUNCTION int numptr_amode(int mode) { - static const int numptr_amode_[AeroConfig::num_modes()] = {12, 17, 25, 29}; - return numptr_amode_[mode]; -} - -// Where lmapcc_val_gas are defined in lmapcc_all -KOKKOS_INLINE_FUNCTION int lmap_gas(int mode) { - static const int lmap_gas_[AeroConfig::num_modes()] = {4, 1}; - return lmap_gas_[mode]; -} -// Where lmapcc_val_aer are defined in lmapcc_all -KOKKOS_INLINE_FUNCTION int lmassptr_amode(int aero_id, int mode) { - static const int lmassptr_amode_[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()] = { - {5, 13, 18, 26}, {6, 14, 19, 27}, {7, 15, 20, 28}, {8, 16, 21, -6}, - {9, -6, 22, -6}, {10, -6, 23, -6}, {11, -6, 24, -6}}; - return lmassptr_amode_[aero_id][mode]; -} - -KOKKOS_INLINE_FUNCTION -void subarea_partition_factors( - const Real - q_int_cell_avg, // in grid cell mean interstitial aerosol mixing ratio - const Real - q_cbn_cell_avg, // in grid cell mean cloud-borne aerosol mixing ratio - const Real fcldy, // in cloudy fraction of the grid cell - const Real fclea, // in clear fraction of the grid cell - Real &part_fac_q_int_clea, // out - Real &part_fac_q_int_cldy) // out -{ - // Calculate mixing ratios of each subarea - - // cloud-borne, cloudy subarea - const Real tmp_q_cbn_cldy = q_cbn_cell_avg / fcldy; - // interstitial, cloudy subarea - const Real tmp_q_int_cldy = - haero::max(0.0, ((q_int_cell_avg + q_cbn_cell_avg) - tmp_q_cbn_cldy)); - // interstitial, clear subarea - const Real tmp_q_int_clea = (q_int_cell_avg - fcldy * tmp_q_int_cldy) / fclea; - - // Calculate the corresponding paritioning factors for interstitial aerosols - // using the above-derived subarea mixing ratios plus the constraint that - // the cloud fraction weighted average of subarea mean need to match grid box - // mean. - - // *** question *** - // use same part_fac_q_int_clea/cldy for everything ? - // use one for number and one for all masses (based on total mass) ? - // use separate ones for everything ? - // maybe one for number and one for all masses is best, - // because number and mass have different activation fractions - // *** question *** - - Real tmp_aa = haero::max(1.e-35, tmp_q_int_clea * fclea) / - haero::max(1.e-35, q_int_cell_avg); - tmp_aa = haero::max(0.0, haero::min(1.0, tmp_aa)); - - part_fac_q_int_clea = tmp_aa / fclea; - part_fac_q_int_cldy = (1.0 - tmp_aa) / fcldy; -} - -KOKKOS_INLINE_FUNCTION -void construct_subareas_1gridcell( - const Real cld, // in - const Real relhumgcm, // in - const Real q_pregaschem[gas_pcnst()], // in q TMRs before - // gas-phase chemistry - const Real q_precldchem[gas_pcnst()], // in q TMRs before - // cloud chemistry - const Real qqcw_precldchem[gas_pcnst()], // in qqcw TMRs before - // cloud chemistry - const Real q[gas_pcnst()], // in current tracer mixing ratios (TMRs) - // *** MUST BE #/kmol-air for number - // *** MUST BE mol/mol-air for mass - const Real qqcw[gas_pcnst()], // in like q but for - // cloud-borner tracers - int &nsubarea, // out - int &ncldy_subarea, // out - int &jclea, // out - int &jcldy, // out - bool iscldy_subarea[maxsubarea()], // out - Real afracsub[maxsubarea()], // out - Real relhumsub[maxsubarea()], // out - Real qsub1[gas_pcnst()][maxsubarea()], // out interstitial - Real qsub2[gas_pcnst()][maxsubarea()], // out interstitial - Real qsub3[gas_pcnst()][maxsubarea()], // out interstitial - Real qqcwsub1[gas_pcnst()][maxsubarea()], // out cloud-borne - Real qqcwsub2[gas_pcnst()][maxsubarea()], // out cloud-borne - Real qqcwsub3[gas_pcnst()][maxsubarea()], // outcloud-borne - Real qaerwatsub3[AeroConfig::num_modes()] - [maxsubarea()], // out aerosol water mixing ratios (mol/mol) - Real qaerwat[AeroConfig::num_modes()] // in aerosol water mixing ratio - // (kg/kg, NOT mol/mol) -) { - static constexpr int num_modes = AeroConfig::num_modes(); - // cloud chemistry is only on when cld(i,k) >= 1.0e-5_wp - // it may be that the macrophysics has a higher threshold that this - const Real fcld_locutoff = 1.0e-5; - const Real fcld_hicutoff = 0.999; - - // qgcmN and qqcwgcmN (N=1:4) are grid-cell mean tracer mixing ratios (TMRs, - // mol/mol or #/kmol) - // N=1 - before gas-phase chemistry - // N=2 - before cloud chemistry - // N=3 - incoming values (before gas-aerosol exchange, newnuc, coag) - // qgcm1, qgcm2, qgcm3 - // qqcwgcm2, qqcwgcm3 - // qaerwatgcm3 ! aerosol water mixing ratios (mol/mol) - - // -------------------------------------------------------------------------------------- - // Determine the number of sub-areas, their fractional areas, and relative - // humidities - // -------------------------------------------------------------------------------------- - // if cloud fraction ~= 0, the grid-cell has a single clear sub-area - // (nsubarea = 1) if cloud fraction ~= 1, the grid-cell has a single cloudy - // sub-area (nsubarea = 1) otherwise, the grid-cell has a - // clear and a cloudy sub-area (nsubarea = 2) - - Real zfcldy = 0; - nsubarea = 0; - ncldy_subarea = 0; - jclea = 0; - jcldy = 0; - - if (cld < fcld_locutoff) { - nsubarea = 1; - jclea = 1; - } else if (cld > fcld_hicutoff) { - zfcldy = 1.0; - nsubarea = 1; - ncldy_subarea = 1; - jcldy = 1; - } else { - zfcldy = cld; - nsubarea = 2; - ncldy_subarea = 1; - jclea = 1; - jcldy = 2; - } - - const Real zfclea = 1.0 - zfcldy; - for (int i = 0; i < maxsubarea(); ++i) - iscldy_subarea[i] = false; - if (jcldy > 0) - iscldy_subarea[jcldy - 1] = true; - for (int i = 0; i < maxsubarea(); ++i) - afracsub[i] = 0.0; - if (jclea > 0) - afracsub[jclea - 1] = zfclea; - if (jcldy > 0) - afracsub[jcldy - 1] = zfcldy; - - // cldy_rh_sameas_clear is just to match mam_refactor. Compiler should - // optimize away. - const int cldy_rh_sameas_clear = 0; - if (ncldy_subarea <= 0) { - for (int i = 0; i < maxsubarea(); ++i) - relhumsub[i] = relhumgcm; - } else if (cldy_rh_sameas_clear > 0) { - for (int i = 0; i < maxsubarea(); ++i) - relhumsub[i] = relhumgcm; - } else { - if (jcldy > 0) { - relhumsub[jcldy - 1] = 1.0; - if (jclea > 0) { - const Real tmpa = - (relhumgcm - afracsub[jcldy - 1]) / afracsub[jclea - 1]; - relhumsub[jclea - 1] = haero::max(0.0, haero::min(1.0, tmpa)); - } - } - } - - // ---------------------------------------------------------------------------- - // Copy grid cell mean mixing ratios. - // These values, together with cloud fraction and a few assumptions, are used - // in the remainder of the subroutine to calculate the sub-area mean mixing - // ratios. - // ---------------------------------------------------------------------------- - // Interstitial aerosols - Real qgcm1[gas_pcnst()], qgcm2[gas_pcnst()], qgcm3[gas_pcnst()]; - for (int i = 0; i < gas_pcnst(); ++i) { - qgcm1[i] = haero::max(0.0, q_pregaschem[i]); - qgcm2[i] = haero::max(0.0, q_precldchem[i]); - qgcm3[i] = haero::max(0.0, q[i]); - } - - // Cloud-borne aerosols - Real qqcwgcm2[gas_pcnst()], qqcwgcm3[gas_pcnst()]; - for (int i = 0; i < gas_pcnst(); ++i) { - qqcwgcm2[i] = haero::max(0.0, qqcw_precldchem[i]); - qqcwgcm3[i] = haero::max(0.0, qqcw[i]); - } - - // aerosol water - Real qaerwatgcm3[num_modes] = {}; - for (int i = 0; i < num_modes; ++i) { - qaerwatgcm3[i] = haero::max(0.0, qaerwat[i]); - } - - // ---------------------------------------------------------------------------- - // Initialize the subarea mean mixing ratios - // ---------------------------------------------------------------------------- - { - const int n = haero::min(maxsubarea(), nsubarea + 1); - for (int i = 0; i < n; ++i) { - for (int j = 0; j < gas_pcnst(); ++j) { - qsub1[j][i] = 0.0; - qsub2[j][i] = 0.0; - qsub3[j][i] = 0.0; - qqcwsub1[j][i] = 0.0; - qqcwsub2[j][i] = 0.0; - qqcwsub3[j][i] = 0.0; - } - for (int j = 0; j < num_modes; ++j) { - qaerwatsub3[j][i] = 0.0; - } - } - } - - // ************************************************************************************************* - // Calculate initial (i.e., before cond/rnam/nnuc/coag) tracer mixing - // ratios within the sub-areas - // - for all-clear or all-cloudy cases, the sub-area TMRs are equal to the - // grid-cell means - // - for partly cloudy case, they are different. This is primarily - // because the - // interstitial aerosol mixing ratios are assumed lower in the cloudy - // sub-area than in the clear sub-area, because much of the aerosol is - // activated in the cloudy sub-area. - // ************************************************************************************************* - // Category I: partly cloudy case - // ************************************************************************************************* - if ((jclea > 0) && (jcldy > 0) && (jclea + jcldy == 3) && (nsubarea == 2)) { - - // --------------------------------------------------------------------- - // Set GAS mixing ratios in sub-areas (for the condensing gases only!!) - // --------------------------------------------------------------------- - for (int lmz = 0; lmz < gas_pcnst(); ++lmz) { - if (lmapcc_all(lmz) == lmapcc_val_gas()) { - - // assume gas in both sub-areas before gas-chem and cloud-chem equal - // grid-cell mean - for (int i = 0; i < nsubarea; ++i) { - qsub1[lmz][i] = qgcm1[lmz]; - qsub2[lmz][i] = qgcm2[lmz]; - } - // assume gas in clear sub-area after cloud-chem equals before - // cloud-chem value - qsub3[lmz][jclea - 1] = qsub2[lmz][jclea - 1]; - // gas in cloud sub-area then determined by grid-cell mean and clear - // values - qsub3[lmz][jcldy - 1] = - (qgcm3[lmz] - zfclea * qsub3[lmz][jclea - 1]) / zfcldy; - - // check that this does not produce a negative value - if (qsub3[lmz][jcldy - 1] < 0.0) { - qsub3[lmz][jcldy - 1] = 0.0; - qsub3[lmz][jclea - 1] = qgcm3[lmz] / zfclea; - } - } - } - // --------------------------------------------------------------------- - // Set CLOUD-BORNE AEROSOL mixing ratios in sub-areas. - // This is straightforward, as the same partitioning factors (0 or 1/f) - // are applied to all mass and number mixing ratios in all modes. - // --------------------------------------------------------------------- - // loop thru log-normal modes - for (int n = 0; n < num_modes; ++n) { - // number - then mass of individual species - of a mode - for (int l2 = -1; l2 < num_species_mode(n); ++l2) { - int lc; - if (l2 == -1) - lc = numptr_amode(n); - else - lc = lmassptr_amode(l2, n); - qqcwsub2[lc][jclea - 1] = 0.0; - qqcwsub2[lc][jcldy - 1] = qqcwgcm2[lc] / zfcldy; - qqcwsub3[lc][jclea - 1] = 0.0; - qqcwsub3[lc][jcldy - 1] = qqcwgcm3[lc] / zfcldy; - } - } - - // --------------------------------------------------------------------- - // Set INTERSTITIAL AEROSOL mixing ratios in sub-areas. - // --------------------------------------------------------------------- - for (int n = 0; n < num_modes; ++n) { - // ------------------------------------- - // Aerosol number - // ------------------------------------- - // grid cell mean, interstitial - Real tmp_q_cellavg_int = qgcm2[numptr_amode(n)]; - // grid cell mean, cloud-borne - Real tmp_q_cellavg_cbn = qqcwgcm2[numptr_amode(n)]; - - Real nmbr_part_fac_clea = 0; - Real nmbr_part_fac_cldy = 0; - subarea_partition_factors(tmp_q_cellavg_int, tmp_q_cellavg_cbn, zfcldy, - zfclea, nmbr_part_fac_clea, nmbr_part_fac_cldy); - - // Apply the partitioning factors to calculate sub-area mean number - // mixing ratios - - const int la = numptr_amode(n); - - qsub2[la][jclea - 1] = qgcm2[la] * nmbr_part_fac_clea; - qsub2[la][jcldy - 1] = qgcm2[la] * nmbr_part_fac_cldy; - qsub3[la][jclea - 1] = qgcm3[la] * nmbr_part_fac_clea; - qsub3[la][jcldy - 1] = qgcm3[la] * nmbr_part_fac_cldy; - - //------------------------------------- - // Aerosol mass - //------------------------------------- - // For aerosol mass, we use the total grid cell mean - // interstitial/cloud-borne mass mixing ratios to come up with the same - // partitioning for all species in the mode. - - // Compute the total mixing ratios by summing up the individual species - - tmp_q_cellavg_int = 0.0; // grid cell mean, interstitial - tmp_q_cellavg_cbn = 0.0; // grid cell mean, cloud-borne - - for (int l2 = 0; l2 < num_species_mode(n); ++l2) { - tmp_q_cellavg_int += qgcm2[lmassptr_amode(l2, n)]; - tmp_q_cellavg_cbn += qqcwgcm2[lmassptr_amode(l2, n)]; - } - Real mass_part_fac_clea = 0; - Real mass_part_fac_cldy = 0; - // Calculate the partitioning factors - subarea_partition_factors(tmp_q_cellavg_int, tmp_q_cellavg_cbn, zfcldy, - zfclea, mass_part_fac_clea, mass_part_fac_cldy); - - // Apply the partitioning factors to calculate sub-area mean mass mixing - // ratios - - for (int l2 = 0; l2 < num_species_mode(n); ++l2) { - const int la = lmassptr_amode(l2, n); - - qsub2[la][jclea - 1] = qgcm2[la] * mass_part_fac_clea; - qsub2[la][jcldy - 1] = qgcm2[la] * mass_part_fac_cldy; - qsub3[la][jclea - 1] = qgcm3[la] * mass_part_fac_clea; - qsub3[la][jcldy - 1] = qgcm3[la] * mass_part_fac_cldy; - } - } - - // ************************************************************************************************* - // Category II: all clear, or cld < 1e-5 - // In this case, zfclea=1 and zfcldy=0 - // ************************************************************************************************* - } else if ((jclea == 1) && (jcldy == 0) && (nsubarea == 1)) { - // - // put all the gases and interstitial aerosols in the clear sub-area - // and set mix-ratios = 0 in cloudy sub-area - // for cloud-borne aerosol, do nothing - // because the grid-cell-mean cloud-borne aerosol will be left - // unchanged (i.e., this routine only changes qqcw when cld >= 1e-5) - // - - for (int lmz = 0; lmz < gas_pcnst(); ++lmz) { - if (0 < lmapcc_all(lmz)) { - qsub1[lmz][jclea - 1] = qgcm1[lmz]; - qsub2[lmz][jclea - 1] = qgcm2[lmz]; - qsub3[lmz][jclea - 1] = qgcm3[lmz]; - qqcwsub2[lmz][jclea - 1] = qqcwgcm2[lmz]; - qqcwsub3[lmz][jclea - 1] = qqcwgcm3[lmz]; - } - } - // ************************************************************************************************* - // Category III: all cloudy, or cld > 0.999 - // in this case, zfcldy= and zfclea=0 - // ************************************************************************************************* - } else if ((jclea == 0) && (jcldy == 1) && (nsubarea == 1)) { - // - // put all the gases and interstitial aerosols in the cloudy sub-area - // and set mix-ratios = 0 in clear sub-area - // - for (int lmz = 0; lmz < gas_pcnst(); ++lmz) { - if (0 < lmapcc_all(lmz)) { - qsub1[lmz][jcldy - 1] = qgcm1[lmz]; - qsub2[lmz][jcldy - 1] = qgcm2[lmz]; - qsub3[lmz][jcldy - 1] = qgcm3[lmz]; - qqcwsub2[lmz][jcldy - 1] = qqcwgcm2[lmz]; - qqcwsub3[lmz][jcldy - 1] = qqcwgcm3[lmz]; - } - } - // ************************************************************************************************* - } else { // this should not happen - EKAT_KERNEL_REQUIRE_MSG(false, "*** modal_aero_amicphys - bad jclea, jcldy, nsubarea!"); - } - // ************************************************************************************************* - - // ------------------------------------------------------------------------------------ - // aerosol water -- how to treat this in sub-areas needs more work/thinking - // currently modal_aero_water_uptake calculates qaerwat using - // the grid-cell mean interstital-aerosol mix-rats and the clear-area rh - for (int jsub = 0; jsub < nsubarea; ++jsub) - for (int i = 0; i < num_modes; ++i) - qaerwatsub3[i][jsub] = qaerwatgcm3[i]; - - // ------------------------------------------------------------------------------------ - if (nsubarea == 1) { - // the j=1 subarea is used for some diagnostics - // but is not used in actual calculations - const int j = 1; - for (int i = 0; i < gas_pcnst(); ++i) { - qsub1[i][j] = 0.0; - qsub2[i][j] = 0.0; - qsub3[i][j] = 0.0; - qqcwsub2[i][j] = 0.0; - qqcwsub3[i][j] = 0.0; - } - } -} - -KOKKOS_INLINE_FUNCTION -void mam_amicphys_1subarea_clear( - const AmicPhysConfig& config, const int nstep, const Real deltat, const int jsub, - const int nsubarea, const bool iscldy_subarea, const Real afracsub, - const Real temp, const Real pmid, const Real pdel, const Real zmid, - const Real pblh, const Real relhum, Real dgn_a[AeroConfig::num_modes()], - Real dgn_awet[AeroConfig::num_modes()], - Real wetdens[AeroConfig::num_modes()], - const Real qgas1[AeroConfig::num_gas_ids()], - const Real qgas3[AeroConfig::num_gas_ids()], - Real qgas4[AeroConfig::num_gas_ids()], - Real qgas_delaa[AeroConfig::num_gas_ids()][nqtendaa()], - const Real qnum3[AeroConfig::num_modes()], - Real qnum4[AeroConfig::num_modes()], - Real qnum_delaa[AeroConfig::num_modes()][nqtendaa()], - const Real qaer3[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()], - Real qaer4[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()], - Real qaer_delaa[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()] - [nqtendaa()], - const Real qwtr3[AeroConfig::num_modes()], - Real qwtr4[AeroConfig::num_modes()]) { - static constexpr int num_gas_ids = AeroConfig::num_gas_ids(); - static constexpr int num_modes = AeroConfig::num_modes(); - static constexpr int num_aerosol_ids = AeroConfig::num_aerosol_ids(); - - static constexpr int igas_h2so4 = static_cast(GasId::H2SO4); - // Turn off nh3 for now. This is a future enhancement. - static constexpr int igas_nh3 = -999888777; // Same as mam_refactor - static constexpr int iaer_so4 = static_cast(AeroId::SO4); - static constexpr int iaer_pom = static_cast(AeroId::POM); - - const AeroId gas_to_aer[num_gas_ids] = {AeroId::SOA, AeroId::SO4, - AeroId::None}; - - const bool l_gas_condense_to_mode[num_gas_ids][num_modes] = { - {true, true, true, true}, - {true, true, true, true}, - {false, false, false, false}}; - enum { NA, ANAL, IMPL }; - const int eqn_and_numerics_category[num_gas_ids] = {IMPL, ANAL, ANAL}; - - // air molar density (kmol/m3) - // const Real r_universal = Constants::r_gas; // [mJ/(K mol)] - const Real r_universal = 8.314467591; // [mJ/(mol)] as in mam_refactor - const Real aircon = pmid / (1000 * r_universal * temp); - const Real alnsg_aer[num_modes] = {0.58778666490211906, 0.47000362924573563, - 0.58778666490211906, 0.47000362924573563}; - const Real uptk_rate_factor[num_gas_ids] = {0.81, 1.0, 1.0}; - // calculates changes to gas and aerosol sub-area TMRs (tracer mixing ratios) - // qgas3, qaer3, qnum3 are the current incoming TMRs - // qgas4, qaer4, qnum4 are the updated outgoing TMRs - // - // this routine calculates changes involving - // gas-aerosol exchange (condensation/evaporation) - // growth from smaller to larger modes (renaming) due to condensation - // new particle nucleation - // coagulation - // transfer of particles from hydrophobic modes to hydrophilic modes - // (aging) - // due to condensation and coagulation - // - // qXXXN (X=gas,aer,wat,num; N=1:4) are sub-area mixing ratios - // XXX=gas - gas species - // XXX=aer - aerosol mass species (excluding water) - // XXX=wat - aerosol water - // XXX=num - aerosol number - // N=1 - before gas-phase chemistry - // N=2 - before cloud chemistry - // N=3 - current incoming values (before gas-aerosol exchange, newnuc, - // coag) N=4 - updated outgoing values (after gas-aerosol exchange, - // newnuc, coag) - // - // qXXX_delaa are TMR changes (not tendencies) - // for different processes, which are used to produce history output - // for a clear sub-area, the processes are condensation/evaporation (and - // associated aging), renaming, coagulation, and nucleation - - Real qgas_cur[num_gas_ids]; - for (int i = 0; i < num_gas_ids; ++i) - qgas_cur[i] = qgas3[i]; - Real qaer_cur[num_aerosol_ids][num_modes]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_cur[i][j] = qaer3[i][j]; - - Real qnum_cur[num_modes]; - for (int j = 0; j < num_modes; ++j) - qnum_cur[j] = qnum3[j]; - Real qwtr_cur[num_modes]; - for (int j = 0; j < num_modes; ++j) - qwtr_cur[j] = qwtr3[j]; - - // qgas_netprod_otrproc = gas net production rate from other processes - // such as gas-phase chemistry and emissions (mol/mol/s) - // this allows the condensation (gasaerexch) routine to apply production and - // condensation loss - // together, which is more accurate numerically - // NOTE - must be >= zero, as numerical method can fail when it is negative - // NOTE - currently only the values for h2so4 and nh3 should be non-zero - Real qgas_netprod_otrproc[num_gas_ids] = {}; - if (config.do_cond && config.gaexch_h2so4_uptake_optaa == 2) { - for (int igas = 0; igas < num_gas_ids; ++igas) { - if (igas == igas_h2so4 || igas == igas_nh3) { - // if config.gaexch_h2so4_uptake_optaa == 2, then - // if qgas increases from pre-gaschem to post-cldchem, - // start from the pre-gaschem mix-ratio and add in the production - // during the integration - // if it decreases, - // start from post-cldchem mix-ratio - // *** currently just do this for h2so4 and nh3 - qgas_netprod_otrproc[igas] = (qgas3[igas] - qgas1[igas]) / deltat; - if (qgas_netprod_otrproc[igas] >= 0.0) - qgas_cur[igas] = qgas1[igas]; - else - qgas_netprod_otrproc[igas] = 0.0; - } - } - } - Real qgas_del_cond[num_gas_ids] = {}; - Real qgas_del_nnuc[num_gas_ids] = {}; - Real qgas_del_cond_only[num_gas_ids] = {}; - Real qaer_del_cond[num_aerosol_ids][num_modes] = {}; - Real qaer_del_rnam[num_aerosol_ids][num_modes] = {}; - Real qaer_del_nnuc[num_aerosol_ids][num_modes] = {}; - Real qaer_del_coag[num_aerosol_ids][num_modes] = {}; - Real qaer_delsub_coag_in[num_aerosol_ids][AeroConfig::max_agepair()] = {}; - Real qaer_delsub_cond[num_aerosol_ids][num_modes] = {}; - Real qaer_delsub_coag[num_aerosol_ids][num_modes] = {}; - Real qaer_del_cond_only[num_aerosol_ids][num_modes] = {}; - Real qnum_del_cond[num_modes] = {}; - Real qnum_del_rnam[num_modes] = {}; - Real qnum_del_nnuc[num_modes] = {}; - Real qnum_del_coag[num_modes] = {}; - Real qnum_delsub_cond[num_modes] = {}; - Real qnum_delsub_coag[num_modes] = {}; - Real qnum_del_cond_only[num_modes] = {}; - Real dnclusterdt = 0.0; - - const int ntsubstep = 1; - Real dtsubstep = deltat; - if (ntsubstep > 1) - dtsubstep = deltat / ntsubstep; - Real del_h2so4_gasprod = - haero::max(qgas3[igas_h2so4] - qgas1[igas_h2so4], 0.0) / ntsubstep; - - // loop over multiple time sub-steps - for (int jtsubstep = 1; jtsubstep <= ntsubstep; ++jtsubstep) { - // gas-aerosol exchange - Real uptkrate_h2so4 = 0.0; - Real del_h2so4_aeruptk = 0.0; - Real qaer_delsub_grow4rnam[num_aerosol_ids][num_modes] = {}; - Real qgas_avg[num_gas_ids] = {}; - Real qnum_sv1[num_modes] = {}; - Real qaer_sv1[num_aerosol_ids][num_modes] = {}; - Real qgas_sv1[num_gas_ids] = {}; - - if (config.do_cond) { - - const bool l_calc_gas_uptake_coeff = jtsubstep == 1; - Real uptkaer[num_gas_ids][num_modes] = {}; - - for (int i = 0; i < num_gas_ids; ++i) - qgas_sv1[i] = qgas_cur[i]; - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_sv1[j][i] = qaer_cur[j][i]; - - // time sub-step - const Real dtsub_soa_fixed = -1.0; - // Integration order - const int nghq = 2; - const int ntot_soamode = 4; - int niter_out = 0; - Real g0_soa_out = 0; - gasaerexch::mam_gasaerexch_1subarea( - nghq, igas_h2so4, igas_nh3, ntot_soamode, gas_to_aer, iaer_so4, - iaer_pom, l_calc_gas_uptake_coeff, l_gas_condense_to_mode, - eqn_and_numerics_category, dtsubstep, dtsub_soa_fixed, temp, pmid, - aircon, num_gas_ids, qgas_cur, qgas_avg, qgas_netprod_otrproc, - qaer_cur, qnum_cur, dgn_awet, alnsg_aer, uptk_rate_factor, uptkaer, - uptkrate_h2so4, niter_out, g0_soa_out); - - if (config.newnuc_h2so4_conc_optaa == 11) - qgas_avg[igas_h2so4] = - 0.5 * (qgas_sv1[igas_h2so4] + qgas_cur[igas_h2so4]); - else if (config.newnuc_h2so4_conc_optaa == 12) - qgas_avg[igas_h2so4] = qgas_cur[igas_h2so4]; - - for (int i = 0; i < num_gas_ids; ++i) - qgas_del_cond[i] += - (qgas_cur[i] - (qgas_sv1[i] + qgas_netprod_otrproc[i] * dtsubstep)); - - for (int i = 0; i < num_modes; ++i) - qnum_delsub_cond[i] = qnum_cur[i] - qnum_sv1[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_delsub_cond[i][j] = qaer_cur[i][j] - qaer_sv1[i][j]; - - // qaer_del_grow4rnam = change in qaer_del_cond during latest condensation - // calculations - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_delsub_grow4rnam[i][j] = qaer_cur[i][j] - qaer_sv1[i][j]; - for (int i = 0; i < num_gas_ids; ++i) - qgas_del_cond_only[i] = qgas_del_cond[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_del_cond_only[i][j] = qaer_delsub_cond[i][j]; - for (int i = 0; i < num_modes; ++i) - qnum_del_cond_only[i] = qnum_delsub_cond[i]; - del_h2so4_aeruptk = - qgas_cur[igas_h2so4] - - (qgas_sv1[igas_h2so4] + qgas_netprod_otrproc[igas_h2so4] * dtsubstep); - } else { - for (int i = 0; i < num_gas_ids; ++i) - qgas_avg[i] = qgas_cur[i]; - } - - // renaming after "continuous growth" - if (config.do_rename) { - constexpr int nmodes = AeroConfig::num_modes(); - constexpr int naerosol_species = AeroConfig::num_aerosol_ids(); - const Real smallest_dryvol_value = 1.0e-25; // BAD_CONSTANT - const int dest_mode_of_mode[nmodes] = {-1, 0, -1, -1}; - - Real qnumcw_cur[num_modes] = {}; - Real qaercw_cur[num_aerosol_ids][num_modes] = {}; - Real qaercw_delsub_grow4rnam[num_aerosol_ids][num_modes] = {}; - Real mean_std_dev[nmodes]; - Real fmode_dist_tail_fac[nmodes]; - Real v2n_lo_rlx[nmodes]; - Real v2n_hi_rlx[nmodes]; - Real ln_diameter_tail_fac[nmodes]; - int num_pairs = 0; - Real diameter_cutoff[nmodes]; - Real ln_dia_cutoff[nmodes]; - Real diameter_threshold[nmodes]; - Real mass_2_vol[naerosol_species] = {0.15, - 6.4971751412429377e-002, - 0.15, - 7.0588235294117650e-003, - 3.0789473684210526e-002, - 5.1923076923076926e-002, - 156.20986883198000}; - - rename::find_renaming_pairs(dest_mode_of_mode, // in - mean_std_dev, // out - fmode_dist_tail_fac, // out - v2n_lo_rlx, // out - v2n_hi_rlx, // out - ln_diameter_tail_fac, // out - num_pairs, // out - diameter_cutoff, // out - ln_dia_cutoff, diameter_threshold); - - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_sv1[j][i] = qaer_cur[j][i]; - Real dgnum_amode[nmodes]; - for (int m = 0; m < nmodes; ++m) { - dgnum_amode[m] = modes(m).nom_diameter; - } - - { - Real qmol_i_cur[num_modes][num_aerosol_ids]; - Real qmol_i_del[num_modes][num_aerosol_ids]; - Real qmol_c_cur[num_modes][num_aerosol_ids]; - Real qmol_c_del[num_modes][num_aerosol_ids]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) { - qmol_i_cur[i][j] = qaer_cur[j][i]; - qmol_i_del[i][j] = qaer_delsub_grow4rnam[j][i]; - qmol_c_cur[i][j] = qaercw_cur[j][i]; - qmol_c_del[i][j] = qaercw_delsub_grow4rnam[j][i]; - } - Rename rename; - rename.mam_rename_1subarea_( - iscldy_subarea, smallest_dryvol_value, dest_mode_of_mode, - mean_std_dev, fmode_dist_tail_fac, v2n_lo_rlx, v2n_hi_rlx, - ln_diameter_tail_fac, num_pairs, diameter_cutoff, ln_dia_cutoff, - diameter_threshold, mass_2_vol, dgnum_amode, qnum_cur, qmol_i_cur, - qmol_i_del, qnumcw_cur, qmol_c_cur, qmol_c_del); - - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) { - qaer_cur[j][i] = qmol_i_cur[i][j]; - qaer_delsub_grow4rnam[j][i] = qmol_i_del[i][j]; - qaercw_cur[j][i] = qmol_c_cur[i][j]; - qaercw_delsub_grow4rnam[j][i] = qmol_c_del[i][j]; - } - } - - for (int i = 0; i < num_modes; ++i) - qnum_del_rnam[i] += qnum_cur[i] - qnum_sv1[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_del_rnam[i][j] += qaer_cur[i][j] - qaer_sv1[i][j]; - } - - // new particle formation (nucleation) - if (config.do_newnuc) { - for (int i = 0; i < num_gas_ids; ++i) - qgas_sv1[i] = qgas_cur[i]; - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - Real qaer_cur_tmp[num_modes][num_aerosol_ids]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) { - qaer_sv1[j][i] = qaer_cur[j][i]; - qaer_cur_tmp[i][j] = qaer_cur[j][i]; - } - Real dnclusterdt_substep = 0; - Real dndt_ait = 0; - Real dmdt_ait = 0; - Real dso4dt_ait = 0; - Real dnh4dt_ait = 0; - Nucleation nucleation; - Nucleation::Config config; - config.dens_so4a_host = 1770; - config.mw_nh4a_host = 115; - config.mw_so4a_host = 115; - config.accom_coef_h2so4 = 0.65; - AeroConfig aero_config; - nucleation.init(aero_config, config); - nucleation.compute_tendencies_( - dtsubstep, temp, pmid, aircon, zmid, pblh, relhum, uptkrate_h2so4, - del_h2so4_gasprod, del_h2so4_aeruptk, qgas_cur, qgas_avg, qnum_cur, - qaer_cur_tmp, qwtr_cur, dndt_ait, dmdt_ait, dso4dt_ait, dnh4dt_ait, - dnclusterdt_substep); - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_cur[j][i] = qaer_cur_tmp[i][j]; - - //! Apply the tendencies to the prognostics. - const int nait = static_cast(ModeIndex::Aitken); - qnum_cur[nait] += dndt_ait * dtsubstep; - - if (dso4dt_ait > 0.0) { - static constexpr int iaer_so4 = static_cast(AeroId::SO4); - static constexpr int igas_h2so4 = static_cast(GasId::H2SO4); - - Real delta_q = dso4dt_ait * dtsubstep; - qaer_cur[iaer_so4][nait] += delta_q; - delta_q = haero::min(delta_q, qgas_cur[igas_h2so4]); - qgas_cur[igas_h2so4] -= delta_q; - } - - if (igas_nh3 > 0 && dnh4dt_ait > 0.0) { - static constexpr int iaer_nh4 = - -9999999; // static_cast(AeroId::NH4); - - Real delta_q = dnh4dt_ait * dtsubstep; - qaer_cur[iaer_nh4][nait] += delta_q; - delta_q = haero::min(delta_q, qgas_cur[igas_nh3]); - qgas_cur[igas_nh3] -= delta_q; - } - for (int i = 0; i < num_gas_ids; ++i) - qgas_del_nnuc[i] += (qgas_cur[i] - qgas_sv1[i]); - for (int i = 0; i < num_modes; ++i) - qnum_del_nnuc[i] += (qnum_cur[i] - qnum_sv1[i]); - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_del_nnuc[j][i] += (qaer_cur[j][i] - qaer_sv1[j][i]); - - dnclusterdt = dnclusterdt + dnclusterdt_substep * (dtsubstep / deltat); - } - - // coagulation part - if (config.do_coag) { - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_sv1[j][i] = qaer_cur[j][i]; - coagulation::mam_coag_1subarea(dtsubstep, temp, pmid, aircon, dgn_a, - dgn_awet, wetdens, qnum_cur, qaer_cur, - qaer_delsub_coag_in); - for (int i = 0; i < num_modes; ++i) - qnum_delsub_coag[i] = qnum_cur[i] - qnum_sv1[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_delsub_coag[j][i] = qaer_cur[j][i] - qaer_sv1[j][i]; - } - - // primary carbon aging - - aging::mam_pcarbon_aging_1subarea( - dgn_a, qnum_cur, qnum_delsub_cond, qnum_delsub_coag, qaer_cur, - qaer_delsub_cond, qaer_delsub_coag, qaer_delsub_coag_in); - - // accumulate sub-step q-dels - if (config.do_coag) { - for (int i = 0; i < num_modes; ++i) - qnum_del_coag[i] += qnum_delsub_coag[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_del_coag[j][i] += qaer_delsub_coag[j][i]; - } - if (config.do_cond) { - for (int i = 0; i < num_modes; ++i) - qnum_del_cond[i] += qnum_delsub_cond[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_del_cond[j][i] += qaer_delsub_cond[j][i]; - } - } - - // final mix ratios - for (int i = 0; i < num_gas_ids; ++i) - qgas4[i] = qgas_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer4[j][i] = qaer_cur[j][i]; - for (int i = 0; i < num_modes; ++i) - qnum4[i] = qnum_cur[i]; - for (int i = 0; i < num_modes; ++i) - qwtr4[i] = qwtr_cur[i]; - - // final mix ratio changes - for (int i = 0; i < num_gas_ids; ++i) { - qgas_delaa[i][iqtend_cond()] = qgas_del_cond[i]; - qgas_delaa[i][iqtend_rnam()] = 0.0; - qgas_delaa[i][iqtend_nnuc()] = qgas_del_nnuc[i]; - qgas_delaa[i][iqtend_coag()] = 0.0; - qgas_delaa[i][iqtend_cond_only()] = qgas_del_cond_only[i]; - } - for (int i = 0; i < num_modes; ++i) { - qnum_delaa[i][iqtend_cond()] = qnum_del_cond[i]; - qnum_delaa[i][iqtend_rnam()] = qnum_del_rnam[i]; - qnum_delaa[i][iqtend_nnuc()] = qnum_del_nnuc[i]; - qnum_delaa[i][iqtend_coag()] = qnum_del_coag[i]; - qnum_delaa[i][iqtend_cond_only()] = qnum_del_cond_only[i]; - } - for (int j = 0; j < num_aerosol_ids; ++j) { - for (int i = 0; i < num_modes; ++i) { - qaer_delaa[j][i][iqtend_cond()] = qaer_del_cond[j][i]; - qaer_delaa[j][i][iqtend_rnam()] = qaer_del_rnam[j][i]; - qaer_delaa[j][i][iqtend_nnuc()] = qaer_del_nnuc[j][i]; - qaer_delaa[j][i][iqtend_coag()] = qaer_del_coag[j][i]; - qaer_delaa[j][i][iqtend_cond_only()] = qaer_del_cond_only[j][i]; - } - } -} - -KOKKOS_INLINE_FUNCTION -void mam_amicphys_1subarea_cloudy( - const AmicPhysConfig& config, const int nstep, const Real deltat, const int jsub, - const int nsubarea, const bool iscldy_subarea, const Real afracsub, - const Real temp, const Real pmid, const Real pdel, const Real zmid, - const Real pblh, const Real relhum, Real dgn_a[AeroConfig::num_modes()], - Real dgn_awet[AeroConfig::num_modes()], - Real wetdens[AeroConfig::num_modes()], - const Real qgas1[AeroConfig::num_gas_ids()], - const Real qgas3[AeroConfig::num_gas_ids()], - Real qgas4[AeroConfig::num_gas_ids()], - Real qgas_delaa[AeroConfig::num_gas_ids()][nqtendaa()], - const Real qnum3[AeroConfig::num_modes()], - Real qnum4[AeroConfig::num_modes()], - Real qnum_delaa[AeroConfig::num_modes()][nqtendaa()], - const Real qaer2[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()], - const Real qaer3[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()], - Real qaer4[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()], - Real qaer_delaa[AeroConfig::num_aerosol_ids()][AeroConfig::num_modes()] - [nqtendaa()], - const Real qwtr3[AeroConfig::num_modes()], - Real qwtr4[AeroConfig::num_modes()], - const Real qnumcw3[AeroConfig::num_modes()], - Real qnumcw4[AeroConfig::num_modes()], - Real qnumcw_delaa[AeroConfig::num_modes()][nqqcwtendaa()], - const Real qaercw2[AeroConfig::num_gas_ids()][AeroConfig::num_modes()], - const Real qaercw3[AeroConfig::num_gas_ids()][AeroConfig::num_modes()], - Real qaercw4[AeroConfig::num_gas_ids()][AeroConfig::num_modes()], - Real qaercw_delaa[AeroConfig::num_gas_ids()][AeroConfig::num_modes()] - [nqqcwtendaa()]) { - - // - // calculates changes to gas and aerosol sub-area TMRs (tracer mixing ratios) - // qgas3, qaer3, qaercw3, qnum3, qnumcw3 are the current incoming TMRs - // qgas4, qaer4, qaercw4, qnum4, qnumcw4 are the updated outgoing TMRs - // - // when config.do_cond = false, this routine only calculates changes involving - // growth from smaller to larger modes (renaming) following cloud chemistry - // so gas TMRs are not changed - // when config.do_cond = true, this routine also calculates changes involving - // gas-aerosol exchange (condensation/evaporation) - // transfer of particles from hydrophobic modes to hydrophilic modes - // (aging) - // due to condensation - // currently this routine does not do - // new particle nucleation - because h2so4 gas conc. should be very low in - // cloudy air coagulation - because cloud-borne aerosol would need to be - // included - // - - // qXXXN (X=gas,aer,wat,num; N=1:4) are sub-area mixing ratios - // XXX=gas - gas species - // XXX=aer - aerosol mass species (excluding water) - // XXX=wat - aerosol water - // XXX=num - aerosol number - // N=1 - before gas-phase chemistry - // N=2 - before cloud chemistry - // N=3 - current incoming values (before gas-aerosol exchange, newnuc, - // coag) N=4 - updated outgoing values (after gas-aerosol exchange, - // newnuc, coag) - // - // qXXX_delaa are TMR changes (not tendencies) - // for different processes, which are used to produce history output - // for a clear sub-area, the processes are condensation/evaporation (and - // associated aging), - // renaming, coagulation, and nucleation - - // qxxx_del_yyyy are mix-ratio changes over full time step (deltat) - // qxxx_delsub_yyyy are mix-ratio changes over time sub-step (dtsubstep) - - static constexpr int num_gas_ids = AeroConfig::num_gas_ids(); - static constexpr int num_modes = AeroConfig::num_modes(); - static constexpr int num_aerosol_ids = AeroConfig::num_aerosol_ids(); - - static constexpr int igas_h2so4 = static_cast(GasId::H2SO4); - // Turn off nh3 for now. This is a future enhancement. - static constexpr int igas_nh3 = -999888777; // Same as mam_refactor - static constexpr int iaer_so4 = static_cast(AeroId::SO4); - static constexpr int iaer_pom = static_cast(AeroId::POM); - - const AeroId gas_to_aer[num_gas_ids] = {AeroId::SOA, AeroId::SO4, - AeroId::None}; - const bool l_gas_condense_to_mode[num_gas_ids][num_modes] = { - {true, true, true, true}, - {true, true, true, true}, - {false, false, false, false}}; - enum { NA, ANAL, IMPL }; - const int eqn_and_numerics_category[num_gas_ids] = {IMPL, ANAL, ANAL}; - // air molar density (kmol/m3) - // In order to try to match the results in mam_refactor - // set r_universal as [mJ/(mol)] as in mam_refactor. - // const Real r_universal = Constants::r_gas; // [mJ/(K mol)] - const Real r_universal = 8.314467591; // [mJ/(mol)] as in mam_refactor - const Real aircon = pmid / (1000 * r_universal * temp); - const Real alnsg_aer[num_modes] = {0.58778666490211906, 0.47000362924573563, - 0.58778666490211906, 0.47000362924573563}; - const Real uptk_rate_factor[num_gas_ids] = {0.81, 1.0, 1.0}; - - Real qgas_cur[num_gas_ids]; - for (int i = 0; i < num_gas_ids; ++i) - qgas_cur[i] = qgas3[i]; - Real qaer_cur[num_aerosol_ids][num_modes]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_cur[i][j] = qaer3[i][j]; - - Real qnum_cur[num_modes]; - for (int j = 0; j < num_modes; ++j) - qnum_cur[j] = qnum3[j]; - Real qwtr_cur[num_modes]; - for (int j = 0; j < num_modes; ++j) - qwtr_cur[j] = qwtr3[j]; - - Real qnumcw_cur[num_modes]; - for (int j = 0; j < num_modes; ++j) - qnumcw_cur[j] = qnumcw3[j]; - - Real qaercw_cur[num_gas_ids][num_modes]; - for (int i = 0; i < num_gas_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw_cur[i][j] = qaercw3[i][j]; - - Real qgas_netprod_otrproc[num_gas_ids] = {}; - if (config.do_cond && config.gaexch_h2so4_uptake_optaa == 2) { - for (int igas = 0; igas < num_gas_ids; ++igas) { - if (igas == igas_h2so4 || igas == igas_nh3) { - // if gaexch_h2so4_uptake_optaa == 2, then - // if qgas increases from pre-gaschem to post-cldchem, - // start from the pre-gaschem mix-ratio and add in the production - // during the integration - // if it decreases, - // start from post-cldchem mix-ratio - // *** currently just do this for h2so4 and nh3 - qgas_netprod_otrproc[igas] = (qgas3[igas] - qgas1[igas]) / deltat; - if (qgas_netprod_otrproc[igas] >= 0.0) - qgas_cur[igas] = qgas1[igas]; - else - qgas_netprod_otrproc[igas] = 0.0; - } - } - } - Real qgas_del_cond[num_gas_ids] = {}; - Real qgas_del_nnuc[num_gas_ids] = {}; - Real qgas_del_cond_only[num_gas_ids] = {}; - Real qaer_del_cond[num_aerosol_ids][num_modes] = {}; - Real qaer_del_rnam[num_aerosol_ids][num_modes] = {}; - Real qaer_del_nnuc[num_aerosol_ids][num_modes] = {}; - Real qaer_del_coag[num_aerosol_ids][num_modes] = {}; - Real qaer_delsub_cond[num_aerosol_ids][num_modes] = {}; - Real qaer_delsub_coag[num_aerosol_ids][num_modes] = {}; - Real qaer_del_cond_only[num_aerosol_ids][num_modes] = {}; - Real qaercw_del_rnam[num_aerosol_ids][num_modes] = {}; - Real qnum_del_cond[num_modes] = {}; - Real qnum_del_rnam[num_modes] = {}; - Real qnum_del_nnuc[num_modes] = {}; - Real qnum_del_coag[num_modes] = {}; - Real qnum_delsub_cond[num_modes] = {}; - Real qnum_delsub_coag[num_modes] = {}; - Real qnum_del_cond_only[num_modes] = {}; - Real qnumcw_del_rnam[num_modes] = {}; - Real qaer_delsub_coag_in[num_aerosol_ids][AeroConfig::max_agepair()] = {}; - - const int ntsubstep = 1; - Real dtsubstep = deltat; - if (ntsubstep > 1) - dtsubstep = deltat / ntsubstep; - - // loop over multiple time sub-steps - - for (int jtsubstep = 1; jtsubstep <= ntsubstep; ++jtsubstep) { - // gas-aerosol exchange - Real uptkrate_h2so4 = 0.0; - Real qgas_avg[num_gas_ids] = {}; - Real qgas_sv1[num_gas_ids] = {}; - Real qnum_sv1[num_modes] = {}; - Real qaer_sv1[num_aerosol_ids][num_modes] = {}; - Real qaer_delsub_grow4rnam[num_aerosol_ids][num_modes] = {}; - - if (config.do_cond) { - - const bool l_calc_gas_uptake_coeff = jtsubstep == 1; - Real uptkaer[num_gas_ids][num_modes] = {}; - - for (int i = 0; i < num_gas_ids; ++i) - qgas_sv1[i] = qgas_cur[i]; - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_sv1[j][i] = qaer_cur[j][i]; - - const int nghq = 2; - const int ntot_soamode = 4; - int niter_out = 0; - Real g0_soa_out = 0; - // time sub-step - const Real dtsub_soa_fixed = -1.0; - gasaerexch::mam_gasaerexch_1subarea( - nghq, igas_h2so4, igas_nh3, ntot_soamode, gas_to_aer, iaer_so4, - iaer_pom, l_calc_gas_uptake_coeff, l_gas_condense_to_mode, - eqn_and_numerics_category, dtsubstep, dtsub_soa_fixed, temp, pmid, - aircon, num_gas_ids, qgas_cur, qgas_avg, qgas_netprod_otrproc, - qaer_cur, qnum_cur, dgn_awet, alnsg_aer, uptk_rate_factor, uptkaer, - uptkrate_h2so4, niter_out, g0_soa_out); - - if (config.newnuc_h2so4_conc_optaa == 11) - qgas_avg[igas_h2so4] = - 0.5 * (qgas_sv1[igas_h2so4] + qgas_cur[igas_h2so4]); - else if (config.newnuc_h2so4_conc_optaa == 12) - qgas_avg[igas_h2so4] = qgas_cur[igas_h2so4]; - - for (int i = 0; i < num_gas_ids; ++i) - qgas_del_cond[i] += - (qgas_cur[i] - (qgas_sv1[i] + qgas_netprod_otrproc[i] * dtsubstep)); - - for (int i = 0; i < num_modes; ++i) - qnum_delsub_cond[i] = qnum_cur[i] - qnum_sv1[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_delsub_cond[i][j] = qaer_cur[i][j] - qaer_sv1[i][j]; - - // qaer_del_grow4rnam = change in qaer_del_cond during latest condensation - // calculations - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_delsub_grow4rnam[i][j] = qaer_cur[i][j] - qaer_sv1[i][j]; - for (int i = 0; i < num_gas_ids; ++i) - qgas_del_cond_only[i] = qgas_del_cond[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_del_cond_only[i][j] = qaer_delsub_cond[i][j]; - for (int i = 0; i < num_modes; ++i) - qnum_del_cond_only[i] = qnum_delsub_cond[i]; - - } else { - for (int i = 0; i < num_gas_ids; ++i) - qgas_avg[i] = qgas_cur[i]; - } - // renaming after "continuous growth" - if (config.do_rename) { - constexpr int nmodes = AeroConfig::num_modes(); - constexpr int naerosol_species = AeroConfig::num_aerosol_ids(); - const Real smallest_dryvol_value = 1.0e-25; // BAD_CONSTANT - const int dest_mode_of_mode[nmodes] = {-1, 0, -1, -1}; - - Real qnumcw_cur[num_modes] = {}; - Real qaercw_cur[num_aerosol_ids][num_modes] = {}; - Real qaercw_delsub_grow4rnam[num_aerosol_ids][num_modes] = {}; - Real mean_std_dev[nmodes]; - Real fmode_dist_tail_fac[nmodes]; - Real v2n_lo_rlx[nmodes]; - Real v2n_hi_rlx[nmodes]; - Real ln_diameter_tail_fac[nmodes]; - int num_pairs = 0; - Real diameter_cutoff[nmodes]; - Real ln_dia_cutoff[nmodes]; - Real diameter_threshold[nmodes]; - Real mass_2_vol[naerosol_species] = {0.15, - 6.4971751412429377e-002, - 0.15, - 7.0588235294117650e-003, - 3.0789473684210526e-002, - 5.1923076923076926e-002, - 156.20986883198000}; - - rename::find_renaming_pairs(dest_mode_of_mode, // in - mean_std_dev, // out - fmode_dist_tail_fac, // out - v2n_lo_rlx, // out - v2n_hi_rlx, // out - ln_diameter_tail_fac, // out - num_pairs, // out - diameter_cutoff, // out - ln_dia_cutoff, diameter_threshold); - - for (int i = 0; i < num_modes; ++i) - qnum_sv1[i] = qnum_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_sv1[j][i] = qaer_cur[j][i]; - Real dgnum_amode[nmodes]; - for (int m = 0; m < nmodes; ++m) { - dgnum_amode[m] = modes(m).nom_diameter; - } - - // qaercw_delsub_grow4rnam = change in qaercw from cloud chemistry - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw_delsub_grow4rnam[i][j] = - (qaercw3[i][j] - qaercw2[i][j]) / ntsubstep; - Real qnumcw_sv1[num_modes]; - for (int i = 0; i < num_modes; ++i) - qnumcw_sv1[i] = qnumcw_cur[i]; - Real qaercw_sv1[num_aerosol_ids][num_modes]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw_sv1[i][j] = qaercw_cur[i][j]; - - { - Real qmol_i_cur[num_modes][num_aerosol_ids]; - Real qmol_i_del[num_modes][num_aerosol_ids]; - Real qmol_c_cur[num_modes][num_aerosol_ids]; - Real qmol_c_del[num_modes][num_aerosol_ids]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) { - qmol_i_cur[i][j] = qaer_cur[j][i]; - qmol_i_del[i][j] = qaer_delsub_grow4rnam[j][i]; - qmol_c_cur[i][j] = qaercw_cur[j][i]; - qmol_c_del[i][j] = qaercw_delsub_grow4rnam[j][i]; - } - - Rename rename; - rename.mam_rename_1subarea_( - iscldy_subarea, smallest_dryvol_value, dest_mode_of_mode, - mean_std_dev, fmode_dist_tail_fac, v2n_lo_rlx, v2n_hi_rlx, - ln_diameter_tail_fac, num_pairs, diameter_cutoff, ln_dia_cutoff, - diameter_threshold, mass_2_vol, dgnum_amode, qnum_cur, qmol_i_cur, - qmol_i_del, qnumcw_cur, qmol_c_cur, qmol_c_del); - - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) { - qaer_cur[j][i] = qmol_i_cur[i][j]; - qaer_delsub_grow4rnam[j][i] = qmol_i_del[i][j]; - qaercw_cur[j][i] = qmol_c_cur[i][j]; - qaercw_delsub_grow4rnam[j][i] = qmol_c_del[i][j]; - } - } - for (int i = 0; i < num_modes; ++i) - qnum_del_rnam[i] += qnum_cur[i] - qnum_sv1[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaer_del_rnam[i][j] += qaer_cur[i][j] - qaer_sv1[i][j]; - for (int i = 0; i < num_modes; ++i) - qnumcw_del_rnam[i] += qnumcw_cur[i] - qnumcw_sv1[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw_del_rnam[i][j] += qaercw_cur[i][j] - qaercw_sv1[i][j]; - } - - // primary carbon aging - if (config.do_cond) { - aging::mam_pcarbon_aging_1subarea( - dgn_a, qnum_cur, qnum_delsub_cond, qnum_delsub_coag, qaer_cur, - qaer_delsub_cond, qaer_delsub_coag, qaer_delsub_coag_in); - } - // accumulate sub-step q-dels - if (config.do_cond) { - for (int i = 0; i < num_modes; ++i) - qnum_del_cond[i] += qnum_delsub_cond[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer_del_cond[j][i] += qaer_delsub_cond[j][i]; - } - } - // final mix ratios - for (int i = 0; i < num_gas_ids; ++i) - qgas4[i] = qgas_cur[i]; - for (int j = 0; j < num_aerosol_ids; ++j) - for (int i = 0; i < num_modes; ++i) - qaer4[j][i] = qaer_cur[j][i]; - for (int i = 0; i < num_modes; ++i) - qnum4[i] = qnum_cur[i]; - for (int i = 0; i < num_modes; ++i) - qwtr4[i] = qwtr_cur[i]; - for (int i = 0; i < num_modes; ++i) - qnumcw4[i] = qnumcw_cur[i]; - for (int i = 0; i < num_gas_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw4[i][j] = qaercw_cur[i][j]; - - // final mix ratio changes - for (int i = 0; i < num_gas_ids; ++i) { - qgas_delaa[i][iqtend_cond()] = qgas_del_cond[i]; - qgas_delaa[i][iqtend_rnam()] = 0.0; - qgas_delaa[i][iqtend_nnuc()] = qgas_del_nnuc[i]; - qgas_delaa[i][iqtend_coag()] = 0.0; - qgas_delaa[i][iqtend_cond_only()] = qgas_del_cond_only[i]; - } - for (int i = 0; i < num_modes; ++i) { - qnum_delaa[i][iqtend_cond()] = qnum_del_cond[i]; - qnum_delaa[i][iqtend_rnam()] = qnum_del_rnam[i]; - qnum_delaa[i][iqtend_nnuc()] = qnum_del_nnuc[i]; - qnum_delaa[i][iqtend_coag()] = qnum_del_coag[i]; - qnum_delaa[i][iqtend_cond_only()] = qnum_del_cond_only[i]; - } - for (int j = 0; j < num_aerosol_ids; ++j) { - for (int i = 0; i < num_modes; ++i) { - qaer_delaa[j][i][iqtend_cond()] = qaer_del_cond[j][i]; - qaer_delaa[j][i][iqtend_rnam()] = qaer_del_rnam[j][i]; - qaer_delaa[j][i][iqtend_nnuc()] = qaer_del_nnuc[j][i]; - qaer_delaa[j][i][iqtend_coag()] = qaer_del_coag[j][i]; - qaer_delaa[j][i][iqtend_cond_only()] = qaer_del_cond_only[j][i]; - } - } - for (int i = 0; i < num_modes; ++i) - qnumcw_delaa[i][iqqcwtend_rnam()] = qnumcw_del_rnam[i]; - for (int i = 0; i < num_aerosol_ids; ++i) - for (int j = 0; j < num_modes; ++j) - qaercw_delaa[i][j][iqqcwtend_rnam()] = qaercw_del_rnam[i][j]; -} - -KOKKOS_INLINE_FUNCTION -void mam_amicphys_1gridcell( - const AmicPhysConfig& config, const int nstep, const Real deltat, const int nsubarea, - const int ncldy_subarea, const bool iscldy_subarea[maxsubarea()], - const Real afracsub[maxsubarea()], const Real temp, const Real pmid, - const Real pdel, const Real zmid, const Real pblh, - const Real relhumsub[maxsubarea()], Real dgn_a[AeroConfig::num_modes()], - Real dgn_awet[AeroConfig::num_modes()], - Real wetdens[AeroConfig::num_modes()], - const Real qsub1[AeroConfig::num_gas_ids()][maxsubarea()], - const Real qsub2[AeroConfig::num_gas_ids()][maxsubarea()], - const Real qqcwsub2[AeroConfig::num_gas_ids()][maxsubarea()], - const Real qsub3[AeroConfig::num_gas_ids()][maxsubarea()], - const Real qqcwsub3[AeroConfig::num_gas_ids()][maxsubarea()], - Real qaerwatsub3[AeroConfig::num_modes()][maxsubarea()], - Real qsub4[AeroConfig::num_gas_ids()][maxsubarea()], - Real qqcwsub4[AeroConfig::num_gas_ids()][maxsubarea()], - Real qaerwatsub4[AeroConfig::num_modes()][maxsubarea()], - Real qsub_tendaa[AeroConfig::num_gas_ids()][nqtendaa()][maxsubarea()], - Real qqcwsub_tendaa[AeroConfig::num_gas_ids()][nqqcwtendaa()][maxsubarea()]) { - - // - // calculates changes to gas and aerosol sub-area TMRs (tracer mixing ratios) - // qsub3 and qqcwsub3 are the incoming current TMRs - // qsub4 and qqcwsub4 are the outgoing updated TMRs - // - // qsubN and qqcwsubN (N=1:4) are tracer mixing ratios (TMRs, mol/mol or - // #/kmol) in sub-areas - // currently there are just clear and cloudy sub-areas - // the N=1:4 have same meanings as for qgcmN - // N=1 - before gas-phase chemistry - // N=2 - before cloud chemistry - // N=3 - incoming values (before gas-aerosol exchange, newnuc, coag) - // N=4 - outgoing values (after gas-aerosol exchange, newnuc, coag) - // qsub_tendaa and qqcwsub_tendaa are TMR tendencies - // for different processes, which are used to produce history output - // the processes are condensation/evaporation (and associated aging), - // renaming, coagulation, and nucleation - - static constexpr int num_gas_ids = AeroConfig::num_gas_ids(); - static constexpr int num_modes = AeroConfig::num_modes(); - static constexpr int num_aerosol_ids = AeroConfig::num_aerosol_ids(); - - // the q--4 values will be equal to q--3 values unless they get changed - for (int i = 0; i < num_gas_ids; ++i) - for (int j = 0; j < maxsubarea(); ++j) { - qsub4[i][j] = qsub3[i][j]; - qqcwsub4[i][j] = qqcwsub3[i][j]; - } - for (int i = 0; i < num_modes; ++i) - for (int j = 0; j < maxsubarea(); ++j) - qaerwatsub4[i][j] = qaerwatsub3[i][j]; - for (int i = 0; i < num_gas_ids; ++i) - for (int j = 0; j < nqtendaa(); ++j) - for (int k = 0; k < maxsubarea(); ++k) - qsub_tendaa[i][j][k] = 0; - for (int i = 0; i < num_gas_ids; ++i) - for (int j = 0; j < nqqcwtendaa(); ++j) - for (int k = 0; k < maxsubarea(); ++k) - qqcwsub_tendaa[i][j][k] = 0.0; - - for (int jsub = 0; jsub < nsubarea; ++jsub) { - AmicPhysConfig sub_config = config; - if (iscldy_subarea[jsub]) { - sub_config.do_cond = config.do_cond; - sub_config.do_rename = config.do_rename; - sub_config.do_newnuc = false; - sub_config.do_coag = false; - } - const bool do_map_gas_sub = sub_config.do_cond || sub_config.do_newnuc; - - // map incoming sub-area mix-ratios to gas/aer/num arrays - Real qgas1[num_gas_ids] = {}; - Real qgas3[num_gas_ids] = {}; - Real qgas4[num_gas_ids] = {}; - if (do_map_gas_sub) { - // for cldy subarea, only do gases if doing gaexch - for (int igas = 0; igas < 2; ++igas) { - const int l = lmap_gas(igas); - qgas1[igas] = qsub1[l][jsub] * fcvt_gas(igas); - qgas3[igas] = qsub3[l][jsub] * fcvt_gas(igas); - qgas4[igas] = qgas3[igas]; - } - } - Real qaer2[num_aerosol_ids][num_modes] = {}; - Real qaer3[num_aerosol_ids][num_modes] = {}; - Real qnum3[num_modes] = {}; - Real qaer4[num_aerosol_ids][num_modes] = {}; - Real qnum4[num_modes] = {}; - Real qwtr3[num_modes] = {}; - Real qwtr4[num_modes] = {}; - for (int n = 0; n < num_modes; ++n) { - qnum3[n] = qsub3[n][jsub] * fcvt_num(); - qnum4[n] = qnum3[n]; - for (int iaer = 0; iaer < num_aerosol_ids; ++iaer) { - qaer2[iaer][n] = qsub2[iaer][jsub] * fcvt_aer(iaer); - qaer3[iaer][n] = qsub3[iaer][jsub] * fcvt_aer(iaer); - qaer4[iaer][n] = qaer3[iaer][n]; - } - qwtr3[n] = qaerwatsub3[n][jsub] * fcvt_wtr(); - qwtr4[n] = qwtr3[n]; - } - Real qaercw2[num_aerosol_ids][num_modes] = {}; - Real qaercw3[num_aerosol_ids][num_modes] = {}; - Real qnumcw3[num_modes] = {}; - Real qaercw4[num_aerosol_ids][num_modes] = {}; - Real qnumcw4[num_modes] = {}; - if (iscldy_subarea[jsub]) { - // only do cloud-borne for cloudy - for (int n = 0; n < num_modes; ++n) { - qnumcw3[n] = qqcwsub3[n][jsub] * fcvt_num(); - qnumcw4[n] = qnumcw3[n]; - for (int iaer = 0; iaer < num_aerosol_ids; ++iaer) { - qaercw2[iaer][n] = qqcwsub2[n][jsub] * fcvt_aer(iaer); - qaercw3[iaer][n] = qqcwsub3[n][jsub] * fcvt_aer(iaer); - qaercw4[iaer][n] = qaercw3[iaer][n]; - } - } - } - - Real qgas_delaa[num_gas_ids][nqtendaa()] = {}; - Real qnum_delaa[num_modes][nqtendaa()] = {}; - Real qnumcw_delaa[num_modes][nqqcwtendaa()] = {}; - Real qaer_delaa[num_aerosol_ids][num_modes][nqtendaa()] = {}; - Real qaercw_delaa[num_aerosol_ids][num_modes][nqqcwtendaa()] = {}; - - if (iscldy_subarea[jsub]) { - mam_amicphys_1subarea_cloudy(sub_config, nstep, deltat, - jsub, nsubarea, iscldy_subarea[jsub], afracsub[jsub], temp, pmid, - pdel, zmid, pblh, relhumsub[jsub], dgn_a, dgn_awet, wetdens, qgas1, - qgas3, qgas4, qgas_delaa, qnum3, qnum4, qnum_delaa, qaer2, qaer3, - qaer4, qaer_delaa, qwtr3, qwtr4, qnumcw3, qnumcw4, qnumcw_delaa, - qaercw2, qaercw3, qaercw4, qaercw_delaa); - } else { - mam_amicphys_1subarea_clear(sub_config, nstep, deltat, - jsub, nsubarea, iscldy_subarea[jsub], afracsub[jsub], temp, pmid, - pdel, zmid, pblh, relhumsub[jsub], dgn_a, dgn_awet, wetdens, qgas1, - qgas3, qgas4, qgas_delaa, qnum3, qnum4, qnum_delaa, qaer3, qaer4, - qaer_delaa, qwtr3, qwtr4); - // map gas/aer/num arrays (mix-ratio and del=change) back to sub-area - // arrays - - if (do_map_gas_sub) { - for (int igas = 0; igas < 2; ++igas) { - const int l = lmap_gas(igas); - qsub4[l][jsub] = qgas4[igas] / fcvt_gas(igas); - for (int i = 0; i < nqtendaa(); ++i) - qsub_tendaa[l][i][jsub] = - qgas_delaa[igas][i] / (fcvt_gas(igas) * deltat); - } - } - for (int n = 0; n < num_modes; ++n) { - qsub4[n][jsub] = qnum4[n] / fcvt_num(); - for (int i = 0; i < nqtendaa(); ++i) - qsub_tendaa[n][i][jsub] = qnum_delaa[n][i] / (fcvt_num() * deltat); - for (int iaer = 0; iaer < num_aerosol_ids; ++iaer) { - qsub4[iaer][jsub] = qaer4[iaer][n] / fcvt_aer(iaer); - for (int i = 0; i < nqtendaa(); ++i) - qsub_tendaa[iaer][i][jsub] = - qaer_delaa[iaer][n][i] / (fcvt_aer(iaer) * deltat); - } - qaerwatsub4[n][jsub] = qwtr4[n] / fcvt_wtr(); - - if (iscldy_subarea[jsub]) { - qqcwsub4[n][jsub] = qnumcw4[n] / fcvt_num(); - for (int i = 0; i < nqqcwtendaa(); ++i) - qqcwsub_tendaa[n][i][jsub] = - qnumcw_delaa[n][i] / (fcvt_num() * deltat); - for (int iaer = 0; iaer < num_aerosol_ids; ++iaer) { - qqcwsub4[iaer][jsub] = qaercw4[iaer][n] / fcvt_aer(iaer); - for (int i = 0; i < nqqcwtendaa(); ++i) - qqcwsub_tendaa[iaer][i][jsub] = - qaercw_delaa[iaer][n][i] / (fcvt_aer(iaer) * deltat); - } - } - } - } - } -} - -} // anonymous namespace - -KOKKOS_INLINE_FUNCTION -void modal_aero_amicphys_intr( - const AmicPhysConfig& config, const int nstep, const Real deltat, const Real t, const Real pmid, const Real pdel, - const Real zm, const Real pblh, const Real qv, const Real cld, - Real q[gas_pcnst()], Real qqcw[gas_pcnst()], const Real q_pregaschem[gas_pcnst()], - const Real q_precldchem[gas_pcnst()], const Real qqcw_precldchem[gas_pcnst()], - Real q_tendbb[gas_pcnst()][nqtendbb()], Real qqcw_tendbb[gas_pcnst()][nqtendbb()], - Real dgncur_a[AeroConfig::num_modes()], - Real dgncur_awet[AeroConfig::num_modes()], - Real wetdens_host[AeroConfig::num_modes()], - Real qaerwat[AeroConfig::num_modes()]) { - - /* - nstep ! model time-step number - nqtendbb ! dimension for q_tendbb - nqqcwtendbb ! dimension f - deltat ! - q(ncol,pver,pcnstxx) ! current tracer mixing ratios (TMRs) - these values are updated (so out /= in) - *** MUST BE #/kmol-air for number - *** MUST BE mol/mol-air for mass - *** NOTE ncol dimension - qqcw(ncol,pver,pcnstxx) - like q but for cloud-borner tracers - these values are updated - q_pregaschem(ncol,pver,pcnstxx) ! q TMRs before gas-phase - chemistry q_precldchem(ncol,pver,pcnstxx) ! q TMRs before cloud - chemistry qqcw_precldchem(ncol,pver,pcnstxx) ! qqcw TMRs before cloud - chemistry q_tendbb(ncol,pver,pcnstxx,nqtendbb()) ! TMR tendencies for - box-model diagnostic output qqcw_tendbb(ncol,pver,pcnstx t(pcols,pver) ! - temperature at model levels (K) pmid(pcols,pver) ! pressure at model - level centers (Pa) pdel(pcols,pver) ! pressure thickness of levels - (Pa) zm(pcols,pver) ! altitude (above ground) at level centers (m) - pblh(pcols) ! planetary boundary layer depth (m) - qv(pcols,pver) ! specific humidity (kg/kg) - cld(ncol,pver) ! cloud fraction (-) *** NOTE ncol dimension - dgncur_a(pcols,pver,ntot_amode) - dgncur_awet(pcols,pver,ntot_amode) - ! dry & wet geo. mean dia. (m) of - number distrib. wetdens_host(pcols,pver,ntot_amode) ! interstitial - aerosol wet density (kg/m3) - - qaerwat(pcols,pver,ntot_amode aerosol water mixing ratio (kg/kg, - NOT mol/mol) - - */ - - // !DESCRIPTION: - // calculates changes to gas and aerosol TMRs (tracer mixing ratios) from - // gas-aerosol exchange (condensation/evaporation) - // growth from smaller to larger modes (renaming) due to both - // condensation and cloud chemistry - // new particle nucleation - // coagulation - // transfer of particles from hydrophobic modes to hydrophilic modes - // (aging) - // due to condensation and coagulation - // - // the incoming mixing ratios (q and qqcw) are updated before output - // - // !REVISION HISTORY: - // RCE 07.04.13: Adapted from earlier version of CAM5 modal aerosol - // routines - // for these processes - // - - static constexpr int num_modes = AeroConfig::num_modes(); - - // qgcmN and qqcwgcmN (N=1:4) are grid-cell mean tracer mixing ratios - // (TMRs, mol/mol or #/kmol) - // N=1 - before gas-phase chemistry - // N=2 - before cloud chemistry - // N=3 - incoming values (before gas-aerosol exchange, newnuc, coag) - // N=4 - outgoing values (after gas-aerosol exchange, newnuc, coag) - - // qsubN and qqcwsubN (N=1:4) are TMRs in sub-areas - // currently there are just clear and cloudy sub-areas - // the N=1:4 have same meanings as for qgcmN - - // q_coltendaa and qqcw_coltendaa are column-integrated tendencies - // for different processes, which are output to history - // the processes are condensation/evaporation (and associated aging), - // renaming, coagulation, and nucleation - - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqtendbb(); ++j) - q_tendbb[i][j] = 0.0, qqcw_tendbb[i][j] = 0.0; - - // get saturation mixing ratio - // call qsat( t(1:ncol,1:pver), pmid(1:ncol,1:pvnner), & - // ev_sat(1:ncol,1:pver), qv_sat(1:ncol,1:pver) ) - const Real epsqs = haero::Constants::weight_ratio_h2o_air; - // Saturation vapor pressure - const Real ev_sat = conversions::vapor_saturation_pressure_magnus(t, pmid); - // Saturation specific humidity - const Real qv_sat = epsqs * ev_sat / (pmid - (1 - epsqs) * ev_sat); - - const Real relhumgcm = haero::max(0.0, haero::min(1.0, qv / qv_sat)); - - // Set up cloudy/clear subareas inside a grid cell - int nsubarea, ncldy_subarea, jclea, jcldy; - bool iscldy_subarea[maxsubarea()]; - Real afracsub[maxsubarea()]; - Real relhumsub[maxsubarea()]; - Real qsub1[gas_pcnst()][maxsubarea()]; - Real qsub2[gas_pcnst()][maxsubarea()]; - Real qsub3[gas_pcnst()][maxsubarea()]; - Real qqcwsub1[gas_pcnst()][maxsubarea()]; - Real qqcwsub2[gas_pcnst()][maxsubarea()]; - Real qqcwsub3[gas_pcnst()][maxsubarea()]; - // aerosol water mixing ratios (mol/mol) - Real qaerwatsub3[AeroConfig::num_modes()][maxsubarea()]; - construct_subareas_1gridcell(cld, relhumgcm, // in - q_pregaschem, q_precldchem, // in - qqcw_precldchem, // in - q, qqcw, // in - nsubarea, ncldy_subarea, jclea, jcldy, // out - iscldy_subarea, afracsub, relhumsub, // out - qsub1, qsub2, qsub3, // out - qqcwsub1, qqcwsub2, qqcwsub3, qaerwatsub3, // out - qaerwat // in - ); - - // Initialize the "after-amicphys" values - Real qsub4[gas_pcnst()][maxsubarea()] = {}; - Real qqcwsub4[gas_pcnst()][maxsubarea()] = {}; - Real qaerwatsub4[AeroConfig::num_modes()][maxsubarea()] = {}; - - // - // start integration - // - Real dgn_a[num_modes], dgn_awet[num_modes], wetdens[num_modes]; - for (int n = 0; n < num_modes; ++n) { - dgn_a[n] = dgncur_a[n]; - dgn_awet[n] = dgncur_awet[n]; - wetdens[n] = haero::max(1000.0, wetdens_host[n]); - } - Real qsub_tendaa[gas_pcnst()][nqtendaa()][maxsubarea()] = {}; - Real qqcwsub_tendaa[gas_pcnst()][nqqcwtendaa()][maxsubarea()] = {}; - mam_amicphys_1gridcell(config, nstep, deltat, - nsubarea, ncldy_subarea, iscldy_subarea, afracsub, t, - pmid, pdel, zm, pblh, relhumsub, dgn_a, dgn_awet, - wetdens, qsub1, qsub2, qqcwsub2, qsub3, qqcwsub3, - qaerwatsub3, qsub4, qqcwsub4, qaerwatsub4, qsub_tendaa, - qqcwsub_tendaa); - - // - // form new grid-mean mix-ratios - Real qgcm4[gas_pcnst()]; - Real qgcm_tendaa[gas_pcnst()][nqtendaa()]; - Real qaerwatgcm4[num_modes]; - if (nsubarea == 1) { - for (int i = 0; i < gas_pcnst(); ++i) - qgcm4[i] = qsub4[i][0]; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqtendaa(); ++j) - qgcm_tendaa[i][j] = qsub_tendaa[i][j][0]; - for (int i = 0; i < num_modes; ++i) - qaerwatgcm4[i] = qaerwatsub4[i][0]; - } else { - for (int i = 0; i < gas_pcnst(); ++i) - qgcm4[i] = 0.0; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqtendaa(); ++j) - qgcm_tendaa[i][j] = 0.0; - for (int n = 0; n < nsubarea; ++n) { - for (int i = 0; i < gas_pcnst(); ++i) - qgcm4[i] += qsub4[i][n] * afracsub[n]; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqtendaa(); ++j) - qgcm_tendaa[i][j] = - qgcm_tendaa[i][j] + qsub_tendaa[i][j][n] * afracsub[n]; - } - for (int i = 0; i < num_modes; ++i) - // for aerosol water use the clear sub-area value - qaerwatgcm4[i] = qaerwatsub4[i][jclea - 1]; - } - Real qqcwgcm4[gas_pcnst()]; - Real qqcwgcm_tendaa[gas_pcnst()][nqqcwtendaa()]; - if (ncldy_subarea <= 0) { - for (int i = 0; i < gas_pcnst(); ++i) - qqcwgcm4[i] = haero::max(0.0, qqcw[i]); - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqqcwtendaa(); ++j) - qqcwgcm_tendaa[i][j] = 0.0; - } else if (nsubarea == 1) { - for (int i = 0; i < gas_pcnst(); ++i) - qqcwgcm4[i] = qqcwsub4[i][0]; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqqcwtendaa(); ++j) - qqcwgcm_tendaa[i][j] = qqcwsub_tendaa[i][j][0]; - } else { - for (int i = 0; i < gas_pcnst(); ++i) - qqcwgcm4[i] = 0.0; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqqcwtendaa(); ++j) - qqcwgcm_tendaa[i][j] = 0.0; - for (int n = 0; n < nsubarea; ++n) { - if (iscldy_subarea[n]) { - for (int i = 0; i < gas_pcnst(); ++i) - qqcwgcm4[i] += qqcwsub4[i][n] * afracsub[n]; - for (int i = 0; i < gas_pcnst(); ++i) - for (int j = 0; j < nqqcwtendaa(); ++j) - qqcwgcm_tendaa[i][j] += qqcwsub_tendaa[i][j][n] * afracsub[n]; - } - } - } - - for (int lmz = 0; lmz < gas_pcnst(); ++lmz) { - if (lmapcc_all(lmz) > 0) { - // HW, to ensure non-negative - q[lmz] = haero::max(qgcm4[lmz], 0.0); - if (lmapcc_all(lmz) >= lmapcc_val_aer()) { - // HW, to ensure non-negative - qqcw[lmz] = haero::max(qqcwgcm4[lmz], 0.0); - } - } - } - for (int i = 0; i < gas_pcnst(); ++i) { - if (iqtend_cond() < nqtendbb()) - q_tendbb[i][iqtend_cond()] = qgcm_tendaa[i][iqtend_cond()]; - if (iqtend_rnam() < nqtendbb()) - q_tendbb[i][iqtend_rnam()] = qgcm_tendaa[i][iqtend_rnam()]; - if (iqtend_nnuc() < nqtendbb()) - q_tendbb[i][iqtend_nnuc()] = qgcm_tendaa[i][iqtend_nnuc()]; - if (iqtend_coag() < nqtendbb()) - q_tendbb[i][iqtend_coag()] = qgcm_tendaa[i][iqtend_coag()]; - if (iqqcwtend_rnam() < nqqcwtendbb()) - qqcw_tendbb[i][iqqcwtend_rnam()] = qqcwgcm_tendaa[i][iqqcwtend_rnam()]; - } - for (int i = 0; i < num_modes; ++i) - qaerwat[i] = qaerwatgcm4[i]; -} - -} // namespace scream::impl diff --git a/components/eamxx/src/physics/mam/mam_coupling.hpp b/components/eamxx/src/physics/mam/mam_coupling.hpp index a3a5f746aa3..68f67f62c5c 100644 --- a/components/eamxx/src/physics/mam/mam_coupling.hpp +++ b/components/eamxx/src/physics/mam/mam_coupling.hpp @@ -647,7 +647,7 @@ void compute_vertical_layer_heights(const Team& team, const auto p_mid = ekat::subview(dry_atm.p_mid, column_index); const auto T_mid = ekat::subview(dry_atm.T_mid, column_index); const auto pseudo_density = ekat::subview(dry_atm.p_del, column_index); - + // NOTE: we are using dry qv. Does calculate_dz require dry or wet? PF::calculate_dz(team, pseudo_density, p_mid, T_mid, qv, // inputs dz);//output @@ -811,196 +811,6 @@ void copy_view_lev_slice(haero::ThreadTeamPolicy team_policy, //inputs }); } -// Because CUDA C++ doesn't allow us to declare and use constants outside of -// KOKKOS_INLINE_FUNCTIONS, we define this macro that allows us to (re)define -// these constants where needed within two such functions so we don't define -// them inconsistently. Yes, it's the 21st century and we're still struggling -// with these basic things. -#define DECLARE_PROG_TRANSFER_CONSTANTS \ - /* mapping of constituent indices to aerosol modes */ \ - const auto Accum = mam4::ModeIndex::Accumulation; \ - const auto Aitken = mam4::ModeIndex::Aitken; \ - const auto Coarse = mam4::ModeIndex::Coarse; \ - const auto PC = mam4::ModeIndex::PrimaryCarbon; \ - const auto NoMode = mam4::ModeIndex::None; \ - static const mam4::ModeIndex mode_for_cnst[gas_pcnst()] = { \ - NoMode, NoMode, NoMode, NoMode, NoMode, NoMode, /* gases (not aerosols) */ \ - Accum, Accum, Accum, Accum, Accum, Accum, Accum, Accum, /* 7 aero species + NMR */ \ - Aitken, Aitken, Aitken, Aitken, Aitken, /* 4 aero species + NMR */ \ - Coarse, Coarse, Coarse, Coarse, Coarse, Coarse, Coarse, Coarse, /* 7 aero species + NMR */ \ - PC, PC, PC, PC, /* 3 aero species + NMR */ \ - }; \ - /* mapping of constituent indices to aerosol species */ \ - const auto SOA = mam4::AeroId::SOA; \ - const auto SO4 = mam4::AeroId::SO4; \ - const auto POM = mam4::AeroId::POM; \ - const auto BC = mam4::AeroId::BC; \ - const auto NaCl = mam4::AeroId::NaCl; \ - const auto DST = mam4::AeroId::DST; \ - const auto MOM = mam4::AeroId::MOM; \ - const auto NoAero = mam4::AeroId::None; \ - static const mam4::AeroId aero_for_cnst[gas_pcnst()] = { \ - NoAero, NoAero, NoAero, NoAero, NoAero, NoAero, /* gases (not aerosols) */ \ - SO4, POM, SOA, BC, DST, NaCl, MOM, NoAero, /* accumulation mode */ \ - SO4, SOA, NaCl, MOM, NoAero, /* aitken mode */ \ - DST, NaCl, SO4, BC, POM, SOA, MOM, NoAero, /* coarse mode */ \ - POM, BC, MOM, NoAero, /* primary carbon mode */ \ - }; \ - /* mapping of constituent indices to gases */ \ - const auto O3 = mam4::GasId::O3; \ - const auto H2O2 = mam4::GasId::H2O2; \ - const auto H2SO4 = mam4::GasId::H2SO4; \ - const auto SO2 = mam4::GasId::SO2; \ - const auto DMS = mam4::GasId::DMS; \ - const auto SOAG = mam4::GasId::SOAG; \ - const auto NoGas = mam4::GasId::None; \ - static const mam4::GasId gas_for_cnst[gas_pcnst()] = { \ - O3, H2O2, H2SO4, SO2, DMS, SOAG, \ - NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, \ - NoGas, NoGas, NoGas, NoGas, NoGas, \ - NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, NoGas, \ - NoGas, NoGas, NoGas, NoGas, \ - }; - -// Given a Prognostics object, transfers data for interstitial aerosols to the -// chemistry work array q, and cloudborne aerosols to the chemistry work array -// qqcw, both at vertical level k. The input and output quantities are stored as -// number/mass mixing ratios. -// NOTE: this mapping is chemistry-mechanism-specific (see mo_sim_dat.F90 -// NOTE: in the relevant preprocessed chemical mechanism) -// NOTE: see mam4xx/aero_modes.hpp to interpret these mode/aerosol/gas -// NOTE: indices -KOKKOS_INLINE_FUNCTION -void transfer_prognostics_to_work_arrays(const mam4::Prognostics &progs, - const int k, - Real q[gas_pcnst()], - Real qqcw[gas_pcnst()]) { - DECLARE_PROG_TRANSFER_CONSTANTS - - // copy number/mass mixing ratios from progs to q and qqcw at level k, - // converting them to VMR - for (int i = 0; i < gas_pcnst(); ++i) { - auto mode_index = mode_for_cnst[i]; - auto aero_id = aero_for_cnst[i]; - auto gas_id = gas_for_cnst[i]; - if (gas_id != NoGas) { // constituent is a gas - int g = static_cast(gas_id); - q[i] = progs.q_gas[g](k); - qqcw[i] = progs.q_gas[g](k); - } else { - int m = static_cast(mode_index); - if (aero_id != NoAero) { // constituent is an aerosol species - int a = aerosol_index_for_mode(mode_index, aero_id); - q[i] = progs.q_aero_i[m][a](k); - qqcw[i] = progs.q_aero_c[m][a](k); - } else { // constituent is a modal number mixing ratio - int m = static_cast(mode_index); - q[i] = progs.n_mode_i[m](k); - qqcw[i] = progs.n_mode_c[m](k); - } - } - } -} - -// converts the quantities in the work arrays q and qqcw from mass/number -// mixing ratios to volume/number mixing ratios -KOKKOS_INLINE_FUNCTION -void convert_work_arrays_to_vmr(const Real q[gas_pcnst()], - const Real qqcw[gas_pcnst()], - Real vmr[gas_pcnst()], - Real vmrcw[gas_pcnst()]) { - DECLARE_PROG_TRANSFER_CONSTANTS - - for (int i = 0; i < gas_pcnst(); ++i) { - auto mode_index = mode_for_cnst[i]; - auto aero_id = aero_for_cnst[i]; - auto gas_id = gas_for_cnst[i]; - if (gas_id != NoGas) { // constituent is a gas - int g = static_cast(gas_id); - const Real mw = mam4::gas_species(g).molecular_weight; - vmr[i] = mam4::conversions::vmr_from_mmr(q[i], mw); - vmrcw[i] = mam4::conversions::vmr_from_mmr(qqcw[i], mw); - } else { - if (aero_id != NoAero) { // constituent is an aerosol species - int a = aerosol_index_for_mode(mode_index, aero_id); - const Real mw = mam4::aero_species(a).molecular_weight; - vmr[i] = mam4::conversions::vmr_from_mmr(q[i], mw); - vmrcw[i] = mam4::conversions::vmr_from_mmr(qqcw[i], mw); - } else { // constituent is a modal number mixing ratio - vmr[i] = q[i]; - vmrcw[i] = qqcw[i]; - } - } - } -} - -// converts the quantities in the work arrays vmrs and vmrscw from mass/number -// mixing ratios to volume/number mixing ratios -KOKKOS_INLINE_FUNCTION -void convert_work_arrays_to_mmr(const Real vmr[gas_pcnst()], - const Real vmrcw[gas_pcnst()], - Real q[gas_pcnst()], - Real qqcw[gas_pcnst()]) { - DECLARE_PROG_TRANSFER_CONSTANTS - - for (int i = 0; i < gas_pcnst(); ++i) { - auto mode_index = mode_for_cnst[i]; - auto aero_id = aero_for_cnst[i]; - auto gas_id = gas_for_cnst[i]; - if (gas_id != NoGas) { // constituent is a gas - int g = static_cast(gas_id); - const Real mw = mam4::gas_species(g).molecular_weight; - q[i] = mam4::conversions::mmr_from_vmr(vmr[i], mw); - qqcw[i] = mam4::conversions::mmr_from_vmr(vmrcw[i], mw); - } else { - if (aero_id != NoAero) { // constituent is an aerosol species - int a = aerosol_index_for_mode(mode_index, aero_id); - const Real mw = mam4::aero_species(a).molecular_weight; - q[i] = mam4::conversions::mmr_from_vmr(vmr[i], mw); - qqcw[i] = mam4::conversions::mmr_from_vmr(vmrcw[i], mw); - } else { // constituent is a modal number mixing ratio - q[i] = vmr[i]; - qqcw[i] = vmrcw[i]; - } - } - } -} - -// Given work arrays with interstitial and cloudborne aerosol data, transfers -// them to the given Prognostics object at the kth vertical level. This is the -// "inverse operator" for transfer_prognostics_to_work_arrays, above. -KOKKOS_INLINE_FUNCTION -void transfer_work_arrays_to_prognostics(const Real q[gas_pcnst()], - const Real qqcw[gas_pcnst()], - mam4::Prognostics &progs, const int k) { - DECLARE_PROG_TRANSFER_CONSTANTS - - // copy number/mass mixing ratios from progs to q and qqcw at level k, - // converting them to VMR - for (int i = 0; i < gas_pcnst(); ++i) { - auto mode_index = mode_for_cnst[i]; - auto aero_id = aero_for_cnst[i]; - auto gas_id = gas_for_cnst[i]; - if (gas_id != NoGas) { // constituent is a gas - int g = static_cast(gas_id); - progs.q_gas[g](k) = q[i]; - } else { - int m = static_cast(mode_index); - if (aero_id != NoAero) { // constituent is an aerosol species - int a = aerosol_index_for_mode(mode_index, aero_id); - progs.q_aero_i[m][a](k) = q[i]; - progs.q_aero_c[m][a](k) = qqcw[i]; - } else { // constituent is a modal number mixing ratio - int m = static_cast(mode_index); - progs.n_mode_i[m](k) = q[i]; - progs.n_mode_c[m](k) = qqcw[i]; - } - } - } -} - -#undef DECLARE_PROG_TRANSFER_CONSTANTS - } // namespace scream::mam_coupling #endif diff --git a/components/eamxx/src/physics/mam/readfiles/photo_table_utils.cpp b/components/eamxx/src/physics/mam/readfiles/photo_table_utils.cpp new file mode 100644 index 00000000000..ab5a0509274 --- /dev/null +++ b/components/eamxx/src/physics/mam/readfiles/photo_table_utils.cpp @@ -0,0 +1,151 @@ +#include + +namespace scream::impl { + +// number of photolysis reactions +using mam4::mo_photo::phtcnt; + +using HostView1D = haero::DeviceType::view_1d::HostMirror; +using HostViewInt1D = haero::DeviceType::view_1d::HostMirror; + +//------------------------------------------------------------------------- +// Reading the photolysis table +//------------------------------------------------------------------------- + +// ON HOST (MPI root only), sets the lng_indexer and pht_alias_mult_1 host views +// according to parameters in our (hardwired) chemical mechanism + +std::vector populate_etfphot_from_e3sm_case() { + // We obtained these values from an e3sm simulations. + // We should only use this function on Host. + std::vector etfphot_data = { + 7.5691227E+11, 8.6525905E+11, 1.0355749E+12, 1.1846288E+12, 2.1524405E+12, + 3.2362584E+12, 3.7289849E+12, 4.4204330E+12, 4.6835350E+12, 6.1217728E+12, + 4.5575051E+12, 5.3491446E+12, 4.7016063E+12, 5.4281722E+12, 4.5023968E+12, + 6.8931981E+12, 6.2012647E+12, 6.1430771E+12, 5.7820385E+12, 7.6770646E+12, + 1.3966509E+13, 1.2105348E+13, 2.8588980E+13, 3.2160821E+13, 2.4978066E+13, + 2.7825401E+13, 2.3276451E+13, 3.6343684E+13, 6.1787886E+13, 7.8009914E+13, + 7.6440824E+13, 7.6291458E+13, 9.4645085E+13, 1.0124628E+14, 1.0354111E+14, + 1.0999650E+14, 1.0889946E+14, 1.1381912E+14, 1.3490042E+14, 1.5941519E+14, + 1.4983265E+14, 1.5184267E+14, 1.5991420E+14, 1.6976697E+14, 1.8771840E+14, + 1.6434367E+14, 1.8371960E+14, 2.1966369E+14, 1.9617879E+14, 2.2399700E+14, + 1.8429912E+14, 2.0129736E+14, 2.0541588E+14, 2.4334962E+14, 3.5077122E+14, + 3.4517894E+14, 3.5749668E+14, 3.6624304E+14, 3.4975113E+14, 3.5566025E+14, + 4.2825273E+14, 4.8406375E+14, 4.9511159E+14, 5.2695368E+14, 5.2401611E+14, + 5.0877746E+14, 4.8780853E+14}; + return etfphot_data; +} + +// This version uses scream_scorpio_interface to read netcdf files. +mam4::mo_photo::PhotoTableData read_photo_table( + const std::string &rsf_file, const std::string &xs_long_file) { + // set up the lng_indexer and pht_alias_mult_1 views based on our + // (hardwired) chemical mechanism + HostViewInt1D lng_indexer_h("lng_indexer", phtcnt); + + int nw, nump, numsza, numcolo3, numalb, nt, np_xs; // table dimensions + scorpio::register_file(rsf_file, scorpio::Read); + // read and broadcast dimension data + nump = scorpio::get_dimlen(rsf_file, "numz"); + numsza = scorpio::get_dimlen(rsf_file, "numsza"); + numalb = scorpio::get_dimlen(rsf_file, "numalb"); + numcolo3 = scorpio::get_dimlen(rsf_file, "numcolo3fact"); + + scorpio::register_file(xs_long_file, scorpio::Read); + nt = scorpio::get_dimlen(xs_long_file, "numtemp"); + nw = scorpio::get_dimlen(xs_long_file, "numwl"); + np_xs = scorpio::get_dimlen(xs_long_file, "numprs"); + + // FIXME: hard-coded for only one photo reaction. + std::string rxt_names[1] = {"jh2o2"}; + int numj = 1; + lng_indexer_h(0) = 0; + // allocate the photolysis table + auto table = mam4::mo_photo::create_photo_table_data( + nw, nt, np_xs, numj, nump, numsza, numcolo3, numalb); + + // allocate host views for table data + auto rsf_tab_h = Kokkos::create_mirror_view(table.rsf_tab); + auto xsqy_h = Kokkos::create_mirror_view(table.xsqy); + auto sza_h = Kokkos::create_mirror_view(table.sza); + auto alb_h = Kokkos::create_mirror_view(table.alb); + auto press_h = Kokkos::create_mirror_view(table.press); + auto colo3_h = Kokkos::create_mirror_view(table.colo3); + auto o3rat_h = Kokkos::create_mirror_view(table.o3rat); + // auto etfphot_h = Kokkos::create_mirror_view(table.etfphot); + auto prs_h = Kokkos::create_mirror_view(table.prs); + + // read file data into our host views + scorpio::read_var(rsf_file, "pm", press_h.data()); + scorpio::read_var(rsf_file, "sza", sza_h.data()); + scorpio::read_var(rsf_file, "alb", alb_h.data()); + scorpio::read_var(rsf_file, "colo3fact", o3rat_h.data()); + scorpio::read_var(rsf_file, "colo3", colo3_h.data()); + // it produces an error. + scorpio::read_var(rsf_file, "RSF", rsf_tab_h.data()); + scorpio::read_var(xs_long_file, "pressure", prs_h.data()); + + // read xsqy data (using lng_indexer_h for the first index) + // FIXME: hard-coded for only one photo reaction. + for(int m = 0; m < phtcnt; ++m) { + auto xsqy_ndx_h = ekat::subview(xsqy_h, m); + scorpio::read_var(xs_long_file, rxt_names[m], xsqy_h.data()); + } + + // populate etfphot by rebinning solar data + HostView1D wc_h("wc", nw), wlintv_h("wlintv", nw), we_h("we", nw + 1); + + scorpio::read_var(rsf_file, "wc", wc_h.data()); + scorpio::read_var(rsf_file, "wlintv", wlintv_h.data()); + for(int i = 0; i < nw; ++i) { + we_h(i) = wc_h(i) - 0.5 * wlintv_h(i); + } + we_h(nw) = wc_h(nw - 1) - 0.5 * wlintv_h(nw - 1); + // populate_etfphot(we_h, etfphot_h); + // FIXME: etfphot_data is hard-coded. + auto etfphot_data = populate_etfphot_from_e3sm_case(); + auto etfphot_h = HostView1D((Real *)etfphot_data.data(), nw); + + scorpio::release_file(rsf_file); + scorpio::release_file(xs_long_file); + + // copy host photolysis table into place on device + Kokkos::deep_copy(table.rsf_tab, rsf_tab_h); + Kokkos::deep_copy(table.xsqy, xsqy_h); + Kokkos::deep_copy(table.sza, sza_h); + Kokkos::deep_copy(table.alb, alb_h); + Kokkos::deep_copy(table.press, press_h); + Kokkos::deep_copy(table.colo3, colo3_h); + Kokkos::deep_copy(table.o3rat, o3rat_h); + Kokkos::deep_copy(table.etfphot, etfphot_h); + Kokkos::deep_copy(table.prs, prs_h); + // set pht_alias_mult_1 to 1 + Kokkos::deep_copy(table.pht_alias_mult_1, 1.0); + Kokkos::deep_copy(table.lng_indexer, lng_indexer_h); + + // compute gradients (on device) + Kokkos::parallel_for( + "del_p", nump - 1, KOKKOS_LAMBDA(int i) { + table.del_p(i) = 1.0 / haero::abs(table.press(i) - table.press(i + 1)); + }); + Kokkos::parallel_for( + "del_sza", numsza - 1, KOKKOS_LAMBDA(int i) { + table.del_sza(i) = 1.0 / (table.sza(i + 1) - table.sza(i)); + }); + Kokkos::parallel_for( + "del_alb", numalb - 1, KOKKOS_LAMBDA(int i) { + table.del_alb(i) = 1.0 / (table.alb(i + 1) - table.alb(i)); + }); + Kokkos::parallel_for( + "del_o3rat", numcolo3 - 1, KOKKOS_LAMBDA(int i) { + table.del_o3rat(i) = 1.0 / (table.o3rat(i + 1) - table.o3rat(i)); + }); + Kokkos::parallel_for( + "dprs", np_xs - 1, KOKKOS_LAMBDA(int i) { + table.dprs(i) = 1.0 / (table.prs(i) - table.prs(i + 1)); + }); + + return table; +} + +} // namespace scream::impl diff --git a/components/eamxx/src/physics/mam/readfiles/tracer_reader_utils.hpp b/components/eamxx/src/physics/mam/readfiles/tracer_reader_utils.hpp new file mode 100644 index 00000000000..2e34db2b496 --- /dev/null +++ b/components/eamxx/src/physics/mam/readfiles/tracer_reader_utils.hpp @@ -0,0 +1,741 @@ +#ifndef EAMXX_MAM_TRACER_READER_UTILS +#define EAMXX_MAM_TRACER_READER_UTILS + +#include +#include + +#include "share/grid/point_grid.hpp" +#include "share/grid/remap/coarsening_remapper.hpp" +#include "share/grid/remap/identity_remapper.hpp" +#include "share/grid/remap/refining_remapper_p2p.hpp" +#include "share/io/scorpio_input.hpp" +#include "share/io/scream_scorpio_interface.hpp" + +namespace scream::mam_coupling { + +using namespace ShortFieldTagsNames; +using view_1d_host = typename KT::view_1d::HostMirror; +using view_2d_host = typename KT::view_2d::HostMirror; + +using ExeSpace = typename KT::ExeSpace; +using ESU = ekat::ExeSpaceUtils; +using C = scream::physics::Constants; +using LIV = ekat::LinInterp; + +// Linoz NetCDF files use levs instead of formula_terms. +// This function allocates a view, so we need to do it during initialization. +// Thus, we assume that source pressure is independent of time, +// which is the case for Linoz files (zonal file). + +inline void compute_p_src_zonal_files(const view_1d &levs, + const view_2d &p_src) { + EKAT_REQUIRE_MSG(p_src.data() != 0, + "Error: p_src has not been allocated. \n"); + EKAT_REQUIRE_MSG(levs.data() != 0, "Error: levs has not been allocated. \n"); + const int ncol = p_src.extent(0); + const int nlevs_data = levs.extent(0); + EKAT_REQUIRE_MSG( + p_src.extent_int(1) == nlevs_data, + "Error: p_src has a different number of levels than the source data. \n"); + Kokkos::parallel_for( + "pressure_computation", + Kokkos::MDRangePolicy >({0, 0}, {ncol, nlevs_data}), + KOKKOS_LAMBDA(const int icol, const int kk) { + // mbar->pascals + // FIXME: Does EAMxx have a better method to + // convert units?" + p_src(icol, kk) = levs(kk) * 100; + }); + Kokkos::fence(); +} + +// We have a similar version in MAM4xx. +// This version was created because the data view cannot be modified +// inside the parallel_for. +// This struct will be used in init while reading nc files. +// The MAM4xx version will be used instead of parallel_for that loops over cols. +struct ForcingHelper { + // This index is in Fortran format. i.e. starts in 1 + int frc_ndx; + // does this file have altitude? + bool file_alt_data; + // number of sectors per forcing + int nsectors; + // offset in output vector from reader + int offset; +}; + +enum TracerFileType { + // file with PS ncol, lev, and time + FORMULA_PS, + // nc zonal file from ncremap + ZONAL, + // vertical emission files + VERT_EMISSION, +}; + +enum TracerDataIndex { BEG = 0, END = 1, OUT = 2 }; + +/* Maximum number of tracers (or fields) that the tracer reader can handle. + Note: We are not allocating memory for MAX_NVARS_TRACER tracers. + Therefore, if a file contains more than this number, it is acceptable to + increase this limit. Currently, Linoz files have 8 fields. */ +constexpr int MAX_NVARS_TRACER = 10; +constexpr int MAX_NUM_VERT_EMISSION_FIELDS = 25; + +// Linoz structures to help manage all of the variables: +struct TracerTimeState { + // 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; +}; // TricerTimeState + +struct TracerData { + TracerData() = default; + TracerData(const int ncol, const int nlev, const int nvars) { + init(ncol, nlev, nvars); + } + void init(const int ncol, const int nlev, const int nvars) { + ncol_ = ncol; + nlev_ = nlev; + nvars_ = nvars; + EKAT_REQUIRE_MSG( + nvars_ <= int(MAX_NVARS_TRACER), + "Error! Number of variables is bigger than NVARS_MAXTRACER. \n"); + } + + int ncol_{-1}; + int nlev_{-1}; + int nvars_{-1}; + int nlevs_data; + int ncols_data; + + // + int offset_time_index_{0}; + + // We cannot use a std::vector + // because we need to access these views from device. + // 0: beg 1: end 3: out + view_2d data[3][MAX_NVARS_TRACER]; + // type of file + TracerFileType file_type; + + // These views are employed in files with the PS variable + // 0: beg 1: end 3: out + view_1d ps[3]; + const_view_1d hyam; + const_view_1d hybm; + view_int_1d work_vert_inter[MAX_NVARS_TRACER]; + + // External forcing file (vertical emission) + // Uses altitude instead of pressure to interpolate data + view_1d altitude_int_; + // + bool has_altitude_{false}; + // pressure source for ZONAL or FORMULA_PS files + view_2d p_src_; + + // only for zonal files + view_1d zonal_levs_; + + void allocate_temporal_views() { + // BEG and OUT data views. + EKAT_REQUIRE_MSG(ncol_ != int(-1), "Error! ncols has not been set. \n"); + EKAT_REQUIRE_MSG(nlev_ != int(-1), "Error! nlevs has not been set. \n"); + EKAT_REQUIRE_MSG(nvars_ != int(-1), "Error! nvars has not been set. \n"); + + for(int ivar = 0; ivar < nvars_; ++ivar) { + data[TracerDataIndex::OUT][ivar] = view_2d("linoz_1_out", ncol_, nlev_); + data[TracerDataIndex::BEG][ivar] = view_2d("linoz_1_out", ncol_, nlev_); + } + + // for vertical interpolation using rebin routine + if(file_type == FORMULA_PS || file_type == ZONAL) { + // we only need work array for FORMULA_PS or ZONAL + for(int ivar = 0; ivar < nvars_; ++ivar) { + work_vert_inter[ivar] = + view_int_1d("allocate_work_vertical_interpolation", ncol_); + } + // we use ncremap and python scripts to convert zonal files to ne4pn4 + // grids. + p_src_ = view_2d("pressure_src_invariant", ncol_, nlev_); + } + + if(file_type == TracerFileType::FORMULA_PS) { + ps[TracerDataIndex::OUT] = view_1d("ps", ncol_); + ps[TracerDataIndex::BEG] = view_1d("ps", ncol_); + } + + if(file_type == TracerFileType::ZONAL) { + // we use ncremap and python scripts to convert zonal files to ne4pn4 + // grids. + compute_p_src_zonal_files(zonal_levs_, p_src_); + } + } +}; + +KOKKOS_INLINE_FUNCTION +Real linear_interp(const Real &x0, const Real &x1, const Real &t) { + return (1 - t) * x0 + t * x1; +} // linear_interp + +// time[3]={year,month, day} +inline util::TimeStamp convert_date(const int date) { + constexpr int ten_thousand = 10000; + constexpr int one_hundred = 100; + + int year = date / ten_thousand; + int month = (date - year * ten_thousand) / one_hundred; + int day = date - year * ten_thousand - month * one_hundred; + return util::TimeStamp(year, month, day, 0, 0, 0); +} +// FIXME: This function is not implemented in eamxx. +// FIXME: Assumes 365 days/year, 30 days/month; +// NOTE: that this assumption is mainly used for plotting. +// NOTE: This is not a direct port from EAM. +//We only use this routine for chlorine. +inline int compute_number_days_from_zero(const util::TimeStamp &ts) { + return ts.get_year() * 365 + ts.get_month() * 30 + ts.get_day(); +} + +inline void create_linoz_chlorine_reader( + const std::string &linoz_chlorine_file, const util::TimeStamp &model_time, + const int chlorine_loading_ymd, // in format YYYYMMDD + std::vector &values, std::vector &time_secs) { + auto time_stamp_beg = convert_date(chlorine_loading_ymd); + + const int offset_time = + compute_number_days_from_zero(time_stamp_beg) - compute_number_days_from_zero(model_time); + scorpio::register_file(linoz_chlorine_file, scorpio::Read); + const int nlevs_time = scorpio::get_time_len(linoz_chlorine_file); + for(int itime = 0; itime < nlevs_time; ++itime) { + int date; + scorpio::read_var(linoz_chlorine_file, "date", &date, itime); + if(date >= chlorine_loading_ymd) { + Real value; + scorpio::read_var(linoz_chlorine_file, "chlorine_loading", &value, itime); + values.push_back(value); + auto time_stamp = convert_date(date); + time_secs.push_back(compute_number_days_from_zero(time_stamp) - offset_time); + } + } // end itime + scorpio::release_file(linoz_chlorine_file); +} + +// Gets the times from the NC file +// Given a date in the format YYYYMMDD, returns its index in the time dimension. +inline void get_time_from_ncfile(const std::string &file_name, + const int cyclical_ymd, // in format YYYYMMDD + int &cyclical_ymd_index, + std::vector &dates) { + // in file_name: name of the NC file + // in cyclical_ymd: date in the format YYYYMMDD + // out cyclical_ymd_index: time index for cyclical_ymd + // out dates: date in YYYYMMDD format + scorpio::register_file(file_name, scorpio::Read); + const int nlevs_time = scorpio::get_time_len(file_name); + cyclical_ymd_index = -1; + for(int itime = 0; itime < nlevs_time; ++itime) { + int date; + scorpio::read_var(file_name, "date", &date, itime); + // std::cout << itime << " date: " << date << "\n"; + if(date >= cyclical_ymd && cyclical_ymd_index == -1) { + cyclical_ymd_index = itime; + } + dates.push_back(date); + } // end itime + + EKAT_REQUIRE_MSG(cyclical_ymd_index >= 0, + "Error! Current model time (" + + std::to_string(cyclical_ymd) + ") is not within " + + "Tracer time period: [" + std::to_string(dates[0]) + + ", " + "(" + std::to_string(dates[nlevs_time - 1]) + + ").\n"); + scorpio::release_file(file_name); +} + +inline Real chlorine_loading_advance(const util::TimeStamp &ts, + std::vector &values, + std::vector &time_secs) { + const int current_time = compute_number_days_from_zero(ts); + int index = 0; + // update index + for(int i = 0; i < int(values.size()); i++) { + if(current_time > time_secs[i]) { + index = i; + break; + } + } // + + const Real delt = (current_time - time_secs[index]) / + (time_secs[index + 1] - time_secs[index]); + return values[index] + delt * (values[index + 1] - values[index]); +} + +// It reads variables that are not time-dependent and independent of columns (no +// MPI involved here). We also obtain the offset_time_index using a date +// (cyclical_ymd) as input. We initialize a few members of tracer_data. +inline void setup_tracer_data(TracerData &tracer_data, // out + const std::string &trace_data_file, // in + const int cyclical_ymd) // in +{ + scorpio::register_file(trace_data_file, scorpio::Read); + // by default, I am assuming a zonal file. + TracerFileType tracer_file_type = ZONAL; + + int nlevs_data = -1; + if(scorpio::has_var(trace_data_file, "lev")) { + nlevs_data = scorpio::get_dimlen(trace_data_file, "lev"); + } + const bool has_altitude = scorpio::has_var(trace_data_file, "altitude"); + + // This type of files use altitude (zi) for vertical interpolation + if(has_altitude) { + nlevs_data = scorpio::get_dimlen(trace_data_file, "altitude"); + tracer_file_type = VERT_EMISSION; + } + EKAT_REQUIRE_MSG( + nlevs_data != -1, + "Error: The file does not contain either lev or altitude. \n"); + + const int ncols_data = scorpio::get_dimlen(trace_data_file, "ncol"); + + // This type of files use model pressure (pmid) for vertical interpolation + if(scorpio::has_var(trace_data_file, "PS")) { + tracer_file_type = FORMULA_PS; + view_1d_host hyam_h("hyam_h", nlevs_data); + view_1d_host hybm_h("hybm_h", nlevs_data); + + scorpio::read_var(trace_data_file, "hyam", hyam_h.data()); + scorpio::read_var(trace_data_file, "hybm", hybm_h.data()); + view_1d hyam("hyam", nlevs_data); + view_1d hybm("hybm", nlevs_data); + Kokkos::deep_copy(hyam, hyam_h); + Kokkos::deep_copy(hybm, hybm_h); + tracer_data.hyam = hyam; + tracer_data.hybm = hybm; + } + + if(tracer_file_type == ZONAL) { + view_1d_host levs_h("levs_h", nlevs_data); + view_1d levs("levs", nlevs_data); + scorpio::read_var(trace_data_file, "lev", levs_h.data()); + Kokkos::deep_copy(levs, levs_h); + tracer_data.zonal_levs_ = levs; + } + + if(tracer_file_type == VERT_EMISSION) { + const int nilevs_data = + scorpio::get_dimlen(trace_data_file, "altitude_int"); + view_1d_host altitude_int_host("altitude_int_host", nilevs_data); + view_1d altitude_int = view_1d("altitude_int", nilevs_data); + scorpio::read_var(trace_data_file, "altitude_int", + altitude_int_host.data()); + Kokkos::deep_copy(altitude_int, altitude_int_host); + tracer_data.altitude_int_ = altitude_int; + } + // time index + { + const int nlevs_time = scorpio::get_dimlen(trace_data_file, "time"); + int cyclical_ymd_index = -1; + for(int itime = 0; itime < nlevs_time; ++itime) { + int date; + scorpio::read_var(trace_data_file, "date", &date, itime); + if(date >= cyclical_ymd) { + cyclical_ymd_index = itime; + break; + } + } // end itime + + EKAT_REQUIRE_MSG(cyclical_ymd_index >= 0, "Error! Current model time (" + + std::to_string(cyclical_ymd) + + ") is not within " + + "Tracer time period.\n"); + + tracer_data.offset_time_index_ = cyclical_ymd_index; + } + + scorpio::release_file(trace_data_file); + tracer_data.file_type = tracer_file_type; + tracer_data.nlevs_data = nlevs_data; + tracer_data.ncols_data = ncols_data; + tracer_data.has_altitude_ = has_altitude; +} +inline std::shared_ptr create_horiz_remapper( + const std::shared_ptr &model_grid, + const std::string &trace_data_file, const std::string &map_file, + const std::vector &var_names, TracerData &tracer_data) { + using namespace ShortFieldTagsNames; + // 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("tracer_horiz_interp_tgt_grid", true); + horiz_interp_tgt_grid->reset_num_vertical_lev(tracer_data.nlevs_data); + + if(tracer_data.file_type == VERT_EMISSION) { + horiz_interp_tgt_grid->reset_field_tag_name(LEV, "altitude"); + horiz_interp_tgt_grid->reset_field_tag_name(ILEV, "altitude_int"); + } + + const int ncols_model = model_grid->get_num_global_dofs(); + std::shared_ptr remapper; + if(tracer_data.ncols_data == ncols_model) { + remapper = std::make_shared( + horiz_interp_tgt_grid, IdentityRemapper::SrcAliasTgt); + } else { + EKAT_REQUIRE_MSG(tracer_data.ncols_data <= ncols_model, + "Error! We do not allow to coarsen tracer external " + "forcing data to fit the " + "model. We only allow\n" + " tracer external forcing 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: tracer external forcing data is on a different grid than the " + "model one,\n" + " but tracer external forcing data remap file is missing from " + "tracer external forcing data 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 layout_3d_mid = tgt_grid->get_3d_scalar_layout(true); + + const auto nondim = ekat::units::Units::nondimensional(); + + for(auto var_name : var_names) { + Field ifield( + FieldIdentifier(var_name, layout_3d_mid, nondim, tgt_grid->name())); + ifield.allocate_view(); + remapper->register_field_from_tgt(ifield); + } + // zonal files do not have the PS variable. + if(tracer_data.file_type == FORMULA_PS) { + Field ps(FieldIdentifier("PS", layout_2d, nondim, tgt_grid->name())); + ps.allocate_view(); + remapper->register_field_from_tgt(ps); + } + remapper->registration_ends(); + return remapper; + +} // create_horiz_remapper + +inline std::shared_ptr create_tracer_data_reader( + const std::shared_ptr &horiz_remapper, + const std::string &tracer_data_file) { + std::vector io_fields; + for(int i = 0; i < horiz_remapper->get_num_fields(); ++i) { + io_fields.push_back(horiz_remapper->get_src_field(i)); + } + const auto io_grid = horiz_remapper->get_src_grid(); + return std::make_shared(tracer_data_file, io_grid, io_fields, + true); +} // create_tracer_data_reader + +inline void update_tracer_data_from_file( + const std::shared_ptr &scorpio_reader, + const int time_index, // zero-based + AbstractRemapper &tracer_horiz_interp, TracerData &tracer_data) { + // 1. read from field + scorpio_reader->read_variables(time_index); + // 2. Run the horiz remapper (it is a do-nothing op if tracer external forcing + // data is on same grid as model) + tracer_horiz_interp.remap(/*forward = */ true); + // + const int nvars = tracer_data.nvars_; + // + for(int i = 0; i < nvars; ++i) { + tracer_data.data[TracerDataIndex::END][i] = + tracer_horiz_interp.get_tgt_field(i).get_view(); + } + + if(tracer_data.file_type == FORMULA_PS) { + // Recall, the fields are registered in the order: tracers, ps + // 3. Copy from the tgt field of the remapper into the spa_data + tracer_data.ps[TracerDataIndex::END] = + tracer_horiz_interp.get_tgt_field(nvars).get_view(); + } + +} // update_tracer_data_from_file +inline void update_tracer_timestate( + const std::shared_ptr &scorpio_reader, + const util::TimeStamp &ts, AbstractRemapper &tracer_horiz_interp, + TracerTimeState &time_state, TracerData &data_tracer) { + // Now we check if we have to update the data that changes monthly + // NOTE: This means that tracer external forcing 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) { + const auto tracer_data = data_tracer.data; + const int nvars = data_tracer.nvars_; + const auto ps = data_tracer.ps; + + // Update the tracer external forcing 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 spa_end'data into spa_beg'data, and read in the new spa_end + for(int ivar = 0; ivar < nvars; ++ivar) { + Kokkos::deep_copy(tracer_data[TracerDataIndex::BEG][ivar], + tracer_data[TracerDataIndex::END][ivar]); + } + + if(data_tracer.file_type == FORMULA_PS) { + Kokkos::deep_copy(ps[TracerDataIndex::BEG], ps[TracerDataIndex::END]); + } + + // Following SPA to time-interpolate data in MAM4xx + // Assume the data is saved monthly and cycles in one year + // Add offset_time_index to support cases where data is saved + // from other periods of time. + // Update the tracer external 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 = + data_tracer.offset_time_index_ + (time_state.current_month + 1) % 12; + update_tracer_data_from_file(scorpio_reader, next_month, + tracer_horiz_interp, data_tracer); + } + +} // END update_tracer_timestate + +// This function is based on the SPA::perform_time_interpolation function. +inline void perform_time_interpolation(const TracerTimeState &time_state, + const TracerData &data_tracer) { + // 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. + // 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; + + // We can ||ize over columns as well as over variables and bands + const auto& data = data_tracer.data; + const auto& file_type = data_tracer.file_type; + + const auto& ps = data_tracer.ps; + const int num_vars = data_tracer.nvars_; + + const int ncol = data_tracer.ncol_; + const int num_vert = data_tracer.nlev_; + + const int outer_iters = ncol * num_vars; + + const auto policy = ESU::get_default_team_policy(outer_iters, num_vert); + + 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"); + + Kokkos::parallel_for( + "linoz_time_interp_loop", policy, KOKKOS_LAMBDA(const Team &team) { + // The policy is over ncols*num_vars, so retrieve icol/ivar + const int icol = team.league_rank() / num_vars; + const int ivar = team.league_rank() % num_vars; + + // Get column of beg/end/out variable + auto var_beg = ekat::subview(data[TracerDataIndex::BEG][ivar], icol); + auto var_end = ekat::subview(data[TracerDataIndex::END][ivar], icol); + auto var_out = ekat::subview(data[TracerDataIndex::OUT][ivar], icol); + + Kokkos::parallel_for( + Kokkos::TeamVectorRange(team, num_vert), [&](const int &k) { + var_out(k) = + linear_interp(var_beg(k), var_end(k), delta_t_fraction); + }); + // linoz files do not have ps variables. + if(ivar == 1 && file_type == FORMULA_PS) { + ps[TracerDataIndex::OUT](icol) = + linear_interp(ps[TracerDataIndex::BEG](icol), + ps[TracerDataIndex::END](icol), delta_t_fraction); + } + }); + Kokkos::fence(); +} // perform_time_interpolation + +inline void compute_source_pressure_levels(const view_1d &ps_src, + const view_2d &p_src, + const const_view_1d &hyam, + const const_view_1d &hybm) { + constexpr auto P0 = C::P0; + const int ncols = ps_src.extent(0); + const int num_vert_packs = p_src.extent(1); + const auto policy = ESU::get_default_team_policy(ncols, num_vert_packs); + + Kokkos::parallel_for( + "tracer_compute_p_src_loop", policy, KOKKOS_LAMBDA(const Team &team) { + const int icol = team.league_rank(); + Kokkos::parallel_for( + Kokkos::TeamVectorRange(team, num_vert_packs), [&](const int k) { + p_src(icol, k) = ps_src(icol) * hybm(k) + P0 * hyam(k); + }); + }); +} // compute_source_pressure_levels + +inline void perform_vertical_interpolation(const view_2d &p_src_c, + const const_view_2d &p_tgt_c, + const TracerData &input, + const view_2d output[]) { + const int ncol = input.ncol_; + const int levsiz = input.nlev_; + const int pver = mam4::nlev; + + const int num_vars = input.nvars_; + // make a local copy of output + view_2d output_local[MAX_NVARS_TRACER]; + EKAT_REQUIRE_MSG( + num_vars <= int(MAX_NVARS_TRACER), + "Error! Number of variables is bigger than NVARS_MAXTRACER. \n"); + for (int ivar = 0; ivar < num_vars; ++ivar) + { + // At this stage, begin/end must have the same horiz dimensions + EKAT_REQUIRE(input.ncol_ == output[ivar].extent_int(0)); + output_local[ivar] = output[ivar]; + } + const int outer_iters = ncol * num_vars; + const auto policy_setup = ESU::get_default_team_policy(outer_iters, pver); + const auto& data = input.data; + + Kokkos::parallel_for( + "vert_interp", policy_setup, + KOKKOS_LAMBDA(typename LIV::MemberType const &team) { + // The policy is over ncols*num_vars, so retrieve icol/ivar + const int icol = team.league_rank() / num_vars; + const int ivar = team.league_rank() % num_vars; + const auto pin_at_icol = ekat::subview(p_src_c, icol); + const auto pmid_at_icol = ekat::subview(p_tgt_c, icol); + const auto& datain = data[TracerDataIndex::OUT][ivar]; + const auto datain_at_icol = ekat::subview(datain, icol); + const auto dataout = output_local[ivar]; + const auto dataout_at_icol = ekat::subview(dataout, icol); + + mam4::vertical_interpolation::vert_interp( + team, levsiz, pver, pin_at_icol, pmid_at_icol, datain_at_icol, + dataout_at_icol); + }); +} + +inline void perform_vertical_interpolation(const const_view_1d &altitude_int, + const const_view_2d &zi, + const TracerData &input, + const view_2d output[]) { + EKAT_REQUIRE_MSG( + input.file_type == VERT_EMISSION, + "Error! vertical interpolation only with altitude variable. \n"); + const int ncols = input.ncol_; + const int num_vars = input.nvars_; + const int num_vertical_lev_target = output[0].extent(1); + const int num_vert_packs = num_vertical_lev_target; + const int outer_iters = ncols * num_vars; + const auto policy_interp = + ESU::get_default_team_policy(outer_iters, num_vert_packs); + // FIXME: Get m2km from emaxx. + const Real m2km = 1e-3; + const auto &src_x = altitude_int; + const int nsrc = input.nlev_; + constexpr int pver = mam4::nlev; + const int pverp = pver + 1; + const auto& data = input.data; + + // make a local copy of output + view_2d output_local[MAX_NVARS_TRACER]; + EKAT_REQUIRE_MSG( + num_vars <= int(MAX_NVARS_TRACER), + "Error! Number of variables is bigger than NVARS_MAXTRACER. \n"); + for (int ivar = 0; ivar < num_vars; ++ivar) + { + // At this stage, begin/end must have the same horiz dimensions + EKAT_REQUIRE(input.ncol_ == output[ivar].extent_int(0)); + output_local[ivar] = output[ivar]; + } + Kokkos::parallel_for( + "tracer_vert_interp_loop", policy_interp, + KOKKOS_LAMBDA(const Team &team) { + const int icol = team.league_rank() / num_vars; + const int ivar = team.league_rank() % num_vars; + + const auto src = ekat::subview(data[TracerDataIndex::OUT][ivar], icol); + const auto trg = ekat::subview(output_local[ivar], icol); + // FIXME: Try to avoid copy of trg_x by modifying rebin + // trg_x + Real trg_x[pver + 1]; + // I am trying to do this: + // model_z(1:pverp) = m2km * state(c)%zi(i,pverp:1:-1) + for(int i = 0; i < pverp; ++i) { + trg_x[pverp - i - 1] = m2km * zi(icol, i); + } + team.team_barrier(); + mam4::vertical_interpolation::rebin(team, nsrc, num_vertical_lev_target, + src_x, trg_x, src, trg); + }); +} + +inline void advance_tracer_data( + const std::shared_ptr &scorpio_reader, // in + AbstractRemapper &tracer_horiz_interp, // out + const util::TimeStamp &ts, // in + TracerTimeState &time_state, TracerData &data_tracer, // out + const const_view_2d &p_tgt, const const_view_2d &zi_tgt, // in + const view_2d output[]) { // out + /* Update the TracerTimeState to reflect the current time, note the addition + * of dt */ + time_state.t_now = ts.frac_of_year_in_days(); + /* Update time state and if the month has changed, update the data.*/ + update_tracer_timestate(scorpio_reader, ts, tracer_horiz_interp, time_state, + data_tracer); + // Step 1. Perform time interpolation + perform_time_interpolation(time_state, data_tracer); + + if(data_tracer.file_type == FORMULA_PS) { + // Step 2. Compute source pressure levels + const auto ps = data_tracer.ps[TracerDataIndex::OUT]; + compute_source_pressure_levels(ps, data_tracer.p_src_, data_tracer.hyam, + data_tracer.hybm); + } + + // Step 3. Perform vertical interpolation + if(data_tracer.file_type == FORMULA_PS || data_tracer.file_type == ZONAL) { + perform_vertical_interpolation(data_tracer.p_src_, p_tgt, data_tracer, + output); + } else if(data_tracer.file_type == VERT_EMISSION) { + perform_vertical_interpolation(data_tracer.altitude_int_, zi_tgt, + data_tracer, output); + } + +} // advance_tracer_data + +} // namespace scream::mam_coupling +#endif // EAMXX_MAM_HELPER_MICRO diff --git a/components/eamxx/src/physics/register_physics.hpp b/components/eamxx/src/physics/register_physics.hpp index d837e96311f..99956bb75f5 100644 --- a/components/eamxx/src/physics/register_physics.hpp +++ b/components/eamxx/src/physics/register_physics.hpp @@ -65,7 +65,7 @@ inline void register_physics () { proc_factory.register_product("Nudging",&create_atmosphere_process); #endif #ifdef EAMXX_HAS_MAM - proc_factory.register_product("mam4_micro",&create_atmosphere_process); + proc_factory.register_product("mam4_aero_microphys",&create_atmosphere_process); proc_factory.register_product("mam4_optics",&create_atmosphere_process); proc_factory.register_product("mam4_drydep",&create_atmosphere_process); proc_factory.register_product("mam4_aci",&create_atmosphere_process); diff --git a/components/eamxx/tests/multi-process/dynamics_physics/mam/homme_shoc_cld_spa_p3_rrtmgp_mam4_wetscav/output.yaml b/components/eamxx/tests/multi-process/dynamics_physics/mam/homme_shoc_cld_spa_p3_rrtmgp_mam4_wetscav/output.yaml index 28c3e31ff19..a48a4eba3c8 100644 --- a/components/eamxx/tests/multi-process/dynamics_physics/mam/homme_shoc_cld_spa_p3_rrtmgp_mam4_wetscav/output.yaml +++ b/components/eamxx/tests/multi-process/dynamics_physics/mam/homme_shoc_cld_spa_p3_rrtmgp_mam4_wetscav/output.yaml @@ -103,7 +103,7 @@ Fields: - aerdepwetis - aerdepwetcw - dgnumwet - - dgncur_a + - dgnum - wetdens - qaerwat - bc_c1 diff --git a/components/eamxx/tests/multi-process/physics_only/mam/p3_mam4_wetscav/output.yaml b/components/eamxx/tests/multi-process/physics_only/mam/p3_mam4_wetscav/output.yaml index c26db8eb681..7d4d355ecc6 100644 --- a/components/eamxx/tests/multi-process/physics_only/mam/p3_mam4_wetscav/output.yaml +++ b/components/eamxx/tests/multi-process/physics_only/mam/p3_mam4_wetscav/output.yaml @@ -57,7 +57,7 @@ Field Names: - aerdepwetis - aerdepwetcw - dgnumwet - - dgncur_a + - dgnum - wetdens - qaerwat output_control: diff --git a/components/eamxx/tests/multi-process/physics_only/mam/shoc_cldfrac_p3_wetscav/output.yaml b/components/eamxx/tests/multi-process/physics_only/mam/shoc_cldfrac_p3_wetscav/output.yaml index 1ae53e6c1c3..da4a83ad7b1 100644 --- a/components/eamxx/tests/multi-process/physics_only/mam/shoc_cldfrac_p3_wetscav/output.yaml +++ b/components/eamxx/tests/multi-process/physics_only/mam/shoc_cldfrac_p3_wetscav/output.yaml @@ -91,7 +91,7 @@ Field Names: - aerdepwetis - aerdepwetcw - dgnumwet - - dgncur_a + - dgnum - wetdens - qaerwat output_control: diff --git a/components/eamxx/tests/single-process/CMakeLists.txt b/components/eamxx/tests/single-process/CMakeLists.txt index 40565d94bf8..2d0ae749c49 100644 --- a/components/eamxx/tests/single-process/CMakeLists.txt +++ b/components/eamxx/tests/single-process/CMakeLists.txt @@ -25,6 +25,7 @@ if (SCREAM_ENABLE_MAM) add_subdirectory(mam/wet_scav) add_subdirectory(mam/emissions) add_subdirectory(mam/constituent_fluxes) + add_subdirectory(mam/aero_microphys) endif() if (SCREAM_TEST_LEVEL GREATER_EQUAL SCREAM_TEST_LEVEL_EXPERIMENTAL) add_subdirectory(zm) diff --git a/components/eamxx/tests/single-process/mam/aero_microphys/CMakeLists.txt b/components/eamxx/tests/single-process/mam/aero_microphys/CMakeLists.txt new file mode 100644 index 00000000000..e1ba85a685f --- /dev/null +++ b/components/eamxx/tests/single-process/mam/aero_microphys/CMakeLists.txt @@ -0,0 +1,64 @@ +include (ScreamUtils) + +set (TEST_BASE_NAME mam4_aero_microphys_standalone) +set (FIXTURES_BASE_NAME ${TEST_BASE_NAME}_generate_output_nc_files) + +# Create the test +CreateADUnitTest(${TEST_BASE_NAME} + LABELS mam4_aero_microphys 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 +set (TEST_INPUT_FILES + scream/init/${EAMxx_tests_IC_FILE_MAM4xx_72lev} + scream/mam4xx/linoz/ne2np4/linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne2np4_c20240724.nc + scream/mam4xx/invariants/ne2np4/oxid_ne2np4_L26_1850-2015_c20240827.nc + scream/mam4xx/linoz/Linoz_Chlorine_Loading_CMIP6_0003-2017_c20171114.nc + scream/mam4xx/photolysis/RSF_GT200nm_v3.0_c080811.nc + scream/mam4xx/photolysis/temp_prs_GT200nm_JPL10_c130206.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_so2_elev_ne2np4_2010_clim_c20240726.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_so4_a1_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_so4_a2_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_pom_a4_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_bc_a4_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_num_a1_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_num_a2_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_num_a4_elev_ne2np4_2010_clim_c20240823.nc + scream/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_soag_elev_ne2np4_2010_clim_c20240823.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_aero_microphys 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/aero_microphys/input.yaml b/components/eamxx/tests/single-process/mam/aero_microphys/input.yaml new file mode 100644 index 00000000000..fa4538bb75a --- /dev/null +++ b/components/eamxx/tests/single-process/mam/aero_microphys/input.yaml @@ -0,0 +1,65 @@ +%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_aero_microphys] + mam4_aero_microphys: + mam4_do_cond: true + mam4_do_newnuc: true + mam4_do_coag: true + mam4_do_rename: true + mam4_o3_tau: 172800.0 + mam4_o3_sfc: 3.0E-008 + mam4_o3_lbl: 4 + mam4_psc_T : 193.0 + mam4_linoz_ymd : 20100101 + mam4_linoz_file_name : ${SCREAM_DATA_DIR}/mam4xx/linoz/ne2np4/linoz1850-2015_2010JPL_CMIP6_10deg_58km_ne2np4_c20240724.nc + mam4_oxid_file_name : ${SCREAM_DATA_DIR}/mam4xx/invariants/ne2np4/oxid_ne2np4_L26_1850-2015_c20240827.nc + mam4_oxid_ymd : 20150101 + mam4_linoz_chlorine_file : ${SCREAM_DATA_DIR}/mam4xx/linoz/Linoz_Chlorine_Loading_CMIP6_0003-2017_c20171114.nc + mam4_chlorine_loading_ymd : 20100101 + mam4_rsf_file : ${SCREAM_DATA_DIR}/mam4xx/photolysis/RSF_GT200nm_v3.0_c080811.nc + mam4_xs_long_file : ${SCREAM_DATA_DIR}/mam4xx/photolysis/temp_prs_GT200nm_JPL10_c130206.nc + verti_emiss_ymd : 20100101 + mam4_so2_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_so2_elev_ne2np4_2010_clim_c20240726.nc + mam4_so4_a1_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_so4_a1_elev_ne2np4_2010_clim_c20240823.nc + mam4_so4_a2_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_so4_a2_elev_ne2np4_2010_clim_c20240823.nc + mam4_pom_a4_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_pom_a4_elev_ne2np4_2010_clim_c20240823.nc + mam4_bc_a4_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated/cmip6_mam4_bc_a4_elev_ne2np4_2010_clim_c20240823.nc + mam4_num_a1_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_num_a1_elev_ne2np4_2010_clim_c20240823.nc + mam4_num_a2_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_num_a2_elev_ne2np4_2010_clim_c20240823.nc + mam4_num_a4_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_num_a4_elev_ne2np4_2010_clim_c20240823.nc + mam4_soag_verti_emiss_file_name : ${SCREAM_DATA_DIR}/mam4xx/emissions/ne2np4/elevated//cmip6_mam4_soag_elev_ne2np4_2010_clim_c20240823.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 + pbl_height: 1.0 + #These should come from the input file + dgnum: [1.246662106183775E-007, 4.081134799487888E-008, 1.103139143795796E-006, 1.000000011686097E-007] + dgnumwet: [2.367209731605067E-007, 6.780643470563889E-008, 3.028011448344027E-006, 1.000000096285154E-007] + wetdens: [1038.67760516297, 1046.20002003441, 1031.74623165457, 1086.79731859184] + +# The parameters for I/O control +Scorpio: + output_yaml_files: ["output.yaml"] +... diff --git a/components/eamxx/tests/single-process/mam/aero_microphys/output.yaml b/components/eamxx/tests/single-process/mam/aero_microphys/output.yaml new file mode 100644 index 00000000000..388781fb8f0 --- /dev/null +++ b/components/eamxx/tests/single-process/mam/aero_microphys/output.yaml @@ -0,0 +1,18 @@ +%YAML 1.1 +--- +filename_prefix: mam4_aero_microphys_standalone_output +Averaging Type: Instant +Fields: + Physics: + Field Names: + - T_mid + # To save these fields make sure to turn on + # OUTPUT_TRACER_FIELDS + #- oxi_fields + #- linoz_fields + #- vertical_emission_fields + +output_control: + Frequency: 1 + frequency_units: nsteps +... diff --git a/components/eamxx/tests/single-process/mam/wet_scav/output.yaml b/components/eamxx/tests/single-process/mam/wet_scav/output.yaml index 9cd0a16bc41..854d18a2e48 100644 --- a/components/eamxx/tests/single-process/mam/wet_scav/output.yaml +++ b/components/eamxx/tests/single-process/mam/wet_scav/output.yaml @@ -8,7 +8,7 @@ Fields: - aerdepwetis - aerdepwetcw - dgnumwet - - dgncur_a + - dgnum - wetdens - qaerwat - bc_c1 diff --git a/externals/mam4xx b/externals/mam4xx index 2eee1988b1e..aae46807bf5 160000 --- a/externals/mam4xx +++ b/externals/mam4xx @@ -1 +1 @@ -Subproject commit 2eee1988b1e6488c904e0a28564bdcadfa190d34 +Subproject commit aae46807bf58d6ffcbc6db620c27568450b2d040