ufs-dev PR#363 PR#339 PR#364

physics/CONV/SAMF/samfdeepcnv.f

@@ -220,10 +220,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 !
 !  parameters for prognostic sigma closure                                                                                                                                                      
       real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
-     &     omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km),
-     &     tentr(im,km)
-      real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
-      parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
+     &     omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km)
+      real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins,
+     &     wc_min, wc_eff 
       logical flag_shallow, flag_mid
 c  physical parameters
 !     parameter(grav=grav,asolfac=0.958)
@@ -347,6 +346,21 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 !************************************************************************
 !
 !
+!   - Initialize parameters related to prognostic closure
+      if (progsigma) then
+         if (progomega) then
+            sigmind  = 0.03
+            sigmins  = 0.03
+            sigminm = 0.03
+            wc_min = 0.2
+         else
+            sigmind  = 0.01
+            sigmins  = 0.03
+            sigminm = 0.03
+            wc_min = 0.2
+         endif
+      endif
+
       km1 = km - 1
 !>  - Initialize column-integrated and other single-value-per-column variable arrays.
 c
@@ -1137,7 +1151,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
         do i=1,im
           if(cnvflg(i) .and.
      &       (k > kbcon(i) .and. k < kmax(i))) then
-              tentr(i,k)=xlamue(i,k)*fent1(i,k)
               tem = cxlamet(i) * frh(i,k) * fent2(i,k)
               xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
               tem1 = cxlamdt(i) * frh(i,k)
@@ -1788,11 +1801,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       if (progomega) then
          call progomega_calc(first_time_step,restart,im,km,
      &        kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
-     &        grav,buo,drag,wush,tentr,bb1,bb2)
+     &        grav,buo,drag,wush,bb1,bb2)
          do k = 1, km
             do i = 1, im
                if (cnvflg(i)) then
-                  if(k > kbcon1(i) .and. k < ktcon(i)) then
+                  if(k >= kbcon1(i) .and. k < ktcon(i)) then
                      omega_u(i,k)=omegaout(i,k)
                      omega_u(i,k)=MAX(omega_u(i,k),-80.)
 !     Convert to m/s for use in convective time-scale:
@@ -1878,7 +1891,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
          do k = 2, km1
             do i = 1, im
                if (cnvflg(i)) then
-                  if(k > kbcon1(i) .and. k < ktcon(i)) then
+                  if(k >= kbcon1(i) .and. k < ktcon(i)) then
                      dp = 1000. * del(i,k)
                      tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
                      omegac(i) = omegac(i) + tem * dp
@@ -2907,27 +2920,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 !  compute convective turn-over time
 !
 !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
-      if(hwrf_samfdeep) then
-       do i= 1, im
-        if(cnvflg(i)) then
-          tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
-          dtconv(i) = tem / wc(i)
-          dtconv(i) = max(dtconv(i),dtmin)
-          dtconv(i) = min(dtconv(i),dtmax)
-        endif
-       enddo
-      else
-       do i= 1, im
-        if(cnvflg(i)) then
-          tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
-          dtconv(i) = tem / wc(i)
-          tfac = 1. + gdx(i) / 75000.
-          dtconv(i) = tfac * dtconv(i)
-          dtconv(i) = max(dtconv(i),dtmin)
-          dtconv(i) = min(dtconv(i),dtmax)
-        endif
-       enddo
-      endif
+      do i = 1, im
+         if (cnvflg(i)) then
+            tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
+            if (progomega) then
+               wc_eff = max(wc(i), wc_min)
+               dtconv(i) = tem / wc_eff
+            else
+               dtconv(i) = tem / wc(i)
+            endif
+         !grid spacing scaling (disabled for HWRF SAMF deep)
+            if (.not. hwrf_samfdeep) then
+               tfac = 1. + gdx(i) / 75000.
+               dtconv(i) = tfac * dtconv(i)
+            endif  
+         !bounds
+            dtconv(i) = max(dtconv(i), dtmin)
+            dtconv(i) = min(dtconv(i), dtmax)
+         endif
+      enddo
 !
 !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
       do i= 1, im

physics/CONV/SAMF/samfshalcnv.f

@@ -171,7 +171,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
      &                     omegac(im),zeta(im,km),dbyo1(im,km),
      &                     sigmab(im),qadv(im,km)
       real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
-     &                     sigminm
+     &                     sigminm,wc_min,wc_eff
       logical flag_shallow,flag_mid
 c  physical parameters
 !     parameter(g=grav,asolfac=0.89)
@@ -205,8 +205,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 !      parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
       parameter(betaw=.03,dxcrtc0=9.e3)
       parameter(h1=0.33333333)
-!  progsigma
-      parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01)
 c  local variables and arrays
       real(kind=kind_phys) pfld(im,km),    to(im,km),     qo(im,km),
      &                     uo(im,km),      vo(im,km),     qeso(im,km),
@@ -254,7 +252,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       real(kind=kind_phys) tf, tcr, tcrf
       parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
 
-
 c-----------------------------------------------------------------------
 !
 ! Initialize CCPP error handling variables
