#include "options.inc"

c
c local defines for isopycnal mixing
c
c


c standard MOM tapering of diff. in case of strong slopes
#define dm_taper
c gerdes koeberle willebrand tapering, however does not work good
c#define gkw_taper

c disable vertical tapering for ahsteep
c does not work for athkdf = ahsteep
#define disable_vertical_tapering_ahsteep
c
c disable entire module for quicker compilation
c
c#define disable_isoneutral_module

c
      module isoneutral_module
c
c-----------------------------------------------------------------------
c     Isoneutral diffusion module for SPFLAME taken from MOM 3
c
c     Original authors:
c           R.C.Pacanowski   rcp@gfdl.gov 
c           S.M.Griffies     smg@gfdl.gov 
c
c     subroutine isoneutral_init is called within routine setup
c
c     subroutine isoneutral is called within routine spflame
c
c     subroutine isoneutralflux is called within tracer
c     see also some other modifications in tracer
c
c     subroutine isoneutral_add_K33 is called from vmixc.F
c
c     subroutine isoneutral_adv is called from diag_snap
c
c     four parameters has to be specified in main namelist 
c     parameters are in spflame_module as well as a switch to
c     enable isoneutral diffusion
c
c                               Jul 2001   c.eden
c-----------------------------------------------------------------------
c

#ifndef disable_isoneutral_module

c      private

#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-20  ! for double precision
#endif

c-----------------------------------------------------------------------
c Arrays and scalars used for all isoneutral mixing schemes.
c-----------------------------------------------------------------------

      real, allocatable :: drho(:,:,:,:)
      real, allocatable :: ddxt(:,:,:,:)
      real, allocatable :: ddyt(:,:,:,:)
      real, allocatable :: ddzt(:,:,:,:)
      real, allocatable :: fzisop(:),agm(:,:)
      real redi_gm(6)
c
c-----------------------------------------------------------------------
c Arrays for either Redi diffusion or Gent McWilliams stirring
c     Ai_ez  = slope*diffusion coefficient on eastern face of T cell
c     Ai_nz  = slope*diffusion coefficient on northern face of T cell
c     Ai_bx  = slope*diffusion coefficient on bottom face of T cell
c     Ai_by  = slope*diffusion coefficient on bottom face of T cell
c-----------------------------------------------------------------------
c
      real, allocatable :: Ai_ez(:,:,:,:,:)
      real, allocatable :: Ai_nz(:,:,:,:,:)
      real, allocatable :: Ai_bx(:,:,:,:,:)
      real, allocatable :: Ai_by(:,:,:,:,:)
c
c-----------------------------------------------------------------------
c Arrays and scalars for just Redi diffusion
c
c Arrays
c     K11,K22,K33 = diagonal components to the Redi diffusion tensor.
c-----------------------------------------------------------------------
c
      real, allocatable :: K11(:,:,:)
      real, allocatable :: K22(:,:,:)
      real, allocatable :: K33(:,:,:)
#ifdef partial_cell
c
c     twice the thickness of quarter cells within T-cells
c     0 is for upper quarter cell and 1 is for lower quarter cell
c
      real, allocatable         :: delqc(:,:,:,:)
#endif
c
#endif

      end module isoneutral_module
c
      subroutine isoneutral_init
c
c=======================================================================
c
c      Initialization for isoneutral mixing schemes.  Each scheme
c      can be turned on alone, or in combination with any of 
c      the others.
c           R.C.Pacanowski   rcp@gfdl.gov 
c           S.M.Griffies     smg@gfdl.gov 
c=======================================================================
c
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_module
      implicit none
      real :: c0=0.,c1=1.0,c2=2.0
      integer :: stdout=6
      integer is,ie,js,je
      integer i,j,k
      real delta1a,delta1b,ft1,delta_iso1
      integer i_delta1,j_delta1,k_delta1

      if (my_pe==0) then
       print*,''
       print*,' Initialization of isoneutral mixing scheme' 
       print*,''
       print*,'   ahisop  = ',ahisop, ' cm^2/s '
       print*,'   athkdf  = ',athkdf, ' cm^2/s '
       print*,'   ahsteep = ',ahsteep,' cm^2/s '
       print*,'   slmx    = ',slmx
       print*,'   del_dm  = ',del_dm
       print*,'   s_dm    = ',s_dm
       print*,''
#if defined gkw_taper
       print*, ' Gerdes, Koberle, and Willebrand taper is used.'
#endif
#if defined dm_taper
       print*,'  Danabasoglu and McWilliams taper is used.'
#endif
#ifdef disable_vertical_tapering_ahsteep
       print*,' vertical tapering disabled for ahsteep '
       if (athkdf == ahisop .and. athkdf == ahsteep) then
        print*,' WARNING: does not work for',
     &    ' ahisop=athkdf = ahsteep'
       endif
#endif
       print*,''
      endif
