gw_drag_b.F90

Go to the documentation of this file.

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.2 (r3024) - 06/17/2009 13:03
!  
!   $Id: gw_drag_b.F90,v 1.2 2013/09/10 19:30:14 jgkim Exp $
MODULE gw_drag_b
  USE mapl_constantsmod, ONLY : mapl_p00, mapl_cp, mapl_grav, &
&  mapl_rgas, mapl_vireps, mapl_kappa
  USE gw_drag_init


  IMPLICIT NONE
!---------------------------------------------------------------------------------
! Purpose:
!
! Module to compute the forcing due to parameterized gravity waves. Both an 
! orographic and an internal source spectrum are considered.
!
! Author: Byron Boville
!         In-Sun Song
!
!---------------------------------------------------------------------------------
!save
! Make default type private to the module_src
  PRIVATE 
!
! PUBLIC: interfaces
!
! interface to actual parameterization
!  PUBLIC gw_intr, gw_sub_intr, gw_drag_prof, gw_bgnd, gw_oro, gw_prof
  PUBLIC gw_main_b
!, gw_drag_prof_b, gw_bgnd_b, gw_drag_prof_bgnd_b

!, gw_intr_fwd, gw_intr_bwd, gw_drag_prof_fwd, gw_drag_prof_bwd, &
!&  gw_oro_fwd, gw_oro_bwd, gw_prof_fwd, gw_prof_bwd
!  PUBLIC gw_main_b, gw_intr_fwd, gw_intr_bwd, gw_drag_prof_fwd, gw_drag_prof_bwd, &
!&  gw_bgnd_fwd, gw_bgnd_bwd, gw_oro_fwd, gw_oro_bwd
!
! PRIVATE: Rest of the data and interfaces are private to this module
!
! effective horizontal wave number for background
  REAL, PARAMETER, PUBLIC :: kwvb=6.28e-5
! effective horizontal wave number for background
  REAL, PARAMETER, PUBLIC :: kwvbeq=6.28e-5/7.
! effective horizontal wave number for orographic
  REAL, PARAMETER, PUBLIC :: kwvo=6.28e-5
! fraction of stress deposited in low level region
  REAL, PARAMETER, PUBLIC :: fracldv=0.0
! max asymmetry between tau(c) and tau(-c)
  REAL, PARAMETER, PUBLIC :: mxasym=0.1
! max range of tau for all c
  REAL, PARAMETER, PUBLIC :: mxrange=0.001
! min value of bouyancy frequency
  REAL, PARAMETER, PUBLIC :: n2min=1.e-8
! critical froude number
  REAL, PARAMETER, PUBLIC :: fcrit2=0.5
! min surface displacment height for orographic waves
  REAL, PARAMETER, PUBLIC :: orohmin=10.
! min wind speed for orographic waves
  REAL, PARAMETER, PUBLIC :: orovmin=2.0
! background source strength (/TAUSCAL)
  REAL, PARAMETER, PUBLIC :: taubgnd=6.4
! minimum (nonzero) stress
  REAL, PARAMETER, PUBLIC :: taumin=1.e-10
! scale factor for background stress source
  REAL, PARAMETER, PUBLIC :: tauscal=0.001
! maximum wind tendency
  REAL, PARAMETER, PUBLIC :: tndmax=500./86400.
! factor to limit tendency to prevent reversing u-c
  REAL, PARAMETER, PUBLIC :: umcfac=0.5
! min (u-c)**2
  REAL, PARAMETER, PUBLIC :: ubmc2mn=0.1
! constant for determining zldv from tau0
  REAL, PARAMETER, PUBLIC :: zldvcon=10.
  REAL, PARAMETER, PUBLIC :: rog=mapl_rgas/mapl_grav
! 1/2 * horizontal wavenumber
  REAL, PARAMETER, PUBLIC :: oroko2=0.5*kwvo
! This is *not* MAPL_PI
  REAL, PARAMETER, PUBLIC :: pi_gwd=4.0*atan(1.0)
! -----------------------------------------------

CONTAINS

SUBROUTINE gw_main_b(pcols, pver, dt, pgwv, effgworo_dev, &
&  effgwbkg_dev, pint_dev, t_dev, u_dev, u_devb, v_dev, v_devb, &
&  sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, &
&  qvt_dev, rog, mapl_vireps_, rlat_dev)
  IMPLICIT NONE
!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! - number of columns                                            
  INTEGER :: pcols
! - number of vertical layers
  INTEGER :: pver
! - time step
  REAL :: dt
! - number of waves allowed                (Default = 4, 0 nullifies)
  INTEGER :: pgwv
! - tendency efficiency for background gwd (Default = 0.125)
  REAL :: effgwbkg_dev(pcols)
! - tendency efficiency for orographic gwd (Default = 0.125)
  REAL :: effgworo_dev(pcols)
! - pressure at the layer edges
  REAL :: pint_dev(pcols, pver+1)
! - temperature at layers
  REAL :: t_dev(pcols, pver)
! - zonal wind at layers
  REAL :: u_dev(pcols, pver)
  REAL :: u_devb(pcols, pver)
! - meridional wind at layers
  REAL :: v_dev(pcols, pver)
  REAL :: v_devb(pcols, pver)
! e standard deviation of orography
  REAL :: sgh_dev(pcols)
! e reference pressure at the layeredges
  REAL :: pref_dev(pver+1)
! c pressure at the layers
  REAL :: pmid_dev(pcols, pver)
! c pressure thickness at the layers
  REAL :: pdel_dev(pcols, pver)
! c 1.0 / pdel
  REAL :: rpdel_dev(pcols, pver)
! c log(pint)
  REAL :: lnpint_dev(pcols, pver+1)
! c height above surface at layers
  REAL :: zm_dev(pcols, pver)
! c height above surface at layers
  REAL :: zi_dev(pcols, pver+1)
! c pressure at the layers
  REAL :: qvt_dev(pcols, pver)
!mg latitude in radian
  REAL :: rlat_dev(pcols)
! zonal wind tendency at layer 
  REAL :: dudt_gwd_dev(pcols, pver)
  REAL :: dudt_gwd_devb(pcols, pver)
! meridional wind tendency at layer 
  REAL :: dvdt_gwd_dev(pcols, pver)
  REAL :: dvdt_gwd_devb(pcols, pver)
! temperature tendency at layer
  REAL :: dtdt_gwd_dev(pcols, pver)
  REAL :: dtdt_gwd_devb(pcols, pver)
! zonal wind tendency at layer due to orography GWD
  REAL :: dudt_org_dev(pcols, pver)
! meridional wind tendency at layer  due to orography GWD
  REAL :: dvdt_org_dev(pcols, pver)
! temperature tendency at layer  due to orography GWD
  REAL :: dtdt_org_dev(pcols, pver)
! zonal      gravity wave surface    stress
  REAL :: taugwdx_dev(pcols)
! meridional gravity wave surface    stress
  REAL :: taugwdy_dev(pcols)
! zonal      orographic gravity wave stress
  REAL :: tauox_dev(pcols, pver+1)
! meridional orographic gravity wave stress
  REAL :: tauoy_dev(pcols, pver+1)
! energy flux of orographic gravity waves
  REAL :: feo_dev(pcols, pver+1)
! pseudoenergy flux of orographic gravity waves
  REAL :: fepo_dev(pcols, pver+1)
! zonal      gravity wave background stress
  REAL :: taubkgx_dev(pcols)
! meridional gravity wave background stress
  REAL :: taubkgy_dev(pcols)
! zonal      background gravity wave stress
  REAL :: taubx_dev(pcols, pver+1)
! meridional background gravity wave stress
  REAL :: tauby_dev(pcols, pver+1)
! energy flux of background gravity waves
  REAL :: feb_dev(pcols, pver+1)
! pseudoenergy flux of background gravity waves
  REAL :: fepb_dev(pcols, pver+1)
! dU/dt below background launch level
  REAL :: utbsrc_dev(pcols, pver)
! dV/dt below background launch level
  REAL :: vtbsrc_dev(pcols, pver)
! dT/dt below background launch level
  REAL :: ttbsrc_dev(pcols, pver)
!---------------------------Local storage-------------------------------
! loop indexes
  INTEGER :: i, k, kc
! launch-level index for orographic
  INTEGER :: kbotoro
! launch-level index for background
  INTEGER :: kbotbg
! top interface of gwd region
  INTEGER :: ktopbg, ktoporo
! top interface of low level stress divergence region
  INTEGER :: kldv
! min value of kldv
  INTEGER :: kldvmn
! index of top interface of source region
  INTEGER :: ksrc
! min value of ksrc
  INTEGER :: ksrcmn
! temperature tendency
  REAL :: ttgw(pver)
! zonal wind tendency
  REAL :: utgw(pver)
! meridional wind tendency
  REAL :: vtgw(pver)
! interface Brunt-Vaisalla frequency
  REAL :: ni(0:pver)
! midpoint Brunt-Vaisalla frequency
  REAL :: nm(pver)
! 1/dp across low level divergence region
  REAL :: rdpldv
! interface density
  REAL :: rhoi(0:pver)
! c=0 sfc. stress (zonal)
  REAL :: tau0x
! c=0 sfc. stress (meridional)
  REAL :: tau0y
! interface temperature
  REAL :: ti(0:pver)
! projection of wind at interfaces
  REAL :: ubi(0:pver)
! projection of wind at midpoints
  REAL :: ubm(pver)
! unit vectors of source wind (x)
  REAL :: xv
! unit vectors of source wind (y)
  REAL :: yv
  REAL :: utosrc(pver)
  REAL :: vtosrc(pver)
  REAL :: ttosrc(pver)
! newtonian cooling coefficients
  REAL :: alpha(0:pver)
! newtonian cooling coefficients
  REAL :: dback(0:pver)
! wave phase speeds
  REAL :: c(-pgwv:pgwv)
! - temperature at layers
  REAL :: att_dev(pcols, pver)
  REAL :: att_devb(pcols, pver)
  REAL :: pwr1
  REAL :: pwx2
  REAL :: pwy2
  REAL :: pwr2
  REAL :: hkl, hkk, tvfac, tv, rog
  REAL :: mapl_vireps_

  isize_s1_tau=-pgwv;isize_s2_tau=0;isize_e1_tau=pgwv;isize_e2_tau=pver

  DO i=1,pcols
    DO k=1,pver
      att_dev(i, k) = t_dev(i, k)*mapl_p00**mapl_kappa*(mapl_p00/&
&        pint_dev(i, k))**(-mapl_kappa)
    END DO
  END DO
  DO i=1,pcols
    zi_dev(i, pver+1) = 0.
    DO k=pver,1,-1
      hkl = lnpint_dev(i, k+1) - lnpint_dev(i, k)
      hkk = 1. - pint_dev(i, k)*hkl*rpdel_dev(i, k)
      tvfac = 1. + mapl_vireps_*qvt_dev(i, k)
      tv = att_dev(i, k)*tvfac
      zm_dev(i, k) = zi_dev(i, k+1) + rog*tv*hkk
      zi_dev(i, k) = zi_dev(i, k+1) + rog*tv*hkl
    END DO
  END DO

  DO i=1,pcols
    If (.not.bypass_bgnd) then
    do_bgnd=.true.
    do_oro=.false.

    CALL gw_intr_fwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&               effgwbkg_dev, pint_dev, att_dev, u_dev, v_dev, sgh_dev, &
&               pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, &
&               zm_dev, rlat_dev, dudt_gwd_dev, dvdt_gwd_dev, &
&               dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, &
&               taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, &
&               taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, &
&               fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    Endif

    If (.not.bypass_oro) then
    do_bgnd=.false.
    do_oro=.true.

    CALL gw_intr_fwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&               effgwbkg_dev, pint_dev, att_dev, u_dev, v_dev, sgh_dev, &
&               pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, &
&               zm_dev, rlat_dev, dudt_gwd_dev, dvdt_gwd_dev, &
&               dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, &
&               taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, &
&               taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, &
&               fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    Endif

    If (.not.bypass_oro) then
    dvdt_gwd_devb = 0.0
    dudt_gwd_devb = 0.0

    do_bgnd=.false.
    do_oro=.true.

    CALL gw_intr_bwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&               effgwbkg_dev, pint_dev, att_dev, u_dev, u_devb&
&               , v_dev, v_devb, sgh_dev, pref_dev, pmid_dev, pdel_dev, &
&               rpdel_dev, lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, &
&               dudt_gwd_devb, dvdt_gwd_dev, dvdt_gwd_devb, dtdt_gwd_dev&
&               , dudt_org_dev, dvdt_org_dev, dtdt_org_dev&
&               , taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev&
&               , taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev&
&               , fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    Endif

    If (.not.bypass_bgnd) then
    dvdt_gwd_devb = 0.0
    dudt_gwd_devb = 0.0

    do_bgnd=.true.
    do_oro=.false.

    CALL gw_intr_bwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&               effgwbkg_dev, pint_dev, att_dev, u_dev, u_devb&
&               , v_dev, v_devb, sgh_dev, pref_dev, pmid_dev, pdel_dev, &
&               rpdel_dev, lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, &
&               dudt_gwd_devb, dvdt_gwd_dev, dvdt_gwd_devb, dtdt_gwd_dev&
&               , dudt_org_dev, dvdt_org_dev, dtdt_org_dev&
&               , taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev&
&               , taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev&
&               , fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    Endif
  END DO

END SUBROUTINE gw_main_b

  SUBROUTINE gw_intr_fwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&    effgwbkg_dev, pint_dev, t_dev, u_dev, v_dev, sgh_dev, pref_dev, &
&    pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, &
&    dudt_gwd_dev, dvdt_gwd_dev, dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev&
&    , dtdt_org_dev, taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, &
&    feo_dev, taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, &
&    fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! - number of columns                                            
    INTEGER, INTENT(IN) :: pcols
! - number of vertical layers
    INTEGER, INTENT(IN) :: pver
! - time step
    REAL, INTENT(IN) :: dt
! - number of waves allowed                (Default = 4, 0 nullifies)
    INTEGER, INTENT(IN) :: pgwv
! - tendency efficiency for background gwd (Default = 0.125)
    REAL :: effgwbkg_dev(pcols)
! - tendency efficiency for orographic gwd (Default = 0.125)
    REAL :: effgworo_dev(pcols)
! - pressure at the layer edges
    REAL :: pint_dev(pcols, pver+1)
! - temperature at layers
    REAL :: t_dev(pcols, pver)
! - zonal wind at layers
    REAL :: u_dev(pcols, pver)
! - meridional wind at layers
    REAL :: v_dev(pcols, pver)
! e standard deviation of orography
    REAL :: sgh_dev(pcols)
! e reference pressure at the layeredges
    REAL :: pref_dev(pver+1)
! c pressure at the layers
    REAL :: pmid_dev(pcols, pver)
! c pressure thickness at the layers
    REAL :: pdel_dev(pcols, pver)
! c 1.0 / pdel
    REAL :: rpdel_dev(pcols, pver)
! c log(pint)
    REAL :: lnpint_dev(pcols, pver+1)
! c height above surface at layers
    REAL :: zm_dev(pcols, pver)
!mg latitude in radian
    REAL :: rlat_dev(pcols)
! zonal wind tendency at layer 
    REAL :: dudt_gwd_dev(pcols, pver)
! meridional wind tendency at layer 
    REAL :: dvdt_gwd_dev(pcols, pver)
! temperature tendency at layer
    REAL :: dtdt_gwd_dev(pcols, pver)
! zonal wind tendency at layer due to orography GWD
    REAL :: dudt_org_dev(pcols, pver)
! meridional wind tendency at layer  due to orography GWD
    REAL :: dvdt_org_dev(pcols, pver)
! temperature tendency at layer  due to orography GWD
    REAL :: dtdt_org_dev(pcols, pver)
! zonal      gravity wave surface    stress
    REAL :: taugwdx_dev(pcols)
! meridional gravity wave surface    stress
    REAL :: taugwdy_dev(pcols)
! zonal      orographic gravity wave stress
    REAL :: tauox_dev(pcols, pver+1)
! meridional orographic gravity wave stress
    REAL :: tauoy_dev(pcols, pver+1)
! energy flux of orographic gravity waves
    REAL :: feo_dev(pcols, pver+1)
! pseudoenergy flux of orographic gravity waves
    REAL :: fepo_dev(pcols, pver+1)
! zonal      gravity wave background stress
    REAL :: taubkgx_dev(pcols)
! meridional gravity wave background stress
    REAL :: taubkgy_dev(pcols)
! zonal      background gravity wave stress
    REAL :: taubx_dev(pcols, pver+1)
! meridional background gravity wave stress
    REAL :: tauby_dev(pcols, pver+1)
! energy flux of background gravity waves
    REAL :: feb_dev(pcols, pver+1)
! pseudoenergy flux of background gravity waves
    REAL :: fepb_dev(pcols, pver+1)
