gw_drag.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

       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 * (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.*atan(1.)
    al0=40.*pi_g/180.
    dlat0=10.*pi_g/180.
!
    call get_rlat_all_p(lchnk, ncol, rlat)
    do i=1,ncol
      flat_gw= 0.5*(1.+tanh((rlat(i)-al0)/dlat0)) + 0.5*(1.+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.)**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,                           &
       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, &
       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(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.
          vt(i,k) = 0.
       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.*ni(i,k)) )
             if (tausat(i,l) .le. taumin) tausat(i,l) = 0.0
             if (ubmc(i,l) / (ubi(i,k+1) - c(l)) .le. 0.0) tausat(i,l) = 0.0
             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. * kwv * ubmc2) * (alpha(k) + ni(i,k)**2/ubmc2 * d(i))
             taudmp = tau(i,l,k+1) * exp(-2.*mi*rog*t(i,k+1)*(piln(i,k+1)-piln(i,k)))
             if (taudmp .le. taumin) taudmp = 0.
             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. - 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 * (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
       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.*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