c
c-----------------------------------------------------------------------
c     Compute the grid factors which set the maximum slopes available
c     for the mixing schemes. 
c-----------------------------------------------------------------------
c
      ft1 = c1/(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
      write(stdout,'(a)')
      write(stdout,'(a,e14.7)')
     &'The diffusion grid factor delta_iso1 =',delta_iso1
      write(stdout,'(a)')
     &'was determined at the grid point'
      write(stdout,'(a,i4,a,e14.7)')
     &'dxt(',i_delta1,') = ',dxt(i_delta1)
      write(stdout,'(a,i4,a,e14.7)')
     &'dyt(',j_delta1,') = ',dyt(j_delta1)
      write(stdout,'(a,i4,a,e14.7)')
     &'dzt(',k_delta1,') = ',dzt(k_delta1)
      write(stdout,'(a,i4,a,e14.7)')
     &'cst(',j_delta1,') = ',cst(j_delta1)
      write(stdout,*)
     &'Without latitudinal filtering, delta_iso1 is the steepest'
      write(stdout,*)
     &'isoneutral slope available for linear stab of Redi and GM.'
      write(stdout,'(/a,e14.7/)')
     & 'Maximum allowable isoneutral slope is specified as slmx = ',slmx
      endif
c
c-----------------------------------------------------------------------
c     vertical structure function for the isoneutral diffusion coeff.
c-----------------------------------------------------------------------
c
      allocate(fzisop(km))
c      fzisop = c1   ! change this for a vertical dependency
      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
c
      if (my_pe==0) then
      write (stdout,'(a)') '  '
      write (stdout,'(a)') ' Vertical structure function "fzisop(k)"='
      write (stdout,'(5(1x,e12.6))') (fzisop(k),k=1,km)
      endif
c
c-----------------------------------------------------------------------
c     initialize arrays
c-----------------------------------------------------------------------
c
      is=max(is_pe,1); ie=min(ie_pe,imt)
      js=max(1,js_pe); je=min(je_pe,jmt)

      allocate( drho(is-1:ie+1,  km,js-1:je+1,2) )
      drho=c0
      allocate( ddxt(is-1:ie+1,  km,js-1:je+1,2) )
      ddxt=c0
      allocate( ddyt(is-1:ie+1,  km,js-1:je+1,2) )
      ddyt=c0
      allocate( ddzt(is-1:ie+1,0:km,js-1:je+1,2) )
      ddzt=c0

      allocate( Ai_ez(is-1:ie,km,js:je,0:1,0:1) )
      Ai_ez=c0
      allocate( K11(is-1:ie,km,js:je) )
      K11=c0

      allocate( Ai_bx(is:ie,km,js:je,0:1,0:1) )
      Ai_bx=c0
      allocate( Ai_by(is:ie,km,js:je,0:1,0:1) )
      Ai_by=c0
      allocate( K33(is:ie,km,js:je) )
      K33=c0

      allocate( Ai_nz(is:ie,km,js-1:je,0:1,0:1) )
      Ai_nz=c0
      allocate( K22(is:ie,km,js-1:je) )
      K22=c0

      redi_gm = c0
c!!!!!!!!!!!!!!!!  here was a bug
      redi_gm(1) = c1
      redi_gm(2) = c1
c
      if(athkdf .eq. ahisop) then
        redi_gm(3) = c0
        redi_gm(4) = c1
      else
        redi_gm(3) = c1
        redi_gm(4) = c0
      endif
      if(athkdf .eq. ahsteep) then
        redi_gm(5) = c0
        redi_gm(6) = c1
      else
        redi_gm(5) = c1
        redi_gm(6) = c0
      endif
      if (my_pe==0) then
      write(stdout,*)' '
      write(stdout,*)'The Redi/GM algorithm coefficients are:'
      do i=1,6
        write(stdout,*)'redi_gm(',i,') = ',redi_gm(i)
      enddo
      endif
c
      allocate( agm(is-1:ie,js-1:je) )
      agm=athkdf   ! change this for horizontal dependency
c
#ifdef partial_cell

      allocate(delqc(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1,0:1) )
c     twice the thickness of quarter cells within T-cells
c     0 is for upper quarter cell and 1 is for lower quarter cell
      delqc=0.

      do j=js_pe,je_pe
         do k=1,km
	  do i=is_pe,ie_pe
            delqc(i,k,j,0) = 2.*fracdz(k,0)*dht(i,k,j)
            delqc(i,k,j,1) = 2.*fracdz(k,1)*dht(i,k,j)
	  enddo
	 enddo
      enddo

      call border_exchg(delqc(:,:,:,0),km,1)
      call border_exchg(delqc(:,:,:,1),km,1)
      if (cyclic) call set_cyclic(delqc(:,:,:,0),km,1)
      if (cyclic) call set_cyclic(delqc(:,:,:,1),km,1)

#endif

      if (my_pe == 0) then
         print*,''
         print*,' done'
         print*,''
      endif
#else
      if (my_pe == 0) then
         print*,''
         print*,' ERROR: isoneutral module disabled'
         print*,''
      endif
      call halt_stop(' in isoneutral_init')
#endif

      end subroutine isoneutral_init



      subroutine isoneutral

#ifndef disable_isoneutral_module
c
c=======================================================================
c
c     Compute isoneutral mixing coefficients.
c
c     input:
c       joff = offset relating row "j" in the MW to latitude "jrow"
c       js   = starting row within the MW for calculations
c       je   = ending row within the MW for calculations
c       is   = starting index longitude within the MW
c       ie   = ending index longitude within the MW
c
c     output:
c       Ai_ez = diffusivity*tapered slope
c               centered on east face of T cells
c       Ai_nz = diffusivity*tapered slope
c               centered on north face of T cells
c       Ai_bx = diffusivity*tapered slope
c               centered on bottom face of T cells
c       Ai_by = diffusivity*tapered slope
c               centered on bottom face of T cells
c            
c=======================================================================
c
c-----------------------------------------------------------------------
c     estimate drho and gradients on sides of T cells 
c-----------------------------------------------------------------------
c
      call isoneutral_elements
c
c-----------------------------------------------------------------------
c     compute Ai_ez and Bi_ez centered on eastern face of T cells
c-----------------------------------------------------------------------
c
      call isoneutral_ai_east 
c
c-----------------------------------------------------------------------
c     compute Ai_nz and Bi_nz centered on the northern face of T cells
c-----------------------------------------------------------------------
c
      call isoneutral_ai_north 
c
c-----------------------------------------------------------------------
c     evaluate Ai_bx, Ai_by, Bi_bx, and Bi_by
c     centered on bottom face of T cells
c-----------------------------------------------------------------------
c
      call isoneutral_ai_bottom
c
#endif
      end subroutine isoneutral



      subroutine isoneutral_elements 
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_module
      implicit none
      integer :: is,ie,js,je
      integer :: i,j,k,n,kp1
      real    :: tprime,sprime
      real T_i,T_j
      integer ip,jp
#ifdef partial_cell
      real dhte,dhtn,dhwe,dhwn
      dhte(i,k,j)   = min(dht(i+1,k,j),dht(i,k,j))
      dhtn(i,k,j)   = min(dht(i,k,j+1),dht(i,k,j))
      dhwe(i,k,j)   = min(dhwt(i+1,k,j),dhwt(i,k,j))
      dhwn(i,k,j)   = min(dhwt(i,k,j+1),dhwt(i,k,j))
      T_i(i,k,j,n,ip) = t(i+ip,max(1,k-1),j,n,taum1) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,taum1) - t(i+ip,k,j,n,taum1))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,taum1) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,taum1) - t(i,k,j+jp,n,taum1))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,taum1)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,taum1)
