Implements CCPPized check_energy in SIMA.

Squashed commit of the following:

commit 8b8a8e0
    Fix merge conflicts

commit bef25fb
    Update atmos_phys external

commit 240ef54
    Set wv_idx in air_composition

    total_hours_in_debugging_one_line = 6

commit 38fdbb2
    Read cp_or_cv_dycore and identify dycore information from CAM snapshot

commit 8fca3a4
    Update vcoord to energy_formula

commit 0d16be9
    New const_get_index logic without cam_ccpp_cap dependency

commit e13e7ea
    Fix build (part 1); update standard names and atmos_phys external

commit c015e28
    Update ncar-physics external

commit 2372f7f
    Fix for b4b to CAM-SIMA (w/ History); add notes on modifications and/or divergence from CAM and/or -SIMA

commit 4ce8427
    Add gmean_mod to src/utils including support infra in physics_grid.F90

commit 7493a1d
    Fix build issues; update factor in get_cp call to 1/dp

commit 041cdfb
    Add is_first_timestep registry

commit 333ad4e
    Initial port of updates to energy budget (port of CAM #761) into CAM-SIMA

commit 8648970
    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

@@ -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"]

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.

src/data/air_composition.F90

@@ -1,18 +1,22 @@
-! 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
    use cam_abortutils,       only: endrun, check_allocate
    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,34 +776,41 @@ 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
 
       ! 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
       integer  :: idx_local(thermodynamic_active_species_num)
       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