#include "dims.h"
      subroutine mgwbgnd (u, v, t, pm, pi, dpm, rdpm, piln, rlat,
     $     kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv,
     $     ngwv, kbot,jlat)
C-----------------------------------------------------------------------
C 
C Driver for multiple gravity wave drag parameterization.
C 
C The parameterization is assumed to operate only where water vapor 
C concentrations are negligible in determining the density.
C 
C-----------------------------------------------------------------------
      implicit none
C-----------------------------------------------------------------------
#include "pmgrid.h"
C-----------------------------------------------------------------------
#include "mgw.h"
C-----------------------------------------------------------------------
C 
C Input variables
C 
      integer
     $     kbot,                        ! index of bottom (source) interface
     $     ngwv,                        ! number of gravity waves to use
     $     jlat,                        ! latitude index
     $     idir(plond)                  ! wind longitude direction index(z=z0)
                                        ! idir=1, westerly, idir, easterly
      real 
     $     u(plond,plev),               ! midpoint zonal wind
     $     v(plond,plev),               ! midpoint meridional wind
     $     t(plond,plev),               ! midpoint temperatures
     $     pm(plond,plev),              ! midpoint pressures
     $     pi(plond,0:plev),            ! interface pressures
     $     dpm(plond,plev),             ! midpoint delta p (pi(k)-pi(k-1))
     $     rdpm(plond,plev),            ! 1. / (pi(k)-pi(k-1))
     $     piln(plond,0:plev),          ! ln(interface pressures)
     $     taueq,                       ! Background tau at equator
     $     taueqw,                      ! taueq width in radians
     $     rlat,                        ! current latitude in radians
     $     orog(plon+1,plat),           ! height of orography in meters
     $     sghdat(plon+1,plat)          ! surface roughness of orography
C 
C Output variables
C 
      integer
     $     kldv(plond),                 ! index of top interface of low level stress 
     $                                  ! divergence region
     $     kldvmn,                      ! min value of kldv
     $     ksrc(plond),                 ! index of top interface of source region
     $     ksrcmn                       ! min value of ksrc

      real
     $     rdpldv(plond),               ! 1/dp across low level divergence region
     $     tau(plond,-pgwv:pgwv,0:plev), ! wave Reynolds stress
     $     ubi(plond,0:plev),           ! projection of wind at interfaces
     $     ubm(plond,plev),             ! projection of wind at midpoints
     $     xv(plond),                   ! unit vectors of source wind (x)
     $     yv(plond)                    ! unit vectors of source wind (y)
C 
C Local workspace
C 
      integer
     $     i,k,l,l1                        ! loop indexes

      real
     $     tauback,                     ! background stress at c=0
     $     tauback1,                    ! background stress at c=0
     $     usrc(plond),                 ! u wind averaged over source region
     $     vsrc(plond),                 ! v wind averaged over source region
     $     al0,                         ! Used in lat dependence of GW spec. 
     $     dlat0,                       ! Used in lat dependence of GW spec.
     $     al1,       		        ! Used in lat dependence of GW spec.
     $	   dlat1,                       ! Used in lat dependence of GW spec.
     $     flat_gw,                     ! The actual lat dependence of GW spec.
     $     pi_g,                        ! 3.14........
     $     fac,                         ! gravity wave forcing
     $     fac1,                        ! gravity wave forcing wu
     $     flt,                         ! local time
     $     xlth,                        ! gravity wave diurnal variation
     $     orogenh,                     ! Mountain enhancement of gravity waves
     $     wugws(36),                   ! Wu/Waters GW distribution
     $     tautmp,
     $     wspd,
     $     whw,
     $     whwe,
     $     wspde
      common/dat/ orog,sghdat
      WUGWS=(/0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.10,
     1 0.35,0.50,0.70,0.85,0.90,0.95,1.0,1.0,1.0,1.0,1.0,
     2 1.0,0.95,0.90,0.85,0.70,0.5,0.35,0.10,0.00,0.00,0.00,
     3 0.00,0.00,0.00,0.00,0.00/)

C 
C---------------------------------------------------------------------------
C Determine the source layer wind and unit vectors, then project winds.
C---------------------------------------------------------------------------
C 
C Just use the source level interface values for the source
c wind speed and direction (unit vector).
C 
      idir = 1
      k = kbot
      do i = 1, plon
         ksrc(i) = k
         kldv(i) = k
	 usrc(i) = 0.5*(u(i,k-1)+u(i,k))
	 vsrc(i) = 0.5*(v(i,k-1)+v(i,k))
