! source file: /Users/oschlies/UVIC/master/testcase/updates/isopyc.F
      subroutine isopi (error, am, ah)

!=======================================================================

!               Initialization for isopycnal mixing scheme

!      -Disopycmix gives either the full or small angle isopycnal
!                  mixing tensor

!      -Disopycmix -Dgent_mcwilliams gives the isopycnal mixing tensor
!       plus the advective velocity parameterization of Gent_McWilliams

!     input:
!       error  = logical to signal problems
!       am     = background mixing coeff for momentum
!       ah     = horizontal mixing coeff for tracers
!       slmx   = max slope of isopycnals
!       ahisop = isopycnal tracer diffusivity(cm**2/sec)
!       athkdf = isopycnal thickness diffusivity (cm**2/sec)

!     output:
!       The above input can be reset via namelist

!     based on code by: R. C. Pacanowski, S. M. Griffies and
!                       G. Danabasoglu
!           Note: Extensive documentation is given in the manual
!=======================================================================

      logical error

      include "param.h"
      include "accel.h"
      include "coord.h"
      include "grdvar.h"
      include "iounit.h"
      include "isopyc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "vmixc.h"

      character(120) :: fname, new_file_name
      integer ib(10), ic(10)
      logical exists, inqvardef
      real tmp(imt,jmt)

      namelist /isopyc/ slmx, ahisop, athkdf, del_dm, s_dm

      write (stdout,'(/,20x,a,/,20x,a,/)')
     & 'I S O P Y C M I X    I N I T I A L I Z A T I O N'
     &,'Isopycnal mixing tensor + GM thickness diffusion'

!-----------------------------------------------------------------------
!     USER INPUT
!     initialize variables (all mixing units are cm**2/sec.)
!-----------------------------------------------------------------------

!     maximum isopycnal slope

      slmx  = 1.0/100.0

!     isopycnal diffusion coefficient

      ahisop = 1.e7

!     isopycnal thickness diffusion coefficient

      athkdf = 1.0e7

!     transition for scaling diffusion coefficient

      del_dm = 4.0/1000.0

!     half width scaling for diffusion coefficient

      s_dm = 1.0/1000.0

!-----------------------------------------------------------------------
!     provide for namelist over-ride of above settings + documentation
!-----------------------------------------------------------------------

      call getunit (io, 'control.in'
     &,               'formatted sequential rewind')
      read (io,isopyc,end=100)
100   continue
      write (stdout,isopyc)
      call relunit (io)

!     set reciprocal of maximum isopycnal slope and scaling half width

      slmxr = c1/slmx
      s_dmr = c1/s_dm

!-----------------------------------------------------------------------
!     check for problems
!-----------------------------------------------------------------------

        write (stdout,*)
     & '==> Note: small angle approximation to isopycnal tensor is'
     &,'             being used.                                     '
        write (stdout,*)
     & '==> Note: Re-scaling of Gerdes et al is being used to        '
     &,'             reduce mixing coeff within areas of steeply     '
     &,'             sloping isopycnals                              '

!-----------------------------------------------------------------------
!     print out tracer diffusion coefficients in isopycnal case
!-----------------------------------------------------------------------

      write(stdout,9102) ahisop, athkdf, ah, am

