Initial implementation of cp_to_cv_dycore/cam_thermo_water_update fro…

…m stash

Needs work to split into cam_thermo_water_update. This is WIP code

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
 
    !===========================================================================

src/data/cam_thermo.F90

@@ -1554,28 +1554,31 @@ 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(:,:)
       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(:)
+      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,78 +1658,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 ',     &
+         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
          write(iulog, *) subname, ' vertical coordinate not supported: ', vcoord
          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
 
    !===========================================================================

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

src/data/registry.xml

@@ -415,6 +415,13 @@
       <dimensions>horizontal_dimension vertical_layer_dimension</dimensions>
       <initial_value>zvir</initial_value>
     </variable>
+    <variable local_name="cp_or_cv_dycore"
+              standard_name="enthalpy_or_internal_energy_scaling_factor_for_energy_consistency?"
+              units="1" type="real" kind="kind_phys"
+              allocatable="allocatable">
+      <long_name>Enthalpy or internal energy scaling factor for energy consistency</long_name>
+      <dimensions>horizontal_dimension vertical_layer_dimension</dimensions>
+    </variable>
     <!-- Constituent Variables -->
     <!-- These are only used to set possible IC file input names, as the constituents object handles allocation. -->
     <variable local_name="q"

src/dynamics/se/dyn_comp.F90

@@ -581,6 +581,7 @@ 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
@@ -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 vertical coordinate
+   vc_dycore = vc_dry_pressure
 
    ! Now allocate and set condenstate vars
    allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers

src/dynamics/tests/dyn_tests_utils.F90

@@ -20,4 +20,10 @@ 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

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 = ()