diff --git a/.gitmodules b/.gitmodules
index fe5de260..5322028a 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -20,7 +20,7 @@
[submodule "ncar-physics"]
path = src/physics/ncar_ccpp
url = https://github.com/ESCOMP/atmospheric_physics
- fxtag = 0ecfcc155ac0387ef9db3304611c6f3ef055ac1d
+ fxtag = e7a599f4bb1533f7cdcd8723b1f864e11578e96c
fxrequired = AlwaysRequired
fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
[submodule "ccs_config"]
diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml
index 5805131c..030c5f87 100644
--- a/cime_config/config_component.xml
+++ b/cime_config/config_component.xml
@@ -161,8 +161,7 @@
-nlev 145 -->
-
+ --physics-suites adiabatic
--physics-suites tj2016 --analytic_ic
--physics-suites kessler --analytic_ic
diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90
index 3daf0168..e3a178f8 100644
--- a/src/control/cam_comp.F90
+++ b/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 ', &
diff --git a/src/control/cam_control_mod.F90 b/src/control/cam_control_mod.F90
index 4dc86bae..d135eea8 100644
--- a/src/control/cam_control_mod.F90
+++ b/src/control/cam_control_mod.F90
@@ -28,10 +28,12 @@ module cam_control_mod
logical, protected :: branch_run ! branch from a previous run; requires a restart file
logical, protected :: post_assim ! We are resuming after a pause
- logical, protected :: aqua_planet ! Flag to run model in "aqua planet" mode
logical, protected :: brnch_retain_casename ! true => branch run may use same caseid as
! the run being branched from
+ !> \section arg_table_cam_control_mod Argument Table
+ !! \htmlinclude arg_table_cam_control_mod.html
+ logical, protected :: aqua_planet ! Flag to run model in "aqua planet" mode
real(r8), protected :: eccen ! Earth's eccentricity factor (unitless) (typically 0 to 0.1)
real(r8), protected :: obliqr ! Earth's obliquity in radians
real(r8), protected :: lambm0 ! Mean longitude of perihelion at the
diff --git a/src/control/cam_control_mod.meta b/src/control/cam_control_mod.meta
new file mode 100644
index 00000000..d7b38dbd
--- /dev/null
+++ b/src/control/cam_control_mod.meta
@@ -0,0 +1,33 @@
+[ccpp-table-properties]
+ name = cam_control_mod
+ type = module
+
+[ccpp-arg-table]
+ name = cam_control_mod
+ type = module
+[ aqua_planet ]
+ standard_name = is_aqua_planet
+ units = flag
+ type = logical
+ dimensions = ()
+[ eccen ]
+ standard_name = planet_orbital_eccentricity_factor
+ units = 1
+ type = real | kind = r8
+ dimensions = ()
+[ obliqr ]
+ standard_name = planet_obliquity
+ long_name = planet's axial tilt (obliquity)
+ units = rad
+ type = real | kind = r8
+ dimensions = ()
+[ lambm0 ]
+ standard_name = mean_longitude_of_perihelion_at_vernal_equinox
+ units = rad
+ type = real | kind = r8
+ dimensions = ()
+[ mvelpp ]
+ standard_name = moving_vernal_equinox_longitude_of_perihelion_plus_pi
+ units = rad
+ type = real | kind = r8
+ dimensions = ()
diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90
index e84fc837..51e7dd6b 100644
--- a/src/data/air_composition.F90
+++ b/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 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
!===========================================================================
!***************************************************************************
@@ -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: '
@@ -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
@@ -707,7 +769,7 @@ 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
@@ -715,14 +777,15 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore)
! 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
@@ -730,11 +793,17 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore)
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 +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
@@ -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
@@ -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
@@ -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
diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90
index 59dd1c83..8330ef64 100644
--- a/src/data/cam_thermo.F90
+++ b/src/data/cam_thermo.F90
@@ -33,8 +33,10 @@ module cam_thermo
! cam_thermo_init: Initialize constituent dependent properties
public :: cam_thermo_init
- ! cam_thermo_update: Update constituent dependent properties
- public :: cam_thermo_update
+ ! cam_thermo_dry_air_update: Update dry air composition dependent properties
+ public :: cam_thermo_dry_air_update
+ ! cam_thermo_water_update: Update water dependent properties
+ public :: cam_thermo_water_update
! get_enthalpy: enthalpy quantity = dp*cp*T
public :: get_enthalpy
! get_virtual_temp: virtual temperature
@@ -77,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 ---------------------
@@ -85,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 --------------------
!
!
@@ -208,51 +211,89 @@ end subroutine cam_thermo_init
!===========================================================================
+ !
!***************************************************************************
!
- ! cam_thermo_update: update species dependent constants for physics
+ ! cam_thermo_dry_air_update: update dry air species dependent constants for physics
!
!***************************************************************************
!
- subroutine cam_thermo_update(mmr, T, ncol, update_thermo_variables, to_moist_factor)
- use air_composition, only: air_composition_update, update_zvirv
- use string_utils, only: to_str
- !-----------------------------------------------------------------------
- ! Update the physics "constants" that vary
- !-------------------------------------------------------------------------
-
- !------------------------------Arguments----------------------------------
+ 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: stringify
- real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array
- real(kind_phys), intent(in) :: T(:,:) ! temperature
- 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), 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(:,:) ! conversion factor to dry if mmr is wet or moist
- real(kind_phys), optional, intent(in) :: to_moist_factor(:,:)
- !
- !---------------------------Local storage-------------------------------
- real(kind_phys):: sponge_factor(SIZE(mmr, 2))
- character(len=*), parameter :: subname = 'cam_thermo_update: '
+ ! Local vars
+ real(kind_phys) :: sponge_factor(SIZE(mmr, 2))
+ character(len=*), parameter :: subname = 'cam_thermo_dry_air_update: '
if (.not. update_thermo_variables) then
return
end if
- if (present(to_moist_factor)) then
- if (SIZE(to_moist_factor, 1) /= ncol) then
- call endrun(subname//'DIM 1 of to_moist_factor is'//to_str(SIZE(to_moist_factor,1))//'but should be'//to_str(ncol))
- end if
+ 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
+
sponge_factor = 1.0_kind_phys
- call air_composition_update(mmr, ncol, to_moist_factor=to_moist_factor)
+ call dry_air_composition_update(mmr, ncol, to_dry_factor=to_dry_factor)
call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:), &
- kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_moist_factor, &
+ kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_dry_factor, &
active_species_idx_dycore=thermodynamic_active_species_idx)
+
+ ! Calculate zvirv for WACCM-X.
call update_zvirv()
- end subroutine cam_thermo_update
+ end subroutine cam_thermo_dry_air_update
+
+ !
+ !***************************************************************************
+ !
+ ! cam_thermo_water_update: update water species dependent constants for physics
+ !
+ !***************************************************************************
+ !
+ 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 (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, cp_or_cv_dycore, to_dry_factor=to_dry_factor)
+ end subroutine cam_thermo_water_update
!===========================================================================
@@ -1554,28 +1595,32 @@ end subroutine cam_thermo_calc_kappav_2hd
!
!***************************************************************************
!
- subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
- vcoord, ps, phis, z_mid, dycore_idx, qidx, te, se, ke, &
- wv, H2O, liq, ice)
+ 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 physconst, only: gravit, latvap, latice
+ use air_composition, only: dry_air_species_num
+ use physconst, only: rga, latvap, latice
! Dummy arguments
! tracer: tracer mixing ratio
+ !
+ ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry
real(kind_phys), intent(in) :: tracer(:,:,:)
+ logical, intent(in) :: moist_mixing_ratio
! pdel: pressure level thickness
- real(kind_phys), intent(in) :: pdel(:,:)
- ! cp_or_cv: dry air heat capacity under constant pressure or
- ! constant volume (depends on vcoord)
+ real(kind_phys), intent(in) :: pdel_in(:,:)
+ ! cp_or_cv: air heat capacity under constant pressure or
+ ! constant volume (depends on energy formula)
real(kind_phys), intent(in) :: cp_or_cv(:,:)
real(kind_phys), intent(in) :: U(:,:)
real(kind_phys), intent(in) :: V(:,:)
real(kind_phys), intent(in) :: T(:,:)
- integer, intent(in) :: vcoord ! vertical coordinate
- real(kind_phys), intent(in), optional :: ps(:)
+ 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(:,:)
! dycore_idx: use dycore index for thermodynamic active species
@@ -1588,8 +1633,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
real(kind_phys), intent(out), optional :: te (:)
! KE: vertically integrated kinetic energy
real(kind_phys), intent(out), optional :: ke (:)
- ! SE: vertically integrated internal+geopotential energy
+ ! SE: vertically integrated enthalpy (pressure coordinate)
+ ! or internal energy (z coordinate)
real(kind_phys), intent(out), optional :: se (:)
+ ! PO: vertically integrated PHIS term (pressure coordinate)
+ ! or potential energy (z coordinate)
+ real(kind_phys), intent(out), optional :: po (:)
! WV: vertically integrated water vapor
real(kind_phys), intent(out), optional :: wv (:)
! liq: vertically integrated liquid
@@ -1599,10 +1648,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
! Local variables
real(kind_phys) :: ke_vint(SIZE(tracer, 1)) ! Vertical integral of KE
- real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of SE
+ real(kind_phys) :: se_vint(SIZE(tracer, 1)) ! Vertical integral of enthalpy or internal energy
+ real(kind_phys) :: po_vint(SIZE(tracer, 1)) ! Vertical integral of PHIS or potential energy
real(kind_phys) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv
real(kind_phys) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq
real(kind_phys) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice
+ real(kind_phys) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) ! moist pressure level thickness
real(kind_phys) :: latsub ! latent heat of sublimation
integer :: ierr
@@ -1633,12 +1684,12 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:)
species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:)
else
- species_idx(:) = thermodynamic_active_species_idx(:)
+ species_idx(:) = thermodynamic_active_species_idx(1:)
species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
end if
else
- species_idx(:) = thermodynamic_active_species_idx(:)
+ species_idx(:) = thermodynamic_active_species_idx(1:)
species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
end if
@@ -1649,78 +1700,96 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
wvidx = wv_idx
end if
+ if (moist_mixing_ratio) then
+ pdel = pdel_in
+ else
+ pdel = pdel_in
+ do qdx = dry_air_species_num+1, thermodynamic_active_species_num
+ pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx))
+ end do
+ end if
+
+ ke_vint = 0._kind_phys
+ se_vint = 0._kind_phys
select case (vcoord)
- case(vc_moist_pressure, vc_dry_pressure)
- if ((.not. present(ps)) .or. (.not. present(phis))) then
- write(iulog, *) subname, ' ps and phis must be present for ', &
- 'moist/dry pressure vertical coordinate'
- call endrun(subname//': ps and phis must be present for '// &
- 'moist/dry pressure vertical coordinate')
+ 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 ', &
+ 'FV/SE energy formula'
+ call endrun(subname//': ptop and phis must be present for '// &
+ 'FV/SE energy formula')
end if
- ke_vint = 0._kind_phys
- se_vint = 0._kind_phys
- wv_vint = 0._kind_phys
+ po_vint = ptop
do kdx = 1, SIZE(tracer, 2)
do idx = 1, SIZE(tracer, 1)
ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * &
- 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit)
+ 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2)) * rga
se_vint(idx) = se_vint(idx) + (T(idx, kdx) * &
- cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit)
- wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * &
- pdel(idx, kdx) / gravit)
+ cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
+ po_vint(idx) = po_vint(idx)+pdel(idx, kdx)
+
end do
end do
do idx = 1, SIZE(tracer, 1)
- se_vint(idx) = se_vint(idx) + (phis(idx) * ps(idx) / gravit)
+ po_vint(idx) = (phis(idx) * po_vint(idx) * rga)
end do
- case(vc_height)
- if (.not. present(z_mid)) then
- write(iulog, *) subname, &
- ' z_mid must be present for height vertical coordinate'
- call endrun(subname//': z_mid must be present for height '// &
- 'vertical coordinate')
+ case(ENERGY_FORMULA_DYCORE_MPAS)
+ if (.not. present(phis)) then
+ write(iulog, *) subname, ' phis must be present for ', &
+ 'MPAS energy formula'
+ call endrun(subname//': phis must be present for '// &
+ 'MPAS energy formula')
end if
- ke_vint = 0._kind_phys
- se_vint = 0._kind_phys
- wv_vint = 0._kind_phys
+ po_vint = 0._kind_phys
do kdx = 1, SIZE(tracer, 2)
do idx = 1, SIZE(tracer, 1)
ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) * &
- 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) / gravit)
+ 0.5_kind_phys * (U(idx, kdx)**2 + V(idx, kdx)**2) * rga)
se_vint(idx) = se_vint(idx) + (T(idx, kdx) * &
- cp_or_cv(idx, kdx) * pdel(idx, kdx) / gravit)
+ cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
! z_mid is height above ground
- se_vint(idx) = se_vint(idx) + (z_mid(idx, kdx) + &
- phis(idx) / gravit) * pdel(idx, kdx)
- wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * &
- pdel(idx, kdx) / gravit)
+ po_vint(idx) = po_vint(idx) + (z_mid(idx, kdx) + &
+ phis(idx) * rga) * pdel(idx, kdx)
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 + ke_vint
+ te = se_vint + po_vint + ke_vint
end if
if (present(se)) then
se = se_vint
end if
+ if (present(po)) then
+ po = po_vint
+ end if
if (present(ke)) then
ke = ke_vint
end if
- if (present(wv)) then
- wv = wv_vint
- end if
!
! vertical integral of total liquid water
!
+ if (.not.moist_mixing_ratio) then
+ pdel = pdel_in! set pseudo density to dry
+ end if
+
+ wv_vint = 0._kind_phys
+ do kdx = 1, SIZE(tracer, 2)
+ do idx = 1, SIZE(tracer, 1)
+ wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) * &
+ pdel(idx, kdx) * rga)
+ end do
+ end do
+ if (present(wv)) wv = wv_vint
+
liq_vint = 0._kind_phys
do qdx = 1, thermodynamic_active_species_liq_num
do kdx = 1, SIZE(tracer, 2)
do idx = 1, SIZE(tracer, 1)
- liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * &
- tracer(idx, kdx, species_liq_idx(qdx)) / gravit)
+ liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) * &
+ tracer(idx, kdx, species_liq_idx(qdx)) * rga)
end do
end do
end do
@@ -1734,7 +1803,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
do kdx = 1, SIZE(tracer, 2)
do idx = 1, SIZE(tracer, 1)
ice_vint(idx) = ice_vint(idx) + (pdel(idx, kdx) * &
- tracer(idx, kdx, species_ice_idx(qdx)) / gravit)
+ tracer(idx, kdx, species_ice_idx(qdx)) * rga)
end do
end do
end do
@@ -1762,7 +1831,6 @@ subroutine get_hydrostatic_energy_1hd(tracer, pdel, cp_or_cv, U, V, T, &
end select
end if
deallocate(species_idx, species_liq_idx, species_ice_idx)
-
end subroutine get_hydrostatic_energy_1hd
!===========================================================================
diff --git a/src/data/cam_thermo_formula.F90 b/src/data/cam_thermo_formula.F90
new file mode 100644
index 00000000..df202c9d
--- /dev/null
+++ b/src/data/cam_thermo_formula.F90
@@ -0,0 +1,39 @@
+module cam_thermo_formula
+
+ use runtime_obj, only: unset_int
+
+ 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 = unset_int
+ ! energy_formula_physics: energy formula used for physics
+ integer, public :: energy_formula_physics = unset_int
+
+ ! Public subroutines
+ public :: cam_thermo_formula_init
+
+contains
+ subroutine cam_thermo_formula_init()
+ use phys_vars_init_check, only: mark_as_initialized
+
+ ! Physics energy formulation is always FV (moist pressure coordinate)
+ energy_formula_physics = ENERGY_FORMULA_DYCORE_FV
+ call mark_as_initialized("total_energy_formula_for_physics")
+
+ end subroutine cam_thermo_formula_init
+
+end module cam_thermo_formula
diff --git a/src/data/cam_thermo_formula.meta b/src/data/cam_thermo_formula.meta
new file mode 100644
index 00000000..f8bf04a1
--- /dev/null
+++ b/src/data/cam_thermo_formula.meta
@@ -0,0 +1,17 @@
+[ccpp-table-properties]
+ name = cam_thermo_formula
+ type = module
+
+[ccpp-arg-table]
+ name = cam_thermo_formula
+ type = module
+[ energy_formula_dycore ]
+ standard_name = total_energy_formula_for_dycore
+ units = 1
+ type = integer
+ dimensions = ()
+[ energy_formula_physics ]
+ standard_name = total_energy_formula_for_physics
+ units = 1
+ type = integer
+ dimensions = ()
diff --git a/src/data/registry.xml b/src/data/registry.xml
index 4a199b9f..91658e20 100644
--- a/src/data/registry.xml
+++ b/src/data/registry.xml
@@ -4,6 +4,7 @@
$SRCROOT/src/utils/spmd_utils.meta
+ $SRCROOT/src/control/cam_control_mod.meta
$SRCROOT/src/control/cam_logfile.meta
$SRCROOT/src/control/camsrfexch.meta
$SRCROOT/src/control/runtime_obj.meta
@@ -13,6 +14,7 @@
$SRCROOT/src/physics/utils/tropopause_climo_read.meta
$SRCROOT/src/data/air_composition.meta
$SRCROOT/src/data/cam_thermo.meta
+ $SRCROOT/src/data/cam_thermo_formula.meta
$SRCROOT/src/data/ref_pres.meta
$SRCROOT/src/dynamics/utils/vert_coord.meta
$SRCROOT/src/dynamics/utils/hycoef.meta
@@ -204,7 +206,7 @@
zi state_zi
@@ -212,7 +214,7 @@
te_ini_phys state_te_ini_phys
@@ -220,7 +222,7 @@
te_cur_phys state_te_cur_phys
@@ -228,7 +230,7 @@
te_ini_dyn state_te_ini_dyn
@@ -236,7 +238,7 @@
te_cur_dyn state_te_cur_dyn
@@ -244,13 +246,34 @@
tw_ini state_tw_ini
horizontal_dimension
tw_cur state_tw_cur
+
+ Total energy using dynamical core formula at the end of physics timestep
+ horizontal_dimension
+ 0.0
+
+
+ flag indicating if dynamical core energy is not consistent with CAM physics and to perform adjustment of temperature and temperature tendency
+ .false.
+ .true.
+
+
+
+ flag indicating if it is the first timestep of an initial run
+
reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure
frontogenesis_function
frontogenesis_angle
- vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula
- vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula
- vertically_integrated_water_vapor_and_condensed_water_of_initial_state
- vertically_integrated_water_vapor_and_condensed_water_of_current_state
+ vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep
+ vertically_integrated_total_energy_using_dycore_energy_formula
+ vertically_integrated_total_water_at_start_of_physics_timestep
+ vertically_integrated_total_water
+ vertically_integrated_total_energy_at_end_of_physics_timestep
tendency_of_air_temperature_due_to_model_physics
@@ -415,6 +439,14 @@
horizontal_dimension vertical_layer_dimension
zvir
+
+ specific heat of air used in the dynamical core (enthalpy for pressure-based dynamical cores and internal energy for z-based dynamical cores)
+ horizontal_dimension vertical_layer_dimension
+ cp_or_cv_dycore
+
dyn_mpas_compute_unit_vector
procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind
procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4
+ procedure, pass, public :: run => dyn_mpas_run
! Accessor subroutines for users to access internal states of MPAS dynamical core.
@@ -2525,19 +2530,21 @@ end subroutine dyn_mpas_compute_unit_vector
!> \date 16 January 2020
!> \details
!> This subroutine computes the edge-normal wind vectors at edge points
- !> (i.e., `u` in MPAS `state` pool) from wind components at cell points
+ !> (i.e., `u` in MPAS `state` pool) from the wind components at cell points
!> (i.e., `uReconstruct{Zonal,Meridional}` in MPAS `diag` pool). In MPAS, the
!> former are PROGNOSTIC variables, while the latter are DIAGNOSTIC variables
!> that are "reconstructed" from the former. This subroutine is essentially the
!> inverse function of that reconstruction. The purpose is to provide an
!> alternative way for MPAS to initialize from zonal and meridional wind
- !> components at cell points.
+ !> components at cell points. If `wind_tendency` is `.true.`, this subroutine
+ !> operates on the wind tendency due to physics instead.
!> \addenda
!> Ported and refactored for CAM-SIMA. (KCW, 2024-05-08)
!
!-------------------------------------------------------------------------------
- subroutine dyn_mpas_compute_edge_wind(self)
+ subroutine dyn_mpas_compute_edge_wind(self, wind_tendency)
class(mpas_dynamical_core_type), intent(in) :: self
+ logical, intent(in) :: wind_tendency
character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_edge_wind'
integer :: cell1, cell2, i
@@ -2561,14 +2568,24 @@ subroutine dyn_mpas_compute_edge_wind(self)
nullify(uedge)
! Make sure halo layers are up-to-date before computation.
- call self % exchange_halo('uReconstructZonal')
- call self % exchange_halo('uReconstructMeridional')
+ if (wind_tendency) then
+ call self % exchange_halo('tend_uzonal')
+ call self % exchange_halo('tend_umerid')
+ else
+ call self % exchange_halo('uReconstructZonal')
+ call self % exchange_halo('uReconstructMeridional')
+ end if
! Input.
call self % get_variable_pointer(nedges, 'dim', 'nEdges')
- call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
- call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
+ if (wind_tendency) then
+ call self % get_variable_pointer(ucellzonal, 'tend_physics', 'tend_uzonal')
+ call self % get_variable_pointer(ucellmeridional, 'tend_physics', 'tend_umerid')
+ else
+ call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
+ call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
+ end if
call self % get_variable_pointer(cellsonedge, 'mesh', 'cellsOnEdge')
call self % get_variable_pointer(east, 'mesh', 'east')
@@ -2576,7 +2593,11 @@ subroutine dyn_mpas_compute_edge_wind(self)
call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors')
! Output.
- call self % get_variable_pointer(uedge, 'state', 'u', time_level=1)
+ if (wind_tendency) then
+ call self % get_variable_pointer(uedge, 'tend_physics', 'tend_ru_physics')
+ else
+ call self % get_variable_pointer(uedge, 'state', 'u', time_level=1)
+ end if
do i = 1, nedges
cell1 = cellsonedge(1, i)
@@ -2607,7 +2628,11 @@ subroutine dyn_mpas_compute_edge_wind(self)
nullify(uedge)
! Make sure halo layers are up-to-date after computation.
- call self % exchange_halo('u')
+ if (wind_tendency) then
+ call self % exchange_halo('tend_ru_physics')
+ else
+ call self % exchange_halo('u')
+ end if
call self % debug_print(subname // ' completed')
end subroutine dyn_mpas_compute_edge_wind
@@ -2640,6 +2665,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
logical, pointer :: config_do_restart
real(rkind), pointer :: config_dt
type(field0dreal), pointer :: field_0d_real
+ type(field2dreal), pointer :: field_2d_real
type(mpas_pool_type), pointer :: mpas_pool
type(mpas_time_type) :: mpas_time
@@ -2651,6 +2677,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
nullify(config_do_restart)
nullify(config_dt)
nullify(field_0d_real)
+ nullify(field_2d_real)
nullify(mpas_pool)
if (coupling_time_interval <= 0) then
@@ -2673,6 +2700,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
end if
self % coupling_time_interval = coupling_time_interval
+ self % number_of_time_steps = 0
call self % debug_print('Coupling time interval is ' // stringify([real(self % coupling_time_interval, rkind)]) // &
' seconds')
@@ -2814,6 +2842,16 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
! Prepare dynamics for time integration.
call mpas_atm_dynamics_init(self % domain_ptr)
+ ! Some additional "scratch" fields are needed for interoperability with CAM-SIMA, but they are not initialized by
+ ! `mpas_atm_dynamics_init`. Initialize them below.
+ call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_uzonal', field_2d_real, timelevel=1)
+ call mpas_allocate_scratch_field(field_2d_real)
+ nullify(field_2d_real)
+
+ call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_umerid', field_2d_real, timelevel=1)
+ call mpas_allocate_scratch_field(field_2d_real)
+ nullify(field_2d_real)
+
call self % debug_print(subname // ' completed')
call self % debug_print('Successful initialization of MPAS dynamical core')
@@ -2867,6 +2905,113 @@ pure function almost_equal(a, b, absolute_tolerance, relative_tolerance)
end function almost_equal
end subroutine dyn_mpas_init_phase4
+ !-------------------------------------------------------------------------------
+ ! subroutine dyn_mpas_run
+ !
+ !> \brief Integrates the dynamical states with time
+ !> \author Michael Duda
+ !> \date 29 February 2020
+ !> \details
+ !> This subroutine calls MPAS dynamical solver in a loop, with each iteration
+ !> of the loop advancing the dynamical states forward by one time step, until
+ !> the coupling time interval is reached.
+ !> Essentially, it closely follows what is done in `atm_core_run`, but without
+ !> any calls to MPAS diagnostics manager or MPAS stream manager.
+ !> \addenda
+ !> Ported and refactored for CAM-SIMA. (KCW, 2024-06-21)
+ !
+ !-------------------------------------------------------------------------------
+ subroutine dyn_mpas_run(self)
+ class(mpas_dynamical_core_type), intent(inout) :: self
+
+ character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_run'
+ character(strkind) :: date_time
+ integer :: ierr
+ real(rkind), pointer :: config_dt
+ type(mpas_pool_type), pointer :: mpas_pool_diag, mpas_pool_mesh, mpas_pool_state
+ type(mpas_time_type) :: mpas_time_end, mpas_time_now ! This derived type is analogous to `ESMF_Time`.
+ type(mpas_timeinterval_type) :: mpas_time_interval ! This derived type is analogous to `ESMF_TimeInterval`.
+
+ call self % debug_print(subname // ' entered')
+
+ nullify(config_dt)
+ nullify(mpas_pool_diag, mpas_pool_mesh, mpas_pool_state)
+
+ call self % get_variable_pointer(config_dt, 'cfg', 'config_dt')
+ call self % get_pool_pointer(mpas_pool_diag, 'diag')
+ call self % get_pool_pointer(mpas_pool_mesh, 'mesh')
+ call self % get_pool_pointer(mpas_pool_state, 'state')
+
+ mpas_time_now = mpas_get_clock_time(self % domain_ptr % clock, mpas_now, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__)
+ end if
+
+ call mpas_get_time(mpas_time_now, datetimestring=date_time, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__)
+ end if
+
+ call self % debug_print('Time integration of MPAS dynamical core begins at ' // trim(adjustl(date_time)))
+
+ call mpas_set_timeinterval(mpas_time_interval, s=self % coupling_time_interval, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to set coupling time interval', subname, __LINE__)
+ end if
+
+ ! The `+` operator is overloaded here.
+ mpas_time_end = mpas_time_now + mpas_time_interval
+
+ ! Integrate until the coupling time interval is reached.
+ ! The `<` operator is overloaded here.
+ do while (mpas_time_now < mpas_time_end)
+ ! Number of time steps that has been completed in this MPAS dynamical core instance.
+ self % number_of_time_steps = self % number_of_time_steps + 1
+
+ ! Advance the dynamical states forward in time by `config_dt` seconds.
+ ! Current states are in time level 1. Upon exit, time level 2 will contain updated states.
+ call atm_do_timestep(self % domain_ptr, config_dt, self % number_of_time_steps)
+
+ ! MPAS `state` pool has two time levels.
+ ! Swap them after advancing a time step.
+ call mpas_pool_shift_time_levels(mpas_pool_state)
+
+ call mpas_advance_clock(self % domain_ptr % clock, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to advance clock', subname, __LINE__)
+ end if
+
+ mpas_time_now = mpas_get_clock_time(self % domain_ptr % clock, mpas_now, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__)
+ end if
+
+ call self % debug_print('Time step ' // stringify([self % number_of_time_steps]) // ' completed')
+ end do
+
+ call mpas_get_time(mpas_time_now, datetimestring=date_time, ierr=ierr)
+
+ if (ierr /= 0) then
+ call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__)
+ end if
+
+ call self % debug_print('Time integration of MPAS dynamical core ends at ' // trim(adjustl(date_time)))
+
+ ! Compute diagnostic variables like `pressure`, `rho` and `theta` from time level 1 of MPAS `state` pool
+ ! by calling upstream MPAS functionality.
+ call atm_compute_output_diagnostics(mpas_pool_state, 1, mpas_pool_diag, mpas_pool_mesh)
+
+ nullify(config_dt)
+ nullify(mpas_pool_diag, mpas_pool_mesh, mpas_pool_state)
+
+ call self % debug_print(subname // ' completed')
+ end subroutine dyn_mpas_run
+
!-------------------------------------------------------------------------------
! function dyn_mpas_get_constituent_name
!
@@ -3148,7 +3293,7 @@ subroutine dyn_mpas_get_pool_pointer(self, pool_pointer, pool_name)
pool_pointer => self % domain_ptr % configs
case ('dim')
pool_pointer => self % domain_ptr % blocklist % dimensions
- case ('diag', 'mesh', 'state', 'tend')
+ case ('diag', 'mesh', 'state', 'tend', 'tend_physics')
call mpas_pool_get_subpool(self % domain_ptr % blocklist % allstructs, trim(adjustl(pool_name)), pool_pointer)
case default
call self % model_error('Unsupported pool name "' // trim(adjustl(pool_name)) // '"', subname, __LINE__)
diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90
index b1892f47..3911af4f 100644
--- a/src/dynamics/mpas/dyn_comp.F90
+++ b/src/dynamics/mpas/dyn_comp.F90
@@ -9,7 +9,7 @@ module dyn_comp
thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, &
thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore
use cam_abortutils, only: check_allocate, endrun
- use cam_constituents, only: const_name, const_is_water_species, num_advected, readtrace
+ use cam_constituents, only: const_name, const_is_dry, const_is_water_species, num_advected, readtrace
use cam_control_mod, only: initial_run
use cam_field_read, only: cam_read_field
use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id
@@ -28,15 +28,16 @@ module dyn_comp
use time_manager, only: get_start_date, get_stop_date, get_step_size, get_run_duration, timemgr_get_calendar_cf
use vert_coord, only: pver, pverp
- ! Modules from CESM Share.
- use shr_file_mod, only: shr_file_getunit
- use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8
- use shr_pio_mod, only: shr_pio_getiosys
-
! Modules from CCPP.
use cam_ccpp_cap, only: cam_constituents_array
use ccpp_kinds, only: kind_phys
use phys_vars_init_check, only: mark_as_initialized, std_name_len
+ use physics_types, only: phys_state
+
+ ! Modules from CESM Share.
+ use shr_file_mod, only: shr_file_getunit
+ use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8
+ use shr_pio_mod, only: shr_pio_getiosys
! Modules from external libraries.
use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open
@@ -49,10 +50,12 @@ module dyn_comp
public :: dyn_export_t
public :: dyn_readnl
public :: dyn_init
- ! public :: dyn_run
+ public :: dyn_run
! public :: dyn_final
public :: dyn_debug_print
+ public :: dyn_exchange_constituent_state
+ public :: reverse
public :: mpas_dynamical_core
public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
public :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max
@@ -184,6 +187,10 @@ 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
+ use phys_vars_init_check, only: mark_as_initialized
+ use physics_types, only: dycore_energy_consistency_adjust
+
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 +207,15 @@ 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
+ call mark_as_initialized('total_energy_formula_for_dycore')
+
+ ! Dynamical core energy is not consistent with CAM physics and requires
+ ! temperature and temperature tendency adjustment at end of physics.
+ dycore_energy_consistency_adjust = .true.
+ call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment')
+
allocate(constituent_name(num_advected), stat=ierr)
call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__)
@@ -251,9 +267,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
! Perform default initialization for all constituents.
! Subsequently, they can be overridden depending on the namelist option (below) and
! the actual availability (checked and handled by MPAS).
- call dyn_debug_print('Calling set_default_constituent')
+ call dyn_debug_print('Calling dyn_exchange_constituent_state')
- call set_default_constituent()
+ call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.false.)
! Namelist option that controls if constituents are to be read from the file.
if (readtrace) then
@@ -357,12 +373,12 @@ subroutine set_analytic_initial_condition()
integer, allocatable :: global_grid_index(:)
real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :)
real(kind_r8), allocatable :: lat_rad(:), lon_rad(:)
- real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces.
- ! Dimension and vertical index orders follow CAM-SIMA convention.
- real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces.
- ! Dimension and vertical index orders follow MPAS convention.
+ real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces.
+ ! Dimension and vertical index orders follow CAM-SIMA convention.
+ real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces.
+ ! Dimension and vertical index orders follow MPAS convention.
- call init_shared_variable()
+ call init_shared_variables()
call set_mpas_state_u()
call set_mpas_state_w()
@@ -370,18 +386,14 @@ subroutine set_analytic_initial_condition()
call set_mpas_state_rho_theta()
call set_mpas_state_rho_base_theta_base()
- deallocate(global_grid_index)
- deallocate(lat_rad, lon_rad)
- deallocate(z_int)
-
- nullify(zgrid)
+ call final_shared_variables()
contains
!> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines.
!> (KCW, 2024-05-13)
- subroutine init_shared_variable()
- character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variable'
+ subroutine init_shared_variables()
+ character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variables'
integer :: ierr
- integer :: k
+ integer :: i
integer, pointer :: indextocellid(:)
real(kind_r8), pointer :: lat_deg(:), lon_deg(:)
@@ -429,17 +441,29 @@ subroutine init_shared_variable()
call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pverp
- z_int(:, k) = zgrid(pverp - k + 1, 1:ncells_solve)
+ do i = 1, ncells_solve
+ z_int(i, :) = reverse(zgrid(:, i))
end do
- end subroutine init_shared_variable
+ end subroutine init_shared_variables
+
+ !> Finalize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines.
+ !> (KCW, 2024-05-13)
+ subroutine final_shared_variables()
+ character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::final_shared_variables'
+
+ deallocate(global_grid_index)
+ deallocate(lat_rad, lon_rad)
+ deallocate(z_int)
+
+ nullify(zgrid)
+ end subroutine final_shared_variables
!> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces).
!> (KCW, 2024-05-13)
subroutine set_mpas_state_u()
character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u'
integer :: ierr
- integer :: k
+ integer :: i
real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :)
call dyn_debug_print('Setting MPAS state "u"')
@@ -457,8 +481,8 @@ subroutine set_mpas_state_u()
call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, u=buffer_2d_real)
! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pver
- ucellzonal(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1)
+ do i = 1, ncells_solve
+ ucellzonal(:, i) = reverse(buffer_2d_real(i, :))
end do
buffer_2d_real(:, :) = 0.0_kind_r8
@@ -466,15 +490,15 @@ subroutine set_mpas_state_u()
call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, v=buffer_2d_real)
! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pver
- ucellmeridional(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1)
+ do i = 1, ncells_solve
+ ucellmeridional(:, i) = reverse(buffer_2d_real(i, :))
end do
deallocate(buffer_2d_real)
nullify(ucellzonal, ucellmeridional)
- call mpas_dynamical_core % compute_edge_wind()
+ call mpas_dynamical_core % compute_edge_wind(.false.)
end subroutine set_mpas_state_u
!> Set MPAS state `w` (i.e., vertical velocity at cell interfaces).
@@ -505,7 +529,7 @@ subroutine set_mpas_state_scalars()
'water_vapor_mixing_ratio_wrt_dry_air'
character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_scalars'
- integer :: i, k
+ integer :: i, j
integer :: ierr
integer, allocatable :: constituent_index(:)
integer, pointer :: index_qv
@@ -531,12 +555,12 @@ subroutine set_mpas_state_scalars()
call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, q=buffer_3d_real, &
m_cnst=constituent_index)
- ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index.
- do i = 1, num_advected
- ! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pver
- scalars(i, k, 1:ncells_solve) = &
- buffer_3d_real(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i))
+ do i = 1, ncells_solve
+ ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index.
+ do j = 1, num_advected
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ scalars(j, :, i) = &
+ reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j)))
end do
end do
@@ -577,6 +601,8 @@ subroutine set_mpas_state_rho_theta()
real(kind_r8), allocatable :: qv_mid_col(:) ! Water vapor mixing ratio (kg/kg) at layer midpoints of each column.
real(kind_r8), allocatable :: t_mid(:, :) ! Temperature (K) at layer midpoints.
real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column.
+ ! Be advised that it is not virtual temperature.
+ ! See doi:10.5065/1DFH-6P97 and doi:10.1175/MWR-D-11-00215.1 for details.
real(kind_r8), pointer :: rho(:, :)
real(kind_r8), pointer :: theta(:, :)
real(kind_r8), pointer :: scalars(:, :, :)
@@ -606,8 +632,8 @@ subroutine set_mpas_state_rho_theta()
call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, t=buffer_2d_real)
! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pver
- t_mid(k, :) = buffer_2d_real(:, pver - k + 1)
+ do i = 1, ncells_solve
+ t_mid(:, i) = reverse(buffer_2d_real(i, :))
end do
deallocate(buffer_2d_real)
@@ -755,19 +781,83 @@ pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0)
end function theta_by_poisson_equation
end subroutine set_analytic_initial_condition
- !> Set default MPAS state `scalars` (i.e., constituents) in accordance with CCPP, which is a component of CAM-SIMA.
- !> (KCW, 2024-07-09)
- subroutine set_default_constituent()
- character(*), parameter :: subname = 'dyn_comp::set_default_constituent'
- integer :: i, k
+ !> Exchange and/or convert constituent states between CAM-SIMA and MPAS.
+ !> If `exchange` is `.true.` and `direction` is "e" or "export", set MPAS state `scalars` from physics state `constituents`.
+ !> If `exchange` is `.true.` and `direction` is "i" or "import", set physics state `constituents` from MPAS state `scalars`.
+ !> Think of it as "exporting/importing constituent states in CAM-SIMA to/from MPAS".
+ !> Otherwise, if `exchange` is `.false.`, no exchange is performed at all.
+ !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratio that has different
+ !> definitions between CAM-SIMA and MPAS (i.e., dry/moist).
+ !> Otherwise, if `conversion` is `.false.`, no conversion is performed at all.
+ !> This subroutine is intentionally designed to have these elaborate controls due to complications in CAM-SIMA.
+ !> Some procedures in CAM-SIMA expect constituent states to be dry, while the others expect them to be moist.
+ !> (KCW, 2024-09-26)
+ subroutine dyn_exchange_constituent_state(direction, exchange, conversion)
+ character(*), intent(in) :: direction
+ logical, intent(in) :: exchange
+ logical, intent(in) :: conversion
+
+ character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_state'
+ integer :: i, j
+ integer :: ierr
+ integer, allocatable :: is_water_species_index(:)
+ logical, allocatable :: is_conversion_needed(:)
+ logical, allocatable :: is_water_species(:)
real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory.
+ real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water species mixing ratios.
real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory.
- call dyn_debug_print('Setting default MPAS state "scalars"')
+ select case (trim(adjustl(direction)))
+ case ('e', 'export')
+ if (exchange) then
+ call dyn_debug_print('Setting MPAS state "scalars" from physics state "constituents"')
+ end if
+
+ if (conversion) then
+ call dyn_debug_print('Converting MPAS state "scalars"')
+ end if
+ case ('i', 'import')
+ if (exchange) then
+ call dyn_debug_print('Setting physics state "constituents" from MPAS state "scalars"')
+ end if
+
+ if (conversion) then
+ call dyn_debug_print('Converting physics state "constituents"')
+ end if
+ case default
+ call endrun('Unsupported exchange direction "' // trim(adjustl(direction)) // '"', subname, __LINE__)
+ end select
nullify(constituents)
nullify(scalars)
+ allocate(is_conversion_needed(num_advected), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'is_conversion_needed(num_advected)', &
+ 'dyn_comp', __LINE__)
+
+ allocate(is_water_species(num_advected), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'is_water_species(num_advected)', &
+ 'dyn_comp', __LINE__)
+
+ do j = 1, num_advected
+ ! All constituent mixing ratios in MPAS are dry.
+ ! Therefore, conversion in between is needed for any constituent mixing ratios that are not dry in CAM-SIMA.
+ is_conversion_needed(j) = .not. const_is_dry(j)
+ is_water_species(j) = const_is_water_species(j)
+ end do
+
+ allocate(is_water_species_index(count(is_water_species)), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'is_water_species_index(count(is_water_species))', &
+ 'dyn_comp', __LINE__)
+
+ allocate(sigma_all_q(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'sigma_all_q(pver)', &
+ 'dyn_comp', __LINE__)
+
constituents => cam_constituents_array()
if (.not. associated(constituents)) then
@@ -776,21 +866,70 @@ subroutine set_default_constituent()
call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
- ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index.
- do i = 1, num_advected
- ! Vertical index order is reversed between CAM-SIMA and MPAS.
- do k = 1, pver
- scalars(i, k, 1:ncells_solve) = &
- constituents(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i))
+ if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then
+ do i = 1, ncells_solve
+ if (conversion .and. any(is_conversion_needed)) then
+ ! The summation term of equation 8 in doi:10.1029/2017MS001257.
+ ! Using equation 7 here is not possible because it requires all constituent mixing ratio to be moist
+ ! on the RHS of it. There is no such guarantee in CAM-SIMA.
+ sigma_all_q(:) = reverse(phys_state % pdel(i, :) / phys_state % pdeldry(i, :))
+ end if
+
+ ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index.
+ do j = 1, num_advected
+ if (exchange) then
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ scalars(j, :, i) = &
+ reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j)))
+ end if
+
+ if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then
+ ! Equation 8 in doi:10.1029/2017MS001257.
+ scalars(j, :, i) = &
+ scalars(j, :, i) * sigma_all_q(:)
+ end if
+ end do
end do
- end do
+ else
+ is_water_species_index(:) = &
+ pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species)
+
+ do i = 1, ncells_solve
+ if (conversion .and. any(is_conversion_needed)) then
+ ! The summation term of equation 8 in doi:10.1029/2017MS001257.
+ sigma_all_q(:) = reverse(1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1))
+ end if
+
+ ! `j` is indexing into `constituents`, so it is regarded as constituent index.
+ do j = 1, num_advected
+ if (exchange) then
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ constituents(i, :, j) = &
+ reverse(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i))
+ end if
+
+ if (conversion .and. is_conversion_needed(j)) then
+ ! Equation 8 in doi:10.1029/2017MS001257.
+ constituents(i, :, j) = &
+ constituents(i, :, j) / sigma_all_q(:)
+ end if
+ end do
+ end do
+ end if
+
+ deallocate(is_conversion_needed)
+ deallocate(is_water_species)
+ deallocate(is_water_species_index)
+ deallocate(sigma_all_q)
nullify(constituents)
nullify(scalars)
- ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
- call mpas_dynamical_core % exchange_halo('scalars')
- end subroutine set_default_constituent
+ if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then
+ ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+ call mpas_dynamical_core % exchange_halo('scalars')
+ end if
+ end subroutine dyn_exchange_constituent_state
!> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized
!> to prevent physics from attempting to read them from a file. These variables are to be exchanged later
@@ -832,13 +971,48 @@ subroutine mark_variable_as_initialized()
do i = 1, num_advected
call mark_as_initialized(trim(adjustl(const_name(i))))
end do
+
+ call mark_as_initialized('specific_heat_of_air_used_in_dycore')
+
+ ! These energy variables are calculated by check_energy_timestep_init
+ ! but need to be marked here
+ call mark_as_initialized('vertically_integrated_total_energy_at_end_of_physics_timestep')
+ call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula')
+ call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep')
+ call mark_as_initialized('vertically_integrated_total_energy_using_physics_energy_formula')
+ call mark_as_initialized('vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep')
+ call mark_as_initialized('vertically_integrated_total_water')
+ call mark_as_initialized('vertically_integrated_total_water_at_start_of_physics_timestep')
end subroutine mark_variable_as_initialized
- ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`.
- ! subroutine dyn_run()
- ! end subroutine dyn_run
+ !> Run MPAS dynamical core to integrate the dynamical states with time.
+ !> (KCW, 2024-07-11)
+ subroutine dyn_run()
+ character(*), parameter :: subname = 'dyn_comp::dyn_run'
+
+ ! MPAS dynamical core will run until the coupling time interval is reached.
+ call mpas_dynamical_core % run()
+ end subroutine dyn_run
! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`.
! subroutine dyn_final()
! end subroutine dyn_final
+
+ !> Helper function for reversing the order of elements in `array`.
+ !> (KCW, 2024-07-17)
+ pure function reverse(array)
+ real(kind_r8), intent(in) :: array(:)
+ real(kind_r8) :: reverse(size(array))
+
+ integer :: n
+
+ n = size(array)
+
+ ! There is nothing to reverse.
+ if (n == 0) then
+ return
+ end if
+
+ reverse(:) = array(n:1:-1)
+ end function reverse
end module dyn_comp
diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90
new file mode 100644
index 00000000..a2066e16
--- /dev/null
+++ b/src/dynamics/mpas/dyn_coupling.F90
@@ -0,0 +1,662 @@
+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_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, &
+ constant_rd => rair, constant_rv => rh2o
+ use runtime_obj, only: cam_runtime_opts
+ use vert_coord, only: pver, pverp
+
+ ! Modules from CCPP.
+ use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+ use ccpp_kinds, only: kind_phys
+ use geopotential_temp, only: geopotential_temp_run
+ 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
+
+ ! Modules from CESM Share.
+ use shr_kind_mod, only: kind_cx => shr_kind_cx, kind_r8 => shr_kind_r8
+
+ implicit none
+
+ private
+ ! Provide APIs required by CAM-SIMA.
+ public :: dynamics_to_physics_coupling
+ public :: physics_to_dynamics_coupling
+contains
+ !> Perform one-way coupling from the dynamics output states to the physics input states.
+ !> The other coupling direction is implemented by its counterpart, `physics_to_dynamics_coupling`.
+ !> (KCW, 2024-07-31)
+ subroutine dynamics_to_physics_coupling()
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling'
+ integer :: column_index
+ integer, allocatable :: is_water_species_index(:)
+ integer, pointer :: index_qv
+ ! Variable name suffixes have the following meanings:
+ ! `*_col`: Variable is of each column.
+ ! `*_int`: Variable is at layer interfaces.
+ ! `*_mid`: Variable is at layer midpoints.
+ real(kind_r8), allocatable :: pd_int_col(:), & ! Dry hydrostatic air pressure (Pa).
+ pd_mid_col(:), & ! Dry hydrostatic air pressure (Pa).
+ p_int_col(:), & ! Full hydrostatic air pressure (Pa).
+ p_mid_col(:), & ! Full hydrostatic air pressure (Pa).
+ z_int_col(:) ! Geometric height (m).
+ real(kind_r8), allocatable :: dpd_col(:), & ! Dry air pressure difference (Pa) between layer interfaces.
+ dp_col(:), & ! Full air pressure difference (Pa) between layer interfaces.
+ dz_col(:) ! Geometric height difference (m) between layer interfaces.
+ real(kind_r8), allocatable :: qv_mid_col(:), & ! Water vapor mixing ratio (kg kg-1).
+ sigma_all_q_mid_col(:) ! Summation of all water species mixing ratios (kg kg-1).
+ real(kind_r8), allocatable :: rhod_mid_col(:), & ! Dry air density (kg m-3).
+ rho_mid_col(:) ! Full air density (kg m-3).
+ real(kind_r8), allocatable :: t_mid_col(:), & ! Temperature (K).
+ tm_mid_col(:), & ! Modified "moist" temperature (K).
+ tv_mid_col(:) ! Virtual temperature (K).
+ real(kind_r8), allocatable :: u_mid_col(:), & ! Eastward wind velocity (m s-1).
+ v_mid_col(:), & ! Northward wind velocity (m s-1).
+ omega_mid_col(:) ! Vertical wind velocity (Pa s-1).
+ real(kind_r8), pointer :: exner(:, :)
+ real(kind_r8), pointer :: rho_zz(:, :)
+ real(kind_r8), pointer :: scalars(:, :, :)
+ real(kind_r8), pointer :: theta_m(:, :)
+ real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :)
+ real(kind_r8), pointer :: zgrid(:, :)
+ real(kind_r8), pointer :: zz(:, :)
+
+ call init_shared_variables()
+
+ call dyn_exchange_constituent_state(direction='i', exchange=.true., conversion=.false.)
+
+ call dyn_debug_print('Setting physics state variables column by column')
+
+ ! Set variables in the `physics_state` derived type column by column.
+ ! This way, peak memory usage can be reduced.
+ do column_index = 1, ncells_solve
+ call update_shared_variables(column_index)
+ call set_physics_state_column(column_index)
+ end do
+
+ call set_physics_state_external()
+
+ call final_shared_variables()
+ contains
+ !> Initialize variables that are shared and repeatedly used by the `update_shared_variables` and
+ !> `set_physics_state_column` internal subroutines.
+ !> (KCW, 2024-07-20)
+ subroutine init_shared_variables()
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variables'
+ integer :: i
+ integer :: ierr
+ logical, allocatable :: is_water_species(:)
+
+ call dyn_debug_print('Preparing for dynamics-physics coupling')
+
+ nullify(index_qv)
+ nullify(exner)
+ nullify(rho_zz)
+ nullify(scalars)
+ nullify(theta_m)
+ nullify(ucellzonal, ucellmeridional, w)
+ nullify(zgrid)
+ nullify(zz)
+
+ allocate(is_water_species(num_advected), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'is_water_species(num_advected)', &
+ 'dyn_coupling', __LINE__)
+
+ do i = 1, num_advected
+ is_water_species(i) = const_is_water_species(i)
+ end do
+
+ allocate(is_water_species_index(count(is_water_species)), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'is_water_species_index(count(is_water_species))', &
+ 'dyn_coupling', __LINE__)
+
+ is_water_species_index(:) = &
+ pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species)
+
+ deallocate(is_water_species)
+
+ allocate(pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(dpd_col(pver), dp_col(pver), dz_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'dpd_col(pver), dp_col(pver), dz_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(qv_mid_col(pver), sigma_all_q_mid_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'qv_mid_col(pver), sigma_all_q_mid_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(rhod_mid_col(pver), rho_mid_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'rhod_mid_col(pver), rho_mid_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(t_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 't_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
+ call mpas_dynamical_core % get_variable_pointer(exner, 'diag', 'exner')
+ call mpas_dynamical_core % get_variable_pointer(rho_zz, 'state', 'rho_zz', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
+ call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
+ call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
+ call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz')
+ end subroutine init_shared_variables
+
+ !> Finalize variables that are shared and repeatedly used by the `update_shared_variables` and
+ !> `set_physics_state_column` internal subroutines.
+ !> (KCW, 2024-07-20)
+ subroutine final_shared_variables()
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::final_shared_variables'
+
+ deallocate(is_water_species_index)
+ deallocate(pd_int_col, pd_mid_col, p_int_col, p_mid_col, z_int_col)
+ deallocate(dpd_col, dp_col, dz_col)
+ deallocate(qv_mid_col, sigma_all_q_mid_col)
+ deallocate(rhod_mid_col, rho_mid_col)
+ deallocate(t_mid_col, tm_mid_col, tv_mid_col)
+ deallocate(u_mid_col, v_mid_col, omega_mid_col)
+
+ nullify(index_qv)
+ nullify(exner)
+ nullify(rho_zz)
+ nullify(scalars)
+ nullify(theta_m)
+ nullify(ucellzonal, ucellmeridional, w)
+ nullify(zgrid)
+ nullify(zz)
+ end subroutine final_shared_variables
+
+ !> Update variables for the specific column, indicated by `i`. This subroutine and `set_physics_state_column`
+ !> should be called in pairs.
+ !> (KCW, 2024-07-30)
+ subroutine update_shared_variables(i)
+ integer, intent(in) :: i
+
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variables'
+ integer :: k
+
+ ! The summation term of equation 5 in doi:10.1029/2017MS001257.
+ sigma_all_q_mid_col(:) = 1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1)
+
+ ! Compute thermodynamic variables.
+
+ ! By definition.
+ z_int_col(:) = zgrid(:, i)
+ dz_col(:) = z_int_col(2:pverp) - z_int_col(1:pver)
+ qv_mid_col(:) = scalars(index_qv, :, i)
+ rhod_mid_col(:) = rho_zz(:, i) * zz(:, i)
+
+ ! Equation 5 in doi:10.1029/2017MS001257.
+ rho_mid_col(:) = rhod_mid_col(:) * sigma_all_q_mid_col(:)
+
+ ! Hydrostatic equation.
+ dpd_col(:) = -rhod_mid_col(:) * constant_g * dz_col(:)
+ dp_col(:) = -rho_mid_col(:) * constant_g * dz_col(:)
+
+ ! By definition of Exner function. Also see below.
+ tm_mid_col(:) = theta_m(:, i) * exner(:, i)
+
+ ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
+ ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
+ t_mid_col(:) = tm_mid_col(:) / &
+ (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:))
+
+ ! Equation 16 in doi:10.1029/2017MS001257.
+ ! The numerator terms are just `tm_mid_col` here (i.e., modified "moist" temperature).
+ tv_mid_col(:) = tm_mid_col(:) / sigma_all_q_mid_col(:)
+
+ ! Equation of state.
+ pd_mid_col(:) = rhod_mid_col(:) * constant_rd * t_mid_col(:)
+ p_mid_col(:) = rho_mid_col(:) * constant_rd * tv_mid_col(:)
+
+ ! By definition.
+ pd_int_col(pverp) = pd_mid_col(pver) + 0.5_kind_r8 * dpd_col(pver)
+ p_int_col(pverp) = p_mid_col(pver) + 0.5_kind_r8 * dp_col(pver)
+
+ ! Integrate downward.
+ do k = pver, 1, -1
+ pd_int_col(k) = pd_int_col(k + 1) - dpd_col(k)
+ p_int_col(k) = p_int_col(k + 1) - dp_col(k)
+ end do
+
+ ! Compute momentum variables.
+
+ ! By definition.
+ u_mid_col(:) = ucellzonal(:, i)
+ v_mid_col(:) = ucellmeridional(:, i)
+ omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * (w(1:pver, i) + w(2:pverp, i))
+ end subroutine update_shared_variables
+
+ !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type.
+ !> This subroutine and `update_shared_variables` should be called in pairs.
+ !> (KCW, 2024-07-30)
+ subroutine set_physics_state_column(i)
+ integer, intent(in) :: i
+
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_column'
+
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ ! Always call `reverse` when assigning anything to/from the `physics_state` derived type.
+
+ phys_state % u(i, :) = reverse(u_mid_col)
+ phys_state % v(i, :) = reverse(v_mid_col)
+ phys_state % omega(i, :) = reverse(omega_mid_col)
+
+ phys_state % psdry(i) = pd_int_col(1)
+ phys_state % pintdry(i, :) = reverse(pd_int_col)
+ phys_state % pmiddry(i, :) = reverse(pd_mid_col)
+ phys_state % pdeldry(i, :) = reverse(-dpd_col)
+ phys_state % lnpintdry(i, :) = log(phys_state % pintdry(i, :))
+ phys_state % lnpmiddry(i, :) = log(phys_state % pmiddry(i, :))
+ phys_state % rpdeldry(i, :) = 1.0_kind_r8 / phys_state % pdeldry(i, :)
+
+ phys_state % ps(i) = p_int_col(1)
+ phys_state % pint(i, :) = reverse(p_int_col)
+ phys_state % pmid(i, :) = reverse(p_mid_col)
+ phys_state % pdel(i, :) = reverse(-dp_col)
+ phys_state % lnpint(i, :) = log(phys_state % pint(i, :))
+ phys_state % lnpmid(i, :) = log(phys_state % pmid(i, :))
+ phys_state % rpdel(i, :) = 1.0_kind_r8 / phys_state % pdel(i, :)
+
+ phys_state % t(i, :) = reverse(t_mid_col)
+
+ phys_state % phis(i) = constant_g * z_int_col(1)
+ end subroutine set_physics_state_column
+
+ !> Set variables in the `physics_state` derived type by calling external procedures.
+ !> (KCW, 2024-07-30)
+ subroutine set_physics_state_external()
+ character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external'
+ character(kind_cx) :: cerr
+ integer :: i
+ integer :: ierr
+ real(kind_phys), allocatable :: minimum_constituents(:)
+ real(kind_phys), pointer :: constituents(:, :, :)
+ type(ccpp_constituent_prop_ptr_t), pointer :: constituent_properties(:)
+
+ call dyn_debug_print('Setting physics state variables externally')
+
+ nullify(constituents)
+ nullify(constituent_properties)
+
+ allocate(minimum_constituents(num_advected), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'minimum_constituents(num_advected)', &
+ 'dyn_coupling', __LINE__)
+
+ do i = 1, num_advected
+ minimum_constituents(i) = const_qmin(i)
+ end do
+
+ constituents => cam_constituents_array()
+
+ if (.not. associated(constituents)) then
+ call endrun('Failed to find variable "constituents"', subname, __LINE__)
+ end if
+
+ constituent_properties => cam_model_const_properties()
+
+ if (.not. associated(constituent_properties)) then
+ call endrun('Failed to find variable "constituent_properties"', subname, __LINE__)
+ end if
+
+ ! 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_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_dry_air_update`.
+ do i = 1, ncells_solve
+ phys_state % exner(i, :) = (phys_state % ps(i) / phys_state % pmid(i, :)) ** cappav(i, :)
+ end do
+
+ ! Note that constituents become moist after this.
+ call dyn_exchange_constituent_state(direction='i', exchange=.false., conversion=.true.)
+
+ ! Impose minimum limits on constituents.
+ call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr)
+
+ if (ierr /= 0) then
+ call endrun('Failed to impose minimum limits on constituents externally' // new_line('') // &
+ 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
+ subname, __LINE__)
+ end if
+
+ ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints).
+ ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_update`.
+ call geopotential_temp_run( &
+ pver, lagrangian_vertical, pver, 1, pverp, 1, num_advected, &
+ phys_state % lnpint, phys_state % pint, phys_state % pmid, phys_state % pdel, phys_state % rpdel, phys_state % t, &
+ constituents(:, :, mpas_dynamical_core % map_constituent_index(index_qv)), constituents, &
+ constituent_properties, rairv, constant_g, zvirv, phys_state % zi, phys_state % zm, ncells_solve, ierr, cerr)
+
+ if (ierr /= 0) then
+ call endrun('Failed to set variable "zi" and "zm" externally' // new_line('') // &
+ 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
+ subname, __LINE__)
+ end if
+
+ ! Set `dse` (i.e., dry static energy).
+ ! Note that `cpairv` is updated externally by `cam_thermo_update`.
+ call update_dry_static_energy_run( &
+ pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr)
+
+ if (ierr /= 0) then
+ call endrun('Failed to set variable "dse" externally' // new_line('') // &
+ 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
+ subname, __LINE__)
+ end if
+
+ deallocate(minimum_constituents)
+
+ nullify(constituents)
+ nullify(constituent_properties)
+ end subroutine set_physics_state_external
+ end subroutine dynamics_to_physics_coupling
+
+ !> Perform one-way coupling from the physics output states to the dynamics input states.
+ !> The other coupling direction is implemented by its counterpart, `dynamics_to_physics_coupling`.
+ !> (KCW, 2024-09-20)
+ subroutine physics_to_dynamics_coupling()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling'
+ integer, pointer :: index_qv
+ real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1)
+ ! before being updated by physics.
+ real(kind_r8), pointer :: rho_zz(:, :)
+ real(kind_r8), pointer :: scalars(:, :, :)
+ real(kind_r8), pointer :: zz(:, :)
+
+ call init_shared_variables()
+
+ call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.true.)
+
+ call set_mpas_physics_tendency_ru()
+ call set_mpas_physics_tendency_rho()
+ call set_mpas_physics_tendency_rtheta()
+
+ call final_shared_variables()
+ contains
+ !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines.
+ !> (KCW, 2024-09-13)
+ subroutine init_shared_variables()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variables'
+ integer :: ierr
+
+ call dyn_debug_print('Preparing for physics-dynamics coupling')
+
+ nullify(index_qv)
+ nullify(rho_zz)
+ nullify(scalars)
+ nullify(zz)
+
+ call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
+ call mpas_dynamical_core % get_variable_pointer(rho_zz, 'state', 'rho_zz', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz')
+
+ allocate(qv_prev(pver, ncells_solve), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'qv_prev(pver, ncells_solve)', &
+ 'dyn_coupling', __LINE__)
+
+ ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta`
+ ! needs it. This must be done before calling `dyn_exchange_constituent_state`.
+ qv_prev(:, :) = scalars(index_qv, :, 1:ncells_solve)
+ end subroutine init_shared_variables
+
+ !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines.
+ !> (KCW, 2024-09-13)
+ subroutine final_shared_variables()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::final_shared_variables'
+
+ deallocate(qv_prev)
+
+ nullify(index_qv)
+ nullify(rho_zz)
+ nullify(scalars)
+ nullify(zz)
+ end subroutine final_shared_variables
+
+ !> Set MPAS physics tendency `tend_ru_physics` (i.e., "coupled" tendency of horizontal velocity at edge interfaces
+ !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
+ !> (KCW, 2024-09-11)
+ subroutine set_mpas_physics_tendency_ru()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru'
+ integer :: i
+ real(kind_r8), pointer :: u_tendency(:, :), v_tendency(:, :)
+
+ call dyn_debug_print('Setting MPAS physics tendency "tend_ru_physics"')
+
+ nullify(u_tendency, v_tendency)
+
+ call mpas_dynamical_core % get_variable_pointer(u_tendency, 'tend_physics', 'tend_uzonal')
+ call mpas_dynamical_core % get_variable_pointer(v_tendency, 'tend_physics', 'tend_umerid')
+
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type.
+ do i = 1, ncells_solve
+ u_tendency(:, i) = reverse(phys_tend % dudt_total(i, :)) * rho_zz(:, i)
+ v_tendency(:, i) = reverse(phys_tend % dvdt_total(i, :)) * rho_zz(:, i)
+ end do
+
+ nullify(u_tendency, v_tendency)
+
+ call mpas_dynamical_core % compute_edge_wind(.true.)
+ end subroutine set_mpas_physics_tendency_ru
+
+ !> Set MPAS physics tendency `tend_rho_physics` (i.e., "coupled" tendency of dry air density due to physics).
+ !> In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
+ !> (KCW, 2024-09-11)
+ subroutine set_mpas_physics_tendency_rho()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho'
+ real(kind_r8), pointer :: rho_tendency(:, :)
+
+ call dyn_debug_print('Setting MPAS physics tendency "tend_rho_physics"')
+
+ nullify(rho_tendency)
+
+ call mpas_dynamical_core % get_variable_pointer(rho_tendency, 'tend_physics', 'tend_rho_physics')
+
+ ! The material derivative of `rho` (i.e., dry air density) is zero for incompressible fluid.
+ rho_tendency(:, 1:ncells_solve) = 0.0_kind_r8
+
+ nullify(rho_tendency)
+
+ ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+ call mpas_dynamical_core % exchange_halo('tend_rho_physics')
+ end subroutine set_mpas_physics_tendency_rho
+
+ !> Set MPAS physics tendency `tend_rtheta_physics` (i.e., "coupled" tendency of modified "moist" potential temperature
+ !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
+ !> (KCW, 2024-09-19)
+ subroutine set_mpas_physics_tendency_rtheta()
+ character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rtheta'
+ integer :: i
+ integer :: ierr
+ ! Variable name suffixes have the following meanings:
+ ! `*_col`: Variable is of each column.
+ ! `*_prev`: Variable is before being updated by physics.
+ ! `*_curr`: Variable is after being updated by physics.
+ real(kind_r8), allocatable :: qv_col_prev(:), qv_col_curr(:) ! Water vapor mixing ratio (kg kg-1).
+ real(kind_r8), allocatable :: rhod_col(:) ! Dry air density (kg m-3).
+ real(kind_r8), allocatable :: t_col_prev(:), t_col_curr(:) ! Temperature (K).
+ real(kind_r8), allocatable :: theta_col_prev(:), theta_col_curr(:) ! Potential temperature (K).
+ real(kind_r8), allocatable :: thetam_col_prev(:), thetam_col_curr(:) ! Modified "moist" potential temperature (K).
+ real(kind_r8), pointer :: theta_m(:, :)
+ real(kind_r8), pointer :: theta_m_tendency(:, :)
+
+ call dyn_debug_print('Setting MPAS physics tendency "tend_rtheta_physics"')
+
+ nullify(theta_m)
+ nullify(theta_m_tendency)
+
+ allocate(qv_col_prev(pver), qv_col_curr(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'qv_col_prev(pver), qv_col_curr(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(rhod_col(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'rhod_col(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(t_col_prev(pver), t_col_curr(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 't_col_prev(pver), t_col_curr(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(theta_col_prev(pver), theta_col_curr(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'theta_col_prev(pver), theta_col_curr(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ allocate(thetam_col_prev(pver), thetam_col_curr(pver), stat=ierr)
+ call check_allocate(ierr, subname, &
+ 'thetam_col_prev(pver), thetam_col_curr(pver)', &
+ 'dyn_coupling', __LINE__)
+
+ call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1)
+ call mpas_dynamical_core % get_variable_pointer(theta_m_tendency, 'tend_physics', 'tend_rtheta_physics')
+
+ ! Set `theta_m_tendency` column by column. This way, peak memory usage can be reduced.
+ do i = 1, ncells_solve
+ qv_col_curr(:) = scalars(index_qv, :, i)
+ qv_col_prev(:) = qv_prev(:, i)
+ rhod_col(:) = rho_zz(:, i) * zz(:, i)
+
+ thetam_col_prev(:) = theta_m(:, i)
+ theta_col_prev(:) = thetam_col_prev(:) / (1.0_kind_r8 + constant_rv / constant_rd * qv_col_prev(:))
+ t_col_prev(:) = t_of_theta_rhod_qv(theta_col_prev, rhod_col, qv_col_prev)
+
+ ! Vertical index order is reversed between CAM-SIMA and MPAS.
+ ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type.
+ t_col_curr(:) = t_col_prev(:) + reverse(phys_tend % dtdt_total(i, :)) * dtime_phys
+ theta_col_curr(:) = theta_of_t_rhod_qv(t_col_curr, rhod_col, qv_col_curr)
+ thetam_col_curr(:) = theta_col_curr(:) * (1.0_kind_r8 + constant_rv / constant_rd * qv_col_curr(:))
+
+ theta_m_tendency(:, i) = (thetam_col_curr(:) - thetam_col_prev(:)) * rho_zz(:, i) / dtime_phys
+ end do
+
+ deallocate(qv_col_prev, qv_col_curr)
+ deallocate(rhod_col)
+ deallocate(t_col_prev, t_col_curr)
+ deallocate(theta_col_prev, theta_col_curr)
+ deallocate(thetam_col_prev, thetam_col_curr)
+
+ nullify(theta_m)
+ nullify(theta_m_tendency)
+
+ ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
+ call mpas_dynamical_core % exchange_halo('tend_rtheta_physics')
+ end subroutine set_mpas_physics_tendency_rtheta
+
+ !> Compute temperature `t` as a function of potential temperature `theta`, dry air density `rhod` and water vapor
+ !> mixing ratio `qv`. The formulation comes from Poisson equation with equation of state plugged in and arranging
+ !> for temperature. This function is the exact inverse of `theta_of_t_rhod_qv`, which means that:
+ !> `t == t_of_theta_rhod_qv(theta_of_t_rhod_qv(t, rhod, qv), rhod, qv)`.
+ !> (KCW, 2024-09-13)
+ pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t)
+ real(kind_r8), intent(in) :: theta, rhod, qv
+ real(kind_r8) :: t
+
+ real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume.
+
+ ! Mayer's relation.
+ constant_cvd = constant_cpd - constant_rd
+
+ ! Poisson equation with equation of state plugged in and arranging for temperature. For equation of state,
+ ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that
+ ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is
+ ! described herein:
+ ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
+ ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
+ !
+ ! In short, solve the below equation set for $T$ in terms of $\theta$, $\rho_d$ and $q_v$:
+ ! \begin{equation*}
+ ! \begin{cases}
+ ! \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
+ ! P &= \rho_d R_d T_m \\
+ ! T_m &= T (1 + \frac{R_v}{R_d} q_v)
+ ! \end{cases}
+ ! \end{equation*}
+ t = (theta ** (constant_cpd / constant_cvd)) * &
+ (((rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv)) / constant_p0) ** &
+ (constant_rd / constant_cvd))
+ end function t_of_theta_rhod_qv
+
+ !> Compute potential temperature `theta` as a function of temperature `t`, dry air density `rhod` and water vapor
+ !> mixing ratio `qv`. The formulation comes from Poisson equation with equation of state plugged in and arranging
+ !> for potential temperature. This function is the exact inverse of `t_of_theta_rhod_qv`, which means that:
+ !> `theta == theta_of_t_rhod_qv(t_of_theta_rhod_qv(theta, rhod, qv), rhod, qv)`.
+ !> (KCW, 2024-09-13)
+ pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta)
+ real(kind_r8), intent(in) :: t, rhod, qv
+ real(kind_r8) :: theta
+
+ real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume.
+
+ ! Mayer's relation.
+ constant_cvd = constant_cpd - constant_rd
+
+ ! Poisson equation with equation of state plugged in and arranging for potential temperature. For equation of state,
+ ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that
+ ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is
+ ! described herein:
+ ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
+ ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
+ !
+ ! In short, solve the below equation set for $\theta$ in terms of $T$, $\rho_d$ and $q_v$:
+ ! \begin{equation*}
+ ! \begin{cases}
+ ! \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
+ ! P &= \rho_d R_d T_m \\
+ ! T_m &= T (1 + \frac{R_v}{R_d} q_v)
+ ! \end{cases}
+ ! \end{equation*}
+ theta = (t ** (constant_cvd / constant_cpd)) * &
+ ((constant_p0 / (rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv))) ** &
+ (constant_rd / constant_cpd))
+ end function theta_of_t_rhod_qv
+ end subroutine physics_to_dynamics_coupling
+end module dyn_coupling
diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90
index 9d53ca43..d089636e 100644
--- a/src/dynamics/mpas/stepon.F90
+++ b/src/dynamics/mpas/stepon.F90
@@ -1,9 +1,14 @@
module stepon
- use camsrfexch, only: cam_out_t
- use dyn_comp, only: dyn_import_t, dyn_export_t
+ ! Modules from CAM-SIMA.
+ use camsrfexch, only: cam_out_t
+ use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run
+ use dyn_coupling, only: dynamics_to_physics_coupling, physics_to_dynamics_coupling
use physics_types, only: physics_state, physics_tend
- use runtime_obj, only: runtime_options
- use shr_kind_mod, only: r8 => shr_kind_r8
+ use runtime_obj, only: runtime_options
+ use time_manager, only: get_step_size
+
+ ! Modules from CCPP.
+ use ccpp_kinds, only: kind_phys
implicit none
@@ -15,48 +20,55 @@ module stepon
public :: stepon_run3
public :: stepon_final
contains
+ ! Called by `cam_init` in `src/control/cam_comp.F90`.
+ subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out)
+ type(runtime_options), intent(in) :: cam_runtime_opts
+ type(dyn_import_t), intent(in) :: dyn_in
+ type(dyn_export_t), intent(in) :: dyn_out
+ end subroutine stepon_init
+
+ ! Called by `cam_timestep_init` in `src/control/cam_comp.F90`.
+ subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out)
+ real(kind_phys), intent(out) :: dtime_phys
+ type(runtime_options), intent(in) :: cam_runtime_opts
+ type(physics_state), intent(in) :: phys_state
+ type(physics_tend), intent(in) :: phys_tend
+ type(dyn_import_t), intent(in) :: dyn_in
+ type(dyn_export_t), intent(in) :: dyn_out
+
+ ! Set timestep for physics.
+ dtime_phys = real(get_step_size(), kind_phys)
+
+ call dynamics_to_physics_coupling()
+ end subroutine stepon_timestep_init
+
+ ! Called by `cam_run2` in `src/control/cam_comp.F90`.
+ subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out)
+ type(runtime_options), intent(in) :: cam_runtime_opts
+ type(physics_state), intent(in) :: phys_state
+ type(physics_tend), intent(in) :: phys_tend
+ type(dyn_import_t), intent(in) :: dyn_in
+ type(dyn_export_t), intent(in) :: dyn_out
+
+ call physics_to_dynamics_coupling()
+ end subroutine stepon_run2
+
+ ! Called by `cam_run3` in `src/control/cam_comp.F90`.
+ subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out)
+ real(kind_phys), intent(in) :: dtime_phys
+ type(runtime_options), intent(in) :: cam_runtime_opts
+ type(cam_out_t), intent(in) :: cam_out
+ type(physics_state), intent(in) :: phys_state
+ type(dyn_import_t), intent(in) :: dyn_in
+ type(dyn_export_t), intent(in) :: dyn_out
-! Called by `cam_init` in `src/control/cam_comp.F90`.
-subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out)
- type(runtime_options), intent(in) :: cam_runtime_opts
- type(dyn_import_t), intent(in) :: dyn_in
- type(dyn_export_t), intent(in) :: dyn_out
-end subroutine stepon_init
-
-! Called by `cam_timestep_init` in `src/control/cam_comp.F90`.
-subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out)
- real(r8), intent(out) :: dtime_phys
- type(runtime_options), intent(in) :: cam_runtime_opts
- type(physics_state), intent(inout) :: phys_state
- type(physics_tend), intent(inout) :: phys_tend
- type(dyn_import_t), intent(inout) :: dyn_in
- type(dyn_export_t), intent(inout) :: dyn_out
-end subroutine stepon_timestep_init
-
-! Called by `cam_run2` in `src/control/cam_comp.F90`.
-subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out)
- type(runtime_options), intent(in) :: cam_runtime_opts
- type(physics_state), intent(inout) :: phys_state
- type(physics_tend), intent(inout) :: phys_tend
- type(dyn_import_t), intent(inout) :: dyn_in
- type(dyn_export_t), intent(inout) :: dyn_out
-end subroutine stepon_run2
-
-! Called by `cam_run3` in `src/control/cam_comp.F90`.
-subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out)
- real(r8), intent(in) :: dtime_phys
- type(runtime_options), intent(in) :: cam_runtime_opts
- type(cam_out_t), intent(inout) :: cam_out
- type(physics_state), intent(inout) :: phys_state
- type(dyn_import_t), intent(inout) :: dyn_in
- type(dyn_export_t), intent(inout) :: dyn_out
-end subroutine stepon_run3
-
-! Called by `cam_final` in `src/control/cam_comp.F90`.
-subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out)
- type(runtime_options), intent(in) :: cam_runtime_opts
- type(dyn_import_t), intent(inout) :: dyn_in
- type(dyn_export_t), intent(inout) :: dyn_out
-end subroutine stepon_final
+ call dyn_run()
+ end subroutine stepon_run3
+ ! Called by `cam_final` in `src/control/cam_comp.F90`.
+ subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out)
+ type(runtime_options), intent(in) :: cam_runtime_opts
+ type(dyn_import_t), intent(in) :: dyn_in
+ type(dyn_export_t), intent(in) :: dyn_out
+ end subroutine stepon_final
end module stepon
diff --git a/src/dynamics/none/dyn_comp.F90 b/src/dynamics/none/dyn_comp.F90
index 968e04e2..9ecb3022 100644
--- a/src/dynamics/none/dyn_comp.F90
+++ b/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
!==============================================================================
diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index bc714e22..ba1bf0ba 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/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
@@ -126,6 +127,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 +161,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 +633,78 @@ 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
+ use pio, only: pio_inq_att, pio_global, 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
+ use physics_types, only: dycore_energy_consistency_adjust
+ use phys_vars_init_check, only: mark_as_initialized
+
+ ! Find which dynamical core is used in and set the energy formulation
+ ! (also called vc_dycore in CAM)
+ !
+ ! This functionality is only used to recognize the originating dynamical core
+ ! from the snapshot file in order to set the energy formulation when running
+ ! with the null dycore. Other dynamical cores set energy_formula_dycore at their
+ ! initialization.
+
+ type(file_desc_t), intent(inout) :: file
+ logical, intent(in) :: grid_is_latlon
+
+ ! Local variables
+ 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
+ dycore_energy_consistency_adjust = .false.
+ if(masterproc) then
+ write(iulog, *) subname, ': Null dycore will use FV dycore energy formula'
+ endif
+ else
+ ! Is SE dycore?
+ ierr = pio_inq_att(file, pio_global, '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
+ dycore_energy_consistency_adjust = .true.
+ 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
+ dycore_energy_consistency_adjust = .true.
+ if(masterproc) then
+ write(iulog, *) subname, ': Null dycore will use MPAS dycore energy formula'
+ endif
+ endif
+ endif
+
+ 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
+ call mark_as_initialized("specific_heat_of_air_used_in_dycore")
+ call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula")
+ call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula")
+ call mark_as_initialized("vertically_integrated_total_water_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_water")
+ call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep")
+
+ end subroutine find_energy_formula
+
end module dyn_grid
diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90
index 8b56e9d9..572f663f 100644
--- a/src/dynamics/se/dp_coupling.F90
+++ b/src/dynamics/se/dp_coupling.F90
@@ -582,18 +582,22 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
use cam_constituents, only: const_qmin
use runtime_obj, only: wv_stdname
use physics_types, only: lagrangian_vertical
- use physconst, only: cpair, gravit, zvir, cappa
- use cam_thermo, only: cam_thermo_update
- use physics_types, only: cpairv, rairv, zvirv
+ use physconst, only: cpair, gravit, zvir
+ use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update
+ use air_composition, only: thermodynamic_active_species_num
+ 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, only: check_energy_timestep_init
use hycoef, only: hyai, ps0
use shr_vmath_mod, only: shr_vmath_log
use shr_kind_mod, only: shr_kind_cx
use dyn_comp, only: ixo, ixo2, ixh, ixh2
+ use cam_thermo_formula,only: ENERGY_FORMULA_DYCORE_SE
! arguments
type(runtime_options), intent(in) :: cam_runtime_opts ! Runtime settings object
@@ -607,7 +611,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
!constituent properties pointer
type(ccpp_constituent_prop_ptr_t), pointer :: const_prop_ptr(:)
- integer :: m, i, k
+ integer :: m, i, k, m_cnst
integer :: ix_q
!Needed for "geopotential_temp" CCPP scheme
@@ -622,6 +626,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
real(r8) :: mmrSum_O_O2_H ! Sum of mass mixing ratios for O, O2, and H
real(r8), parameter :: mmrMin=1.e-20_r8 ! lower limit of o2, o, and h mixing ratios
real(r8), parameter :: N2mmrMin=1.e-6_r8 ! lower limit of o2, o, and h mixing ratios
+ real(r8), parameter :: H2lim=6.e-5_r8 ! H2 limiter: 10x global H2 MMR (Roble, 1995)
!----------------------------------------------------------------------------
! Nullify pointers
@@ -683,14 +688,23 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
end do
! wet pressure variables (should be removed from physics!)
+ factor_array(:,:) = 1.0_kind_phys
+ !$omp parallel do num_threads(horz_num_threads) private (k, i, m_cnst)
+ do m_cnst = dry_air_species_num + 1, thermodynamic_active_species_num
+ ! include all water species in the factor array.
+ m = thermodynamic_active_species_idx(m_cnst)
+ do k = 1, nlev
+ do i = 1, pcols
+ ! at this point all q's are dry
+ factor_array(i,k) = factor_array(i,k) + const_data_ptr(i,k,m)
+ end do
+ end do
+ end do
!$omp parallel do num_threads(horz_num_threads) private (k, i)
- do k=1,nlev
- do i=1, pcols
- ! to be consistent with total energy formula in physic's check_energy module only
- ! include water vapor in moist dp
- factor_array(i,k) = 1._kind_phys+const_data_ptr(i,k,ix_q)
- phys_state%pdel(i,k) = phys_state%pdeldry(i,k)*factor_array(i,k)
+ do k = 1, nlev
+ do i = 1, pcols
+ phys_state%pdel(i,k) = phys_state%pdeldry(i,k) * factor_array(i,k)
end do
end do
@@ -721,31 +735,12 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
do k = 1, nlev
do i = 1, pcols
phys_state%rpdel(i,k) = 1._kind_phys/phys_state%pdel(i,k)
- phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa
- end do
- end do
-
- ! all tracers (including moisture) are in dry mixing ratio units
- ! physics expect water variables moist
- factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev)
-
- !$omp parallel do num_threads(horz_num_threads) private (m, k, i)
- do m=1, num_advected
- do k = 1, nlev
- do i=1, pcols
- !This should ideally check if a constituent is a wet
- !mixing ratio or not, but until that is working properly
- !in the CCPP framework just check for the water species status
- !instead, which is all that CAM physics requires:
- if (const_is_water_species(m)) then
- const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m)
- end if
- end do
end do
end do
!------------------------------------------------------------
- ! Ensure O2 + O + H (N2) mmr greater than one.
+ ! Apply limiters to mixing ratios of major species (WACCMX):
+ ! Ensure N2 = 1 - (O2 + O + H) mmr is greater than 0
! Check for unusually large H2 values and set to lower value.
!------------------------------------------------------------
if (cam_runtime_opts%waccmx_option() == 'ionosphere' .or. &
@@ -769,8 +764,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
endif
- if(const_data_ptr(i,k,ixh2) .gt. 6.e-5_r8) then
- const_data_ptr(i,k,ixh2) = 6.e-5_r8
+ if(const_data_ptr(i,k,ixh2) > H2lim) then
+ const_data_ptr(i,k,ixh2) = H2lim
endif
end do
@@ -789,11 +784,61 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
! Compute molecular viscosity(kmvis) and conductivity(kmcnd).
! Update zvirv registry variable; calculated for WACCM-X.
!-----------------------------------------------------------------------------
+ if (dry_air_species_num > 0) then
+ call cam_thermo_dry_air_update( &
+ 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
- call cam_thermo_update(const_data_ptr, phys_state%t, pcols, &
- 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 = const_data_ptr, & ! dry MMR
+ ncol = pcols, &
+ pver = pver, &
+ energy_formula = ENERGY_FORMULA_DYCORE_SE, &
+ cp_or_cv_dycore = cp_or_cv_dycore &
+ )
- !Call geopotential_temp CCPP scheme:
+ !$omp parallel do num_threads(horz_num_threads) private (k, i)
+ do k = 1, nlev
+ do i = 1, pcols
+ phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k)
+ end do
+ end do
+
+ ! ========= Q is dry ^^^ ---- Q is moist vvv ========= !
+
+ !
+ ! CAM physics expects that: water tracers (including moisture) are moist; the rest dry mixing ratio
+ ! at this point Q is converted to moist.
+ !
+ factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev)
+
+ !$omp parallel do num_threads(horz_num_threads) private (m, k, i)
+ do m = 1, num_advected
+ do k = 1, nlev
+ do i = 1, pcols
+ ! This should ideally check if a constituent is a wet
+ ! mixing ratio or not, but until that is working properly
+ ! in the CCPP framework just check for the water species status
+ ! instead, which is all that CAM physics requires:
+ if (const_is_water_species(m)) then
+ const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m)
+ end if
+ end do
+ end do
+ end do
+
+ ! Call geopotential_temp CCPP scheme:
call geopotential_temp_run(pver, lagrangian_vertical, pver, 1, &
pverp, 1, num_advected, phys_state%lnpint, &
phys_state%pint, phys_state%pmid, phys_state%pdel, &
@@ -807,12 +852,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_timestep_init(phys_state, phys_tend, pbuf_chnk)
-#endif
-
end subroutine derived_phys_dry
!=========================================================================================
diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90
index 1341b9f4..cf88d7f9 100644
--- a/src/dynamics/se/dycore/prim_advance_mod.F90
+++ b/src/dynamics/se/dycore/prim_advance_mod.F90
@@ -1674,8 +1674,12 @@ subroutine calc_tot_energy_dynamics(elem,fvm,nets,nete,tl,tl_qdp,outfld_name_suf
call get_dp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),MASS_MIXING_RATIO,thermodynamic_active_species_idx_dycore,&
elem(ie)%state%dp3d(:,:,:,tl),pdel,ps=ps,ptop=hyai(1)*ps0)
call get_cp(elem(ie)%state%Qdp(:,:,:,1:qsize,tl_qdp),&
- .false.,cp,dp_dry=elem(ie)%state%dp3d(:,:,:,tl),&
+ .false.,cp,factor=1.0_r8/elem(ie)%state%dp3d(:,:,:,tl),&
active_species_idx_dycore=thermodynamic_active_species_idx_dycore)
+
+ ! TODO: need to port cam6_3_109 changes to total energy using get_hydrostatic_energy
+ ! https://github.com/ESCOMP/CAM/pull/761/files#diff-946bde17289e2f42e43e64413610aa11d102deda8b5199ddaa5b71e67e5d517a
+
do k = 1, nlev
do j=1,np
do i = 1, np
diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90
index ab52d91c..4f70ae2c 100644
--- a/src/dynamics/se/dyn_comp.F90
+++ b/src/dynamics/se/dyn_comp.F90
@@ -566,6 +566,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
use dyn_thermo, only: get_molecular_diff_coef_reference
!use cam_history, only: addfld, add_default, horiz_only, register_vector_field
use gravity_waves_sources, only: gws_init
+ use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_SE
+ use physics_types, only: dycore_energy_consistency_adjust
+ use phys_vars_init_check, only: mark_as_initialized
!SE dycore:
use prim_advance_mod, only: prim_advance_init
@@ -642,6 +645,14 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out)
real(r8) :: tau0, krange, otau0, scale
real(r8) :: km_sponge_factor_local(nlev+1)
!----------------------------------------------------------------------------
+ ! Set dynamical core energy formula for use in cam_thermo.
+ energy_formula_dycore = ENERGY_FORMULA_DYCORE_SE
+ call mark_as_initialized("total_energy_formula_for_dycore")
+
+ ! Dynamical core energy is not consistent with CAM physics and requires
+ ! temperature and temperature tendency adjustment at end of physics.
+ dycore_energy_consistency_adjust = .true.
+ call mark_as_initialized("flag_for_dycore_energy_consistency_adjustment")
! Now allocate and set condenstate vars
allocate(cnst_name_gll(qsize), stat=iret) ! constituent names for gll tracers
@@ -1852,6 +1863,17 @@ subroutine read_inidat(dyn_in)
call mark_as_initialized("tendency_of_air_temperature_due_to_model_physics")
call mark_as_initialized("tendency_of_eastward_wind_due_to_model_physics")
call mark_as_initialized("tendency_of_northward_wind_due_to_model_physics")
+ call mark_as_initialized("specific_heat_of_air_used_in_dycore")
+
+ ! These energy variables are calculated by check_energy_timestep_init
+ ! but need to be marked here
+ call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_energy_using_physics_energy_formula")
+ call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_energy_using_dycore_energy_formula")
+ call mark_as_initialized("vertically_integrated_total_water_at_start_of_physics_timestep")
+ call mark_as_initialized("vertically_integrated_total_water")
+ call mark_as_initialized("vertically_integrated_total_energy_at_end_of_physics_timestep")
end subroutine read_inidat
diff --git a/src/dynamics/utils/dyn_thermo.F90 b/src/dynamics/utils/dyn_thermo.F90
index c4c4723c..9b465031 100644
--- a/src/dynamics/utils/dyn_thermo.F90
+++ b/src/dynamics/utils/dyn_thermo.F90
@@ -61,7 +61,7 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore)
!Declare local variables:
real(kind_phys), allocatable :: tracer_phys(:,:,:,:)
real(kind_phys), allocatable :: cp_phys(:,:,:)
- real(kind_phys), allocatable :: dp_dry_phys(:,:,:)
+ real(kind_phys), allocatable :: factor_phys(:,:,:)
!check_allocate variables:
integer :: iret !allocate status integer
@@ -70,11 +70,16 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore)
!Check if kinds are different:
if (kind_phys == kind_dyn) then
- !The dynamics and physics kind is the same, so just call the physics
- !routine directly:
- call get_cp_phys(tracer,inv_cp,cp, &
- dp_dry=dp_dry, &
- active_species_idx_dycore=active_species_idx_dycore)
+ ! The dynamics and physics kind is the same, so just call the physics
+ ! routine directly:
+ if(present(dp_dry)) then
+ call get_cp_phys(tracer,inv_cp,cp, &
+ factor=1.0_kind_phys/dp_dry, &
+ active_species_idx_dycore=active_species_idx_dycore)
+ else
+ call get_cp_phys(tracer,inv_cp,cp, &
+ active_species_idx_dycore=active_species_idx_dycore)
+ endif
else
@@ -95,18 +100,18 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore)
!Allocate and set optional variables:
if (present(dp_dry)) then
- allocate(dp_dry_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret)
+ allocate(factor_phys(size(dp_dry, 1), size(dp_dry, 2), size(dp_dry,3)), stat=iret)
call check_allocate(iret, subname, &
- 'dp_dry_phys', &
+ 'factor_phys', &
file=__FILE__, line=__LINE__)
!Set optional local variable:
- dp_dry_phys = real(dp_dry, kind_phys)
+ factor_phys = 1.0_kind_phys/real(dp_dry, kind_phys)
end if
- !Call physics routine using local vriables with matching kinds:
+ !Call physics routine using local variables with matching kinds:
call get_cp_phys(tracer_phys,inv_cp,cp_phys, &
- dp_dry=dp_dry_phys, &
+ factor=factor_phys, &
active_species_idx_dycore=active_species_idx_dycore)
!Set output variables back to dynamics kind:
@@ -116,8 +121,8 @@ subroutine get_cp(tracer,inv_cp,cp,dp_dry,active_species_idx_dycore)
deallocate(tracer_phys)
deallocate(cp_phys)
- if (allocated(dp_dry_phys)) then
- deallocate(dp_dry_phys)
+ if (allocated(factor_phys)) then
+ deallocate(factor_phys)
end if
@@ -957,7 +962,7 @@ subroutine get_rho_dry(tracer,temp,ptop,dp_dry,tracer_mass,&
end if
- !Call physics routine using local vriables with matching kinds:
+ !Call physics routine using local variables with matching kinds:
call get_rho_dry_phys(tracer_phys,temp_phys, &
ptop_phys, dp_dry_phys,tracer_mass, &
rho_dry=rho_dry_phys, &
diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp
index 0ecfcc15..e7a599f4 160000
--- a/src/physics/ncar_ccpp
+++ b/src/physics/ncar_ccpp
@@ -1 +1 @@
-Subproject commit 0ecfcc155ac0387ef9db3304611c6f3ef055ac1d
+Subproject commit e7a599f4bb1533f7cdcd8723b1f864e11578e96c
diff --git a/src/physics/utils/phys_comp.F90 b/src/physics/utils/phys_comp.F90
index 5dbfc20a..aa1feef6 100644
--- a/src/physics/utils/phys_comp.F90
+++ b/src/physics/utils/phys_comp.F90
@@ -174,10 +174,12 @@ subroutine phys_init()
use physics_grid, only: columns_on_task
use vert_coord, only: pver, pverp
use cam_thermo, only: cam_thermo_init
+ use cam_thermo_formula, only: cam_thermo_formula_init
use physics_types, only: allocate_physics_types_fields
use cam_ccpp_cap, only: cam_ccpp_physics_initialize
call cam_thermo_init(columns_on_task, pver, pverp)
+ call cam_thermo_formula_init()
call allocate_physics_types_fields(columns_on_task, pver, pverp, &
set_init_val_in=.true., reallocate_in=.false.)
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index 39fc0f99..8dd3dca3 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -39,8 +39,10 @@ module physics_grid
public :: get_rlat_p ! latitude of a physics column in radians
public :: get_rlon_p ! longitude of a physics column in radians
public :: get_area_p ! area of a physics column in radians squared
+ public :: get_wght_p ! weight of a physics column in radians squared
public :: get_rlat_all_p ! latitudes of physics cols on task (radians)
public :: get_rlon_all_p ! longitudes of physics cols on task (radians)
+ public :: get_wght_all_p ! weights of physics cols on task
public :: get_dyn_col_p ! dynamics local blk number and blk offset(s)
public :: global_index_p ! global column index of a physics column
public :: local_index_p ! local column index of a physics column
@@ -376,8 +378,6 @@ end subroutine phys_grid_init
!========================================================================
real(r8) function get_dlat_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! latitude of a physics column in degrees
! Dummy argument
@@ -396,8 +396,6 @@ end function get_dlat_p
!========================================================================
real(r8) function get_dlon_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! longitude of a physics column in degrees
! Dummy argument
@@ -416,8 +414,6 @@ end function get_dlon_p
!========================================================================
real(r8) function get_rlat_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! latitude of a physics column in radians
! Dummy argument
@@ -436,8 +432,6 @@ end function get_rlat_p
!========================================================================
real(r8) function get_rlon_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! longitude of a physics column in radians
! Dummy argument
@@ -456,8 +450,6 @@ end function get_rlon_p
!========================================================================
real(r8) function get_area_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! area of a physics column in radians squared
! Dummy argument
@@ -475,9 +467,25 @@ end function get_area_p
!========================================================================
+ real(r8) function get_wght_p(index)
+ ! weight of a physics column in radians squared
+
+ ! Dummy argument
+ integer, intent(in) :: index
+ ! Local variables
+ character(len=128) :: errmsg
+ character(len=*), parameter :: subname = 'get_wght_p'
+
+ ! Check that input is valid:
+ call check_phys_input(subname, index)
+
+ get_wght_p = phys_columns(index)%weight
+
+ end function get_wght_p
+
+ !========================================================================
+
subroutine get_rlat_all_p(rlatdim, rlats)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
!-----------------------------------------------------------------------
!
! get_rlat_all_p: Return all latitudes (in radians) on task.
@@ -506,8 +514,6 @@ end subroutine get_rlat_all_p
!========================================================================
subroutine get_rlon_all_p(rlondim, rlons)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
!-----------------------------------------------------------------------
!
! get_rlon_all_p: Return all longitudes (in radians) on task.
@@ -535,8 +541,35 @@ end subroutine get_rlon_all_p
!========================================================================
+ subroutine get_wght_all_p(wghtdim, wghts)
+ !-----------------------------------------------------------------------
+ !
+ ! get_wght_all_p: Return all weights on task.
+ !
+ !-----------------------------------------------------------------------
+ ! Dummy Arguments
+ integer, intent(in) :: wghtdim ! declared size of output array
+ real(r8), intent(out) :: wghts(wghtdim) ! array of weights
+
+ ! Local variables
+ integer :: index ! loop index
+ character(len=128) :: errmsg
+ character(len=*), parameter :: subname = 'get_wght_all_p: '
+
+ !-----------------------------------------------------------------------
+
+ ! Check that input is valid:
+ call check_phys_input(subname, wghtdim)
+
+ do index = 1, wghtdim
+ wghts(index) = phys_columns(index)%weight
+ end do
+
+ end subroutine get_wght_all_p
+
+ !========================================================================
+
subroutine get_dyn_col_p(index, blk_num, blk_ind)
- use cam_logfile, only: iulog
use cam_abortutils, only: endrun
! Return the dynamics local block number and block offset(s) for
! the physics column indicated by .
@@ -568,8 +601,6 @@ end subroutine get_dyn_col_p
!========================================================================
integer function global_index_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! global column index of a physics column
! Dummy argument
@@ -586,8 +617,6 @@ integer function global_index_p(index)
end function global_index_p
integer function local_index_p(index)
- use cam_logfile, only: iulog
- use cam_abortutils, only: endrun
! global column index of a physics column
! Dummy argument
diff --git a/src/physics/utils/physics_grid.meta b/src/physics/utils/physics_grid.meta
index fcf934b5..fc4798d6 100644
--- a/src/physics/utils/physics_grid.meta
+++ b/src/physics/utils/physics_grid.meta
@@ -38,13 +38,13 @@
protected = True
[ lat_rad ]
standard_name = latitude
- units = radians
+ units = rad
type = real | kind = kind_phys
dimensions = (horizontal_dimension)
protected = True
[ lon_rad ]
standard_name = longitude
- units = radians
+ units = rad
type = real | kind = kind_phys
dimensions = (horizontal_dimension)
protected = True
@@ -62,7 +62,7 @@
protected = True
[ area ]
standard_name = cell_angular_area
- units = steradian
+ units = sr
type = real | kind = kind_phys
dimensions = (horizontal_dimension)
protected = True
diff --git a/src/utils/gmean_mod.F90 b/src/utils/gmean_mod.F90
new file mode 100644
index 00000000..7959ebf0
--- /dev/null
+++ b/src/utils/gmean_mod.F90
@@ -0,0 +1,263 @@
+module gmean_mod
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Perform global mean calculations for energy conservation and other checks.
+ !
+ ! Method:
+ ! Reproducible (scalable):
+ ! Convert to fixed point (integer representation) to enable
+ ! reproducibility when using MPI collectives.
+ ! If error checking is on (via setting reprosum_diffmax > 0 and
+ ! reprosum_recompute = .true. in user_nl_cpl), shr_reprosum_calc will
+ ! check the accuracy of its computation with a fast but
+ ! non-reproducible algorithm. If any error is reported, report
+ ! the difference and the expected sum and abort run (call endrun)
+ !
+ ! gmean_mod in to_be_ccppized is different from the CAM version and
+ ! has chunk support removed.
+ !
+ !
+ !-----------------------------------------------------------------------
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use physics_grid, only: pcols => columns_on_task
+ use perf_mod, only: t_startf, t_stopf
+ use cam_logfile, only: iulog
+
+ implicit none
+ private
+
+ public :: gmean ! compute global mean of 2D fields on physics decomposition
+
+ interface gmean
+ module procedure gmean_arr
+ module procedure gmean_scl
+ end interface gmean
+
+ private :: gmean_fixed_repro
+ private :: gmean_float_norepro
+
+ ! Set do_gmean_tests to .true. to run a gmean challenge test
+ logical, private :: do_gmean_tests = .false.
+
+CONTAINS
+
+ !
+ !========================================================================
+ !
+
+ subroutine gmean_arr (arr, arr_gmean, nflds)
+ use shr_strconvert_mod, only: toString
+ use spmd_utils, only: masterproc
+ use cam_abortutils, only: endrun
+ use shr_reprosum_mod, only: shr_reprosum_reldiffmax, shr_reprosum_recompute, shr_reprosum_tolExceeded
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Compute the global mean of each field in "arr" in the physics grid
+ !
+ ! Method is to call shr_reprosum_calc (called from gmean_fixed_repro)
+ !-----------------------------------------------------------------------
+ !
+ ! Arguments
+ !
+ integer, intent(in) :: nflds ! number of fields
+ real(r8), intent(in) :: arr(pcols, nflds)
+ real(r8), intent(out) :: arr_gmean(nflds) ! global means
+ !
+ ! Local workspace
+ !
+ real(r8) :: rel_diff(2, nflds)
+ integer :: ifld ! field index
+ integer :: num_err
+ logical :: write_warning
+ !
+ !-----------------------------------------------------------------------
+ !
+ call t_startf('gmean_arr')
+ call t_startf ('gmean_fixed_repro')
+ call gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds)
+ call t_stopf ('gmean_fixed_repro')
+
+ ! check that "fast" reproducible sum is accurate enough. If not, calculate
+ ! using old method
+ write_warning = masterproc
+ num_err = 0
+ if (shr_reprosum_tolExceeded('gmean', nflds, write_warning, &
+ iulog, rel_diff)) then
+ if (shr_reprosum_recompute) then
+ do ifld = 1, nflds
+ if (rel_diff(1, ifld) > shr_reprosum_reldiffmax) then
+ call gmean_float_norepro(arr(:,ifld), arr_gmean(ifld), ifld)
+ num_err = num_err + 1
+ end if
+ end do
+ end if
+ end if
+ call t_stopf('gmean_arr')
+ if (num_err > 0) then
+ call endrun('gmean: '//toString(num_err)//' reprosum errors found')
+ end if
+
+ end subroutine gmean_arr
+
+ !
+ !========================================================================
+ !
+
+ subroutine gmean_scl (arr, gmean)
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Compute the global mean of a single field in "arr" in the physics grid
+ !
+ !-----------------------------------------------------------------------
+ !
+ ! Arguments
+ !
+ real(r8), intent(in) :: arr(pcols)
+ ! Input array
+ real(r8), intent(out) :: gmean ! global means
+ !
+ ! Local workspace
+ !
+ integer, parameter :: nflds = 1
+ real(r8) :: gmean_array(nflds)
+ real(r8) :: array(pcols, nflds)
+ integer :: ncols, lchnk
+
+ array(:ncols, 1) = arr(:ncols)
+ call gmean_arr(array, gmean_array, nflds)
+ gmean = gmean_array(1)
+
+ end subroutine gmean_scl
+
+ !
+ !========================================================================
+ !
+
+ subroutine gmean_float_norepro(arr, repro_sum, index)
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Compute the global mean of in the physics chunked
+ ! decomposition using a fast but non-reproducible algorithm.
+ ! Log that value along with the value computed by
+ ! shr_reprosum_calc ()
+ !
+ !-----------------------------------------------------------------------
+
+ use physconst, only: pi
+ use spmd_utils, only: masterproc, masterprocid, mpicom
+ use mpi, only: mpi_real8, mpi_sum
+ use physics_grid, only: get_wght_p
+ !
+ ! Arguments
+ !
+ real(r8), intent(in) :: arr(pcols)
+ real(r8), intent(in) :: repro_sum ! Value computed by reprosum
+ integer, intent(in) :: index ! Index of field in original call
+ !
+ ! Local workspace
+ !
+ integer :: icol
+ integer :: ierr
+ real(r8) :: wght
+ real(r8) :: check
+ real(r8) :: check_sum
+ real(r8), parameter :: pi4 = 4.0_r8 * pi
+
+ !
+ !-----------------------------------------------------------------------
+ !
+ ! Calculate and print out non-reproducible value
+ check = 0.0_r8
+ do icol = 1, pcols
+ wght = get_wght_p(icol)
+ check = check + arr(icol) * wght
+ end do
+ call MPI_reduce(check, check_sum, 1, mpi_real8, mpi_sum, &
+ masterprocid, mpicom, ierr)
+
+ ! normalization
+ check_sum = check_sum / pi4
+
+ if (masterproc) then
+ write(iulog, '(a,i0,2(a,e20.13e2))') 'gmean(', index, ') = ', &
+ check_sum, ', reprosum reported ', repro_sum
+ end if
+
+ end subroutine gmean_float_norepro
+
+ !
+ !========================================================================
+ !
+ subroutine gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds)
+ !-----------------------------------------------------------------------
+ !
+ ! Purpose:
+ ! Compute the global mean of each field in "arr" in the physics grid
+ ! with a reproducible yet scalable implementation
+ ! based on a fixed-point algorithm.
+ !
+ !-----------------------------------------------------------------------
+ use spmd_utils, only: mpicom
+ use physics_grid, only: get_wght_all_p
+ use physics_grid, only: ngcols_p => num_global_phys_cols
+ use physconst, only: pi
+ use shr_reprosum_mod, only: shr_reprosum_calc
+ use cam_abortutils, only: check_allocate
+ !
+ ! Arguments
+ !
+ integer, intent(in) :: nflds ! number of fields
+ real(r8), intent(in) :: arr(pcols,nflds)
+ ! arr_gmean: output global sums
+ real(r8), intent(out) :: arr_gmean(nflds)
+ ! rel_diff: relative and absolute differences from shr_reprosum_calc
+ real(r8), intent(out) :: rel_diff(2, nflds)
+ !
+ ! Local workspace
+ !
+ real(r8), parameter :: pi4 = 4.0_r8 * pi
+ character(len=*), parameter :: subname = 'gmean_fixed_repro: '
+
+ integer :: icol, ifld ! column, field indices
+ integer :: errflg
+
+ real(r8) :: wght(pcols) ! integration weights
+ real(r8), allocatable :: xfld(:,:) ! weighted summands
+
+ errflg = 0
+
+ allocate(xfld(pcols, nflds), stat=errflg)
+ call check_allocate(errflg, subname, 'xfld(pcols, nflds)', &
+ file=__FILE__, line=__LINE__)
+
+ ! pre-weight summands
+ call get_wght_all_p(pcols, wght)
+
+ do ifld = 1, nflds
+ do icol = 1, pcols
+ xfld(icol, ifld) = arr(icol, ifld) * wght(icol)
+ end do
+ end do
+
+ ! call fixed-point algorithm
+ call shr_reprosum_calc ( &
+ arr = xfld, &
+ arr_gsum = arr_gmean, &
+ nsummands = pcols, & ! # of local summands
+ dsummands = pcols, & ! declared first dimension of arr.
+ nflds = nflds, &
+ commid = mpicom, &
+ rel_diff = rel_diff &
+ )
+
+ deallocate(xfld)
+ ! final normalization
+ arr_gmean(:) = arr_gmean(:) / pi4
+
+ end subroutine gmean_fixed_repro
+
+end module gmean_mod
diff --git a/test/existing-test-failures.txt b/test/existing-test-failures.txt
index 253b5633..05e0a145 100644
--- a/test/existing-test-failures.txt
+++ b/test/existing-test-failures.txt
@@ -2,12 +2,7 @@ SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_intel.cam-outfrq_kessler_mpas_derecho
SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_gnu.cam-outfrq_kessler_mpas_derecho (Overall: FAIL)
- will fail until MPAS is fully integrated
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FTJ16.derecho_intel.cam-outfrq_se_cslam (Overall: FAIL)
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FKESSLER.derecho_intel.cam-outfrq_se_cslam (Overall: FAIL)
-SMS_Ln2.ne3pg3_ne3pg3_mg37.FPHYStest.derecho_intel.cam-outfrq_kessler_derecho (Overall: FAIL)
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_intel.cam-outfrq_se_cslam_analy_ic (Overall: FAIL)
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FTJ16.derecho_gnu.cam-outfrq_se_cslam (Overall: FAIL)
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FKESSLER.derecho_gnu.cam-outfrq_se_cslam (Overall: FAIL)
-SMS_Ln2.ne3pg3_ne3pg3_mg37.FPHYStest.derecho_gnu.cam-outfrq_kessler_derecho (Overall: FAIL)
-SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_gnu.cam-outfrq_se_cslam_analy_ic (Overall: FAIL)
- - will fail until https://github.com/ESCOMP/CAM-SIMA/pull/316 is merged
+SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_intel.cam-outfrq_se_cslam_analy_ic (Overall: PEND) details:
+SMS_Ln9.ne5pg3_ne5pg3_mg37.FCAM7.derecho_gnu.cam-outfrq_se_cslam_analy_ic (Overall: PEND) details:
+ - build failure due to dadadj_apply_qv_tendency removal, needs to be updated to use constituent tendency updater atmospheric_physics#179
+ - also expected runtime failure CAM-SIMA#335
diff --git a/test/unit/sample_files/phys_types_dup_section.meta b/test/unit/sample_files/phys_types_dup_section.meta
index 65127e8e..23db9b12 100644
--- a/test/unit/sample_files/phys_types_dup_section.meta
+++ b/test/unit/sample_files/phys_types_dup_section.meta
@@ -13,13 +13,13 @@
protected = True
[ latitude ]
standard_name = latitude
- units = radians
+ units = rad
type = real | kind = kind_phys
dimensions = (horizontal_dimension)
protected = True
[ longitude ]
standard_name = longitude
- units = radians
+ units = rad
type = real | kind = kind_phys
dimensions = (horizontal_dimension)
protected = True
diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml
index 0ac72a21..dbd84b99 100644
--- a/tools/stdnames_to_inputnames_dictionary.xml
+++ b/tools/stdnames_to_inputnames_dictionary.xml
@@ -160,42 +160,47 @@
state_zi
-
+
te_ini_phys
state_te_ini_phys
-
+
te_cur_phys
state_te_cur_phys
-
+
te_ini_dyn
state_te_ini_dyn
-
+
te_cur_dyn
state_te_cur_dyn
-
+
tw_ini
state_tw_ini
-
+
tw_cur
state_tw_cur
+
+
+ cp_or_cv_dycore
+
+
RHO