From 37fdbfbf58cdebbb982a8c23857267e837c39e47 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 16 Dec 2024 11:09:31 -0500 Subject: [PATCH 1/2] Implements cam_thermo_water_update and CCPPized check_energy (#316) Originator(s): jimmielin Description (include the issue title, and the keyword ['closes', 'fixes', 'resolves'] followed by the issue number): All changes are bit-for-bit, except those noted: Implements `cam_thermo_water_update`: - updates `cp_or_cv_dycore` (`specific_heat_for_air_used_in_dycore`) in `air_composition.F90` for SE and MPAS dynamical cores - read `cp_or_cv_dycore` from CAM snapshot (refer to companion CAM PR) - energy formula (that has to be matching dycore) is recognized and set by null-dycore in `dyn_grid.F90` by looking at global attributes of the initial file; see `find_energy_formula` - added `is_first_timestep` logical state flag Ports `cam_thermo` and related updates in `air_composition` and `dp_coupling` from CAM 6.3.109 (https://github.com/ESCOMP/CAM/pull/761/files) - update to hydrostatic energy calculation - changes `get_cp`, `get_R` in `air_composition.F90` to use moist mixing ratios - **answer-changing:** update to moist-to-dry (for physics) conversion in `dp_coupling::derived_phys_dry` to account for all water tracers instead of just Q - **answer-changing:** update to not-really-"exner" calculation to use composition-dependent `cappav` instead of `cappa` Changes `vcoord` in `dyn_tests_utils` (old CAM) to `energy_formula` now in `cam_thermo_formula` (separated out into a different file to avoid dependency issues) - `vc_moist_pressure` is now `ENERGY_FORMULA_DYCORE_FV`; `vc_dry_pressure` is `_SE`; `vc_height` is `_MPAS` - these are just integer flags (0,1,2) and values are kept consistent with old CAM and their use in dynamics tests Ports global mean utility module (`gmean_mod.F90`), de-chunkized from CAM: - Implements `get_wght` in `physics_grid` for weighted sum calculation Imports `check_energy_chng` and `check_energy_fix` from `atmospheric_physics` Describe any changes made to build system: N/A Describe any changes made to the namelist: contained within ncar-physics List any changes to the defaults for the input datasets (e.g. boundary datasets): - Added `cp_or_cv_dycore` in CAM snapshots List all files eliminated and why: N/A List all files added and what they do: ``` - energy_formula A src/data/cam_thermo_formula.F90 A src/data/cam_thermo_formula.meta - gmean A src/utils/gmean_mod.F90 ``` List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) ``` - ncar-physics update M .gitmodules M src/physics/ncar_ccpp - `is_first_timestep` M src/control/cam_comp.F90 - cam_thermo_water_update for cp_or_cv_dycore (include in registry; read from ic) M src/data/air_composition.F90 M src/data/cam_thermo.F90 M src/data/registry.xml M tools/stdnames_to_inputnames_dictionary.xml M src/dynamics/se/dp_coupling.F90 M src/dynamics/se/dycore/prim_advance_mod.F90 M src/dynamics/se/dyn_comp.F90 M src/dynamics/utils/dyn_thermo.F90 - energy_formula M src/physics/utils/phys_comp.F90 M src/dynamics/mpas/dyn_comp.F90 M src/dynamics/none/dyn_comp.F90 M src/dynamics/none/dyn_grid.F90 - gmean M src/physics/utils/physics_grid.F90 ``` Note: bit-for-bit in check_energy with CAM is tricky to validate without dycore updates to SE; may need to merge #301 first --------- Co-authored-by: Kuan-Chih Wang --- cime_config/config_component.xml | 3 +- src/control/cam_comp.F90 | 23 +- src/data/air_composition.F90 | 144 ++++++++--- src/data/cam_thermo.F90 | 232 +++++++++++------ src/data/cam_thermo_formula.F90 | 39 +++ src/data/cam_thermo_formula.meta | 17 ++ src/data/registry.xml | 51 +++- src/dynamics/mpas/dyn_comp.F90 | 25 ++ src/dynamics/mpas/dyn_coupling.F90 | 22 +- src/dynamics/none/dyn_comp.F90 | 3 + src/dynamics/none/dyn_grid.F90 | 81 ++++++ src/dynamics/se/dp_coupling.F90 | 125 ++++++---- src/dynamics/se/dycore/prim_advance_mod.F90 | 6 +- src/dynamics/se/dyn_comp.F90 | 22 ++ src/dynamics/utils/dyn_thermo.F90 | 33 +-- src/physics/utils/phys_comp.F90 | 2 + src/physics/utils/physics_grid.F90 | 67 +++-- src/utils/gmean_mod.F90 | 263 ++++++++++++++++++++ test/existing-test-failures.txt | 13 +- tools/stdnames_to_inputnames_dictionary.xml | 17 +- 20 files changed, 955 insertions(+), 233 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/cime_config/config_component.xml b/cime_config/config_component.xml index 5805131c..030c5f87 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -161,8 +161,7 @@ -nlev 145 --> - + --physics-suites adiabatic --physics-suites tj2016 --analytic_ic --physics-suites kessler --analytic_ic diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 3daf0168..e3a178f8 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, nstep use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -149,9 +150,6 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & character(len=cx) :: errmsg !----------------------------------------------------------------------- - dtime_phys = 0.0_r8 - call mark_as_initialized('timestep_for_physics') - call init_pio_subsystem() ! Initializations using data passed from coupler. @@ -167,12 +165,20 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & call cam_ctrl_set_orbit(eccen, obliqr, lambm0, mvelpp) - call timemgr_init( & dtime, calendar, start_ymd, start_tod, ref_ymd, & ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod, & perpetual_run, perpetual_ymd, initial_run_in) + dtime_phys = 0.0_r8 + call mark_as_initialized('timestep_for_physics') + + is_first_timestep = .true. + call mark_as_initialized('is_first_timestep') + + nstep = get_nstep() + call mark_as_initialized('current_timestep_number') + ! Get current fractional calendar day. Needs to be updated at every timestep. calday = get_curr_calday() call mark_as_initialized('fractional_calendar_days_on_end_of_current_timestep') @@ -268,6 +274,10 @@ 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() + nstep = get_nstep() + !---------------------------------------------------------- ! First phase of dynamics (at least couple from dynamics to physics) ! Return time-step for physics from dynamics. @@ -514,10 +524,6 @@ subroutine cam_final(cam_out, cam_in) type(cam_out_t), pointer :: cam_out ! Output from CAM to surface type(cam_in_t), pointer :: cam_in ! Input from merged surface to CAM - ! - ! Local variable - ! - integer :: nstep ! Current timestep number. !----------------------------------------------------------------------- call phys_final() @@ -540,7 +546,6 @@ subroutine cam_final(cam_out, cam_in) call shr_sys_flush( iulog ) ! Flush all output to the CAM log file if (masterproc) then - nstep = get_nstep() write(iulog,9300) nstep-1,nstep 9300 format (//'Number of completed timesteps:',i6,/,'Time step ',i6, & ' partially done to provide convectively adjusted and ', & diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index e84fc837..51e7dd6b 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 @@ -12,7 +13,9 @@ module air_composition 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 +228,7 @@ subroutine air_composition_init() ! !************************************************************************ ! - ! add prognostic components of dry air + ! add prognostic components of air ! !************************************************************************ ! @@ -309,6 +312,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 +514,68 @@ 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 + ! (enthalpy for pressure-based dynamical cores and internal energy for z-based dynamical cores) + !--------------------------------------------------------------------------- + !=========================================================================== + subroutine water_composition_update(mmr, ncol, energy_formula, cp_or_cv_dycore, to_dry_factor) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + use string_utils, only: stringify + + 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), intent(out) :: cp_or_cv_dycore(:,:) ! enthalpy or heat capacity, dycore dependent [J K-1 kg-1] + 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 + 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 + ! Note: species index subset to 1: because SIMA currently uses index 0. See GitHub issue #334 in ESCOMP/CAM-SIMA. + 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 + ! Note: species index subset to 1: because SIMA currently uses index 0. See GitHub issue #334 in ESCOMP/CAM-SIMA. + 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; doi: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 (value = '//stringify((/energy_formula/))//') not supported') + end if + end subroutine water_composition_update !=========================================================================== !*************************************************************************** @@ -639,27 +685,34 @@ 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 + ! + ! if factor not present then tracer must be a dry mixing ratio + ! if factor present tracer*factor must be a 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(:,:) + ! factor: to convert tracer to dry mixing ratio + ! if provided, then tracer is not a dry mass 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 +728,37 @@ 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 + 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 +769,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 +777,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 +793,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 +912,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 +926,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 +945,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 +1012,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..8330ef64 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,89 @@ 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 - use string_utils, only: to_str - !----------------------------------------------------------------------- - ! Update the physics "constants" that vary - !------------------------------------------------------------------------- - - !------------------------------Arguments---------------------------------- + subroutine cam_thermo_dry_air_update(mmr, T, ncol, pver, update_thermo_variables, to_dry_factor) + use air_composition, only: dry_air_composition_update + use air_composition, only: update_zvirv + use string_utils, only: stringify - 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), 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) :: pver ! number of vertical levels + 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(:,:) ! conversion factor to dry if mmr is wet or moist - 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 '//stringify((/SIZE(to_dry_factor,1)/))//' but should be '//stringify((/ncol/))) + end if + if (SIZE(to_dry_factor, 2) /= pver) then + call endrun(subname//'DIM 2 of to_dry_factor is '//stringify((/SIZE(to_dry_factor,2)/))//' but should be '//stringify((/pver/))) + 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, pver, energy_formula, cp_or_cv_dycore, to_dry_factor) + use air_composition, only: water_composition_update + use string_utils, only: stringify + !----------------------------------------------------------------------- + ! Update the physics "constants" that vary + !------------------------------------------------------------------------- + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: pver ! number of vertical levels + integer, intent(in) :: energy_formula + real(kind_phys), intent(out) :: cp_or_cv_dycore(:,:) ! enthalpy or heat capacity, dycore dependent [J K-1 kg-1] + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + character(len=*), parameter :: subname = 'cam_thermo_water_update: ' + + if (present(to_dry_factor)) then + if (SIZE(to_dry_factor, 1) /= ncol) then + call endrun(subname//'DIM 1 of to_dry_factor is '//stringify((/SIZE(to_dry_factor,1)/))//' but should be '//stringify((/ncol/))) + end if + if (SIZE(to_dry_factor, 2) /= pver) then + call endrun(subname//'DIM 2 of to_dry_factor is '//stringify((/SIZE(to_dry_factor,2)/))//' but should be '//stringify((/pver/))) + end if + end if + + call water_composition_update(mmr, ncol, energy_formula, cp_or_cv_dycore, to_dry_factor=to_dry_factor) + end subroutine cam_thermo_water_update !=========================================================================== @@ -1554,28 +1595,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(:,:) - ! cp_or_cv: dry air heat capacity under constant pressure or - ! constant volume (depends on vcoord) + real(kind_phys), intent(in) :: pdel_in(:,:) + ! cp_or_cv: air heat capacity under constant pressure or + ! constant volume (depends on energy formula) 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 +1633,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 +1648,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 @@ -1633,12 +1684,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:) species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:) else - species_idx(:) = thermodynamic_active_species_idx(:) + species_idx(:) = thermodynamic_active_species_idx(1:) species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) end if else - species_idx(:) = thermodynamic_active_species_idx(:) + species_idx(:) = thermodynamic_active_species_idx(1:) species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) end if @@ -1649,78 +1700,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 +1803,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 +1831,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..df202c9d --- /dev/null +++ b/src/data/cam_thermo_formula.F90 @@ -0,0 +1,39 @@ +module cam_thermo_formula + + use runtime_obj, only: unset_int + + 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 = unset_int + ! energy_formula_physics: energy formula used for physics + integer, public :: energy_formula_physics = unset_int + + ! 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 0d7bbcf9..91658e20 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -14,6 +14,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 @@ -205,7 +206,7 @@ zi state_zi @@ -213,7 +214,7 @@ te_ini_phys state_te_ini_phys @@ -221,7 +222,7 @@ te_cur_phys state_te_cur_phys @@ -229,7 +230,7 @@ te_ini_dyn state_te_ini_dyn @@ -237,7 +238,7 @@ te_cur_dyn state_te_cur_dyn @@ -245,13 +246,34 @@ tw_ini state_tw_ini horizontal_dimension tw_cur state_tw_cur + + Total energy using dynamical core formula at the end of physics timestep + horizontal_dimension + 0.0 + + + flag indicating if dynamical core energy is not consistent with CAM physics and to perform adjustment of temperature and temperature tendency + .false. + .true. + + + + flag indicating if it is the first timestep of an initial run + reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure frontogenesis_function frontogenesis_angle - vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula - vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula - vertically_integrated_water_vapor_and_condensed_water_of_initial_state - vertically_integrated_water_vapor_and_condensed_water_of_current_state + vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + vertically_integrated_total_energy_using_dycore_energy_formula + vertically_integrated_total_water_at_start_of_physics_timestep + vertically_integrated_total_water + vertically_integrated_total_energy_at_end_of_physics_timestep tendency_of_air_temperature_due_to_model_physics @@ -416,6 +439,14 @@ horizontal_dimension vertical_layer_dimension zvir + + specific heat of air used in the dynamical core (enthalpy for pressure-based dynamical cores and internal energy for z-based dynamical cores) + horizontal_dimension vertical_layer_dimension + cp_or_cv_dycore + Run MPAS dynamical core to integrate the dynamical states with time. diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 6757158f..a2066e16 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -2,7 +2,8 @@ module dyn_coupling ! Modules from CAM-SIMA. use cam_abortutils, only: check_allocate, endrun use cam_constituents, only: const_is_water_species, const_qmin, num_advected - use cam_thermo, only: cam_thermo_update + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_MPAS use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, reverse, mpas_dynamical_core, & ncells_solve use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, & @@ -18,6 +19,7 @@ module dyn_coupling use physics_types, only: cappav, cpairv, rairv, zvirv, & dtime_phys, lagrangian_vertical, & phys_state, phys_tend + use physics_types, only: cp_or_cv_dycore use qneg, only: qneg_run use static_energy, only: update_dry_static_energy_run use string_utils, only: stringify @@ -326,15 +328,25 @@ subroutine set_physics_state_external() call endrun('Failed to find variable "constituent_properties"', subname, __LINE__) end if - ! Update `cappav`, `cpairv`, `rairv`, `zvirv`, etc. as needed by calling `cam_thermo_update`. + ! Update `cappav`, `cpairv`, `rairv`, `zvirv`, etc. as needed by calling `cam_thermo_dry_air_update`. ! Note that this subroutine expects constituents to be dry. - call cam_thermo_update( & - constituents, phys_state % t, ncells_solve, cam_runtime_opts % update_thermodynamic_variables()) + call cam_thermo_dry_air_update( & + constituents, phys_state % t, ncells_solve, pver, cam_runtime_opts % update_thermodynamic_variables()) + + ! update cp_or_cv_dycore in SIMA state. + ! (note: at this point q is dry) + call cam_thermo_water_update( & + mmr = constituents, & ! dry MMR + ncol = ncells_solve, & + pver = pver, & + energy_formula = ENERGY_FORMULA_DYCORE_MPAS, & + cp_or_cv_dycore = cp_or_cv_dycore & + ) ! This variable name is really misleading. It actually represents the reciprocal of Exner function ! with respect to surface pressure. This definition is sometimes used for boundary layer work. See ! the paragraph below equation 1.5.1c in doi:10.1007/978-94-009-3027-8. - ! Also note that `cappav` is updated externally by `cam_thermo_update`. + ! Also note that `cappav` is updated externally by `cam_thermo_dry_air_update`. do i = 1, ncells_solve phys_state % exner(i, :) = (phys_state % ps(i) / phys_state % pmid(i, :)) ** cappav(i, :) end do diff --git a/src/dynamics/none/dyn_comp.F90 b/src/dynamics/none/dyn_comp.F90 index 968e04e2..9ecb3022 100644 --- a/src/dynamics/none/dyn_comp.F90 +++ b/src/dynamics/none/dyn_comp.F90 @@ -60,6 +60,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) type(dyn_import_t), intent(out) :: dyn_in type(dyn_export_t), intent(out) :: dyn_out + ! Note: dynamical core energy formula is set in dyn_grid based on dynamical core + ! that provided the initial conditions file + end subroutine dyn_init !============================================================================== diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index bc714e22..ba1bf0ba 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -49,6 +49,7 @@ module dyn_grid ! Private module routines private :: find_units private :: find_dimension + private :: find_energy_formula !============================================================================== CONTAINS @@ -126,6 +127,7 @@ subroutine model_grid_init() ! We will handle errors for this routine call pio_seterrorhandling(fh_ini, PIO_BCAST_ERROR, oldmethod=err_handling) + ! Find the latitude variable and dimension(s) call cam_pio_find_var(fh_ini, (/ 'lat ', 'lat_d ', 'latitude' /), lat_name, & lat_vardesc, var_found) @@ -159,6 +161,11 @@ subroutine model_grid_init() write(iulog, *) subname, ': Grid is unstructured' end if end if + + ! Find the dynamical core from which snapshot was saved to populate energy formula used + ! Some information about the grid is needed to determine this. + call find_energy_formula(fh_ini, grid_is_latlon) + ! Find the longitude variable and dimension(s) call cam_pio_find_var(fh_ini, (/ 'lon ', 'lon_d ', 'longitude' /), lon_name, & lon_vardesc, var_found) @@ -626,4 +633,78 @@ subroutine find_dimension(file, dim_names, found_name, dim_len) end if end subroutine find_dimension + !=========================================================================== + + subroutine find_energy_formula(file, grid_is_latlon) + use pio, only: file_desc_t + use pio, only: pio_inq_att, pio_global, PIO_NOERR + use cam_thermo_formula, only: energy_formula_physics, energy_formula_dycore + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_MPAS + use physics_types, only: dycore_energy_consistency_adjust + use phys_vars_init_check, only: mark_as_initialized + + ! Find which dynamical core is used in and set the energy formulation + ! (also called vc_dycore in CAM) + ! + ! This functionality is only used to recognize the originating dynamical core + ! from the snapshot file in order to set the energy formulation when running + ! with the null dycore. Other dynamical cores set energy_formula_dycore at their + ! initialization. + + 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 + dycore_energy_consistency_adjust = .false. + 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 + dycore_energy_consistency_adjust = .true. + 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 + dycore_energy_consistency_adjust = .true. + 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 + call mark_as_initialized("flag_for_dycore_energy_consistency_adjustment") + + ! Mark other energy variables calculated by check_energy_timestep_init + ! here since it will always run when required + call mark_as_initialized("specific_heat_of_air_used_in_dycore") + call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_total_water_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_water") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep") + + end subroutine find_energy_formula + end module dyn_grid diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 8b56e9d9..572f663f 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_types, only: cp_or_cv_dycore 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 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,61 @@ 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. !----------------------------------------------------------------------------- + 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, & + pver = pver, & + update_thermo_variables = cam_runtime_opts%update_thermodynamic_variables() & + ) + else + zvirv(:,:) = zvir + end if - call cam_thermo_update(const_data_ptr, phys_state%t, pcols, & - cam_runtime_opts%update_thermodynamic_variables()) + ! + ! update cp_or_cv_dycore in SIMA state. + ! (note: at this point q is dry) + ! + call cam_thermo_water_update( & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + pver = pver, & + energy_formula = ENERGY_FORMULA_DYCORE_SE, & + cp_or_cv_dycore = cp_or_cv_dycore & + ) - !Call geopotential_temp CCPP scheme: + !$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_run(pver, lagrangian_vertical, pver, 1, & pverp, 1, num_advected, phys_state%lnpint, & phys_state%pint, phys_state%pmid, phys_state%pdel, & @@ -807,12 +852,6 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) phys_state%phis, phys_state%dse, cpairv, & errflg, errmsg) -!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) -#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..4f70ae2c 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -566,6 +566,9 @@ 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 + use physics_types, only: dycore_energy_consistency_adjust + use phys_vars_init_check, only: mark_as_initialized !SE dycore: use prim_advance_mod, only: prim_advance_init @@ -642,6 +645,14 @@ 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 + call mark_as_initialized("total_energy_formula_for_dycore") + + ! Dynamical core energy is not consistent with CAM physics and requires + ! temperature and temperature tendency adjustment at end of physics. + dycore_energy_consistency_adjust = .true. + call mark_as_initialized("flag_for_dycore_energy_consistency_adjustment") ! Now allocate and set condenstate vars allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers @@ -1852,6 +1863,17 @@ subroutine read_inidat(dyn_in) call mark_as_initialized("tendency_of_air_temperature_due_to_model_physics") call mark_as_initialized("tendency_of_eastward_wind_due_to_model_physics") call mark_as_initialized("tendency_of_northward_wind_due_to_model_physics") + call mark_as_initialized("specific_heat_of_air_used_in_dycore") + + ! These energy variables are calculated by check_energy_timestep_init + ! but need to be marked here + call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_total_water_at_start_of_physics_timestep") + call mark_as_initialized("vertically_integrated_total_water") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep") end subroutine read_inidat 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/utils/phys_comp.F90 b/src/physics/utils/phys_comp.F90 index 5dbfc20a..aa1feef6 100644 --- a/src/physics/utils/phys_comp.F90 +++ b/src/physics/utils/phys_comp.F90 @@ -174,10 +174,12 @@ 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 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..8dd3dca3 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 @@ -376,8 +378,6 @@ end subroutine phys_grid_init !======================================================================== real(r8) function get_dlat_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! latitude of a physics column in degrees ! Dummy argument @@ -396,8 +396,6 @@ end function get_dlat_p !======================================================================== real(r8) function get_dlon_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! longitude of a physics column in degrees ! Dummy argument @@ -416,8 +414,6 @@ end function get_dlon_p !======================================================================== real(r8) function get_rlat_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! latitude of a physics column in radians ! Dummy argument @@ -436,8 +432,6 @@ end function get_rlat_p !======================================================================== real(r8) function get_rlon_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! longitude of a physics column in radians ! Dummy argument @@ -456,8 +450,6 @@ end function get_rlon_p !======================================================================== real(r8) function get_area_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! area of a physics column in radians squared ! Dummy argument @@ -475,9 +467,25 @@ end function get_area_p !======================================================================== + real(r8) function get_wght_p(index) + ! 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 !----------------------------------------------------------------------- ! ! get_rlat_all_p: Return all latitudes (in radians) on task. @@ -506,8 +514,6 @@ end subroutine get_rlat_all_p !======================================================================== subroutine get_rlon_all_p(rlondim, rlons) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun !----------------------------------------------------------------------- ! ! get_rlon_all_p: Return all longitudes (in radians) on task. @@ -535,8 +541,35 @@ end subroutine get_rlon_all_p !======================================================================== + subroutine get_wght_all_p(wghtdim, wghts) + !----------------------------------------------------------------------- + ! + ! 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 ! Return the dynamics local block number and block offset(s) for ! the physics column indicated by . @@ -568,8 +601,6 @@ end subroutine get_dyn_col_p !======================================================================== integer function global_index_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! global column index of a physics column ! Dummy argument @@ -586,8 +617,6 @@ integer function global_index_p(index) end function global_index_p integer function local_index_p(index) - use cam_logfile, only: iulog - use cam_abortutils, only: endrun ! global column index of a physics column ! Dummy argument diff --git a/src/utils/gmean_mod.F90 b/src/utils/gmean_mod.F90 new file mode 100644 index 00000000..7959ebf0 --- /dev/null +++ b/src/utils/gmean_mod.F90 @@ -0,0 +1,263 @@ +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 + 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 + use cam_abortutils, only: check_allocate + ! + ! 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 + ! + real(r8), parameter :: pi4 = 4.0_r8 * pi + character(len=*), parameter :: subname = 'gmean_fixed_repro: ' + + integer :: icol, ifld ! column, field indices + integer :: errflg + + real(r8) :: wght(pcols) ! integration weights + real(r8), allocatable :: xfld(:,:) ! weighted summands + + errflg = 0 + + allocate(xfld(pcols, nflds), stat=errflg) + call check_allocate(errflg, subname, 'xfld(pcols, nflds)', & + file=__FILE__, line=__LINE__) + + ! 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(:) / pi4 + + end subroutine gmean_fixed_repro + +end module gmean_mod diff --git a/test/existing-test-failures.txt b/test/existing-test-failures.txt index 253b5633..05e0a145 100644 --- a/test/existing-test-failures.txt +++ b/test/existing-test-failures.txt @@ -2,12 +2,7 @@ SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_intel.cam-outfrq_kessler_mpas_derecho SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_gnu.cam-outfrq_kessler_mpas_derecho (Overall: FAIL) - will fail until MPAS is fully integrated -SMS_Ln9.ne5pg3_ne5pg3_mg37.FTJ16.derecho_intel.cam-outfrq_se_cslam (Overall: FAIL) -SMS_Ln9.ne5pg3_ne5pg3_mg37.FKESSLER.derecho_intel.cam-outfrq_se_cslam (Overall: FAIL) -SMS_Ln2.ne3pg3_ne3pg3_mg37.FPHYStest.derecho_intel.cam-outfrq_kessler_derecho (Overall: FAIL) -SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_intel.cam-outfrq_se_cslam_analy_ic (Overall: FAIL) -SMS_Ln9.ne5pg3_ne5pg3_mg37.FTJ16.derecho_gnu.cam-outfrq_se_cslam (Overall: FAIL) -SMS_Ln9.ne5pg3_ne5pg3_mg37.FKESSLER.derecho_gnu.cam-outfrq_se_cslam (Overall: FAIL) -SMS_Ln2.ne3pg3_ne3pg3_mg37.FPHYStest.derecho_gnu.cam-outfrq_kessler_derecho (Overall: FAIL) -SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_gnu.cam-outfrq_se_cslam_analy_ic (Overall: FAIL) - - will fail until https://github.com/ESCOMP/CAM-SIMA/pull/316 is merged +SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_intel.cam-outfrq_se_cslam_analy_ic (Overall: PEND) details: +SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_gnu.cam-outfrq_se_cslam_analy_ic (Overall: PEND) details: + - build failure due to dadadj_apply_qv_tendency removal, needs to be updated to use constituent tendency updater atmospheric_physics#179 + - also expected runtime failure CAM-SIMA#335 diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 0ac72a21..dbd84b99 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -160,42 +160,47 @@ state_zi - + te_ini_phys state_te_ini_phys - + te_cur_phys state_te_cur_phys - + te_ini_dyn state_te_ini_dyn - + te_cur_dyn state_te_cur_dyn - + tw_ini state_tw_ini - + tw_cur state_tw_cur + + + cp_or_cv_dycore + + RHO From 2ed783aacea3046afc0c0e97f3b6252aff41473b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 16 Dec 2024 12:08:04 -0500 Subject: [PATCH 2/2] Update standard names for tropopause_find (#329) Tag name (required for release branches): Originator(s): @jimmielin Description (include the issue title, and the keyword ['closes', 'fixes', 'resolves'] followed by the issue number): Fixes #308 by updating tropopause_find standard names. Describe any changes made to build system: N/A Describe any changes made to the namelist: N/A List any changes to the defaults for the input datasets (e.g. boundary datasets): N/A List all files eliminated and why: N/A List all files added and what they do: N/A List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) Updates standard names for tropopause_find M src/physics/utils/tropopause_climo_read.F90 M src/physics/utils/tropopause_climo_read.meta If there are new failures (compared to the `test/existing-test-failures.txt` file), have them OK'd by the gatekeeper, note them here, and add them to the file. If there are baseline differences, include the test and the reason for the diff. What is the nature of the change? Roundoff? derecho/intel/aux_sima: derecho/gnu/aux_sima: If this changes climate describe any run(s) done to evaluate the new climate in enough detail that it(they) could be reproduced: CAM-SIMA date used for the baseline comparison tests if different than latest: --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- src/physics/utils/tropopause_climo_read.F90 | 4 ++-- src/physics/utils/tropopause_climo_read.meta | 10 +++++----- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/.gitmodules b/.gitmodules index 5322028a..0fdc3af0 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/ESCOMP/atmospheric_physics - fxtag = e7a599f4bb1533f7cdcd8723b1f864e11578e96c + fxtag = 491e56247815ef23bfd8dba65d1e3c3b78ba164a fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index e7a599f4..491e5624 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit e7a599f4bb1533f7cdcd8723b1f864e11578e96c +Subproject commit 491e56247815ef23bfd8dba65d1e3c3b78ba164a diff --git a/src/physics/utils/tropopause_climo_read.F90 b/src/physics/utils/tropopause_climo_read.F90 index 7033d91b..c9369142 100644 --- a/src/physics/utils/tropopause_climo_read.F90 +++ b/src/physics/utils/tropopause_climo_read.F90 @@ -233,8 +233,8 @@ subroutine tropopause_climo_read_file() !-------------------------------------------------------- ! Mark variables as initialized so they are not read from initial conditions !-------------------------------------------------------- - call mark_as_initialized('tropopause_air_pressure_from_climatology_dataset') - call mark_as_initialized('tropopause_calendar_days_from_climatology') + call mark_as_initialized('tropopause_air_pressure_from_tropopause_climatology_dataset') + call mark_as_initialized('tropopause_calendar_days_from_tropopause_climatology') end subroutine tropopause_climo_read_file end module tropopause_climo_read diff --git a/src/physics/utils/tropopause_climo_read.meta b/src/physics/utils/tropopause_climo_read.meta index 6d4d8538..e42cc9b0 100644 --- a/src/physics/utils/tropopause_climo_read.meta +++ b/src/physics/utils/tropopause_climo_read.meta @@ -6,18 +6,18 @@ name = tropopause_climo_read type = module [ tropp_slices ] - standard_name = number_of_months_in_year + standard_name = number_of_time_slices_in_tropopause_climatology_dataset units = 1 type = integer dimensions = () [ tropp_p_loc ] - standard_name = tropopause_air_pressure_from_climatology_dataset + standard_name = tropopause_air_pressure_from_tropopause_climatology_dataset units = Pa type = real | kind = kind_phys - dimensions = (horizontal_dimension, number_of_months_in_year) + dimensions = (horizontal_dimension, number_of_time_slices_in_tropopause_climatology_dataset) [ tropp_days ] - standard_name = tropopause_calendar_days_from_climatology + standard_name = tropopause_calendar_days_from_tropopause_climatology long_name = Climatological tropopause calendar day indices from file units = 1 type = real | kind = kind_phys - dimensions = (number_of_months_in_year) + dimensions = (number_of_time_slices_in_tropopause_climatology_dataset)