Skip to content

Commit

Permalink
Implements cam_thermo_water_update and CCPPized check_energy (#316)
Browse files Browse the repository at this point in the history 
Originator(s): jimmielin

Description (include the issue title, and the keyword ['closes',
'fixes', 'resolves'] followed by the issue number):

All changes are bit-for-bit, except those noted:

Implements `cam_thermo_water_update`:
- updates `cp_or_cv_dycore` (`specific_heat_for_air_used_in_dycore`) in
`air_composition.F90` for SE and MPAS dynamical cores
- read `cp_or_cv_dycore` from CAM snapshot (refer to companion CAM PR)
- energy formula (that has to be matching dycore) is recognized and set
by null-dycore in `dyn_grid.F90` by looking at global attributes of the
initial file; see `find_energy_formula`
- added `is_first_timestep` logical state flag

Ports `cam_thermo` and related updates in `air_composition` and
`dp_coupling` from CAM 6.3.109
(https://github.com/ESCOMP/CAM/pull/761/files)
- update to hydrostatic energy calculation
- changes `get_cp`, `get_R` in `air_composition.F90` to use moist mixing
ratios
- **answer-changing:** update to moist-to-dry (for physics) conversion
in `dp_coupling::derived_phys_dry` to account for all water tracers
instead of just Q
- **answer-changing:** update to not-really-"exner" calculation to use
composition-dependent `cappav` instead of `cappa`

Changes `vcoord` in `dyn_tests_utils` (old CAM) to `energy_formula` now
in `cam_thermo_formula` (separated out into a different file to avoid
dependency issues)
- `vc_moist_pressure` is now `ENERGY_FORMULA_DYCORE_FV`;
`vc_dry_pressure` is `_SE`; `vc_height` is `_MPAS`
- these are just integer flags (0,1,2) and values are kept consistent
with old CAM and their use in dynamics tests

Ports global mean utility module (`gmean_mod.F90`), de-chunkized from
CAM:
- Implements `get_wght` in `physics_grid` for weighted sum calculation

Imports `check_energy_chng` and `check_energy_fix` from
`atmospheric_physics`

Describe any changes made to build system: N/A

Describe any changes made to the namelist: contained within ncar-physics

List any changes to the defaults for the input datasets (e.g. boundary
datasets):
- Added `cp_or_cv_dycore` in CAM snapshots

List all files eliminated and why: N/A

List all files added and what they do:

```
- energy_formula
A       src/data/cam_thermo_formula.F90
A       src/data/cam_thermo_formula.meta

- gmean
A       src/utils/gmean_mod.F90
```

List all existing files that have been modified, and describe the
changes:
(Helpful git command: `git diff --name-status
development...<your_branch_name>`)

```
- ncar-physics update
M       .gitmodules
M       src/physics/ncar_ccpp

- `is_first_timestep`
M       src/control/cam_comp.F90

- cam_thermo_water_update for cp_or_cv_dycore (include in registry; read from ic)
M       src/data/air_composition.F90
M       src/data/cam_thermo.F90
M       src/data/registry.xml
M       tools/stdnames_to_inputnames_dictionary.xml
M       src/dynamics/se/dp_coupling.F90
M       src/dynamics/se/dycore/prim_advance_mod.F90
M       src/dynamics/se/dyn_comp.F90
M       src/dynamics/utils/dyn_thermo.F90

- energy_formula
M       src/physics/utils/phys_comp.F90
M       src/dynamics/mpas/dyn_comp.F90
M       src/dynamics/none/dyn_comp.F90
M       src/dynamics/none/dyn_grid.F90

- gmean
M       src/physics/utils/physics_grid.F90
```

Note: bit-for-bit in check_energy with CAM is tricky to validate without
dycore updates to SE; may need to merge #301 first

---------

Co-authored-by: Kuan-Chih Wang <[email protected]>
  • Loading branch information
@jimmielin @kuanchihwang
jimmielin and kuanchihwang authored Dec 16, 2024
1 parent edf84cd commit 37fdbfb
Show file tree
Hide file tree
Showing 20 changed files with 955 additions and 233 deletions.
3 changes: 1 addition & 2 deletions cime_config/config_component.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
23 changes: 14 additions & 9 deletions src/control/cam_comp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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')
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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()
Expand All @@ -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 ', &
Expand Down
144 changes: 111 additions & 33 deletions src/data/air_composition.F90
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand Down Expand Up @@ -225,7 +228,7 @@ subroutine air_composition_init()
!
!************************************************************************
!
! add prognostic components of dry air
! add prognostic components of air
!
!************************************************************************
!
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 GitHub issue #334 in ESCOMP/CAM-SIMA.
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 GitHub issue #334 in ESCOMP/CAM-SIMA.
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

!===========================================================================
!***************************************************************************
Expand Down Expand Up @@ -639,27 +685,34 @@ 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(:,:)
! factor: to convert tracer to dry mixing ratio
! if provided, then tracer is not a dry 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: '

Expand All @@ -675,28 +728,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
Expand All @@ -707,34 +769,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
Expand Down Expand Up @@ -843,9 +912,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
Expand All @@ -856,6 +926,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
Expand All @@ -874,12 +945,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
Expand Down Expand Up @@ -934,7 +1012,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
Expand Down
Loading

0 comments on commit 37fdbfb

Please sign in to comment.