@@ -276,8 +273,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 
       if (progsigma) then
          dxcrt=10.e3
+         dxcrtas=500.e3
       else
          dxcrt=15.e3
+         dxcrtas=500.e3
       endif
 
 c-----------------------------------------------------------------------
@@ -295,7 +294,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       prsl = prslp * 0.001
       del  = delp  * 0.001
 !************************************************************************
-!
+!   - Initialize parameters related to prognostic closure
+      if (progsigma) then
+ 	 if (progomega) then
+            sigmind  = 0.03
+            sigmins  = 0.03
+            sigminm = 0.03
+            wc_min = 0.2
+         else
+   	    sigmind  = 0.01
+            sigmins  = 0.03
+            sigminm = 0.03
+            wc_min = 0.2
+ 	 endif
+      endif
+!     
       km1 = km - 1
 c
 c  initialize arrays
@@ -1520,11 +1533,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       if (progomega) then
          call progomega_calc(first_time_step,restart,im,km,
      &        kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
-     &        grav,buo,drag,wush,xlamue,bb1,bb2)
+     &        grav,buo,drag,wush,bb1,bb2)
          do k = 1, km
             do i = 1, im
                if (cnvflg(i)) then
-                  if(k > kbcon1(i) .and. k < ktcon(i)) then
+                  if(k >= kbcon1(i) .and. k < ktcon(i)) then
                      omega_u(i,k)=omegaout(i,k)
                      omega_u(i,k)=MAX(omega_u(i,k),-80.)
 !     Convert to m/s for use in convective time-scale:
@@ -1611,7 +1624,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
          do k = 2, km1
             do i = 1, im
                if (cnvflg(i)) then
-                  if(k > kbcon1(i) .and. k < ktcon(i)) then
+                  if(k >= kbcon1(i) .and. k < ktcon(i)) then
                      dp = 1000. * del(i,k)
                      tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
                      omegac(i) = omegac(i) + tem * dp
@@ -1636,7 +1649,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
          do k = 2, km1
             do i = 1, im
                if (cnvflg(i)) then
-                  if(k > kbcon1(i) .and. k < ktcon(i)) then
+                  if(k >= kbcon1(i) .and. k < ktcon(i)) then
                      if(omega_u(i,k) .ne. 0.)then
                         zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
                      else
@@ -1955,21 +1968,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 !  compute convective turn-over time
 !
 !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
-      do i= 1, im
-        if(cnvflg(i)) then
-          tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
-          dtconv(i) = tem / wc(i)
-          if (.not.hwrf_samfshal) then
-            tfac = 1. + gdx(i) / 75000.
-            dtconv(i) = tfac * dtconv(i)
-          endif
-          dtconv(i) = max(dtconv(i),dtmin)
-          dtconv(i) = max(dtconv(i),dt2)
-          dtconv(i) = min(dtconv(i),dtmax)
-        endif
+      do i = 1, im
+         if (cnvflg(i)) then
+            tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
+            if (progomega) then
+               wc_eff = max(wc(i), wc_min)
+               dtconv(i) = tem / wc_eff
+            else
+               dtconv(i) = tem / wc(i)
+            endif
+! - grid spacing scaling (disabled for HWRF shallow option)
+            if (.not. hwrf_samfshal) then
+               tfac = 1. + gdx(i) / 75000.
+               dtconv(i) = tfac * dtconv(i)
+            endif
+! - limits
+            dtconv(i) = max(dtconv(i), dtmin)
+            dtconv(i) = max(dtconv(i),dt2)
+            dtconv(i) = min(dtconv(i), dtmax)
+         endif
       enddo
-!
-!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
+!     
+!     > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
       do i= 1, im
         if(cnvflg(i)) then
           sumx(i) = 0.

physics/CONV/progomega_calc.f90

@@ -20,7 +20,7 @@ module progomega
 !!\section gen_progomega progomega_calc General Algorithm
 
    subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, &
-        zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2)
+        zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2)
 
      use machine,  only : kind_phys
      use funcphys, only : fpvs  
@@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
      integer, intent(in)  :: kbcon1(im),ktcon(im)
      real(kind=kind_phys), intent(in)  :: delt,grav,bb1,bb2
      real(kind=kind_phys), intent(in)  :: omegain(im,km), del(im,km),zi(im,km)
-     real(kind=kind_phys), intent(in)  :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km)
+     real(kind=kind_phys), intent(in)  :: drag(im,km),buo(im,km),wush(im,km)
      real(kind=kind_phys), intent(inout) :: omegaout(im,km)
      logical, intent(in)               :: cnvflg(im),first_time_step,flag_restart
      real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km)
      real(kind=kind_phys) :: RHS(im,km),Kd(im,km)
-     real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3
+     real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3
      integer              :: i,k
 
-     entrn = 0.8E-4 !0.5E-4 !m^-1
-     Kdn   = 0.5E-4 !2.9E-4 !m^-1
-     lbb1  = 0.5 !1.0 
-     lbb2  = 3.2 !3.0
-     lbb3  = 0.5 !0.5
-
+     lbb1  = 1.5
+     lbb2  = 0.6
+     lbb3  = 1.2
 
      !Initialization 2D
      do k = 1,km