#endif
c
c=======================================================================
c     Estimate drho and normal gradients on faces of T cells
c=======================================================================
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     drho_dt and drho_ds at centers of T cells
c-----------------------------------------------------------------------
c
      do j=js-1,je+1
        do k=1,km
         call drhodt(t(is-1,k,j,1,taum1),t(is-1,k,j,2,taum1),
     &             drho(is-1,k,j,1),k,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
         call drhods(t(is-1,k,j,1,taum1),t(is-1,k,j,2,taum1),
     &             drho(is-1,k,j,2),k,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
        enddo
      enddo
      call set_cyclic(drho(:,:,:,1),km,1)
      call set_cyclic(drho(:,:,:,2),km,1)
c
c-----------------------------------------------------------------------
c     gradients at bottom face of T cells
c-----------------------------------------------------------------------
c
      do n=1,2
        do j=js-1,je+1
          do k=1,km
            kp1 = min(k+1,km)
            do i=is-1,ie+1
	      ddzt(i,k,j,n) = tmask(i,kp1,j)*
     &                        (t(i,k,j,n,taum1) - t(i,kp1,j,n,taum1))
#ifdef partial_cell
     &                      /dhwt(i,k,j)
#else
     &                      *dzwr(k)
#endif
            enddo
          enddo
          do i=is,ie
	    ddzt(i,0,j,n) = 0.
	  enddo
        enddo
        call set_cyclic(ddzt(:,:,:,n),km+1,1)
      enddo
c
c-----------------------------------------------------------------------
c     gradients at eastern face of T cells
c-----------------------------------------------------------------------
c
      do n=1,2
        do j=js-1,je+1
          do k=1,km
            do i=is-1,ie
	      ddxt(i,k,j,n) = tmask(i,k,j)*tmask(i+1,k,j)*cstr(j)*
     &              dxur(i)*(t_i(i,k,j,n,1) - t_i(i,k,j,n,0))
            enddo
          enddo
	enddo
        call set_cyclic(ddxt(:,:,:,n),km,1)
      enddo
c
c-----------------------------------------------------------------------
c     gradients at northern face of T cells
c-----------------------------------------------------------------------
c
      do n=1,2
        do j=js-1,je
          do k=1,km
            kp1 = min(k+1,km)
            do i=is-1,ie+1
	      ddyt(i,k,j,n) = tmask(i,k,j)*tmask(i,k,j+1)*dyur(j)*
     &                 (t_j(i,k,j,n,1) - t_j(i,k,j,n,0))  
            enddo
          enddo
	enddo
        call set_cyclic(ddyt(:,:,:,n),km,1)
      enddo
#endif
      end subroutine isoneutral_elements


      subroutine isoneutral_ai_east 
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_module
      implicit none
      integer is,ie,js,je
      integer i,j,k,kp2
      integer ip,kr
c
      real drodxe
      drodxe(i,k,j,ip) =    drho(i+ip,k,j,1)*ddxt(i,k,j,1) + 
     &                      drho(i+ip,k,j,2)*ddxt(i,k,j,2) 
      real drodze
      drodze(i,k,j,ip,kr) = drho(i+ip,k,j,1)*ddzt(i+ip,k-1+kr,j,1) + 
     &                      drho(i+ip,k,j,2)*ddzt(i+ip,k-1+kr,j,2)

      real Ai11,Ai13,Asteep,sumz,sxe,taper,fact1,fact2
c
c=======================================================================
c     Compute Ai_ez, Bi_ez, and K11 on center of east face of T cell.
c     Note re-scaling factor which reduces diffusivity where abs slope
c     "sxe"  exceeds the critical slope "sc" for the small slope approx.
c=======================================================================
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,je
        do k=1,km
          kp2 = min(km,k+2) 
          do i=is-1,ie
            Ai11 = (redi_gm(3)*ahisop+redi_gm(4)*agm(i,j))
            if(athkdf.eq.ahisop .and. athkdf.eq.ahsteep) then
    	      K11(i,k,j) = Ai11*tmask(i,k,j)*tmask(i+1,k,j)*fzisop(k)
            else
#ifdef disable_vertical_tapering_ahsteep
              Asteep = redi_gm(5)*ahsteep+redi_gm(6)*agm(i,j)
     &                  *fzisop(k)
#else
              Asteep = (redi_gm(5)*ahsteep+redi_gm(6)*agm(i,j))
     &                  *fzisop(k)
#endif
              Ai13   = (redi_gm(1)*Ai11-redi_gm(2)*agm(i,j))
     &                  *fzisop(k)
	      sumz = 0.
              do kr=0,1
                do ip=0,1
                  sxe    = -drodxe(i,k,j,ip)/(drodze(i,k,j,ip,kr)-epsln)
#ifdef dm_taper 
                  fact1 = (abs(sxe)-del_dm)/s_dm   
                  fact2 = sign(1.0,fact1)*min(19.0,abs(fact1)) 
                  taper = 0.5*(1.-tanh(fact2))   
#endif
#ifdef gkw_taper
                  if (abs(sxe) .gt. slmx) then
                    taper = (slmx/(abs(sxe) + epsln))**2
                  else
                    taper = 1.
                  endif
#endif
	          sumz = sumz +
#ifdef partial_cell
     &                 min(delqc(i,k,j,kr),delqc(i+1,k,j,kr))
#else          
     &                 dzw(k-1+kr)
#endif         
     &                 *tmask(i,k,j)*tmask(i+1,k,j)
     &                 *max(Asteep,Ai11*taper)
	          Ai_ez(i,k,j,ip,kr) = 
     &                 Ai13*taper*sxe*tmask(i,k,j)*tmask(i+1,k,j)
	        enddo
	      enddo
#ifdef partial_cell
              K11(i,k,j) = 0.25*sumz
#else      
              K11(i,k,j) = .5*dzt2r(k)*sumz
#endif     
            endif
          enddo
        enddo
c        call setbcx (Ai_ez(1,1,j,0,0), imt, km)
c        call setbcx (Ai_ez(1,1,j,1,0), imt, km)
c        call setbcx (Ai_ez(1,1,j,0,1), imt, km)
c        call setbcx (Ai_ez(1,1,j,1,1), imt, km)
c        call setbcx (K11(1,1,j), imt, km)
      enddo
c     
#endif
      end subroutine isoneutral_ai_east


      subroutine isoneutral_ai_north
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_module
      implicit none
      integer is,ie,js,je
      integer i,j,k,kp2
      integer kr,jq
      real drodyn
      drodyn(i,k,j,jq) =    drho(i,k,j+jq,1)*ddyt(i,k,j,1) + 
     &                      drho(i,k,j+jq,2)*ddyt(i,k,j,2) 
      real drodzn
      drodzn(i,k,j,jq,kr) = drho(i,k,j+jq,1)*ddzt(i,k-1+kr,j+jq,1) + 
     &                      drho(i,k,j+jq,2)*ddzt(i,k-1+kr,j+jq,2)

      real Ai22,Asteep,Ai23,sumz,taper,syn,fact1,fact2
c     
c=======================================================================
c     Compute "Ai_nz" and "Bi_nz" on center of north face of T cell.
c     Note re-scaling factor which reduces diffusivity where abs slope
c     "syn"  exceeds the critical slope "sc" for the small slope approx.
c=======================================================================
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-1,je

        do k=1,km
          kp2 = min(km,k+2) 
          do i=is,ie
            Ai22 = (redi_gm(3)*ahisop+redi_gm(4)*agm(i,j))
            if(athkdf.eq.ahisop .and. athkdf.eq.ahsteep) then
    	      K22(i,k,j) = Ai22*tmask(i,k,j)*tmask(i,k,j+1)*fzisop(k)
            else
#ifdef disable_vertical_tapering_ahsteep
              Asteep = redi_gm(5)*ahsteep+redi_gm(6)*agm(i,j)
     &                  *fzisop(k)
#else
              Asteep = (redi_gm(5)*ahsteep+redi_gm(6)*agm(i,j))
     &                  *fzisop(k)
#endif
              Ai23   = (redi_gm(1)*Ai22-redi_gm(2)*agm(i,j))
     &                  *fzisop(k)
  	      sumz = 0.
              do kr=0,1
                do jq=0,1
                  syn    = -drodyn(i,k,j,jq)/(drodzn(i,k,j,jq,kr)-epsln)
#ifdef dm_taper 
                  fact1 = (abs(syn)-del_dm)/s_dm   
                  fact2 = sign(1.0,fact1)*min(19.0,abs(fact1)) 
                  taper = .5*(1.-tanh(fact2))   
#endif
#ifdef gkw_taper
                  if (abs(syn) .gt. slmx) then
                    taper = (slmx/(abs(syn) + epsln))**2
                  else
                    taper = 1.
                  endif
#endif
  	          sumz = sumz +
#ifdef partial_cell
     &                 min(delqc(i,k,j,kr),delqc(i,k,j+1,kr))
#else           
     &                 dzw(k-1+kr)
#endif          
     &                 *tmask(i,k,j)*tmask(i,k,j+1)
     &                 *max(Asteep,Ai22*taper)
	          Ai_nz(i,k,j,jq,kr) = 
     &                 Ai23*taper*syn*tmask(i,k,j)*tmask(i,k,j+1)
	        enddo
	      enddo
#ifdef partial_cell
              K22(i,k,j) = 0.25*sumz
#else      
              K22(i,k,j) = .5*dzt2r(k)*sumz
#endif     
            endif
          enddo
        enddo
c        call setbcx (Ai_nz(1,1,j,0,0), imt, km)
c        call setbcx (Ai_nz(1,1,j,1,0), imt, km)
c        call setbcx (Ai_nz(1,1,j,0,1), imt, km)
c        call setbcx (Ai_nz(1,1,j,1,1), imt, km)
c        call setbcx (K22(1,1,j), imt, km)
      enddo
#endif     
      end subroutine isoneutral_ai_north


      subroutine isoneutral_ai_bottom 
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_module
      implicit none
      integer is,ie,js,je
      integer i,j,k,kp2,kp1
      integer kr,jq,ip
      real drodxb
      drodxb(i,k,j,ip,kr) = drho(i,k+kr,j,1)*ddxt(i-1+ip,k+kr,j,1) + 
     &                      drho(i,k+kr,j,2)*ddxt(i-1+ip,k+kr,j,2) 
      real drodyb
      drodyb(i,k,j,jq,kr) = drho(i,k+kr,j,1)*ddyt(i,k+kr,j-1+jq,1) + 
     &                      drho(i,k+kr,j,2)*ddyt(i,k+kr,j-1+jq,2) 
      real drodzb
      drodzb(i,k,j,kr) =    drho(i,k+kr,j,1)*ddzt(i,k,j,1) + 
     &                      drho(i,k+kr,j,2)*ddzt(i,k,j,2)
      real fzisopb
      real Ai31,Ai32,Ai33,sumx,sxb,taper,facty,syb,sumy
      real fact1,fact2
c
c=======================================================================
c     compute Ai_bx, Ai_by, Bi_bx, Bi_by, and K33 at the center
c     of the bottom face of T cells.
c=======================================================================
c
c-----------------------------------------------------------------------
c     compute "Ai_bx", "Ai_by", "Bi_bx", "Bi_by", & K33 on bottom face
c     of T cell. Note re-scaling factor to reduce mixing coefficient
c     "Ai" where slopes excede the critical slope "sc" for the small
c     slope approx, gm_skew, or biharmonic_rm.  
c-----------------------------------------------------------------------
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=1,km-1
	  kp1   = min(k+1,km)
	  kp2   = min(k+2,km)
	  fzisopb = .5*(fzisop(k+1)+fzisop(k))
          do i=is,ie
            Ai31  = fzisopb
     &       *(  redi_gm(1)*(redi_gm(3)*ahisop+redi_gm(4)*agm(i,j))
     &         + redi_gm(2)*agm(i,j))
            Ai32  = Ai31
            Ai33  = (redi_gm(3)*ahisop+redi_gm(4)*agm(i,j))*fzisopb
c           eastward slopes at the base of T cells
            sumx  = 0.
            do ip=0,1
              do kr=0,1
                sxb = -drodxb(i,k,j,ip,kr)/(drodzb(i,k,j,kr)-epsln)
#ifdef dm_taper
                fact1 = (abs(sxb)-del_dm)/s_dm 
                fact2 = sign(1.0,fact1)*min(19.0,abs(fact1))
                taper = .5*(1.-tanh(fact2))   
#endif
#ifdef gkw_taper
                if (abs(sxb) .gt. slmx) then
                  taper = (slmx/(abs(sxb) + epsln))**2
                else
                  taper = 1.
                endif
#endif
	        sumx = sumx
     &          + dxu(i-1+ip)*tmask(i,k+1,j)*Ai33*taper*sxb**2
#ifdef partial_cell
     &          *min(delqc(i-1+ip,k+kr,j,1-kr),delqc(i+ip,k+kr,j,1-kr))
#endif
	        Ai_bx(i,k,j,ip,kr) =
     &               Ai31*taper*sxb*tmask(i,k+1,j)
c
              enddo
            enddo
c           northward slopes at the base of T cells
            sumy  = 0.
            do jq=0,1
	      facty = csu(j-1+jq)*dyu(j-1+jq)
              do kr=0,1
                syb = -drodyb(i,k,j,jq,kr)/(drodzb(i,k,j,kr)-epsln)
#ifdef dm_taper
                fact1 = (abs(syb)-del_dm)/s_dm
                fact2 = sign(1.0,fact1)*min(19.0,abs(fact1))
                taper =  .5*(1.-tanh(fact2))  
#endif
#ifdef gkw_taper
                if (abs(syb) .gt. slmx) then
                  taper = (slmx/(abs(syb) + epsln))**2
                else
                  taper = 1.
                endif
#endif
	        sumy = sumy
     &          + facty*tmask(i,k+1,j)*Ai33*taper*syb**2
#ifdef partial_cell
     &          *min(delqc(i,k+kr,j-1+jq,1-kr),delqc(i,k+kr,j+jq,1-kr))
#endif
	        Ai_by(i,k,j,jq,kr) =
     &               Ai32*taper*syb*tmask(i,k+1,j)
              enddo
            enddo
            K33(i,k,j) = (dxt4r(i)*sumx + dyt4r(j)*cstr(j)*sumy)
#ifdef partial_cell
     &                   /dhwt(i,k,j)
#endif
          enddo
        enddo
c        call setbcx (Ai_bx(1,1,j,1,0), imt, km)
c        call setbcx (Ai_bx(1,1,j,0,0), imt, km)
c        call setbcx (Ai_bx(1,1,j,1,1), imt, km)
c        call setbcx (Ai_bx(1,1,j,0,1), imt, km)
c        call setbcx (Ai_by(1,1,j,1,0), imt, km)
c        call setbcx (Ai_by(1,1,j,0,0), imt, km)
c        call setbcx (Ai_by(1,1,j,1,1), imt, km)
c        call setbcx (Ai_by(1,1,j,0,1), imt, km)
c        call setbcx (K33(1,1,j), imt, km)
      enddo
#endif
      end subroutine isoneutral_ai_bottom 



      subroutine isoneutral_flux(n,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
c     K11(is_pe-1:ie_pe,km,js_pe:je_pe)
c     Ai_ez(is_pe-1:ie_pe,km,js_pe:je_pe,0:1,0:1)
c
c     K22(is_pe:ie_pe,km,js_pe-1:je_pe)
c     Ai_nz(is_pe:ie_pe,km,js_pe-1:je_pe,0:1,0:1)
c
      use spflame_module
      use isoneutral_module
      implicit none
      integer :: n
      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)
#ifndef disable_isoneutral_module
      integer is,ie,js,je
      integer i,j,k
      integer km1,km1kr,kpkr,jm1,kr,ip,jq,jp
      real sumz,sumx,sumy

      real T_i,T_j
#ifdef partial_cell
      real dhte,dhtn,dhwe,dhwn
      dhte(i,k,j)   = min(dht(i+1,k,j),dht(i,k,j))
      dhtn(i,k,j)   = min(dht(i,k,j+1),dht(i,k,j))
      dhwe(i,k,j)   = min(dhwt(i+1,k,j),dhwt(i,k,j))
      dhwn(i,k,j)   = min(dhwt(i,k,j+1),dhwt(i,k,j))
      T_i(i,k,j,n,ip) = t(i+ip,max(1,k-1),j,n,taum1) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,taum1) - t(i+ip,k,j,n,taum1))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,taum1) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,taum1) - t(i,k,j+jp,n,taum1))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,taum1)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,taum1)
