From c0ec619536bd29740627cb2dc1da106b61dd435c Mon Sep 17 00:00:00 2001 From: Michael Toy Date: Tue, 19 Sep 2023 02:36:14 +0000 Subject: [PATCH] Added tendency limiter for mesosphere and horizontal wave number filter for orographic gravity wave drag in UGWP -- Issue #95 --- physics/drag_suite.F90 | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 22f122e71..ff68f4216 100644 --- a/physics/drag_suite.F90 +++ b/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)