************************************************************************
      subroutine ice_adv (nx,ny,npt,uice,vice,cice,hice,thice,
     *                 isk,tp,fcice,fhice,fthice,dx2i,dy2i,ycos,inx,iny)
************************************************************************
      implicit real(a-h,o-z),integer(i-n)
      include 'diffiso.h'
      dimension isk(npt),inx(nx),iny(ny)

      dimension uice(npt),vice(npt),cice(npt),hice(npt),thice(npt),
     *     tp(npt),fcice(npt),fhice(npt),fthice(npt),
     *          dx2i(1), dy2i(1), ycos(1)
c
      
c    input:  
c         uice, vice           ice velocity = ocean velocity
c         cice, hice, thice    ice tracer quantities
c         fx,fy                from ice_force

c    output: 
c         fcice, fhice, fthice     
c           
c  d(hu)/dt = - f x (hu) - h*grad(P_top) + [fx,fy]

c  d(h*cice)/dt  = advection
c  d(h*hice)/dt  = advection
c  d(h*thice)/dt = advection

c  make the forcings for the ice momentum (uice,vice) advection:
c     this is a crude approximation, until ice_force is finished
      
c  make the forcings for the tracer (cice,hice,thice) advection:
      
      call divergence_tvd(nx,ny,npt,npt,cice,uice,vice,
     *        dx2i,dy2i,ycos,isk,inx,iny,tp,hupi_tr)
      
      do i = 1, npt
         fcice(i) = fcice(i) - tp(i)
      enddo
      
      call divergence_tvd(nx,ny,npt,npt,hice,uice,vice,
     *        dx2i,dy2i,ycos,isk,inx,iny,tp,hupi_tr)

      do i = 1, npt
         fhice(i) = fhice(i) - tp(i)
      enddo
      
      call divergence_tvd(nx,ny,npt,npt,thice,uice,vice,
     *        dx2i,dy2i,ycos,isk,inx,iny,tp,hupi_tr)

      do i = 1, npt
         fthice(i) = fthice(i) - tp(i)
      enddo
      
      return
      end
      
      subroutine ice_diff(npt,nz,cice,hice,thice,fcice,fhice,fthice,
     *       lxxk,lyyk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek,
     *       saxmk,saxpk,saymk,saypk)

c  d(hu)/dt =  "advection stuff from ice_move" + div(k*grad(hu)) 
      include 'comm_para.h'
      include 'comm_diff.h'
      include 'diffiso.h'
      include 'icedyn.h'

      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +        ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      dimension cice(npt),hice(npt),thice(npt), tp(npt,2)
      dimension fcice(npt),fhice(npt),fthice(npt)
      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)
      
      coef_diff = dyice_diff

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

c  will need to also compute ice momentum diffusion here

      call dfdx_a2c(cice,tp,npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,dxp)
      call dfdx_c2a(tp,tp(1,2),npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,saxmk,saxpk)
      
      do i = 1, npt
         fcice(i) = fcice(i) + coef_diff*tp(i,2)
      enddo

      call dfdy_a2c(cice,tp,npt,npk,nyk,nck,lyyk,snyk,isyk,dyp)
      call dfdy_c2a(tp,tp(1,2),npt,npk,nyk,nck,lyyk,
     *     snyk,isyk,csy,csyc,saymk,saypk)
      
      do i = 1, npt
         fcice(i) = fcice(i) + coef_diff*tp(i,2)
      enddo
      
      call dfdx_a2c(hice,tp,npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,dxp)
      call dfdx_c2a(tp,tp(1,2),npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,saxmk,saxpk)
      
      do i = 1, npt
         fhice(i) = fhice(i) + coef_diff*tp(i,2)
      enddo
      
      call dfdy_a2c(hice,tp,npt,npk,nyk,nck,lyyk,snyk,isyk,dyp)
      call dfdy_c2a(tp,tp(1,2),npt,npk,nyk,nck,lyyk,
     *     snyk,isyk,csy,csyc,saymk,saypk)
      
      do i = 1, npt
         fhice(i) = fhice(i) + coef_diff*tp(i,2)
      enddo
      
      call dfdx_a2c(thice,tp,npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,dxp)
      call dfdx_c2a(tp,tp(1,2),npt,npk,nxk,nck,lxxk,
     *     snxk,npbk,lpbcwk,lpbcek,isk,saxmk,saxpk)
      
      do i = 1, npt
         fthice(i) = fthice(i) + coef_diff*tp(i,2)
      enddo

      call dfdy_a2c(thice,tp,npt,npk,nyk,nck,lyyk,snyk,isyk,dyp)
      call dfdy_c2a(tp,tp(1,2),npt,npk,nyk,nck,lyyk,
     *     snyk,isyk,csy,csyc,saymk,saypk)
      
      do i = 1, npt
         fthice(i) = fthice(i) + coef_diff*tp(i,2)
      enddo
      
      return
      end
      
      subroutine ice_updat(npt,binv,abinv,h,cice,hice,thice,
     *                                      fcice,fhice,fthice)

      implicit real(a-h,o-z),integer(i-n)
      dimension cice(npt),hice(npt),thice(npt)
      dimension fcice(npt),fhice(npt),fthice(npt)

c  d(hu)/dt =  advection + diffusion

      do i = 1, npt
         cice(i)  = cice(i) + binv*fcice(i)
         thice(i) = thice(i)+ binv*fthice(i)
         hice(i)  = hice(i) + binv*fhice(i)

         fcice(i) = abinv*fcice(i)
         fhice(i) = abinv*fhice(i)
         fthice(i)= abinv*fthice(i)
      enddo
      
      
      return
      end