#endif

c
c=======================================================================
c     isoneutral diffusive tracer fluxes are computed.
c=======================================================================
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     construct total isoneutral tracer flux at east face of "T" cells 
c-----------------------------------------------------------------------
c
      if (athkdf .ne. ahisop) then     
        do j=js,je
          do k=1,km
            km1 = max(k-1,1)
            do i=is-1,ie
	      sumz    = 0.
              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)
#ifdef partial_cell
     &                *min(delqc(i,k,j,kr),delqc(i+1,k,j,kr))
     &                /dhwt(i+ip,km1kr,j)
c      dz_wtr(i,k,j) = c1/dhwt(i,k,j)
#endif
     &             *(t(i+ip,km1kr,j,n,taum1)-t(i+ip,kpkr,j,n,taum1))
                enddo
	      enddo
              diff_fe(i,k,j) = diff_fe(i,k,j)+ 
#ifdef partial_cell
     &                                     0.25*sumz
#else
     &                                     sumz*.5*dzt2r(k)
#endif
     &         + cstr(j)*dxur(i)*(T_i(i,k,j,n,1) - T_i(i,k,j,n,0))
     &                          *K11(i,k,j)
            enddo
          enddo
c          call setbcx (diff_fe(1,1,j), imt, km)
        enddo
      else 
        do j=js,je
          do k=1,km
            do i=is-1,ie
              diff_fe(i,k,j) = diff_fe(i,k,j)+ 
     &              K11(i,k,j)*cstr(j)*dxur(i)
     &             *(T_i(i,k,j,n,1) - T_i(i,k,j,n,0))
