c  dyn_diff.f is a collection of fortran routines that will
c  allow to diffuse tracers along isopyncal surfaces and also
c  reduce the available potential energy by adding eddy induced
c  velocities to the eulerian mean velocities
c  first the density slope needs to be determined CALL SLOPE
c  then the isopycnal diffusion  can be computed  CALL DIFF_ISO
c  finally the eddy induced velocities are obtained CALL ADV_ISO
c
c  you need to add those to the Eulerian u,v,w BEFORE the advection
c  is called.
c
c  PARAMETER:
c    sigzmin : minimum vertical density gradient  (1e-6 [kg/m^4])
c    coef_diff_adv(z) : transfer coef for eddy induced velocity (1000 [m^2/s])
c    slmax : maximum slope for isopycnal diffusion (1e-2)
c    alpha  :  1 isopycnal diffusion, 0 horizontal diffusion 
c    coef_diff(z) : isopycnal diffusion coeff (1000 [m^2/s])
c    eps : ratio between isopycnal and diapycnal mixing (0)
c    
c   version 1.0 Aug 1996
c   Martin Visbeck
c
      subroutine diff_init(npt,nz,iglob,mgrid)

      common/grid/nxp,nyp,nxyc,nzl,nbx,nby,ncs,land,nlo,npbc

      include 'comm_data.h'
      include 'comm_diff.h'

      parameter (TORAD = 3.14159265/180., REARTH = 6378000.)

      do k = 1, npt
         j = (iox(k)-1)/nxp + 1
         i = iox(k) - (j-1)*nxp
         if (j.lt.nyp) then
            dyp(k) = TORAD * REARTH * (ym(j+1)-ym(j))
         else 
            dyp(k) = TORAD * REARTH * (ym(nyp)-ym(nyp-1))
         endif

         if (mgrid.eq.2) then
            csy(k) = cos(TORAD * ym(j))
            csyc(k) = cos(TORAD * ym(j) + dyp(k)/2./REARTH )
         else
            csy(k) = cos(TORAD * (ym(nyp)+ym(1))/2.)
            csyc(k)= csy(k)
         endif
         
         deg2met  = TORAD * REARTH * csy(k)

         
         if (i.lt.nxp) then
            dxp(k) = deg2met * (xm(i+1) - xm(i))
         else 
            dxp(k) = deg2met * (xm(nxp) - xm(nxp-1))
         endif
      enddo
      
      if (iglob.eq.1) then
         do k = 1, npt
            deg2met  = TORAD * REARTH * csy(k)
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            if (i.eq.nxp) then
               dxp(k)=deg2met*((xm(2)-xm(1))+(xm(nxp)-xm(nxp-1)))/2.
            endif
         enddo
      endif

      return
      end

c     ------------------------------------------------------------------
      subroutine diff_iso(alpha,coef_diff,npt,nz,nzi,h,tr,ftr,
     *       lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek,
     *       saxmk,saxpk,saymk,saypk)
c     ------------------------------------------------------------------

c     input:
c        h   - layer depths
c        tr  - tracer quantity (times h)
c       ftr  - forcing term to which diffusion will be added
c       slx,sly  - isopycnal slopes
c       coef_diff - diffusion coeficient
c       eps  - ratio of diapyncal / isopycnal
c       alpha - mixing slope (0-horizontal, 1-isopycnal)
c     output:
c       ftr  - forcing term to which diffusion will be added
c     diagnostics:
c       gtr  - diffusion tensor quantity  
c       trx,try,trz - tracer gradients
c       

c     subroutine that calculates diffusion tensor term in the
c         tracer equations
 
      include 'comm_para.h'
      include 'comm_diff.h'
      include 'diffiso.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +        ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension h(npt,nz),tr(npt,nz),ftr(npt,nz),tp(npt)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), nzi(npt),
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)

 
      eps     = diffiso_eps
      slmax   = diffiso_slmax
      n_order = diffiso_order

      if (alpha.eq.0.) then

         do k = 1, nz   
            npk = nptk(k)
            nxk = nbxk(k)
            nyk = nbyk(k)
            nck = ncsk(k)
            npbk = npbck(k)
            
            do j = 1, npk
               i = isk(j,k)
               tp(i) = tr(i,k)
            enddo

            do iorder = 1, n_order - 1
               call dfdx_a2c(tp,gtr,npt,npk,nxk,nck,lxxk(1,k),
     *              snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp)
               call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),
     *              snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *              saxmk(1,k),saxpk(1,k))
            enddo

            call dfdxh_a2c(tp,gtr,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp,h(1,k))
            call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
            
            do j = 1, npk
               i = isk(j,k)
               ftr(i,k) = ftr(i,k) + coef_diff*tp(i)
               tp(i) = tr(i,k)
            enddo
            
            do iorder = 1, n_order - 1
               call dfdy_a2c(tp,gtr,npt,npk,nyk,nck,lyyk(1,k),
     *              snyk(1,k),isyk(1,k),dyp,csyc)
               call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),
     *              snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
            enddo

            call dfdyh_a2c(tp,gtr,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),dyp,csyc,h(1,k))
            call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
            
            do j = 1, npk
               i = isk(j,k)
               ftr(i,k) = ftr(i,k) + coef_diff*tp(i)
            enddo
         enddo
         
      else
         