c         usrc(i) = u(i,k)
c         vsrc(i) = v(i,k)
         ubi(i,kbot) = sqrt (usrc(i)**2 + vsrc(i)**2)
         xv(i) = usrc(i) / ubi(i,k)
         yv(i) = vsrc(i) / ubi(i,k)
         if(xv(i).le.-0.1) idir(i) = -1
      end do
C 
C Project the local wind at midpoints onto the source wind.
C 
      do k = 1, kbot
         do i = 1, plon
            ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
         end do
      end do

C 
C Compute the interface wind projection by averaging the midpoint winds.
C Use the top level wind at the top interface.
C 
      do i = 1, plon
         ubi(i,0) = ubm(i,1)
      end do
      do k = 1, kbot-1
         do i = 1, plon
            ubi(i,k) = 0.5 * (ubm(i,k) + ubm(i,k+1))
         end do
      end do
c +
c-----------------------------------------------------------------------
c Gravity wave sources
c-----------------------------------------------------------------------
c -
c +
c Determine the background stress at c=0
c -
c     fac = 0.120
      fac = faca
c
c     fac1= faca*0.5
      fac1= faca*0.75
c     fac1= faca
c     fac1= faca/2.
c     fac1= 0.00
c
      tauback = taubgnd * tauscal*fac
c
! 01/04/01: taueqw and taueq as per kibo12:
      whwe = 30.
      wspde = 0.
      taueqw = 15.*3.14159/180.
      taueq = taubgnd * tauscal*fac1*exp(-(rlat/taueqw)**2)
c     taueq = 0.
c
! 12/11/00: eliminate tauback1 as per kibo12
!     tauback1 = taubgnd * tauscal*fac1
c++
c	Include dependence on latitude:
c       The lat function was obtained from RR Garcia as
c	currently used in his 2D model
c       [Added by F. Sassi on May 30, 1997]
c--
        pi_g=4.*atan(1.)
        al0=20.*pi_g/180.
        dlat0=10.*pi_g/180.
        al1=55.*pi_g/180.
        dlat1=70.*pi_g/180.
! 12/11/00: flat_gw commented out as per kibo12
!       flat_gw=
!    $  (0.5*(1.+tanh((rlat-al0)/dlat0))*facnh
!    $  +0.5*(1.+tanh(-(rlat+al0)/dlat0))*facsh)
!    $       *abs(sin(2.*rlat))
c       flat_gw=amax1(flat_gw,0.10)
c       flat_gw=1.0
c    ****
c    ****  cos(2*rlat) distribution
c    ****
! 12/11/00: flat_gw as per kibo12
c       flat_gw = abs(2.*sin(rlat)*cos(rlat))
c       flat_gw=amax1(flat_gw,1.e-5)
        flat_gw = abs(2.*sin(rlat)*cos(rlat))
        if(rlat.gt.0.) flat_gw = flat_gw*facnh
        if(rlat.le.0.) flat_gw = flat_gw*facsh
        flat_gw=amax1(flat_gw,1.e-5)
	if(rlat.gt.0.) wspd = 20.
	if(rlat.le.0.) wspd = 20.
	if(rlat.gt.0.) whw = 30.
	if(rlat.le.0.) whw  =30.

c
c Set the phase speeds and wave numbers in the direction of the source wind.
c Set the source stress magnitude (positive only, note that the sign of the
c stress is the same as (c-u)).
c
	tauback=tauback*flat_gw
c
c       tauback=tauback*flat_gw+tauback1*wugws(jlat)
c       tauback=tauback*wugws(jlat)
c
c  **** equinox waves
      do l = 1, ngwv
	 do i = 1, plon
	    l1 = l*idir(i)
	    xlth = uth+(-180.+(i-1)*5.)/15.
! 12/11/00: flt as per kibo12:
c           flt = 1.0+0.5*sin((xlth-10.)*pi_g/12.)
	    flt = 1.
	    tauback1 = tauback*flt
c           tauback1 = tauback
! 12/11/00: taueq and tau as per kibo12:
            taueq = taueq*flt
	 tau(i,l1,kbot) = (exp(-((cphs(l1)-idir(i)*wspd)/whw)**2)
     $                       *tauback1
     $                     +exp(-((cphs(l1)-idir(i)*wspde)/whwe)**2)
     $                       *taueq)*abs(xv(i))
     $                      +tauback1*exp(-(cphs(l1)/whw)**2)
     $                       *abs(yv(i))*0.5
