Skip to content

Commit

Permalink
Fix wet-to-dry conversion and Held-Suarez ICs.
Browse files Browse the repository at this point in the history
  • Loading branch information
@nusbaume
nusbaume committed Sep 6, 2024
1 parent 63ca8be commit 5c54cb0
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 44 deletions.
2 changes: 1 addition & 1 deletion .gitmodules
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
[submodule "history"]
path = src/history/buffers
url = https://github.com/peverwhee/history_output
fxtag = sima-history
fxtag = 7cd5a736016be3a970601e65c9d9dbefa31ff9e9
fxrequired = AlwaysRequired
fxDONOTUSEurl = https://github.com/peverwhee/history_output
[submodule "mpas"]
Expand Down
63 changes: 43 additions & 20 deletions src/dynamics/se/dp_coupling.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,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 @@ -389,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 @@ -551,11 +573,14 @@ 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: const_get_index, const_qmin
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
Expand Down Expand Up @@ -583,7 +608,7 @@ 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
Expand All @@ -604,11 +629,7 @@ 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()
Expand Down Expand Up @@ -667,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 @@ -708,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
84 changes: 61 additions & 23 deletions src/dynamics/tests/initial_conditions/ic_held_suarez.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,12 @@ module ic_held_suarez

subroutine hs94_set_ic(latvals, lonvals, U, V, T, PS, PHIS, &
Q, m_cnst, mask, verbose)
!use const_init, only: cnst_init_default
use cam_constituents, only: const_get_index
use shr_kind_mod, only: cx => shr_kind_cx
use ccpp_kinds, only: kind_phys
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use cam_ccpp_cap, only: cam_model_const_properties
use cam_constituents, only: const_get_index, const_qmin
use runtime_obj, only: wv_stdname

!-----------------------------------------------------------------------
!
Expand All @@ -49,14 +53,21 @@ subroutine hs94_set_ic(latvals, lonvals, U, V, T, PS, PHIS, &
! Local variables
logical, allocatable :: mask_use(:)
logical :: verbose_use
logical :: const_has_default
integer :: i, k, m
integer :: ix_cld_liq, ix_rain
integer :: ix_q
integer :: ncol
integer :: nlev
integer :: ncnst
integer :: iret
real(kind_phys) :: const_default_value !Constituent default value
real(kind_phys) :: const_qmin_value !Constituent minimum value
character(len=cx) :: errmsg !CCPP error message
character(len=*), parameter :: subname = 'HS94_SET_IC'

!Private array of constituent properties (for property interface functions)
type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:)

allocate(mask_use(size(latvals)), stat=iret)
call check_allocate(iret, subname, 'mask_use(size(latvals))', &
file=__FILE__, line=__LINE__)
Expand All @@ -76,12 +87,6 @@ subroutine hs94_set_ic(latvals, lonvals, U, V, T, PS, PHIS, &
verbose_use = .true.
end if

!set constituent indices
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.)

ncol = size(latvals, 1)
nlev = -1
if (present(U)) then
Expand Down Expand Up @@ -137,32 +142,65 @@ subroutine hs94_set_ic(latvals, lonvals, U, V, T, PS, PHIS, &
end if

if (present(Q)) then
!Get water vapor constituent index:
call const_get_index(wv_stdname, ix_q)

!Get constituent properties object:
const_props => cam_model_const_properties()

!Determine array sizes:
nlev = size(Q, 2)
ncnst = size(m_cnst, 1)

!Loop over all constituents:
do m = 1, ncnst
if (m_cnst(m) == 1) then
if (m_cnst(m) == ix_q) then
! No water vapor in Held-Suarez
do k = 1, nlev
where(mask_use)
Q(:,k,m_cnst(m)) = 0.0_r8
Q(:,k,m) = 0.0_r8
end where
end do
!Un-comment once constituents are working in CAMDEN -JN:
#if 0

if(masterproc .and. verbose_use) then
write(iulog,*) ' ', trim(cnst_name(m_cnst(m))), ' initialized by "',subname,'"'
write(iulog,*) ' ', wv_stdname, ' initialized by "',subname,'"'
end if

else
call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),&
mask=mask_use, verbose=verbose_use, notfound=.false.)
#else
else
!Initialize cloud liquid and rain until constituent routines are enabled:
Q(:,:,ix_cld_liq) = 0.0_r8
Q(:,:,ix_rain) = 0.0_r8
#endif
end if
!Extract constituent minimum value:
const_qmin_value = const_qmin(m_cnst(m))

!Initialize constituent to its minimum value:
Q(:,:,m) = real(const_qmin_value, r8)

!Check for default value in constituent properties object:
call const_props(m_cnst(m))%has_default(const_has_default, &
iret, &
errmsg)
if (iret /= 0) then
call endrun(errmsg, file=__FILE__, line=__LINE__)
end if

if (const_has_default) then

!If default value exists, then extract default value
!from constituent properties object:
call const_props(m_cnst(m))%default_value(const_default_value, &
iret, &
errmsg)
if (iret /= 0) then
call endrun(errmsg, file=__FILE__, line=__LINE__)
end if

!Set constituent to default value in masked region:
do k=1,nlev
where(mask_use)
Q(:,k,m) = real(const_default_value, r8)
end where
end do

end if !has_default
end if !water vapor
end do
end if

Expand Down

0 comments on commit 5c54cb0

Please sign in to comment.