c.....compute vertical tracer gradient 
      call dfdz_ff(tr,trz,npt,nz,nzi,h)

c.....compute horizontal tracer gradient 
      do k = 1, nz  
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

c         call grad (nx,ny,npt,tr(1,k),dx2i,dy2i,ycos,inx,iny,isk(1,k),
c     *              trx(1,k),try(1,k),nxk,nyk,lxxk(1,k),lyyk(1,k))
         call dfdxk(tr(1,k),trx(1,k),npt,npk,0,nxk,nyk,nck,lxxk(1,k),
     *        lyxk(1,k),snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         call dfdyk(tr(1,k),try(1,k),npt,npk,0,nyk,nxk,nck,lyyk(1,k),
     *         lxyk(1,k),snyk(1,k),isyk(1,k),saymk(1,k),saypk(1,k))
      enddo

c do diffusion in x-direction
      do k = 1, nz   
c.....multiply gradients with slopes and diffusion factors  
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         call dfdx_a2c(tr(1,k),tp,npt,npk,nxk,nck,lxxk(1,k),
     *       snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp)
         do j = 1, npk - 1
            i = isk(j,k)
            ip= i + 1
            h_ave = (h(i,k)+h(ip,k))/2.
            ax= alpha*(slx(i,k)+slx(ip,k))/2.
            ay= alpha*(sly(i,k)+sly(ip,k))/2.
            ax2 = ax*ax
            ay2 = ay*ay
            axy = ax*ay
            sl = sqrt(ax2 + ay2)
            slfac =  max(0.,1. - sl/slmax * slred)
            eps1 = (1-eps)
            fac = h_ave*coef_diff/(1+ax2+ay2)
            tryc = (try(i,k)+try(ip,k))/2.
            trzc = (trz(i,k)+trz(ip,k))/2.
         
            gtr(i) = fac*(
     *           (1+eps*ax2+ay2)*tp(i) 
     *              + (-axy*eps1*tryc + ax*eps1*trzc)*slfac)

         enddo                     
         call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),snxk(1,k),npbk,
     *        lpbcwk(1,k),lpbcek(1,k),isk(1,k),saxmk(1,k),saxpk(1,k))

         do j = 1, npk
            i = isk(j,k)
            ftr(i,k) = ftr(i,k) + tp(i)
         enddo
      enddo

c do diffusion in y-direction
      do k = 1, nz   
c.....multiply gradients with slopes and diffusion factors  
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         call dfdy_a2c(tr(1,k),tp,npt,npk,nyk,nck,lyyk(1,k),
     *       snyk(1,k),isyk(1,k),dyp,csyc)
         do j = 1, npk - 1
            i = isyk(j,k)
            ip= isyk(j+1,k)
            h_ave = (h(i,k)+h(ip,k))/2.
            ax= alpha*(slx(i,k)+slx(ip,k))/2.
            ay= alpha*(sly(i,k)+sly(ip,k))/2.
            ax2 = ax*ax
            ay2 = ay*ay
            axy = ax*ay
            sl = sqrt(ax2 + ay2)
            slfac =  max(0.,1. - sl/slmax * slred)
            eps1 = csyc(i)*(1-eps)
            fac = h_ave*coef_diff/(1+ax2+ay2)
            trxc = (trx(i,k)+trx(ip,k))/2.
            trzc = (trz(i,k)+trz(ip,k))/2.
         
            gtr(i) = fac*(
     *             -axy*eps1*trxc*slfac + (1+ax2+eps*ay2)*tp(i) 
     *            + ay*eps1*trzc*slfac)
         enddo                     

         call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),snyk(1,k),isyk(1,k),
     *                 csy,saymk(1,k),saypk(1,k))

         do j = 1, npk
            i = isk(j,k)
            ftr(i,k) = ftr(i,k) + tp(i)
         enddo
      enddo

