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

!=======================================================================
!     set viscosity coefficient on bottom face of "u" cells
!     set diffusion coefficient on bottom face of "t" cells

!     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

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

      include "param.h"
      include "coord.h"
      include "mw.h"
      include "switch.h"
      include "vmixc.h"
      include "isopyc.h"
      include "tidal_kv.h"
      include "diag.h"
      include "grdvar.h"
      include "levind.h"
      include "timeavgs.h"

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

      if (js .gt. je) return

!-----------------------------------------------------------------------
!     limit the longitude and latitude indices
!-----------------------------------------------------------------------

      istrt = max(2,is)
      iend  = min(imt-1,ie)
      jstrt = max(2,js-1)
      jend  = je-1

!-----------------------------------------------------------------------
!     constant vertical mixing coefficients
!-----------------------------------------------------------------------

      do j=jstrt,jend
        jrow = j + joff
        do i=istrt,iend
          do k=1,km
            visc_cbu(i,k,j) = kappa_m

!           calculate N^2 = -g/rho drhodz on bottom of cell face
!           (where K33 and diff_cbt = kappa_h are defined). Note that
!           N2 is not guaranteed to be positive. If instability occurs,
!           convective adjustment will eliminate it.
!           drodzb is defined in isopyc.h

!           ZN2 is defined on T-cell bottom (zw pt)
            ZN2 = -gravrho0r*drodzb(i,k,j,0)

!           height above bottom
            if (kmt(i,jrow) .ne. 0.0) then
              hab = zw(k) - zw(kmt(i,jrow) - 1)
            else
              hab = 0.0
            endif
            if (Zn2 .ne. 0.0) then
              zkappa = 0.33*ogamma*edr(i,jrow)*exp(hab*zetar)/
     &                 (ZN2*(1-exp(-zetar*zw(kmt(i,jrow)))))
            else
              zkappa = 0.0
            endif
! Andreas: high mixing in Southern Ocean below 500 m as in observations
            if (tlat(i,j+joff) .lt. -40. .and.
     &           zw(kmt(i,j+joff)).gt.50000.) zkappa = max(zkappa,1.00)
!           limit diff_cbt
            diff_cbt(i,k,j) = max(kappa_h, min(100., zkappa + kappa_h))
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     accumulate time average diapycnal (without K33) diffusivity
!-----------------------------------------------------------------------

      if (timavgperts .and. .not. euler2) then
        do j=jstrt,jend
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              ta_diff_cbt(i,k,jrow) = ta_diff_cbt(i,k,jrow) +
     &                              diff_cbt(i,k,j)
            enddo
          enddo
        enddo
      endif

!-----------------------------------------------------------------------
!     Add K33 component to vertical diffusion coefficient
!-----------------------------------------------------------------------

      do j=jstrt,jend
        do i=istrt,iend
          do k=1,km
            diff_cbt(i,k,j) = diff_cbt(i,k,j) + K33(i,k,j)
          enddo
        enddo
      enddo

      return
      end