From 28d6d7b2c9ff50e545c64a2f8055eea5cc7b173d Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Wed, 4 Sep 2024 09:59:49 -0600 Subject: [PATCH 1/6] Implement constituents in FKESSLER analytic ICs, and fix cpair bug in SE dycore. --- src/control/cam_comp.F90 | 63 ++++----- src/control/runtime_obj.F90 | 3 + src/data/air_composition.F90 | 16 ++- src/dynamics/se/dp_coupling.F90 | 37 ++++-- src/dynamics/se/dycore/prim_state_mod.F90 | 12 +- src/dynamics/se/dyn_comp.F90 | 14 +- .../initial_conditions/ic_baroclinic.F90 | 124 +++++++++++------- 7 files changed, 160 insertions(+), 109 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index b2fba9a2..09b96d78 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -16,7 +16,7 @@ module cam_comp use spmd_utils, only: masterproc, mpicom use cam_control_mod, only: cam_ctrl_init, cam_ctrl_set_orbit - use cam_control_mod, only: cam_ctrl_set_physics_type + use cam_physics_control, only: cam_ctrl_set_physics_type use cam_control_mod, only: caseid, ctitle use runtime_opts, only: read_namelist use runtime_obj, only: cam_runtime_opts @@ -25,8 +25,7 @@ module cam_comp use time_manager, only: is_first_step, is_first_restart_step use camsrfexch, only: cam_out_t, cam_in_t - use physics_types, only: phys_state, phys_tend - use physics_types, only: dtime_phys + use physics_types, only: phys_state, phys_tend, dtime_phys use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -59,7 +58,6 @@ module cam_comp ! Private interface (here to avoid circular dependency) private :: cam_register_constituents - !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------- @@ -87,7 +85,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use dyn_comp, only: dyn_init ! use cam_restart, only: cam_read_restart use camsrfexch, only: hub2atm_alloc, atm2hub_alloc -! use cam_history, only: hist_init_files + use cam_history, only: history_init_files ! use history_scam, only: scm_intht use cam_pio_utils, only: init_pio_subsystem use cam_instance, only: inst_suffix @@ -235,9 +233,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! if (single_column) then ! call scm_intht() ! end if -!!XXgoldyXX: v need to import this -! call hist_init_files(model_doi_url, caseid, ctitle) -!!XXgoldyXX: ^ need to import this + call history_init_files(model_doi_url, caseid, ctitle) end subroutine cam_init @@ -396,10 +392,8 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & ! file output. ! !----------------------------------------------------------------------- -! use cam_history, only: wshist, wrapup ! use cam_restart, only: cam_write_restart ! use qneg_module, only: qneg_print_summary - use time_manager, only: is_last_step type(cam_out_t), intent(inout) :: cam_out ! Output from CAM to surface type(cam_in_t), intent(inout) :: cam_in ! Input from surface to CAM @@ -410,18 +404,6 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & integer, intent(in), optional :: day_spec ! Simulation day integer, intent(in), optional :: sec_spec ! Secs in current simulation day - !---------------------------------------------------------- - ! History and restart logic: Write and/or dispose history - ! tapes if required - !---------------------------------------------------------- - ! -!!XXgoldyXX: v need to import this -! call t_barrierf('sync_wshist', mpicom) -! call t_startf('wshist') -! call wshist() -! call t_stopf('wshist') -!!XXgoldyXX: ^ need to import this - ! ! Write restart files ! @@ -441,37 +423,39 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & call t_stopf('cam_write_restart') end if -!!XXgoldyXX: v need to import this -! call t_startf ('cam_run4_wrapup') -! call wrapup(rstwr, nlend) -! call t_stopf ('cam_run4_wrapup') -!!XXgoldyXX: ^ need to import this - - call shr_sys_flush(iulog) - end subroutine cam_run4 ! !----------------------------------------------------------------------- ! - subroutine cam_timestep_final(do_ncdata_check) + subroutine cam_timestep_final(rstwr, nlend, do_ncdata_check) !----------------------------------------------------------------------- ! ! Purpose: Timestep final runs at the end of each timestep ! !----------------------------------------------------------------------- - use phys_comp, only: phys_timestep_final - + use phys_comp, only: phys_timestep_final + use cam_history, only: history_write_files + use cam_history, only: history_wrap_up + logical, intent(in) :: rstwr ! write restart file + logical, intent(in) :: nlend ! this is final timestep !Flag for whether a snapshot (ncdata) check should be run or not - logical, intent(in) :: do_ncdata_check + ! - flag is true if this is not the first or last step + logical, intent(in) :: do_ncdata_check + if (do_ncdata_check .or. get_nstep() == 0) then + call history_write_files() + ! peverwhee - todo: handle restarts + call history_wrap_up(rstwr, nlend) + end if ! !---------------------------------------------------------- ! PHYS_TIMESTEP_FINAL Call the Physics package !---------------------------------------------------------- ! call phys_timestep_final(do_ncdata_check) + call shr_sys_flush(iulog) end subroutine cam_timestep_final @@ -551,6 +535,7 @@ subroutine cam_register_constituents(cam_runtime_opts) use cam_ccpp_cap, only: cam_ccpp_number_constituents use cam_ccpp_cap, only: cam_model_const_properties use cam_ccpp_cap, only: cam_ccpp_is_scheme_constituent + use runtime_obj, only: wv_stdname ! Dummy arguments type(runtime_options), intent(in) :: cam_runtime_opts @@ -561,6 +546,7 @@ subroutine cam_register_constituents(cam_runtime_opts) integer :: errflg character(len=512) :: errmsg type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:) + type(ccpp_constituent_properties_t), allocatable, target :: dynamic_constituents(:) character(len=*), parameter :: subname = 'cam_register_constituents: ' ! Initalize error flag and message: @@ -569,9 +555,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 +573,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 +623,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..f30cb2ab 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 + ! 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..d4056a5d 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 !========================================================================================= @@ -553,7 +555,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 + use cam_constituents, only: const_get_index, const_qmin use physics_types, only: lagrangian_vertical use physconst, only: cpair, gravit, zvir, cappa use cam_thermo, only: cam_thermo_update @@ -561,7 +563,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 @@ -587,6 +589,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) integer :: errflg character(len=shr_kind_cx) :: errmsg + character(len=*), parameter :: subname = 'derived_phys_dry' + !-------------------------------------------- ! Variables needed for WACCM-X !-------------------------------------------- @@ -610,6 +614,20 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) 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 @@ -736,6 +754,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 +766,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 +784,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..e391dd6f 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -9,7 +9,8 @@ module dyn_comp use cam_constituents, only: const_get_index, const_is_wet, const_qmin use cam_constituents, only: readtrace use air_composition, only: const_is_water_species -use cam_control_mod, only: initial_run, simple_phys +use cam_control_mod, only: initial_run +use cam_physics_control, only: simple_phys use cam_initfiles, only: initial_file_get_id, topo_file_get_id, pertlim use dyn_grid, only: ini_grid_name, timelevel, hvcoord, edgebuf @@ -1307,7 +1308,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 +1319,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 +1351,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_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 From fa0b5e1c0c22855f8813a7c70d26415a6447f14a Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 6 Sep 2024 13:59:04 -0600 Subject: [PATCH 2/6] Properly implement constituents in Held-Suarez ICs, and fix wet-to-dry conversions in dp_coupling. --- src/dynamics/se/dp_coupling.F90 | 63 +++++++++----- .../initial_conditions/ic_held_suarez.F90 | 84 ++++++++++++++----- 2 files changed, 104 insertions(+), 43 deletions(-) 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 From 0d7919c524ed3afba37c1ecc407ada8870fa73bd Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 10 Sep 2024 15:08:50 -0600 Subject: [PATCH 3/6] Add constituent hooks to remaining analytic IC options. --- .../initial_conditions/ic_baro_dry_jw06.F90 | 98 ++++++++++++------- .../initial_conditions/ic_us_standard_atm.F90 | 95 ++++++++++++------ .../tests/namelist_definition_analy_ic.xml | 4 +- 3 files changed, 130 insertions(+), 67 deletions(-) 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..7f894168 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,63 @@ 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 + ncnst = size(m_cnst) 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) + + 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 + end if !l3d_vars deallocate(mask_use) 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..5bb19cbe 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,58 @@ 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) + + ncnst = size(m_cnst) 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 - 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.) + write(iulog,*) ' ', wv_stdname, ' initialized by '//subname 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) + + 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) From 894c3532f9c28ebcb8b2435ad51095475e2994f5 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 10 Sep 2024 15:38:01 -0600 Subject: [PATCH 4/6] Remove code added by accident during testing. --- src/control/cam_comp.F90 | 55 +++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 20 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 09b96d78..0782738a 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -16,7 +16,7 @@ module cam_comp use spmd_utils, only: masterproc, mpicom use cam_control_mod, only: cam_ctrl_init, cam_ctrl_set_orbit - use cam_physics_control, only: cam_ctrl_set_physics_type + use cam_control_mod, only: cam_ctrl_set_physics_type use cam_control_mod, only: caseid, ctitle use runtime_opts, only: read_namelist use runtime_obj, only: cam_runtime_opts @@ -25,7 +25,8 @@ module cam_comp use time_manager, only: is_first_step, is_first_restart_step use camsrfexch, only: cam_out_t, cam_in_t - use physics_types, only: phys_state, phys_tend, dtime_phys + use physics_types, only: phys_state, phys_tend + use physics_types, only: dtime_phys use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -58,6 +59,7 @@ module cam_comp ! Private interface (here to avoid circular dependency) private :: cam_register_constituents + !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------- @@ -85,7 +87,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use dyn_comp, only: dyn_init ! use cam_restart, only: cam_read_restart use camsrfexch, only: hub2atm_alloc, atm2hub_alloc - use cam_history, only: history_init_files +! use cam_history, only: hist_init_files ! use history_scam, only: scm_intht use cam_pio_utils, only: init_pio_subsystem use cam_instance, only: inst_suffix @@ -233,7 +235,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! if (single_column) then ! call scm_intht() ! end if - call history_init_files(model_doi_url, caseid, ctitle) +!!XXgoldyXX: v need to import this +! call hist_init_files(model_doi_url, caseid, ctitle) +!!XXgoldyXX: ^ need to import this end subroutine cam_init @@ -392,8 +396,10 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & ! file output. ! !----------------------------------------------------------------------- +! use cam_history, only: wshist, wrapup ! use cam_restart, only: cam_write_restart ! use qneg_module, only: qneg_print_summary + use time_manager, only: is_last_step type(cam_out_t), intent(inout) :: cam_out ! Output from CAM to surface type(cam_in_t), intent(inout) :: cam_in ! Input from surface to CAM @@ -404,6 +410,18 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & integer, intent(in), optional :: day_spec ! Simulation day integer, intent(in), optional :: sec_spec ! Secs in current simulation day + !---------------------------------------------------------- + ! History and restart logic: Write and/or dispose history + ! tapes if required + !---------------------------------------------------------- + ! +!!XXgoldyXX: v need to import this +! call t_barrierf('sync_wshist', mpicom) +! call t_startf('wshist') +! call wshist() +! call t_stopf('wshist') +!!XXgoldyXX: ^ need to import this + ! ! Write restart files ! @@ -423,39 +441,37 @@ subroutine cam_run4(cam_out, cam_in, rstwr, nlend, & call t_stopf('cam_write_restart') end if +!!XXgoldyXX: v need to import this +! call t_startf ('cam_run4_wrapup') +! call wrapup(rstwr, nlend) +! call t_stopf ('cam_run4_wrapup') +!!XXgoldyXX: ^ need to import this + + call shr_sys_flush(iulog) + end subroutine cam_run4 ! !----------------------------------------------------------------------- ! - subroutine cam_timestep_final(rstwr, nlend, do_ncdata_check) + subroutine cam_timestep_final(do_ncdata_check) !----------------------------------------------------------------------- ! ! Purpose: Timestep final runs at the end of each timestep ! !----------------------------------------------------------------------- - use phys_comp, only: phys_timestep_final - use cam_history, only: history_write_files - use cam_history, only: history_wrap_up - logical, intent(in) :: rstwr ! write restart file - logical, intent(in) :: nlend ! this is final timestep + use phys_comp, only: phys_timestep_final + !Flag for whether a snapshot (ncdata) check should be run or not - ! - flag is true if this is not the first or last step - logical, intent(in) :: do_ncdata_check + logical, intent(in) :: do_ncdata_check - if (do_ncdata_check .or. get_nstep() == 0) then - call history_write_files() - ! peverwhee - todo: handle restarts - call history_wrap_up(rstwr, nlend) - end if ! !---------------------------------------------------------- ! PHYS_TIMESTEP_FINAL Call the Physics package !---------------------------------------------------------- ! call phys_timestep_final(do_ncdata_check) - call shr_sys_flush(iulog) end subroutine cam_timestep_final @@ -526,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 @@ -535,7 +552,6 @@ subroutine cam_register_constituents(cam_runtime_opts) use cam_ccpp_cap, only: cam_ccpp_number_constituents use cam_ccpp_cap, only: cam_model_const_properties use cam_ccpp_cap, only: cam_ccpp_is_scheme_constituent - use runtime_obj, only: wv_stdname ! Dummy arguments type(runtime_options), intent(in) :: cam_runtime_opts @@ -546,7 +562,6 @@ subroutine cam_register_constituents(cam_runtime_opts) integer :: errflg character(len=512) :: errmsg type(ccpp_constituent_prop_ptr_t), pointer :: const_props(:) - type(ccpp_constituent_properties_t), allocatable, target :: dynamic_constituents(:) character(len=*), parameter :: subname = 'cam_register_constituents: ' ! Initalize error flag and message: From e39a48a4544266808bb426cca87eec0d11ce6311 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 10 Sep 2024 15:56:34 -0600 Subject: [PATCH 5/6] Fix comment. --- src/data/air_composition.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index f30cb2ab..f202ae97 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -459,7 +459,7 @@ subroutine air_composition_init() !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 - ! instead of just N2 in this configuration + ! 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 From a28e7be8d951fc845b6ea0736490f1abaa4fe639 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 13 Sep 2024 14:56:10 -0600 Subject: [PATCH 6/6] Implement code review suggestions. --- src/dynamics/se/dyn_comp.F90 | 3 +-- .../initial_conditions/ic_baro_dry_jw06.F90 | 21 +++++++++++++------ .../initial_conditions/ic_us_standard_atm.F90 | 11 +++++++++- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index e391dd6f..8d74b66e 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -9,8 +9,7 @@ module dyn_comp use cam_constituents, only: const_get_index, const_is_wet, const_qmin use cam_constituents, only: readtrace use air_composition, only: const_is_water_species -use cam_control_mod, only: initial_run -use cam_physics_control, only: simple_phys +use cam_control_mod, only: initial_run, simple_phys use cam_initfiles, only: initial_file_get_id, topo_file_get_id, pertlim use dyn_grid, only: ini_grid_name, timelevel, hvcoord, edgebuf 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 7f894168..87bdf868 100644 --- a/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 +++ b/src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90 @@ -243,11 +243,11 @@ 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) + !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 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) @@ -261,7 +261,12 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & end if 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 do m = 1, ncnst @@ -274,6 +279,10 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & !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 @@ -298,8 +307,8 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, & end if !has_default end do !m_cnst - end if !lq - end if !l3d_vars + end if !vcoord + end if !lq deallocate(mask_use) 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 5bb19cbe..41a13a21 100644 --- a/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 +++ b/src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90 @@ -222,7 +222,12 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, & !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) == m_cnst_ix_q) then ! No water vapor in profile @@ -242,8 +247,12 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, & !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__) + call endrun(errmsg, file=__FILE__, line=__LINE__) end if if (const_has_default) then