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