! source file: /Users/oschlies/UVIC/master/source/mom/adv_vel.F
      subroutine adv_vel (joff, js, je, is, ie)

!=======================================================================
!     calculate advection velocities for momentum and tracer equations

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW

!     output:
!       adv_vet = advection velocity on east face of "t" cell
!       adv_vnt = advection velocity on north face of "t" cell
!       adv_vbt = advection velocity on bottom face of "t" cell
!       adv_veu = advection velocity on east face of "u" cell
!       adv_vnu = advection velocity on north face of "u" cell
!       adv_vbu = advection velocity on bottom face of "u" cell

!     based on code by: R. C. Pacanowski
!=======================================================================

      include "param.h"
      include "coord.h"
      include "grdvar.h"
      include "levind.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     bail out if starting row exceeds ending row
!-----------------------------------------------------------------------

      if (js .gt. je) return

!-----------------------------------------------------------------------
!     limit the longitude indices
!-----------------------------------------------------------------------

      istrt = max(2,is)
      iend  = min(imt-1,ie)

!-----------------------------------------------------------------------
!     advection velocity on northern face of "T" cells. Note the
!     embedded cosine.
!     adv_vnt = WT_AVG_X(u(1,1,1,2,tau))
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,km
          do i=istrt,iend
            adv_vnt(i,k,j) = (u(i,k,j,2,tau)*dxu(i) +
     &                     u(i-1,k,j,2,tau)*dxu(i-1))*csu(jrow)*dxt2r(i)
          enddo
        enddo
        call setbcx (adv_vnt(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     advection velocity on the eastern face of "T" cells
!     adv_vnt = WT_AVG_Y(u(1,1,1,1,tau))
!-----------------------------------------------------------------------

      jstbe = max(js,jsmw)
      do j=jstbe,je
        jrow = j + joff
        do k=1,km
          do i=istrt-1,iend+1
            adv_vet(i,k,j) = (u(i,k,j,1,tau)*dyu(jrow) +
     &                     u(i,k,j-1,1,tau)*dyu(jrow-1))*dyt2r(jrow)
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     construct vertical velocity on the bottom face of "T" cells
!-----------------------------------------------------------------------

      do j=jstbe,je
        jrow = j + joff

!       set "adv_vbt" at surface to 0.0 (rigid-lid) or dh/dt (free surf)

        do i=istrt,iend
          adv_vbt(i,0,j)   = c0
        enddo

!       construct divergence of advection velocity * level thickness

        do k=1,km
          do i=istrt,iend
            adv_vbt(i,k,j) =
     &                   ((adv_vet(i,k,j) - adv_vet(i-1,k,j))*dxtr(i)
     &                   +(adv_vnt(i,k,j) - adv_vnt(i,k,j-1))*dytr(jrow)
     &                   )*cstr(jrow)*dzt(k)
          enddo
        enddo

!       integrate downward to define "adv_vbt" at the bottom of levels

        do k=1,km
          do i=istrt,iend
            adv_vbt(i,k,j) = adv_vbt(i,k,j) + adv_vbt(i,k-1,j)
          enddo
        enddo

        call setbcx (adv_vbt(1,0,j), imt, km+1)

      enddo

!-----------------------------------------------------------------------
!     construct advection velocity on the northern face of "u" cells by
!     averaging advection velocity on northern face of "t" cells
!     note: je-1 is used instead of jemw to account for possible non
!           integral number of MW`s in jmt
!     adv_vnu = LINEAR_INTRP_Y(WT_AVG_X(adv_vnt))
!-----------------------------------------------------------------------

      jsun = max(js,jsmw)-1
      do j=jsun,je-1
        jrow = j + joff
        dyr  = dytr(jrow+1)
        do k=1,km
          do i=istrt,iend
            adv_vnu(i,k,j) = ((adv_vnt(i,k,j)*duw(i)
     &                       + adv_vnt(i+1,k,j)*due(i)
     &                        )*dus(jrow+1) +
     &                        (adv_vnt(i,k,j+1)*duw(i)
     &                       + adv_vnt(i+1,k,j+1)*due(i)
     &                        )*dun(jrow))*dyr*dxur(i)
          enddo
        enddo
        call setbcx (adv_vnu(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     construct advection velocity on the eastern face of "u" cells by
!     averaging advection velocity on eastern face of "t" cells
!     note: take special care of zonal b.c. on this term.
!     adv_veu = LINEAR_INTRP_X(WT_AVG_Y(adv_vet))
!-----------------------------------------------------------------------

      jsube = max(js-1,jsmw)
      do j=jsube,je-1
        jrow = j + joff
        dyr  = dyur(jrow)
        do k=1,km
          do i=istrt-1,iend
            adv_veu(i,k,j) = ((adv_vet(i,k,j)*dus(jrow)
     &                       + adv_vet(i,k,j+1)*dun(jrow)
     &                        )*duw(i+1) +
     &                        (adv_vet(i+1,k,j)*dus(jrow)
     &                       + adv_vet(i+1,k,j+1)*dun(jrow)
     &                        )*due(i))*dyr*dxtr(i+1)
          enddo
        enddo
        call setbcx (adv_veu(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     construct advection velocity on the bottom face of "u" cells by
!     averaging advection velocity on bottom face of "t" cells
!-----------------------------------------------------------------------

      do j=jsube,je-1
        jrow = j + joff
        dyn  = dun(jrow)*cst(jrow+1)
        dys  = dus(jrow)*cst(jrow)
        dyr  = dyur(jrow)*csur(jrow)
        do k=0,km
          do i=istrt,iend
            asw = duw(i)*dys
            anw = duw(i)*dyn
            ase = due(i)*dys
            ane = due(i)*dyn
            adv_vbu(i,k,j) = dyr*dxur(i)*(
     &                    adv_vbt(i,k,j)*asw + adv_vbt(i+1,k,j)*ase
     &                  + adv_vbt(i,k,j+1)*anw + adv_vbt(i+1,k,j+1)*ane)
          enddo
        enddo

        call setbcx (adv_vbu(1,0,j), imt, km+1)

      enddo

!#define debug_adv_vel
      return
      end