Skip to content

Update check_flux_state_consistency() computation #2525

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Update check_flux_state_consistency() computation #2525

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1fd2fea
Update check_flux_state_consistency() computation
tcclevenger Sep 1, 2023
840e717
Use original surf_evap in conservation checks
tcclevenger Sep 28, 2023
ade1a8e
Various name changes and indent fixes
tcclevenger Sep 28, 2023
0cdcd75
Adjust current computation to include excess_mass
tcclevenger Oct 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
173 changes: 97 additions & 76 deletions components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,19 +64,19 @@ void SHOCMacrophysics::set_grids(const std::shared_ptr<const GridsManager> grids
const auto s2 = s*s;

// These variables are needed by the interface, but not actually passed to shoc_main.
add_field<Required>("omega", scalar3d_layout_mid, Pa/s, grid_name, ps);
add_field<Required>("surf_sens_flux", scalar2d_layout_col, W/m2, grid_name);
add_field<Required>("surf_mom_flux", surf_mom_flux_layout, N/m2, grid_name);

add_field<Updated>("surf_evap", scalar2d_layout_col, kg/m2/s, grid_name);
add_field<Updated> ("T_mid", scalar3d_layout_mid, K, grid_name, ps);
add_field<Updated> ("qv", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps);
add_field<Required>("omega", scalar3d_layout_mid, Pa/s, grid_name, ps);
add_field<Required>("surf_sens_flux", scalar2d_layout_col, W/m2, grid_name);
add_field<Required>("surf_evap", scalar2d_layout_col, kg/m2/s, grid_name);
add_field<Required>("surf_mom_flux", surf_mom_flux_layout, N/m2, grid_name);

// If TMS is a process, add surface drag coefficient to required fields
if (m_params.get<bool>("apply_tms", false)) {
add_field<Required>("surf_drag_coeff_tms", scalar2d_layout_col, kg/s/m2, grid_name);
}

add_field<Updated> ("T_mid", scalar3d_layout_mid, K, grid_name, ps);
add_field<Updated> ("qv", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps);

// Input variables
add_field<Required>("p_mid", scalar3d_layout_mid, Pa, grid_name, ps);
add_field<Required>("p_int", scalar3d_layout_int, Pa, grid_name, ps);
Expand All @@ -85,11 +85,11 @@ void SHOCMacrophysics::set_grids(const std::shared_ptr<const GridsManager> grids

// Input/Output variables
add_field<Updated>("tke", scalar3d_layout_mid, m2/s2, grid_name, "tracers", ps);
add_field<Updated>("horiz_winds", horiz_wind_layout, m/s, grid_name, ps);
add_field<Updated>("sgs_buoy_flux", scalar3d_layout_mid, K*(m/s), grid_name, ps);
add_field<Updated>("eddy_diff_mom", scalar3d_layout_mid, m2/s, grid_name, ps);
add_field<Updated>("horiz_winds", horiz_wind_layout, m/s, grid_name, ps);
add_field<Updated>("sgs_buoy_flux", scalar3d_layout_mid, K*(m/s), grid_name, ps);
add_field<Updated>("eddy_diff_mom", scalar3d_layout_mid, m2/s, grid_name, ps);
add_field<Updated>("qc", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps);
add_field<Updated>("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name, ps);
add_field<Updated>("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name, ps);

// Output variables
add_field<Computed>("pbl_height", scalar2d_layout_col, m, grid_name);
Expand Down Expand Up @@ -159,7 +159,8 @@ void SHOCMacrophysics::init_buffers(const ATMBufferManager &buffer_manager)
// 1d scalar views
using scalar_view_t = decltype(m_buffer.cell_length);
scalar_view_t* _1d_scalar_view_ptrs[Buffer::num_1d_scalar_ncol] =
{&m_buffer.cell_length, &m_buffer.wpthlp_sfc, &m_buffer.wprtp_sfc, &m_buffer.upwp_sfc, &m_buffer.vpwp_sfc
{&m_buffer.cell_length, &m_buffer.wpthlp_sfc, &m_buffer.wprtp_sfc,
&m_buffer.upwp_sfc, &m_buffer.vpwp_sfc, &m_buffer.modified_surf_evap
#ifdef SCREAM_SMALL_KERNELS
, &m_buffer.se_b, &m_buffer.ke_b, &m_buffer.wv_b, &m_buffer.wl_b
, &m_buffer.se_a, &m_buffer.ke_a, &m_buffer.wv_a, &m_buffer.wl_a
Expand Down Expand Up @@ -233,47 +234,48 @@ void SHOCMacrophysics::initialize_impl (const RunType run_type)
// Initialize all of the structures that are passed to shoc_main in run_impl.
// Note: Some variables in the structures are not stored in the field manager. For these
// variables a local view is constructed.
const auto& T_mid = get_field_out("T_mid").get_view<Spack**>();
const auto& p_mid = get_field_in("p_mid").get_view<const Spack**>();
const auto& p_int = get_field_in("p_int").get_view<const Spack**>();
const auto& pseudo_density = get_field_in("pseudo_density").get_view<const Spack**>();
const auto& omega = get_field_in("omega").get_view<const Spack**>();
const auto& surf_sens_flux = get_field_in("surf_sens_flux").get_view<const Real*>();
const auto& surf_evap = get_field_in("surf_evap").get_view<const Real*>();
Copy link
Contributor

@oksanaguba oksanaguba Sep 28, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

needs to be grabbed back

const auto& surf_mom_flux = get_field_in("surf_mom_flux").get_view<const Real**>();
const auto& qtracers = get_group_out("tracers").m_bundle->get_view<Spack***>();
const auto& qc = get_field_out("qc").get_view<Spack**>();
const auto& qv = get_field_out("qv").get_view<Spack**>();
const auto& tke = get_field_out("tke").get_view<Spack**>();
const auto& cldfrac_liq = get_field_out("cldfrac_liq").get_view<Spack**>();
const auto& sgs_buoy_flux = get_field_out("sgs_buoy_flux").get_view<Spack**>();
const auto& tk = get_field_out("eddy_diff_mom").get_view<Spack**>();
const auto& inv_qc_relvar = get_field_out("inv_qc_relvar").get_view<Spack**>();
const auto& phis = get_field_in("phis").get_view<const Real*>();
const auto& T_mid = get_field_out("T_mid").get_view<Spack**>();
const auto& p_mid = get_field_in("p_mid").get_view<const Spack**>();
const auto& p_int = get_field_in("p_int").get_view<const Spack**>();
const auto& pseudo_density = get_field_in("pseudo_density").get_view<const Spack**>();
const auto& omega = get_field_in("omega").get_view<const Spack**>();
const auto& surf_sens_flux = get_field_in("surf_sens_flux").get_view<const Real*>();
const auto& surf_evap = get_field_in("surf_evap").get_view<const Real*>();
const auto& surf_mom_flux = get_field_in("surf_mom_flux").get_view<const Real**>();
const auto& qtracers = get_group_out("tracers").m_bundle->get_view<Spack***>();
const auto& qc = get_field_out("qc").get_view<Spack**>();
const auto& qv = get_field_out("qv").get_view<Spack**>();
const auto& tke = get_field_out("tke").get_view<Spack**>();
const auto& cldfrac_liq = get_field_out("cldfrac_liq").get_view<Spack**>();
const auto& sgs_buoy_flux = get_field_out("sgs_buoy_flux").get_view<Spack**>();
const auto& tk = get_field_out("eddy_diff_mom").get_view<Spack**>();
const auto& inv_qc_relvar = get_field_out("inv_qc_relvar").get_view<Spack**>();
const auto& phis = get_field_in("phis").get_view<const Real*>();

// Alias local variables from temporary buffer
auto z_mid = m_buffer.z_mid;
auto z_int = m_buffer.z_int;
auto cell_length = m_buffer.cell_length;
auto wpthlp_sfc = m_buffer.wpthlp_sfc;
auto wprtp_sfc = m_buffer.wprtp_sfc;
auto upwp_sfc = m_buffer.upwp_sfc;
auto vpwp_sfc = m_buffer.vpwp_sfc;
auto rrho = m_buffer.rrho;
auto rrho_i = m_buffer.rrho_i;
auto thv = m_buffer.thv;
auto dz = m_buffer.dz;
auto zt_grid = m_buffer.zt_grid;
auto zi_grid = m_buffer.zi_grid;
auto wtracer_sfc = m_buffer.wtracer_sfc;
auto wm_zt = m_buffer.wm_zt;
auto inv_exner = m_buffer.inv_exner;
auto thlm = m_buffer.thlm;
auto qw = m_buffer.qw;
auto dse = m_buffer.dse;
auto tke_copy = m_buffer.tke_copy;
auto qc_copy = m_buffer.qc_copy;
auto shoc_ql2 = m_buffer.shoc_ql2;
auto z_mid = m_buffer.z_mid;
auto z_int = m_buffer.z_int;
auto cell_length = m_buffer.cell_length;
auto wpthlp_sfc = m_buffer.wpthlp_sfc;
auto wprtp_sfc = m_buffer.wprtp_sfc;
auto upwp_sfc = m_buffer.upwp_sfc;
auto vpwp_sfc = m_buffer.vpwp_sfc;
auto rrho = m_buffer.rrho;
auto rrho_i = m_buffer.rrho_i;
auto thv = m_buffer.thv;
auto dz = m_buffer.dz;
auto zt_grid = m_buffer.zt_grid;
auto zi_grid = m_buffer.zi_grid;
auto wtracer_sfc = m_buffer.wtracer_sfc;
auto wm_zt = m_buffer.wm_zt;
auto inv_exner = m_buffer.inv_exner;
auto thlm = m_buffer.thlm;
auto qw = m_buffer.qw;
auto dse = m_buffer.dse;
auto tke_copy = m_buffer.tke_copy;
auto qc_copy = m_buffer.qc_copy;
auto shoc_ql2 = m_buffer.shoc_ql2;
auto modified_surf_evap = m_buffer.modified_surf_evap;

// For now, set z_int(i,nlevs) = z_surf = 0
const Real z_surf = 0.0;
Expand All @@ -288,7 +290,7 @@ void SHOCMacrophysics::initialize_impl (const RunType run_type)
}

shoc_preprocess.set_variables(m_num_cols,m_num_levs,m_num_tracers,z_surf,m_cell_area,m_cell_lat,
T_mid,p_mid,p_int,pseudo_density,omega,phis,surf_sens_flux,surf_evap,
T_mid,p_mid,p_int,pseudo_density,omega,phis,surf_sens_flux,modified_surf_evap,
surf_mom_flux,qtracers,qv,qc,qc_copy,tke,tke_copy,z_mid,z_int,cell_length,
dse,rrho,rrho_i,thv,dz,zt_grid,zi_grid,wpthlp_sfc,wprtp_sfc,upwp_sfc,vpwp_sfc,
wtracer_sfc,wm_zt,inv_exner,thlm,qw);
Expand Down Expand Up @@ -431,6 +433,16 @@ void SHOCMacrophysics::run_impl (const double dt)
const auto scan_policy = ekat::ExeSpaceUtils<KT::ExeSpace>::get_thread_range_parallel_scan_team_policy(m_num_cols, nlev_packs);
const auto default_policy = ekat::ExeSpaceUtils<KT::ExeSpace>::get_default_team_policy(m_num_cols, nlev_packs);

// The global surf_evap field should not change between iterations, but
// within an iteration, if the param check_flux_state_consistency=true,
// surf_evap used for shoc preprocessor may need to be updated. The local
// modified_surf_evap view is used for input in preprocessor for both param cases.
Kokkos::deep_copy(m_buffer.modified_surf_evap,
get_field_in("surf_evap").get_view<const Real*>());
if (m_params.get<bool>("check_flux_state_consistency", false)) {
check_flux_state_consistency(dt);
}

// Preprocessing of SHOC inputs. Kernel contains a parallel_scan,
// so a special TeamPolicy is required.
Kokkos::parallel_for("shoc_preprocess",
Expand All @@ -442,15 +454,11 @@ void SHOCMacrophysics::run_impl (const double dt)
apply_turbulent_mountain_stress();
}

if (m_params.get<bool>("check_flux_state_consistency", false)) {
check_flux_state_consistency(dt);
}

// For now set the host timestep to the shoc timestep. This forces
// number of SHOC timesteps (nadv) to be 1.
// TODO: input parameter?
hdtime = dt;
m_nadv = std::max(static_cast<int>(round(hdtime/dt)),1);
m_hdtime = dt;
m_nadv = std::max(static_cast<int>(round(m_hdtime/dt)),1);

// Reset internal WSM variables.
workspace_mgr.reset_internals();
Expand Down Expand Up @@ -499,11 +507,10 @@ void SHOCMacrophysics::check_flux_state_consistency(const double dt)
{
using PC = scream::physics::Constants<Real>;
const Real gravit = PC::gravit;
const Real qmin = 1e-12; // minimum permitted constituent concentration (kg/kg)
Copy link
Contributor

@oksanaguba oksanaguba Sep 28, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bogensch @AaronDonahue @PeterCaldwell this lower bound was taken from EAM setup (iirc EAM won't run if this bount is set to 0.0). Do we need a certain nontrivial bound in scream or can we always assume lower acceptable bounds for tracers are zero. In this case it is only for vapor. If a certain bound is needed, Conrad will put it in physics constants.

Copy link
Contributor

@oksanaguba oksanaguba Sep 28, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe it is a question for @brhillman instead, since in EAM the comment mentions radiation

    real(r8),intent(in)    :: qminc  ! minimum value of mass mixing ratio (kg/kg)
                                     ! normally 0., except water 1.E-12, for radiation.


const auto& pseudo_density = get_field_in ("pseudo_density").get_view<const Spack**>();
const auto& surf_evap = get_field_out("surf_evap").get_view<Real*>();
const auto& qv = get_field_out("qv").get_view<Spack**>();
const auto& pseudo_density = get_field_in ("pseudo_density").get_view<const Spack**>();
const auto& modified_surf_evap = m_buffer.modified_surf_evap; // Use/update the local modified_surf_evap view
const auto& qv = get_field_out("qv").get_view<Spack**>();

const auto nlevs = m_num_levs;
const auto nlev_packs = ekat::npack<Spack>(nlevs);
Expand All @@ -518,28 +525,42 @@ void SHOCMacrophysics::check_flux_state_consistency(const double dt)
const auto& pseudo_density_i = ekat::subview(pseudo_density, i);
const auto& qv_i = ekat::subview(qv, i);

// reciprocal of pseudo_density at the bottom layer
const auto rpdel = 1.0/pseudo_density_i(last_pack_idx)[last_pack_entry];
// Define lambda to get column vapor mass on level pack k
auto column_vapor_mass_on_level = [&](const int k) {
return qv_i(k)*pseudo_density_i(k);
};

// Define column_vapor_mass on surface
const auto column_vapor_mass_sfc = column_vapor_mass_on_level(last_pack_idx)[last_pack_entry];

// Define surface evaporative mass
const auto surf_evap_mass = modified_surf_evap(i)*dt*gravit;

// Check if the negative surface latent heat flux can exhaust
// the moisture in the lowest model level. If so, apply fixer.
const auto condition = surf_evap(i) - (qmin - qv_i(last_pack_idx)[last_pack_entry])/(dt*gravit*rpdel);
if (condition < 0) {
const auto cc = abs(surf_evap(i)*dt*gravit);
const auto excess_mass = surf_evap_mass + column_vapor_mass_sfc;
if (excess_mass < 0) {
Copy link
Contributor

@oksanaguba oksanaguba Sep 28, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

qmin?

// Calculate the total column vapor mass
Real total_column_vap_mass = ekat::ExeSpaceUtils<KT::ExeSpace>::view_reduction(team, 0, nlevs, column_vapor_mass_on_level);
EKAT_KERNEL_ASSERT_MSG(total_column_vap_mass >= std::abs(surf_evap_mass),
"Error! Total mass of column vapor should be greater "
"than the magnitude of mass of surf_evap.\n");

auto tracer_mass = [&](const int k) {
return qv_i(k)*pseudo_density_i(k);
};
Real mm = ekat::ExeSpaceUtils<KT::ExeSpace>::view_reduction(team, 0, nlevs, tracer_mass);
// Set modified_surf_evap to exactly exhaust the moisture at the bottom layer
modified_surf_evap(i) = -1.0*column_vapor_mass_sfc/(dt*gravit);

EKAT_KERNEL_ASSERT_MSG(mm >= cc, "Error! Total mass of column vapor should be greater than mass of surf_evap.\n");
// Set total_column_vap_mass to only be from levels 0:nlev-2 (exclude sfc value)
total_column_vap_mass -= column_vapor_mass_sfc;

// Redistribute excess mass
Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_packs), [&](const int& k) {
const auto adjust = cc*qv_i(k)*pseudo_density_i(k)/mm;
qv_i(k) = (qv_i(k)*pseudo_density_i(k) - adjust)/pseudo_density_i(k);
});
// Define mask for accessing all but surface level
const auto range_pack = ekat::range<IntSmallPack>(k*Spack::n);
const Smask index_above_sfc = range_pack < nlevs-1;

surf_evap(i) = 0;
const auto adjust = excess_mass*column_vapor_mass_on_level(k)/total_column_vap_mass;
qv_i(k).set(index_above_sfc, (column_vapor_mass_on_level(k) + adjust)/pseudo_density_i(k));
});
}
});
}
Expand Down
9 changes: 5 additions & 4 deletions components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,7 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
} // set_variables

