From c99c034afdf0bddb46f4e38cc5ef86af3ae1bc32 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 13 Sep 2024 15:35:24 -0600 Subject: [PATCH] Implement constituents infrastructure in analytic IC routines (#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 #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: --- src/control/cam_comp.F90 | 10 +- src/control/runtime_obj.F90 | 3 + src/data/air_composition.F90 | 16 ++- src/dynamics/se/dp_coupling.F90 | 96 ++++++++++---- src/dynamics/se/dycore/prim_state_mod.F90 | 12 +- src/dynamics/se/dyn_comp.F90 | 11 +- .../initial_conditions/ic_baro_dry_jw06.F90 | 107 ++++++++++----- .../initial_conditions/ic_baroclinic.F90 | 124 +++++++++++------- .../initial_conditions/ic_held_suarez.F90 | 84 ++++++++---- .../initial_conditions/ic_us_standard_atm.F90 | 104 ++++++++++----- .../tests/namelist_definition_analy_ic.xml | 4 +- 11 files changed, 389 insertions(+), 182 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index b2fba9a2..0782738a 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -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 @@ -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__) @@ -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, & @@ -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) diff --git a/src/control/runtime_obj.F90 b/src/control/runtime_obj.F90 index 0c84c3f5..157545a8 100644 --- a/src/control/runtime_obj.F90 +++ b/src/control/runtime_obj.F90 @@ -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 diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index ea653020..f202ae97 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -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 @@ -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 @@ -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 diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 290b97a7..8b56e9d9 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -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 !========================================================================================= @@ -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 @@ -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 @@ -549,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: 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 @@ -561,7 +588,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 @@ -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 !-------------------------------------------- @@ -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 @@ -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 @@ -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 @@ -736,6 +777,11 @@ 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 @@ -743,6 +789,7 @@ 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. !----------------------------------------------------------------------------- + call cam_thermo_update(const_data_ptr, phys_state%t, pcols, & cam_runtime_opts%update_thermodynamic_variables()) @@ -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 diff --git a/src/dynamics/se/dycore/prim_state_mod.F90 b/src/dynamics/se/dycore/prim_state_mod.F90 index acd746d9..6395c169 100644 --- a/src/dynamics/se/dycore/prim_state_mod.F90 +++ b/src/dynamics/se/dycore/prim_state_mod.F90 @@ -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 @@ -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 ! @@ -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? @@ -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 diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index f04a3b26..8d74b66e 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -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 @@ -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 @@ -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 diff --git a/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 b/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 index bab2239b..87bdf868 100644 --- a/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 +++ b/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 @@ -50,10 +50,13 @@ module ic_baro_dry_jw06 subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & Q, m_cnst, mask, verbose) - use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height - use cam_constituents, only: const_get_index - !use constituents, only: cnst_name - !use const_init, only: cnst_init_default + use shr_kind_mod, only: cx => shr_kind_cx + use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height + use runtime_obj, only: wv_stdname + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use cam_ccpp_cap, only: cam_model_const_properties + use ccpp_kinds, only: kind_phys + use cam_constituents, only: const_get_index, const_qmin !----------------------------------------------------------------------- ! @@ -79,13 +82,15 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & logical, allocatable :: mask_use(:) logical :: verbose_use logical :: lu,lv,lt,lq,l3d_vars + logical :: const_has_default integer :: i, k, m integer :: ncol integer :: nlev integer :: ncnst integer :: iret - integer :: ix_rain, ix_cld_liq + integer :: ix_q, m_cnst_ix_q character(len=*), parameter :: subname = 'BC_DRY_JW06_SET_IC' + character(len=cx) :: errmsg !CCPP error message real(r8) :: tmp real(r8) :: r(size(latvals)) real(r8) :: eta @@ -93,9 +98,15 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & real(r8) :: perturb_lon, perturb_lat real(r8) :: phi_vertical real(r8) :: u_wind(size(latvals)) + real(kind_phys) :: const_default_value !Constituent default value + real(kind_phys) :: const_qmin_value !Constituent minimum value - a_omega = rearth*omega - exponent = rair*gamma/gravit + !Private array of constituent properties (for property interface functions) + type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:) + + !Set local constants: + a_omega = rearth*omega + exponent = rair*gamma/gravit allocate(mask_use(size(latvals)), stat=iret) call check_allocate(iret, subname, 'mask_use(size(latvals))', & @@ -116,10 +127,6 @@ subroutine bc_dry_jw06_set_ic(vcoord, 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) - call const_get_index('rain_mixing_ratio_wrt_moist_air_and_condensed_water', ix_rain) - ncol = size(latvals, 1) nlev = -1 @@ -236,40 +243,72 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & end if end if if (lq) then + !Get water vapor constituent index: + call const_get_index(wv_stdname, ix_q) + + !Determine which "Q" variable index matches water vapor: + m_cnst_ix_q = findloc(m_cnst, ix_q, dim=1) + do k = 1, nlev where(mask_use) - Q(:,k,1) = 0.0_r8 + Q(:,k,m_cnst_ix_q) = 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(1))), ' initialized by "',subname,'"' + write(iulog,*) ' ', wv_stdname, ' initialized by "',subname,'"' end if -#endif end if end if -!Un-comment once constituents are working in CAMDEN -JN: -#if 0 - if (lq) then - ncnst = size(m_cnst, 1) - if ((vcoord == vc_moist_pressure) .or. (vcoord == vc_dry_pressure)) then - do m = 2, ncnst - call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),& - mask=mask_use, verbose=verbose_use, notfound=.false.) - end do - end if - end if -#else if (lq) then + !Determine total number of constituents: + ncnst = size(m_cnst) + + !Extract constituent properties from CCPP constituents object: + const_props => cam_model_const_properties() + if ((vcoord == vc_moist_pressure) .or. (vcoord == vc_dry_pressure)) then - !Initialize cloud liquid and rain until constituent routines are enabled: - Q(:,:,ix_cld_liq) = 0.0_r8 - Q(:,:,ix_rain) = 0.0_r8 - end if - end if -#endif + do m = 1, ncnst + + !Skip water vapor, as it was aleady set above: + if (m == m_cnst_ix_q) cycle + + !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 do !m_cnst + end if !vcoord + end if !lq deallocate(mask_use) diff --git a/src/dynamics/tests/initial_conditions/ic_baroclinic.F90 b/src/dynamics/tests/initial_conditions/ic_baroclinic.F90 index 082aeca1..e2f76c06 100644 --- a/src/dynamics/tests/initial_conditions/ic_baroclinic.F90 +++ b/src/dynamics/tests/initial_conditions/ic_baroclinic.F90 @@ -11,9 +11,8 @@ module ic_baroclinic use cam_abortutils, only: endrun use spmd_utils, only: masterproc - use physconst, only : rair, gravit, rearth, pi, omega, epsilo - use hycoef, only : hyai, hybi, hyam, hybm, ps0 - use cam_constituents, only: const_get_index + use physconst, only: rair, gravit, rearth, pi, omega, epsilo + use hycoef, only: hyai, hybi, hyam, hybm, ps0 implicit none private @@ -76,11 +75,15 @@ module ic_baroclinic subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & Q, m_cnst, mask, verbose) - use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height - !use constituents, only: cnst_name - !use const_init, only: cnst_init_default - use inic_analytic_utils, only: analytic_ic_is_moist - use cam_abortutils, only: check_allocate + use shr_kind_mod, only: cx => shr_kind_cx + use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height + use runtime_obj, only: wv_stdname + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use cam_ccpp_cap, only: cam_model_const_properties + use ccpp_kinds, only: kind_phys + use cam_constituents, only: const_get_index, const_qmin + use inic_analytic_utils, only: analytic_ic_is_moist + use cam_abortutils, only: check_allocate !----------------------------------------------------------------------- ! @@ -107,22 +110,27 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & ! Local variables logical, allocatable :: mask_use(:) logical :: verbose_use - integer :: ix_cld_liq, ix_rain + integer :: ix_q, m_cnst_ix_q integer :: i, k, m integer :: ncol integer :: nlev integer :: ncnst integer :: iret character(len=*), parameter :: subname = 'BC_WAV_SET_IC' + character(len=cx) :: errmsg !CCPP error message real(r8) :: ztop,ptop real(r8) :: uk,vk,Tvk,qk,pk !mid-level state real(r8) :: psurface real(r8) :: wvp,qdry - logical :: lU, lV, lT, lQ, l3d_vars - logical :: cnst1_is_moisture + real(kind_phys) :: const_default_value !Constituent default value + real(kind_phys) :: const_qmin_value !Constituent minimum value + logical :: lU, lV, lT, lQ, l3d_vars, const_has_default real(r8), allocatable :: pdry_half(:), pwet_half(:),zdry_half(:),zk(:) real(r8), allocatable :: zmid(:,:) ! layer midpoint heights for test tracer initialization + !Private array of constituent properties (for property interface functions) + type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:) + if ((vcoord == vc_moist_pressure) .or. (vcoord == vc_dry_pressure)) then ! ! pressure-based vertical coordinate @@ -159,9 +167,6 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & mask_use = .true. end if - call const_get_index('cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', ix_cld_liq) - call const_get_index('rain_mixing_ratio_wrt_moist_air_and_condensed_water', ix_rain) - if (present(verbose)) then verbose_use = verbose else @@ -235,13 +240,17 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & if (lt) nlev = size(T, 2) if (lq) then + !Get water vapor constituent index: + call const_get_index(wv_stdname, ix_q) + + !Determine which "Q" variable index matches water vapor: + m_cnst_ix_q = findloc(m_cnst, ix_q, dim=1) + + !Determine vertical levels for constituents: nlev = size(Q, 2) - ! check whether first constituent in Q is water vapor. - cnst1_is_moisture = m_cnst(1) == 1 allocate(zmid(size(Q, 1),nlev), stat=iret) call check_allocate(iret, subname, 'zmid(size(Q, 1),nlev)', & file=__FILE__, line=__LINE__) - end if allocate(zk(nlev), stat=iret) @@ -306,12 +315,14 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & pk = moist_pressure_given_z(zk(k),latvals(i)) qk = qv_given_moist_pressure(pk,latvals(i)) else - qk = 0.d0 + qk = 0._r8 + end if + if (lq) then + Q(i,k,m_cnst_ix_q) = qk end if - if (lq .and. cnst1_is_moisture) Q(i,k,1) = qk if (lt) then tvk = Tv_given_z(zk(k),latvals(i)) - T(i,k) = tvk / (1.d0 + Mvap * qk) + T(i,k) = tvk / (1._r8 + Mvap * qk) end if end if end do @@ -339,8 +350,8 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & else qdry = 0.0_r8 end if - if (lq .and. cnst1_is_moisture) then - Q(i,k,1) = qdry + if (lq) then + Q(i,k,m_cnst_ix_q) = qdry end if if (lt) then ! @@ -356,35 +367,58 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, & if(lu .and. masterproc.and. verbose_use) write(iulog,*) ' U initialized by "',subname,'"' if(lv .and. masterproc.and. verbose_use) write(iulog,*) ' V initialized by "',subname,'"' if(lt .and. masterproc.and. verbose_use) write(iulog,*) ' T initialized by "',subname,'"' -!Un-comment once constituents are working in CAMDEN -JN: -#if 0 - if(lq .and. cnst1_is_moisture .and. masterproc.and. verbose_use) write(iulog,*) & - ' ', trim(cnst_name(m_cnst(1))), ' initialized by "',subname,'"' -#endif - end if + if(lq .and. masterproc.and. verbose_use) then + write(iulog,*) ' ', wv_stdname, ' initialized by "',subname,'"' + end if + + end if !l3d_vars -!Un-comment once constituents are working in CAMDEN -JN: -#if 0 if (lq) then - ncnst = size(m_cnst, 1) - do m = 1, ncnst - ! water vapor already done above - if (m_cnst(m) == 1) cycle + !Get constituent properties object: + const_props => cam_model_const_properties() - call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m),& - mask=mask_use, verbose=verbose_use, notfound=.false.,& - z=zmid) + do m = 1, size(m_cnst) - end do + !Skip water vapor, as it was aleady set above: + if (m == m_cnst_ix_q) cycle + + !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 do !m_cnst end if ! lq -#else - if (lq) then - !Initialize cloud liquid and rain until constituent routines are enabled: - Q(:,:,ix_cld_liq) = 0.0_r8 - Q(:,:,ix_rain) = 0.0_r8 - end if -#endif deallocate(mask_use) if (l3d_vars) then 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 diff --git a/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 b/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 index 0f91f131..41a13a21 100644 --- a/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 +++ b/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 @@ -7,19 +7,6 @@ module ic_us_standard_atmosphere ! !------------------------------------------------------------------------------- -use shr_kind_mod, only: r8 => shr_kind_r8 -use spmd_utils, only: masterproc - -use hycoef, only: ps0, hyam, hybm -use physconst, only: gravit -!use constituents, only: cnst_name -!use const_init, only: cnst_init_default - -use std_atm_profile, only: std_atm_pres, std_atm_height, std_atm_temp - -use cam_logfile, only: iulog -use cam_abortutils, only: endrun, check_allocate - implicit none private save @@ -33,6 +20,19 @@ module ic_us_standard_atmosphere subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, & PHIS_OUT, Q, m_cnst, mask, verbose) + use shr_kind_mod, only: r8 => shr_kind_r8, cx => shr_kind_cx + use ccpp_kinds, only: kind_phys + use spmd_utils, only: masterproc + use hycoef, only: ps0, hyam, hybm + use physconst, only: gravit + use std_atm_profile, only: std_atm_pres, std_atm_height, std_atm_temp + use cam_logfile, only: iulog + use cam_abortutils, only: endrun, check_allocate + use runtime_obj, only: wv_stdname + 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 + !---------------------------------------------------------------------------- ! ! Set initial values for static atmosphere with vertical profile from US @@ -58,14 +58,22 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, & ! Local variables logical, allocatable :: mask_use(:) logical :: verbose_use + logical :: const_has_default integer :: i, k, m integer :: ncol integer :: nlev, nlevp integer :: ncnst integer :: iret + integer :: ix_q, m_cnst_ix_q character(len=*), parameter :: subname = 'us_std_atm_set_ic' + character(len=cx) :: errmsg !CCPP error message real(r8) :: psurf(1) real(r8), allocatable :: pmid(:), zmid(:), zmid2d(:,:) + real(kind_phys) :: const_default_value !Constituent default value + real(kind_phys) :: const_qmin_value !Constituent minimum value + + !Private array of constituent properties (for property interface functions) + type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:) !---------------------------------------------------------------------------- ! check input consistency @@ -208,35 +216,67 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, & zmid2d = 0.5_r8*(zint(:,1:nlev) + zint(:,2:nlev+1)) end if - ncnst = size(m_cnst, 1) + !Get water vapor constituent index: + call const_get_index(wv_stdname, ix_q) + + !Determine which "Q" variable index matches water vapor: + m_cnst_ix_q = findloc(m_cnst, ix_q, dim=1) + + !Determine total number of constituents: + ncnst = size(m_cnst) + + !Extract constituent properties from CCPP constituents object: + const_props => cam_model_const_properties() + do m = 1, ncnst - if (m_cnst(m) == 1) then + if (m_cnst(m) == m_cnst_ix_q) then ! No water vapor in profile do k = 1, nlev where(mask_use) Q(:,k,m_cnst(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 - if (present(zint)) then - call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),& - mask=mask_use, verbose=verbose_use, notfound=.false., z=zmid2d) - else - call cnst_init_default(m_cnst(m), latvals, lonvals, Q(:,:,m_cnst(m)),& - mask=mask_use, verbose=verbose_use, notfound=.false.) - end if -#else - else - !Initialize cloud liquid and rain until constituent routines are enabled: - Q(:,:,m_cnst(m)) = 0.0_r8 -#endif - end if - end do + + !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 !ncnst if (allocated(zmid2d)) deallocate(zmid2d) diff --git a/src/dynamics/tests/namelist_definition_analy_ic.xml b/src/dynamics/tests/namelist_definition_analy_ic.xml index 0c17db2f..f71490bc 100644 --- a/src/dynamics/tests/namelist_definition_analy_ic.xml +++ b/src/dynamics/tests/namelist_definition_analy_ic.xml @@ -8,7 +8,9 @@ char*80 dyn_test analytic_ic_nl - none,held_suarez_1994,moist_baroclinic_wave_dcmip2016,dry_baroclinic_wave_dcmip2016,dry_baroclinic_wave_jw2006 + + none,held_suarez_1994,moist_baroclinic_wave_dcmip2016,dry_baroclinic_wave_dcmip2016,dry_baroclinic_wave_jw2006,us_standard_atmosphere + Specify the type of analytic initial conditions for an initial run. held_suarez_1994: Initial conditions specified in Held and Suarez (1994)