From 86489702461192ccb13ae7cff9d13f0db004812b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 12 Sep 2024 16:24:55 -0400 Subject: [PATCH] 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 = ()