#ifdef partial_cell
ce    I just try to debug here without knowing what I do
c     &             *dhte(i,k,j)
#endif
            enddo
          enddo
c          call setbcx (diff_fe(1,1,j), imt, km)
        enddo
      endif
c
c-----------------------------------------------------------------------
c     construct total isoneutral tracer flux at north face of "T" cells 
c-----------------------------------------------------------------------
c    
      if (athkdf .ne. ahisop) then     
        do j=js-1,je
          jm1 = max(1,j-1)
          do k=1,km
            do i=is,ie
	      sumz    = 0.
              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)
#ifdef partial_cell
     &                *min(delqc(i,k,j,kr),delqc(i,k,j+1,kr))
     &                /dhwt(i,km1kr,j+jq)
#endif
     &                 *(t(i,km1kr,j+jq,n,taum1)-t(i,kpkr,j+jq,n,taum1))
                enddo
	      enddo
             diff_fn(i,k,j) = diff_fn(i,k,j)+
#ifdef partial_cell
     &                 csu(j)*0.25*sumz
#else
     &                sumz*csu(j)*0.5*dzt2r(k)
#endif
     &        + (T_j(i,k,j,n,1)-T_j(i,k,j,n,0))*csu_dyur(j)
     &          *K22(i,k,j)
            enddo
          enddo
