      subroutine new_topo (nxp,nyp,nz,npt,zin,dzin,hin,nsig,sigma,dept,h,nptk,nzi)
      dimension zin(nz+1),dzin(nz+1),dept(npt),nptk(nz),h(npt,nz)
     *         ,nzi(npt),hin(nz),sigma(nz)
      include 'comm_new.h'
      
      do k = 1, nz
         do i = 1, npt
c            h(i,k) = -98765432.
            h(i,k) = 0.
         enddo
      enddo
      
      do i = 1, npt
         dep = dept(i)
         n = 1
         do k = 2, nz
            if (dep.lt.(zin(k)+zin(k-1))/2. + depth_fac*hin(k)) goto 10
            n = k
         enddo
   10    nzi(i) = n
      enddo
      
      if (initbt.eq.0) then
         do i = 1, npt
            mz = nzi(i)
            do k = 1, mz
               h(i,k) = dzin(k+1) + dzin(k)
            enddo
         enddo
      elseif (initbt.eq.3) then
         mmz = 0
         do i = 1, npt
            mmz = max(nzi(i),mmz)
         enddo
         do i = 1, npt
            mz = nzi(i)
            do k = 1, mz - 1
               h(i,k) = dzin(k+1) + dzin(k)
            enddo
            if (mz.eq.mmz) then
               h(i,mz) = dzin(mz) + (dept(i)-zin(mz))
            else
               h(i,mz) = dzin(mz+1) + dzin(mz)
            endif
         enddo
      else
         do i = 1, npt
            dep = dept(i)
            mz = nzi(i)
            do k = 1, mz - 1
               h(i,k) = dzin(k+1) + dzin(k)
            enddo
            dzp = dep - zin(mz)
            if (initbt.eq.1) dzp = min(dzin(mz+1),dzp)
            h(i,mz) = dzin(mz) + dzp
         enddo
      endif
      
      return
      end

c     ---------------------------------------------------------------------
      subroutine dfdxk (f,df,npt,npk,nbu,nbx,nby,ncs,lxx,lyx,snx
     *           ,npbc,lpbcw,lpbce,isk,saxm,saxp)
c     ---------------------------------------------------------------------
c............ a version with Periodic B.C. (Senya Basin)
      implicit real(a-h,o-z),integer(i-n)
      logical use_hi
      common /order/ use_hi
      dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4), lyx(nby)
      dimension lpbcw(npbc), lpbce(npbc),isk(npk),saxm(npt),saxp(npt)
      
c     compute centered first derivative in flux form for entire grid.
      do i = 2, npk - 1
         j = isk(i)
c            df(j) = 0.5* ( saxp(j)*(f(j+1)+f(j)) - saxm(j)*(f(j)+f(j-1)) )
         df(j) = 0.5* saxp(j)*(f(j+1)-f(j-1))
      enddo
      
c....................... periodic B.C.
      do i = 1, npbc
         i2 = lpbce(i)
         f4 = f(i2)
         f3 = f(i2-1)
         i1 = lpbcw(i)
         f5 = f(i1)
         f6 = f(i1+1)
c            df(i1) = 0.5* (saxp(i1)*(f6+f5) - saxm(i1)*(f5+f4))
c            df(i2) = 0.5* (saxp(i2)*(f4+f5) - saxm(i2)*(f3+f4))
         df(i1) = 0.5* saxp(i1)*(f6 - f4)
         df(i2) = 0.5* saxp(i2)*(f5 - f3)
      enddo
      
      if (nbu .eq. -1 .or. nbu .eq. 2) then
         do i = 1, nbx
            i1 = lxx(i,1)
            df(i1) = 0.            
         enddo
      else
c  if mbox = 1, this is linear extrapolation to boundary point 
c                         (N.B. f1 is the midcell value)
         do i = 1, nbx
            i1 = lxx(i,1)
            i2 = lxx(i,2)
            f1 = f(i1)
            f2 = (f1 + f(i2))/2.
            df(i1) = snx(i)*saxp(i1)*(f2 - f1)
         enddo
      endif
      
      
      if (nbu .eq. -1 .or. nbu.eq.1) then