c do diffusion in z-direction
      call dfdz_a2c(tr,trz,npt,nz,nzi,h)

      do i = 1, npt
         nzb = nzi(i)
         do k = 1, nzb-1
            kp = k + 1
            ax= alpha*(slx(i,k)+slx(i,kp))/2.
            ay= alpha*(sly(i,k)+sly(i,kp))/2.
            ax2 = ax*ax
            ay2 = ay*ay
            axy = ax*ay
            sl = sqrt(ax2 + ay2)
            slfac =  max(0.,1. - sl/slmax * slred)
            eps1 = (1-eps)
            fac = coef_diff/(1+ax2+ay2) * slfac
            trxc = (trx(i,k)+trx(i,kp))/2.
            tryc = (try(i,k)+try(i,kp))/2.
            
            gtrz(i,k) = fac * ( 
     *              ax*eps1*trxc + ay*eps1*tryc + (eps+ax2+ay2)*trz(i,k))

          enddo
       enddo
 
      call dfdzh_c2a(gtrz,trz,npt,nz,nzi)

      do k = 1, nz    
         npk = nptk(k)
         do j = 1, npk
            i = isk(j,k)
            ftr(i,k) = ftr(i,k) + trz(i,k)
         enddo
      enddo
      endif

      return
      end

c     ------------------------------------------------------------------
      subroutine diff_hor(coef_diff,coef_dz,npt,nz,nzi,h,tr,ftr,
     *     lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek,
     *     saxmk,saxpk,saymk,saypk)
c     ------------------------------------------------------------------

c     input:
c        h   - layer depths
c        tr  - tracer quantity (times h)
c       ftr  - forcing term to which diffusion will be added
c       coef_diff - diffusion coeficient
c     output:
c       ftr  - forcing term to which diffusion will be added
c     diagnostics:
c       gtr  - diffusion tensor quantity  
c       trx,try,trz - tracer gradients
c       

c     subroutine that calculates diffusion tensor term in the
c         tracer equations
      
      include 'comm_para.h'
      include 'comm_diff.h'
      include 'diffiso.h'

      common /vert/  zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1),
     *               bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ),
     *               facz(MAXNZ), ttr_fac(MAXNZ)
      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +     ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension h(npt,nz),tr(npt,nz),ftr(npt,nz),tp(npt,2)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), nzi(npt),
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)

      
      n_order = diffiso_order

      do k = 1, nz   
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         coef = coef_diff * exp(-(zin(k)-zin(1))/coef_dz)
         
         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
            tp(i,2) = Kmap(i)*h(i,k)
         enddo

         do iorder = 1, n_order - 1
            call dfdy_a2c(tp,gtr,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),dyp,csyc)
            call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         enddo

         call dfdyh_a2c(tp,gtr,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),dyp,csyc,tp(1,2))
         call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            ftr(i,k) = ftr(i,k) + coef*tp(i,1)
            tp(i,1) = tr(i,k)
            tp(i,2) = Kaspect(i)*tp(i,2)
         enddo

         do iorder = 1, n_order - 1
            call dfdx_a2c(tp,gtr,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp)
            call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
         enddo

         call dfdxh_a2c(tp,gtr,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp,tp(1,2))
         call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            ftr(i,k) = ftr(i,k) + coef*tp(i,1)
         enddo
         
      enddo
      
      return
      end

c     ------------------------------------------------------------------
      subroutine diff_isogm(alpha,coef_diff,npt,nz,nzi,h,tr,ftr,
     *     lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek,
     *     saxmk,saxpk,saymk,saypk,ft_hd,ft_vd,flor,xmld)
c     ------------------------------------------------------------------

c     input:
c        h   - layer depths
c        tr  - tracer quantity (times h)
c       ftr  - forcing term to which diffusion will be added
c       slx,sly  - isopycnal slopes
c       coef_diff - diffusion coeficient
c       eps  - ratio of diapyncal / isopycnal
c       alpha - mixing slope (0-horizontal, 1-isopycnal)
c       Kmap - the variable weighting 
c     output:
c       ftr  - forcing term to which diffusion will be added
c     diagnostics:
c       gtr  - diffusion tensor quantity  
c       trx,try,trz - tracer gradients
c       

c     subroutine that calculates diffusion tensor term in the
c         tracer equations

c     assumes small slope approximation and diffiso_cadv = coef_diff
      
      include 'comm_para.h'
      include 'comm_diff.h'
      include 'diffiso.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +     ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension h(npt,nz),tr(npt,nz),ftr(npt,nz),tp(npt,2)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), nzi(npt),
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)
     *     ,ft_hd(npt,nz) ,ft_vd(npt,nz), xmld(npt)

      
      slmax = diffiso_slmax
      eps   = diffiso_eps
      eps1 = 1 - eps
      f0 = f0 + tscl*(1.-f0)

