Merge pull request #1010 from grantfirl/ufs-dev-PR57

UFS-dev PR#57

physics/drag_suite.F90

@@ -342,7 +342,7 @@ subroutine drag_suite_run(                                           &
    real(kind=kind_phys), intent(inout) ::                        &
      &                   dudt(:,:),dvdt(:,:),                &
      &                   dtdt(:,:)
-   real(kind=kind_phys), intent(out) :: rdxzb(:)
+   real(kind=kind_phys), intent(inout) :: rdxzb(:)
    real(kind=kind_phys), intent(in) ::                           &
      &                            u1(:,:),v1(:,:),           &
      &                            t1(:,:),q1(:,:),           &
@@ -605,7 +605,6 @@ subroutine drag_suite_run(                                           &
       else
          xland(i)=2.0
       endif
-      RDXZB(i)  = 0.0
    enddo
 
 !--- calculate scale-aware tapering factors
@@ -818,6 +817,8 @@ subroutine drag_suite_run(                                           &
 
    do i=its,im
 
+      RDXZB(i)  = 0.0
+
       if ( ls_taper(i).GT.1.E-02 ) then
 
 !

physics/drag_suite.meta

@@ -536,7 +536,7 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
-  intent = out
+  intent = inout
 [dx]
   standard_name = characteristic_grid_lengthscale
   long_name = size of the grid cell

physics/progsigma_calc.f90

@@ -13,8 +13,8 @@
 !!\section gen_progsigma progsigma_calc General Algorithm 
       subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
            flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,    &
-           delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout,           &
-           sigmab,errmsg,errflg)
+           delt,qadv,kbcon1,ktcon,cnvflg,sigmain,sigmaout,           &
+           sigmab)
 !                                                           
 !                                                                                                                                             
       use machine,  only : kind_phys
@@ -25,7 +25,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 !     intent in
       integer, intent(in)  :: im,km,kbcon1(im),ktcon(im)
       real(kind=kind_phys), intent(in)  :: hvap,delt
-      real(kind=kind_phys), intent(in)  :: prevsq(im,km), q(im,km),del(im,km),    &
+      real(kind=kind_phys), intent(in)  :: qadv(im,km),del(im,km),    &
            qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km),           &
            omega_u(im,km),zeta(im,km)
       logical, intent(in)  :: flag_init,flag_restart,cnvflg(im),flag_shallow
@@ -34,14 +34,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 !     intent out
       real(kind=kind_phys), intent(out) :: sigmaout(im,km)
       real(kind=kind_phys), intent(out) :: sigmab(im)
-      character(len=*),     intent(out) :: errmsg
-      integer,              intent(out) :: errflg
+
 
 !     Local variables
       integer              :: i,k,km1
       real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im)
       real(kind=kind_phys) :: mcons(im),fdqa(im),form(im,km),              &
-           qadv(im,km),dp(im,km),inbu(im,km)                         
+           dp(im,km),inbu(im,km)                         
 
 
       real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2,   &
@@ -77,21 +76,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
          mcons(i)=0.
       enddo
 
-      !Initial computations, dynamic q-tendency
-      if(flag_init .and. .not.flag_restart)then
-         do k = 1,km
-            do i = 1,im
-               qadv(i,k)=0.
-            enddo
-         enddo
-      else
-         do k = 1,km
-            do i = 1,im
-               qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
-            enddo
-         enddo
-      endif
-
       do k = 2,km1
           do i = 1,im
              if(cnvflg(i))then
@@ -133,7 +117,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
                 mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k))
                 buy2 = termD(i)+mcon+mcons(i)
 !               Do the integral over buoyant layers with positive mcon acc from surface
-                if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then
+                if(dbyo1(i,k)>0 .and. buy2 > 0.)then
                    inbu(i,k)=1.
                 endif
                 inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k))

physics/samfdeepcnv.f

@@ -102,7 +102,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       real(kind=kind_phys), intent(in) :: nthresh
       real(kind=kind_phys), intent(in) :: ca_deep(:)
       real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:),     &