c..................set the derivative along the boundary equal to zero.
         do i = 1, nby
            df(lyx(i)) = 0.
         enddo
      endif
      
      return
c     end of dfdxk.
      end

c     ------------------------------------------------------------------
      subroutine dfdyk(f,df,npt,npk,nbv,nby,nbx,ncs,lyy,lxy,sny,isy,
     *                     saym,sayp)
c     ------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      logical use_hi
      common /order/ use_hi
      dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4),lxy(nbx)
     *          ,isy(npk),sayp(npt),saym(npt)
c   note, isy: k-y-comp  -> x-comp, 
c         lyy: k-y-comp-bound -> x-comp , etc.
c
      do i = 2, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         jm  = isy(i-1)
c            df(j)= 0.5*( sayp(j)*(f(jp)+f(j)) - saym(j)*(f(j)+f(jm)) )
         df(j)= 0.5*sayp(j)*(f(jp)- f(jm))
      enddo
      
      do i = 1, nby
         i1 = lyy(i,1)
         i2 = lyy(i,2)
         f1 = f(i1)
         f2 = (f1 + f(i2))/2.
         df(i1) = sny(i)*sayp(i1)*(f2 - f1)
      enddo
      
      if(nbv.eq.1) then
c.....................set the derivative along the x boundary equal to zero.
         do i = 1, nbx
            df(lxy(i)) = 0.
         enddo
      endif

      return
c     end of dfdyk.
      end

c     ---------------------------------------------------------------------
      subroutine dfdx_dp (f,df,npt,npk,nbu,nbx,nby,ncs,lxx,lyx,snx
     *           ,npbc,lpbcw,lpbce,isk,saxm,saxp)
c     ---------------------------------------------------------------------
c     double precision version of dfdxk - only in differences of pres

      real*8 f(npt)
      integer lpbcw(npbc), lpbce(npbc),isk(npk),lxx(nbx+ncs,4),lyx(nby)
      real saxm(npt),saxp(npt),snx(nbx),df(npt)
      
      do i = 2, npk - 1
         j = isk(i)
c            df(j) = 0.5* ( saxp(j)*(f(j+1)+f(j)) - saxm(j)*(f(j)+f(j-1)) )
         df(j) = 0.5* saxp(j)*(f(j+1)-f(j-1))
      enddo
      
c....................... periodic B.C.
      do i = 1, npbc
         i2 = lpbce(i)
         f4 = f(i2)
         f3 = f(i2-1)
         i1 = lpbcw(i)
         f5 = f(i1)
         f6 = f(i1+1)
         df(i1) = 0.5* saxp(i1)*(f6 - f4)
         df(i2) = 0.5* saxp(i2)*(f5 - f3)
      enddo
      
      do i = 1, nbx
         i1 = lxx(i,1)
         i2 = lxx(i,2)
         f1 = f(i1)
         f2 = (f1 + f(i2))/2.
         df(i1) = snx(i)*saxp(i1)*(f2 - f1)
      enddo
      
      return
c     end of dfdxk.
      end

c     ------------------------------------------------------------------
      subroutine dfdy_dp(f,df,npt,npk,nbv,nby,nbx,ncs,lyy,lxy,sny,isy,
     *                     saym,sayp)
c     ------------------------------------------------------------------
      real*8 f(npt)
      real sny(nby),sayp(npt),saym(npt),df(npt)
      integer lyy(nby+ncs,4),lxy(nbx),isy(npk)
 
      do i = 2, npk-1
         j   = isy(i)
         jp  = isy(i+1)
         jm  = isy(i-1)
c            df(j)= 0.5*( sayp(j)*(f(jp)+f(j)) - saym(j)*(f(j)+f(jm)) )
         df(j)= 0.5*sayp(j)*(f(jp)- f(jm))
      enddo
      
      do i = 1, nby
         i1 = lyy(i,1)
         i2 = lyy(i,2)
         f1 = f(i1)
         f2 = (f1 + f(i2))/2.
         df(i1) = sny(i)*sayp(i1)*(f2 - f1)
      enddo
      
      return
c     end of dfdyk.
      end

c     ---------------------------------------------------------------------
      subroutine divergence (nx,ny,npt,npk,uc,vc,dx2i,dy2i,ycos,
     *                       isk,inx,iny,div)