c          call setbcx (diff_fn(1,1,j), imt, km)
        enddo
      else 
        do j=js-1,je
          do k=1,km
            do i=is,ie
              diff_fn(i,k,j) = diff_fn(i,k,j)+
     &         K22(i,k,j)*csu_dyur(j)
     &        *(T_j(i,k,j,n,1)-T_j(i,k,j,n,0))
#ifdef partial_cell
ce    I just try to debug here without knowing what I do
c     &        *dhtn(i,k,j)
#endif
            enddo
          enddo
c          call setbcx (diff_fn(1,1,j), imt, km)
        enddo
      endif

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 implicitly. Note that there are some
c     cancellations of dxu(i-1+ip) and dyu(jrow-1+jq) 
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
	    sumx    = 0.
            do ip=0,1
              do kr=0,1
	        sumx = sumx
     &            + Ai_bx(i,k,j,ip,kr)*cstr(j)*
#ifdef partial_cell
     &          min(delqc(i-1+ip,k+kr,j,1-kr),delqc(i+ip,k+kr,j,1-kr))*
#endif
     &             (T_i(i-1+ip,k+kr,j,n,1) - T_i(i-1+ip,k+kr,j,n,0))
              enddo
	    enddo
	    sumy    = 0.
            do jq=0,1
              do kr=0,1
	        sumy    = sumy 
     &           + Ai_by(i,k,j,jq,kr)*csu(j-1+jq)*
