Merge pull request #111 from mdtoyNOAA/ufs/dev_meso_fix_new_ksmax

Fixed instability in upper atmosphere due to OGWD

physics/drag_suite.F90

@@ -460,6 +460,8 @@ subroutine drag_suite_run(                                           &
    real(kind=kind_phys), parameter       ::  ce      = 0.8
    real(kind=kind_phys), parameter       ::  cg      = 0.5
    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
 
 !
@@ -472,7 +474,7 @@ subroutine drag_suite_run(                                           &
                             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,ksmax,dtfac_meso
 !
    logical              ::  ldrag(im),icrilv(im),                 &
                             flag(im),kloop1(im)
@@ -887,6 +889,14 @@ subroutine drag_suite_run(                                           &
          ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
          ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
          ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
+!  Check if mesoscale gravity waves will propagate vertically or be evanescent
+!  and not impart a drag force -- consider the maximum sub-grid horizontal
+!  topographic wavelength to be one-half the horizontal grid spacing -- calculate
+!  ksmax accordingly
+         ksmax = 4.0*pi/dx(i)   ! based on wavelength = 0.5*dx(i)
+         if ( bnv2(i,1).gt.0.0 ) then
+            ldrag(i) = ldrag(i) .or. sqrt(bnv2(i,1))*rulow(i).lt.ksmax
+         endif
 !
 !  set all ri low level values to the low level value
 !
@@ -1106,7 +1116,19 @@ subroutine drag_suite_run(                                           &
          enddo
 !
          do k = kts,km
-            taud_ms(i,k)  = taud_ms(i,k)*dtfac(i)* 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_ms(i,k).ne.0.)                                  &
+                  dtfac_meso = min(dtfac_meso,facmeso*abs(velco(i,k)    &
+                     /(deltim*rcs*taud_ms(i,k))))
+            end if
+
+            taud_ms(i,k)  = taud_ms(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_ms(i,k) * xn(i)

physics/module_sf_noahmplsm.F90

@@ -2116,7 +2116,7 @@ subroutine energy (parameters,ice    ,vegtyp ,ist    ,nsnow  ,nsoil  , & !in
 ! thermal properties of soil, snow, lake, and frozen soil
 
   call thermoprop (parameters,nsoil   ,nsnow   ,isnow   ,ist     ,dzsnso  , & !in
-                   dt      ,snowh   ,snice   ,snliq   , & !in
+                   dt      ,snowh   ,snice   ,snliq   , shdfac, & !in
                    smc     ,sh2o    ,tg      ,stc     ,ur      , & !in
                    lat     ,z0m     ,zlvl    ,vegtyp  ,  & !in
                    df      ,hcpct   ,snicev  ,snliqv  ,epore   , & !out
@@ -2463,7 +2463,7 @@ end subroutine energy
 
 !>\ingroup NoahMP_LSM
   subroutine thermoprop (parameters,nsoil   ,nsnow   ,isnow   ,ist     ,dzsnso  , & !in
-                         dt      ,snowh   ,snice   ,snliq   , & !in
+                         dt      ,snowh   ,snice   ,snliq   , shdfac,   & !in
                          smc     ,sh2o    ,tg      ,stc     ,ur      , & !in
                          lat     ,z0m     ,zlvl    ,vegtyp  , & !in
                          df      ,hcpct   ,snicev  ,snliqv  ,epore   , & !out
@@ -2480,6 +2480,7 @@ subroutine thermoprop (parameters,nsoil   ,nsnow   ,isnow   ,ist     ,dzsnso  ,
   real (kind=kind_phys)                           , intent(in)  :: dt      !< time step [s]
   real (kind=kind_phys), dimension(-nsnow+1:    0), intent(in)  :: snice   !< snow ice mass (kg/m2)
   real (kind=kind_phys), dimension(-nsnow+1:    0), intent(in)  :: snliq   !< snow liq mass (kg/m2)
+  real (kind=kind_phys)                           , intent(in)  :: shdfac !< green vegetation fraction [0.0-1.0]
   real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in)  :: dzsnso  !< thickness of snow/soil layers [m]
   real (kind=kind_phys), dimension(       1:nsoil), intent(in)  :: smc     !< soil moisture (ice + liq.) [m3/m3]
   real (kind=kind_phys), dimension(       1:nsoil), intent(in)  :: sh2o    !< liquid soil moisture [m3/m3]
@@ -2539,6 +2540,7 @@ subroutine thermoprop (parameters,nsoil   ,nsnow   ,isnow   ,ist     ,dzsnso  ,
 ! not in use because of the separation of the canopy layer from the ground.
 ! but this may represent the effects of leaf litter (niu comments)
 !       df1 = df1 * exp (sbeta * shdfac)
+        df(1) = df(1) * exp (sbeta * shdfac)
 
 ! compute lake thermal properties 
 ! (no consideration of turbulent mixing for this version)
@@ -4888,7 +4890,7 @@ subroutine bare_flux (parameters,nsnow   ,nsoil   ,isnow   ,dt      ,sag     , &
      end if
     endif ! 4
 
-! use sfc_diag to calculate t2mv and q2v for opt_sfc=1&3
+! use sfc_diag to calculate t2mb and q2b for opt_sfc=1&3
     if(opt_diag ==3) then
      if(opt_sfc == 1 .or. opt_sfc == 3) then
 
@@ -5823,7 +5825,8 @@ subroutine thermalz0(parameters,    fveg,          z0m, z0mg,       zlvl,
 
       elseif (opt_trs == chen09) then
 
-        z0m_out = exp(fveg * log(z0m)      + (1.0 - fveg) * log(z0mg))
+!       z0m_out = exp(fveg * log(z0m)      + (1.0 - fveg) * log(z0mg))
+        z0m_out = fveg * z0m      + (1.0 - fveg) * z0mg
         czil    = 10.0 ** (- 0.4 * parameters%hvt)
 
         reyn = ustarx*z0m_out/viscosity                      ! Blumel99 eqn 36c
@@ -5873,15 +5876,15 @@ subroutine thermalz0(parameters,    fveg,          z0m, z0mg,       zlvl,
 
         z0h_out = z0m_out
 
-      elseif (opt_trs == tessel) then
+      elseif (opt_trs == chen09 .or. opt_trs == tessel) then
 
         if (vegtyp <= 5) then
           z0h_out = z0m_out
         else
           z0h_out = z0m_out * 0.01
         endif
 
-      elseif (opt_trs == blumel99 .or. opt_trs == chen09) then
+      elseif (opt_trs == blumel99) then
 
         reyn = ustarx*z0m_out/viscosity                      ! Blumel99 eqn 36c
         if (reyn > 2.0) then

physics/noahmpdrv.F90

@@ -450,7 +450,7 @@ subroutine noahmpdrv_run                                       &
   integer    :: iopt_pedo = 1 ! option for pedotransfer function
   integer    :: iopt_crop = 0 ! option for crop model
   integer    :: iopt_gla  = 2 ! option for glacier treatment
-  integer    :: iopt_z0m  = 2 ! option for z0m treatment
+  integer    :: iopt_z0m  = 1 ! option for z0m treatment
 
 !
 !  ---  local inputs to noah-mp and glacier subroutines; listed in order in noah-mp call

physics/sfc_diag_post.F90

@@ -14,16 +14,17 @@ module sfc_diag_post
 !!
 #endif
       subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, opt_diag, dry, lssav, dtf, con_eps, con_epsm1, pgr,&
-                 t2mmp,q2mp, t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax,                  &
+                 vegtype,t2mmp,q2mp, t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax,                  &
                          wind10mmax, u10mmax, v10mmax, dpt2m, errmsg, errflg)
 
         use machine,               only: kind_phys, kind_dbl_prec
 
         implicit none
 
-        integer,                              intent(in) :: im, lsm, lsm_noahmp,opt_diag
-        logical,                              intent(in) :: lssav
-        real(kind=kind_phys),                 intent(in) :: dtf, con_eps, con_epsm1
+        integer,                             intent(in) :: im, lsm, lsm_noahmp,opt_diag
+        integer,              dimension(:),  intent(in) :: vegtype    !  vegetation type (integer index)
+        logical,                             intent(in) :: lssav
+        real(kind=kind_phys),                intent(in) :: dtf, con_eps, con_epsm1
         logical             , dimension(:),  intent(in) :: dry
         real(kind=kind_phys), dimension(:),  intent(in) :: pgr, u10m, v10m
         real(kind=kind_phys), dimension(:),  intent(inout) :: t2m, q2m, tmpmin, tmpmax, spfhmin, spfhmax
@@ -41,12 +42,23 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, opt_diag, dry, lssav, dtf, co
         errflg = 0
 
         if (lsm == lsm_noahmp) then
-         if (opt_diag == 2 .or. opt_diag == 3)then
+!     over shrublands use opt_diag=2
+          do i=1, im
+           if(dry(i)) then
+             if (vegtype(i) == 6 .or. vegtype(i) == 7  & 
+                .or. vegtype(i) == 16) then
+              t2m(i) = t2mmp(i)
+              q2m(i) = q2mp(i)
+             endif
+           endif
+          enddo
+
+         if (opt_diag == 2 .or. opt_diag == 3) then
           do i=1,im
             if(dry(i)) then
               t2m(i) = t2mmp(i)
               q2m(i) = q2mp(i)
-            endif
+             endif
           enddo
          endif
         endif

physics/sfc_diag_post.meta

@@ -81,6 +81,13 @@
   type = real
   kind = kind_phys
   intent = in
+[vegtype]
+  standard_name = vegetation_type_classification
+  long_name = vegetation type at each grid cell
+  units = index
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent= in
 [t2mmp]
   standard_name = temperature_at_2m_from_noahmp
   long_name = 2 meter temperature from noahmp

physics/sfcsub.F

@@ -7491,9 +7491,6 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw,     &
           endif
           call abort
         endif
-!
-!  soil type
-        print *,'in FIXREAD fnsotc =',fnsotc 
 !
         if(fnsotc(1:8).ne.'        ') then
           if ( index(fnsotc, "tileX.nc") == 0) then ! grib file