Modified 'dtfac' OGWD tendency limiter in drag_suite_psl to match that in drag_suite_run

Author: mdtoyNOAA

physics/GWD/drag_suite.F90

@@ -1674,10 +1674,9 @@ subroutine drag_suite_psl(                                           &
    real(kind=kind_phys), parameter       ::  hmt_min = 50.
    real(kind=kind_phys), parameter       ::  oc_min  = 1.0
    real(kind=kind_phys), parameter       ::  oc_max  = 10.0
-! 7.5 mb -- 33 km ... 0.01 kgm-3 reduce gwd drag above cutoff level
-   real(kind=kind_phys), parameter       ::  pcutoff  = 7.5e2  
-! 0.76 mb -- 50 km ...0.001 kgm-3  --- 0.1 mb 65 km 0.0001 kgm-3
-   real(kind=kind_phys), parameter       ::  pcutoff_den  = 0.01  !  
+   real(kind=kind_phys), parameter       ::  sgmalolev  = 0.5  ! max sigma lvl for dtfac
+   real(kind=kind_phys), parameter       ::  plolevmeso = 70.0 ! pres lvl for mesosphere OGWD reduction (Pa)
+   real(kind=kind_phys), parameter       ::  facmeso    = 0.5  ! fractional velocity reduction for OGWD
 
    integer,parameter    ::  kpblmin = 2
 
@@ -1691,7 +1690,7 @@ subroutine drag_suite_psl(                                           &
                             rcsks,wdir,ti,rdz,tem2,dw2,shr2,      &
                             bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro,  &
                             rim,temc,tem1,efact,temv,dtaux,dtauy, &
-                            dtauxb,dtauyb,eng0,eng1
+                            dtauxb,dtauyb,eng0,eng1,dtfac_meso
    real(kind=kind_phys) ::  denfac
 !
    logical              ::  ldrag(im),icrilv(im),                 &
@@ -1706,7 +1705,7 @@ subroutine drag_suite_psl(                                           &
                             rulow(im),bnv(im),                    &
                             oa(im),ol(im),oc(im),                 &
                             oass(im),olss(im),                    &
-                            roll(im),dtfac(im,km),                   &
+                            roll(im),dtfac(im),                   &
                             brvf(im),xlinv(im),                   &
                             delks(im),delks1(im),                 &
                             bnv2(im,km),usqj(im,km),              &
@@ -1861,6 +1860,7 @@ subroutine drag_suite_psl(                                           &
      oass(i)       = 0.0
      olss(i)       = 0.0
      ulow (i)      = 0.0
+     dtfac(i)      = 1.0
      rstoch(i)     = 0.0
      ldrag(i)      = .false.
      icrilv(i)     = .false.
@@ -1877,7 +1877,6 @@ subroutine drag_suite_psl(                                           &
        taud_bl(i,k) = 0.0
        dtaux2d(i,k) = 0.0
        dtauy2d(i,k) = 0.0
-       dtfac(i,k)   = 1.0
      enddo
    enddo
 !
@@ -2536,30 +2535,38 @@ subroutine drag_suite_psl(                                           &
             taud_bl(i,klcap) = taud_bl(i,klcap) * factop
          enddo
 !
-!  if the gravity wave drag would force a critical line
-!  in the lower ksmm1 layers during the next deltim timestep,
-!  then only apply drag until that critical line is reached.
+!  if the gravity wave drag + blocking would force a critical line
+!  in the layers below pressure-based 'sigma' level = sgmalolev during the next deltim
+!  timestep, then only apply drag until that critical line is reached, i.e.,
+!  reduce drag to limit resulting wind components to zero
+!  Note: 'sigma' = prsi(k)/prsi(k=1), where prsi(k=1) is the surface pressure
 !
          do k = kts,kpblmax-1
-           if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.)                  then
-               dtfac(i,k) = min(dtfac(i,k),abs(velco(i,k)                 &
-                       /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
-           endif
-         enddo
-! apply limiter to mesosphere drag, reduce the drag by density factor 10-3
-! prevent wind reversal...
-!
-         do k = kpblmax,km-1
-           if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff) then
-               denfac = min(ro(i,k)/pcutoff_den,1.)
-               dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k)          &
+            if (prsi(i,k).ge.sgmalolev*prsi(i,1)) then 
+               if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.)                   &
+                  dtfac(i) = min(dtfac(i),abs(velco(i,k)                &
                        /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
-           endif
+            else
+               exit
+            endif
          enddo
 !
          do k = kts,km
-            taud_ls(i,k)  = taud_ls(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
-            taud_bl(i,k)  = taud_bl(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
+
+            ! Check if well into mesosphere -- if so, perform similar reduction of
+            ! velocity tendency due to mesoscale GWD to prevent sudden reversal of
+            ! wind direction (similar to above)
+            dtfac_meso = 1.0
+            if (prsl(i,k).le.plolevmeso) then 
+               if (taud_ls(i,k).ne.0.)                                  &    
+                  dtfac_meso = min(dtfac_meso,facmeso*abs(velco(i,k)    &    
+                     /(deltim*rcs*taud_ls(i,k))))
+            end if
+
+            taud_ls(i,k)  = taud_ls(i,k)*dtfac(i)*dtfac_meso*           &    
+                               ls_taper(i) *(1.-rstoch(i))
+            taud_bl(i,k)  = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i))
+
             dtaux  = taud_ls(i,k) * xn(i)
             dtauy  = taud_ls(i,k) * yn(i)
             dtauxb = taud_bl(i,k) * xn(i)