! dU/dt below background launch level
    REAL :: utbsrc_dev(pcols, pver)
! dV/dt below background launch level
    REAL :: vtbsrc_dev(pcols, pver)
! dT/dt below background launch level
    REAL :: ttbsrc_dev(pcols, pver)
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: i, ii, k, kc
! launch-level index for orographic
    INTEGER :: kbotoro
! launch-level index for background
    INTEGER :: kbotbg
! top interface of gwd region
    INTEGER :: ktopbg, ktoporo
! top interface of low level stress divergence region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! temperature tendency
    REAL :: ttgw(pver)
! zonal wind tendency
    REAL :: utgw(pver)
! meridional wind tendency
    REAL :: vtgw(pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! interface density
    REAL :: rhoi(0:pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
! interface temperature
    REAL :: ti(0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
    REAL :: utosrc(pver)
    REAL :: vtosrc(pver)
    REAL :: ttosrc(pver)
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! Assign wave phase speeds
! ------------------------
    DO kc=-pgwv,pgwv
      c(kc) = 10.0*kc
    END DO
! Determine the bounds of the background and orographic stress regions
    ktoporo = 0
    ktopbg = 0
    kbotoro = pver
    DO k=0,pver
      IF (pref_dev(k+1) .LT. p_src) kbotbg = k
! spectrum source at 400 mb
    END DO
    DO k=0,pver
! Profiles of background state variables
      CALL gw_prof(i, k, pcols, pver, u_dev, v_dev, t_dev, pmid_dev, &
&             pint_dev, rhoi, ni, ti, nm)
    END DO
!-----------------------------------------------------------------------------
! Non-orographic backgound gravity wave spectrum
!-----------------------------------------------------------------------------
    IF (pgwv .GT. 0 .AND. do_bgnd) THEN
! Determine the wave source for a background spectrum at ~400 mb
      CALL gw_bgnd_fwd(i, pcols, pver, c, u_dev, v_dev, t_dev, pmid_dev&
&                 , pint_dev, pdel_dev, rpdel_dev, lnpint_dev, rlat_dev, &
&                 kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, &
&                 yv, pgwv, kbotbg)
      CALL pushreal4array(vtosrc, pver)
      CALL pushreal4array(utosrc, pver)
      CALL pushreal4array(vtgw, pver)
      CALL pushreal4array(utgw, pver)
      CALL pushreal4array(tau, (2*pgwv+1)*(pver+1))
! Solve for the drag profile
      CALL gw_drag_prof_bgnd(i, pcols, pver, pgwv, pgwv, kbotbg, ktopbg&
&                       , c, u_dev, v_dev, t_dev, pint_dev, pdel_dev, &
&                       rpdel_dev, lnpint_dev, rlat_dev, rhoi, ni, ti, nm&
&                       , dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, &
&                       rdpldv, tau, ubi, ubm, xv, yv, utgw, vtgw, ttgw, &
&                       taubx_dev, tauby_dev, feb_dev, fepb_dev, utosrc, &
&                       vtosrc, ttosrc, tau0x, tau0y, effgwbkg_dev)
! Add the momentum tendencies to the output tendency arrays
      DO k=1,pver
        dudt_gwd_dev(i, k) = utgw(k) + utosrc(k)
        dvdt_gwd_dev(i, k) = vtgw(k) + vtosrc(k)
      END DO
      CALL pushinteger4(0)
    ELSE
!             dudt_gwd_dev(i,k) = utgw(k) + utosrc(k)
!             dvdt_gwd_dev(i,k) = vtgw(k) + vtosrc(k)
!             dtdt_gwd_dev(i,k) = ttgw(k) + ttosrc(k)
! zero net tendencies if no spectrum computed
      DO k=1,pver
        dudt_gwd_dev(i, k) = 0.
        dvdt_gwd_dev(i, k) = 0.
      END DO
      CALL pushinteger4(1)
    END IF
!            dtdt_gwd_dev(i,k) = 0.
!-----------------------------------------------------------------------------
! Orographic stationary gravity wave
!-----------------------------------------------------------------------------
! Determine the orographic wave source
    IF (do_oro) THEN
      CALL gw_oro_fwd(i, pcols, pver, pgwv, u_dev, v_dev, t_dev, sgh_dev&
&                , pmid_dev, pint_dev, pdel_dev, zm_dev, nm, kldv, kldvmn&
&                , ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, kbotoro, &
&                rlat_dev)
! Solve for the drag profile
      CALL gw_drag_prof_fwd(i, pcols, pver, pgwv, 0, kbotoro, ktoporo, c&
&                      , u_dev, v_dev, t_dev, pint_dev, pdel_dev, &
&                      rpdel_dev, lnpint_dev, rlat_dev, rhoi, ni, ti, nm&
&                      , dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, &
&                      rdpldv, tau, ubi, ubm, xv, yv, utgw, vtgw, ttgw, &
&                      tauox_dev, tauoy_dev, feo_dev, fepo_dev, utosrc, &
&                      vtosrc, ttosrc, tau0x, tau0y, effgworo_dev)
! Add the orographic tendencies to the spectrum tendencies
! Compute the temperature tendency from energy conservation (includes spectrum).
      DO k=1,pver
        dudt_gwd_dev(i, k) = dudt_gwd_dev(i, k) + utgw(k)
        dvdt_gwd_dev(i, k) = dvdt_gwd_dev(i, k) + vtgw(k)
      END DO
      CALL pushinteger4(1)
    ELSE
      CALL pushinteger4(0)
    END IF
!          dtdt_gwd_dev(i,k) = dtdt_gwd_dev(i,k) + ttgw(k)
    DO k=1,pver
      CALL pushreal4(u_dev(i, k))
      u_dev(i, k) = u_dev(i, k) + dudt_gwd_dev(i, k)*dt
      CALL pushreal4(v_dev(i, k))
      v_dev(i, k) = v_dev(i, k) + dvdt_gwd_dev(i, k)*dt
    END DO
    CALL pushinteger4(ktoporo)
    CALL pushreal4array(ni, pver + 1)
    CALL pushreal4array(vtgw, pver)
    CALL pushreal4array(utgw, pver)
    CALL pushinteger4(kbotbg)
    CALL pushreal4(xv)
    CALL pushinteger4(kbotoro)
    CALL pushreal4array(rhoi, pver + 1)
    CALL pushreal4(yv)
    CALL pushinteger4(ktopbg)
    CALL pushreal4array(ubi, pver + 1)
    CALL pushreal4array(c, 2*pgwv + 1)
  END SUBROUTINE gw_intr_fwd
!  Differentiation of gw_intr in reverse (adjoint) mode (backward sweep):
!   gradient, with respect to input variables: dvdt_gwd_dev u_dev
!                dudt_gwd_dev v_dev
!   of linear combination of output variables: dvdt_gwd_dev u_dev
!                dudt_gwd_dev v_dev
  SUBROUTINE gw_intr_bwd(i, pcols, pver, dt, pgwv, effgworo_dev, &
&    effgwbkg_dev, pint_dev, t_dev, u_dev, u_devb, v_dev, v_devb, sgh_dev&
&    , pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, &
&    rlat_dev, dudt_gwd_dev, dudt_gwd_devb, dvdt_gwd_dev, dvdt_gwd_devb, &
&    dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev&
&    , taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, &
&    taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, fepb_dev, &
&    utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: pcols
    INTEGER, INTENT(IN) :: pver
    REAL, INTENT(IN) :: dt
    INTEGER, INTENT(IN) :: pgwv
    REAL :: effgwbkg_dev(pcols)
    REAL :: effgworo_dev(pcols)
    REAL :: pint_dev(pcols, pver+1)
    REAL :: t_dev(pcols, pver)
    REAL :: u_dev(pcols, pver)
    REAL :: u_devb(pcols, pver)
    REAL :: v_dev(pcols, pver)
    REAL :: v_devb(pcols, pver)
    REAL :: sgh_dev(pcols)
    REAL :: pref_dev(pver+1)
    REAL :: pmid_dev(pcols, pver)
    REAL :: pdel_dev(pcols, pver)
    REAL :: rpdel_dev(pcols, pver)
    REAL :: lnpint_dev(pcols, pver+1)
    REAL :: zm_dev(pcols, pver)
    REAL :: rlat_dev(pcols)
    REAL :: dudt_gwd_dev(pcols, pver)
    REAL :: dudt_gwd_devb(pcols, pver)
    REAL :: dvdt_gwd_dev(pcols, pver)
    REAL :: dvdt_gwd_devb(pcols, pver)
    REAL :: dtdt_gwd_dev(pcols, pver)
    REAL :: dudt_org_dev(pcols, pver)
    REAL :: dvdt_org_dev(pcols, pver)
    REAL :: dtdt_org_dev(pcols, pver)
    REAL :: taugwdx_dev(pcols)
    REAL :: taugwdy_dev(pcols)
    REAL :: tauox_dev(pcols, pver+1)
    REAL :: tauoy_dev(pcols, pver+1)
    REAL :: feo_dev(pcols, pver+1)
    REAL :: fepo_dev(pcols, pver+1)
    REAL :: taubkgx_dev(pcols)
    REAL :: taubkgy_dev(pcols)
    REAL :: taubx_dev(pcols, pver+1)
    REAL :: tauby_dev(pcols, pver+1)
    REAL :: feb_dev(pcols, pver+1)
    REAL :: fepb_dev(pcols, pver+1)
    REAL :: utbsrc_dev(pcols, pver)
    REAL :: vtbsrc_dev(pcols, pver)
    REAL :: ttbsrc_dev(pcols, pver)
    INTEGER :: i, ii, k, kc
    INTEGER :: kbotoro
    INTEGER :: kbotbg
    INTEGER :: ktopbg, ktoporo
    INTEGER :: kldv
    INTEGER :: kldvmn
    INTEGER :: ksrc
    INTEGER :: ksrcmn
    REAL :: ttgw(pver)
    REAL :: utgw(pver)
    REAL :: utgwb(pver)
    REAL :: vtgw(pver)
    REAL :: vtgwb(pver)
    REAL :: ni(0:pver)
    REAL :: nm(pver)
    REAL :: rdpldv
    REAL :: rhoi(0:pver)
    REAL :: tau0x
    REAL :: tau0y
    REAL :: ti(0:pver)
    REAL :: ubi(0:pver)
    REAL :: ubib(0:pver)
    REAL :: ubm(pver)
    REAL :: ubmb(pver)
    REAL :: xv
    REAL :: xvb
    REAL :: yv
    REAL :: yvb
    REAL :: utosrc(pver)
    REAL :: utosrcb(pver)
    REAL :: vtosrc(pver)
    REAL :: vtosrcb(pver)
    REAL :: ttosrc(pver)
    REAL :: alpha(0:pver)
    REAL :: dback(0:pver)
    REAL :: c(-pgwv:pgwv)
    REAL :: tau(-pgwv:pgwv, 0:pver)
    REAL :: taub(-pgwv:pgwv, 0:pver)
    INTEGER :: branch
    CALL popreal4array(c, 2*pgwv + 1)
    CALL popreal4array(ubi, pver + 1)
    CALL popinteger4(ktopbg)
    CALL popreal4(yv)
    CALL popreal4array(rhoi, pver + 1)
    CALL popinteger4(kbotoro)
    CALL popreal4(xv)
    CALL popinteger4(kbotbg)
    CALL popreal4array(utgw, pver)
    CALL popreal4array(vtgw, pver)
    CALL popreal4array(ni, pver + 1)
    CALL popinteger4(ktoporo)
    DO k=pver,1,-1
      CALL popreal4(v_dev(i, k))
      dvdt_gwd_devb(i, k) = dvdt_gwd_devb(i, k) + dt*v_devb(i, k)
      CALL popreal4(u_dev(i, k))
      dudt_gwd_devb(i, k) = dudt_gwd_devb(i, k) + dt*u_devb(i, k)
    END DO
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      taub(:, :) = 0.0_8
      vtgwb(:) = 0.0_8
      utgwb(:) = 0.0_8
      ubib(:) = 0.0_8
      ubmb(:) = 0.0_8
    ELSE
      vtgwb(:) = 0.0_8
      utgwb(:) = 0.0_8
      DO k=pver,1,-1
        vtgwb(k) = vtgwb(k) + dvdt_gwd_devb(i, k)
        utgwb(k) = utgwb(k) + dudt_gwd_devb(i, k)
      END DO
      CALL gw_drag_prof_bwd(i, pcols, pver, pgwv, 0, kbotoro, ktoporo, c&
&                      , u_dev, v_dev, t_dev, pint_dev, pdel_dev, &
&                      rpdel_dev, lnpint_dev, rlat_dev, rhoi, ni, ti, nm&
&                      , dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, &
&                      rdpldv, tau, taub, ubi, ubib, ubm, xv, xvb, yv, &
&                      yvb, utgw, utgwb, vtgw, vtgwb, ttgw, tauox_dev, &
&                      tauoy_dev, feo_dev, fepo_dev, utosrc, vtosrc, &
&                      ttosrc, tau0x, tau0y, effgworo_dev)
      CALL gw_oro_bwd(i, pcols, pver, pgwv, u_dev, u_devb, v_dev, v_devb&
&                , t_dev, sgh_dev, pmid_dev, pint_dev, pdel_dev, zm_dev, &
&                nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taub, ubi, &
&                ubib, ubm, ubmb, xv, xvb, yv, yvb, kbotoro, rlat_dev)
    END IF
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      vtosrcb(:) = 0.0_8
      utosrcb(:) = 0.0_8
      DO k=pver,1,-1
        vtgwb(k) = vtgwb(k) + dvdt_gwd_devb(i, k)
        vtosrcb(k) = vtosrcb(k) + dvdt_gwd_devb(i, k)
        dvdt_gwd_devb(i, k) = 0.0_8
        utgwb(k) = utgwb(k) + dudt_gwd_devb(i, k)
        utosrcb(k) = utosrcb(k) + dudt_gwd_devb(i, k)
        dudt_gwd_devb(i, k) = 0.0_8
      END DO
      CALL popreal4array(tau, (2*pgwv+1)*(pver+1))
      CALL popreal4array(utgw, pver)
      CALL popreal4array(vtgw, pver)
      CALL popreal4array(utosrc, pver)
      CALL popreal4array(vtosrc, pver)
      CALL gw_drag_prof_bgnd_b(i, pcols, pver, pgwv, pgwv, kbotbg, &
&                         ktopbg, c, u_dev, v_dev, t_dev, pint_dev, &
&                         pdel_dev, rpdel_dev, lnpint_dev, rlat_dev, rhoi&
&                         , ni, ti, nm, dt, alpha, dback, kldv, kldvmn, &
&                         ksrc, ksrcmn, rdpldv, tau, taub, ubi, ubib, ubm&
&                         , xv, xvb, yv, yvb, utgw, utgwb, vtgw, vtgwb, &
&                         ttgw, taubx_dev, tauby_dev, feb_dev, fepb_dev, &
&                         utosrc, utosrcb, vtosrc, vtosrcb, ttosrc, tau0x&
&                         , tau0y, effgwbkg_dev)
      CALL gw_bgnd_bwd(i, pcols, pver, c, u_dev, u_devb, v_dev, v_devb, &
&                 t_dev, pmid_dev, pint_dev, pdel_dev, rpdel_dev, &
&                 lnpint_dev, rlat_dev, kldv, kldvmn, ksrc, ksrcmn, &
&                 rdpldv, tau, ubi, ubib, ubm, ubmb, xv, xvb, yv, yvb, &
&                 pgwv, kbotbg)
    ELSE
      DO k=pver,1,-1
        dvdt_gwd_devb(i, k) = 0.0_8
        dudt_gwd_devb(i, k) = 0.0_8
      END DO
    END IF
  END SUBROUTINE gw_intr_bwd
  SUBROUTINE gw_intr(i, pcols, pver, dt, pgwv, effgworo_dev, &
&    effgwbkg_dev, pint_dev, t_dev, u_dev, v_dev, sgh_dev, pref_dev, &
&    pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, &
&    dudt_gwd_dev, dvdt_gwd_dev, dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev&
&    , dtdt_org_dev, taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, &
&    feo_dev, taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, &
&    fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! - number of columns                                            
    INTEGER, INTENT(IN) :: pcols
! - number of vertical layers
    INTEGER, INTENT(IN) :: pver
! - time step
    REAL, INTENT(IN) :: dt
! - number of waves allowed                (Default = 4, 0 nullifies)
    INTEGER, INTENT(IN) :: pgwv
! - tendency efficiency for background gwd (Default = 0.125)
    REAL :: effgwbkg_dev(pcols)
! - tendency efficiency for orographic gwd (Default = 0.125)
    REAL :: effgworo_dev(pcols)
! - pressure at the layer edges
    REAL :: pint_dev(pcols, pver+1)
! - temperature at layers
    REAL :: t_dev(pcols, pver)
! - zonal wind at layers
    REAL :: u_dev(pcols, pver)
! - meridional wind at layers
    REAL :: v_dev(pcols, pver)
! e standard deviation of orography
    REAL :: sgh_dev(pcols)
! e reference pressure at the layeredges
    REAL :: pref_dev(pver+1)
! c pressure at the layers
    REAL :: pmid_dev(pcols, pver)
! c pressure thickness at the layers
    REAL :: pdel_dev(pcols, pver)
! c 1.0 / pdel
    REAL :: rpdel_dev(pcols, pver)
! c log(pint)
    REAL :: lnpint_dev(pcols, pver+1)
! c height above surface at layers
    REAL :: zm_dev(pcols, pver)
!mg latitude in radian
    REAL :: rlat_dev(pcols)
! zonal wind tendency at layer 
    REAL :: dudt_gwd_dev(pcols, pver)
! meridional wind tendency at layer 
    REAL :: dvdt_gwd_dev(pcols, pver)
! temperature tendency at layer
    REAL :: dtdt_gwd_dev(pcols, pver)
! zonal wind tendency at layer due to orography GWD
    REAL :: dudt_org_dev(pcols, pver)
! meridional wind tendency at layer  due to orography GWD
    REAL :: dvdt_org_dev(pcols, pver)
! temperature tendency at layer  due to orography GWD
    REAL :: dtdt_org_dev(pcols, pver)
! zonal      gravity wave surface    stress
    REAL :: taugwdx_dev(pcols)
! meridional gravity wave surface    stress
    REAL :: taugwdy_dev(pcols)
! zonal      orographic gravity wave stress
    REAL :: tauox_dev(pcols, pver+1)
! meridional orographic gravity wave stress
    REAL :: tauoy_dev(pcols, pver+1)
! energy flux of orographic gravity waves
    REAL :: feo_dev(pcols, pver+1)
! pseudoenergy flux of orographic gravity waves
    REAL :: fepo_dev(pcols, pver+1)
! zonal      gravity wave background stress
    REAL :: taubkgx_dev(pcols)
! meridional gravity wave background stress
    REAL :: taubkgy_dev(pcols)
! zonal      background gravity wave stress
    REAL :: taubx_dev(pcols, pver+1)
! meridional background gravity wave stress
    REAL :: tauby_dev(pcols, pver+1)
! energy flux of background gravity waves
    REAL :: feb_dev(pcols, pver+1)
! pseudoenergy flux of background gravity waves
    REAL :: fepb_dev(pcols, pver+1)
! dU/dt below background launch level
    REAL :: utbsrc_dev(pcols, pver)
! dV/dt below background launch level
    REAL :: vtbsrc_dev(pcols, pver)
! dT/dt below background launch level
    REAL :: ttbsrc_dev(pcols, pver)
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: i, ii, k, kc
! launch-level index for orographic
    INTEGER :: kbotoro
! launch-level index for background
    INTEGER :: kbotbg
! top interface of gwd region
    INTEGER :: ktopbg, ktoporo
! top interface of low level stress divergence region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! temperature tendency
    REAL :: ttgw(pver)
! zonal wind tendency
    REAL :: utgw(pver)
! meridional wind tendency
    REAL :: vtgw(pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! interface density
    REAL :: rhoi(0:pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
! interface temperature
    REAL :: ti(0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
    REAL :: utosrc(pver)
    REAL :: vtosrc(pver)
    REAL :: ttosrc(pver)
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
!-----------------------------------------------------------------------------
! Assign newtonian cooling coefficients
! -------------------------------------
    DO k=0,pver
      alpha(k) = 0.0
      dback(k) = 0.0
    END DO
! Assign wave phase speeds
! ------------------------
    DO kc=-pgwv,pgwv
      c(kc) = 10.0*kc
    END DO
! Determine the bounds of the background and orographic stress regions
    ktoporo = 0
    ktopbg = 0
    kbotoro = pver
    DO k=0,pver
      IF (pref_dev(k+1) .LT. p_src) kbotbg = k
! spectrum source at 400 mb
    END DO
    DO k=0,pver
! Profiles of background state variables
      CALL gw_prof(i, k, pcols, pver, u_dev, v_dev, t_dev, pmid_dev, &
&             pint_dev, rhoi, ni, ti, nm)
    END DO
!-----------------------------------------------------------------------------
! Non-orographic backgound gravity wave spectrum
!-----------------------------------------------------------------------------
    IF (pgwv .GT. 0 .AND. do_bgnd) THEN
! Determine the wave source for a background spectrum at ~400 mb
      CALL gw_bgnd(i, pcols, pver, c, u_dev, v_dev, t_dev, pmid_dev, &
&             pint_dev, pdel_dev, rpdel_dev, lnpint_dev, rlat_dev, kldv, &
&             kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, pgwv, &
&             kbotbg)
! Solve for the drag profile
      CALL gw_drag_prof_bgnd(i, pcols, pver, pgwv, pgwv, kbotbg, ktopbg&
&                       , c, u_dev, v_dev, t_dev, pint_dev, pdel_dev, &
&                       rpdel_dev, lnpint_dev, rlat_dev, rhoi, ni, ti, nm&
&                       , dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, &
&                       rdpldv, tau, ubi, ubm, xv, yv, utgw, vtgw, ttgw, &
&                       taubx_dev, tauby_dev, feb_dev, fepb_dev, utosrc, &
&                       vtosrc, ttosrc, tau0x, tau0y, effgwbkg_dev)
! Add the momentum tendencies to the output tendency arrays
      DO k=1,pver
        dudt_gwd_dev(i, k) = utgw(k) + utosrc(k)
        dvdt_gwd_dev(i, k) = vtgw(k) + vtosrc(k)
      END DO
    ELSE
!             dudt_gwd_dev(i,k) = utgw(k) + utosrc(k)
!             dvdt_gwd_dev(i,k) = vtgw(k) + vtosrc(k)
!             dtdt_gwd_dev(i,k) = ttgw(k) + ttosrc(k)
! zero net tendencies if no spectrum computed
      DO k=1,pver
        dudt_gwd_dev(i, k) = 0.
        dvdt_gwd_dev(i, k) = 0.
      END DO
    END IF
!            dtdt_gwd_dev(i,k) = 0.
!-----------------------------------------------------------------------------
! Orographic stationary gravity wave
!-----------------------------------------------------------------------------
! Determine the orographic wave source
    IF (do_oro) THEN
      CALL gw_oro(i, pcols, pver, pgwv, u_dev, v_dev, t_dev, sgh_dev, &
&            pmid_dev, pint_dev, pdel_dev, zm_dev, nm, kldv, kldvmn, ksrc&
&            , ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, kbotoro, rlat_dev)
! Solve for the drag profile
      CALL gw_drag_prof(i, pcols, pver, pgwv, 0, kbotoro, ktoporo, c, &
&                  u_dev, v_dev, t_dev, pint_dev, pdel_dev, rpdel_dev, &
&                  lnpint_dev, rlat_dev, rhoi, ni, ti, nm, dt, alpha, &
&                  dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, &
&                  ubm, xv, yv, utgw, vtgw, ttgw, tauox_dev, tauoy_dev, &
&                  feo_dev, fepo_dev, utosrc, vtosrc, ttosrc, tau0x, &
&                  tau0y, effgworo_dev)
! Add the orographic tendencies to the spectrum tendencies
! Compute the temperature tendency from energy conservation (includes spectrum).
      DO k=1,pver
        dudt_gwd_dev(i, k) = dudt_gwd_dev(i, k) + utgw(k)
        dvdt_gwd_dev(i, k) = dvdt_gwd_dev(i, k) + vtgw(k)
      END DO
    END IF
!          dtdt_gwd_dev(i,k) = dtdt_gwd_dev(i,k) + ttgw(k)
    DO k=1,pver
      u_dev(i, k) = u_dev(i, k) + dudt_gwd_dev(i, k)*dt
      v_dev(i, k) = v_dev(i, k) + dvdt_gwd_dev(i, k)*dt
    END DO
!       t_dev(i,k) = t_dev(i,k) + dtdt_gwd_dev(i,k)*dt
    RETURN
  END SUBROUTINE gw_intr
!================================================================================
  SUBROUTINE gw_prof(i, k, pcols, pver, u, v, t, pm, pi, rhoi, ni, ti, &
&    nm)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! current atmospheric column
    INTEGER, INTENT(IN) :: i
! current atmospheric layer
    INTEGER, INTENT(IN) :: k
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of vertical layers
    INTEGER, INTENT(IN) :: pver
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! midpoint pressures
    REAL :: pm(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! interface density
    REAL :: rhoi(0:pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! interface temperature
    REAL :: ti(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
!---------------------------Local storage-------------------------------
    REAL :: dtdp
! Brunt-Vaisalla frequency squared
    REAL :: n2
! Brunt-Vaisalla frequency squared
    REAL :: t_new, t_new_
    INTRINSIC max
    INTRINSIC sqrt
    REAL :: max1
!-----------------------------------------------------------------------------
! Determine the interface densities and Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------
! The top interface values are calculated assuming an isothermal atmosphere 
! above the top level.
    IF (k .EQ. 0) THEN
      ti(k) = t(i, k+1)
      CALL get_ti(t_new, ti(k))
      rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
      ni(k) = sqrt(mapl_grav*mapl_grav/(mapl_cp*t_new))
! Interior points use centered differences
    ELSE IF (k .GT. 0 .AND. k .LT. pver) THEN
      ti(k) = 0.5*(t(i, k)+t(i, k+1))
      CALL get_ti(t_new, ti(k))
      rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
      dtdp = (t(i, k+1)-t(i, k))/(pm(i, k+1)-pm(i, k))
      CALL get_ti(t_new_, dtdp)
      n2 = mapl_grav*mapl_grav/t_new*(1./mapl_cp-rhoi(k)*t_new_)
      IF (n2min .LT. n2) THEN
        max1 = n2
      ELSE
        max1 = n2min
      END IF
      ni(k) = sqrt(max1)
! Bottom interface uses bottom level temperature, density; next interface
! B-V frequency.
    ELSE IF (k .EQ. pver) THEN
      ti(k) = t(i, k)
      CALL get_ti(t_new, ti(k))
      rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
      ni(k) = ni(k-1)
    END IF
!-----------------------------------------------------------------------------
! Determine the midpoint Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------
    IF (k .GT. 0) nm(k) = 0.5*(ni(k-1)+ni(k))
    RETURN
  END SUBROUTINE gw_prof
!  Differentiation of gw_oro in reverse (adjoint) mode (forward sweep):
!   gradient, with respect to input variables: tau u v ubi ubm
!   of linear combination of output variables: tau u v xv yv ubi
!================================================================================
  SUBROUTINE gw_oro_fwd(i, pcols, pver, pgwv, u, v, t, sgh, pm, pi, dpm&
&    , zm, nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv&
&    , kbot, rlat)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! number of waves allowed
    INTEGER, INTENT(IN) :: pgwv
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! standard deviation of orography
    REAL :: sgh(pcols)
! midpoint pressures
    REAL :: pm(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! midpoint heights
    REAL :: zm(pcols, pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! top interface of low level stress div region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! 1/dp across low level divergence region
    REAL :: rdpldv
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
    INTEGER :: kbot
    REAL :: rlat(pcols)
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k
! Source-layer basic-state wind
    REAL :: ubsrc
! surface streamline displacment height (2*sgh)
    REAL :: hdsp
! max orographic sdv to use
    REAL :: sghmax
! c=0 stress from orography
    REAL :: tauoro
!    real    :: zldv                              ! top of the low level stress divergence region
! b-f frequency averaged over source region
    REAL :: nsrc
! interface pressure at top of source region
    REAL :: psrc
! density averaged over source region
    REAL :: rsrc
! u wind averaged over source region
    REAL :: usrc
! v wind averaged over source region
    REAL :: vsrc
    REAL :: u_new, v_new, t_new
    REAL :: min1
    INTRINSIC min
    INTRINSIC sqrt
! Begins
!---------------------------------------------------------------------------
! 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.
!---------------------------------------------------------------------------
    ksrc = pver - 1
    psrc = pi(i, pver-1)
    CALL get_ti(t_new, t(i, pver))
    rsrc = pm(i, pver)/(mapl_rgas*t_new)*dpm(i, pver)
    CALL get_uv(u_new, u(i, pver))
    CALL get_uv(v_new, v(i, pver))
    usrc = u_new*dpm(i, pver)
    vsrc = v_new*dpm(i, pver)
    nsrc = nm(pver)*dpm(i, pver)
    hdsp = 2.0*sgh(i)
    DO k=pver-1,pver/2,-1
      IF (hdsp .GT. sqrt(zm(i, k)*zm(i, k+1))) THEN
        ksrc = k - 1
        psrc = pi(i, k-1)
        CALL get_ti(t_new, t(i, k))
        rsrc = rsrc + pm(i, k)/(mapl_rgas*t_new)*dpm(i, k)
        CALL get_uv(u_new, u(i, k))
        CALL get_uv(v_new, v(i, k))
        usrc = usrc + u_new*dpm(i, k)
        vsrc = vsrc + v_new*dpm(i, k)
        nsrc = nsrc + nm(k)*dpm(i, k)
        CALL pushinteger4(2)
      ELSE
        CALL pushinteger4(1)
      END IF
    END DO
    rsrc = rsrc/(pi(i, pver)-psrc)
    usrc = usrc/(pi(i, pver)-psrc)
    vsrc = vsrc/(pi(i, pver)-psrc)
    nsrc = nsrc/(pi(i, pver)-psrc)
    IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
      ubsrc = sqrt(ubmc2mn)
      CALL pushreal4(xv)
      xv = 1.
      CALL pushreal4(yv)
      yv = 0.
      CALL pushinteger4(0)
    ELSE
      ubsrc = sqrt(usrc**2 + vsrc**2)
      CALL pushreal4(xv)
      xv = usrc/ubsrc
      CALL pushreal4(yv)
      yv = vsrc/ubsrc
      CALL pushinteger4(1)
    END IF
! Project the local wind at midpoints onto the source wind.
    DO k=1,pver
      ubm(k) = u(i, k)*xv + v(i, k)*yv
    END DO
    CALL pushreal4(ubi(0))
! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
    ubi(0) = ubm(1)
    DO k=1,pver
      CALL pushreal4(ubi(k))
      ubi(k) = ubm(k)
    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.
    IF (ubsrc .GT. orovmin .AND. hdsp .GT. orohmin) THEN
      sghmax = fcrit2*(ubsrc/nsrc)**2
      IF (hdsp**2 .GT. sghmax) THEN
        min1 = sghmax
        CALL pushinteger4(1)
      ELSE
        min1 = hdsp**2
        CALL pushinteger4(0)
      END IF
      tauoro = oroko2*min1*rsrc*nsrc*ubsrc
      CALL pushinteger4(0)
    ELSE
      tauoro = 0.
      ksrc = pver
      CALL pushinteger4(1)
    END IF
! tauoro is nontrivial when ubsrc is positive. However, if ubi(ksrc) is negative
! [note that the sign of ubsrc is irrelevant to the sign of ubi(ksrc)], orographic
! GWs can propagative upward passing through the negative basic-state wind.
! The following is to prevent this physically unjustified simulation.
! In addition, even if ubi(ksrc) > 0, if ubm(ksrc) < 0 .and. ubi(ksrc-1) < 0,
! negative wave stress leads to the acceleration of the negative ubm(ksrc).
! This result is physically inconsistent. However, if ubm(ksrc) > 0 .and.
! ubi(ksrc-1) < 0, GWs are filtered in physically consistent way, and
! decelerate the positive ubm(ksrc). Therefore, GWs are also assumed not to be
! launched when ubm(ksrc) < 0.
    IF (ubi(kbot) .LT. 0. .OR. ubm(kbot) .LT. 0.) THEN
      tauoro = 0.
      ksrc = pver
      CALL pushinteger4(1)
    ELSE
      CALL pushinteger4(0)
    END IF
! Sets kbot equal to ksrc
    kbot = ksrc
! 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).
    tau(0, kbot) = tauoro
! 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.
    CALL pushreal4(rsrc)
    CALL pushreal4(ubsrc)
    CALL pushreal4(psrc)
    CALL pushreal4(nsrc)
    CALL pushreal4(min1)
    CALL pushreal4(vsrc)
    CALL pushreal4(usrc)
  END SUBROUTINE gw_oro_fwd
!  Differentiation of gw_oro in reverse (adjoint) mode (backward sweep):
!   gradient, with respect to input variables: tau u v ubi ubm
!   of linear combination of output variables: tau u v xv yv ubi
!================================================================================
  SUBROUTINE gw_oro_bwd(i, pcols, pver, pgwv, u, ub, v, vb, t, sgh, pm, &
&    pi, dpm, zm, nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taub, ubi&
&    , ubib, ubm, ubmb, xv, xvb, yv, yvb, kbot, rlat)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i
    INTEGER, INTENT(IN) :: pcols
    INTEGER, INTENT(IN) :: pver
    INTEGER, INTENT(IN) :: pgwv
    REAL :: u(pcols, pver)
    REAL :: ub(pcols, pver)
    REAL :: v(pcols, pver)
    REAL :: vb(pcols, pver)
    REAL :: t(pcols, pver)
    REAL :: sgh(pcols)
    REAL :: pm(pcols, pver)
    REAL :: pi(pcols, 0:pver)
    REAL :: dpm(pcols, pver)
    REAL :: zm(pcols, pver)
    REAL :: nm(pver)
    INTEGER :: kldv
    INTEGER :: kldvmn
    INTEGER :: ksrc
    INTEGER :: ksrcmn
    REAL :: rdpldv
    REAL :: tau(-pgwv:pgwv, 0:pver)
    REAL :: taub(-pgwv:pgwv, 0:pver)
    REAL :: ubi(0:pver)
    REAL :: ubib(0:pver)
    REAL :: ubm(pver)
    REAL :: ubmb(pver)
    REAL :: xv
    REAL :: xvb
    REAL :: yv
    REAL :: yvb
    INTEGER :: kbot
    REAL :: rlat(pcols)
    INTEGER :: k
    REAL :: ubsrc
    REAL :: ubsrcb
    REAL :: hdsp
    REAL :: sghmax
    REAL :: sghmaxb
    REAL :: tauoro
    REAL :: tauorob
    REAL :: nsrc
    REAL :: psrc
    REAL :: rsrc
    REAL :: usrc
    REAL :: usrcb
    REAL :: vsrc
    REAL :: vsrcb
    REAL :: u_new, v_new, t_new
    REAL :: u_newb, v_newb
    INTEGER :: branch
    REAL :: min1
    REAL :: tempb0
    REAL :: min1b
    REAL :: tempb
    INTRINSIC min
    INTRINSIC sqrt
    CALL popreal4(usrc)
    CALL popreal4(vsrc)
    CALL popreal4(min1)
    CALL popreal4(nsrc)
    CALL popreal4(psrc)
    CALL popreal4(ubsrc)
    CALL popreal4(rsrc)
    tauorob = taub(0, kbot)
    taub(0, kbot) = 0.0_8
    CALL popinteger4(branch)
    IF (.NOT.branch .LT. 1) tauorob = 0.0_8
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      tempb0 = oroko2*rsrc*nsrc*tauorob
      min1b = ubsrc*tempb0
      ubsrcb = min1*tempb0
      CALL popinteger4(branch)
      IF (branch .LT. 1) THEN
        sghmaxb = 0.0_8
      ELSE
        sghmaxb = min1b
      END IF
      ubsrcb = ubsrcb + fcrit2*2*ubsrc*sghmaxb/nsrc**2
    ELSE
      ubsrcb = 0.0_8
    END IF
    ubmb(:) = 0.0_8
    DO k=pver,1,-1
      CALL popreal4(ubi(k))
      ubmb(k) = ubmb(k) + ubib(k)
      ubib(k) = 0.0_8
    END DO
    CALL popreal4(ubi(0))
    ubmb(1) = ubmb(1) + ubib(0)
    ubib(0) = 0.0_8
    DO k=pver,1,-1
      ub(i, k) = ub(i, k) + xv*ubmb(k)
      xvb = xvb + u(i, k)*ubmb(k)
      vb(i, k) = vb(i, k) + yv*ubmb(k)
      yvb = yvb + v(i, k)*ubmb(k)
      ubmb(k) = 0.0_8
    END DO
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      CALL popreal4(yv)
      CALL popreal4(xv)
      vsrcb = 0.0_8
      usrcb = 0.0_8
    ELSE
      CALL popreal4(yv)
      ubsrcb = ubsrcb - usrc*xvb/ubsrc**2 - vsrc*yvb/ubsrc**2
      IF (usrc**2 + vsrc**2 .EQ. 0.0) THEN
        tempb = 0.0
      ELSE
        tempb = ubsrcb/(2.0*sqrt(usrc**2+vsrc**2))
      END IF
      vsrcb = 2*vsrc*tempb + yvb/ubsrc
      CALL popreal4(xv)
      usrcb = 2*usrc*tempb + xvb/ubsrc
    END IF
    vsrcb = vsrcb/(pi(i, pver)-psrc)
    usrcb = usrcb/(pi(i, pver)-psrc)
    DO k=pver/2,pver-1,1
      CALL popinteger4(branch)
      IF (.NOT.branch .LT. 2) THEN
        v_newb = dpm(i, k)*vsrcb
        u_newb = dpm(i, k)*usrcb
        CALL get_uv_b(v_new, v_newb, v(i, k), vb(i, k))
        CALL get_uv_b(u_new, u_newb, u(i, k), ub(i, k))
      END IF
    END DO
    v_newb = dpm(i, pver)*vsrcb
    u_newb = dpm(i, pver)*usrcb
    CALL get_uv_b(v_new, v_newb, v(i, pver), vb(i, pver))
    CALL get_uv_b(u_new, u_newb, u(i, pver), ub(i, pver))
  END SUBROUTINE gw_oro_bwd
!================================================================================
  SUBROUTINE gw_oro(i, pcols, pver, pgwv, u, v, t, sgh, pm, pi, dpm, zm&
&    , nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
&    kbot, rlat)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! number of waves allowed
    INTEGER, INTENT(IN) :: pgwv
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! standard deviation of orography
    REAL :: sgh(pcols)
! midpoint pressures
    REAL :: pm(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! midpoint heights
    REAL :: zm(pcols, pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! top interface of low level stress div region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! 1/dp across low level divergence region
    REAL :: rdpldv
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
    INTEGER :: kbot
    REAL :: rlat(pcols)
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k
! Source-layer basic-state wind
    REAL :: ubsrc
! surface streamline displacment height (2*sgh)
    REAL :: hdsp
! max orographic sdv to use
    REAL :: sghmax
! c=0 stress from orography
    REAL :: tauoro
!    real    :: zldv                              ! top of the low level stress divergence region
! b-f frequency averaged over source region
    REAL :: nsrc
! interface pressure at top of source region
    REAL :: psrc
! density averaged over source region
    REAL :: rsrc
! u wind averaged over source region
    REAL :: usrc
! v wind averaged over source region
    REAL :: vsrc
    REAL :: u_new, v_new, t_new
    REAL :: min1
    INTRINSIC min
    INTRINSIC sqrt
! Begins
!---------------------------------------------------------------------------
! 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.
!---------------------------------------------------------------------------
    ksrc = pver - 1
    kldv = pver - 1
    psrc = pi(i, pver-1)
    CALL get_ti(t_new, t(i, pver))
    rsrc = pm(i, pver)/(mapl_rgas*t_new)*dpm(i, pver)
    CALL get_uv(u_new, u(i, pver))
    CALL get_uv(v_new, v(i, pver))
    usrc = u_new*dpm(i, pver)
    vsrc = v_new*dpm(i, pver)
    nsrc = nm(pver)*dpm(i, pver)
    hdsp = 2.0*sgh(i)
    DO k=pver-1,pver/2,-1
      IF (hdsp .GT. sqrt(zm(i, k)*zm(i, k+1))) THEN
        ksrc = k - 1
        kldv = k - 1
        psrc = pi(i, k-1)
        CALL get_ti(t_new, t(i, k))
        rsrc = rsrc + pm(i, k)/(mapl_rgas*t_new)*dpm(i, k)
        CALL get_uv(u_new, u(i, k))
        CALL get_uv(v_new, v(i, k))
        usrc = usrc + u_new*dpm(i, k)
        vsrc = vsrc + v_new*dpm(i, k)
        nsrc = nsrc + nm(k)*dpm(i, k)
      END IF
    END DO
    rsrc = rsrc/(pi(i, pver)-psrc)
    usrc = usrc/(pi(i, pver)-psrc)
    vsrc = vsrc/(pi(i, pver)-psrc)
    nsrc = nsrc/(pi(i, pver)-psrc)
    IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
      ubsrc = sqrt(ubmc2mn)
      xv = 1.
      yv = 0.
    ELSE
      ubsrc = sqrt(usrc**2 + vsrc**2)
      xv = usrc/ubsrc
      yv = vsrc/ubsrc
    END IF
! Project the local wind at midpoints onto the source wind.
    DO k=1,pver
      ubm(k) = u(i, k)*xv + v(i, k)*yv
    END DO
! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
    ubi(0) = ubm(1)
    DO k=1,pver
      ubi(k) = ubm(k)
    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.
    IF (ubsrc .GT. orovmin .AND. hdsp .GT. orohmin) THEN
      sghmax = fcrit2*(ubsrc/nsrc)**2
      IF (hdsp**2 .GT. sghmax) THEN
        min1 = sghmax
      ELSE
        min1 = hdsp**2
      END IF
      tauoro = oroko2*min1*rsrc*nsrc*ubsrc
    ELSE
      tauoro = 0.
      ksrc = pver
      kldv = pver
    END IF
! tauoro is nontrivial when ubsrc is positive. However, if ubi(ksrc) is negative
! [note that the sign of ubsrc is irrelevant to the sign of ubi(ksrc)], orographic
! GWs can propagative upward passing through the negative basic-state wind.
! The following is to prevent this physically unjustified simulation.
! In addition, even if ubi(ksrc) > 0, if ubm(ksrc) < 0 .and. ubi(ksrc-1) < 0,
! negative wave stress leads to the acceleration of the negative ubm(ksrc).
! This result is physically inconsistent. However, if ubm(ksrc) > 0 .and.
! ubi(ksrc-1) < 0, GWs are filtered in physically consistent way, and
! decelerate the positive ubm(ksrc). Therefore, GWs are also assumed not to be
! launched when ubm(ksrc) < 0.
    IF (ubi(kbot) .LT. 0. .OR. ubm(kbot) .LT. 0.) THEN
      tauoro = 0.
      ksrc = pver
      kldv = pver
    END IF
! Sets kbot equal to ksrc
    kbot = ksrc
! 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).
    tau(0, kbot) = tauoro
! 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
    IF (ksrcmn .GT. ksrc) THEN
      ksrcmn = ksrc
    ELSE
      ksrcmn = ksrcmn
    END IF
    IF (kldvmn .GT. kldv) THEN
      kldvmn = kldv
    ELSE
      kldvmn = kldvmn
    END IF
    IF (kldv .NE. pver) rdpldv = 1./(pi(i, kldv)-pi(i, pver))
! kldvmn is always pver because FRACLDV == 0.
    IF (fracldv .LE. 0.) kldvmn = pver
    RETURN
  END SUBROUTINE gw_oro
!  Differentiation of gw_bgnd in reverse (adjoint) mode (forward sweep):
!   gradient, with respect to input variables: u v
!   of linear combination of output variables: u v xv yv ubi ubm
!===============================================================================
  SUBROUTINE gw_bgnd_fwd(i, pcols, pver, c, u, v, t, pm, pi, dpm, rdpm, &
&    piln, rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, &
&    yv, ngwv, kbot)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
! 
! The parameterization is assumed to operate only where water vapor 
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! number of current column
    INTEGER :: i
! number of atmospheric columns
    INTEGER :: pcols
! number of atmospheric columns
    INTEGER :: pver
! index of bottom (source) interface
    INTEGER :: kbot
! number of gravity waves to use
    INTEGER :: ngwv
! wave phase speeds
    REAL :: c(-ngwv:ngwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! midpoint pressures
    REAL :: pm(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
! latitude in radians for columns
    REAL :: rlat(pcols)
! top interface of low level stress divergence region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! 1/dp across low level divergence region
    REAL :: rdpldv
! wave Reynolds stress
    REAL :: tau(-ngwv:ngwv, 0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! background stress at c=0
    REAL :: tauback
! u wind averaged over source region
    REAL :: usrc
! v wind averaged over source region
    REAL :: vsrc
    REAL :: ubsrc
! Used in lat dependence of GW spec. 
    REAL :: al0
! Used in lat dependence of GW spec.
    REAL :: dlat0
    REAL :: latdeg
! The actual lat dependence of GW spec.
    REAL :: flat_gw
    REAL :: u_new, v_new, tau_new
    REAL :: tmp
    INTRINSIC exp
    INTRINSIC max
    INTRINSIC abs
    REAL :: x1
    INTRINSIC dble
    REAL :: abs6
    REAL :: abs5
    REAL :: abs4
    REAL :: abs3
    REAL :: abs2
    REAL :: abs1
    INTRINSIC sqrt
    REAL :: y1
!---------------------------------------------------------------------------
! 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).
    ubm = 0.
    ubi = 0.
    tau = 0.
    usrc = 0.5*(u(i, kbot+1)+u(i, kbot))
    vsrc = 0.5*(v(i, kbot+1)+v(i, kbot))
    x1 = sqrt(usrc**2 + vsrc**2)
    y1 = sqrt(ubmc2mn)
    IF (x1 .LT. y1) THEN
      ubsrc = y1
      CALL pushinteger4(1)
    ELSE
      ubsrc = x1
      CALL pushinteger4(0)
    END IF
    IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
      xv = 1.0
      yv = 0.0
      CALL pushinteger4(1)
    ELSE
      xv = usrc/ubsrc
      yv = vsrc/ubsrc
      CALL pushinteger4(0)
    END IF
! Project the local wind at midpoints onto the source wind.
    DO k=1,pver
      ubm(k) = u(i, k)*xv + v(i, k)*yv
    END DO
! Compute the bottom interface wind projection using the midpoint winds.
    ubi(0) = ubm(1)
    DO k=1,pver
      ubi(k) = ubm(k)
    END DO
    CALL pushreal4(ubsrc)
    CALL pushreal4(vsrc)
    CALL pushreal4(usrc)
!-----------------------------------------------------------------------
! Gravity wave sources
!-----------------------------------------------------------------------
! Determine the background stress at c=0
    tauback = taubgnd*tauscal
! Include dependence on latitude:
    latdeg = rlat(i)*180./pi_gwd
!
    IF (-15.3 .LT. latdeg .AND. latdeg .LT. 15.3) THEN
      IF (latdeg .GE. 0.) THEN
        abs1 = latdeg
      ELSE
        abs1 = -latdeg
      END IF
      flat_gw = 1.2*exp(-(dble((abs1-3.)/8.0)**2))
      IF (latdeg .GE. 0.) THEN
        abs2 = latdeg
      ELSE
        abs2 = -latdeg
      END IF
      IF (flat_gw .LT. 1.2 .AND. abs2 .LE. 3.) flat_gw = 1.2
    ELSE IF (latdeg .GT. -31. .AND. latdeg .LE. -15.3) THEN
      flat_gw = 0.10
    ELSE IF (latdeg .LT. 31. .AND. latdeg .GE. 15.3) THEN
      flat_gw = 0.10
    ELSE IF (latdeg .GT. -60. .AND. latdeg .LE. -31.) THEN
      IF (latdeg .GE. 0.) THEN
        abs3 = latdeg
      ELSE
        abs3 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs3-60.)/23.)**2))
    ELSE IF (latdeg .LT. 60. .AND. latdeg .GE. 31.) THEN
      IF (latdeg .GE. 0.) THEN
        abs4 = latdeg
      ELSE
        abs4 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs4-60.)/23.)**2))
    ELSE IF (latdeg .LE. -60.) THEN
      IF (latdeg .GE. 0.) THEN
        abs5 = latdeg
      ELSE
        abs5 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs5-60.)/70.)**2))
    ELSE IF (latdeg .GE. 60.) THEN
      IF (latdeg .GE. 0.) THEN
        abs6 = latdeg
      ELSE
        abs6 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs6-60.)/70.)**2))
    END IF
    tauback = tauback*flat_gw
! 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
      tau(l, kbot) = tauback*exp(-((c(l)/30.)**2))
      tmp = tau(l, kbot)
      tau(-l, kbot) = tmp
    END DO
    tau(0, kbot) = tauback
! 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.
  END SUBROUTINE gw_bgnd_fwd
!  Differentiation of gw_bgnd in reverse (adjoint) mode (backward sweep):
!   gradient, with respect to input variables: u v
!   of linear combination of output variables: u v xv yv ubi ubm
!===============================================================================
  SUBROUTINE gw_bgnd_bwd(i, pcols, pver, c, u, ub, v, vb, t, pm, pi, dpm&
&    , rdpm, piln, rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, &
&    ubib, ubm, ubmb, xv, xvb, yv, yvb, ngwv, kbot)
    IMPLICIT NONE
    INTEGER :: i
    INTEGER :: pcols
    INTEGER :: pver
    INTEGER :: kbot
    INTEGER :: ngwv
    REAL :: c(-ngwv:ngwv)
    REAL :: u(pcols, pver)
    REAL :: ub(pcols, pver)
    REAL :: v(pcols, pver)
    REAL :: vb(pcols, pver)
    REAL :: t(pcols, pver)
    REAL :: pm(pcols, pver)
    REAL :: pi(pcols, 0:pver)
    REAL :: dpm(pcols, pver)
    REAL :: rdpm(pcols, pver)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
    INTEGER :: kldv
    INTEGER :: kldvmn
    INTEGER :: ksrc
    INTEGER :: ksrcmn
    REAL :: rdpldv
    REAL :: tau(-ngwv:ngwv, 0:pver)
    REAL :: ubi(0:pver)
    REAL :: ubib(0:pver)
    REAL :: ubm(pver)
    REAL :: ubmb(pver)
    REAL :: xv
    REAL :: xvb
    REAL :: yv
    REAL :: yvb
    INTEGER :: k, l
    REAL :: tauback
    REAL :: usrc
    REAL :: usrcb
    REAL :: vsrc
    REAL :: vsrcb
    REAL :: ubsrc
    REAL :: ubsrcb
    REAL :: al0
    REAL :: dlat0
    REAL :: latdeg
    REAL :: flat_gw
    REAL :: u_new, v_new, tau_new
    INTEGER :: branch
    INTRINSIC exp
    INTRINSIC max
    INTRINSIC abs
    REAL :: x1
    REAL :: tempb
    REAL :: x1b
    INTRINSIC dble
    REAL :: abs6
    REAL :: abs5
    REAL :: abs4
    REAL :: abs3
    REAL :: abs2
    REAL :: abs1
    INTRINSIC sqrt
    REAL :: y1
    CALL popreal4(usrc)
    CALL popreal4(vsrc)
    CALL popreal4(ubsrc)
    DO k=pver,1,-1
      ubmb(k) = ubmb(k) + ubib(k)
      ubib(k) = 0.0_8
    END DO
    ubmb(1) = ubmb(1) + ubib(0)
    DO k=pver,1,-1
      ub(i, k) = ub(i, k) + xv*ubmb(k)
      xvb = xvb + u(i, k)*ubmb(k)
      vb(i, k) = vb(i, k) + yv*ubmb(k)
      yvb = yvb + v(i, k)*ubmb(k)
      ubmb(k) = 0.0_8
    END DO
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      vsrcb = yvb/ubsrc
      ubsrcb = -(usrc*xvb/ubsrc**2) - vsrc*yvb/ubsrc**2
      usrcb = xvb/ubsrc
    ELSE
      ubsrcb = 0.0_8
      vsrcb = 0.0_8
      usrcb = 0.0_8
    END IF
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      x1b = ubsrcb
    ELSE
      x1b = 0.0_8
    END IF
    IF (usrc**2 + vsrc**2 .EQ. 0.0) THEN
      tempb = 0.0
    ELSE
      tempb = x1b/(2.0*sqrt(usrc**2+vsrc**2))
    END IF
    usrcb = usrcb + 2*usrc*tempb
    vsrcb = vsrcb + 2*vsrc*tempb
    vb(i, kbot+1) = vb(i, kbot+1) + 0.5*vsrcb
    vb(i, kbot) = vb(i, kbot) + 0.5*vsrcb
    ub(i, kbot+1) = ub(i, kbot+1) + 0.5*usrcb
    ub(i, kbot) = ub(i, kbot) + 0.5*usrcb
  END SUBROUTINE gw_bgnd_bwd
