Add a namelist option of newshear to control the use of the updated

vertical wind shear impacts on GFS scale-aware TKE-based EDMF PBL and
saSAS CP. (From @WeiguoWang-NOAA)

physics/CONV/SAMF/samfdeepcnv.f

@@ -84,7 +84,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &    clam,c0s,c1,betal,betas,evef,pgcon,asolfac,                   &
      &    do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep,      &
      &    rainevap,sigmain,sigmaout,betadcu,betamcu,betascu,            &
-     &    maxMF, do_mynnedmf,errmsg,errflg)
+     &    maxMF, do_mynnedmf,newshear,errmsg,errflg)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -117,6 +117,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &   cnvw(:,:),  cnvc(:,:)
 
       integer, intent(out) :: kbot(:), ktop(:)
+      integer, intent(in)  ::  newshear
       real(kind=kind_phys), intent(out) :: cldwrk(:),                   &
      &   rn(:),                                                         &
      &   ud_mf(:,:),dd_mf(:,:), dt_mf(:,:)
@@ -242,8 +243,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       parameter(cinpcrmx=180.,cinpcrmn=120.)
 !     parameter(cinacrmx=-120.,cinacrmn=-120.)
       parameter(cinacrmx=-120.,cinacrmn=-80.)
-      parameter(bb1=4.0,bb2=0.8,csmf=0.2)
-      parameter(tkcrt=2.,cmxfac=15.)
+!      parameter(bb1=4.0,bb2=0.8,csmf=0.2)
+!      parameter(tkcrt=2.,cmxfac=15.)
+      parameter(bb1=4.0,bb2=0.8)
+      parameter(cmxfac=15.)
       parameter(betaw=.03)
 
 !
@@ -307,6 +310,14 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       real(kind=kind_phys) tf, tcr, tcrf
       parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
 
+       if (newshear == 0) then
+          csmf=0.0
+          tkcrt=99999.0
+       else
+          csmf=0.2
+          tkcrt=2.0
+       endif
+
       ! Initialize CCPP error handling variables
       errmsg = ''
       errflg = 0

physics/CONV/SAMF/samfdeepcnv.meta

@@ -698,6 +698,13 @@
   type = real
   kind = kind_phys
   intent = out
+[newshear]
+  standard_name = choice_of_using_updated_shear_effects_in_TKE_EDMF_and_SAS
+  long_name = choice of using updated shear effects in TKE EDMF and SAS
+  units = none
+  dimensions = ()
+  type = integer
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/CONV/SAMF/samfshalcnv.f

@@ -57,7 +57,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
      &     rn,kbot,ktop,kcnv,islimsk,garea,                             &
      &     dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,                       &
      &     clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal,                & 
-     &     sigmain,sigmaout,betadcu,betamcu,betascu,errmsg,errflg)
+     &     sigmain,sigmaout,betadcu,betamcu,betascu,newshear,           &
+     &     errmsg,errflg)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -92,6 +93,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
      &     restart,progsigma
       character(len=*), intent(out) :: errmsg
       integer,          intent(out) :: errflg
+      integer,          intent(in)  :: newshear
 !
 !  local variables
       integer              i,j,indx, k, kk, km1, n
@@ -190,8 +192,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 !  shevf is an enhancing evaporation factor for shallow convection
       parameter(cinacrmx=-120.,shevf=2.0)
       parameter(dtmax=10800.,dtmin=600.)
-      parameter(bb1=4.0,bb2=0.8,csmf=0.2)
-      parameter(tkcrt=2.,cmxfac=10.)
+!      parameter(bb1=4.0,bb2=0.8,csmf=0.2)
+      parameter(bb1=4.0,bb2=0.8)
+      parameter(cmxfac=10.)
+!      parameter(tkcrt=2.,cmxfac=10.)
 !      parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
       parameter(betaw=.03,dxcrtc0=9.e3)
       parameter(h1=0.33333333)
@@ -246,6 +250,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 
 
 c-----------------------------------------------------------------------
+      if (newshear == 0) then
+        csmf=0.0
+        tkcrt=99999.0
+      else
+        csmf=0.2
+        tkcrt=2.0
+      endif
 !
 ! Initialize CCPP error handling variables
       errmsg = ''

physics/CONV/SAMF/samfshalcnv.meta

@@ -505,6 +505,13 @@
   dimensions = ()
   type = real
   intent = in
+[newshear]
+  standard_name = choice_of_using_updated_shear_effects_in_TKE_EDMF_and_SAS
+  long_name = choice of using updated shear effects in TKE EDMF and SAS
+  units = none
+  dimensions = ()
+  type = integer
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/PBL/SATMEDMF/mfscuq.f

@@ -13,7 +13,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
      &   cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,
      &   thlx,thvx,thlvx,gdx,thetae,
      &   krad,mrad,radmin,buo,wush,tkemean,vez0fun,xmfd,
-     &   tcdo,qcdo,ucdo,vcdo,xlamdeq,a1)
+     &   tcdo,qcdo,ucdo,vcdo,xlamdeq,a1,newshear)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -26,6 +26,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
 !
       integer            im, ix,  km, kmscu, ntcw, ntrac1
 !    &,                  me
+      integer            newshear
       integer   krad(im), mrad(im)
 !
       logical cnvflg(im)
@@ -83,7 +84,8 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
       parameter(gocp=g/cp)
       parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
       parameter(ce0=0.4,cm=1.0,cq=1.0,pgcon=0.55)