9102  format(/' ahisop = ',e12.6,' along isopyncal tracer mixing '
     &,  '(cm**2/sec) '/' athkdf = ',e12.6,' isopycnal thickness '
     &,'diffusion (cm**2/sec) '/' ah = ',e12.6,' cm**2/sec for '
     &,' am = ',e12.6,' cm**2/sec for horizontal mixing of momentum'/)
      ft = c1/(4.0*ahisop*dtts)
      delta_iso = dxt(1)*cst(jmt/2)*dzt(1)*ft
      i_delta = 1
      j_delta = 1
      k_delta = 1
      do jrow=2,jmt-1
        do i=2,imt-1
          do k=1,km
            delta1 = dxt(i)*cst(jrow)*dzt(k)*ft
            delta2 = dyt(jrow)*dzt(k)*ft
            if (delta_iso .ge. delta1 .or. delta_iso .ge. delta2) then
              i_delta = i
              j_delta = jrow
              k_delta = k
              delta_iso = min(delta1,delta2)
            endif
          enddo
        enddo
      enddo
      if (delta_iso .lt. p5) then
        s_minus = (c1 - sqrt(c1 - 4.0*delta_iso**2))/(c2*delta_iso)
        s_plus  = c1/s_minus
      else
        s_minus = c0
        s_plus  = c0
      endif

      write(stdout,'(a)')
     &'------------------------------------------------'
      write(stdout,'(a,e14.7)')
     &'The isopycnal diffusion grid factor delta_iso =',delta_iso
      write(stdout,'(a)')
     &'was determined at the grid point'
      write(stdout,'(a,i4,a,e14.7)')
     &'dxt(',i_delta,') = ',dxt(i_delta)
      write(stdout,'(a,i4,a,e14.7)')
     &'dyt(',j_delta,') = ',dyt(j_delta)
      write(stdout,'(a,i4,a,e14.7)')
     &'dzt(',k_delta,') = ',dzt(k_delta)
      write(stdout,'(a,i4,a,e14.7)')
     &'cst(',j_delta,') = ',cst(j_delta)
      write(stdout,'(a,e14.7)')
     &'dtts             = ',dtts
      write(stdout,'(a,e14.7)')
     &'ahisop           = ',ahisop

      write(stdout,'(/a,e14.7/)')
     & 'Maximum allowable isopycnal slope is specified as slmx = ',slmx

      write(stdout,'(a)')
     &'------------------------------------------------'

!-----------------------------------------------------------------------
!     store the square root of the tracer timestep acceleration values
!     into variable "dtxsqr" for use in isopycnal mixing
!-----------------------------------------------------------------------

      do k=1,km
       dtxsqr(k) = sqrt(dtxcel(k))
      enddo

      write (stdout,'(a)') ' Acceleration with depth "dtxsqr(k)"='
      write (stdout,'(5(1x,e12.6))') (dtxsqr(k),k=1,km)

!-----------------------------------------------------------------------
!     read ocean diffusion factor
!-----------------------------------------------------------------------

      fisop(:,:,:) = 1.
      fname = new_file_name ("odf.nc")
      inquire (file=trim(fname), exist=exists)
      if (exists) then
        ib(:) = 1
        ic(:) = imt
        ic(2) = jmt
        ic(3) = km
        call openfile (fname, iou)
        exists = inqvardef('odf', iou)
        if (exists) then
          write(*,*) "reading ocean diffusion factor from: ",fname
          call getvara ('odf', iou, imt*jmt*km, ib, ic, fisop, c1, c0)
        endif
        call closefile (iou)
      endif

!-----------------------------------------------------------------------
!     initialize arrays
!-----------------------------------------------------------------------
      do j=1,jmw
        do k=1,km
          do i=1,imt
            alphai(i,k,j) = c0
            betai(i,k,j)  = c0
          enddo
        enddo
      enddo

      do j=jsmw,jemw
        do n=1,2
          do k=1,km
            do i=1,imt
              ddxt(i,k,j,n) = c0
            enddo
          enddo
        enddo
      enddo

      do j=1,jemw
        do n=1,2
          do k=1,km
            do i=1,imt
              ddyt(i,k,j,n) = c0
            enddo
          enddo
        enddo
      enddo
      do j=1,jmw
        do n=1,2
          do k=0,km
            do i=1,imt
              ddzt(i,k,j,n) = c0
            enddo
          enddo
        enddo
      enddo

      do j=jsmw,jemw
        do k=1,km
          do i=1,imt
            do kr=0,1
              do ip=0,1
                Ai_ez(i,k,j,ip,kr)  = c0
              enddo
            enddo
            do kr=0,1
              do ip=0,1
                Ai_bx(i,k,j,ip,kr) = c0
              enddo
            enddo
            do kr=0,1
              do jq=0,1
                Ai_by(i,k,j,jq,kr) = c0
              enddo
            enddo
            K33(i,k,j) = c0
            K11(i,k,j) = c0
            adv_vetiso(i,k,j) = c0
          enddo
        enddo
      enddo

      do j=1,jemw
        do k=1,km
          do i=1,imt
            do kr=0,1
              do jq=0,1
                Ai_nz(i,k,j,jq,kr)  = c0
              enddo
            enddo
            K22(i,k,j)        = c0
            adv_vntiso(i,k,j) = c0
          enddo
        enddo
      enddo

      do j=jsmw,jemw
        do k=0,km
          do i=1,imt
            adv_vbtiso(i,k,j) = c0
            adv_fbiso(i,k,j)  = c0
          enddo
        enddo
      enddo

      return
      end

      subroutine elements (joff, js, je, is, ie)

!=======================================================================
!     Estimate alpha, beta, and normal gradients on faces of T cells
!=======================================================================

      include "param.h"
      include "accel.h"
      include "grdvar.h"
      include "isopyc.h"
      include "levind.h"
      include "mw.h"
      include "state.h"
      include "dens.h"

!-----------------------------------------------------------------------
!     alpha and beta at centers of T cells
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,km
          do i=is,ie
            tprime = t(i,k,j,1,taum1)-to(k)
            sprime = t(i,k,j,2,taum1)-so(k)
            alphai(i,k,j) = drodt (tprime, sprime, k)
            betai(i,k,j)  = drods (tprime, sprime, k)
          enddo
        enddo
        call setbcx (alphai(1,1,j), imt, km)
        call setbcx (betai(1,1,j),  imt, km)
      enddo

!-----------------------------------------------------------------------
!     gradients at bottom face of T cells
!-----------------------------------------------------------------------

      do j=js,je
        do n=1,2
          do k=1,km
            kp1 = min(k+1,km)
            do i=is,ie
              ddzt(i,k,j,n) = tmask(i,kp1,j)*dzwr(k)*
     &                        (t(i,k,j,n,taum1) - t(i,kp1,j,n,taum1))
            enddo
          enddo
          do i=is,ie
            ddzt(i,0,j,n) = c0
          enddo
          call setbcx (ddzt(1,0,j,n), imt, km+1)
        enddo
      enddo