c compute horizontal diffusion 

      do k = 1, nz   
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         do j = 1, npk
            i = isk(j,k)
            tp(i,2) = h(i,k)*Kmap(i)
         enddo

         call dfdyh_a2c(tr(1,k),gtr,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),dyp,csyc,tp(1,2))
         call dfdy_c2a(gtr,tp,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            diff = f0*coef_diff*tp(i,1)
            ftr(i,k) = ftr(i,k) + diff
            ft_hd(i,k) = ft_hd(i,k) + flor*diff
            tp(i,2) = tp(i,2)*Kaspect(i)
         enddo
         
         call dfdxh_a2c(tr(1,k),gtr,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),dxp,tp(1,2))
         call dfdx_c2a(gtr,tp,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            diff =  f0*coef_diff*tp(i,1)
            ftr(i,k) = ftr(i,k) + diff
            ft_hd(i,k) = ft_hd(i,k) + flor*diff
         enddo
         
      enddo
      

      if (alpha.ne.0.) then
c compute vertical diffusion

         do k = 1, nz   
            npk = nptk(k)
            nxk = nbxk(k)
            nyk = nbyk(k)
            nck = ncsk(k)
            npbk = npbck(k)
            
            call dfdxk(tr(1,k),trx(1,k),npt,npk,0,nxk,nyk,nck,lxxk(1,k),
     *           lyxk(1,k),snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
            call dfdyk(tr(1,k),try(1,k),npt,npk,0,nyk,nxk,nck,lyyk(1,k),
     *           lxyk(1,k),snyk(1,k),isyk(1,k),saymk(1,k),saypk(1,k))
         enddo
         
         call dfdz_a2c(tr,trz,npt,nz,nzi,h)
         
         do i = 1, npt
            nzb = nzi(i)
            do k = 1, nzb-1
               kp = k + 1
               ax= alpha*(slx(i,k)+slx(i,kp))/2.
               ay= alpha*(sly(i,k)+sly(i,kp))/2.
               ax2 = ax*ax
               ay2 = ay*ay
               axy = ax*ay
               sl = sqrt(ax2 + ay2)
               slfac =  max(0.,1. - sl/slmax * slred)
               fac = f0*coef_diff*Kmap(i)/(1+ax2+ay2) * slfac
               facml = min(float(k),xmld(i))/xmld(i)
               facx = Kaspect(i)
               trxc = (trx(i,k)+trx(i,kp))/2.
               tryc = (try(i,k)+try(i,kp))/2.
               
               gtrz(i,k) = facml* fac * ( 
     *       2.*eps1*(facx*ax*trxc + ay*tryc) + (eps + ax2+ay2)*trz(i,k))
               
            enddo
         enddo
         
         call dfdzh_c2a(gtrz,trz,npt,nz,nzi)
         
         do k = 1, nz    
            npk = nptk(k)
            do j = 1, npk
               i = isk(j,k)
               diff = trz(i,k)
               ftr(i,k) = ftr(i,k) + diff
               ft_vd(i,k) = ft_vd(i,k) + flor*diff
            enddo
         enddo
      endif
      
      return
      end
      
c     ------------------------------------------------------------------
      subroutine slope(npt,sig,sigl,nz,nzi,h,lxxk,lyyk,lxyk,lyxk,snxk,snyk,
     *       isyk,isk,lok,tp,lpbcwk,lpbcek,saxmk,saxpk,saymk,saypk,xmld)
c     ------------------------------------------------------------------

c     input:
c       sig   - potential density, reference to surface
c       sigl  - locally referenced potential density
c       slmax - maximum allowed slope for slx,sly
c     output:
c       sigx,sigy - horizontal density gradients (central)
c       sigz    - vertical density gradient (plain)
c       slx  - zonal slope of density surfaces (A-grid)
c       sly  - meridional slope of density surfaces (A-grid)

c     subroutine that calculates slopes of density surfaces
 
      include 'comm_para.h'
      include 'diffiso.h'
      include 'comm_diff.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     *    ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension sig(npt,nz),nzi(nz),h(npt,nz),sigl(npt,nz)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz),xmld(npt),
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)
 

      slmax = diffiso_slmax

c..... compute vertical density gradients, centered, for use here only
      call dfdz_ff(sig,sigz,npt,nz,nzi,h)

      do i=1,npt
        xmld(i) = 1.
      enddo
      do k = 1, nz  
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)
         nok = nlok(k)

c.....compute horizontal gradients
         call dfdxk(sigl(1,k),sigx(1,k),npt,npk,0,nxk,nyk,nck,lxxk(1,k),
     *        lyxk(1,k),snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         call dfdyk(sigl(1,k),sigy(1,k),npt,npk,0,nyk,nxk,nck,lyyk(1,k),
     *         lxyk(1,k),snyk(1,k),isyk(1,k),saymk(1,k),saypk(1,k))


c..... compute slopes
c     when we go to the GM as an antisymmetric diffusion tensor, then
c        we will not need sigx or sigy or sigz to be kept around

         do j = 1, npk
            i = isk(j,k)
            sigz(i,k) = min(sigz(i,k),0.)
            slx(i,k) = - sigx(i,k)/
     *                 min(sigz(i,k),sigzmin-abs(sigx(i,k)/slmax))
            sigx(i,k) = - slx(i,k)*(sigz(i,k)+sigzmin)
            sly(i,k) = - sigy(i,k)/
     *                 min(sigz(i,k),sigzmin-abs(sigy(i,k)/slmax))
            sigy(i,k) = - sly(i,k)*(sigz(i,k)+sigzmin)
            dsig = sig(i,k) - sig(i,1)

         enddo

      enddo


      return
      end

c     ---------------------------------------------------------------------
      subroutine dfdx_a2c (f,df,npt,npk,nbx,ncs,lxx,snx
     *           ,npbc,lpbcw,lpbce,isk,dxp)
c     ---------------------------------------------------------------------
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk),dxp(npt)

      do i = 1, npk - 1
         j = isk(i)
         df(j) = (f(j+1)-f(j))/dxp(j)
      enddo