@@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
         enddo
      enddo
 
+     do k = 1,km
+        do i = 1,im
+           if(cnvflg(i))then
+              if(omega(i,k) < 1.0E-5) then
+                 omega(i,k) = 0.
+              endif
+           endif
+        enddo
+     enddo
+
      if(first_time_step .and. .not. flag_restart)then
         do k = 1,km
            do i = 1,im
@@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
      do k = 2, km
         do i = 1, im
            if (cnvflg(i)) then
-              if (k > kbcon1(i) .and. k < ktcon(i)) then
-
-                 ! Aerodynamic drag parameter
-                 Kd(i,k) = (tentr(i,k)/entrn)*Kdn
+              if (k >= kbcon1(i) .and. k < ktcon(i)) then
 
                  ! Scale by dp/dz to have equation in Pa/s
                  !(dp/dz > 0)
@@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
                  !termC - Adds buoyancy forcing 
 
                  !Coefficients for the quadratic equation
-                 termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz)))
-                 termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz
+                 termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz))) 
+                 termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz
                  termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) &
                       - delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp
                  !Compute the discriminant

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_common.F90

@@ -14,15 +14,15 @@ module GFS_PBL_generic_common
       subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6,          &
                                           imp_physics_thompson, ltaerosol,mraerosol,   &
                                           imp_physics_mg, ntgl, imp_physics_gfdl, &
-                                          imp_physics_zhao_carr, imp_physics_nssl,&
+                                          imp_physics_nssl,                       &
                                           nssl_hail_on, nssl_ccn_on, nssl_3moment, kk, &
                                           errmsg, errflg)
       implicit none
       !
       integer, intent(in )          :: imp_physics, imp_physics_wsm6,          &
                                        imp_physics_thompson,                   &
                                        imp_physics_mg, ntgl, imp_physics_gfdl, &
-                                       imp_physics_zhao_carr,imp_physics_nssl
+                                       imp_physics_nssl
       logical, intent(in )          :: ltaerosol, mraerosol, nssl_hail_on, nssl_ccn_on, nssl_3moment
       integer, intent(out)          :: kk
       character(len=*), intent(out) :: errmsg
@@ -53,9 +53,6 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6,          &
       elseif (imp_physics == imp_physics_gfdl) then
 ! GFDL MP
         kk = 7
-      elseif (imp_physics == imp_physics_zhao_carr) then
-! Zhao/Carr/Sundqvist
-        kk = 3
       elseif (imp_physics == imp_physics_nssl) then
         IF ( nssl_hail_on ) THEN
           kk = 16

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90

@@ -11,7 +11,7 @@ module GFS_PBL_generic_post
       subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,                                                            &
         ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev,nqrimef,          &
         trans_aero, ntchs, ntchm, ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz,                                            &
-        imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, imp_physics_mg,          &
+        imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_mg,                                 &
         imp_physics_fer_hires, imp_physics_nssl, nssl_ccn_on, ltaerosol, mraerosol, nssl_hail_on, nssl_3moment,                &
         cplflx, cplaqm, cplchm, lssav, flag_for_pbl_generic_tend, ldiag3d, lsidea, hybedmf, do_shoc, satmedmf,                 &
         shinhong, do_ysu, dvdftra, dusfc1, dvsfc1, dtsfc1, dqsfc1, dtf, dudt, dvdt, dtdt, htrsw, htrlw, xmu,                   &
@@ -33,7 +33,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       integer, intent(in) :: ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz
       logical, intent(in) :: trans_aero
       integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6
-      integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires
+      integer, intent(in) :: imp_physics_mg, imp_physics_fer_hires
       integer, intent(in) :: imp_physics_nssl
       logical, intent(in) :: nssl_ccn_on, nssl_hail_on, nssl_3moment
       logical, intent(in) :: ltaerosol, cplflx, cplaqm, cplchm, lssav, ldiag3d, lsidea, use_med_flux, mraerosol
@@ -109,7 +109,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
           call set_aerosol_tracer_index(imp_physics, imp_physics_wsm6,          &
                                         imp_physics_thompson, ltaerosol,mraerosol,   &
                                         imp_physics_mg, ntgl, imp_physics_gfdl, &
-                                        imp_physics_zhao_carr, imp_physics_nssl,&
+                                        imp_physics_nssl,&
                                         nssl_hail_on, nssl_ccn_on, nssl_3moment, kk, &
                                         errmsg, errflg)
           if (errflg /= 0) return
@@ -244,14 +244,6 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
               dqdt(i,k,ntoz) = dvdftra(i,k,7)
             enddo
           enddo
-        elseif (imp_physics == imp_physics_zhao_carr) then
-          do k=1,levs
-            do i=1,im
-              dqdt(i,k,1)    = dvdftra(i,k,1)
-              dqdt(i,k,ntcw) = dvdftra(i,k,2)
-              dqdt(i,k,ntoz) = dvdftra(i,k,3)
-            enddo
-          enddo
         elseif (imp_physics == imp_physics_nssl ) then
   ! nssl
             IF ( nssl_hail_on ) THEN