Skip to content

Commit

Permalink
Implement constituents infrastructure in analytic IC routines (ESCOMP…
Browse files Browse the repository at this point in the history 
…#299)

Originator(s): nusbaume, adamrher

Summary:

Adds in the subroutine/function hooks needed for the analytic ICs to
properly interact with the CCPP constituents object. This PR also
contains some bug fixes for the air composition and SE dycore
dynamics-physics coupling routines that were found via testing with the
FKESSLER and FHS94 compsets. Finally, this PR also includes some slight
cleanup of the use of the water vapor constituent standard name within
the core CAM-SIMA host code.

Fixes ESCOMP#287

cod reviewed by:  adamrher, peverwhee

Describe any changes made to build system:  N/A

Describe any changes made to the namelist:

M       src/dynamics/tests/namelist_definition_analy_ic.xml
    - Added the US standard atmosphere option to the analytic IC types.

List any changes to the defaults for the input datasets (e.g. boundary
datasets): N/A

List all files eliminated and why: N/A

List all files added and what they do: N/A

List all existing files that have been modified, and describe the
changes:

M       src/control/cam_comp.F90
M       src/control/runtime_obj.F90
- Added new "wv_stdname" parameter to cleanup water vapor standard name
usage.

M       src/data/air_composition.F90
- Fix cpair bug when no dry air species are listed, and use new
"wv_stdname" parameter.

M       src/dynamics/se/dp_coupling.F90
- Implement qneg and remaining constituents infrastructure, and add
missing wet-to-dry conversion step.

M       src/dynamics/se/dycore/prim_state_mod.F90
    -  Whitespace cleanup.

M       src/dynamics/se/dyn_comp.F90
M       src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90
M       src/dynamics/tests/initial_conditions/ic_baroclinic.F90
M       src/dynamics/tests/initial_conditions/ic_held_suarez.F90
M       src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90
    - Implement new constituents infrastructure in analytic IC routines.

If there are new failures (compare to the existing-test-failures.txt
file),
have them OK'd by the gatekeeper, note them here, and add them to the
file.
If there are baseline differences, include the test and the reason for
the
diff. What is the nature of the change? Roundoff?

derecho/intel/aux_sima:  All Pass

derecho/gnu/aux_sima:  All Pass

CAM-SIMA date used for the baseline comparison tests if different than
latest:
  • Loading branch information
@nusbaume
nusbaume authored Sep 13, 2024
1 parent 6b643f6 commit c99c034
Show file tree
Hide file tree
Showing 11 changed files with 389 additions and 182 deletions.
10 changes: 4 additions & 6 deletions src/control/cam_comp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -542,6 +542,7 @@ subroutine cam_register_constituents(cam_runtime_opts)
! physics suite being invoked during this run.
use cam_abortutils, only: endrun, check_allocate
use runtime_obj, only: runtime_options
use runtime_obj, only: wv_stdname
use phys_comp, only: phys_suite_name
use cam_constituents, only: cam_constituents_init
use cam_constituents, only: const_set_qmin, const_get_index
Expand Down Expand Up @@ -569,9 +570,7 @@ subroutine cam_register_constituents(cam_runtime_opts)

! Check if water vapor is already marked as a constituent by the
! physics:
call cam_ccpp_is_scheme_constituent( &
"water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water", &
is_constituent, errflg, errmsg)
call cam_ccpp_is_scheme_constituent(wv_stdname, is_constituent, errflg, errmsg)

if (errflg /= 0) then
call endrun(subname//trim(errmsg), file=__FILE__, line=__LINE__)
Expand All @@ -589,7 +588,7 @@ subroutine cam_register_constituents(cam_runtime_opts)

! Register the constituents so they can be advected:
call host_constituents(1)%instantiate( &
std_name="water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water", &
std_name=wv_stdname, &
long_name="water vapor mixing ratio w.r.t moist air and condensed_water", &
units="kg kg-1", &
default_value=0._kind_phys, &
Expand Down Expand Up @@ -639,8 +638,7 @@ subroutine cam_register_constituents(cam_runtime_opts)
if (phys_suite_name /= 'held_suarez_1994') then !Held-Suarez is "dry" physics

! Get constituent index for water vapor:
call const_get_index("water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water", &
const_idx)
call const_get_index(wv_stdname, const_idx)

! Set new minimum value:
call const_set_qmin(const_idx, 1.E-12_kind_phys)
Expand Down
3 changes: 3 additions & 0 deletions src/control/runtime_obj.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@ module runtime_obj
integer, public, parameter :: unset_int = huge(1)
real(r8), public, parameter :: unset_real = huge(1.0_r8)

! Water vapor constituent standard name
character(len=*), public, parameter :: wv_stdname = 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water'

! Public interfaces and data

!> \section arg_table_runtime_options Argument Table
Expand Down
16 changes: 13 additions & 3 deletions src/data/air_composition.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,10 @@ subroutine air_composition_init()
use cam_logfile, only: iulog
use physconst, only: r_universal, cpwv
use physconst, only: rh2o, cpliq, cpice
use physconst, only: cpair, rair
use physics_grid, only: pcols => columns_on_task
use vert_coord, only: pver
use runtime_obj, only: wv_stdname
use cam_constituents, only: const_name, num_advected
use cam_constituents, only: const_set_thermo_active
use cam_constituents, only: const_set_water_species
Expand Down Expand Up @@ -305,9 +307,8 @@ subroutine air_composition_init()
!
! Q
!
case('water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water')
call air_species_info('water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', &
ix, mw)
case(wv_stdname) !water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water
call air_species_info(wv_stdname, ix, mw)
thermodynamic_active_species_idx(icnst) = ix
thermodynamic_active_species_cp (icnst) = cpwv
thermodynamic_active_species_cv (icnst) = cv3 / mw
Expand Down Expand Up @@ -455,6 +456,15 @@ subroutine air_composition_init()
has_ice = .false.
end do

!Set dry air thermodynamic properities if no dry air species provided:
if (dry_species_num == 0) then
!Note: The zeroeth index is used to represent all of dry
! air instead of just N2 in this configuration
thermodynamic_active_species_cp(0) = cpair
thermodynamic_active_species_cv(0) = cpair - rair
thermodynamic_active_species_R(0) = rair
end if

water_species_in_air_num = water_species_num
dry_air_species_num = dry_species_num
thermodynamic_active_species_num = water_species_num + dry_species_num
Expand Down
96 changes: 67 additions & 29 deletions src/dynamics/se/dp_coupling.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ module dp_coupling

real(r8), allocatable :: q_prev(:,:,:) ! Previous Q for computing tendencies

real(kind_phys), allocatable :: qmin_vals(:) !Consitutent minimum values array

!=========================================================================================
CONTAINS
!=========================================================================================
Expand Down Expand Up @@ -320,6 +322,8 @@ subroutine p_d_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_in, tl_f, t
use test_fvm_mapping, only: test_mapping_overwrite_tendencies
use test_fvm_mapping, only: test_mapping_output_mapped_tendencies
use cam_ccpp_cap, only: cam_constituents_array
use cam_constituents, only: num_advected
use cam_constituents, only: const_is_water_species

! SE dycore:
use bndry_mod, only: bndry_exchange
Expand Down Expand Up @@ -387,9 +391,29 @@ subroutine p_d_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_in, tl_f, t
uv_tmp = 0.0_r8
dq_tmp = 0.0_r8

! Grab pointer to constituent array
!Grab pointer to constituent array
const_data_ptr => cam_constituents_array()

!Convert wet mixing ratios to dry, which for CAM
!configurations is only the water species:
!$omp parallel do num_threads(max_num_threads) private (k, i, m)
do ilyr = 1, nlev
do icol=1, pcols
!Determine wet to dry adjustment factor:
factor = phys_state%pdel(icol,ilyr)/phys_state%pdeldry(icol,ilyr)

!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 configurations require:
do m=1, num_advected
if (const_is_water_species(m)) then
const_data_ptr(icol,ilyr,m) = factor*const_data_ptr(icol,ilyr,m)
end if
end do
end do
end do

if (.not. allocated(q_prev)) then
call endrun('p_d_coupling: q_prev not allocated')
end if
Expand Down Expand Up @@ -549,19 +573,22 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
! mixing ratios are converted to a wet basis. Initialize geopotential heights.
! Finally compute energy and water column integrals of the physics input state.

! use constituents, only: qmin
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use cam_ccpp_cap, only: cam_constituents_array
use cam_ccpp_cap, only: cam_model_const_properties
use cam_constituents, only: num_advected
use cam_constituents, only: const_is_water_species
use cam_constituents, only: const_get_index
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 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 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
Expand All @@ -581,12 +608,14 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
type(ccpp_constituent_prop_ptr_t), pointer :: const_prop_ptr(:)

integer :: m, i, k
integer :: ix_q, ix_cld_liq, ix_rain
integer :: ix_q

!Needed for "geopotential_temp" CCPP scheme
integer :: errflg
character(len=shr_kind_cx) :: errmsg

character(len=*), parameter :: subname = 'derived_phys_dry'

!--------------------------------------------
! Variables needed for WACCM-X
!--------------------------------------------
Expand All @@ -600,16 +629,26 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
nullify(const_prop_ptr)

! Set constituent indices
call const_get_index('water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', ix_q)
call const_get_index('cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', &
ix_cld_liq, abort=.false.)
call const_get_index('rain_mixing_ratio_wrt_moist_air_and_condensed_water', &
ix_rain, abort=.false.)
call const_get_index(wv_stdname, ix_q)

! Grab pointer to constituent and properties arrays
const_data_ptr => cam_constituents_array()
const_prop_ptr => cam_model_const_properties()

! Create qmin array (if not already generated):
if (.not.allocated(qmin_vals)) then
allocate(qmin_vals(size(const_prop_ptr)), stat=errflg)
call check_allocate(errflg, subname, &
'qmin_vals(size(cam_model_const_properties))', &
file=__FILE__, line=__LINE__)


! Set relevant minimum values for each constituent:
do m = 1, size(qmin_vals)
qmin_vals(m) = const_qmin(m)
end do
end if

! Evaluate derived quantities

! dry pressure variables
Expand Down Expand Up @@ -649,7 +688,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
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 in moist dp
! 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)
end do
Expand Down Expand Up @@ -690,16 +729,18 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
! 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 (k, i)
do k = 1, nlev
do i=1, pcols
const_data_ptr(i,k,ix_q) = factor_array(i,k)*const_data_ptr(i,k,ix_q)
if (ix_cld_liq /= -1) then
const_data_ptr(i,k,ix_cld_liq) = factor_array(i,k)*const_data_ptr(i,k,ix_cld_liq)
end if
if (ix_rain /= -1) then
const_data_ptr(i,k,ix_rain) = factor_array(i,k)*const_data_ptr(i,k,ix_rain)
end if
!$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

Expand Down Expand Up @@ -736,13 +777,19 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
end do
endif

! Ensure tracers are all greater than or equal to their
! minimum-allowed value:
call qneg_run('D_P_COUPLING', columns_on_task, pver, &
qmin_vals, const_data_ptr, errflg, errmsg)

!-----------------------------------------------------------------------------
! Call cam_thermo_update. If cam_runtime_opts%update_thermodynamic_variables()
! returns .true., cam_thermo_update will compute cpairv, rairv, mbarv, and cappav as
! constituent dependent variables. It will also:
! Compute molecular viscosity(kmvis) and conductivity(kmcnd).
! Update zvirv registry variable; calculated for WACCM-X.
!-----------------------------------------------------------------------------

call cam_thermo_update(const_data_ptr, phys_state%t, pcols, &
cam_runtime_opts%update_thermodynamic_variables())

Expand All @@ -760,15 +807,6 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
phys_state%phis, phys_state%dse, cpairv, &
errflg, errmsg)

! Ensure tracers are all positive
! Please note this cannot be used until the 'qmin'
! array is publicly provided by the CCPP constituent object. -JN
#if 0
call qneg_run('D_P_COUPLING', columns_on_task, pver, &
qmin, const_data_ptr, &
errflg, errmsg)
#endif

!Remove once check_energy scheme exists in CAMDEN:
#if 0
! Compute energy and water integrals of input state
Expand Down
12 changes: 6 additions & 6 deletions src/dynamics/se/dycore/prim_state_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -418,11 +418,11 @@ subroutine adjust_nsplit(elem, tl,hybrid,nets,nete, fvm, omega_cn)
use cam_abortutils, only: endrun
use control_mod, only: nu_top
!
type (element_t), intent(inout) :: elem(:)
type (element_t), intent(inout) :: elem(:)
type (TimeLevel_t), target, intent(in) :: tl
type (hybrid_t), intent(in) :: hybrid
integer, intent(in) :: nets,nete
type(fvm_struct), intent(inout) :: fvm(:)
type(fvm_struct), intent(inout) :: fvm(:)
real (kind=r8), intent(in) :: omega_cn(2,nets:nete)
! Local variables...
integer :: k,ie
Expand Down Expand Up @@ -457,7 +457,7 @@ subroutine adjust_nsplit(elem, tl,hybrid,nets,nete, fvm, omega_cn)
nsplit=2*nsplit_baseline
fvm_supercycling = rsplit
fvm_supercycling_jet = rsplit
nu_top=2.0_r8*nu_top
nu_top=2.0_r8*nu_top
!
! write diagnostics to log file
!
Expand All @@ -470,7 +470,7 @@ subroutine adjust_nsplit(elem, tl,hybrid,nets,nete, fvm, omega_cn)
end if
dtime = get_step_size()
tstep = dtime / real(nsplit*qsplit*rsplit, r8)

else if (nsplit.ne.nsplit_baseline.and.max_o(1)<0.4_r8*threshold) then
!
! should nsplit be reduced again?
Expand All @@ -480,9 +480,9 @@ subroutine adjust_nsplit(elem, tl,hybrid,nets,nete, fvm, omega_cn)
fvm_supercycling = rsplit
fvm_supercycling_jet = rsplit
nu_top=nu_top/2.0_r8

! nu_div_scale_top(:) = 1.0_r8

dtime = get_step_size()
tstep = dtime / real(nsplit*qsplit*rsplit, r8)
if(hybrid%masterthread) then
Expand Down
11 changes: 8 additions & 3 deletions src/dynamics/se/dyn_comp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1307,7 +1307,7 @@ subroutine read_inidat(dyn_in)
file=__FILE__, line=__LINE__)

