diff --git a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.cpp b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.cpp index 5031a75b5c6..d20f64ec0f6 100644 --- a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.cpp +++ b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.cpp @@ -64,19 +64,19 @@ void SHOCMacrophysics::set_grids(const std::shared_ptr grids const auto s2 = s*s; // These variables are needed by the interface, but not actually passed to shoc_main. - add_field("omega", scalar3d_layout_mid, Pa/s, grid_name, ps); - add_field("surf_sens_flux", scalar2d_layout_col, W/m2, grid_name); - add_field("surf_mom_flux", surf_mom_flux_layout, N/m2, grid_name); - - add_field("surf_evap", scalar2d_layout_col, kg/m2/s, grid_name); - add_field ("T_mid", scalar3d_layout_mid, K, grid_name, ps); - add_field ("qv", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps); + add_field("omega", scalar3d_layout_mid, Pa/s, grid_name, ps); + add_field("surf_sens_flux", scalar2d_layout_col, W/m2, grid_name); + add_field("surf_evap", scalar2d_layout_col, kg/m2/s, grid_name); + add_field("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("apply_tms", false)) { add_field("surf_drag_coeff_tms", scalar2d_layout_col, kg/s/m2, grid_name); } + add_field ("T_mid", scalar3d_layout_mid, K, grid_name, ps); + add_field ("qv", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps); + // Input variables add_field("p_mid", scalar3d_layout_mid, Pa, grid_name, ps); add_field("p_int", scalar3d_layout_int, Pa, grid_name, ps); @@ -85,11 +85,11 @@ void SHOCMacrophysics::set_grids(const std::shared_ptr grids // Input/Output variables add_field("tke", scalar3d_layout_mid, m2/s2, grid_name, "tracers", ps); - add_field("horiz_winds", horiz_wind_layout, m/s, grid_name, ps); - add_field("sgs_buoy_flux", scalar3d_layout_mid, K*(m/s), grid_name, ps); - add_field("eddy_diff_mom", scalar3d_layout_mid, m2/s, grid_name, ps); + add_field("horiz_winds", horiz_wind_layout, m/s, grid_name, ps); + add_field("sgs_buoy_flux", scalar3d_layout_mid, K*(m/s), grid_name, ps); + add_field("eddy_diff_mom", scalar3d_layout_mid, m2/s, grid_name, ps); add_field("qc", scalar3d_layout_mid, Qunit, grid_name, "tracers", ps); - add_field("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name, ps); + add_field("cldfrac_liq", scalar3d_layout_mid, nondim, grid_name, ps); // Output variables add_field("pbl_height", scalar2d_layout_col, m, grid_name); @@ -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 @@ -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(); - const auto& p_mid = get_field_in("p_mid").get_view(); - const auto& p_int = get_field_in("p_int").get_view(); - const auto& pseudo_density = get_field_in("pseudo_density").get_view(); - const auto& omega = get_field_in("omega").get_view(); - const auto& surf_sens_flux = get_field_in("surf_sens_flux").get_view(); - const auto& surf_evap = get_field_in("surf_evap").get_view(); - const auto& surf_mom_flux = get_field_in("surf_mom_flux").get_view(); - const auto& qtracers = get_group_out("tracers").m_bundle->get_view(); - const auto& qc = get_field_out("qc").get_view(); - const auto& qv = get_field_out("qv").get_view(); - const auto& tke = get_field_out("tke").get_view(); - const auto& cldfrac_liq = get_field_out("cldfrac_liq").get_view(); - const auto& sgs_buoy_flux = get_field_out("sgs_buoy_flux").get_view(); - const auto& tk = get_field_out("eddy_diff_mom").get_view(); - const auto& inv_qc_relvar = get_field_out("inv_qc_relvar").get_view(); - const auto& phis = get_field_in("phis").get_view(); + const auto& T_mid = get_field_out("T_mid").get_view(); + const auto& p_mid = get_field_in("p_mid").get_view(); + const auto& p_int = get_field_in("p_int").get_view(); + const auto& pseudo_density = get_field_in("pseudo_density").get_view(); + const auto& omega = get_field_in("omega").get_view(); + const auto& surf_sens_flux = get_field_in("surf_sens_flux").get_view(); + const auto& surf_evap = get_field_in("surf_evap").get_view(); + const auto& surf_mom_flux = get_field_in("surf_mom_flux").get_view(); + const auto& qtracers = get_group_out("tracers").m_bundle->get_view(); + const auto& qc = get_field_out("qc").get_view(); + const auto& qv = get_field_out("qv").get_view(); + const auto& tke = get_field_out("tke").get_view(); + const auto& cldfrac_liq = get_field_out("cldfrac_liq").get_view(); + const auto& sgs_buoy_flux = get_field_out("sgs_buoy_flux").get_view(); + const auto& tk = get_field_out("eddy_diff_mom").get_view(); + const auto& inv_qc_relvar = get_field_out("inv_qc_relvar").get_view(); + const auto& phis = get_field_in("phis").get_view(); // 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; @@ -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); @@ -431,6 +433,16 @@ void SHOCMacrophysics::run_impl (const double dt) const auto scan_policy = ekat::ExeSpaceUtils::get_thread_range_parallel_scan_team_policy(m_num_cols, nlev_packs); const auto default_policy = ekat::ExeSpaceUtils::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()); + if (m_params.get("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", @@ -442,15 +454,11 @@ void SHOCMacrophysics::run_impl (const double dt) apply_turbulent_mountain_stress(); } - if (m_params.get("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(round(hdtime/dt)),1); + m_hdtime = dt; + m_nadv = std::max(static_cast(round(m_hdtime/dt)),1); // Reset internal WSM variables. workspace_mgr.reset_internals(); @@ -499,11 +507,10 @@ void SHOCMacrophysics::check_flux_state_consistency(const double dt) { using PC = scream::physics::Constants; const Real gravit = PC::gravit; - const Real qmin = 1e-12; // minimum permitted constituent concentration (kg/kg) - const auto& pseudo_density = get_field_in ("pseudo_density").get_view(); - const auto& surf_evap = get_field_out("surf_evap").get_view(); - const auto& qv = get_field_out("qv").get_view(); + const auto& pseudo_density = get_field_in ("pseudo_density").get_view(); + 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(); const auto nlevs = m_num_levs; const auto nlev_packs = ekat::npack(nlevs); @@ -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) { + // Calculate the total column vapor mass + Real total_column_vap_mass = ekat::ExeSpaceUtils::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::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(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)); + }); } }); } diff --git a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp index 6a6298755d6..a12072ceffe 100644 --- a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp +++ b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp @@ -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; @@ -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 @@ -408,6 +408,7 @@ class SHOCMacrophysics : public scream::AtmosphereProcess uview_1d wprtp_sfc; uview_1d upwp_sfc; uview_1d vpwp_sfc; + uview_1d modified_surf_evap; #ifdef SCREAM_SMALL_KERNELS uview_1d se_b; uview_1d ke_b; @@ -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::view_1d m_cell_area; KokkosTypes::view_1d m_cell_lat;