From 2772bb4aa280d4c4dc4d5253fa2fc672dde09b48 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Sep 2024 20:04:06 -0400 Subject: [PATCH 01/20] Port check_energy_chng, check_energy_timestep_init to CCPP-ized subroutines --- bld/configure | 1 + src/dynamics/eul/dp_coupling.F90 | 4 +- src/dynamics/fv/dp_coupling.F90 | 4 +- src/dynamics/mpas/dp_coupling.F90 | 4 +- src/dynamics/se/dp_coupling.F90 | 4 +- src/physics/cam/check_energy.F90 | 302 +-------------------------- src/physics/cam/check_energy_cam.F90 | 211 +++++++++++++++++++ src/physics/cam/physpkg.F90 | 5 +- src/physics/cam/zm_conv_intr.F90 | 1 - src/physics/cam7/physpkg.F90 | 6 +- 10 files changed, 232 insertions(+), 310 deletions(-) create mode 100644 src/physics/cam/check_energy_cam.F90 diff --git a/bld/configure b/bld/configure index 86c6c4e65b..4b898b9d88 100755 --- a/bld/configure +++ b/bld/configure @@ -2340,6 +2340,7 @@ sub write_filepath print $fh "$camsrcdir/src/atmos_phys/schemes/tropopause_find\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/zhang_mcfarlane\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/dry_adiabatic_adjust\n"; + print $fh "$camsrcdir/src/atmos_phys/schemes/check_energy\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/utilities\n"; # Dynamics package and test utilities diff --git a/src/dynamics/eul/dp_coupling.F90 b/src/dynamics/eul/dp_coupling.F90 index bc900e2d0e..8e3eb64a21 100644 --- a/src/dynamics/eul/dp_coupling.F90 +++ b/src/dynamics/eul/dp_coupling.F90 @@ -14,7 +14,7 @@ module dp_coupling use physconst, only: cpair, gravit, rair, zvir use air_composition, only: rairv use geopotential, only: geopotential_t - use check_energy, only: check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_timestep_init #if (defined SPMD) use spmd_dyn, only: buf1, buf1win, buf2, buf2win, & spmdbuf_siz, local_dp_map, & @@ -291,7 +291,7 @@ subroutine d_p_coupling(ps, t3, u3, v3, q3, & ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk ) + call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk ) end do diff --git a/src/dynamics/fv/dp_coupling.F90 b/src/dynamics/fv/dp_coupling.F90 index fc02821471..af2d13b12e 100644 --- a/src/dynamics/fv/dp_coupling.F90 +++ b/src/dynamics/fv/dp_coupling.F90 @@ -12,7 +12,7 @@ module dp_coupling use physconst, only: gravit, zvir use air_composition, only: cpairv, rairv use geopotential, only: geopotential_t - use check_energy, only: check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_timestep_init use dynamics_vars, only: T_FVDYCORE_GRID, t_fvdycore_state use dyn_internal_state,only: get_dyn_state use dyn_comp, only: dyn_import_t, dyn_export_t, fv_print_dpcoup_warn @@ -621,7 +621,7 @@ subroutine d_p_coupling(grid, phys_state, phys_tend, pbuf2d, dyn_out) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do call t_stopf('derived_fields') diff --git a/src/dynamics/mpas/dp_coupling.F90 b/src/dynamics/mpas/dp_coupling.F90 index 10d75b4b8c..c312f1e2f9 100644 --- a/src/dynamics/mpas/dp_coupling.F90 +++ b/src/dynamics/mpas/dp_coupling.F90 @@ -413,7 +413,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d) ! MPAS prognostic fields. use geopotential, only: geopotential_t - use check_energy, only: check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_timestep_init use shr_vmath_mod, only: shr_vmath_log use phys_control, only: waccmx_is use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update @@ -564,7 +564,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 919b7f3510..9de31132a2 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -546,7 +546,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) 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 check_energy_cam,only: check_energy_cam_timestep_init use hycoef, only: hyai, ps0 use shr_vmath_mod, only: shr_vmath_log use qneg_module, only: qneg3 @@ -707,7 +707,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) cpairv(1:ncol,:,lchnk), errflg, errmsg) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do ! lchnk diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 9c569387e0..5e1df3a967 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -26,10 +26,9 @@ module check_energy use gmean_mod, only: gmean use physconst, only: gravit, rga, latvap, latice, cpair, rair - use air_composition, only: cpairv, rairv, cp_or_cv_dycore + use air_composition, only: cpairv, cp_or_cv_dycore use physics_types, only: physics_state, physics_tend, physics_ptend, physics_ptend_init use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind - use time_manager, only: is_first_step use cam_logfile, only: iulog use scamMod, only: single_column, use_camiop, heat_glob_scm use cam_history, only: outfld, write_camiop @@ -45,8 +44,6 @@ module check_energy public :: check_energy_register ! register fields in physics buffer public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean public :: check_energy_init ! initialization of module - public :: check_energy_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes - public :: check_energy_chng ! check changes in integrals against cumulative boundary fluxes public :: check_energy_gmean ! global means of physics input and output total energy public :: check_energy_fix ! add global mean energy difference as a heating public :: check_tracers_init ! initialize tracer integrals and cumulative boundary fluxes @@ -67,11 +64,11 @@ module check_energy ! Physics buffer indices - integer :: teout_idx = 0 ! teout index in physics buffer - integer :: dtcore_idx = 0 ! dtcore index in physics buffer - integer :: dqcore_idx = 0 ! dqcore index in physics buffer - integer :: ducore_idx = 0 ! ducore index in physics buffer - integer :: dvcore_idx = 0 ! dvcore index in physics buffer + integer, public :: teout_idx = 0 ! teout index in physics buffer + integer, public :: dtcore_idx = 0 ! dtcore index in physics buffer + integer, public :: dqcore_idx = 0 ! dqcore index in physics buffer + integer, public :: ducore_idx = 0 ! ducore index in physics buffer + integer, public :: dvcore_idx = 0 ! dvcore index in physics buffer type check_tracers_data real(r8) :: tracer(pcols,pcnst) ! initial vertically integrated total (kinetic + static) energy @@ -219,293 +216,6 @@ end subroutine check_energy_init !=============================================================================== - subroutine check_energy_timestep_init(state, tend, pbuf, col_type) - use cam_thermo, only: get_hydrostatic_energy - use physics_buffer, only: physics_buffer_desc, pbuf_set_field - use cam_abortutils, only: endrun - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure - use physics_types, only: phys_te_idx, dyn_te_idx -!----------------------------------------------------------------------- -! Compute initial values of energy and water integrals, -! zero cumulative tendencies -!----------------------------------------------------------------------- -!------------------------------Arguments-------------------------------- - - type(physics_state), intent(inout) :: state - type(physics_tend ), intent(inout) :: tend - type(physics_buffer_desc), pointer :: pbuf(:) - integer, optional :: col_type ! Flag inidicating whether using grid or subcolumns -!---------------------------Local storage------------------------------- - real(r8) :: cp_or_cv(state%psetcols,pver) - integer lchnk ! chunk identifier - integer ncol ! number of atmospheric columns -!----------------------------------------------------------------------- - - lchnk = state%lchnk - ncol = state%ncol - - ! cp_or_cv needs to be allocated to a size which matches state and ptend - ! If psetcols == pcols, cpairv is the correct size and just copy into cp_or_cv - ! If psetcols > pcols and all cpairv match cpair, then assign the constant cpair - - if (state%psetcols == pcols) then - cp_or_cv(:,:) = cpairv(:,:,lchnk) - else if (state%psetcols > pcols .and. all(cpairv(:,:,lchnk) == cpair)) then - cp_or_cv(1:ncol,:) = cpair - else - call endrun('check_energy_timestep_init: cpairv is not allowed to vary when subcolumns are turned on') - end if - ! - ! CAM physics total energy - ! - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), & - vc_physics, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol),& - te = state%te_ini(1:ncol,phys_te_idx), H2O = state%tw_ini(1:ncol)) - ! - ! Dynamical core total energy - ! - state%temp_ini(:ncol,:) = state%T(:ncol,:) - state%z_ini(:ncol,:) = state%zm(:ncol,:) - if (vc_dycore == vc_height) then - ! - ! MPAS specific hydrostatic energy computation (internal energy) - ! - if (state%psetcols == pcols) then - cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - else - cp_or_cv(:ncol,:) = cpair-rair - endif - - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), & - vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & - z_mid = state%z_ini(1:ncol,:), & - te = state%te_ini(1:ncol,dyn_te_idx), H2O = state%tw_ini(1:ncol)) - else if (vc_dycore == vc_dry_pressure) then - ! - ! SE specific hydrostatic energy (enthalpy) - ! - if (state%psetcols == pcols) then - cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - else - cp_or_cv(:ncol,:) = cpair - endif - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), & - vc_dry_pressure, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & - te = state%te_ini(1:ncol,dyn_te_idx), H2O = state%tw_ini(1:ncol)) - else - ! - ! dycore energy is the same as physics - ! - state%te_ini(1:ncol,dyn_te_idx) = state%te_ini(1:ncol,phys_te_idx) - end if - state%te_cur(:ncol,:) = state%te_ini(:ncol,:) - state%tw_cur(:ncol) = state%tw_ini(:ncol) - -! zero cummulative boundary fluxes - tend%te_tnd(:ncol) = 0._r8 - tend%tw_tnd(:ncol) = 0._r8 - - state%count = 0 - -! initialize physics buffer - if (is_first_step()) then - call pbuf_set_field(pbuf, teout_idx, state%te_ini(:,dyn_te_idx), col_type=col_type) - end if - - end subroutine check_energy_timestep_init - -!=============================================================================== - - subroutine check_energy_chng(state, tend, name, nstep, ztodt, & - flx_vap, flx_cnd, flx_ice, flx_sen) - use cam_thermo, only: get_hydrostatic_energy - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure - use cam_abortutils, only: endrun - use physics_types, only: phys_te_idx, dyn_te_idx -!----------------------------------------------------------------------- -! Check that the energy and water change matches the boundary fluxes -!----------------------------------------------------------------------- -!------------------------------Arguments-------------------------------- - - type(physics_state) , intent(inout) :: state - type(physics_tend ) , intent(inout) :: tend - character*(*),intent(in) :: name ! parameterization name for fluxes - integer , intent(in ) :: nstep ! current timestep number - real(r8), intent(in ) :: ztodt ! 2 delta t (model time increment) - real(r8), intent(in ) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) - real(r8), intent(in ) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?) - real(r8), intent(in ) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?) - real(r8), intent(in ) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2) - -!******************** BAB ****************************************************** -!******* Note that the precip and ice fluxes are in precip units (m/s). ******** -!******* I would prefer to have kg/m2/s. ******** -!******* I would also prefer liquid (not total) and ice fluxes ******** -!******************************************************************************* - -!---------------------------Local storage------------------------------- - - real(r8) :: te_xpd(state%ncol) ! expected value (f0 + dt*boundary_flux) - real(r8) :: te_dif(state%ncol) ! energy of input state - original energy - real(r8) :: te_tnd(state%ncol) ! tendency from last process - real(r8) :: te_rer(state%ncol) ! relative error in energy column - - real(r8) :: tw_xpd(state%ncol) ! expected value (w0 + dt*boundary_flux) - real(r8) :: tw_dif(state%ncol) ! tw_inp - original water - real(r8) :: tw_tnd(state%ncol) ! tendency from last process - real(r8) :: tw_rer(state%ncol) ! relative error in water column - - real(r8) :: te(state%ncol) ! vertical integral of total energy - real(r8) :: tw(state%ncol) ! vertical integral of total water - real(r8) :: cp_or_cv(state%psetcols,pver) ! cp or cv depending on vcoord - real(r8) :: scaling(state%psetcols,pver) ! scaling for conversion of temperature increment - real(r8) :: temp(state%ncol,pver) ! temperature - - real(r8) :: se(state%ncol) ! enthalpy or internal energy (J/m2) - real(r8) :: po(state%ncol) ! surface potential or potential energy (J/m2) - real(r8) :: ke(state%ncol) ! kinetic energy (J/m2) - real(r8) :: wv(state%ncol) ! column integrated vapor (kg/m2) - real(r8) :: liq(state%ncol) ! column integrated liquid (kg/m2) - real(r8) :: ice(state%ncol) ! column integrated ice (kg/m2) - - integer lchnk ! chunk identifier - integer ncol ! number of atmospheric columns - integer i ! column index -!----------------------------------------------------------------------- - - lchnk = state%lchnk - ncol = state%ncol - - ! If psetcols == pcols, cpairv is the correct size and just copy into cp_or_cv - ! If psetcols > pcols and all cpairv match cpair, then assign the constant cpair - - if (state%psetcols == pcols) then - cp_or_cv(:,:) = cpairv(:,:,lchnk) - else if (state%psetcols > pcols .and. all(cpairv(:,:,:) == cpair)) then - cp_or_cv(:,:) = cpair - else - call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') - end if - - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver), & - vc_physics, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & - te = te(1:ncol), H2O = tw(1:ncol), se=se(1:ncol),po=po(1:ncol), & - ke=ke(1:ncol),wv=wv(1:ncol),liq=liq(1:ncol),ice=ice(1:ncol)) - ! compute expected values and tendencies - do i = 1, ncol - ! change in static energy and total water - te_dif(i) = te(i) - state%te_cur(i,phys_te_idx) - tw_dif(i) = tw(i) - state%tw_cur(i) - - ! expected tendencies from boundary fluxes for last process - te_tnd(i) = flx_vap(i)*(latvap+latice) - (flx_cnd(i) - flx_ice(i))*1000._r8*latice + flx_sen(i) - tw_tnd(i) = flx_vap(i) - flx_cnd(i) *1000._r8 - - ! cummulative tendencies from boundary fluxes - tend%te_tnd(i) = tend%te_tnd(i) + te_tnd(i) - tend%tw_tnd(i) = tend%tw_tnd(i) + tw_tnd(i) - - ! expected new values from previous state plus boundary fluxes - te_xpd(i) = state%te_cur(i,phys_te_idx) + te_tnd(i)*ztodt - tw_xpd(i) = state%tw_cur(i) + tw_tnd(i)*ztodt - - ! relative error, expected value - input state / previous state - te_rer(i) = (te_xpd(i) - te(i)) / state%te_cur(i,phys_te_idx) - end do - - ! relative error for total water (allow for dry atmosphere) - tw_rer = 0._r8 - where (state%tw_cur(:ncol) > 0._r8) - tw_rer(:ncol) = (tw_xpd(:ncol) - tw(:ncol)) / state%tw_cur(:ncol) - end where - - ! error checking - if (print_energy_errors) then - if (any(abs(te_rer(1:ncol)) > 1.E-14_r8 .or. abs(tw_rer(1:ncol)) > 1.E-10_r8)) then - do i = 1, ncol - ! the relative error threshold for the water budget has been reduced to 1.e-10 - ! to avoid messages generated by QNEG3 calls - ! PJR- change to identify if error in energy or water - if (abs(te_rer(i)) > 1.E-14_r8 ) then - state%count = state%count + 1 - write(iulog,*) "significant energy conservation error after ", name, & - " count", state%count, " nstep", nstep, "chunk", lchnk, "col", i - write(iulog,*) te(i),te_xpd(i),te_dif(i),tend%te_tnd(i)*ztodt, & - te_tnd(i)*ztodt,te_rer(i) - endif - if ( abs(tw_rer(i)) > 1.E-10_r8) then - state%count = state%count + 1 - write(iulog,*) "significant water conservation error after ", name, & - " count", state%count, " nstep", nstep, "chunk", lchnk, "col", i - write(iulog,*) tw(i),tw_xpd(i),tw_dif(i),tend%tw_tnd(i)*ztodt, & - tw_tnd(i)*ztodt,tw_rer(i) - end if - end do - end if - end if - - ! copy new value to state - - do i = 1, ncol - state%te_cur(i,phys_te_idx) = te(i) - state%tw_cur(i) = tw(i) - end do - - ! - ! Dynamical core total energy - ! - if (vc_dycore == vc_height) then - ! - ! compute cv if vertical coordinate is height: cv = cp - R - ! - ! Note: cp_or_cv set above for pressure coordinate - if (state%psetcols == pcols) then - cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - else - cp_or_cv(:ncol,:) = cpair-rair - endif - scaling(:,:) = cpairv(:,:,lchnk)/cp_or_cv(:,:) !cp/cv scaling - temp(1:ncol,:) = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:)) - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver), & - vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & - z_mid = state%z_ini(1:ncol,:), & - te = state%te_cur(1:ncol,dyn_te_idx), H2O = state%tw_cur(1:ncol)) - else if (vc_dycore == vc_dry_pressure) then - ! - ! SE specific hydrostatic energy - ! - if (state%psetcols == pcols) then - cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - scaling(:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk) - else - cp_or_cv(:ncol,:) = cpair - scaling(:ncol,:) = 1.0_r8 - endif - ! - ! enthalpy scaling for energy consistency - ! - temp(1:ncol,:) = state%temp_ini(1:ncol,:)+scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:)) - call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & - state%pdel(1:ncol,1:pver), cp_or_cv(1:ncol,1:pver), & - state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), temp(1:ncol,1:pver), & - vc_dry_pressure, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & - te = state%te_cur(1:ncol,dyn_te_idx), H2O = state%tw_cur(1:ncol)) - else - state%te_cur(1:ncol,dyn_te_idx) = te(1:ncol) - end if - end subroutine check_energy_chng - - subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 new file mode 100644 index 0000000000..2bc48d5327 --- /dev/null +++ b/src/physics/cam/check_energy_cam.F90 @@ -0,0 +1,211 @@ +! Shim for CCPP-ized check_energy routines in CAM +module check_energy_cam + use shr_kind_mod, only: r8 => shr_kind_r8 + use ppgrid, only: pcols, pver + use constituents, only: pcnst + use physics_types, only: physics_state, physics_tend + + implicit none + private + + public :: check_energy_cam_chng ! check changes in integrals against cumulative boundary fluxes + public :: check_energy_cam_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes + +contains + ! Compute initial values of energy and water integrals, + ! zero cumulative tendencies + subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) + use physics_buffer, only: physics_buffer_desc, pbuf_set_field + use cam_abortutils, only: endrun + use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure + use physics_types, only: phys_te_idx, dyn_te_idx + use time_manager, only: is_first_step + use physconst, only: cpair, rair + use air_composition, only: cpairv, cp_or_cv_dycore + + ! To remove, pbuf indices in check_energy + use check_energy, only: teout_idx + + ! CCPP-ized subroutine + use check_energy_chng, only: check_energy_chng_timestep_init + + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + type(physics_buffer_desc), pointer :: pbuf(:) + integer, optional :: col_type ! Flag indicating whether using grid or subcolumns + + real(r8) :: local_cp_phys(state%psetcols,pver) + real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + character(len=512) :: errmsg + integer :: errflg + + lchnk = state%lchnk + ncol = state%ncol + + ! The code below is split into not-subcolumns and subcolumns code, as there is different handling of the + ! cp passed into the hydrostatic energy call. CAM-SIMA does not support subcolumns, so we keep this special + ! handling inside this shim module. (hplin, 9/9/24) + if(state%psetcols == pcols) then + ! No subcolumns + local_cp_phys(:,:) = cpairv(:,:,lchnk) + local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) + else if (state%psetcols > pcols) then + ! Subcolumns code + ! Subcolumns specific error handling + if(.not. all(cpairv(:,:,lchnk) == cpair)) then + call endrun('check_energy_timestep_init: cpairv is not allowed to vary when subcolumns are turned on') + endif + + local_cp_phys(1:ncol,:) = cpair + + if (vc_dycore == vc_height) then + ! MPAS specific hydrostatic energy computation (internal energy) + local_cp_or_cv_dycore(:ncol,:) = cpair-rair + else if(vc_dycore == vc_dry_pressure) then + ! SE specific hydrostatic energy (enthalpy) + local_cp_or_cv_dycore(:ncol,:) = cpair + else + ! cp_or_cv is not used in the underlying subroutine, zero it out to be sure + local_cp_or_cv_dycore(:ncol,:) = 0.0_r8 + endif + end if + + ! Call CCPP-ized underlying subroutine. + call check_energy_chng_timestep_init( & + ncol = ncol, & + pver = pver, & + pcnst = pcnst, & + q = state%q(1:ncol,1:pver,1:pcnst), & + pdel = state%pdel(1:ncol,1:pver), & + u = state%u(1:ncol,1:pver), & + v = state%v(1:ncol,1:pver), & + T = state%T(1:ncol,1:pver), & + pintdry = state%pintdry(1:ncol,1:pver), & + phis = state%phis(1:ncol), & + zm = state%zm(1:ncol,:), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & + cp_phys = local_cp_phys(1:ncol,:), & + cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & + te_ini_phys = state%te_ini(1:ncol,phys_te_idx), & + te_ini_dyn = state%te_ini(1:ncol,dyn_te_idx), & + tw_ini = state%tw_ini(1:ncol), & + te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & + te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & + tw_cur = state%tw_cur(1:ncol), & + tend_te_tnd = tend%te_tnd(1:ncol), & + tend_tw_tnd = tend%tw_tnd(1:ncol), & + count = state%count, & + vc_physics = vc_physics, & + vc_dycore = vc_dycore, & + errmsg = errmsg, & + errflg = errflg & + ) + + ! initialize physics buffer + if (is_first_step()) then + call pbuf_set_field(pbuf, teout_idx, state%te_ini(:,dyn_te_idx), col_type=col_type) + end if + + end subroutine check_energy_cam_timestep_init + + ! Check that the energy and water change matches the boundary fluxes + subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & + flx_vap, flx_cnd, flx_ice, flx_sen) + use cam_thermo, only: get_hydrostatic_energy + use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure + use cam_abortutils, only: endrun + use physics_types, only: phys_te_idx, dyn_te_idx + use physconst, only: cpair, rair + use air_composition, only: cpairv, cp_or_cv_dycore + + ! CCPP-ized subroutine + use check_energy_chng, only: check_energy_chng_run + + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + character*(*),intent(in) :: name ! parameterization name for fluxes + integer , intent(in) :: nstep ! current timestep number + real(r8), intent(in) :: ztodt ! 2 delta t (model time increment) + real(r8), intent(in) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) + real(r8), intent(in) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?) + real(r8), intent(in) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?) + real(r8), intent(in) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2) + + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: i ! column index + real(r8) :: local_cp_phys(state%psetcols,pver) + real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) + real(r8) :: scaling_dycore(state%ncol,pver) + character(len=512) :: errmsg + integer :: errflg + + lchnk = state%lchnk + ncol = state%ncol + + if(state%psetcols == pcols) then + ! No subcolumns + local_cp_phys(:,:) = cpairv(:,:,lchnk) + local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) + + scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + elseif(state%psetcols > pcols) then + ! Subcolumns + if(.not. all(cpairv(:,:,:) == cpair)) then + call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') + endif + + local_cp_phys(:,:) = cpair + + ! Note: cp_or_cv set above for pressure coordinate + if (vc_dycore == vc_height) then + ! compute cv if vertical coordinate is height: cv = cp - R + local_cp_or_cv_dycore(:ncol,:) = cpair-rair + scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + else if (vc_dycore == vc_dry_pressure) then + ! SE specific hydrostatic energy + local_cp_or_cv_dycore(:ncol,:) = cpair + scaling_dycore(:ncol,:) = 1.0_r8 + else + ! Moist pressure... use phys formula + local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:,lchnk) + scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + end if + endif + + ! Call CCPP-ized underlying subroutine. + call check_energy_chng_run( & + ncol = ncol, & + pver = pver, & + pcnst = pcnst, & + q = state%q(1:ncol,1:pver,1:pcnst), & + pdel = state%pdel(1:ncol,1:pver), & + u = state%u(1:ncol,1:pver), & + v = state%v(1:ncol,1:pver), & + T = state%T(1:ncol,1:pver), & + pintdry = state%pintdry(1:ncol,1:pver), & + phis = state%phis(1:ncol), & + zm = state%zm(1:ncol,:), & + cp_phys = local_cp_phys(1:ncol,:), & + cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & + scaling_dycore = scaling_dycore(1:ncol,:), & + te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & + te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & + tw_cur = state%tw_cur(1:ncol), & + tend_te_tnd = tend%te_tnd(1:ncol), & + tend_tw_tnd = tend%tw_tnd(1:ncol), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & + count = state%count, & + ztodt = ztodt, & + vc_physics = vc_physics, & + vc_dycore = vc_dycore, & + errmsg = errmsg, & + errflg = errflg & + ) + + end subroutine check_energy_cam_chng +end module check_energy_cam \ No newline at end of file diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index ba36670ce8..227274355e 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -2044,7 +2044,8 @@ subroutine tphysbc (ztodt, state, & use convect_deep, only: convect_deep_tend, convect_deep_tend_2, deep_scheme_does_scav_trans use time_manager, only: is_first_step, get_nstep use convect_shallow, only: convect_shallow_tend - use check_energy, only: check_energy_chng, check_energy_fix, check_energy_timestep_init + use check_energy_cam,only: check_energy_cam_timestep_init, check_energy_cam_chng + use check_energy, only: check_energy_fix use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use dycore, only: dycore_is @@ -2630,7 +2631,7 @@ subroutine tphysbc (ztodt, state, & call subcol_gen(state, tend, state_sc, tend_sc, pbuf) !Initialize check energy for subcolumns - call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) + call check_energy_cam_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) end if if (trim(cam_take_snapshot_before) == "microp_section") then diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 4113c33a4b..5db6d1bc03 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -393,7 +393,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & use time_manager, only: get_nstep, is_first_step use physics_buffer, only : pbuf_get_field, physics_buffer_desc, pbuf_old_tim_idx use constituents, only: pcnst, cnst_get_ind, cnst_is_convtran1 - use check_energy, only: check_energy_chng use physconst, only: gravit, latice, latvap, tmelt, cpwv, cpliq, rh2o use phys_control, only: cam_physpkg_is diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index 83d03c46d1..ad5bc9913e 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -1366,7 +1366,8 @@ subroutine tphysac (ztodt, cam_in, & use aoa_tracers, only: aoa_tracers_timestep_tend use physconst, only: rhoh2o use aero_model, only: aero_model_drydep - use check_energy, only: check_energy_chng, tot_energy_phys + use check_energy_cam, only: check_energy_cam_timestep_init, check_energy_cam_chng + use check_energy, only: tot_energy_phys use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use time_manager, only: get_nstep use cam_abortutils, only: endrun @@ -1405,7 +1406,6 @@ subroutine tphysac (ztodt, cam_in, & use aero_model, only: aero_model_wetdep use aero_wetdep_cam, only: wetdep_lq use physics_buffer, only: col_type_subcol - use check_energy, only: check_energy_timestep_init use carma_intr, only: carma_wetdep_tend, carma_timestep_tend, carma_emission_tend use carma_flags_mod, only: carma_do_aerosol, carma_do_emission, carma_do_detrain use carma_flags_mod, only: carma_do_cldice, carma_do_cldliq, carma_do_wetdep @@ -1775,7 +1775,7 @@ subroutine tphysac (ztodt, cam_in, & call subcol_gen(state, tend, state_sc, tend_sc, pbuf) !Initialize check energy for subcolumns - call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) + call check_energy_cam_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) end if if (trim(cam_take_snapshot_before) == "microp_section") then From 50afaa6917eb7a40f4f344a8943be30a71c0240f Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Sep 2024 23:56:53 -0400 Subject: [PATCH 02/20] Update physical constants to pass into subroutine --- src/physics/cam/check_energy_cam.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 index 2bc48d5327..4f686b598d 100644 --- a/src/physics/cam/check_energy_cam.F90 +++ b/src/physics/cam/check_energy_cam.F90 @@ -118,7 +118,7 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure use cam_abortutils, only: endrun use physics_types, only: phys_te_idx, dyn_te_idx - use physconst, only: cpair, rair + use physconst, only: cpair, rair, latice, latvap use air_composition, only: cpairv, cp_or_cv_dycore ! CCPP-ized subroutine @@ -201,6 +201,8 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & z_ini = state%z_ini(:ncol,:), & count = state%count, & ztodt = ztodt, & + latice = latice, & + latvap = latvap, & vc_physics = vc_physics, & vc_dycore = vc_dycore, & errmsg = errmsg, & From 80d5b86fb3d2018c978642c21902b31af6c56716 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 10 Sep 2024 00:11:10 -0400 Subject: [PATCH 03/20] Build updates; now use check_energy_cam_chng --- src/physics/cam/check_energy_cam.F90 | 7 ++++- src/physics/cam/physpkg.F90 | 42 ++++++++++++++-------------- src/physics/cam7/physpkg.F90 | 39 +++++++++++++------------- src/physics/simple/physpkg.F90 | 11 ++++---- src/physics/spcam/crm_physics.F90 | 4 +-- src/physics/spcam/spcam_drivers.F90 | 7 +++-- 6 files changed, 59 insertions(+), 51 deletions(-) diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 index 4f686b598d..d5ae76ca48 100644 --- a/src/physics/cam/check_energy_cam.F90 +++ b/src/physics/cam/check_energy_cam.F90 @@ -171,7 +171,7 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & scaling_dycore(:ncol,:) = 1.0_r8 else ! Moist pressure... use phys formula - local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:,lchnk) + local_cp_or_cv_dycore(:,:) = local_cp_phys(:,:) scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling end if endif @@ -205,6 +205,11 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & latvap = latvap, & vc_physics = vc_physics, & vc_dycore = vc_dycore, & + name = name, & + flx_vap = flx_vap, & + flx_cnd = flx_cnd, & + flx_ice = flx_ice, & + flx_sen = flx_sen, & errmsg = errmsg, & errflg = errflg & ) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 227274355e..21e693fa17 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1390,7 +1390,7 @@ subroutine tphysac (ztodt, cam_in, & use aero_model, only: aero_model_drydep use carma_intr, only: carma_emission_tend, carma_timestep_tend use carma_flags_mod, only: carma_do_aerosol, carma_do_emission - use check_energy, only: check_energy_chng, tot_energy_phys + use check_energy, only: tot_energy_phys use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use time_manager, only: get_nstep use cam_abortutils, only: endrun @@ -1615,7 +1615,7 @@ subroutine tphysac (ztodt, cam_in, & call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& fh2o, surfric, obklen, flx_heat) end if - call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) + call check_energy_cam_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, & cam_in%cflx) end if @@ -1677,9 +1677,9 @@ subroutine tphysac (ztodt, cam_in, & call t_stopf('rayleigh_friction') if (do_clubb_sgs) then - call check_energy_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero) else - call check_energy_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, & + call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, & zero, cam_in%shf) endif @@ -1719,7 +1719,7 @@ subroutine tphysac (ztodt, cam_in, & call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, obklen=obklen, ustar=surfric) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero) call t_stopf('carma_timestep_tend') end if @@ -1759,7 +1759,7 @@ subroutine tphysac (ztodt, cam_in, & end if ! Check energy integrals - call check_energy_chng(state, tend, "gwdrag", nstep, ztodt, zero, & + call check_energy_cam_chng(state, tend, "gwdrag", nstep, ztodt, zero, & zero, zero, flx_heat) call t_stopf('gw_tend') @@ -1789,7 +1789,7 @@ subroutine tphysac (ztodt, cam_in, & end if ! Check energy integrals - call check_energy_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) ! Lunar tides call lunar_tides_tend( state, ptend ) @@ -1801,7 +1801,7 @@ subroutine tphysac (ztodt, cam_in, & end if call physics_update(state, ptend, ztodt, tend) ! Check energy integrals - call check_energy_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero) ! Ion drag calculation call t_startf ( 'iondrag' ) @@ -1851,7 +1851,7 @@ subroutine tphysac (ztodt, cam_in, & endif ! Check energy integrals - call check_energy_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) call t_stopf ( 'iondrag' ) @@ -1866,7 +1866,7 @@ subroutine tphysac (ztodt, cam_in, & call outfld( 'VTEND_NDG', ptend%v, pcols, lchnk) end if call physics_update(state,ptend,ztodt,tend) - call check_energy_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) endif !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv @@ -2257,7 +2257,7 @@ subroutine tphysbc (ztodt, state, & if (.not.dycore_is('EUL')) then call check_energy_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) + call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) end if call tot_energy_phys(state, 'phBP') @@ -2385,7 +2385,7 @@ subroutine tphysbc (ztodt, state, & ! Check energy integrals, including "reserved liquid" flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol) snow_dp(:ncol) = snow_dp(:ncol) + rice(:ncol) - call check_energy_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero) + call check_energy_cam_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero) snow_dp(:ncol) = snow_dp(:ncol) - rice(:ncol) ! @@ -2428,7 +2428,7 @@ subroutine tphysbc (ztodt, state, & end if flx_cnd(:ncol) = prec_sh(:ncol) + rliq2(:ncol) - call check_energy_chng(state, tend, "convect_shallow", nstep, ztodt, zero, flx_cnd, snow_sh, zero) + call check_energy_cam_chng(state, tend, "convect_shallow", nstep, ztodt, zero, flx_cnd, snow_sh, zero) call check_tracers_chng(state, tracerint, "convect_shallow", nstep, ztodt, zero_tracers) @@ -2460,9 +2460,9 @@ subroutine tphysbc (ztodt, state, & ! Before the detrainment, the reserved condensate is all liquid, but if CARMA is doing ! detrainment, then the reserved condensate is snow. if (carma_do_detrain) then - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero) else - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero) end if end if @@ -2484,7 +2484,7 @@ subroutine tphysbc (ztodt, state, & cam_in%ts, cam_in%sst, zdu) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "cldwat_tend", nstep, ztodt, zero, prec_str, snow_str, zero) + call check_energy_cam_chng(state, tend, "cldwat_tend", nstep, ztodt, zero, prec_str, snow_str, zero) call t_stopf('rk_stratiform_tend') @@ -2556,7 +2556,7 @@ subroutine tphysbc (ztodt, state, & flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "macrop_tend", nstep, ztodt, & + call check_energy_cam_chng(state, tend, "macrop_tend", nstep, ztodt, & zero, flx_cnd(:ncol)/cld_macmic_num_steps, & det_ice(:ncol)/cld_macmic_num_steps, & flx_heat(:ncol)/cld_macmic_num_steps) @@ -2603,7 +2603,7 @@ subroutine tphysbc (ztodt, state, & end if ! Use actual qflux (not lhf/latvap) for consistency with surface fluxes and revised code - call check_energy_chng(state, tend, "clubb_tend", nstep, ztodt, & + call check_energy_cam_chng(state, tend, "clubb_tend", nstep, ztodt, & cam_in%cflx(:ncol,1)/cld_macmic_num_steps, & flx_cnd(:ncol)/cld_macmic_num_steps, & det_ice(:ncol)/cld_macmic_num_steps, & @@ -2704,7 +2704,7 @@ subroutine tphysbc (ztodt, state, & flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) end if - call check_energy_chng(state_sc, tend_sc, "microp_tend_subcol", & + call check_energy_cam_chng(state_sc, tend_sc, "microp_tend_subcol", & nstep, ztodt, zero_sc, & prec_str_sc(:state_sc%ncol)/cld_macmic_num_steps, & snow_str_sc(:state_sc%ncol)/cld_macmic_num_steps, zero_sc) @@ -2736,7 +2736,7 @@ subroutine tphysbc (ztodt, state, & flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "microp_tend", nstep, ztodt, & + call check_energy_cam_chng(state, tend, "microp_tend", nstep, ztodt, & zero, prec_str(:ncol)/cld_macmic_num_steps, & snow_str(:ncol)/cld_macmic_num_steps, zero) @@ -2889,7 +2889,7 @@ subroutine tphysbc (ztodt, state, & flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) + call check_energy_cam_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) call t_stopf('radiation') diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index ad5bc9913e..b31a3fcd04 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -1644,7 +1644,7 @@ subroutine tphysac (ztodt, cam_in, & call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "clubb_emissions_tend", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "clubb_emissions_tend", nstep, ztodt, zero, zero, zero, zero) call t_stopf('clubb_emissions_tend') @@ -1672,9 +1672,9 @@ subroutine tphysac (ztodt, cam_in, & ! CARMA is doing ! detrainment, then the reserved condensate is snow. if (carma_do_detrain) then - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero) else - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero) end if end if @@ -1749,7 +1749,7 @@ subroutine tphysac (ztodt, cam_in, & end if ! Use actual qflux (not lhf/latvap) for consistency with surface fluxes and revised code - call check_energy_chng(state, tend, "clubb_tend", nstep, ztodt, & + call check_energy_cam_chng(state, tend, "clubb_tend", nstep, ztodt, & cam_in%cflx(:ncol,1)/cld_macmic_num_steps, & flx_cnd(:ncol)/cld_macmic_num_steps, & det_ice(:ncol)/cld_macmic_num_steps, & @@ -1848,7 +1848,7 @@ subroutine tphysac (ztodt, cam_in, & fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) end if - call check_energy_chng(state_sc, tend_sc, "microp_tend_subcol", & + call check_energy_cam_chng(state_sc, tend_sc, "microp_tend_subcol", & nstep, ztodt, zero_sc, & prec_str_sc(:state_sc%ncol)/cld_macmic_num_steps, & snow_str_sc(:state_sc%ncol)/cld_macmic_num_steps, zero_sc) @@ -1880,7 +1880,7 @@ subroutine tphysac (ztodt, cam_in, & fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "microp_tend", nstep, ztodt, & + call check_energy_cam_chng(state, tend, "microp_tend", nstep, ztodt, & zero, prec_str(:ncol)/cld_macmic_num_steps, & snow_str(:ncol)/cld_macmic_num_steps, zero) @@ -2036,7 +2036,7 @@ subroutine tphysac (ztodt, cam_in, & fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) + call check_energy_cam_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) call t_stopf('radiation') @@ -2113,7 +2113,7 @@ subroutine tphysac (ztodt, cam_in, & call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) end if - call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) + call check_energy_cam_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, & cam_in%cflx) end if @@ -2175,9 +2175,9 @@ subroutine tphysac (ztodt, cam_in, & call t_stopf('rayleigh_friction') if (do_clubb_sgs) then - call check_energy_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero) else - call check_energy_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, & + call check_energy_cam_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, & zero, cam_in%shf) endif @@ -2221,7 +2221,7 @@ subroutine tphysac (ztodt, cam_in, & call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, obklen=obklen, ustar=surfric) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero) call t_stopf('carma_timestep_tend') end if @@ -2260,7 +2260,7 @@ subroutine tphysac (ztodt, cam_in, & end if ! Check energy integrals - call check_energy_chng(state, tend, "gwdrag", nstep, ztodt, zero, & + call check_energy_cam_chng(state, tend, "gwdrag", nstep, ztodt, zero, & zero, zero, flx_heat) call t_stopf('gw_tend') @@ -2290,7 +2290,7 @@ subroutine tphysac (ztodt, cam_in, & end if ! Check energy integrals - call check_energy_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) ! Lunar tides call lunar_tides_tend( state, ptend ) @@ -2302,7 +2302,7 @@ subroutine tphysac (ztodt, cam_in, & end if call physics_update(state, ptend, ztodt, tend) ! Check energy integrals - call check_energy_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero) ! Ion drag calculation call t_startf ( 'iondrag' ) @@ -2351,7 +2351,7 @@ subroutine tphysac (ztodt, cam_in, & endif ! Check energy integrals - call check_energy_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) call t_stopf ( 'iondrag' ) @@ -2366,7 +2366,7 @@ subroutine tphysac (ztodt, cam_in, & call outfld( 'VTEND_NDG', ptend%v, pcols, lchnk) end if call physics_update(state,ptend,ztodt,tend) - call check_energy_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) endif !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv @@ -2517,7 +2517,8 @@ subroutine tphysbc (ztodt, state, & use convect_deep, only: convect_deep_tend use time_manager, only: is_first_step, get_nstep use convect_diagnostics,only: convect_diagnostics_calc - use check_energy, only: check_energy_chng, check_energy_fix + use check_energy_cam,only: check_energy_cam_chng + use check_energy, only: check_energy_fix use check_energy, only: check_tracers_data, check_tracers_init use check_energy, only: tot_energy_phys use dycore, only: dycore_is @@ -2701,7 +2702,7 @@ subroutine tphysbc (ztodt, state, & if (.not.dycore_is('EUL')) then call check_energy_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) + call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) end if @@ -2831,7 +2832,7 @@ subroutine tphysbc (ztodt, state, & ! Check energy integrals, including "reserved liquid" flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol) snow_dp(:ncol) = snow_dp(:ncol) + rice(:ncol) - call check_energy_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero) + call check_energy_cam_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero) snow_dp(:ncol) = snow_dp(:ncol) - rice(:ncol) !=================================================== diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 8c9c1586ef..22f1171351 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -503,7 +503,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) use air_composition, only: cpairv, cp_or_cv_dycore use time_manager, only: get_nstep use nudging, only: Nudge_Model, Nudge_ON, nudging_timestep_tend - use check_energy, only: check_energy_chng + use check_energy_cam,only: check_energy_cam_chng ! Arguments ! @@ -595,7 +595,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) if (Nudge_Model .and. Nudge_ON) then call nudging_timestep_tend(state,ptend) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) endif call tot_energy_phys(state, 'phAP') @@ -728,7 +728,8 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) use cam_diagnostics, only: diag_conv_tend_ini, diag_conv, diag_export use cam_history, only: outfld use time_manager, only: get_nstep - use check_energy, only: check_energy_chng, check_energy_fix, check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_chng + use check_energy, only: check_energy_fix, check_energy_timestep_init use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use chemistry, only: chem_is_active, chem_timestep_tend @@ -833,7 +834,7 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) if (adiabatic .and. (.not. dycore_is('EUL'))) then call check_energy_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) + call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) end if @@ -968,7 +969,7 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) ! surface flux is computed and supplied as an argument to ! check_energy_chng to account for how the simplified physics forcings are ! changing the total exnergy. - call check_energy_chng(state, tend, "tphysidl", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "tphysidl", nstep, ztodt, zero, zero, zero, zero) if (chem_is_active()) then call t_startf('simple_chem') diff --git a/src/physics/spcam/crm_physics.F90 b/src/physics/spcam/crm_physics.F90 index a1d9d2560d..735d6c3630 100644 --- a/src/physics/spcam/crm_physics.F90 +++ b/src/physics/spcam/crm_physics.F90 @@ -750,7 +750,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf, cam_in) use crmx_crm_module, only: crm use crmx_microphysics, only: nmicro_fields use physconst, only: latvap - use check_energy, only: check_energy_chng + use check_energy_cam, only: check_energy_cam_chng use phys_grid, only: get_rlat_all_p, get_rlon_all_p, get_lon_all_p, get_lat_all_p use modal_aero_calcsize, only: modal_aero_calcsize_sub use micro_pumas_utils, only: size_dist_param_liq, mg_liq_props, mincld, qsmall @@ -2152,7 +2152,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf, cam_in) end if ! check water and energy conservation - call check_energy_chng(state_loc, tend_loc, "crm_tend", nstep, ztodt, zero, & + call check_energy_cam_chng(state_loc, tend_loc, "crm_tend", nstep, ztodt, zero, & prec_dp(:ncol)+(qli_hydro(:ncol,2)-qli_hydro(:ncol,1))/ztodt/1000._r8, & snow_dp(:ncol)+(qi_hydro(:ncol,2)-qi_hydro(:ncol,1))/ztodt/1000._r8, radflux) diff --git a/src/physics/spcam/spcam_drivers.F90 b/src/physics/spcam/spcam_drivers.F90 index ebe6507607..00ab2448bd 100644 --- a/src/physics/spcam/spcam_drivers.F90 +++ b/src/physics/spcam/spcam_drivers.F90 @@ -300,7 +300,8 @@ subroutine tphysbc_spcam (ztodt, state, & use cam_history, only: outfld use constituents, only: pcnst, qmin, cnst_get_ind use time_manager, only: get_nstep - use check_energy, only: check_energy_chng, check_energy_fix + use check_energy_cam,only: check_energy_cam_chng + use check_energy, only: check_energy_fix use check_energy, only: check_tracers_data, check_tracers_init use dycore, only: dycore_is use radiation, only: radiation_tend @@ -443,7 +444,7 @@ subroutine tphysbc_spcam (ztodt, state, & if (dycore_is('LR') .or. dycore_is('SE')) then call check_energy_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) + call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld('EFIX', flx_heat, pcols,lchnk) end if ! Save state for convective tendency calculations. @@ -557,7 +558,7 @@ subroutine tphysbc_spcam (ztodt, state, & call physics_update(state, ptend, ztodt, tend) - call check_energy_chng(state, tend, "spradheat", nstep, ztodt, zero, zero, zero, zero) + call check_energy_cam_chng(state, tend, "spradheat", nstep, ztodt, zero, zero, zero, zero) call t_stopf('radiation') From 03209f0918ce253ea326963a9ee56eab1e2bd89c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 10 Sep 2024 00:16:43 -0400 Subject: [PATCH 04/20] Add missing USE in physbc --- src/physics/cam/physpkg.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 21e693fa17..bab29327b2 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1392,6 +1392,7 @@ subroutine tphysac (ztodt, cam_in, & use carma_flags_mod, only: carma_do_aerosol, carma_do_emission use check_energy, only: tot_energy_phys use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng + use check_energy_cam, only: check_energy_cam_chng use time_manager, only: get_nstep use cam_abortutils, only: endrun use dycore, only: dycore_is From 3cc3c3e33cc026cf3148a29a988b59a6d868cf81 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 10 Sep 2024 19:34:41 -0400 Subject: [PATCH 05/20] Fix array oob pcols --- src/physics/cam/check_energy_cam.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 index d5ae76ca48..7e1c289b96 100644 --- a/src/physics/cam/check_energy_cam.F90 +++ b/src/physics/cam/check_energy_cam.F90 @@ -49,7 +49,7 @@ subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) ! handling inside this shim module. (hplin, 9/9/24) if(state%psetcols == pcols) then ! No subcolumns - local_cp_phys(:,:) = cpairv(:,:,lchnk) + local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) else if (state%psetcols > pcols) then ! Subcolumns code @@ -148,10 +148,10 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & if(state%psetcols == pcols) then ! No subcolumns - local_cp_phys(:,:) = cpairv(:,:,lchnk) + local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling elseif(state%psetcols > pcols) then ! Subcolumns if(.not. all(cpairv(:,:,:) == cpair)) then @@ -164,14 +164,14 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & if (vc_dycore == vc_height) then ! compute cv if vertical coordinate is height: cv = cp - R local_cp_or_cv_dycore(:ncol,:) = cpair-rair - scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling else if (vc_dycore == vc_dry_pressure) then ! SE specific hydrostatic energy local_cp_or_cv_dycore(:ncol,:) = cpair scaling_dycore(:ncol,:) = 1.0_r8 else ! Moist pressure... use phys formula - local_cp_or_cv_dycore(:,:) = local_cp_phys(:,:) + local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:) scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling end if endif From f7f4139a6f1333eb80375d6adcfb8a34247b22e3 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 7 Oct 2024 10:56:11 -0600 Subject: [PATCH 06/20] Update atmos_phys to check_energy version --- .gitmodules | 4 +- src/atmos_phys | 2 +- src/physics/cam/cam_diagnostics.F90 | 2 +- src/physics/cam/check_energy.F90 | 164 ++------------------------- src/physics/cam/check_energy_cam.F90 | 164 ++++++++++++++++++++++++++- src/physics/cam/physpkg.F90 | 6 +- src/physics/cam7/physpkg.F90 | 7 +- src/physics/simple/physpkg.F90 | 8 +- src/physics/spcam/spcam_drivers.F90 | 5 +- 9 files changed, 184 insertions(+), 178 deletions(-) diff --git a/.gitmodules b/.gitmodules index 224ee10052..e295d8fec1 100644 --- a/.gitmodules +++ b/.gitmodules @@ -35,8 +35,8 @@ [submodule "atmos_phys"] path = src/atmos_phys - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_05_001 + url = https://github.com/jimmielin/atmospheric_physics + fxtag = df80b9c1 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index f8ce60bf40..df80b9c1e7 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit f8ce60bf40f800623f8eb3065021ec5dfa9e6b45 +Subproject commit df80b9c1e73064f9a0196282a92b66a2041f6a57 diff --git a/src/physics/cam/cam_diagnostics.F90 b/src/physics/cam/cam_diagnostics.F90 index 97dad2ba01..00d6859a84 100644 --- a/src/physics/cam/cam_diagnostics.F90 +++ b/src/physics/cam/cam_diagnostics.F90 @@ -2019,7 +2019,7 @@ subroutine diag_phys_tend_writeout_dry(state, pbuf, tend, ztodt) ! !--------------------------------------------------------------- - use check_energy, only: check_energy_get_integrals + use check_energy_cam,only: check_energy_get_integrals use physconst, only: cpair ! Arguments diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 5e1df3a967..d71f3df8f6 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -24,13 +24,11 @@ module check_energy use ppgrid, only: pcols, pver, begchunk, endchunk use spmd_utils, only: masterproc - use gmean_mod, only: gmean use physconst, only: gravit, rga, latvap, latice, cpair, rair use air_composition, only: cpairv, cp_or_cv_dycore - use physics_types, only: physics_state, physics_tend, physics_ptend, physics_ptend_init + use physics_types, only: physics_state use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind use cam_logfile, only: iulog - use scamMod, only: single_column, use_camiop, heat_glob_scm use cam_history, only: outfld, write_camiop implicit none @@ -42,10 +40,7 @@ module check_energy ! Public methods public :: check_energy_readnl ! read namelist values public :: check_energy_register ! register fields in physics buffer - public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean public :: check_energy_init ! initialization of module - public :: check_energy_gmean ! global means of physics input and output total energy - public :: check_energy_fix ! add global mean energy difference as a heating public :: check_tracers_init ! initialize tracer integrals and cumulative boundary fluxes public :: check_tracers_chng ! check changes in integrals against cumulative boundary fluxes @@ -55,13 +50,6 @@ module check_energy logical :: print_energy_errors = .false. - real(r8) :: teout_glob ! global mean energy of output state - real(r8) :: teinp_glob ! global mean energy of input state - real(r8) :: tedif_glob ! global mean energy difference - real(r8) :: psurf_glob ! global mean surface pressure - real(r8) :: ptopb_glob ! global mean top boundary pressure - real(r8) :: heat_glob ! global mean heating rate - ! Physics buffer indices integer, public :: teout_idx = 0 ! teout index in physics buffer @@ -88,6 +76,9 @@ subroutine check_energy_readnl(nlfile) use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_logical use cam_abortutils, only: endrun + ! update the CCPP-ized namelist option + use check_energy_chng, only: check_energy_chng_init + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input ! Local variables @@ -120,6 +111,9 @@ subroutine check_energy_readnl(nlfile) write(iulog,*) ' print_energy_errors =', print_energy_errors end if + ! update the CCPP-ized namelist option + call check_energy_chng_init(print_energy_errors_in=print_energy_errors) + end subroutine check_energy_readnl !=============================================================================== @@ -154,28 +148,6 @@ subroutine check_energy_register() end subroutine check_energy_register -!=============================================================================== - -subroutine check_energy_get_integrals( tedif_glob_out, heat_glob_out ) - -!----------------------------------------------------------------------- -! Purpose: Return energy integrals -!----------------------------------------------------------------------- - - real(r8), intent(out), optional :: tedif_glob_out - real(r8), intent(out), optional :: heat_glob_out - -!----------------------------------------------------------------------- - - if ( present(tedif_glob_out) ) then - tedif_glob_out = tedif_glob - endif - if ( present(heat_glob_out) ) then - heat_glob_out = heat_glob - endif - -end subroutine check_energy_get_integrals - !================================================================================================ subroutine check_energy_init() @@ -214,128 +186,6 @@ subroutine check_energy_init() end subroutine check_energy_init -!=============================================================================== - - subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) - - use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk - use physics_types, only: dyn_te_idx - use cam_history, only: write_camiop -!----------------------------------------------------------------------- -! Compute global mean total energy of physics input and output states -! computed consistently with dynamical core vertical coordinate -! (under hydrostatic assumption) -!----------------------------------------------------------------------- -!------------------------------Arguments-------------------------------- - - type(physics_state), intent(in ), dimension(begchunk:endchunk) :: state - type(physics_buffer_desc), pointer :: pbuf2d(:,:) - - real(r8), intent(in) :: dtime ! physics time step - integer , intent(in) :: nstep ! current timestep number - -!---------------------------Local storage------------------------------- - integer :: ncol ! number of active columns - integer :: lchnk ! chunk index - - real(r8) :: te(pcols,begchunk:endchunk,4) - ! total energy of input/output states (copy) - real(r8) :: te_glob(4) ! global means of total energy - real(r8), pointer :: teout(:) -!----------------------------------------------------------------------- - - ! Copy total energy out of input and output states - do lchnk = begchunk, endchunk - ncol = state(lchnk)%ncol - ! input energy using dynamical core energy formula - te(:ncol,lchnk,1) = state(lchnk)%te_ini(:ncol,dyn_te_idx) - ! output energy - call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout) - - te(:ncol,lchnk,2) = teout(1:ncol) - ! surface pressure for heating rate - te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1) - ! model top pressure for heating rate (not constant for z-based vertical coordinate!) - te(:ncol,lchnk,4) = state(lchnk)%pint(:ncol,1) - end do - - ! Compute global means of input and output energies and of - ! surface pressure for heating rate (assume uniform ptop) - call gmean(te, te_glob, 4) - - if (begchunk .le. endchunk) then - teinp_glob = te_glob(1) - teout_glob = te_glob(2) - psurf_glob = te_glob(3) - ptopb_glob = te_glob(4) - - ! Global mean total energy difference - tedif_glob = teinp_glob - teout_glob - heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob) - if (masterproc) then - write(iulog,'(1x,a9,1x,i8,5(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, & - heat_glob, psurf_glob, ptopb_glob - end if - else - heat_glob = 0._r8 - end if ! (begchunk .le. endchunk) - - end subroutine check_energy_gmean - -!=============================================================================== - subroutine check_energy_fix(state, ptend, nstep, eshflx) - -!----------------------------------------------------------------------- -! Add heating rate required for global mean total energy conservation -!----------------------------------------------------------------------- -!------------------------------Arguments-------------------------------- - - type(physics_state), intent(in ) :: state - type(physics_ptend), intent(out) :: ptend - - integer , intent(in ) :: nstep ! time step number - real(r8), intent(out ) :: eshflx(pcols) ! effective sensible heat flux - -!---------------------------Local storage------------------------------- - integer :: i ! column - integer :: ncol ! number of atmospheric columns in chunk - integer :: lchnk ! chunk number - real(r8) :: heat_out(pcols) -!----------------------------------------------------------------------- - lchnk = state%lchnk - ncol = state%ncol - - call physics_ptend_init(ptend, state%psetcols, 'chkenergyfix', ls=.true.) - -#if ( defined OFFLINE_DYN ) - ! disable the energy fix for offline driver - heat_glob = 0._r8 -#endif - - ! Special handling of energy fix for SCAM - supplied via CAMIOP - zero's for normal IOPs - if (single_column) then - if ( use_camiop) then - heat_glob = heat_glob_scm(1) - else - heat_glob = 0._r8 - endif - endif - ptend%s(:ncol,:pver) = heat_glob - - if (nstep > 0 .and. write_camiop) then - heat_out(:ncol) = heat_glob - call outfld('heat_glob', heat_out(:ncol), pcols, lchnk) - endif - -! compute effective sensible heat flux - do i = 1, ncol - eshflx(i) = heat_glob * (state%pint(i,pver+1) - state%pint(i,1)) * rga - end do - - return - end subroutine check_energy_fix - - !=============================================================================== subroutine check_tracers_init(state, tracerint) diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 index 7e1c289b96..6931cfb3ee 100644 --- a/src/physics/cam/check_energy_cam.F90 +++ b/src/physics/cam/check_energy_cam.F90 @@ -1,4 +1,4 @@ -! Shim for CCPP-ized check_energy routines in CAM +! CAM interface for CCPP-ized check_energy routines in CAM module check_energy_cam use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver @@ -8,9 +8,23 @@ module check_energy_cam implicit none private - public :: check_energy_cam_chng ! check changes in integrals against cumulative boundary fluxes + public :: check_energy_cam_chng ! check changes in integrals against cumulative boundary fluxes public :: check_energy_cam_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes + public :: check_energy_cam_fix ! add heating rate required for global mean total energy conservation + + ! These subroutines cannot be CCPP-ized and replicated from check_energy + public :: check_energy_gmean ! global means of physics input and output total energy + public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean + + ! Module variables used for check_energy_gmean + real(r8) :: teout_glob ! global mean energy of output state + real(r8) :: teinp_glob ! global mean energy of input state + real(r8) :: tedif_glob ! global mean energy difference + real(r8) :: psurf_glob ! global mean surface pressure + real(r8) :: ptopb_glob ! global mean top boundary pressure + real(r8) :: heat_glob ! global mean heating rate + contains ! Compute initial values of energy and water integrals, ! zero cumulative tendencies @@ -215,4 +229,148 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & ) end subroutine check_energy_cam_chng -end module check_energy_cam \ No newline at end of file + + ! Add heating rate required for global mean total energy conservation + subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) + use physics_types, only: physics_ptend, physics_ptend_init + use physconst, only: rga + + ! SCAM support + use scamMod, only: single_column, use_camiop, heat_glob_scm + + ! CCPP-ized subroutine + use check_energy_fix, only: check_energy_fix_run + + type(physics_state), intent(in) :: state + type(physics_ptend), intent(out) :: ptend + + integer , intent(in) :: nstep ! time step number + real(r8), intent(out) :: eshflx(pcols) ! effective sensible heat flux + + integer :: ncol ! number of atmospheric columns in chunk + integer :: lchnk ! chunk number + real(r8) :: heat_out(pcols) + + lchnk = state%lchnk + ncol = state%ncol + + call physics_ptend_init(ptend, state%psetcols, 'chkenergyfix', ls=.true.) + +#if ( defined OFFLINE_DYN ) + ! disable the energy fix for offline driver + heat_glob = 0._r8 +#endif + + ! Special handling of energy fix for SCAM - supplied via CAMIOP - zero's for normal IOPs + if (single_column) then + if (use_camiop) then + heat_glob = heat_glob_scm(1) + else + heat_glob = 0._r8 + endif + endif + + if (nstep > 0 .and. write_camiop) then + heat_out(:ncol) = heat_glob + call outfld('heat_glob', heat_out(:ncol), pcols, lchnk) + endif + + ! Call the CCPP-ized subroutine (for non-SCAM) + ! to compute the effective sensible heat flux and save to ptend%s + call check_energy_fix_run( & + ncol = ncol, & + pver = pver, & + pint = state%pint(:ncol,:), & + rga = rga, & + heat_glob = heat_glob. & + ptend_s = ptend%s(:ncol,:), & + eshflx = eshflx(:ncol) & + ) + + end subroutine check_energy_cam_fix + + ! Compute global mean total energy of physics input and output states + ! computed consistently with dynamical core vertical coordinate + ! (under hydrostatic assumption) + ! + ! This subroutine cannot use the CCPP-ized equivalent because + ! it is dependent on chunks. + subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) + use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_chunk + use physics_types, only: dyn_te_idx + use cam_history, only: write_camiop + use ppgrid, only: begchunk, endchunk + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use gmean_mod, only: gmean + use physconst, only: gravit + + ! To remove, pbuf indices in check_energy + use check_energy, only: teout_idx + + type(physics_state), intent(in), dimension(begchunk:endchunk) :: state + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + real(r8), intent(in) :: dtime ! physics time step + integer , intent(in) :: nstep ! current timestep number + + integer :: ncol ! number of active columns + integer :: lchnk ! chunk index + + real(r8) :: te(pcols,begchunk:endchunk,4) + ! total energy of input/output states (copy) + real(r8) :: te_glob(4) ! global means of total energy + real(r8), pointer :: teout(:) + + ! Copy total energy out of input and output states + do lchnk = begchunk, endchunk + ncol = state(lchnk)%ncol + ! input energy using dynamical core energy formula + te(:ncol,lchnk,1) = state(lchnk)%te_ini(:ncol,dyn_te_idx) + ! output energy + call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout) + + te(:ncol,lchnk,2) = teout(1:ncol) + ! surface pressure for heating rate + te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1) + ! model top pressure for heating rate (not constant for z-based vertical coordinate!) + te(:ncol,lchnk,4) = state(lchnk)%pint(:ncol,1) + end do + + ! Compute global means of input and output energies and of + ! surface pressure for heating rate (assume uniform ptop) + call gmean(te, te_glob, 4) + + if (begchunk .le. endchunk) then + teinp_glob = te_glob(1) + teout_glob = te_glob(2) + psurf_glob = te_glob(3) + ptopb_glob = te_glob(4) + + ! Global mean total energy difference + tedif_glob = teinp_glob - teout_glob + heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob) + if (masterproc) then + write(iulog,'(1x,a9,1x,i8,5(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, & + heat_glob, psurf_glob, ptopb_glob + end if + else + heat_glob = 0._r8 + end if ! (begchunk .le. endchunk) + + end subroutine check_energy_gmean + + ! Return energy integrals (module variables) + subroutine check_energy_get_integrals(tedif_glob_out, heat_glob_out) + real(r8), intent(out), optional :: tedif_glob_out + real(r8), intent(out), optional :: heat_glob_out + + if ( present(tedif_glob_out) ) then + tedif_glob_out = tedif_glob + endif + + if ( present(heat_glob_out) ) then + heat_glob_out = heat_glob + endif + end subroutine check_energy_get_integrals +end module check_energy_cam diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index bab29327b2..25f46c3a45 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1065,7 +1065,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics,only: diag_allocate, diag_physvar_ic - use check_energy, only: check_energy_gmean + use check_energy_cam, only: check_energy_gmean use phys_control, only: phys_getopts use spcam_drivers, only: tphysbc_spcam use spmd_utils, only: mpicom @@ -2046,7 +2046,7 @@ subroutine tphysbc (ztodt, state, & use time_manager, only: is_first_step, get_nstep use convect_shallow, only: convect_shallow_tend use check_energy_cam,only: check_energy_cam_timestep_init, check_energy_cam_chng - use check_energy, only: check_energy_fix + use check_energy_cam,only: check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use dycore, only: dycore_is @@ -2256,7 +2256,7 @@ subroutine tphysbc (ztodt, state, & call tot_energy_phys(state, 'phBF') call tot_energy_phys(state, 'dyBF',vc=vc_dycore) if (.not.dycore_is('EUL')) then - call check_energy_fix(state, ptend, nstep, flx_heat) + call check_energy_cam_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index b31a3fcd04..1ce28d3a93 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -1056,7 +1056,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics,only: diag_allocate, diag_physvar_ic - use check_energy, only: check_energy_gmean + use check_energy_cam, only: check_energy_gmean use spmd_utils, only: mpicom use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate use cam_history, only: outfld, write_camiop @@ -2517,8 +2517,7 @@ subroutine tphysbc (ztodt, state, & use convect_deep, only: convect_deep_tend use time_manager, only: is_first_step, get_nstep use convect_diagnostics,only: convect_diagnostics_calc - use check_energy_cam,only: check_energy_cam_chng - use check_energy, only: check_energy_fix + use check_energy_cam,only: check_energy_cam_chng, check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init use check_energy, only: tot_energy_phys use dycore, only: dycore_is @@ -2700,7 +2699,7 @@ subroutine tphysbc (ztodt, state, & call tot_energy_phys(state, 'dyBF',vc=vc_dycore) if (.not.dycore_is('EUL')) then - call check_energy_fix(state, ptend, nstep, flx_heat) + call check_energy_cam_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 22f1171351..3e594d9073 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -305,7 +305,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics, only: diag_allocate, diag_physvar_ic - use check_energy, only: check_energy_gmean + use check_energy_cam, only: check_energy_gmean use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate ! @@ -728,8 +728,8 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) use cam_diagnostics, only: diag_conv_tend_ini, diag_conv, diag_export use cam_history, only: outfld use time_manager, only: get_nstep - use check_energy_cam, only: check_energy_cam_chng - use check_energy, only: check_energy_fix, check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_chng, check_energy_cam_fix + use check_energy, only: check_energy_timestep_init use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use chemistry, only: chem_is_active, chem_timestep_tend @@ -832,7 +832,7 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) call t_startf('energy_fixer') if (adiabatic .and. (.not. dycore_is('EUL'))) then - call check_energy_fix(state, ptend, nstep, flx_heat) + call check_energy_cam_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) diff --git a/src/physics/spcam/spcam_drivers.F90 b/src/physics/spcam/spcam_drivers.F90 index 00ab2448bd..cd7d3f6370 100644 --- a/src/physics/spcam/spcam_drivers.F90 +++ b/src/physics/spcam/spcam_drivers.F90 @@ -300,8 +300,7 @@ subroutine tphysbc_spcam (ztodt, state, & use cam_history, only: outfld use constituents, only: pcnst, qmin, cnst_get_ind use time_manager, only: get_nstep - use check_energy_cam,only: check_energy_cam_chng - use check_energy, only: check_energy_fix + use check_energy_cam,only: check_energy_cam_chng, check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init use dycore, only: dycore_is use radiation, only: radiation_tend @@ -442,7 +441,7 @@ subroutine tphysbc_spcam (ztodt, state, & call t_startf('energy_fixer') if (dycore_is('LR') .or. dycore_is('SE')) then - call check_energy_fix(state, ptend, nstep, flx_heat) + call check_energy_cam_fix(state, ptend, nstep, flx_heat) call physics_update(state, ptend, ztodt, tend) call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld('EFIX', flx_heat, pcols,lchnk) From f6f464bb6b73aa52c2850cdb4595ab4148774fc4 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 10 Oct 2024 16:26:34 -0600 Subject: [PATCH 07/20] Fix build errors -- does not include bld/configure atmos_phys path updates --- src/physics/cam/check_energy.F90 | 2 +- src/physics/cam/check_energy_cam.F90 | 12 ++++++++---- src/physics/simple/physpkg.F90 | 2 +- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index d71f3df8f6..9f9faa10fe 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -29,7 +29,7 @@ module check_energy use physics_types, only: physics_state use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind use cam_logfile, only: iulog - use cam_history, only: outfld, write_camiop + use cam_history, only: outfld implicit none private diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 index 6931cfb3ee..e1f143c0c2 100644 --- a/src/physics/cam/check_energy_cam.F90 +++ b/src/physics/cam/check_energy_cam.F90 @@ -50,6 +50,7 @@ subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) real(r8) :: local_cp_phys(state%psetcols,pver) real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) + real(r8) :: teout(state%ncol) ! dummy teout argument integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns character(len=512) :: errmsg @@ -91,6 +92,7 @@ subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) ncol = ncol, & pver = pver, & pcnst = pcnst, & + is_first_timestep = is_first_step(), & q = state%q(1:ncol,1:pver,1:pcnst), & pdel = state%pdel(1:ncol,1:pver), & u = state%u(1:ncol,1:pver), & @@ -99,8 +101,6 @@ subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) pintdry = state%pintdry(1:ncol,1:pver), & phis = state%phis(1:ncol), & zm = state%zm(1:ncol,:), & - temp_ini = state%temp_ini(:ncol,:), & - z_ini = state%z_ini(:ncol,:), & cp_phys = local_cp_phys(1:ncol,:), & cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & te_ini_phys = state%te_ini(1:ncol,phys_te_idx), & @@ -111,7 +111,10 @@ subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) tw_cur = state%tw_cur(1:ncol), & tend_te_tnd = tend%te_tnd(1:ncol), & tend_tw_tnd = tend%tw_tnd(1:ncol), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & count = state%count, & + teout = teout(1:ncol), & ! dummy argument - actual teout written to pbuf directly below vc_physics = vc_physics, & vc_dycore = vc_dycore, & errmsg = errmsg, & @@ -237,6 +240,8 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) ! SCAM support use scamMod, only: single_column, use_camiop, heat_glob_scm + use cam_history, only: write_camiop + use cam_history, only: outfld ! CCPP-ized subroutine use check_energy_fix, only: check_energy_fix_run @@ -282,7 +287,7 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) pver = pver, & pint = state%pint(:ncol,:), & rga = rga, & - heat_glob = heat_glob. & + heat_glob = heat_glob, & ptend_s = ptend%s(:ncol,:), & eshflx = eshflx(:ncol) & ) @@ -298,7 +303,6 @@ end subroutine check_energy_cam_fix subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_chunk use physics_types, only: dyn_te_idx - use cam_history, only: write_camiop use ppgrid, only: begchunk, endchunk use spmd_utils, only: masterproc use cam_logfile, only: iulog diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 3e594d9073..eaea0ff1dd 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -729,7 +729,7 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) use cam_history, only: outfld use time_manager, only: get_nstep use check_energy_cam, only: check_energy_cam_chng, check_energy_cam_fix - use check_energy, only: check_energy_timestep_init + use check_energy_cam, only: check_energy_cam_timestep_init use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use chemistry, only: chem_is_active, chem_timestep_tend From 2ffa522d8403f73aba0f306a2c35dd50319c33a5 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 14 Oct 2024 12:16:26 -0400 Subject: [PATCH 08/20] Move back to check_energy.F90 for FV3 external compatibility --- src/dynamics/eul/dp_coupling.F90 | 4 +- src/dynamics/fv/dp_coupling.F90 | 4 +- src/dynamics/mpas/dp_coupling.F90 | 4 +- src/dynamics/se/dp_coupling.F90 | 4 +- src/physics/cam/cam_diagnostics.F90 | 2 +- src/physics/cam/check_energy.F90 | 390 ++++++++++++++++++++++++++- src/physics/cam/check_energy_cam.F90 | 380 -------------------------- src/physics/cam/physpkg.F90 | 12 +- src/physics/cam7/physpkg.F90 | 8 +- src/physics/simple/physpkg.F90 | 8 +- src/physics/spcam/crm_physics.F90 | 2 +- src/physics/spcam/spcam_drivers.F90 | 2 +- 12 files changed, 403 insertions(+), 417 deletions(-) delete mode 100644 src/physics/cam/check_energy_cam.F90 diff --git a/src/dynamics/eul/dp_coupling.F90 b/src/dynamics/eul/dp_coupling.F90 index 8e3eb64a21..bc900e2d0e 100644 --- a/src/dynamics/eul/dp_coupling.F90 +++ b/src/dynamics/eul/dp_coupling.F90 @@ -14,7 +14,7 @@ module dp_coupling use physconst, only: cpair, gravit, rair, zvir use air_composition, only: rairv use geopotential, only: geopotential_t - use check_energy_cam, only: check_energy_cam_timestep_init + use check_energy, only: check_energy_timestep_init #if (defined SPMD) use spmd_dyn, only: buf1, buf1win, buf2, buf2win, & spmdbuf_siz, local_dp_map, & @@ -291,7 +291,7 @@ subroutine d_p_coupling(ps, t3, u3, v3, q3, & ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk ) + call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk ) end do diff --git a/src/dynamics/fv/dp_coupling.F90 b/src/dynamics/fv/dp_coupling.F90 index af2d13b12e..fc02821471 100644 --- a/src/dynamics/fv/dp_coupling.F90 +++ b/src/dynamics/fv/dp_coupling.F90 @@ -12,7 +12,7 @@ module dp_coupling use physconst, only: gravit, zvir use air_composition, only: cpairv, rairv use geopotential, only: geopotential_t - use check_energy_cam, only: check_energy_cam_timestep_init + use check_energy, only: check_energy_timestep_init use dynamics_vars, only: T_FVDYCORE_GRID, t_fvdycore_state use dyn_internal_state,only: get_dyn_state use dyn_comp, only: dyn_import_t, dyn_export_t, fv_print_dpcoup_warn @@ -621,7 +621,7 @@ subroutine d_p_coupling(grid, phys_state, phys_tend, pbuf2d, dyn_out) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do call t_stopf('derived_fields') diff --git a/src/dynamics/mpas/dp_coupling.F90 b/src/dynamics/mpas/dp_coupling.F90 index c312f1e2f9..10d75b4b8c 100644 --- a/src/dynamics/mpas/dp_coupling.F90 +++ b/src/dynamics/mpas/dp_coupling.F90 @@ -413,7 +413,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d) ! MPAS prognostic fields. use geopotential, only: geopotential_t - use check_energy_cam, only: check_energy_cam_timestep_init + use check_energy, only: check_energy_timestep_init use shr_vmath_mod, only: shr_vmath_log use phys_control, only: waccmx_is use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update @@ -564,7 +564,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 9de31132a2..919b7f3510 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -546,7 +546,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) use phys_control, only: waccmx_is use geopotential, only: geopotential_t use static_energy, only: update_dry_static_energy_run - use check_energy_cam,only: check_energy_cam_timestep_init + 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 @@ -707,7 +707,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) cpairv(1:ncol,:,lchnk), errflg, errmsg) ! Compute energy and water integrals of input state pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) - call check_energy_cam_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) + call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk) end do ! lchnk diff --git a/src/physics/cam/cam_diagnostics.F90 b/src/physics/cam/cam_diagnostics.F90 index 00d6859a84..97dad2ba01 100644 --- a/src/physics/cam/cam_diagnostics.F90 +++ b/src/physics/cam/cam_diagnostics.F90 @@ -2019,7 +2019,7 @@ subroutine diag_phys_tend_writeout_dry(state, pbuf, tend, ztodt) ! !--------------------------------------------------------------- - use check_energy_cam,only: check_energy_get_integrals + use check_energy, only: check_energy_get_integrals use physconst, only: cpair ! Arguments diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 9f9faa10fe..f318cb51e1 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -21,37 +21,52 @@ module check_energy !--------------------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: pcols, pver, begchunk, endchunk + use ppgrid, only: pcols, pver use spmd_utils, only: masterproc - use physconst, only: gravit, rga, latvap, latice, cpair, rair + use physconst, only: rga use air_composition, only: cpairv, cp_or_cv_dycore use physics_types, only: physics_state use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind use cam_logfile, only: iulog - use cam_history, only: outfld implicit none private -! Public types: + ! Public types: public check_tracers_data -! Public methods - public :: check_energy_readnl ! read namelist values - public :: check_energy_register ! register fields in physics buffer - public :: check_energy_init ! initialization of module - public :: check_tracers_init ! initialize tracer integrals and cumulative boundary fluxes - public :: check_tracers_chng ! check changes in integrals against cumulative boundary fluxes + ! Public methods - not CCPP-ized + public :: check_tracers_init ! initialize tracer integrals and cumulative boundary fluxes + public :: check_tracers_chng ! check changes in integrals against cumulative boundary fluxes + public :: tot_energy_phys ! calculate and output total energy and axial angular momentum diagnostics - public :: tot_energy_phys ! calculate and output total energy and axial angular momentum diagnostics + ! These subroutines cannot be CCPP-ized + public :: check_energy_readnl ! read namelist values + public :: check_energy_register ! register fields in physics buffer + public :: check_energy_init ! initialization of module + public :: check_energy_gmean ! global means of physics input and output total energy + public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean -! Private module data + ! Public methods - CAM interfaces to CCPP version: + public :: check_energy_cam_chng ! check changes in integrals against cumulative boundary fluxes + public :: check_energy_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes + ! name is retained for FV3 compatibility + public :: check_energy_cam_fix ! add heating rate required for global mean total energy conservation + + ! Private module data logical :: print_energy_errors = .false. -! Physics buffer indices + ! used for check_energy_gmean + real(r8) :: teout_glob ! global mean energy of output state + real(r8) :: teinp_glob ! global mean energy of input state + real(r8) :: tedif_glob ! global mean energy difference + real(r8) :: psurf_glob ! global mean surface pressure + real(r8) :: ptopb_glob ! global mean top boundary pressure + real(r8) :: heat_glob ! global mean heating rate + ! Physics buffer indices integer, public :: teout_idx = 0 ! teout index in physics buffer integer, public :: dtcore_idx = 0 ! dtcore index in physics buffer integer, public :: dqcore_idx = 0 ! dqcore index in physics buffer @@ -542,5 +557,354 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) end subroutine tot_energy_phys + ! Compute global mean total energy of physics input and output states + ! computed consistently with dynamical core vertical coordinate + ! (under hydrostatic assumption) + ! + ! This subroutine cannot use the CCPP-ized equivalent because + ! it is dependent on chunks. + subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) + use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_chunk + use physics_types, only: dyn_te_idx + use ppgrid, only: begchunk, endchunk + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use gmean_mod, only: gmean + use physconst, only: gravit + + type(physics_state), intent(in), dimension(begchunk:endchunk) :: state + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + real(r8), intent(in) :: dtime ! physics time step + integer , intent(in) :: nstep ! current timestep number + + integer :: ncol ! number of active columns + integer :: lchnk ! chunk index + + real(r8) :: te(pcols,begchunk:endchunk,4) + ! total energy of input/output states (copy) + real(r8) :: te_glob(4) ! global means of total energy + real(r8), pointer :: teout(:) + + ! Copy total energy out of input and output states + do lchnk = begchunk, endchunk + ncol = state(lchnk)%ncol + ! input energy using dynamical core energy formula + te(:ncol,lchnk,1) = state(lchnk)%te_ini(:ncol,dyn_te_idx) + ! output energy + call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout) + + te(:ncol,lchnk,2) = teout(1:ncol) + ! surface pressure for heating rate + te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1) + ! model top pressure for heating rate (not constant for z-based vertical coordinate!) + te(:ncol,lchnk,4) = state(lchnk)%pint(:ncol,1) + end do + + ! Compute global means of input and output energies and of + ! surface pressure for heating rate (assume uniform ptop) + call gmean(te, te_glob, 4) + + if (begchunk .le. endchunk) then + teinp_glob = te_glob(1) + teout_glob = te_glob(2) + psurf_glob = te_glob(3) + ptopb_glob = te_glob(4) + + ! Global mean total energy difference + tedif_glob = teinp_glob - teout_glob + heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob) + if (masterproc) then + write(iulog,'(1x,a9,1x,i8,5(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, & + heat_glob, psurf_glob, ptopb_glob + end if + else + heat_glob = 0._r8 + end if ! (begchunk .le. endchunk) + + end subroutine check_energy_gmean + + ! Return energy integrals (module variables) + subroutine check_energy_get_integrals(tedif_glob_out, heat_glob_out) + real(r8), intent(out), optional :: tedif_glob_out + real(r8), intent(out), optional :: heat_glob_out + + if ( present(tedif_glob_out) ) then + tedif_glob_out = tedif_glob + endif + + if ( present(heat_glob_out) ) then + heat_glob_out = heat_glob + endif + end subroutine check_energy_get_integrals + + ! Compute initial values of energy and water integrals, + ! zero cumulative tendencies + subroutine check_energy_timestep_init(state, tend, pbuf, col_type) + use physics_buffer, only: physics_buffer_desc, pbuf_set_field + use cam_abortutils, only: endrun + use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure + use physics_types, only: physics_tend + use physics_types, only: phys_te_idx, dyn_te_idx + use time_manager, only: is_first_step + use physconst, only: cpair, rair + use air_composition, only: cpairv, cp_or_cv_dycore + + ! CCPP-ized subroutine + use check_energy_chng, only: check_energy_chng_timestep_init + + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + type(physics_buffer_desc), pointer :: pbuf(:) + integer, optional :: col_type ! Flag indicating whether using grid or subcolumns + + real(r8) :: local_cp_phys(state%psetcols,pver) + real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) + real(r8) :: teout(state%ncol) ! dummy teout argument + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + character(len=512) :: errmsg + integer :: errflg + + lchnk = state%lchnk + ncol = state%ncol + + ! The code below is split into not-subcolumns and subcolumns code, as there is different handling of the + ! cp passed into the hydrostatic energy call. CAM-SIMA does not support subcolumns, so we keep this special + ! handling inside this shim module. (hplin, 9/9/24) + if(state%psetcols == pcols) then + ! No subcolumns + local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) + local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) + else if (state%psetcols > pcols) then + ! Subcolumns code + ! Subcolumns specific error handling + if(.not. all(cpairv(:,:,lchnk) == cpair)) then + call endrun('check_energy_timestep_init: cpairv is not allowed to vary when subcolumns are turned on') + endif + + local_cp_phys(1:ncol,:) = cpair + + if (vc_dycore == vc_height) then + ! MPAS specific hydrostatic energy computation (internal energy) + local_cp_or_cv_dycore(:ncol,:) = cpair-rair + else if(vc_dycore == vc_dry_pressure) then + ! SE specific hydrostatic energy (enthalpy) + local_cp_or_cv_dycore(:ncol,:) = cpair + else + ! cp_or_cv is not used in the underlying subroutine, zero it out to be sure + local_cp_or_cv_dycore(:ncol,:) = 0.0_r8 + endif + end if + + ! Call CCPP-ized underlying subroutine. + call check_energy_chng_timestep_init( & + ncol = ncol, & + pver = pver, & + pcnst = pcnst, & + is_first_timestep = is_first_step(), & + q = state%q(1:ncol,1:pver,1:pcnst), & + pdel = state%pdel(1:ncol,1:pver), & + u = state%u(1:ncol,1:pver), & + v = state%v(1:ncol,1:pver), & + T = state%T(1:ncol,1:pver), & + pintdry = state%pintdry(1:ncol,1:pver), & + phis = state%phis(1:ncol), & + zm = state%zm(1:ncol,:), & + cp_phys = local_cp_phys(1:ncol,:), & + cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & + te_ini_phys = state%te_ini(1:ncol,phys_te_idx), & + te_ini_dyn = state%te_ini(1:ncol,dyn_te_idx), & + tw_ini = state%tw_ini(1:ncol), & + te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & + te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & + tw_cur = state%tw_cur(1:ncol), & + tend_te_tnd = tend%te_tnd(1:ncol), & + tend_tw_tnd = tend%tw_tnd(1:ncol), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & + count = state%count, & + teout = teout(1:ncol), & ! dummy argument - actual teout written to pbuf directly below + vc_physics = vc_physics, & + vc_dycore = vc_dycore, & + errmsg = errmsg, & + errflg = errflg & + ) + + ! initialize physics buffer + if (is_first_step()) then + call pbuf_set_field(pbuf, teout_idx, state%te_ini(:,dyn_te_idx), col_type=col_type) + end if + + end subroutine check_energy_timestep_init + + ! Check that the energy and water change matches the boundary fluxes + subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & + flx_vap, flx_cnd, flx_ice, flx_sen) + use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure + use cam_abortutils, only: endrun + use physics_types, only: phys_te_idx, dyn_te_idx + use physics_types, only: physics_tend + use physconst, only: cpair, rair, latice, latvap + use air_composition, only: cpairv, cp_or_cv_dycore + + ! CCPP-ized subroutine + use check_energy_chng, only: check_energy_chng_run + + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + character*(*),intent(in) :: name ! parameterization name for fluxes + integer , intent(in) :: nstep ! current timestep number + real(r8), intent(in) :: ztodt ! 2 delta t (model time increment) + real(r8), intent(in) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) + real(r8), intent(in) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?) + real(r8), intent(in) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?) + real(r8), intent(in) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2) + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + real(r8) :: local_cp_phys(state%psetcols,pver) + real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) + real(r8) :: scaling_dycore(state%ncol,pver) + character(len=512) :: errmsg + integer :: errflg + + lchnk = state%lchnk + ncol = state%ncol + + if(state%psetcols == pcols) then + ! No subcolumns + local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) + + ! Only if vertical coordinate is vc_height or vc_dry_pressure, cp_or_cv_dycore + ! is nonzero. + if(vc_dycore == vc_height .or. vc_dycore == vc_dry_pressure) then + local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) + + scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + endif + elseif(state%psetcols > pcols) then + ! Subcolumns + if(.not. all(cpairv(:,:,:) == cpair)) then + call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') + endif + + local_cp_phys(:,:) = cpair + + ! Note: cp_or_cv set above for pressure coordinate + if (vc_dycore == vc_height) then + ! compute cv if vertical coordinate is height: cv = cp - R + local_cp_or_cv_dycore(:ncol,:) = cpair-rair + scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + else if (vc_dycore == vc_dry_pressure) then + ! SE specific hydrostatic energy + local_cp_or_cv_dycore(:ncol,:) = cpair + scaling_dycore(:ncol,:) = 1.0_r8 + else + ! Moist pressure... use phys formula + local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:) + scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + end if + endif + + ! Call CCPP-ized underlying subroutine. + call check_energy_chng_run( & + ncol = ncol, & + pver = pver, & + pcnst = pcnst, & + q = state%q(1:ncol,1:pver,1:pcnst), & + pdel = state%pdel(1:ncol,1:pver), & + u = state%u(1:ncol,1:pver), & + v = state%v(1:ncol,1:pver), & + T = state%T(1:ncol,1:pver), & + pintdry = state%pintdry(1:ncol,1:pver), & + phis = state%phis(1:ncol), & + zm = state%zm(1:ncol,:), & + cp_phys = local_cp_phys(1:ncol,:), & + cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & + scaling_dycore = scaling_dycore(1:ncol,:), & + te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & + te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & + tw_cur = state%tw_cur(1:ncol), & + tend_te_tnd = tend%te_tnd(1:ncol), & + tend_tw_tnd = tend%tw_tnd(1:ncol), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & + count = state%count, & + ztodt = ztodt, & + latice = latice, & + latvap = latvap, & + vc_physics = vc_physics, & + vc_dycore = vc_dycore, & + name = name, & + flx_vap = flx_vap, & + flx_cnd = flx_cnd, & + flx_ice = flx_ice, & + flx_sen = flx_sen, & + errmsg = errmsg, & + errflg = errflg & + ) + + end subroutine check_energy_cam_chng + + ! Add heating rate required for global mean total energy conservation + subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) + use physics_types, only: physics_ptend, physics_ptend_init + use physconst, only: rga + + ! SCAM support + use scamMod, only: single_column, use_camiop, heat_glob_scm + use cam_history, only: write_camiop + use cam_history, only: outfld + + ! CCPP-ized subroutine + use check_energy_fix, only: check_energy_fix_run + + type(physics_state), intent(in) :: state + type(physics_ptend), intent(out) :: ptend + + integer , intent(in) :: nstep ! time step number + real(r8), intent(out) :: eshflx(pcols) ! effective sensible heat flux + + integer :: ncol ! number of atmospheric columns in chunk + integer :: lchnk ! chunk number + real(r8) :: heat_out(pcols) + + lchnk = state%lchnk + ncol = state%ncol + + call physics_ptend_init(ptend, state%psetcols, 'chkenergyfix', ls=.true.) + +#if ( defined OFFLINE_DYN ) + ! disable the energy fix for offline driver + heat_glob = 0._r8 +#endif + + ! Special handling of energy fix for SCAM - supplied via CAMIOP - zero's for normal IOPs + if (single_column) then + if (use_camiop) then + heat_glob = heat_glob_scm(1) + else + heat_glob = 0._r8 + endif + endif + + if (nstep > 0 .and. write_camiop) then + heat_out(:ncol) = heat_glob + call outfld('heat_glob', heat_out(:ncol), pcols, lchnk) + endif + + ! Call the CCPP-ized subroutine (for non-SCAM) + ! to compute the effective sensible heat flux and save to ptend%s + call check_energy_fix_run( & + ncol = ncol, & + pver = pver, & + pint = state%pint(:ncol,:), & + rga = rga, & + heat_glob = heat_glob, & + ptend_s = ptend%s(:ncol,:), & + eshflx = eshflx(:ncol) & + ) + + end subroutine check_energy_cam_fix end module check_energy diff --git a/src/physics/cam/check_energy_cam.F90 b/src/physics/cam/check_energy_cam.F90 deleted file mode 100644 index e1f143c0c2..0000000000 --- a/src/physics/cam/check_energy_cam.F90 +++ /dev/null @@ -1,380 +0,0 @@ -! CAM interface for CCPP-ized check_energy routines in CAM -module check_energy_cam - use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: pcols, pver - use constituents, only: pcnst - use physics_types, only: physics_state, physics_tend - - implicit none - private - - public :: check_energy_cam_chng ! check changes in integrals against cumulative boundary fluxes - public :: check_energy_cam_timestep_init ! timestep initialization of energy integrals and cumulative boundary fluxes - - public :: check_energy_cam_fix ! add heating rate required for global mean total energy conservation - - ! These subroutines cannot be CCPP-ized and replicated from check_energy - public :: check_energy_gmean ! global means of physics input and output total energy - public :: check_energy_get_integrals ! get energy integrals computed in check_energy_gmean - - ! Module variables used for check_energy_gmean - real(r8) :: teout_glob ! global mean energy of output state - real(r8) :: teinp_glob ! global mean energy of input state - real(r8) :: tedif_glob ! global mean energy difference - real(r8) :: psurf_glob ! global mean surface pressure - real(r8) :: ptopb_glob ! global mean top boundary pressure - real(r8) :: heat_glob ! global mean heating rate - -contains - ! Compute initial values of energy and water integrals, - ! zero cumulative tendencies - subroutine check_energy_cam_timestep_init(state, tend, pbuf, col_type) - use physics_buffer, only: physics_buffer_desc, pbuf_set_field - use cam_abortutils, only: endrun - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure - use physics_types, only: phys_te_idx, dyn_te_idx - use time_manager, only: is_first_step - use physconst, only: cpair, rair - use air_composition, only: cpairv, cp_or_cv_dycore - - ! To remove, pbuf indices in check_energy - use check_energy, only: teout_idx - - ! CCPP-ized subroutine - use check_energy_chng, only: check_energy_chng_timestep_init - - type(physics_state), intent(inout) :: state - type(physics_tend ), intent(inout) :: tend - type(physics_buffer_desc), pointer :: pbuf(:) - integer, optional :: col_type ! Flag indicating whether using grid or subcolumns - - real(r8) :: local_cp_phys(state%psetcols,pver) - real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) - real(r8) :: teout(state%ncol) ! dummy teout argument - integer :: lchnk ! chunk identifier - integer :: ncol ! number of atmospheric columns - character(len=512) :: errmsg - integer :: errflg - - lchnk = state%lchnk - ncol = state%ncol - - ! The code below is split into not-subcolumns and subcolumns code, as there is different handling of the - ! cp passed into the hydrostatic energy call. CAM-SIMA does not support subcolumns, so we keep this special - ! handling inside this shim module. (hplin, 9/9/24) - if(state%psetcols == pcols) then - ! No subcolumns - local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) - local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - else if (state%psetcols > pcols) then - ! Subcolumns code - ! Subcolumns specific error handling - if(.not. all(cpairv(:,:,lchnk) == cpair)) then - call endrun('check_energy_timestep_init: cpairv is not allowed to vary when subcolumns are turned on') - endif - - local_cp_phys(1:ncol,:) = cpair - - if (vc_dycore == vc_height) then - ! MPAS specific hydrostatic energy computation (internal energy) - local_cp_or_cv_dycore(:ncol,:) = cpair-rair - else if(vc_dycore == vc_dry_pressure) then - ! SE specific hydrostatic energy (enthalpy) - local_cp_or_cv_dycore(:ncol,:) = cpair - else - ! cp_or_cv is not used in the underlying subroutine, zero it out to be sure - local_cp_or_cv_dycore(:ncol,:) = 0.0_r8 - endif - end if - - ! Call CCPP-ized underlying subroutine. - call check_energy_chng_timestep_init( & - ncol = ncol, & - pver = pver, & - pcnst = pcnst, & - is_first_timestep = is_first_step(), & - q = state%q(1:ncol,1:pver,1:pcnst), & - pdel = state%pdel(1:ncol,1:pver), & - u = state%u(1:ncol,1:pver), & - v = state%v(1:ncol,1:pver), & - T = state%T(1:ncol,1:pver), & - pintdry = state%pintdry(1:ncol,1:pver), & - phis = state%phis(1:ncol), & - zm = state%zm(1:ncol,:), & - cp_phys = local_cp_phys(1:ncol,:), & - cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & - te_ini_phys = state%te_ini(1:ncol,phys_te_idx), & - te_ini_dyn = state%te_ini(1:ncol,dyn_te_idx), & - tw_ini = state%tw_ini(1:ncol), & - te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & - te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & - tw_cur = state%tw_cur(1:ncol), & - tend_te_tnd = tend%te_tnd(1:ncol), & - tend_tw_tnd = tend%tw_tnd(1:ncol), & - temp_ini = state%temp_ini(:ncol,:), & - z_ini = state%z_ini(:ncol,:), & - count = state%count, & - teout = teout(1:ncol), & ! dummy argument - actual teout written to pbuf directly below - vc_physics = vc_physics, & - vc_dycore = vc_dycore, & - errmsg = errmsg, & - errflg = errflg & - ) - - ! initialize physics buffer - if (is_first_step()) then - call pbuf_set_field(pbuf, teout_idx, state%te_ini(:,dyn_te_idx), col_type=col_type) - end if - - end subroutine check_energy_cam_timestep_init - - ! Check that the energy and water change matches the boundary fluxes - subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & - flx_vap, flx_cnd, flx_ice, flx_sen) - use cam_thermo, only: get_hydrostatic_energy - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure - use cam_abortutils, only: endrun - use physics_types, only: phys_te_idx, dyn_te_idx - use physconst, only: cpair, rair, latice, latvap - use air_composition, only: cpairv, cp_or_cv_dycore - - ! CCPP-ized subroutine - use check_energy_chng, only: check_energy_chng_run - - type(physics_state), intent(inout) :: state - type(physics_tend ), intent(inout) :: tend - character*(*),intent(in) :: name ! parameterization name for fluxes - integer , intent(in) :: nstep ! current timestep number - real(r8), intent(in) :: ztodt ! 2 delta t (model time increment) - real(r8), intent(in) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) - real(r8), intent(in) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?) - real(r8), intent(in) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?) - real(r8), intent(in) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2) - - integer :: lchnk ! chunk identifier - integer :: ncol ! number of atmospheric columns - integer :: i ! column index - real(r8) :: local_cp_phys(state%psetcols,pver) - real(r8) :: local_cp_or_cv_dycore(state%psetcols,pver) - real(r8) :: scaling_dycore(state%ncol,pver) - character(len=512) :: errmsg - integer :: errflg - - lchnk = state%lchnk - ncol = state%ncol - - if(state%psetcols == pcols) then - ! No subcolumns - local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) - local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) - - scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling - elseif(state%psetcols > pcols) then - ! Subcolumns - if(.not. all(cpairv(:,:,:) == cpair)) then - call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') - endif - - local_cp_phys(:,:) = cpair - - ! Note: cp_or_cv set above for pressure coordinate - if (vc_dycore == vc_height) then - ! compute cv if vertical coordinate is height: cv = cp - R - local_cp_or_cv_dycore(:ncol,:) = cpair-rair - scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling - else if (vc_dycore == vc_dry_pressure) then - ! SE specific hydrostatic energy - local_cp_or_cv_dycore(:ncol,:) = cpair - scaling_dycore(:ncol,:) = 1.0_r8 - else - ! Moist pressure... use phys formula - local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:) - scaling_dycore(:ncol,:) = cpairv(:,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling - end if - endif - - ! Call CCPP-ized underlying subroutine. - call check_energy_chng_run( & - ncol = ncol, & - pver = pver, & - pcnst = pcnst, & - q = state%q(1:ncol,1:pver,1:pcnst), & - pdel = state%pdel(1:ncol,1:pver), & - u = state%u(1:ncol,1:pver), & - v = state%v(1:ncol,1:pver), & - T = state%T(1:ncol,1:pver), & - pintdry = state%pintdry(1:ncol,1:pver), & - phis = state%phis(1:ncol), & - zm = state%zm(1:ncol,:), & - cp_phys = local_cp_phys(1:ncol,:), & - cp_or_cv_dycore = local_cp_or_cv_dycore(1:ncol,:), & - scaling_dycore = scaling_dycore(1:ncol,:), & - te_cur_phys = state%te_cur(1:ncol,phys_te_idx), & - te_cur_dyn = state%te_cur(1:ncol,dyn_te_idx), & - tw_cur = state%tw_cur(1:ncol), & - tend_te_tnd = tend%te_tnd(1:ncol), & - tend_tw_tnd = tend%tw_tnd(1:ncol), & - temp_ini = state%temp_ini(:ncol,:), & - z_ini = state%z_ini(:ncol,:), & - count = state%count, & - ztodt = ztodt, & - latice = latice, & - latvap = latvap, & - vc_physics = vc_physics, & - vc_dycore = vc_dycore, & - name = name, & - flx_vap = flx_vap, & - flx_cnd = flx_cnd, & - flx_ice = flx_ice, & - flx_sen = flx_sen, & - errmsg = errmsg, & - errflg = errflg & - ) - - end subroutine check_energy_cam_chng - - ! Add heating rate required for global mean total energy conservation - subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) - use physics_types, only: physics_ptend, physics_ptend_init - use physconst, only: rga - - ! SCAM support - use scamMod, only: single_column, use_camiop, heat_glob_scm - use cam_history, only: write_camiop - use cam_history, only: outfld - - ! CCPP-ized subroutine - use check_energy_fix, only: check_energy_fix_run - - type(physics_state), intent(in) :: state - type(physics_ptend), intent(out) :: ptend - - integer , intent(in) :: nstep ! time step number - real(r8), intent(out) :: eshflx(pcols) ! effective sensible heat flux - - integer :: ncol ! number of atmospheric columns in chunk - integer :: lchnk ! chunk number - real(r8) :: heat_out(pcols) - - lchnk = state%lchnk - ncol = state%ncol - - call physics_ptend_init(ptend, state%psetcols, 'chkenergyfix', ls=.true.) - -#if ( defined OFFLINE_DYN ) - ! disable the energy fix for offline driver - heat_glob = 0._r8 -#endif - - ! Special handling of energy fix for SCAM - supplied via CAMIOP - zero's for normal IOPs - if (single_column) then - if (use_camiop) then - heat_glob = heat_glob_scm(1) - else - heat_glob = 0._r8 - endif - endif - - if (nstep > 0 .and. write_camiop) then - heat_out(:ncol) = heat_glob - call outfld('heat_glob', heat_out(:ncol), pcols, lchnk) - endif - - ! Call the CCPP-ized subroutine (for non-SCAM) - ! to compute the effective sensible heat flux and save to ptend%s - call check_energy_fix_run( & - ncol = ncol, & - pver = pver, & - pint = state%pint(:ncol,:), & - rga = rga, & - heat_glob = heat_glob, & - ptend_s = ptend%s(:ncol,:), & - eshflx = eshflx(:ncol) & - ) - - end subroutine check_energy_cam_fix - - ! Compute global mean total energy of physics input and output states - ! computed consistently with dynamical core vertical coordinate - ! (under hydrostatic assumption) - ! - ! This subroutine cannot use the CCPP-ized equivalent because - ! it is dependent on chunks. - subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) - use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_chunk - use physics_types, only: dyn_te_idx - use ppgrid, only: begchunk, endchunk - use spmd_utils, only: masterproc - use cam_logfile, only: iulog - use gmean_mod, only: gmean - use physconst, only: gravit - - ! To remove, pbuf indices in check_energy - use check_energy, only: teout_idx - - type(physics_state), intent(in), dimension(begchunk:endchunk) :: state - type(physics_buffer_desc), pointer :: pbuf2d(:,:) - - real(r8), intent(in) :: dtime ! physics time step - integer , intent(in) :: nstep ! current timestep number - - integer :: ncol ! number of active columns - integer :: lchnk ! chunk index - - real(r8) :: te(pcols,begchunk:endchunk,4) - ! total energy of input/output states (copy) - real(r8) :: te_glob(4) ! global means of total energy - real(r8), pointer :: teout(:) - - ! Copy total energy out of input and output states - do lchnk = begchunk, endchunk - ncol = state(lchnk)%ncol - ! input energy using dynamical core energy formula - te(:ncol,lchnk,1) = state(lchnk)%te_ini(:ncol,dyn_te_idx) - ! output energy - call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout) - - te(:ncol,lchnk,2) = teout(1:ncol) - ! surface pressure for heating rate - te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1) - ! model top pressure for heating rate (not constant for z-based vertical coordinate!) - te(:ncol,lchnk,4) = state(lchnk)%pint(:ncol,1) - end do - - ! Compute global means of input and output energies and of - ! surface pressure for heating rate (assume uniform ptop) - call gmean(te, te_glob, 4) - - if (begchunk .le. endchunk) then - teinp_glob = te_glob(1) - teout_glob = te_glob(2) - psurf_glob = te_glob(3) - ptopb_glob = te_glob(4) - - ! Global mean total energy difference - tedif_glob = teinp_glob - teout_glob - heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob) - if (masterproc) then - write(iulog,'(1x,a9,1x,i8,5(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, & - heat_glob, psurf_glob, ptopb_glob - end if - else - heat_glob = 0._r8 - end if ! (begchunk .le. endchunk) - - end subroutine check_energy_gmean - - ! Return energy integrals (module variables) - subroutine check_energy_get_integrals(tedif_glob_out, heat_glob_out) - real(r8), intent(out), optional :: tedif_glob_out - real(r8), intent(out), optional :: heat_glob_out - - if ( present(tedif_glob_out) ) then - tedif_glob_out = tedif_glob - endif - - if ( present(heat_glob_out) ) then - heat_glob_out = heat_glob - endif - end subroutine check_energy_get_integrals -end module check_energy_cam diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 25f46c3a45..1d431b2dbb 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1065,7 +1065,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics,only: diag_allocate, diag_physvar_ic - use check_energy_cam, only: check_energy_gmean + use check_energy, only: check_energy_gmean use phys_control, only: phys_getopts use spcam_drivers, only: tphysbc_spcam use spmd_utils, only: mpicom @@ -1392,7 +1392,7 @@ subroutine tphysac (ztodt, cam_in, & use carma_flags_mod, only: carma_do_aerosol, carma_do_emission use check_energy, only: tot_energy_phys use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng - use check_energy_cam, only: check_energy_cam_chng + use check_energy, only: check_energy_cam_chng use time_manager, only: get_nstep use cam_abortutils, only: endrun use dycore, only: dycore_is @@ -2045,8 +2045,8 @@ subroutine tphysbc (ztodt, state, & use convect_deep, only: convect_deep_tend, convect_deep_tend_2, deep_scheme_does_scav_trans use time_manager, only: is_first_step, get_nstep use convect_shallow, only: convect_shallow_tend - use check_energy_cam,only: check_energy_cam_timestep_init, check_energy_cam_chng - use check_energy_cam,only: check_energy_cam_fix + use check_energy, only: check_energy_timestep_init, check_energy_cam_chng + use check_energy, only: check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use dycore, only: dycore_is @@ -2257,7 +2257,9 @@ subroutine tphysbc (ztodt, state, & call tot_energy_phys(state, 'dyBF',vc=vc_dycore) if (.not.dycore_is('EUL')) then call check_energy_cam_fix(state, ptend, nstep, flx_heat) + call physics_update(state, ptend, ztodt, tend) + call check_energy_cam_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) call outfld( 'EFIX', flx_heat , pcols, lchnk ) end if @@ -2632,7 +2634,7 @@ subroutine tphysbc (ztodt, state, & call subcol_gen(state, tend, state_sc, tend_sc, pbuf) !Initialize check energy for subcolumns - call check_energy_cam_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) + call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) end if if (trim(cam_take_snapshot_before) == "microp_section") then diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index 1ce28d3a93..ead4935939 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -1056,7 +1056,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics,only: diag_allocate, diag_physvar_ic - use check_energy_cam, only: check_energy_gmean + use check_energy, only: check_energy_gmean use spmd_utils, only: mpicom use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate use cam_history, only: outfld, write_camiop @@ -1366,7 +1366,7 @@ subroutine tphysac (ztodt, cam_in, & use aoa_tracers, only: aoa_tracers_timestep_tend use physconst, only: rhoh2o use aero_model, only: aero_model_drydep - use check_energy_cam, only: check_energy_cam_timestep_init, check_energy_cam_chng + use check_energy, only: check_energy_timestep_init, check_energy_cam_chng use check_energy, only: tot_energy_phys use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use time_manager, only: get_nstep @@ -1775,7 +1775,7 @@ subroutine tphysac (ztodt, cam_in, & call subcol_gen(state, tend, state_sc, tend_sc, pbuf) !Initialize check energy for subcolumns - call check_energy_cam_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) + call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) end if if (trim(cam_take_snapshot_before) == "microp_section") then @@ -2517,7 +2517,7 @@ subroutine tphysbc (ztodt, state, & use convect_deep, only: convect_deep_tend use time_manager, only: is_first_step, get_nstep use convect_diagnostics,only: convect_diagnostics_calc - use check_energy_cam,only: check_energy_cam_chng, check_energy_cam_fix + use check_energy, only: check_energy_cam_chng, check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init use check_energy, only: tot_energy_phys use dycore, only: dycore_is diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index eaea0ff1dd..c236f7f021 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -305,7 +305,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) !----------------------------------------------------------------------- use time_manager, only: get_nstep use cam_diagnostics, only: diag_allocate, diag_physvar_ic - use check_energy_cam, only: check_energy_gmean + use check_energy, only: check_energy_gmean use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate ! @@ -503,7 +503,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) use air_composition, only: cpairv, cp_or_cv_dycore use time_manager, only: get_nstep use nudging, only: Nudge_Model, Nudge_ON, nudging_timestep_tend - use check_energy_cam,only: check_energy_cam_chng + use check_energy, only: check_energy_cam_chng ! Arguments ! @@ -728,8 +728,8 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) use cam_diagnostics, only: diag_conv_tend_ini, diag_conv, diag_export use cam_history, only: outfld use time_manager, only: get_nstep - use check_energy_cam, only: check_energy_cam_chng, check_energy_cam_fix - use check_energy_cam, only: check_energy_cam_timestep_init + use check_energy, only: check_energy_cam_chng, check_energy_cam_fix + use check_energy, only: check_energy_timestep_init use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng use check_energy, only: tot_energy_phys use chemistry, only: chem_is_active, chem_timestep_tend diff --git a/src/physics/spcam/crm_physics.F90 b/src/physics/spcam/crm_physics.F90 index 735d6c3630..8812d2be72 100644 --- a/src/physics/spcam/crm_physics.F90 +++ b/src/physics/spcam/crm_physics.F90 @@ -750,7 +750,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf, cam_in) use crmx_crm_module, only: crm use crmx_microphysics, only: nmicro_fields use physconst, only: latvap - use check_energy_cam, only: check_energy_cam_chng + use check_energy, only: check_energy_cam_chng use phys_grid, only: get_rlat_all_p, get_rlon_all_p, get_lon_all_p, get_lat_all_p use modal_aero_calcsize, only: modal_aero_calcsize_sub use micro_pumas_utils, only: size_dist_param_liq, mg_liq_props, mincld, qsmall diff --git a/src/physics/spcam/spcam_drivers.F90 b/src/physics/spcam/spcam_drivers.F90 index cd7d3f6370..b9f7a596cc 100644 --- a/src/physics/spcam/spcam_drivers.F90 +++ b/src/physics/spcam/spcam_drivers.F90 @@ -300,7 +300,7 @@ subroutine tphysbc_spcam (ztodt, state, & use cam_history, only: outfld use constituents, only: pcnst, qmin, cnst_get_ind use time_manager, only: get_nstep - use check_energy_cam,only: check_energy_cam_chng, check_energy_cam_fix + use check_energy, only: check_energy_cam_chng, check_energy_cam_fix use check_energy, only: check_tracers_data, check_tracers_init use dycore, only: dycore_is use radiation, only: radiation_tend From 4b6de4306cce342c6eab46333004e63d4afcb989 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 15 Oct 2024 12:47:16 -0400 Subject: [PATCH 09/20] Fix subcol test with FV dycore --- src/physics/cam/check_energy.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index f318cb51e1..998fbb9e18 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -801,9 +801,9 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & local_cp_or_cv_dycore(:ncol,:) = cpair scaling_dycore(:ncol,:) = 1.0_r8 else - ! Moist pressure... use phys formula - local_cp_or_cv_dycore(:ncol,:) = local_cp_phys(:ncol,:) - scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling + ! Moist pressure... use phys formula, cp_or_cv_dycore is unused. Reset for safety + local_cp_or_cv_dycore(:ncol,:) = 0.0_r8 + scaling_dycore(:ncol,:) = 0.0_r8 end if endif From d924e4e434929aa47a89d875e80b7974ce1cdee5 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 23 Oct 2024 13:19:24 -0400 Subject: [PATCH 10/20] Update atmos_phys; move energy formula to cam_thermo_formula; write cp_or_cv_dycore to snapshot --- .gitmodules | 2 +- src/atmos_phys | 2 +- src/control/cam_snapshot_common.F90 | 9 ++++- src/physics/cam/check_energy.F90 | 56 +++++++++++++++-------------- src/utils/cam_thermo_formula.F90 | 14 ++++++++ 5 files changed, 53 insertions(+), 30 deletions(-) create mode 100644 src/utils/cam_thermo_formula.F90 diff --git a/.gitmodules b/.gitmodules index e295d8fec1..0edbb78972 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/jimmielin/atmospheric_physics - fxtag = df80b9c1 + fxtag = 952ebdd fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index df80b9c1e7..952ebddfa2 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit df80b9c1e73064f9a0196282a92b66a2041f6a57 +Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 diff --git a/src/control/cam_snapshot_common.F90 b/src/control/cam_snapshot_common.F90 index ffae561370..fa558b3ad2 100644 --- a/src/control/cam_snapshot_common.F90 +++ b/src/control/cam_snapshot_common.F90 @@ -81,7 +81,7 @@ module cam_snapshot_common integer :: cam_snapshot_before_num, cam_snapshot_after_num ! Note the maximum number of variables for each type -type (snapshot_type) :: state_snapshot(29) +type (snapshot_type) :: state_snapshot(30) type (snapshot_type) :: cnst_snapshot(pcnst) type (snapshot_type) :: tend_snapshot(6) type (snapshot_type) :: cam_in_snapshot(30) @@ -283,6 +283,9 @@ subroutine cam_state_snapshot_init(cam_snapshot_before_num_in, cam_snapshot_afte call snapshot_addfld( nstate_var, state_snapshot, cam_snapshot_before_num, cam_snapshot_after_num, & 'state%te_cur_dyn', 'state_te_cur_dyn', 'unset', horiz_only) + call snapshot_addfld( nstate_var, state_snapshot, cam_snapshot_before_num, cam_snapshot_after_num, & + 'air_composition_cp_or_cv_dycore', 'cp_or_cv_dycore', 'J kg-1 K-1', 'lev') + end subroutine cam_state_snapshot_init subroutine cam_cnst_snapshot_init(cam_snapshot_before_num, cam_snapshot_after_num) @@ -741,6 +744,7 @@ end subroutine snapshot_addfld subroutine state_snapshot_all_outfld(lchnk, file_num, state) use physics_types, only: phys_te_idx, dyn_te_idx + use air_composition, only: cp_or_cv_dycore integer, intent(in) :: lchnk integer, intent(in) :: file_num @@ -843,6 +847,9 @@ subroutine state_snapshot_all_outfld(lchnk, file_num, state) case ('state%te_cur_dyn') call outfld(state_snapshot(i)%standard_name, state%te_cur(:, dyn_te_idx), pcols, lchnk) + case ('air_composition_cp_or_cv_dycore') + call outfld(state_snapshot(i)%standard_name, cp_or_cv_dycore(:,:,lchnk), pcols, lchnk) + case default call endrun('ERROR in state_snapshot_all_outfld: no match found for '//trim(state_snapshot(i)%ddt_string)) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 998fbb9e18..774764c563 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -415,7 +415,8 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) use cam_thermo, only: get_hydrostatic_energy,thermo_budget_num_vars,thermo_budget_vars, & wvidx,wlidx,wiidx,seidx,poidx,keidx,moidx,mridx,ttidx,teidx use cam_history, only: outfld - use dyn_tests_utils, only: vc_physics, vc_height, vc_dry_pressure + use dyn_tests_utils, only: vc_physics + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS use cam_abortutils, only: endrun use cam_history_support, only: max_fieldname_len @@ -467,7 +468,7 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) end if if (state%psetcols == pcols) then - if (vc_loc == vc_height .or. vc_loc == vc_dry_pressure) then + if (vc_loc == ENERGY_FORMULA_DYCORE_MPAS .or. vc_loc == ENERGY_FORMULA_DYCORE_SE) then cp_or_cv(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) else cp_or_cv(:ncol,:) = cpairv(:ncol,:,lchnk) @@ -476,7 +477,7 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) call endrun('tot_energy_phys: energy diagnostics not implemented/tested for subcolumns') end if - if (vc_loc == vc_height .or. vc_loc == vc_dry_pressure) then + if (vc_loc == ENERGY_FORMULA_DYCORE_MPAS .or. vc_loc == ENERGY_FORMULA_DYCORE_SE) then scaling(:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv(:ncol,:)!scaling for energy consistency else scaling(:ncol,:) = 1.0_r8 !internal energy / enthalpy same as CAM physics @@ -643,7 +644,8 @@ end subroutine check_energy_get_integrals subroutine check_energy_timestep_init(state, tend, pbuf, col_type) use physics_buffer, only: physics_buffer_desc, pbuf_set_field use cam_abortutils, only: endrun - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure + use dyn_tests_utils, only: vc_physics, vc_dycore + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS use physics_types, only: physics_tend use physics_types, only: phys_te_idx, dyn_te_idx use time_manager, only: is_first_step @@ -671,7 +673,7 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type) ! The code below is split into not-subcolumns and subcolumns code, as there is different handling of the ! cp passed into the hydrostatic energy call. CAM-SIMA does not support subcolumns, so we keep this special - ! handling inside this shim module. (hplin, 9/9/24) + ! handling inside this CAM interface. (hplin, 9/9/24) if(state%psetcols == pcols) then ! No subcolumns local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) @@ -685,10 +687,10 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type) local_cp_phys(1:ncol,:) = cpair - if (vc_dycore == vc_height) then + if (vc_dycore == ENERGY_FORMULA_DYCORE_MPAS) then ! MPAS specific hydrostatic energy computation (internal energy) local_cp_or_cv_dycore(:ncol,:) = cpair-rair - else if(vc_dycore == vc_dry_pressure) then + else if(vc_dycore == ENERGY_FORMULA_DYCORE_SE) then ! SE specific hydrostatic energy (enthalpy) local_cp_or_cv_dycore(:ncol,:) = cpair else @@ -721,12 +723,12 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type) tw_cur = state%tw_cur(1:ncol), & tend_te_tnd = tend%te_tnd(1:ncol), & tend_tw_tnd = tend%tw_tnd(1:ncol), & - temp_ini = state%temp_ini(:ncol,:), & - z_ini = state%z_ini(:ncol,:), & + temp_ini = state%temp_ini(:ncol,:), & + z_ini = state%z_ini(:ncol,:), & count = state%count, & teout = teout(1:ncol), & ! dummy argument - actual teout written to pbuf directly below - vc_physics = vc_physics, & - vc_dycore = vc_dycore, & + energy_formula_physics = vc_physics, & + energy_formula_dycore = vc_dycore, & errmsg = errmsg, & errflg = errflg & ) @@ -741,15 +743,16 @@ end subroutine check_energy_timestep_init ! Check that the energy and water change matches the boundary fluxes subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & flx_vap, flx_cnd, flx_ice, flx_sen) - use dyn_tests_utils, only: vc_physics, vc_dycore, vc_height, vc_dry_pressure - use cam_abortutils, only: endrun - use physics_types, only: phys_te_idx, dyn_te_idx - use physics_types, only: physics_tend - use physconst, only: cpair, rair, latice, latvap - use air_composition, only: cpairv, cp_or_cv_dycore + use dyn_tests_utils, only: vc_physics, vc_dycore + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + use cam_abortutils, only: endrun + use physics_types, only: phys_te_idx, dyn_te_idx + use physics_types, only: physics_tend + use physconst, only: cpair, rair, latice, latvap + use air_composition, only: cpairv, cp_or_cv_dycore ! CCPP-ized subroutine - use check_energy_chng, only: check_energy_chng_run + use check_energy_chng, only: check_energy_chng_run type(physics_state), intent(inout) :: state type(physics_tend ), intent(inout) :: tend @@ -776,9 +779,8 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & ! No subcolumns local_cp_phys(:ncol,:) = cpairv(:ncol,:,lchnk) - ! Only if vertical coordinate is vc_height or vc_dry_pressure, cp_or_cv_dycore - ! is nonzero. - if(vc_dycore == vc_height .or. vc_dycore == vc_dry_pressure) then + ! Only if using MPAS or SE energy formula cp_or_cv_dycore is nonzero. + if(vc_dycore == ENERGY_FORMULA_DYCORE_MPAS .or. vc_dycore == ENERGY_FORMULA_DYCORE_SE) then local_cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:,lchnk) scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling @@ -792,11 +794,11 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & local_cp_phys(:,:) = cpair ! Note: cp_or_cv set above for pressure coordinate - if (vc_dycore == vc_height) then + if (vc_dycore == ENERGY_FORMULA_DYCORE_MPAS) then ! compute cv if vertical coordinate is height: cv = cp - R local_cp_or_cv_dycore(:ncol,:) = cpair-rair scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling - else if (vc_dycore == vc_dry_pressure) then + else if (vc_dycore == ENERGY_FORMULA_DYCORE_SE) then ! SE specific hydrostatic energy local_cp_or_cv_dycore(:ncol,:) = cpair scaling_dycore(:ncol,:) = 1.0_r8 @@ -832,10 +834,10 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & z_ini = state%z_ini(:ncol,:), & count = state%count, & ztodt = ztodt, & - latice = latice, & - latvap = latvap, & - vc_physics = vc_physics, & - vc_dycore = vc_dycore, & + latice = latice, & + latvap = latvap, & + energy_formula_physics = vc_physics, & + energy_formula_dycore = vc_dycore, & name = name, & flx_vap = flx_vap, & flx_cnd = flx_cnd, & diff --git a/src/utils/cam_thermo_formula.F90 b/src/utils/cam_thermo_formula.F90 new file mode 100644 index 0000000000..7781e9da9c --- /dev/null +++ b/src/utils/cam_thermo_formula.F90 @@ -0,0 +1,14 @@ +module cam_thermo_formula + + implicit none + private + save + + ! energy_formula options for use by CCPPized check_energy + integer, public, parameter :: ENERGY_FORMULA_DYCORE_FV = 0 ! vc_moist_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_SE = 1 ! vc_dry_pressure + integer, public, parameter :: ENERGY_FORMULA_DYCORE_MPAS = 2 ! vc_height + + !REMOVECAM: in CAM, energy_formula_physics and energy_formula_dycore still uses vc_physics + ! and vc_dycore in dyn_tests_utils. The values are the same. +end module cam_thermo_formula From 46b82279c76644b410979b327dabb813afee01a7 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 28 Oct 2024 11:09:15 -0400 Subject: [PATCH 11/20] Update external to fix bug; ChangeLog draft --- .gitmodules | 2 +- doc/ChangeLog | 72 +++++++++++++++++++++++++++++++++++++++++++++++--- src/atmos_phys | 2 +- 3 files changed, 71 insertions(+), 5 deletions(-) diff --git a/.gitmodules b/.gitmodules index 0edbb78972..e74baf681e 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/jimmielin/atmospheric_physics - fxtag = 952ebdd + fxtag = 44d724d fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/doc/ChangeLog b/doc/ChangeLog index 8294c223f2..f9aded72d5 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,71 @@ =============================================================== +Tag name: +Originator(s): jimmielin +Date: 25 Oct 2024 +One-line Summary: Implement CCPPized check_energy_chng and check_energy_fix +Github PR URL: https://github.com/ESCOMP/CAM/pull/1180 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): +- CCPPize check_energy_chng and check_energy_fix (https://github.com/ESCOMP/CAM/issues/1138) + CAM interfaces have been kept in check_energy.F90 instead of a new module because check_tracers is not yet CCPPized, there are other non-CCPPized routines in module, and for compatibility with FV3 external calls +- Save air_composition cp_or_cv_dycore into state snapshot +- Separate out "energy_formula_physics"/"energy_formula_dycore" definitions used in get_hydrostatic_energy from dyn_tests_utils "vcoord" for an eventual change in SIMA. + +Describe any changes made to build system: N/A + +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: + +List all files eliminated: N/A + +List all files added and what they do: +- separate out energy formula definitions used in get_hydrostatic_energy from dyn_tests_utils "vcoord" +A src/utils/cam_thermo_formula.F90 + +List all existing files that have been modified, and describe the changes: + +Update atmos_phys external with CCPPized check_energy +M .gitmodules +M bld/configure +M src/atmos_phys + +Save cp_or_cv_dycore in state snapshot +M src/control/cam_snapshot_common.F90 + +New CAM interfaces to CCPPized routines, when avaialble +M src/physics/cam/check_energy.F90 + +Update calls to CCPPized routines, when available +M src/physics/cam/physpkg.F90 +M src/physics/cam/zm_conv_intr.F90 +M src/physics/cam7/physpkg.F90 +M src/physics/simple/physpkg.F90 +M src/physics/spcam/crm_physics.F90 +M src/physics/spcam/spcam_drivers.F90 + +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. + +derecho/intel/aux_cam: + +derecho/nvhpc/aux_cam: + +izumi/nag/aux_cam: + +izumi/gnu/aux_cam: + +Summarize any changes to answers: B4B + +=============================================================== + Tag name: cam6_4_043 Originator(s): eaton Date: 25 Oct 2024 @@ -31,12 +97,12 @@ Describe any changes made to the namelist: . change default value of seasalt_emis_scale to 0.75 for cam7 (both lt and mt) This is a cam7 tuning mod from issue #1143 -. update ubc_file_path for cam7 (lt only) to +. update ubc_file_path for cam7 (lt only) to atm/cam/chem/ubc/b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensAvg123.cam.h0zm.H2O.1849-2014_c240604.nc List any changes to the defaults for the boundary datasets: -. update ubc_file_path for cam7 (lt only) to +. update ubc_file_path for cam7 (lt only) to atm/cam/chem/ubc/b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensAvg123.cam.h0zm.H2O.1849-2014_c240604.nc Describe any substantial timing or memory changes: none @@ -64,7 +130,7 @@ bld/namelist_files/namelist_defaults_cam.xml This is a cam7 tuning mod from issue #1143 bld/namelist_files/use_cases/1850_cam_lt.xml -. update ubc_file_path to +. update ubc_file_path to atm/cam/chem/ubc/b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensAvg123.cam.h0zm.H2O.1849-2014_c240604.nc cime_config/testdefs/testlist_cam.xml diff --git a/src/atmos_phys b/src/atmos_phys index 952ebddfa2..44d724d8e8 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit 952ebddfa22dd796578fba8d73db6128c8db88c1 +Subproject commit 44d724d8e88ab0b910543a9fdecd24fea9ca62f7 From 7244285a9702f8ec99478bc5dc9d6991a634826a Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 18 Nov 2024 11:26:46 -0500 Subject: [PATCH 12/20] Update build issues; atmos_phys external to head of PR --- .gitmodules | 2 +- src/atmos_phys | 2 +- src/physics/cam/check_energy.F90 | 5 +++-- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.gitmodules b/.gitmodules index 195d003bc2..c478c0c686 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/jimmielin/atmospheric_physics - fxtag = 7b188e35 + fxtag = b1e2c2a9 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index 7b188e35fe..b1e2c2a9d1 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit 7b188e35fe26247289d3bdbbbdd902cde382765a +Subproject commit b1e2c2a9d188172f888b0b5c903c5a12a741f30e diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 774764c563..9fbc2d592c 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -814,6 +814,7 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & ncol = ncol, & pver = pver, & pcnst = pcnst, & + iulog = iulog, & q = state%q(1:ncol,1:pver,1:pcnst), & pdel = state%pdel(1:ncol,1:pver), & u = state%u(1:ncol,1:pver), & @@ -852,7 +853,7 @@ end subroutine check_energy_cam_chng ! Add heating rate required for global mean total energy conservation subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) use physics_types, only: physics_ptend, physics_ptend_init - use physconst, only: rga + use physconst, only: gravit ! SCAM support use scamMod, only: single_column, use_camiop, heat_glob_scm @@ -902,7 +903,7 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) ncol = ncol, & pver = pver, & pint = state%pint(:ncol,:), & - rga = rga, & + gravit = gravit, & heat_glob = heat_glob, & ptend_s = ptend%s(:ncol,:), & eshflx = eshflx(:ncol) & From e3e6a131208feaf72360d5fa5767ddf638ecc3c9 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 18 Nov 2024 11:52:16 -0500 Subject: [PATCH 13/20] Update check_energy_cam_fix interface --- src/physics/cam/check_energy.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 9fbc2d592c..592faeab00 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -872,6 +872,7 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) integer :: ncol ! number of atmospheric columns in chunk integer :: lchnk ! chunk number real(r8) :: heat_out(pcols) + character(len=64) :: dummy_scheme_name ! dumy scheme name for CCPP-ized scheme lchnk = state%lchnk ncol = state%ncol @@ -906,7 +907,8 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) gravit = gravit, & heat_glob = heat_glob, & ptend_s = ptend%s(:ncol,:), & - eshflx = eshflx(:ncol) & + eshflx = eshflx(:ncol), & + scheme_name = dummy_scheme_name ) end subroutine check_energy_cam_fix From a10a9557acc06679839e64682c9c3f5f4b96663b Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 18 Nov 2024 13:15:13 -0500 Subject: [PATCH 14/20] Update check_energy_cam_fix interface (2) --- src/physics/cam/check_energy.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 592faeab00..3f203bae2c 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -872,7 +872,7 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) integer :: ncol ! number of atmospheric columns in chunk integer :: lchnk ! chunk number real(r8) :: heat_out(pcols) - character(len=64) :: dummy_scheme_name ! dumy scheme name for CCPP-ized scheme + character(len=64) :: dummy_scheme_name ! dummy scheme name for CCPP-ized scheme lchnk = state%lchnk ncol = state%ncol @@ -908,7 +908,7 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) heat_glob = heat_glob, & ptend_s = ptend%s(:ncol,:), & eshflx = eshflx(:ncol), & - scheme_name = dummy_scheme_name + scheme_name = dummy_scheme_name & ) end subroutine check_energy_cam_fix From 902c495ebb7ec85b107ed7cba8bf003136b9ca04 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 19 Nov 2024 12:26:51 -0500 Subject: [PATCH 15/20] Update atmos_phys to atmos_phys0_07_000 --- .gitmodules | 4 ++-- src/atmos_phys | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index c478c0c686..40e189dab5 100644 --- a/.gitmodules +++ b/.gitmodules @@ -35,8 +35,8 @@ [submodule "atmos_phys"] path = src/atmos_phys - url = https://github.com/jimmielin/atmospheric_physics - fxtag = b1e2c2a9 + url = https://github.com/ESCOMP/atmospheric_physics + fxtag = atmos_phys0_07_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index b1e2c2a9d1..e10f81151e 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit b1e2c2a9d188172f888b0b5c903c5a12a741f30e +Subproject commit e10f81151ecef5416462472ef0e52810dc82e412 From bc5428f4eaaf55641176b61034a1a71363a943cd Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 20 Nov 2024 21:41:42 -0500 Subject: [PATCH 16/20] Address review comments --- src/control/cam_snapshot_common.F90 | 2 ++ src/physics/cam/check_energy.F90 | 15 +++++++++------ 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/control/cam_snapshot_common.F90 b/src/control/cam_snapshot_common.F90 index fa558b3ad2..0f49d27a02 100644 --- a/src/control/cam_snapshot_common.F90 +++ b/src/control/cam_snapshot_common.F90 @@ -848,6 +848,8 @@ subroutine state_snapshot_all_outfld(lchnk, file_num, state) call outfld(state_snapshot(i)%standard_name, state%te_cur(:, dyn_te_idx), pcols, lchnk) case ('air_composition_cp_or_cv_dycore') + ! this field is not part of physics state (it is in air_composition) + ! but describes the atmospheric thermodynamic state and thus saved within the snapshot call outfld(state_snapshot(i)%standard_name, cp_or_cv_dycore(:,:,lchnk), pcols, lchnk) case default diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 3f203bae2c..5350371cef 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -461,6 +461,9 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) lchnk = state%lchnk ncol = state%ncol + ! The "vertical coordinate" parameter is equivalent to the dynamical core + ! energy formula parameter, which controls the dycore energy formula used + ! by get_hydrostatic_energy. if (present(vc)) then vc_loc = vc else @@ -758,11 +761,11 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & type(physics_tend ), intent(inout) :: tend character*(*),intent(in) :: name ! parameterization name for fluxes integer , intent(in) :: nstep ! current timestep number - real(r8), intent(in) :: ztodt ! 2 delta t (model time increment) - real(r8), intent(in) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) - real(r8), intent(in) :: flx_cnd(:) ! (pcols) -boundary flux of liquid+ice (m/s) (precip?) - real(r8), intent(in) :: flx_ice(:) ! (pcols) -boundary flux of ice (m/s) (snow?) - real(r8), intent(in) :: flx_sen(:) ! (pcols) -boundary flux of sensible heat (w/m2) + real(r8), intent(in) :: ztodt ! physics timestep (s) + real(r8), intent(in) :: flx_vap(:) ! (pcols) - boundary flux of vapor (kg/m2/s) + real(r8), intent(in) :: flx_cnd(:) ! (pcols) - boundary flux of lwe liquid+ice (m/s) + real(r8), intent(in) :: flx_ice(:) ! (pcols) - boundary flux of lwe ice (m/s) + real(r8), intent(in) :: flx_sen(:) ! (pcols) - boundary flux of sensible heat (W/m2) integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns @@ -785,7 +788,7 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & scaling_dycore(:ncol,:) = cpairv(:ncol,:,lchnk)/local_cp_or_cv_dycore(:ncol,:) ! cp/cv scaling endif - elseif(state%psetcols > pcols) then + else if(state%psetcols > pcols) then ! Subcolumns if(.not. all(cpairv(:,:,:) == cpair)) then call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') From c0b9d8bc47d2a7b658d1b45a31e723f5ca7055a3 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 21 Nov 2024 08:51:05 -0500 Subject: [PATCH 17/20] Update comment for tot_energy_phys vc parameter --- src/physics/cam/check_energy.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 5350371cef..7af3a32bda 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -425,7 +425,7 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc) type(physics_state), intent(inout) :: state character(len=*), intent(in) :: outfld_name_suffix ! suffix for "outfld" - integer, optional, intent(in) :: vc ! vertical coordinate + integer, optional, intent(in) :: vc ! vertical coordinate (controls energy formula to use) !---------------------------Local storage------------------------------- real(r8) :: se(pcols) ! Dry Static energy (J/m2) From bdb3b768da79bebb4154a8e1a4597883de700261 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 30 Dec 2024 23:54:15 -0500 Subject: [PATCH 18/20] Update to head of atmos_phys with check_energy_fix_run arguments --- .gitmodules | 2 +- src/atmos_phys | 2 +- src/physics/cam/check_energy.F90 | 21 +++++++++++++-------- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/.gitmodules b/.gitmodules index 91690771e1..d07d5f8ac5 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_07_000 + fxtag = 8e3d0bdff98ff70a4e815160ebd13a150c662a12 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index e10f81151e..8e3d0bdff9 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit e10f81151ecef5416462472ef0e52810dc82e412 +Subproject commit 8e3d0bdff98ff70a4e815160ebd13a150c662a12 diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 7af3a32bda..d1d59e173f 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -877,6 +877,9 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) real(r8) :: heat_out(pcols) character(len=64) :: dummy_scheme_name ! dummy scheme name for CCPP-ized scheme + integer :: errflg + character(len=512) :: errmsg + lchnk = state%lchnk ncol = state%ncol @@ -904,14 +907,16 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) ! Call the CCPP-ized subroutine (for non-SCAM) ! to compute the effective sensible heat flux and save to ptend%s call check_energy_fix_run( & - ncol = ncol, & - pver = pver, & - pint = state%pint(:ncol,:), & - gravit = gravit, & - heat_glob = heat_glob, & - ptend_s = ptend%s(:ncol,:), & - eshflx = eshflx(:ncol), & - scheme_name = dummy_scheme_name & + ncol = ncol, & + pver = pver, & + pint = state%pint(:ncol,:), & + gravit = gravit, & + heat_glob = heat_glob, & + ptend_s = ptend%s(:ncol,:), & + eshflx = eshflx(:ncol), & + scheme_name = dummy_scheme_name, & + errmsg = errmsg, & + errflg = errflg & ) end subroutine check_energy_cam_fix From 0881a22bde02cf059f3d383f0de13cf35a492108 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 31 Dec 2024 09:57:33 -0500 Subject: [PATCH 19/20] Update to atmos_phys0_07_001 --- .gitmodules | 2 +- src/atmos_phys | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 29c2e7bdac..1809a94e55 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/ESCOMP/atmospheric_physics - fxtag = 8e3d0bdff98ff70a4e815160ebd13a150c662a12 + fxtag = atmos_phys0_07_001 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/src/atmos_phys b/src/atmos_phys index 8e3d0bdff9..c3de8468f7 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit 8e3d0bdff98ff70a4e815160ebd13a150c662a12 +Subproject commit c3de8468f7b245a939448f4ca6d3ef386584e92d From b7b4c3cfad202facbb6064e33f607c26f1988f03 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 31 Dec 2024 14:02:57 -0500 Subject: [PATCH 20/20] Finalize ChangeLog for cam6_4_050. --- doc/ChangeLog | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/doc/ChangeLog b/doc/ChangeLog index dcc4560393..8ae657fda9 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,7 +1,7 @@ =============================================================== -Tag name: +Tag name: cam6_4_050 Originator(s): jimmielin -Date: 30 Dec 2024 +Date: 31 Dec 2024 One-line Summary: Implement CCPPized check_energy_chng and check_energy_fix Github PR URL: https://github.com/ESCOMP/CAM/pull/1180 @@ -37,7 +37,7 @@ M src/atmos_phys Save cp_or_cv_dycore in state snapshot M src/control/cam_snapshot_common.F90 -New CAM interfaces to CCPPized routines, when avaialble +New CAM interfaces to CCPPized routines, when available M src/physics/cam/check_energy.F90 Update calls to CCPPized routines, when available @@ -55,11 +55,22 @@ appropriate machine below. All failed tests must be justified. derecho/intel/aux_cam: -derecho/nvhpc/aux_cam: +SMS_D_Ln9.f19_f19_mg17.FXHIST.derecho_intel.cam-outfrq9s_amie (Overall: FAIL) +SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL) + - pre-existing failures due to build-namelist error requiring CLM/CTSM external update. + +ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL) +SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: PASS) + - pre-existing failure due to HEMCO not having reproducible results issues #1018 and #856 + +derecho/nvhpc/aux_cam: All PASS izumi/nag/aux_cam: -izumi/gnu/aux_cam: +DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) + - pre-existing failure -- issue #670 + +izumi/gnu/aux_cam: All PASS Summarize any changes to answers: B4B