!-----------------------------------------------------------------------
!     gradients at eastern face of T cells
!-----------------------------------------------------------------------

      do j=max(js-1,2), je-1
        jrow = j + joff
        do n=1,2
          do k=1,km
            do i=is,ie
              ddxt(i,k,j,n) = tmask(i,k,j)*tmask(i+1,k,j)*cstr(jrow)*
     &                   dxur(i)*(t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
          call setbcx (ddxt(1,1,j,n), imt, km)
        enddo
      enddo

!-----------------------------------------------------------------------
!     gradients at northern face of T cells
!-----------------------------------------------------------------------

      do j=max(js-1,1), je-1
        jrow = j + joff
        do n=1,2
          do k=1,km
            kp1 = min(k+1,km)
            do i=is,ie
              ddyt(i,k,j,n) = tmask(i,k,j)*tmask(i,k,j+1)*dyur(jrow)*
     &                 (t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
          call setbcx (ddyt(1,1,j,n), imt, km)
        enddo
      enddo

      return
      end

      subroutine isopyc (joff, js, je, is, ie)

!=======================================================================

!     Compute isopycnal diffusion coefficients Ai_ez, Ai_nz, Ai_bx,
!     Ai_by, and the Gent/McWilliams advection velocities.

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

!     output:
!       Ai_ez = diffusion coefficient centered on east face of T cells
!       Ai_nz = diffusion coefficient centered on north face of T cells
!       Ai_bx = diffusion coefficient centered on bottom face of T cells
!       Ai_by = diffusion coefficient centered on bottom face of T cells

!       adv_vetiso = isopycnal advective vel on east face of "T" cell
!       adv_vntiso = isopycnal advective vel on north face of "T" cell
!               (Note: this includes the cosine factor as in "adv_vnt")
!       adv_vbtiso = isopycnal advective vel on bottom face of "T" cell

!=======================================================================

      include "param.h"

!-----------------------------------------------------------------------
!     set local constants
!-----------------------------------------------------------------------

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

!-----------------------------------------------------------------------
!     estimate alpha, beta, and gradients on sides of T cells
!-----------------------------------------------------------------------

      call elements (joff, js, je, is, ie)

!-----------------------------------------------------------------------
!     compute Ai_ez centered on eastern face of T cells
!     1 <= MW < last: calculate for rows jsmw..jemw
!     last MW       : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------

      call ai_east (joff, max(js-1,2), je-1, is, ie)

!-----------------------------------------------------------------------
!     compute Ai_nz centered on the northern face of T cells
!     first MW     : calculate for rows 1..jemw
!     1 < MW < last: calculate for rows jsmw..jemw
!     last MW      : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------

      call ai_north (joff, max(js-1,1), je-1, is, ie)

!-----------------------------------------------------------------------
!     evaluate Ai_bx & Ai_by centered on bottom face of T cells
!     1 <= MW < last: calculate for rows jsmw..jemw
!     last MW       : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------

      call ai_bottom (joff, max(js-1,2), je-1, is, ie)

!-----------------------------------------------------------------------
!     compute isopycnal advective velocities for tracers
!     first MW     : calculate (viso)           for rows 1..jemw
!                              (uiso,wiso)      for rows 2..jemw
!     1 < MW < last: calculate (uiso,viso,wiso) for rows jsmw..jemw
!     last MW      : calculate (uiso,viso,wiso) for rows jsmw..<=jemw
!-----------------------------------------------------------------------

      call isopyc_adv (joff, max(js-1,1), je-1, is, ie)

      return
      end

      subroutine ai_east (joff, js, je, is, ie)

!=======================================================================
!     compute "Ai_ez" & "K11" at the center of eastern face of T cells.
!=======================================================================

      include "param.h"
      include "accel.h"
      include "coord.h"
      include "grdvar.h"
      include "hmixc.h"
      include "isopyc.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     compute "Ai_ez" on east face of T cell. Note re-scaling factor
!     which reduces mixing coefficient "Ai" where slope "sxe"
!     exceeds the critical slope "sc" for the small slope approx. For
!     the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,km
          sc = c1/(slmxr*dtxsqr(k))
          dzt4r = p5*dzt2r(k)
          do i=2,imtm1
            Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k))
            sumz = c0
            do kr=0,1
              do ip=0,1
                sxe = abs(drodxe(i,k,j,ip)/(drodze(i,k,j,ip,kr)+epsln))
                if (sxe .gt. sc) then
                  Ai_ez(i,k,j,ip,kr)  =  Ai0*tmask(i,k,j)*tmask(i+1,k,j)
     &                                   *(sc/(sxe + epsln))**2
                else
                  Ai_ez(i,k,j,ip,kr)  =  Ai0*tmask(i,k,j)*tmask(i+1,k,j)
                endif
                sumz = sumz + dzw(k-1+kr)*Ai_ez(i,k,j,ip,kr)
              enddo
            enddo
            K11(i,k,j) = dzt4r*sumz
          enddo
        enddo
        call setbcx (Ai_ez(1,1,j,0,0), imt, km)
        call setbcx (Ai_ez(1,1,j,1,0), imt, km)
        call setbcx (Ai_ez(1,1,j,0,1), imt, km)
        call setbcx (Ai_ez(1,1,j,1,1), imt, km)
        call setbcx (K11(1,1,j), imt, km)
      enddo

      return
      end

      subroutine ai_north (joff, js, je, is, ie)

!=======================================================================
!     compute "Ai_nz" & "K22" at the center of northern face of T cells.
!=======================================================================

      include "param.h"
      include "accel.h"
      include "coord.h"
      include "grdvar.h"
      include "hmixc.h"
      include "isopyc.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     compute "Ai_nz" on north face of T cell. Note re-scaling factor