!===============================================================================
  SUBROUTINE gw_bgnd(i, pcols, pver, c, u, v, t, pm, pi, dpm, rdpm, piln&
&    , rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
&    ngwv, kbot)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
! 
! The parameterization is assumed to operate only where water vapor 
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
! number of current column
    INTEGER :: i
! number of atmospheric columns
    INTEGER :: pcols
! number of atmospheric columns
    INTEGER :: pver
! index of bottom (source) interface
    INTEGER :: kbot
! number of gravity waves to use
    INTEGER :: ngwv
! wave phase speeds
    REAL :: c(-ngwv:ngwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! midpoint pressures
    REAL :: pm(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
! latitude in radians for columns
    REAL :: rlat(pcols)
! top interface of low level stress divergence region
    INTEGER :: kldv
! min value of kldv
    INTEGER :: kldvmn
! index of top interface of source region
    INTEGER :: ksrc
! min value of ksrc
    INTEGER :: ksrcmn
! 1/dp across low level divergence region
    REAL :: rdpldv
! wave Reynolds stress
    REAL :: tau(-ngwv:ngwv, 0:pver)
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! background stress at c=0
    REAL :: tauback
! u wind averaged over source region
    REAL :: usrc
! v wind averaged over source region
    REAL :: vsrc
    REAL :: ubsrc
! Used in lat dependence of GW spec. 
    REAL :: al0
! Used in lat dependence of GW spec.
    REAL :: dlat0
    REAL :: latdeg
! The actual lat dependence of GW spec.
    REAL :: flat_gw
    REAL :: u_new, v_new, tau_new
    INTRINSIC exp
    INTRINSIC max
    INTRINSIC abs
    REAL :: x1
    INTRINSIC dble
    REAL :: abs6
    REAL :: abs5
    REAL :: abs4
    REAL :: abs3
    REAL :: abs2
    REAL :: abs1
    INTRINSIC sqrt
    REAL :: y1
!---------------------------------------------------------------------------
! 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).
    ubm = 0.
    ubi = 0.
    xv = 0.
    yv = 0.
    tau = 0.
    ksrc = kbot
    kldv = kbot
    usrc = 0.5*(u(i, kbot+1)+u(i, kbot))
    vsrc = 0.5*(v(i, kbot+1)+v(i, kbot))
    x1 = sqrt(usrc**2 + vsrc**2)
    y1 = sqrt(ubmc2mn)
    IF (x1 .LT. y1) THEN
      ubsrc = y1
    ELSE
      ubsrc = x1
    END IF
    IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
      xv = 1.0
      yv = 0.0
    ELSE
      xv = usrc/ubsrc
      yv = vsrc/ubsrc
    END IF
! Project the local wind at midpoints onto the source wind.
    DO k=1,pver
      ubm(k) = u(i, k)*xv + v(i, k)*yv
    END DO
! Compute the bottom interface wind projection using the midpoint winds.
    ubi(0) = ubm(1)
    DO k=1,pver
      ubi(k) = ubm(k)
    END DO
!-----------------------------------------------------------------------
! Gravity wave sources
!-----------------------------------------------------------------------
! Determine the background stress at c=0
    tauback = taubgnd*tauscal
! Include dependence on latitude:
    latdeg = rlat(i)*180./pi_gwd
!
    IF (-15.3 .LT. latdeg .AND. latdeg .LT. 15.3) THEN
      IF (latdeg .GE. 0.) THEN
        abs1 = latdeg
      ELSE
        abs1 = -latdeg
      END IF
      flat_gw = 1.2*exp(-(dble((abs1-3.)/8.0)**2))
      IF (latdeg .GE. 0.) THEN
        abs2 = latdeg
      ELSE
        abs2 = -latdeg
      END IF
      IF (flat_gw .LT. 1.2 .AND. abs2 .LE. 3.) flat_gw = 1.2
    ELSE IF (latdeg .GT. -31. .AND. latdeg .LE. -15.3) THEN
      flat_gw = 0.10
    ELSE IF (latdeg .LT. 31. .AND. latdeg .GE. 15.3) THEN
      flat_gw = 0.10
    ELSE IF (latdeg .GT. -60. .AND. latdeg .LE. -31.) THEN
      IF (latdeg .GE. 0.) THEN
        abs3 = latdeg
      ELSE
        abs3 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs3-60.)/23.)**2))
    ELSE IF (latdeg .LT. 60. .AND. latdeg .GE. 31.) THEN
      IF (latdeg .GE. 0.) THEN
        abs4 = latdeg
      ELSE
        abs4 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs4-60.)/23.)**2))
    ELSE IF (latdeg .LE. -60.) THEN
      IF (latdeg .GE. 0.) THEN
        abs5 = latdeg
      ELSE
        abs5 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs5-60.)/70.)**2))
    ELSE IF (latdeg .GE. 60.) THEN
      IF (latdeg .GE. 0.) THEN
        abs6 = latdeg
      ELSE
        abs6 = -latdeg
      END IF
      flat_gw = 0.50*exp(-(dble((abs6-60.)/70.)**2))
    END IF
    tauback = tauback*flat_gw
