Implements cam_thermo_water_update and CCPPized check_energy

cime_config/config_component.xml

@@ -161,8 +161,7 @@
       <value compset="_CAM\d0%WX.*%SDYN">-nlev 145</value> -->
 
       <!-- Simple models -->
-<!--      <value compset="_CAM%ADIAB">-phys adiabatic</value>
-      <value compset="_CAM%DABIP04">-phys adiabatic</value> -->
+      <value compset="_CAM%ADIAB">--physics-suites 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/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, nstep
    use dyn_comp,                  only: dyn_import_t, dyn_export_t
 
    use perf_mod,                  only: t_barrierf, t_startf, t_stopf
@@ -149,9 +150,6 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
       character(len=cx)        :: errmsg
       !-----------------------------------------------------------------------
 
-      dtime_phys = 0.0_r8
-      call mark_as_initialized('timestep_for_physics')
-
       call init_pio_subsystem()
 
       ! Initializations using data passed from coupler.
@@ -167,12 +165,20 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
 
       call cam_ctrl_set_orbit(eccen, obliqr, lambm0, mvelpp)
 
-
       call timemgr_init(                                                      &
            dtime, calendar, start_ymd, start_tod, ref_ymd,                    &
            ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod,                   &
            perpetual_run, perpetual_ymd, initial_run_in)
 
+      dtime_phys = 0.0_r8
+      call mark_as_initialized('timestep_for_physics')
+
+      is_first_timestep = .true.
+      call mark_as_initialized('is_first_timestep')
+
+      nstep = get_nstep()
+      call mark_as_initialized('current_timestep_number')
+
       ! Get current fractional calendar day. Needs to be updated at every timestep.
       calday = get_curr_calday()
       call mark_as_initialized('fractional_calendar_days_on_end_of_current_timestep')
@@ -268,6 +274,10 @@ 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()
+      nstep = get_nstep()
+
       !----------------------------------------------------------
       ! First phase of dynamics (at least couple from dynamics to physics)
       ! Return time-step for physics from dynamics.
@@ -514,10 +524,6 @@ subroutine cam_final(cam_out, cam_in)
       type(cam_out_t), pointer :: cam_out ! Output from CAM to surface
       type(cam_in_t),  pointer :: cam_in  ! Input from merged surface to CAM
 
-      !
-      ! Local variable
-      !
-      integer                  :: nstep   ! Current timestep number.
       !-----------------------------------------------------------------------
 
       call phys_final()
@@ -540,7 +546,6 @@ subroutine cam_final(cam_out, cam_in)
       call shr_sys_flush( iulog )   ! Flush all output to the CAM log file
 
       if (masterproc) then
-         nstep = get_nstep()
          write(iulog,9300) nstep-1,nstep
 9300     format (//'Number of completed timesteps:',i6,/,'Time step ',i6,     &
               ' partially done to provide convectively adjusted and ',        &

src/data/air_composition.F90

@@ -1,4 +1,5 @@
-! 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
@@ -12,7 +13,9 @@ module air_composition
    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 +228,7 @@ subroutine air_composition_init()
       !
       !************************************************************************
       !
-      ! add prognostic components of dry air
+      ! add prognostic components of air
       !
       !************************************************************************
       !
@@ -309,6 +312,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 +514,68 @@ 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
+   ! (enthalpy for pressure-based dynamical cores and internal energy for z-based dynamical cores)
+   !---------------------------------------------------------------------------
+   !===========================================================================
+   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),           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
+      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
+         ! 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; 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 (value = '//stringify((/energy_formula/))//') not supported')
+      end if
+   end subroutine water_composition_update
 
    !===========================================================================
    !***************************************************************************
@@ -639,27 +685,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
+      !
+      ! 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(:,:,:)
-      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(:,:)
+      ! 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
       !    (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 +727,37 @@ 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
+         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 +768,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 +911,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 +925,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 +944,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 +1011,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