Implements cam_thermo_water_update and CCPPized check_energy

cime_config/config_component.xml

@@ -162,7 +162,6 @@
 
       <!-- Simple models -->
       <value compset="_CAM%ADIAB">--physics-suites adiabatic</value>
-<!--      <value compset="_CAM%DABIP04">-phys adiabatic</value> -->
       <value compset="_CAM%TJ16">--physics-suites tj2016 --analytic_ic</value>
 <!--      <value compset="_CAM%KESSLER">-phys kessler -chem terminator -analytic_ic</value> -->
       <value compset="_CAM%KESSLER">--physics-suites kessler --analytic_ic</value>

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,45 +539,42 @@ 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
          ! 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 +692,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

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: 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) :: 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 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))
@@ -240,7 +241,10 @@ 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
+        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
 
@@ -263,18 +267,32 @@ 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, 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
-      integer,                   intent(in) :: ncol       ! number of columns
-      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: '
+
+      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)
+      call water_composition_update(mmr, ncol, energy_formula, cp_or_cv_dycore, to_dry_factor=to_dry_factor)
     end subroutine cam_thermo_water_update
 
    !===========================================================================

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

src/data/registry.xml

@@ -264,15 +264,15 @@
     <variable local_name="dycore_energy_consistency_adjust"
               standard_name="flag_for_dycore_energy_consistency_adjustment"
               units="flag" type="logical">
-      <long_name>flag indicating if dynamical core energy not consistent with CAM physics and to perform adjustment of temperature and temperature tendency</long_name>
+      <long_name>flag indicating if dynamical core energy is not consistent with CAM physics and to perform adjustment of temperature and temperature tendency</long_name>
       <initial_value>.false.</initial_value>
       <initial_value dycore="SE,MPAS">.true.</initial_value>
     </variable>
     <!-- Timestep properties -->
     <variable local_name="is_first_timestep"
               standard_name="is_first_timestep"
               units="flag" type="logical">
-      <long_name>flag indicating if is first timestep of initial run</long_name>
+      <long_name>flag indicating if it is the first timestep of an initial run</long_name>
     </variable>
     <!-- State tendencies -->
     <variable local_name="dTdt_total"

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

src/dynamics/none/dyn_grid.F90

@@ -62,7 +62,6 @@ 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
@@ -688,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

src/dynamics/se/dp_coupling.F90

@@ -588,11 +588,11 @@ 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
    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
@@ -793,20 +793,23 @@ 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
       zvirv(:,:) = zvir
    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( &
-        mmr            = const_data_ptr,          & ! dry MMR
-        ncol           = pcols,                   &
-        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)
@@ -853,12 +856,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
 
 !=========================================================================================