subroutine reinit(phi,nx,ny)
*******************************************************************
*     reinitialize the level set function so it is the distance
*     function from the zero level set.
*     Given phi on the grid, calculate the time it takes for a
*     propagating front with propagating speed 1, which starts as 
*     the zero level set of phi, to arrive at the grid point (i,j). 
*     The time for this particular grid point is stored in tarv(i,j).
*     Two sides of the zero level set are dealt with separately in
*     two sweeps (k=1,2 in loop 2).
*
*         phi:    level set function;
*         fo,fn:  local arrays used to contain the propagating
*                 information of the "front", fo (initial) and
*                 fn (final) states.
*******************************************************************
      integer nx,ny,nxm,nym,mxym,i,j,k,ndone,kdel,kt,kk
      include 'gridsz'
      parameter(kdel=10)
      logical done(0:nxm+1,0:nym+1)
      double precision t,dt,del,h,hinv,sign,dth
      double precision phi(0:nx+1,0:ny+1)
      double precision fo(0:nxm+1,0:nym+1),fn(0:nxm+1,0:nym+1),
     +                 tarv(0:nxm+1,0:nym+1)
      common/meshsz/h,hinv
      del = dble(kdel)*h
      dt = 0.5D0*h
      dth = 0.5D0
      kt = 2*kdel
c
      do 2 kk=1,2
      sign = (-1.D0)**(kk+1)
      ndone = 0
      t = 0.D0
      do 8 j=0,ny+1
      do 8 i=0,nx+1
         fo(i,j) = sign*phi(i,j)
 8    continue
c     mark all the points on the nonnegative side:
      do 10 j=0,ny+1
      do 10 i=0,nx+1
         done(i,j) = fo(i,j) .ge. 0.D0
         if(done(i,j)) ndone = ndone + 1
 10   continue
*
      do 5 k=1,kt
*
      do 50 j=1,ny
      do 50 i=1,nx
         fn(i,j) = fo(i,j) + dth*dsqrt(
     +             dmin1(fo(i,j)-fo(i-1,j),0.D0)**2 +
     +             dmax1(fo(i+1,j)-fo(i,j),0.D0)**2 +
     +             dmin1(fo(i,j)-fo(i,j-1),0.D0)**2 +
     +             dmax1(fo(i,j+1)-fo(i,j),0.D0)**2 )
 50   continue
      do 60 j=1,ny
         fn(0,j) = fn(1,j) -h
         fn(nx+1,j) = fn(nx,j) - h
 60   continue
      do 65 i=1,nx
         fn(i,0) = fn(i,1) - h
         fn(i,ny+1) = fn(i,ny) - h
 65   continue
      fn(0,0) = fn(1,1) - h
      fn(nx+1,0) = fn(nx,1) - h
      fn(0,ny+1) = fn(1,ny) - h
      fn(nx+1,ny+1) = fn(nx,ny) - h
c
      t = t + dt
c
      do 20 j=0,ny+1
      do 20 i=0,nx+1
         if( .not. done(i,j) .and. fn(i,j) .ge. 0.D0 ) then
            tarv(i,j) = t - dt*fn(i,j)/(fn(i,j)-fo(i,j))
            ndone = ndone + 1
            done(i,j) = .true.
         endif
 20   continue
c
      if(ndone .eq. (nx+2)*(ny+2)) goto 2
c
      do 25 j=0,ny+1
      do 25 i=0,nx+1
         fo(i,j) = fn(i,j)
 25   continue
c
 5    continue
*
      do 30 j=0,ny+1
      do 30 i=0,nx+1
         if( .not. done(i,j) ) tarv(i,j) = del
 30   continue

 2    continue

      do 70 j=0,ny+1
      do 70 i=0,nx+1
         phi(i,j) = dsign(tarv(i,j),phi(i,j))
 70   continue
      end