! 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
      tau(l, kbot) = tauback*exp(-((c(l)/30.)**2))
      tau(-l, kbot) = tau(l, kbot)
    END DO
    tau(0, kbot) = tauback
! 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
!  Differentiation of gw_drag_prof in reverse (adjoint) mode (forward sweep):
!   gradient, with respect to input variables: tau xv ut yv vt
!                ubi
!   of linear combination of output variables: ut vt
!===============================================================================
  SUBROUTINE gw_drag_prof_fwd(i, pcols, pver, pgwv, ngwv, kbot, ktop, c&
&    , u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, &
&    dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut&
&    , vt, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, &
&    tau0y, effgw)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! index of bottom (source) interface
    INTEGER, INTENT(IN) :: kbot
! index of top interface of gwd region
    INTEGER, INTENT(IN) :: ktop
! number of gravity waves possible
    INTEGER, INTENT(IN) :: pgwv
! number of gravity waves to use
    INTEGER, INTENT(IN) :: ngwv
! top interface of low level stress  divergence region
    INTEGER, INTENT(IN) :: kldv
! min value of kldv
    INTEGER, INTENT(IN) :: kldvmn
! index of top interface of source region
    INTEGER, INTENT(IN) :: ksrc
! min value of ksrc
    INTEGER, INTENT(IN) :: ksrcmn
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
! interface density
    REAL :: rhoi(0:pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! interface temperature
    REAL :: ti(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! time step
    REAL :: dt
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
! tendency efficiency for gwd 
    REAL :: effgw(pcols)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! zonal wind tendency
    REAL :: ut(pver)
! meridional wind tendency
    REAL :: vt(pver)
! temperature tendency
    REAL :: tt(pver)
! Total zonal GW momentum flux
    REAL :: taugwx(pcols, 0:pver)
! Total meridional GW momentum flux
    REAL :: taugwy(pcols, 0:pver)
! Total GW energy flux
    REAL :: fegw(pcols, 0:pver)
! Total GW pseudo energy flux
    REAL :: fepgw(pcols, 0:pver)
! Total U tendency below launch level
    REAL :: dusrc(pver)
! Total V tendency below launch level
    REAL :: dvsrc(pver)
! Total V tendency below launch level
    REAL :: dtsrc(pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! fraction of dsat to use
    REAL :: dscal
! imaginary part of vertical wavenumber
    REAL :: mi
! stress after damping
    REAL :: taudmp
!real    :: taumax(pcols)                     ! max(tau) for any l
! saturation stress
    REAL :: tausat
! (ub-c)
    REAL :: ubmc
! (ub-c)**2
    REAL :: ubmc2
! ubar tendency
    REAL :: ubt
! tbar tendency
    REAL :: tbt
! ubar tendency from wave l
    REAL :: ubtl
! saturation tendency
    REAL :: ubtlsat
! layer pressure
    REAL :: pm
! layer density
    REAL :: rhom
! launch level height
    REAL :: zlb
! c-u
    REAL :: cmu
    REAL :: dzm, hscal, tautmp
    REAL :: utl
    REAL :: ttl
! zonal pseudomomentum flux spectrum
    REAL :: fpmx
! meridional pseudomomentum flux spectrum
    REAL :: fpmy
! energy flux (p'w') spectrum
    REAL :: fe
! pseudoenergy flux (p'w'+U rho u'w') spectrum
    REAL :: fpe
    REAL :: fpml
    REAL :: fpmt
    REAL :: fpel
    REAL :: fpet
    REAL :: dusrcl
    REAL :: dvsrcl
    REAL :: dtsrcl
    REAL :: zi
    REAL :: effkwvmap
    REAL :: zfac
    REAL :: uhtmax
    REAL :: utfac
    REAL :: tau_min, cmu_, t_new
    INTRINSIC max
    INTRINSIC sign
    INTRINSIC abs
    REAL :: x2
    REAL :: x1
    INTRINSIC sqrt
! Initialize gravity wave drag tendencies to zero
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) THEN
          tau(l, k) = 0.
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
    END DO
    DO k=1,pver
      ut(k) = 0.
      vt(k) = 0.
    END DO
! Initialize total momentum and energy fluxes to zero
    DO k=0,pver
      taugwx(i, k) = 0.
      taugwy(i, k) = 0.
      fegw(i, k) = 0.
      fepgw(i, k) = 0.
    END DO
!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------
! 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.
! Loop from bottom to top to get stress profiles
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) THEN
          CALL pushreal4(ubmc)
!!!!!!!!jk             d = dback(k)
          ubmc = ubi(k) - c(l)
          IF (ngwv .GT. 0) THEN
            CALL pushreal4(effkwvmap)
            CALL get_effkwvmap_1(effkwvmap, rlat(i))
            IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
&                .LT. 15.0) THEN
              CALL pushinteger4(0)
              effkwvmap = fcrit2*kwvbeq
            ELSE
              CALL pushinteger4(0)
            END IF
          ELSE
            IF (pi(i, k) .LT. 1000.0) THEN
              zfac = (pi(i, k)/1000.0)**3
            ELSE
              zfac = 1.0
            END IF
            IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
              CALL pushreal4(effkwvmap)
              CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
              CALL pushinteger4(1)
            ELSE
              CALL pushreal4(effkwvmap)
              CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
              CALL pushinteger4(2)
            END IF
          END IF
          x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
          IF (x1 .GE. 0.) THEN
            tausat = x1
            CALL pushinteger4(1)
          ELSE
            tausat = -x1
            CALL pushinteger4(0)
          END IF
          IF (tausat .LE. taumin) THEN
            tausat = 0.0
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) THEN
            tausat = 0.0
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop) THEN
            tausat = 0.
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
!
          IF (k .EQ. ktop + 1) THEN
            tausat = tausat*0.02
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 2) THEN
            tausat = tausat*0.05
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 3) THEN
            tausat = tausat*0.10
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 4) THEN
            tausat = tausat*0.20
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 5) THEN
            tausat = tausat*0.50
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
!
          tau_min = tausat
          IF (tau_min .GT. tau(l, k+1)) THEN
            tau_min = tau(l, k+1)
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          tau(l, k) = tau_min
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
    END DO
!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------
! Accumulate the mean wind tendency over wavenumber.
! Loop over levels from top to bottom
    DO k=ktop+1,pver
      CALL pushreal4(ubt)
      ubt = 0.0
      DO l=-ngwv,ngwv
        IF (k .LE. kbot) THEN
          CALL pushreal4(ubtl)
! Determine the wind tendency including excess stress carried down from above.
          ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
! Calculate the sign of wind tendency
          utl = sign(ubtl, c(l) - ubi(k))
! Accumulate the mean wind tendency over wavenumber.
          ubt = ubt + utl
! Calculate irreversible temperature tendency associated with gravity wave breaking.
! Adding frictional heating associated with the GW momentum forcing
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
! Project the mean wind tendency onto the components and scale by "efficiency".
      IF (k .LE. kbot) THEN
        ut(k) = ubt*xv*effgw(i)
        vt(k) = ubt*yv*effgw(i)
        CALL pushinteger4(2)
      ELSE
        CALL pushinteger4(1)
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates energy and momentum flux profiles
!-----------------------------------------------------------------------
    DO l=-ngwv,ngwv
      DO k=ktop,pver
        IF (k .LE. kbot) THEN
          cmu = c(l) - ubi(k)
          CALL get_cmu(cmu_, cmu)
          fpmx = cmu_*tau(l, k)*xv*effgw(i)
          fpmy = cmu_*tau(l, k)*yv*effgw(i)
          fe = cmu*cmu_*tau(l, k)*effgw(i)
          fpe = c(l)*cmu_*tau(l, k)*effgw(i)
! Record outputs for GW fluxes
          taugwx(i, k) = taugwx(i, k) + fpmx
          taugwy(i, k) = taugwy(i, k) + fpmy
          fegw(i, k) = fegw(i, k) + fe
          fepgw(i, k) = fepgw(i, k) + fpe
        END IF
      END DO
    END DO
!-----------------------------------------------------------------------
! Adjust efficiency factor to prevent unrealistically strong forcing
!-----------------------------------------------------------------------
    uhtmax = 0.0
    utfac = 1.0
    DO k=1,pver
      x2 = sqrt(ut(k)**2 + vt(k)**2)
      IF (x2 .LT. uhtmax) THEN
        uhtmax = uhtmax
        CALL pushinteger4(2)
      ELSE
        uhtmax = x2
        CALL pushinteger4(1)
      END IF
    END DO
    IF (uhtmax .GT. tndmax) THEN
      utfac = tndmax/uhtmax
      CALL pushinteger4(1)
    ELSE
      CALL pushinteger4(0)
    END IF