-      parameter(tkcrt=2.,cmxfac=5.)
+!      parameter(tkcrt=2.,cmxfac=5.)
+      parameter(cmxfac=5.)
       parameter(qmin=1.e-8,qlmin=1.e-12)
       parameter(b1=0.45,f1=0.15)
       parameter(a2=0.5)
@@ -93,6 +95,12 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt,
 !
 !************************************************************************
 !!
+      if (newshear == 1 ) then
+          tkcrt=2.0
+      else
+          tkcrt=99999.0
+      endif
+
       totflg = .true.
       do i=1,im
         totflg = totflg .and. (.not. cnvflg(i))

physics/PBL/SATMEDMF/satmedmfvdifq.F

@@ -82,7 +82,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
      &     prsi,del,prsl,prslk,phii,phil,delt,                          &
      &     dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,                &
      &     kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr,          &
-     &     rlmx,elmx,sfc_rlm,tc_pbl,                                    &
+     &     rlmx,elmx,sfc_rlm,tc_pbl,newshear,                           &
      &     ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind,       &
      &     index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d,       &
      &     errmsg,errflg)
@@ -97,6 +97,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
      &                        ntke, ntqv
       integer, intent(in)  :: sfc_rlm
       integer, intent(in)  :: tc_pbl
+      integer, intent(in)  :: newshear
       integer, intent(in)  :: kinver(:)
       integer, intent(out) :: kpbl(:)
       logical, intent(in)  :: gen_tend,ldiag3d
@@ -279,9 +280,16 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
       parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0)
       parameter(vc0=1.0,zc0=1.0)
       parameter(ck1=0.15,ch1=0.15)
-      parameter(cs0=0.4,csmf=0.5)
+!      parameter(cs0=0.4,csmf=0.5)
       parameter(rchck=1.5,ndt=20)
 
+      cs0=0.4
+      csmf=0.5
+      if(newshear == 0) then
+       cs0=0.2
+       csmf=0.0
+      endif
+
       if (tc_pbl == 0) then
         ck0 = 0.4
         ch0 = 0.4
@@ -457,6 +465,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
         ptem = (z0(i) - z0lo) / (z0up - z0lo)
         ptem = min(max(ptem, 0.), 1.)
         vez0fun(i) = (1. + vc0 * tem1) * (1. + zc0 * ptem)
+        if ( newshear == 0) vez0fun(i)=1.0
       enddo
 !
 !>  - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae),
@@ -1014,14 +1023,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
       call mfpbltq(im,im,km,kmpbl,ntcw,ntrac1,dt2,
      &    pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
      &    gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf,
-     &    tcko,qcko,ucko,vcko,xlamue,bl_upfr)
+     &    tcko,qcko,ucko,vcko,xlamue,bl_upfr,newshear)
 !>  - Call mfscuq(), which is a new mass-flux parameterization for 
 !!    stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq
       call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2,
      &    scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
      &    thlx,thvx,thlvx,gdx,thetae,
      &    krad,mrad,radmin,buod,wush,tkemean,vez0fun,xmfd,
-     &    tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr)
+     &    tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr,newshear)
 
       if (tc_pbl == 1) then
 !>  - unify mass fluxes with Cu

physics/PBL/SATMEDMF/satmedmfvdifq.meta

@@ -611,6 +611,13 @@
   dimensions = ()
   type = integer
   intent = in
+[newshear]
+  standard_name = choice_of_using_updated_shear_effects_in_TKE_EDMF_and_SAS
+  long_name = choice of using updated shear effects in TKE EDMF and SAS
+  units = none
+  dimensions = ()
+  type = integer
+  intent = in
 [ntqv]
   standard_name = index_of_specific_humidity_in_tracer_concentration_array
   long_name = tracer index for water vapor (specific humidity)

physics/PBL/mfpbltq.f

@@ -12,7 +12,7 @@ module mfpbltq_mod
       subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
      &   cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
      &   gdx,hpbl,kpbl,vpert,buo,wush,tkemean,vez0fun,xmf,
-     &   tcko,qcko,ucko,vcko,xlamueq,a1)
+     &   tcko,qcko,ucko,vcko,xlamueq,a1,newshear)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -25,6 +25,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
 !
       integer              im, ix, km, kmpbl, ntcw, ntrac1
 !    &,                    me
+      integer              newshear
       integer              kpbl(im)
       logical              cnvflg(im)
       real(kind=kind_phys) delt
@@ -74,13 +75,21 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt,
       parameter(g=grav)
       parameter(gocp=g/cp)
       parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
-      parameter(ce0=0.4,cm=1.0,cq=1.0,tkcrt=2.,cmxfac=5.)
+!      parameter(ce0=0.4,cm=1.0,cq=1.0,tkcrt=2.,cmxfac=5.)
+      parameter(ce0=0.4,cm=1.0,cq=1.0,cmxfac=5.)
       parameter(qmin=1.e-8,qlmin=1.e-12)
       parameter(alp=1.5,vpertmax=3.0,pgcon=0.55)
       parameter(b1=0.5,f1=0.15)
 !
 !************************************************************************
 !!
+      if (newshear == 1 ) then
+        tkcrt=2.0
+      else
+        tkcrt=99999.0
+      endif
+
+
       totflg = .true.
       do i=1,im
         totflg = totflg .and. (.not. cnvflg(i))