c     ---------------------------------------------------------------------
c............ computes the divergence of (u v 0):

c            div(u,v) = d(f*uc)/dx + d(f*vc)/dy

c............ on the whole uncompressed domain, with f, uc and vc on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension uc(npt),vc(npt),div(npt),
     *          dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *          isk(npt),inx(npt),iny(npt)

      do j = 1, ny+1
         do i = 1, nx+1
            t1(i,j) = 0.
            t2(i,j) = 0.
            t3(i,j) = 0.
         enddo
      enddo
      
      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = uc(ipt)
         t2(i,j) = vc(ipt)
      enddo

      do j = 1, ny
         t1(nx+1,j) = t1(1,j)
      enddo

      do i = 1, nx
         ip = i+1
         do j = 1, ny
            up = t1(ip,j) + t1(i,j)
            t3(i,j) = t3(i,j)  + up* dx2i(i ,j)
            t3(ip,j)= t3(ip,j) - up* dx2i(ip,j)
         enddo
      enddo
      do j = 1, ny
         t3(1,j)  = t3(1,j)  + t3(nx+1,j)
      enddo

      do j = 1, ny-1
         jp = j+1
         yc = ycos(j)
         yp = ycos(jp)
         do i = 1, nx
            vp = yp*t2(i,jp) + yc*t2(i,j)
            t3(i,j) = t3(i,j)  + vp* dy2i(i,j) / yc
            t3(i,jp)= t3(i,jp) - vp* dy2i(i,jp)/ yp
         enddo
      enddo

      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         div(ipt) = t3(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine curlz (nx,ny,npt,npk,u,v,dx2i,dy2i,ycos,
     *                       isk,inx,iny,curl,nbx,nby,lxx,lyy)
c     ---------------------------------------------------------------------
c............ computes vertical component of curl of (u v 0)
c                 using Neumann boundary conditions

c                         curlz = dv/dx  - du/dy

c............ on the whole uncompressed domain, with u and v on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension u(npt),v(npt),curl(npt),lxx(nbx),lyy(nby),
     *          dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *          isk(npt),inx(npt),iny(npt)

      do j = 0, ny+1
         do i = 0, nx+1
            t1(i,j) = 0.
            t2(i,j) = 0.
         enddo
      enddo
       
      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = v(ipt)
         t2(i,j) = u(ipt)
      enddo

      if (iperx) then
         do j = 1, ny
            t1(0,j)    = t1(nx,j)
            t1(nx+1,j) = t1(1,j)
         enddo
      endif

      do j = 1, ny
         do i = 1, nx
            ip = i+1
            im = i-1
            t3(i,j) = (t1(ip,j) - t1(im,j))* dx2i(i,j)
         enddo
      enddo

      do k = 1, nbx
         ipt = lxx(k)
         i = inx(ipt)
         j = iny(ipt)
         t3(i,j) = 0.
      enddo

      do j = 1, ny
         do i = 1, nx
            jp = j+1
            jm = j-1
            t1(i,j)=(ycos(jp)*t2(i,jp)- ycos(jm)*t2(i,jm))* dy2i(i,j)/ycos(j)
         enddo
      enddo

      do k = 1, nby
         ipt = lyy(k)
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = 0.
      enddo


      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         curl(ipt) = t3(i,j) - t1(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine curlH (nx,ny,npt,psi,dept,dx2i,dy2i,ycos,
     *                       inx,iny,u,v)
c     ---------------------------------------------------------------------
c............ computes components of curl of (0 0 psi)

c            curl = ( -d(psi)/dy * e_1 + d(psi)/dx * e_2 )/ total_depth

c............ on the whole uncompressed domain, with f, uc and vc on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension psi(0:nx+1,0:ny+1),u(npt),v(npt),dx2i(0:nx+1,ny),dy2i(nx,ny),dept(npt),
     *          ycos(0:ny+1)
      dimension inx(npt),iny(npt)

      if (iperx) then
         do j = 1, ny
            psi(0,j)    = psi(nx,j)
            psi(nx+1,j) = psi(1,j)
         enddo
      endif

      do j = 1, ny
         do i = 1, nx
            ip = i+1
            im = i-1
            t2(i,j)= (psi(ip,j) - psi(im,j))* dx2i(i,j)
            jp = j+1
            jm = j-1
            t3(i,j)= (psi(i,jp)- psi(i,jm))* dy2i(i,j)
         enddo
      enddo

      do ipt = 1, npt
         i = inx(ipt)
         j = iny(ipt)
         depi = 1./dept(ipt)
         u(ipt) = - depi*t3(i,j)
         v(ipt) =   depi*t2(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine grad (nx,ny,npt,f,dx2i,dy2i,ycos,inx,iny,isk,u,v,
     *                 nbx,nby,lxx,lyy)
c     ---------------------------------------------------------------------
c............ computes horizontal components of grad of f, Neumann boundary
c                                                              conditions

c            gradH =  df/dx * e_1 + df/dy * e_2
      include 'comm_topo.h'


      dimension f(npt),u(npt),v(npt),dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *  lxx(nbx),lyy(nby)
      dimension isk(npt),inx(npt),iny(npt)

      do j = 0, ny+1
         do i = 0, nx+1
            t1(i,j) = 0.
         enddo
      enddo
       
      do ipt = 1, npt
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = f(ipt)
      enddo

      if (iperx) then
         do j = 1, ny
            t1(0,j)    = t1(nx,j)
            t1(nx+1,j) = t1(1,j)
         enddo
      endif

      do j = 1, ny
         do i = 1, nx
            ip = i+1
            im = i-1
            t2(i,j)= (t1(ip,j) - t1(im,j))* dx2i(i,j)
            jp = j+1
            jm = j-1
            t3(i,j)= (t1(i,jp)- t1(i,jm))* dy2i(i,j)
         enddo
      enddo

      do k = 1, nbx
         ipt = lxx(k)
         i = inx(ipt)
         j = iny(ipt)
         t2(i,j) = 0.
      enddo
      do k = 1, nby
         ipt = lyy(k)
         i = inx(ipt)
         j = iny(ipt)
         t3(i,j) = 0.
      enddo

      do ipt = 1, npt
         i = inx(ipt)
         j = iny(ipt)
         u(ipt) = t2(i,j)
         v(ipt) = t3(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine divergence_up (nx,ny,npt,npk,f,uc,vc,dx2i,dy2i,ycos,
     *                       isk,inx,iny,div,hupi)
c     ---------------------------------------------------------------------
c............ computes the divergence of (u v 0):

c            div(u,v) = d(f*uc)/dx + d(f*vc)/dy

c............ on the whole uncompressed domain, with f, uc and vc on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension f(npt),uc(npt),vc(npt),div(npt),
     *          dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *          isk(npt),inx(npt),iny(npt)

      do j = 1, ny+1
         do i = 1, nx+1
            t1(i,j) = 0.
            t2(i,j) = 0.
            t3(i,j) = 0.
            t4(i,j) = 0.
         enddo
      enddo
      
      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = uc(ipt)
         t2(i,j) = vc(ipt)
         t3(i,j) = f(ipt)
      enddo

      do j = 1, ny
         t1(nx+1,j) = t1(1,j)
         t3(nx+1,j) = t3(1,j)
      enddo

      do j = 1, ny
         do i = 1, nx
            ip = i+1
            fc = t3(i,j)
            fp = t3(ip,j)

            up = t1(ip,j) + t1(i,j)
            cp = hupi
            dp = 1.-cp
            if (up.lt.0) then
c  flow into cell
               cp = dp
               dp = hupi
            endif

            fp = cp*fp + dp*fc
            fdx = fp*up

            t4(i,j) = t4(i,j)  + fdx* dx2i(i ,j)
            t4(ip,j)= t4(ip,j) - fdx* dx2i(ip,j)
         enddo
         t4(1,j)  = t4(1,j)  + t4(nx+1,j)
      enddo

      do i = 1, nx
         do j = 1, ny-1
            jp = j+1
            yc = ycos(j)
            yp = ycos(jp)
            fc = t3(i,j)
            fp = t3(i,jp)

            vp = yp*t2(i,jp) + yc*t2(i,j)
            cp = hupi
            dp = 1.-cp
            if (vp.lt.0) then
c  flow into cell
               cp = dp
               dp = hupi
            endif

            fp = cp*fp + dp*fc
            fdy = fp*vp

            t4(i,j) = t4(i,j)  + fdy* dy2i(i,j) / yc
            t4(i,jp)= t4(i,jp) - fdy* dy2i(i,jp)/ yp
         enddo
      enddo

      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         div(ipt) = t4(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine divergence_mul(nx,ny,npt,npk,f,uc,vc,dx2i,dy2i,ycos,
     *                       isk,inx,iny,div)
c     ---------------------------------------------------------------------
c............ computes the divergence of (u v 0):

c            div(u,v) = d(f*uc)/dx + d(f*vc)/dy

c............ on the whole uncompressed domain, with f, uc and vc on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension f(npt),uc(npt),vc(npt),div(npt),
     *          dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *          isk(npt),inx(npt),iny(npt)

      do j = 1, ny+1
         do i = 1, nx+1
            t1(i,j) = 0.
            t2(i,j) = 0.
            t4(i,j) = 0.
         enddo
      enddo
      
      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         fi = f(ipt)
         t1(i,j) = fi*uc(ipt)
         t2(i,j) = fi*vc(ipt)
      enddo

      do j = 1, ny
         t1(nx+1,j) = t1(1,j)
      enddo

      do i = 1, nx
         ip = i+1
         do j = 1, ny
            up = t1(ip,j) + t1(i,j)
            t4(i,j) = t4(i,j)  + up* dx2i(i ,j)
            t4(ip,j)= t4(ip,j) - up* dx2i(ip,j)
         enddo
      enddo
      do j = 1, ny
         t4(1,j)  = t4(1,j)  + t4(nx+1,j)
      enddo

      do j = 1, ny-1
         jp = j+1
         yc = ycos(j)
         yp = ycos(jp)
         do i = 1, nx
            vp = yp*t2(i,jp) + yc*t2(i,j)
            t4(i,j) = t4(i,j)  + vp* dy2i(i,j) / yc
            t4(i,jp)= t4(i,jp) - vp* dy2i(i,jp)/ yp
         enddo
      enddo

      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         div(ipt) = t4(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine make_big (nx,ny,npt,psi,psin,inx,iny,isk,mask)
c     ---------------------------------------------------------------------
      include 'comm_topo.h'

      dimension psi(npt),psin(0:nx+1,0:ny+1)
      dimension isk(npt),inx(npt),iny(npt),mask(nx,ny)

      do j = 0, ny+1
         do i = 0, nx+1
            psin(i,j) = 0.
         enddo
      enddo
       
      do ipt = 1, npt
         i = inx(ipt)
         j = iny(ipt)
         psin(i,j) = psi(ipt)
      enddo

      do j = 1, ny
         do i = 1, nx
            if (mask(i,j).eq.0) then
               if (i+1.le.nx.and.mask(i+1,j).ne.0) then
                  psin(i,j) = psin(i+1,j)
               elseif (i-1.ge.1.and.mask(i-1,j).ne.0) then
                  psin(i,j) = psin(i-1,j)
               elseif (j+1.le.ny.and.mask(i,j+1).ne.0) then
                  psin(i,j) = psin(i,j+1)
               elseif (j-1.ge.1.and.mask(i,j-1).ne.0) then
                  psin(i,j) = psin(i,j-1)
               endif
            endif
         enddo
      enddo

      if (iperx) then
         do j = 1, ny
            psin(0,j)    = psin(nx,j)
            psin(nx+1,j) = psin(1,j)
         enddo
      endif

      return
      end

c     ---------------------------------------------------------------------
      subroutine make_small (nx,ny,npt,psi,psin,inx,iny,isk)
c     ---------------------------------------------------------------------
      include 'comm_topo.h'

      dimension psi(npt),psin(0:nx+1,0:ny+1)
      dimension isk(npt),inx(npt),iny(npt)

      do ipt = 1, npt
         i = inx(ipt)
         j = iny(ipt)
         psi(ipt) = psin(i,j)
      enddo

      return
      end

c     ---------------------------------------------------------------------
      subroutine divergence_tvd (nx,ny,npt,npk,f,uc,vc,dx2i,dy2i,ycos,
     *                       isk,inx,iny,div,hupi)
c     ---------------------------------------------------------------------
c............ computes the divergence of (u v 0):

c            div(u,v) = d(f*uc)/dx + d(f*vc)/dy

c............ on the whole uncompressed domain, with f, uc and vc on
c.............the compressed ocean grid, level k
      include 'comm_topo.h'

      dimension f(npt),uc(npt),vc(npt),div(npt),
     *          dx2i(0:nx+1,ny),dy2i(nx,ny),ycos(0:ny+1),
     *          isk(npt),inx(npt),iny(npt)

      do j = 1, ny+1
         do i = 1, nx+1
            t1(i,j) = 0.
            t2(i,j) = 0.
            t3(i,j) = 0.
            t4(i,j) = 0.
         enddo
      enddo
      
      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         t1(i,j) = uc(ipt)
         t2(i,j) = vc(ipt)
         t3(i,j) = f(ipt)
      enddo

      do j = 1, ny
         t1(nx+1,j) = t1(1,j)
         t3(nx+1,j) = t3(1,j)
      enddo

      do j = 1, ny
         do i = 1, nx+1
            ftp(i,j) = 0.
            ftm(i,j) = 0.
         enddo

         do i = 1, nx
            ip = i+1
            tk  = t3(i,j)
            tkp = t3(ip,j)

            wik = t1(ip,j) + t1(i,j)

            wp = 0.
            if (wik.gt.0.)wp = wik
            wm = 0.
            if (wik.lt.0.)wm = wik

            flux(i,j) = wm * tkp + wp * tk
            dt = tkp - tk
            ftp(ip,j) = wp * dt
            ftm(i ,j) = wm * dt
         enddo

         do i = 1, nx
            ip = i+1
            
            rp = 0
            if (ftp(i,j).ne.0) rp = ftp(ip,j)/ftp(i,j)
            rp = max(0.,rp)
            rm = 0
            if (ftm(ip,j).ne.0) rm = ftm(i,j)/ftm(ip,j)
            rm = max(0.,rm)
            sp = rp/(1.+rp) 
            sm = rm/(1.+rm) 
            
            flx = flux(i,j) + sp*ftp(i,j)- sm*ftm(ip,j)

            t4(i,j) = t4(i,j)  + flx* dx2i(i ,j)
            t4(ip,j)= t4(ip,j) - flx* dx2i(ip,j)
         enddo
         t4(1,j)  = t4(1,j)  + t4(nx+1,j)
      enddo

      do i = 1, nx
         do j = 1, ny
            ftp(i,j) = 0.
            ftm(i,j) = 0.
         enddo

         do j = 1, ny-1
            jp = j+1
            yc = ycos(j)
            yp = ycos(jp)

            tk  = t3(i,j)
            tkp = t3(i,jp)

            wik = yp*t2(i,jp) + yc*t2(i,j)
            wp = 0.
            if (wik.gt.0.)wp = wik
            wm = 0.
            if (wik.lt.0.)wm = wik

            flux(i,j) = wm * tkp + wp * tk
            dt = tkp - tk
            ftp(i,jp)  = wp * dt
            ftm(i,j ) = wm * dt
         enddo

         do j = 1, ny-1
            jp = j+1
            yc = ycos(j)
            yp = ycos(jp)

            rp = 0
            if (ftp(i,j).ne.0) rp = ftp(i,jp)/ftp(i,j)
            rp = max(0.,rp)
            rm = 0
            if (ftm(i,jp).ne.0) rm = ftm(i,j)/ftm(i,jp)
            rm = max(0.,rm)
            sp = rp/(1.+rp) 
            sm = rm/(1.+rm) 
            
            flx = flux(i,j) + sp*ftp(i,j)- sm*ftm(i,jp)

            t4(i,j) = t4(i,j)  + flx* dy2i(i,j) / yc
            t4(i,jp)= t4(i,jp) - flx* dy2i(i,jp)/ yp
         enddo
      enddo

      do k = 1, npk
         ipt = isk(k)
         i = inx(ipt)
         j = iny(ipt)
         div(ipt) = t4(i,j)
      enddo

      return
      end