void set_mass_and_energy_fluxes (const view_1d_const& surf_evap_, const view_1d_const surf_sens_flux_,
const view_1d& vapor_flux_, const view_1d& water_flux_,
const view_1d& vapor_flux_, const view_1d& water_flux_,
const view_1d& ice_flux_, const view_1d& heat_flux_)
{
compute_mass_and_energy_fluxes = true;
Expand All @@ -389,9 +389,9 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
// Structure for storing local variables initialized using the ATMBufferManager
struct Buffer {
#ifndef SCREAM_SMALL_KERNELS
static constexpr int num_1d_scalar_ncol = 5;
static constexpr int num_1d_scalar_ncol = 6;
#else
static constexpr int num_1d_scalar_ncol = 18;
static constexpr int num_1d_scalar_ncol = 19;
#endif
static constexpr int num_1d_scalar_nlev = 1;
#ifndef SCREAM_SMALL_KERNELS
Expand All @@ -408,6 +408,7 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
uview_1d<Real> wprtp_sfc;
uview_1d<Real> upwp_sfc;
uview_1d<Real> vpwp_sfc;
uview_1d<Real> modified_surf_evap;
#ifdef SCREAM_SMALL_KERNELS
uview_1d<Real> se_b;
uview_1d<Real> ke_b;
Expand Down Expand Up @@ -504,7 +505,7 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
Int m_npbl;
Int m_nadv;
Int m_num_tracers;
Int hdtime;
Int m_hdtime;

KokkosTypes<DefaultDevice>::view_1d<const Real> m_cell_area;
KokkosTypes<DefaultDevice>::view_1d<const Real> m_cell_lat;
Expand Down