!jkim    effgw(i) = effgw(i)*utfac
    DO k=1,pver
      CALL pushreal4(ut(k))
      ut(k) = ut(k)*utfac
      CALL pushreal4(vt(k))
      vt(k) = vt(k)*utfac
    END DO
    CALL pushreal4(ubmc)
    CALL pushreal4(effkwvmap)
    CALL pushreal4(ubtl)
    CALL pushreal4(utfac)
    CALL pushreal4(uhtmax)
    CALL pushreal4(ubt)
  END SUBROUTINE gw_drag_prof_fwd
!  Differentiation of gw_drag_prof in reverse (adjoint) mode (backward sweep):
!   gradient, with respect to input variables: tau xv ut yv vt
!                ubi
!   of linear combination of output variables: ut vt
!===============================================================================
  SUBROUTINE gw_drag_prof_bwd(i, pcols, pver, pgwv, ngwv, kbot, ktop, c&
&    , u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, &
&    dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taub, ubi, ubib, ubm&
&    , xv, xvb, yv, yvb, ut, utb, vt, vtb, tt, taugwx, taugwy, fegw, &
&    fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y, effgw)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i
    INTEGER, INTENT(IN) :: pcols
    INTEGER, INTENT(IN) :: pver
    INTEGER, INTENT(IN) :: kbot
    INTEGER, INTENT(IN) :: ktop
    INTEGER, INTENT(IN) :: pgwv
    INTEGER, INTENT(IN) :: ngwv
    INTEGER, INTENT(IN) :: kldv
    INTEGER, INTENT(IN) :: kldvmn
    INTEGER, INTENT(IN) :: ksrc
    INTEGER, INTENT(IN) :: ksrcmn
    REAL :: c(-pgwv:pgwv)
    REAL :: u(pcols, pver)
    REAL :: v(pcols, pver)
    REAL :: t(pcols, pver)
    REAL :: pi(pcols, 0:pver)
    REAL :: dpm(pcols, pver)
    REAL :: rdpm(pcols, pver)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
    REAL :: rhoi(0:pver)
    REAL :: ni(0:pver)
    REAL :: ti(0:pver)
    REAL :: nm(pver)
    REAL :: dt
    REAL :: alpha(0:pver)
    REAL :: dback(0:pver)
    REAL :: rdpldv
    REAL :: ubi(0:pver)
    REAL :: ubib(0:pver)
    REAL :: ubm(pver)
    REAL :: xv
    REAL :: xvb
    REAL :: yv
    REAL :: yvb
    REAL :: effgw(pcols)
    REAL :: tau(-pgwv:pgwv, 0:pver)
    REAL :: taub(-pgwv:pgwv, 0:pver)
    REAL :: ut(pver)
    REAL :: utb(pver)
    REAL :: vt(pver)
    REAL :: vtb(pver)
    REAL :: tt(pver)
    REAL :: taugwx(pcols, 0:pver)
    REAL :: taugwy(pcols, 0:pver)
    REAL :: fegw(pcols, 0:pver)
    REAL :: fepgw(pcols, 0:pver)
    REAL :: dusrc(pver)
    REAL :: dvsrc(pver)
    REAL :: dtsrc(pver)
    REAL :: tau0x
    REAL :: tau0y
    INTEGER :: k, l
    REAL :: dscal
    REAL :: mi
    REAL :: taudmp
    REAL :: tausat
    REAL :: tausatb
    REAL :: ubmc
    REAL :: ubmcb
    REAL :: ubmc2
    REAL :: ubt
    REAL :: ubtb
    REAL :: tbt
    REAL :: ubtl
    REAL :: ubtlb
    REAL :: ubtlsat
    REAL :: pm
    REAL :: rhom
    REAL :: zlb
    REAL :: cmu
    REAL :: dzm, hscal, tautmp
    REAL :: utl
    REAL :: utlb
    REAL :: ttl
    REAL :: fpmx
    REAL :: fpmy
    REAL :: fe
    REAL :: fpe
    REAL :: fpml
    REAL :: fpmt
    REAL :: fpel
    REAL :: fpet
    REAL :: dusrcl
    REAL :: dvsrcl
    REAL :: dtsrcl
    REAL :: zi
    REAL :: effkwvmap
    REAL :: zfac
    REAL :: uhtmax
    REAL :: uhtmaxb
    REAL :: utfac
    REAL :: utfacb
    REAL :: tau_min, cmu_, t_new
    REAL :: tau_minb
    INTEGER :: branch
    REAL :: tempb0
    INTRINSIC max
    INTRINSIC sign
    INTRINSIC abs
    REAL :: x2
    REAL :: x1
    REAL :: x2b
    REAL :: tempb
    REAL :: x1b
    INTRINSIC sqrt
    CALL popreal4(ubt)
    CALL popreal4(uhtmax)
    CALL popreal4(utfac)
    CALL popreal4(ubtl)
    CALL popreal4(effkwvmap)
    CALL popreal4(ubmc)
    utfacb = 0.0_8
    DO k=pver,1,-1
      CALL popreal4(vt(k))
      CALL popreal4(ut(k))
      utfacb = utfacb + ut(k)*utb(k) + vt(k)*vtb(k)
      vtb(k) = utfac*vtb(k)
      utb(k) = utfac*utb(k)
    END DO
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      uhtmaxb = 0.0_8
    ELSE
      uhtmaxb = -(tndmax*utfacb/uhtmax**2)
    END IF
    DO k=pver,1,-1
      CALL popinteger4(branch)
      IF (branch .LT. 2) THEN
        x2b = uhtmaxb
        uhtmaxb = 0.0_8
      ELSE
        x2b = 0.0_8
      END IF
      IF (ut(k)**2 + vt(k)**2 .EQ. 0.0) THEN
        tempb0 = 0.0
      ELSE
        tempb0 = x2b/(2.0*sqrt(ut(k)**2+vt(k)**2))
      END IF
      utb(k) = utb(k) + 2*ut(k)*tempb0
      vtb(k) = vtb(k) + 2*vt(k)*tempb0
    END DO
    taub(:, :) = 0.0_8
    xvb = 0.0_8
    yvb = 0.0_8
    DO k=pver,ktop+1,-1
      CALL popinteger4(branch)
      IF (branch .LT. 2) THEN
        ubtb = 0.0_8
      ELSE
        ubtb = effgw(i)*xv*utb(k) + effgw(i)*yv*vtb(k)
        yvb = yvb + effgw(i)*ubt*vtb(k)
        vtb(k) = 0.0_8
        xvb = xvb + effgw(i)*ubt*utb(k)
        utb(k) = 0.0_8
      END IF
      DO l=ngwv,-ngwv,-1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) THEN
          utlb = ubtb
          ubtlb = sign(1.d0, ubtl*(c(l)-ubi(k)))*utlb
          CALL popreal4(ubtl)
          tempb = mapl_grav*rdpm(i, k)*ubtlb
          taub(l, k) = taub(l, k) + tempb
          taub(l, k-1) = taub(l, k-1) - tempb
        END IF
      END DO
      CALL popreal4(ubt)
    END DO
    ubib(:) = 0.0_8
    DO l=ngwv,-ngwv,-1
      DO k=ktop,pver-1,1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) THEN
          tau_minb = taub(l, k)
          taub(l, k) = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) THEN
            taub(l, k+1) = taub(l, k+1) + tau_minb
            tau_minb = 0.0_8
          END IF
          tausatb = tau_minb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.50*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.20*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.10*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.05*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.02*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (branch .LT. 1) THEN
            x1b = -tausatb
          ELSE
            x1b = tausatb
          END IF
          ubmcb = effkwvmap*rhoi(k)*3*ubmc**2*x1b/(2.*ni(k))
          CALL popinteger4(branch)
          IF (branch .LT. 2) THEN
            IF (branch .LT. 1) THEN
              CALL popreal4(effkwvmap)
            ELSE
              CALL popreal4(effkwvmap)
            END IF
          ELSE
            CALL popreal4(effkwvmap)
          END IF
          CALL popreal4(ubmc)
          ubib(k) = ubib(k) + ubmcb
        END IF
      END DO
    END DO
    DO k=pver,1,-1
      vtb(k) = 0.0_8
      utb(k) = 0.0_8
    END DO
    DO l=ngwv,-ngwv,-1
      DO k=ktop,pver-1,1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) taub(l, k) = 0.0_8
      END DO
    END DO
  END SUBROUTINE gw_drag_prof_bwd
!===============================================================================
  SUBROUTINE gw_drag_prof(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, &
&    v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback&
&    , kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut, vt&
&    , tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y&
&    , effgw)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! index of bottom (source) interface
    INTEGER, INTENT(IN) :: kbot
! index of top interface of gwd region
    INTEGER, INTENT(IN) :: ktop
! number of gravity waves possible
    INTEGER, INTENT(IN) :: pgwv
! number of gravity waves to use
    INTEGER, INTENT(IN) :: ngwv
! top interface of low level stress  divergence region
    INTEGER, INTENT(IN) :: kldv
! min value of kldv
    INTEGER, INTENT(IN) :: kldvmn
! index of top interface of source region
    INTEGER, INTENT(IN) :: ksrc
! min value of ksrc
    INTEGER, INTENT(IN) :: ksrcmn
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
! interface density
    REAL :: rhoi(0:pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! interface temperature
    REAL :: ti(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! time step
    REAL :: dt
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
! tendency efficiency for gwd 
    REAL :: effgw(pcols)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! zonal wind tendency
    REAL :: ut(pver)
! meridional wind tendency
    REAL :: vt(pver)
! temperature tendency
    REAL :: tt(pver)
! Total zonal GW momentum flux
    REAL :: taugwx(pcols, 0:pver)
! Total meridional GW momentum flux
    REAL :: taugwy(pcols, 0:pver)
! Total GW energy flux
    REAL :: fegw(pcols, 0:pver)
! Total GW pseudo energy flux
    REAL :: fepgw(pcols, 0:pver)
! Total U tendency below launch level
    REAL :: dusrc(pver)
! Total V tendency below launch level
    REAL :: dvsrc(pver)
! Total V tendency below launch level
    REAL :: dtsrc(pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! fraction of dsat to use
    REAL :: dscal
! imaginary part of vertical wavenumber
    REAL :: mi
! stress after damping
    REAL :: taudmp
!real    :: taumax(pcols)                     ! max(tau) for any l
! saturation stress
    REAL :: tausat
! (ub-c)
    REAL :: ubmc
! (ub-c)**2
    REAL :: ubmc2
! ubar tendency
    REAL :: ubt
! tbar tendency
    REAL :: tbt
! ubar tendency from wave l
    REAL :: ubtl
! saturation tendency
    REAL :: ubtlsat
! layer pressure
    REAL :: pm
! layer density
    REAL :: rhom
! launch level height
    REAL :: zlb
! c-u
    REAL :: cmu
    REAL :: dzm, hscal, tautmp
    REAL :: utl
    REAL :: ttl
! zonal pseudomomentum flux spectrum
    REAL :: fpmx
! meridional pseudomomentum flux spectrum
    REAL :: fpmy
! energy flux (p'w') spectrum
    REAL :: fe
! pseudoenergy flux (p'w'+U rho u'w') spectrum
    REAL :: fpe
    REAL :: fpml
    REAL :: fpmt
    REAL :: fpel
    REAL :: fpet
    REAL :: dusrcl
    REAL :: dvsrcl
    REAL :: dtsrcl
    REAL :: zi
    REAL :: effkwvmap
    REAL :: zfac
    REAL :: uhtmax
    REAL :: utfac
    REAL :: tau_min, cmu_, t_new
    INTRINSIC max
    INTRINSIC sign
    INTRINSIC abs
    REAL :: x2
    REAL :: x1
    INTRINSIC sqrt
! Initialize gravity wave drag tendencies to zero
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) tau(l, k) = 0.
      END DO
    END DO
    DO k=1,pver
      ut(k) = 0.
      vt(k) = 0.
      tt(k) = 0.
      dusrc(k) = 0.
      dvsrc(k) = 0.
      dtsrc(k) = 0.
    END DO
! Initialize total momentum and energy fluxes to zero
    DO k=0,pver
      taugwx(i, k) = 0.
      taugwy(i, k) = 0.
      fegw(i, k) = 0.
      fepgw(i, k) = 0.
    END DO
! Initialize surface wave stress at c = 0 to zero
    tau0x = 0.
    tau0y = 0.
!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------
! 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.
! Loop from bottom to top to get stress profiles
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) THEN
!!!!!!!!jk             d = dback(k)
          ubmc = ubi(k) - c(l)
          IF (ngwv .GT. 0) THEN
            CALL get_effkwvmap_1(effkwvmap, rlat(i))
            IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
&                .LT. 15.0) effkwvmap = fcrit2*kwvbeq
          ELSE
            IF (pi(i, k) .LT. 1000.0) THEN
              zfac = (pi(i, k)/1000.0)**3
            ELSE
              zfac = 1.0
            END IF
            IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
              CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
            ELSE
              CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
            END IF
          END IF
          x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
          IF (x1 .GE. 0.) THEN
            tausat = x1
          ELSE
            tausat = -x1
          END IF
          IF (tausat .LE. taumin) tausat = 0.0
          IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) tausat = 0.0
          IF (k .EQ. ktop) tausat = 0.
!
          IF (k .EQ. ktop + 1) tausat = tausat*0.02
          IF (k .EQ. ktop + 2) tausat = tausat*0.05
          IF (k .EQ. ktop + 3) tausat = tausat*0.10
          IF (k .EQ. ktop + 4) tausat = tausat*0.20
          IF (k .EQ. ktop + 5) tausat = tausat*0.50
!
          tau_min = tausat
          IF (tau_min .GT. tau(l, k+1)) tau_min = tau(l, k+1)
          tau(l, k) = tau_min
        END IF
      END DO
    END DO
!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------
! Accumulate the mean wind tendency over wavenumber.
! Loop over levels from top to bottom
    DO k=ktop+1,pver
      ubt = 0.0
      tbt = 0.0
      DO l=-ngwv,ngwv
        IF (k .LE. kbot) THEN
! Determine the wind tendency including excess stress carried down from above.
          ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
! Calculate the sign of wind tendency
          utl = sign(ubtl, c(l) - ubi(k))
! Accumulate the mean wind tendency over wavenumber.
          ubt = ubt + utl
! Calculate irreversible temperature tendency associated with gravity wave breaking.
! Adding frictional heating associated with the GW momentum forcing
          ttl = (c(l)-ubm(k))*utl/mapl_cp
! Adding frictional heating associated with the GW momentum forcing
          tbt = tbt + ttl
        END IF
      END DO
! Project the mean wind tendency onto the components and scale by "efficiency".
      IF (k .LE. kbot) THEN
        ut(k) = ubt*xv*effgw(i)
        vt(k) = ubt*yv*effgw(i)
        tt(k) = tbt*effgw(i)
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates wind and temperature tendencies below launch level for
! energy and momentum conservation (does nothing for orographic GWs).
!-----------------------------------------------------------------------
!----kimmmmm here
! Calculate launch level height
    zlb = 0.
    DO k=ktop+1,pver
      IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
        pm = (pi(i, k-1)+pi(i, k))*0.5
        CALL get_ti(t_new, t(i, k))
        rhom = pm/(mapl_rgas*t_new)
        zlb = zlb + dpm(i, k)/mapl_grav/rhom
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates energy and momentum flux profiles
!-----------------------------------------------------------------------
    DO l=-ngwv,ngwv
      DO k=ktop,pver
        IF (k .LE. kbot) THEN
          cmu = c(l) - ubi(k)
          CALL get_cmu(cmu_, cmu)
          fpmx = cmu_*tau(l, k)*xv*effgw(i)
          fpmy = cmu_*tau(l, k)*yv*effgw(i)
          fe = cmu*cmu_*tau(l, k)*effgw(i)
          fpe = c(l)*cmu_*tau(l, k)*effgw(i)
          IF (k .EQ. kbot) THEN
            fpml = fpmx*xv + fpmy*yv
            fpel = fpe
          END IF
          IF (k .EQ. ktop) THEN
            fpmt = fpmx*xv + fpmy*yv
            fpet = fpe
          END IF
! Record outputs for GW fluxes
          taugwx(i, k) = taugwx(i, k) + fpmx
          taugwy(i, k) = taugwy(i, k) + fpmy
          fegw(i, k) = fegw(i, k) + fe
          fepgw(i, k) = fepgw(i, k) + fpe
        END IF
      END DO
      DO k=ktop+1,pver
        IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
! Add sub-source wind and temperature tendencies
          pm = (pi(i, k-1)+pi(i, k))*0.5
          CALL get_ti(t_new, t(i, k))
          rhom = pm/(mapl_rgas*t_new)
          dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
          dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
          dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
! Add sub-source wind and temperature tendencies
          dusrc(k) = dusrc(k) + dusrcl
          dvsrc(k) = dvsrc(k) + dvsrcl
          dtsrc(k) = dtsrc(k) + dtsrcl
        END IF
      END DO
    END DO
! For orographic waves, sub-source tendencies are set equal to zero.
    DO k=1,pver
      IF (ngwv .EQ. 0) THEN
        dusrc(k) = 0.0
        dvsrc(k) = 0.0
        dtsrc(k) = 0.0
      END IF
    END DO
