diff --git a/Externals_CAM.cfg b/Externals_CAM.cfg index bfba41d0db..f58df387ea 100644 --- a/Externals_CAM.cfg +++ b/Externals_CAM.cfg @@ -55,9 +55,9 @@ tag = ALI_ARMS_v1.0.1 required = True [atmos_phys] -tag = atmos_phys0_00_011 +tag = atmos_phys0_01_000 protocol = git -repo_url = https://github.com/NCAR/atmospheric_physics +repo_url = https://github.com/ESCOMP/atmospheric_physics required = True local_path = src/atmos_phys diff --git a/bld/configure b/bld/configure index 262ac38e6d..d1f86cb03e 100755 --- a/bld/configure +++ b/bld/configure @@ -2269,6 +2269,7 @@ sub write_filepath print $fh "$camsrcdir/src/cpl/$cpl\n"; print $fh "$camsrcdir/src/control\n"; print $fh "$camsrcdir/src/utils\n"; + print $fh "$camsrcdir/src/utils/cam_ccpp\n"; print $fh "$camsrcdir/src/atmos_phys/utilities\n"; diff --git a/doc/ChangeLog b/doc/ChangeLog index 4a2c43d13f..7e45ac811b 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,133 @@ =============================================================== +Tag name: cam6_3_134 +Originator(s): nusbaume, jimmielin +Date: 31 Oct 2023 +One-line Summary: Update atmospheric_physics external +Github PR URL: https://github.com/ESCOMP/CAM/pull/891 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + +This PR updates the atmospheric_physics external in CAM, +which requires some source modifications on the CAM-side. +There has also been some replacement of CAM code with CCPP +physics routines in order to reduce duplication. + +This PR also fixes a bug found in Kessler where it was using +the wrong pressure when converting to/from potential temperature +via the exner function. + +Finally, this change also brings in a rename of the SE +dycore's 'time_mod' module to 'se_dyn_time_mod' in order +to avoid name collision with GEOS-Chem code. + +Fixes #752 -> Update atmospheric_physics external +Fixes #802 -> Kessler physics using inconsistent reference pressures + +Closes #904 -> Rename SE dycore time_mod to dyn_time_mod (PR) + +Describe any changes made to build system: + +Added the "src/utils/cam_ccpp" and "src/atmos_phys/utilities" directory +to the list of directories used during compilation. + +Describe any changes made to the namelist: N/A + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: cacraigucar, PeterHjortLauritzen, jtruesdal, peverwhee + +List all files eliminated: + +R src/dynamics/se/dycore/time_mod.F90 + - Renamed to 'se_dyn_time_mod' +R src/utils/ccpp_kinds.F90 + - Moved to 'src/utils/cam_ccpp/ccpp_kinds.F90' + +List all files added and what they do: + +A src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 + - Creates a "CCPP Constituent Properties" object for CAM. + +List all existing files that have been modified, and describe the changes: + +M Externals_CAM.cfg + - Update 'atmospheric_physics' external + +M bld/configure + - Add 'src/utils/cam_ccpp' and 'src/atmos_phys/utilities' to build file paths + +M src/dynamics/se/dp_coupling.F90 + - Use CCPP-ized 'update_dry_static_energy_run' subroutine + +M src/physics/cam/geopotential.F90 + - Use CCPP-ized 'geopotential_temp_run' subroutine + +M src/physics/cam/ref_pres.F90 + - Update comment in order to avoid future confusion + +M src/physics/cam/physpkg.F90 +M src/physics/cam_dev/physpkg.F90 +M src/physics/simple/held_suarez_cam.F90 +M src/physics/simple/physpkg.F90 + - Update physics calls and add new const. prop. DDT to work with new + atmospheric_physics external + +M src/physics/simple/kessler_cam.F90 +- Update physics calls and add new const. prop. DDT to work with new + atmospheric_physics external. Also replace use of "standard pressure" + with "reference pressure" in order to avoid unphysical mismatches. + +M src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 +M src/dynamics/se/dycore/prim_advance_mod.F90 +M src/dynamics/se/dycore/prim_advection_mod.F90 +M src/dynamics/se/dycore/prim_driver_mod.F90 +M src/dynamics/se/dycore/prim_init.F90 +M src/dynamics/se/dycore/prim_state_mod.F90 +M src/dynamics/se/stepon.F90 + - Replace 'time_mod' with 'se_dyn_time_mod' + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +cheyenne/intel/aux_cam: + + All tests have NLCOMP failures due to changes in input data paths on Cheyenne. + + ERP_Ln9_Vnuopc.C96_C96_mg17.F2000climo.cheyenne_intel.cam-outfrq9s_mg3 (Overall: FAIL) details: + FAIL ERP_Ln9_Vnuopc.C96_C96_mg17.F2000climo.cheyenne_intel.cam-outfrq9s_mg3 NLCOMP + FAIL ERP_Ln9_Vnuopc.C96_C96_mg17.F2000climo.cheyenne_intel.cam-outfrq9s_mg3 MODEL_BUILD time=3 + ERP_Ln9_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq9s (Overall: FAIL) details: + FAIL ERP_Ln9_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq9s NLCOMP + FAIL ERP_Ln9_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq9s COMPARE_base_rest + SMS_Lh12_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq3h (Overall: DIFF) details: + FAIL SMS_Lh12_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq3h NLCOMP + FAIL SMS_Lh12_Vnuopc.f09_f09_mg17.FCSD_HCO.cheyenne_intel.cam-outfrq3h BASELINE /glade/p/cesm/amwg/cesm_baselines/cam6_3_133: DIFF + - pre-existing failure + +izumi/nag/aux_cam: + + DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details: + FAIL DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae RUN time=9 + - pre-existing failure + + ERS_Ln27_Vnuopc.ne5pg3_ne5pg3_mg37.FKESSLER.izumi_nag.cam-outfrq9s (Overall: DIFF) details: + FAIL ERS_Ln27_Vnuopc.ne5pg3_ne5pg3_mg37.FKESSLER.izumi_nag.cam-outfrq9s BASELINE /fs/cgd/csm/models/atm/cam/pretag_bl/cam6_3_133_nag: DIFF + - expected failure + +izumi/gnu/aux_cam: All PASS + +Summarize any changes to answers: + All compsets which use Kessler idealized physics + will have larger-than-roundoff answer changes due + to reference pressure bug fix. + +=============================================================== + Tag name: cam6_3_133 Originator(s): fvitt Date: 19 Oct 2023 diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 056eb45f93..61b7fc54e9 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -54,7 +54,7 @@ subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out) use phys_control, only: use_gw_front, use_gw_front_igw use hycoef, only: hyai, ps0 use fvm_mapping, only: dyn2phys_vector, dyn2phys_all_vars - use time_mod, only: timelevel_qdp + use se_dyn_time_mod, only: timelevel_qdp use control_mod, only: qsplit use test_fvm_mapping, only: test_mapping_overwrite_dyn_state, test_mapping_output_phys_state use prim_advance_mod, only: tot_energy_dyn @@ -546,15 +546,17 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) use shr_const_mod, only: shr_const_rwv use phys_control, only: waccmx_is use geopotential, only: geopotential_t + use static_energy, only: update_dry_static_energy_run use check_energy, only: check_energy_timestep_init use hycoef, only: hyai, ps0 use shr_vmath_mod, only: shr_vmath_log use qneg_module, only: qneg3 use dyn_tests_utils, only: vc_dry_pressure + use shr_kind_mod, only: shr_kind_cx ! arguments type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend - type(physics_buffer_desc), pointer :: pbuf2d(:,:) + type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! local variables integer :: lchnk @@ -564,6 +566,10 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) integer :: m, i, k, ncol, m_cnst type(physics_buffer_desc), pointer :: pbuf_chnk(:) + + !Needed for "update_dry_static_energy" CCPP scheme + integer :: errflg + character(len=shr_kind_cx) :: errmsg !---------------------------------------------------------------------------- ! Evaluate derived quantities @@ -686,18 +692,17 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) end do ! Compute initial geopotential heights - based on full pressure - call geopotential_t (phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint , & - phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , & - phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv , & - phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol ) + call geopotential_t(phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint, & + phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , & + phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv , & + phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol) ! Compute initial dry static energy, include surface geopotential - do k = 1, pver - do i = 1, ncol - phys_state(lchnk)%s(i,k) = cpairv(i,k,lchnk)*phys_state(lchnk)%t(i,k) & - + gravit*phys_state(lchnk)%zm(i,k) + phys_state(lchnk)%phis(i) - end do - end do + call update_dry_static_energy_run(pver, gravit, phys_state(lchnk)%t(1:ncol,:), & + phys_state(lchnk)%zm(1:ncol,:), & + phys_state(lchnk)%phis(1:ncol), & + phys_state(lchnk)%s(1:ncol,:), & + cpairv(1:ncol,:,lchnk), errflg, errmsg) ! Ensure tracers are all positive call qneg3('D_P_COUPLING',lchnk ,ncol ,pcols ,pver , & diff --git a/src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 b/src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 index b5bcebd845..5da18d76b8 100644 --- a/src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 +++ b/src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 @@ -6,7 +6,7 @@ module fvm_consistent_se_cslam use cam_abortutils, only: endrun use cam_logfile, only: iulog - use time_mod, only: timelevel_t + use se_dyn_time_mod, only: timelevel_t use element_mod, only: element_t use fvm_control_volume_mod, only: fvm_struct use hybrid_mod, only: hybrid_t, config_thread_region, get_loop_ranges, threadOwnsVertLevel diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90 index 1246b77e34..ed7a627ec4 100644 --- a/src/dynamics/se/dycore/prim_advance_mod.F90 +++ b/src/dynamics/se/dycore/prim_advance_mod.F90 @@ -53,7 +53,7 @@ subroutine prim_advance_exp(elem, fvm, deriv, hvcoord, hybrid,dt, tl, nets, net use element_mod, only: element_t use hybvcoord_mod, only: hvcoord_t use hybrid_mod, only: hybrid_t - use time_mod, only: TimeLevel_t, timelevel_qdp, tevolve + use se_dyn_time_mod, only: TimeLevel_t, timelevel_qdp, tevolve use fvm_control_volume_mod, only: fvm_struct use cam_thermo, only: get_kappa_dry use air_composition, only: thermodynamic_active_species_num diff --git a/src/dynamics/se/dycore/prim_advection_mod.F90 b/src/dynamics/se/dycore/prim_advection_mod.F90 index 7c54abc2cd..17ad85ba61 100644 --- a/src/dynamics/se/dycore/prim_advection_mod.F90 +++ b/src/dynamics/se/dycore/prim_advection_mod.F90 @@ -23,7 +23,7 @@ module prim_advection_mod use element_mod, only: element_t use fvm_control_volume_mod, only: fvm_struct use hybvcoord_mod, only: hvcoord_t - use time_mod, only: TimeLevel_t, TimeLevel_Qdp + use se_dyn_time_mod, only: TimeLevel_t, TimeLevel_Qdp use control_mod, only: nu_q, nu_p, limiter_option, hypervis_subcycle_q, rsplit use edge_mod, only: edgevpack, edgevunpack, initedgebuffer, initedgesbuffer diff --git a/src/dynamics/se/dycore/prim_driver_mod.F90 b/src/dynamics/se/dycore/prim_driver_mod.F90 index af22869f24..6cfb52e356 100644 --- a/src/dynamics/se/dycore/prim_driver_mod.F90 +++ b/src/dynamics/se/dycore/prim_driver_mod.F90 @@ -28,8 +28,8 @@ subroutine prim_init2(elem, fvm, hybrid, nets, nete, tl, hvcoord) use dimensions_mod, only: irecons_tracer, fvm_supercycling use dimensions_mod, only: fv_nphys, nc use parallel_mod, only: syncmp - use time_mod, only: timelevel_t, tstep, phys_tscale, nsplit, TimeLevel_Qdp - use time_mod, only: nsplit_baseline,rsplit_baseline + use se_dyn_time_mod, only: timelevel_t, tstep, phys_tscale, nsplit, TimeLevel_Qdp + use se_dyn_time_mod, only: nsplit_baseline,rsplit_baseline use prim_state_mod, only: prim_printstate use control_mod, only: runtype, topology, rsplit, qsplit, rk_stage_user, & nu, nu_q, nu_div, hypervis_subcycle, hypervis_subcycle_q, & @@ -218,7 +218,7 @@ subroutine prim_run_subcycle(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,nsubst ! ! use hybvcoord_mod, only : hvcoord_t - use time_mod, only: TimeLevel_t, timelevel_update, timelevel_qdp, nsplit + use se_dyn_time_mod, only: TimeLevel_t, timelevel_update, timelevel_qdp, nsplit use control_mod, only: statefreq,qsplit, rsplit, variable_nsplit use prim_advance_mod, only: applycamforcing use prim_advance_mod, only: tot_energy_dyn,compute_omega @@ -407,7 +407,7 @@ subroutine prim_step(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord, rstep) ! tl%n0 time t + dt_q ! use hybvcoord_mod, only: hvcoord_t - use time_mod, only: TimeLevel_t, timelevel_update + use se_dyn_time_mod, only: TimeLevel_t, timelevel_update use control_mod, only: statefreq, qsplit, nu_p use thread_mod, only: omp_get_thread_num use prim_advance_mod, only: prim_advance_exp diff --git a/src/dynamics/se/dycore/prim_init.F90 b/src/dynamics/se/dycore/prim_init.F90 index afbd94869e..42a336f65c 100644 --- a/src/dynamics/se/dycore/prim_init.F90 +++ b/src/dynamics/se/dycore/prim_init.F90 @@ -28,7 +28,7 @@ subroutine prim_init1(elem, fvm, par, Tl) use element_mod, only: element_t, allocate_element_desc use fvm_mod, only: fvm_init1 use mesh_mod, only: MeshUseMeshFile - use time_mod, only: timelevel_init, timelevel_t + use se_dyn_time_mod, only: timelevel_init, timelevel_t use mass_matrix_mod, only: mass_matrix use derivative_mod, only: allocate_subcell_integration_matrix_cslam use derivative_mod, only: allocate_subcell_integration_matrix_physgrid diff --git a/src/dynamics/se/dycore/prim_state_mod.F90 b/src/dynamics/se/dycore/prim_state_mod.F90 index 2f4bcbb2db..3075c0e125 100644 --- a/src/dynamics/se/dycore/prim_state_mod.F90 +++ b/src/dynamics/se/dycore/prim_state_mod.F90 @@ -4,7 +4,7 @@ module prim_state_mod use dimensions_mod, only: nlev, np, nc, qsize_d, ntrac_d use parallel_mod, only: ordered use hybrid_mod, only: hybrid_t - use time_mod, only: timelevel_t, TimeLevel_Qdp, time_at + use se_dyn_time_mod, only: timelevel_t, TimeLevel_Qdp, time_at use control_mod, only: qsplit, statediag_numtrac use global_norms_mod, only: global_integrals_general use element_mod, only: element_t @@ -24,7 +24,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn) use air_composition, only: thermodynamic_active_species_idx_dycore, dry_air_species_num use air_composition, only: thermodynamic_active_species_num,thermodynamic_active_species_idx use cam_control_mod, only: initial_run - use time_mod, only: tstep + use se_dyn_time_mod, only: tstep use control_mod, only: rsplit, qsplit use perf_mod, only: t_startf, t_stopf type (element_t), intent(inout) :: elem(:) @@ -345,10 +345,10 @@ end subroutine prim_printstate_cslam_gamma subroutine adjust_nsplit(elem, tl,hybrid,nets,nete, fvm, omega_cn) use dimensions_mod, only: ksponge_end use dimensions_mod, only: fvm_supercycling, fvm_supercycling_jet - use time_mod, only: tstep + use se_dyn_time_mod, only: tstep use control_mod, only: rsplit, qsplit use perf_mod, only: t_startf, t_stopf - use time_mod, only: nsplit, nsplit_baseline,rsplit_baseline + use se_dyn_time_mod, only: nsplit, nsplit_baseline,rsplit_baseline use control_mod, only: qsplit, rsplit use time_manager, only: get_step_size use cam_abortutils, only: endrun diff --git a/src/dynamics/se/dycore/time_mod.F90 b/src/dynamics/se/dycore/se_dyn_time_mod.F90 similarity index 98% rename from src/dynamics/se/dycore/time_mod.F90 rename to src/dynamics/se/dycore/se_dyn_time_mod.F90 index fdd68af06a..4dfd981661 100644 --- a/src/dynamics/se/dycore/time_mod.F90 +++ b/src/dynamics/se/dycore/se_dyn_time_mod.F90 @@ -1,4 +1,4 @@ -module time_mod +module se_dyn_time_mod !------------------ use shr_kind_mod, only: r8=>shr_kind_r8 !------------------ @@ -132,4 +132,4 @@ subroutine TimeLevel_update(tl,uptype) !$OMP BARRIER end subroutine TimeLevel_update -end module time_mod +end module se_dyn_time_mod diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index ac46b8cc46..312349eb44 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -42,7 +42,7 @@ module dyn_comp use dimensions_mod, only: qsize, use_cslam use element_mod, only: element_t, elem_state_t use fvm_control_volume_mod, only: fvm_struct -use time_mod, only: nsplit +use se_dyn_time_mod, only: nsplit use edge_mod, only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer use edgetype_mod, only: EdgeBuffer_t use bndry_mod, only: bndry_exchange @@ -963,11 +963,11 @@ subroutine dyn_run(dyn_state) use air_composition, only: thermodynamic_active_species_idx_dycore use prim_driver_mod, only: prim_run_subcycle use dimensions_mod, only: cnst_name_gll - use time_mod, only: tstep, nsplit, timelevel_qdp + use se_dyn_time_mod, only: tstep, nsplit, timelevel_qdp use hybrid_mod, only: config_thread_region, get_loop_ranges use control_mod, only: qsplit, rsplit, ftype_conserve use thread_mod, only: horz_num_threads - use time_mod, only: tevolve + use se_dyn_time_mod, only: tevolve type(dyn_export_t), intent(inout) :: dyn_state diff --git a/src/dynamics/se/dyn_grid.F90 b/src/dynamics/se/dyn_grid.F90 index 766fb893d7..293f7402dd 100644 --- a/src/dynamics/se/dyn_grid.F90 +++ b/src/dynamics/se/dyn_grid.F90 @@ -48,7 +48,7 @@ module dyn_grid use prim_init, only: prim_init1 use edge_mod, only: initEdgeBuffer use edgetype_mod, only: EdgeBuffer_t -use time_mod, only: TimeLevel_t +use se_dyn_time_mod, only: TimeLevel_t use dof_mod, only: UniqueCoords, UniquePoints implicit none @@ -133,7 +133,7 @@ subroutine dyn_grid_init() use hybrid_mod, only: hybrid_t, init_loop_ranges, & get_loop_ranges, config_thread_region use control_mod, only: qsplit, rsplit - use time_mod, only: tstep, nsplit + use se_dyn_time_mod, only: tstep, nsplit use fvm_mod, only: fvm_init2, fvm_init3, fvm_pg_init use dimensions_mod, only: irecons_tracer use comp_gll_ctr_vol, only: gll_grid_write diff --git a/src/dynamics/se/restart_dynamics.F90 b/src/dynamics/se/restart_dynamics.F90 index f5b3c6a0d8..0c630c9336 100644 --- a/src/dynamics/se/restart_dynamics.F90 +++ b/src/dynamics/se/restart_dynamics.F90 @@ -46,7 +46,7 @@ module restart_dynamics use dimensions_mod, only: np, npsq, ne, nlev, qsize, nelemd, nc, ntrac, use_cslam use dof_mod, only: UniquePoints use element_mod, only: element_t -use time_mod, only: tstep, TimeLevel_Qdp +use se_dyn_time_mod, only: tstep, TimeLevel_Qdp use edge_mod, only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer use edgetype_mod, only: EdgeBuffer_t diff --git a/src/dynamics/se/stepon.F90 b/src/dynamics/se/stepon.F90 index 88dda66c3d..60e73a0ece 100644 --- a/src/dynamics/se/stepon.F90 +++ b/src/dynamics/se/stepon.F90 @@ -99,7 +99,7 @@ subroutine stepon_run1( dtime_out, phys_state, phys_tend, & use dp_coupling, only: d_p_coupling use physics_buffer, only: physics_buffer_desc - use time_mod, only: tstep ! dynamics timestep + use se_dyn_time_mod,only: tstep ! dynamics timestep real(r8), intent(out) :: dtime_out ! Time-step type(physics_state), intent(inout) :: phys_state(begchunk:endchunk) @@ -152,7 +152,7 @@ subroutine stepon_run2(phys_state, phys_tend, dyn_in, dyn_out) use dp_coupling, only: p_d_coupling use dyn_grid, only: TimeLevel - use time_mod, only: TimeLevel_Qdp + use se_dyn_time_mod, only: TimeLevel_Qdp use control_mod, only: qsplit use prim_advance_mod, only: tot_energy_dyn @@ -207,7 +207,7 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out) use dyn_comp, only: dyn_run use advect_tend, only: compute_adv_tends_xyz use dyn_grid, only: TimeLevel - use time_mod, only: TimeLevel_Qdp + use se_dyn_time_mod,only: TimeLevel_Qdp use control_mod, only: qsplit ! arguments real(r8), intent(in) :: dtime ! Time-step @@ -254,7 +254,7 @@ subroutine diag_dynvar_ic(elem, fvm) use cam_history, only: write_inithist, outfld, hist_fld_active, fieldname_len use dyn_grid, only: TimeLevel - use time_mod, only: TimeLevel_Qdp ! dynamics typestep + use se_dyn_time_mod, only: TimeLevel_Qdp ! dynamics typestep use control_mod, only: qsplit use hybrid_mod, only: config_thread_region, get_loop_ranges use hybrid_mod, only: hybrid_t diff --git a/src/physics/cam/geopotential.F90 b/src/physics/cam/geopotential.F90 index 7672fc92e0..ad49e470c4 100644 --- a/src/physics/cam/geopotential.F90 +++ b/src/physics/cam/geopotential.F90 @@ -7,7 +7,7 @@ module geopotential ! The hydrostatic matrix elements must be consistent with the dynamics algorithm. ! The diagonal element is the itegration weight from interface k+1 to midpoint k. ! The offdiagonal element is the weight between interfaces. -! +! ! Author: B.Boville, Feb 2001 from earlier code by Boville and S.J. Lin !--------------------------------------------------------------------------------- @@ -29,17 +29,18 @@ subroutine geopotential_t( & t , q , rair , gravit , zvir , & zi , zm , ncol ) -!----------------------------------------------------------------------- -! -! Purpose: -! Compute the geopotential height (above the surface) at the midpoints and +!----------------------------------------------------------------------- +! +! Purpose: +! Compute the geopotential height (above the surface) at the midpoints and ! interfaces using the input temperatures and pressures. ! !----------------------------------------------------------------------- -use ppgrid, only : pcols -use air_composition, only: thermodynamic_active_species_num -use air_composition, only: thermodynamic_active_species_idx, dry_air_species_num +use ppgrid, only: pcols +use constituents, only: pcnst, cnst_get_ind +use ccpp_constituent_prop_mod, only: ccpp_const_props !CCPP constituent properties array (CAM version) +use geopotential_temp, only: geopotential_temp_run !CCPP version !------------------------------Arguments-------------------------------- ! ! Input arguments @@ -65,6 +66,8 @@ subroutine geopotential_t( & ! !---------------------------Local variables----------------------------- ! + logical :: lagrang ! Lagrangian vertical coordinate flag + integer :: ixq ! state constituent array index for water vapor integer :: i,k,idx ! Lon, level indices, water species index real(r8) :: hkk(ncol) ! diagonal element of hydrostatic matrix real(r8) :: hkl(ncol) ! off-diagonal element @@ -74,28 +77,35 @@ subroutine geopotential_t( & real(r8) :: qfac(ncol,pver) ! factor to convert from wet to dry mixing ratio real(r8) :: sum_dry_mixing_ratio(ncol,pver)! sum of dry water mixing ratios + !CCPP-required variables (not used): + integer :: errflg + character(len=512) :: errmsg + ! !----------------------------------------------------------------------- ! - rog(:ncol,:) = rair(:ncol,:) / gravit - -! The surface height is zero by definition. - - do i = 1,ncol - zi(i,pverp) = 0.0_r8 - end do - - ! Compute zi, zm from bottom up. - ! Note, zi(i,k) is the interface above zm(i,k) + !Determine index for water vapor mass mixing ratio + call cnst_get_ind('Q', ixq) ! ! original code for backwards compatability with FV and EUL ! if (.not.(dycore_is('MPAS') .or. dycore_is('SE'))) then + + !dry air gas constant over gravity + rog(:ncol,:) = rair(:ncol,:) / gravit + + ! The surface height is zero by definition. + do i = 1,ncol + zi(i,pverp) = 0.0_r8 + end do + + ! Compute zi, zm from bottom up. + ! Note, zi(i,k) is the interface above zm(i,k) do k = pver, 1, -1 - + ! First set hydrostatic elements consistent with dynamics - + if ((dycore_is('LR') .or. dycore_is('FV3'))) then do i = 1,ncol hkl(i) = piln(i,k+1) - piln(i,k) @@ -107,83 +117,37 @@ subroutine geopotential_t( & hkk(i) = 0.5_r8 * hkl(i) end do end if - - ! Now compute tv, zm, zi - - do i = 1,ncol - tvfac = 1._r8 + zvir(i,k) * q(i,k,1) - tv = t(i,k) * tvfac - - zm(i,k) = zi(i,k+1) + rog(i,k) * tv * hkk(i) - zi(i,k) = zi(i,k+1) + rog(i,k) * tv * hkl(i) - end do - end do - else - ! - ! For the computation of generalized virtual temperature (equation 16 - ! in Lauritzen et al. (2018); https://doi.org/10.1029/2017MS001257) - ! - ! Compute factor for converting wet to dry mixing ratio (eq.7) - ! - qfac = 1.0_r8 - do idx = dry_air_species_num + 1,thermodynamic_active_species_num - do k = 1,pver - do i = 1,ncol - qfac(i,k) = qfac(i,k)-q(i,k,thermodynamic_active_species_idx(idx)) - end do - end do - end do - qfac = 1.0_r8/qfac - - ! Compute sum of dry water mixing ratios - sum_dry_mixing_ratio = 1.0_r8 - do idx = dry_air_species_num + 1,thermodynamic_active_species_num - do k = 1,pver - do i = 1,ncol - sum_dry_mixing_ratio(i,k) = sum_dry_mixing_ratio(i,k)& - +q(i,k,thermodynamic_active_species_idx(idx))*qfac(i,k) - end do - end do - end do - sum_dry_mixing_ratio(:,:) = 1.0_r8/sum_dry_mixing_ratio(:,:) - ! Compute zi, zm from bottom up. - ! Note, zi(i,k) is the interface above zm(i,k) - do k = pver, 1, -1 - - ! First set hydrostatic elements consistent with dynamics - - ! - ! the outcommented code is left for when/if we will support - ! FV3 and/or FV with condensate loading - ! - -! if ((dycore_is('LR') .or. dycore_is('FV3'))) then -! do i = 1,ncol -! hkl(i) = piln(i,k+1) - piln(i,k) -! hkk(i) = 1._r8 - pint(i,k) * hkl(i) * rpdel(i,k) -! end do -! else!MPAS, SE or EUL - ! - ! For SE : pmid = 0.5*(pint(k+1)+pint(k)) - ! For MPAS : pmid is computed from theta_m, rhodry, etc. - ! - do i = 1,ncol - hkl(i) = pdel(i,k) / pmid(i,k) - hkk(i) = 0.5_r8 * hkl(i) - end do -! end if - + ! Now compute tv, zm, zi - + do i = 1,ncol - tvfac = (1._r8 + (zvir(i,k)+1.0_r8) * q(i,k,1)*qfac(i,k))*sum_dry_mixing_ratio(i,k) + tvfac = 1._r8 + zvir(i,k) * q(i,k,ixq) tv = t(i,k) * tvfac - + zm(i,k) = zi(i,k+1) + rog(i,k) * tv * hkk(i) zi(i,k) = zi(i,k+1) + rog(i,k) * tv * hkl(i) end do end do + else !Using MPAS or SE dycore + + !Determine vertical coordinate type, + !NOTE: Currently the FV (LR) or FV3 dycores + ! do not allow for condensate loading, + ! so for now 'lagrang' will always be FALSE. + if ((dycore_is('LR') .or. dycore_is('FV3'))) then + lagrang = .true. + else + lagrang = .false. + end if + + !Use CCPP version of geopotential_t: + call geopotential_temp_run(pver, lagrang, pver, 1, pverp, 1, & + pcnst, piln(1:ncol,:), pint(1:ncol,:), pmid(1:ncol,:), & + pdel(1:ncol,:), rpdel(1:ncol,:), t(1:ncol,:), & + q(1:ncol,:,ixq), q(1:ncol,:,:), ccpp_const_props, & + rair(1:ncol,:), gravit, zvir(1:ncol,:), zi(1:ncol,:), & + zm(1:ncol,:), ncol, errflg, errmsg) + end if - return end subroutine geopotential_t end module geopotential diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 706b9dcdee..c7d5b2524b 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -786,6 +786,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use phys_grid_ctem, only: phys_grid_ctem_init use cam_budget, only: cam_budget_init + use ccpp_constituent_prop_mod, only: ccpp_const_props_init + ! Input/output arguments type(physics_state), pointer :: phys_state(:) type(physics_tend ), pointer :: phys_tend(:) @@ -978,6 +980,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) end if + ! Initialize CAM CCPP constituent properties array + ! for use in CCPP-ized physics schemes: + call ccpp_const_props_init() + ! Initialize qneg3 and qneg4 call qneg_init() @@ -989,7 +995,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Initialize the budget capability call cam_budget_init() - + ! addfld calls for U, V tendency budget variables that are output in ! tphysac, tphysbc call addfld ( 'UTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by deep convection') @@ -1938,7 +1944,7 @@ subroutine tphysac (ztodt, cam_in, & else ! ! for moist-mixing ratio based dycores - ! + ! ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call ! call set_dry_to_wet(state) @@ -1960,7 +1966,7 @@ subroutine tphysac (ztodt, cam_in, & if (vc_dycore == vc_height.or.vc_dycore == vc_dry_pressure) then ! ! MPAS and SE specific scaling of temperature for enforcing energy consistency - ! (and to make sure that temperature dependent diagnostic tendencies + ! (and to make sure that temperature dependent diagnostic tendencies ! are computed correctly; e.g. dtcore) ! scaling(1:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk) diff --git a/src/physics/cam/ref_pres.F90 b/src/physics/cam/ref_pres.F90 index 742652db11..f0d5994b81 100644 --- a/src/physics/cam/ref_pres.F90 +++ b/src/physics/cam/ref_pres.F90 @@ -1,14 +1,14 @@ module ref_pres !-------------------------------------------------------------------------- -! +! ! Provides access to reference pressures for use by the physics ! parameterizations. The pressures are provided by the dynamical core ! since it determines the grid used by the physics. -! +! ! Note that the init method for this module is called before the init ! method in physpkg; therefore, most physics modules can use these ! reference pressures during their init phases. -! +! !-------------------------------------------------------------------------- use shr_kind_mod, only: r8=>shr_kind_r8 @@ -25,7 +25,7 @@ module ref_pres ! surface pressure ('eta' coordinate) real(r8), protected :: ptop_ref ! Top of model -real(r8), protected :: psurf_ref ! Surface pressure +real(r8), protected :: psurf_ref ! reference pressure ! Number of top levels using pure pressure representation integer, protected :: num_pr_lev diff --git a/src/physics/cam_dev/physpkg.F90 b/src/physics/cam_dev/physpkg.F90 index b288a17177..1c461c9a1c 100644 --- a/src/physics/cam_dev/physpkg.F90 +++ b/src/physics/cam_dev/physpkg.F90 @@ -769,6 +769,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cam_budget, only: cam_budget_init use phys_grid_ctem, only: phys_grid_ctem_init + use ccpp_constituent_prop_mod, only: ccpp_const_props_init + ! Input/output arguments type(physics_state), pointer :: phys_state(:) type(physics_tend ), pointer :: phys_tend(:) @@ -947,6 +949,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) end if + ! Initialize CAM CCPP constituent properties array + ! for use in CCPP-ized physics schemes: + call ccpp_const_props_init() + ! Initialize qneg3 and qneg4 call qneg_init() diff --git a/src/physics/simple/held_suarez_cam.F90 b/src/physics/simple/held_suarez_cam.F90 index 68df115d18..78def4b354 100644 --- a/src/physics/simple/held_suarez_cam.F90 +++ b/src/physics/simple/held_suarez_cam.F90 @@ -1,13 +1,13 @@ module held_suarez_cam - !----------------------------------------------------------------------- - ! + !----------------------------------------------------------------------- + ! ! Purpose: Implement idealized Held-Suarez forcings ! Held, I. M., and M. J. Suarez, 1994: 'A proposal for the ! intercomparison of the dynamical cores of atmospheric general ! circulation models.' ! Bulletin of the Amer. Meteor. Soc., vol. 75, pp. 1825-1830. - ! + ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 @@ -31,24 +31,22 @@ module held_suarez_cam real(r8), parameter :: ka = 1._r8/(86400._r8 * efolda) ! 1./efolding_time for temperature diss. real(r8), parameter :: ks = 1._r8/(86400._r8 * efolds) -!======================================================================= +!======================================================================= contains -!======================================================================= +!======================================================================= - subroutine held_suarez_init(pbuf2d) + subroutine held_suarez_init() use physics_buffer, only: physics_buffer_desc use cam_history, only: addfld, add_default - use physconst, only: cappa, cpair - use ref_pres, only: pref_mid_norm, psurf_ref + use ref_pres, only: psurf_ref use held_suarez_1994, only: held_suarez_1994_init - type(physics_buffer_desc), pointer :: pbuf2d(:,:) - - character(len=512) :: errmsg - integer :: errflg + ! Local variables + character(len=512) :: errmsg + integer :: errflg ! Set model constant values - call held_suarez_1994_init(pver, cappa, cpair, psurf_ref, pref_mid_norm, errmsg, errflg) + call held_suarez_1994_init(psurf_ref, errmsg, errflg) ! This field is added by radiation when full physics is used call addfld('QRS', (/ 'lev' /), 'A', 'K/s', & @@ -57,7 +55,8 @@ subroutine held_suarez_init(pbuf2d) end subroutine held_suarez_init subroutine held_suarez_tend(state, ptend, ztodt) - use air_composition, only: cpairv + use air_composition, only: cappav, cpairv + use ref_pres, only: pref_mid_norm use phys_grid, only: get_rlat_all_p use physics_types, only: physics_state, physics_ptend use physics_types, only: physics_ptend_init @@ -84,8 +83,9 @@ subroutine held_suarez_tend(state, ptend, ztodt) real(r8) :: pmid(pcols,pver) ! mid-point pressure integer :: i, k ! Longitude, level indices - character(len=512) :: errmsg - integer :: errflg + character(len=64) :: scheme_name ! CCPP-required variables (not used in CAM) + character(len=512) :: errmsg + integer :: errflg ! !----------------------------------------------------------------------- @@ -104,8 +104,11 @@ subroutine held_suarez_tend(state, ptend, ztodt) ! initialize individual parameterization tendencies call physics_ptend_init(ptend, state%psetcols, 'held_suarez', ls=.true., lu=.true., lv=.true.) - call held_suarez_1994_run(pver, ncol, clat, state%pmid, & - state%u, state%v, state%t, ptend%u, ptend%v, ptend%s, errmsg, errflg) + call held_suarez_1994_run(pver, ncol, pref_mid_norm, clat, cappav(1:ncol,:,lchnk), & + cpairv(1:ncol,:,lchnk), state%pmid(1:ncol,:), & + state%u(1:ncol,:), state%v(1:ncol,:), state%t(1:ncol,:), & + ptend%u(1:ncol,:), ptend%v(1:ncol,:), ptend%s(1:ncol,:), & + scheme_name, errmsg, errflg) ! Note, we assume that there are no subcolumns in simple physics pmid(:ncol,:) = ptend%s(:ncol, :)/cpairv(:ncol,:,lchnk) diff --git a/src/physics/simple/kessler_cam.F90 b/src/physics/simple/kessler_cam.F90 index d6319962f1..93e961cd97 100644 --- a/src/physics/simple/kessler_cam.F90 +++ b/src/physics/simple/kessler_cam.F90 @@ -39,22 +39,26 @@ end subroutine kessler_register !======================================================================================== - subroutine kessler_cam_init(pbuf2d) + subroutine kessler_cam_init() - use physconst, only: cpair, latvap,pstd, rair, rhoh2o + use physconst, only: latvap, rhoh2o + use ref_pres, only: psurf_ref use constituents, only: cnst_name, cnst_longname, bpcnst, apcnst - use cam_history, only: addfld, add_default, horiz_only + use cam_history, only: addfld, add_default use cam_abortutils, only: endrun - use state_converters, only: pres_to_density_dry_init - use kessler, only: kessler_init - use state_converters, only: calc_exner_init + use kessler, only: kessler_init - integer :: errflg - character(len=512) :: errmsg + ! + !---------------------------Local workspace----------------------------- + ! + integer :: errflg + character(len=512) :: errmsg - type(physics_buffer_desc), pointer :: pbuf2d(:,:) + ! + !----------------------------------------------------------------------- + ! errflg = 0 @@ -68,52 +72,43 @@ subroutine kessler_cam_init(pbuf2d) call add_default(cnst_name(ixrain), 1, ' ') ! Initialize Kessler with CAM physical constants - - if (errflg == 0) then - call kessler_init(cpair, latvap, pstd, rair, rhoh2o, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_cam_init error: Error returned from kessler_init: '//trim(errmsg)) - end if - end if - if (errflg == 0) then - call pres_to_density_dry_init(cpair, rair, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_cam_init error: Error returned from pres_to_density_dry_init: '//trim(errmsg)) - end if - end if - if (errflg == 0) then - call calc_exner_init(errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_cam_init error: Error returned from calc_exner_init: '//trim(errmsg)) - end if - end if + call kessler_init(latvap, psurf_ref, rhoh2o, errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_cam_init error: Error returned from kessler_init: '//trim(errmsg)) + end if end subroutine kessler_cam_init !======================================================================================== subroutine kessler_tend(state, ptend, ztodt, pbuf) - !----------------------------------------------------------------------- - ! + !----------------------------------------------------------------------- + ! ! Purpose: Run Kessler physics (see kessler.F90) - ! + ! !----------------------------------------------------------------------- - use shr_kind_mod, only: SHR_KIND_CM - use physconst, only: cpair, rair, zvir - use physics_types, only: physics_state, physics_ptend - use physics_types, only: physics_ptend_init - use constituents, only: pcnst, cnst_name, cnst_type + use shr_kind_mod, only: SHR_KIND_CM + use ref_pres, only: psurf_ref + use physconst, only: rh2o + use air_composition, only: cpairv, rairv + use physics_types, only: physics_state, physics_ptend + use physics_types, only: physics_ptend_init + use constituents, only: pcnst, cnst_name, cnst_type use kessler, only: kessler_run - use state_converters, only: wet_to_dry_run - use state_converters, only: dry_to_wet_run - use state_converters, only: pres_to_density_dry_run + use state_converters, only: wet_to_dry_water_vapor_run + use state_converters, only: wet_to_dry_cloud_liquid_water_run + use state_converters, only: wet_to_dry_rain_run + use state_converters, only: dry_to_wet_water_vapor_run + use state_converters, only: dry_to_wet_cloud_liquid_water_run + use state_converters, only: dry_to_wet_rain_run use state_converters, only: temp_to_potential_temp_run use state_converters, only: calc_exner_run + use state_converters, only: calc_dry_air_ideal_gas_density_run use state_converters, only: potential_temp_to_temp_run - use cam_abortutils, only: endrun - use cam_history, only: outfld + use cam_abortutils, only: endrun + use cam_history, only: outfld ! arguments @@ -122,28 +117,36 @@ subroutine kessler_tend(state, ptend, ztodt, pbuf) type(physics_ptend), intent(out) :: ptend ! Package tendencies type(physics_buffer_desc), pointer :: pbuf(:) + ! !---------------------------Local workspace----------------------------- ! - integer :: lchnk ! chunk identifier - integer :: ncol ! number of atmospheric columns + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns integer :: lyr_surf integer :: lyr_toa - real(r8) :: rho(pcols,pver) ! Dry air density - real(r8) :: pk(pcols,pver) ! exner func. - real(r8) :: th(pcols,pver) ! Potential temp. - real(r8) :: temp(pcols,pver) ! temperature - real(r8) :: qv(pcols,pver) ! Water vapor - real(r8) :: qc(pcols,pver) ! Cloud water - real(r8) :: qr(pcols,pver) ! Rain water - integer :: k,rk ! vert. indices - logical :: lq(pcnst) ! Calc tendencies? - character(len=SHR_KIND_CM) :: errmsg + real(r8) :: zvirv(pcols,pver) ! ratio of water vapor to dry air constants - 1 + real(r8) :: rho(pcols,pver) ! Dry air density + real(r8) :: pk(pcols,pver) ! exner func. + real(r8) :: th(pcols,pver) ! Potential temp. + real(r8) :: temp(pcols,pver) ! temperature + real(r8) :: qv(pcols,pver) ! Water vapor mixing ratio wrt moist air + real(r8) :: qc(pcols,pver) ! Cloud water mixing ratio wrt moist air + real(r8) :: qr(pcols,pver) ! Rain mixing ratio wrt moist air + real(r8) :: qv_dry(pcols,pver) ! Water vapor mixing ratio wrt dry air + real(r8) :: qc_dry(pcols,pver) ! Cloud water mixing ratio wrt dry air + real(r8) :: qr_dry(pcols,pver) ! Rain mixing ratio wrt dry air - integer :: errflg + integer :: k,rk ! vert. indices + logical :: lq(pcnst) ! Calc tendencies? + character(len=SHR_KIND_CM) :: errmsg ! CCPP physics scheme error message + + integer :: errflg ! CCPP physics scheme error flag - real(r8), pointer :: prec_sed(:) ! total precip from cloud sedimentation - real(r8), pointer :: relhum(:,:) ! relative humidity + character(len=64) :: scheme_name ! CCPP physics scheme name (not used in CAM) + + real(r8), pointer :: prec_sed(:) ! total precip from cloud sedimentation + real(r8), pointer :: relhum(:,:) ! relative humidity integer :: i @@ -169,63 +172,110 @@ subroutine kessler_tend(state, ptend, ztodt, pbuf) call pbuf_get_field(pbuf, relhum_idx, relhum) do k = 1, pver - ! Create temporaries for state variables changed by Kessler routine - temp(:ncol,k) = state%t(:ncol,k) - qv(:ncol,k) = state%q(:ncol,k,1) - qc(:ncol,k) = state%q(:ncol,k,ixcldliq) - qr(:ncol,k) = state%q(:ncol,k,ixrain) + ! Create temporaries for state variables changed by Kessler routine + temp(:ncol,k) = state%t(:ncol,k) + qv(:ncol,k) = state%q(:ncol,k,1) + qc(:ncol,k) = state%q(:ncol,k,ixcldliq) + qr(:ncol,k) = state%q(:ncol,k,ixrain) + + !Also calculate gas constant ratio: + zvirv(:ncol,k) = rh2o/rairv(:ncol,k, lchnk) - 1._r8 end do - - if (errflg == 0) then - call calc_exner_run(ncol, pver, cpair, rair, state%pmid, pk, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from calc_exner_run: '//trim(errmsg)) - end if + + ! Calculate Exner function: + call calc_exner_run(ncol, pver, cpairv(1:ncol,:,lchnk), rairv(1:ncol,:,lchnk), & + psurf_ref, state%pmid(1:ncol,:), pk(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from calc_exner_run: '//trim(errmsg)) + end if + + ! Calculate potential temperature: + call temp_to_potential_temp_run(ncol, pver, temp(1:ncol,:), pk(1:ncol,:), & + th(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from temp_to_potential_temp_run: '//trim(errmsg)) end if - if (errflg == 0) then - call temp_to_potential_temp_run(ncol, pver, temp, pk, th, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from temp_to_potential_temp_run: '//trim(errmsg)) - end if + + ! Calculate density using ideal gas law: + call calc_dry_air_ideal_gas_density_run(ncol, pver, rairv(1:ncol,:,lchnk), & + state%pmiddry(1:ncol,:), temp(1:ncol,:), & + rho(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from pres_to_density_dry_run: '//trim(errmsg)) end if - if (errflg == 0) then - call pres_to_density_dry_run(ncol, pver, state%pmiddry, temp, rho, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from pres_to_density_dry_run: '//trim(errmsg)) - end if + + ! Convert moist air mixing ratios to dry air mixing ratios: + !--------------------------------------------------------- + call wet_to_dry_water_vapor_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), qv(1:ncol,:), & + qv_dry(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from wet_to_dry_water_vapor_run: '//trim(errmsg)) end if - if (errflg == 0) then - call wet_to_dry_run(ncol, pver, state%pdel, state%pdeldry, qv, qc, qr, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from wet_to_dry_run: '//trim(errmsg)) - end if + + call wet_to_dry_cloud_liquid_water_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), qc(1:ncol,:), & + qc_dry(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from wet_to_dry_cloud_liquid_water_run: '//trim(errmsg)) end if - if (errflg == 0) then - call kessler_run(ncol, pver, ztodt, lyr_surf, lyr_toa, & - rho, state%zm, pk, th, qv, qc, qr, prec_sed, relhum, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from kessler_run: '//trim(errmsg)) - end if + + call wet_to_dry_rain_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), qr(1:ncol,:), & + qr_dry(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from wet_to_dry_rain_run: '//trim(errmsg)) end if - if (errflg == 0) then - call potential_temp_to_temp_run(ncol, pver, th, pk, temp, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from potential_temp_to_temp_run: '//trim(errmsg)) - end if + !--------------------------------------------------------- + + ! Run Kessler physics scheme: + call kessler_run(ncol, pver, ztodt, lyr_surf, lyr_toa, cpairv(1:ncol,:,lchnk), & + rairv(1:ncol,:,lchnk), rho(1:ncol,:), state%zm(1:ncol,:), & + pk(1:ncol,:), th(1:ncol,:), qv_dry(1:ncol,:), qc_dry(1:ncol,:), & + qr_dry(1:ncol,:), prec_sed(1:ncol), relhum(1:ncol,:), & + scheme_name, errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from kessler_run: '//trim(errmsg)) end if - if (errflg == 0) then - call dry_to_wet_run(ncol, pver, state%pdel, state%pdeldry, qv, qc, qr, errmsg, errflg) - if (errflg /=0) then - call endrun('kessler_tend error: Error returned from dry_to_wet_run: '//trim(errmsg)) - end if + + ! Calculate air temperature from potential temperature: + call potential_temp_to_temp_run(ncol, pver, th(1:ncol,:), pk(1:ncol,:), & + temp(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from potential_temp_to_temp_run: '//trim(errmsg)) + end if + + ! Convert dry air mixing ratios to moist air mixing ratios: + !--------------------------------------------------------- + call dry_to_wet_water_vapor_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), qv_dry(1:ncol,:), & + qv(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from dry_to_wet_water_vapor_run: '//trim(errmsg)) + end if + + call dry_to_wet_cloud_liquid_water_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), & + qc_dry(1:ncol,:), qc(1:ncol,:), errmsg, & + errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from dry_to_wet_cloud_liquid_water_run: '//trim(errmsg)) + end if + + call dry_to_wet_rain_run(ncol, pver, state%pdel(1:ncol,:), & + state%pdeldry(1:ncol,:), qr_dry(1:ncol,:), & + qr(1:ncol,:), errmsg, errflg) + if (errflg /=0) then + call endrun('kessler_tend error: Error returned from dry_to_wet_rain_run: '//trim(errmsg)) end if + !--------------------------------------------------------- ! Back out tendencies from updated fields do k = 1, pver - ptend%s(:ncol,k) = (th(:ncol,k)*pk(:ncol,k) - state%t(:ncol,k)) * cpair / ztodt - ptend%q(:ncol,k,1) = (qv(:ncol,k) - state%q(:ncol,k,1)) / ztodt - ptend%q(:ncol,k,ixcldliq) = (qc(:ncol,k) - state%q(:ncol,k,ixcldliq)) / ztodt - ptend%q(:ncol,k,ixrain) = (qr(:ncol,k) - state%q(:ncol,k,ixrain)) / ztodt + ptend%s(:ncol,k) = (th(:ncol,k)*pk(:ncol,k) - state%t(:ncol,k)) * cpairv(:ncol,k,lchnk) / ztodt + ptend%q(:ncol,k,1) = (qv(:ncol,k) - state%q(:ncol,k,1)) / ztodt + ptend%q(:ncol,k,ixcldliq) = (qc(:ncol,k) - state%q(:ncol,k,ixcldliq)) / ztodt + ptend%q(:ncol,k,ixrain) = (qr(:ncol,k) - state%q(:ncol,k,ixrain)) / ztodt end do ! Output liquid tracers diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 31c41def58..a296fd2fdb 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -209,6 +209,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cam_snapshot, only: cam_snapshot_init use cam_budget, only: cam_budget_init + use ccpp_constituent_prop_mod, only: ccpp_const_props_init + ! Input/output arguments type(physics_state), pointer :: phys_state(:) type(physics_tend ), pointer :: phys_tend(:) @@ -259,9 +261,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) end if if (ideal_phys) then - call held_suarez_init(pbuf2d) + call held_suarez_init() else if (kessler_phys) then - call kessler_cam_init(pbuf2d) + call kessler_cam_init() else if (tj2016_phys) then call thatcher_jablonowski_init(pbuf2d) else if (frierson_phys) then @@ -277,6 +279,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call chem_init(phys_state,pbuf2d) end if + ! Initialize CAM CCPP constituent properties array + ! for use in CCPP-ized physics schemes: + call ccpp_const_props_init() + ! Initialize qneg3 and qneg4 call qneg_init() diff --git a/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 new file mode 100644 index 0000000000..edbc4d59e9 --- /dev/null +++ b/src/utils/cam_ccpp/ccpp_constituent_prop_mod.F90 @@ -0,0 +1,190 @@ +! This module is the CAM version of the CCPP generated module of the same name +module ccpp_constituent_prop_mod + + implicit none + private + + ! Define CAM version of constituent properties mod + type, public :: ccpp_constituent_prop_ptr_t + logical, private :: thermo_active = .false. + logical, private :: water_species = .false. + contains + procedure :: standard_name => ccp_get_standard_name + procedure :: is_thermo_active => ccp_is_thermo_active + procedure :: is_water_species => ccp_is_water_species + procedure :: set_thermo_active => ccp_set_thermo_active + procedure :: set_water_species => ccp_set_water_species + end type ccpp_constituent_prop_ptr_t + + ! CCPP properties init routine + public :: ccpp_const_props_init + + ! Public properties DDT variable: + type(ccpp_constituent_prop_ptr_t), allocatable, public :: ccpp_const_props(:) + +contains + +!+++++++++++++++++++++++++++++++++++++++ +!CCPP constituent properties DDT methods +!+++++++++++++++++++++++++++++++++++++++ + + subroutine ccp_get_standard_name(this, std_name, errcode, errmsg) + ! Return this constituent's standard name + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(in) :: this + character(len=*), intent(out) :: std_name + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + std_name = 'Not Used!' + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = 'Still Not Used!' + end if + + end subroutine ccp_get_standard_name + + !------ + + subroutine ccp_is_thermo_active(this, val_out, errcode, errmsg) + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(in) :: this + logical, intent(out) :: val_out + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Pass back thermo active property: + val_out = this%thermo_active + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_is_thermo_active + + !------ + + subroutine ccp_is_water_species(this, val_out, errcode, errmsg) + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(in) :: this + logical, intent(out) :: val_out + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Pass back water species property: + val_out = this%water_species + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_is_water_species + + !------ + + subroutine ccp_set_thermo_active(this, thermo_flag, errcode, errmsg) + ! Set whether this constituent is thermodynamically active, which + ! means that certain physics schemes will use this constitutent + ! when calculating thermodynamic quantities (e.g. enthalpy). + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(inout) :: this + logical, intent(in) :: thermo_flag + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Set thermodynamically active flag for this constituent: + this%thermo_active = thermo_flag + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_set_thermo_active + + !------ + + subroutine ccp_set_water_species(this, water_flag, errcode, errmsg) + ! Set whether this constituent is a water species, which means + ! that this constituent represents a particular phase or type + ! of water in the atmosphere. + + ! Dummy arguments + class(ccpp_constituent_prop_ptr_t), intent(inout) :: this + logical, intent(in) :: water_flag + integer, optional, intent(out) :: errcode + character(len=*), optional, intent(out) :: errmsg + + ! Set thermodynamically active flag for this constituent: + this%water_species = water_flag + + ! Provide err values if requested: + if(present(errcode)) then + errcode = 0 + end if + if(present(errmsg)) then + errmsg = '' + end if + + end subroutine ccp_set_water_species + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!CAM-equivalent CCPP constituents initialization routine +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +subroutine ccpp_const_props_init() + + ! Use statements: + use constituents, only: pcnst + use cam_abortutils, only: handle_allocate_error + use air_composition, only: dry_air_species_num + use air_composition, only: thermodynamic_active_species_idx + + ! Local variables: + integer :: ierr + integer :: m + + character(len=*), parameter :: subname = 'ccpp_const_prop_init:' + + ! Allocate constituents object: + allocate(ccpp_const_props(pcnst), stat=ierr) + + ! Check if allocation succeeded: + call handle_allocate_error(ierr, subname, 'ccpp_const_props(pcnst)') + + ! Set "thermo_active" property: + do m = 1,pcnst + if(any(thermodynamic_active_species_idx == m)) then + call ccpp_const_props(m)%set_thermo_active(.true.) + end if + end do + + ! Set "water_species" property: + do m=1,pcnst + if(any(thermodynamic_active_species_idx(dry_air_species_num+1:) == m)) then + call ccpp_const_props(m)%set_water_species(.true.) + end if + end do + +end subroutine ccpp_const_props_init + +end module ccpp_constituent_prop_mod diff --git a/src/utils/ccpp_kinds.F90 b/src/utils/cam_ccpp/ccpp_kinds.F90 similarity index 65% rename from src/utils/ccpp_kinds.F90 rename to src/utils/cam_ccpp/ccpp_kinds.F90 index 505001a625..c95c2b08c3 100644 --- a/src/utils/ccpp_kinds.F90 +++ b/src/utils/cam_ccpp/ccpp_kinds.F90 @@ -1,4 +1,4 @@ -! This module is a placeholder for the CCPP generated module of the same name +! This module is the CAM version of the CCPP generated module of the same name module ccpp_kinds use shr_kind_mod, only: kind_phys => shr_kind_r8