-     &     tmf(:,:),q(:,:), prevsq(:,:)
+     &     tmf(:,:,:),q(:,:), prevsq(:,:)
       real(kind=kind_phys), intent(out) :: rainevap(:)
       real(kind=kind_phys), intent(out) :: sigmaout(:,:)
       logical, intent(in)  :: do_ca,ca_closure,ca_entr,ca_trigger
@@ -209,9 +209,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &                     bb1,     bb2,     wucb
 !
 !  parameters for prognostic sigma closure                                                                                                                                                      
-      real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km),
-     &     omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im)
-      real(kind=kind_phys) gravinv
+      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)
+      real(kind=kind_phys) gravinv,invdelt
       logical flag_shallow
 c  physical parameters
 !     parameter(grav=grav,asolfac=0.958)
@@ -306,6 +306,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       errflg = 0
 
       gravinv = 1./grav
+      invdelt = 1./delt
 
       elocp = hvap/cp
       el2orc = hvap*hvap/(rv*cp)
@@ -585,7 +586,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
          do i = 1, im
             dbyo1(i,k)=0.
             zdqca(i,k)=0.
-            qlks(i,k)=0.
             omega_u(i,k)=0.
             zeta(i,k)=1.0
          enddo
@@ -1515,7 +1515,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
                 pwavo(i) = pwavo(i) + pwo(i,k)
 !               cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
                 cnvwt(i,k) = etah * qlk * grav / dp
-                qlks(i,k)=qlk
+                zdqca(i,k)=dq/eta(i,k)
               endif
 !
 !  compute buoyancy and drag for updraft velocity
@@ -1690,7 +1690,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
                 pwavo(i) = pwavo(i) + pwo(i,k)
 !               cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
                 cnvwt(i,k) = etah * qlk * grav / dp
-                qlks(i,k)=qlk
+                zdqca(i,k)=dq/eta(i,k)
               endif
             endif
           endif
@@ -1860,28 +1860,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
           if(dq > 0.) then
             qlko_ktcon(i) = dq
             qcko(i,k) = qrch
-            qlks(i,k) = qlko_ktcon(i)
+            zdqca(i,k) = dq
           endif
         endif
       enddo
       endif
 c
 
-c store term needed for "termC" in prognostic area fraction closure
-      if(progsigma)then
-         do k = 2, km1
-            do i = 1, im
-               dp = 1000. * del(i,k)
-               if (cnvflg(i)) then
-                  if(k > kbcon(i) .and. k < ktcon(i)) then
-                     zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) +
-     &                    pwo(i,k)+dellal(i,k))
-                  endif
-               endif
-            enddo
-         enddo
-      endif
-
 ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then
 ccccc   print *, ' aa1(i) before dwndrft =', aa1(i)
 ccccc endif
@@ -2885,11 +2870,33 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 
 !> - From Bengtsson et al. (2022) \cite Bengtsson_2022 prognostic closure scheme, equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
       if(progsigma)then
+
+!Initial computations, dynamic q-tendency                                                                                                                                               
+         if(first_time_step .and. .not.restart)then
+            do k = 1,km
+               do i = 1,im
+                  qadv(i,k)=0.
+               enddo
+            enddo
+         else
+            do k = 1,km
+               do i = 1,im
+                  qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
+               enddo
+            enddo
+         endif
+
+         do k = 1,km
+            do i = 1,im
+               tmfq(i,k)=tmf(i,k,1)
+            enddo
+         enddo
+
          flag_shallow = .false.
          call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
-     &        del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
-     &        prevsq,q,kbcon1,ktcon,cnvflg,
-     &        sigmain,sigmaout,sigmab,errmsg,errflg)
+     &        del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
+     &        qadv,kbcon1,ktcon,cnvflg,
+     &        sigmain,sigmaout,sigmab)
       endif
 
 !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer.

physics/samfdeepcnv.meta

@@ -70,10 +70,10 @@
   type = logical
   intent = in
 [tmf]
-  standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
-  long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
+  standard_name = tendency_of_vertically_diffused_tracer_concentration
+  long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme
   units = kg kg-1 s-1
-  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers)
   type = real
   kind = kind_phys
   intent = in

physics/samfshalcnv.f

