!
module gw_drag 5,8
!---------------------------------------------------------------------------------
! Purpose:
!
! Module to compute the forcing due to parameterized gravity waves. Both an
! orographic and an internal source spectrum are considered.
!
! Author: Byron Boville
!
!---------------------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver
use physics_types, only: physics_state, physics_ptend
use cam_history, only: outfld
use scamMod, only: single_column
use cam_logfile, only: iulog
use abortutils, only: endrun
implicit none
save
private ! Make default type private to the module
!
! PUBLIC: interfaces
!
public gw_drag_readnl ! Read namelist
public gw_inti ! Initialization
public gw_intr ! interface to actual parameterization
!
! PRIVATE: Rest of the data and interfaces are private to this module
!
real(r8), parameter :: unset_r8 = huge(1._r8)
! fcrit2 has been made a namelist variable to facilitate backwards compatibility
! with the CAM3 version of this parameterization. In CAM3 fcrit2=0.5
real(r8) :: fcrit2 = unset_r8 ! critical froude number squared
integer, parameter :: pgwv = 0 ! number of waves allowed
integer :: kbotbg, kbotoro ! interface of gwd source
integer :: ktopbg, ktoporo ! top interface of gwd region
real(r8) :: alpha(0:pver) ! newtonian cooling coefficients
real(r8) :: c(-pgwv:pgwv) ! list of wave phase speeds
real(r8) :: cpair ! specific heat of dry air (constant p)
real(r8) :: dback ! background diffusivity
real(r8) :: effkwv ! effective wavenumber (fcrit2*kwv)
real(r8) :: effgw_oro ! tendency efficiency for orographic gw
real(r8) :: effgw_spec=.125_r8 ! tendency efficiency for internal gw
real(r8) :: fracldv ! fraction of stress deposited in low level region
real(r8) :: g ! acceleration of gravity
real(r8) :: kwv ! effective horizontal wave number
real(r8) :: lzmax ! maximum vertical wavelength at source
real(r8) :: mxasym ! max asymmetry between tau(c) and tau(-c)
real(r8) :: mxrange ! max range of tau for all c
real(r8) :: n2min ! min value of bouyancy frequency
real(r8) :: oroko2 ! 1/2 * horizontal wavenumber
real(r8) :: orohmin ! min surface displacment height for orographic waves
real(r8) :: orovmin ! min wind speed for orographic waves
real(r8) :: r ! gas constant for dry air
real(r8) :: rog ! r / g
real(r8) :: taubgnd ! background source strength (/tauscal)
real(r8) :: taumin ! minimum (nonzero) stress
real(r8) :: tauscal ! scale factor for background stress source
real(r8) :: tndmax ! maximum wind tendency
real(r8) :: umcfac ! factor to limit tendency to prevent reversing u-c
real(r8) :: ubmc2mn ! min (u-c)**2
real(r8) :: zldvcon ! constant for determining zldv from tau0
!===============================================================================
contains
!===============================================================================
subroutine gw_drag_readnl(nlfile) 1,9
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'gw_drag_readnl'
namelist /gw_drag_nl/ fcrit2
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'gw_drag_nl', status=ierr)
if (ierr == 0) then
read(unitn, gw_drag_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(fcrit2, 1, mpir8, 0, mpicom)
#endif
! Error checking:
! check if fcrit2 was set.
if (fcrit2 == unset_r8) then
call endrun('gw_dray_readnl: fcrit2 must be set via the namelist')
end if
end subroutine gw_drag_readnl
!================================================================================
subroutine gw_inti (cpairx, cpwv, gx, rx, hypi) 1,16
!-----------------------------------------------------------------------
! Time independent initialization for multiple gravity wave parameterization.
!-----------------------------------------------------------------------
use cam_history, only: addfld, add_default, phys_decomp
use dycore, only: dycore_is
use phys_control,only: phys_getopts
!------------------------------Arguments--------------------------------
real(r8), intent(in) :: cpairx ! specific heat of dry air (constant p)
real(r8), intent(in) :: cpwv ! specific heat of water vapor (constant p)
real(r8), intent(in) :: gx ! acceleration of gravity
real(r8), intent(in) :: rx ! gas constant for dry air
real(r8), intent(in) :: hypi(pver+1) ! reference interface pressures
!---------------------------Local storage-------------------------------
integer :: k
logical :: history_budget ! output tendencies and state variables for CAM4
! temperature, water vapor, cloud ice and cloud
! liquid budgets.
!-----------------------------------------------------------------------
! Copy model constants
cpair = cpairx
g = gx
r = rx
! Set MGWD constants
kwv = 6.28e-5_r8 ! 100 km wave length
dback = 0.05_r8 ! background diffusivity
tauscal= 0.001_r8 ! scale factor for background stress
taubgnd= 6.4_r8 ! background stress amplitude
zldvcon= 10._r8 ! constant for determining zldv
lzmax = 7.E3_r8 ! maximum vertical wavelength at source (m)
if ( dycore_is ('LR') ) then
effgw_oro = 0.125_r8
fracldv= 0.0_r8
! Restore these to work with turulent mountain stress
! effgw_oro = 1.0
! fracldv= 0.7 ! fraction of tau0 diverged in low level region
else
effgw_oro = 0.125_r8
fracldv= 0.0_r8
endif
! Set phase speeds
do k = -pgwv, pgwv
c(k) = 10._r8 * k ! 0, +/- 10, +/- 20, ... m/s
end do
if (masterproc) then
write(iulog,*) ' '
write(iulog,*) 'GW_INTI: pgwv = ', pgwv
write(iulog,*) 'GW_INTI: c(l) = ', c
write(iulog,*) 'GW_INTI: effgw_oro = ', effgw_oro
write(iulog,*) 'GW_INTI: kwv = ', kwv
write(iulog,*) 'GW_INTI: fcrit2 = ', fcrit2
write(iulog,*) ' '
end if
! Set radiative damping times
do k = 0, pver
alpha(k) = 1.e-6_r8 ! about 10 days.
end do
! Min and max values to keep things reasonable
mxasym = 0.1_r8 ! max factor of 10 from |tau(c)| to |tau(-c)|
mxrange= 0.001_r8 ! factor of 100 from max to min |tau(c)|
n2min = 1.e-8_r8 ! min value of Brunt-Vaisalla freq squared
orohmin= 10._r8 ! min surface displacement for orographic wave drag
orovmin= 2._r8 ! min wind speed for orographic wave drag
taumin = 1.e-10_r8 ! min stress considered > 0
tndmax = 500._r8 / 86400._r8 ! max permitted tendency (500 m/s/day)
umcfac = 0.5_r8 ! max permitted reduction in u-c
ubmc2mn= 0.01_r8 ! min value of (u-c)^2
! Determine other derived constants
oroko2 = 0.5_r8 * kwv
effkwv = fcrit2 * kwv
rog = r/g
! Determine the bounds of the background and orographic stress regions
ktopbg = 0
kbotoro = pver
do k = 0, pver
if (hypi(k+1) .lt. 10000._r8) kbotbg = k ! spectrum source at 100 mb
!!$ if (hypi(k+1) .lt. 3000.) ktoporo = k
end do
ktoporo = 0
if (masterproc) then
write(iulog,*) 'KTOPBG =',ktopbg
write(iulog,*) 'KBOTBG =',kbotbg
write(iulog,*) 'KTOPORO =',ktoporo
write(iulog,*) 'KBOTORO =',kbotoro
end if
! Declare history variables for orgraphic term
call addfld ('TTGWORO ','K/s ',pver, 'A','T tendency - orographic gravity wave drag',phys_decomp)
call addfld ('UTGWORO ','m/s2 ',pver, 'A','U tendency - orographic gravity wave drag',phys_decomp)
call addfld ('VTGWORO ','m/s2 ',pver, 'A','V tendency - orographic gravity wave drag',phys_decomp)
call addfld ('TAUGWX ','N/m2 ',1, 'A','Zonal gravity wave surface stress', phys_decomp)
call addfld ('TAUGWY ','N/m2 ',1, 'A','Meridional gravity wave surface stress', phys_decomp)
call phys_getopts(history_budget_out = history_budget)
if ( history_budget ) then
call add_default ('TTGWORO', 1, ' ')
end if
! call add_default ('UTGWORO ', 1, ' ')
! call add_default ('VTGWORO ', 1, ' ')
! call add_default ('TAUGWX ', 1, ' ')
! call add_default ('TAUGWY ', 1, ' ')
! Declare history variables for spectrum
if (pgwv > 0) then
call addfld ('TTGWSPEC','K/s ',pver, 'A','T tendency - gravity wave spectrum', phys_decomp)
call addfld ('UTGWSPEC','m/s2 ',pver, 'A','U tendency - gravity wave spectrum', phys_decomp)
call addfld ('VTGWSPEC','m/s2 ',pver, 'A','V tendency - gravity wave spectrum', phys_decomp)
call add_default ('UTGWSPEC', 1, ' ')
call add_default ('VTGWSPEC', 1, ' ')
end if
return
end subroutine gw_inti
!===============================================================================
subroutine gw_intr (state, sgh, pblh, dt, ptend, landfrac, kvh) 1,13
!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography
real(r8), intent(in) :: pblh(pcols) ! planetary boundary layer height
real(r8), intent(in) :: dt ! time step
real(r8), intent(in) :: landfrac(pcols) ! Land fraction
real(r8), intent(in) :: kvh(pcols,0:pver) ! molecular thermal diffusivity
type(physics_state), intent(in) :: state ! physics state structure
type(physics_ptend), intent(inout):: ptend ! parameterization tendency structure
!---------------------------Local storage-------------------------------
integer :: lchnk ! chunk identifier
integer :: ncol ! number of atmospheric columns
integer :: i,k ! loop indexes
integer :: kldv(pcols) ! top interface of low level stress divergence region
integer :: kldvmn ! min value of kldv
integer :: ksrc(pcols) ! index of top interface of source region
integer :: ksrcmn ! min value of ksrc
real(r8) :: ttgw(pcols,pver) ! temperature tendency
real(r8) :: utgw(pcols,pver) ! zonal wind tendency
real(r8) :: vtgw(pcols,pver) ! meridional wind tendency
real(r8) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency
real(r8) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency
real(r8) :: rdpldv(pcols) ! 1/dp across low level divergence region
real(r8) :: rhoi(pcols,0:pver) ! interface density
real(r8) :: tau(pcols,-pgwv:pgwv,0:pver) ! wave Reynolds stress
real(r8) :: tau0x(pcols) ! c=0 sfc. stress (zonal)
real(r8) :: tau0y(pcols) ! c=0 sfc. stress (meridional)
real(r8) :: ti(pcols,0:pver) ! interface temperature
real(r8) :: ubi(pcols,0:pver) ! projection of wind at interfaces
real(r8) :: ubm(pcols,pver) ! projection of wind at midpoints
real(r8) :: xv(pcols) ! unit vectors of source wind (x)
real(r8) :: yv(pcols) ! unit vectors of source wind (y)
!-----------------------------------------------------------------------------
lchnk = state%lchnk
ncol = state%ncol
! Profiles of background state variables
call gw_prof(lchnk, ncol, &
state%u , state%v , state%t , state%pmid , state%pint, &
rhoi , ni , ti , nm)
!-----------------------------------------------------------------------------
! Non-orographic backgound gravity wave spectrum
!-----------------------------------------------------------------------------
if (pgwv >0) then
! Determine the wave source for a background spectrum at ~100 mb
call gw_bgnd (lchnk , ncol , &
state%u , state%v , state%t , state%pmid , state%pint , &
state%pdel , state%rpdel, state%lnpint,kldv , kldvmn , &
ksrc , ksrcmn , rdpldv , tau , ubi , &
ubm , xv , yv , PGWV , kbotbg )
! Solve for the drag profile
call gw_drag_prof (lchnk , ncol , &
PGWV , kbotbg , ktopbg , state%u , state%v , &
state%t , state%pint , state%pdel , state%rpdel, state%lnpint,&
rhoi , ni , ti , nm , dt , &
kldv , kldvmn , ksrc , ksrcmn , rdpldv , &
tau , ubi , ubm , xv , yv , &
effgw_spec , utgw , vtgw , tau0x , tau0y )
! Add the momentum tendencies to the output tendency arrays
do k = 1, pver
do i = 1, ncol
ptend%u(i,k) = utgw(i,k)
ptend%v(i,k) = vtgw(i,k)
end do
end do
! Write output fields to history file
call outfld ('UTGWSPEC', utgw, pcols, lchnk)
call outfld ('VTGWSPEC', vtgw, pcols, lchnk)
! zero net tendencies if no spectrum computed
else
ptend%u = 0._r8
ptend%v = 0._r8
end if
!-----------------------------------------------------------------------------
! Orographic stationary gravity wave
!-----------------------------------------------------------------------------
! Determine the orographic wave source
call gw_oro (lchnk, ncol, &
state%u , state%v , state%t , sgh , state%pmid , &
state%pint , state%pdel , state%zm , nm , pblh , &
kldv , kldvmn , ksrc , ksrcmn , rdpldv , &
tau , ubi , ubm , xv , yv )
! Solve for the drag profile
call gw_drag_prof (lchnk, ncol, &
0 , kbotoro , ktoporo , state%u , state%v , &
state%t , state%pint , state%pdel , state%rpdel, state%lnpint,&
rhoi , ni , ti , nm , dt , &
kldv , kldvmn , ksrc , ksrcmn , rdpldv , &
tau , ubi , ubm , xv , yv , &
effgw_oro , utgw , vtgw , tau0x , tau0y )
! Add the orographic tendencies to the spectrum tendencies
! Compute the temperature tendency from energy conservation (includes spectrum).
do k = 1, pver
do i = 1, ncol
ptend%u(i,k) = ptend%u(i,k) + utgw(i,k) * landfrac(i)
ptend%v(i,k) = ptend%v(i,k) + vtgw(i,k) * landfrac(i)
ptend%s(i,k) = -(ptend%u(i,k) * (state%u(i,k) + ptend%u(i,k)*0.5_r8*dt) &
+ptend%v(i,k) * (state%v(i,k) + ptend%v(i,k)*0.5_r8*dt))
ttgw(i,k) = ptend%s(i,k) / cpair
utgw(i,k) = ptend%u(i,k)
vtgw(i,k) = ptend%v(i,k)
end do
end do
! Set flags for nonzero tendencies, q not yet affected by gwd
ptend%name = "vertical diffusion"
ptend%lq(:) = .FALSE.
ptend%ls = .TRUE.
ptend%lu = .TRUE.
ptend%lv = .TRUE.
! Write output fields to history file
call outfld ('UTGWORO', utgw, pcols, lchnk)
call outfld ('VTGWORO', vtgw, pcols, lchnk)
call outfld ('TTGWORO', ttgw, pcols, lchnk)
call outfld ('TAUGWX', tau0x, pcols, lchnk)
call outfld ('TAUGWY', tau0y, pcols, lchnk)
call outfld ('SGH ', sgh, pcols, lchnk)
return
end subroutine gw_intr
!===============================================================================
subroutine gw_prof (lchnk, ncol, u, v, t, pm, pi, rhoi, ni, ti, nm) 1
!-----------------------------------------------------------------------
! Compute profiles of background state quantities for the multiple
! gravity wave drag parameterization.
!
! The parameterization is assumed to operate only where water vapor
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind
real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind
real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures
real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures
real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures
real(r8), intent(out) :: rhoi(pcols,0:pver) ! interface density
real(r8), intent(out) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency
real(r8), intent(out) :: ti(pcols,0:pver) ! interface temperature
real(r8), intent(out) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency
!---------------------------Local storage-------------------------------
integer :: i,k ! loop indexes
real(r8) :: dtdp
real(r8) :: n2 ! Brunt-Vaisalla frequency squared
!-----------------------------------------------------------------------------
! Determine the interface densities and Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------
! The top interface values are calculated assuming an isothermal atmosphere
! above the top level.
k = 0
do i = 1, ncol
ti(i,k) = t(i,k+1)
rhoi(i,k) = pi(i,k) / (r*ti(i,k))
ni(i,k) = sqrt (g*g / (cpair*ti(i,k)))
end do
! Interior points use centered differences
do k = 1, pver-1
do i = 1, ncol
ti(i,k) = 0.5_r8 * (t(i,k) + t(i,k+1))
rhoi(i,k) = pi(i,k) / (r*ti(i,k))
dtdp = (t(i,k+1)-t(i,k)) / (pm(i,k+1)-pm(i,k))
n2 = g*g/ti(i,k) * (1._r8/cpair - rhoi(i,k)*dtdp)
ni(i,k) = sqrt (max (n2min, n2))
end do
end do
! Bottom interface uses bottom level temperature, density; next interface
! B-V frequency.
k = pver
do i = 1, ncol
ti(i,k) = t(i,k)
rhoi(i,k) = pi(i,k) / (r*ti(i,k))
ni(i,k) = ni(i,k-1)
end do
!-----------------------------------------------------------------------------
! Determine the midpoint Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------
do k=1,pver
do i=1,ncol
nm(i,k) = 0.5_r8 * (ni(i,k-1) + ni(i,k))
end do
end do
return
end subroutine gw_prof
!===============================================================================
subroutine gw_oro (lchnk, ncol, & 1
u, v, t, sgh, pm, pi, dpm, zm, nm, pblh, &
kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv)
!-----------------------------------------------------------------------
! Orographic source for multiple gravity wave drag parameterization.
!
! The stress is returned for a single wave with c=0, over orography.
! For points where the orographic variance is small (including ocean),
! the returned stress is zero.
!------------------------------Arguments--------------------------------
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind
real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind
real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures
real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography
real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures
real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures
real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1))
real(r8), intent(in) :: zm(pcols,pver) ! midpoint heights
real(r8), intent(in) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency
real(r8), intent(in) :: pblh(pcols) ! planetary boundary layer height
integer, intent(out) :: kldv(pcols) ! top interface of low level stress div region
integer, intent(out) :: kldvmn ! min value of kldv
integer, intent(out) :: ksrc(pcols) ! index of top interface of source region
integer, intent(out) :: ksrcmn ! min value of ksrc
real(r8), intent(out) :: rdpldv(pcols) ! 1/dp across low level divergence region
real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress
real(r8), intent(out) :: ubi(pcols,0:pver) ! projection of wind at interfaces
real(r8), intent(out) :: ubm(pcols,pver) ! projection of wind at midpoints
real(r8), intent(out) :: xv(pcols) ! unit vectors of source wind (x)
real(r8), intent(out) :: yv(pcols) ! unit vectors of source wind (y)
!---------------------------Local storage-------------------------------
integer :: i,k ! loop indexes
real(r8) :: pil ! don't we have pi somewhere?
real(r8) :: lzsrc ! vertical wavelength at source
real(r8) :: hdsp(pcols) ! surface streamline displacment height (2*sgh)
real(r8) :: sghmax ! max orographic sdv to use
real(r8) :: tauoro(pcols) ! c=0 stress from orography
real(r8) :: zldv(pcols) ! top of the low level stress divergence region
real(r8) :: nsrc(pcols) ! b-f frequency averaged over source region
real(r8) :: psrc(pcols) ! interface pressure at top of source region
real(r8) :: rsrc(pcols) ! density averaged over source region
real(r8) :: usrc(pcols) ! u wind averaged over source region
real(r8) :: vsrc(pcols) ! v wind averaged over source region
!---------------------------------------------------------------------------
! Average the basic state variables for the wave source over the depth of
! the orographic standard deviation. Here we assume that the apropiate
! values of wind, stability, etc. for determining the wave source are
! averages over the depth of the atmosphere pentrated by the typical mountain.
! Reduces to the bottom midpoint values when sgh=0, such as over ocean.
!---------------------------------------------------------------------------
k = pver
do i = 1, ncol
ksrc(i) = k-1
psrc(i) = pi(i,k-1)
rsrc(i) = pm(i,k)/(r*t(i,k)) * dpm(i,k)
usrc(i) = u(i,k) * dpm(i,k)
vsrc(i) = v(i,k) * dpm(i,k)
nsrc(i) = nm(i,k)* dpm(i,k)
hdsp(i) = 2.0_r8 * sgh(i)
end do
do k = pver-1, pver/2, -1
do i = 1, ncol
if (hdsp(i) > sqrt(zm(i,k)*zm(i,k+1))) then
ksrc(i) = k-1
psrc(i) = pi(i,k-1)
rsrc(i) = rsrc(i) + pm(i,k) / (r*t(i,k))* dpm(i,k)
usrc(i) = usrc(i) + u(i,k) * dpm(i,k)
vsrc(i) = vsrc(i) + v(i,k) * dpm(i,k)
nsrc(i) = nsrc(i) + nm(i,k)* dpm(i,k)
end if
end do
end do
do i = 1, ncol
rsrc(i) = rsrc(i) / (pi(i,pver) - psrc(i))
usrc(i) = usrc(i) / (pi(i,pver) - psrc(i))
vsrc(i) = vsrc(i) / (pi(i,pver) - psrc(i))
nsrc(i) = nsrc(i) / (pi(i,pver) - psrc(i))
if (single_column) then
! needed the following fix when winds are identically 0
! orig -> ubi(i,pver) = sqrt (usrc(i)**2 + vsrc(i)**2)
ubi(i,pver) = max(sqrt (usrc(i)**2 + vsrc(i)**2),orovmin)
else
ubi(i,pver) = sqrt (usrc(i)**2 + vsrc(i)**2)
#if (defined BFB_CAM_SCAM_IOP )
ubi(i,pver) = max(sqrt (usrc(i)**2 + vsrc(i)**2),orovmin)
#endif
endif
xv(i) = usrc(i) / ubi(i,pver)
yv(i) = vsrc(i) / ubi(i,pver)
end do
! Project the local wind at midpoints onto the source wind.
do k = 1, pver
do i = 1, ncol
ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
end do
end do
! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
do i = 1, ncol
ubi(i,0) = ubm(i,1)
end do
do k = 1, pver-1
do i = 1, ncol
ubi(i,k) = 0.5_r8 * (ubm(i,k) + ubm(i,k+1))
end do
end do
!---------------------------------------------------------------------------
! Determine the depth of the low level stress divergence region, as
! the max of the source region depth and 1/2 the vertical wavelength at the
! source.
!---------------------------------------------------------------------------
pil = acos(-1._r8)
do i = 1, ncol
lzsrc = min(2._r8 * pil * usrc(i) / nsrc(i), lzmax)
zldv(i) = max(hdsp(i), 0.5_r8 * lzsrc)
end do
! find the index of the top of the low level divergence region
kldv(:) = pver-1
do k = pver-1, pver/2, -1
do i = 1, ncol
if (zldv(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
kldv(i) = k-1
end if
end do
end do
! Determine the orographic c=0 source term following McFarlane (1987).
! Set the source top interface index to pver, if the orographic term is zero.
do i = 1, ncol
if ((ubi(i,pver) .gt. orovmin) .and. (hdsp(i) .gt. orohmin)) then
sghmax = fcrit2 * (ubi(i,pver) / nsrc(i))**2
tauoro(i) = oroko2 * min(hdsp(i)**2, sghmax) * rsrc(i) * nsrc(i) * ubi(i,pver)
else
tauoro(i) = 0._r8
ksrc(i) = pver
kldv(i) = pver
end if
end do
! Set the phase speeds and wave numbers in the direction of the source wind.
! Set the source stress magnitude (positive only, note that the sign of the
! stress is the same as (c-u).
do i = 1, ncol
tau(i,0,pver) = tauoro(i)
end do
! Determine the min value of kldv and ksrc for limiting later loops
! and the pressure at the top interface of the low level stress divergence
! region.
ksrcmn = pver
kldvmn = pver
do i = 1, ncol
ksrcmn = min(ksrcmn, ksrc(i))
kldvmn = min(kldvmn, kldv(i))
if (kldv(i) .ne. pver) then
rdpldv(i) = 1._r8 / (pi(i,kldv(i)) - pi(i,pver))
end if
end do
if (fracldv .le. 0._r8) kldvmn = pver
return
end subroutine gw_oro
!===============================================================================
subroutine gw_bgnd (lchnk, ncol, & 1,2
u, v, t, pm, pi, dpm, rdpm, piln, &
kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
ngwv, kbot)
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
!
! The parameterization is assumed to operate only where water vapor
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------
use phys_grid, only: get_rlat_all_p
!------------------------------Arguments--------------------------------
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: kbot ! index of bottom (source) interface
integer, intent(in) :: ngwv ! number of gravity waves to use
real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind
real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind
real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures
real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures
real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures
real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1))
real(r8), intent(in) :: rdpm(pcols,pver) ! 1. / (pi(k)-pi(k-1))
real(r8), intent(in) :: piln(pcols,0:pver) ! ln(interface pressures)
integer, intent(out) :: kldv(pcols) ! top interface of low level stress divergence region
integer, intent(out) :: kldvmn ! min value of kldv
integer, intent(out) :: ksrc(pcols) ! index of top interface of source region
integer, intent(out) :: ksrcmn ! min value of ksrc
real(r8), intent(in) :: rdpldv(pcols) ! 1/dp across low level divergence region
real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress
real(r8), intent(out) :: ubi(pcols,0:pver) ! projection of wind at interfaces
real(r8), intent(out) :: ubm(pcols,pver) ! projection of wind at midpoints
real(r8), intent(out) :: xv(pcols) ! unit vectors of source wind (x)
real(r8), intent(out) :: yv(pcols) ! unit vectors of source wind (y)
!---------------------------Local storage-------------------------------
integer :: i,k,l ! loop indexes
real(r8) :: rlat(pcols) ! latitude in radians for columns
real(r8) :: tauback(pcols) ! background stress at c=0
real(r8) :: usrc(pcols) ! u wind averaged over source region
real(r8) :: vsrc(pcols) ! v wind averaged over source region
real(r8) :: al0 ! Used in lat dependence of GW spec.
real(r8) :: dlat0 ! Used in lat dependence of GW spec.
real(r8) :: flat_gw ! The actual lat dependence of GW spec.
real(r8) :: pi_g ! 3.14........
!---------------------------------------------------------------------------
! Determine the source layer wind and unit vectors, then project winds.
!---------------------------------------------------------------------------
! Just use the source level interface values for the source
! wind speed and direction (unit vector).
k = kbot
do i = 1, ncol
ksrc(i) = k
kldv(i) = k
usrc(i) = 0.5_r8*(u(i,k+1)+u(i,k))
vsrc(i) = 0.5_r8*(v(i,k+1)+v(i,k))
ubi(i,kbot) = sqrt (usrc(i)**2 + vsrc(i)**2)
xv(i) = usrc(i) / ubi(i,k)
yv(i) = vsrc(i) / ubi(i,k)
end do
! Project the local wind at midpoints onto the source wind.
do k = 1, kbot
do i = 1, ncol
ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
end do
end do
! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
do i = 1, ncol
ubi(i,0) = ubm(i,1)
end do
do k = 1, kbot-1
do i = 1, ncol
ubi(i,k) = 0.5_r8 * (ubm(i,k) + ubm(i,k+1))
end do
end do
!-----------------------------------------------------------------------
! Gravity wave sources
!-----------------------------------------------------------------------
! Determine the background stress at c=0
do i=1,ncol
tauback(i) = taubgnd * tauscal
enddo
! Include dependence on latitude:
! The lat function was obtained by RR Garcia as
! currently used in his 2D model
! [Added by F. Sassi on May 30, 1997]
pi_g=4._r8*atan(1._r8)
al0=40._r8*pi_g/180._r8
dlat0=10._r8*pi_g/180._r8
!
call get_rlat_all_p(lchnk, ncol, rlat)
do i=1,ncol
flat_gw= 0.5_r8*(1._r8+tanh((rlat(i)-al0)/dlat0)) + 0.5_r8*(1._r8+tanh(-(rlat(i)+al0)/dlat0))
flat_gw=max(flat_gw,0.2_r8)
tauback(i)=tauback(i)*flat_gw
enddo
! Set the phase speeds and wave numbers in the direction of the source wind.
! Set the source stress magnitude (positive only, note that the sign of the
! stress is the same as (c-u).
do l = 1, ngwv
do i = 1, ncol
tau(i, l,kbot) = tauback(i) * exp(-(c(l)/30._r8)**2)
tau(i,-l,kbot) = tau(i, l,kbot)
end do
end do
do i = 1, ncol
tau(i,0,kbot) = tauback(i)
end do
! Determine the min value of kldv and ksrc for limiting later loops
! and the pressure at the top interface of the low level stress divergence
! region.
ksrcmn = pver
kldvmn = pver
return
end subroutine gw_bgnd
!===============================================================================
subroutine gw_drag_prof (lchnk, ncol, & 2
ngwv, kbot, ktop, u, v, t, &
pi, dpm, rdpm, piln, rhoi, ni, ti, nm, dt, &
kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
effgw, ut, vt, tau0x, tau0y)
!-----------------------------------------------------------------------
! Solve for the drag profile from the multiple gravity wave drag
! parameterization.
! 1. scan up from the wave source to determine the stress profile
! 2. scan down the stress profile to determine the tendencies
! => apply bounds to the tendency
! a. from wkb solution
! b. from computational stability constraints
! => adjust stress on interface below to reflect actual bounded tendency
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: kbot ! index of bottom (source) interface
integer, intent(in) :: ktop ! index of top interface of gwd region
integer, intent(in) :: ngwv ! number of gravity waves to use
integer, intent(in) :: kldv(pcols) ! top interface of low level stress divergence region
integer, intent(in) :: kldvmn ! min value of kldv
integer, intent(in) :: ksrc(pcols) ! index of top interface of source region
integer, intent(in) :: ksrcmn ! min value of ksrc
real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind
real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind
real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures
real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures
real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1))
real(r8), intent(in) :: rdpm(pcols,pver) ! 1. / (pi(k)-pi(k-1))
real(r8), intent(in) :: piln(pcols,0:pver) ! ln(interface pressures)
real(r8), intent(in) :: rhoi(pcols,0:pver) ! interface density
real(r8), intent(in) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency
real(r8), intent(in) :: ti(pcols,0:pver) ! interface temperature
real(r8), intent(in) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency
real(r8), intent(in) :: dt ! time step
real(r8), intent(in) :: rdpldv(pcols) ! 1/dp across low level divergence region
real(r8), intent(in) :: ubi(pcols,0:pver) ! projection of wind at interfaces
real(r8), intent(in) :: ubm(pcols,pver) ! projection of wind at midpoints
real(r8), intent(in) :: xv(pcols) ! unit vectors of source wind (x)
real(r8), intent(in) :: yv(pcols) ! unit vectors of source wind (y)
real(r8), intent(in) :: effgw ! tendency efficiency
real(r8), intent(inout) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress
real(r8), intent(out) :: ut(pcols,pver) ! zonal wind tendency
real(r8), intent(out) :: vt(pcols,pver) ! meridional wind tendency
real(r8), intent(out) :: tau0x(pcols) ! c=0 sfc. stress (zonal)
real(r8), intent(out) :: tau0y(pcols) ! c=0 sfc. stress (meridional)
!---------------------------Local storage-------------------------------
integer :: i,k,l ! loop indexes
real(r8) :: d(pcols) ! "total" diffusivity
real(r8) :: dsat(pcols,-pgwv:pgwv) ! saturation diffusivity
real(r8) :: dscal ! fraction of dsat to use
real(r8) :: mi ! imaginary part of vertical wavenumber
real(r8) :: taudmp ! stress after damping
real(r8) :: taumax(pcols) ! max(tau) for any l
real(r8) :: tausat(pcols,-pgwv:pgwv) ! saturation stress
real(r8) :: ubmc(pcols,-pgwv:pgwv) ! (ub-c)
real(r8) :: ubmc2 ! (ub-c)**2
real(r8) :: ubt(pcols,pver) ! ubar tendency
real(r8) :: ubtl ! ubar tendency from wave l
real(r8) :: ubtlsat ! saturation tendency
! Initialize gravity wave drag tendencies to zero
do k=1,pver
do i=1,pcols
ut(i,k) = 0._r8
vt(i,k) = 0._r8
end do
end do
!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------
! Loop from bottom to top to get stress profiles
do k = kbot-1, ktop, -1
! Determine the absolute value of the saturation stress and the diffusivity
! for each wave.
! Define critical levels where the sign of (u-c) changes between interfaces.
d(:ncol) = dback
do l = -ngwv, ngwv
do i = 1, ncol
ubmc(i,l) = ubi(i,k) - c(l)
tausat(i,l) = abs (effkwv * rhoi(i,k) * ubmc(i,l)**3 / (2._r8*ni(i,k)) )
if (tausat(i,l) .le. taumin) tausat(i,l) = 0.0_r8
if (single_column) then
! needed the following fix when winds are identically 0
! ie. thermal equilibrium case
if ((ubi(i,k+1) .ne. c(l))) then
if (ubmc(i,l) / (ubi(i,k+1) - c(l)) .le. 0.0_r8) tausat(i,l) = 0.0_r8
else
tausat(i,l) = 0.0_r8
end if
else
if (ubmc(i,l) / (ubi(i,k+1) - c(l)) .le. 0.0_r8) tausat(i,l) = 0.0_r8
end if
dsat(i,l) = (ubmc(i,l) / ni(i,k))**2 * &
(effkwv * ubmc(i,l)**2 / (rog * ti(i,k) * ni(i,k)) - alpha(k))
dscal = min (1.0_r8, tau(i,l,k+1) / (tausat(i,l)+taumin))
d(i) = max( d(i), dscal * dsat(i,l))
end do
end do
! Compute stress for each wave. The stress at this level is the min of
! the saturation stress and the stress at the level below reduced by damping.
! The sign of the stress must be the same as at the level below.
do l = -ngwv, ngwv
do i = 1, ncol
ubmc2 = max(ubmc(i,l)**2, ubmc2mn)
mi = ni(i,k) / (2._r8 * kwv * ubmc2) * (alpha(k) + ni(i,k)**2/ubmc2 * d(i))
taudmp = tau(i,l,k+1) * exp(-2._r8*mi*rog*t(i,k+1)*(piln(i,k+1)-piln(i,k)))
if (taudmp .le. taumin) taudmp = 0._r8
tau(i,l,k) = min (taudmp, tausat(i,l))
end do
end do
! The orographic stress term does not change across the source region
! Note that k ge ksrcmn cannot occur without an orographic source term
if (k .ge. ksrcmn) then
do i = 1, ncol
if (k .ge. ksrc(i)) then
tau(i,0,k) = tau(i,0,pver)
end if
end do
end if
! Require that the orographic term decrease linearly (with pressure)
! within the low level stress divergence region. This supersedes the above
! requirment of constant stress within the source region.
! Note that k ge kldvmn cannot occur without an orographic source term, since
! kldvmn=pver then and k<=pver-1
if (k .ge. kldvmn) then
do i = 1, ncol
if (k .ge. kldv(i)) then
tau(i,0,k) = min (tau(i,0,k), tau(i,0,pver) * &
(1._r8 - fracldv * (pi(i,k)-pi(i,pver)) * rdpldv(i)))
end if
end do
end if
! Apply lower bounds to the stress if ngwv > 0.
if (ngwv .ge. 1) then
! Determine the max value of tau for any l
do i = 1, ncol
taumax(i) = tau(i,-ngwv,k)
end do
do l = -ngwv+1, ngwv
do i = 1, ncol
taumax(i) = max(taumax(i), tau(i,l,k))
end do
end do
do i = 1, ncol
taumax(i) = mxrange * taumax(i)
end do
! Set the min value of tau for each wave to the max of mxrange*taumax or
! mxasym*tau(-c)
do l = 1, ngwv
do i = 1, ncol
tau(i, l,k) = max(tau(i, l,k), taumax(i))
tau(i, l,k) = max(tau(i, l,k), mxasym*tau(i,-l,k))
tau(i,-l,k) = max(tau(i,-l,k), taumax(i))
tau(i,-l,k) = max(tau(i,-l,k), mxasym*tau(i, l,k))
end do
end do
l = 0
do i = 1, ncol
tau(i,l,k) = max(tau(i,l,k), mxasym * 0.5_r8 * (tau(i,l-1,k) + tau(i,l+1,k)))
end do
end if
end do
! Put an upper bound on the stress at the top interface to force some stress
! divergence in the top layer. This prevents the easterlies from running
! away in the summer mesosphere, since most of the gravity wave activity
! will pass through a top interface at 75--80 km under strong easterly
! conditions.
!++BAB fix to match ccm3.10
!!$ do l = -ngwv, ngwv
!!$ do i = 1, ncol
!!$ tau(i,l,0) = min(tau(i,l,0), 0.5*tau(i,l,1))
!!$ end do
!!$ end do
!--BAB fix to match ccm3.10
!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------
! Loop over levels from top to bottom
do k = ktop+1, kbot
! Accumulate the mean wind tendency over wavenumber.
do i = 1, ncol
ubt (i,k) = 0.0_r8
end do
do l = -ngwv, ngwv
do i = 1, ncol
! Determine the wind tendency including excess stress carried down from above.
ubtl = g * (tau(i,l,k)-tau(i,l,k-1)) * rdpm(i,k)
! Require that the tendency be no larger than the analytic solution for
! a saturated region [proportional to (u-c)^3].
ubtlsat = effkwv * abs((c(l)-ubm(i,k))**3) /(2._r8*rog*t(i,k)*nm(i,k))
ubtl = min(ubtl, ubtlsat)
! Apply tendency limits to maintain numerical stability.
! 1. du/dt < |c-u|/dt so u-c cannot change sign (u^n+1 = u^n + du/dt * dt)
! 2. du/dt < tndmax so that ridicuously large tendency are not permitted
ubtl = min(ubtl, umcfac * abs(c(l)-ubm(i,k)) / dt)
ubtl = min(ubtl, tndmax)
! Accumulate the mean wind tendency over wavenumber.
ubt (i,k) = ubt (i,k) + sign(ubtl, c(l)-ubm(i,k))
! Redetermine the effective stress on the interface below from the wind
! tendency. If the wind tendency was limited above, then the new stress
! will be small than the old stress and will cause stress divergence in
! the next layer down. This has the effect of smoothing large stress
! divergences downward while conserving total stress.
tau(i,l,k) = tau(i,l,k-1) + ubtl * dpm(i,k) / g
end do
end do
! Project the mean wind tendency onto the components and scale by "efficiency".
do i = 1, ncol
ut(i,k) = ubt(i,k) * xv(i) * effgw
vt(i,k) = ubt(i,k) * yv(i) * effgw
end do
! End of level loop
end do
!-----------------------------------------------------------------------
! Project the c=0 stress (scaled) in the direction of the source wind for recording
! on the output file.
!-----------------------------------------------------------------------
do i = 1, ncol
tau0x(i) = tau(i,0,kbot) * xv(i) * effgw
tau0y(i) = tau(i,0,kbot) * yv(i) * effgw
end do
return
end subroutine gw_drag_prof
end module gw_drag