! source file: /Users/csomes/Research/Models/UVic_ESCM/2.9/sim/nonredfield_dom/clean_up_code/seaice_wind_bug/afs_high/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
!=======================================================================

      implicit none

      integer i, k, j, ip, kr, jq, js, je, istrt, is, iend, ie, jstrt
      integer jend, jrow, joff, ks

      real zn2, hab, zkappa

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.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

            if (tlat(i,j+joff) .lt. -40.) zkappa = max(zkappa,1.00)
cao below is the old code used in "my" standard version and in Schmittner et al. (2009). a.o., 27.12.10
!           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(10., zkappa + kappa_h))
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     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

!-----------------------------------------------------------------------
!     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

      return
      end