@@ -70,7 +70,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       real(kind=kind_phys), intent(in) ::  delt
       real(kind=kind_phys), intent(in) :: psp(:), delp(:,:),            &
      &   prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:),            &
-     &   qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:)
+     &   qmicro(:,:),tmf(:,:,:),prevsq(:,:),q(:,:)
 
       real(kind=kind_phys), intent(in) :: sigmain(:,:)
 !
@@ -156,10 +156,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 cc
 
 !  parameters for prognostic sigma closure
-      real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km),
+      real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
      &                     omegac(im),zeta(im,km),dbyo1(im,km),
-     &                     sigmab(im)
-      real(kind=kind_phys) gravinv,dxcrtas
+     &                     sigmab(im),qadv(im,km)
+      real(kind=kind_phys) gravinv,dxcrtas,invdelt
       logical flag_shallow
 c  physical parameters
 !     parameter(g=grav,asolfac=0.89)
@@ -247,6 +247,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       errflg = 0
 
       gravinv = 1./grav
+      invdelt = 1./delt
 
       elocp = hvap/cp
       el2orc = hvap*hvap/(rv*cp)
@@ -524,7 +525,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
          do i = 1, im
             dbyo1(i,k)=0.
             zdqca(i,k)=0.
-            qlks(i,k)=0.
             omega_u(i,k)=0.
             zeta(i,k)=1.0
          enddo
@@ -1270,7 +1270,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
                 qcko(i,k)= qlk + qrch
                 pwo(i,k) = etah * c0t(i,k) * dz * qlk
                 cnvwt(i,k) = etah * qlk * grav / dp
-                qlks(i,k)=qlk
+                zdqca(i,k)=dq/eta(i,k)
               endif
 !
 !  compute buoyancy and drag for updraft velocity
@@ -1435,7 +1435,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
                 qcko(i,k) = qlk + qrch
                 pwo(i,k) = etah * c0t(i,k) * dz * qlk
                 cnvwt(i,k) = etah * qlk * grav / dp
-                qlks(i,k)=qlk
+                zdqca(i,k)=dq/eta(i,k)
               endif
             endif
           endif
@@ -1601,24 +1601,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
           if(dq > 0.) then
             qlko_ktcon(i) = dq
             qcko(i,k) = qrch
-            qlks(i,k) = qlko_ktcon(i)
+            zdqca(i,k) = dq
           endif
         endif
       enddo
       endif
 c
 
-       do k = 2, km1
-        do i = 1, im
-           if (cnvflg(i)) then
-              if(k > kbcon(i) .and. k < ktcon(i)) then
-                 zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) +
-     &                pwo(i,k)+dellal(i,k))
-              endif
-           endif
-        enddo
-      enddo
-
 c--- compute precipitation efficiency in terms of windshear
 c
 !! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
@@ -1935,11 +1924,32 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 c
 !> - From Bengtsson et al. (2022) \cite Bengtsson_2022 prognostic closure scheme, equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
       if(progsigma)then
+!     Initial computations, dynamic q-tendency
+         if(first_time_step .and. .not.restart)then
+            do k = 1,km
+               do i = 1,im
+                  qadv(i,k)=0.
+               enddo
+            enddo
+         else
+            do k = 1,km
+               do i = 1,im
+                  qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
+               enddo
+            enddo
+         endif
+
+         do k = 1,km
+            do i = 1,im
+               tmfq(i,k)=tmf(i,k,1)
+            enddo
+         enddo
+
          flag_shallow = .true.
          call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
-     &        del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
-     &        prevsq,q,kbcon1,ktcon,cnvflg,
-     &        sigmain,sigmaout,sigmab,errmsg,errflg)
+     &        del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
+     &        qadv,kbcon1,ktcon,cnvflg,
+     &        sigmain,sigmaout,sigmab)
       endif
 
 !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity.

physics/samfshalcnv.meta

@@ -70,10 +70,10 @@
   type = logical
   intent = in
 [tmf]
-  standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
-  long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
+  standard_name = tendency_of_vertically_diffused_tracer_concentration
+  long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme
   units = kg kg-1 s-1
-  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers)
   type = real
   kind = kind_phys
   intent = in