Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Implement constituents infrastructure in analytic IC routines #299

Merged
merged 6 commits into from
Sep 13, 2024
Merged

Implement constituents infrastructure in analytic IC routines #299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
28d6d7b
Implement constituents in FKESSLER analytic ICs, and fix cpair bug in…
nusbaume Sep 4, 2024
fa0b5e1
Properly implement constituents in Held-Suarez ICs, and fix wet-to-dr…
nusbaume Sep 6, 2024
0d7919c
Add constituent hooks to remaining analytic IC options.
nusbaume Sep 10, 2024
894c353
Remove code added by accident during testing.
nusbaume Sep 10, 2024
e39a48a
Fix comment.
nusbaume Sep 10, 2024
a28e7be
Implement code review suggestions.
nusbaume Sep 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 4 additions & 6 deletions src/control/cam_comp.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading