subroutine index(m,n,n1,a,b,c,d,ix0,ix1,iy0,iy1,
     1     x,y,x1,y1,cox1,coy1,hs1,hs2,ijump,xyjump)
C
C*******************************************************************************
C								                
C   index indices all the grid in the following way:                            
C        ijump(i,j,1) = 1,  regular grid inside the interface              
C        ijump(i,j,1) = 2,  irregular grid inside the interface
C        ijump(i,j,1) = 3,  irregular grid on the interface
C        ijump(i,j,1) = 4,  irregular grid outside the interface
C        ijump(i,j,1) = 5,  regular grid outside the interface
C   Also index gives information of a point on the interface, either on
C   x = x(i) or y = y(j) line depending on which point has shorter distance
C   to (x_i, y_j).
C        ijump(i,j,2), the index for the spline interval, say (s_k, s_{k+1})
C	 xyjump(i,j,8), the corresponding parameter s for such point
C
C  a,b,c,d: On entry, the rectangular region, a <= x <= b, c <= y <= d.
C  m,n: On entry, the number of the grid in x and y direction respectively.
C  x(m), y(n): On entry, double precision, the x-grid and y-grid.
C	x(i) = a + (i-1)*(b-a)/m, y(j) = c + (j-1)*(d-c)/n.
C  n1: On entry, the number of the point on the interface.
C  x1(n1+1), y1(n1+1): On entry, the given points on the interface.
C  cox1(n1+1), coy1(n1+1): On entry, the moments for the splines x= x(s) and
C			   y = y(s), got from subroutine splcl.
C  hs1(n1+1), hs2(n1+1): On entry, the parameter for the splines x= x(s) and
C                        y = y(s), got from subroutine splcl.
C  xmin(n1+1), xmax(n1+1), ymin(n1+1), ymax(n1+1): Output, The xmin(i) and 
C              xmax(i) are min {x(s)} and max {x(s)} over the interval 
C	       (s_i, s_{i+1}) respectively. 
C  ix0, ix1, iy0, iy1, Output, the low and up bound for the indeces which
C        contain the irregular grid.
C
C*******************************************************************************
C
C	implicit double precision (a-h,o-z)
	implicit real*8 (a-h,o-z)
	parameter(nn = 20000)
        dimension x(m),y(n),cox1(n1+1),coy1(n1+1),x1(n1+1),y1(n1+1),
     1 	    hs1(n1+1),hs2(n1+1),xmin(nn),xmax(nn),ymin(nn),
     2      ymax(nn),w1(nn),w2(nn),ind(nn),ijump(m,n,3),
     3	    xyjump(m,n,9),ind2(nn),w3(nn)

	h = (b -a)/m
	tol = 1.0d-8

C----- First find the extreme values of x and y in each spline interval. -------

	xmaxt = x1(1)
	xmint = x1(1)
	ymint = y1(1)
	ymaxt = y1(1)

	do 10 i=1,n1

	      xt = x1(i+1)
	      a1 = (xt-x1(i))/hs2(i+1) - (2.0d0*cox1(i) +
     1             cox1(i+1) )*hs2(i+1)/6.0d0
              a2 = cox1(i)/2.0d0
              a3 = (cox1(i+1) - cox1(i) )/(hs2(i+1)*6.0d0)
C		  The spline coeffcients on the interval (s(i),s(i+1))

	   xmin(i) = min(x1(i),xt)
	   xmax(i) = max(x1(i),xt)

	   call rootp2(a1,2.0d0*a2,3.0d0*a3,r1,r2,info)
	     if(info .eq. 2) then
	        r1 = r1 + hs1(i)
	        r2 = r2 + hs1(i)
	        if(hs1(i) .le. r1 .and. r1 .le. hs1(i+1)) then
	          t = r1 - hs1(i)
	          r1 = x1(i) + t*(a1 + t*(a2 + a3*t))
	          xmin(i) = min(xmin(i),r1)
	          xmax(i) = max(xmax(i),r1)
	        endif

		if(hs1(i) .le. r2 .and. r2 .le. hs1(i+1)) then
		  t = r2 - hs1(i)
                  r2 = x1(i) + t*(a1 + t*(a2 + a3*t))
                  xmin(i) = min(xmin(i),r2)
                  xmax(i) = max(xmax(i),r2)
                endif

	        if(xmint .gt. xmin(i)) xmint = xmin(i)
	        if(xmaxt .lt. xmax(i)) xmaxt = xmax(i)
	     endif

              xt = y1(i+1)

              a1 = (xt-y1(i))/hs2(i+1) - (2.0d0*coy1(i) +
     1             coy1(i+1) )*hs2(i+1)/6.0d0
              a2 = coy1(i)/2.0d0
              a3 = (coy1(i+1) - coy1(i) )/(hs2(i+1)*6.0d0)

           ymin(i) = min(y1(i),xt)
           ymax(i) = max(y1(i),xt)

           call rootp2(a1,2.0d0*a2,3.0d0*a3,r1,r2,info)
             if(info .eq. 2) then
	     r1 = r1 + hs1(i)
             r2 = r2 + hs1(i)
                if(hs1(i) .le. r1 .and. r1 .le. hs1(i+1)) then
                  t = r1 - hs1(i)
                  r1 = y1(i) + t*(a1 + t*(a2 + a3*t))
                  ymin(i) = min(ymin(i),r1)
                  ymax(i) = max(ymax(i),r1)
                endif

                if(hs1(i) .le. r2 .and. r2 .le. hs1(i+1)) then
                  t = r2 - hs1(i)
                  r2 = y1(i) + t*(a1 + t*(a2 + a3*t))
                  ymin(i) = min(ymin(i),r2)
                  ymax(i) = max(ymax(i),r2)
                endif

                if(ymint .gt. ymin(i)) ymint = ymin(i)
                if(ymaxt .lt. ymax(i)) ymaxt = ymax(i)
             endif

 10	continue

C------------------- End of the computation of extreme values ------------------
C
        t = (xmaxt - a)/h
        if ( int(t) .eq. t) then
           ix1= int(t)
        else
           ix1= int(t) + 1
        endif

           ix0= int( (xmint - a)/h )

        iy0= int((ymint - c)/h )
        t = (ymaxt - c)/h
        if ( int(t) .eq. t) then
           iy1= int(t)
        else
           iy1= int(t) + 1
        endif

C
C  Having determined the low and up bound of the indeces which contain the
C  irregular grid point. Now we gre going to find all the intersection 
C  between x = x(i) and the interface x= x(s).
                    
	do 300 i=ix0-1,ix1+1
	   do 310 j=iy0-1,iy1+1
		xyjump(i,j,7) = 2.0d0
 310	   continue
 300	continue

	do 20 i=ix0,ix1
	  nr = 0
	  do 30 k=1,n1

	    if( xmin(k) .le. x(i) .and. x(i) .le. xmax(k)) then
	      a0 = x1(k) - x(i)
	      xt = x1(k+1)
	      a1 = (xt-x1(k))/hs2(k+1) - (2.0d0*cox1(k) +
     1             cox1(k+1) )*hs2(k+1)/6.0d0
              a2 = cox1(k)/2.0d0
              a3 = (cox1(k+1) - cox1(k) )/(hs2(k+1)*6.0d0)

	      call rootp3(a0,a1,a2,a3,r1,r2,r3,tol,info)

C   In case the one of two ends x1(k), x1(k+1) is the intersection and the
C   root finding progam failed to find it due to the round-off errors.

                if( abs(x1(k) - x(i)) .le. 1.0d-6 .or.
     1                  abs(xt - x(i)) .le. 1.0d-6 ) then

                  if(info .eq. 1) then
                    info = 2
             if( abs(x1(k) - x(i)) .le. abs(xt - x(i))) then
                    r2 = 0.0d0
             else
                    r2 = hs2(k+1)
             endif
                  else
                    if(info .eq. 2) then
                        info = 3
             if( abs(x1(k) - x(i)) .le. abs(xt - x(i))) then
                    r3 = 0.0d0
             else
                    r3 = hs2(k+1)
             endif
                    else
                        info = 4
             if( abs(x1(k) - x(i)) .le. abs(xt - x(i))) then
                    r4 = 0.0d0
             else
                    r4 = hs2(k+1)
             endif
                    endif
                  endif
                endif

C ----- If the curve is: x(s) = x(i), y(s) = ymin(k) + (ymax(k)-ymin(k))
C       *(s- hs1(k))/hs2(k+1), then all the points between ymin(i) and
C       ymax(i) are on the curve.   --------------------------------------------

	      if(info .eq. 110) then
		t = ( ymin(k) - c )/h
	        if(t.eq.int(t)) then
		   j0 = int(t)
		else
		   j0 = int(t) + 1
		endif
		j1 = int(( ymax(k) - c)/h )
		do 40 j=j0,j1
		   ijump(i,j,1) = 3
		   ijump(i,j,2) = k
		   xyjump(i,j,8) = (y(j)-ymin(k))*hs2(k+1)/(ymax(k)
     1		         -ymin(k)) + hs1(k)
 40		continue
		goto 30
C-------------goto next spline interval ----------------------------------------
	      endif

C --Determine the number of zeros, the curve parameters and the y coordinates --
	     yt = y1(k+1)

	  do 50 k1 = 1, info
	    if(k1.eq.1) then
	      rt = r1 +hs1(k)
	      if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
	 	nr = nr + 1
		w1(nr) = rt
		call splval2(rt,y1(k),yt,coy1(k),coy1(k+1),hs1(k),
     1			    hs2(k+1),w2(nr))
		ind(nr) = k
	      endif
	    endif
	    if(k1 .eq. 2) then
              rt = r2 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
		nr = nr + 1
                w1(nr) = rt 
                call splval2(rt,y1(k),yt,coy1(k),coy1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
		ind(nr) = k
	      endif
            endif
            if(k1 .eq. 3) then
              rt = r3 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
                nr = nr + 1
                w1(nr) = rt 
                call splval2(rt,y1(k),yt,coy1(k),coy1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
                ind(nr) = k
	      endif
            endif
            if(k1 .eq. 4) then
              rt = r4 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
                nr = nr + 1
                w1(nr) = rt
                call splval2(rt,y1(k),yt,coy1(k),coy1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
                ind(nr) = k
	call splval(n,rt,x1,y1,cox1,coy1,hs1,hs2,x4,y4,info)

              endif
            endif
 50	  continue
	endif
 30	continue

	if(nr .eq. 0 ) goto 20
	do 60 j=iy0-1,iy1+1
	  if(ijump(i,j,1) .eq. 3) goto 60
C-------- No zeros or only have the case when x(s) = x(i) ----------------------

	do 55 k8=1,nr
  	   ind2(k8) = k8
 	   w3(k8) = dabs(y(j)-w2(k8))
 55	continue

C  Now ordering the | y(j) - y( s_p) | in descending.

	  do 65 k2=1,nr
	     do 70 k1=k2+1,nr
		if(w3(k1) .lt. w3(k2) ) then
		   t = w3(k2)
		   w3(k2) = w3(k1)
		   w3(k1) = t

		   it = ind2(k2)
		   ind2(k2) = ind2(k1)
		   ind2(k1) = it
		endif
 70	     continue
 65	  continue

	  do 80 k0=1,nr
	     k3 = ind2(k0) 

	t = y(j)-w2(k3)
	if( abs(t) .le. 1.0d-6) then
	   ijump(i,j,1) = 3
	   xyjump(i,j,7) = abs(t)
	   ijump(i,j,2) = ind(k3)
	   xyjump(i,j,8) = w1(k3)
	   goto 60
	endif
	
C------------------- Compute dx ------------------------------------------------
	
	kk = ind(k3)
	xt = x1(kk+1)

	      a1 = (xt-x1(kk))/hs2(kk+1) - (2.0d0*cox1(kk) +
     1             cox1(kk+1) )*hs2(kk+1)/6.0d0
              a2 = cox1(kk)/2.0d0
              a3 = (cox1(kk+1) - cox1(kk) )/(hs2(kk+1)*6.0d0)
              t = w1(k3) - hs1(kk)
	 dx = a1 + t*(2.0d0*a2 + 3.0d0*a3*t)

C----------- Finish the computation dx. Now determine the index ----------------

	 if( abs(dx).le.1.0d-4) goto 80

	 t = y(j) - w2(k3)
	 if ( t*dx .gt. 0.0d0) then
	    if(abs(t) .ge. h) then
		ijump(i,j,1) = 1
	        xyjump(i,j,7) = abs(t)
	    else
	        ijump(i,j,1) = 2
		ijump(i,j,2) = kk
		xyjump(i,j,8) = w1(k3)
	        xyjump(i,j,7) = dabs(t)
	    endif
	else
            if(abs(t) .gt. h) then
                ijump(i,j,1) = 5
		xyjump(i,j,7) = abs(t)
            else
                ijump(i,j,1) = 4
                ijump(i,j,2) = kk
		xyjump(i,j,8) = w1(k3)
                xyjump(i,j,7) = abs(t)
            endif
	endif
	goto 60
 80	continue
 60	continue

 20	continue

c       do 2000 j=n-1,1,-1
c         write(*,1000) (ijump(i,j,1),i=1,m-1)
c2000     continue
c1000    format(39I2)

C------------------ Finishing indexing in x-direction. Now in y-direction.------

	do 120 j=iy0,iy1
	  nr = 0
	  do 130 k=1,n1

	    if( ymin(k) .le. y(j) .and. y(j) .le. ymax(k)) then
	      a0 = y1(k) - y(j)
	      xt = y1(k+1)
	      a1 = (xt-y1(k))/hs2(k+1) - (2.0d0*coy1(k) +
     1             coy1(k+1) )*hs2(k+1)/6.0d0
              a2 = coy1(k)/2.0d0
              a3 = (coy1(k+1) - coy1(k) )/(hs2(k+1)*6.0d0)
	      call rootp3(a0,a1,a2,a3,r1,r2,r3,tol,info)

		if( abs(y1(k) - y(j)) .le. 1.0d-6 .or.
     1			abs(xt - y(j)) .le. 1.0d-6 ) then
		  if(info .eq. 1) then
		    info = 2
	     if( abs(y1(k) - y(j)) .le. abs(xt - y(j))) then
		    r2 = 0.0d0
	     else
		    r2 = hs2(k+1)
	     endif
		  else
		    if(info .eq. 2) then
			info = 3
             if( abs(y1(k) - y(j)) .le. abs(xt - y(j))) then
                    r3 = 0.0d0
             else
                    r3 = hs2(k+1)
             endif
		    else
 			info = 4
             if( abs(y1(k) - y(j)) .le. abs(xt - y(j))) then
                    r4 = 0.0d0
             else
                    r4 = hs2(k+1)
             endif
		    endif
		  endif
		endif

C ----- If the curve is: y(s) = y(j), x(s) = xmin(k) + (xmax(k)-xmin(k))
C       *(s- hs1(k))/hs2(k+1), then all the points between xmin(i) and
C       xmax(i) are on the curve.   --------------------------------------------

	      if(info .eq. 10) then
		t = ( xmin(k) - a )/h
	        if(t.eq.int(t)) then
		   j0 = int(t)
		else
		   j0 = int(t) + 1
		endif
		j1 = int(( xmax(k) - a)/h )
		do 140 i=j0,j1
		   ijump(i,j,1) = 3
		   xyjump(i,j,7) = 0.0d0
		   ijump(i,j,2) = k
		   xyjump(i,j,8) = (x(i)-xmin(k))*hs2(k+1)/(xmax(k)
     1		         -xmin(k)) + hs1(k)
 140		continue
		goto 130
C-------------goto next spline interval ----------------------------------------
	      endif

C --Determine the number of zeros, the curve parameters and the x coordinates --
	     yt = x1(k+1)

	  do 150 k1 = 1, info
	    if(k1.eq.1) then
	      rt = r1 +hs1(k)
	      if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
	 	nr = nr + 1
		w1(nr) = rt
		call splval2(rt,x1(k),yt,cox1(k),cox1(k+1),hs1(k),
     1			    hs2(k+1),w2(nr))
		ind(nr) = k
	      endif
	    endif
	    if(k1 .eq. 2) then
              rt = r2 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
		nr = nr + 1
                w1(nr) = rt 
                call splval2(rt,x1(k),yt,cox1(k),cox1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
		ind(nr) = k
	      endif
            endif
            if(k1 .eq. 3) then
              rt = r3 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
                nr = nr + 1
                w1(nr) = rt 
                call splval2(rt,x1(k),yt,cox1(k),cox1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
                ind(nr) = k
	      endif
            endif
            if(k1 .eq. 4) then
              rt = r4 + hs1(k)
              if( hs1(k) .le. rt .and. rt .le. hs1(k+1)) then
                nr = nr + 1
                w1(nr) = rt
                call splval2(rt,x1(k),yt,cox1(k),cox1(k+1),hs1(k),
     1                      hs2(k+1),w2(nr))
                ind(nr) = k
              endif
            endif
 150	  continue
	endif
 130	continue

	if(nr .eq. 0 ) goto 120

	do 160 i=ix0-1,ix1+1
	  if(ijump(i,j,1) .eq. 3) goto 160
C-------- No zeros or only have the case when y(s) = y(j) ----------------------

          do 155 k1=1,nr
             ind2(k1) = k1
             w3(k1) = abs(x(i)-w2(k1))
 155       continue

          do 165 k2=1,nr
             do 170 k1=k2+1,nr
                if(w3(k1) .lt. w3(k2) ) then
                   t = w3(k2)
                   w3(k2) = w3(k1)
                   w3(k1) = t

                   it = ind2(k2)
                   ind2(k2) = ind2(k1)
                   ind2(k1) = it
                endif
 170          continue
 165       continue

          do 180 k0=1,nr
             k3 = ind2(k0)
	t = x(i)-w2(k3)
	if( abs(t) .le. 1.0d-6) then
	   ijump(i,j,1) = 3
	   xyjump(i,j,7) = abs(t)
	   ijump(i,j,2) = ind(k3)
	   xyjump(i,j,8) = w1(k3)
	   goto 160
	endif
	
C------------------- Compute dy ------------------------------------------------

	kk = ind(k3)
	xt = y1(kk+1)

	      a1 = (xt-y1(kk))/hs2(kk+1) - (2.0d0*coy1(kk) +
     1             coy1(kk+1) )*hs2(kk+1)/6.0d0
              a2 = coy1(kk)/2.0d0
              a3 = (coy1(kk+1) - coy1(kk) )/(hs2(kk+1)*6.0d0)
              t = w1(k3) - hs1(kk)
	 dy = a1 + t*(2.0d0*a2 + 3.0d0*a3*t)

C----------- Finish the computation dy. Now determine the index ----------------

	 if( abs(dy).le.1.0d-4) goto 180
	
	 t = x(i) - w2(k3)
	 if ( t*dy .gt. 0.0d0) then
	    if(abs(t) .gt. h) then
		if(ijump(i,j,1) .eq. 0) then
		   ijump(i,j,1) = 5
		   xyjump(i,j,7) = abs(t)
		endif
	    else
	      if( abs(x(i)-w2(k3)) .lt. xyjump(i,j,7)) then
		ijump(i,j,1) = 4
		ijump(i,j,2) = kk
		xyjump(i,j,7) = abs(x(i)-w2(k3))
		xyjump(i,j,8) = w1(k3)
	      endif
	    endif
	else
            if(abs(t) .ge. h) then
                if(ijump(i,j,1) .eq. 0) then
		   ijump(i,j,1) = 1
		   xyjump(i,j,7) = abs(t)
		endif	
            else
 	      if( abs(x(i)-w2(k3)) .lt. xyjump(i,j,7)) then
                ijump(i,j,1) = 2
                ijump(i,j,2) = kk
		xyjump(i,j,8) = w1(k3)
		xyjump(i,j,7) = abs(x(i)-w2(k3))
	      endif
            endif
	endif
	
	goto 160
 180	continue
 160	continue
 120    continue

	do 200 i=ix0,ix1
	   do 210 j=iy0,iy1
	      if(ijump(i,j,1) .eq. 3) then
C                 labeling the grid point next to the grid point which the 
C		  interface going through. Only some regular grid points
C		  outside of the interface may not be labled.

	if(ijump(i-1,j,1) .eq. 0 .or. ijump(i-1,j,1) .eq. 5) then
      		    ijump(i-1,j,1) = 4
		    ijump(i-1,j,2 ) = ijump(i,j,2)
		    xyjump(i-1,j,8) = xyjump(i,j,8)
                    xyjump(i-1,j,7) = h
	endif

        if(ijump(i+1,j,1) .eq. 0 .or. ijump(i+1,j,1) .eq. 5) then
                    ijump(i+1,j,1) = 4
		    ijump(i+1,j,2 ) = ijump(i,j,2)
                    xyjump(i+1,j,8) = xyjump(i,j,8)
                    xyjump(i+1,j,7) = h
        endif

        if(ijump(i,j-1,1) .eq. 0 .or. ijump(i,j-1,1) .eq. 5) then
                    ijump(i,j-1,1) = 4
                    ijump(i,j-1,2 ) = ijump(i,j,2)
                    xyjump(i,j-1,8) = xyjump(i,j,8)
                    xyjump(i,j-1,7) = h
        endif

        if(ijump(i,j+1,1) .eq. 0 .or. ijump(i,j+1,1) .eq. 5) then
                    ijump(i,j+1,1) = 4
                    ijump(i,j+1,2 ) = ijump(i,j,2)
                    xyjump(i,j+1,8) = xyjump(i,j,8)
                    xyjump(i,j+1,7) = h
        endif

	      endif
 210	   continue
 200	continue

C--------------- Lable those points left as regular points outside of the 
C                Interface  --------------------------------------------

	do 380 i=1,m
	  do 390 j=1,n
		if(ijump(i,j,1) .eq. 0) ijump(i,j,1) = 5
 390	  continue
 380	continue

C-------------------------------------------------------------------------------

C ------ Store some information at irregular grid points for later use

        do 410 i=ix0-1, ix1+1
           do 420 j=iy0-1, iy1+1

             if(ijump(i,j,1) .ne. 1 .and. ijump(i,j,1) .ne. 5) then

        ks = ijump(i,j,2)
        s = xyjump(i,j,8)

        call splval4(s,x1(ks),x1(ks+1),cox1(ks),cox1(ks+1)
     1                    ,hs1(ks),hs2(ks+1),x3,dxs,dxss)
        call splval4(s,y1(ks),y1(ks+1),coy1(ks),coy1(ks+1)
     1                     ,hs1(ks),hs2(ks+1),y3,dys,dyss)

        tlen = dsqrt(dxs*dxs + dys*dys)
        coc = dys/tlen
        sic = -dxs/tlen
        dsn = 1.0d0/(coc*dys - sic*dxs)
        dxx = (dxss*coc + dyss*sic)*dsn*dsn

c	if(ks.eq.1) write(*,*)'dxx ',dxx,ks,s,dxs,
c    1  	dys,dsn,coc,sic,dxss,dyss,x1(ks),y1(ks)

c       if(ks.eq.1) write(*,*)'dxx2 ', y1(ks),y1(ks+1),coy1(ks)
c    1    ,coy1(ks+1),hs1(ks),hs2(ks+1),y3,dys,dyss

	dsnn = -(dyss*coc - dxss*sic)*dsn*dsn*dsn
        xyjump(i,j,1) = coc
        xyjump(i,j,2) = sic
        xyjump(i,j,3) = dsn
        xyjump(i,j,4) = dxx
        xyjump(i,j,5) = dsnn

	xyjump(i,j,6) = x3 
	xyjump(i,j,7) = y3 
	    endif
 420       continue
 410     continue

        return
        end