c....................... periodic B.C.
      do i = 1, npbc
         ie = lpbce(i)
         iw = lpbcw(i)
         df(ie)   = (f(iw) - f(ie))/dxp(ie)
      enddo
      
      return
      end

c     ---------------------------------------------------------------------
      subroutine dfdxh_a2c (f,df,npt,npk,nbx,ncs,lxx,snx
     *           ,npbc,lpbcw,lpbce,isk,dxp,h)
c     ---------------------------------------------------------------------
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk),dxp(npt),h(npt)

      do i = 1, npk - 1
         j = isk(i)
         df(j) = 0.5*(h(j+1)+h(j))*(f(j+1)-f(j))/dxp(j)
      enddo
c....................... periodic B.C.
      do i = 1, npbc
         ie = lpbce(i)
         iw = lpbcw(i)
         df(ie)   = 0.5*(h(iw)+h(ie))*(f(iw) - f(ie))/dxp(ie)
      enddo
c      do i = 1, nbx
c         i1 = lxx(i,1)
c         i2 = lxx(i,2)
c         if (snx(i).gt.0) df(i1) = 0.
c         if (snx(i).lt.0) df(i2) = 0.
c      enddo
      
      return
      end

c     ---------------------------------------------------------------------
      subroutine dfdx_c2a (f,df,npt,npk,nbx,ncs,lxx,snx
     *           ,npbc,lpbcw,lpbce,isk,saxm,saxp)
c     ---------------------------------------------------------------------
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk),saxm(npt),saxp(npt)

      do i = 2, npk
         j = isk(i)
         df(j) = saxp(j)*f(j) - saxm(j)*f(j-1)
      enddo
c....................... periodic B.C.
      do i = 1, npbc
         ie = lpbce(i)
         iw = lpbcw(i)
         df(iw)   = saxp(iw)*f(iw) - saxm(iw)*f(ie) 
      enddo
      
      do i = 1, nbx
         i1 = lxx(i,1)
         i2 = lxx(i,2)
         f1 = saxp(i1)*f(i1)
         if (snx(i).lt.0) f1 = -saxm(i1)*f(i2)
         df(i1) = f1
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine sdfdx_a2c (f,df,npt,npk,nbx,ncs,lxx,snx
     *           ,npbc,lpbcw,lpbce,isk,sxp)
c     ---------------------------------------------------------------------
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4),sxp(npt)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk)

c      -0.25*del_x^2(df/dx) from the a-grid to the c-grid

      do i = 1, npk - 1
         j = isk(i)
         df(j) = (f(j)-f(j+1))*sxp(j)
      enddo
c....................... periodic B.C.
      do i = 1, npbc
         ie = lpbce(i)
         iw = lpbcw(i)
         df(ie)  = (f(ie) - f(iw))*sxp(ie)
      enddo
      
      return
      end

c     ---------------------------------------------------------------------
      subroutine sdfdxh_a2c (f,df,npt,npk,nbx,ncs,lxx,snx
     *           ,npbc,lpbcw,lpbce,isk,h,sxp)
c     ---------------------------------------------------------------------
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4),sxp(npt),h(npt)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk)

c      -0.25*h*del_x^2(df/dx) from the a-grid to the c-grid

      do i = 1, npk - 1
         j = isk(i)
         df(j) = (f(j)-f(j+1))*sxp(j)*0.5*(h(j+1)+h(j))
      enddo
c....................... periodic B.C.
      do i = 1, npbc
         ie = lpbce(i)
         iw = lpbcw(i)
         df(ie)   = (f(ie) - f(iw))*sxp(ie)*0.5*(h(iw)+h(ie))
      enddo
      
      return
      end

c     ------------------------------------------------------------------
      subroutine dfdy_a2c(f,df,npt,npk,nby,ncs,lyy,sny,isy,dyp,csyc)
c     ------------------------------------------------------------------
c      dfdy in flux form
      implicit real(a-h,o-z),integer(i-n)
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4)
     *          ,isy(npk),dyp(npt),csyc(npt)
 
      do i = 1, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         df(j)=(f(jp)-f(j))/dyp(j)* csyc(j)
      enddo
      
      return
      end

c     ------------------------------------------------------------------
      subroutine dfdyh_a2c(f,df,npt,npk,nby,ncs,lyy,sny,isy,dyp,csyc,h)