!     which reduces mixing coefficient "Ai" where slope "syn"
!     exceeds the critical slope "sc" for the small slope approx. For
!     the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,km
          sc = c1/(slmxr*dtxsqr(k))
          dzt4r = p5*dzt2r(k)
          do i=2,imtm1
            Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k))
            sumz = c0
            do kr=0,1
              do jq=0,1
                syn =abs(drodyn(i,k,j,jq)/(drodzn(i,k,j,jq,kr)+epsln))
                if (syn .gt. sc) then
                  Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
     &                                   *(sc/(syn + epsln))**2
                else
                  Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
                endif
                sumz = sumz + dzw(k-1+kr)*Ai_nz(i,k,j,jq,kr)
              enddo
            enddo
            K22(i,k,j) = dzt4r*sumz
          enddo
        enddo
        call setbcx (Ai_nz(1,1,j,0,0), imt, km)
        call setbcx (Ai_nz(1,1,j,1,0), imt, km)
        call setbcx (Ai_nz(1,1,j,0,1), imt, km)
        call setbcx (Ai_nz(1,1,j,1,1), imt, km)
        call setbcx (K22(1,1,j), imt, km)
      enddo

      return
      end

      subroutine ai_bottom (joff, js, je, is, ie)

!=======================================================================
!     compute Ai_bx, Ai_by, and K33 at the center of the bottom face of
!     T cells.
!=======================================================================

      include "param.h"
      include "accel.h"
      include "grdvar.h"
      include "hmixc.h"
      include "isopyc.h"
      include "mw.h"

!-----------------------------------------------------------------------
!     compute "Ai_bx", "Ai_by", & K33 on bottom face of T cell. Note
!     re-scaling factor to reduce mixing coefficient "Ai" where slopes
!     "sxb" and "syb" exceeds the critical slope "sc" for the small slope approx. For
!     the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,kmm1
          sc = c1/(slmxr*dtxsqr(k))
          kp1   = min(k+1,km)
          do i=2,imtm1
            Ai0 = ahisop*p5*(fisop(i,jrow,k+1) + fisop(i,jrow,k))

!           eastward slopes at the base of T cells

            sumx = c0
            do ip=0,1
              do kr=0,1
                sxb = abs(drodxb(i,k,j,ip,kr)/(drodzb(i,k,j,kr)+epsln))
                if (sxb .gt. sc) then
                  Ai_bx(i,k,j,ip,kr)  =  Ai0*tmask(i,k+1,j)
     &                                  *(sc/(sxb + epsln))**2
                else
                  Ai_bx(i,k,j,ip,kr)  =  Ai0*tmask(i,k+1,j)
                endif
                sumx = sumx + dxu(i-1+ip)*Ai_bx(i,k,j,ip,kr)*sxb**2
              enddo
            enddo

!           northward slopes at the base of T cells

            sumy = c0
            do jq=0,1
              facty = csu(jrow-1+jq)*dyu(jrow-1+jq)
              do kr=0,1
                syb = abs(drodyb(i,k,j,jq,kr)/(drodzb(i,k,j,kr)+epsln))
                if (syb .gt. sc) then
                  Ai_by(i,k,j,jq,kr)  =  Ai0*tmask(i,k+1,j)
     &                                  *(sc/(syb + epsln))**2
                else
                  Ai_by(i,k,j,jq,kr)  =  Ai0*tmask(i,k+1,j)
                endif
                sumy = sumy + facty*Ai_by(i,k,j,jq,kr)*syb**2
              enddo
            enddo

            K33(i,k,j) = dxt4r(i)*sumx + dyt4r(jrow)*cstr(jrow)*sumy
          enddo
        enddo
        call setbcx (Ai_bx(1,1,j,1,0), imt, km)
        call setbcx (Ai_bx(1,1,j,0,0), imt, km)
        call setbcx (Ai_bx(1,1,j,1,1), imt, km)
        call setbcx (Ai_bx(1,1,j,0,1), imt, km)
        call setbcx (Ai_by(1,1,j,1,0), imt, km)
        call setbcx (Ai_by(1,1,j,0,0), imt, km)
        call setbcx (Ai_by(1,1,j,1,1), imt, km)
        call setbcx (Ai_by(1,1,j,0,1), imt, km)
        call setbcx (K33(1,1,j), imt, km)
      enddo

      return
      end

      subroutine isoflux (joff, js, je, is, ie, n)

