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
15 changes: 9 additions & 6 deletions physics/module_sf_noahmplsm.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion physics/noahmpdrv.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
24 changes: 18 additions & 6 deletions physics/sfc_diag_post.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
7 changes: 7 additions & 0 deletions physics/sfc_diag_post.meta
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 0 additions & 3 deletions physics/sfcsub.F
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading