From 86489702461192ccb13ae7cff9d13f0db004812b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Sep 2024 16:24:55 -0400 Subject: [PATCH 01/41] Initial implementation of cp_to_cv_dycore/cam_thermo_water_update from stash Needs work to split into cam_thermo_water_update. This is WIP code --- src/data/air_composition.F90 | 31 ++++++- src/data/cam_thermo.F90 | 108 +++++++++++++++--------- src/data/physconst.F90 | 4 + src/data/registry.xml | 7 ++ src/dynamics/se/dyn_comp.F90 | 3 + src/dynamics/tests/dyn_tests_utils.F90 | 6 ++ src/dynamics/tests/dyn_tests_utils.meta | 17 ++++ 7 files changed, 134 insertions(+), 42 deletions(-) create mode 100644 src/dynamics/tests/dyn_tests_utils.meta diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index f202ae97..d69685a9 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -6,6 +6,7 @@ module air_composition use runtime_obj, only: unset_real, unset_int use phys_vars_init_check, only: std_name_len use physics_types, only: cpairv, rairv, cappav, mbarv, zvirv + use physics_types, only: cp_or_cv_dycore implicit none private @@ -514,12 +515,16 @@ end subroutine air_composition_init !------------------------------------------------------------------------- !=========================================================================== - subroutine air_composition_update(mmr, ncol, to_moist_factor) + subroutine air_composition_update(mmr, ncol, vcoord, to_moist_factor) + use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) + character(len=*), parameter :: subname = 'air_composition_update' + call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, & rairv(:ncol, :), fact=to_moist_factor) call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx, & @@ -529,6 +534,30 @@ subroutine air_composition_update(mmr, ncol, to_moist_factor) cappav(:ncol,:) = rairv(:ncol,:) / cpairv(:ncol,:) + ! update enthalpy or internal energy scaling factor for energy consistency with CAM physics + ! FIXME hplin 9/5/24: why is air_composition_update in CAM-SIMA to_moist_factor + ! whereas water_composition_update in CAM to_dry_factor. + ! get_cp_dry implementation appears to be the same?, but + ! get_cp implementation appears to take an 1/factor. + ! Nothing in CAM-SIMA passes to_moist_factor to cam_thermo_update for now, so it is not possible + ! to know which is correct. In CAM, physpkg.F90 calls cam_thermo_water_update with to_dry_factor=pdel/pdeldry. + if (vcoord == vc_moist_pressure) then + ! FV: moist pressure vertical coordinate does not need update. + else if (vcoord == vc_dry_pressure) then + ! SE + call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & + dp_dry=to_moist_factor, active_species_idx_dycore=thermodynamic_active_species_idx) + else if (vcoord == vc_height) then + ! MPAS + ! Note hplin 9/5/24: the below code is "unused and untested" and reproduced verbatim (together with this warning) + ! from CAM air_composition.F90 + call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, cp_or_cv_dycore(:ncol,:), & + fact=to_moist_factor, R_dry=rairv(:ncol,:)) + cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) + else + call endrun(subname//': dycore vcoord not supported') + end if + end subroutine air_composition_update !=========================================================================== diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 59dd1c83..233f92c7 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -1554,20 +1554,23 @@ 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_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 physconst, only: rga, latvap, latice ! Dummy arguments ! tracer: tracer mixing ratio + ! + ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry real(kind_phys), intent(in) :: tracer(:,:,:) + logical, intent(in) :: moist_mixing_ratio ! pdel: pressure level thickness - real(kind_phys), intent(in) :: pdel(:,:) + real(kind_phys), intent(in) :: pdel_in(:,:) ! cp_or_cv: dry air heat capacity under constant pressure or ! constant volume (depends on vcoord) real(kind_phys), intent(in) :: cp_or_cv(:,:) @@ -1575,7 +1578,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & 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(:) + 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 +1591,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 +1606,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & ! Local variables real(kind_phys) :: ke_vint(SIZE(tracer, 1)) ! Vertical integral of KE - real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of SE + real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of enthalpy or internal energy + real(kind_phys) :: po_vint(SIZE(tracer, 1)) ! Vertical integral of PHIS or potential energy real(kind_phys) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv real(kind_phys) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq real(kind_phys) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice + real(kind_phys) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) ! moist pressure level thickness real(kind_phys) :: latsub ! latent heat of sublimation integer :: ierr @@ -1649,51 +1658,56 @@ 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 ', & + if (.not. present(ptop).or. (.not. present(phis))) then + write(iulog, *) subname, ' ptop and phis must be present for ', & 'moist/dry pressure vertical coordinate' - call endrun(subname//': ps and phis must be present for '// & + call endrun(subname//': ptop and phis must be present for '// & 'moist/dry pressure vertical coordinate') 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') + if (.not. present(phis)) then + write(iulog, *) subname, ' phis must be present for ', & + 'heigt-based vertical coordinate' + call endrun(subname//': phis must be present for '// & + 'height-based vertical coordinate') 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 @@ -1701,26 +1715,39 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & call endrun(subname//': vertical coordinate 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 +1761,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 +1789,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/physconst.F90 b/src/data/physconst.F90 index 8a528342..774b8f3d 100644 --- a/src/data/physconst.F90 +++ b/src/data/physconst.F90 @@ -114,6 +114,7 @@ subroutine physconst_readnl(nlfile) use mpi, only: mpi_real8 use cam_logfile, only: iulog use runtime_obj, only: unset_real + use dyn_tests_utils, only: vc_physics, vc_moist_pressure ! Dummy argument: filepath for file containing namelist input character(len=*), intent(in) :: nlfile @@ -287,6 +288,9 @@ subroutine physconst_readnl(nlfile) ez = omega / sqrt(0.375_kind_phys) + ! set physics vertical coordinate info + vc_physics = vc_moist_pressure + end subroutine physconst_readnl end module physconst diff --git a/src/data/registry.xml b/src/data/registry.xml index 41b31416..5ae88e03 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -415,6 +415,13 @@ horizontal_dimension vertical_layer_dimension zvir + + Enthalpy or internal energy scaling factor for energy consistency + horizontal_dimension vertical_layer_dimension + \section arg_table_dyn_tests_utils Argument Table +!! \htmlinclude dyn_tests_utils.html + integer :: vc_dycore ! vertical coordinate of dynamical core - set in dyn_comp.F90 + integer :: vc_physics ! vertical coordinate of physics - set in physconst.F90 + public :: vc_dycore, vc_physics + end module dyn_tests_utils diff --git a/src/dynamics/tests/dyn_tests_utils.meta b/src/dynamics/tests/dyn_tests_utils.meta new file mode 100644 index 00000000..91f30f6d --- /dev/null +++ b/src/dynamics/tests/dyn_tests_utils.meta @@ -0,0 +1,17 @@ +[ccpp-table-properties] + name = dyn_tests_utils + type = module + +[ccpp-arg-table] + name = dyn_tests_utils + type = module +[ vc_dycore ] + standard_name = vertical_coordinate_for_dynamical_core? + units = none + type = integer + dimensions = () +[ vc_physics ] + standard_name = vertical_coordinate_for_physics? + units = none + type = integer + dimensions = () From 333ad4eed8520f92c61a6731bfd589d2f873f097 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 16 Sep 2024 11:04:29 -0400 Subject: [PATCH 02/41] Initial port of updates to energy budget (port of CAM #761) into CAM-SIMA --- src/data/air_composition.F90 | 140 +++++++++++++++++++++----------- src/data/cam_thermo.F90 | 80 +++++++++++------- src/dynamics/se/dp_coupling.F90 | 111 +++++++++++++++++-------- 3 files changed, 220 insertions(+), 111 deletions(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index d69685a9..0da0a5b2 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 @@ -13,7 +14,8 @@ 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 @@ -511,54 +513,66 @@ 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, vcoord, to_moist_factor) - use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + 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 - integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - - character(len=*), parameter :: subname = 'air_composition_update' + 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 + + !=========================================================================== + !--------------------------------------------------------------------------- + ! water_composition_update: Update generalized cp or cv depending on dycore + !--------------------------------------------------------------------------- + !=========================================================================== + subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) + use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore + 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 - ! FIXME hplin 9/5/24: why is air_composition_update in CAM-SIMA to_moist_factor - ! whereas water_composition_update in CAM to_dry_factor. - ! get_cp_dry implementation appears to be the same?, but - ! get_cp implementation appears to take an 1/factor. - ! Nothing in CAM-SIMA passes to_moist_factor to cam_thermo_update for now, so it is not possible - ! to know which is correct. In CAM, physpkg.F90 calls cam_thermo_water_update with to_dry_factor=pdel/pdeldry. + ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module + if (vcoord == vc_moist_pressure) then ! FV: moist pressure vertical coordinate does not need update. else if (vcoord == vc_dry_pressure) then ! SE call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & - dp_dry=to_moist_factor, active_species_idx_dycore=thermodynamic_active_species_idx) + factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx, & + cpdry=cpairv(:ncol,:)) else if (vcoord == vc_height) then ! MPAS - ! Note hplin 9/5/24: the below code is "unused and untested" and reproduced verbatim (together with this warning) - ! from CAM air_composition.F90 - call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, cp_or_cv_dycore(:ncol,:), & - fact=to_moist_factor, R_dry=rairv(:ncol,:)) + call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, & + cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) + + ! internal energy coefficient for MPAS + ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) else call endrun(subname//': dycore vcoord not supported') end if - end subroutine air_composition_update + end subroutine water_composition_update !=========================================================================== !*************************************************************************** @@ -668,27 +682,33 @@ end subroutine get_cp_dry_2hd ! !*************************************************************************** ! - subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) use cam_abortutils, only: endrun use string_utils, only: to_str ! Dummy arguments ! tracer: Tracer array + ! + ! factor not present then tracer must be dry mixing ratio + ! if factor present tracer*factor must be dry mixing ratio + ! real(kind_phys), intent(in) :: tracer(:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) + ! dp: if provided then tracer is mass not mixing ratio + real(kind_phys), optional, intent(in) :: factor(:,:) ! active_species_idx_dycore: array of indices for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:) ! Local variables integer :: qdx, itrac real(kind_phys) :: sum_species(SIZE(cp, 1), SIZE(cp, 2)) real(kind_phys) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2)) - real(kind_phys) :: factor(SIZE(cp, 1), SIZE(cp, 2)) + real(kind_phys) :: factor_local(SIZE(cp, 1), SIZE(cp, 2)) integer :: idx_local(thermodynamic_active_species_num) character(len=*), parameter :: subname = 'get_cp_1hd: ' @@ -704,28 +724,38 @@ 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 + ! FIXME hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + 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 @@ -736,7 +766,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 @@ -744,14 +774,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 @@ -759,11 +790,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 @@ -872,9 +909,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 @@ -885,6 +923,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 @@ -903,12 +942,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 @@ -963,7 +1009,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 233f92c7..2b0af969 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 @@ -208,51 +210,67 @@ 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 - !------------------------------------------------------------------------- + subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) + use air_composition, only: dry_air_composition_update + use air_composition, only: update_zvirv + use string_utils, only: int2str - !------------------------------Arguments---------------------------------- + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) + real(kind_phys), intent(in) :: T(:,:) ! temperature + integer, intent(in) :: ncol ! number of columns + logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables + ! false: do not calculate composition-dependent thermo variables + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr moist convert - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array - real(kind_phys), intent(in) :: T(:,:) ! temperature - integer, intent(in) :: ncol ! number of columns - logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables - ! false: do not calculate composition-dependent thermo variables - - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - ! - !---------------------------Local storage------------------------------- - real(kind_phys):: sponge_factor(SIZE(mmr, 2)) - character(len=*), parameter :: subname = 'cam_thermo_update: ' + ! Local vars + real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) + character(len=*), parameter :: subname = 'cam_thermo_dry_air_update: ' - if (.not. update_thermo_variables) then - return + if (present(to_dry_factor)) then + if (SIZE(to_dry_factor, 1) /= ncol) then + call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) + end if end if - 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 - 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, vcoord, to_dry_factor) + use air_composition, only: water_composition_update + !----------------------------------------------------------------------- + ! Update the physics "constants" that vary + !------------------------------------------------------------------------- + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + call water_composition_update(mmr, ncol, vcoord, to_dry_factor=to_dry_factor) + end subroutine cam_thermo_water_update !=========================================================================== diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 8b56e9d9..64efb9b2 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -583,7 +583,10 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 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 use physics_grid, only: columns_on_task use geopotential_temp, only: geopotential_temp_run @@ -607,7 +610,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 +625,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 +687,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 + m = thermodynamic_active_species_idx(m_cnst) + ! for now m is just == ix_q (one entry) + 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 +734,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 +763,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,9 +783,60 @@ 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, & + 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 module air_composition. + ! (note: at this point q is dry) + ! + call cam_thermo_water_update( & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + vcoord = vc_dry_pressure & + ) + + !$omp parallel do num_threads(horz_num_threads) private (k, i) + do k = 1, nlev + do i = 1, pcols + ! CAM version: + ! phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) + + ! CAM-SIMA version (uses constant cappa): + phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa + 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, & From 041cdfb250dc0088f87d6c1ecb391c015e6b4b57 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 16 Sep 2024 13:17:21 -0400 Subject: [PATCH 03/41] Add is_first_timestep registry --- .gitmodules | 4 ++-- src/control/cam_comp.F90 | 7 +++++++ src/data/registry.xml | 6 ++++++ src/dynamics/se/dp_coupling.F90 | 5 +++-- src/physics/ncar_ccpp | 2 +- 5 files changed, 19 insertions(+), 5 deletions(-) diff --git a/.gitmodules b/.gitmodules index 94610494..f35d1e2e 100644 --- a/.gitmodules +++ b/.gitmodules @@ -19,8 +19,8 @@ fxDONOTUSEurl = https://github.com/MPAS-Dev/MPAS-Model.git [submodule "ncar-physics"] path = src/physics/ncar_ccpp - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = e95c172d7a5a0ebf054f420b08416228e211baa3 + url = https://github.com/jimmielin/atmospheric_physics + fxtag = 996a6d5 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 538bd641..7667b584 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -29,6 +29,7 @@ module cam_comp use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys use physics_types, only: calday + use physics_types, only: is_first_timestep use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -151,6 +152,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & dtime_phys = 0.0_r8 call mark_as_initialized('timestep_for_physics') + is_first_timestep = is_first_step() + call mark_as_initialized('is_first_timestep') + call init_pio_subsystem() ! Initializations using data passed from coupler. @@ -263,6 +267,9 @@ subroutine cam_timestep_init() use phys_comp, only: phys_timestep_init use stepon, only: stepon_timestep_init + ! Update timestep flags in physics state + is_first_timestep = is_first_step() + !---------------------------------------------------------- ! First phase of dynamics (at least couple from dynamics to physics) ! Return time-step for physics from dynamics. diff --git a/src/data/registry.xml b/src/data/registry.xml index 5ae88e03..e55e8a4c 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -251,6 +251,12 @@ horizontal_dimension tw_cur state_tw_cur + + + flag indicating if is first timestep of initial run + Date: Mon, 16 Sep 2024 13:41:33 -0400 Subject: [PATCH 04/41] Fix build issues; update factor in get_cp call to 1/dp --- src/control/cam_comp.F90 | 2 +- src/data/air_composition.F90 | 14 ++++++++--- src/data/cam_thermo.F90 | 3 ++- src/data/registry.xml | 7 +++--- src/dynamics/tests/dyn_tests_utils.meta | 8 +++--- src/dynamics/utils/dyn_thermo.F90 | 33 ++++++++++++++----------- 6 files changed, 41 insertions(+), 26 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 7667b584..27fd57a7 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -152,7 +152,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & dtime_phys = 0.0_r8 call mark_as_initialized('timestep_for_physics') - is_first_timestep = is_first_step() + is_first_timestep = .true. call mark_as_initialized('is_first_timestep') call init_pio_subsystem() diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 0da0a5b2..43ce33d0 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -16,6 +16,7 @@ module air_composition public :: air_composition_init 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 @@ -533,7 +534,7 @@ subroutine dry_air_composition_update(mmr, ncol, to_dry_factor) cappav(:ncol,:) = rairv(:ncol,:) / cpairv(:ncol,:) - end subroutine air_composition_update + end subroutine dry_air_composition_update !=========================================================================== !--------------------------------------------------------------------------- @@ -557,12 +558,19 @@ subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) ! FV: moist pressure vertical coordinate does not need update. else if (vcoord == vc_dry_pressure) then ! SE + + ! TODO hplin 9/17/24: for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) + ! (whereas CAM only allocates 1-indexed) I subset it here. But from the meaning of the code arguments + ! it is unknown where it was meant to be thermodynamic_active_species_idx_dycore. + ! This should be verified during code review. call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & - factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx, & + factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx(1:), & cpdry=cpairv(:ncol,:)) else if (vcoord == vc_height) then ! MPAS - call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, & + + ! TODO hplin 9/17/24 same here. + 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 diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 2b0af969..2c7f62a5 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -220,7 +220,7 @@ end subroutine cam_thermo_init subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) use air_composition, only: dry_air_composition_update use air_composition, only: update_zvirv - use string_utils, only: int2str + use string_utils, only: to_str 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 @@ -1579,6 +1579,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & 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 air_composition, only: dry_air_species_num use physconst, only: rga, latvap, latice ! Dummy arguments diff --git a/src/data/registry.xml b/src/data/registry.xml index e55e8a4c..de0c7439 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -16,6 +16,7 @@ $SRCROOT/src/data/ref_pres.meta $SRCROOT/src/dynamics/utils/vert_coord.meta $SRCROOT/src/dynamics/utils/hycoef.meta + $SRCROOT/src/dynamics/tests/dyn_tests_utils.meta @@ -254,7 +255,7 @@ + units="flag" type="logical"> flag indicating if is first timestep of initial run @@ -422,8 +423,8 @@ zvir Enthalpy or internal energy scaling factor for energy consistency horizontal_dimension vertical_layer_dimension diff --git a/src/dynamics/tests/dyn_tests_utils.meta b/src/dynamics/tests/dyn_tests_utils.meta index 91f30f6d..d236ec6f 100644 --- a/src/dynamics/tests/dyn_tests_utils.meta +++ b/src/dynamics/tests/dyn_tests_utils.meta @@ -6,12 +6,12 @@ name = dyn_tests_utils type = module [ vc_dycore ] - standard_name = vertical_coordinate_for_dynamical_core? - units = none + standard_name = vertical_coordinate_for_dynamical_core_tbd + units = 1 type = integer dimensions = () [ vc_physics ] - standard_name = vertical_coordinate_for_physics? - units = none + standard_name = vertical_coordinate_for_physics_tbd + units = 1 type = integer dimensions = () 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, & From 4ce842786b8b64374057058dd22a21dc3d2576eb Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 20 Sep 2024 14:39:16 -0400 Subject: [PATCH 05/41] Add gmean_mod to src/utils including support infra in physics_grid.F90 --- src/physics/utils/physics_grid.F90 | 52 ++++++ src/utils/gmean_mod.F90 | 256 +++++++++++++++++++++++++++++ 2 files changed, 308 insertions(+) create mode 100644 src/utils/gmean_mod.F90 diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90 index 39fc0f99..cf28f98d 100644 --- a/src/physics/utils/physics_grid.F90 +++ b/src/physics/utils/physics_grid.F90 @@ -39,8 +39,10 @@ module physics_grid public :: get_rlat_p ! latitude of a physics column in radians public :: get_rlon_p ! longitude of a physics column in radians public :: get_area_p ! area of a physics column in radians squared + public :: get_wght_p ! weight of a physics column in radians squared public :: get_rlat_all_p ! latitudes of physics cols on task (radians) public :: get_rlon_all_p ! longitudes of physics cols on task (radians) + public :: get_wght_all_p ! weights of physics cols on task public :: get_dyn_col_p ! dynamics local blk number and blk offset(s) public :: global_index_p ! global column index of a physics column public :: local_index_p ! local column index of a physics column @@ -475,6 +477,26 @@ end function get_area_p !======================================================================== + real(r8) function get_wght_p(index) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + ! weight of a physics column in radians squared + + ! Dummy argument + integer, intent(in) :: index + ! Local variables + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_p' + + ! Check that input is valid: + call check_phys_input(subname, index) + + get_wght_p = phys_columns(index)%weight + + end function get_wght_p + + !======================================================================== + subroutine get_rlat_all_p(rlatdim, rlats) use cam_logfile, only: iulog use cam_abortutils, only: endrun @@ -535,6 +557,36 @@ end subroutine get_rlon_all_p !======================================================================== + subroutine get_wght_all_p(wghtdim, wghts) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + !----------------------------------------------------------------------- + ! + ! get_wght_all_p: Return all weights on task. + ! + !----------------------------------------------------------------------- + ! Dummy Arguments + integer, intent(in) :: wghtdim ! declared size of output array + real(r8), intent(out) :: wghts(wghtdim) ! array of weights + + ! Local variables + integer :: index ! loop index + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_all_p: ' + + !----------------------------------------------------------------------- + + ! Check that input is valid: + call check_phys_input(subname, wghtdim) + + do index = 1, wghtdim + wghts(index) = phys_columns(index)%weight + end do + + end subroutine get_wght_all_p + + !======================================================================== + subroutine get_dyn_col_p(index, blk_num, blk_ind) use cam_logfile, only: iulog use cam_abortutils, only: endrun diff --git a/src/utils/gmean_mod.F90 b/src/utils/gmean_mod.F90 new file mode 100644 index 00000000..ca6dd6d4 --- /dev/null +++ b/src/utils/gmean_mod.F90 @@ -0,0 +1,256 @@ +module gmean_mod + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Perform global mean calculations for energy conservation and other checks. + ! + ! Method: + ! Reproducible (scalable): + ! Convert to fixed point (integer representation) to enable + ! reproducibility when using MPI collectives. + ! If error checking is on (via setting reprosum_diffmax > 0 and + ! reprosum_recompute = .true. in user_nl_cpl), shr_reprosum_calc will + ! check the accuracy of its computation with a fast but + ! non-reproducible algorithm. If any error is reported, report + ! the difference and the expected sum and abort run (call endrun) + ! + ! gmean_mod in to_be_ccppized is different from the CAM version and + ! has chunk support removed. + ! + ! + !----------------------------------------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use physics_grid, only: pcols => columns_on_task + use perf_mod, only: t_startf, t_stopf + use cam_logfile, only: iulog + + implicit none + private + + public :: gmean ! compute global mean of 2D fields on physics decomposition + + interface gmean + module procedure gmean_arr + module procedure gmean_scl + end interface gmean + + private :: gmean_fixed_repro + private :: gmean_float_norepro + + ! Set do_gmean_tests to .true. to run a gmean challenge test + logical, private :: do_gmean_tests = .false. + +CONTAINS + + ! + !======================================================================== + ! + + subroutine gmean_arr (arr, arr_gmean, nflds) + use shr_strconvert_mod, only: toString + use spmd_utils, only: masterproc + use cam_abortutils, only: endrun + use shr_reprosum_mod, only: shr_reprosum_reldiffmax, shr_reprosum_recompute, shr_reprosum_tolExceeded + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! + ! Method is to call shr_reprosum_calc (called from gmean_fixed_repro) + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols, nflds) + real(r8), intent(out) :: arr_gmean(nflds) ! global means + ! + ! Local workspace + ! + real(r8) :: rel_diff(2, nflds) + integer :: ifld ! field index + integer :: num_err + logical :: write_warning + ! + !----------------------------------------------------------------------- + ! + call t_startf('gmean_arr') + call t_startf ('gmean_fixed_repro') + call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds) + call t_stopf ('gmean_fixed_repro') + + ! check that "fast" reproducible sum is accurate enough. If not, calculate + ! using old method + write_warning = masterproc + num_err = 0 + if (shr_reprosum_tolExceeded('gmean', nflds, write_warning, & + iulog, rel_diff)) then + if (shr_reprosum_recompute) then + do ifld = 1, nflds + if (rel_diff(1, ifld) > shr_reprosum_reldiffmax) then + call gmean_float_norepro(arr(:,ifld), arr_gmean(ifld), ifld) + num_err = num_err + 1 + end if + end do + end if + end if + call t_stopf('gmean_arr') + if (num_err > 0) then + call endrun('gmean: '//toString(num_err)//' reprosum errors found') + end if + + end subroutine gmean_arr + + ! + !======================================================================== + ! + + subroutine gmean_scl (arr, gmean) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of a single field in "arr" in the physics grid + ! + !----------------------------------------------------------------------- + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + ! Input array, chunked + real(r8), intent(out):: gmean ! global means + ! + ! Local workspace + ! + integer, parameter :: nflds = 1 + real(r8) :: gmean_array(nflds) + real(r8) :: array(pcols, nflds) + integer :: ncols, lchnk + + array(:ncols, 1) = arr(:ncols) + call gmean_arr(array, gmean_array, nflds) + gmean = gmean_array(1) + + end subroutine gmean_scl + + ! + !======================================================================== + ! + + subroutine gmean_float_norepro(arr, repro_sum, index) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of in the physics chunked + ! decomposition using a fast but non-reproducible algorithm. + ! Log that value along with the value computed by + ! shr_reprosum_calc () + ! + !----------------------------------------------------------------------- + + use physconst, only: pi + use spmd_utils, only: masterproc, masterprocid, mpicom + use mpi, only: mpi_real8, mpi_sum + use physics_grid, only: get_wght_p + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + real(r8), intent(in) :: repro_sum ! Value computed by reprosum + integer, intent(in) :: index ! Index of field in original call + ! + ! Local workspace + ! + integer :: icol + integer :: ierr + real(r8) :: wght + real(r8) :: check + real(r8) :: check_sum + real(r8), parameter :: pi4 = 4.0_r8 * pi + + ! + !----------------------------------------------------------------------- + ! + ! Calculate and print out non-reproducible value + check = 0.0_r8 + do icol = 1, pcols + wght = get_wght_p(icol) + check = check + arr(icol) * wght + end do + call MPI_reduce(check, check_sum, 1, mpi_real8, mpi_sum, & + masterprocid, mpicom, ierr) + + ! normalization + check_sum = check_sum / pi4 + + if (masterproc) then + write(iulog, '(a,i0,2(a,e20.13e2))') 'gmean(', index, ') = ', & + check_sum, ', reprosum reported ', repro_sum + end if + + end subroutine gmean_float_norepro + + ! + !======================================================================== + ! + subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! with a reproducible yet scalable implementation + ! based on a fixed-point algorithm. + ! + !----------------------------------------------------------------------- + use spmd_utils, only: mpicom + use physics_grid, only: get_wght_all_p + use physics_grid, only: ngcols_p => num_global_phys_cols + use physconst, only: pi + use shr_reprosum_mod, only: shr_reprosum_calc + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols,nflds) + ! arr_gmean: output global sums + real(r8), intent(out) :: arr_gmean(nflds) + ! rel_diff: relative and absolute differences from shr_reprosum_calc + real(r8), intent(out) :: rel_diff(2, nflds) + ! + ! Local workspace + ! + integer :: icol, ifld ! column, field indices + + real(r8) :: wght(pcols) ! integration weights + real(r8), allocatable :: xfld(:,:) ! weighted summands + ! + !----------------------------------------------------------------------- + ! + allocate(xfld(pcols, nflds)) + + ! pre-weight summands + call get_wght_all_p(pcols, wght) + + do ifld = 1, nflds + do icol = 1, pcols + xfld(icol, ifld) = arr(icol, ifld) * wght(icol) + end do + end do + + ! call fixed-point algorithm + call shr_reprosum_calc ( & + arr = xfld, & + arr_gsum = arr_gmean, & + nsummands = pcols, & ! # of local summands + dsummands = pcols, & ! declared first dimension of arr. + nflds = nflds, & + commid = mpicom, & + rel_diff = rel_diff & + ) + + deallocate(xfld) + ! final normalization + arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi) + + end subroutine gmean_fixed_repro + +end module gmean_mod From 2372f7fc2fc1c71900f57ce1131c2b323f8803c3 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 2 Oct 2024 17:52:57 -0400 Subject: [PATCH 06/41] Fix for b4b to CAM-SIMA (w/ History); add notes on modifications and/or divergence from CAM and/or -SIMA --- src/data/air_composition.F90 | 15 ++++++++------- src/data/cam_thermo.F90 | 4 ++++ src/dynamics/se/dp_coupling.F90 | 16 ++++++++++++++-- src/dynamics/se/dycore/prim_advance_mod.F90 | 6 +++++- 4 files changed, 31 insertions(+), 10 deletions(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 43ce33d0..443565eb 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -559,17 +559,15 @@ subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) else if (vcoord == vc_dry_pressure) then ! SE - ! TODO hplin 9/17/24: for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) - ! (whereas CAM only allocates 1-indexed) I subset it here. But from the meaning of the code arguments - ! it is unknown where it was meant to be thermodynamic_active_species_idx_dycore. - ! This should be verified during code review. + ! **TEMP** TODO hplin 9/17/24: + ! for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) + ! (whereas CAM only allocates 1-indexed) I subset it here. This should be verified during code + ! review. call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx(1:), & cpdry=cpairv(:ncol,:)) else if (vcoord == vc_height) then ! MPAS - - ! TODO hplin 9/17/24 same here. call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), & cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) @@ -745,7 +743,10 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpd end do if (dry_air_species_num == 0) then - ! FIXME hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + ! **TMP** TODO CHECK hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + ! and dry_air_species_num is == 0 in CAM-SIMA as of 10/1/24 which means this clause will + ! change answers (?) -- previously in CAM-SIMA there was only a call to get_cp_dry + ! and not the two IF-clauses here. sum_cp = thermodynamic_active_species_cp(0) else if (present(cpdry)) then ! diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 2c7f62a5..a6de3d2b 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -233,6 +233,10 @@ subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_d 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_dry_factor)) then if (SIZE(to_dry_factor, 1) /= ncol) then call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index ea251163..1eca0982 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -690,9 +690,16 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! 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) + ! **TEMP** TODO CHECK hplin: check indices to use here + ! in cam6_3_109 after fix in 6_3_127: loop from dry_air_species_num + 1, thermodynamic_active_species_num + ! in cam-sima: factor_array only contains m = ix_q + ! hplin -- to get same answers as CAM-SIMA ix_q would need to be used, + ! but the CAM version appears to be correct. this should be science checked + ! // **TEMP** do m_cnst = dry_air_species_num + 1, thermodynamic_active_species_num m = thermodynamic_active_species_idx(m_cnst) - ! for now m is just == ix_q (one entry) + ! write(6,*) "hplin dp_coupling m, m_cnst", m, m_cnst + ! m = ix_q do k = 1, nlev do i = 1, pcols ! at this point all q's are dry @@ -784,6 +791,10 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! Compute molecular viscosity(kmvis) and conductivity(kmcnd). ! Update zvirv registry variable; calculated for WACCM-X. !----------------------------------------------------------------------------- + ! **TEMP** TODO CHECK hplin: CAM has this if-clause for dry_air_species_num > 0 + ! or otherwise uses zvirv = zvir. CAM-SIMA previously did not have this, and + ! instead has a switch for update_thermodynamic_variables. Check if we still want + ! this if-clause or change it to something else. if (dry_air_species_num > 0) then call cam_thermo_dry_air_update( & mmr = const_data_ptr, & ! dry MMR @@ -808,6 +819,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !$omp parallel do num_threads(horz_num_threads) private (k, i) do k = 1, nlev do i = 1, pcols + ! **TEMP** TODO CHECK hplin: CAM and CAM-SIMA version differ ! CAM version: ! phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) @@ -839,7 +851,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end do end do - !Call geopotential_temp CCPP scheme: + ! Call geopotential_temp CCPP scheme: call geopotential_temp_run(pver, lagrangian_vertical, pver, 1, & pverp, 1, num_advected, phys_state%lnpint, & phys_state%pint, phys_state%pmid, phys_state%pdel, & 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 From c015e28ab9ea5e590230c05df7c5411d7adf8b6a Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 15 Oct 2024 15:04:39 -0400 Subject: [PATCH 07/41] Update ncar-physics external --- src/data/cam_thermo.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index a6de3d2b..a776cece 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -1717,7 +1717,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & case(vc_height) if (.not. present(phis)) then write(iulog, *) subname, ' phis must be present for ', & - 'heigt-based vertical coordinate' + 'height-based vertical coordinate' call endrun(subname//': phis must be present for '// & 'height-based vertical coordinate') end if From e13e7ea0de0592bd2517cfcac32ea0bf3fabc6a6 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 15 Oct 2024 17:44:53 -0400 Subject: [PATCH 08/41] Fix build (part 1); update standard names and atmos_phys external --- src/data/registry.xml | 2 +- src/dynamics/tests/dyn_tests_utils.meta | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index de0c7439..daf04535 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -423,7 +423,7 @@ zvir Enthalpy or internal energy scaling factor for energy consistency diff --git a/src/dynamics/tests/dyn_tests_utils.meta b/src/dynamics/tests/dyn_tests_utils.meta index d236ec6f..7bfcbf1c 100644 --- a/src/dynamics/tests/dyn_tests_utils.meta +++ b/src/dynamics/tests/dyn_tests_utils.meta @@ -6,12 +6,12 @@ name = dyn_tests_utils type = module [ vc_dycore ] - standard_name = vertical_coordinate_for_dynamical_core_tbd + standard_name = total_energy_formula_for_dycore units = 1 type = integer dimensions = () [ vc_physics ] - standard_name = vertical_coordinate_for_physics_tbd + standard_name = total_energy_formula_for_physics units = 1 type = integer dimensions = () From 0d16be936ddfbb1e0c3a223671fadf15e8c955cd Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 15 Oct 2024 17:45:09 -0400 Subject: [PATCH 09/41] New const_get_index logic without cam_ccpp_cap dependency --- cime_config/cam_autogen.py | 10 ++++++++++ src/physics/utils/cam_constituents.F90 | 8 ++++++++ 2 files changed, 18 insertions(+) diff --git a/cime_config/cam_autogen.py b/cime_config/cam_autogen.py index 3d5ed312..f5467a6a 100644 --- a/cime_config/cam_autogen.py +++ b/cime_config/cam_autogen.py @@ -614,6 +614,7 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, # Copy to_be_ccppized utility modules to the build directory, # as SIMA cam_constituents depends on them. +<<<<<<< HEAD # Note: to_be_ccppized utility modules to be removed once functionality is migrated # to SIMA or CCPPized in atmospheric_physics. atm_phys_to_be_ccppized_dir = os.path.join(atm_phys_top_dir, "to_be_ccppized") @@ -630,6 +631,15 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, for to_be_ccppized_file in atm_phys_to_be_ccppized_files: shutil.copy(to_be_ccppized_file, physics_blddir) # end for +======= + atm_phys_to_be_ccppized_dir = os.path.join(atm_phys_top_dir, "to_be_ccppized") + + # Check that the directory exists and copy to build directory + if os.path.isdir(atm_phys_to_be_ccppized_dir): + atm_phys_to_be_ccppized_files = glob.glob(os.path.join(atm_phys_to_be_ccppized_dir, "*.F90")) + for to_be_ccppized_file in atm_phys_to_be_ccppized_files: + shutil.copy(to_be_ccppized_file, physics_blddir) +>>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) # end if if do_gen_ccpp or do_gen_nl: diff --git a/src/physics/utils/cam_constituents.F90 b/src/physics/utils/cam_constituents.F90 index 19329153..3883aba5 100644 --- a/src/physics/utils/cam_constituents.F90 +++ b/src/physics/utils/cam_constituents.F90 @@ -292,7 +292,11 @@ subroutine const_get_index(name, cindex, abort, warning, caller) use cam_abortutils, only: endrun use cam_logfile, only: iulog use phys_vars_init_check, only: std_name_len +<<<<<<< HEAD use string_utils, only: stringify +======= + use string_utils, only: to_str +>>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) ! Get the index of a constituent with standard name, . ! Setting optional argument to .false. returns control to @@ -319,7 +323,11 @@ subroutine const_get_index(name, cindex, abort, warning, caller) call ccpp_const_get_idx(const_props, name, cindex, errmsg, errcode) if (errcode /= 0) then +<<<<<<< HEAD call endrun(subname//"Error "//stringify((/errcode/))//": "// & +======= + call endrun(subname//"Error "//to_str(errcode)//": "// & +>>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) trim(errmsg), file=__FILE__, line=__LINE__) endif From 8fca3a42b6db828e185ec4df51665e1dcae8eae4 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 22 Oct 2024 17:56:43 -0400 Subject: [PATCH 10/41] Update vcoord to energy_formula --- src/data/air_composition.F90 | 18 +++--- src/data/cam_thermo.F90 | 29 +++++----- src/data/cam_thermo_formula.F90 | 24 ++++++++ .../cam_thermo_formula.meta} | 8 +-- src/data/physconst.F90 | 4 -- src/data/registry.xml | 2 +- src/dynamics/mpas/dyn_comp.F90 | 5 ++ src/dynamics/none/dyn_comp.F90 | 3 + src/dynamics/none/dyn_grid.F90 | 57 +++++++++++++++++++ src/dynamics/se/dp_coupling.F90 | 30 +++------- src/dynamics/se/dyn_comp.F90 | 6 +- src/dynamics/tests/dyn_tests_utils.F90 | 6 -- 12 files changed, 130 insertions(+), 62 deletions(-) create mode 100644 src/data/cam_thermo_formula.F90 rename src/{dynamics/tests/dyn_tests_utils.meta => data/cam_thermo_formula.meta} (71%) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 443565eb..0a014c9a 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -541,12 +541,12 @@ end subroutine dry_air_composition_update ! water_composition_update: Update generalized cp or cv depending on dycore !--------------------------------------------------------------------------- !=========================================================================== - subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) - use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array - integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: energy_formula ! energy formula for dynamical core real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) character(len=*), parameter :: subname = 'water_composition_update' @@ -554,9 +554,9 @@ subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) ! update enthalpy or internal energy scaling factor for energy consistency with CAM physics ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module - if (vcoord == vc_moist_pressure) then + if (energy_formula == ENERGY_FORMULA_DYCORE_FV) then ! FV: moist pressure vertical coordinate does not need update. - else if (vcoord == vc_dry_pressure) then + else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then ! SE ! **TEMP** TODO hplin 9/17/24: @@ -566,7 +566,7 @@ subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) 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 (vcoord == vc_height) then + else if (energy_formula == ENERGY_FORMULA_DYCORE_MPAS) then ! MPAS call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), & cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) @@ -575,7 +575,7 @@ subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) else - call endrun(subname//': dycore vcoord not supported') + call endrun(subname//': dycore energy formula not supported') end if end subroutine water_composition_update diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index a776cece..09dff1f1 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -79,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 --------------------- @@ -87,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 -------------------- ! ! @@ -262,7 +263,7 @@ end subroutine cam_thermo_dry_air_update ! !*************************************************************************** ! - subroutine cam_thermo_water_update(mmr, ncol, vcoord, to_dry_factor) + subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor) use air_composition, only: water_composition_update !----------------------------------------------------------------------- ! Update the physics "constants" that vary @@ -270,10 +271,10 @@ subroutine cam_thermo_water_update(mmr, ncol, vcoord, to_dry_factor) real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: vcoord + integer, intent(in) :: energy_formula real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) - call water_composition_update(mmr, ncol, vcoord, to_dry_factor=to_dry_factor) + call water_composition_update(mmr, ncol, energy_formula, to_dry_factor=to_dry_factor) end subroutine cam_thermo_water_update !=========================================================================== @@ -1580,8 +1581,8 @@ 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 air_composition, only: dry_air_species_num use physconst, only: rga, latvap, latice @@ -1600,7 +1601,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & real(kind_phys), intent(in) :: U(:,:) real(kind_phys), intent(in) :: V(:,:) real(kind_phys), intent(in) :: T(:,:) - integer, intent(in) :: vcoord ! vertical coordinate + 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(:,:) @@ -1693,12 +1694,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & ke_vint = 0._kind_phys se_vint = 0._kind_phys select case (vcoord) - case(vc_moist_pressure, vc_dry_pressure) + 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 ', & - 'moist/dry pressure vertical coordinate' + 'FV/SE energy formula' call endrun(subname//': ptop and phis must be present for '// & - 'moist/dry pressure vertical coordinate') + 'FV/SE energy formula') end if po_vint = ptop do kdx = 1, SIZE(tracer, 2) @@ -1714,12 +1715,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & do idx = 1, SIZE(tracer, 1) po_vint(idx) = (phis(idx) * po_vint(idx) * rga) end do - case(vc_height) + case(ENERGY_FORMULA_DYCORE_MPAS) if (.not. present(phis)) then write(iulog, *) subname, ' phis must be present for ', & - 'height-based vertical coordinate' + 'MPAS energy formula' call endrun(subname//': phis must be present for '// & - 'height-based vertical coordinate') + 'MPAS energy formula') end if po_vint = 0._kind_phys do kdx = 1, SIZE(tracer, 2) @@ -1734,8 +1735,8 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & 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 + po_vint + ke_vint diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90 new file mode 100644 index 00000000..56a56cf7 --- /dev/null +++ b/src/data/cam_thermo_formula.F90 @@ -0,0 +1,24 @@ +module cam_thermo_formula + + implicit none + private + save + + ! saves energy formula to use for physics and dynamical core + ! for use in cam_thermo, air_composition and other modules + ! separated into cam_thermo_formula module for clean dependencies + + ! energy_formula options (was vcoord in CAM and stored in dyn_tests_utils) + integer, public, parameter :: ENERGY_FORMULA_DYCORE_FV = 0 ! vc_moist_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_SE = 1 ! vc_dry_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_MPAS = 2 ! vc_height + + !> \section arg_table_cam_thermo_formula Argument Table + !! \htmlinclude cam_thermo_formula.html + ! energy_formula_dycore: energy formula used for dynamical core + ! written by the dynamical core + integer, public :: energy_formula_dycore + ! energy_formula_physics: energy formula used for physics + integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + +end module cam_thermo_formula diff --git a/src/dynamics/tests/dyn_tests_utils.meta b/src/data/cam_thermo_formula.meta similarity index 71% rename from src/dynamics/tests/dyn_tests_utils.meta rename to src/data/cam_thermo_formula.meta index 7bfcbf1c..f8bf04a1 100644 --- a/src/dynamics/tests/dyn_tests_utils.meta +++ b/src/data/cam_thermo_formula.meta @@ -1,16 +1,16 @@ [ccpp-table-properties] - name = dyn_tests_utils + name = cam_thermo_formula type = module [ccpp-arg-table] - name = dyn_tests_utils + name = cam_thermo_formula type = module -[ vc_dycore ] +[ energy_formula_dycore ] standard_name = total_energy_formula_for_dycore units = 1 type = integer dimensions = () -[ vc_physics ] +[ energy_formula_physics ] standard_name = total_energy_formula_for_physics units = 1 type = integer diff --git a/src/data/physconst.F90 b/src/data/physconst.F90 index 774b8f3d..8a528342 100644 --- a/src/data/physconst.F90 +++ b/src/data/physconst.F90 @@ -114,7 +114,6 @@ subroutine physconst_readnl(nlfile) use mpi, only: mpi_real8 use cam_logfile, only: iulog use runtime_obj, only: unset_real - use dyn_tests_utils, only: vc_physics, vc_moist_pressure ! Dummy argument: filepath for file containing namelist input character(len=*), intent(in) :: nlfile @@ -288,9 +287,6 @@ subroutine physconst_readnl(nlfile) ez = omega / sqrt(0.375_kind_phys) - ! set physics vertical coordinate info - vc_physics = vc_moist_pressure - end subroutine physconst_readnl end module physconst diff --git a/src/data/registry.xml b/src/data/registry.xml index daf04535..bda3800d 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -13,10 +13,10 @@ $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 - $SRCROOT/src/dynamics/tests/dyn_tests_utils.meta diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 30e4b2a5..666a6eff 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -184,6 +184,8 @@ end subroutine dyn_readnl ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) + use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS + type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out @@ -200,6 +202,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) nullify(pio_init_file) nullify(pio_topo_file) + ! Set dynamical core energy formula for use in cam_thermo. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS + allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) 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..8e5a2a69 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 @@ -61,6 +62,7 @@ subroutine model_grid_init() use pio, only: PIO_BCAST_ERROR, pio_seterrorhandling use pio, only: pio_get_var, pio_freedecomp use pio, only: pio_read_darray + use pio, only: pio_inq_att use spmd_utils, only: npes, iam use cam_pio_utils, only: cam_pio_handle_error, cam_pio_find_var use cam_pio_utils, only: cam_pio_var_info, pio_subsystem @@ -126,6 +128,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 +162,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 +634,53 @@ 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, var_desc_t + use pio, only: pio_inq_att, 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 + + ! Find which dynamical core is used in and set the energy formulation + ! (also called vc_dycore in CAM) + + type(file_desc_t), intent(inout) :: file + logical, intent(in) :: grid_is_latlon + + ! Local variables + type(var_desc_t) :: vardesc + integer :: ierr, xtype + character(len=*), parameter :: subname = 'find_energy_formula' + + energy_formula_dycore = -1 + + ! Is FV dycore? (has lat lon dimension) + if(grid_is_latlon) then + energy_formula_dycore = ENERGY_FORMULA_DYCORE_FV + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use FV dycore energy formula' + endif + else + ! Is SE dycore? + ierr = pio_inq_att(file, vardesc, 'ne', xtype) + if (ierr == PIO_NOERR) then + ! Has ne property - is SE dycore. + ! if has fv_nphys then is physics grid (ne..pg..), but the energy formulation is the same. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use SE dycore energy formula' + endif + else + ! Is unstructured and is MPAS dycore + ! there are no global attributes to identify MPAS dycore, so this has to do for now. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use MPAS dycore energy formula' + endif + endif + endif + + end subroutine + end module dyn_grid diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 1eca0982..22f8d6b8 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -582,12 +582,12 @@ 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 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 + use physics_types, only: cpairv, rairv, zvirv, cappav use physics_grid, only: columns_on_task use geopotential_temp, only: geopotential_temp_run use static_energy, only: update_dry_static_energy_run @@ -597,7 +597,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 dyn_tests_utils, only: vc_dry_pressure + use cam_thermo_formula,only: ENERGY_FORMULA_DYCORE_SE ! arguments type(runtime_options), intent(in) :: cam_runtime_opts ! Runtime settings object @@ -690,16 +690,9 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! 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) - ! **TEMP** TODO CHECK hplin: check indices to use here - ! in cam6_3_109 after fix in 6_3_127: loop from dry_air_species_num + 1, thermodynamic_active_species_num - ! in cam-sima: factor_array only contains m = ix_q - ! hplin -- to get same answers as CAM-SIMA ix_q would need to be used, - ! but the CAM version appears to be correct. this should be science checked - ! // **TEMP** 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) - ! write(6,*) "hplin dp_coupling m, m_cnst", m, m_cnst - ! m = ix_q do k = 1, nlev do i = 1, pcols ! at this point all q's are dry @@ -746,7 +739,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end do !------------------------------------------------------------ - ! Apply limiters to mixing ratios of major species (waccmx): + ! 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. !------------------------------------------------------------ @@ -811,20 +804,15 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! (note: at this point q is dry) ! call cam_thermo_water_update( & - mmr = const_data_ptr, & ! dry MMR - ncol = pcols, & - vcoord = vc_dry_pressure & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + energy_formula = ENERGY_FORMULA_DYCORE_SE & ) !$omp parallel do num_threads(horz_num_threads) private (k, i) do k = 1, nlev do i = 1, pcols - ! **TEMP** TODO CHECK hplin: CAM and CAM-SIMA version differ - ! CAM version: - ! phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) - - ! CAM-SIMA version (uses constant cappa): - phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa + phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) end do end do diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index c9906b1a..5f8dd2e6 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -566,6 +566,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use dyn_thermo, only: get_molecular_diff_coef_reference !use cam_history, only: addfld, add_default, horiz_only, register_vector_field use gravity_waves_sources, only: gws_init + use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_SE !SE dycore: use prim_advance_mod, only: prim_advance_init @@ -581,7 +582,6 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use control_mod, only: runtype, raytau0, raykrange, rayk0, molecular_diff, nu_top use test_fvm_mapping, only: test_mapping_addfld use control_mod, only: vert_remap_uvTq_alg, vert_remap_tracer_alg - use dyn_tests_utils, only: vc_dycore, vc_dry_pressure ! Dummy arguments: type(runtime_options), intent(in) :: cam_runtime_opts @@ -643,8 +643,8 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) real(r8) :: tau0, krange, otau0, scale real(r8) :: km_sponge_factor_local(nlev+1) !---------------------------------------------------------------------------- - ! Set dynamical core vertical coordinate - vc_dycore = vc_dry_pressure + ! Set dynamical core energy formula for use in cam_thermo. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE ! Now allocate and set condenstate vars allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers diff --git a/src/dynamics/tests/dyn_tests_utils.F90 b/src/dynamics/tests/dyn_tests_utils.F90 index 9ac2a2f1..3a3596b0 100644 --- a/src/dynamics/tests/dyn_tests_utils.F90 +++ b/src/dynamics/tests/dyn_tests_utils.F90 @@ -20,10 +20,4 @@ module dyn_tests_utils integer, parameter :: vc_height = 2 ! Height vertical coord public :: vc_moist_pressure, vc_dry_pressure, vc_height -!> \section arg_table_dyn_tests_utils Argument Table -!! \htmlinclude dyn_tests_utils.html - integer :: vc_dycore ! vertical coordinate of dynamical core - set in dyn_comp.F90 - integer :: vc_physics ! vertical coordinate of physics - set in physconst.F90 - public :: vc_dycore, vc_physics - end module dyn_tests_utils From 38fdbb21c4da0c625f9c964d503b12ad909c3775 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 24 Oct 2024 14:52:30 -0400 Subject: [PATCH 11/41] Read cp_or_cv_dycore and identify dycore information from CAM snapshot --- src/data/cam_thermo_formula.F90 | 13 +++++++++++++ src/data/registry.xml | 1 + src/dynamics/none/dyn_grid.F90 | 16 ++++++++++------ src/physics/utils/phys_comp.F90 | 2 ++ tools/stdnames_to_inputnames_dictionary.xml | 5 +++++ 5 files changed, 31 insertions(+), 6 deletions(-) diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90 index 56a56cf7..c3e2c15b 100644 --- a/src/data/cam_thermo_formula.F90 +++ b/src/data/cam_thermo_formula.F90 @@ -21,4 +21,17 @@ module cam_thermo_formula ! energy_formula_physics: energy formula used for physics integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + ! Public subroutines + public :: cam_thermo_formula_init + +contains + subroutine cam_thermo_formula_init() + use phys_vars_init_check, only: mark_as_initialized + + ! Physics energy formulation is always FV (moist pressure coordinate) + energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + call mark_as_initialized("total_energy_formula_for_physics") + + end subroutine cam_thermo_formula_init + end module cam_thermo_formula diff --git a/src/data/registry.xml b/src/data/registry.xml index bda3800d..d2445ab9 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -428,6 +428,7 @@ allocatable="allocatable"> Enthalpy or internal energy scaling factor for energy consistency horizontal_dimension vertical_layer_dimension + cp_or_cv_dycore diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index 8e5a2a69..d18cbc07 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -637,10 +637,11 @@ end subroutine find_dimension !=========================================================================== subroutine find_energy_formula(file, grid_is_latlon) - use pio, only: file_desc_t, var_desc_t - use pio, only: pio_inq_att, 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 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 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) @@ -649,7 +650,6 @@ subroutine find_energy_formula(file, grid_is_latlon) logical, intent(in) :: grid_is_latlon ! Local variables - type(var_desc_t) :: vardesc integer :: ierr, xtype character(len=*), parameter :: subname = 'find_energy_formula' @@ -663,7 +663,7 @@ subroutine find_energy_formula(file, grid_is_latlon) endif else ! Is SE dycore? - ierr = pio_inq_att(file, vardesc, 'ne', xtype) + 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. @@ -681,6 +681,10 @@ subroutine find_energy_formula(file, grid_is_latlon) endif endif + if(energy_formula_dycore /= -1) then + call mark_as_initialized("total_energy_formula_for_dycore") + endif + end subroutine end module dyn_grid diff --git a/src/physics/utils/phys_comp.F90 b/src/physics/utils/phys_comp.F90 index abe0428a..59ec39ab 100644 --- a/src/physics/utils/phys_comp.F90 +++ b/src/physics/utils/phys_comp.F90 @@ -134,6 +134,7 @@ subroutine phys_init() use physics_grid, only: columns_on_task use vert_coord, only: pver, pverp use cam_thermo, only: cam_thermo_init + use cam_thermo_formula, only: cam_thermo_formula_init use physics_types, only: allocate_physics_types_fields use cam_ccpp_cap, only: cam_ccpp_physics_initialize use cam_ccpp_cap, only: ccpp_physics_suite_part_list @@ -142,6 +143,7 @@ subroutine phys_init() integer :: i_group call cam_thermo_init(columns_on_task, pver, pverp) + call cam_thermo_formula_init() call allocate_physics_types_fields(columns_on_task, pver, pverp, & set_init_val_in=.true., reallocate_in=.false.) diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 0ac72a21..78dc69f2 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -196,6 +196,11 @@ state_tw_cur + + + cp_or_cv_dycore + + RHO From 240ef548bc3573097a739c67a5f3cac671022677 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 24 Oct 2024 19:32:07 -0400 Subject: [PATCH 12/41] Set wv_idx in air_composition total_hours_in_debugging_one_line = 6 --- src/data/air_composition.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 0a014c9a..f450c182 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -229,7 +229,7 @@ subroutine air_composition_init() ! !************************************************************************ ! - ! add prognostic components of dry air + ! add prognostic components of air ! !************************************************************************ ! @@ -313,6 +313,7 @@ subroutine air_composition_init() ! case(wv_stdname) !water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water call air_species_info(wv_stdname, ix, mw) + wv_idx = ix ! set water species index for use in get_hydrostatic_energy thermodynamic_active_species_idx(icnst) = ix thermodynamic_active_species_cp (icnst) = cpwv thermodynamic_active_species_cv (icnst) = cv3 / mw From bef25fbc46642a65a6ec390e83ffdb537424db67 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 25 Oct 2024 12:03:33 -0400 Subject: [PATCH 13/41] Update atmos_phys external --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index f35d1e2e..5f4d639c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 996a6d5 + fxtag = 952ebdd 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 996a6d51..952ebddf 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 996a6d517128496285cd3bcae12aed3be11da622 +Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 From 8b8a8e089b578c7c336098c2831f462ec24811c9 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 25 Oct 2024 12:09:18 -0400 Subject: [PATCH 14/41] Fix merge conflicts --- cime_config/cam_autogen.py | 10 ---------- src/physics/utils/cam_constituents.F90 | 8 -------- 2 files changed, 18 deletions(-) diff --git a/cime_config/cam_autogen.py b/cime_config/cam_autogen.py index f5467a6a..3d5ed312 100644 --- a/cime_config/cam_autogen.py +++ b/cime_config/cam_autogen.py @@ -614,7 +614,6 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, # Copy to_be_ccppized utility modules to the build directory, # as SIMA cam_constituents depends on them. -<<<<<<< HEAD # Note: to_be_ccppized utility modules to be removed once functionality is migrated # to SIMA or CCPPized in atmospheric_physics. atm_phys_to_be_ccppized_dir = os.path.join(atm_phys_top_dir, "to_be_ccppized") @@ -631,15 +630,6 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, for to_be_ccppized_file in atm_phys_to_be_ccppized_files: shutil.copy(to_be_ccppized_file, physics_blddir) # end for -======= - atm_phys_to_be_ccppized_dir = os.path.join(atm_phys_top_dir, "to_be_ccppized") - - # Check that the directory exists and copy to build directory - if os.path.isdir(atm_phys_to_be_ccppized_dir): - atm_phys_to_be_ccppized_files = glob.glob(os.path.join(atm_phys_to_be_ccppized_dir, "*.F90")) - for to_be_ccppized_file in atm_phys_to_be_ccppized_files: - shutil.copy(to_be_ccppized_file, physics_blddir) ->>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) # end if if do_gen_ccpp or do_gen_nl: diff --git a/src/physics/utils/cam_constituents.F90 b/src/physics/utils/cam_constituents.F90 index 3883aba5..19329153 100644 --- a/src/physics/utils/cam_constituents.F90 +++ b/src/physics/utils/cam_constituents.F90 @@ -292,11 +292,7 @@ subroutine const_get_index(name, cindex, abort, warning, caller) use cam_abortutils, only: endrun use cam_logfile, only: iulog use phys_vars_init_check, only: std_name_len -<<<<<<< HEAD use string_utils, only: stringify -======= - use string_utils, only: to_str ->>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) ! Get the index of a constituent with standard name, . ! Setting optional argument to .false. returns control to @@ -323,11 +319,7 @@ subroutine const_get_index(name, cindex, abort, warning, caller) call ccpp_const_get_idx(const_props, name, cindex, errmsg, errcode) if (errcode /= 0) then -<<<<<<< HEAD call endrun(subname//"Error "//stringify((/errcode/))//": "// & -======= - call endrun(subname//"Error "//to_str(errcode)//": "// & ->>>>>>> aa86287 (New const_get_index logic without cam_ccpp_cap dependency) trim(errmsg), file=__FILE__, line=__LINE__) endif From 53cba060f460bef4636249d701401927d5e5059b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 25 Oct 2024 12:11:15 -0400 Subject: [PATCH 15/41] Implements CCPPized check_energy in SIMA. Squashed commit of the following: commit 8b8a8e089b578c7c336098c2831f462ec24811c9 Fix merge conflicts commit bef25fbc46642a65a6ec390e83ffdb537424db67 Update atmos_phys external commit 240ef548bc3573097a739c67a5f3cac671022677 Set wv_idx in air_composition total_hours_in_debugging_one_line = 6 commit 38fdbb21c4da0c625f9c964d503b12ad909c3775 Read cp_or_cv_dycore and identify dycore information from CAM snapshot commit 8fca3a42b6db828e185ec4df51665e1dcae8eae4 Update vcoord to energy_formula commit 0d16be936ddfbb1e0c3a223671fadf15e8c955cd New const_get_index logic without cam_ccpp_cap dependency commit e13e7ea0de0592bd2517cfcac32ea0bf3fabc6a6 Fix build (part 1); update standard names and atmos_phys external commit c015e28ab9ea5e590230c05df7c5411d7adf8b6a Update ncar-physics external commit 2372f7fc2fc1c71900f57ce1131c2b323f8803c3 Fix for b4b to CAM-SIMA (w/ History); add notes on modifications and/or divergence from CAM and/or -SIMA commit 4ce842786b8b64374057058dd22a21dc3d2576eb Add gmean_mod to src/utils including support infra in physics_grid.F90 commit 7493a1dfa70ecd6e837308e087a0cff4dc9cd2c2 Fix build issues; update factor in get_cp call to 1/dp commit 041cdfb250dc0088f87d6c1ecb391c015e6b4b57 Add is_first_timestep registry commit 333ad4eed8520f92c61a6731bfd589d2f873f097 Initial port of updates to energy budget (port of CAM #761) into CAM-SIMA commit 86489702461192ccb13ae7cff9d13f0db004812b Initial implementation of cp_to_cv_dycore/cam_thermo_water_update from stash Needs work to split into cam_thermo_water_update. This is WIP code --- .gitmodules | 4 +- src/control/cam_comp.F90 | 7 + src/data/air_composition.F90 | 151 +++++++++--- src/data/cam_thermo.F90 | 204 ++++++++++------ src/data/cam_thermo_formula.F90 | 37 +++ src/data/cam_thermo_formula.meta | 17 ++ src/data/registry.xml | 15 ++ src/dynamics/mpas/dyn_comp.F90 | 5 + src/dynamics/none/dyn_comp.F90 | 3 + src/dynamics/none/dyn_grid.F90 | 61 +++++ src/dynamics/se/dp_coupling.F90 | 122 +++++++--- src/dynamics/se/dycore/prim_advance_mod.F90 | 6 +- src/dynamics/se/dyn_comp.F90 | 3 + src/dynamics/utils/dyn_thermo.F90 | 33 +-- src/physics/ncar_ccpp | 2 +- src/physics/utils/phys_comp.F90 | 2 + src/physics/utils/physics_grid.F90 | 52 ++++ src/utils/gmean_mod.F90 | 256 ++++++++++++++++++++ tools/stdnames_to_inputnames_dictionary.xml | 5 + 19 files changed, 819 insertions(+), 166 deletions(-) create mode 100644 src/data/cam_thermo_formula.F90 create mode 100644 src/data/cam_thermo_formula.meta create mode 100644 src/utils/gmean_mod.F90 diff --git a/.gitmodules b/.gitmodules index 94610494..5f4d639c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -19,8 +19,8 @@ fxDONOTUSEurl = https://github.com/MPAS-Dev/MPAS-Model.git [submodule "ncar-physics"] path = src/physics/ncar_ccpp - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = e95c172d7a5a0ebf054f420b08416228e211baa3 + url = https://github.com/jimmielin/atmospheric_physics + fxtag = 952ebdd fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 538bd641..27fd57a7 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -29,6 +29,7 @@ module cam_comp use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys use physics_types, only: calday + use physics_types, only: is_first_timestep use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -151,6 +152,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & dtime_phys = 0.0_r8 call mark_as_initialized('timestep_for_physics') + is_first_timestep = .true. + call mark_as_initialized('is_first_timestep') + call init_pio_subsystem() ! Initializations using data passed from coupler. @@ -263,6 +267,9 @@ subroutine cam_timestep_init() use phys_comp, only: phys_timestep_init use stepon, only: stepon_timestep_init + ! Update timestep flags in physics state + is_first_timestep = is_first_step() + !---------------------------------------------------------- ! First phase of dynamics (at least couple from dynamics to physics) ! Return time-step for physics from dynamics. diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index f202ae97..f450c182 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -1,4 +1,5 @@ -! air_composition module defines major species of the atmosphere and manages the physical properties that are dependent on the composition of air +! air_composition module defines major species of the atmosphere and manages +! the physical properties that are dependent on the composition of air module air_composition use ccpp_kinds, only: kind_phys @@ -6,13 +7,16 @@ module air_composition use runtime_obj, only: unset_real, unset_int use phys_vars_init_check, only: std_name_len use physics_types, only: cpairv, rairv, cappav, mbarv, zvirv + use physics_types, only: cp_or_cv_dycore implicit none private save public :: air_composition_init - public :: air_composition_update + public :: dry_air_composition_update + public :: water_composition_update + ! get_cp_dry: (generalized) heat capacity for dry air public :: get_cp_dry ! get_cp: (generalized) heat capacity @@ -225,7 +229,7 @@ subroutine air_composition_init() ! !************************************************************************ ! - ! add prognostic components of dry air + ! add prognostic components of air ! !************************************************************************ ! @@ -309,6 +313,7 @@ subroutine air_composition_init() ! case(wv_stdname) !water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water call air_species_info(wv_stdname, ix, mw) + wv_idx = ix ! set water species index for use in get_hydrostatic_energy thermodynamic_active_species_idx(icnst) = ix thermodynamic_active_species_cp (icnst) = cpwv thermodynamic_active_species_cv (icnst) = cv3 / mw @@ -510,26 +515,71 @@ end subroutine air_composition_init !=========================================================================== !----------------------------------------------------------------------- - ! air_composition_update: Update the physics "constants" that vary + ! dry_air_composition_update: Update the physics "constants" that vary !------------------------------------------------------------------------- !=========================================================================== - subroutine air_composition_update(mmr, ncol, to_moist_factor) + subroutine dry_air_composition_update(mmr, ncol, to_dry_factor) - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + !(mmr = dry mixing ratio, if not, use to_dry_factor to convert!) + real(kind_phys), intent(in) :: mmr(:,:,:) ! mixing ratios for species dependent dry air integer, intent(in) :: ncol ! number of columns - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, & - rairv(:ncol, :), fact=to_moist_factor) + rairv(:ncol, :), fact=to_dry_factor) call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - cpairv(:ncol,:), fact=to_moist_factor) + cpairv(:ncol,:), fact=to_dry_factor) call get_mbarv(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - mbarv(:ncol,:), fact=to_moist_factor) + mbarv(:ncol,:), fact=to_dry_factor) cappav(:ncol,:) = rairv(:ncol,:) / cpairv(:ncol,:) - end subroutine air_composition_update + end subroutine dry_air_composition_update + + !=========================================================================== + !--------------------------------------------------------------------------- + ! water_composition_update: Update generalized cp or cv depending on dycore + !--------------------------------------------------------------------------- + !=========================================================================== + subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: energy_formula ! energy formula for dynamical core + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + character(len=*), parameter :: subname = 'water_composition_update' + + ! update enthalpy or internal energy scaling factor for energy consistency with CAM physics + ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module + + if (energy_formula == ENERGY_FORMULA_DYCORE_FV) then + ! FV: moist pressure vertical coordinate does not need update. + else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then + ! SE + + ! **TEMP** TODO hplin 9/17/24: + ! for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) + ! (whereas CAM only allocates 1-indexed) I subset it here. This should be verified during code + ! review. + call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & + factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx(1:), & + cpdry=cpairv(:ncol,:)) + else if (energy_formula == ENERGY_FORMULA_DYCORE_MPAS) then + ! MPAS + call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), & + cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) + + ! internal energy coefficient for MPAS + ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) + cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) + else + call endrun(subname//': dycore energy formula not supported') + end if + + end subroutine water_composition_update !=========================================================================== !*************************************************************************** @@ -639,27 +689,33 @@ end subroutine get_cp_dry_2hd ! !*************************************************************************** ! - subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) use cam_abortutils, only: endrun use string_utils, only: to_str ! Dummy arguments ! tracer: Tracer array + ! + ! factor not present then tracer must be dry mixing ratio + ! if factor present tracer*factor must be dry mixing ratio + ! real(kind_phys), intent(in) :: tracer(:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) + ! dp: if provided then tracer is mass not mixing ratio + real(kind_phys), optional, intent(in) :: factor(:,:) ! active_species_idx_dycore: array of indices for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:) ! Local variables integer :: qdx, itrac real(kind_phys) :: sum_species(SIZE(cp, 1), SIZE(cp, 2)) real(kind_phys) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2)) - real(kind_phys) :: factor(SIZE(cp, 1), SIZE(cp, 2)) + real(kind_phys) :: factor_local(SIZE(cp, 1), SIZE(cp, 2)) integer :: idx_local(thermodynamic_active_species_num) character(len=*), parameter :: subname = 'get_cp_1hd: ' @@ -675,28 +731,41 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) idx_local = thermodynamic_active_species_idx end if - if (present(dp_dry)) then - factor = 1.0_kind_phys / dp_dry + if (present(factor)) then + factor_local = factor else - factor = 1.0_kind_phys + factor_local = 1.0_kind_phys end if + sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) - sum_species(:,:) = sum_species(:,:) + & - (tracer(:,:,itrac) * factor(:,:)) + sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:)) end do - ! Get heat capacity at constant pressure (Cp) for dry air: - call get_cp_dry(tracer, idx_local, sum_cp, fact=factor) + if (dry_air_species_num == 0) then + ! **TMP** TODO CHECK hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + ! and dry_air_species_num is == 0 in CAM-SIMA as of 10/1/24 which means this clause will + ! change answers (?) -- previously in CAM-SIMA there was only a call to get_cp_dry + ! and not the two IF-clauses here. + sum_cp = thermodynamic_active_species_cp(0) + else if (present(cpdry)) then + ! + ! if cpdry is known don't recompute + ! + sum_cp = cpdry + else + ! Get heat capacity at constant pressure (Cp) for dry air: + call get_cp_dry(tracer, idx_local, sum_cp, fact=factor_local) + end if ! Add water species to Cp: do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) sum_cp(:,:) = sum_cp(:,:) + & - (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * & - factor(:,:)) + (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * factor_local(:,:)) end do + if (inv_cp) then cp = sum_species / sum_cp else @@ -707,7 +776,7 @@ end subroutine get_cp_1hd !=========================================================================== - subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_2hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) ! Version of get_cp for arrays that have a second horizontal index use cam_abortutils, only: endrun use string_utils, only: to_str @@ -715,14 +784,15 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) ! Dummy arguments ! tracer: Tracer array real(kind_phys), intent(in) :: tracer(:,:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:,:) + real(kind_phys), optional, intent(in) :: factor(:,:,:) ! active_species_idx_dycore: array of indicies for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:,:) ! Local variables integer :: jdx @@ -730,11 +800,17 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) character(len=*), parameter :: subname = 'get_cp_2hd: ' do jdx = 1, SIZE(cp, 2) - if (present(dp_dry)) then - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & - dp_dry=dp_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + if (present(factor).and.present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) + else if (present(factor)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + else if (present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) else - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& active_species_idx_dycore=active_species_idx_dycore) end if end do @@ -843,9 +919,10 @@ end subroutine get_R_dry_2hd ! !*************************************************************************** ! - subroutine get_R_1hd(tracer, active_species_idx, R, fact) + subroutine get_R_1hd(tracer, active_species_idx, R, fact, Rdry) use cam_abortutils, only: endrun use string_utils, only: to_str + use physconst, only: rair ! Dummy arguments ! tracer: !tracer array @@ -856,6 +933,7 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) real(kind_phys), intent(out) :: R(:, :) ! fact: optional factor for converting tracer to dry mixing ratio real(kind_phys), optional, intent(in) :: fact(:, :) + real(kind_phys), optional, intent(in) :: Rdry(:, :) ! Local variables integer :: qdx, itrac @@ -874,12 +952,19 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) call endrun(subname//"SIZE mismatch in dimension 2 "// & to_str(SIZE(fact, 2))//' /= '//to_str(SIZE(factor, 2))) end if - call get_R_dry(tracer, active_species_idx, R, fact=fact) factor = fact(:,:) else - call get_R_dry(tracer, active_species_idx, R) factor = 1.0_kind_phys end if + + if (dry_air_species_num == 0) then + R = rair + else if (present(Rdry)) then + R = Rdry + else + call get_R_dry(tracer, active_species_idx, R, fact=factor) + end if + idx_local = active_species_idx sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num @@ -934,7 +1019,7 @@ end subroutine get_R_2hd !************************************************************************************************************************* ! subroutine get_mbarv_1hd(tracer, active_species_idx, mbarv_in, fact) - use physconst, only: mwdry, rair, cpair + use physconst, only: mwdry real(kind_phys), intent(in) :: tracer(:,:,:) !tracer array integer, intent(in) :: active_species_idx(:) !index of active species in tracer real(kind_phys), intent(out) :: mbarv_in(:,:) !molecular weight of dry air diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 59dd1c83..09dff1f1 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -33,8 +33,10 @@ module cam_thermo ! cam_thermo_init: Initialize constituent dependent properties public :: cam_thermo_init - ! cam_thermo_update: Update constituent dependent properties - public :: cam_thermo_update + ! cam_thermo_dry_air_update: Update dry air composition dependent properties + public :: cam_thermo_dry_air_update + ! cam_thermo_water_update: Update water dependent properties + public :: cam_thermo_water_update ! get_enthalpy: enthalpy quantity = dp*cp*T public :: get_enthalpy ! get_virtual_temp: virtual temperature @@ -77,6 +79,7 @@ module cam_thermo ! mixing_ratio options integer, public, parameter :: DRY_MIXING_RATIO = 1 integer, public, parameter :: MASS_MIXING_RATIO = 2 + !> \section arg_table_cam_thermo Argument Table !! \htmlinclude cam_thermo.html !--------------- Variables below here are for WACCM-X --------------------- @@ -85,7 +88,7 @@ module cam_thermo ! kmcnd: molecular conductivity J m-1 s-1 K-1 real(kind_phys), public, protected, allocatable :: kmcnd(:,:) - !------------- Variables for consistent themodynamics -------------------- + !------------- Variables for consistent thermodynamics -------------------- ! ! @@ -208,51 +211,71 @@ end subroutine cam_thermo_init !=========================================================================== + ! !*************************************************************************** ! - ! cam_thermo_update: update species dependent constants for physics + ! cam_thermo_dry_air_update: update dry air species dependent constants for physics ! !*************************************************************************** ! - subroutine cam_thermo_update(mmr, T, ncol, update_thermo_variables, to_moist_factor) - use air_composition, only: air_composition_update, update_zvirv + subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) + use air_composition, only: dry_air_composition_update + use air_composition, only: update_zvirv use string_utils, only: to_str - !----------------------------------------------------------------------- - ! Update the physics "constants" that vary - !------------------------------------------------------------------------- - !------------------------------Arguments---------------------------------- + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) + real(kind_phys), intent(in) :: T(:,:) ! temperature + integer, intent(in) :: ncol ! number of columns + logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables + ! false: do not calculate composition-dependent thermo variables + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr moist convert - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array - real(kind_phys), intent(in) :: T(:,:) ! temperature - integer, intent(in) :: ncol ! number of columns - logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables - ! false: do not calculate composition-dependent thermo variables - - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - ! - !---------------------------Local storage------------------------------- - real(kind_phys):: sponge_factor(SIZE(mmr, 2)) - character(len=*), parameter :: subname = 'cam_thermo_update: ' + ! Local vars + real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) + character(len=*), parameter :: subname = 'cam_thermo_dry_air_update: ' if (.not. update_thermo_variables) then return end if - if (present(to_moist_factor)) then - if (SIZE(to_moist_factor, 1) /= ncol) then - call endrun(subname//'DIM 1 of to_moist_factor is'//to_str(SIZE(to_moist_factor,1))//'but should be'//to_str(ncol)) - end if + if (present(to_dry_factor)) then + if (SIZE(to_dry_factor, 1) /= ncol) then + call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) + end if end if + sponge_factor = 1.0_kind_phys - call air_composition_update(mmr, ncol, to_moist_factor=to_moist_factor) + call dry_air_composition_update(mmr, ncol, to_dry_factor=to_dry_factor) call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:), & - kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_moist_factor, & + kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_dry_factor, & active_species_idx_dycore=thermodynamic_active_species_idx) + + ! Calculate zvirv for WACCM-X. call update_zvirv() - end subroutine cam_thermo_update + end subroutine cam_thermo_dry_air_update + + ! + !*************************************************************************** + ! + ! cam_thermo_water_update: update water species dependent constants for physics + ! + !*************************************************************************** + ! + subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor) + use air_composition, only: water_composition_update + !----------------------------------------------------------------------- + ! Update the physics "constants" that vary + !------------------------------------------------------------------------- + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: energy_formula + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + call water_composition_update(mmr, ncol, energy_formula, to_dry_factor=to_dry_factor) + end subroutine cam_thermo_water_update !=========================================================================== @@ -1554,28 +1577,32 @@ end subroutine cam_thermo_calc_kappav_2hd ! !*************************************************************************** ! - subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & - vcoord, ps, phis, z_mid, dycore_idx, qidx, te, se, ke, & - wv, H2O, liq, ice) + subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & + cp_or_cv, U, V, T, vcoord, ptop, phis, z_mid, dycore_idx, qidx, & + te, se, po, ke, wv, H2O, liq, ice) + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS use cam_logfile, only: iulog - use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure use air_composition, only: wv_idx - use physconst, only: gravit, latvap, latice + use air_composition, only: dry_air_species_num + use physconst, only: rga, latvap, latice ! Dummy arguments ! tracer: tracer mixing ratio + ! + ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry real(kind_phys), intent(in) :: tracer(:,:,:) + logical, intent(in) :: moist_mixing_ratio ! pdel: pressure level thickness - real(kind_phys), intent(in) :: pdel(:,:) + real(kind_phys), intent(in) :: pdel_in(:,:) ! cp_or_cv: dry air heat capacity under constant pressure or ! constant volume (depends on vcoord) real(kind_phys), intent(in) :: cp_or_cv(:,:) real(kind_phys), intent(in) :: U(:,:) real(kind_phys), intent(in) :: V(:,:) real(kind_phys), intent(in) :: T(:,:) - integer, intent(in) :: vcoord ! vertical coordinate - real(kind_phys), intent(in), optional :: ps(:) + integer, intent(in) :: vcoord !REMOVECAM - vcoord or energy formula to use + real(kind_phys), intent(in), optional :: ptop(:) real(kind_phys), intent(in), optional :: phis(:) real(kind_phys), intent(in), optional :: z_mid(:,:) ! dycore_idx: use dycore index for thermodynamic active species @@ -1588,8 +1615,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & real(kind_phys), intent(out), optional :: te (:) ! KE: vertically integrated kinetic energy real(kind_phys), intent(out), optional :: ke (:) - ! SE: vertically integrated internal+geopotential energy + ! SE: vertically integrated enthalpy (pressure coordinate) + ! or internal energy (z coordinate) real(kind_phys), intent(out), optional :: se (:) + ! PO: vertically integrated PHIS term (pressure coordinate) + ! or potential energy (z coordinate) + real(kind_phys), intent(out), optional :: po (:) ! WV: vertically integrated water vapor real(kind_phys), intent(out), optional :: wv (:) ! liq: vertically integrated liquid @@ -1599,10 +1630,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & ! Local variables real(kind_phys) :: ke_vint(SIZE(tracer, 1)) ! Vertical integral of KE - real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of SE + real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of enthalpy or internal energy + real(kind_phys) :: po_vint(SIZE(tracer, 1)) ! Vertical integral of PHIS or potential energy real(kind_phys) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv real(kind_phys) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq real(kind_phys) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice + real(kind_phys) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) ! moist pressure level thickness real(kind_phys) :: latsub ! latent heat of sublimation integer :: ierr @@ -1649,78 +1682,96 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & wvidx = wv_idx end if + if (moist_mixing_ratio) then + pdel = pdel_in + else + pdel = pdel_in + do qdx = dry_air_species_num+1, thermodynamic_active_species_num + pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx)) + end do + end if + + ke_vint = 0._kind_phys + se_vint = 0._kind_phys select case (vcoord) - case(vc_moist_pressure, vc_dry_pressure) - if ((.not. present(ps)) .or. (.not. present(phis))) then - write(iulog, *) subname, ' ps and phis must be present for ', & - 'moist/dry pressure vertical coordinate' - call endrun(subname//': ps and phis must be present for '// & - 'moist/dry pressure vertical coordinate') + case(ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE) + if (.not. present(ptop).or. (.not. present(phis))) then + write(iulog, *) subname, ' ptop and phis must be present for ', & + 'FV/SE energy formula' + call endrun(subname//': ptop and phis must be present for '// & + 'FV/SE energy formula') end if - ke_vint = 0._kind_phys - se_vint = 0._kind_phys - wv_vint = 0._kind_phys + po_vint = ptop do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * & - 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit) + 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2)) * rga se_vint(idx) = se_vint(idx) + (T(idx, kdx) * & - cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit) - wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & - pdel(idx, kdx) / gravit) + cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga) + po_vint(idx) = po_vint(idx)+pdel(idx, kdx) + end do end do do idx = 1, SIZE(tracer, 1) - se_vint(idx) = se_vint(idx) + (phis(idx) * ps(idx) / gravit) + po_vint(idx) = (phis(idx) * po_vint(idx) * rga) end do - case(vc_height) - if (.not. present(z_mid)) then - write(iulog, *) subname, & - ' z_mid must be present for height vertical coordinate' - call endrun(subname//': z_mid must be present for height '// & - 'vertical coordinate') + case(ENERGY_FORMULA_DYCORE_MPAS) + if (.not. present(phis)) then + write(iulog, *) subname, ' phis must be present for ', & + 'MPAS energy formula' + call endrun(subname//': phis must be present for '// & + 'MPAS energy formula') end if - ke_vint = 0._kind_phys - se_vint = 0._kind_phys - wv_vint = 0._kind_phys + po_vint = 0._kind_phys do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * & - 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit) + 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) * rga) se_vint(idx) = se_vint(idx) + (T(idx, kdx) * & - cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit) + cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga) ! z_mid is height above ground - se_vint(idx) = se_vint(idx) + (z_mid(idx, kdx) + & - phis(idx) / gravit) * pdel(idx, kdx) - wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & - pdel(idx, kdx) / gravit) + po_vint(idx) = po_vint(idx) + (z_mid(idx, kdx) + & + phis(idx) * rga) * pdel(idx, kdx) end do end do case default - write(iulog, *) subname, ' vertical coordinate not supported: ', vcoord - call endrun(subname//': vertical coordinate not supported') + write(iulog, *) subname, ' energy formula not supported: ', vcoord + call endrun(subname//': energy formula not supported') end select if (present(te)) then - te = se_vint + ke_vint + te = se_vint + po_vint + ke_vint end if if (present(se)) then se = se_vint end if + if (present(po)) then + po = po_vint + end if if (present(ke)) then ke = ke_vint end if - if (present(wv)) then - wv = wv_vint - end if ! ! vertical integral of total liquid water ! + if (.not.moist_mixing_ratio) then + pdel = pdel_in! set pseudo density to dry + end if + + wv_vint = 0._kind_phys + do kdx = 1, SIZE(tracer, 2) + do idx = 1, SIZE(tracer, 1) + wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * & + pdel(idx, kdx) * rga) + end do + end do + if (present(wv)) wv = wv_vint + liq_vint = 0._kind_phys do qdx = 1, thermodynamic_active_species_liq_num do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) - liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * & - tracer(idx, kdx, species_liq_idx(qdx)) / gravit) + liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * & + tracer(idx, kdx, species_liq_idx(qdx)) * rga) end do end do end do @@ -1734,7 +1785,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & do kdx = 1, SIZE(tracer, 2) do idx = 1, SIZE(tracer, 1) ice_vint(idx) = ice_vint(idx) + (pdel(idx, kdx) * & - tracer(idx, kdx, species_ice_idx(qdx)) / gravit) + tracer(idx, kdx, species_ice_idx(qdx)) * rga) end do end do end do @@ -1762,7 +1813,6 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, & end select end if deallocate(species_idx, species_liq_idx, species_ice_idx) - end subroutine get_hydrostatic_energy_1hd !=========================================================================== diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90 new file mode 100644 index 00000000..c3e2c15b --- /dev/null +++ b/src/data/cam_thermo_formula.F90 @@ -0,0 +1,37 @@ +module cam_thermo_formula + + implicit none + private + save + + ! saves energy formula to use for physics and dynamical core + ! for use in cam_thermo, air_composition and other modules + ! separated into cam_thermo_formula module for clean dependencies + + ! energy_formula options (was vcoord in CAM and stored in dyn_tests_utils) + integer, public, parameter :: ENERGY_FORMULA_DYCORE_FV = 0 ! vc_moist_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_SE = 1 ! vc_dry_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_MPAS = 2 ! vc_height + + !> \section arg_table_cam_thermo_formula Argument Table + !! \htmlinclude cam_thermo_formula.html + ! energy_formula_dycore: energy formula used for dynamical core + ! written by the dynamical core + integer, public :: energy_formula_dycore + ! energy_formula_physics: energy formula used for physics + integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + + ! Public subroutines + public :: cam_thermo_formula_init + +contains + subroutine cam_thermo_formula_init() + use phys_vars_init_check, only: mark_as_initialized + + ! Physics energy formulation is always FV (moist pressure coordinate) + energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + call mark_as_initialized("total_energy_formula_for_physics") + + end subroutine cam_thermo_formula_init + +end module cam_thermo_formula diff --git a/src/data/cam_thermo_formula.meta b/src/data/cam_thermo_formula.meta new file mode 100644 index 00000000..f8bf04a1 --- /dev/null +++ b/src/data/cam_thermo_formula.meta @@ -0,0 +1,17 @@ +[ccpp-table-properties] + name = cam_thermo_formula + type = module + +[ccpp-arg-table] + name = cam_thermo_formula + type = module +[ energy_formula_dycore ] + standard_name = total_energy_formula_for_dycore + units = 1 + type = integer + dimensions = () +[ energy_formula_physics ] + standard_name = total_energy_formula_for_physics + units = 1 + type = integer + dimensions = () diff --git a/src/data/registry.xml b/src/data/registry.xml index 41b31416..d2445ab9 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -13,6 +13,7 @@ $SRCROOT/src/physics/utils/tropopause_climo_read.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta + $SRCROOT/src/data/cam_thermo_formula.meta $SRCROOT/src/data/ref_pres.meta $SRCROOT/src/dynamics/utils/vert_coord.meta $SRCROOT/src/dynamics/utils/hycoef.meta @@ -251,6 +252,12 @@ horizontal_dimension tw_cur state_tw_cur + + + flag indicating if is first timestep of initial run + horizontal_dimension vertical_layer_dimension zvir + + Enthalpy or internal energy scaling factor for energy consistency + horizontal_dimension vertical_layer_dimension + cp_or_cv_dycore + and set the energy formulation + ! (also called vc_dycore in CAM) + + type(file_desc_t), intent(inout) :: file + logical, intent(in) :: grid_is_latlon + + ! Local variables + integer :: ierr, xtype + character(len=*), parameter :: subname = 'find_energy_formula' + + energy_formula_dycore = -1 + + ! Is FV dycore? (has lat lon dimension) + if(grid_is_latlon) then + energy_formula_dycore = ENERGY_FORMULA_DYCORE_FV + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use FV dycore energy formula' + endif + else + ! Is SE dycore? + ierr = pio_inq_att(file, pio_global, 'ne', xtype) + if (ierr == PIO_NOERR) then + ! Has ne property - is SE dycore. + ! if has fv_nphys then is physics grid (ne..pg..), but the energy formulation is the same. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use SE dycore energy formula' + endif + else + ! Is unstructured and is MPAS dycore + ! there are no global attributes to identify MPAS dycore, so this has to do for now. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS + if(masterproc) then + write(iulog, *) subname, ': Null dycore will use MPAS dycore energy formula' + endif + endif + endif + + if(energy_formula_dycore /= -1) then + call mark_as_initialized("total_energy_formula_for_dycore") + endif + + end subroutine + end module dyn_grid diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 8b56e9d9..22f8d6b8 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -582,18 +582,22 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) use cam_constituents, only: const_qmin use runtime_obj, only: wv_stdname use physics_types, only: lagrangian_vertical - use physconst, only: cpair, gravit, zvir, cappa - use cam_thermo, only: cam_thermo_update - use physics_types, only: cpairv, rairv, zvirv + use physconst, only: cpair, gravit, zvir + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use air_composition, only: thermodynamic_active_species_num + use air_composition, only: thermodynamic_active_species_idx + use air_composition, only: dry_air_species_num + use physics_types, only: cpairv, rairv, zvirv, cappav use physics_grid, only: columns_on_task use geopotential_temp, only: geopotential_temp_run use static_energy, only: update_dry_static_energy_run use qneg, only: qneg_run -! use check_energy, only: check_energy_timestep_init +! use check_energy_chng, only: check_energy_chng_timestep_init use hycoef, only: hyai, ps0 use shr_vmath_mod, only: shr_vmath_log use shr_kind_mod, only: shr_kind_cx use dyn_comp, only: ixo, ixo2, ixh, ixh2 + use cam_thermo_formula,only: ENERGY_FORMULA_DYCORE_SE ! arguments type(runtime_options), intent(in) :: cam_runtime_opts ! Runtime settings object @@ -607,7 +611,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !constituent properties pointer type(ccpp_constituent_prop_ptr_t), pointer :: const_prop_ptr(:) - integer :: m, i, k + integer :: m, i, k, m_cnst integer :: ix_q !Needed for "geopotential_temp" CCPP scheme @@ -622,6 +626,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) real(r8) :: mmrSum_O_O2_H ! Sum of mass mixing ratios for O, O2, and H real(r8), parameter :: mmrMin=1.e-20_r8 ! lower limit of o2, o, and h mixing ratios real(r8), parameter :: N2mmrMin=1.e-6_r8 ! lower limit of o2, o, and h mixing ratios + real(r8), parameter :: H2lim=6.e-5_r8 ! H2 limiter: 10x global H2 MMR (Roble, 1995) !---------------------------------------------------------------------------- ! Nullify pointers @@ -683,14 +688,23 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end do ! wet pressure variables (should be removed from physics!) + factor_array(:,:) = 1.0_kind_phys + !$omp parallel do num_threads(horz_num_threads) private (k, i, m_cnst) + do m_cnst = dry_air_species_num + 1, thermodynamic_active_species_num + ! include all water species in the factor array. + m = thermodynamic_active_species_idx(m_cnst) + do k = 1, nlev + do i = 1, pcols + ! at this point all q's are dry + factor_array(i,k) = factor_array(i,k) + const_data_ptr(i,k,m) + end do + end do + end do !$omp parallel do num_threads(horz_num_threads) private (k, i) - do k=1,nlev - do i=1, pcols - ! to be consistent with total energy formula in physic's check_energy module only - ! include water vapor in moist dp - factor_array(i,k) = 1._kind_phys+const_data_ptr(i,k,ix_q) - phys_state%pdel(i,k) = phys_state%pdeldry(i,k)*factor_array(i,k) + do k = 1, nlev + do i = 1, pcols + phys_state%pdel(i,k) = phys_state%pdeldry(i,k) * factor_array(i,k) end do end do @@ -721,31 +735,12 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) do k = 1, nlev do i = 1, pcols phys_state%rpdel(i,k) = 1._kind_phys/phys_state%pdel(i,k) - phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa - end do - end do - - ! all tracers (including moisture) are in dry mixing ratio units - ! physics expect water variables moist - factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) - - !$omp parallel do num_threads(horz_num_threads) private (m, k, i) - do m=1, num_advected - do k = 1, nlev - do i=1, pcols - !This should ideally check if a constituent is a wet - !mixing ratio or not, but until that is working properly - !in the CCPP framework just check for the water species status - !instead, which is all that CAM physics requires: - if (const_is_water_species(m)) then - const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) - end if - end do end do end do !------------------------------------------------------------ - ! Ensure O2 + O + H (N2) mmr greater than one. + ! Apply limiters to mixing ratios of major species (WACCMX): + ! Ensure N2 = 1 - (O2 + O + H) mmr is greater than 0 ! Check for unusually large H2 values and set to lower value. !------------------------------------------------------------ if (cam_runtime_opts%waccmx_option() == 'ionosphere' .or. & @@ -769,8 +764,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) endif - if(const_data_ptr(i,k,ixh2) .gt. 6.e-5_r8) then - const_data_ptr(i,k,ixh2) = 6.e-5_r8 + if(const_data_ptr(i,k,ixh2) > H2lim) then + const_data_ptr(i,k,ixh2) = H2lim endif end do @@ -789,11 +784,62 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! Compute molecular viscosity(kmvis) and conductivity(kmcnd). ! Update zvirv registry variable; calculated for WACCM-X. !----------------------------------------------------------------------------- + ! **TEMP** TODO CHECK hplin: CAM has this if-clause for dry_air_species_num > 0 + ! or otherwise uses zvirv = zvir. CAM-SIMA previously did not have this, and + ! instead has a switch for update_thermodynamic_variables. Check if we still want + ! this if-clause or change it to something else. + if (dry_air_species_num > 0) then + call cam_thermo_dry_air_update( & + mmr = const_data_ptr, & ! dry MMR + T = phys_state%t, & + ncol = pcols, & + update_thermo_variables = cam_runtime_opts%update_thermodynamic_variables() & + ) + else + zvirv(:,:) = zvir + end if + + ! + ! update cp_or_cv_dycore in module air_composition. + ! (note: at this point q is dry) + ! + call cam_thermo_water_update( & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + energy_formula = ENERGY_FORMULA_DYCORE_SE & + ) - call cam_thermo_update(const_data_ptr, phys_state%t, pcols, & - cam_runtime_opts%update_thermodynamic_variables()) + !$omp parallel do num_threads(horz_num_threads) private (k, i) + do k = 1, nlev + do i = 1, pcols + phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) + end do + end do + + ! ========= Q is dry ^^^ ---- Q is moist vvv ========= ! + + ! + ! CAM physics expects that: water tracers (including moisture) are moist; the rest dry mixing ratio + ! at this point Q is converted to moist. + ! + factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) + + !$omp parallel do num_threads(horz_num_threads) private (m, k, i) + do m = 1, num_advected + do k = 1, nlev + do i = 1, pcols + ! This should ideally check if a constituent is a wet + ! mixing ratio or not, but until that is working properly + ! in the CCPP framework just check for the water species status + ! instead, which is all that CAM physics requires: + if (const_is_water_species(m)) then + const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) + end if + end do + end do + end do - !Call geopotential_temp CCPP scheme: + ! Call geopotential_temp CCPP scheme: call geopotential_temp_run(pver, lagrangian_vertical, pver, 1, & pverp, 1, num_advected, phys_state%lnpint, & phys_state%pint, phys_state%pmid, phys_state%pdel, & @@ -810,7 +856,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !Remove once check_energy scheme exists in CAMDEN: #if 0 ! Compute energy and water integrals of input state - call check_energy_timestep_init(phys_state, phys_tend, pbuf_chnk) + ! call check_energy_chng_timestep_init(...) #endif end subroutine derived_phys_dry diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90 index 1341b9f4..cf88d7f9 100644 --- a/src/dynamics/se/dycore/prim_advance_mod.F90 +++ b/src/dynamics/se/dycore/prim_advance_mod.F90 @@ -1674,8 +1674,12 @@ subroutine calc_tot_energy_dynamics(elem,fvm,nets,nete,tl,tl_qdp,outfld_name_suf call get_dp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),MASS_MIXING_RATIO,thermodynamic_active_species_idx_dycore,& elem(ie)%state%dp3d(:,:,:,tl),pdel,ps=ps,ptop=hyai(1)*ps0) call get_cp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),& - .false.,cp,dp_dry=elem(ie)%state%dp3d(:,:,:,tl),& + .false.,cp,factor=1.0_r8/elem(ie)%state%dp3d(:,:,:,tl),& active_species_idx_dycore=thermodynamic_active_species_idx_dycore) + + ! TODO: need to port cam6_3_109 changes to total energy using get_hydrostatic_energy + ! https://github.com/ESCOMP/CAM/pull/761/files#diff-946bde17289e2f42e43e64413610aa11d102deda8b5199ddaa5b71e67e5d517a + do k = 1, nlev do j=1,np do i = 1, np diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index ab52d91c..5f8dd2e6 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -566,6 +566,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use dyn_thermo, only: get_molecular_diff_coef_reference !use cam_history, only: addfld, add_default, horiz_only, register_vector_field use gravity_waves_sources, only: gws_init + use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_SE !SE dycore: use prim_advance_mod, only: prim_advance_init @@ -642,6 +643,8 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) real(r8) :: tau0, krange, otau0, scale real(r8) :: km_sponge_factor_local(nlev+1) !---------------------------------------------------------------------------- + ! Set dynamical core energy formula for use in cam_thermo. + energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE ! Now allocate and set condenstate vars allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers diff --git a/src/dynamics/utils/dyn_thermo.F90 b/src/dynamics/utils/dyn_thermo.F90 index c4c4723c..9b465031 100644 --- a/src/dynamics/utils/dyn_thermo.F90 +++ b/src/dynamics/utils/dyn_thermo.F90 @@ -61,7 +61,7 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Declare local variables: real(kind_phys), allocatable :: tracer_phys(:,:,:,:) real(kind_phys), allocatable :: cp_phys(:,:,:) - real(kind_phys), allocatable :: dp_dry_phys(:,:,:) + real(kind_phys), allocatable :: factor_phys(:,:,:) !check_allocate variables: integer :: iret !allocate status integer @@ -70,11 +70,16 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Check if kinds are different: if (kind_phys == kind_dyn) then - !The dynamics and physics kind is the same, so just call the physics - !routine directly: - call get_cp_phys(tracer,inv_cp,cp, & - dp_dry=dp_dry, & - active_species_idx_dycore=active_species_idx_dycore) + ! The dynamics and physics kind is the same, so just call the physics + ! routine directly: + if(present(dp_dry)) then + call get_cp_phys(tracer,inv_cp,cp, & + factor=1.0_kind_phys/dp_dry, & + active_species_idx_dycore=active_species_idx_dycore) + else + call get_cp_phys(tracer,inv_cp,cp, & + active_species_idx_dycore=active_species_idx_dycore) + endif else @@ -95,18 +100,18 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) !Allocate and set optional variables: if (present(dp_dry)) then - allocate(dp_dry_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret) + allocate(factor_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret) call check_allocate(iret, subname, & - 'dp_dry_phys', & + 'factor_phys', & file=__FILE__, line=__LINE__) !Set optional local variable: - dp_dry_phys = real(dp_dry, kind_phys) + factor_phys = 1.0_kind_phys/real(dp_dry, kind_phys) end if - !Call physics routine using local vriables with matching kinds: + !Call physics routine using local variables with matching kinds: call get_cp_phys(tracer_phys,inv_cp,cp_phys, & - dp_dry=dp_dry_phys, & + factor=factor_phys, & active_species_idx_dycore=active_species_idx_dycore) !Set output variables back to dynamics kind: @@ -116,8 +121,8 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore) deallocate(tracer_phys) deallocate(cp_phys) - if (allocated(dp_dry_phys)) then - deallocate(dp_dry_phys) + if (allocated(factor_phys)) then + deallocate(factor_phys) end if @@ -957,7 +962,7 @@ subroutine get_rho_dry(tracer,temp,ptop,dp_dry,tracer_mass,& end if - !Call physics routine using local vriables with matching kinds: + !Call physics routine using local variables with matching kinds: call get_rho_dry_phys(tracer_phys,temp_phys, & ptop_phys, dp_dry_phys,tracer_mass, & rho_dry=rho_dry_phys, & diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index e95c172d..952ebddf 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit e95c172d7a5a0ebf054f420b08416228e211baa3 +Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 diff --git a/src/physics/utils/phys_comp.F90 b/src/physics/utils/phys_comp.F90 index abe0428a..59ec39ab 100644 --- a/src/physics/utils/phys_comp.F90 +++ b/src/physics/utils/phys_comp.F90 @@ -134,6 +134,7 @@ subroutine phys_init() use physics_grid, only: columns_on_task use vert_coord, only: pver, pverp use cam_thermo, only: cam_thermo_init + use cam_thermo_formula, only: cam_thermo_formula_init use physics_types, only: allocate_physics_types_fields use cam_ccpp_cap, only: cam_ccpp_physics_initialize use cam_ccpp_cap, only: ccpp_physics_suite_part_list @@ -142,6 +143,7 @@ subroutine phys_init() integer :: i_group call cam_thermo_init(columns_on_task, pver, pverp) + call cam_thermo_formula_init() call allocate_physics_types_fields(columns_on_task, pver, pverp, & set_init_val_in=.true., reallocate_in=.false.) diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90 index 39fc0f99..cf28f98d 100644 --- a/src/physics/utils/physics_grid.F90 +++ b/src/physics/utils/physics_grid.F90 @@ -39,8 +39,10 @@ module physics_grid public :: get_rlat_p ! latitude of a physics column in radians public :: get_rlon_p ! longitude of a physics column in radians public :: get_area_p ! area of a physics column in radians squared + public :: get_wght_p ! weight of a physics column in radians squared public :: get_rlat_all_p ! latitudes of physics cols on task (radians) public :: get_rlon_all_p ! longitudes of physics cols on task (radians) + public :: get_wght_all_p ! weights of physics cols on task public :: get_dyn_col_p ! dynamics local blk number and blk offset(s) public :: global_index_p ! global column index of a physics column public :: local_index_p ! local column index of a physics column @@ -475,6 +477,26 @@ end function get_area_p !======================================================================== + real(r8) function get_wght_p(index) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + ! weight of a physics column in radians squared + + ! Dummy argument + integer, intent(in) :: index + ! Local variables + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_p' + + ! Check that input is valid: + call check_phys_input(subname, index) + + get_wght_p = phys_columns(index)%weight + + end function get_wght_p + + !======================================================================== + subroutine get_rlat_all_p(rlatdim, rlats) use cam_logfile, only: iulog use cam_abortutils, only: endrun @@ -535,6 +557,36 @@ end subroutine get_rlon_all_p !======================================================================== + subroutine get_wght_all_p(wghtdim, wghts) + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + !----------------------------------------------------------------------- + ! + ! get_wght_all_p: Return all weights on task. + ! + !----------------------------------------------------------------------- + ! Dummy Arguments + integer, intent(in) :: wghtdim ! declared size of output array + real(r8), intent(out) :: wghts(wghtdim) ! array of weights + + ! Local variables + integer :: index ! loop index + character(len=128) :: errmsg + character(len=*), parameter :: subname = 'get_wght_all_p: ' + + !----------------------------------------------------------------------- + + ! Check that input is valid: + call check_phys_input(subname, wghtdim) + + do index = 1, wghtdim + wghts(index) = phys_columns(index)%weight + end do + + end subroutine get_wght_all_p + + !======================================================================== + subroutine get_dyn_col_p(index, blk_num, blk_ind) use cam_logfile, only: iulog use cam_abortutils, only: endrun diff --git a/src/utils/gmean_mod.F90 b/src/utils/gmean_mod.F90 new file mode 100644 index 00000000..ca6dd6d4 --- /dev/null +++ b/src/utils/gmean_mod.F90 @@ -0,0 +1,256 @@ +module gmean_mod + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Perform global mean calculations for energy conservation and other checks. + ! + ! Method: + ! Reproducible (scalable): + ! Convert to fixed point (integer representation) to enable + ! reproducibility when using MPI collectives. + ! If error checking is on (via setting reprosum_diffmax > 0 and + ! reprosum_recompute = .true. in user_nl_cpl), shr_reprosum_calc will + ! check the accuracy of its computation with a fast but + ! non-reproducible algorithm. If any error is reported, report + ! the difference and the expected sum and abort run (call endrun) + ! + ! gmean_mod in to_be_ccppized is different from the CAM version and + ! has chunk support removed. + ! + ! + !----------------------------------------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use physics_grid, only: pcols => columns_on_task + use perf_mod, only: t_startf, t_stopf + use cam_logfile, only: iulog + + implicit none + private + + public :: gmean ! compute global mean of 2D fields on physics decomposition + + interface gmean + module procedure gmean_arr + module procedure gmean_scl + end interface gmean + + private :: gmean_fixed_repro + private :: gmean_float_norepro + + ! Set do_gmean_tests to .true. to run a gmean challenge test + logical, private :: do_gmean_tests = .false. + +CONTAINS + + ! + !======================================================================== + ! + + subroutine gmean_arr (arr, arr_gmean, nflds) + use shr_strconvert_mod, only: toString + use spmd_utils, only: masterproc + use cam_abortutils, only: endrun + use shr_reprosum_mod, only: shr_reprosum_reldiffmax, shr_reprosum_recompute, shr_reprosum_tolExceeded + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! + ! Method is to call shr_reprosum_calc (called from gmean_fixed_repro) + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols, nflds) + real(r8), intent(out) :: arr_gmean(nflds) ! global means + ! + ! Local workspace + ! + real(r8) :: rel_diff(2, nflds) + integer :: ifld ! field index + integer :: num_err + logical :: write_warning + ! + !----------------------------------------------------------------------- + ! + call t_startf('gmean_arr') + call t_startf ('gmean_fixed_repro') + call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds) + call t_stopf ('gmean_fixed_repro') + + ! check that "fast" reproducible sum is accurate enough. If not, calculate + ! using old method + write_warning = masterproc + num_err = 0 + if (shr_reprosum_tolExceeded('gmean', nflds, write_warning, & + iulog, rel_diff)) then + if (shr_reprosum_recompute) then + do ifld = 1, nflds + if (rel_diff(1, ifld) > shr_reprosum_reldiffmax) then + call gmean_float_norepro(arr(:,ifld), arr_gmean(ifld), ifld) + num_err = num_err + 1 + end if + end do + end if + end if + call t_stopf('gmean_arr') + if (num_err > 0) then + call endrun('gmean: '//toString(num_err)//' reprosum errors found') + end if + + end subroutine gmean_arr + + ! + !======================================================================== + ! + + subroutine gmean_scl (arr, gmean) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of a single field in "arr" in the physics grid + ! + !----------------------------------------------------------------------- + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + ! Input array, chunked + real(r8), intent(out):: gmean ! global means + ! + ! Local workspace + ! + integer, parameter :: nflds = 1 + real(r8) :: gmean_array(nflds) + real(r8) :: array(pcols, nflds) + integer :: ncols, lchnk + + array(:ncols, 1) = arr(:ncols) + call gmean_arr(array, gmean_array, nflds) + gmean = gmean_array(1) + + end subroutine gmean_scl + + ! + !======================================================================== + ! + + subroutine gmean_float_norepro(arr, repro_sum, index) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of in the physics chunked + ! decomposition using a fast but non-reproducible algorithm. + ! Log that value along with the value computed by + ! shr_reprosum_calc () + ! + !----------------------------------------------------------------------- + + use physconst, only: pi + use spmd_utils, only: masterproc, masterprocid, mpicom + use mpi, only: mpi_real8, mpi_sum + use physics_grid, only: get_wght_p + ! + ! Arguments + ! + real(r8), intent(in) :: arr(pcols) + real(r8), intent(in) :: repro_sum ! Value computed by reprosum + integer, intent(in) :: index ! Index of field in original call + ! + ! Local workspace + ! + integer :: icol + integer :: ierr + real(r8) :: wght + real(r8) :: check + real(r8) :: check_sum + real(r8), parameter :: pi4 = 4.0_r8 * pi + + ! + !----------------------------------------------------------------------- + ! + ! Calculate and print out non-reproducible value + check = 0.0_r8 + do icol = 1, pcols + wght = get_wght_p(icol) + check = check + arr(icol) * wght + end do + call MPI_reduce(check, check_sum, 1, mpi_real8, mpi_sum, & + masterprocid, mpicom, ierr) + + ! normalization + check_sum = check_sum / pi4 + + if (masterproc) then + write(iulog, '(a,i0,2(a,e20.13e2))') 'gmean(', index, ') = ', & + check_sum, ', reprosum reported ', repro_sum + end if + + end subroutine gmean_float_norepro + + ! + !======================================================================== + ! + subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Compute the global mean of each field in "arr" in the physics grid + ! with a reproducible yet scalable implementation + ! based on a fixed-point algorithm. + ! + !----------------------------------------------------------------------- + use spmd_utils, only: mpicom + use physics_grid, only: get_wght_all_p + use physics_grid, only: ngcols_p => num_global_phys_cols + use physconst, only: pi + use shr_reprosum_mod, only: shr_reprosum_calc + ! + ! Arguments + ! + integer, intent(in) :: nflds ! number of fields + real(r8), intent(in) :: arr(pcols,nflds) + ! arr_gmean: output global sums + real(r8), intent(out) :: arr_gmean(nflds) + ! rel_diff: relative and absolute differences from shr_reprosum_calc + real(r8), intent(out) :: rel_diff(2, nflds) + ! + ! Local workspace + ! + integer :: icol, ifld ! column, field indices + + real(r8) :: wght(pcols) ! integration weights + real(r8), allocatable :: xfld(:,:) ! weighted summands + ! + !----------------------------------------------------------------------- + ! + allocate(xfld(pcols, nflds)) + + ! pre-weight summands + call get_wght_all_p(pcols, wght) + + do ifld = 1, nflds + do icol = 1, pcols + xfld(icol, ifld) = arr(icol, ifld) * wght(icol) + end do + end do + + ! call fixed-point algorithm + call shr_reprosum_calc ( & + arr = xfld, & + arr_gsum = arr_gmean, & + nsummands = pcols, & ! # of local summands + dsummands = pcols, & ! declared first dimension of arr. + nflds = nflds, & + commid = mpicom, & + rel_diff = rel_diff & + ) + + deallocate(xfld) + ! final normalization + arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi) + + end subroutine gmean_fixed_repro + +end module gmean_mod diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 0ac72a21..78dc69f2 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -196,6 +196,11 @@ state_tw_cur + + + cp_or_cv_dycore + + RHO From a4d3866a9c9f1c45436cd9d524973339eb3c6e4e Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 4 Nov 2024 14:34:47 -0500 Subject: [PATCH 16/41] Add adiabatic scheme and check_energy_fix --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 5f4d639c..577a92de 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 952ebdd + fxtag = 69621f8 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 952ebddf..69621f82 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 +Subproject commit 69621f822fac2b181b8c77d4ec933a613f875e18 From e0450f986507b8ac46359781b527b4e7450e4bdc Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 4 Nov 2024 15:22:35 -0500 Subject: [PATCH 17/41] Fix index underflow in get_hydrostatic_energy; update ncar-physics external to latest --- .gitmodules | 2 +- src/data/cam_thermo.F90 | 4 ++-- src/physics/ncar_ccpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.gitmodules b/.gitmodules index 577a92de..d0c0cbb1 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 69621f8 + fxtag = 87a45418 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 09dff1f1..3b1dfbfd 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -1666,12 +1666,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & 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 diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 69621f82..87a45418 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 69621f822fac2b181b8c77d4ec933a613f875e18 +Subproject commit 87a45418a647a08fc9626b1f04c8d6c0e460bfc0 From 6e5873e0a30add62a47b0b890a301153e73bb6a0 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 5 Nov 2024 11:54:06 -0500 Subject: [PATCH 18/41] Update with ne3np4 inic files --- cime_config/namelist_definition_cam.xml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index ba8fbe40..0d8b77dc 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -104,6 +104,10 @@ ${DIN_LOC_ROOT}/atm/cam/inic/gaus/cami_0000-01-01_8x16_L26_c030228.nc ${DIN_LOC_ROOT}/atm/cam/inic/gaus/cami_0000-09-01_8x16_L26_c030918.nc ${DIN_LOC_ROOT}/atm/cam/inic/gaus/cami_0000-01-01_8x16_L30_c090102.nc + ${DIN_LOC_ROOT}/atm/cam/inic/se/cam6_QPC6_topo_ne3pg3_mg37_L32_01-01-31_c221214.nc + ${DIN_LOC_ROOT}/atm/cam/inic/se/CESM2.F2000climo.ne3np4.cam.i.0003-09-01-00000.nc + ${DIN_LOC_ROOT}/atm/cam/inic/se/cami_0000-01-01_ne3np4_L30_c120315.nc + ${DIN_LOC_ROOT}/atm/cam/inic/se/cami_0000-01-01_ne3np4_L26_c120525.nc ${DIN_LOC_ROOT}/atm/cam/inic/homme/cami-mam3_0000-01_ne5np4_L30.140707.nc ${DIN_LOC_ROOT}/atm/cam/inic/se/ape_topo_cam4_ne16np4_L26_c171020.nc ${DIN_LOC_ROOT}/atm/cam/inic/se/ape_topo_cam4_ne16np4_L30_c171020.nc From 768040ee438b4ba5a839ce6bdb69981db201887e Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 5 Nov 2024 11:57:02 -0500 Subject: [PATCH 19/41] Now mark energy consistency variables in registry as initialized in dycore initialization --- src/data/registry.xml | 7 ------- src/dynamics/mpas/dyn_comp.F90 | 4 +++- src/dynamics/se/dyn_comp.F90 | 12 ++++++++++++ 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index d2445ab9..176e3272 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -210,7 +210,6 @@ type="real" kind="kind_phys" allocatable="pointer"> horizontal_dimension - te_ini_phys state_te_ini_phys horizontal_dimension - te_cur_phys state_te_cur_phys horizontal_dimension - te_ini_dyn state_te_ini_dyn horizontal_dimension - te_cur_dyn state_te_cur_dyn horizontal_dimension - tw_ini state_tw_ini horizontal_dimension - tw_cur state_tw_cur Enthalpy or internal energy scaling factor for energy consistency horizontal_dimension vertical_layer_dimension - cp_or_cv_dycore diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 666a6eff..a484d7ec 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -184,7 +184,8 @@ end subroutine dyn_readnl ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) - use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS + use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS + use phys_vars_init_check, only: mark_as_initialized type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in @@ -204,6 +205,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Set dynamical core energy formula for use in cam_thermo. energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS + call mark_as_initialized("total_energy_formula_for_dycore") allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 5f8dd2e6..59c92b55 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -567,6 +567,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) !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 phys_vars_init_check, only: mark_as_initialized !SE dycore: use prim_advance_mod, only: prim_advance_init @@ -645,6 +646,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) !---------------------------------------------------------------------------- ! 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") ! Now allocate and set condenstate vars allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers @@ -1855,6 +1857,16 @@ 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_for_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_of_initial_state_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") + call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") end subroutine read_inidat From 03100c1591689c2223321d6a359bc380763f6951 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 5 Nov 2024 18:09:29 -0500 Subject: [PATCH 20/41] Fix remaining issues in adiabatic scheme; add teout to registry and initialize mark --- .gitmodules | 2 +- src/data/registry.xml | 16 ++++++++++++++++ src/dynamics/none/dyn_grid.F90 | 11 +++++++++++ src/dynamics/se/dyn_comp.F90 | 1 + src/physics/ncar_ccpp | 2 +- 5 files changed, 30 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index d0c0cbb1..45bebfd8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 87a45418 + fxtag = 003ece7 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/data/registry.xml b/src/data/registry.xml index 176e3272..c476248c 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -210,6 +210,7 @@ type="real" kind="kind_phys" allocatable="pointer"> horizontal_dimension + te_ini_phys state_te_ini_phys horizontal_dimension + te_cur_phys state_te_cur_phys horizontal_dimension + te_ini_dyn state_te_ini_dyn horizontal_dimension + te_cur_dyn state_te_cur_dyn horizontal_dimension + 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 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_at_end_of_physics tendency_of_air_temperature_due_to_model_physics @@ -422,6 +437,7 @@ allocatable="allocatable"> Enthalpy or internal energy scaling factor for energy consistency horizontal_dimension vertical_layer_dimension + cp_or_cv_dycore diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index d18cbc07..07729c30 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -685,6 +685,17 @@ subroutine find_energy_formula(file, grid_is_latlon) call mark_as_initialized("total_energy_formula_for_dycore") endif + ! Mark other energy variables calculated by check_energy_timestep_init + ! here since it will always run when required + call mark_as_initialized("specific_heat_for_air_used_in_dycore") + call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_physics_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") + call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") + call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics") + end subroutine end module dyn_grid diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 59c92b55..0701ae15 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -1867,6 +1867,7 @@ subroutine read_inidat(dyn_in) call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics") end subroutine read_inidat diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 87a45418..003ece7f 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 87a45418a647a08fc9626b1f04c8d6c0e460bfc0 +Subproject commit 003ece7f3a88cef889fb72392600e6b0702cfcaf From 3a1a5b915f3060622251ede675684a46ebd36580 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 8 Nov 2024 11:18:37 -0500 Subject: [PATCH 21/41] Update ncar-physics external for dycore_energy_consistency_adjust --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 45bebfd8..96271203 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 003ece7 + fxtag = e87c44d 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 003ece7f..e87c44d1 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 003ece7f3a88cef889fb72392600e6b0702cfcaf +Subproject commit e87c44d14ceb7c2fdda788af60f5eb53ff047f83 From 7a706dae77e6cfc7d20f5fc0922c9c57c8a25e2a Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 Nov 2024 12:42:13 -0500 Subject: [PATCH 22/41] Update atmospheric_physics external; standard names to address review comments --- .gitmodules | 2 +- src/data/registry.xml | 6 +++--- src/dynamics/mpas/dyn_comp.F90 | 14 ++++++++++++++ src/dynamics/none/dyn_grid.F90 | 4 ++-- src/dynamics/se/dyn_comp.F90 | 4 ++-- src/physics/ncar_ccpp | 2 +- tools/stdnames_to_inputnames_dictionary.xml | 2 +- 7 files changed, 24 insertions(+), 10 deletions(-) diff --git a/.gitmodules b/.gitmodules index 96271203..e2e358c0 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = e87c44d + fxtag = 7b46000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/data/registry.xml b/src/data/registry.xml index c476248c..2c91606a 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -253,7 +253,7 @@ tw_cur state_tw_cur Total energy using dynamical core formula at the end of physics timestep @@ -333,7 +333,7 @@ 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_at_end_of_physics + vertically_integrated_total_energy_at_end_of_physics_timestep tendency_of_air_temperature_due_to_model_physics @@ -432,7 +432,7 @@ zvir Enthalpy or internal energy scaling factor for energy consistency diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index a484d7ec..b55224b7 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -839,6 +839,20 @@ subroutine mark_variable_as_initialized() do i = 1, num_advected call mark_as_initialized(trim(adjustl(const_name(i)))) end do + + call mark_as_initialized('specific_heat_of_dry_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_of_initial_state_using_physics_energy_formula') + call mark_as_initialized('vertically_integrated_total_energy_of_current_state_using_physics_energy_formula') + call mark_as_initialized('vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula') + call mark_as_initialized('vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula') + call mark_as_initialized('vertically_integrated_water_vapor_and_condensed_water_of_initial_state') + call mark_as_initialized('vertically_integrated_water_vapor_and_condensed_water_of_current_state') + call mark_as_initialized('vertically_integrated_total_energy_at_end_of_physics_timestep') + + end subroutine mark_variable_as_initialized ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`. diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index a01d3e0f..587769eb 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -687,14 +687,14 @@ subroutine find_energy_formula(file, grid_is_latlon) ! Mark other energy variables calculated by check_energy_timestep_init ! here since it will always run when required - call mark_as_initialized("specific_heat_for_air_used_in_dycore") + call mark_as_initialized("specific_heat_of_dry_air_used_in_dycore") call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_physics_energy_formula") call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_physics_energy_formula") call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula") call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") - call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep") end subroutine find_energy_formula diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 0701ae15..35633239 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -1857,7 +1857,7 @@ 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_for_air_used_in_dycore") + call mark_as_initialized("specific_heat_of_dry_air_used_in_dycore") ! These energy variables are calculated by check_energy_timestep_init ! but need to be marked here @@ -1867,7 +1867,7 @@ subroutine read_inidat(dyn_in) call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") - call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics") + call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep") end subroutine read_inidat diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index e87c44d1..7b460000 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit e87c44d14ceb7c2fdda788af60f5eb53ff047f83 +Subproject commit 7b46000093aa3e4d2ab90b9bfadf8a06a1d87957 diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 78dc69f2..ba5161e3 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -196,7 +196,7 @@ state_tw_cur - + cp_or_cv_dycore From 11e9a200506f29bde2c54985a88b1edaffb6d947 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 Nov 2024 13:27:45 -0500 Subject: [PATCH 23/41] Update nstep from model state; add FADIAB physics-suite --- .gitmodules | 2 +- cime_config/config_component.xml | 4 ++-- src/control/cam_comp.F90 | 24 +++++++++++------------- src/physics/ncar_ccpp | 2 +- 4 files changed, 15 insertions(+), 17 deletions(-) diff --git a/.gitmodules b/.gitmodules index e2e358c0..fe9ec66b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 7b46000 + fxtag = 1003f81 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index c93a9de8..5551d778 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -159,8 +159,8 @@ -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 27fd57a7..a656e18c 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -29,7 +29,7 @@ module cam_comp use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys use physics_types, only: calday - use physics_types, only: is_first_timestep + use 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,12 +149,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') - - is_first_timestep = .true. - call mark_as_initialized('is_first_timestep') - call init_pio_subsystem() ! Initializations using data passed from coupler. @@ -170,12 +164,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') @@ -269,6 +271,7 @@ subroutine cam_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) @@ -516,10 +519,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() @@ -542,7 +541,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/physics/ncar_ccpp b/src/physics/ncar_ccpp index 7b460000..1003f81c 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 7b46000093aa3e4d2ab90b9bfadf8a06a1d87957 +Subproject commit 1003f81c12f722ba3c68bcc89f9dbf9bf0beca4c From 800bf8b591b3e79254c2c851bb92c7998fc99f49 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 Nov 2024 21:30:48 -0500 Subject: [PATCH 24/41] Update standard names from review comments; add flag for energy consistency adjustment --- .gitmodules | 2 +- src/data/cam_thermo.F90 | 4 +-- src/data/registry.xml | 29 +++++++++++++-------- src/dynamics/mpas/dyn_comp.F90 | 20 +++++++++----- src/dynamics/none/dyn_grid.F90 | 18 ++++++++----- src/dynamics/se/dyn_comp.F90 | 20 +++++++++----- src/physics/ncar_ccpp | 2 +- tools/stdnames_to_inputnames_dictionary.xml | 14 +++++----- 8 files changed, 66 insertions(+), 43 deletions(-) diff --git a/.gitmodules b/.gitmodules index fe9ec66b..ca7b944f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 1003f81 + fxtag = a53ff9b fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 3b1dfbfd..b6f0f568 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -1595,8 +1595,8 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & logical, intent(in) :: moist_mixing_ratio ! pdel: pressure level thickness real(kind_phys), intent(in) :: pdel_in(:,:) - ! cp_or_cv: dry air heat capacity under constant pressure or - ! constant volume (depends on vcoord) + ! 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(:,:) diff --git a/src/data/registry.xml b/src/data/registry.xml index 2c91606a..b9955806 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -205,7 +205,7 @@ zi state_zi @@ -213,7 +213,7 @@ te_ini_phys state_te_ini_phys @@ -221,7 +221,7 @@ te_cur_phys state_te_cur_phys @@ -229,7 +229,7 @@ te_ini_dyn state_te_ini_dyn @@ -237,7 +237,7 @@ te_cur_dyn state_te_cur_dyn @@ -245,7 +245,7 @@ tw_ini state_tw_ini @@ -260,6 +260,13 @@ horizontal_dimension 0.0 + + flag indicating if dynamical core energy not consistent with CAM physics and to perform adjustment of temperature and temperature tendency + .false. + .true. + 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 @@ -432,7 +439,7 @@ zvir Enthalpy or internal energy scaling factor for energy consistency diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index b55224b7..660dc496 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -185,6 +185,7 @@ end subroutine dyn_readnl ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS + use physics_types, only: dycore_energy_consistency_adjust use phys_vars_init_check, only: mark_as_initialized type(runtime_options), intent(in) :: cam_runtime_opts @@ -207,6 +208,11 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS 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") + allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) @@ -840,16 +846,16 @@ subroutine mark_variable_as_initialized() call mark_as_initialized(trim(adjustl(const_name(i)))) end do - call mark_as_initialized('specific_heat_of_dry_air_used_in_dycore') + 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_of_initial_state_using_physics_energy_formula') - call mark_as_initialized('vertically_integrated_total_energy_of_current_state_using_physics_energy_formula') - call mark_as_initialized('vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula') - call mark_as_initialized('vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula') - call mark_as_initialized('vertically_integrated_water_vapor_and_condensed_water_of_initial_state') - call mark_as_initialized('vertically_integrated_water_vapor_and_condensed_water_of_current_state') + 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') diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index 587769eb..a859fb16 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -641,6 +641,7 @@ subroutine find_energy_formula(file, grid_is_latlon) 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 @@ -658,6 +659,7 @@ subroutine find_energy_formula(file, grid_is_latlon) ! 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 @@ -668,6 +670,7 @@ subroutine find_energy_formula(file, grid_is_latlon) ! 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 @@ -675,6 +678,7 @@ subroutine find_energy_formula(file, grid_is_latlon) ! 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 @@ -687,13 +691,13 @@ subroutine find_energy_formula(file, grid_is_latlon) ! 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_dry_air_used_in_dycore") - call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_physics_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_physics_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") - call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") - call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") + 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 diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 35633239..4f70ae2c 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -567,6 +567,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) !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: @@ -648,6 +649,11 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) 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 call check_allocate(iret, subname, 'cnst_name_gll(qsize)', & @@ -1857,16 +1863,16 @@ 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_dry_air_used_in_dycore") + 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_of_initial_state_using_physics_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_physics_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula") - call mark_as_initialized("vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula") - call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_initial_state") - call mark_as_initialized("vertically_integrated_water_vapor_and_condensed_water_of_current_state") + 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/physics/ncar_ccpp b/src/physics/ncar_ccpp index 1003f81c..a53ff9bd 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 1003f81c12f722ba3c68bcc89f9dbf9bf0beca4c +Subproject commit a53ff9bd853c7519488219972e14ac4ee469297a diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index ba5161e3..dbd84b99 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -160,43 +160,43 @@ 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 From 916cb9f07bf74a2921418db831d2b164484bda34 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 Nov 2024 21:43:38 -0500 Subject: [PATCH 25/41] rga -> gravit --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index ca7b944f..bff58e03 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = a53ff9b + fxtag = 5e52537 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 a53ff9bd..5e525373 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit a53ff9bd853c7519488219972e14ac4ee469297a +Subproject commit 5e52537390e9fc65b898e76d28ae8148f171dba3 From b7abae1e8ac05009e9c69ca625e3568f7c0d5ef8 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 Nov 2024 22:26:56 -0500 Subject: [PATCH 26/41] gmean diagnostic updates --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index bff58e03..0cdd8f76 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 5e52537 + fxtag = cef429f 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 5e525373..cef429fe 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 5e52537390e9fc65b898e76d28ae8148f171dba3 +Subproject commit cef429feb6e15569270576edbcbf088111f4bb16 From 6105ca8687c6f84bfa9785a6a27cad0e311322bd Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 15 Nov 2024 12:54:18 -0500 Subject: [PATCH 27/41] Address review comments (update external) --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 0cdd8f76..227d8951 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = cef429f + fxtag = 7b188e35 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 cef429fe..7b188e35 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit cef429feb6e15569270576edbcbf088111f4bb16 +Subproject commit 7b188e35fe26247289d3bdbbbdd902cde382765a From f9997077cc3fe9ab7ac3e7675b7482a77804ef88 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 10 Dec 2024 15:28:33 -0500 Subject: [PATCH 28/41] Address review comments --- src/data/air_composition.F90 | 4 ---- src/data/registry.xml | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 4e87d2fc..833b8e4b 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -744,10 +744,6 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpd end do if (dry_air_species_num == 0) then - ! **TMP** TODO CHECK hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA - ! and dry_air_species_num is == 0 in CAM-SIMA as of 10/1/24 which means this clause will - ! change answers (?) -- previously in CAM-SIMA there was only a call to get_cp_dry - ! and not the two IF-clauses here. sum_cp = thermodynamic_active_species_cp(0) else if (present(cpdry)) then ! diff --git a/src/data/registry.xml b/src/data/registry.xml index 06a5758f..d449bf09 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -443,7 +443,7 @@ standard_name="specific_heat_of_air_used_in_dycore" units="J kg-1 K-1" type="real" kind="kind_phys" allocatable="allocatable"> - Enthalpy or internal energy scaling factor for energy consistency + 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 From 3d2fd910b6a980399b6c88e5638e57e13901a025 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 11 Dec 2024 19:38:16 -0500 Subject: [PATCH 29/41] Address review comments --- cime_config/config_component.xml | 1 - src/data/air_composition.F90 | 18 +++++++-------- src/data/cam_thermo.F90 | 20 +++++++++++++---- src/data/cam_thermo_formula.F90 | 6 +++-- src/data/registry.xml | 4 ++-- src/dynamics/none/dyn_grid.F90 | 1 - src/dynamics/se/dp_coupling.F90 | 1 + src/physics/utils/physics_grid.F90 | 23 -------------------- src/utils/gmean_mod.F90 | 35 ++++++++++++++++++------------ 9 files changed, 52 insertions(+), 57 deletions(-) diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index 1b1a7215..030c5f87 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -162,7 +162,6 @@ --physics-suites adiabatic - --physics-suites tj2016 --analytic_ic --physics-suites kessler --analytic_ic diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 833b8e4b..a27f45b1 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -544,6 +544,7 @@ end subroutine dry_air_composition_update !=========================================================================== subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + use string_utils, only: stringify real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array integer, intent(in) :: ncol ! number of columns @@ -559,24 +560,21 @@ subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) ! FV: moist pressure vertical coordinate does not need update. else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then ! SE - - ! **TEMP** TODO hplin 9/17/24: - ! for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected) - ! (whereas CAM only allocates 1-indexed) I subset it here. This should be verified during code - ! review. + ! Note: species index subset to 1: because SIMA currently uses index 0. See #334. 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 #334. call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), & cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) ! internal energy coefficient for MPAS - ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) + ! (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 not supported') + call endrun(subname//': dycore energy formula (value = '//stringify((/energy_formula/))//') not supported') end if end subroutine water_composition_update @@ -696,14 +694,14 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpd ! Dummy arguments ! tracer: Tracer array ! - ! factor not present then tracer must be dry mixing ratio - ! if factor present tracer*factor must be dry mixing ratio + ! 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(:,:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) - ! dp: if provided then tracer is mass not mixing ratio + ! if provided then tracer is not a 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 diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index b6f0f568..89fc6d44 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -221,14 +221,14 @@ end subroutine cam_thermo_init subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) use air_composition, only: dry_air_composition_update use air_composition, only: update_zvirv - use string_utils, only: to_str + use string_utils, only: stringify real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) real(kind_phys), intent(in) :: T(:,:) ! temperature integer, intent(in) :: ncol ! number of columns logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables ! false: do not calculate composition-dependent thermo variables - real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr moist convert + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr wet or moist convert ! Local vars real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) @@ -240,7 +240,7 @@ subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_d if (present(to_dry_factor)) then if (SIZE(to_dry_factor, 1) /= ncol) then - call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) + call endrun(subname//'DIM 1 of to_dry_factor is'//stringify((/SIZE(to_dry_factor,1)/))//'but should be'//stringify((/ncol/))) end if end if @@ -263,7 +263,7 @@ end subroutine cam_thermo_dry_air_update ! !*************************************************************************** ! - subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor) + subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, to_dry_factor) use air_composition, only: water_composition_update !----------------------------------------------------------------------- ! Update the physics "constants" that vary @@ -271,9 +271,21 @@ subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor) real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: pver ! number of vertical levels integer, intent(in) :: energy_formula 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, to_dry_factor=to_dry_factor) end subroutine cam_thermo_water_update diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90 index c3e2c15b..df202c9d 100644 --- a/src/data/cam_thermo_formula.F90 +++ b/src/data/cam_thermo_formula.F90 @@ -1,5 +1,7 @@ module cam_thermo_formula + use runtime_obj, only: unset_int + implicit none private save @@ -17,9 +19,9 @@ module cam_thermo_formula !! \htmlinclude cam_thermo_formula.html ! energy_formula_dycore: energy formula used for dynamical core ! written by the dynamical core - integer, public :: energy_formula_dycore + integer, public :: energy_formula_dycore = unset_int ! energy_formula_physics: energy formula used for physics - integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV + integer, public :: energy_formula_physics = unset_int ! Public subroutines public :: cam_thermo_formula_init diff --git a/src/data/registry.xml b/src/data/registry.xml index d449bf09..91658e20 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -264,7 +264,7 @@ - flag indicating if dynamical core energy not consistent with CAM physics and to perform adjustment of temperature and temperature tendency + flag indicating if dynamical core energy is not consistent with CAM physics and to perform adjustment of temperature and temperature tendency .false. .true. @@ -272,7 +272,7 @@ - flag indicating if is first timestep of initial run + flag indicating if it is the first timestep of an initial run . @@ -620,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 @@ -638,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 index ca6dd6d4..7959ebf0 100644 --- a/src/utils/gmean_mod.F90 +++ b/src/utils/gmean_mod.F90 @@ -115,16 +115,16 @@ subroutine gmean_scl (arr, gmean) ! ! Arguments ! - real(r8), intent(in) :: arr(pcols) - ! Input array, chunked - real(r8), intent(out):: gmean ! global means + 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 + 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) @@ -192,7 +192,7 @@ end subroutine gmean_float_norepro ! !======================================================================== ! - subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) + subroutine gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds) !----------------------------------------------------------------------- ! ! Purpose: @@ -206,6 +206,7 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) 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 ! @@ -218,14 +219,20 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) ! ! Local workspace ! - integer :: icol, ifld ! column, field indices + 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 - ! - !----------------------------------------------------------------------- - ! - allocate(xfld(pcols, nflds)) + + 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) @@ -249,7 +256,7 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds) deallocate(xfld) ! final normalization - arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi) + arr_gmean(:) = arr_gmean(:) / pi4 end subroutine gmean_fixed_repro From fc1789b50ed1a605e176ff2b46c3cb1911f15847 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 11 Dec 2024 20:11:33 -0500 Subject: [PATCH 30/41] Remove mention of check_energy_timestep_init in dp_coupling since it is initialized by CCPP --- src/dynamics/se/dp_coupling.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 9198b84c..450a691c 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -592,7 +592,6 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) use geopotential_temp, only: geopotential_temp_run use static_energy, only: update_dry_static_energy_run use qneg, only: qneg_run -! use check_energy_chng, only: check_energy_chng_timestep_init use hycoef, only: hyai, ps0 use shr_vmath_mod, only: shr_vmath_log use shr_kind_mod, only: shr_kind_cx @@ -854,12 +853,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_chng_timestep_init(...) -#endif - end subroutine derived_phys_dry !========================================================================================= From 58a151ad026c03090c80b26be4979c252c65125c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 12:57:34 -0500 Subject: [PATCH 31/41] Fix build --- src/data/cam_thermo.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 89fc6d44..8d6715f3 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -265,6 +265,7 @@ end subroutine cam_thermo_dry_air_update ! subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, to_dry_factor) use air_composition, only: water_composition_update + use string_utils, only: stringify !----------------------------------------------------------------------- ! Update the physics "constants" that vary !------------------------------------------------------------------------- From 63dff8fd831df75d3cd2c6ef1d31ab19ee59b77c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 15:36:54 -0500 Subject: [PATCH 32/41] Add pver dim check in cam_thermo_dry_air_update --- src/data/cam_thermo.F90 | 10 +++++++--- src/dynamics/se/dp_coupling.F90 | 1 + 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 8d6715f3..b638747b 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -218,17 +218,18 @@ end subroutine cam_thermo_init ! !*************************************************************************** ! - subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) + 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 (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(:,:) ! if mmr wet or moist convert + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! conversion factor to dry if mmr is wet or moist ! Local vars real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) @@ -242,6 +243,9 @@ subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_d 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 @@ -270,7 +274,7 @@ subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, to_dry_facto ! Update the physics "constants" that vary !------------------------------------------------------------------------- - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + 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 diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 450a691c..acdc2371 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -792,6 +792,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) mmr = const_data_ptr, & ! dry MMR T = phys_state%t, & ncol = pcols, & + pver = pver, & update_thermo_variables = cam_runtime_opts%update_thermodynamic_variables() & ) else From 55ea5b333df123d2e305580ef0f7cfc9f01e2dff Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 15:50:24 -0500 Subject: [PATCH 33/41] Pass cp_or_cv_dycore explicitly to cam_thermo_water_update in preparation for CCPP scheme --- src/data/air_composition.F90 | 16 +++++++--------- src/data/cam_thermo.F90 | 15 ++++++++------- src/dynamics/se/dp_coupling.F90 | 10 ++++++---- 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index a27f45b1..9140794c 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -7,7 +7,6 @@ module air_composition use runtime_obj, only: unset_real, unset_int use phys_vars_init_check, only: std_name_len use physics_types, only: cpairv, rairv, cappav, mbarv, zvirv - use physics_types, only: cp_or_cv_dycore implicit none private @@ -540,22 +539,22 @@ 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, to_dry_factor) + 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), optional, intent(in) :: to_dry_factor(:,:) + 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 - ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module - if (energy_formula == ENERGY_FORMULA_DYCORE_FV) then ! FV: moist pressure vertical coordinate does not need update. else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then @@ -576,7 +575,6 @@ subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor) else call endrun(subname//': dycore energy formula (value = '//stringify((/energy_formula/))//') not supported') end if - end subroutine water_composition_update !=========================================================================== diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index b638747b..aeab3d28 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -267,18 +267,19 @@ end subroutine cam_thermo_dry_air_update ! !*************************************************************************** ! - subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, to_dry_factor) + 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), optional, intent(in) :: to_dry_factor(:,:) + 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: ' @@ -291,7 +292,7 @@ subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, to_dry_facto end if end if - call water_composition_update(mmr, ncol, energy_formula, to_dry_factor=to_dry_factor) + call water_composition_update(mmr, ncol, energy_formula, cp_or_cv_dycore, to_dry_factor=to_dry_factor) end subroutine cam_thermo_water_update !=========================================================================== diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index acdc2371..11001da6 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -588,6 +588,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 @@ -804,10 +805,11 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! (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 & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + pver = pver, & + energy_formula = ENERGY_FORMULA_DYCORE_SE, & + cp_or_cv_dycore = cp_or_cv_dycore & ) !$omp parallel do num_threads(horz_num_threads) private (k, i) From a8b977b0c12a19652ef2085be8f84057493f805a Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 16:17:19 -0500 Subject: [PATCH 34/41] Fixes for MPAS dynamical core to support cam_thermo_water_update --- src/dynamics/mpas/dyn_coupling.F90 | 22 +++++++++++++++++----- src/dynamics/se/dp_coupling.F90 | 2 +- 2 files changed, 18 insertions(+), 6 deletions(-) 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/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 11001da6..4122a218 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -801,7 +801,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end if ! - ! update cp_or_cv_dycore in module air_composition. + ! update cp_or_cv_dycore in SIMA state. ! (note: at this point q is dry) ! call cam_thermo_water_update( & From 7448d737d9f6ec769788c505bdd9d304ec9a371c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 16:19:30 -0500 Subject: [PATCH 35/41] Fix FPHYStest test failure --- src/dynamics/none/dyn_grid.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index 465fb703..54c98f02 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -687,6 +687,7 @@ subroutine find_energy_formula(file, grid_is_latlon) 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 From 4bbb9200be29bc170bde0690aa6d77105c08dba5 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Dec 2024 18:09:27 -0500 Subject: [PATCH 36/41] Address review comments --- src/data/air_composition.F90 | 7 ++++--- src/data/cam_thermo.F90 | 8 ++++---- src/dynamics/none/dyn_grid.F90 | 5 +++++ src/dynamics/se/dp_coupling.F90 | 4 ---- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index 9140794c..51e7dd6b 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -559,13 +559,13 @@ subroutine water_composition_update(mmr, ncol, energy_formula, cp_or_cv_dycore, ! 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 #334. + ! 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 #334. + ! 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,:)) @@ -699,7 +699,8 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpd ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) - ! if provided then tracer is not a mass mixing ratio + ! 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 diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index aeab3d28..8330ef64 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -241,10 +241,10 @@ subroutine cam_thermo_dry_air_update(mmr, T, ncol, pver, update_thermo_variables 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/))) + 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/))) + 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 @@ -285,10 +285,10 @@ subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, cp_or_cv_dyc 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/))) + 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/))) + 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 diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90 index 54c98f02..ba1bf0ba 100644 --- a/src/dynamics/none/dyn_grid.F90 +++ b/src/dynamics/none/dyn_grid.F90 @@ -645,6 +645,11 @@ subroutine find_energy_formula(file, grid_is_latlon) ! 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 diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 4122a218..572f663f 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -784,10 +784,6 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! Compute molecular viscosity(kmvis) and conductivity(kmcnd). ! Update zvirv registry variable; calculated for WACCM-X. !----------------------------------------------------------------------------- - ! **TEMP** TODO CHECK hplin: CAM has this if-clause for dry_air_species_num > 0 - ! or otherwise uses zvirv = zvir. CAM-SIMA previously did not have this, and - ! instead has a switch for update_thermodynamic_variables. Check if we still want - ! this if-clause or change it to something else. if (dry_air_species_num > 0) then call cam_thermo_dry_air_update( & mmr = const_data_ptr, & ! dry MMR From e28e9b9f771301c361e0fe082c7ea2b341c7f3a9 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Dec 2024 15:47:29 -0500 Subject: [PATCH 37/41] Update src/dynamics/mpas/dyn_comp.F90 Co-authored-by: Kuan-Chih Wang --- src/dynamics/mpas/dyn_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 2cb9c482..b8d881f8 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -188,8 +188,8 @@ end subroutine dyn_readnl ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS - use physics_types, only: dycore_energy_consistency_adjust use phys_vars_init_check, only: mark_as_initialized + use physics_types, only: dycore_energy_consistency_adjust type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in From 2769a0d71ee102735e4543f066d65a938bff5f15 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Dec 2024 15:47:38 -0500 Subject: [PATCH 38/41] Update src/dynamics/mpas/dyn_comp.F90 Co-authored-by: Kuan-Chih Wang --- src/dynamics/mpas/dyn_comp.F90 | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index b8d881f8..d5e124f3 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -976,15 +976,13 @@ subroutine mark_variable_as_initialized() ! 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_at_end_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_energy_using_dycore_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_physics_energy_formula_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') - - + call mark_as_initialized('vertically_integrated_total_water_at_start_of_physics_timestep') end subroutine mark_variable_as_initialized !> Run MPAS dynamical core to integrate the dynamical states with time. From 835bfbea994ce6f47deeb25f294ee0a78ca1a987 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Dec 2024 15:48:20 -0500 Subject: [PATCH 39/41] Update src/dynamics/mpas/dyn_comp.F90 Co-authored-by: Kuan-Chih Wang --- src/dynamics/mpas/dyn_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index d5e124f3..32d7f551 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -209,7 +209,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Set dynamical core energy formula for use in cam_thermo. energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS - call mark_as_initialized("total_energy_formula_for_dycore") + 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. From 1df32517fc191791a7433ea21dbc47bd4e1f0eb9 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Dec 2024 15:48:29 -0500 Subject: [PATCH 40/41] Update src/dynamics/mpas/dyn_comp.F90 Co-authored-by: Kuan-Chih Wang --- src/dynamics/mpas/dyn_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 32d7f551..3911af4f 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -214,7 +214,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! 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") + call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment') allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) From dc50feacbe3acf83c79e27472746ae203658d263 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Dec 2024 21:27:04 -0500 Subject: [PATCH 41/41] Update list of existing test failures --- test/existing-test-failures.txt | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) 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