c     ------------------------------------------------------------------
c      dfdy in flux form
      implicit real(a-h,o-z),integer(i-n)
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4)
     *          ,isy(npk),dyp(npt),h(npt),csyc(npt)
 
      do i = 1, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         df(j)=0.5*(h(jp)+h(j))*(f(jp)-f(j))/dyp(j)* csyc(j)
      enddo
c      do i = 1, nby
c         i1 = lyy(i,1)
c         i2 = lyy(i,2)
c         if (sny(i).gt.0) df(i1) = 0.
c         if (sny(i).lt.0) df(i2) = 0.
c      enddo
      
      return
      end


c     ------------------------------------------------------------------
      subroutine dfdy_c2a(f,df,npt,npk,nby,ncs,lyy,sny,isy,csy,saym,sayp)
c     ------------------------------------------------------------------
c      dfdy in flux form, from c-grid to a-grid
      implicit real(a-h,o-z),integer(i-n)
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4)
     *          ,isy(npk),saym(npt),sayp(npt),csy(npt)
 
      do i = 2, npk
         j   = isy(i)
         jm  = isy(i-1)
         df(j)=(sayp(j)*f(j)-saym(j)*f(jm))/csy(j)
      enddo
      
      do i = 1, nby
         i1 = lyy(i,1)
         f1 = sayp(i1)*f(i1)
         i2 = lyy(i,2)
         if (sny(i).lt.0) f1 = -saym(i1)*f(i2)
         df(i1) = f1/csy(i1)
      enddo

      return
      end

c     ------------------------------------------------------------------
      subroutine sdfdy_a2c(f,df,npt,npk,nby,ncs,lyy,sny,isy,csyc,syp)
c     ------------------------------------------------------------------
c      dfdy in flux form
      implicit real(a-h,o-z),integer(i-n)
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4),syp(npt)
     *          ,isy(npk),csyc(npt)

c      -0.25*del_y^2(d(f)/dy) from the a-grid to the c-grid
 
      do i = 1, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         df(j)=(f(j)-f(jp))*syp(j)*csyc(j)
      enddo
      
      return
      end

c     ------------------------------------------------------------------
      subroutine sdfdyh_a2c(f,df,npt,npk,nby,ncs,lyy,sny,isy,csyc,h,syp)
c     ------------------------------------------------------------------
c      dfdy in flux form
      implicit real(a-h,o-z),integer(i-n)
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4),h(npt),syp(npt)
     *          ,isy(npk),csyc(npt)

c      -0.25*h*del_y^2(df/dy) from the a-grid to the c-grid
 
      do i = 1, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         df(j)=(f(j)-f(jp))*syp(j)*0.5*(h(jp)+h(j))*csyc(j)
      enddo
      
      return
      end

      subroutine dfdz_ff(f,df,npt,nz,nzi,h)
c-------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension h(npt,nz), f(npt,nz), df(npt,nz)
      dimension nzi(npt)

c.....    nzi  : number of layers
c.....    h    : layer depth
c.....    f    : quantity to be differentiated
c.....    df   : (central difference) vertical derivative of tracer 
c
      do i = 1, npt
         nzb = nzi(i)
         do k=  2, nzb -1
            kp = k + 1
            km = k - 1
            fp = f(i,kp)   
            fm = f(i,km)   
            df(i,k) = (fm - fp) / (2.*h(i,k))
         enddo
         df(i,1  ) = (f(i,1) - f(i,2)) / h(i,1)
         df(i,nzb) = (f(i,nzb-1) - f(i,nzb)) / h(i,nzb)
      enddo
      
      return
      end

      subroutine dfdz_a2c(f,df,npt,nz,nzi,h)
c-------------------------------------------------------------
c  a-grid means layer centers
c  c-grid means layer interfaces
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension h(npt,nz), f(npt,nz), df(npt,nz)
      dimension nzi(npt)

c.....    nzi  : number of layers
c.....    h    : layer depth
c.....    f    : quantity to be differentiated
c.....    df   : vertical derivative of tracer 
c
      do i = 1, npt
         nzb = nzi(i)
         dz = h(i,1)
         do k=  1, nzb -1
            kp = k + 1
            dz = 2.*h(i,k) - dz
            df(i,k) = (f(i,k) - f(i,kp)) / dz
         enddo
      enddo
      
      return
      end


      subroutine dfdz_c2a(f,df,npt,nz,nzi,h)
c-------------------------------------------------------------
c  a-grid means layer centers
c  c-grid means layer interfaces
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension h(npt,nz), f(npt,nz), df(npt,nz)
      dimension nzi(npt)

