#include "options.inc"

c#define weird_option

      module isopycnic_module
c
c-----------------------------------------------------------------------
c     Isopycnic diffusion module for SPFLAME taken from MOM 2
c     ( the old Redi/Cox code)
c
c     Original authors:
c             R.C.Pacanowski  (implemented Gokhans version for MOM 2)
c                              rcp@gfdl.gov
c             Gokhan Danabasoglu (pre MOM 2 version including
c                                 Gent/Mcwilliams transport velocities) 
c                                 gokhan@isis.cgd.ucar.edu
c             K. Dixon        (implemented Cox version in MOM 1)
c                              kd@gfdl.gov
c             M. Cox          (original isopycnal code in COX model)
c
c     subroutine isopycnic_init is called within routine setup
c
c     subroutine isopycnic is called within routine spflame
c
c     subroutine isopycnic_flux is called within tracer
c     see also some other modifications in tracer
c
c     subroutine isopycnic_add_K33 is called from vmixc.F
c
c     four parameters has to be specified in main namelist 
c
c
c     parameters are in spflame_module as well as a switch to
c     enable isopycnic diffusion
c
c                               aug 2001   c.eden
c-----------------------------------------------------------------------
c
      implicit none

#if defined  IRIX_host || defined ALPHA_host || defined LINUX_host || defined SUN_host 
      real,parameter :: epsln   = 1.0e-15  !  for single precision
#else
      real,parameter :: epsln   = 1.0e-25  !  for double precision
#endif

      integer, parameter :: nrpl = 11
      real :: dptlim(nrpl+1)=(/0.0e2, 500.0e2, 1000.0e2, 1500.0e2, 
     &                       2000.0e2, 2500.0e2, 3000.0e2, 3500.0e2, 
     &                       4000.0e2, 4500.0e2, 5000.0e2, 5500.0e2/)

      integer :: krplin(nrpl)=0
      integer, allocatable  :: kisrpl(:)

      real, allocatable :: fzisop(:)

      real, allocatable :: K1(:,:,:,:) 
      real, allocatable :: K2(:,:,:,:) 
      real, allocatable :: K3(:,:,:,:) 

      real, allocatable        :: adv_vetiso(:,:,:) 
      real, allocatable        :: adv_vntiso(:,:,:) 
      real, allocatable        :: adv_vbtiso(:,:,:) 

      real, allocatable :: dciso1(:,:,:) 
      real, allocatable :: dciso2(:,:,:) 

      real, allocatable :: rhoi(:,:,:,:) 
      real :: dslope=0.0008, slopec=0.007

      real, allocatable :: krit(:,:) 

      end module isopycnic_module

      subroutine isopycnic_init 
      use spflame_module
      use isopycnic_module
      implicit none
      integer :: k,m
      real ::  dptmid,t1,t2
      integer is,ie,js,je
      real :: ft1, delta_iso1  ,delta1a,delta1b
      integer :: i_delta1, j_delta1, k_delta1,i,j

      slopec=del_dm
      dslope=s_dm

      if (my_pe==0) then
       print*,''
       print*,' Initialization of isopycnal mixing scheme' 
       print*,'  (the old Redi/Cox code)'
       print*,''
       print*,'   ahisop  = ',ahisop, ' cm^2/s '
       print*,'   athkdf  = ',athkdf, ' cm^2/s '
       print*,'   slmx    = ',slmx
       print*,'   del_dm  = ',slopec
       print*,'   s_dm    = ',dslope
       print*,''
       if (enable_isopycnic_horvar) then
          print*,' horizontal varying thickness diffusion coeff.'
          print*,' and isopycnic diffusion coeff.'
          print*,' depending on local Rossby radius'
          print*,' (meant for eddy-permitting cases only)'
       endif
       if (enable_no_isopyc_medwater) then
        print*,''
        print*,
     &    ' decreasing isopycnal diffusivity in medwater outflow-area'
        print*,''
       endif
      endif
c
c   checks
c
      if (enable_rotated_grid .and.enable_isopycnic_horvar) then
         print*,' ERROR in isopycnic: horvar thicknes diff'
         print*,' does not work with rotated grid'
         call halt_stop('in isopycnic')
      endif
#ifdef partial_cell
         print*,' ERROR in isopycnic: partial cells '
         print*,' do not work with old Redi/Cox scheme'
         call halt_stop('in isopycnic')
#endif

c
c-----------------------------------------------------------------------
c     Compute the grid factors which set the maximum slopes available
c     for the mixing schemes. 
c-----------------------------------------------------------------------
c
      ft1 = 1.0/(4.0*max(ahisop,athkdf)*dt + epsln)
      delta_iso1  = dzt(1)*ft1*dxt(1)*cst(jmt/2)
      i_delta1 = 1
      j_delta1 = 1
      k_delta1 = 1
c
      do j=2,jmt-1
        do i=2,imt-1
	  do k=1,km
	    delta1a = dxt(i)*cst(j)*dzt(k)*ft1
	    delta1b = dyt(j)*dzt(k)*ft1
	    if ( delta_iso1 .ge. delta1a
     &      .or. delta_iso1 .ge. delta1b) then
              i_delta1 = i
              j_delta1 = j
              k_delta1 = k
	      delta_iso1  = min(delta1a,delta1b)
	    endif 
	  enddo
        enddo
      enddo
c
      if (my_pe==0) then
       print'(a)',' '
       print'(a,e14.7)',
     &'The diffusion grid factor delta_iso1 =',delta_iso1
       print'(a)', 'was determined at the grid point'
       print'(a,i4,a,e14.7)', 'dxt(',i_delta1,') = ',dxt(i_delta1)
       print'(a,i4,a,e14.7)', 'dyt(',j_delta1,') = ',dyt(j_delta1)
       print'(a,i4,a,e14.7)', 'dzt(',k_delta1,') = ',dzt(k_delta1)
       print'(a,i4,a,e14.7)', 'cst(',j_delta1,') = ',cst(j_delta1)
       print*,
     &'Without latitudinal filtering, delta_iso1 is the steepest'
       print*,
     &'isoneutral slope available for linear stab of Redi and GM.'
       print'(/a,e14.7/)',
     &'Maximum allowable isoneutral slope is specified as slmx = ',slmx
      endif