!-----------------------------------------------------------------------
! Adjust efficiency factor to prevent unrealistically strong forcing
!-----------------------------------------------------------------------
    uhtmax = 0.0
    utfac = 1.0
    DO k=1,pver
      x2 = sqrt(ut(k)**2 + vt(k)**2)
      IF (x2 .LT. uhtmax) THEN
        uhtmax = uhtmax
      ELSE
        uhtmax = x2
      END IF
    END DO
    IF (uhtmax .GT. tndmax) utfac = tndmax/uhtmax
!jkim    effgw(i) = effgw(i)*utfac
    DO k=1,pver
      ut(k) = ut(k)*utfac
      vt(k) = vt(k)*utfac
      tt(k) = tt(k)*utfac
      dusrc(k) = dusrc(k)*utfac
      dvsrc(k) = dvsrc(k)*utfac
      dtsrc(k) = dtsrc(k)*utfac
    END DO
!-----------------------------------------------------------------------
! Project the c=0 stress (scaled) in the direction of the source wind
! for recording on the output file.
!-----------------------------------------------------------------------
    tau0x = tau(0, kbot)*xv*effgw(i)*utfac
    tau0y = tau(0, kbot)*yv*effgw(i)*utfac
    RETURN
  END SUBROUTINE gw_drag_prof
!  Differentiation of gw_drag_prof_bgnd in reverse (adjoint) mode:
!   gradient, with respect to input variables: xv yv ubi
!   of linear combination of output variables: tau ut vt ubi dvsrc
!                dusrc
  SUBROUTINE gw_drag_prof_bgnd_b(i, pcols, pver, pgwv, ngwv, kbot, ktop&
&    , c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha&
&    , dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taub, ubi, ubib, &
&    ubm, xv, xvb, yv, yvb, ut, utb, vt, vtb, tt, taugwx, taugwy, fegw, &
&    fepgw, dusrc, dusrcb, dvsrc, dvsrcb, dtsrc, tau0x, tau0y, effgw)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! index of bottom (source) interface
    INTEGER, INTENT(IN) :: kbot
! index of top interface of gwd region
    INTEGER, INTENT(IN) :: ktop
! number of gravity waves possible
    INTEGER, INTENT(IN) :: pgwv
! number of gravity waves to use
    INTEGER, INTENT(IN) :: ngwv
! top interface of low level stress  divergence region
    INTEGER, INTENT(IN) :: kldv
! min value of kldv
    INTEGER, INTENT(IN) :: kldvmn
! index of top interface of source region
    INTEGER, INTENT(IN) :: ksrc
! min value of ksrc
    INTEGER, INTENT(IN) :: ksrcmn
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
! interface density
    REAL :: rhoi(0:pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! interface temperature
    REAL :: ti(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! time step
    REAL :: dt
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! projection of wind at interfaces
    REAL :: ubi(0:pver)
    REAL :: ubib(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
    REAL :: xvb
! unit vectors of source wind (y)
    REAL :: yv
    REAL :: yvb
! tendency efficiency for gwd 
    REAL :: effgw(pcols)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
    REAL :: taub(-pgwv:pgwv, 0:pver)
! zonal wind tendency
    REAL :: ut(pver)
    REAL :: utb(pver)
! meridional wind tendency
    REAL :: vt(pver)
    REAL :: vtb(pver)
! temperature tendency
    REAL :: tt(pver)
! Total zonal GW momentum flux
    REAL :: taugwx(pcols, 0:pver)
! Total meridional GW momentum flux
    REAL :: taugwy(pcols, 0:pver)
! Total GW energy flux
    REAL :: fegw(pcols, 0:pver)
! Total GW pseudo energy flux
    REAL :: fepgw(pcols, 0:pver)
! Total U tendency below launch level
    REAL :: dusrc(pver)
    REAL :: dusrcb(pver)
! Total V tendency below launch level
    REAL :: dvsrc(pver)
    REAL :: dvsrcb(pver)
! Total V tendency below launch level
    REAL :: dtsrc(pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! fraction of dsat to use
    REAL :: dscal
! imaginary part of vertical wavenumber
    REAL :: mi
! stress after damping
    REAL :: taudmp
!real    :: taumax(pcols)                     ! max(tau) for any l
! saturation stress
    REAL :: tausat
    REAL :: tausatb
! (ub-c)
    REAL :: ubmc
    REAL :: ubmcb
! (ub-c)**2
    REAL :: ubmc2
! ubar tendency
    REAL :: ubt
    REAL :: ubtb
! tbar tendency
    REAL :: tbt
! ubar tendency from wave l
    REAL :: ubtl
    REAL :: ubtlb
! saturation tendency
    REAL :: ubtlsat
! layer pressure
    REAL :: pm
! layer density
    REAL :: rhom
! launch level height
    REAL :: zlb
! c-u
    REAL :: cmu
    REAL :: dzm, hscal, tautmp
    REAL :: utl
    REAL :: utlb
    REAL :: ttl
! zonal pseudomomentum flux spectrum
    REAL :: fpmx
    REAL :: fpmxb
! meridional pseudomomentum flux spectrum
    REAL :: fpmy
    REAL :: fpmyb
! energy flux (p'w') spectrum
    REAL :: fe
! pseudoenergy flux (p'w'+U rho u'w') spectrum
    REAL :: fpe
    REAL :: fpml
    REAL :: fpmlb
    REAL :: fpmt
    REAL :: fpmtb
    REAL :: fpel
    REAL :: fpet
    REAL :: dusrcl
    REAL :: dusrclb
    REAL :: dvsrcl
    REAL :: dvsrclb
    REAL :: dtsrcl
    REAL :: zi
    REAL :: effkwvmap
    REAL :: zfac
    REAL :: uhtmax
    REAL :: uhtmaxb
    REAL :: utfac
    REAL :: utfacb
    REAL :: tau_min, cmu_, t_new
    REAL :: tau_minb
    INTEGER :: branch
    REAL :: tempb4
    REAL :: tempb3
    REAL :: tempb2
    REAL :: tempb1
    REAL :: tempb0
    INTRINSIC max
    INTRINSIC sign
    INTRINSIC abs
    REAL :: x2
    REAL :: x1
    REAL :: x2b
    REAL :: tempb
    REAL :: x1b
    INTRINSIC sqrt
! Initialize gravity wave drag tendencies to zero
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) tau(l, k) = 0.
      END DO
    END DO
    DO k=1,pver
      ut(k) = 0.
      vt(k) = 0.
      dusrc(k) = 0.
      dvsrc(k) = 0.
    END DO
!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------
! 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.
! Loop from bottom to top to get stress profiles
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) THEN
          CALL pushreal4(ubmc)
!!!!!!!!jk             d = dback(k)
          ubmc = ubi(k) - c(l)
          IF (ngwv .GT. 0) THEN
            CALL pushreal4(effkwvmap)
            CALL get_effkwvmap_1(effkwvmap, rlat(i))
            IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
&                .LT. 15.0) THEN
              CALL pushinteger4(0)
              effkwvmap = fcrit2*kwvbeq
            ELSE
              CALL pushinteger4(0)
            END IF
          ELSE
            IF (pi(i, k) .LT. 1000.0) THEN
              zfac = (pi(i, k)/1000.0)**3
            ELSE
              zfac = 1.0
            END IF
            IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
              CALL pushreal4(effkwvmap)
              CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
              CALL pushinteger4(1)
            ELSE
              CALL pushreal4(effkwvmap)
              CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
              CALL pushinteger4(2)
            END IF
          END IF
          x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
          IF (x1 .GE. 0.) THEN
            tausat = x1
            CALL pushinteger4(1)
          ELSE
            tausat = -x1
            CALL pushinteger4(0)
          END IF
          IF (tausat .LE. taumin) THEN
            tausat = 0.0
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) THEN
            tausat = 0.0
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop) THEN
            tausat = 0.
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
!
          IF (k .EQ. ktop + 1) THEN
            tausat = tausat*0.02
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 2) THEN
            tausat = tausat*0.05
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 3) THEN
            tausat = tausat*0.10
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 4) THEN
            tausat = tausat*0.20
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop + 5) THEN
            tausat = tausat*0.50
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
!
          tau_min = tausat
          IF (tau_min .GT. tau(l, k+1)) THEN
            tau_min = tau(l, k+1)
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          tau(l, k) = tau_min
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
    END DO
!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------
! Accumulate the mean wind tendency over wavenumber.
! Loop over levels from top to bottom
    DO k=ktop+1,pver
      CALL pushreal4(ubt)
      ubt = 0.0
      DO l=-ngwv,ngwv
        IF (k .LE. kbot) THEN
          CALL pushreal4(ubtl)
! Determine the wind tendency including excess stress carried down from above.
          ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
! Calculate the sign of wind tendency
          utl = sign(ubtl, c(l) - ubi(k))
! Accumulate the mean wind tendency over wavenumber.
          ubt = ubt + utl
! Calculate irreversible temperature tendency associated with gravity wave breaking.
! Adding frictional heating associated with the GW momentum forcing
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
! Project the mean wind tendency onto the components and scale by "efficiency".
      IF (k .LE. kbot) THEN
        ut(k) = ubt*xv*effgw(i)
        vt(k) = ubt*yv*effgw(i)
        CALL pushinteger4(2)
      ELSE
        CALL pushinteger4(1)
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates wind and temperature tendencies below launch level for
! energy and momentum conservation (does nothing for orographic GWs).
!-----------------------------------------------------------------------
!----kimmmmm here
! Calculate launch level height
    zlb = 0.
    DO k=ktop+1,pver
      IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
        pm = (pi(i, k-1)+pi(i, k))*0.5
        CALL get_ti(t_new, t(i, k))
        rhom = pm/(mapl_rgas*t_new)
        zlb = zlb + dpm(i, k)/mapl_grav/rhom
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates energy and momentum flux profiles
!-----------------------------------------------------------------------
    DO l=-ngwv,ngwv
      DO k=ktop,pver
        IF (k .LE. kbot) THEN
          cmu = c(l) - ubi(k)
          CALL pushreal4(cmu_)
          CALL get_cmu(cmu_, cmu)
          CALL pushreal4(fpmx)
          fpmx = cmu_*tau(l, k)*xv*effgw(i)
          CALL pushreal4(fpmy)
          fpmy = cmu_*tau(l, k)*yv*effgw(i)
          IF (k .EQ. kbot) THEN
            CALL pushreal4(fpml)
            fpml = fpmx*xv + fpmy*yv
            CALL pushinteger4(1)
          ELSE
            CALL pushinteger4(0)
          END IF
          IF (k .EQ. ktop) THEN
            CALL pushreal4(fpmt)
            fpmt = fpmx*xv + fpmy*yv
            CALL pushinteger4(3)
          ELSE
            CALL pushinteger4(2)
          END IF
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
      DO k=ktop+1,pver
        IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
          pm = (pi(i, k-1)+pi(i, k))*0.5
          CALL get_ti(t_new, t(i, k))
          CALL pushreal4(rhom)
          rhom = pm/(mapl_rgas*t_new)
          dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
          dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
! Add sub-source wind and temperature tendencies
          dusrc(k) = dusrc(k) + dusrcl
          dvsrc(k) = dvsrc(k) + dvsrcl
          CALL pushinteger4(2)
        ELSE
          CALL pushinteger4(1)
        END IF
      END DO
    END DO
! For orographic waves, sub-source tendencies are set equal to zero.
    DO k=1,pver
      IF (ngwv .EQ. 0) THEN
        dusrc(k) = 0.0
        dvsrc(k) = 0.0
        CALL pushinteger4(2)
      ELSE
        CALL pushinteger4(1)
      END IF
    END DO
!-----------------------------------------------------------------------
! Adjust efficiency factor to prevent unrealistically strong forcing
!-----------------------------------------------------------------------
    uhtmax = 0.0
    utfac = 1.0
    DO k=1,pver
      x2 = sqrt(ut(k)**2 + vt(k)**2)
      IF (x2 .LT. uhtmax) THEN
        uhtmax = uhtmax
        CALL pushinteger4(2)
      ELSE
        uhtmax = x2
        CALL pushinteger4(1)
      END IF
    END DO
    IF (uhtmax .GT. tndmax) THEN
      utfac = tndmax/uhtmax
      CALL pushinteger4(1)
    ELSE
      CALL pushinteger4(0)
    END IF
    utfacb = 0.0_8
    DO k=pver,1,-1
      utfacb = utfacb + dusrc(k)*dusrcb(k) + ut(k)*utb(k) + vt(k)*vtb(k)&
&        + dvsrc(k)*dvsrcb(k)
      dvsrcb(k) = utfac*dvsrcb(k)
      dusrcb(k) = utfac*dusrcb(k)
      vtb(k) = utfac*vtb(k)
      utb(k) = utfac*utb(k)
    END DO
    CALL popinteger4(branch)
    IF (branch .LT. 1) THEN
      uhtmaxb = 0.0_8
    ELSE
      uhtmaxb = -(tndmax*utfacb/uhtmax**2)
    END IF
    DO k=pver,1,-1
      CALL popinteger4(branch)
      IF (branch .LT. 2) THEN
        x2b = uhtmaxb
        uhtmaxb = 0.0_8
      ELSE
        x2b = 0.0_8
      END IF
      IF (ut(k)**2 + vt(k)**2 .EQ. 0.0) THEN
        tempb4 = 0.0
      ELSE
        tempb4 = x2b/(2.0*sqrt(ut(k)**2+vt(k)**2))
      END IF
      utb(k) = utb(k) + 2*ut(k)*tempb4
      vtb(k) = vtb(k) + 2*vt(k)*tempb4
    END DO
    DO k=pver,1,-1
      CALL popinteger4(branch)
      IF (.NOT.branch .LT. 2) THEN
        dvsrcb(k) = 0.0_8
        dusrcb(k) = 0.0_8
      END IF
    END DO
    xvb = 0.0_8
    yvb = 0.0_8
    fpmlb = 0.0_8
    fpmtb = 0.0_8
    DO l=ngwv,-ngwv,-1
      DO k=pver,ktop+1,-1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) THEN
          dvsrclb = dvsrcb(k)
          dusrclb = dusrcb(k)
          tempb2 = -(dvsrclb/(rhom*zlb))
          tempb3 = -(dusrclb/(rhom*zlb))
          fpmlb = fpmlb + xv*tempb3 + yv*tempb2
          fpmtb = fpmtb - xv*tempb3 - yv*tempb2
          yvb = yvb + (fpml-fpmt)*tempb2
          xvb = xvb + (fpml-fpmt)*tempb3
          CALL popreal4(rhom)
        END IF
      END DO
      DO k=pver,ktop,-1
        CALL popinteger4(branch)
        IF (branch .LT. 3) THEN
          IF (branch .LT. 2) THEN
            GOTO 100
          ELSE
            fpmxb = 0.0_8
            fpmyb = 0.0_8
          END IF
        ELSE
          CALL popreal4(fpmt)
          fpmxb = xv*fpmtb
          xvb = xvb + fpmx*fpmtb
          fpmyb = yv*fpmtb
          yvb = yvb + fpmy*fpmtb
          fpmtb = 0.0_8
        END IF
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 1) THEN
          CALL popreal4(fpml)
          fpmxb = fpmxb + xv*fpmlb
          xvb = xvb + fpmx*fpmlb
          fpmyb = fpmyb + yv*fpmlb
          yvb = yvb + fpmy*fpmlb
          fpmlb = 0.0_8
        END IF
        CALL popreal4(fpmy)
        tempb0 = cmu_*effgw(i)*fpmyb
        tempb1 = cmu_*effgw(i)*fpmxb
        taub(l, k) = taub(l, k) + xv*tempb1 + yv*tempb0
        yvb = yvb + tau(l, k)*tempb0
        CALL popreal4(fpmx)
        xvb = xvb + tau(l, k)*tempb1
        CALL popreal4(cmu_)
 100    CONTINUE
      END DO
    END DO
    DO k=pver,ktop+1,-1
      CALL popinteger4(branch)
      IF (branch .LT. 2) THEN
        ubtb = 0.0_8
      ELSE
        ubtb = effgw(i)*xv*utb(k) + effgw(i)*yv*vtb(k)
        yvb = yvb + effgw(i)*ubt*vtb(k)
        vtb(k) = 0.0_8
        xvb = xvb + effgw(i)*ubt*utb(k)
        utb(k) = 0.0_8
      END IF
      DO l=ngwv,-ngwv,-1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) THEN
          utlb = ubtb
          ubtlb = sign(1.d0, ubtl*(c(l)-ubi(k)))*utlb
          CALL popreal4(ubtl)
          tempb = mapl_grav*rdpm(i, k)*ubtlb
          taub(l, k) = taub(l, k) + tempb
          taub(l, k-1) = taub(l, k-1) - tempb
        END IF
      END DO
      CALL popreal4(ubt)
    END DO
    DO l=ngwv,-ngwv,-1
      DO k=ktop,pver-1,1
        CALL popinteger4(branch)
        IF (.NOT.branch .LT. 2) THEN
          tau_minb = taub(l, k)
          taub(l, k) = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) THEN
            taub(l, k+1) = taub(l, k+1) + tau_minb
            tau_minb = 0.0_8
          END IF
          tausatb = tau_minb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.50*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.20*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.10*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.05*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.02*tausatb
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (.NOT.branch .LT. 1) tausatb = 0.0_8
          CALL popinteger4(branch)
          IF (branch .LT. 1) THEN
            x1b = -tausatb
          ELSE
            x1b = tausatb
          END IF
          ubmcb = effkwvmap*rhoi(k)*3*ubmc**2*x1b/(2.*ni(k))
          CALL popinteger4(branch)
          IF (branch .LT. 2) THEN
            IF (branch .LT. 1) THEN
              CALL popreal4(effkwvmap)
            ELSE
              CALL popreal4(effkwvmap)
            END IF
          ELSE
            CALL popreal4(effkwvmap)
          END IF
          CALL popreal4(ubmc)
          ubib(k) = ubib(k) + ubmcb
        END IF
      END DO
    END DO
  END SUBROUTINE gw_drag_prof_bgnd_b
  SUBROUTINE gw_drag_prof_bgnd(i, pcols, pver, pgwv, ngwv, kbot, ktop, c&
&    , u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, &
&    dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut&
&    , vt, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, &
&    tau0y, effgw)
    IMPLICIT NONE
