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"
#if defined held_larichev
# include "hmixc.h"
#endif
#include "mw.h"
#include "switch.h"
#include "vmixc.h"
#if defined isopycmix || defined redi_diffusion
# include "isopyc.h"
# if defined uvic_tidal_kv
#  include "tidal_kv.h"
#  include "diag.h"
#  include "grdvar.h"
#  include "levind.h"
# endif
#endif
#if defined uvic_kv_out
# include "timeavgs.h"
#endif
#if defined constvmix && defined implicitvmix
      dimension temp(imt,km,jsmw:jmw)
#endif

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

#if defined constvmix

!-----------------------------------------------------------------------
!     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
# if defined uvic_tidal_kv

!           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 
!            if (tlat(i,j+joff) .lt. -40.) zkappa = max(zkappa,1.00)
!           limit diff_cbt
#  if defined bryan_lewis_vertical
            diff_cbt(i,k,j) = max(Ahv(k), min(100., zkappa + Ahv(k)))
#  else
            diff_cbt(i,k,j) = max(kappa_h, min(100., zkappa + kappa_h))
#  endif
# elif defined bryan_lewis_vertical
            diff_cbt(i,k,j) = Ahv(k)
# else
            diff_cbt(i,k,j) = kappa_h
# endif
          enddo
        enddo
      enddo

# if defined implicitvmix
      do ks=1,2

!       find density

        call statec (t(1,1,1,1,taum1), t(1,1,1,2,taum1), temp(1,1,jsmw)
     &,                jstrt, jend, istrt, iend, ks)

!       set large diffusion coefficient between unstable layers
!       (note: viscosity coefficient is not limited but could be here)

        do j=jstrt,jend
          do k=ks,kmm1,2
            do i=istrt,iend
              if (temp(i,k,j) .gt. temp(i,k+1,j)) then
                diff_cbt(i,k,j) = diff_cbt_limit*tmask(i,k+1,j)
              endif
            enddo
          enddo
        enddo
      enddo
# endif
#endif

#if defined ppvmix

!-----------------------------------------------------------------------
!     for momentum and tracers based on the pacanowski & philander
!     richardson mixing scheme (JPO vol 11, #11, 1981).
!-----------------------------------------------------------------------

      call ppmix (joff, js, je, istrt, iend)
#endif
#if defined uvic_kv_out
!-----------------------------------------------------------------------
!     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
#endif
#if defined isopycmix || defined redi_diffusion

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

#if defined trace_indices
      write (stdout,'(2x,5(a,i4))')
     & "=> In vmixc: js=",js," je=",je," joff=",joff
     &," jstrt=",jstrt," jend=",jend
#endif
      return
      end