#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 tau0contains!===============================================================================  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.