From 53cba060f460bef4636249d701401927d5e5059b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 25 Oct 2024 12:11:15 -0400 Subject: [PATCH] Implements CCPPized check_energy in SIMA. Squashed commit of the following: commit 8b8a8e089b578c7c336098c2831f462ec24811c9 Fix merge conflicts commit bef25fbc46642a65a6ec390e83ffdb537424db67 Update atmos_phys external commit 240ef548bc3573097a739c67a5f3cac671022677 Set wv_idx in air_composition total_hours_in_debugging_one_line = 6 commit 38fdbb21c4da0c625f9c964d503b12ad909c3775 Read cp_or_cv_dycore and identify dycore information from CAM snapshot commit 8fca3a42b6db828e185ec4df51665e1dcae8eae4 Update vcoord to energy_formula commit 0d16be936ddfbb1e0c3a223671fadf15e8c955cd New const_get_index logic without cam_ccpp_cap dependency commit e13e7ea0de0592bd2517cfcac32ea0bf3fabc6a6 Fix build (part 1); update standard names and atmos_phys external commit c015e28ab9ea5e590230c05df7c5411d7adf8b6a Update ncar-physics external commit 2372f7fc2fc1c71900f57ce1131c2b323f8803c3 Fix for b4b to CAM-SIMA (w/ History); add notes on modifications and/or divergence from CAM and/or -SIMA commit 4ce842786b8b64374057058dd22a21dc3d2576eb Add gmean_mod to src/utils including support infra in physics_grid.F90 commit 7493a1dfa70ecd6e837308e087a0cff4dc9cd2c2 Fix build issues; update factor in get_cp call to 1/dp commit 041cdfb250dc0088f87d6c1ecb391c015e6b4b57 Add is_first_timestep registry commit 333ad4eed8520f92c61a6731bfd589d2f873f097 Initial port of updates to energy budget (port of CAM #761) into CAM-SIMA commit 86489702461192ccb13ae7cff9d13f0db004812b Initial implementation of cp_to_cv_dycore/cam_thermo_water_update from stash Needs work to split into cam_thermo_water_update. This is WIP code --- .gitmodules | 4 +- src/control/cam_comp.F90 | 7 + src/data/air_composition.F90 | 151 +++++++++--- src/data/cam_thermo.F90 | 204 ++++++++++------ src/data/cam_thermo_formula.F90 | 37 +++ src/data/cam_thermo_formula.meta | 17 ++ src/data/registry.xml | 15 ++ src/dynamics/mpas/dyn_comp.F90 | 5 + src/dynamics/none/dyn_comp.F90 | 3 + src/dynamics/none/dyn_grid.F90 | 61 +++++ src/dynamics/se/dp_coupling.F90 | 122 +++++++--- src/dynamics/se/dycore/prim_advance_mod.F90 | 6 +- src/dynamics/se/dyn_comp.F90 | 3 + src/dynamics/utils/dyn_thermo.F90 | 33 +-- src/physics/ncar_ccpp | 2 +- src/physics/utils/phys_comp.F90 | 2 + src/physics/utils/physics_grid.F90 | 52 ++++ src/utils/gmean_mod.F90 | 256 ++++++++++++++++++++ tools/stdnames_to_inputnames_dictionary.xml | 5 + 19 files changed, 819 insertions(+), 166 deletions(-) create mode 100644 src/data/cam_thermo_formula.F90 create mode 100644 src/data/cam_thermo_formula.meta create mode 100644 src/utils/gmean_mod.F90 diff --git a/.gitmodules b/.gitmodules index 94610494..5f4d639c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -19,8 +19,8 @@ fxDONOTUSEurl = https://github.com/MPAS-Dev/MPAS-Model.git [submodule "ncar-physics"] path = src/physics/ncar_ccpp - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = e95c172d7a5a0ebf054f420b08416228e211baa3 + url = https://github.com/jimmielin/atmospheric_physics + fxtag = 952ebdd fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 538bd641..27fd57a7 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -29,6 +29,7 @@ module cam_comp use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys use physics_types, only: calday + use physics_types, only: is_first_timestep use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -151,6 +152,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & dtime_phys = 0.0_r8 call mark_as_initialized('timestep_for_physics') + is_first_timestep = .true. + call mark_as_initialized('is_first_timestep') + call init_pio_subsystem() ! Initializations using data passed from coupler. @@ -263,6 +267,9 @@ subroutine cam_timestep_init() use phys_comp, only: phys_timestep_init use stepon, only: stepon_timestep_init + ! Update timestep flags in physics state + is_first_timestep = is_first_step() + !---------------------------------------------------------- ! First phase of dynamics (at least couple from dynamics to physics) ! Return time-step for physics from dynamics. diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index f202ae97..f450c182 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -1,4 +1,5 @@ -! air_composition module defines major species of the atmosphere and manages the physical properties that are dependent on the composition of air +! air_composition module defines major species of the atmosphere and manages +! the physical properties that are dependent on the composition of air module air_composition use ccpp_kinds, only: kind_phys @@ -6,13 +7,16 @@ module air_composition use runtime_obj, only: unset_real, unset_int use phys_vars_init_check, only: std_name_len use physics_types, only: cpairv, rairv, cappav, mbarv, zvirv + use physics_types, only: cp_or_cv_dycore implicit none private save public :: air_composition_init - public :: air_composition_update + public :: dry_air_composition_update + public :: water_composition_update + ! get_cp_dry: (generalized) heat capacity for dry air public :: get_cp_dry ! get_cp: (generalized) heat capacity @@ -225,7 +229,7 @@ subroutine air_composition_init() ! !************************************************************************ ! - ! add prognostic components of dry air + ! add prognostic components of air ! !************************************************************************ ! @@ -309,6 +313,7 @@ subroutine air_composition_init() ! case(wv_stdname) !water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water call air_species_info(wv_stdname, ix, mw) + wv_idx = ix ! set water species index for use in get_hydrostatic_energy thermodynamic_active_species_idx(icnst) = ix thermodynamic_active_species_cp (icnst) = cpwv thermodynamic_active_species_cv (icnst) = cv3 / mw @@ -510,26 +515,71 @@ end subroutine air_composition_init !=========================================================================== !----------------------------------------------------------------------- - ! air_composition_update: Update the physics "constants" that vary + ! dry_air_composition_update: Update the physics "constants" that vary !------------------------------------------------------------------------- !=========================================================================== - subroutine air_composition_update(mmr, ncol, to_moist_factor) + subroutine dry_air_composition_update(mmr, ncol, to_dry_factor) - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + !(mmr = dry mixing ratio, if not, use to_dry_factor to convert!) + real(kind_phys), intent(in) :: mmr(:,:,:) ! mixing ratios for species dependent dry air integer, intent(in) :: ncol ! number of columns - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, & - rairv(:ncol, :), fact=to_moist_factor) + rairv(:ncol, :), fact=to_dry_factor) call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - cpairv(:ncol,:), fact=to_moist_factor) + cpairv(:ncol,:), fact=to_dry_factor) call get_mbarv(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - mbarv(:ncol,:), fact=to_moist_factor) + mbarv(:ncol,:), fact=to_dry_factor) cappav(:ncol,:) = rairv(:ncol,:) / cpairv(:ncol,:) - end subroutine air_composition_update + end subroutine dry_air_composition_update + + !=========================================================================== + !--------------------------------------------------------------------------- + ! water_composition_update: Update generalized cp or cv depending on dycore + !--------------------------------------------------------------------------- + !=========================================================================== + subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: energy_formula ! energy formula for dynamical core + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + character(len=*), parameter :: subname = 'water_composition_update' + + ! update enthalpy or internal energy scaling factor for energy consistency with CAM physics + ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module + + if (energy_formula == ENERGY_FORMULA_DYCORE_FV) then + ! FV: moist pressure vertical coordinate does not need update. + else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then + ! SE + + ! **TEMP** TODO hplin 9/17/24: + ! for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) + ! (whereas CAM only allocates 1-indexed) I subset it here. This should be verified during code + ! review. + call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & + factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx(1:), & + cpdry=cpairv(:ncol,:)) + else if (energy_formula == ENERGY_FORMULA_DYCORE_MPAS) then + ! MPAS + call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), & + cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) + + ! internal energy coefficient for MPAS + ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) + cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) + else + call endrun(subname//': dycore energy formula not supported') + end if + + end subroutine water_composition_update !=========================================================================== !*************************************************************************** @@ -639,27 +689,33 @@ end subroutine get_cp_dry_2hd ! !*************************************************************************** ! - subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) use cam_abortutils, only: endrun use string_utils, only: to_str ! Dummy arguments ! tracer: Tracer array + ! + ! factor not present then tracer must be dry mixing ratio + ! if factor present tracer*factor must be dry mixing ratio + ! real(kind_phys), intent(in) :: tracer(:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) + ! dp: if provided then tracer is mass not mixing ratio + real(kind_phys), optional, intent(in) :: factor(:,:) ! active_species_idx_dycore: array of indices for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:) ! Local variables integer :: qdx, itrac real(kind_phys) :: sum_species(SIZE(cp, 1), SIZE(cp, 2)) real(kind_phys) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2)) - real(kind_phys) :: factor(SIZE(cp, 1), SIZE(cp, 2)) + real(kind_phys) :: factor_local(SIZE(cp, 1), SIZE(cp, 2)) integer :: idx_local(thermodynamic_active_species_num) character(len=*), parameter :: subname = 'get_cp_1hd: ' @@ -675,28 +731,41 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) idx_local = thermodynamic_active_species_idx end if - if (present(dp_dry)) then - factor = 1.0_kind_phys / dp_dry + if (present(factor)) then + factor_local = factor else - factor = 1.0_kind_phys + factor_local = 1.0_kind_phys end if + sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) - sum_species(:,:) = sum_species(:,:) + & - (tracer(:,:,itrac) * factor(:,:)) + sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:)) end do - ! Get heat capacity at constant pressure (Cp) for dry air: - call get_cp_dry(tracer, idx_local, sum_cp, fact=factor) + if (dry_air_species_num == 0) then + ! **TMP** TODO CHECK hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + ! and dry_air_species_num is == 0 in CAM-SIMA as of 10/1/24 which means this clause will + ! change answers (?) -- previously in CAM-SIMA there was only a call to get_cp_dry + ! and not the two IF-clauses here. + sum_cp = thermodynamic_active_species_cp(0) + else if (present(cpdry)) then + ! + ! if cpdry is known don't recompute + ! + sum_cp = cpdry + else + ! Get heat capacity at constant pressure (Cp) for dry air: + call get_cp_dry(tracer, idx_local, sum_cp, fact=factor_local) + end if ! Add water species to Cp: do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) sum_cp(:,:) = sum_cp(:,:) + & - (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * & - factor(:,:)) + (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * factor_local(:,:)) end do + if (inv_cp) then cp = sum_species / sum_cp else @@ -707,7 +776,7 @@ end subroutine get_cp_1hd !=========================================================================== - subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_2hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) ! Version of get_cp for arrays that have a second horizontal index use cam_abortutils, only: endrun use string_utils, only: to_str @@ -715,14 +784,15 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) ! Dummy arguments ! tracer: Tracer array real(kind_phys), intent(in) :: tracer(:,:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:,:) + real(kind_phys), optional, intent(in) :: factor(:,:,:) ! active_species_idx_dycore: array of indicies for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:,:) ! Local variables integer :: jdx @@ -730,11 +800,17 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) character(len=*), parameter :: subname = 'get_cp_2hd: ' do jdx = 1, SIZE(cp, 2) - if (present(dp_dry)) then - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & - dp_dry=dp_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + if (present(factor).and.present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) + else if (present(factor)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + else if (present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) else - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& active_species_idx_dycore=active_species_idx_dycore) end if end do @@ -843,9 +919,10 @@ end subroutine get_R_dry_2hd ! !*************************************************************************** ! - subroutine get_R_1hd(tracer, active_species_idx, R, fact) + subroutine get_R_1hd(tracer, active_species_idx, R, fact, Rdry) use cam_abortutils, only: endrun use string_utils, only: to_str + use physconst, only: rair ! Dummy arguments ! tracer: !tracer array @@ -856,6 +933,7 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) real(kind_phys), intent(out) :: R(:, :) ! fact: optional factor for converting tracer to dry mixing ratio real(kind_phys), optional, intent(in) :: fact(:, :) + real(kind_phys), optional, intent(in) :: Rdry(:, :) ! Local variables integer :: qdx, itrac @@ -874,12 +952,19 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) call endrun(subname//"SIZE mismatch in dimension 2 "// & to_str(SIZE(fact, 2))//' /= '//to_str(SIZE(factor, 2))) end if - call get_R_dry(tracer, active_species_idx, R, fact=fact) factor = fact(:,:) else - call get_R_dry(tracer, active_species_idx, R) factor = 1.0_kind_phys end if + + if (dry_air_species_num == 0) then + R = rair + else if (present(Rdry)) then + R = Rdry + else + call get_R_dry(tracer, active_species_idx, R, fact=factor) + end if + idx_local = active_species_idx sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num @@ -934,7 +1019,7 @@ end subroutine get_R_2hd !************************************************************************************************************************* ! subroutine get_mbarv_1hd(tracer, active_species_idx, mbarv_in, fact) - use physconst, only: mwdry, rair, cpair + use physconst, only: mwdry real(kind_phys), intent(in) :: tracer(:,:,:) !tracer array integer, intent(in) :: active_species_idx(:) !index of active species in tracer real(kind_phys), intent(out) :: mbarv_in(:,:) !molecular weight of dry air diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 59dd1c83..09dff1f1 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -33,8 +33,10 @@ module cam_thermo ! cam_thermo_init: Initialize constituent dependent properties public :: cam_thermo_init - ! cam_thermo_update: Update constituent dependent properties - public :: cam_thermo_update + ! cam_thermo_dry_air_update: Update dry air composition dependent properties + public :: cam_thermo_dry_air_update + ! cam_thermo_water_update: Update water dependent properties + public :: cam_thermo_water_update ! get_enthalpy: enthalpy quantity = dp*cp*T public :: get_enthalpy ! get_virtual_temp: virtual temperature @@ -77,6 +79,7 @@ module cam_thermo ! mixing_ratio options integer, public, parameter :: DRY_MIXING_RATIO = 1 integer, public, parameter :: MASS_MIXING_RATIO = 2 + !> \section arg_table_cam_thermo Argument Table !! \htmlinclude cam_thermo.html !--------------- Variables below here are for WACCM-X --------------------- @@ -85,7 +88,7 @@ module cam_thermo ! kmcnd: molecular conductivity J m-1 s-1 K-1 real(kind_phys), public, protected, allocatable :: kmcnd(:,:) - !------------- Variables for consistent themodynamics -------------------- + !------------- Variables for consistent thermodynamics -------------------- ! ! @@ -208,51 +211,71 @@ end subroutine cam_thermo_init !=========================================================================== + ! !*************************************************************************** ! - ! cam_thermo_update: update species dependent constants for physics + ! cam_thermo_dry_air_update: update dry air species dependent constants for physics ! !*************************************************************************** ! - subroutine cam_thermo_update(mmr, T, ncol, update_thermo_variables, to_moist_factor) - use air_composition, only: air_composition_update, update_zvirv + subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) + use air_composition, only: dry_air_composition_update + use air_composition, only: update_zvirv use string_utils, only: to_str - !----------------------------------------------------------------------- - ! Update the physics "constants" that vary - !------------------------------------------------------------------------- - !------------------------------Arguments---------------------------------- + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) + real(kind_phys), intent(in) :: T(:,:) ! temperature + integer, intent(in) :: ncol ! number of columns + logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables + ! false: do not calculate composition-dependent thermo variables + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr moist convert - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array - real(kind_phys), intent(in) :: T(:,:) ! temperature - integer, intent(in) :: ncol ! number of columns - logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables - ! false: do not calculate composition-dependent thermo variables - - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - ! - !---------------------------Local storage------------------------------- - real(kind_phys):: sponge_factor(SIZE(mmr, 2)) - character(len=*), parameter :: subname = 'cam_thermo_update: ' + ! Local vars + real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) + character(len=*), parameter :: subname = 'cam_thermo_dry_air_update: ' if (.not. update_thermo_variables) then return end if - if (present(to_moist_factor)) then - if (SIZE(to_moist_factor, 1) /= ncol) then - call endrun(subname//'DIM 1 of to_moist_factor is'//to_str(SIZE(to_moist_factor,1))//'but should be'//to_str(ncol)) - end if + if (present(to_dry_factor)) then + if (SIZE(to_dry_factor, 1) /= ncol) then + call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) + end if end if + sponge_factor = 1.0_kind_phys - call air_composition_update(mmr, ncol, to_moist_factor=to_moist_factor) + call dry_air_composition_update(mmr, ncol, to_dry_factor=to_dry_factor) call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:), & - kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_moist_factor, & + kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_dry_factor, & active_species_idx_dycore=thermodynamic_active_species_idx) + + ! Calculate zvirv for WACCM-X. call update_zvirv() - end subroutine cam_thermo_update + end subroutine cam_thermo_dry_air_update + + ! + !*************************************************************************** + ! + ! cam_thermo_water_update: update water species dependent constants for physics + ! + !*************************************************************************** + ! + subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor) + use air_composition, only: water_composition_update + !----------------------------------------------------------------------- + ! Update the physics "constants" that vary + !------------------------------------------------------------------------- + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: energy_formula + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + call water_composition_update(mmr, ncol, energy_formula, to_dry_factor=to_dry_factor) + end subroutine cam_thermo_water_update !=========================================================================== @@ -1554,28 +1577,32 @@ end subroutine cam_thermo_calc_kappav_2hd ! !*************************************************************************** ! - subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & - vcoord, ps, phis, z_mid, dycore_idx, qidx, te, se, ke, & - wv, H2O, liq, ice) + subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & + cp_or_cv, U, V, T, vcoord, ptop, phis, z_mid, dycore_idx, qidx, & + te, se, po, ke, wv, H2O, liq, ice) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS use cam_logfile, only: iulog - use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure use air_composition, only: wv_idx - use physconst, only: gravit, latvap, latice + use air_composition, only: dry_air_species_num + use physconst, only: rga, latvap, latice ! Dummy arguments ! tracer: tracer mixing ratio + ! + ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry real(kind_phys), intent(in) :: tracer(:,:,:) + logical, intent(in) :: moist_mixing_ratio ! pdel: pressure level thickness - real(kind_phys), intent(in) :: pdel(:,:) + real(kind_phys), intent(in) :: pdel_in(:,:) ! cp_or_cv: dry air heat capacity under constant pressure or ! constant volume (depends on vcoord) real(kind_phys), intent(in) :: cp_or_cv(:,:) real(kind_phys), intent(in) :: U(:,:) real(kind_phys), intent(in) :: V(:,:) real(kind_phys), intent(in) :: T(:,:) - integer, intent(in) :: vcoord ! vertical coordinate - real(kind_phys), intent(in), optional :: ps(:) + integer, intent(in) :: vcoord !REMOVECAM - vcoord or energy formula to use + real(kind_phys), intent(in), optional :: ptop(:) real(kind_phys), intent(in), optional :: phis(:) real(kind_phys), intent(in), optional :: z_mid(:,:) ! dycore_idx: use dycore index for thermodynamic active species @@ -1588,8 +1615,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & real(kind_phys), intent(out), optional :: te (:) ! KE: vertically integrated kinetic energy real(kind_phys), intent(out), optional :: ke (:) - ! SE: vertically integrated internal+geopotential energy + ! SE: vertically integrated enthalpy (pressure coordinate) + ! or internal energy (z coordinate) real(kind_phys), intent(out), optional :: se (:) + ! PO: vertically integrated PHIS term (pressure coordinate) + ! or potential energy (z coordinate) + real(kind_phys), intent(out), optional :: po (:) ! WV: vertically integrated water vapor real(kind_phys), intent(out), optional :: wv (:) ! liq: vertically integrated liquid @@ -1599,10 +1630,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & ! Local variables real(kind_phys) :: ke_vint(SIZE(tracer, 1)) ! Vertical integral of KE - real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of SE + real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of enthalpy or internal energy + real(kind_phys) :: po_vint(SIZE(tracer, 1)) ! Vertical integral of PHIS or potential energy real(kind_phys) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv real(kind_phys) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq real(kind_phys) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice + real(kind_phys) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) ! moist pressure level thickness real(kind_phys) :: latsub ! latent heat of sublimation integer :: ierr @@ -1649,78 +1682,96 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & wvidx = wv_idx end if + if (moist_mixing_ratio) then + pdel = pdel_in + else + pdel = pdel_in + do qdx = dry_air_species_num+1, thermodynamic_active_species_num + pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx)) + end do + end if + + ke_vint = 0._kind_phys + se_vint = 0._kind_phys select case (vcoord) - case(vc_moist_pressure, vc_dry_pressure) - if ((.not. present(ps)) .or. (.not. present(phis))) then - write(iulog, *) subname, ' ps and phis must be present for ', & - 'moist/dry pressure vertical coordinate' - call endrun(subname//': ps and phis must be present for '// & - 'moist/dry pressure vertical coordinate') + case(ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE) + if (.not. present(ptop).or. (.not. present(phis))) then + write(iulog, *) subname, ' ptop and phis must be present for ', & + 'FV/SE energy formula' + call endrun(subname//': ptop and phis must be present for '// & + 'FV/SE energy formula') end if - ke_vint = 0._kind_phys - se_vint = 0._kind_phys - wv_vint = 0._kind_phys + po_vint = ptop do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * & - 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit) + 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2)) * rga se_vint(idx) = se_vint(idx) + (T(idx, kdx) * & - cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit) - wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & - pdel(idx, kdx) / gravit) + cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga) + po_vint(idx) = po_vint(idx)+pdel(idx, kdx) + end do end do do idx = 1, SIZE(tracer, 1) - se_vint(idx) = se_vint(idx) + (phis(idx) * ps(idx) / gravit) + po_vint(idx) = (phis(idx) * po_vint(idx) * rga) end do - case(vc_height) - if (.not. present(z_mid)) then - write(iulog, *) subname, & - ' z_mid must be present for height vertical coordinate' - call endrun(subname//': z_mid must be present for height '// & - 'vertical coordinate') + case(ENERGY_FORMULA_DYCORE_MPAS) + if (.not. present(phis)) then + write(iulog, *) subname, ' phis must be present for ', & + 'MPAS energy formula' + call endrun(subname//': phis must be present for '// & + 'MPAS energy formula') end if - ke_vint = 0._kind_phys - se_vint = 0._kind_phys - wv_vint = 0._kind_phys + po_vint = 0._kind_phys do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * & - 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit) + 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) * rga) se_vint(idx) = se_vint(idx) + (T(idx, kdx) * & - cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit) + cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga) ! z_mid is height above ground - se_vint(idx) = se_vint(idx) + (z_mid(idx, kdx) + & - phis(idx) / gravit) * pdel(idx, kdx) - wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & - pdel(idx, kdx) / gravit) + po_vint(idx) = po_vint(idx) + (z_mid(idx, kdx) + & + phis(idx) * rga) * pdel(idx, kdx) end do end do case default - write(iulog, *) subname, ' vertical coordinate not supported: ', vcoord - call endrun(subname//': vertical coordinate not supported') + write(iulog, *) subname, ' energy formula not supported: ', vcoord + call endrun(subname//': energy formula not supported') end select if (present(te)) then - te = se_vint + ke_vint + te = se_vint + po_vint + ke_vint end if if (present(se)) then se = se_vint end if + if (present(po)) then + po = po_vint + end if if (present(ke)) then ke = ke_vint end if - if (present(wv)) then - wv = wv_vint - end if ! ! vertical integral of total liquid water ! + if (.not.moist_mixing_ratio) then + pdel = pdel_in! set pseudo density to dry + end if + + wv_vint = 0._kind_phys + do kdx = 1, SIZE(tracer, 2) + do idx = 1, SIZE(tracer, 1) + wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & + pdel(idx, kdx) * rga) + end do + end do + if (present(wv)) wv = wv_vint + liq_vint = 0._kind_phys do qdx = 1, thermodynamic_active_species_liq_num do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) - liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * & - tracer(idx, kdx, species_liq_idx(qdx)) / gravit) + liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * & + tracer(idx, kdx, species_liq_idx(qdx)) * rga) end do end do end do @@ -1734,7 +1785,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ice_vint(idx) = ice_vint(idx) + (pdel(idx, kdx) * & - tracer(idx, kdx, species_ice_idx(qdx)) / gravit) + tracer(idx, kdx, species_ice_idx(qdx)) * rga) end do end do end do @@ -1762,7 +1813,6 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & end select end if deallocate(species_idx, species_liq_idx, species_ice_idx) - end subroutine get_hydrostatic_energy_1hd !=========================================================================== diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90 new file mode 100644 index 00000000..c3e2c15b --- /dev/null +++ b/src/data/cam_thermo_formula.F90 @@ -0,0 +1,37 @@ +module cam_thermo_formula + + implicit none + private + save + + ! saves energy formula to use for physics and dynamical core + ! for use in cam_thermo, air_composition and other modules + ! separated into cam_thermo_formula module for clean dependencies + + ! energy_formula options (was vcoord in CAM and stored in dyn_tests_utils) + integer, public, parameter :: ENERGY_FORMULA_DYCORE_FV = 0 ! vc_moist_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_SE = 1 ! vc_dry_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_MPAS = 2 ! vc_height + + !> \section arg_table_cam_thermo_formula Argument Table + !! \htmlinclude cam_thermo_formula.html + ! energy_formula_dycore: energy formula used for dynamical core + ! written by the dynamical core + integer, public :: energy_formula_dycore + ! energy_formula_physics: energy formula used for physics + integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + + ! Public subroutines + public :: cam_thermo_formula_init + +contains + subroutine cam_thermo_formula_init() + use phys_vars_init_check, only: mark_as_initialized + + ! Physics energy formulation is always FV (moist pressure coordinate) + energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + call mark_as_initialized("total_energy_formula_for_physics") + + end subroutine cam_thermo_formula_init + +end module cam_thermo_formula diff --git a/src/data/cam_thermo_formula.meta b/src/data/cam_thermo_formula.meta new file mode 100644 index 00000000..f8bf04a1 --- /dev/null +++ b/src/data/cam_thermo_formula.meta @@ -0,0 +1,17 @@ +[ccpp-table-properties] + name = cam_thermo_formula + type = module + +[ccpp-arg-table] + name = cam_thermo_formula + type = module +[ energy_formula_dycore ] + standard_name = total_energy_formula_for_dycore + units = 1 + type = integer + dimensions = () +[ energy_formula_physics ] + standard_name = total_energy_formula_for_physics + units = 1 + type = integer + dimensions = () diff --git a/src/data/registry.xml b/src/data/registry.xml index 41b31416..d2445ab9 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -13,6 +13,7 @@ $SRCROOT/src/physics/utils/tropopause_climo_read.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta + $SRCROOT/src/data/cam_thermo_formula.meta $SRCROOT/src/data/ref_pres.meta $SRCROOT/src/dynamics/utils/vert_coord.meta $SRCROOT/src/dynamics/utils/hycoef.meta @@ -251,6 +252,12 @@ horizontal_dimension tw_cur state_tw_cur + + + flag indicating if is first timestep of initial run + horizontal_dimension vertical_layer_dimension zvir + + Enthalpy or internal energy scaling factor for energy consistency + horizontal_dimension vertical_layer_dimension + cp_or_cv_dycore + and set the energy formulation + ! (also called vc_dycore in CAM) + + type(file_desc_t), intent(inout) :: file + logical, intent(in) :: grid_is_latlon + + ! Local variables + integer :: ierr, xtype + character(len=*), parameter :: subname = 'find_energy_formula' + + energy_formula_dycore = -1 + + ! Is FV dycore? (has lat lon dimension) + if(grid_is_latlon) then + energy_formula_dycore = ENERGY_FORMULA_DYCORE_FV + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use FV dycore energy formula' + endif + else + ! Is SE dycore? + ierr = pio_inq_att(file, pio_global, 'ne', xtype) + if (ierr == PIO_NOERR) then + ! Has ne property - is SE dycore. + ! if has fv_nphys then is physics grid (ne..pg..), but the energy formulation is the same. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use SE dycore energy formula' + endif + else + ! Is unstructured and is MPAS dycore + ! there are no global attributes to identify MPAS dycore, so this has to do for now. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use MPAS dycore energy formula' + endif + endif + endif + + if(energy_formula_dycore /= -1) then + call mark_as_initialized("total_energy_formula_for_dycore") + endif + + end subroutine + end module dyn_grid diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 8b56e9d9..22f8d6b8 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -582,18 +582,22 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) use cam_constituents, only: const_qmin use runtime_obj, only: wv_stdname use physics_types, only: lagrangian_vertical - use physconst, only: cpair, gravit, zvir, cappa - use cam_thermo, only: cam_thermo_update - use physics_types, only: cpairv, rairv, zvirv + use physconst, only: cpair, gravit, zvir + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use air_composition, only: thermodynamic_active_species_num + use air_composition, only: thermodynamic_active_species_idx + use air_composition, only: dry_air_species_num + use physics_types, only: cpairv, rairv, zvirv, cappav use physics_grid, only: columns_on_task use geopotential_temp, only: geopotential_temp_run use static_energy, only: update_dry_static_energy_run use qneg, only: qneg_run -! use check_energy, only: check_energy_timestep_init +! use check_energy_chng, only: check_energy_chng_timestep_init use hycoef, only: hyai, ps0 use shr_vmath_mod, only: shr_vmath_log use shr_kind_mod, only: shr_kind_cx use dyn_comp, only: ixo, ixo2, ixh, ixh2 + use cam_thermo_formula,only: ENERGY_FORMULA_DYCORE_SE ! arguments type(runtime_options), intent(in) :: cam_runtime_opts ! Runtime settings object @@ -607,7 +611,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !constituent properties pointer type(ccpp_constituent_prop_ptr_t), pointer :: const_prop_ptr(:) - integer :: m, i, k + integer :: m, i, k, m_cnst integer :: ix_q !Needed for "geopotential_temp" CCPP scheme @@ -622,6 +626,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) real(r8) :: mmrSum_O_O2_H ! Sum of mass mixing ratios for O, O2, and H real(r8), parameter :: mmrMin=1.e-20_r8 ! lower limit of o2, o, and h mixing ratios real(r8), parameter :: N2mmrMin=1.e-6_r8 ! lower limit of o2, o, and h mixing ratios + real(r8), parameter :: H2lim=6.e-5_r8 ! H2 limiter: 10x global H2 MMR (Roble, 1995) !---------------------------------------------------------------------------- ! Nullify pointers @@ -683,14 +688,23 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end do ! wet pressure variables (should be removed from physics!) + factor_array(:,:) = 1.0_kind_phys + !$omp parallel do num_threads(horz_num_threads) private (k, i, m_cnst) + do m_cnst = dry_air_species_num + 1, thermodynamic_active_species_num + ! include all water species in the factor array. + m = thermodynamic_active_species_idx(m_cnst) + do k = 1, nlev + do i = 1, pcols + ! at this point all q's are dry + factor_array(i,k) = factor_array(i,k) + const_data_ptr(i,k,m) + end do + end do + end do !$omp parallel do num_threads(horz_num_threads) private (k, i) - do k=1,nlev - do i=1, pcols - ! to be consistent with total energy formula in physic's check_energy module only - ! include water vapor in moist dp - factor_array(i,k) = 1._kind_phys+const_data_ptr(i,k,ix_q) - phys_state%pdel(i,k) = phys_state%pdeldry(i,k)*factor_array(i,k) + do k = 1, nlev + do i = 1, pcols + phys_state%pdel(i,k) = phys_state%pdeldry(i,k) * factor_array(i,k) end do end do @@ -721,31 +735,12 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) do k = 1, nlev do i = 1, pcols phys_state%rpdel(i,k) = 1._kind_phys/phys_state%pdel(i,k) - phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa - end do - end do - - ! all tracers (including moisture) are in dry mixing ratio units - ! physics expect water variables moist - factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) - - !$omp parallel do num_threads(horz_num_threads) private (m, k, i) - do m=1, num_advected - do k = 1, nlev - do i=1, pcols - !This should ideally check if a constituent is a wet - !mixing ratio or not, but until that is working properly - !in the CCPP framework just check for the water species status - !instead, which is all that CAM physics requires: - if (const_is_water_species(m)) then - const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) - end if - end do end do end do !------------------------------------------------------------ - ! Ensure O2 + O + H (N2) mmr greater than one. + ! Apply limiters to mixing ratios of major species (WACCMX): + ! Ensure N2 = 1 - (O2 + O + H) mmr is greater than 0 ! Check for unusually large H2 values and set to lower value. !------------------------------------------------------------ if (cam_runtime_opts%waccmx_option() == 'ionosphere' .or. & @@ -769,8 +764,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) endif - if(const_data_ptr(i,k,ixh2) .gt. 6.e-5_r8) then - const_data_ptr(i,k,ixh2) = 6.e-5_r8 + if(const_data_ptr(i,k,ixh2) > H2lim) then + const_data_ptr(i,k,ixh2) = H2lim endif end do @@ -789,11 +784,62 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! Compute molecular viscosity(kmvis) and conductivity(kmcnd). ! Update zvirv registry variable; calculated for WACCM-X. !----------------------------------------------------------------------------- + ! **TEMP** TODO CHECK hplin: CAM has this if-clause for dry_air_species_num > 0 + ! or otherwise uses zvirv = zvir. CAM-SIMA previously did not have this, and + ! instead has a switch for update_thermodynamic_variables. Check if we still want + ! this if-clause or change it to something else. + if (dry_air_species_num > 0) then + call cam_thermo_dry_air_update( & + mmr = const_data_ptr, & ! dry MMR + T = phys_state%t, & + ncol = pcols, & + update_thermo_variables = cam_runtime_opts%update_thermodynamic_variables() & + ) + else + zvirv(:,:) = zvir + end if + + ! + ! update cp_or_cv_dycore in module air_composition. + ! (note: at this point q is dry) + ! + call cam_thermo_water_update( & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + energy_formula = ENERGY_FORMULA_DYCORE_SE & + ) - call cam_thermo_update(const_data_ptr, phys_state%t, pcols, & - cam_runtime_opts%update_thermodynamic_variables()) + !$omp parallel do num_threads(horz_num_threads) private (k, i) + do k = 1, nlev + do i = 1, pcols + phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) + end do + end do + + ! ========= Q is dry ^^^ ---- Q is moist vvv ========= ! + + ! + ! CAM physics expects that: water tracers (including moisture) are moist; the rest dry mixing ratio + ! at this point Q is converted to moist. + ! + factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) + + !$omp parallel do num_threads(horz_num_threads) private (m, k, i) + do m = 1, num_advected + do k = 1, nlev + do i = 1, pcols + ! This should ideally check if a constituent is a wet + ! mixing ratio or not, but until that is working properly + ! in the CCPP framework just check for the water species status + ! instead, which is all that CAM physics requires: + if (const_is_water_species(m)) then + const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) + end if + end do + end do + end do - !Call geopotential_temp CCPP scheme: + ! Call geopotential_temp CCPP scheme: call geopotential_temp_run(pver, lagrangian_vertical, pver, 1, & pverp, 1, num_advected, phys_state%lnpint, & phys_state%pint, phys_state%pmid, phys_state%pdel, & @@ -810,7 +856,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !Remove once check_energy scheme exists in CAMDEN: #if 0 ! Compute energy and water integrals of input state - call check_energy_timestep_init(phys_state, phys_tend, pbuf_chnk) + ! call check_energy_chng_timestep_init(...) #endif end subroutine derived_phys_dry diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90 index 1341b9f4..cf88d7f9 100644 --- a/src/dynamics/se/dycore/prim_advance_mod.F90 +++ b/src/dynamics/se/dycore/prim_advance_mod.F90 @@ -1674,8 +1674,12 @@ subroutine calc_tot_energy_dynamics(elem,fvm,nets,nete,tl,tl_qdp,outfld_name_suf call get_dp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),MASS_MIXING_RATIO,thermodynamic_active_species_idx_dycore,& elem(ie)%state%dp3d(:,:,:,tl),pdel,ps=ps,ptop=hyai(1)*ps0) call get_cp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),& - .false.,cp,dp_dry=elem(ie)%state%dp3d(:,:,:,tl),& + .false.,cp,factor=1.0_r8/elem(ie)%state%dp3d(:,:,:,tl),& active_species_idx_dycore=thermodynamic_active_species_idx_dycore) + + ! TODO: need to port cam6_3_109 changes to total energy using get_hydrostatic_energy + ! https://github.com/ESCOMP/CAM/pull/761/files#diff-946bde17289e2f42e43e64413610aa11d102deda8b5199ddaa5b71e67e5d517a + do k = 1, nlev do j=1,np do i = 1, np diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index ab52d91c..5f8dd2e6 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -566,6 +566,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use dyn_thermo, only: get_molecular_diff_coef_reference !use cam_history, only: addfld, add_default, horiz_only, register_vector_field use gravity_waves_sources, only: gws_init + use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_SE !SE dycore: use prim_advance_mod, only: prim_advance_init @@ -642,6 +643,8 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) real(r8) :: tau0, krange, otau0, scale real(r8) :: km_sponge_factor_local(nlev+1) !---------------------------------------------------------------------------- + ! Set dynamical core energy formula for use in cam_thermo. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE ! Now allocate and set condenstate vars allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers diff --git a/src/dynamics/utils/dyn_thermo.F90 b/src/dynamics/utils/dyn_thermo.F90 index c4c4723c..9b465031 100644 --- a/src/dynamics/utils/dyn_thermo.F90 +++ b/src/dynamics/utils/dyn_thermo.F90 @@ -61,7 +61,7 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Declare local variables: real(kind_phys), allocatable :: tracer_phys(:,:,:,:) real(kind_phys), allocatable :: cp_phys(:,:,:) - real(kind_phys), allocatable :: dp_dry_phys(:,:,:) + real(kind_phys), allocatable :: factor_phys(:,:,:) !check_allocate variables: integer :: iret !allocate status integer @@ -70,11 +70,16 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Check if kinds are different: if (kind_phys == kind_dyn) then - !The dynamics and physics kind is the same, so just call the physics - !routine directly: - call get_cp_phys(tracer,inv_cp,cp, & - dp_dry=dp_dry, & - active_species_idx_dycore=active_species_idx_dycore) + ! The dynamics and physics kind is the same, so just call the physics + ! routine directly: + if(present(dp_dry)) then + call get_cp_phys(tracer,inv_cp,cp, & + factor=1.0_kind_phys/dp_dry, & + active_species_idx_dycore=active_species_idx_dycore) + else + call get_cp_phys(tracer,inv_cp,cp, & + active_species_idx_dycore=active_species_idx_dycore) + endif else @@ -95,18 +100,18 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Allocate and set optional variables: if (present(dp_dry)) then - allocate(dp_dry_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret) + allocate(factor_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret) call check_allocate(iret, subname, & - 'dp_dry_phys', & + 'factor_phys', & file=__FILE__, line=__LINE__) !Set optional local variable: - dp_dry_phys = real(dp_dry, kind_phys) + factor_phys = 1.0_kind_phys/real(dp_dry, kind_phys) end if - !Call physics routine using local vriables with matching kinds: + !Call physics routine using local variables with matching kinds: call get_cp_phys(tracer_phys,inv_cp,cp_phys, & - dp_dry=dp_dry_phys, & + factor=factor_phys, & active_species_idx_dycore=active_species_idx_dycore) !Set output variables back to dynamics kind: @@ -116,8 +121,8 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) deallocate(tracer_phys) deallocate(cp_phys) - if (allocated(dp_dry_phys)) then - deallocate(dp_dry_phys) + if (allocated(factor_phys)) then + deallocate(factor_phys) end if @@ -957,7 +962,7 @@ subroutine get_rho_dry(tracer,temp,ptop,dp_dry,tracer_mass,& end if - !Call physics routine using local vriables with matching kinds: + !Call physics routine using local variables with matching kinds: call get_rho_dry_phys(tracer_phys,temp_phys, & ptop_phys, dp_dry_phys,tracer_mass, & rho_dry=rho_dry_phys, & diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index e95c172d..952ebddf 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit e95c172d7a5a0ebf054f420b08416228e211baa3 +Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 diff --git a/src/physics/utils/phys_comp.F90 b/src/physics/utils/phys_comp.F90 index abe0428a..59ec39ab 100644 --- a/src/physics/utils/phys_comp.F90 +++ b/src/physics/utils/phys_comp.F90 @@ -134,6 +134,7 @@ subroutine phys_init() use physics_grid, only: columns_on_task use vert_coord, only: pver, pverp use cam_thermo, only: cam_thermo_init + use cam_thermo_formula, only: cam_thermo_formula_init use physics_types, only: allocate_physics_types_fields use cam_ccpp_cap, only: cam_ccpp_physics_initialize use cam_ccpp_cap, only: ccpp_physics_suite_part_list @@ -142,6 +143,7 @@ subroutine phys_init() integer :: i_group call cam_thermo_init(columns_on_task, pver, pverp) + call cam_thermo_formula_init() call allocate_physics_types_fields(columns_on_task, pver, pverp, & set_init_val_in=.true., reallocate_in=.false.) diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90 index 39fc0f99..cf28f98d 100644 --- a/src/physics/utils/physics_grid.F90 +++ b/src/physics/utils/physics_grid.F90 @@ -39,8 +39,10 @@ module physics_grid public :: get_rlat_p ! latitude of a physics column in radians public :: get_rlon_p ! longitude of a physics column in radians public :: get_area_p ! area of a physics column in radians squared + public :: get_wght_p ! weight of a physics column in radians squared public :: get_rlat_all_p ! latitudes of physics cols on task (radians) public :: get_rlon_all_p ! longitudes of physics cols on task (radians) + public :: get_wght_all_p ! weights of physics cols on task public :: get_dyn_col_p ! dynamics local blk number and blk offset(s) public :: global_index_p ! global column index of a physics column public :: local_index_p ! local column index of a physics column @@ -475,6 +477,26 @@ end function get_area_p !======================================================================== + real(r8) function get_wght_p(index) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + ! weight of a physics column in radians squared + + ! Dummy argument + integer, intent(in) :: index + ! Local variables + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_p' + + ! Check that input is valid: + call check_phys_input(subname, index) + + get_wght_p = phys_columns(index)%weight + + end function get_wght_p + + !======================================================================== + subroutine get_rlat_all_p(rlatdim, rlats) use cam_logfile, only: iulog use cam_abortutils, only: endrun @@ -535,6 +557,36 @@ end subroutine get_rlon_all_p !======================================================================== + subroutine get_wght_all_p(wghtdim, wghts) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + !----------------------------------------------------------------------- + ! + ! get_wght_all_p: Return all weights on task. + ! + !----------------------------------------------------------------------- + ! Dummy Arguments + integer, intent(in) :: wghtdim ! declared size of output array + real(r8), intent(out) :: wghts(wghtdim) ! array of weights + + ! Local variables + integer :: index ! loop index + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_all_p: ' + + !----------------------------------------------------------------------- + + ! Check that input is valid: + call check_phys_input(subname, wghtdim) + + do index = 1, wghtdim + wghts(index) = phys_columns(index)%weight + end do + + end subroutine get_wght_all_p + + !======================================================================== + subroutine get_dyn_col_p(index, blk_num, blk_ind) use cam_logfile, only: iulog use cam_abortutils, only: endrun diff --git a/src/utils/gmean_mod.F90 b/src/utils/gmean_mod.F90 new file mode 100644 index 00000000..ca6dd6d4 --- /dev/null +++ b/src/utils/gmean_mod.F90 @@ -0,0 +1,256 @@ +module gmean_mod + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Perform global mean calculations for energy conservation and other checks. + ! + ! Method: + ! Reproducible (scalable): + ! Convert to fixed point (integer representation) to enable + ! reproducibility when using MPI collectives. + ! If error checking is on (via setting reprosum_diffmax > 0 and + ! reprosum_recompute = .true. in user_nl_cpl), shr_reprosum_calc will + ! check the accuracy of its computation with a fast but + ! non-reproducible algorithm. If any error is reported, report + ! the difference and the expected sum and abort run (call endrun) + ! + ! gmean_mod in to_be_ccppized is different from the CAM version and + ! has chunk support removed. + ! + ! + !----------------------------------------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use physics_grid, only: pcols => columns_on_task + use perf_mod, only: t_startf, t_stopf + use cam_logfile, only: iulog + + implicit none + private + + public :: gmean ! compute global mean of 2D fields on physics decomposition + + interface gmean + module procedure gmean_arr + module procedure gmean_scl + end interface gmean + + private :: gmean_fixed_repro + private :: gmean_float_norepro + + ! Set do_gmean_tests to .true. to run a gmean challenge test + logical, private :: do_gmean_tests = .false. + +CONTAINS + + ! + !======================================================================== + ! + + subroutine gmean_arr (arr, arr_gmean, nflds) + use shr_strconvert_mod, only: toString + use spmd_utils, only: masterproc + use cam_abortutils, only: endrun + use shr_reprosum_mod, only: shr_reprosum_reldiffmax, shr_reprosum_recompute, shr_reprosum_tolExceeded + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! + ! Method is to call shr_reprosum_calc (called from gmean_fixed_repro) + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols, nflds) + real(r8), intent(out) :: arr_gmean(nflds) ! global means + ! + ! Local workspace + ! + real(r8) :: rel_diff(2, nflds) + integer :: ifld ! field index + integer :: num_err + logical :: write_warning + ! + !----------------------------------------------------------------------- + ! + call t_startf('gmean_arr') + call t_startf ('gmean_fixed_repro') + call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds) + call t_stopf ('gmean_fixed_repro') + + ! check that "fast" reproducible sum is accurate enough. If not, calculate + ! using old method + write_warning = masterproc + num_err = 0 + if (shr_reprosum_tolExceeded('gmean', nflds, write_warning, & + iulog, rel_diff)) then + if (shr_reprosum_recompute) then + do ifld = 1, nflds + if (rel_diff(1, ifld) > shr_reprosum_reldiffmax) then + call gmean_float_norepro(arr(:,ifld), arr_gmean(ifld), ifld) + num_err = num_err + 1 + end if + end do + end if + end if + call t_stopf('gmean_arr') + if (num_err > 0) then + call endrun('gmean: '//toString(num_err)//' reprosum errors found') + end if + + end subroutine gmean_arr + + ! + !======================================================================== + ! + + subroutine gmean_scl (arr, gmean) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of a single field in "arr" in the physics grid + ! + !----------------------------------------------------------------------- + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + ! Input array, chunked + real(r8), intent(out):: gmean ! global means + ! + ! Local workspace + ! + integer, parameter :: nflds = 1 + real(r8) :: gmean_array(nflds) + real(r8) :: array(pcols, nflds) + integer :: ncols, lchnk + + array(:ncols, 1) = arr(:ncols) + call gmean_arr(array, gmean_array, nflds) + gmean = gmean_array(1) + + end subroutine gmean_scl + + ! + !======================================================================== + ! + + subroutine gmean_float_norepro(arr, repro_sum, index) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of in the physics chunked + ! decomposition using a fast but non-reproducible algorithm. + ! Log that value along with the value computed by + ! shr_reprosum_calc () + ! + !----------------------------------------------------------------------- + + use physconst, only: pi + use spmd_utils, only: masterproc, masterprocid, mpicom + use mpi, only: mpi_real8, mpi_sum + use physics_grid, only: get_wght_p + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + real(r8), intent(in) :: repro_sum ! Value computed by reprosum + integer, intent(in) :: index ! Index of field in original call + ! + ! Local workspace + ! + integer :: icol + integer :: ierr + real(r8) :: wght + real(r8) :: check + real(r8) :: check_sum + real(r8), parameter :: pi4 = 4.0_r8 * pi + + ! + !----------------------------------------------------------------------- + ! + ! Calculate and print out non-reproducible value + check = 0.0_r8 + do icol = 1, pcols + wght = get_wght_p(icol) + check = check + arr(icol) * wght + end do + call MPI_reduce(check, check_sum, 1, mpi_real8, mpi_sum, & + masterprocid, mpicom, ierr) + + ! normalization + check_sum = check_sum / pi4 + + if (masterproc) then + write(iulog, '(a,i0,2(a,e20.13e2))') 'gmean(', index, ') = ', & + check_sum, ', reprosum reported ', repro_sum + end if + + end subroutine gmean_float_norepro + + ! + !======================================================================== + ! + subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! with a reproducible yet scalable implementation + ! based on a fixed-point algorithm. + ! + !----------------------------------------------------------------------- + use spmd_utils, only: mpicom + use physics_grid, only: get_wght_all_p + use physics_grid, only: ngcols_p => num_global_phys_cols + use physconst, only: pi + use shr_reprosum_mod, only: shr_reprosum_calc + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols,nflds) + ! arr_gmean: output global sums + real(r8), intent(out) :: arr_gmean(nflds) + ! rel_diff: relative and absolute differences from shr_reprosum_calc + real(r8), intent(out) :: rel_diff(2, nflds) + ! + ! Local workspace + ! + integer :: icol, ifld ! column, field indices + + real(r8) :: wght(pcols) ! integration weights + real(r8), allocatable :: xfld(:,:) ! weighted summands + ! + !----------------------------------------------------------------------- + ! + allocate(xfld(pcols, nflds)) + + ! pre-weight summands + call get_wght_all_p(pcols, wght) + + do ifld = 1, nflds + do icol = 1, pcols + xfld(icol, ifld) = arr(icol, ifld) * wght(icol) + end do + end do + + ! call fixed-point algorithm + call shr_reprosum_calc ( & + arr = xfld, & + arr_gsum = arr_gmean, & + nsummands = pcols, & ! # of local summands + dsummands = pcols, & ! declared first dimension of arr. + nflds = nflds, & + commid = mpicom, & + rel_diff = rel_diff & + ) + + deallocate(xfld) + ! final normalization + arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi) + + end subroutine gmean_fixed_repro + +end module gmean_mod diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 0ac72a21..78dc69f2 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -196,6 +196,11 @@ state_tw_cur + + + cp_or_cv_dycore + + RHO