!-----------------------------------------------------------------------
! 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--------------------------------
! number of current column
    INTEGER, INTENT(IN) :: i
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pcols
! number of atmospheric columns
    INTEGER, INTENT(IN) :: pver
! index of bottom (source) interface
    INTEGER, INTENT(IN) :: kbot
! index of top interface of gwd region
    INTEGER, INTENT(IN) :: ktop
! number of gravity waves possible
    INTEGER, INTENT(IN) :: pgwv
! number of gravity waves to use
    INTEGER, INTENT(IN) :: ngwv
! top interface of low level stress  divergence region
    INTEGER, INTENT(IN) :: kldv
! min value of kldv
    INTEGER, INTENT(IN) :: kldvmn
! index of top interface of source region
    INTEGER, INTENT(IN) :: ksrc
! min value of ksrc
    INTEGER, INTENT(IN) :: ksrcmn
! wave phase speeds
    REAL :: c(-pgwv:pgwv)
! midpoint zonal wind
    REAL :: u(pcols, pver)
! midpoint meridional wind
    REAL :: v(pcols, pver)
! midpoint temperatures
    REAL :: t(pcols, pver)
! interface pressures
    REAL :: pi(pcols, 0:pver)
! midpoint delta p (pi(k)-pi(k-1))
    REAL :: dpm(pcols, pver)
! 1. / (pi(k)-pi(k-1))
    REAL :: rdpm(pcols, pver)
! ln(interface pressures)
    REAL :: piln(pcols, 0:pver)
    REAL :: rlat(pcols)
! interface density
    REAL :: rhoi(0:pver)
! interface Brunt-Vaisalla frequency
    REAL :: ni(0:pver)
! interface temperature
    REAL :: ti(0:pver)
! midpoint Brunt-Vaisalla frequency
    REAL :: nm(pver)
! time step
    REAL :: dt
! newtonian cooling coefficients
    REAL :: alpha(0:pver)
! newtonian cooling coefficients
    REAL :: dback(0:pver)
! 1/dp across low level divergence region
    REAL :: rdpldv
! projection of wind at interfaces
    REAL :: ubi(0:pver)
! projection of wind at midpoints
    REAL :: ubm(pver)
! unit vectors of source wind (x)
    REAL :: xv
! unit vectors of source wind (y)
    REAL :: yv
! tendency efficiency for gwd 
    REAL :: effgw(pcols)
! wave Reynolds stress
    REAL :: tau(-pgwv:pgwv, 0:pver)
! zonal wind tendency
    REAL :: ut(pver)
! meridional wind tendency
    REAL :: vt(pver)
! temperature tendency
    REAL :: tt(pver)
! Total zonal GW momentum flux
    REAL :: taugwx(pcols, 0:pver)
! Total meridional GW momentum flux
    REAL :: taugwy(pcols, 0:pver)
! Total GW energy flux
    REAL :: fegw(pcols, 0:pver)
! Total GW pseudo energy flux
    REAL :: fepgw(pcols, 0:pver)
! Total U tendency below launch level
    REAL :: dusrc(pver)
! Total V tendency below launch level
    REAL :: dvsrc(pver)
! Total V tendency below launch level
    REAL :: dtsrc(pver)
! c=0 sfc. stress (zonal)
    REAL :: tau0x
! c=0 sfc. stress (meridional)
    REAL :: tau0y
!---------------------------Local storage-------------------------------
! loop indexes
    INTEGER :: k, l
! fraction of dsat to use
    REAL :: dscal
! imaginary part of vertical wavenumber
    REAL :: mi
! stress after damping
    REAL :: taudmp
!real    :: taumax(pcols)                     ! max(tau) for any l
! saturation stress
    REAL :: tausat
! (ub-c)
    REAL :: ubmc
! (ub-c)**2
    REAL :: ubmc2
! ubar tendency
    REAL :: ubt
! tbar tendency
    REAL :: tbt
! ubar tendency from wave l
    REAL :: ubtl
! saturation tendency
    REAL :: ubtlsat
! layer pressure
    REAL :: pm
! layer density
    REAL :: rhom
! launch level height
    REAL :: zlb
! c-u
    REAL :: cmu
    REAL :: dzm, hscal, tautmp
    REAL :: utl
    REAL :: ttl
! zonal pseudomomentum flux spectrum
    REAL :: fpmx
! meridional pseudomomentum flux spectrum
    REAL :: fpmy
! energy flux (p'w') spectrum
    REAL :: fe
! pseudoenergy flux (p'w'+U rho u'w') spectrum
    REAL :: fpe
    REAL :: fpml
    REAL :: fpmt
    REAL :: fpel
    REAL :: fpet
    REAL :: dusrcl
    REAL :: dvsrcl
    REAL :: dtsrcl
    REAL :: zi
    REAL :: effkwvmap
    REAL :: zfac
    REAL :: uhtmax
    REAL :: utfac
    REAL :: tau_min, cmu_, t_new
    INTRINSIC max
    INTRINSIC sign
    INTRINSIC abs
    REAL :: x2
    REAL :: x1
    INTRINSIC sqrt
! Initialize gravity wave drag tendencies to zero
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) tau(l, k) = 0.
      END DO
    END DO
    DO k=1,pver
      ut(k) = 0.
      vt(k) = 0.
      tt(k) = 0.
      dusrc(k) = 0.
      dvsrc(k) = 0.
      dtsrc(k) = 0.
    END DO
! Initialize total momentum and energy fluxes to zero
    DO k=0,pver
      taugwx(i, k) = 0.
      taugwy(i, k) = 0.
      fegw(i, k) = 0.
      fepgw(i, k) = 0.
    END DO
! Initialize surface wave stress at c = 0 to zero
    tau0x = 0.
    tau0y = 0.
!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------
! 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.
! Loop from bottom to top to get stress profiles
    DO l=-ngwv,ngwv
      DO k=pver-1,ktop,-1
        IF (k .LE. kbot - 1) THEN
!!!!!!!!jk             d = dback(k)
          ubmc = ubi(k) - c(l)
          IF (ngwv .GT. 0) THEN
            CALL get_effkwvmap_1(effkwvmap, rlat(i))
            IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
&                .LT. 15.0) effkwvmap = fcrit2*kwvbeq
          ELSE
            IF (pi(i, k) .LT. 1000.0) THEN
              zfac = (pi(i, k)/1000.0)**3
            ELSE
              zfac = 1.0
            END IF
            IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
              CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
            ELSE
              CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
            END IF
          END IF
          x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
          IF (x1 .GE. 0.) THEN
            tausat = x1
          ELSE
            tausat = -x1
          END IF
          IF (tausat .LE. taumin) tausat = 0.0
          IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) tausat = 0.0
          IF (k .EQ. ktop) tausat = 0.
!
          IF (k .EQ. ktop + 1) tausat = tausat*0.02
          IF (k .EQ. ktop + 2) tausat = tausat*0.05
          IF (k .EQ. ktop + 3) tausat = tausat*0.10
          IF (k .EQ. ktop + 4) tausat = tausat*0.20
          IF (k .EQ. ktop + 5) tausat = tausat*0.50
!
          tau_min = tausat
          IF (tau_min .GT. tau(l, k+1)) tau_min = tau(l, k+1)
          tau(l, k) = tau_min
        END IF
      END DO
    END DO
!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------
! Accumulate the mean wind tendency over wavenumber.
! Loop over levels from top to bottom
    DO k=ktop+1,pver
      ubt = 0.0
      tbt = 0.0
      DO l=-ngwv,ngwv
        IF (k .LE. kbot) THEN
! Determine the wind tendency including excess stress carried down from above.
          ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
! Calculate the sign of wind tendency
          utl = sign(ubtl, c(l) - ubi(k))
! Accumulate the mean wind tendency over wavenumber.
          ubt = ubt + utl
! Calculate irreversible temperature tendency associated with gravity wave breaking.
! Adding frictional heating associated with the GW momentum forcing
          ttl = (c(l)-ubm(k))*utl/mapl_cp
! Adding frictional heating associated with the GW momentum forcing
          tbt = tbt + ttl
        END IF
      END DO
! Project the mean wind tendency onto the components and scale by "efficiency".
      IF (k .LE. kbot) THEN
        ut(k) = ubt*xv*effgw(i)
        vt(k) = ubt*yv*effgw(i)
        tt(k) = tbt*effgw(i)
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates wind and temperature tendencies below launch level for
! energy and momentum conservation (does nothing for orographic GWs).
!-----------------------------------------------------------------------
!----kimmmmm here
! Calculate launch level height
    zlb = 0.
    DO k=ktop+1,pver
      IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
        pm = (pi(i, k-1)+pi(i, k))*0.5
        CALL get_ti(t_new, t(i, k))
        rhom = pm/(mapl_rgas*t_new)
        zlb = zlb + dpm(i, k)/mapl_grav/rhom
      END IF
    END DO
!-----------------------------------------------------------------------
! Calculates energy and momentum flux profiles
!-----------------------------------------------------------------------
    DO l=-ngwv,ngwv
      DO k=ktop,pver
        IF (k .LE. kbot) THEN
          cmu = c(l) - ubi(k)
          CALL get_cmu(cmu_, cmu)
          fpmx = cmu_*tau(l, k)*xv*effgw(i)
          fpmy = cmu_*tau(l, k)*yv*effgw(i)
          fe = cmu*cmu_*tau(l, k)*effgw(i)
          fpe = c(l)*cmu_*tau(l, k)*effgw(i)
          IF (k .EQ. kbot) THEN
            fpml = fpmx*xv + fpmy*yv
            fpel = fpe
          END IF
          IF (k .EQ. ktop) THEN
            fpmt = fpmx*xv + fpmy*yv
            fpet = fpe
          END IF
! Record outputs for GW fluxes
          taugwx(i, k) = taugwx(i, k) + fpmx
          taugwy(i, k) = taugwy(i, k) + fpmy
          fegw(i, k) = fegw(i, k) + fe
          fepgw(i, k) = fepgw(i, k) + fpe
        END IF
      END DO
      DO k=ktop+1,pver
        IF (k .GE. kbot + 1) THEN
! Define layer pressure and density
          pm = (pi(i, k-1)+pi(i, k))*0.5
          CALL get_ti(t_new, t(i, k))
          rhom = pm/(mapl_rgas*t_new)
          dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
          dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
          dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
! Add sub-source wind and temperature tendencies
          dusrc(k) = dusrc(k) + dusrcl
          dvsrc(k) = dvsrc(k) + dvsrcl
          dtsrc(k) = dtsrc(k) + dtsrcl
        END IF
      END DO
    END DO
! For orographic waves, sub-source tendencies are set equal to zero.
    DO k=1,pver
      IF (ngwv .EQ. 0) THEN
        dusrc(k) = 0.0
        dvsrc(k) = 0.0
        dtsrc(k) = 0.0
      END IF
    END DO
!-----------------------------------------------------------------------
! Adjust efficiency factor to prevent unrealistically strong forcing
!-----------------------------------------------------------------------
    uhtmax = 0.0
    utfac = 1.0
    DO k=1,pver
      x2 = sqrt(ut(k)**2 + vt(k)**2)
      IF (x2 .LT. uhtmax) THEN
        uhtmax = uhtmax
      ELSE
        uhtmax = x2
      END IF
    END DO
    IF (uhtmax .GT. tndmax) utfac = tndmax/uhtmax
!jkim    effgw(i) = effgw(i)*utfac
    DO k=1,pver
      ut(k) = ut(k)*utfac
      vt(k) = vt(k)*utfac
      tt(k) = tt(k)*utfac
      dusrc(k) = dusrc(k)*utfac
      dvsrc(k) = dvsrc(k)*utfac
      dtsrc(k) = dtsrc(k)*utfac
    END DO
!-----------------------------------------------------------------------
! Project the c=0 stress (scaled) in the direction of the source wind
! for recording on the output file.
!-----------------------------------------------------------------------
    tau0x = tau(0, kbot)*xv*effgw(i)*utfac
    tau0y = tau(0, kbot)*yv*effgw(i)*utfac
    RETURN
  END SUBROUTINE gw_drag_prof_bgnd
!  Differentiation of get_uv in reverse (adjoint) mode:
!   gradient, with respect to input variables: uv_in
!   of linear combination of output variables: uv_out uv_in
  SUBROUTINE get_uv_b(uv_out, uv_outb, uv_in, uv_inb)
    IMPLICIT NONE
    REAL :: uv_out, uv_in
    REAL :: uv_outb, uv_inb
    LOGICAL :: src_bypass
    IF (src_bypass) THEN
      uv_inb = uv_inb + uv_outb
    ELSE
      uv_inb = uv_inb + uv_outb
    END IF
  END SUBROUTINE get_uv_b
  SUBROUTINE get_uv(uv_out, uv_in)
    IMPLICIT NONE
    REAL :: uv_out, uv_in
    LOGICAL :: src_bypass
    IF (src_bypass) THEN
      uv_out = uv_in
    ELSE
      uv_out = uv_in
    END IF
  END SUBROUTINE get_uv
  SUBROUTINE get_effkwvmap_1(effkwvmap_, rlat_)
    IMPLICIT NONE
    REAL :: tanh1, tanh2, effkwvmap_, rlat_
    INTRINSIC tanh
    tanh1 = (rlat_*180./pi_gwd-20.0)/6.0
    tanh2 = -((rlat_*180./pi_gwd+20.0)/6.0)
    tanh1 = tanh(tanh1)
    tanh2 = tanh(tanh2)
    effkwvmap_ = fcrit2*kwvb*(0.5*(1.0+tanh1)+0.5*(1.0+tanh2))
  END SUBROUTINE get_effkwvmap_1
  SUBROUTINE get_effkwvmap_2(effkwvmap_, rlat_, zfac)
    IMPLICIT NONE
    REAL :: effkwvmap_, rlat_, zfac
    INTRINSIC exp
    INTRINSIC abs
    REAL :: abs1
    IF (rlat_*180./pi_gwd + 20. .GE. 0.) THEN
      abs1 = rlat_*180./pi_gwd + 20.
    ELSE
      abs1 = -(rlat_*180./pi_gwd+20.)
    END IF
    effkwvmap_ = fcrit2*kwvo*zfac*(2.0-0.65*exp(-(abs1/5.)))
  END SUBROUTINE get_effkwvmap_2
  SUBROUTINE get_effkwvmap_3(effkwvmap_, rlat_, zfac)
    IMPLICIT NONE
    REAL :: effkwvmap_, rlat_, zfac
    INTRINSIC exp
    INTRINSIC abs
    REAL :: abs1
    IF (rlat_*180./pi_gwd + 20. .GE. 0.) THEN
      abs1 = rlat_*180./pi_gwd + 20.
    ELSE
      abs1 = -(rlat_*180./pi_gwd+20.)
    END IF
    effkwvmap_ = fcrit2*kwvo*zfac*(0.7+0.65*exp(-(abs1/5.)))
  END SUBROUTINE get_effkwvmap_3
  SUBROUTINE get_cmu(cmu_, cmu)
    IMPLICIT NONE
    REAL :: cmu_, cmu
    INTRINSIC sign
    cmu_ = sign(1.0, cmu)
  END SUBROUTINE get_cmu
  SUBROUTINE get_ti(t_out, t_in)
    IMPLICIT NONE
    REAL :: t_out, t_in
    LOGICAL :: src_bypass
    t_out = t_in
  END SUBROUTINE get_ti
END MODULE gw_drag_b