!=======================================================================
!     isopycnal diffusive tracer fluxes are computed.
!=======================================================================

      include "param.h"
      include "coord.h"
      include "grdvar.h"
      include "hmixc.h"
      include "isopyc.h"
      include "levind.h"
      include "mw.h"
      include "vmixc.h"

!-----------------------------------------------------------------------
!     construct total isopycnal tracer flux at east face of "T" cells
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,km
          dzt4r = p5*dzt2r(k)
          do i=2,imtm1
            sumz = c0
            do kr=0,1
              km1kr = max(k-1+kr,1)
              kpkr = min(k+kr,km)
              do ip=0,1
                sumz = sumz - Ai_ez(i,k,j,ip,kr)
     &                *(t(i+ip,km1kr,j,n,taum1)-t(i+ip,kpkr,j,n,taum1))
     &                *drodxe(i,k,j,ip)
     &               /(drodze(i,k,j,ip,kr) + epsln)
              enddo
            enddo
            flux_x = dzt4r*sumz
            diff_fe(i,k,j) = diff_fe(i,k,j) + K11(i,k,j)*cstdxur(i,j)*
     &                       (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
     &                     + flux_x
          enddo
        enddo
        call setbcx (diff_fe(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     construct total isopycnal tracer flux at north face of "T" cells
!-----------------------------------------------------------------------

      do j=js-1,je
        jrow = j + joff
        do k=1,km
          csu_dzt4r = csu(jrow)*p5*dzt2r(k)
          do i=2,imtm1
            sumz = c0
            do kr=0,1
              km1kr = max(k-1+kr,1)
              kpkr = min(k+kr,km)
              do jq=0,1
                sumz = sumz - Ai_nz(i,k,j,jq,kr)
     &                 *(t(i,km1kr,j+jq,n,taum1)-t(i,kpkr,j+jq,n,taum1))
     &                 *drodyn(i,k,j,jq)
     &               /(drodzn(i,k,j,jq,kr) + epsln)
              enddo
            enddo
            flux_y = csu_dzt4r*sumz
            diff_fn(i,k,j) = diff_fn(i,k,j)+K22(i,k,j)*csu_dyur(jrow)*
     &                       (t(i,k,j+1,n,taum1)-t(i,k,j,n,taum1))
     &                    + flux_y
          enddo
        enddo
        call setbcx (diff_fn(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     compute the vertical tracer flux "diff_fbiso" containing the K31
!     and K32 components which are to be solved explicitly. The K33
!     component will be treated implicitly. Note that there are some
!     cancellations of dxu(i-1+ip) and dyu(jrow-1+jq)
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,kmm1
          do i=2,imtm1
            sumx = c0
            do ip=0,1
              do kr=0,1
                sumx = sumx
     &            - Ai_bx(i,k,j,ip,kr)*cstr(jrow)*
     &             (t(i+ip,k+kr,j,n,taum1) - t(i-1+ip,k+kr,j,n,taum1))
     &             *drodxb(i,k,j,ip,kr)
     &              /(drodzb(i,k,j,kr) + epsln)
              enddo
            enddo
            sumy = c0
            do jq=0,1
              do kr=0,1
                sumy = sumy
     &           - Ai_by(i,k,j,jq,kr)*csu(jrow-1+jq)*
     &            (t(i,k+kr,j+jq,n,taum1)-t(i,k+kr,j-1+jq,n,taum1))
     &               *drodyb(i,k,j,jq,kr)
     &              /(drodzb(i,k,j,kr) + epsln)
              enddo
            enddo
            diff_fbiso(i,k,j) = dxt4r(i)*sumx
     &                          + dyt4r(jrow)*cstr(jrow)*sumy
          enddo
        enddo
        do i=2,imtm1
          diff_fbiso(i,0,j)  = c0
          diff_fbiso(i,km,j) = c0
        enddo
        call setbcx (diff_fbiso(1,0,j), imt, km+1)
      enddo

!-----------------------------------------------------------------------
!     compute advective tracer flux at the center of the bottom face of
!     the "T" cells
!-----------------------------------------------------------------------

      do j=js,je
        do k=1,kmm1
          do i=2,imt-1
            adv_fbiso(i,k,j) = adv_vbtiso(i,k,j)*
     &                         (t(i,k,j,n,taum1) + t(i,k+1,j,n,taum1))
          enddo
        enddo
      enddo

!     now consider the top and bottom boundaries

      do j=js,je
        do i=2,imt-1
          adv_fbiso(i,0,j)  = c0
          adv_fbiso(i,km,j) = c0
        enddo
      enddo

      return
      end

      subroutine isopyc_adv (joff, js, je, is, ie)

!=======================================================================
!     compute isopycnal transport velocities.

!     based on code by:  G. Danabasoglu
!=======================================================================

      include "param.h"
      include "accel.h"
      include "coord.h"
      include "grdvar.h"
      include "isopyc.h"
      include "levind.h"
      include "mw.h"
      include "switch.h"
      include "timeavgs.h"

      dimension top_bc(km), bot_bc(km)

      do k=1,km
        top_bc(k) = c1
        bot_bc(k) = c1
      enddo
      top_bc(1)  = c0
      bot_bc(km) = c0

!-----------------------------------------------------------------------
!     compute the meridional component of the isopycnal mixing velocity
!     at the center of the northern face of the "T" cells.
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff
        do k=1,km

          sc = c1/(slmxr*dtxsqr(k))
          km1 = max(k-1,1)
          kp1 = min(k+1,km)
          do i=1,imt
            Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k))

            at =     (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,km1,j)
     &              + alphai(i,km1,j+1))
            bt =     (betai(i,k,j) + betai(i,k,j+1) + betai(i,km1,j)
     &              + betai(i,km1,j+1))
            stn = -(at*(ddyt(i,k,j,1) + ddyt(i,km1,j,1))
     &               + bt*(ddyt(i,k,j,2) + ddyt(i,km1,j,2))) /
     &                (at*(ddzt(i,km1,j,1) + ddzt(i,km1,j+1,1))
     &               + bt*(ddzt(i,km1,j,2) + ddzt(i,km1,j+1,2))+epsln)

            ab =     (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,kp1,j)
     &              + alphai(i,kp1,j+1))
            bb =     (betai(i,k,j) + betai(i,k,j+1) + betai(i,kp1,j)
     &              + betai(i,kp1,j+1))
            sbn = -(ab*(ddyt(i,k,j,1) + ddyt(i,kp1,j,1))
     &               + bb*(ddyt(i,k,j,2) + ddyt(i,kp1,j,2))) /
     &                (ab*(ddzt(i,k,j,1) + ddzt(i,k,j+1,1))
     &               + bb*(ddzt(i,k,j,2) + ddzt(i,k,j+1,2))+epsln)
            absstn = abs(stn)
            abssbn = abs(sbn)
            if (absstn .gt. sc) then
              ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
     &              *(sc/(absstn + epsln))**2
            else
              ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
            endif
            if (abssbn .gt. sc) then
              ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
     &              *(sc/(abssbn + epsln))**2
            else
              ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
            endif
            adv_vntiso(i,k,j) = -(ath_t*stn*top_bc(k) -
     &                            ath_b*sbn*bot_bc(k))*dztr(k)*csu(jrow)
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     compute the zonal component of the isopycnal mixing velocity
!     at the center of the eastern face of the "T" grid cell.
!-----------------------------------------------------------------------

      jstrt = max(js,jsmw)
      do j=jstrt,je
        jrow = j + joff
        do k=1,km
          sc = c1/(slmxr*dtxsqr(k))
          km1 = max(k-1,1)
          kp1 = min(k+1,km)
          do i=1,imtm1
            Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k))

            at =     (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,km1,j)
     &              + alphai(i+1,km1,j))
            bt =     (betai(i,k,j) + betai(i+1,k,j) + betai(i,km1,j)
     &              + betai(i+1,km1,j))
            ste =-(at*(ddxt(i,k,j,1) + ddxt(i,km1,j,1))
     &           + bt*(ddxt(i,k,j,2) + ddxt(i,km1,j,2)))
     &          / (at*(ddzt(i,km1,j,1) + ddzt(i+1,km1,j,1))
     &           + bt*(ddzt(i,km1,j,2) + ddzt(i+1,km1,j,2))+epsln)

            ab =     (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,kp1,j)
     &              + alphai(i+1,kp1,j))
            bb =     (betai(i,k,j) + betai(i+1,k,j) + betai(i,kp1,j)
     &              + betai(i+1,kp1,j))
            sbe =-(ab*(ddxt(i,k,j,1) + ddxt(i,kp1,j,1))
     &           + bb*(ddxt(i,k,j,2) + ddxt(i,kp1,j,2))) /
     &            (ab*(ddzt(i,k,j,1) + ddzt(i+1,k,j,1))
     &           + bb*(ddzt(i,k,j,2) + ddzt(i+1,k,j,2))+epsln)

            absste = abs(ste)
            abssbe = abs(sbe)
            if (absste .gt. sc) then
              ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
     &              *(sc/(absste + epsln))**2
            else
              ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
            endif
            if (abssbe .gt. sc) then
              ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
     &              *(sc/(abssbe + epsln))**2
            else
              ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
            endif
            adv_vetiso(i,k,j) = -(ath_t*ste*top_bc(k) -
     &                            ath_b*sbe*bot_bc(k))*dztr(k)
          enddo
        enddo
      enddo

