#include "dims.h"
      subroutine mgworo (u, v, t, sgh, pm, pi, dpm, zm, nm, pblh,rlat,
     $     kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv)
C-----------------------------------------------------------------------
C 
C Orographic source for multiple gravity wave drag parameterization.
c 
c The stress is returned for a single wave with c=0, over orography.
c For points where the orographic variance is small (including ocean),
c the returned stress is zero. 
C 
C-----------------------------------------------------------------------
      implicit none
C-----------------------------------------------------------------------
#include "pmgrid.h"
C-----------------------------------------------------------------------
#include "mgw.h"
C-----------------------------------------------------------------------
C 
C Input variables
C 
      real 
     $     u(plond,plev),               ! midpoint zonal wind
     $     v(plond,plev),               ! midpoint meridional wind
     $     t(plond,plev),               ! midpoint temperatures
     $     sgh(plond),                  ! standard deviation of orography
     $     pm(plond,plev),              ! midpoint pressures
     $     pi(plond,0:plev),            ! interface pressures
     $     dpm(plond,plev),             ! midpoint delta p (pi(k)-pi(k-1))
     $     zm(plond,plev),              ! midpoint heights
     $     nm(plond,plev),              ! midpoint Brunt-Vaisalla frequency
     $     pblh(plond),                 ! planetary boundary layer height
     $     disp(plond),                 ! twice standard deviation
     $     tautmp
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)
     $     rlat
C 
C Local workspace
C 
      integer
     $     i,k                          ! loop indexes

      real
     $     sghmax,                      ! max orographic sdv to use
     $     tauoro(plond),               ! c=0 stress from orography
     $     zldv(plond),                 ! top of the low level stress divergence region
     $     nsrc(plond),                 ! b-f frequency averaged over source region
     $     psrc(plond),                 ! interface pressure at top of source region
     $     rsrc(plond),                 ! density averaged over source region
     $     usrc(plond),                 ! u wind averaged over source region
     $     vsrc(plond)                  ! v wind averaged over source region
C 
C---------------------------------------------------------------------------
c +
c Average the basic state variables for the wave source over the depth of
c the orographic standard deviation. Here we assume that the apropiate
c values of wind, stability, etc. for determining the wave source are 
c averages over the depth of the atmosphere pentrated by the typical mountain.
c Reduces to the bottom midpoint values when sgh=0, such as over ocean.
c 
c Also determine the depth of the low level stress divergence region, as
c the max of the boundary layer depth and the source region depth. This
c can be done here if the stress magnitude does not determine the depth,
c otherwise it must be done below.
c -
      disp = 2.*sgh
      k = plev
      do i = 1, plon
         ksrc(i) = k-1
         kldv(i) = k-1
         psrc(i) = pi(i,k-1)
         rsrc(i) = pm(i,k)/(r*t(i,k)) * dpm(i,k)
         usrc(i) = u(i,k) * dpm(i,k)
         vsrc(i) = v(i,k) * dpm(i,k)
         nsrc(i) = nm(i,k)* dpm(i,k)
      end do
      do k = plev-1, plev/2, -1
         do i = 1, plon
C            if (0.5*sgh(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
            if (disp(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
               ksrc(i) = k-1
               kldv(i) = k-1
               psrc(i) = pi(i,k-1)
               rsrc(i) = rsrc(i)
     $              + pm(i,k) / (r*t(i,k))* dpm(i,k)
               usrc(i) = usrc(i) + u(i,k) * dpm(i,k)
               vsrc(i) = vsrc(i) + v(i,k) * dpm(i,k)
               nsrc(i) = nsrc(i) + nm(i,k)* dpm(i,k)
            elseif (pblh(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
               kldv(i) = k-1
            end if
         end do
      end do
      do i = 1, plon
         rsrc(i) = rsrc(i) / (pi(i,plev) - psrc(i))
         usrc(i) = usrc(i) / (pi(i,plev) - psrc(i))
         vsrc(i) = vsrc(i) / (pi(i,plev) - psrc(i))
         nsrc(i) = nsrc(i) / (pi(i,plev) - psrc(i))

         ubi(i,plev) = sqrt (usrc(i)**2 + vsrc(i)**2)
         xv(i) = usrc(i) / ubi(i,plev)
         yv(i) = vsrc(i) / ubi(i,plev)
      end do
C 
C Project the local wind at midpoints onto the source wind.
C 
      do k = 1, plev
         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, plev-1
         do i = 1, plon
            ubi(i,k) = 0.5 * (ubm(i,k) + ubm(i,k+1))
         end do
      end do
c +
c Determine the orographic c=0 source term following McFarlane (1987).
c Set the source top interface index to plev, if the orographic term is zero.
c -
      tautmp = 1.e-4
c      write(*,*)'tauoro'
c      write(*,*)(tauoro(i),i=1,plon)
c      write(*,*)'ubi'
c      write(*,*)(ubi(i,plev),i=1,plon)
c      write(*,*)'nsrc'
c      write(*,*)(nsrc(i),i=1,plon)
c      write(*,*)'disp'
c      write(*,*)(disp(i),i=1,plon)
c      write(*,*)'rsrc'
c      write(*,*)(rsrc(i),i=1,plon)
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 
      if(rlat.gt.0.6981.and.rlat.lt.1.57) tautmp = 1.4e-3
     $     *exp(-(rlat-1.134)**2/0.2684**2)
      if(rlat.gt.-1.396.and.rlat.lt.-1.047) tautmp = 
     $     amax1(5.e-4*(sin((rlat+1.047)/
     $     0.1745*1.5708))**2,1.e-4)
      tauoro = tautmp
      do i = 1, plon
         tau(i,0,kbotoro) = tauoro(i)
      end do
c +
c Find the top interface of the low level stress divergence region according
c to the maximum depth of three criterion.
c 1. source region depth
c 2. planetary boundary layer depth
c 3. 10 * (u_*) / N where u_* is defined from the gravity wave stresss
c = sqrt(tau/rho) using source region values
c -
c$$$if (kbot .lt. plev) then
c$$$do i = 1, plon
c$$$kldv(i) = kbot
c$$$end do
c$$$else
c$$$do i = 1, plon
c$$$zldv(i) = max (pblh(i), sgh(i)
c$$$zldv(i) = max (zdlv(i),
c$$$$           zldvcon * sqrt(tau(i,0,k)/rsrc(i)) / nsrc(i))
c$$$end do
c$$$kldv(i) = plev-1
c$$$do k = plev-1, plev/2, -1
c$$$do i = 1, plon
c$$$if (zldv(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
c$$$kldv(i)  = k-1
c$$$end do
c$$$end do
c$$$end if
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
      do i = 1, plon
         ksrcmn = min(ksrcmn, ksrc(i))
         kldvmn = min(kldvmn, kldv(i))
         if (kldv(i) .ne. plev) then
            rdpldv(i) = 1. / (pi(i,kldv(i)) - pi(i,plev))
         end if
      end do
      if (fracldv .le. 0.) kldvmn = plev
      return
      end