Update vcoord to energy_formula

src/data/air_composition.F90

@@ -541,22 +541,22 @@ 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'
 
       ! 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

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,18 +263,18 @@ 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
       !-------------------------------------------------------------------------
 
       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

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

src/dynamics/tests/dyn_tests_utils.meta → 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

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

src/data/registry.xml

@@ -13,10 +13,10 @@
   <metadata_file>$SRCROOT/src/physics/utils/tropopause_climo_read.meta</metadata_file>
   <metadata_file>$SRCROOT/src/data/air_composition.meta</metadata_file>
   <metadata_file>$SRCROOT/src/data/cam_thermo.meta</metadata_file>
+  <metadata_file>$SRCROOT/src/data/cam_thermo_formula.meta</metadata_file>
   <metadata_file>$SRCROOT/src/data/ref_pres.meta</metadata_file>
   <metadata_file>$SRCROOT/src/dynamics/utils/vert_coord.meta</metadata_file>
   <metadata_file>$SRCROOT/src/dynamics/utils/hycoef.meta</metadata_file>
-  <metadata_file>$SRCROOT/src/dynamics/tests/dyn_tests_utils.meta</metadata_file>
   <!-- Variables registered to physics_types -->
   <file name="physics_types" type="module">
     <use module="ccpp_kinds" reference="kind_phys"/>

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__)
 

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

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 <file> 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