Skip to content

Commit

Permalink
EAMxx: Fixes after rebase and fixes a possible race condition in sethet
Browse files Browse the repository at this point in the history
  • Loading branch information
singhbalwinder committed Nov 23, 2024
1 parent 8b32ee5 commit 3e55254
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 60 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,8 @@
#include <physics/mam/eamxx_mam_microphysics_process_interface.hpp>
// impl namespace for some driver level functions for microphysics

#include "readfiles/photo_table_utils.cpp"
#include "readfiles/find_season_index_utils.hpp"
#include "physics/rrtmgp/shr_orb_mod_c2f.hpp"
#include "readfiles/find_season_index_utils.hpp"
#include "readfiles/photo_table_utils.cpp"

namespace scream {
Expand Down Expand Up @@ -136,6 +135,13 @@ void MAMMicrophysics::set_grids(
// Total cloud fraction [fraction]
add_field<Required>("cldfrac_liq", scalar3d_mid, nondim, grid_name);

// Evaporation from stratiform rain [kg/kg/s]
add_field<Required>("nevapr", scalar3d_mid, kg / kg / s, grid_name);

// Stratiform rain production rate [kg/kg/s]
add_field<Required>("precip_total_tend", scalar3d_mid, kg / kg / s,
grid_name);

// precipitation liquid mass [kg/m2]
add_field<Required>("precip_liq_surf_mass", scalar3d_mid, kg / m2, grid_name);

Expand Down Expand Up @@ -248,7 +254,7 @@ void MAMMicrophysics::set_grids(
LinozHorizInterp_, linoz_file_name_);

// linoz reader
const auto io_grid_linoz = LinozHorizInterp_->get_src_grid();
const auto io_grid_linoz = LinozHorizInterp_->get_tgt_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 =
Expand Down Expand Up @@ -276,7 +282,7 @@ void MAMMicrophysics::set_grids(
TracerHorizInterp_, oxid_file_name_);

const int nvars = int(var_names.size());
const auto io_grid = TracerHorizInterp_->get_src_grid();
const auto io_grid = TracerHorizInterp_->get_tgt_grid();
const int num_cols_io =
io_grid->get_num_local_dofs(); // Number of columns on this rank
const int num_levs_io =
Expand All @@ -301,78 +307,90 @@ void MAMMicrophysics::set_grids(
"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";
std::string item_name = "mam4_" + var_name + "_elevated_emiss_file_name";
const auto file_name = m_params.get<std::string>(item_name);
vert_emis_file_name_[var_name] = file_name;
elevated_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"] = {
elevated_emis_var_names_["so2"] = {"BB", "ENE_ELEV", "IND_ELEV",
"contvolc"};
elevated_emis_var_names_["so4_a1"] = {"BB", "ENE_ELEV", "IND_ELEV",
"contvolc"};
elevated_emis_var_names_["so4_a2"] = {"contvolc"};
elevated_emis_var_names_["pom_a4"] = {"BB"};
elevated_emis_var_names_["bc_a4"] = {"BB"};
elevated_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"};
elevated_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"};
elevated_emis_var_names_["num_a4"] = {"num_a1_BC_ELEV_BB",
"num_a1_POM_ELEV_BB"};
elevated_emis_var_names_["soag"] = {"SOAbb_src", "SOAbg_src", "SOAff_src"};

int verti_emiss_cyclical_ymd = m_params.get<int>("verti_emiss_ymd");
int elevated_emiss_cyclical_ymd = m_params.get<int>("elevated_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];
const auto file_name = elevated_emis_file_name_[var_name];
const auto var_names = elevated_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);
elevated_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

auto file_reader = scream::mam_coupling::create_tracer_data_reader(
hor_rem, file_name, data_tracer.file_type);
ElevatedEmissionsHorizInterp_.push_back(hor_rem);
ElevatedEmissionsDataReader_.push_back(file_reader);
elevated_emis_data_.push_back(data_tracer);
} // var_name elevated 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 auto file_name = elevated_emis_file_name_[var_name];
const auto var_names = elevated_emis_var_names_[var_name];
const int nvars = static_cast<int>(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();
forcings_[i].frc_ndx = i + 1;
const auto io_grid_emis =
ElevatedEmissionsHorizInterp_[i]->get_tgt_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_;
elevated_emis_data_[i].init(num_cols_io_emis, num_levs_io_emis, nvars);
elevated_emis_data_[i].allocate_temporal_views();
forcings_[i].file_alt_data = elevated_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_);
elevated_emis_output_[isp + offset_emis_ver] =
view_2d("elevated_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),
offset_emis_ver <= int(mam_coupling::MAX_NUM_ELEVATED_EMISSIONS_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");
"MAX_NUM_ELEVATED_EMISSIONS_FIELDS. Increase the "
"MAX_NUM_ELEVATED_EMISSIONS_FIELDS in tracer_reader_utils.hpp \n");

} // Tracer external forcing data

{
const std::string season_wes_file =
m_params.get<std::string>("mam4_season_wes_file");
const auto &clat = col_latitudes_;
mam_coupling::find_season_index_reader(season_wes_file, clat,
index_season_lai_);
}
} // set_grids

// ================================================================
Expand Down Expand Up @@ -544,6 +562,9 @@ void MAMMicrophysics::initialize_impl(const RunType run_type) {

const int photo_table_len = get_photo_table_work_len(photo_table_);
work_photo_table_ = view_2d("work_photo_table", ncol_, photo_table_len);
const int sethet_work_len = mam4::mo_sethet::get_total_work_len_sethet();
work_set_het_ = view_2d("work_set_het_array", ncol_, sethet_work_len);
cmfdqr_ = view_1d("cmfdqr_", nlev_);

// here's where we store per-column photolysis rates
photo_rates_ = view_3d("photo_rates", ncol_, nlev_, mam4::mo_photo::phtcnt);
Expand All @@ -561,8 +582,8 @@ void MAMMicrophysics::initialize_impl(const RunType run_type) {

for(int i = 0; i < static_cast<int>(extfrc_lst_.size()); ++i) {
scream::mam_coupling::update_tracer_data_from_file(
VertEmissionsDataReader_[i], curr_month, *VertEmissionsHorizInterp_[i],
vert_emis_data_[i]);
ElevatedEmissionsDataReader_[i], curr_month,
*ElevatedEmissionsHorizInterp_[i], elevated_emis_data_[i]);
}

invariants_ = view_3d("invarians", ncol_, nlev_, mam4::gas_chemistry::nfs);
Expand Down Expand Up @@ -599,6 +620,15 @@ void MAMMicrophysics::run_impl(const double dt) {
Kokkos::parallel_for("preprocess", scan_policy, preprocess_);
Kokkos::fence();

//----------- Variables from microphysics scheme -------------

// Evaporation from stratiform rain [kg/kg/s]
const auto &nevapr = get_field_in("nevapr").get_view<const Real **>();

// Stratiform rain production rate [kg/kg/s]
const auto &prain =
get_field_in("precip_total_tend").get_view<const Real **>();

const auto wet_geometric_mean_diameter_i =
get_field_in("dgnumwet").get_view<const Real ***>();
const auto dry_geometric_mean_diameter_i =
Expand Down Expand Up @@ -701,20 +731,21 @@ void MAMMicrophysics::run_impl(const double dt) {
linoz_output); // out
Kokkos::fence();

vert_emiss_time_state_.t_now = ts.frac_of_year_in_days();
int i = 0;
elevated_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 auto file_name = elevated_emis_file_name_[var_name];
const auto var_names = elevated_emis_var_names_[var_name];
const int nsectors = int(var_names.size());
view_2d vert_emis_output[nsectors];
view_2d elevated_emis_output[nsectors];
for(int isp = 0; isp < nsectors; ++isp) {
vert_emis_output[isp] = vert_emis_output_[isp + forcings_[i].offset];
elevated_emis_output[isp] =
elevated_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);
ElevatedEmissionsDataReader_[i], *ElevatedEmissionsHorizInterp_[i], ts,
elevated_emiss_time_state_, elevated_emis_data_[i], dry_atm_.p_mid,
dry_atm_.z_iface, elevated_emis_output);
i++;
Kokkos::fence();
}
Expand Down Expand Up @@ -789,10 +820,10 @@ void MAMMicrophysics::run_impl(const double dt) {
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 auto elevated_emis_output = elevated_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];
Expand All @@ -807,6 +838,8 @@ void MAMMicrophysics::run_impl(const double dt) {
clsmap_4[i] = mam4::gas_chemistry::clsmap_4[i];
permute_4[i] = mam4::gas_chemistry::permute_4[i];
}
const auto &cmfdqr = cmfdqr_;
const auto &work_set_het = work_set_het_;
const mam4::seq_drydep::Data drydep_data =
mam4::seq_drydep::set_gas_drydep_data();
const auto qv = wet_atm_.qv;
Expand Down Expand Up @@ -848,7 +881,7 @@ void MAMMicrophysics::run_impl(const double dt) {
// 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];
const auto field = elevated_emis_output[isec + forcings[i].offset];
forcings_in[i].fields_data[isec] = ekat::subview(field, icol);
}
} // extcnt for loop
Expand Down Expand Up @@ -879,6 +912,9 @@ void MAMMicrophysics::run_impl(const double dt) {
ekat::subview(linoz_dPmL_dO3col, icol);
const auto linoz_cariolle_pscs_icol =
ekat::subview(linoz_cariolle_pscs, icol);
const auto nevapr_icol = ekat::subview(nevapr, icol);
const auto prain_icol = ekat::subview(prain, icol);
const auto work_set_het_icol = ekat::subview(work_set_het, icol);

// Surface temperature
const Real sfc_air_temp = atm.temperature(surface_lev);
Expand Down Expand Up @@ -940,9 +976,14 @@ void MAMMicrophysics::run_impl(const double dt) {

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, drydep_data, dvel, dflx, progs);

// update constituent fluxes with gas drydep fluxes (dflx)
wet_diameter_icol, wetdens_icol, dry_atm.phis(icol), cmfdqr,
prain_icol, nevapr_icol, work_set_het_icol, drydep_data, dvel, dflx,
progs);

