diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist index 65dcb193626c..823531c565bf 100755 --- a/components/eam/bld/build-namelist +++ b/components/eam/bld/build-namelist @@ -4049,8 +4049,8 @@ add_default($nl, 'zmconv_accr_fac'); add_default($nl, 'zmconv_micro_dcs'); add_default($nl, 'zmconv_MCSP_heat_coeff'); add_default($nl, 'zmconv_MCSP_moisture_coeff'); -add_default($nl, 'zmconv_MCSP_uwind_coeff'); -add_default($nl, 'zmconv_MCSP_vwind_coeff'); +add_default($nl, 'zmconv_MCSP_mom_coeff'); +add_default($nl, 'zmconv_MCSP_use_full_shear'); # # moist convection rainwater coefficients # These no more needed for EAM, including uwshcu_rpen. No need to reset for 'cam5' diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index b42f03b8564f..a56f8a4c7a7c 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -1270,8 +1270,8 @@ 0.0D0 0.0D0 -0.0D0 -0.0D0 +0.0D0 +.false. diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index 258c3ef00e27..2f3a91816be7 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -3075,15 +3075,17 @@ Multiscale Coherent System Parameterization (MCSP), moisture coefficient Default: set by build-namelist - -Multiscale Coherent System Parameterization (MCSP), zonal wind coefficient +Multiscale Coherent System Parameterization (MCSP), momentum coefficient +(dimensionless fraction of the storm-relative shear) Default: set by build-namelist - -Multiscale Coherent System Parameterization (MCSP), meridional wind coefficient +Multiscale Coherent System Parameterization (MCSP), use the full (u,v) shear +vector instead of zonal-only shear. Default (.false.) preserves E3SMv3 behavior. Default: set by build-namelist diff --git a/components/eam/docs/user-guide/namelist_parameters.md b/components/eam/docs/user-guide/namelist_parameters.md index 5864fcad6383..78a50e23c774 100644 --- a/components/eam/docs/user-guide/namelist_parameters.md +++ b/components/eam/docs/user-guide/namelist_parameters.md @@ -135,20 +135,20 @@ | Conv. micro. Parameters | Description | Default value | | ------------------------- | ----------------------------------------------------------------- | ---------------------- | | `zmconv_microp` | Convective microphysics option in ZM convection scheme | `true` | -| `zmconv_auto_fac` | Cloud droplet-rain autoconversion enhancement factor in the convective microphysics scheme | `7.0` | -| `zmconv_accr_fac` | Cloud droplet-rain accretion enhancement factor in the convective microphysics scheme | `1.5` | -| `zmconv_micro_dcs` | Autoconversion size threshold for cloud ice to snow (m) | `150.E-6` | +| `zmconv_auto_fac` | Cloud droplet-rain autoconversion factor in the ZM microphysics | `7.0` | +| `zmconv_accr_fac` | Cloud droplet-rain accretion factor in the ZM microphysics | `1.5` | +| `zmconv_micro_dcs` | Autoconversion size threshold for cloud ice to snow (m) | `150.E-6` | | Mass flux adj. Parameters | Description | Default value | | ------------------------- | ----------------------------------------------------------------- | ---------------------- | -| `zmconv_clos_dyn_adj` | Apply mass flux adjustment to ZM convection scheme | `true` | - -| MCSP Parameters | Description | Default value | -| ---------------------------- | ----------------------------------------------------------------- | ---------------------- | -| `zmconv_mcsp_heat_coeff` | MCSP heating coefficient | `0.3` | -| `zmconv_mcsp_moisture_coeff` | MCSP moisture coefficient | `0.0` | -| `zmconv_mcsp_uwind_coeff` | MCSP zonal wind coefficient | `0.0` | -| `zmconv_mcsp_vwind_coeff` | MCSP meridional wind coefficient | `0.0` | +| `zmconv_clos_dyn_adj` | Apply mass flux adjustment to ZM convection scheme | `true` | + +| MCSP Parameters | Description | Default value | +| ---------------------------- | ---------------------------------------------------------------| ---------------------- | +| `zmconv_mcsp_heat_coeff` | MCSP heating coefficient | `0.3` | +| `zmconv_mcsp_moisture_coeff` | MCSP moisture coefficient | `0.0` | +| `zmconv_mcsp_mom_coeff` | MCSP momentum coefficient (fraction of storm-relative shear) | `0.0` | +| `zmconv_mcsp_use_full_shear` | MCSP use full (u,v) shear vector instead of zonal-only shear | `.false.` | ## Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package diff --git a/components/eam/src/physics/cam/zm/zm_conv.F90 b/components/eam/src/physics/cam/zm/zm_conv.F90 index 11a5de6d0ac4..4002e593c572 100644 --- a/components/eam/src/physics/cam/zm/zm_conv.F90 +++ b/components/eam/src/physics/cam/zm/zm_conv.F90 @@ -50,7 +50,8 @@ module zm_conv public zm_conv_main_init ! ZM scheme initialization public zm_conv_main ! ZM scheme calculations public zm_conv_evap ! ZM scheme evaporation of precip -#ifdef SCREAM_CONFIG_IS_CMAKE ! Public for testing only +#ifdef SCREAM_CONFIG_IS_CMAKE + ! Public for testing only public zm_calc_fractional_entrainment public zm_downdraft_properties public zm_cloud_properties @@ -110,7 +111,7 @@ end subroutine zm_conv_main_init !=================================================================================================== -subroutine zm_get_gather_index(pcols, ncol, pver, pverp, is_first_step, cape, dcape, & +subroutine zm_get_gather_index(pcols, ncol, is_first_step, cape, dcape, & cape_threshold_loc, gather_index, lengath) !---------------------------------------------------------------------------- ! Purpose: determine length of gathered arrays @@ -118,8 +119,6 @@ subroutine zm_get_gather_index(pcols, ncol, pver, pverp, is_first_step, cape, dc ! Arguments integer, intent(in ) :: pcols ! maximum number of columns integer, intent(in ) :: ncol ! actual number of columns - integer, intent(in ) :: pver ! number of mid-point vertical levels - integer, intent(in ) :: pverp ! number of interface vertical levels logical(btype), intent(in ) :: is_first_step ! flag for first step of run real(r8), dimension(pcols), intent(in ) :: cape ! conv. avail. potential energy [J] real(r8), dimension(pcols), intent(in ) :: dcape ! CAPE generated by dycor (dCAPE) [J] @@ -128,7 +127,7 @@ subroutine zm_get_gather_index(pcols, ncol, pver, pverp, is_first_step, cape, dc integer, intent( out) :: lengath ! number of columns for gathered arrays !---------------------------------------------------------------------------- ! Local variables - integer :: i, ii + integer :: i !---------------------------------------------------------------------------- ! set local threshold to be used for zm_closure() if ( zm_param%trig_dcape .and. (.not.is_first_step) ) then @@ -239,7 +238,6 @@ subroutine zm_conv_main(pcols, ncol, pver, pverp, is_first_step, time_step, & real(r8), dimension(pcols) :: t_pcl_lcl ! parcel temperature at lcl [K] integer, dimension(pcols) :: lcl ! base level index of deep cumulus convection integer, dimension(pcols) :: lel ! index of highest theoretical convective plume - integer, dimension(pcols) :: lon ! index of onset level for deep convection integer, dimension(pcols) :: msemax_klev ! index of level with largest moist static energy real(r8), dimension(pcols,pver) :: t_pcl_m1 ! time n-1 parcel temperatures @@ -247,7 +245,6 @@ subroutine zm_conv_main(pcols, ncol, pver, pverp, is_first_step, time_step, & real(r8), dimension(pcols) :: t_pcl_lcl_m1 ! time n-1 parcel Temperature at LCL integer, dimension(pcols) :: lcl_m1 ! time n-1 base level index of deep cumulus convection integer, dimension(pcols) :: lel_m1 ! time n-1 index of highest theoretical convective plume - integer, dimension(pcols) :: lon_m1 ! time n-1 index of onset level for deep convection integer, dimension(pcols) :: msemax_klev_m1 ! time n-1 index of level with largest moist static energy real(r8), dimension(pcols) :: cape_m1 ! time n-1 CAPE @@ -295,7 +292,7 @@ subroutine zm_conv_main(pcols, ncol, pver, pverp, is_first_step, time_step, & type(zm_microp_st) :: loc_microp_st ! local (gathered) convective microphysics state and tendencies - integer i, ii, k, kk ! loop iterators + integer i, k ! loop iterators integer msg ! number of missing moisture levels at the top of model real(r8) qdifr @@ -440,7 +437,7 @@ subroutine zm_conv_main(pcols, ncol, pver, pverp, is_first_step, time_step, & !---------------------------------------------------------------------------- ! determine whether active columns for gathering - call zm_get_gather_index(pcols, ncol, pver, pverp, is_first_step, cape, dcape, & + call zm_get_gather_index(pcols, ncol, is_first_step, cape, dcape, & cape_threshold_loc, gather_index, lengath) if (lengath.eq.0) then ! Deallocate local microphysics arrays before returning diff --git a/components/eam/src/physics/cam/zm/zm_conv_intr.F90 b/components/eam/src/physics/cam/zm/zm_conv_intr.F90 index aa1be347e52d..1fbdb5e3cb6f 100644 --- a/components/eam/src/physics/cam/zm/zm_conv_intr.F90 +++ b/components/eam/src/physics/cam/zm/zm_conv_intr.F90 @@ -109,17 +109,17 @@ subroutine zm_conv_readnl(nlfile) real(r8) :: zmconv_micro_dcs = unset_r8 real(r8) :: zmconv_MCSP_heat_coeff = 0._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 + real(r8) :: zmconv_MCSP_mom_coeff = 0._r8 + logical :: zmconv_MCSP_use_full_shear = .false. !---------------------------------------------------------------------------- namelist /zmconv_nl/ zmconv_tau, zmconv_alfa, zmconv_ke, zmconv_dmpdz, & zmconv_tpert_fix, zmconv_tp_fac, zmconv_tiedke_add, & - zmconv_c0_lnd, zmconv_c0_ocn, & + zmconv_c0_lnd, zmconv_c0_ocn, & zmconv_cape_cin, zmconv_mx_bot_lyr_adj, & zmconv_trig_dcape, zmconv_trig_ull, zmconv_clos_dyn_adj, & zmconv_microp, zmconv_auto_fac, zmconv_accr_fac, zmconv_micro_dcs, & zmconv_MCSP_heat_coeff, zmconv_MCSP_moisture_coeff, & - zmconv_MCSP_uwind_coeff, zmconv_MCSP_vwind_coeff + zmconv_MCSP_mom_coeff, zmconv_MCSP_use_full_shear !---------------------------------------------------------------------------- if (masterproc) then @@ -157,12 +157,12 @@ subroutine zm_conv_readnl(nlfile) zm_param%micro_dcs = zmconv_micro_dcs ! mesoscale coherent structure parameterization (MCSP) parameters - zm_param%mcsp_t_coeff = zmconv_MCSP_heat_coeff - zm_param%mcsp_q_coeff = zmconv_MCSP_moisture_coeff - zm_param%mcsp_u_coeff = zmconv_MCSP_uwind_coeff - zm_param%mcsp_v_coeff = zmconv_MCSP_vwind_coeff + zm_param%mcsp_t_coeff = zmconv_MCSP_heat_coeff + zm_param%mcsp_q_coeff = zmconv_MCSP_moisture_coeff + zm_param%mcsp_mom_coeff = zmconv_MCSP_mom_coeff + zm_param%mcsp_use_full_shear = zmconv_MCSP_use_full_shear if ( abs(zm_param%mcsp_t_coeff)>0._r8 .or. abs(zm_param%mcsp_q_coeff)>0._r8 .or. & - abs(zm_param%mcsp_u_coeff)>0._r8 .or. abs(zm_param%mcsp_v_coeff)>0._r8 ) then + abs(zm_param%mcsp_mom_coeff)>0._r8 ) then zm_param%mcsp_enabled = .true. else zm_param%mcsp_enabled = .false. @@ -483,8 +483,7 @@ subroutine zm_conv_tend(pblh, mcon, cme, tpert, dlftot, pflx, zdu, & ! 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. + logical :: do_mcsp_mom = .false. ! MCSP history output variables real(r8), dimension(pcols,pver) :: mcsp_dt_out ! MCSP tendency for DSE @@ -674,11 +673,10 @@ subroutine zm_conv_tend(pblh, mcon, cme, tpert, dlftot, pflx, zdu, & ! 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. + do_mcsp_mom = .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) + ls=do_mcsp_t, lq=do_mcsp_q, lu=do_mcsp_mom, lv=do_mcsp_mom) call zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & ztodt, jctop, zm_const, zm_param, & diff --git a/components/eam/src/physics/cam/zm/zm_conv_mcsp.F90 b/components/eam/src/physics/cam/zm/zm_conv_mcsp.F90 index d7f38a7cb047..3319b45d64a0 100644 --- a/components/eam/src/physics/cam/zm/zm_conv_mcsp.F90 +++ b/components/eam/src/physics/cam/zm/zm_conv_mcsp.F90 @@ -44,7 +44,9 @@ module zm_conv_mcsp public :: zm_conv_mcsp_init ! Initialize MCSP output fields public :: zm_conv_mcsp_tend ! Perform MCSP tendency calculations +#ifndef SCREAM_CONFIG_IS_CMAKE public :: zm_conv_mcsp_hist ! Write diagnostic quantities to history files +#endif #ifdef SCREAM_CONFIG_IS_CMAKE public :: zm_conv_mcsp_calculate_shear ! just for testing #endif @@ -80,9 +82,11 @@ end subroutine zm_conv_mcsp_init !=================================================================================================== -subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear) +subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, shear_u, shear_v) !---------------------------------------------------------------------------- - ! Purpose: calculate shear for MCSP + ! Purpose: calculate storm-relative shear vector for MCSP. The zonal and + ! meridional components are returned separately; the caller decides whether to + ! gate/scale on the zonal component alone (legacy) or the full shear magnitude. !---------------------------------------------------------------------------- #ifdef SCREAM_CONFIG_IS_CMAKE use zm_eamxx_bridge_methods, only: vertinterp @@ -97,23 +101,29 @@ subroutine zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, 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 + real(r8), dimension(pcols), intent( out) :: shear_u ! zonal component of storm-relative shear + real(r8), dimension(pcols), intent( out) :: shear_v ! meridional component of storm-relative 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 !---------------------------------------------------------------------------- - ! Interpolate wind to pressure level specified by MCSP_storm_speed_pref + ! Interpolate winds 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 + ! calculate low-level shear components. The -999 sentinel (when the surface is + ! above the storm reference level) fails the caller's shear-magnitude gate. do i = 1,ncol if (state_pmid(i,pver).gt.MCSP_storm_speed_pref) then - mcsp_shear(i) = storm_u(i)-state_u(i,pver) + shear_u(i) = storm_u(i)-state_u(i,pver) + shear_v(i) = storm_v(i)-state_v(i,pver) else - mcsp_shear(i) = -999 + shear_u(i) = -999 + shear_v(i) = 0 end if end do @@ -185,10 +195,13 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & 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 + real(r8), dimension(pcols) :: shear_u ! zonal component of storm-relative shear + real(r8), dimension(pcols) :: shear_v ! meridional component of storm-relative shear + real(r8) :: cos_struct ! vertical structure factor for momentum tendencies + 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 + logical :: do_mcsp_mom = .false. ! internal flag to enable momentum tendency calculations !---------------------------------------------------------------------------- @@ -197,10 +210,12 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & !---------------------------------------------------------------------------- ! initialize variables - if (zm_param%mcsp_t_coeff>0) do_mcsp_t = .true. - if (zm_param%mcsp_q_coeff>0) do_mcsp_q = .true. - if (zm_param%mcsp_u_coeff>0) do_mcsp_u = .true. - if (zm_param%mcsp_v_coeff>0) do_mcsp_v = .true. + if (zm_param%mcsp_t_coeff>0) do_mcsp_t = .true. + if (zm_param%mcsp_q_coeff>0) do_mcsp_q = .true. + ! momentum is active for any non-zero coefficient; a negative coefficient + ! flips the leading sign in the tendency, switching up-gradient forcing + ! (default) to down-gradient forcing. + if (abs(zm_param%mcsp_mom_coeff)>0) do_mcsp_mom = .true. zm_avg_tend_s(1:ncol) = 0 zm_avg_tend_q(1:ncol) = 0 @@ -219,7 +234,7 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & !---------------------------------------------------------------------------- ! calculate shear - call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear ) + call zm_conv_mcsp_calculate_shear( pcols, ncol, pver, state_pmid, state_u, state_v, shear_u, shear_v ) !---------------------------------------------------------------------------- ! calculate mass weighted column average tendencies from ZM @@ -259,6 +274,17 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & do i = 1,ncol + ! The momentum forcing always follows the (u,v) shear vector (du~shear_u, + ! dv~shear_v). The mcsp_use_full_shear option only controls the scalar shear + ! used for the activation threshold and diagnostic: the full shear magnitude, + ! or the legacy signed zonal shear (which reproduces E3SMv3, since the same + ! shear also gates the T/q tendencies). + if (zm_param%mcsp_use_full_shear) then + mcsp_shear(i) = sqrt( shear_u(i)*shear_u(i) + shear_v(i)*shear_v(i) ) + else + mcsp_shear(i) = shear_u(i) ! legacy signed zonal shear diagnostic + end if + ! 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 @@ -275,8 +301,18 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & ! specify the assumed vertical structure if (do_mcsp_t) mcsp_tend_s(i,k) = -1*zm_param%mcsp_t_coeff * bfb_sin(2.0_r8*zm_const%pi*(pdepth_mid_k/pdepth_total)) if (do_mcsp_q) mcsp_tend_q(i,k) = -1*zm_param%mcsp_q_coeff * bfb_sin(2.0_r8*zm_const%pi*(pdepth_mid_k/pdepth_total)) - if (do_mcsp_u) mcsp_tend_u(i,k) = zm_param%mcsp_u_coeff * (bfb_cos(zm_const%pi*(pdepth_mid_k/pdepth_total))) - if (do_mcsp_v) mcsp_tend_v(i,k) = zm_param%mcsp_v_coeff * (bfb_cos(zm_const%pi*(pdepth_mid_k/pdepth_total))) + ! Momentum tendencies are scaled by the storm-relative shear vector so + ! mcsp_mom_coeff is a dimensionless O(0.01-0.1) fraction rather than a + ! raw acceleration. The implied wind increment over the step is + ! du = -mcsp_mom_coeff * shear_u * cos(...) -> tend_u = du / ztodt + ! The leading minus sign makes the forcing up-gradient (amplifying the + ! shear), consistent with organized-convection momentum transport + ! (Moncrieff), and bounds the KE correction. + if (do_mcsp_mom) then + cos_struct = bfb_cos(zm_const%pi*(pdepth_mid_k/pdepth_total)) + mcsp_tend_u(i,k) = -1*zm_param%mcsp_mom_coeff * shear_u(i) * cos_struct / ztodt + mcsp_tend_v(i,k) = -1*zm_param%mcsp_mom_coeff * shear_v(i) * cos_struct / ztodt + end if ! 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) @@ -287,7 +323,7 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & 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 + if (do_mcsp_mom) 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) @@ -326,15 +362,15 @@ subroutine zm_conv_mcsp_tend( pcols, ncol, pver, pverp, & ! 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 + if (do_mcsp_mom) 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) + 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_mom) ptend_u(i,k) = ptend_u(i,k) + mcsp_du_out(i,k) + if (do_mcsp_mom) 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)/zm_const%cpair @@ -349,15 +385,14 @@ end subroutine zm_conv_mcsp_tend !=================================================================================================== +#ifndef SCREAM_CONFIG_IS_CMAKE 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 ) !---------------------------------------------------------------------------- ! 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 @@ -371,7 +406,6 @@ subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & 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 !---------------------------------------------------------------------------- -#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 ) @@ -380,11 +414,11 @@ subroutine zm_conv_mcsp_hist( lchnk, pcols, pver, & 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 ) -#endif !---------------------------------------------------------------------------- return end subroutine zm_conv_mcsp_hist +#endif !=================================================================================================== diff --git a/components/eam/src/physics/cam/zm/zm_conv_types.F90 b/components/eam/src/physics/cam/zm/zm_conv_types.F90 index 9c3c53a1dea8..ddae1f303f76 100644 --- a/components/eam/src/physics/cam/zm/zm_conv_types.F90 +++ b/components/eam/src/physics/cam/zm/zm_conv_types.F90 @@ -77,11 +77,11 @@ module zm_conv_types real(r8) :: accr_fac = unset_r8 ! ZM microphysics enhancement factor for droplet-rain accretion real(r8) :: micro_dcs = unset_r8 ! ZM microphysics size threshold for cloud ice to snow autoconversion [m] ! MCSP parameters - logical :: mcsp_enabled = .false. ! flag for mesoscale coherent structure parameterization (MSCP) - real(r8) :: mcsp_t_coeff = 0 ! MCSP coefficient for temperature tendencies - real(r8) :: mcsp_q_coeff = 0 ! MCSP coefficient for specific humidity tendencies - real(r8) :: mcsp_u_coeff = 0 ! MCSP coefficient for zonal momentum tendencies - real(r8) :: mcsp_v_coeff = 0 ! MCSP coefficient for meridional momentum tendencies + logical :: mcsp_enabled = .false. ! flag for mesoscale coherent structure parameterization (MSCP) + real(r8) :: mcsp_t_coeff = 0 ! MCSP coefficient for temperature tendencies + real(r8) :: mcsp_q_coeff = 0 ! MCSP coefficient for specific humidity tendencies + real(r8) :: mcsp_mom_coeff = 0 ! MCSP momentum coefficient (dimensionless fraction of the shear) + logical :: mcsp_use_full_shear = .false. ! use full (u,v) shear vector instead of zonal-only shear (default preserves E3SMv3) end type zm_param_t !=================================================================================================== @@ -213,11 +213,14 @@ subroutine zm_param_set_for_testing(zm_param) zm_param%accr_fac = 1.5D0 zm_param%micro_dcs = 150.E-6 ! MCSP parameters - zm_param%mcsp_enabled = .true. - zm_param%mcsp_t_coeff = 0.3 - zm_param%mcsp_q_coeff = 0 - zm_param%mcsp_u_coeff = 0 - zm_param%mcsp_v_coeff = 0 + zm_param%mcsp_enabled = .true. + zm_param%mcsp_t_coeff = 0.3 + zm_param%mcsp_q_coeff = 0 + ! Use a negative momentum coefficient so the unit tests exercise the + ! do_mcsp_mom path AND confirm the guards admit a negative coefficient + ! (down-gradient forcing). Must match the C++ zm_opts_init test setup. + zm_param%mcsp_mom_coeff = -0.05 + zm_param%mcsp_use_full_shear = .true. end subroutine zm_param_set_for_testing !=================================================================================================== @@ -256,8 +259,8 @@ subroutine zm_param_mpi_broadcast(zm_param) call mpibcast(zm_param%mcsp_enabled, 1, mpilog, 0, mpicom) ! MCSP parameters call mpibcast(zm_param%mcsp_t_coeff, 1, mpir8, 0, mpicom) call mpibcast(zm_param%mcsp_q_coeff, 1, mpir8, 0, mpicom) - call mpibcast(zm_param%mcsp_u_coeff, 1, mpir8, 0, mpicom) - call mpibcast(zm_param%mcsp_v_coeff, 1, mpir8, 0, mpicom) + call mpibcast(zm_param%mcsp_mom_coeff, 1, mpir8, 0, mpicom) + call mpibcast(zm_param%mcsp_use_full_shear, 1, mpilog, 0, mpicom) #endif end subroutine zm_param_mpi_broadcast @@ -300,8 +303,8 @@ subroutine zm_param_print(zm_param) write(iulog,*) indent,'mcsp_enabled : ',zm_param%mcsp_enabled write(iulog,*) indent,'mcsp_t_coeff : ',zm_param%mcsp_t_coeff write(iulog,*) indent,'mcsp_q_coeff : ',zm_param%mcsp_q_coeff - write(iulog,*) indent,'mcsp_u_coeff : ',zm_param%mcsp_u_coeff - write(iulog,*) indent,'mcsp_v_coeff : ',zm_param%mcsp_v_coeff + write(iulog,*) indent,'mcsp_mom_coeff : ',zm_param%mcsp_mom_coeff + write(iulog,*) indent,'mcsp_use_full_shear : ',zm_param%mcsp_use_full_shear write(iulog,*) '' call shr_sys_flush(iulog) end if ! masterproc diff --git a/components/eamxx/cime_config/namelist_defaults_eamxx.xml b/components/eamxx/cime_config/namelist_defaults_eamxx.xml index 29c61a6260fe..054afdf29b73 100644 --- a/components/eamxx/cime_config/namelist_defaults_eamxx.xml +++ b/components/eamxx/cime_config/namelist_defaults_eamxx.xml @@ -289,11 +289,11 @@ be lost if SCREAM_HACK_XML is not enabled. 0.002 1 - true - 0.3 - 0.0 - 0.0 - 0.0 + true + 0.3 + 0.0 + 0.0 + false diff --git a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp b/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp index 9cf88d10f326..b5ddc420ec06 100644 --- a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp +++ b/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp @@ -138,10 +138,10 @@ void ZMDeepConvection::initialize_impl (const RunType) zm_output.init_host_mirrors(m_ncol, m_nlev); // initialize variables on the fortran side zm::zm_eamxx_bridge_init( m_nlev, ZMF::s_zm_opts.limcnv+1, - ZMF::s_zm_opts.trig_dcape, - ZMF::s_zm_opts.trig_ull, - ZMF::s_zm_opts.clos_dyn_adj, - ZMF::s_zm_opts.mcsp_enabled); + ZMF::s_zm_opts.trig_dcape, ZMF::s_zm_opts.trig_ull, + ZMF::s_zm_opts.clos_dyn_adj, ZMF::s_zm_opts.mcsp_enabled, + ZMF::s_zm_opts.mcsp_t_coeff, ZMF::s_zm_opts.mcsp_q_coeff, + ZMF::s_zm_opts.mcsp_mom_coeff, ZMF::s_zm_opts.mcsp_use_full_shear ); } // if use_fortran_bridge //---------------------------------------------------------------------------- } // ZMDeepConvection::initialize_impl diff --git a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.cpp b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.cpp index 8312423be48c..eb104d6ec04e 100644 --- a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.cpp +++ b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.cpp @@ -9,7 +9,9 @@ extern "C" { void zm_eamxx_bridge_init_c( Int pver_in, Int limcnv_in, bool trig_dcape_in, bool trig_ull_in, - bool clos_dyn_adj_in, bool mcsp_enabled_in ); + bool clos_dyn_adj_in, bool mcsp_enabled_in, + Real mcsp_t_coeff_in, Real mcsp_q_coeff_in, + Real mcsp_mom_coeff_in, bool mcsp_use_full_shear_in ); void zm_eamxx_bridge_run_c( Int ncol, // 01 Real dtime, // 02 @@ -63,10 +65,14 @@ namespace zm { void zm_eamxx_bridge_init( Int pver, Int limcnv_in, bool trig_dcape_in, bool trig_ull_in, - bool clos_dyn_adj_in, bool mcsp_enabled_in ){ + bool clos_dyn_adj_in, bool mcsp_enabled_in, + Real mcsp_t_coeff_in, Real mcsp_q_coeff_in, + Real mcsp_mom_coeff_in, bool mcsp_use_full_shear_in ){ zm_eamxx_bridge_init_c( pver, limcnv_in, trig_dcape_in, trig_ull_in, - clos_dyn_adj_in, mcsp_enabled_in ); + clos_dyn_adj_in, mcsp_enabled_in, + mcsp_t_coeff_in, mcsp_q_coeff_in, + mcsp_mom_coeff_in, mcsp_use_full_shear_in ); } void zm_eamxx_bridge_run( Int ncol, Int pver, Real dtime, bool is_first_step, diff --git a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.hpp b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.hpp index c5270fae6a10..782506fbc110 100644 --- a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.hpp +++ b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge.hpp @@ -16,7 +16,9 @@ using ZMF = zm::Functions; // Glue functions to call fortran from from C++ with the Data struct void zm_eamxx_bridge_init( Int pver, Int limcnv_in, bool trig_dcape_in, bool trig_ull_in, - bool clos_dyn_adj_in, bool mcsp_enabled_in ); + bool clos_dyn_adj_in, bool mcsp_enabled_in, + Real mcsp_t_coeff_in, Real mcsp_q_coeff_in, + Real mcsp_mom_coeff_in, bool mcsp_use_full_shear_in ); void zm_eamxx_bridge_run( Int ncol, Int pver, Real dtime, bool is_first_step, ZMF::ZmInputState& zm_input, ZMF::ZmOutputTend& zm_output, diff --git a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge_main.F90 b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge_main.F90 index a1f90149cbec..2889387f2459 100644 --- a/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge_main.F90 +++ b/components/eamxx/src/physics/zm/fortran_bridge/zm_eamxx_bridge_main.F90 @@ -28,7 +28,9 @@ module zm_eamxx_bridge_main subroutine zm_eamxx_bridge_init_c( pver_in, limcnv_in, & trig_dcape_in, trig_ull_in, & - clos_dyn_adj_in, mcsp_enabled_in ) bind(C) + clos_dyn_adj_in, mcsp_enabled_in, & + mcsp_t_coeff_in, mcsp_q_coeff_in, & + mcsp_mom_coeff_in, mcsp_use_full_shear_in ) bind(C) use mpi use zm_conv, only: zm_const, zm_param use zm_conv_types, only: zm_const_set_for_testing, zm_param_set_for_testing @@ -42,6 +44,10 @@ subroutine zm_eamxx_bridge_init_c( pver_in, limcnv_in, & logical(kind=c_bool),value, intent(in) :: trig_ull_in logical(kind=c_bool),value, intent(in) :: clos_dyn_adj_in logical(kind=c_bool),value, intent(in) :: mcsp_enabled_in + real(kind=c_real), value, intent(in) :: mcsp_t_coeff_in + real(kind=c_real), value, intent(in) :: mcsp_q_coeff_in + real(kind=c_real), value, intent(in) :: mcsp_mom_coeff_in + logical(kind=c_bool),value, intent(in) :: mcsp_use_full_shear_in !----------------------------------------------------------------------------- ! Local variables integer :: mpi_rank, ierror @@ -67,6 +73,10 @@ subroutine zm_eamxx_bridge_init_c( pver_in, limcnv_in, & zm_param%trig_ull = trig_ull_in zm_param%clos_dyn_adj = clos_dyn_adj_in zm_param%mcsp_enabled = mcsp_enabled_in + zm_param%mcsp_t_coeff = mcsp_t_coeff_in + zm_param%mcsp_q_coeff = mcsp_q_coeff_in + zm_param%mcsp_mom_coeff = mcsp_mom_coeff_in + zm_param%mcsp_use_full_shear = mcsp_use_full_shear_in call zm_param_print(zm_param) !----------------------------------------------------------------------------- call wv_sat_init() diff --git a/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_calculate_shear_impl.hpp b/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_calculate_shear_impl.hpp index d296b4f6d22e..79da5f7ce859 100644 --- a/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_calculate_shear_impl.hpp +++ b/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_calculate_shear_impl.hpp @@ -19,15 +19,20 @@ void Functions::zm_conv_mcsp_calculate_shear( const Int& pver, // number of mid-point vertical levels const uview_1d& state_pmid, // physics state mid-point pressure const uview_1d& state_u, // physics state u momentum + const uview_1d& state_v, // physics state v momentum // Outputs - Real& mcsp_shear) + Real& shear_u, // zonal component of storm-relative shear + Real& shear_v) // meridional component of storm-relative shear { //---------------------------------------------------------------------------- - // Purpose: calculate shear for MCSP + // Purpose: calculate storm-relative shear vector for MCSP. The zonal and + // meridional components are returned separately; the caller decides whether to + // gate/scale on the zonal component alone (legacy) or the full shear magnitude. //---------------------------------------------------------------------------- // Local variables Real storm_u = 0; // u wind at storm reference level set by MCSP_storm_speed_pref + Real storm_v = 0; // v wind at storm reference level set by MCSP_storm_speed_pref //---------------------------------------------------------------------------- // Interpolate wind to pressure level specified by MCSP_storm_speed_pref @@ -47,29 +52,35 @@ void Functions::zm_conv_mcsp_calculate_shear( // Linear interpolation - computed redundantly on all team threads so the // result is consistent across the team (no single/broadcast needed). Using - // Kokkos::single here would update mcsp_shear on only one thread, leaving the + // Kokkos::single here would update the outputs on only one thread, leaving the // other threads with stale values and causing a divergent branch (and GPU - // deadlock) where mcsp_shear later gates a team-collective. + // deadlock) where the shear later gates a team-collective. // (removing the Kokkos::single fixed a hang issue) if (state_pmid(0) >= ZMC::MCSP_storm_speed_pref) { storm_u = state_u(0); + storm_v = state_v(0); } else if (state_pmid(pver - 1) < ZMC::MCSP_storm_speed_pref) { storm_u = state_u(pver - 1); + storm_v = state_v(pver - 1); } else { EKAT_KERNEL_ASSERT(k_below < pver-1); const Real dpu = ZMC::MCSP_storm_speed_pref - state_pmid(k_below); const Real dpl = state_pmid(k_below+1) - ZMC::MCSP_storm_speed_pref; storm_u = (state_u(k_below)*dpl + state_u(k_below+1)*dpu) / (dpl + dpu); + storm_v = (state_v(k_below)*dpl + state_v(k_below+1)*dpu) / (dpl + dpu); } //---------------------------------------------------------------------------- - // calculate low-level shear + // calculate low-level shear components. The -999 sentinel (when the surface is + // above the storm reference level) fails the caller's shear-magnitude gate. if (state_pmid(pver - 1) > ZMC::MCSP_storm_speed_pref) { - mcsp_shear = storm_u - state_u(pver - 1); + shear_u = storm_u - state_u(pver - 1); + shear_v = storm_v - state_v(pver - 1); } else { - mcsp_shear = -999; + shear_u = -999; + shear_v = 0; } } diff --git a/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_tend_impl.hpp b/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_tend_impl.hpp index 4ff8c0b35bb9..96e005651f8e 100644 --- a/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_tend_impl.hpp +++ b/components/eamxx/src/physics/zm/impl/zm_conv_mcsp_tend_impl.hpp @@ -54,10 +54,12 @@ void Functions::zm_conv_mcsp_tend( //---------------------------------------------------------------------------- // initialize variables - const bool do_mcsp_t = (runtime_opt.mcsp_t_coeff > 0); - const bool do_mcsp_q = (runtime_opt.mcsp_q_coeff > 0); - const bool do_mcsp_u = (runtime_opt.mcsp_u_coeff > 0); - const bool do_mcsp_v = (runtime_opt.mcsp_v_coeff > 0); + const bool do_mcsp_t = (runtime_opt.mcsp_t_coeff > 0); + const bool do_mcsp_q = (runtime_opt.mcsp_q_coeff > 0); + // momentum is active for any non-zero coefficient; a negative coefficient + // flips the leading sign in the tendency, switching up-gradient forcing + // (default) to down-gradient forcing. + const bool do_mcsp_mom = (runtime_opt.mcsp_mom_coeff != 0); // Allocate temporary arrays uview_1d mcsp_tend_s, mcsp_tend_q, mcsp_tend_u, mcsp_tend_v; @@ -86,9 +88,21 @@ void Functions::zm_conv_mcsp_tend( team.team_barrier(); //---------------------------------------------------------------------------- - // calculate shear - - zm_conv_mcsp_calculate_shear(team, pver, state_pmid, state_u, mcsp_shear); + // calculate shear vector and derive the scalar metric used for gating/scaling + + Real shear_u = 0; // zonal component of storm-relative shear + Real shear_v = 0; // meridional component of storm-relative shear + zm_conv_mcsp_calculate_shear(team, pver, state_pmid, state_u, state_v, shear_u, shear_v); + + // The momentum forcing always follows the (u,v) shear vector (du~shear_u, + // dv~shear_v). The mcsp_use_full_shear option only controls the scalar shear + // used for the activation threshold and diagnostic: the full shear magnitude, + // or the legacy signed zonal shear (which reproduces E3SMv3, since the same + // shear also gates the T/q tendencies). The gate below uses abs(mcsp_shear), + // which equals the magnitude either way. + mcsp_shear = runtime_opt.mcsp_use_full_shear + ? std::sqrt(shear_u * shear_u + shear_v * shear_v) + : shear_u; //---------------------------------------------------------------------------- // calculate mass weighted column average tendencies from ZM @@ -143,8 +157,19 @@ void Functions::zm_conv_mcsp_tend( // specify the assumed vertical structure if (do_mcsp_t) mcsp_tend_s(k) = -1 * runtime_opt.mcsp_t_coeff * std::sin(2 * PC::Pi * (pdepth_mid_k / pdepth_total)); if (do_mcsp_q) mcsp_tend_q(k) = -1 * runtime_opt.mcsp_q_coeff * std::sin(2 * PC::Pi * (pdepth_mid_k / pdepth_total)); - if (do_mcsp_u) mcsp_tend_u(k) = runtime_opt.mcsp_u_coeff * (std::cos(PC::Pi * (pdepth_mid_k / pdepth_total))); - if (do_mcsp_v) mcsp_tend_v(k) = runtime_opt.mcsp_v_coeff * (std::cos(PC::Pi * (pdepth_mid_k / pdepth_total))); + // Momentum tendencies are scaled by the storm-relative shear vector so + // mcsp_mom_coeff is a dimensionless O(0.01-0.1) fraction rather than a + // raw acceleration. The implied wind increment over the step is + // du = -mcsp_mom_coeff * shear_u * cos(...) -> tend_u = du / ztodt + // The leading minus sign makes the forcing up-gradient (amplifying the + // shear), consistent with organized-convection momentum transport + // (Moncrieff). Using the shear components keeps the tendency aligned + // with the shear vector and bounds the KE correction. + if (do_mcsp_mom) { + const Real cos_struct = std::cos(PC::Pi * (pdepth_mid_k / pdepth_total)); + mcsp_tend_u(k) = -1 * runtime_opt.mcsp_mom_coeff * shear_u * cos_struct / ztodt; + mcsp_tend_v(k) = -1 * runtime_opt.mcsp_mom_coeff * shear_v * cos_struct / ztodt; + } // scale the vertical structure by the ZM heating/drying tendencies if (do_mcsp_t) mcsp_tend_s(k) = zm_avg_tend_s * mcsp_tend_s(k); @@ -155,9 +180,9 @@ void Functions::zm_conv_mcsp_tend( if (do_mcsp_q) avg_q += mcsp_tend_q(k) * state_pdel(k) / pdel_sum; // integrate the change in kinetic energy (KE) for energy fixer - if (do_mcsp_u || do_mcsp_v) { + if (do_mcsp_mom) { const Real tend_k = (2 * mcsp_tend_u(k) * ztodt * state_u(k) + mcsp_tend_u(k) * mcsp_tend_u(k) * ztodt * ztodt - + 2 * mcsp_tend_v(k) * ztodt * state_v(k) + mcsp_tend_v(k) * mcsp_tend_v(k) * ztodt * ztodt) / 2 / ztodt; + +2 * mcsp_tend_v(k) * ztodt * state_v(k) + mcsp_tend_v(k) * mcsp_tend_v(k) * ztodt * ztodt) / 2 / ztodt; avg_k += tend_k * state_pdel(k) / pdel_sum; } }, @@ -184,15 +209,15 @@ void Functions::zm_conv_mcsp_tend( // make sure kinetic energy correction is added to DSE tendency // to conserve total energy whenever momentum tendencies are calculated - if (do_mcsp_u || do_mcsp_v) { + if (do_mcsp_mom) { mcsp_ds_out(k) = mcsp_ds_out(k) - mcsp_avg_tend_k; } // update output tendencies - if (do_mcsp_t) ptend_s(k) = ptend_s(k) + mcsp_ds_out(k); - if (do_mcsp_q) ptend_q(k) = ptend_q(k) + mcsp_dq_out(k); - if (do_mcsp_u) ptend_u(k) = ptend_u(k) + mcsp_du_out(k); - if (do_mcsp_v) ptend_v(k) = ptend_v(k) + mcsp_dv_out(k); + if (do_mcsp_t) ptend_s(k) = ptend_s(k) + mcsp_ds_out(k); + if (do_mcsp_q) ptend_q(k) = ptend_q(k) + mcsp_dq_out(k); + if (do_mcsp_mom) ptend_u(k) = ptend_u(k) + mcsp_du_out(k); + if (do_mcsp_mom) ptend_v(k) = ptend_v(k) + mcsp_dv_out(k); // adjust units for diagnostic outputs if (do_mcsp_t) mcsp_ds_out(k) = mcsp_ds_out(k) / PC::Cpair.value; diff --git a/components/eamxx/src/physics/zm/impl/zm_opts_impl.hpp b/components/eamxx/src/physics/zm/impl/zm_opts_impl.hpp index e940a4caed54..ad8c3f2d9d83 100644 --- a/components/eamxx/src/physics/zm/impl/zm_opts_impl.hpp +++ b/components/eamxx/src/physics/zm/impl/zm_opts_impl.hpp @@ -41,8 +41,13 @@ void Functions::zm_opts_init() s_zm_opts.mcsp_enabled = true; s_zm_opts.mcsp_t_coeff = ZMC::MCSP_t_coeff_default; s_zm_opts.mcsp_q_coeff = ZMC::MCSP_q_coeff_default; - s_zm_opts.mcsp_u_coeff = ZMC::MCSP_u_coeff_default; - s_zm_opts.mcsp_v_coeff = ZMC::MCSP_v_coeff_default; + // Use a negative momentum coefficient so the unit tests exercise the + // do_mcsp_mom path AND confirm the guards admit a negative coefficient + // (down-gradient forcing). Must match the Fortran zm_param_set_for_testing + // setup. (This routine is only used by the standalone test harness; the real + // eamxx run overwrites s_zm_opts via load_runtime_options from XML.) + s_zm_opts.mcsp_mom_coeff = -0.05; + s_zm_opts.mcsp_use_full_shear = true; // // set up table values of saturation vapor pressure diff --git a/components/eamxx/src/physics/zm/tests/infra/zm_c2f_bridge.f90 b/components/eamxx/src/physics/zm/tests/infra/zm_c2f_bridge.f90 index fede6c75eed0..1c28b4d1c2f6 100644 --- a/components/eamxx/src/physics/zm/tests/infra/zm_c2f_bridge.f90 +++ b/components/eamxx/src/physics/zm/tests/infra/zm_c2f_bridge.f90 @@ -171,14 +171,14 @@ subroutine compute_cape_from_parcel_bridge_f(pcols, ncol, pver, pverp, num_cin, call compute_cape_from_parcel(pcols, ncol, pver, pverp, num_cin, num_msg, temperature, tv, sp_humidity, pint, msemax_klev, lcl_pmid, lcl_klev, zm_const, zm_param, parcel_qsat, parcel_temp, parcel_vtemp, eql_klev, cape) end subroutine compute_cape_from_parcel_bridge_f -subroutine zm_conv_mcsp_calculate_shear_bridge_f(pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear) bind(C) +subroutine zm_conv_mcsp_calculate_shear_bridge_f(pcols, ncol, pver, state_pmid, state_u, state_v, shear_u, shear_v) bind(C) use zm_conv_mcsp, only : zm_conv_mcsp_calculate_shear integer(kind=c_int) , value, intent(in) :: pcols, ncol, pver real(kind=c_real) , intent(in), dimension(pcols, pver) :: state_pmid, state_u, state_v - real(kind=c_real) , intent(out), dimension(pcols) :: mcsp_shear + real(kind=c_real) , intent(out), dimension(pcols) :: shear_u, shear_v - call zm_conv_mcsp_calculate_shear(pcols, ncol, pver, state_pmid, state_u, state_v, mcsp_shear) + call zm_conv_mcsp_calculate_shear(pcols, ncol, pver, state_pmid, state_u, state_v, shear_u, shear_v) end subroutine zm_conv_mcsp_calculate_shear_bridge_f subroutine zm_conv_mcsp_tend_bridge_f(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_ds_out, mcsp_dq_out, mcsp_du_out, mcsp_dv_out, mcsp_freq, mcsp_shear, zm_depth) bind(C) diff --git a/components/eamxx/src/physics/zm/tests/infra/zm_test_data.cpp b/components/eamxx/src/physics/zm/tests/infra/zm_test_data.cpp index e064e0d410df..1f813bd51832 100644 --- a/components/eamxx/src/physics/zm/tests/infra/zm_test_data.cpp +++ b/components/eamxx/src/physics/zm/tests/infra/zm_test_data.cpp @@ -54,7 +54,7 @@ void compute_dilute_parcel_bridge_f(Int pcols, Int ncol, Int pver, Int num_msg, void compute_cape_from_parcel_bridge_f(Int pcols, Int ncol, Int pver, Int pverp, Int num_cin, Int num_msg, Real* temperature, Real* tv, Real* sp_humidity, Real* pint, Int* msemax_klev, Real* lcl_pmid, Int* lcl_klev, Real* parcel_qsat, Real* parcel_temp, Real* parcel_vtemp, Int* eql_klev, Real* cape); -void zm_conv_mcsp_calculate_shear_bridge_f(Int pcols, Int ncol, Int pver, Real* state_pmid, Real* state_u, Real* state_v, Real* mcsp_shear); +void zm_conv_mcsp_calculate_shear_bridge_f(Int pcols, Int ncol, Int pver, Real* state_pmid, Real* state_u, Real* state_v, Real* shear_u, Real* shear_v); void zm_conv_mcsp_tend_bridge_f(Int pcols, Int ncol, Int pver, Int pverp, Real ztodt, Int* jctop, Real* state_pmid, Real* state_pint, Real* state_pdel, Real* state_s, Real* state_q, Real* state_u, Real* state_v, Real* ptend_zm_s, Real* ptend_zm_q, Real* ptend_s, Real* ptend_q, Real* ptend_u, Real* ptend_v, Real* mcsp_ds_out, Real* mcsp_dq_out, Real* mcsp_du_out, Real* mcsp_dv_out, Real* mcsp_freq, Real* mcsp_shear, Real* zm_depth); @@ -864,7 +864,7 @@ void zm_conv_mcsp_calculate_shear_f(ZmConvMcspCalculateShearData& d) { d.transition(); zm_opts_init_f(); - zm_conv_mcsp_calculate_shear_bridge_f(d.pcols, d.ncol, d.pver, d.state_pmid, d.state_u, d.state_v, d.mcsp_shear); + zm_conv_mcsp_calculate_shear_bridge_f(d.pcols, d.ncol, d.pver, d.state_pmid, d.state_u, d.state_v, d.shear_u, d.shear_v); zm_opts_finalize_f(); d.transition(); } @@ -874,14 +874,15 @@ void zm_conv_mcsp_calculate_shear(ZmConvMcspCalculateShearData& d) zm_opts_init(); // create device views and copy - std::vector vec1dr_in(1); - ekat::host_to_device({d.mcsp_shear}, d.pcols, vec1dr_in); + std::vector vec1dr_in(2); + ekat::host_to_device({d.shear_u, d.shear_v}, d.pcols, vec1dr_in); std::vector vec2dr_in(3); ekat::host_to_device({d.state_pmid, d.state_u, d.state_v}, d.pcols, d.pver, vec2dr_in); view1dr_d - mcsp_shear_d(vec1dr_in[0]); + shear_u_d(vec1dr_in[0]), + shear_v_d(vec1dr_in[1]); view2dr_d state_pmid_d(vec2dr_in[0]), @@ -900,18 +901,21 @@ void zm_conv_mcsp_calculate_shear(ZmConvMcspCalculateShearData& d) // after this. const auto state_pmid_c = ekat::subview(state_pmid_d, i); const auto state_u_c = ekat::subview(state_u_d, i); + const auto state_v_c = ekat::subview(state_v_d, i); ZMF::zm_conv_mcsp_calculate_shear( team, pver, state_pmid_c, state_u_c, - mcsp_shear_d(i)); + state_v_c, + shear_u_d(i), + shear_v_d(i)); }); // Now get arrays - std::vector vec1dr_out = {mcsp_shear_d}; - ekat::device_to_host({d.mcsp_shear}, d.pcols, vec1dr_out); + std::vector vec1dr_out = {shear_u_d, shear_v_d}; + ekat::device_to_host({d.shear_u, d.shear_v}, d.pcols, vec1dr_out); zm_finalize_cxx(); } diff --git a/components/eamxx/src/physics/zm/tests/infra/zm_test_data.hpp b/components/eamxx/src/physics/zm/tests/infra/zm_test_data.hpp index 73ec02123945..36e2562833d0 100644 --- a/components/eamxx/src/physics/zm/tests/infra/zm_test_data.hpp +++ b/components/eamxx/src/physics/zm/tests/infra/zm_test_data.hpp @@ -345,7 +345,7 @@ struct ZmConvMcspCalculateShearData : public PhysicsTestData { Real *state_pmid, *state_u, *state_v; // Outputs - Real *mcsp_shear; + Real *shear_u, *shear_v; ZmConvMcspCalculateShearData(Int pcols_, Int ncol_, Int pver_) : PhysicsTestData({ @@ -354,7 +354,7 @@ struct ZmConvMcspCalculateShearData : public PhysicsTestData { }, { {&state_pmid, &state_u, &state_v}, - {&mcsp_shear} + {&shear_u, &shear_v} }), pcols(pcols_), ncol(ncol_), pver(pver_) {} diff --git a/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_calculate_shear_tests.cpp b/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_calculate_shear_tests.cpp index 5bc3adbb9412..14ac2d5ca783 100644 --- a/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_calculate_shear_tests.cpp +++ b/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_calculate_shear_tests.cpp @@ -65,9 +65,13 @@ struct UnitWrap::UnitTest::TestZmConvMcspCalculateShear : public UnitWrap::Un for (Int i = 0; i < num_runs; ++i) { ZmConvMcspCalculateShearData& d_baseline = baseline_data[i]; ZmConvMcspCalculateShearData& d_test = test_data[i]; - REQUIRE(d_baseline.total(d_baseline.mcsp_shear) == d_test.total(d_test.mcsp_shear)); - for (Int k = 0; k < d_baseline.total(d_baseline.mcsp_shear); ++k) { - REQUIRE(d_baseline.mcsp_shear[k] == d_test.mcsp_shear[k]); + REQUIRE(d_baseline.total(d_baseline.shear_u) == d_test.total(d_test.shear_u)); + for (Int k = 0; k < d_baseline.total(d_baseline.shear_u); ++k) { + REQUIRE(d_baseline.shear_u[k] == d_test.shear_u[k]); + } + REQUIRE(d_baseline.total(d_baseline.shear_v) == d_test.total(d_test.shear_v)); + for (Int k = 0; k < d_baseline.total(d_baseline.shear_v); ++k) { + REQUIRE(d_baseline.shear_v[k] == d_test.shear_v[k]); } } } diff --git a/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_tend_tests.cpp b/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_tend_tests.cpp index 6c370eb29033..80f0d4644880 100644 --- a/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_tend_tests.cpp +++ b/components/eamxx/src/physics/zm/tests/zm_conv_mcsp_tend_tests.cpp @@ -79,14 +79,17 @@ struct UnitWrap::UnitTest::TestZmConvMcspTend : public UnitWrap::UnitTest: REQUIRE(d_baseline.total(d_baseline.ptend_s) == d_test.total(d_test.mcsp_du_out)); REQUIRE(d_baseline.total(d_baseline.ptend_s) == d_test.total(d_test.mcsp_dv_out)); for (Int k = 0; k < d_baseline.total(d_baseline.ptend_s); ++k) { + // Momentum (u/v) fields use Approx() because the momentum tendency is + // built from cos(), which is std::cos in C++ but bfb_cos in Fortran; + // the two are not bit-for-bit. DSE (s) fields likewise use Approx(). REQUIRE(d_baseline.ptend_s[k] == Approx(d_test.ptend_s[k]).margin(margin)); REQUIRE(d_baseline.ptend_q[k] == d_test.ptend_q[k]); - REQUIRE(d_baseline.ptend_u[k] == d_test.ptend_u[k]); - REQUIRE(d_baseline.ptend_v[k] == d_test.ptend_v[k]); + REQUIRE(d_baseline.ptend_u[k] == Approx(d_test.ptend_u[k]).margin(margin)); + REQUIRE(d_baseline.ptend_v[k] == Approx(d_test.ptend_v[k]).margin(margin)); REQUIRE(d_baseline.mcsp_ds_out[k] == Approx(d_test.mcsp_ds_out[k]).margin(margin)); REQUIRE(d_baseline.mcsp_dq_out[k] == d_test.mcsp_dq_out[k]); - REQUIRE(d_baseline.mcsp_du_out[k] == d_test.mcsp_du_out[k]); - REQUIRE(d_baseline.mcsp_dv_out[k] == d_test.mcsp_dv_out[k]); + REQUIRE(d_baseline.mcsp_du_out[k] == Approx(d_test.mcsp_du_out[k]).margin(margin)); + REQUIRE(d_baseline.mcsp_dv_out[k] == Approx(d_test.mcsp_dv_out[k]).margin(margin)); } REQUIRE(d_baseline.total(d_baseline.mcsp_freq) == d_test.total(d_test.mcsp_freq)); REQUIRE(d_baseline.total(d_baseline.mcsp_freq) == d_test.total(d_test.mcsp_shear)); diff --git a/components/eamxx/src/physics/zm/zm_functions.hpp b/components/eamxx/src/physics/zm/zm_functions.hpp index c4a3d748243c..74a791a57dbe 100644 --- a/components/eamxx/src/physics/zm/zm_functions.hpp +++ b/components/eamxx/src/physics/zm/zm_functions.hpp @@ -122,8 +122,7 @@ struct Functions { static inline constexpr Real MCSP_shear_max = 200.0; // max shear value for MCSP to be active static inline constexpr Real MCSP_t_coeff_default = 0.3; // default MCSP temperature coefficient static inline constexpr Real MCSP_q_coeff_default = 0.0; // default MCSP sp. humidity coefficient - static inline constexpr Real MCSP_u_coeff_default = 0.0; // default MCSP U-wind coefficient - static inline constexpr Real MCSP_v_coeff_default = 0.0; // default MCSP V-wind coefficient + static inline constexpr Real MCSP_mom_coeff_default = 0.0; // default MCSP momentum coefficient (fraction of shear) // Default values for ZmRuntimeOpt fields static inline constexpr Real alfa = 0.14; // default downdraft proportionality factor static inline constexpr Real ke = 2.5E-6; // default evaporation efficiency @@ -175,8 +174,8 @@ struct Functions { mcsp_enabled = params.get("mcsp_enabled", true); mcsp_t_coeff = params.get("mcsp_t_coeff", ZMC::MCSP_t_coeff_default); mcsp_q_coeff = params.get("mcsp_q_coeff", ZMC::MCSP_q_coeff_default); - mcsp_u_coeff = params.get("mcsp_u_coeff", ZMC::MCSP_u_coeff_default); - mcsp_v_coeff = params.get("mcsp_v_coeff", ZMC::MCSP_v_coeff_default); + mcsp_mom_coeff = params.get("mcsp_mom_coeff", ZMC::MCSP_mom_coeff_default); + mcsp_use_full_shear = params.get("mcsp_use_full_shear", false); // determine SVP table size (add two to make the table slightly too big, just in case) plenest = static_cast(ZMC::tmax-ZMC::tmin) + 3; @@ -251,8 +250,8 @@ struct Functions { os << indent << "mcsp_enabled : " << mcsp_enabled << "\n"; os << indent << "mcsp_t_coeff : " << mcsp_t_coeff << "\n"; os << indent << "mcsp_q_coeff : " << mcsp_q_coeff << "\n"; - os << indent << "mcsp_u_coeff : " << mcsp_u_coeff << "\n"; - os << indent << "mcsp_v_coeff : " << mcsp_v_coeff << "\n"; + os << indent << "mcsp_mom_coeff : " << mcsp_mom_coeff << "\n"; + os << indent << "mcsp_use_full_shear : " << mcsp_use_full_shear << "\n"; os << std::endl; os.flags(saved_flags); } @@ -283,8 +282,8 @@ struct Functions { bool mcsp_enabled; // flag for mesoscale coherent structure parameterization (MSCP) Real mcsp_t_coeff; // MCSP coefficient for temperature tendencies Real mcsp_q_coeff; // MCSP coefficient for specific humidity tendencies - Real mcsp_u_coeff; // MCSP coefficient for zonal momentum tendencies - Real mcsp_v_coeff; // MCSP coefficient for meridional momentum tendencies + Real mcsp_mom_coeff; // MCSP momentum coefficient (dimensionless fraction of the shear) + bool mcsp_use_full_shear; // use full (u,v) shear vector instead of zonal-only shear // saturation vapor pressure (svp) table Int plenest; // saturation vapor pressure table size view_1d estbl; // saturation vapor pressure table values @@ -751,8 +750,10 @@ struct Functions { const Int& pver, // number of mid-point vertical levels const uview_1d& state_pmid, // physics state mid-point pressure const uview_1d& state_u, // physics state u momentum + const uview_1d& state_v, // physics state v momentum // Outputs - Real& mcsp_shear); + Real& shear_u, // zonal component of storm-relative shear + Real& shear_v); // meridional component of storm-relative shear KOKKOS_FUNCTION static void zm_conv_mcsp_tend(