c
c-----------------------------------------------------------------------
c     determine the isopycnal reference pressure levels for the "t"
c     grid point levels, using the depths at the "t" grid points as the
c     reference depth (pressure)
c-----------------------------------------------------------------------
c
      dptlim(1)    = 0.
      dptlim(nrpl+1) = zw(km)
      allocate( kisrpl(km) )

      do k=1,km
        do m=2,nrpl+1
          if (zt(k).gt.dptlim(m-1) .and. zt(k).le.dptlim(m)) then
            kisrpl(k) = m-1
            go to 101
          endif
        enddo
101     continue
        if (kisrpl(k) .lt. 1 .or. kisrpl(k) .gt. nrpl) then
          print*,' =>Error: kisrpl is ',kisrpl(k),' at k ',k
          call halt_stop(' in isopycnic_init')
        endif
      enddo
c
c-----------------------------------------------------------------------
c     the indices used in isopycnal mixing indicating the location of
c     the reference pressure levels in the 20-level table of polynomial
c     expansion variables are computed
c
c     REMARK: because the polynomial expansion coefficients are
c             functions of the reference potential temperature and
c             salinity profiles, at the reference pressure level
c             the corresponding potential temperature and salinity
c             values will be used.
c-----------------------------------------------------------------------
c
c
      do m=2,nrpl+1
        dptmid = .5*(dptlim(m-1)+dptlim(m))
        if (dptmid .le. zt(1)) then
          krplin(m-1) = 1
        elseif (dptmid .gt. zt(km)) then
          krplin(m-1) = km
        elseif (dptmid.gt.zt(1) .and. dptmid.le.zt(km)) then
          do k=2,km
            if (zt(k) .ge. dptmid) then
              t1 = zt(k)-dptmid
              t2 = dptmid-zt(k-1)
              if (t1 .gt. t2) then
                krplin(m-1) = k-1
              else
                krplin(m-1) = k
              endif
              go to 102 
            endif
          enddo
102       continue
        endif
        if (krplin(m-1) .lt. 1 .or. krplin(m-1) .gt. km) then
          print*,' Error: krplin is ', krplin(m-1),' at m ',m-1
          call halt_stop(' in isopycnic_init')
        endif
      enddo
c
  96  format (/,' isopycnal reference pressure levels (kisrpl) = ',
     &        20(1x,i4))
  97  format (/,' reference pressure level indices (krplin) = ',
     &        20(1x,i4))

      if (my_pe==0) then
       print 96, (kisrpl(k),k=1,km)
       print 97, (krplin(m),m=1,nrpl)
      endif
c
c
      allocate( fzisop(km) ); fzisop=1.0
      do k=1,km
        fzisop(k)=1.0*(1.-exp(-zt(k)/100.e2))*
     &       ( 0.25+0.75*(1.-tanh((zt(k)-2000.e2)/1000.e2))/2. )
      enddo

      if (my_pe==0) then
       print'(a)', '  '
       print'(a)', ' Vertical structure function "fzisop(k)"='
       print'(5(1x,e12.6))', (fzisop(k),k=1,km)
      endif

c-----------------------------------------------------------------------
c     initialize arrays
c-----------------------------------------------------------------------
      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      allocate( rhoi(is-1:ie+1,km,js-1:je+1,nrpl));   rhoi=0.

      allocate( K1(is-1:ie,km,js  :je,3:3));      K1=0.
      allocate( K2(is  :ie,km,js-1:je,3:3));      K2=0.
      allocate( K3(is  :ie,km,js  :je,1:3));      K3=0.
      allocate( dciso1(is-1:ie,km,js  :je));      dciso1=0.
      allocate( dciso2(is  :ie,km,js-1:je));      dciso2=0.

      allocate( adv_vetiso(is-1:ie,  km,js  :je));adv_vetiso=0.
      allocate( adv_vntiso(is  :ie,  km,js-1:je));adv_vntiso=0.
      allocate( adv_vbtiso(is  :ie,0:km,js  :je));adv_vbtiso=0.
c
      allocate(krit(is_pe-1:ie_pe+1,js_pe-1:je_pe+1) ); krit=1.0
c
      if (enable_isopycnic_horvar)  then
c        call isopycnic_write_krit_init
      endif

      if (my_pe == 0) then
        print*,''
        print*,' done'
        print*,''
      endif

      end subroutine isopycnic_init




      subroutine isopycnic ()
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je,m,i,j,k
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     compute normalized densities for each isopycnal reference pressure
c     level using a 3rd order polynomial fit to the equation of state.
c     for each isopycnal reference pressure level, the same reference
c     potential temperature, reference salinity and expansion coeff
c     values are used at all of the vertical levels.
c-----------------------------------------------------------------------
c
      do m=1,nrpl
	do j=js-1,je+1
          do k=1,km
            call model_dens(t(is-1,k,j,1,taum1), t(is-1,k,j,2,taum1),
     &            rhoi(is-1,k,j,m), krplin(m) ,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,krplin(m),j)
#endif
     &                       )
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
      if (enable_isopycnic_horvar)  then
        call isopycnic_krit
c        if (snapshot_time_step) call isopycnic_write_krit
      endif
c-----------------------------------------------------------------------
c
      if (enable_no_isopyc_medwater) then
        call isopycnic_krit_med
      endif 
c
c-----------------------------------------------------------------------
c     evaluate K2(,,3) centered on the northern face of "T" cells
c-----------------------------------------------------------------------
c
      call isopycnic_k2_3 
c
c-----------------------------------------------------------------------
c     evaluate K1(,,3) centered on eastern face of "T" cells
c-----------------------------------------------------------------------
c
      call isopycnic_k1_3 
c
c-----------------------------------------------------------------------
c     evaluate K3(,,1..3) centered on bottom face of "T" cells
c-----------------------------------------------------------------------
c
      call isopycnic_k3_123