// Update constituent fluxes with gas drydep fluxes (dflx)
// FIXME: Possible units mismatch (dflx is in kg/cm2/s but
// constituent_fluxes is kg/m2/s) (Following mimics Fortran code
// behavior but we should look into it)
for(int ispc = offset_aerosol; ispc < mam4::pcnst; ++ispc) {
constituent_fluxes(icol, ispc) = dflx[ispc - offset_aerosol];
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ class MAMMicrophysics final : public scream::AtmosphereProcess {

using view_1d_host = typename KT::view_1d<Real>::HostMirror;

using view_int_2d = typename KT::template view_2d<int>;
using view_int_2d = typename KT::template view_2d<int>;

// a thread team dispatched to a single vertical column
using ThreadTeam = mam4::ThreadTeam;
Expand Down Expand Up @@ -239,7 +239,8 @@ class MAMMicrophysics final : public scream::AtmosphereProcess {
std::vector<mam_coupling::TracerData> elevated_emis_data_;
std::map<std::string, std::string> elevated_emis_file_name_;
std::map<std::string, std::vector<std::string>> elevated_emis_var_names_;
view_2d elevated_emis_output_[mam_coupling::MAX_NUM_ELEVATED_EMISSIONS_FIELDS];
view_2d
elevated_emis_output_[mam_coupling::MAX_NUM_ELEVATED_EMISSIONS_FIELDS];
view_3d extfrc_;
mam_coupling::ForcingHelper forcings_[mam4::gas_chemistry::extcnt];

Expand Down

0 comments on commit 3e55254

Please sign in to comment.