#ifdef partial_cell
     &          min(delqc(i,k+kr,j-1+jq,1-kr),delqc(i,k+kr,j+jq,1-kr))*
#endif
     &            (T_j(i,k+kr,j-1+jq,n,1)-T_j(i,k+kr,j-1+jq,n,0))
              enddo
	    enddo
            diff_fb(i,k,j) = (dxt4r(i)*sumx +dyt4r(j)*cstr(j)*sumy)
#ifdef partial_cell
     &                          /dhwt(i,k,j)
#endif
          enddo
        enddo
c        call setbcx (diff_fb(1,0,j), imt, km+1)
      enddo
#endif
      end subroutine isoneutral_flux



      subroutine isoneutral_add_K33
#ifndef disable_isoneutral_module
      use spflame_module
      use isoneutral_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) = 
     &            (diff_cbt(i,k,j) + K33(i,k,j))*tmask(i,k+1,j)
          enddo
        enddo
      enddo
#endif
      end subroutine isoneutral_add_K33




      subroutine isoneutral_adv(adv_vetiso,adv_vntiso,adv_vbtiso)
c=======================================================================
c     compute GM eddy-induced transport velocities.
c     author:  R.C. Pacanowski e-mail => rcp@gfdl.gov 
c=======================================================================
      use spflame_module
      use isoneutral_module
      implicit none

      real :: adv_vetiso(is_pe-1:ie_pe,  km,js_pe  :je_pe) 
      real :: adv_vntiso(is_pe  :ie_pe,  km,js_pe-1:je_pe) 
      real :: adv_vbtiso(is_pe  :ie_pe,0:km,js_pe  :je_pe) 

#ifndef disable_isoneutral_module
      integer is,ie,js,je
      integer i,j,k,km1,kp1
      real ath0,at,bt,stn,ab,bb,sbn,ath_t,ath_b
      real ste,sbe,fact_t1,fact_b1,fact_t,fact_b
