gw_drag.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

module gw_drag

!---------------------------------------------------------------------------------
! 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 precision
  use ppgrid,         only: pcols, pver
  use physics_types,  only: physics_state, physics_ptend
  use pmgrid, only:         masterproc
  use history, only:        outfld

  implicit none

  save
  private                         ! Make default type private to the module
!
! PUBLIC: interfaces
!
  public gw_inti                  ! Initialization
  public gw_intr                  ! interface to actual parameterization
!
! PRIVATE: Rest of the data and interfaces are private to this module
!
  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               ! tendency efficiency
  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) :: 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) :: fcrit2              ! critical froude number
  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_inti (cpairx, cpwv, gx, rx, hypi)
!-----------------------------------------------------------------------
! Time independent initialization for multiple gravity wave parameterization.
!-----------------------------------------------------------------------
    use history,    only: addfld, add_default, phys_decomp

!------------------------------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

!-----------------------------------------------------------------------
! Copy model constants
    cpair  = cpairx
    g      = gx
    r      = rx

! Set MGWD constants
    effgw  = 0.125            ! efficiency of the tendencies
    kwv    = 6.28e-5          ! 100 km wave length
    dback  = 0.05             ! background diffusivity
    fcrit2 = 0.5              ! critical froude number squared
    tauscal= 0.001            ! scale factor for background stress
    taubgnd= 6.4              ! background stress amplitude

    fracldv= 0.0              ! fraction of tau0 diverged in low level region
    zldvcon= 10.              ! constant for determining zldv

! Set phase speeds 
    do k = -pgwv, pgwv
       c(k)   = 10. * k       ! 0, +/- 10, +/- 20, ... m/s
    end do

    if (masterproc) then
       write(6,*) ' '
       write(6,*) 'GW_INTI: pgwv = ', pgwv
       write(6,*) 'GW_INTI: c(l) = ', c
       write(6,*) ' '
    end if

! Set radiative damping times
    do k = 0, pver
       alpha(k) = 1.e-6       ! about 10 days.
    end do

! Min and max values to keep things reasonable
    mxasym = 0.1              ! max factor of 10 from |tau(c)| to |tau(-c)|
    mxrange= 0.001            ! factor of 100 from max to min |tau(c)|
    n2min  = 1.e-8            ! min value of Brunt-Vaisalla freq squared
    orohmin= 10.              ! min surface displacement for orographic wave drag
    orovmin=  2.              ! min wind speed for orographic wave drag
    taumin = 1.e-10           ! min stress considered > 0
    tndmax = 500. / 86400.    ! max permitted tendency (500 m/s/day)
    umcfac = 0.5              ! max permitted reduction in u-c
    ubmc2mn= 0.01             ! min value of (u-c)^2

! Determine other derived constants
    oroko2 = 0.5 * 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.) kbotbg  = k    ! spectrum source at 100 mb
!!$       if (hypi(k+1) .lt.  3000.) ktoporo = k
    end do
    ktoporo = 0

    if (masterproc) then
       write (6,*) 'KTOPBG  =',ktopbg
       write (6,*) 'KBOTBG  =',kbotbg
       write (6,*) 'KTOPORO =',ktoporo
       write (6,*) '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 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)

!-----------------------------------------------------------------------
! 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

    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         , &
            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.
       ptend%v = 0.
    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         , &
         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)
          ptend%v(i,k) = ptend%v(i,k) + vtgw(i,k)
          ptend%s(i,k) = -(ptend%u(i,k) * (state%u(i,k) + ptend%u(i,k)*0.5*dt) &
                          +ptend%v(i,k) * (state%v(i,k) + ptend%v(i,k)*0.5*dt))
          ttgw(i,k) = ptend%s(i,k) / cpair
       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)
!-----------------------------------------------------------------------
! 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.