gw_drag.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!-----------------------------------------------------------------------
!------------------------------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 * (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./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 * (ni(i,k-1) + ni(i,k))
       end do
    end do

    return
  end subroutine gw_prof

!===============================================================================

  subroutine gw_oro (lchnk, ncol,                &
       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) :: 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.
! 
! Also determine the depth of the low level stress divergence region, as
! the max of the boundary layer depth and the source region depth. This
! can be done here if the stress magnitude does not determine the depth,
! otherwise it must be done below.
!---------------------------------------------------------------------------
    k = pver
    do i = 1, ncol
       ksrc(i) = k-1
       kldv(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 * 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
             kldv(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)
          elseif (pblh(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
             kldv(i) = k-1
          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))

       ubi(i,pver) = sqrt (usrc(i)**2 + vsrc(i)**2)
       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 * (ubm(i,k) + ubm(i,k+1))
       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.
          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
! +
! Find the top interface of the low level stress divergence region according
! to the maximum depth of three criterion.
! 1. source region depth
! 2. planetary boundary layer depth
! 3. 10 * (u_*) / N where u_* is defined from the gravity wave stresss
! = sqrt(tau/rho) using source region values
! -
!   if (kbot .lt. pver) then
!   do i = 1, ncol
!   kldv(i) = kbot
!   end do
!   else
!   do i = 1, ncol
!   zldv(i) = max (pblh(i), sgh(i)
!   zldv(i) = max (zdlv(i),
!   $           zldvcon * sqrt(tau(i,0,k)/rsrc(i)) / nsrc(i))
!   end do
!   kldv(i) = 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 do
!   end do
!   end if

! 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. / (pi(i,kldv(i)) - pi(i,pver))
       end if
    end do
    if (fracldv .le. 0.) kldvmn = pver

    return
  end subroutine gw_oro

!===============================================================================
  subroutine gw_bgnd (lchnk, ncol, &
       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*(u(i,k+1)+u(i,k))
       vsrc(i) = 0.5*(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)