c
c-----------------------------------------------------------------------
c Arrays for GM.  Note, the advection velocities are still needed 
c with the default gm_skew approach, since generally wish to map 
c these velocities in snapshots.
c
c     adv_vetiso = zonal isopycnal mixing velocity computed at the 
c                  center of the eastern face of the "t" cells
c     adv_vntiso = meridional isopycnal mixing velocity computed at
c                  the center of the northern face of the "t" cells
c                  (Note: this includes the cosine as in "adv_vnt")
c     adv_vbtiso = vertical isopycnal mixing velocity computed at the
c                  center of the top face of the "t" cells
c     adv_fbiso  = "adv_vbtiso" * (tracer) evaluated at the center of
c                  the bottom face of the "t" cells
c-----------------------------------------------------------------------
c
      real top_bc(km), bot_bc(km)

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
      top_bc = 1.; bot_bc = 1.
      top_bc(1)  = 0.; bot_bc(km) = 0.
c
c-----------------------------------------------------------------------
c     compute the meridional component of the isoneutral mixing velocity
c     at the center of the northern face of the "T" cells.
c-----------------------------------------------------------------------
c
      do j=js-1,je
        do k=1,km
          km1 = max(k-1,1)
          kp1 = min(k+1,km)
          do i=is,ie
            Ath0 = agm(i,j)*fzisop(k)
            at =     (drho(i,k,j,1) + drho(i,k,j+1,1) + drho(i,km1,j,1)
     &              + drho(i,km1,j+1,1))
            bt =     (drho(i,k,j,2) + drho(i,k,j+1,2) + drho(i,km1,j,2)
     &              + drho(i,km1,j+1,2))
            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 =     (drho(i,k,j,1) + drho(i,k,j+1,1) + drho(i,kp1,j,1)
     &              + drho(i,kp1,j+1,1))
            bb =     (drho(i,k,j,2) + drho(i,k,j+1,2) + drho(i,kp1,j,2)
     &              + drho(i,kp1,j+1,2))
            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)
#ifdef dm_taper
            fact_t1 =  (abs(stn)-del_dm)/s_dm
            fact_b1 =  (abs(sbn)-del_dm)/s_dm
            fact_t =  sign(1.0,fact_t1)*min(19.0,abs(fact_t1))
            fact_b =  sign(1.0,fact_b1)*min(19.0,abs(fact_b1))
	    ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
     &             *.5*(1.-tanh(fact_t))
	    ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
     &             *.5*(1.-tanh(fact_b))
#endif
#ifdef gkw_taper
            if (abs(stn) .gt. slmx) then
	      ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
     &              *(slmx/(abs(stn) + epsln))**2
	    else
	      ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
	    endif
            if (abs(sbn) .gt. slmx) then
	      ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
     &              *(slmx/(abs(sbn) + epsln))**2
	    else
	      ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
	    endif
#endif
            adv_vntiso(i,k,j) = -(ath_t*stn*top_bc(k) -
     &                            ath_b*sbn*bot_bc(k))*dztr(k)*csu(j)
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute the zonal component of the isoneutral mixing velocity
c     at the center of the eastern face of the "T" grid cell.
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          km1 = max(k-1,1)
          kp1 = min(k+1,km)
          do i=is-1,ie
            Ath0 = agm(i,j)*fzisop(k)
            at =     (drho(i,k,j,1) + drho(i+1,k,j,1) + drho(i,km1,j,1)
     &              + drho(i+1,km1,j,1))
            bt =     (drho(i,k,j,2) + drho(i+1,k,j,2) + drho(i,km1,j,2)
     &              + drho(i+1,km1,j,2))
            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 =     (drho(i,k,j,1) + drho(i+1,k,j,1) + drho(i,kp1,j,1)
     &              + drho(i+1,kp1,j,1))
            bb =     (drho(i,k,j,2) + drho(i+1,k,j,2) + drho(i,kp1,j,2)
     &              + drho(i+1,kp1,j,2))
            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)

#ifdef dm_taper
            fact_t1 =  (abs(ste)-del_dm)/s_dm
            fact_b1 =  (abs(sbe)-del_dm)/s_dm
            fact_t =  sign(1.0,fact_t1)*min(19.0,abs(fact_t1))
            fact_b =  sign(1.0,fact_b1)*min(19.0,abs(fact_b1))
	    ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
     &             *.5*(1.-tanh(fact_t))
	    ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
     &             *.5*(1.-tanh(fact_b))
#endif
#ifdef gkw_taper
            if (abs(ste) .gt. slmx) then
	      ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
     &              *(slmx/(abs(ste) + epsln))**2
	    else
	      ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
	    endif
            if (abs(sbe) .gt. slmx) then
	      ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
     &              *(slmx/(abs(sbe) + epsln))**2
	    else
	      ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
	    endif
#endif
            adv_vetiso(i,k,j) = -(ath_t*ste*top_bc(k) - 
     &                            ath_b*sbe*bot_bc(k))*dztr(k)
          enddo
        enddo
      enddo
c
c-----------------------------------------------------------------------
c     compute the vertical component of the isoneutral mixing velocity
c     at the center of the bottom face of the "T" cells, using the
c     continuity equation for the GM velocities
c-----------------------------------------------------------------------
c
      do j=js,je
        do i=is,ie
          adv_vbtiso(i,0,j) = 0.
        enddo
      enddo
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(i,j),j) = 0.
        enddo
      enddo

#endif
      end subroutine isoneutral_adv 


