From 09d3f295ea16d2eb329e4e54df8558656aee0af4 Mon Sep 17 00:00:00 2001 From: Kai Huang Date: Wed, 11 Mar 2026 15:17:23 -0600 Subject: [PATCH 1/4] Kai's ZM_MCSP modifications --- bld/namelist_files/namelist_defaults_cam.xml | 6 + bld/namelist_files/namelist_definition.xml | 29 +- src/control/runtime_opts.F90 | 8 +- src/physics/cam/zm_conv_MCSP.F90 | 464 +++++++++++++++++++ src/physics/cam/zm_conv_intr.F90 | 72 ++- 5 files changed, 575 insertions(+), 4 deletions(-) create mode 100644 src/physics/cam/zm_conv_MCSP.F90 diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 75acdede44..e86e93a0dc 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -2891,6 +2891,12 @@ .false. 0.5 + + 0.0D0 + 0.0D0 + 0.0D0 + 0.0D0 + 1.0D0 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index d767872618..7e041b8a74 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3508,8 +3508,6 @@ Number of subcolumns in VAMP Generator Default: 10 - - + + + +MCSP heating coefficient controlling the tendency intensity added. +Default: 0.0 + + + +MCSP moistening coefficient controlling the tendency intensity added. +Default: 0.0 + + + +MCSP uwnd tendency coefficient controlling the tendency intensity added. +Default: 0.0 + + + +MCSP vwnd tendency coefficient controlling the tendency intensity added. +Default: 0.0 + + + shr_kind_r8 + use cam_abortutils, only: endrun + use cam_logfile, only: iulog +#endif + use physconst, only: cpair, pi + + + implicit none + private + + public :: zm_conv_mcsp_init ! Initialize MCSP output fields +!KH:++++++++++++++++ + public :: zm_conv_mcsp_readnl ! Read MCSP namelist +!KH:++++++++++++++++ + public :: zm_conv_mcsp_tend ! Perform MCSP tendency calculations + public :: zm_conv_mcsp_hist ! Write diagnostic quantities to history files + + real(r8), parameter :: MCSP_storm_speed_pref = 600e2_r8 ! pressure level for winds in MCSP calculation [Pa] + real(r8), parameter :: MCSP_conv_depth_min = 500e2_r8 ! pressure thickness of convective heating [Pa] + real(r8), parameter :: mcsp_shear_min = 3.0_r8 ! min shear value for MCSP to be active + real(r8), parameter :: mcsp_shear_max = 200.0_r8 ! max shear value for MCSP to be active +!KH:--------- +! real(r8) :: zmconv_MCSP_heat_coeff = 0.3_r8 +! real(r8) :: zmconv_MCSP_moisture_coeff = 0._r8 +! real(r8) :: zmconv_MCSP_uwind_coeff = 0._r8 +! real(r8) :: zmconv_MCSP_vwind_coeff = 0._r8 +!KH:--------- + +!KH:+++++++++ + real(r8) :: zmconv_MCSP_heat_coeff = 0.0_r8 + real(r8) :: zmconv_MCSP_moisture_coeff = 0.0_r8 + real(r8) :: zmconv_MCSP_uwind_coeff = 0.0_r8 + real(r8) :: zmconv_MCSP_vwind_coeff = 0.0_r8 +!KH:+++++++++ + + +!=================================================================================================== +contains +!=================================================================================================== + +subroutine zm_conv_mcsp_init() + !---------------------------------------------------------------------------- + ! Purpose: initialize MCSP output fields + !---------------------------------------------------------------------------- +#ifndef SCREAM_CONFIG_IS_CMAKE + use cam_history, only: addfld, horiz_only + use mpishorthand + !---------------------------------------------------------------------------- + !++ + call addfld('MCSP_DT_max', horiz_only, 'A', 'K/s', 'MCSP max T tendency') + !-- + call addfld('MCSP_DT', (/'lev'/), 'A', 'K/s', 'MCSP T tendency') + call addfld('MCSP_DQ', (/'lev'/), 'A', 'kg/kg/s', 'MCSP qv tendency') + call addfld('MCSP_DU', (/'lev'/), 'A', 'm/s/day', 'MCSP U wind tendency') + call addfld('MCSP_DV', (/'lev'/), 'A', 'm/s/day', 'MCSP V wind tendency') + + call addfld('MCSP_freq', horiz_only, 'A', '1', 'MCSP frequency of activation') + call addfld('MCSP_shear', horiz_only, 'A', 'm/s', 'MCSP vertical shear of zonal wind') + call addfld('MCSP_zm_depth', horiz_only, 'A', 'Pa', 'ZM convection depth for MCSP') +#endif /* SCREAM_CONFIG_IS_CMAKE */ + +end subroutine zm_conv_mcsp_init + +!KH:+++++++++++++++++++++++++++ +!========================================================================================= + +subroutine zm_conv_mcsp_readnl(nlfile) + + use spmd_utils, only: mpicom, masterproc, masterprocid, mpi_real8, mpi_integer, mpi_logical + use namelist_utils, only: find_group_name + + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! Local variables + integer :: unitn, ierr + character(len=*), parameter :: subname = 'zm_conv_mcsp_readnl' + + namelist /zmconv_mcsp_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & + zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff + !----------------------------------------------------------------------------- + + if (masterproc) then + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'zmconv_mcsp_nl', status=ierr) + if (ierr == 0) then + read(unitn, zmconv_mcsp_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist') + end if + end if + close(unitn) + + end if + + ! Broadcast namelist variables + call mpi_bcast(zmconv_MCSP_heat_coeff, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") + call mpi_bcast(zmconv_MCSP_moisture_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") + call mpi_bcast(zmconv_MCSP_uwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") + call mpi_bcast(zmconv_MCSP_vwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") + +end subroutine zm_conv_mcsp_readnl +!KH:+++++++++++++++++++++++++++++++ + +!=================================================================================================== + +subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear) + !---------------------------------------------------------------------------- + ! Purpose: calculate shear for MCSP + !---------------------------------------------------------------------------- +#ifdef SCREAM_CONFIG_IS_CMAKE + use zm_eamxx_bridge_methods, only: vertinterp +#else + use interpolate_data, only: vertinterp +#endif + !---------------------------------------------------------------------------- + ! Arguments + integer, intent(in ) :: pcols ! number of atmospheric columns (max) + integer, intent(in ) :: ncol ! number of atmospheric columns (actual) + integer, intent(in ) :: pver ! number of mid-point vertical levels + real(r8), dimension(pcols,pver), intent(in ) :: state_pmid ! physics state mid-point pressure + real(r8), dimension(pcols,pver), intent(in ) :: state_u ! physics state u momentum + real(r8), dimension(pcols,pver), intent(in ) :: state_v ! physics state v momentum + real(r8), dimension(pcols), intent( out) :: mcsp_shear + !---------------------------------------------------------------------------- + ! Local variables + integer :: i + real(r8), dimension(pcols) :: storm_u ! u wind at storm reference level set by MCSP_storm_speed_pref + real(r8), dimension(pcols) :: storm_v ! v wind at storm reference level set by MCSP_storm_speed_pref + real(r8), dimension(pcols) :: storm_u_shear ! u shear at storm reference level set by MCSP_storm_speed_pref + real(r8), dimension(pcols) :: storm_v_shear ! v shear at storm reference level set by MCSP_storm_speed_pref + !---------------------------------------------------------------------------- + ! Interpolate wind to pressure level specified by MCSP_storm_speed_pref + call vertinterp( ncol, pcols, pver, state_pmid, MCSP_storm_speed_pref, state_u, storm_u ) + call vertinterp( ncol, pcols, pver, state_pmid, MCSP_storm_speed_pref, state_v, storm_v ) + + !---------------------------------------------------------------------------- + ! calculate low-level shear + do i = 1,ncol + if (state_pmid(i,pver).gt.MCSP_storm_speed_pref) then + storm_u_shear(i) = storm_u(i)-state_u(i,pver) + storm_v_shear(i) = storm_v(i)-state_v(i,pver) + else + storm_u_shear(i) = -999 + storm_v_shear(i) = -999 + end if + mcsp_shear(i) = storm_u_shear(i) + end do + + !---------------------------------------------------------------------------- + return + +end subroutine zm_conv_mcsp_calculate_shear + +!=================================================================================================== + +subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & + ztodt, jctop, & + state_pmid, state_pint, state_pdel, & + state_s, state_q, state_u, state_v, & + ptend_zm_s, ptend_zm_q, & + ptend_s, ptend_q, ptend_u, ptend_v, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) + !---------------------------------------------------------------------------- + ! Purpose: perform MCSP tendency calculations + !---------------------------------------------------------------------------- + ! Arguments + integer, intent(in ) :: pcols ! number of atmospheric columns (max) + integer, intent(in ) :: ncol ! number of atmospheric columns (actual) + integer, intent(in ) :: pver ! number of mid-point vertical levels + integer, intent(in ) :: pverp ! number of interface vertical levels + real(r8), intent(in ) :: ztodt ! 2x physics time step + integer, dimension(pcols), intent(in ) :: jctop ! cloud top level indices + real(r8), dimension(pcols,pver), intent(in ) :: state_pmid ! physics state mid-point pressure + real(r8), dimension(pcols,pverp), intent(in ) :: state_pint ! physics state interface pressure + real(r8), dimension(pcols,pver), intent(in ) :: state_pdel ! physics state pressure thickness + real(r8), dimension(pcols,pver), intent(in ) :: state_s ! physics state dry static energy + real(r8), dimension(pcols,pver), intent(in ) :: state_q ! physics state specific humidity + real(r8), dimension(pcols,pver), intent(in ) :: state_u ! physics state u momentum + real(r8), dimension(pcols,pver), intent(in ) :: state_v ! physics state v momentum + real(r8), dimension(pcols,pver), intent(in ) :: ptend_zm_s ! input ZM tendency for dry static energy (DSE) + real(r8), dimension(pcols,pver), intent(in ) :: ptend_zm_q ! input ZM tendency for specific humidity (qv) + real(r8), dimension(pcols,pver), intent(inout) :: ptend_s ! output tendency of DSE + real(r8), dimension(pcols,pver), intent(inout) :: ptend_q ! output tendency of qv + real(r8), dimension(pcols,pver), intent(inout) :: ptend_u ! output tendency of u-wind + real(r8), dimension(pcols,pver), intent(inout) :: ptend_v ! output tendency of v-wind + real(r8), dimension(pcols,pver), intent( out) :: mcsp_dt_out! final MCSP tendency for DSE + real(r8), dimension(pcols,pver), intent( out) :: mcsp_dq_out! final MCSP tendency for qv + real(r8), dimension(pcols,pver), intent( out) :: mcsp_du_out! final MCSP tendency for u wind + real(r8), dimension(pcols,pver), intent( out) :: mcsp_dv_out! final MCSP tendency for v wind + real(r8), dimension(pcols), intent( out) :: mcsp_freq ! MSCP frequency for output + real(r8), dimension(pcols), intent( out) :: mcsp_shear ! shear used to check against threshold + real(r8), dimension(pcols), intent( out) :: zm_depth ! pressure depth of ZM heating + !++ + real(r8), dimension(pcols), intent( out) :: mcsp_tend_s_max ! max MCSP heating tendency + !-- + !---------------------------------------------------------------------------- + ! Local variables + integer :: i, k + + real(r8) :: tend_k ! temporary variable for kinetic energy tendency + real(r8) :: pdepth_mid_k ! temporary pressure depth used for vertical structure + real(r8) :: pdepth_total ! temporary pressure depth used for vertical structure + + real(r8), dimension(pcols) :: zm_avg_tend_s ! mass weighted column average DSE tendency from ZM + real(r8), dimension(pcols) :: zm_avg_tend_q ! mass weighted column average qv tendency from ZM + real(r8), dimension(pcols) :: pdel_sum ! column integrated pressure thickness + + real(r8), dimension(pcols,pver) :: mcsp_tend_s ! MCSP tendency before energy fixer for DSE + real(r8), dimension(pcols,pver) :: mcsp_tend_q ! MCSP tendency before energy fixer for qv + real(r8), dimension(pcols,pver) :: mcsp_tend_u ! MCSP tendency before energy fixer for u wind + real(r8), dimension(pcols,pver) :: mcsp_tend_v ! MCSP tendency before energy fixer for v wind + + real(r8), dimension(pcols) :: mcsp_avg_tend_s ! mass weighted column average MCSP tendency of DSE + real(r8), dimension(pcols) :: mcsp_avg_tend_q ! mass weighted column average MCSP tendency of qv + real(r8), dimension(pcols) :: mcsp_avg_tend_k ! mass weighted column average MCSP tendency of kinetic energy + + logical :: do_mcsp_t = .false. ! internal flag to enable tendency calculations + logical :: do_mcsp_q = .false. ! internal flag to enable tendency calculations + logical :: do_mcsp_u = .false. ! internal flag to enable tendency calculations + logical :: do_mcsp_v = .false. ! internal flag to enable tendency calculations + + !---------------------------------------------------------------------------- + ! initialize variables + + if (zmconv_MCSP_heat_coeff>0) do_mcsp_t = .true. + if (zmconv_MCSP_moisture_coeff>0) do_mcsp_q = .true. + if (zmconv_MCSP_uwind_coeff>0) do_mcsp_u = .true. + if (zmconv_MCSP_vwind_coeff>0) do_mcsp_v = .true. + + zm_avg_tend_s(1:ncol) = 0 + zm_avg_tend_q(1:ncol) = 0 + + pdel_sum(1:ncol) = 0 + + mcsp_avg_tend_s(1:ncol) = 0 + mcsp_avg_tend_q(1:ncol) = 0 + mcsp_avg_tend_k(1:ncol) = 0 + + mcsp_tend_s(1:ncol,1:pver) = 0 + mcsp_tend_q(1:ncol,1:pver) = 0 + mcsp_tend_u(1:ncol,1:pver) = 0 + mcsp_tend_v(1:ncol,1:pver) = 0 + !++ + mcsp_tend_s_max = 0._r8 + !-- + !---------------------------------------------------------------------------- + ! calculate shear + + call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear ) + + !---------------------------------------------------------------------------- + ! calculate mass weighted column average tendencies from ZM + + zm_depth(1:ncol) = 0 + do i = 1,ncol + if (jctop(i) .ne. pver) then + ! integrate pressure and ZM tendencies over column + do k = jctop(i),pver + zm_avg_tend_s(i) = zm_avg_tend_s(i) + ptend_zm_s(i,k) * state_pdel(i,k) + zm_avg_tend_q(i) = zm_avg_tend_q(i) + ptend_zm_q(i,k) * state_pdel(i,k) + pdel_sum(i) = pdel_sum(i) + state_pdel(i,k) + end do + ! normalize integrated ZM tendencies by total mass + zm_avg_tend_s(i) = zm_avg_tend_s(i) / pdel_sum(i) + zm_avg_tend_q(i) = zm_avg_tend_q(i) / pdel_sum(i) + ! calculate diagnostic zm_depth + zm_depth(i) = state_pint(i,pver+1) - state_pmid(i,jctop(i)) + else + zm_avg_tend_s(i) = 0 + zm_avg_tend_q(i) = 0 + zm_depth(i) = 0 + end if + end do + + !---------------------------------------------------------------------------- + ! Note: To conserve total energy we need to account for the kinteic energy tendency + ! which we can obtain from the velocity tendencies based on the following: + ! KE_new = (u_new^2 + v_new^2)/2 + ! = [ (u_old+du)^2 + (v_old+dv)^2 ]/2 + ! = [ ( u_old^2 + 2*u_old*du + du^2 ) + ( v_old^2 + 2*v_old*dv + dv^2 ) ]/2 + ! = ( u_old^2 + v_old^2 )/2 + ( 2*u_old*du + du^2 + 2*v_old*dv + dv^2 )/2 + ! = KE_old + [ 2*u_old*du + du^2 + 2*v_old*dv + dv^2 ] /2 + + !---------------------------------------------------------------------------- + ! calculate MCSP tendencies + + do i = 1,ncol + + ! check that ZM produced tendencies over a depth that exceeds the threshold + if ( zm_depth(i) >= MCSP_conv_depth_min ) then + ! check that ZM provided a non-zero column total heating + if ( zm_avg_tend_s(i) > 0 ) then + ! check that there is sufficient wind shear to justify coherent organization + if ( abs(mcsp_shear(i)).ge.mcsp_shear_min .and. & + abs(mcsp_shear(i)).lt.mcsp_shear_max ) then + do k = jctop(i),pver + + ! See eq 7-8 of Moncrieff et al. (2017) - also eq (5) of Moncrieff & Liu (2006) + pdepth_mid_k = state_pint(i,pver+1) - state_pmid(i,k) + pdepth_total = state_pint(i,pver+1) - state_pmid(i,jctop(i)) + + ! specify the assumed vertical structure + if (do_mcsp_t) mcsp_tend_s(i,k) = -1*zmconv_MCSP_heat_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) + if (do_mcsp_q) mcsp_tend_q(i,k) = -1*zmconv_MCSP_moisture_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) + if (do_mcsp_u) mcsp_tend_u(i,k) = zmconv_MCSP_uwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) + if (do_mcsp_v) mcsp_tend_v(i,k) = zmconv_MCSP_vwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) + + ! scale the vertical structure by the ZM heating/drying tendencies + if (do_mcsp_t) mcsp_tend_s(i,k) = zm_avg_tend_s(i) * mcsp_tend_s(i,k) + !++ + if (do_mcsp_t) mcsp_tend_s_max(i) = zm_avg_tend_s(i) * zmconv_MCSP_heat_coeff / cpair + !-- + if (do_mcsp_q) mcsp_tend_q(i,k) = zm_avg_tend_q(i) * mcsp_tend_q(i,k) + + ! integrate the DSE/qv tendencies for energy/mass fixer + if (do_mcsp_t) mcsp_avg_tend_s(i) = mcsp_avg_tend_s(i) + mcsp_tend_s(i,k) * state_pdel(i,k) / pdel_sum(i) + if (do_mcsp_q) mcsp_avg_tend_q(i) = mcsp_avg_tend_q(i) + mcsp_tend_q(i,k) * state_pdel(i,k) / pdel_sum(i) + + ! integrate the change in kinetic energy (KE) for energy fixer + if (do_mcsp_u.or.do_mcsp_v) then + tend_k = ( 2.0_r8*mcsp_tend_u(i,k)*ztodt*state_u(i,k) + mcsp_tend_u(i,k)*mcsp_tend_u(i,k)*ztodt*ztodt & + +2.0_r8*mcsp_tend_v(i,k)*ztodt*state_v(i,k) + mcsp_tend_v(i,k)*mcsp_tend_v(i,k)*ztodt*ztodt )/2.0_r8/ztodt + mcsp_avg_tend_k(i) = mcsp_avg_tend_k(i) + tend_k*state_pdel(i,k) / pdel_sum(i) + end if + + end do ! k = jctop(i),pver + end if ! shear threshold + end if ! zm_avg_tend_s(i) > 0 + end if ! zm_depth(i) >= MCSP_conv_depth_min + end do + + !---------------------------------------------------------------------------- + ! calculate final output tendencies + + mcsp_dt_out(1:ncol,1:pver) = 0 + mcsp_dq_out(1:ncol,1:pver) = 0 + mcsp_du_out(1:ncol,1:pver) = 0 + mcsp_dv_out(1:ncol,1:pver) = 0 + + mcsp_freq(1:ncol) = 0 + + do i = 1,ncol + do k = jctop(i),pver + + ! update frequency if MCSP contributes any tendency in the column + if ( abs(mcsp_tend_s(i,k)).gt.0 .or. abs(mcsp_tend_q(i,k)).gt.0 .or.& + abs(mcsp_tend_u(i,k)).gt.0 .or. abs(mcsp_tend_v(i,k)).gt.0 ) then + mcsp_freq(i) = 1 + end if + + ! subtract mass weighted average tendencies for energy/mass conservation + mcsp_dt_out(i,k) = mcsp_tend_s(i,k) - mcsp_avg_tend_s(i) + mcsp_dq_out(i,k) = mcsp_tend_q(i,k) - mcsp_avg_tend_q(i) + mcsp_du_out(i,k) = mcsp_tend_u(i,k) + mcsp_dv_out(i,k) = mcsp_tend_v(i,k) + + ! make sure kinetic energy correction is added to DSE tendency + ! to conserve total energy whenever momentum tendencies are calculated + if (do_mcsp_u.or.do_mcsp_v) then + mcsp_dt_out(i,k) = mcsp_dt_out(i,k) - mcsp_avg_tend_k(i) + end if + + ! update output tendencies + if (do_mcsp_t) ptend_s(i,k) = ptend_s(i,k) + mcsp_dt_out(i,k) + if (do_mcsp_q) ptend_q(i,k) = ptend_q(i,k) + mcsp_dq_out(i,k) + if (do_mcsp_u) ptend_u(i,k) = ptend_u(i,k) + mcsp_du_out(i,k) + if (do_mcsp_v) ptend_v(i,k) = ptend_v(i,k) + mcsp_dv_out(i,k) + + ! adjust units for diagnostic outputs + if (do_mcsp_t) mcsp_dt_out(i,k) = mcsp_dt_out(i,k)/cpair + + end do + end do + + !---------------------------------------------------------------------------- + return + +end subroutine zm_conv_mcsp_tend + +!=================================================================================================== + +subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) + !---------------------------------------------------------------------------- + ! Purpose: write diagnostic quantities to history files + !---------------------------------------------------------------------------- +#ifndef SCREAM_CONFIG_IS_CMAKE + use cam_history, only: outfld +#endif + !---------------------------------------------------------------------------- + ! Arguments + integer, intent(in) :: lchnk ! chunk identifier + integer, intent(in) :: pcols ! number of atmospheric columns (max) + integer, intent(in) :: pver ! number of mid-point vertical levels + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dt_out ! final MCSP tendency for DSE + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dq_out ! final MCSP tendency for qv + real(r8), dimension(pcols,pver), intent(in) :: mcsp_du_out ! final MCSP tendency for u wind + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dv_out ! final MCSP tendency for v wind + real(r8), dimension(pcols), intent(in) :: mcsp_freq ! MSCP frequency for output + real(r8), dimension(pcols), intent(in) :: mcsp_shear ! shear used to check against threshold + real(r8), dimension(pcols), intent(in) :: zm_depth ! pressure depth of ZM heating + !++ + real(r8), dimension(pcols), intent(in) :: mcsp_tend_s_max ! max MCSP heating tendency + !-- + !---------------------------------------------------------------------------- +#ifndef SCREAM_CONFIG_IS_CMAKE + ! write out MCSP diagnostic history fields + call outfld('MCSP_DT', mcsp_dt_out, pcols, lchnk ) + call outfld('MCSP_DQ', mcsp_dq_out, pcols, lchnk ) + call outfld('MCSP_DU', mcsp_du_out, pcols, lchnk ) + call outfld('MCSP_DV', mcsp_dv_out, pcols, lchnk ) + call outfld('MCSP_freq', mcsp_freq, pcols, lchnk ) + call outfld('MCSP_shear', mcsp_shear, pcols, lchnk ) + call outfld('MCSP_zm_depth', zm_depth, pcols, lchnk ) + !++ + call outfld('MCSP_DT_max', mcsp_tend_s_max, pcols, lchnk) + !-- +#endif + !---------------------------------------------------------------------------- + return + +end subroutine zm_conv_mcsp_hist + +!=================================================================================================== + +end module zm_conv_mcsp diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 6a4f9d8cfa..5f44907e28 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -78,7 +78,6 @@ module zm_conv_intr real(r8) :: zmconv_parcel_hscale = unset_r8 ! Fraction of PBL depth over which to mix initial parcel real(r8) :: zmconv_tau = unset_r8 ! Timescale for convection - ! indices for fields in the physics buffer integer :: cld_idx = 0 integer :: icwmrdp_idx = 0 @@ -225,6 +224,9 @@ subroutine zm_conv_init(pref_edge) use spmd_utils, only: masterproc use phys_control, only: phys_deepconv_pbl, phys_getopts, cam_physpkg_is use physics_buffer, only: pbuf_get_index +!++ MCSP + use zm_conv_mcsp, only: zm_conv_mcsp_init +!-- MCSP implicit none @@ -308,6 +310,10 @@ subroutine zm_conv_init(pref_edge) call add_default('ZMMTT ', history_budget_histfile_num, ' ') end if +!++ MCSP + call zm_conv_mcsp_init() +!-- MCSP + ! ! Limit deep convection to regions below 40 mb ! Note this calculation is repeated in the shallow convection interface @@ -383,6 +389,9 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & use phys_control, only: cam_physpkg_is use ccpp_constituent_prop_mod, only: ccpp_const_props +!++ MCSP + use zm_conv_mcsp, only: zm_conv_mcsp_tend, zm_conv_mcsp_hist +!-- MCSP ! Arguments @@ -427,6 +436,28 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & type(physics_state) :: state1 ! locally modify for evaporation to use, not returned type(physics_ptend),target :: ptend_loc ! package tendencies +!++ MCSP + type(physics_ptend),target :: ptend_mcsp ! MCSP output tendencies + + ! flags for MCSP tendencies + logical :: do_mcsp_t = .false. + logical :: do_mcsp_q(pcnst) = .false. + logical :: do_mcsp_u = .false. + logical :: do_mcsp_v = .false. + + ! MCSP history output variables + real(r8), dimension(pcols,pver) :: mcsp_dt_out ! MCSP tendency for DSE + real(r8), dimension(pcols,pver) :: mcsp_dq_out ! MCSP tendency for qv + real(r8), dimension(pcols,pver) :: mcsp_du_out ! MCSP tendency for u wind + real(r8), dimension(pcols,pver) :: mcsp_dv_out ! MCSP tendency for v wind + real(r8), dimension(pcols) :: mcsp_freq ! MSCP frequency for output + real(r8), dimension(pcols) :: mcsp_shear ! shear used to check against threshold + real(r8), dimension(pcols) :: zm_depth ! pressure depth of ZM heating + !++ MCSP + real(r8), dimension(pcols) :: mcsp_tend_s_max ! max MCSP heating tendency + !-- MCSP +!-- MCSP + ! physics buffer fields real(r8), pointer, dimension(:) :: prec ! total precipitation real(r8), pointer, dimension(:) :: snow ! snow from ZM convection @@ -593,6 +624,45 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & end do call outfld('CAPE', cape, pcols, lchnk) ! RBN - CAPE output + +!++ MCSP + !---------------------------------------------------------------------------- + ! mesoscale coherent structure parameterization (MCSP) + ! Note that this modifies the tendencies produced by zm_convr(), such that + ! history variables like ZMDT will include the effects of MCSP + + ! Set these all to True so that the tedency variables get allocated. The internal flags to + ! calculate each MCSP tedency will behave the same, but having them always allocated avoids + ! a problem with bridging to this routine from C++ for porting ZM to EAMxx + do_mcsp_t = .true. + do_mcsp_q(1) = .true. + do_mcsp_u = .true. + do_mcsp_v = .true. + + call physics_ptend_init( ptend_mcsp, state%psetcols, 'zm_conv_mcsp_tend', & + ls=do_mcsp_t, lq=do_mcsp_q, lu=do_mcsp_u, lv=do_mcsp_v) + + call zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & + ztodt, int(jctop), & + state%pmid, state%pint, state%pdel, & + state%s, state%q, state%u, state%v, & + ptend_loc%s, ptend_loc%q(:,:,1), & + ptend_mcsp%s(:,:), ptend_mcsp%q(:,:,1), & + ptend_mcsp%u(:,:), ptend_mcsp%v(:,:), & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) + + call zm_conv_mcsp_hist( lchnk, pcols, pver, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) + + ! add MCSP tendencies to ZM convective tendencies + call physics_ptend_sum( ptend_mcsp, ptend_loc, ncol) + call physics_ptend_dealloc(ptend_mcsp) + +!-- MCSP + + ! ! Output fractional occurance of ZM convection ! From d04e02cce6ea98011e10a3a534e1ada46fe4f150 Mon Sep 17 00:00:00 2001 From: Kai Huang Date: Thu, 26 Mar 2026 15:27:10 -0600 Subject: [PATCH 2/4] KH comments and SCREAM lines removed. If statement added to zm_conv_mcsp_tend so that no calculations will be done if mcsp coefficients are all zero. Test runs are underway. --- bld/namelist_files/namelist_defaults_cam.xml | 2 +- bld/namelist_files/namelist_definition.xml | 10 +- src/control/runtime_opts.F90 | 4 - src/physics/cam/zm_conv_MCSP.F90 | 222 +++++++++---------- src/physics/cam/zm_conv_intr.F90 | 2 +- 5 files changed, 107 insertions(+), 133 deletions(-) diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index e86e93a0dc..0dc50f2f6b 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -2891,7 +2891,7 @@ .false. 0.5 - + 0.0D0 0.0D0 0.0D0 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 7e041b8a74..befebb3772 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3676,28 +3676,28 @@ the zonally averaged zonal wind in the sponge layer. Default: 0.0 - + + group="zmconv_nl" valid_values="" > MCSP heating coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_nl" valid_values="" > MCSP moistening coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_nl" valid_values="" > MCSP uwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_nl" valid_values="" > MCSP vwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 5f425e1380..6b56bf7df8 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -50,9 +50,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use cldfrc2m, only: cldfrc2m_readnl use rk_stratiform_cam, only: rk_stratiform_cam_readnl use zm_conv_intr, only: zm_conv_readnl -!KH:++++++++++++++++ use zm_conv_mcsp, only: zm_conv_mcsp_readnl -!KH:++++++++++++++++ use hk_conv, only: hkconv_readnl use uwshcu, only: uwshcu_readnl use pkg_cld_sediment, only: cld_sediment_readnl @@ -156,9 +154,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call cldfrc_readnl(nlfilename) call cldfrc2m_readnl(nlfilename) call zm_conv_readnl(nlfilename) -!KH:+++++++++++++++++ call zm_conv_mcsp_readnl(nlfilename) -!KH:+++++++++++++++++ call rk_stratiform_cam_readnl(nlfilename) call hkconv_readnl(nlfilename) call uwshcu_readnl(nlfilename) diff --git a/src/physics/cam/zm_conv_MCSP.F90 b/src/physics/cam/zm_conv_MCSP.F90 index b0e0b5450f..1658993d65 100644 --- a/src/physics/cam/zm_conv_MCSP.F90 +++ b/src/physics/cam/zm_conv_MCSP.F90 @@ -27,13 +27,10 @@ module zm_conv_mcsp ! https://doi.org/10.1029/2019GL085316 ! !---------------------------------------------------------------------------- -#ifdef SCREAM_CONFIG_IS_CMAKE - use zm_eamxx_bridge_params, only: r8 -#else + use shr_kind_mod, only: r8=>shr_kind_r8 use cam_abortutils, only: endrun use cam_logfile, only: iulog -#endif use physconst, only: cpair, pi @@ -41,9 +38,7 @@ module zm_conv_mcsp private public :: zm_conv_mcsp_init ! Initialize MCSP output fields -!KH:++++++++++++++++ public :: zm_conv_mcsp_readnl ! Read MCSP namelist -!KH:++++++++++++++++ public :: zm_conv_mcsp_tend ! Perform MCSP tendency calculations public :: zm_conv_mcsp_hist ! Write diagnostic quantities to history files @@ -51,20 +46,11 @@ module zm_conv_mcsp real(r8), parameter :: MCSP_conv_depth_min = 500e2_r8 ! pressure thickness of convective heating [Pa] real(r8), parameter :: mcsp_shear_min = 3.0_r8 ! min shear value for MCSP to be active real(r8), parameter :: mcsp_shear_max = 200.0_r8 ! max shear value for MCSP to be active -!KH:--------- -! real(r8) :: zmconv_MCSP_heat_coeff = 0.3_r8 -! real(r8) :: zmconv_MCSP_moisture_coeff = 0._r8 -! real(r8) :: zmconv_MCSP_uwind_coeff = 0._r8 -! real(r8) :: zmconv_MCSP_vwind_coeff = 0._r8 -!KH:--------- - -!KH:+++++++++ - real(r8) :: zmconv_MCSP_heat_coeff = 0.0_r8 - real(r8) :: zmconv_MCSP_moisture_coeff = 0.0_r8 - real(r8) :: zmconv_MCSP_uwind_coeff = 0.0_r8 - real(r8) :: zmconv_MCSP_vwind_coeff = 0.0_r8 -!KH:+++++++++ + real(r8) :: zmconv_MCSP_heat_coeff + real(r8) :: zmconv_MCSP_moisture_coeff + real(r8) :: zmconv_MCSP_uwind_coeff + real(r8) :: zmconv_MCSP_vwind_coeff !=================================================================================================== contains @@ -74,7 +60,6 @@ subroutine zm_conv_mcsp_init() !---------------------------------------------------------------------------- ! Purpose: initialize MCSP output fields !---------------------------------------------------------------------------- -#ifndef SCREAM_CONFIG_IS_CMAKE use cam_history, only: addfld, horiz_only use mpishorthand !---------------------------------------------------------------------------- @@ -89,11 +74,10 @@ subroutine zm_conv_mcsp_init() call addfld('MCSP_freq', horiz_only, 'A', '1', 'MCSP frequency of activation') call addfld('MCSP_shear', horiz_only, 'A', 'm/s', 'MCSP vertical shear of zonal wind') call addfld('MCSP_zm_depth', horiz_only, 'A', 'Pa', 'ZM convection depth for MCSP') -#endif /* SCREAM_CONFIG_IS_CMAKE */ end subroutine zm_conv_mcsp_init -!KH:+++++++++++++++++++++++++++ + !========================================================================================= subroutine zm_conv_mcsp_readnl(nlfile) @@ -107,15 +91,15 @@ subroutine zm_conv_mcsp_readnl(nlfile) integer :: unitn, ierr character(len=*), parameter :: subname = 'zm_conv_mcsp_readnl' - namelist /zmconv_mcsp_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & - zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff + namelist /zmconv_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & + zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff !----------------------------------------------------------------------------- if (masterproc) then open( newunit=unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'zmconv_mcsp_nl', status=ierr) + call find_group_name(unitn, 'zmconv_nl', status=ierr) if (ierr == 0) then - read(unitn, zmconv_mcsp_nl, iostat=ierr) + read(unitn, zmconv_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if @@ -126,16 +110,15 @@ subroutine zm_conv_mcsp_readnl(nlfile) ! Broadcast namelist variables call mpi_bcast(zmconv_MCSP_heat_coeff, 1, mpi_integer, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") + if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") call mpi_bcast(zmconv_MCSP_moisture_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") + if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") call mpi_bcast(zmconv_MCSP_uwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") + if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") call mpi_bcast(zmconv_MCSP_vwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") + if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") end subroutine zm_conv_mcsp_readnl -!KH:+++++++++++++++++++++++++++++++ !=================================================================================================== @@ -143,11 +126,7 @@ subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, !---------------------------------------------------------------------------- ! Purpose: calculate shear for MCSP !---------------------------------------------------------------------------- -#ifdef SCREAM_CONFIG_IS_CMAKE - use zm_eamxx_bridge_methods, only: vertinterp -#else use interpolate_data, only: vertinterp -#endif !---------------------------------------------------------------------------- ! Arguments integer, intent(in ) :: pcols ! number of atmospheric columns (max) @@ -280,91 +259,94 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & !++ mcsp_tend_s_max = 0._r8 !-- - !---------------------------------------------------------------------------- - ! calculate shear - - call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear ) - - !---------------------------------------------------------------------------- - ! calculate mass weighted column average tendencies from ZM - - zm_depth(1:ncol) = 0 - do i = 1,ncol - if (jctop(i) .ne. pver) then - ! integrate pressure and ZM tendencies over column - do k = jctop(i),pver - zm_avg_tend_s(i) = zm_avg_tend_s(i) + ptend_zm_s(i,k) * state_pdel(i,k) - zm_avg_tend_q(i) = zm_avg_tend_q(i) + ptend_zm_q(i,k) * state_pdel(i,k) - pdel_sum(i) = pdel_sum(i) + state_pdel(i,k) - end do - ! normalize integrated ZM tendencies by total mass - zm_avg_tend_s(i) = zm_avg_tend_s(i) / pdel_sum(i) - zm_avg_tend_q(i) = zm_avg_tend_q(i) / pdel_sum(i) - ! calculate diagnostic zm_depth - zm_depth(i) = state_pint(i,pver+1) - state_pmid(i,jctop(i)) - else - zm_avg_tend_s(i) = 0 - zm_avg_tend_q(i) = 0 - zm_depth(i) = 0 - end if - end do - - !---------------------------------------------------------------------------- - ! Note: To conserve total energy we need to account for the kinteic energy tendency - ! which we can obtain from the velocity tendencies based on the following: - ! KE_new = (u_new^2 + v_new^2)/2 - ! = [ (u_old+du)^2 + (v_old+dv)^2 ]/2 - ! = [ ( u_old^2 + 2*u_old*du + du^2 ) + ( v_old^2 + 2*v_old*dv + dv^2 ) ]/2 - ! = ( u_old^2 + v_old^2 )/2 + ( 2*u_old*du + du^2 + 2*v_old*dv + dv^2 )/2 - ! = KE_old + [ 2*u_old*du + du^2 + 2*v_old*dv + dv^2 ] /2 - - !---------------------------------------------------------------------------- - ! calculate MCSP tendencies + if (zmconv_MCSP_heat_coeff>0 .OR. zmconv_MCSP_moisture_coeff>0 .OR. zmconv_MCSP_uwind_coeff>0 .OR. zmconv_MCSP_vwind_coeff>0) then + !---------------------------------------------------------------------------- + ! calculate shear + + call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear ) + + !---------------------------------------------------------------------------- + ! calculate mass weighted column average tendencies from ZM + + zm_depth(1:ncol) = 0 + do i = 1,ncol + if (jctop(i) .ne. pver) then + ! integrate pressure and ZM tendencies over column + do k = jctop(i),pver + zm_avg_tend_s(i) = zm_avg_tend_s(i) + ptend_zm_s(i,k) * state_pdel(i,k) + zm_avg_tend_q(i) = zm_avg_tend_q(i) + ptend_zm_q(i,k) * state_pdel(i,k) + pdel_sum(i) = pdel_sum(i) + state_pdel(i,k) + end do + ! normalize integrated ZM tendencies by total mass + zm_avg_tend_s(i) = zm_avg_tend_s(i) / pdel_sum(i) + zm_avg_tend_q(i) = zm_avg_tend_q(i) / pdel_sum(i) + ! calculate diagnostic zm_depth + zm_depth(i) = state_pint(i,pver+1) - state_pmid(i,jctop(i)) + else + zm_avg_tend_s(i) = 0 + zm_avg_tend_q(i) = 0 + zm_depth(i) = 0 + end if + end do - do i = 1,ncol + !---------------------------------------------------------------------------- + ! Note: To conserve total energy we need to account for the kinteic energy tendency + ! which we can obtain from the velocity tendencies based on the following: + ! KE_new = (u_new^2 + v_new^2)/2 + ! = [ (u_old+du)^2 + (v_old+dv)^2 ]/2 + ! = [ ( u_old^2 + 2*u_old*du + du^2 ) + ( v_old^2 + 2*v_old*dv + dv^2 ) ]/2 + ! = ( u_old^2 + v_old^2 )/2 + ( 2*u_old*du + du^2 + 2*v_old*dv + dv^2 )/2 + ! = KE_old + [ 2*u_old*du + du^2 + 2*v_old*dv + dv^2 ] /2 + + !---------------------------------------------------------------------------- + ! calculate MCSP tendencies + + do i = 1,ncol + + ! check that ZM produced tendencies over a depth that exceeds the threshold + if ( zm_depth(i) >= MCSP_conv_depth_min ) then + ! check that ZM provided a non-zero column total heating + if ( zm_avg_tend_s(i) > 0 ) then + ! check that there is sufficient wind shear to justify coherent organization + if ( abs(mcsp_shear(i)).ge.mcsp_shear_min .and. & + abs(mcsp_shear(i)).lt.mcsp_shear_max ) then + do k = jctop(i),pver + + ! See eq 7-8 of Moncrieff et al. (2017) - also eq (5) of Moncrieff & Liu (2006) + pdepth_mid_k = state_pint(i,pver+1) - state_pmid(i,k) + pdepth_total = state_pint(i,pver+1) - state_pmid(i,jctop(i)) + + ! specify the assumed vertical structure + if (do_mcsp_t) mcsp_tend_s(i,k) = -1*zmconv_MCSP_heat_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) + if (do_mcsp_q) mcsp_tend_q(i,k) = -1*zmconv_MCSP_moisture_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) + if (do_mcsp_u) mcsp_tend_u(i,k) = zmconv_MCSP_uwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) + if (do_mcsp_v) mcsp_tend_v(i,k) = zmconv_MCSP_vwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) + + ! scale the vertical structure by the ZM heating/drying tendencies + if (do_mcsp_t) mcsp_tend_s(i,k) = zm_avg_tend_s(i) * mcsp_tend_s(i,k) + !++ + if (do_mcsp_t) mcsp_tend_s_max(i) = zm_avg_tend_s(i) * zmconv_MCSP_heat_coeff / cpair + !-- + if (do_mcsp_q) mcsp_tend_q(i,k) = zm_avg_tend_q(i) * mcsp_tend_q(i,k) + + ! integrate the DSE/qv tendencies for energy/mass fixer + if (do_mcsp_t) mcsp_avg_tend_s(i) = mcsp_avg_tend_s(i) + mcsp_tend_s(i,k) * state_pdel(i,k) / pdel_sum(i) + if (do_mcsp_q) mcsp_avg_tend_q(i) = mcsp_avg_tend_q(i) + mcsp_tend_q(i,k) * state_pdel(i,k) / pdel_sum(i) + + ! integrate the change in kinetic energy (KE) for energy fixer + if (do_mcsp_u.or.do_mcsp_v) then + tend_k = ( 2.0_r8*mcsp_tend_u(i,k)*ztodt*state_u(i,k) + mcsp_tend_u(i,k)*mcsp_tend_u(i,k)*ztodt*ztodt & + +2.0_r8*mcsp_tend_v(i,k)*ztodt*state_v(i,k) + mcsp_tend_v(i,k)*mcsp_tend_v(i,k)*ztodt*ztodt )/2.0_r8/ztodt + mcsp_avg_tend_k(i) = mcsp_avg_tend_k(i) + tend_k*state_pdel(i,k) / pdel_sum(i) + end if + + end do ! k = jctop(i),pver + end if ! shear threshold + end if ! zm_avg_tend_s(i) > 0 + end if ! zm_depth(i) >= MCSP_conv_depth_min + end do - ! check that ZM produced tendencies over a depth that exceeds the threshold - if ( zm_depth(i) >= MCSP_conv_depth_min ) then - ! check that ZM provided a non-zero column total heating - if ( zm_avg_tend_s(i) > 0 ) then - ! check that there is sufficient wind shear to justify coherent organization - if ( abs(mcsp_shear(i)).ge.mcsp_shear_min .and. & - abs(mcsp_shear(i)).lt.mcsp_shear_max ) then - do k = jctop(i),pver - - ! See eq 7-8 of Moncrieff et al. (2017) - also eq (5) of Moncrieff & Liu (2006) - pdepth_mid_k = state_pint(i,pver+1) - state_pmid(i,k) - pdepth_total = state_pint(i,pver+1) - state_pmid(i,jctop(i)) - - ! specify the assumed vertical structure - if (do_mcsp_t) mcsp_tend_s(i,k) = -1*zmconv_MCSP_heat_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) - if (do_mcsp_q) mcsp_tend_q(i,k) = -1*zmconv_MCSP_moisture_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) - if (do_mcsp_u) mcsp_tend_u(i,k) = zmconv_MCSP_uwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) - if (do_mcsp_v) mcsp_tend_v(i,k) = zmconv_MCSP_vwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) - - ! scale the vertical structure by the ZM heating/drying tendencies - if (do_mcsp_t) mcsp_tend_s(i,k) = zm_avg_tend_s(i) * mcsp_tend_s(i,k) - !++ - if (do_mcsp_t) mcsp_tend_s_max(i) = zm_avg_tend_s(i) * zmconv_MCSP_heat_coeff / cpair - !-- - if (do_mcsp_q) mcsp_tend_q(i,k) = zm_avg_tend_q(i) * mcsp_tend_q(i,k) - - ! integrate the DSE/qv tendencies for energy/mass fixer - if (do_mcsp_t) mcsp_avg_tend_s(i) = mcsp_avg_tend_s(i) + mcsp_tend_s(i,k) * state_pdel(i,k) / pdel_sum(i) - if (do_mcsp_q) mcsp_avg_tend_q(i) = mcsp_avg_tend_q(i) + mcsp_tend_q(i,k) * state_pdel(i,k) / pdel_sum(i) - - ! integrate the change in kinetic energy (KE) for energy fixer - if (do_mcsp_u.or.do_mcsp_v) then - tend_k = ( 2.0_r8*mcsp_tend_u(i,k)*ztodt*state_u(i,k) + mcsp_tend_u(i,k)*mcsp_tend_u(i,k)*ztodt*ztodt & - +2.0_r8*mcsp_tend_v(i,k)*ztodt*state_v(i,k) + mcsp_tend_v(i,k)*mcsp_tend_v(i,k)*ztodt*ztodt )/2.0_r8/ztodt - mcsp_avg_tend_k(i) = mcsp_avg_tend_k(i) + tend_k*state_pdel(i,k) / pdel_sum(i) - end if - - end do ! k = jctop(i),pver - end if ! shear threshold - end if ! zm_avg_tend_s(i) > 0 - end if ! zm_depth(i) >= MCSP_conv_depth_min - end do + endif !---------------------------------------------------------------------------- ! calculate final output tendencies @@ -422,9 +404,7 @@ subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & !---------------------------------------------------------------------------- ! Purpose: write diagnostic quantities to history files !---------------------------------------------------------------------------- -#ifndef SCREAM_CONFIG_IS_CMAKE use cam_history, only: outfld -#endif !---------------------------------------------------------------------------- ! Arguments integer, intent(in) :: lchnk ! chunk identifier @@ -441,7 +421,6 @@ subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & real(r8), dimension(pcols), intent(in) :: mcsp_tend_s_max ! max MCSP heating tendency !-- !---------------------------------------------------------------------------- -#ifndef SCREAM_CONFIG_IS_CMAKE ! write out MCSP diagnostic history fields call outfld('MCSP_DT', mcsp_dt_out, pcols, lchnk ) call outfld('MCSP_DQ', mcsp_dq_out, pcols, lchnk ) @@ -453,7 +432,6 @@ subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & !++ call outfld('MCSP_DT_max', mcsp_tend_s_max, pcols, lchnk) !-- -#endif !---------------------------------------------------------------------------- return diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 5f44907e28..18e1b48890 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -633,7 +633,7 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & ! Set these all to True so that the tedency variables get allocated. The internal flags to ! calculate each MCSP tedency will behave the same, but having them always allocated avoids - ! a problem with bridging to this routine from C++ for porting ZM to EAMxx + do_mcsp_t = .true. do_mcsp_q(1) = .true. do_mcsp_u = .true. From a806aba794399409906d9d1562d1c1d253578407 Mon Sep 17 00:00:00 2001 From: Kai Huang Date: Tue, 31 Mar 2026 15:38:56 -0600 Subject: [PATCH 3/4] version passing the test runs with updates --- bld/namelist_files/namelist_definition.xml | 8 ++++---- src/physics/cam/zm_conv_MCSP.F90 | 18 ++++++++++-------- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index befebb3772..6bcc68520c 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3679,25 +3679,25 @@ Default: 0.0 + group="zmconv_mcsp_nl" valid_values="" > MCSP heating coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_mcsp_nl" valid_values="" > MCSP moistening coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_mcsp_nl" valid_values="" > MCSP uwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 + group="zmconv_mcsp_nl" valid_values="" > MCSP vwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 diff --git a/src/physics/cam/zm_conv_MCSP.F90 b/src/physics/cam/zm_conv_MCSP.F90 index 1658993d65..803d91eeb3 100644 --- a/src/physics/cam/zm_conv_MCSP.F90 +++ b/src/physics/cam/zm_conv_MCSP.F90 @@ -91,15 +91,15 @@ subroutine zm_conv_mcsp_readnl(nlfile) integer :: unitn, ierr character(len=*), parameter :: subname = 'zm_conv_mcsp_readnl' - namelist /zmconv_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & + namelist /zmconv_mcsp_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff !----------------------------------------------------------------------------- if (masterproc) then open( newunit=unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'zmconv_nl', status=ierr) + call find_group_name(unitn, 'zmconv_mcsp_nl', status=ierr) if (ierr == 0) then - read(unitn, zmconv_nl, iostat=ierr) + read(unitn, zmconv_mcsp_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if @@ -110,13 +110,13 @@ subroutine zm_conv_mcsp_readnl(nlfile) ! Broadcast namelist variables call mpi_bcast(zmconv_MCSP_heat_coeff, 1, mpi_integer, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") call mpi_bcast(zmconv_MCSP_moisture_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") call mpi_bcast(zmconv_MCSP_uwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") call mpi_bcast(zmconv_MCSP_vwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") + if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") end subroutine zm_conv_mcsp_readnl @@ -259,7 +259,9 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & !++ mcsp_tend_s_max = 0._r8 !-- + if (zmconv_MCSP_heat_coeff>0 .OR. zmconv_MCSP_moisture_coeff>0 .OR. zmconv_MCSP_uwind_coeff>0 .OR. zmconv_MCSP_vwind_coeff>0) then + !---------------------------------------------------------------------------- ! calculate shear @@ -346,7 +348,7 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & end if ! zm_depth(i) >= MCSP_conv_depth_min end do - endif + end if !---------------------------------------------------------------------------- ! calculate final output tendencies From 5ac2fe335438ed6236d8c22b5e44845d46df1716 Mon Sep 17 00:00:00 2001 From: Kai Huang Date: Fri, 8 May 2026 09:53:37 -0600 Subject: [PATCH 4/4] updated version with CCPP-compliance --- bld/build-namelist | 65 ++- bld/namelist_files/namelist_defaults_cam.xml | 11 +- bld/namelist_files/namelist_definition.xml | 36 +- src/control/runtime_opts.F90 | 8 +- src/physics/cam/convect_deep.F90 | 14 +- src/physics/cam/mcsp_intr.F90 | 237 ++++++++++ src/physics/cam/physpkg.F90 | 80 +++- src/physics/cam/zm_conv_MCSP.F90 | 444 ------------------- src/physics/cam/zm_conv_intr.F90 | 80 +--- src/physics/cam7/physpkg.F90 | 72 ++- 10 files changed, 463 insertions(+), 584 deletions(-) create mode 100644 src/physics/cam/mcsp_intr.F90 delete mode 100644 src/physics/cam/zm_conv_MCSP.F90 diff --git a/bld/build-namelist b/bld/build-namelist index 29777bc449..1aeb2e0f4a 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -583,16 +583,13 @@ if ( ($chem ne 'none') or ( $prog_species ) ){ my $chem_proc_src = $cfg->get('chem_proc_src'); my $chem_src_dir = $cfg->get('chem_src_dir'); - my ( $gas_wetdep_list, $gas_wetdep_ice_uptake_list, $aer_wetdep_list, - $aer_sol_facti, $aer_sol_factb, $aer_scav_coef, $gas_drydep_list, $aer_drydep_list ) = + my ( $gas_wetdep_list, $aer_wetdep_list, $aer_sol_facti, $aer_sol_factb, $aer_scav_coef, + $aer_drydep_list, $gas_drydep_list ) = set_dep_lists( $chem, $cfgdir, $chem_proc_src, $chem_src_dir, $nl, $print ); if (length($gas_wetdep_list)>2){ add_default($nl, 'gas_wetdep_method' ); add_default($nl, 'gas_wetdep_list', 'val'=>$gas_wetdep_list ); - if (($chem =~ /_slh/) and (length($gas_wetdep_ice_uptake_list)>2)) { - add_default($nl, 'gas_wetdep_ice_uptake_list', 'val'=>$gas_wetdep_ice_uptake_list ); - } } if (length($aer_wetdep_list)>2){ @@ -1940,11 +1937,11 @@ if (($chem =~ /trop_mozart/ or $chem =~ /trop_strat/ or $chem =~ /waccm_tsmlt/) 'ISOP -> ' => 'isop_emis_file', 'TOLUENE -> ' => 'toluene_emis_file', ); - if (!($chem =~ /_vbs/) and !($chem =~ /_tsmlt/) and !($chem =~ /_t1s/)) { + if (!($chem =~ /_vbs/) and !($chem =~ /_tsmlt/)) { %species = (%species, 'NH3 -> ' => 'nh3_emis_file'); } - if (!($chem =~ /_vbs/) and !($chem =~ /_tsmlt/) and !($chem =~ /_t1s/)) { + if (!($chem =~ /_vbs/) and !($chem =~ /_tsmlt/)) { %species = (%species, 'C10H16 -> ' => 'c10h16_emis_file', ); } @@ -2360,25 +2357,6 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { 'soag_bb_srf_file' => 'SOAG' ); } - # for short-lived halogen (SLH) chemistry - if ($chem =~ /_slh/) { - %species = (%species, - 'c2cl4_slh_emis_file' => 'C2CL4', - 'ch2br2_slh_emis_file' => 'CH2BR2', - 'ch2brcl_slh_emis_file' => 'CH2BRCL', - 'ch2cl2_slh_emis_file' => 'CH2CL2', - 'ch2i2_slh_emis_file' => 'CH2I2', - 'ch2ibr_slh_emis_file' => 'CH2IBR', - 'ch2icl_slh_emis_file' => 'CH2ICL', - 'ch3i_slh_emis_file' => 'CH3I', - 'chbr2cl_slh_emis_file' => 'CHBR2CL', - 'chbr3_slh_emis_file' => 'CHBR3', - 'chbrcl2_slh_emis_file' => 'CHBRCL2', - 'hoi_slh_emis_file' => 'HOI', - 'i2_slh_emis_file' => 'I2' ); - } - - # for mid-atmos gas-phase chemistry if ($chem =~ /trop_strat/ or $chem =~ /_tsmlt/ or $chem =~ /waccm_ma/ or $chem =~ /waccm_t4ma/) { %species = (%species, @@ -2392,7 +2370,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { 'CH2O_bb_srf_file' => 'CH2O' ); } - # for troposphere gas-phase chemistry + # for troposphere gas-phase chemistry if ($chem =~ /trop_strat/ or $chem =~ /_tsmlt/ or $chem =~ /_t4ma/) { %species = (%species, 'BIGALK_an_srf_file' => 'BIGALK', @@ -2429,7 +2407,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { %species = (%species, 'E90_srf_file' => 'E90' ); } - if ($chem !~ /_t4s/ and $chem !~ /_t4ma/) { + if ($chem !~ /_ts4/ and $chem !~ /_t4ma/) { %species = (%species, 'BENZENE_an_srf_file' => 'BENZENE', 'BENZENE_bb_srf_file' => 'BENZENE', @@ -2450,10 +2428,10 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { 'XYLENES_an_srf_file' => 'XYLENES', 'XYLENES_bb_srf_file' => 'XYLENES' ) ; } - if ($chem =~ /trop_strat_mam4_ts2/ or $chem =~ /trop_strat_mam5_t2s1/) { + if ($chem =~ /trop_strat_mam4_ts2/ or $chem =~ /trop_strat_mam5_ts2/) { %species = (%species, 'MTERP_bb_srf_file' => 'APIN') ; - } elsif ($chem =~ /_t4s/ or $chem =~ /_t4ma/) { + } elsif ($chem =~ /_ts4/ or $chem =~ /_t4ma/) { %species = (%species, 'MTERP_bb_srf_file' => 'TERP') ; } else { @@ -2478,7 +2456,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { 'IVOC_bb_srf_file' => 'IVOCbb', 'SVOC_an_srf_file' => 'SVOCff', 'SVOC_bb_srf_file' => 'SVOCbb' ); - } elsif ($chem !~ /_t4s/ and $chem !~ /_t4ma/) { + } elsif ($chem !~ /_ts4/ and $chem !~ /_t4ma/) { %species = (%species, 'IVOC_an_srf_file' => 'IVOC', 'IVOC_bb_srf_file' => 'IVOC', @@ -2503,7 +2481,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { $first = 0; } } - if ($chem eq 'trop_mam4' or $chem eq 'waccm_sc_mam4' or $chem eq 'ghg_mam4' or $chem =~ /_t4s/ or $chem =~ /_t4ma/) { + if ($chem eq 'trop_mam4' or $chem eq 'waccm_sc_mam4' or $chem eq 'ghg_mam4' or $chem =~ /_ts4/ or $chem =~ /_t4ma/) { # SOA yields (used for the interactive emissions) have been calculated based on the VBS yields in CAM-chem. # Duseong S. Jo, et al. to be submitted to GMD, 2023 -- see https://github.com/ESCOMP/CAM/pull/727 discussion for additional detail. my %soae_fctrs = ('BENZENE_an_srf_file' => '2.5592D0', @@ -2679,7 +2657,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { add_default($nl, 'megan_factors_file'); add_default($nl, 'megan_mapped_emisfctrs', 'val'=>'.false.'); } - if ($chem =~ /trop_strat_mam4_vbs/ or $chem =~ /trop_strat_mam5_vbsext/ or $chem =~ /trop_strat_mam5_t1s/ or $chem =~ /trop_strat_mam5_slh/) { + if ($chem =~ /trop_strat_mam4_vbs/ or $chem =~ /trop_strat_mam5_vbs/) { my $val = "'ISOP = isoprene'," . "'MTERP = carene_3 + pinene_a + thujene_a + bornene + terpineol_4 + terpineol_a + terpinyl_ACT_a " . "+ myrtenal + sabinene + pinene_b + camphene + fenchene_a + limonene + phellandrene_a + terpinene_a " @@ -2718,7 +2696,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { add_default($nl, 'megan_factors_file'); add_default($nl, 'megan_mapped_emisfctrs', 'val'=>'.false.'); } - if ($chem =~ /trop_strat_mam5_t4s2/ or $chem =~ /_t4ma/) { + if ($chem =~ /trop_strat_mam5_ts4/ or $chem =~ /_t4ma/) { my $val = "'ISOP = isoprene'," . "'TERP = carene_3 + pinene_a + thujene_a + bornene + terpineol_4 + terpineol_a + terpinyl_ACT_a +'," . "' myrtenal + sabinene + pinene_b + camphene + fenchene_a + limonene + phellandrene_a + terpinene_a +'," @@ -2750,7 +2728,7 @@ if ($phys =~ /cam6/ or $phys =~ /cam7/) { add_default($nl, 'megan_factors_file'); add_default($nl, 'megan_mapped_emisfctrs', 'val'=>'.false.'); } - if ($chem =~ /trop_strat_mam4_ts2/ or $chem =~ /trop_strat_mam5_t2s1/) { + if ($chem =~ /trop_strat_mam4_ts2/ or $chem =~ /trop_strat_mam5_ts2/) { my $val = "'ISOP = isoprene'," . "'APIN = pinene_a + myrtenal'," . "'BPIN = carene_3 + thujene_a + bornene + fenchene_a + pinene_b + sabinene + camphene + terpineol_4 + terpineol_a + terpinyl_ACT_a'," @@ -3725,11 +3703,6 @@ if (!$simple_phys) { add_default($nl, 'do_iss'); } -# Sponge layer vertical diffusion -if (!$simple_phys) { - add_default($nl, 'diff_sponge_fac'); -} - # Convective water in radiation if (!$simple_phys) { add_default($nl, 'conv_water_in_rad'); @@ -3782,6 +3755,18 @@ if (!$simple_phys) { add_default($nl, 'zmconv_parcel_hscale'); } +# MCSP +if (!$simple_phys) { + add_default($nl, 'MCSP_heat_coeff'); + add_default($nl, 'MCSP_moisture_coeff'); + add_default($nl, 'MCSP_uwind_coeff'); + add_default($nl, 'MCSP_vwind_coeff'); + add_default($nl, 'MCSP_storm_speed_pref'); + add_default($nl, 'MCSP_conv_depth_min'); + add_default($nl, 'mcsp_shear_min'); +} + + # moist convection rainwater coefficients my $shallow_scheme = $nl->get_value('shallow_scheme'); $shallow_scheme =~ s/['"]//g; # strip quotes "' diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 0dc50f2f6b..72455ff0ca 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -2892,10 +2892,13 @@ 0.5 - 0.0D0 - 0.0D0 - 0.0D0 - 0.0D0 + 0.0D0 + 0.0D0 + 0.0D0 + 0.0D0 +60000.0 +50000.0 + 3.0 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 6bcc68520c..6eb9509b21 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -3678,30 +3678,50 @@ Default: 0.0 - + MCSP heating coefficient controlling the tendency intensity added. Default: 0.0 - + MCSP moistening coefficient controlling the tendency intensity added. Default: 0.0 - + MCSP uwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 - + MCSP vwnd tendency coefficient controlling the tendency intensity added. Default: 0.0 + +Reference pressure level in Pa for low-level zonal wind shear calculation for MCSP. +Default: 60000.0 + + + +Minimum convection depth for MCSP to be activated in Pa. +Default: 50000.0 + + + +Minimum uwind difference/shear between reference pressure level and surface for MCSP to be activated. +Default: 3.0 + + + + diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 6b56bf7df8..118cbc3c32 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -50,7 +50,9 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use cldfrc2m, only: cldfrc2m_readnl use rk_stratiform_cam, only: rk_stratiform_cam_readnl use zm_conv_intr, only: zm_conv_readnl - use zm_conv_mcsp, only: zm_conv_mcsp_readnl + !++ MCSP + use mcsp_intr, only: mcsp_readnl + !-- MCSP use hk_conv, only: hkconv_readnl use uwshcu, only: uwshcu_readnl use pkg_cld_sediment, only: cld_sediment_readnl @@ -154,7 +156,9 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call cldfrc_readnl(nlfilename) call cldfrc2m_readnl(nlfilename) call zm_conv_readnl(nlfilename) - call zm_conv_mcsp_readnl(nlfilename) + !++ MCSP + call mcsp_readnl(nlfilename) + !-- MCSP call rk_stratiform_cam_readnl(nlfilename) call hkconv_readnl(nlfilename) call uwshcu_readnl(nlfilename) diff --git a/src/physics/cam/convect_deep.F90 b/src/physics/cam/convect_deep.F90 index 920200b8b4..5328daac15 100644 --- a/src/physics/cam/convect_deep.F90 +++ b/src/physics/cam/convect_deep.F90 @@ -29,6 +29,12 @@ module convect_deep convect_deep_tend_2, &! return tendencies deep_scheme_does_scav_trans ! = .t. if scheme does scavenging and conv. transport +!++ MCSP + public :: jctop1 + + integer :: jctop1(pcols) +!-- MCSP + ! Private module data character(len=16) :: deep_scheme ! default set in phys_control.F90, use namelist to change ! Physics buffer indices @@ -163,7 +169,7 @@ end subroutine convect_deep_init !subroutine convect_deep_tend(state, ptend, tdt, pbuf) subroutine convect_deep_tend( & - mcon ,cme , & + mcon ,cme , & zdu , & rliq ,rice , & ztodt , & @@ -259,13 +265,17 @@ subroutine convect_deep_tend( & call pbuf_get_field(pbuf, pblh_idx, pblh) call pbuf_get_field(pbuf, tpert_idx, tpert) - call zm_conv_tend( pblh ,mcon ,cme , & + call zm_conv_tend( pblh ,mcon ,cme , & tpert ,zdu , & rliq ,rice , & ztodt , & jctop, jcbot , & state ,ptend ,landfrac, pbuf) + !++ MCSP + jctop1 = int(jctop(:pcols)) + !-- MCSP + end select ! If we added temperature tendency to pbuf, set it now. diff --git a/src/physics/cam/mcsp_intr.F90 b/src/physics/cam/mcsp_intr.F90 new file mode 100644 index 0000000000..10ba43c137 --- /dev/null +++ b/src/physics/cam/mcsp_intr.F90 @@ -0,0 +1,237 @@ +module mcsp_intr + !---------------------------------------------------------------------------- + ! Purpose: methods for mesoscale coherent structure parameterization (MCSP) + ! This scheme essentially redistributes the ZM heating and drying + ! tendencies vertically in order to mimic the effects of mesoscale + ! organization. It can also add momentum tendencies, although this + ! capability has not been extensively tested + !---------------------------------------------------------------------------- + ! References: + ! + ! Moncrieff, M. W., & Liu, C. (2006). Representing convective organization in + ! prediction models by a hybrid strategy. J. Atmos. Sci., 63, 3404–3420. + ! https://doi.org/10.1175/JAS3812.1 + ! + ! Chen, C.-C., Richter, J. H., Liu, C., Moncrieff, M. W., Tang, Q., Lin, W., + ! et al. (2021). Effects of organized convection parameterization on the MJO + ! and precipitation in E3SMv1. Part I: Mesoscale heating. J. Adv. + ! Mod. Earth Sys., 13, e2020MS002401, https://doi.org/10.1029/2020MS002401 + ! + ! Moncrieff, M. W., C. Liu, and P. Bogenschutz, 2017: Simulation, Modeling, + ! and Dynamically Based Parameterization of Organized Tropical Convection + ! for Global Climate Models. J. Atmos. Sci., 74, 1363–1380, + ! https://doi.org/10.1175/JAS-D-16-0166.1. + ! + ! Moncrieff, M. W. (2019). Toward a Dynamical Foundation for Organized Convection + ! Parameterization in GCMs. Geophys. Res. Lett., 46, 14103–14108. + ! https://doi.org/10.1029/2019GL085316 + ! + !---------------------------------------------------------------------------- + use shr_kind_mod, only: r8=>shr_kind_r8 + use cam_abortutils, only: endrun + use cam_logfile, only: iulog + + + implicit none + private + + public :: mcsp_register ! Register fields wit the output buffer + public :: mcsp_readnl ! Read MCSP namelist + public :: mcsp_intr_init ! initialize MCSP namelist variables + public :: mcsp_tend ! Perform MCSP tendency calculations + + real(r8) :: MCSP_storm_speed_pref ! pressure level for winds in MCSP calculation [Pa] + real(r8) :: MCSP_conv_depth_min ! pressure thickness of convective heating [Pa] + real(r8) :: mcsp_shear_min ! min shear value for MCSP to be active + + real(r8) :: MCSP_heat_coeff ! heating coefficient for MCSP + real(r8) :: MCSP_moisture_coeff ! moisture coefficient for MCSP + real(r8) :: MCSP_uwind_coeff ! uwind coefficient for MCSP + real(r8) :: MCSP_vwind_coeff ! vwind coefficient for MCSP + +!========================================================================================= +contains +!========================================================================================= +subroutine mcsp_register() + !---------------------------------------------------------------------------- + ! Purpose: register MCSP output fields + !---------------------------------------------------------------------------- + use cam_history, only: addfld, horiz_only + !---------------------------------------------------------------------------- + call addfld('MCSP_DT_max', horiz_only, 'A', 'K/s', 'MCSP max T tendency') + call addfld('MCSP_DT', (/'lev'/), 'A', 'K/s', 'MCSP T tendency') + call addfld('MCSP_DQ', (/'lev'/), 'A', 'kg/kg/s', 'MCSP qv tendency') + call addfld('MCSP_DU', (/'lev'/), 'A', 'm/s/day', 'MCSP U wind tendency') + call addfld('MCSP_DV', (/'lev'/), 'A', 'm/s/day', 'MCSP V wind tendency') + + call addfld('MCSP_freq', horiz_only, 'A', '1', 'MCSP frequency of activation') + call addfld('MCSP_shear', horiz_only, 'A', 'm/s', 'MCSP vertical shear of zonal wind') + call addfld('MCSP_conv_depth', horiz_only, 'A', 'Pa', 'ZM convection depth for MCSP') + +end subroutine mcsp_register + +!========================================================================================= + +subroutine mcsp_readnl(nlfile) + + use spmd_utils, only: mpicom, masterproc, masterprocid, mpi_real8, mpi_integer, mpi_logical + use namelist_utils, only: find_group_name + + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! Local variables + integer :: unitn, ierr + character(len=*), parameter :: subname = 'mcsp_readnl' + + namelist /mcsp_nl/ MCSP_heat_coeff, MCSP_moisture_coeff, & + MCSP_uwind_coeff, MCSP_vwind_coeff, & + MCSP_storm_speed_pref, MCSP_conv_depth_min, mcsp_shear_min + !----------------------------------------------------------------------------- + + if (masterproc) then + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'mcsp_nl', status=ierr) + if (ierr == 0) then + read(unitn, mcsp_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist') + end if + end if + close(unitn) + + end if + + ! Broadcast namelist variables + call mpi_bcast(MCSP_heat_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_heat_coeff") + call mpi_bcast(MCSP_moisture_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_moisture_coeff") + call mpi_bcast(MCSP_uwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_uwind_coeff") + call mpi_bcast(MCSP_vwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_vwind_coeff") + call mpi_bcast(MCSP_storm_speed_pref, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_storm_speed_pref") + call mpi_bcast(MCSP_conv_depth_min, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: MCSP_conv_depth_min") + call mpi_bcast(mcsp_shear_min, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun("mcsp_readnl: FATAL: mpi_bcast: mcsp_shear_min") + +end subroutine mcsp_readnl + +!=================================================================================================== + +subroutine mcsp_intr_init() + use MCSP, only: MCSP_init + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + + implicit none + + call MCSP_init(MCSP_heat_coeff, MCSP_moisture_coeff, MCSP_uwind_coeff, MCSP_vwind_coeff, & + MCSP_storm_speed_pref, MCSP_conv_depth_min, mcsp_shear_min, masterproc, iulog) + +end subroutine mcsp_intr_init + +!=================================================================================================== + +subroutine mcsp_tend( state, ptend, ztodt, jctop, tend_t, tend_q ) + + use MCSP, only: MCSP_run + use physconst, only: cpair + use physics_types, only: physics_state, physics_ptend, physics_ptend_init + use ppgrid, only: pver, pverp, pcols + use physconst, only: cpair, pi + use constituents, only: pcnst + + ! Arguments + type(physics_state), intent( in) :: state ! Physics state variables + type(physics_ptend), intent(out) :: ptend ! individual parameterization tendencies + real(r8), intent(in ) :: ztodt ! 2x physics time step + integer, dimension(pcols), intent(in ) :: jctop ! cloud top level indices + real(r8), dimension(pcols,pver), intent(in ) :: tend_t ! input deep convective temperature tendency + real(r8), dimension(pcols,pver), intent(in ) :: tend_q ! input deep convective tendency for specific humidity (qv) + + ! local variables + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns (actual) + real(r8), dimension(pcols,pver) :: mcsp_dt_out! final MCSP tendency for DSE + real(r8), dimension(pcols,pver) :: mcsp_dq_out! final MCSP tendency for qv + real(r8), dimension(pcols,pver) :: mcsp_du_out! final MCSP tendency for u wind + real(r8), dimension(pcols,pver) :: mcsp_dv_out! final MCSP tendency for v wind + real(r8), dimension(pcols) :: mcsp_freq ! MSCP frequency for output + real(r8), dimension(pcols) :: mcsp_shear ! shear used to check against threshold + real(r8), dimension(pcols) :: conv_depth ! pressure depth of deep convection + real(r8), dimension(pcols) :: mcsp_tend_s_max ! max MCSP heating tendency + logical :: lq(pcnst) + + !---------------------------------------------------------------------------- + ! Perform MCSP tendency calculations + !---------------------------------------------------------------------------- + lchnk = state%lchnk + ncol = state%ncol + lq(:) = .false. + lq(1) = .true. + + call physics_ptend_init(ptend, state%psetcols, 'MCSP', ls = .true., lq = lq, lu = .true., lv = .true.) + + call MCSP_run( pcols, ncol, pver, pverp, cpair, pi, & + ztodt, jctop, & + state%pmid, state%pint, state%pdel, & + state%s, state%q(:,:,1), state%u, state%v, & + tend_t, tend_q, & + ptend%s, ptend%q(:,:,1), ptend%u, ptend%v, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, conv_depth, mcsp_tend_s_max ) + + !---------------------------------------------------------------------------- + ! outpuf MCSP variables + !---------------------------------------------------------------------------- + call mcsp_hist( lchnk, pcols, pver, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, conv_depth, mcsp_tend_s_max ) + +end subroutine mcsp_tend + +!=================================================================================================== + +subroutine mcsp_hist( lchnk, pcols, pver, & + mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & + mcsp_freq, mcsp_shear, conv_depth, mcsp_tend_s_max ) + !---------------------------------------------------------------------------- + ! Purpose: write diagnostic quantities to history files + !---------------------------------------------------------------------------- + use cam_history, only: outfld + !---------------------------------------------------------------------------- + ! Arguments + integer, intent(in) :: lchnk ! chunk identifier + integer, intent(in) :: pcols ! number of atmospheric columns (max) + integer, intent(in) :: pver ! number of mid-point vertical levels + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dt_out ! final MCSP tendency for DSE + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dq_out ! final MCSP tendency for qv + real(r8), dimension(pcols,pver), intent(in) :: mcsp_du_out ! final MCSP tendency for u wind + real(r8), dimension(pcols,pver), intent(in) :: mcsp_dv_out ! final MCSP tendency for v wind + real(r8), dimension(pcols), intent(in) :: mcsp_freq ! MSCP frequency for output + real(r8), dimension(pcols), intent(in) :: mcsp_shear ! shear used to check against threshold + real(r8), dimension(pcols), intent(in) :: conv_depth ! pressure depth of ZM heating + real(r8), dimension(pcols), intent(in) :: mcsp_tend_s_max ! max MCSP heating tendency + !---------------------------------------------------------------------------- + ! write out MCSP diagnostic history fields + call outfld('MCSP_DT', mcsp_dt_out, pcols, lchnk ) + call outfld('MCSP_DQ', mcsp_dq_out, pcols, lchnk ) + call outfld('MCSP_DU', mcsp_du_out, pcols, lchnk ) + call outfld('MCSP_DV', mcsp_dv_out, pcols, lchnk ) + call outfld('MCSP_freq', mcsp_freq, pcols, lchnk ) + call outfld('MCSP_shear', mcsp_shear, pcols, lchnk ) + call outfld('MCSP_conv_depth', conv_depth, pcols, lchnk ) + call outfld('MCSP_DT_max', mcsp_tend_s_max, pcols, lchnk) + !---------------------------------------------------------------------------- + return + +end subroutine mcsp_hist + +!=================================================================================================== + +end module mcsp_intr + + diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index ad0837ea9b..c08d77962a 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -137,6 +137,9 @@ subroutine phys_register use ghg_data, only: ghg_data_register use vertical_diffusion, only: vd_register use convect_deep, only: convect_deep_register + !++ MCSP + use mcsp_intr, only: mcsp_register + !-- MCSP use convect_shallow, only: convect_shallow_register use radiation, only: radiation_register use co2_cycle, only: co2_register @@ -315,6 +318,9 @@ subroutine phys_register ! deep convection call convect_deep_register + !MCSP + call mcsp_register() + ! shallow convection call convect_shallow_register @@ -733,6 +739,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cldfrc2m, only: cldfrc2m_init use co2_cycle, only: co2_init, co2_transport use convect_deep, only: convect_deep_init + !++ MCSP + use mcsp_intr, only: mcsp_intr_init + !-- MCSP use convect_shallow, only: convect_shallow_init use constituents, only: cnst_get_ind use cam_diagnostics, only: diag_init @@ -921,6 +930,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call convect_deep_init(pref_edge) + !++ MCSP + call mcsp_intr_init() + !-- MCSP + if( microp_scheme == 'RK' ) then call rk_stratiform_cam_init() elseif( microp_scheme == 'MG' ) then @@ -1411,7 +1424,7 @@ subroutine tphysac (ztodt, cam_in, & use dycore, only: dycore_is use cam_control_mod, only: aqua_planet use mo_gas_phase_chemdr,only: map2chm - use clybryiy_fam, only: clybryiy_fam_set + use clybry_fam, only: clybry_fam_set use charge_neutrality, only: charge_balance use qbo, only: qbo_relax use iondrag, only: iondrag_calc, do_waccm_ions @@ -2073,7 +2086,7 @@ subroutine tphysac (ztodt, cam_in, & call diag_phys_tend_writeout (state, pbuf, tend, ztodt, qini, cldliqini, cldiceini) - call clybryiy_fam_set( ncol, lchnk, map2chm, state%q, pbuf ) + call clybry_fam_set( ncol, lchnk, map2chm, state%q, pbuf ) ! clean CARMA diagnostics object if (associated(carma_diags_obj)) then @@ -2158,7 +2171,7 @@ subroutine tphysbc (ztodt, state, & use cloud_diagnostics, only: cloud_diagnostics_calc use perf_mod use mo_gas_phase_chemdr,only: map2chm - use clybryiy_fam, only: clybryiy_fam_adj + use clybry_fam, only: clybry_fam_adj use clubb_intr, only: clubb_tend_cam use sslt_rebin, only: sslt_rebin_adv use tropopause, only: tropopause_output @@ -2176,6 +2189,14 @@ subroutine tphysbc (ztodt, state, & use dyn_tests_utils, only: vc_dycore use surface_emissions_mod,only: surface_emissions_set use elevated_emissions_mod,only: elevated_emissions_set + !++ MCSP + use mcsp_intr, only: mcsp_tend + use save_ttend_from_convect_deep, only : save_ttend_from_convect_deep_timestep_init, save_ttend_from_convect_deep_run + use save_qtend_from_convect_deep, only : save_qtend_from_convect_deep_timestep_init, save_qtend_from_convect_deep_run + use zm_conv_intr, only: ttend_s + use convect_deep, only: jctop1 + use physconst, only: cpair + !-- MCSP ! Arguments @@ -2213,6 +2234,12 @@ subroutine tphysbc (ztodt, state, & real(r8) dlf(pcols,pver) ! Detraining cld H20 from shallow + deep convections real(r8) dlf2(pcols,pver) ! Detraining cld H20 from shallow convections real(r8) rtdt ! 1./ztodt + !++ MCSP + real(r8) ttend_dp(pcols,pver) ! temperature tendency from deep convection + real(r8) qtend_dp(pcols,pver) ! water vapor from deep convection + character(len=512) :: errmsg + integer :: errflg + !-- MCSP integer lchnk ! chunk identifier integer ncol ! number of atmospheric columns @@ -2337,16 +2364,16 @@ subroutine tphysbc (ztodt, state, & if (state_debug_checks) & call physics_state_check(state, name="before tphysbc (dycore?)") - call clybryiy_fam_adj( ncol, lchnk, map2chm, state%q, pbuf ) + call clybry_fam_adj( ncol, lchnk, map2chm, state%q, pbuf ) - ! Since clybryiy_fam_adj operates directly on the tracers, and has no + ! Since clybry_fam_adj operates directly on the tracers, and has no ! physics_update call, re-run qneg3. call qneg3('TPHYSBCc',lchnk ,ncol ,pcols ,pver , & 1, pcnst, qmin ,state%q ) - ! Validate output of clybryiy_fam_adj. + ! Validate output of clybry_fam_adj. if (state_debug_checks) & - call physics_state_check(state, name="clybryiy_fam_adj") + call physics_state_check(state, name="clybry_fam_adj") ! ! Dump out "before physics" state @@ -2454,6 +2481,11 @@ subroutine tphysbc (ztodt, state, & flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) end if + !++ MCSP + call save_ttend_from_convect_deep_timestep_init(ncol, pver, ttend_dp, errmsg, errflg) + call save_qtend_from_convect_deep_timestep_init(ncol, pver, qtend_dp, errmsg, errflg) + !-- MCSP + call convect_deep_tend( & cmfmc, cmfcme, & zdu, & @@ -2461,6 +2493,11 @@ subroutine tphysbc (ztodt, state, & ztodt, & state, ptend, cam_in%landfrac, pbuf) + !++ MCSP + call save_ttend_from_convect_deep_run(ncol, pver, ttend_s, cpair, ttend_dp, errmsg, errflg) + call save_qtend_from_convect_deep_run(ncol, pver, ptend%q(:,:,1), qtend_dp, errmsg, errflg) + !-- MCSP + if ( (trim(cam_take_snapshot_after) == "convect_deep_tend") .and. & (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then call cam_snapshot_ptend_outfld(ptend, lchnk) @@ -2481,6 +2518,26 @@ subroutine tphysbc (ztodt, state, & call t_stopf('convect_deep_tend') + !++ MCSP + call t_startf ('MCSP_tend') + + if (trim(cam_take_snapshot_before) == "mcsp_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) + end if + + call mcsp_tend( state, ptend, ztodt, jctop1, ttend_dp, qtend_dp) + + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "mcsp_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + flx_heat, cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, det_s, det_ice, net_flx) + end if + + call t_stopf('MCSP_tend') + !-- MCSP + call pbuf_get_field(pbuf, prec_dp_idx, prec_dp ) call pbuf_get_field(pbuf, snow_dp_idx, snow_dp ) call pbuf_get_field(pbuf, prec_sh_idx, prec_sh ) @@ -3102,6 +3159,15 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) use phys_grid_ctem, only: phys_grid_ctem_diags use surface_emissions_mod,only: surface_emissions_adv use elevated_emissions_mod,only: elevated_emissions_adv + !++ MCSP + use mcsp_intr, only: mcsp_tend + use save_ttend_from_convect_deep, only : save_ttend_from_convect_deep_timestep_init, save_ttend_from_convect_deep_run + use save_qtend_from_convect_deep, only : save_qtend_from_convect_deep_timestep_init, save_qtend_from_convect_deep_run + use zm_conv_intr, only: ttend_s + use convect_deep, only: jctop1 + use physconst, only: cpair + !-- MCSP + implicit none diff --git a/src/physics/cam/zm_conv_MCSP.F90 b/src/physics/cam/zm_conv_MCSP.F90 deleted file mode 100644 index 803d91eeb3..0000000000 --- a/src/physics/cam/zm_conv_MCSP.F90 +++ /dev/null @@ -1,444 +0,0 @@ -module zm_conv_mcsp - !---------------------------------------------------------------------------- - ! Purpose: methods for mesoscale coherent structure parameterization (MCSP) - ! This scheme essentially redistributes the ZM heating and drying - ! tendencies vertically in order to mimic the effects of mesoscale - ! organization. It can also add momentum tendencies, although this - ! capability has not been extensively tested - !---------------------------------------------------------------------------- - ! References: - ! - ! Moncrieff, M. W., & Liu, C. (2006). Representing convective organization in - ! prediction models by a hybrid strategy. J. Atmos. Sci., 63, 3404–3420. - ! https://doi.org/10.1175/JAS3812.1 - ! - ! Chen, C.-C., Richter, J. H., Liu, C., Moncrieff, M. W., Tang, Q., Lin, W., - ! et al. (2021). Effects of organized convection parameterization on the MJO - ! and precipitation in E3SMv1. Part I: Mesoscale heating. J. Adv. - ! Mod. Earth Sys., 13, e2020MS002401, https://doi.org/10.1029/2020MS002401 - ! - ! Moncrieff, M. W., C. Liu, and P. Bogenschutz, 2017: Simulation, Modeling, - ! and Dynamically Based Parameterization of Organized Tropical Convection - ! for Global Climate Models. J. Atmos. Sci., 74, 1363–1380, - ! https://doi.org/10.1175/JAS-D-16-0166.1. - ! - ! Moncrieff, M. W. (2019). Toward a Dynamical Foundation for Organized Convection - ! Parameterization in GCMs. Geophys. Res. Lett., 46, 14103–14108. - ! https://doi.org/10.1029/2019GL085316 - ! - !---------------------------------------------------------------------------- - - use shr_kind_mod, only: r8=>shr_kind_r8 - use cam_abortutils, only: endrun - use cam_logfile, only: iulog - use physconst, only: cpair, pi - - - implicit none - private - - public :: zm_conv_mcsp_init ! Initialize MCSP output fields - public :: zm_conv_mcsp_readnl ! Read MCSP namelist - public :: zm_conv_mcsp_tend ! Perform MCSP tendency calculations - public :: zm_conv_mcsp_hist ! Write diagnostic quantities to history files - - real(r8), parameter :: MCSP_storm_speed_pref = 600e2_r8 ! pressure level for winds in MCSP calculation [Pa] - real(r8), parameter :: MCSP_conv_depth_min = 500e2_r8 ! pressure thickness of convective heating [Pa] - real(r8), parameter :: mcsp_shear_min = 3.0_r8 ! min shear value for MCSP to be active - real(r8), parameter :: mcsp_shear_max = 200.0_r8 ! max shear value for MCSP to be active - - real(r8) :: zmconv_MCSP_heat_coeff - real(r8) :: zmconv_MCSP_moisture_coeff - real(r8) :: zmconv_MCSP_uwind_coeff - real(r8) :: zmconv_MCSP_vwind_coeff - -!=================================================================================================== -contains -!=================================================================================================== - -subroutine zm_conv_mcsp_init() - !---------------------------------------------------------------------------- - ! Purpose: initialize MCSP output fields - !---------------------------------------------------------------------------- - use cam_history, only: addfld, horiz_only - use mpishorthand - !---------------------------------------------------------------------------- - !++ - call addfld('MCSP_DT_max', horiz_only, 'A', 'K/s', 'MCSP max T tendency') - !-- - call addfld('MCSP_DT', (/'lev'/), 'A', 'K/s', 'MCSP T tendency') - call addfld('MCSP_DQ', (/'lev'/), 'A', 'kg/kg/s', 'MCSP qv tendency') - call addfld('MCSP_DU', (/'lev'/), 'A', 'm/s/day', 'MCSP U wind tendency') - call addfld('MCSP_DV', (/'lev'/), 'A', 'm/s/day', 'MCSP V wind tendency') - - call addfld('MCSP_freq', horiz_only, 'A', '1', 'MCSP frequency of activation') - call addfld('MCSP_shear', horiz_only, 'A', 'm/s', 'MCSP vertical shear of zonal wind') - call addfld('MCSP_zm_depth', horiz_only, 'A', 'Pa', 'ZM convection depth for MCSP') - -end subroutine zm_conv_mcsp_init - - -!========================================================================================= - -subroutine zm_conv_mcsp_readnl(nlfile) - - use spmd_utils, only: mpicom, masterproc, masterprocid, mpi_real8, mpi_integer, mpi_logical - use namelist_utils, only: find_group_name - - character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - - ! Local variables - integer :: unitn, ierr - character(len=*), parameter :: subname = 'zm_conv_mcsp_readnl' - - namelist /zmconv_mcsp_nl/ zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & - zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff - !----------------------------------------------------------------------------- - - if (masterproc) then - open( newunit=unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'zmconv_mcsp_nl', status=ierr) - if (ierr == 0) then - read(unitn, zmconv_mcsp_nl, iostat=ierr) - if (ierr /= 0) then - call endrun(subname // ':: ERROR reading namelist') - end if - end if - close(unitn) - - end if - - ! Broadcast namelist variables - call mpi_bcast(zmconv_MCSP_heat_coeff, 1, mpi_integer, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_heat_coeff") - call mpi_bcast(zmconv_MCSP_moisture_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_moisture_coeff") - call mpi_bcast(zmconv_MCSP_uwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_uwind_coeff") - call mpi_bcast(zmconv_MCSP_vwind_coeff, 1, mpi_real8, masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun("zm_conv_mcsp_readnl: FATAL: mpi_bcast: zmconv_MCSP_vwind_coeff") - -end subroutine zm_conv_mcsp_readnl - -!=================================================================================================== - -subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear) - !---------------------------------------------------------------------------- - ! Purpose: calculate shear for MCSP - !---------------------------------------------------------------------------- - use interpolate_data, only: vertinterp - !---------------------------------------------------------------------------- - ! Arguments - integer, intent(in ) :: pcols ! number of atmospheric columns (max) - integer, intent(in ) :: ncol ! number of atmospheric columns (actual) - integer, intent(in ) :: pver ! number of mid-point vertical levels - real(r8), dimension(pcols,pver), intent(in ) :: state_pmid ! physics state mid-point pressure - real(r8), dimension(pcols,pver), intent(in ) :: state_u ! physics state u momentum - real(r8), dimension(pcols,pver), intent(in ) :: state_v ! physics state v momentum - real(r8), dimension(pcols), intent( out) :: mcsp_shear - !---------------------------------------------------------------------------- - ! Local variables - integer :: i - real(r8), dimension(pcols) :: storm_u ! u wind at storm reference level set by MCSP_storm_speed_pref - real(r8), dimension(pcols) :: storm_v ! v wind at storm reference level set by MCSP_storm_speed_pref - real(r8), dimension(pcols) :: storm_u_shear ! u shear at storm reference level set by MCSP_storm_speed_pref - real(r8), dimension(pcols) :: storm_v_shear ! v shear at storm reference level set by MCSP_storm_speed_pref - !---------------------------------------------------------------------------- - ! Interpolate wind to pressure level specified by MCSP_storm_speed_pref - call vertinterp( ncol, pcols, pver, state_pmid, MCSP_storm_speed_pref, state_u, storm_u ) - call vertinterp( ncol, pcols, pver, state_pmid, MCSP_storm_speed_pref, state_v, storm_v ) - - !---------------------------------------------------------------------------- - ! calculate low-level shear - do i = 1,ncol - if (state_pmid(i,pver).gt.MCSP_storm_speed_pref) then - storm_u_shear(i) = storm_u(i)-state_u(i,pver) - storm_v_shear(i) = storm_v(i)-state_v(i,pver) - else - storm_u_shear(i) = -999 - storm_v_shear(i) = -999 - end if - mcsp_shear(i) = storm_u_shear(i) - end do - - !---------------------------------------------------------------------------- - return - -end subroutine zm_conv_mcsp_calculate_shear - -!=================================================================================================== - -subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & - ztodt, jctop, & - state_pmid, state_pint, state_pdel, & - state_s, state_q, state_u, state_v, & - ptend_zm_s, ptend_zm_q, & - ptend_s, ptend_q, ptend_u, ptend_v, & - mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & - mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) - !---------------------------------------------------------------------------- - ! Purpose: perform MCSP tendency calculations - !---------------------------------------------------------------------------- - ! Arguments - integer, intent(in ) :: pcols ! number of atmospheric columns (max) - integer, intent(in ) :: ncol ! number of atmospheric columns (actual) - integer, intent(in ) :: pver ! number of mid-point vertical levels - integer, intent(in ) :: pverp ! number of interface vertical levels - real(r8), intent(in ) :: ztodt ! 2x physics time step - integer, dimension(pcols), intent(in ) :: jctop ! cloud top level indices - real(r8), dimension(pcols,pver), intent(in ) :: state_pmid ! physics state mid-point pressure - real(r8), dimension(pcols,pverp), intent(in ) :: state_pint ! physics state interface pressure - real(r8), dimension(pcols,pver), intent(in ) :: state_pdel ! physics state pressure thickness - real(r8), dimension(pcols,pver), intent(in ) :: state_s ! physics state dry static energy - real(r8), dimension(pcols,pver), intent(in ) :: state_q ! physics state specific humidity - real(r8), dimension(pcols,pver), intent(in ) :: state_u ! physics state u momentum - real(r8), dimension(pcols,pver), intent(in ) :: state_v ! physics state v momentum - real(r8), dimension(pcols,pver), intent(in ) :: ptend_zm_s ! input ZM tendency for dry static energy (DSE) - real(r8), dimension(pcols,pver), intent(in ) :: ptend_zm_q ! input ZM tendency for specific humidity (qv) - real(r8), dimension(pcols,pver), intent(inout) :: ptend_s ! output tendency of DSE - real(r8), dimension(pcols,pver), intent(inout) :: ptend_q ! output tendency of qv - real(r8), dimension(pcols,pver), intent(inout) :: ptend_u ! output tendency of u-wind - real(r8), dimension(pcols,pver), intent(inout) :: ptend_v ! output tendency of v-wind - real(r8), dimension(pcols,pver), intent( out) :: mcsp_dt_out! final MCSP tendency for DSE - real(r8), dimension(pcols,pver), intent( out) :: mcsp_dq_out! final MCSP tendency for qv - real(r8), dimension(pcols,pver), intent( out) :: mcsp_du_out! final MCSP tendency for u wind - real(r8), dimension(pcols,pver), intent( out) :: mcsp_dv_out! final MCSP tendency for v wind - real(r8), dimension(pcols), intent( out) :: mcsp_freq ! MSCP frequency for output - real(r8), dimension(pcols), intent( out) :: mcsp_shear ! shear used to check against threshold - real(r8), dimension(pcols), intent( out) :: zm_depth ! pressure depth of ZM heating - !++ - real(r8), dimension(pcols), intent( out) :: mcsp_tend_s_max ! max MCSP heating tendency - !-- - !---------------------------------------------------------------------------- - ! Local variables - integer :: i, k - - real(r8) :: tend_k ! temporary variable for kinetic energy tendency - real(r8) :: pdepth_mid_k ! temporary pressure depth used for vertical structure - real(r8) :: pdepth_total ! temporary pressure depth used for vertical structure - - real(r8), dimension(pcols) :: zm_avg_tend_s ! mass weighted column average DSE tendency from ZM - real(r8), dimension(pcols) :: zm_avg_tend_q ! mass weighted column average qv tendency from ZM - real(r8), dimension(pcols) :: pdel_sum ! column integrated pressure thickness - - real(r8), dimension(pcols,pver) :: mcsp_tend_s ! MCSP tendency before energy fixer for DSE - real(r8), dimension(pcols,pver) :: mcsp_tend_q ! MCSP tendency before energy fixer for qv - real(r8), dimension(pcols,pver) :: mcsp_tend_u ! MCSP tendency before energy fixer for u wind - real(r8), dimension(pcols,pver) :: mcsp_tend_v ! MCSP tendency before energy fixer for v wind - - real(r8), dimension(pcols) :: mcsp_avg_tend_s ! mass weighted column average MCSP tendency of DSE - real(r8), dimension(pcols) :: mcsp_avg_tend_q ! mass weighted column average MCSP tendency of qv - real(r8), dimension(pcols) :: mcsp_avg_tend_k ! mass weighted column average MCSP tendency of kinetic energy - - logical :: do_mcsp_t = .false. ! internal flag to enable tendency calculations - logical :: do_mcsp_q = .false. ! internal flag to enable tendency calculations - logical :: do_mcsp_u = .false. ! internal flag to enable tendency calculations - logical :: do_mcsp_v = .false. ! internal flag to enable tendency calculations - - !---------------------------------------------------------------------------- - ! initialize variables - - if (zmconv_MCSP_heat_coeff>0) do_mcsp_t = .true. - if (zmconv_MCSP_moisture_coeff>0) do_mcsp_q = .true. - if (zmconv_MCSP_uwind_coeff>0) do_mcsp_u = .true. - if (zmconv_MCSP_vwind_coeff>0) do_mcsp_v = .true. - - zm_avg_tend_s(1:ncol) = 0 - zm_avg_tend_q(1:ncol) = 0 - - pdel_sum(1:ncol) = 0 - - mcsp_avg_tend_s(1:ncol) = 0 - mcsp_avg_tend_q(1:ncol) = 0 - mcsp_avg_tend_k(1:ncol) = 0 - - mcsp_tend_s(1:ncol,1:pver) = 0 - mcsp_tend_q(1:ncol,1:pver) = 0 - mcsp_tend_u(1:ncol,1:pver) = 0 - mcsp_tend_v(1:ncol,1:pver) = 0 - !++ - mcsp_tend_s_max = 0._r8 - !-- - - if (zmconv_MCSP_heat_coeff>0 .OR. zmconv_MCSP_moisture_coeff>0 .OR. zmconv_MCSP_uwind_coeff>0 .OR. zmconv_MCSP_vwind_coeff>0) then - - !---------------------------------------------------------------------------- - ! calculate shear - - call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear ) - - !---------------------------------------------------------------------------- - ! calculate mass weighted column average tendencies from ZM - - zm_depth(1:ncol) = 0 - do i = 1,ncol - if (jctop(i) .ne. pver) then - ! integrate pressure and ZM tendencies over column - do k = jctop(i),pver - zm_avg_tend_s(i) = zm_avg_tend_s(i) + ptend_zm_s(i,k) * state_pdel(i,k) - zm_avg_tend_q(i) = zm_avg_tend_q(i) + ptend_zm_q(i,k) * state_pdel(i,k) - pdel_sum(i) = pdel_sum(i) + state_pdel(i,k) - end do - ! normalize integrated ZM tendencies by total mass - zm_avg_tend_s(i) = zm_avg_tend_s(i) / pdel_sum(i) - zm_avg_tend_q(i) = zm_avg_tend_q(i) / pdel_sum(i) - ! calculate diagnostic zm_depth - zm_depth(i) = state_pint(i,pver+1) - state_pmid(i,jctop(i)) - else - zm_avg_tend_s(i) = 0 - zm_avg_tend_q(i) = 0 - zm_depth(i) = 0 - end if - end do - - !---------------------------------------------------------------------------- - ! Note: To conserve total energy we need to account for the kinteic energy tendency - ! which we can obtain from the velocity tendencies based on the following: - ! KE_new = (u_new^2 + v_new^2)/2 - ! = [ (u_old+du)^2 + (v_old+dv)^2 ]/2 - ! = [ ( u_old^2 + 2*u_old*du + du^2 ) + ( v_old^2 + 2*v_old*dv + dv^2 ) ]/2 - ! = ( u_old^2 + v_old^2 )/2 + ( 2*u_old*du + du^2 + 2*v_old*dv + dv^2 )/2 - ! = KE_old + [ 2*u_old*du + du^2 + 2*v_old*dv + dv^2 ] /2 - - !---------------------------------------------------------------------------- - ! calculate MCSP tendencies - - do i = 1,ncol - - ! check that ZM produced tendencies over a depth that exceeds the threshold - if ( zm_depth(i) >= MCSP_conv_depth_min ) then - ! check that ZM provided a non-zero column total heating - if ( zm_avg_tend_s(i) > 0 ) then - ! check that there is sufficient wind shear to justify coherent organization - if ( abs(mcsp_shear(i)).ge.mcsp_shear_min .and. & - abs(mcsp_shear(i)).lt.mcsp_shear_max ) then - do k = jctop(i),pver - - ! See eq 7-8 of Moncrieff et al. (2017) - also eq (5) of Moncrieff & Liu (2006) - pdepth_mid_k = state_pint(i,pver+1) - state_pmid(i,k) - pdepth_total = state_pint(i,pver+1) - state_pmid(i,jctop(i)) - - ! specify the assumed vertical structure - if (do_mcsp_t) mcsp_tend_s(i,k) = -1*zmconv_MCSP_heat_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) - if (do_mcsp_q) mcsp_tend_q(i,k) = -1*zmconv_MCSP_moisture_coeff * sin(2.0_r8*pi*(pdepth_mid_k/pdepth_total)) - if (do_mcsp_u) mcsp_tend_u(i,k) = zmconv_MCSP_uwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) - if (do_mcsp_v) mcsp_tend_v(i,k) = zmconv_MCSP_vwind_coeff * (cos(pi*(pdepth_mid_k/pdepth_total))) - - ! scale the vertical structure by the ZM heating/drying tendencies - if (do_mcsp_t) mcsp_tend_s(i,k) = zm_avg_tend_s(i) * mcsp_tend_s(i,k) - !++ - if (do_mcsp_t) mcsp_tend_s_max(i) = zm_avg_tend_s(i) * zmconv_MCSP_heat_coeff / cpair - !-- - if (do_mcsp_q) mcsp_tend_q(i,k) = zm_avg_tend_q(i) * mcsp_tend_q(i,k) - - ! integrate the DSE/qv tendencies for energy/mass fixer - if (do_mcsp_t) mcsp_avg_tend_s(i) = mcsp_avg_tend_s(i) + mcsp_tend_s(i,k) * state_pdel(i,k) / pdel_sum(i) - if (do_mcsp_q) mcsp_avg_tend_q(i) = mcsp_avg_tend_q(i) + mcsp_tend_q(i,k) * state_pdel(i,k) / pdel_sum(i) - - ! integrate the change in kinetic energy (KE) for energy fixer - if (do_mcsp_u.or.do_mcsp_v) then - tend_k = ( 2.0_r8*mcsp_tend_u(i,k)*ztodt*state_u(i,k) + mcsp_tend_u(i,k)*mcsp_tend_u(i,k)*ztodt*ztodt & - +2.0_r8*mcsp_tend_v(i,k)*ztodt*state_v(i,k) + mcsp_tend_v(i,k)*mcsp_tend_v(i,k)*ztodt*ztodt )/2.0_r8/ztodt - mcsp_avg_tend_k(i) = mcsp_avg_tend_k(i) + tend_k*state_pdel(i,k) / pdel_sum(i) - end if - - end do ! k = jctop(i),pver - end if ! shear threshold - end if ! zm_avg_tend_s(i) > 0 - end if ! zm_depth(i) >= MCSP_conv_depth_min - end do - - end if - - !---------------------------------------------------------------------------- - ! calculate final output tendencies - - mcsp_dt_out(1:ncol,1:pver) = 0 - mcsp_dq_out(1:ncol,1:pver) = 0 - mcsp_du_out(1:ncol,1:pver) = 0 - mcsp_dv_out(1:ncol,1:pver) = 0 - - mcsp_freq(1:ncol) = 0 - - do i = 1,ncol - do k = jctop(i),pver - - ! update frequency if MCSP contributes any tendency in the column - if ( abs(mcsp_tend_s(i,k)).gt.0 .or. abs(mcsp_tend_q(i,k)).gt.0 .or.& - abs(mcsp_tend_u(i,k)).gt.0 .or. abs(mcsp_tend_v(i,k)).gt.0 ) then - mcsp_freq(i) = 1 - end if - - ! subtract mass weighted average tendencies for energy/mass conservation - mcsp_dt_out(i,k) = mcsp_tend_s(i,k) - mcsp_avg_tend_s(i) - mcsp_dq_out(i,k) = mcsp_tend_q(i,k) - mcsp_avg_tend_q(i) - mcsp_du_out(i,k) = mcsp_tend_u(i,k) - mcsp_dv_out(i,k) = mcsp_tend_v(i,k) - - ! make sure kinetic energy correction is added to DSE tendency - ! to conserve total energy whenever momentum tendencies are calculated - if (do_mcsp_u.or.do_mcsp_v) then - mcsp_dt_out(i,k) = mcsp_dt_out(i,k) - mcsp_avg_tend_k(i) - end if - - ! update output tendencies - if (do_mcsp_t) ptend_s(i,k) = ptend_s(i,k) + mcsp_dt_out(i,k) - if (do_mcsp_q) ptend_q(i,k) = ptend_q(i,k) + mcsp_dq_out(i,k) - if (do_mcsp_u) ptend_u(i,k) = ptend_u(i,k) + mcsp_du_out(i,k) - if (do_mcsp_v) ptend_v(i,k) = ptend_v(i,k) + mcsp_dv_out(i,k) - - ! adjust units for diagnostic outputs - if (do_mcsp_t) mcsp_dt_out(i,k) = mcsp_dt_out(i,k)/cpair - - end do - end do - - !---------------------------------------------------------------------------- - return - -end subroutine zm_conv_mcsp_tend - -!=================================================================================================== - -subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & - mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & - mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) - !---------------------------------------------------------------------------- - ! Purpose: write diagnostic quantities to history files - !---------------------------------------------------------------------------- - use cam_history, only: outfld - !---------------------------------------------------------------------------- - ! Arguments - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: pcols ! number of atmospheric columns (max) - integer, intent(in) :: pver ! number of mid-point vertical levels - real(r8), dimension(pcols,pver), intent(in) :: mcsp_dt_out ! final MCSP tendency for DSE - real(r8), dimension(pcols,pver), intent(in) :: mcsp_dq_out ! final MCSP tendency for qv - real(r8), dimension(pcols,pver), intent(in) :: mcsp_du_out ! final MCSP tendency for u wind - real(r8), dimension(pcols,pver), intent(in) :: mcsp_dv_out ! final MCSP tendency for v wind - real(r8), dimension(pcols), intent(in) :: mcsp_freq ! MSCP frequency for output - real(r8), dimension(pcols), intent(in) :: mcsp_shear ! shear used to check against threshold - real(r8), dimension(pcols), intent(in) :: zm_depth ! pressure depth of ZM heating - !++ - real(r8), dimension(pcols), intent(in) :: mcsp_tend_s_max ! max MCSP heating tendency - !-- - !---------------------------------------------------------------------------- - ! write out MCSP diagnostic history fields - call outfld('MCSP_DT', mcsp_dt_out, pcols, lchnk ) - call outfld('MCSP_DQ', mcsp_dq_out, pcols, lchnk ) - call outfld('MCSP_DU', mcsp_du_out, pcols, lchnk ) - call outfld('MCSP_DV', mcsp_dv_out, pcols, lchnk ) - call outfld('MCSP_freq', mcsp_freq, pcols, lchnk ) - call outfld('MCSP_shear', mcsp_shear, pcols, lchnk ) - call outfld('MCSP_zm_depth', zm_depth, pcols, lchnk ) - !++ - call outfld('MCSP_DT_max', mcsp_tend_s_max, pcols, lchnk) - !-- - !---------------------------------------------------------------------------- - return - -end subroutine zm_conv_mcsp_hist - -!=================================================================================================== - -end module zm_conv_mcsp diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 18e1b48890..6f7789d000 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -40,6 +40,12 @@ module zm_conv_intr zm_conv_tend, &! return tendencies zm_conv_tend_2 ! return tendencies + !++ MCSP + public :: ttend_s + + real(r8) :: ttend_s(pcols,pver) + !-- MCSP + public zmconv_ke, zmconv_ke_lnd ! needed by convect_shallow integer ::& ! indices for fields in the physics buffer @@ -78,6 +84,7 @@ module zm_conv_intr real(r8) :: zmconv_parcel_hscale = unset_r8 ! Fraction of PBL depth over which to mix initial parcel real(r8) :: zmconv_tau = unset_r8 ! Timescale for convection + ! indices for fields in the physics buffer integer :: cld_idx = 0 integer :: icwmrdp_idx = 0 @@ -224,9 +231,6 @@ subroutine zm_conv_init(pref_edge) use spmd_utils, only: masterproc use phys_control, only: phys_deepconv_pbl, phys_getopts, cam_physpkg_is use physics_buffer, only: pbuf_get_index -!++ MCSP - use zm_conv_mcsp, only: zm_conv_mcsp_init -!-- MCSP implicit none @@ -310,10 +314,6 @@ subroutine zm_conv_init(pref_edge) call add_default('ZMMTT ', history_budget_histfile_num, ' ') end if -!++ MCSP - call zm_conv_mcsp_init() -!-- MCSP - ! ! Limit deep convection to regions below 40 mb ! Note this calculation is repeated in the shallow convection interface @@ -389,9 +389,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & use phys_control, only: cam_physpkg_is use ccpp_constituent_prop_mod, only: ccpp_const_props -!++ MCSP - use zm_conv_mcsp, only: zm_conv_mcsp_tend, zm_conv_mcsp_hist -!-- MCSP ! Arguments @@ -436,28 +433,6 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & type(physics_state) :: state1 ! locally modify for evaporation to use, not returned type(physics_ptend),target :: ptend_loc ! package tendencies -!++ MCSP - type(physics_ptend),target :: ptend_mcsp ! MCSP output tendencies - - ! flags for MCSP tendencies - logical :: do_mcsp_t = .false. - logical :: do_mcsp_q(pcnst) = .false. - logical :: do_mcsp_u = .false. - logical :: do_mcsp_v = .false. - - ! MCSP history output variables - real(r8), dimension(pcols,pver) :: mcsp_dt_out ! MCSP tendency for DSE - real(r8), dimension(pcols,pver) :: mcsp_dq_out ! MCSP tendency for qv - real(r8), dimension(pcols,pver) :: mcsp_du_out ! MCSP tendency for u wind - real(r8), dimension(pcols,pver) :: mcsp_dv_out ! MCSP tendency for v wind - real(r8), dimension(pcols) :: mcsp_freq ! MSCP frequency for output - real(r8), dimension(pcols) :: mcsp_shear ! shear used to check against threshold - real(r8), dimension(pcols) :: zm_depth ! pressure depth of ZM heating - !++ MCSP - real(r8), dimension(pcols) :: mcsp_tend_s_max ! max MCSP heating tendency - !-- MCSP -!-- MCSP - ! physics buffer fields real(r8), pointer, dimension(:) :: prec ! total precipitation real(r8), pointer, dimension(:) :: snow ! snow from ZM convection @@ -625,44 +600,9 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & call outfld('CAPE', cape, pcols, lchnk) ! RBN - CAPE output -!++ MCSP - !---------------------------------------------------------------------------- - ! mesoscale coherent structure parameterization (MCSP) - ! Note that this modifies the tendencies produced by zm_convr(), such that - ! history variables like ZMDT will include the effects of MCSP - - ! Set these all to True so that the tedency variables get allocated. The internal flags to - ! calculate each MCSP tedency will behave the same, but having them always allocated avoids - - do_mcsp_t = .true. - do_mcsp_q(1) = .true. - do_mcsp_u = .true. - do_mcsp_v = .true. - - call physics_ptend_init( ptend_mcsp, state%psetcols, 'zm_conv_mcsp_tend', & - ls=do_mcsp_t, lq=do_mcsp_q, lu=do_mcsp_u, lv=do_mcsp_v) - - call zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & - ztodt, int(jctop), & - state%pmid, state%pint, state%pdel, & - state%s, state%q, state%u, state%v, & - ptend_loc%s, ptend_loc%q(:,:,1), & - ptend_mcsp%s(:,:), ptend_mcsp%q(:,:,1), & - ptend_mcsp%u(:,:), ptend_mcsp%v(:,:), & - mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & - mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) - - call zm_conv_mcsp_hist( lchnk, pcols, pver, & - mcsp_dt_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, & - mcsp_freq, mcsp_shear, zm_depth, mcsp_tend_s_max ) - - ! add MCSP tendencies to ZM convective tendencies - call physics_ptend_sum( ptend_mcsp, ptend_loc, ncol) - call physics_ptend_dealloc(ptend_mcsp) - -!-- MCSP - - + !++ MCSP + ttend_s = ptend_loc%s(:pcols,:) + !-- MCSP ! ! Output fractional occurance of ZM convection ! diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index a3f0c49fe6..06daae3b82 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -133,6 +133,9 @@ subroutine phys_register use ghg_data, only: ghg_data_register use vertical_diffusion, only: vd_register use convect_deep, only: convect_deep_register + !++ MCSP + use mcsp_intr, only: mcsp_register + !-- MCSP use convect_diagnostics,only: convect_diagnostics_register use radiation, only: radiation_register use co2_cycle, only: co2_register @@ -305,6 +308,9 @@ subroutine phys_register ! deep convection call convect_deep_register + !MCSP + call mcsp_register() + ! convection diagnostics call convect_diagnostics_register @@ -731,6 +737,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cldfrc2m, only: cldfrc2m_init use co2_cycle, only: co2_init, co2_transport use convect_deep, only: convect_deep_init + !++ MCSP + use mcsp_intr, only: mcsp_intr_init + !-- MCSP use convect_diagnostics,only: convect_diagnostics_init use cam_diagnostics, only: diag_init use gw_drag_cam, only: gw_drag_cam_init @@ -917,6 +926,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call convect_deep_init(pref_edge) + !++ MCSP + call mcsp_intr_init() + !-- MCSP + if (.not. do_clubb_sgs) call macrop_driver_init(pbuf2d) call microp_aero_init(phys_state,pbuf2d) call microp_driver_init(pbuf2d) @@ -1402,7 +1415,7 @@ subroutine tphysac (ztodt, cam_in, & use dycore, only: dycore_is use cam_control_mod, only: aqua_planet use mo_gas_phase_chemdr,only: map2chm - use clybryiy_fam, only: clybryiy_fam_set + use clybry_fam, only: clybry_fam_set use charge_neutrality, only: charge_balance use qbo, only: qbo_relax use iondrag, only: iondrag_calc, do_waccm_ions @@ -1487,6 +1500,7 @@ subroutine tphysac (ztodt, cam_in, & real(r8) dlf(pcols,pver) ! Detraining cld H20 from shallow + deep convections real(r8) rtdt ! 1./ztodt + real(r8) :: rliq(pcols) ! vertical integral of liquid not yet in q(ixcldliq) real(r8) :: det_s (pcols) ! vertical integral of detrained static energy from ice real(r8) :: det_ice(pcols) ! vertical integral of detrained ice @@ -2556,7 +2570,7 @@ subroutine tphysac (ztodt, cam_in, & call diag_phys_tend_writeout (state, pbuf, tend, ztodt, qini, cldliqini, cldiceini) - call clybryiy_fam_set( ncol, lchnk, map2chm, state%q, pbuf ) + call clybry_fam_set( ncol, lchnk, map2chm, state%q, pbuf ) ! output these here -- after updates by chem_timestep_tend or export_fields within the current time step if (associated(cam_out%nhx_nitrogen_flx)) then @@ -2620,7 +2634,7 @@ subroutine tphysbc (ztodt, state, & use radiation, only: radiation_tend use perf_mod use mo_gas_phase_chemdr,only: map2chm - use clybryiy_fam, only: clybryiy_fam_adj + use clybry_fam, only: clybry_fam_adj use cam_abortutils, only: endrun use subcol_utils, only: is_subcol_on use qneg_module, only: qneg3 @@ -2629,6 +2643,14 @@ subroutine tphysbc (ztodt, state, & use dyn_tests_utils, only: vc_dycore use surface_emissions_mod,only: surface_emissions_set use elevated_emissions_mod,only: elevated_emissions_set + !++ MCSP + use mcsp_intr, only: mcsp_tend + use save_ttend_from_convect_deep, only : save_ttend_from_convect_deep_timestep_init, save_ttend_from_convect_deep_run + use save_qtend_from_convect_deep, only : save_qtend_from_convect_deep_timestep_init, save_qtend_from_convect_deep_run + use zm_conv_intr, only: ttend_s + use convect_deep, only: jctop1 + use physconst, only: cpair + !-- MCSP ! Arguments @@ -2660,6 +2682,12 @@ subroutine tphysbc (ztodt, state, & real(r8) dlf(pcols,pver) ! Detraining cld H20 from shallow + deep convections real(r8) dlf2(pcols,pver) ! Detraining cld H20 from shallow convections real(r8) rtdt ! 1./ztodt + !++ MCSP + real(r8) ttend_dp(pcols,pver) ! temperature tendency from deep convection + real(r8) qtend_dp(pcols,pver) ! water vapor from deep convection + character(len=512) :: errmsg + integer :: errflg + !-- MCSP integer lchnk ! chunk identifier integer ncol ! number of atmospheric columns @@ -2767,16 +2795,16 @@ subroutine tphysbc (ztodt, state, & call physics_state_check(state, name="before tphysbc (dycore?)") end if - call clybryiy_fam_adj( ncol, lchnk, map2chm, state%q, pbuf ) + call clybry_fam_adj( ncol, lchnk, map2chm, state%q, pbuf ) - ! Since clybryiy_fam_adj operates directly on the tracers, and has no + ! Since clybry_fam_adj operates directly on the tracers, and has no ! physics_update call, re-run qneg3. call qneg3('TPHYSBCc',lchnk ,ncol ,pcols ,pver , & 1, pcnst, qmin ,state%q ) - ! Validate output of clybryiy_fam_adj. + ! Validate output of clybry_fam_adj. if (state_debug_checks) then - call physics_state_check(state, name="clybryiy_fam_adj") + call physics_state_check(state, name="clybry_fam_adj") end if ! ! Dump out "before physics" state @@ -2884,6 +2912,11 @@ subroutine tphysbc (ztodt, state, & cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) end if + !++ MCSP + call save_ttend_from_convect_deep_timestep_init(ncol, pver, ttend_dp, errmsg, errflg) + call save_qtend_from_convect_deep_timestep_init(ncol, pver, qtend_dp, errmsg, errflg) + !-- MCSP + call convect_deep_tend( & cmfmc, cmfcme, & zdu, & @@ -2891,6 +2924,11 @@ subroutine tphysbc (ztodt, state, & ztodt, & state, ptend, cam_in%landfrac, pbuf) + !++ MCSP + call save_ttend_from_convect_deep_run(ncol, pver, ttend_s, cpair, ttend_dp, errmsg, errflg) + call save_qtend_from_convect_deep_run(ncol, pver, ptend%q(:,:,1), qtend_dp, errmsg, errflg) + !-- MCSP + if ( (trim(cam_take_snapshot_after) == "convect_deep_tend") .and. & (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then call cam_snapshot_ptend_outfld(ptend, lchnk) @@ -2911,6 +2949,26 @@ subroutine tphysbc (ztodt, state, & call t_stopf('convect_deep_tend') + !++ MCSP + call t_startf ('MCSP_tend') + if (trim(cam_take_snapshot_before) == "mcsp_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call mcsp_tend( state, ptend, ztodt, jctop1, ttend_dp, qtend_dp) + + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "mcsp_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call t_stopf('MCSP_tend') + !-- MCSP + + call pbuf_get_field(pbuf, prec_dp_idx, prec_dp ) call pbuf_get_field(pbuf, snow_dp_idx, snow_dp ) call pbuf_get_field(pbuf, prec_sh_idx, prec_sh )