c    $                       *abs(yv(i))
	 tau(i,-l1,kbot) = (exp(-((cphs(l1)+idir(i)*wspd)/whw)**2)
     $                       *tauback1
     $                     +exp(-((cphs(l1)+idir(i)*wspde)/whwe)**2)
     $                       *taueq)*abs(xv(i))
     $                      +tauback1*exp(-(cphs(l1)/whw)**2)
     $                       *abs(yv(i))*0.5
c    $                       *abs(yv(i))
c
c    waltershied spectra
c
c        tau(i,l1,kbot) = (1./(cphs(l1)-idir(i)*20.)**2*15.
c    $                       *tauback1*abs(xv(i))
c    $                      +tauback1*1./cphs(l1)**2*15.
c    $                       *abs(yv(i))*0.5)
c    $                    *l*3.
c        tau(i,-l1,kbot) = (1./(cphs(l1)+idir(i)*20.)**2*15.
c    $                       *tauback1*abs(xv(i))
c    $                      +tauback1*1./cphs(l1)**2*15.
c    $                       *abs(yv(i))*0.5)
c    $                    *l*3.
c
c        tau(i,l1,kbot) = (1./(cphs(l1)-idir(i)*20.)**2*10.
c    $                       *tauback1*abs(xv(i))
c    $                      +tauback1*1./cphs(l1)**2*10.
c    $                       *abs(yv(i))*0.5)
c    $                    *l*3.
c        tau(i,-l1,kbot) = (1./(cphs(l1)+idir(i)*20.)**2*10.
c    $                       *tauback1*abs(xv(i))
c    $                      +tauback1*1./cphs(l1)**2*10.
c    $                       *abs(yv(i))*0.5)
c    $                    *l*3.
c
	 end do
      end do
c
c  **** solstice waves
c
c     do l = 1, ngwv
c        do i = 1, plon
c           l1 = l*idir(i)
c           xlth = uth+(-180.+(i-1)*5.)/15.
c           flt = 1.0+0.5*sin((xlth-10.)*pi_g/12.)
c           tauback1 = tauback*flt
c        tau(i,l1,kbot) = (exp(-((cphs(l1)-idir(i)*5.)/30.)**2)
c    $                       *tauback1
c    $                     +exp(-(cphs(l1)/10.)**2)
c    $                       *taueq)*abs(xv(i))
c    $                      +tauback1*exp(-(cphs(l1)/30.)**2)
c    $                       *abs(yv(i))*0.5
c        tau(i,-l1,kbot) = (exp(-((cphs(l1)+idir(i)*5.)/30.)**2)
c    $                       *tauback1
c    $                    +exp(-(cphs(l1)/10.)**2)
c    $                       *taueq)*abs(xv(i))
c    $                      +tauback1*exp(-(cphs(l1)/30.)**2)
c    $                       *abs(yv(i))*0.5
c        end do
c     end do
c
! 12/11/00: tautmp as per kibo12:
c     tautmp = 1.e-4
      tautmp = 1.e-6
c     if(rlat.gt.0.6981.and.rlat.lt.1.57) tautmp = 1.4e-3
c    $     *exp(-(rlat-1.134)**2/0.2684**2)*facoron
c     if(rlat.gt.-1.57.and.rlat.lt.-0.6981) tautmp = 1.4e-3
c    $     *exp(-(rlat+1.134)**2/0.2684**2)*facoros
c
      do i = 1, plon
	    xlth = uth+(-180.+(i-1)*5.)/15.
! 12/11/00: flt and tau as per kibo12:
c           flt = 1.0+0.5*sin((xlth-10.)*pi_g/12.)
            flt = 1.0
	 tau(i,0,kbotoro) = tautmp+1.0e-4*sghdat(i,jlat)/1240.
     $                      *flt
c        tau(i,0,kbotoro) = tautmp+1.0e-3*sghdat(i,jlat)/1240.
c    $                      *flt
c        tau(i,0,kbotoro) = tautmp+2.0e-3*sghdat(i,jlat)/1240.
c    $                      *flt
c        tau(i,0,kbotoro) = tautmp+4.0e-3*sghdat(i,jlat)/1240.
c    $                      *flt
      end do
c
c +
c Determine the min value of kldv and ksrc for limiting later loops
c and the pressure at the top interface of the low level stress divergence
c region.
c -
      ksrcmn = plev
      kldvmn = plev

      return
      end