From 5c54cb0fa1feb993bd6040343a412d59de337ccf Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 6 Sep 2024 14:04:06 -0600 Subject: [PATCH] Fix wet-to-dry conversion and Held-Suarez ICs. --- .gitmodules | 2 +- src/dynamics/se/dp_coupling.F90 | 63 +++++++++----- .../initial_conditions/ic_held_suarez.F90 | 84 ++++++++++++++----- 3 files changed, 105 insertions(+), 44 deletions(-) diff --git a/.gitmodules b/.gitmodules index c0b72575..17da03b7 100644 --- a/.gitmodules +++ b/.gitmodules @@ -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"] diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index d4056a5d..8b56e9d9 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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() @@ -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 @@ -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 diff --git a/src/dynamics/tests/initial_conditions/ic_held_suarez.F90 b/src/dynamics/tests/initial_conditions/ic_held_suarez.F90 index 18a5d06d..75dd0c1a 100644 --- a/src/dynamics/tests/initial_conditions/ic_held_suarez.F90 +++ b/src/dynamics/tests/initial_conditions/ic_held_suarez.F90 @@ -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 !----------------------------------------------------------------------- ! @@ -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__) @@ -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 @@ -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