!     set the boundary conditions

      do j=jstrt,je
        call setbcx (adv_vetiso(1,1,j), imt, km)
      enddo

!-----------------------------------------------------------------------
!     compute the vertical component of the isopycnal mixing velocity
!     at the center of the bottom face of the "T" cells, using the
!     continuity equation for the isopycnal mixing velocities
!-----------------------------------------------------------------------

      do j=jstrt,je
        do i=1,imt
          adv_vbtiso(i,0,j) = c0
        enddo
      enddo

      do j=jstrt,je
        jrow = j + joff
        do k=1,kmm1
          do i=2,imt
            adv_vbtiso(i,k,j) = dzt(k)*cstr(jrow)*(
     &      (adv_vetiso(i,k,j) - adv_vetiso(i-1,k,j))*dxtr(i) +
     &      (adv_vntiso(i,k,j) - adv_vntiso(i,k,j-1))*dytr(jrow))
          enddo
        enddo
      enddo

      do j=jstrt,je
        do k=1,kmm1
          do i=2,imt
            adv_vbtiso(i,k,j) = adv_vbtiso(i,k,j) + adv_vbtiso(i,k-1,j)
          enddo
        enddo
      enddo

      do j=jstrt,je
        jrow = j + joff
        do i=2,imt
          adv_vbtiso(i,kmt(i,jrow),j) = c0
        enddo
      enddo

!     set the boundary conditions

      do j=jstrt,je
        call setbcx (adv_vbtiso(1,0,j), imt, km+1)
      enddo

!-----------------------------------------------------------------------
!     accumulate time average gm velocities
!-----------------------------------------------------------------------

      if (timavgperts .and. .not. euler2) then
        do j=jstrt,je
          jrow = j + joff
          do k=1,kmm1
            do i=2,imt
              ta_vetiso(i,k,jrow) = ta_vetiso(i,k,jrow) +
     &                              adv_vetiso(i,k,j)
              ta_vntiso(i,k,jrow) = ta_vntiso(i,k,jrow) +
     &                              adv_vntiso(i,k,j)
              ta_vbtiso(i,k,jrow) = ta_vbtiso(i,k,jrow) +
     &                              adv_vbtiso(i,k,j)
            enddo
          enddo
        enddo
      endif
!#  define debug_adv_vel_iso
      return
      end