Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed instability in upper atmosphere due to OGWD #111

Merged
merged 9 commits into from
Oct 23, 2023
26 changes: 24 additions & 2 deletions physics/drag_suite.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

!
Expand All @@ -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)
Expand Down Expand Up @@ -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
!
Expand Down Expand Up @@ -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)
Expand Down
Loading