do m_cnst = 1, qsize
m_ind(m_cnst) = m_cnst
m_ind(m_cnst) = thermodynamic_active_species_idx(m_cnst)
end do

! Init tracers on the GLL grid. Note that analytic_ic_set_ic makes
Expand All @@ -1318,9 +1318,11 @@ subroutine read_inidat(dyn_in)
V=dbuf4(:,:,:,(qsize+3)), T=dbuf4(:,:,:,(qsize+4)), &
Q=dbuf4(:,:,:,1:qsize), m_cnst=m_ind, mask=pmask(:), &
PHIS_IN=PHIS_tmp)
deallocate(m_ind)

! Deallocate variables that are no longer used:
deallocate(glob_ind)
deallocate(phis_tmp)

do ie = 1, nelemd
indx = 1
do j = 1, np
Expand Down Expand Up @@ -1348,13 +1350,16 @@ subroutine read_inidat(dyn_in)
do i = 1, np
! Set qtmp at the unique columns only
if (pmask(((ie - 1) * npsq) + indx)) then
qtmp(i,j,:,ie,m_cnst) = dbuf4(indx, :, ie, m_cnst)
qtmp(i,j,:,ie,m_ind(m_cnst)) = dbuf4(indx, :, ie, m_cnst)
end if
indx = indx + 1
end do
end do
end do
end do

! Deallocate variables that are not longer used:
deallocate(m_ind)
deallocate(dbuf4)

else
Expand Down
Loading

0 comments on commit c99c034

Please sign in to comment.