c
c-----------------------------------------------------------------------
c     compute isopycnal advective velocities for tracers
c-----------------------------------------------------------------------
c
      call isopycnic_adv 
c
      end subroutine isopycnic


      subroutine isopycnic_krit
c
c=======================================================================
c     compute isopycnal transport velocities 
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none

      integer is,ie,js,je,i,j,k
      real :: fxa,lm,cm,beta,Kr,dy
      real, allocatable :: ro1(:),ro2(:),bvfrqs(:,:,:)
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
      allocate(  ro1(is_pe-1:ie_pe+1),ro2(is_pe-1:ie_pe+1) )
      allocate( bvfrqs(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      ro1=0.;ro2=0.; bvfrqs=0.;
c
c    compute N2
c
      do j=js-1,je
        do k=1,km-1
         call model_dens(t(is-1,k  ,j,1,tau),t(is-1,k  ,j,2,tau),
     &                  ro1(is-1),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
         call model_dens(t(is-1,k+1,j,1,tau),t(is-1,k+1,j,2,tau),
     &                  ro2(is-1),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
         do i=is-1,ie
#ifdef partial_cell
          fxa = -grav/(rho0*dhwt(i,k,j))
#else
          fxa = -grav/rho0*dzwr(k)
#endif
          bvfrqs(i,k,j) =fxa*(ro1(i)-ro2(i))*tmask(i,k+1,j)
         enddo
        enddo
      enddo
       bvfrqs(:,km,:) = 0.  ! lower boundary condition
c
c    now N
c
      do j=js-1,je
        do k=1,km
         do i=is-1,ie
           if (bvfrqs(i,k,j) >=0.) then
             bvfrqs(i,k,j)=sqrt(bvfrqs(i,k,j))
           else
             bvfrqs(i,k,j)=0.
           endif
         enddo
        enddo
      enddo
c
c    now integrate N
c
      krit=0.0
      do j=js-1,je
        do k=1,km
         do i=is-1,ie
           krit(i,j)=krit(i,j)+bvfrqs(i,k,j)*dzt(k)
         enddo
        enddo
      enddo
c
c    now compute local Rossby radius approximated as in 
c    Chelton, et al. 1998, JPO 28 and compare with grid 
c
      do j=js,je
        beta  = 2.*omega*csu(j)/radius
        do i=is,ie
         CM = krit(i,j)/PI
         lm = MIN( CM/(ABS(cori(i,j,1))+epsln), sqrt(CM/2.0/BETA) )
         dy=max( dxt(i)*cst(j),dyt(j) )
         Kr = (LM/DY-isohv_a1)*isohv_a2 
         krit(i,j)=1.-( (exp(kr)-exp(-kr))/(exp(kr)+exp(-kr))+1.0)/2.0
        enddo
      enddo
      deallocate(ro1,ro2,bvfrqs)

      call border_exchg(krit,1,1)
      call set_cyclic(krit,1,1)

      end subroutine isopycnic_krit


      subroutine isopycnic_krit_med
c
c=======================================================================
c    special for HS4 model and Gibraltar region
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je,i,j,k
      real distance,dstnce_isopyc
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js-1,je
         do i=is-1,ie
c           distance = dstnce(yt(j),xt(i),yt(682),xt(112))
          distance = dstnce_isopyc(yt(j),xt(i),35.8585,269.2083)
           if (distance < 500.0) then
             if (distance <= 300.0) then
               krit(i,j)=krit(i,j)*0.00001
             else
               krit(i,j)=krit(i,j)*(200.0-(500.0-distance)*0.005)
             endif
           endif           
         enddo
       enddo
       end subroutine isopycnic_krit_med


       real function dstnce_isopyc(phi1,rla1,phi2,rla2)
c returns distance (in km) between two geographical
c points
c        written by   p.herrmann
c      phi1 : first latitude
c      rla1 : first longitude
c      phi2 : second latitude
c      rla2 : second longitude
       implicit none
       real :: p1,p2,rl1,rl2,x,xx, phi1,phi2,rla1,rla2,d
       real, parameter :: s=0.0174533 
       p1=phi1*s; p2=phi2*s; rl1=rla1*s; rl2=rla2*s
       x=sin(p1)*sin(p2)+cos(p1)*cos(p2)*cos(rl2-rl1)
       xx=abs(x)
       if(xx.gt.1.0) x=1.0
       d=atan(sqrt((1-x)/(1+x)))*222.24/s; dstnce_isopyc=d
      end function dstnce_isopyc


      subroutine isopycnic_k1_3 
c
c=======================================================================
c     
c     compute "K1(,,3)" at the center of the eastern face of "T" cells
c     use "c1e10" to keep the exponents in range.
c
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je,m,i,j,k
      real :: c1e10 = 1.0e10, fxd, fxe, olmask
      real :: fxa, fxb, fxc, slope
      real :: e(is_pe-1:ie_pe,km,js_pe:je_pe,3)
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     d(rho_barx_barz)/dz centered on eastern face of "t" cells
c     Note: values involving ocean surface and ocean bottom are
c           estimated afterwards using a linear extrapolation
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=2,km-1
          m = kisrpl(k)
          fxd = c1e10*0.5*dzt2r(k)
          do i=is-1,ie
            e(i,k,j,3) = fxd*(rhoi(i  ,k-1,j,m) - rhoi(i  ,k+1,j,m)
     &                       +rhoi(i+1,k-1,j,m) - rhoi(i+1,k+1,j,m))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     linearly extrapolate densities to ocean surface for calculation
c     of d(rho_barx_barz)/dz involving level 1.
c
c     REMARK: requires min(kmt(i,jrow)) = 2 cells in ocean.
c-----------------------------------------------------------------------
c
      k   = 1
      fxd = c1e10*dztr(k)
      fxe = dzw(k-1)+dzw(k)
      m   = kisrpl(k)
      do j=js,je
        do i=is-1,ie
          fxa        = 0.5*(rhoi(i,k+1,j,m) + rhoi(i+1,k+1,j,m))
          fxb        = 0.5*(rhoi(i,k,j,m) + rhoi(i+1,k,j,m))
          fxc        = dzwr(k)*(fxb*fxe - fxa*dzw(k-1))
          e(i,k,j,3) = fxd*(fxc - 0.5*(fxa+fxb))
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     linearly extrapolate densities to ocean bottom for calculation
c     of d(rho_barx_barz)/dz involving bottom level.
c-----------------------------------------------------------------------
c
      e(:,km,:,3) = 0.
c
      do j=js,je
        do i=is-1,ie
          k = min(kmt_big(i,j),kmt_big(i+1,j))
          if (k .ne. 0) then
            fxe        = dzw(k-1)+dzw(k)
            m          = kisrpl(k)
            fxa        = 0.5*(rhoi(i,k-1,j,m) + rhoi(i+1,k-1,j,m))
            fxb        = 0.5*(rhoi(i,k  ,j,m) + rhoi(i+1,k  ,j,m))
            fxc        = dzwr(k-1)*(fxb*fxe - fxa*dzw(k))
            e(i,k,j,3) = dztr(k)*c1e10*(0.5*(fxa+fxb) - fxc)
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     "e(,,,1)" = d(rho)/dx centered on east face of "T" cells
c     "e(,,,2)" = d(rho_barx_bary)/dy centered on east face of "T" cells
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          m = kisrpl(k)
          do i=is-1,ie
            e(i,k,j,1) = tmask(i,k,j)*tmask(i+1,k,j)*cstr(j)*dxur(i)
     &                     *c1e10*(rhoi(i+1,k,j,m) - rhoi(i,k,j,m))
            e(i,k,j,2) = dyt4r(j)*c1e10*(
     &                     rhoi(i  ,k,j+1,m) - rhoi(i  ,k,j-1,m)
     &                   + rhoi(i+1,k,j+1,m) - rhoi(i+1,k,j-1,m))
           enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     if any one of the 4 neighboring corner grid points is a land point,
c     set "e(i,k,j,2)" to zero. note that "e(i,k,j,2)" will be used
c     only in the slope check.
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
            olmask = tmask(i,k,j-1)*tmask(i,k,j+1)*tmask(i+1,k,j-1)
     &   	    *tmask(i+1,k,j+1)
            if (olmask .eq. 0.0)  e(i,k,j,2) = 0.0
           enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute "K1", using "slmxr" to limit vertical slope of isopycnal
c     to guard against numerical instabilities. 
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          do i=is-1,ie
            fxa = 0.
            fxb = sign(1.,e(i,k,j,3))/(abs(e(i,k,j,3))+epsln)
            slope = fxb*sqrt(e(i,k,j,1)**2+e(i,k,j,2)**2)
            if (slope .le. 0. .and. slope .ge. (-slmx)) then
              fxa = 0.5*(1.+tanh((slope+slopec)/dslope))
            endif
            dciso1(i,k,j) = fxa*fzisop(k)
            K1(i,k,j,3)  = -fxb*e(i,k,j,1)*dciso1(i,k,j)
          enddo
        enddo
      enddo
      end subroutine  isopycnic_k1_3



      subroutine isopycnic_k2_3 
c     
c=======================================================================
c     compute "K2(,,3)" at the center of the northern face of "T" cells
c     use "c1e10" to keep the exponents in range.
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je,m,i,j,k
      real :: c1e10 = 1.0e10, fxd, fxe, olmask
      real :: fxa, fxb, fxc, slope
      real :: e(is_pe:ie_pe,km,js_pe-1:je_pe,3)
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     d(rho_bary_barz)/dz centered on northern face of "T" cells
c     Note: values involving ocean surface and ocean bottom are
c           estimated afterwards using a linear extrapolation
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=2,km-1
          m = kisrpl(k)
          fxd = c1e10*0.5*dzt2r(k)
          do i=is,ie
            e(i,k,j,3) = fxd*(rhoi(i,k-1,j  ,m) - rhoi(i,k+1,j  ,m)
     &                       +rhoi(i,k-1,j+1,m) - rhoi(i,k+1,j+1,m))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     linearly extrapolate densities to ocean surface for calculation
c     of d(rho_bary_barz)/dz involving level 1.
c
c     REMARK: requires min(kmt(i,jrow)) = 2 cells in ocean.
c-----------------------------------------------------------------------
c
      k   = 1
      fxd = c1e10*dztr(k) 
      fxe = dzw(k-1)+dzw(k)
      m   = kisrpl(k)
      do j=js-1,je
        do i=is,ie
          fxa        = 0.5*(rhoi(i,k+1,j,m) + rhoi(i,k+1,j+1,m))
          fxb        = 0.5*(rhoi(i,k  ,j,m) + rhoi(i,k  ,j+1,m))
          fxc        = dzwr(k)*(fxb*fxe - fxa*dzw(k-1))
          e(i,k,j,3) = fxd*(fxc - 0.5*(fxa+fxb))
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     linearly extrapolate densities to ocean bottom for calculation
c     of d(rho_bary_barz)/dz involving bottom level.
c-----------------------------------------------------------------------
c
      e(:,km,:,3) = 0.
c
      do j=js-1,je
        do i=is,ie
          k = min(kmt_big(i,j),kmt_big(i,j+1))
          if (k .ne. 0) then
            fxe        = dzw(k-1)+dzw(k)
            m          = kisrpl(k)
            fxa        = 0.5*(rhoi(i,k-1,j,m) + rhoi(i,k-1,j+1,m))
            fxb        = 0.5*(rhoi(i,k  ,j,m) + rhoi(i,k  ,j+1,m))
            fxc        = dzwr(k-1)*(fxb*fxe - fxa*dzw(k))
            e(i,k,j,3) = dztr(k)*c1e10*(0.5*(fxa+fxb) - fxc)
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     "e(,,,1)" = d(rho_barx_bary)/dx centered on north face of "T" cells
c     "e(,,,2)" = d(rho)/dy on north face of "T" cells
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          m = kisrpl(k)
          do i=is,ie
            e(i,k,j,1) = csur(j)*dxt4r(i)*c1e10*(
     &                    rhoi(i+1,k,j+1,m) - rhoi(i-1,k,j+1,m)
     &                  + rhoi(i+1,k,j  ,m) - rhoi(i-1,k,j  ,m))
            e(i,k,j,2) = tmask(i,k,j)*tmask(i,k,j+1)*dyur(j)*c1e10
     &                  *(rhoi(i,k,j+1,m) - rhoi(i,k,j,m))   
           enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     if any one of the 4 neighboring corner grid points is a land point,
c     set "e(i,k,j,1)" to zero. note that "e(i,k,j,1)" will be used
c     only in the slope check.
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
            olmask = tmask(i-1,k,j+1)*tmask(i+1,k,j+1)*tmask(i-1,k,j)
     &              *tmask(i+1,k,j)
            if (olmask == 0.)  e(i,k,j,1) = 0.
           enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute "K2", using "slmxr" to limit vertical slope of isopycnal
c     to guard against numerical instabilities. 
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
            fxa = 0.
            fxb = sign(1.,e(i,k,j,3))/(abs(e(i,k,j,3))+epsln)
            slope = fxb*sqrt(e(i,k,j,1)**2+e(i,k,j,2)**2)
            if (slope .le. 0. .and. slope .ge. (-slmx)) then
              fxa = 0.5*(1.+tanh((slope+slopec)/dslope))
            endif
            dciso2(i,k,j) = fxa *fzisop(k)
            K2(i,k,j,3)  = -fxb*e(i,k,j,2)*dciso2(i,k,j)
          enddo
        enddo
      enddo
      end subroutine isopycnic_k2_3







      subroutine isopycnic_k3_123 
c
c=======================================================================
c     compute K2(,,,1:3) at the center of the bottom face of "T" cells
c     use "c1e10" to keep the exponents in range.
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je,m,i,j,k
      real :: c1e10 = 1.0e10, fxd, fxe, olmask
      real :: fxa, fxb, fxc, slope
      real :: e(is_pe:ie_pe,km,js_pe:je_pe,3)

c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
      do j=js,je
        do k=2,km
          m = kisrpl(k)
          do i=is,ie
            e(i,k-1,j,1) = cstr(j)*dxt4r(i)*c1e10
     &               *(tmask(i-1,k-1,j)*tmask(i,k-1,j)*(rhoi(i,k-1,j,m)
     &               -rhoi(i-1,k-1,j,m))
     &               +tmask(i,k-1,j)*tmask(i+1,k-1,j)*(rhoi(i+1,k-1,j,m)
     &               -rhoi(i,k-1,j,m))
     &               +tmask(i-1,k,j)*tmask(i,k,j)*(rhoi(i,k,j,m)
     &               -rhoi(i-1,k,j,m))
     &               +tmask(i,k,j)*tmask(i+1,k,j)*(rhoi(i+1,k,j,m)
     &               -rhoi(i,k,j,m)))
            e(i,k-1,j,2) = dyt4r(j)*c1e10
     &               *(tmask(i,k-1,j-1)*tmask(i,k-1,j)*(rhoi(i,k-1,j,m)
     &               -rhoi(i,k-1,j-1,m))
     &               +tmask(i,k-1,j)*tmask(i,k-1,j+1)*(rhoi(i,k-1,j+1,m)
     &               -rhoi(i,k-1,j,m))
     &               +tmask(i,k,j-1)*tmask(i,k,j)*(rhoi(i,k,j,m)
     &               -rhoi(i,k,j-1,m))
     &               +tmask(i,k,j)*tmask(i,k,j+1)*(rhoi(i,k,j+1,m)
     &               -rhoi(i,k,j,m)))
            e(i,k-1,j,3) = dzwr(k-1)*tmask(i,k-1,j)*tmask(i,k,j)*c1e10
     &               *(rhoi(i,k-1,j,m) - rhoi(i,k,j,m))
           enddo
        enddo
        k = km
	e(:,k,j,1) = 0.
	e(:,k,j,2) = 0.
        e(:,k,j,3) = 0.
      enddo
c
c-----------------------------------------------------------------------
c     compute "K3", using "slmxr" to limit vertical slope of isopycnal
c     to guard against numerical instabilities.  
c-----------------------------------------------------------------------
c

      do j=js,je
        do k=1,km
          do i=is,ie
            fxa = 0.
            fxb = sign(1.,e(i,k,j,3))/(abs(e(i,k,j,3))+epsln)
            slope = fxb*sqrt(e(i,k,j,1)**2+e(i,k,j,2)**2)
            if (slope .le. 0. .and. slope .ge. (-slmx)) then
              fxa = 0.5*(1.+tanh((slope+slopec)/dslope))
            endif
            fxc = fxb*fxa
ce            scaling with fzisop
ce            same form as without isopycmixspatialvar
ce            correct ??
     &         *0.5*(fzisop(min(km,k+1))+fzisop(k))
            K3(i,k,j,1) = -fxc*e(i,k,j,1)
            K3(i,k,j,2) = -fxc*e(i,k,j,2)
            K3(i,k,j,3) = fxb*fxb*(e(i,k,j,1)**2+e(i,k,j,2)**2)
     &                    *fxa
ce             fxc is not used, so scale fxa again with fzisop
     &         *0.5*(fzisop(min(km,k+1))+fzisop(k))
          enddo
        enddo
      enddo
c
      end subroutine isopycnic_K3_123




      subroutine isopycnic_adv
c
c=======================================================================
c     compute isopycnal transport velocities 
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none

      integer is,ie,js,je,m,i,j,k
      real :: fxa,fxb
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     compute the meridional component of the isopycnal mixing velocity
c     at the center of the northern face of the "t" cells.
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=2,km-1
          fxa = -dzt2r(k)*athkdf*csu(j)
          do i=is,ie
            adv_vntiso(i,k,j) = fxa*tmask(i,k,j)*tmask(i,k,j+1)*(
     &                          K2(i,k-1,j,3) - K2(i,k+1,j,3))
     &                          *krit(i,j)
          enddo
        enddo
      enddo
c
c     consider the top and bottom levels. "K2" is assumed to be zero
c     at the ocean top and bottom.
c
      k = 1
      fxa = -dzt2r(k)*athkdf
      do j=js-1,je
        do i=is,ie
          adv_vntiso(i,k,j) = -fxa*tmask(i,k,j)*tmask(i,k,j+1)*csu(j)
     &                        *(K2(i,k,j,3) + K2(i,k+1,j,3))
     &                        *krit(i,j)
        enddo
      enddo
c
      adv_vntiso(:,km,:) = 0.
c
      do j=js-1,je
        do i=is,ie
          k = min(kmt_big(i,j),kmt_big(i,j+1))
          if (k .ne. 0) then
            adv_vntiso(i,k,j) = -dzt2r(k)*athkdf*csu(j)*tmask(i,k,j)
     &                                          *krit(i,j)
     &                   *tmask(i,k,j+1)*(K2(i,k,j,3) + K2(i,k-1,j,3))
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute the zonal component of the isopycnal mixing velocity
c     at the center of the eastern face of the "t" grid box.
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=2,km-1
          fxa = -dzt2r(k)*athkdf
          do i=is-1,ie
            adv_vetiso(i,k,j) = fxa*tmask(i,k,j)*tmask(i+1,k,j)
     &                          *krit(i,j)
     &                          *(K1(i,k-1,j,3) - K1(i,k+1,j,3))
          enddo
        enddo
      enddo
c
c     consider the top and bottom levels. "K1" is assumed to be zero
c     at the ocean top and bottom.
c
      k = 1
      fxa = -dzt2r(k)*athkdf
      do j=js,je
        do i=is-1,ie
          adv_vetiso(i,k,j) = -fxa*tmask(i,k,j)*tmask(i+1,k,j)
     &                        *(K1(i,k,j,3)+K1(i,k+1,j,3))
     &                        *krit(i,j)
        enddo
      enddo
c
      adv_vetiso(:,km,:) = 0.
c
      do j=js,je
        do i=is-1,ie
          k = min(kmt_big(i,j),kmt_big(i+1,j))
          if (k .ne. 0) then
            adv_vetiso(i,k,j) = -dzt2r(k)*athkdf*tmask(i,k,j)
     &                     *tmask(i+1,k,j)*(K1(i,k,j,3)+K1(i,k-1,j,3))
     &                     *krit(i,j)
          endif
        enddo
      enddo

#ifdef weird_option
      do j=js,je  ! note the error in boundaries
       do k=1,km
        do i=is,ie
            fxa = adv_vntiso(i,k,j) 
            fxb = adv_vetiso(i,k,j) 
            adv_vntiso(i,k,j) = fxa + fxb
            adv_vetiso(i,k,j) = fxb - fxa
        enddo
       enddo
      enddo
#endif
c
c
c-----------------------------------------------------------------------
c     compute the vertical component of the isopycnal mixing velocity
c     at the center of the bottom face of the "t" cells, using the
c     continuity equation for the isopycnal mixing velocities
c-----------------------------------------------------------------------
c
      adv_vbtiso(:,0,:) = 0.
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            adv_vbtiso(i,k,j) = dzt(k)*cstr(j)*(
     &      (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(j)) 
          enddo
        enddo
      enddo
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            adv_vbtiso(i,k,j) = adv_vbtiso(i,k,j) + adv_vbtiso(i,k-1,j)
          enddo
        enddo
      enddo
c
      do j=js,je
        do i=is,ie
          adv_vbtiso(i,kmt_big(i,j),j) = 0.
        enddo
      enddo
c
      end subroutine isopycnic_adv



      subroutine isopycnic_flux (m,diff_fe,diff_fn,diff_fb)
c
c=======================================================================
c     isopycnal diffusive tracer fluxes are computed.
c     horizontal fluxes are added (!) to diff_fe and diff_fn
c     while diff_fb is set (!) to K13 and K23 components which
c     are solved explicitly. 
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer, intent(in) :: m
      real diff_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real diff_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real diff_fb(is_pe:ie_pe,0:km,js_pe:je_pe)

      integer is,ie,js,je,i,j,k
      real :: fxa,fxb,fxc,fxe
      real temp(is_pe-1:ie_pe,km,js_pe-1:je_pe)
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     first compute the vertical tracer flux "temp" at the northern
c     face of "t" cells.
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=2,km-1
          do i=is,ie
            temp(i,k,j) = 0.5*dzt2r(k)*(
     &                      t(i,k-1,j+1,m,taum1) - t(i,k+1,j+1,m,taum1)
     &                     +t(i,k-1,j  ,m,taum1) - t(i,k+1,j  ,m,taum1))
     &                     
           enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     now consider the top level, assuming that the surface tracer 
c     values are the same as the ones at "k"=1
c-----------------------------------------------------------------------
c
      k = 1
      do j=js-1,je
        do i=is,ie
          temp(i,k,j) = 0.5*dzt2r(k)*(
     &                     t(i,k,j+1,m,taum1) - t(i,k+1,j+1,m,taum1)
     &                    +t(i,k,j  ,m,taum1) - t(i,k+1,j  ,m,taum1))
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     finally, consider the bottom level. the extrapolative estimator
c     is used to compute the tracer values at the ocean bottom.
c-----------------------------------------------------------------------
c
      temp(:,km,:) = 0.
c
      do j=js-1,je
        do i=is,ie
          k = min(kmt_big(i,j),kmt_big(i,j+1))
          if (k .ne. 0) then
            fxe = dzw(k-1)+dzw(k)
            fxa = 0.5*(t(i,k-1,j+1,m,taum1) + t(i,k-1,j,m,taum1))
            fxb = 0.5*(t(i,k  ,j+1,m,taum1) + t(i,k  ,j,m,taum1))
            fxc = dzwr(k-1)*(fxb*fxe-fxa*dzw(k))
            temp(i,k,j) = dztr(k)*(0.5*(fxa+fxb) - fxc)
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute of meridional tracer flux at northern face of "T" cells.
c     first calculate the effects of purely horizontal diffusion, using
c     the sum of the along isopycnal and background horizontal diffusion
c     coefficients as the total diffusion coefficient.     
c     add in the effects of the along isopycnal diffusion computed
c     using "K2" component of the tensor and apply land/sea masks
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          do i=is,ie
            fxa = ahisop*dciso2(i,k,j)
            fxa = fxa*krit(i,j)
            diff_fn(i,k,j) = diff_fn(i,k,j)+(fxa*dyur(j)*(
     &                        t(i,k,j+1,m,taum1)-t(i,k,j,m,taum1))
     &                        + ahisop*K2(i,k,j,3)*temp(i,k,j)
     &                          *krit(i,j)
     &                       )*csu(j)*tmask(i,k,j+1)*tmask(i,k,j)
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute the vertical tracer flux "temp" at the eastern
c     face of "t" cells.
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=2,km-1
          do i=is-1,ie
            temp(i,k,j) = 0.5*dzt2r(k)*(
     &                    t(i+1,k-1,j,m,taum1) - t(i+1,k+1,j,m,taum1)
     &                   +t(i  ,k-1,j,m,taum1) - t(i  ,k+1,j,m,taum1))
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     now consider the top level, assuming that the surface tracer 
c     values are the same as the ones at "k"=1
c-----------------------------------------------------------------------
c
      k = 1
      do j=js,je
        do i=is-1,ie
          temp(i,k,j) = 0.5*dzt2r(k)*(
     &                  t(i+1,k,j,m,taum1) -t(i+1,k+1,j,m,taum1)
     &                 +t(i  ,k,j,m,taum1) -t(i  ,k+1,j,m,taum1))
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     finally, consider the bottom level. the extrapolative estimator
c     is used to compute the tracer values at the ocean bottom.
c-----------------------------------------------------------------------
c
      temp(:,km,:) = 0.
c
      do j=js,je
        do i=is-1,ie
          k = min(kmt_big(i,j),kmt_big(i+1,j))
          if (k .ne. 0) then
            fxe          = dzw(k-1)+dzw(k)
            fxa          = 0.5*(t(i,k-1,j,m,taum1)+t(i+1,k-1,j,m,taum1))
            fxb          = 0.5*(t(i,k,j,m,taum1)+t(i+1,k,j,m,taum1))
            fxc          = dzwr(k-1)*(fxb*fxe - fxa*dzw(k))
            temp(i,k,j) = dztr(k)*(0.5*(fxa+fxb)-fxc)
          endif
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute of zonal tracer flux at eastern face of "T" cells.
c     first calculate the effects of purely horizontal diffusion, using
c     the sum of the along isopycnal and background horizontal diffusion
c     coefficients as the total diffusion coefficient.     
c     add in the effects of the along isopycnal diffusion computed
c     using "K1" component of the tensor and apply land/sea masks      
c-----------------------------------------------------------------------
c      
      do j=js,je
        do k=1,km
          do i=is-1,ie
            fxa = ahisop*dciso1(i,k,j)
            fxa = fxa*krit(i,j)
            diff_fe(i,k,j) = diff_fe(i,k,j)+(fxa*cstr(j)*dxur(i)*(
     &                        t(i+1,k,j,m,taum1) - t(i,k,j,m,taum1))
     &                        + ahisop*K1(i,k,j,3)*temp(i,k,j)
     &                          *krit(i,j)
     &                       )*tmask(i+1,k,j)*tmask(i,k,j)
          enddo
        enddo
      enddo

c
c-----------------------------------------------------------------------
c     compute the vertical tracer flux "diff_fbiso" containing the K31
c     and K32 components which are to be solved explicitly. The K33
c     component will be treated semi-implicitly by diff_fb.
c-----------------------------------------------------------------------
c

      do j=js,je
        do k=2,km
          do i=is,ie
            diff_fb(i,k-1,j) = 
     &                 tmask(i,k-1,j)*tmask(i,k,j)
     &                  *(ahisop*
     &                           krit(i,j)*
     &                           K3(i,k-1,j,1)*cstr(j)*dxt4r(i)*(
     &      tmask(i-1,k,j)*(t(i,k,j,m,taum1)-t(i-1,k,j,m,taum1))
     &     +tmask(i-1,k-1,j)*(t(i,k-1,j,m,taum1)-t(i-1,k-1,j,m,taum1))
     &     +tmask(i+1,k,j)*(t(i+1,k,j,m,taum1)-t(i,k,j,m,taum1))
     &     +tmask(i+1,k-1,j)*(t(i+1,k-1,j,m,taum1)-t(i,k-1,j,m,taum1)))
     &                   +ahisop*
     &                           krit(i,j)*
     &                           K3(i,k-1,j,2)*dyt4r(j)*(
     &      tmask(i,k,j-1)*(t(i,k,j,m,taum1)-t(i,k,j-1,m,taum1))
     &     +tmask(i,k-1,j-1)*(t(i,k-1,j,m,taum1)-t(i,k-1,j-1,m,taum1))
     &     +tmask(i,k,j+1)*(t(i,k,j+1,m,taum1)-t(i,k,j,m,taum1))
     &     +tmask(i,k-1,j+1)*(t(i,k-1,j+1,m,taum1)-t(i,k-1,j,m,taum1)))
     &                   )
          enddo
        enddo
      enddo
      end subroutine isopycnic_flux




      subroutine isopycnic_add_K33
c
c-----------------------------------------------------------------------
c     add K33 component of isopycnal slope to vertical diffusion coeff
c-----------------------------------------------------------------------
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer is,ie,js,je
      integer i,j,k

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
        do k=1,km-1
          do i=is,ie
            diff_cbt(i,k,j) = (ahisop*K3(i,k,j,3)
     &                          *krit(i,j)
     &                         + diff_cbt(i,k,j) )*tmask(i,k+1,j)


          enddo
        enddo
      enddo
      end subroutine isopycnic_add_K33



      subroutine isopycnic_add_gm(n,adv_fe,adv_fn,adv_fb)
c
c=======================================================================
c     add isopycnal transport velocities to advective fluxes
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer, intent(in) :: n
      real :: adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real :: adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real :: adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
      integer :: is,ie,js,je,i,j,k
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate 2*advective flux across eastern face of "T" cells.
c-----------------------------------------------------------------------
c
       do j=js,je
        do k=1,km
         do i=is-1,ie
	   adv_fe(i,k,j) = adv_fe(i,k,j)+
     &                     adv_vetiso(i,k,j)*(t(i,  k,j,n,taum1) + 
     &                                        t(i+1,k,j,n,taum1))
         enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c    calculate 2*advective flux across northern face of "T" cells 
c-----------------------------------------------------------------------
c
      do j=js-1,je
       do k=1,km
        do i=is,ie
	 adv_fn(i,k,j) = adv_fn(i,k,j)+
     &                   adv_vntiso(i,k,j)*(t(i,k,j  ,n,taum1) + 
     &                                      t(i,k,j+1,n,taum1))
        enddo
       enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute advective tracer flux at the center of the bottom face of
c     the "T" cells
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            adv_fb(i,k,j) = adv_fb(i,k,j)+adv_vbtiso(i,k,j)*
     &                         (t(i,k,j,n,taum1) + t(i,k+1,j,n,taum1))
          enddo
        enddo
      enddo
      end subroutine isopycnic_add_gm


      subroutine isopycnic_add_div_gm(n,is,ie,js,je)
c
c=======================================================================
c     add divergence of bolus transport fluxes
c=======================================================================
c
      use spflame_module
      use isopycnic_module
      implicit none
      integer, intent(in) :: n,is,ie,js,je
      integer i,j,k
      real adv_fbiso(is:ie,0:km,js:je),xh(is:ie)
      real ADV_Tx,ADV_Ty,ADV_Tz

      ADV_Tx(i,k,j) = cstdxt2r(i,j)*
     &                   (adv_vetiso(i,k,j)*
     &                    (t(i+1,k,j,n,taum1) + t(i,k,j,n,taum1))
     &                   -adv_vetiso(i-1,k,j)*
     &                    (t(i,k,j,n,taum1) + t(i-1,k,j,n,taum1)))
      ADV_Ty(i,k,j) = cstdyt2r(j)*
     &                   (adv_vntiso(i,k,j)*
     &                   (t(i,k,j+1,n,taum1) + t(i,k,j,n,taum1))
     &                  -adv_vntiso(i,k,j-1)*
     &                   (t(i,k,j,n,taum1) + t(i,k,j-1,n,taum1)))   
      ADV_Tz(i,k,j) = dzt2r(k)*(adv_fbiso(i,k-1,j)-adv_fbiso(i,k,j))
c
c-----------------------------------------------------------------------
c     compute advective tracer flux at the center of the bottom face of
c     the "T" cells
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            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
c
c     now consider the top and bottom boundaries
c
      adv_fbiso(:,0,:)  = 0.
      adv_fbiso(:,km,:) = 0.
c
      do j=js,je
       do k=1,km
        do i=is, ie
         xh(i) = -ADV_Tx(i,k,j) - ADV_Ty(i,k,j)
        enddo
        do i=is, ie
         xh(i) = xh(i) - ADV_Tz(i,k,j)
        enddo
        do i=is, ie
         t(i,k,j,n,taup1)=t(i,k,j,n,taup1)+xh(i)*c2dt*tmask(i,k,j)
        enddo
       enddo
      enddo
      end subroutine isopycnic_add_div_gm





#ifdef notdef

      subroutine isopycnic_write_krit_init
      use spflame_module
      use isopycnic_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: ncid,iret,lat_tdim,lon_tdim,itimedim
      integer :: dims(4),lon_tid,lat_tid,itimeid,kritid
      real :: spval=-9.9e12
      character name*24, unit*16, text*80

      if (my_pe==0) then
        print*,' Writing krit to  NetCDF output file '
        ncid = nccre ('krit.cdf', NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
c       dimensions
        lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
        Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
        iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c       grid variables
        dims(1)  = Lon_tdim
        Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
        dims(1)  = Lat_tdim
        Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
        dims(1)  = iTimedim
        iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
c       attributes of the grid
        name = 'Longitude on T grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on T grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Time                    '; unit = 'days            '
        call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
        call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c       variables
        dims=(/lon_tdim,lat_tdim,itimedim,1/)
        kritid=ncvdef (ncid,'krit', NCFLOAT,3,dims,iret)
        call dvcdf(ncid,kritid,'Kriterium',9,' ',1,spval)
        call ncendf(ncid, iret)
        call ncvpt(ncid, Lon_tid, 1, imt,xt, iret)
        call ncvpt(ncid, Lat_tid, 1, jmt,yt, iret)
        call ncclos (ncid, iret)
      endif
#endif
      end subroutine isopycnic_write_krit_init

      subroutine isopycnic_write_krit
      use spflame_module
      use isopycnic_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: ncid,iret,npe
      integer :: kritid,itdimid,m,start(4),count(4),itimeid
      real :: spval=-9.9e12,fxa

      do npe=0,n_pes
       call barrier
       if (my_pe==npe) then

        iret=nf_open('krit.cdf',NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_varid(ncid,'krit',kritid)
        iret=nf_inq_varid(ncid,'Time',itimeid)

        iret=nf_inq_dimid(ncid,'Time',itdimid)
        iret=nf_inq_dimlen(ncid, itdimid,m)
        if (my_pe==0) then
         m=m+1
         call read_stamp(current_stamp,fxa)
         print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',fxa,')',
     &          ' (time steps in file : ',m,')'
         iret= nf_put_vara_real (ncid,itimeid,m,1,fxa)
        endif
        start=(/is_pe,js_pe,m,1/)
        count=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret=nf_put_vara_real(ncid,kritid,start,count,
     &           krit(is_pe:ie_pe,js_pe:je_pe))
        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
#endif
      end subroutine isopycnic_write_krit

#endif