c.....    nzi  : number of layers
c.....    h    : layer depth
c.....    f    : quantity to be differentiated
c.....    df   : vertical derivative of tracer 
c
      do i = 1, npt
         nzb = nzi(i)
         do k=  2, nzb -1
            df(i,k) = (f(i,k-1) - f(i,k)) / h(i,k)
         enddo
         df(i,1)   = (0. - f(i,1)) / h(i,1)
         df(i,nzb) = (f(i,nzb-1) - 0.) / h(i,nzb)
      enddo
      
      return
      end

      subroutine dfdzh_c2a(f,df,npt,nz,nzi)
c-------------------------------------------------------------
c  a-grid means layer centers
c  c-grid means layer interfaces
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension f(npt,nz), df(npt,nz)
      dimension nzi(npt)

c.....    nzi  : number of layers
c.....    f    : quantity to be differentiated
c.....    df   : vertical derivative of tracer 
c
      do i = 1, npt
         nzb = nzi(i)
         do k=  2, nzb -1
            df(i,k) = f(i,k-1) - f(i,k)
         enddo
         df(i,1)   = - f(i,1)
         df(i,nzb) = f(i,nzb-1)
      enddo
      
      return
      end

c     ------------------------------------------------------------------
      subroutine cons_shap_scl(nstep,npt,nz,nz1,nz2,h,tr,tp,
     *       lxxk,lyyk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,
     *       saxmk,saxpk,saymk,saypk,sxp,syp)
c     ------------------------------------------------------------------
      include 'comm_para.h'
      include 'comm_diff.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +        ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension h(npt,nz),tr(npt,nz),tp(npt,4),sxp(npt,nz),syp(npt,nz)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), 
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)

      common /new_shap/nordu,nshapu,mshapu,dshapu, nordh,nshaph,mshaph,dshaph
     *               , nordp,nshapp,mshapp,dshapp
     *               , nordd,nshapd,mshapd,dshapd,dshapv

      if (mshaph .ge. 0) return

      if (nstep.eq.0 .or. (nshaph.ne.0 .and. mod(nstep,nshaph).eq.0))  then
 
      do k = nz1, nz2
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordh - 1
            call sdfdx_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),sxp(1,k))
            call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
         enddo
         
         call sdfdxh_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),h(1,k),sxp(1,k))
         call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshaph* tp(i,1)/ h(i,k)
         enddo
         
         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordh - 1
            call sdfdy_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csyc,syp(1,k))
            call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         enddo

         call sdfdyh_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csyc,h(1,k),syp(1,k))
         call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshaph* tp(i,1)/ h(i,k)
         enddo
      enddo

      endif
      
      return
      end
      
c     ------------------------------------------------------------------
      subroutine cons_shap_psi(nstep,npt,nz,tr,tp,
     *       lxxk,lyyk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,
     *       saxmk,saxpk,saymk,saypk,sxp,syp)
c     ------------------------------------------------------------------
      include 'comm_para.h'
      include 'comm_diff.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +        ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension tr(npt,nz),tp(npt,4),sxp(npt,nz),syp(npt,nz)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), 
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)

      common /new_shap/nordu,nshapu,mshapu,dshapu, nordh,nshaph,mshaph,dshaph
     *               , nordp,nshapp,mshapp,dshapp
     *               , nordd,nshapd,mshapd,dshapd,dshapv

      if (mshapp .ge. 0) return

      if (nstep.eq.0 .or. (nshapp.ne.0 .and. mod(nstep,nshapp).eq.0))  then
 
         k = 1
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordp - 1
            call sdfdx_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),sxp(1,k))
            call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
         enddo
         
         call sdfdx_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),sxp(1,k))
         call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshapp* tp(i,1)
         enddo
         
         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordp - 1
            call sdfdy_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csyc,syp(1,k))
            call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         enddo

         call sdfdy_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csyc,syp(1,k))
         call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshapp* tp(i,1)
         enddo

      endif
      
      return
      end
      
c     ------------------------------------------------------------------
      subroutine cons_shap_dens(nstep,npt,nz,nz1,nz2,h,tr,tp,
     *       lxxk,lyyk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,
     *       saxmk,saxpk,saymk,saypk,sxp,syp)
