diff --git a/bld/build-namelist b/bld/build-namelist
index bd91edc9f1..78e3c8c23c 100755
--- a/bld/build-namelist
+++ b/bld/build-namelist
@@ -3795,6 +3795,7 @@ my $do_gw_convect_sh = ($nl->get_value('use_gw_convect_sh') =~ /$TRUE/io);
my $do_gw_movmtn_pbl = ($nl->get_value('use_gw_movmtn_pbl') =~ /$TRUE/io);
my $do_gw_rdg_beta = ($nl->get_value('use_gw_rdg_beta') =~ /$TRUE/io);
my $do_gw_rdg_gamma = ($nl->get_value('use_gw_rdg_gamma') =~ /$TRUE/io);
+my $do_gw_rdg_resid = ($nl->get_value('use_gw_rdg_resid') =~ /$TRUE/io);
my $do_divstream = ($nl->get_value('gw_rdg_do_divstream') =~ /$TRUE/io);
@@ -3872,6 +3873,10 @@ if ($do_gw_rdg_beta) {
add_default($nl, 'gw_prndl');
}
+if ($do_gw_rdg_resid) {
+ add_default($nl, 'effgw_rdg_resid');
+}
+
if ($do_gw_rdg_gamma) {
add_default($nl, 'n_rdg_gamma', 'val'=>'-1');
add_default($nl, 'effgw_rdg_gamma', 'val'=>'1.0D0');
diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml
index 308e128304..3374ad09b4 100644
--- a/bld/namelist_files/namelist_defaults_cam.xml
+++ b/bld/namelist_files/namelist_defaults_cam.xml
@@ -836,6 +836,10 @@
0.5D0
0.5D0
+JULIO
+JULIO
+JULIO
+
0.03D0
diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml
index e0c4b46778..239c19f986 100644
--- a/bld/namelist_files/namelist_definition.xml
+++ b/bld/namelist_files/namelist_definition.xml
@@ -1338,6 +1338,12 @@ Whether or not to enable gravity waves from PBL moving mountains source.
Default: .false.
+
+JULIO - need to add comment
+Default: .false.
+
+
Gravity wave spectrum dimension (wave numbers are from -pgwv to pgwv).
@@ -1441,6 +1447,12 @@ Max efficiency associated with anisotropic OGW.
Default: 1.0
+
+JULIO - need to add comment and appropriate default
+Default: ?????? JULIO
+
+
Drag coefficient for obstacles in low-level flow.
diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90
index aeab27a5c6..3d50877cbe 100644
--- a/src/physics/cam/gw_drag.F90
+++ b/src/physics/cam/gw_drag.F90
@@ -110,6 +110,11 @@ module gw_drag
! Beres (shallow convection).
real(r8) :: effgw_beres_sh = unset_r8
+
+ ! JULIO - Please put in appropriate comment
+ logical :: use_gw_rdg_resid = false
+ read(r8) :: effgw_rdg_resid = unset_r8
+
! Horzontal wavelengths [m].
real(r8), parameter :: wavelength_mid = 1.e5_r8
real(r8), parameter :: wavelength_long = 3.e5_r8
@@ -168,7 +173,9 @@ module gw_drag
integer, parameter :: prdg = 16
real(r8), allocatable, dimension(:,:), target :: &
- rdg_gbxar
+ rdg_gbxar, &
+ rdg_isovar, &
+ rdg_isowgt
! Meso Beta
real(r8), allocatable, dimension(:,:,:), target :: &
@@ -250,7 +257,9 @@ subroutine gw_drag_readnl(nlfile)
rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, bnd_rdggm, &
gw_oro_south_fac, gw_limit_tau_without_eff, &
gw_lndscl_sgh, gw_prndl, gw_apply_tndmax, gw_qbo_hdepth_scaling, &
- gw_top_taper, front_gaussian_width, alpha_gw_movmtn
+ gw_top_taper, front_gaussian_width, alpha_gw_movmtn, use_gw_rdg_resid, &
+ effgw_rdg_resid
+
!----------------------------------------------------------------------
if (use_simple_phys) return
@@ -359,6 +368,11 @@ subroutine gw_drag_readnl(nlfile)
call mpi_bcast(alpha_gw_movmtn, 1, mpi_real8, mstrid, mpicom, ierr)
if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: alpha_gw_movmtn")
+ call mpi_bcast(use_gw_rdg_resid, 1, mpi_logical, mstrid, mpicom, ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_resid")
+ call mpi_bcast(effgw_rdg_resid, 1, mpi_real8, mstrid, mpicom, ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_resid")
+
! Check if fcrit2 was set.
call shr_assert(fcrit2 /= unset_r8, &
"gw_drag_readnl: fcrit2 must be set via the namelist."// &
@@ -648,6 +662,8 @@ subroutine gw_init()
! Get beta ridge data
allocate( &
rdg_gbxar(pcols,begchunk:endchunk), &
+ rdg_isovar(pcols,begchunk:endchunk), &
+ rdg_isowgt(pcols,begchunk:endchunk), &
rdg_hwdth(pcols,prdg,begchunk:endchunk), &
rdg_clngt(pcols,prdg,begchunk:endchunk), &
rdg_mxdis(pcols,prdg,begchunk:endchunk), &
@@ -659,6 +675,14 @@ subroutine gw_init()
if (.not. found) call endrun(sub//': ERROR: GBXAR not found on topo file')
rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2
+ call infld('ISOVAR', fh_topo, dim1name, dim2name, 1, pcols, &
+ begchunk, endchunk, rdg_isovar, found, gridname='physgrid')
+ if (.not. found) call endrun(sub//': ERROR: ISOVAR not found on topo file')
+
+ call infld('ISOWGT', fh_topo, dim1name, dim2name, 1, pcols, &
+ begchunk, endchunk, rdg_isowgt, found, gridname='physgrid')
+ if (.not. found) call endrun(sub//': ERROR: ISOWGT not found on topo file')
+
call infld('HWDTH', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
1, prdg, begchunk, endchunk, rdg_hwdth, found, gridname='physgrid')
if (.not. found) call endrun(sub//': ERROR: HWDTH not found on topo file')
@@ -734,15 +758,39 @@ subroutine gw_init()
call addfld('ZMGW', (/ 'lev' /) , 'A' ,'m' , &
'midlayer geopotential heights in GW code ' )
+
+ call addfld('NIEGW', (/ 'ilev' /) , 'I' ,'1/s' , &
+ 'interface BV freq in GW code ' )
+ call addfld('NMEGW', (/ 'lev' /) , 'I' ,'1/s' , &
+ 'midlayer BV freq in GW code ' )
+ call addfld('RHOIEGW', (/ 'ilev' /) , 'I' ,'kg/m^3' , &
+ 'interface density in GW code ' )
+ call addfld('PINTEGW', (/ 'ilev' /) , 'I' ,'Pa' , &
+ 'interface density in GW code ' )
+
call addfld('TAUM1_DIAG' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')
call addfld('TAU1RDGBETAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
'Ridge based momentum flux profile')
- call addfld('UBM1BETA', (/ 'lev' /) , 'A' ,'s-1' , &
+ call addfld('UBM1BETA', (/ 'lev' /) , 'A' ,'m s-1' , &
'On-ridge wind profile ' )
- call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I' ,'m s-1' , &
+ call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I' ,'m s-2' , &
'On-ridge wind tendency from ridge 1 ')
+ call addfld('TAURESIDBETAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
+ 'Ridge based momentum flux profile')
+ call addfld('UBMRESIDBETA', (/ 'lev' /) , 'I' ,'m s-1' , &
+ 'On-ridge wind profile ' )
+ call addfld('UBIRESIDBETA', (/ 'ilev' /) , 'I' ,'m s-1' , &
+ 'On-ridge wind profile (interface) ' )
+ call addfld('SRC_LEVEL_RESIDBETA', horiz_only , 'I' ,'1' , &
+ 'src level index for ridge residual ' )
+ call addfld('TAUORO_RESID', horiz_only , 'I' ,'N m-2' , &
+ 'Surface mom flux from ridge reisdual ' )
+ call addfld('TAUDIAG_RESID' , (/ 'ilev' /) , 'I' ,'N m-2' , &
+ 'Ridge based momentum flux profile')
+
+
do i = 1, 6
write(cn, '(i1)') i
call addfld('TAU'//cn//'RDGBETAY' , (/ 'ilev' /), 'I', 'N m-2', &
@@ -1580,6 +1628,12 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
real(r8), pointer :: angll(:,:)
! anisotropy of ridges.
real(r8), pointer :: anixy(:,:)
+ ! sqrt(residual variance) not repr by ridges (assumed isotropic).
+ real(r8), pointer :: isovar(:)
+ ! area fraction of res variance
+ real(r8), pointer :: isowgt(:)
+
+
! Gamma ridges
! width of ridges.
@@ -2257,6 +2311,8 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
mxdis => rdg_mxdis(:ncol,:,lchnk)
angll => rdg_angll(:ncol,:,lchnk)
anixy => rdg_anixy(:ncol,:,lchnk)
+ isovar => rdg_isovar(:ncol,lchnk)
+ isowgt => rdg_isowgt(:ncol,lchnk)
where(mxdis < 0._r8)
mxdis = 0._r8
@@ -2276,6 +2332,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
nm, ni, rhoi, kvtt, q, dse, &
effgw_rdg_beta, effgw_rdg_beta_max, &
hwdth, clngt, gbxar, mxdis, angll, anixy, &
+ isovar, isowgt, &
rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
ptend, flx_heat)
@@ -2306,6 +2363,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
nm, ni, rhoi, kvtt, q, dse, &
effgw_rdg_gamma, effgw_rdg_gamma_max, &
hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, &
+ isovar, isowgt, &
rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, &
ptend, flx_heat)
@@ -2347,11 +2405,12 @@ subroutine gw_rdg_calc( &
effgw_rdg, effgw_rdg_max, &
hwdth, clngt, gbxar, &
mxdis, angll, anixy, &
+ isovar, isowgt, &
rdg_cd_llb, trpd_leewv, &
ptend, flx_heat)
use coords_1d, only: Coords1D
- use gw_rdg, only: gw_rdg_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
+ use gw_rdg, only: gw_rdg_src, gw_rdg_resid_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
use gw_common, only: gw_drag_prof, energy_change
character(len=5), intent(in) :: type ! BETA or GAMMA
@@ -2385,6 +2444,9 @@ subroutine gw_rdg_calc( &
real(r8), intent(in) :: angll(ncol,prdg) ! orientation of ridges.
real(r8), intent(in) :: anixy(ncol,prdg) ! Anisotropy parameter.
+ real(r8), intent(in) :: isovar(ncol) ! sqrt of residual variance
+ real(r8), intent(in) :: isowgt(ncol) ! area frac of residual variance
+
real(r8), intent(in) :: rdg_cd_llb ! Drag coefficient for low-level flow
logical, intent(in) :: trpd_leewv
@@ -2604,6 +2666,58 @@ subroutine gw_rdg_calc( &
end do ! end of loop over multiple ridges
+ ! Add additional GW from residual variance. Assumed isotropic
+ !kwvrdg = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!!
+ kwvrdg = 0.001_r8 / ( 100._r8 )
+ effgw = 1.0_r8 * isowgt
+ tauoro = 0._r8
+
+ call gw_rdg_resid_src(ncol, band_oro, p, &
+ u, v, t, isovar, kwvrdg, zi, nm, &
+ src_level, tend_level, tau, ubm, ubi, xv, yv, &
+ ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, phase_speeds, tauoro )
+
+ call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
+ t, vramp, &
+ piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
+ effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
+ ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
+ kwvrdg=kwvrdg, &
+ satfac_in = 1._r8, lapply_vdiff=gw_rdg_do_vdiff , tau_diag=tau_diag )
+
+ ! Add the tendencies from isotropic residual to the totals.
+ do k = 1, pver
+ ! diagnostics
+ utrdg(:,k) = utrdg(:,k) + utgw(:,k)
+ vtrdg(:,k) = vtrdg(:,k) + vtgw(:,k)
+ ttrdg(:,k) = ttrdg(:,k) + ttgw(:,k)
+ ! physics tendencies
+ ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
+ ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
+ ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
+ end do
+
+ do m = 1, pcnst
+ do k = 1, pver
+ ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
+ end do
+ end do
+
+ do k = 1, pver+1
+ taurx0(:,k) = tau(:,0,k)*xv
+ taury0(:,k) = tau(:,0,k)*yv
+ taurx(:,k) = taurx(:,k) + taurx0(:,k)
+ taury(:,k) = taury(:,k) + taury0(:,k)
+ end do
+
+ call outfld('TAUDIAG_RESID', tau_diag, ncol, lchnk)
+ call outfld('TAUORO_RESID', tauoro , ncol, lchnk)
+ call outfld('TAURESID'//trim(type)//'M', tau(:,0,:), ncol, lchnk)
+ call outfld('UBMRESID'//trim(type), ubm, ncol, lchnk)
+ call outfld('UBIRESID'//trim(type), ubi, ncol, lchnk)
+ call outfld('SRC_LEVEL_RESID'//trim(type), 1._r8*src_level , ncol, lchnk)
+ ! end of residual variance calc
+
! Calculate energy change for output to CAM's energy checker.
call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
ptend%v(:ncol,:), ptend%s(:ncol,:), de)
diff --git a/src/physics/cam/gw_rdg.F90 b/src/physics/cam/gw_rdg.F90
index b5a2a2137f..4c3e1ef745 100644
--- a/src/physics/cam/gw_rdg.F90
+++ b/src/physics/cam/gw_rdg.F90
@@ -19,6 +19,7 @@ module gw_rdg
! Public interface
public :: gw_rdg_readnl
public :: gw_rdg_src
+public :: gw_rdg_resid_src
public :: gw_rdg_belowpeak
public :: gw_rdg_break_trap
public :: gw_rdg_do_vdiff
@@ -175,6 +176,211 @@ subroutine gw_rdg_readnl(nlfile)
end subroutine gw_rdg_readnl
+!==========================================================================
+subroutine gw_rdg_resid_src(ncol, band, p, &
+ u, v, t, mxdis, kwvrdg, zi, nm, &
+ src_level, tend_level, tau, ubm, ubi, xv, yv, &
+ ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, c, tauoro )
+ use gw_common, only: rair, GWBand
+ use gw_utils, only: dot_2d, midpoint_interp, get_unit_vector
+ !-----------------------------------------------------------------------
+ ! 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--------------------------------
+ ! Column dimension.
+ integer, intent(in) :: ncol
+
+ ! Band to emit orographic waves in.
+ ! Regardless, we will only ever emit into l = 0.
+ type(GWBand), intent(in) :: band
+ ! Pressure coordinates.
+ type(Coords1D), intent(in) :: p
+
+
+ ! Midpoint zonal/meridional winds. ( m s-1)
+ real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
+ ! Midpoint temperatures. (K)
+ real(r8), intent(in) :: t(ncol,pver)
+ ! Height estimate for ridge (m) [anisotropic orography].
+ real(r8), intent(in) :: mxdis(ncol)
+ ! horiz wavenumber [anisotropic orography].
+ real(r8), intent(in) :: kwvrdg(ncol)
+ ! Interface altitudes above ground (m).
+ real(r8), intent(in) :: zi(ncol,pver+1)
+ ! Midpoint Brunt-Vaisalla frequencies (s-1).
+ real(r8), intent(in) :: nm(ncol,pver)
+
+ ! Indices of top gravity wave source level and lowest level where wind
+ ! tendencies are allowed.
+ integer, intent(out) :: src_level(ncol)
+ integer, intent(out) :: tend_level(ncol)
+
+ ! Averages over source region.
+ real(r8), intent(out) :: nsrc(ncol) ! B-V frequency.
+ real(r8), intent(out) :: rsrc(ncol) ! Density.
+ real(r8), intent(out) :: usrc(ncol) ! Zonal wind.
+ real(r8), intent(out) :: vsrc(ncol) ! Meridional wind.
+ real(r8), intent(out) :: ubmsrc(ncol) ! On-obstacle wind.
+ ! normalized wavenumber
+ real(r8), intent(out) :: m2src(ncol)
+
+
+ ! Wave Reynolds stress.
+ real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
+ ! Projection of wind at midpoints and interfaces.
+ real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
+ ! Unit vectors of source wind (zonal and meridional components).
+ real(r8), intent(out) :: xv(ncol), yv(ncol)
+ ! Phase speeds.
+ real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
+ ! source level mom. flux
+ real(r8), intent(out) :: tauoro(ncol)
+
+ !---------------------------Local Storage-------------------------------
+ ! Column and level indices.
+ integer :: i, k
+
+ ! Surface streamline displacement height (2*sgh).
+ real(r8) :: hdsp(ncol)
+
+ ! Difference in interface pressure across source region.
+ real(r8) :: dpsrc(ncol)
+ ! Thickness of downslope wind region.
+ real(r8) :: ddw(ncol)
+ ! Thickness of linear wave region.
+ real(r8) :: dwv(ncol)
+ ! Wind speed in source region.
+ real(r8) :: wmsrc(ncol)
+ ! source level mom. flux
+ !real(r8) :: tauoro(ncol)
+
+ real(r8) :: ragl(ncol)
+ real(r8) :: Fcrit_res,sghmax
+
+!--------------------------------------------------------------------------
+! Check that ngwav is equal to zero, otherwise end the job
+!--------------------------------------------------------------------------
+ if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not')
+
+!--------------------------------------------------------------------------
+! Average the basic state variables for the wave source over the depth of
+! the orographic standard deviation. Here we assume that the appropiate
+! values of wind, stability, etc. for determining the wave source are
+! averages over the depth of the atmosphere penterated by the typical
+! mountain.
+! Reduces to the bottom midpoint values when mxdis=0, such as over ocean.
+!--------------------------------------------------------------------------
+
+ Fcrit_res = 1.0_r8
+ hdsp = mxdis ! no longer multipied by 2
+ where(hdsp < 10._r8)
+ hdsp = 0._r8
+ end where
+
+ src_level = pver+1
+
+ tau(:,0,:) = 0.0_r8
+
+ ! Find depth of "source layer" for mountain waves
+ ! i.e., between ground and mountain top
+ do k = pver, 1, -1
+ do i = 1, ncol
+ ! Need to have h >= z(k+1) here or code will bomb when h=0.
+ if ( (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k)) ) then
+ src_level(i) = k
+ end if
+ end do
+ end do
+
+ rsrc = 0._r8
+ usrc = 0._r8
+ vsrc = 0._r8
+ nsrc = 0._r8
+ do i = 1, ncol
+ do k = pver, src_level(i), -1
+ rsrc(i) = rsrc(i) + p%mid(i,k) / (rair*t(i,k))* p%del(i,k)
+ usrc(i) = usrc(i) + u(i,k) * p%del(i,k)
+ vsrc(i) = vsrc(i) + v(i,k) * p%del(i,k)
+ nsrc(i) = nsrc(i) + nm(i,k)* p%del(i,k)
+ end do
+ end do
+
+
+ do i = 1, ncol
+ dpsrc(i) = p%ifc(i,pver+1) - p%ifc(i,src_level(i))
+ end do
+
+ rsrc = rsrc / dpsrc
+ usrc = usrc / dpsrc
+ vsrc = vsrc / dpsrc
+ nsrc = nsrc / dpsrc
+
+ ! Get the unit vector components and magnitude at the surface.
+ call get_unit_vector(usrc, vsrc, xv, yv, wmsrc )
+
+ ubmsrc = wmsrc
+
+ ! Project the local wind at midpoints onto the source wind.
+ do k = 1, pver
+ ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
+ end do
+
+ ! Compute the interface wind projection by averaging the midpoint winds.
+ ! Use the top level wind at the top interface.
+ ubi(:,1) = ubm(:,1)
+
+ ubi(:,2:pver) = midpoint_interp(ubm)
+
+ ! The minimum stratification allowing GW behavior
+ ! should really depend on horizontal scale since
+ !
+ ! m^2 ~ (N/U)^2 - k^2
+ !
+
+ m2src = ( (nsrc/(ubmsrc+0.01_r8))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01_r8))**2)
+
+ ! Compute the interface wind projection by averaging the midpoint winds.
+ ! Use the top level wind at the top interface.
+ ubi(:,1) = ubm(:,1)
+ ubi(:,2:pver) = midpoint_interp(ubm)
+ ubi(:,pver+1) = ubm(:,pver)
+
+
+
+ ! 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 ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
+ sghmax = Fcrit_res * (ubmsrc(i) / nsrc(i))**2
+ tauoro(i) = 0.5_r8 * kwvrdg(i) * min(hdsp(i)**2, sghmax) * &
+ rsrc(i) * nsrc(i) * ubmsrc(i)
+ else
+ tauoro(i) = 0._r8
+ end if
+ end do
+
+ do i = 1, ncol
+ do k=src_level(i),pver+1
+ tau(i,0,k) = tauoro(i)
+ end do
+ end do
+
+
+ ! Allow wind tendencies all the way to the model bottom.
+ tend_level = pver
+
+ ! No spectrum; phase speed is just 0.
+ c = 0._r8
+
+end subroutine gw_rdg_resid_src
+
+
+!==========================================================================
+
subroutine gw_rdg_src(ncol, band, p, &
u, v, t, mxdis, angxy, anixy, kwvrdg, iso, zi, nm, &
src_level, tend_level, bwv_level ,tlb_level , tau, ubm, ubi, xv, yv, &