c     ------------------------------------------------------------------
      include 'comm_para.h'
      include 'comm_diff.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +        ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension h(npt,nz),tr(npt,nz),tp(npt,4),sxp(npt,nz),syp(npt,nz)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), 
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)

      common /new_shap/nordu,nshapu,mshapu,dshapu, nordh,nshaph,mshaph,dshaph
     *               , nordp,nshapp,mshapp,dshapp
     *               , nordd,nshapd,mshapd,dshapd,dshapv

      if (mshapd .ge. 0) return

      if (nstep.eq.0 .or. (nshapd.ne.0 .and. mod(nstep,nshapd).eq.0))  then
 
      do k = nz1, nz2
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordd - 1
            call sdfdx_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),sxp(1,k))
            call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *           snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *           saxmk(1,k),saxpk(1,k))
         enddo
         
         call sdfdxh_a2c(tp,tp(1,2),npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),h(1,k),sxp(1,k))
         call dfdx_c2a(tp(1,2),tp,npt,npk,nxk,nck,lxxk(1,k),
     *        snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *        saxmk(1,k),saxpk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshapd* tp(i,1)/ h(i,k)
         enddo
         
         do j = 1, npk
            i = isk(j,k)
            tp(i,1) = tr(i,k)
         enddo

         do ns = 1, nordd - 1
            call sdfdy_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csyc,syp(1,k))
            call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *           snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         enddo

         call sdfdyh_a2c(tp,tp(1,2),npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csyc,h(1,k),syp(1,k))
         call dfdy_c2a(tp(1,2),tp,npt,npk,nyk,nck,lyyk(1,k),
     *        snyk(1,k),isyk(1,k),csy,saymk(1,k),saypk(1,k))
         
         do j = 1, npk
            i = isk(j,k)
            tr(i,k) = tr(i,k) - dshapd* tp(i,1)/ h(i,k)
         enddo
      enddo

      endif
      
      return
      end
      
c     ------------------------------------------------------------------
      subroutine diff_div(coef_diff,npt,nz,nzi,fhd,fu,fv,
     *     lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek,
     *     saxmk,saxpk,saymk,saypk,nx,ny,iox,mask,isteer)
c     ------------------------------------------------------------------
c     subroutine to add divergence damping to momentum equations

c     input:
c        fhd      - layer divergence
c        fu,fv    - momentum forcing terms
c       coef_diff - diffusion coeficient
c     output:
c        fu,fv    - momentum forcing terms + divergence damping
c     ------------------------------------------------------------------

      include 'comm_para.h'
      include 'comm_diff.h'
      include 'diffiso.h'

      common /vert/  zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1),
     *               bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ),
     *               facz(MAXNZ), ttr_fac(MAXNZ)
      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +     ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension fhd(npt,nz), fu(npt,nz),tp(npt,2)
      dimension fv(npt,nz)
      dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz),
     *     snxk(MAXNB,nz),snyk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *     lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz), nzi(npt),
     *     saxmk(npt,nz), saxpk(npt,nz), saymk(npt,nz), saypk(npt,nz)
      dimension iox(npt), mask(nx*ny,nz)

      do k = 1, nz   
         npk = nptk(k)
         nxk = nbxk(k)
         nyk = nbyk(k)
         nck = ncsk(k)
         npbk = npbck(k)

         call dfdxk(fhd(1,k),tp(1,1),npt,npk,0,nxk,nyk,nck,lxxk(1,k),
     *         lyxk(1,k),snxk(1,k),npbk,lpbcwk(1,k),lpbcek(1,k),isk(1,k),
     *         saxmk(1,k),saxpk(1,k))
         call dfdyk(fhd(1,k),tp(1,2),npt,npk,0,nyk,nxk,nck,lyyk(1,k),
     *         lxyk(1,k),snyk(1,k),isyk(1,k),
     *         saymk(1,k),saypk(1,k))
         
         if (isteer.eq.1) then
            nxy = nx*ny
            do i = 1, nbxk(k)
               ik = lxxk(i,k)
               gridfac = saypk(ik,k)/saxpk(ik,k)
               ixy = iox(ik)
               i_N = 0
               i_S = 0
               if (ixy+nx.le.nxy) i_N = mask(ixy + nx,k)
               if (ixy-nx.ge.1)   i_S = mask(ixy - nx,k)
               if (i_N.ne.0) tp(i_N,2) = tp(i_N,2) - 0.5* tp(ik,1)* gridfac
               if (i_S.ne.0) tp(i_S,2) = tp(i_S,2) + 0.5* tp(ik,1)* gridfac
            enddo
            do i = 1, nbyk(k)
               ik = lyyk(i,k)
               gridfac = saxpk(ik,k)/saypk(ik,k)
               ixy = iox(ik)
               i_E = 0
               i_W = 0
               if (ixy+1.le.nxy) i_E = mask(ixy + 1,k)
               if (ixy-1.ge.1)   i_W = mask(ixy - 1,k)
               if (i_E.ne.0) tp(i_E,1) = tp(i_E,1) - 0.5* tp(ik,2)* gridfac
               if (i_W.ne.0) tp(i_W,1) = tp(i_W,1) + 0.5* tp(ik,2)* gridfac
            enddo
            do i = 1, nbxk(k)
               ik = lxxk(i,k)
               tp(ik,1) = 0.
            enddo
            do i = 1, nbyk(k)
               ik = lyyk(i,k)
               tp(ik,2) = 0.
            enddo
         endif

         do j = 1, npk
            i = isk(j,k)
            fu(i,k) = fu(i,k) + coef_diff*tp(i,1)
            fv(i,k) = fv(i,k) + coef_diff*tp(i,2)
         enddo

      enddo
      
      return
      end