c-----------------------------------------------------------------------
c	masks.f
c	Written by Mijan Huq Fri Aug 25 15:54:02 CDT 1995
c	Routines to 
c		1. initialize the masked region given a mask center and a mask radius
c		2. update the masked region given an apparent horizon surface.
c	Added routines for updating points that have been added to computatational
c	domain that were originally undefined.
c-----------------------------------------------------------------------
c
		subroutine initmask(msk, grx, gry, grz, mskrad, mskcntr, 
     .				hx,hy,hz, min,max,NX, NY, NZ)
			implicit none
			
			integer	NX, NY, NZ

			real	msk(NX,NY,NZ), grx(NX), gry(NY), grz(NZ), 
     .			mskrad,mskcntr(3), hx,hy,hz, min,max

c		Local variables

			integer	i, j, k

			real		rr

			hx = (max - min)/float(NX - 1)
			hy = (max - min)/float(NY - 1)
			hz = (max - min)/float(NZ - 1)

			do i = 1, NX
				grx(i) = hx*float(i-1) + min
			end do

			do i = 1, NY
				gry(i) = hy*float(i-1) + min
			end do

			do i = 1, NZ
				grz(i) = hz*float(i-1) + min
			end do

			do i = 1, NX
				do j = 1, NY
					do k = 1, NZ
						rr = sqrt((grx(i) - mskcntr(1))**2 + 
     .				(gry(j) - mskcntr(2))**2 + 
     .				(grz(k) - mskcntr(3))**2)
						if(rr .lt. mskrad)then
							msk(i,j,k) = 0.0
						else
							msk(i,j,k) = 1.0
						end if
					end do
				end do
			end do	

			write(0,*)'initialized data'
			write(0,*)'hx = ',hx
			write(0,*)'hy = ',hy
			write(0,*)'hz = ',hz

			write(0,*)'NX,NY,NZ=>',NX,NY,NZ
			write(0,*)'min max =>',min,max
			
			write(0,*)'mask radius = ',mskrad
			write(0,*)'mask center = ',mskcntr(1),mskcntr(2),mskcntr(3)

			return

			end
c-----------------------------------------------------------------------
c
c	Routine to redefine center and regrid 2-surface from original surface
c
		subroutine nsrfcntr(S,Sf,Sc,center,dtheta,dphi,NST,NSP)
			implicit none

			integer NST, NSP

			real	S(NST,NSP,3), Sc(NST,NSP,3), Sf(NST,NSP), center(3),
     .			dtheta, dphi

c		Local variables

			integer	i, j, k, nt, np

			real	nc(3), ssq(3), Pi, th, ph, rn, tn,pn

			Pi = 4.0 * atan(1.0)



c		first calculate approximate new center

			ssq(1) = 0.0
			ssq(2) = 0.0
			ssq(3) = 0.0

			do i = 1, NST
				do j = 1, NSP
					ssq(1) = ssq(1) + Sc(i,j,1)
					ssq(2) = ssq(2) + Sc(i,j,2)
					ssq(3) = ssq(3) + Sc(i,j,3)
				end do
			end do

			nc(1) = ssq(1)/float(NST*NSP)
			nc(2) = ssq(2)/float(NST*NSP)
			nc(3) = ssq(3)/float(NST*NSP)

c		Now calculate new spherical coordinate points

			do i = 1, NST
				do j = 1, NSP
					S(i,j,1) = -99999.0
				end do
			end do

			do i = 1, NST
				do j = 1, NSP
					rn = sqrt((Sc(i,j,1) - nc(1))**2 + 
     .		(Sc(i,j,2) - nc(2))**2 + (Sc(i,j,3) - nc(3))**2)

					tn = acos((Sc(i,j,3) - nc(3))/rn)

					pn = 2.0*Pi + atan2((Sc(i,j,2)-nc(2)),
     .					(Sc(i,j,1)-nc(1)))

					if(pn .ge. 2.0 * Pi)then
						pn = pn - 2.0*Pi
					end if

c			Find nearest new grid point
					nt = int(tn/dtheta) + 1
					np = int(pn/dphi) + 1

					write(0,*)nt,np

c			For now do a simple copy

					S(nt,np,1) = rn

				end do
			end do

			do i = 1, NST
				do j =1 ,NSP
					if(S(i,j,1) .eq. -99999.0)then
						write(0,*)'ERROR in nsrfcntr : missed point',i,j
					else

						Sf(i,j) = S(i,j,1)

						Sc(i,j,1) = S(i,j,1)*sin(S(i,j,2))*cos(S(i,j,3))
						Sc(i,j,2) = S(i,j,1)*sin(S(i,j,2))*sin(S(i,j,3))
						Sc(i,j,3) = S(i,j,1)*cos(S(i,j,2))

					end if
				end do
			end do
			return
		end
				
c-----------------------------------------------------------------------
c	Routine to update the mask on the next slice

		subroutine updmsk(S,Sc,grx,gry,grz,msk,mskcntr,center,
     .				NST,NSP,NX,NY,NZ)
			implicit none

			integer	NST, NSP, NX, NY, NZ

			real	S(NST,NSP,3), Sc(NST,NSP,3), center(3),
     .			grx(NX), gry(NY), grz(NZ),
     .			msk(NX,NY,NZ), mskcntr(3)

c		Local variables

			integer i, j, k, l, m

			real	rmax, rmin, rr, th, ph, srr, hx

			real	intpsf

			hx = grx(2) - grx(1)


c		First thing is to find the rmax and rmin 

			rmin = 99999.0
			rmax = -99999.0
			do i = 1, NST
				do j = 1, NSP
						if(S(i,j,1) .gt. rmax)then
							rmax = S(i,j,1)
						end if
						if(S(i,j,1) .lt. rmin)then
							rmin = S(i,j,1)
						end if
				end do
			end do
	
			rmin = rmin - hx
			rmax = rmax + hx

			write(0,*)'min max -> ',rmin,rmax


c		Next locate points inside the surface and tag

			do i = 1, NX
				do j = 1, NY
					do k = 1, NZ

						rr = sqrt((grx(i) - center(1))**2 + 
     .		(gry(j) - center(2))**2 + (grz(k) - center(3))**2)

c			If outside of outer boundary of neighborhood of surface
						if( rr .ge. rmax)then
							msk(i,j,k) = 1.0
						end if

c			If inside of inner boundary of neighborhood of surface
c				also if rr = 0.0 relative to surface center
						if(rr .le. rmin .or. rr .eq. 0.0)then
							msk(i,j,k) = 0.0
						end if

c			If inside neighborhood of surface
						if( rr .le. rmax .and. rr .ge. rmin)then
c							Calculate theta and phi of point
							call calcsph(grx(i),gry(j), grz(k),th,ph,center)

c							Carry out linear interpolation to get surface r
							 srr = intpsf(th,ph,S,NST,NSP)

							if( rr .le. srr)then
								msk(i,j,k) = 0.0
							else
								msk(i,j,k) = 1.0
							end if

							
						end if

					end do
				end do
			end do

			mskcntr(1) = center(1)
			mskcntr(2) = center(2)
			mskcntr(3) = center(3)

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		subroutine calcsph(x,y,z,th,ph,c)
			implicit none
			
			real	x, y, z, th, ph, c(3), tr, Pi

			Pi = 4.0 * atan(1.0)

			tr = sqrt((x-c(1))**2 + (y-c(2))**2 + (z-c(3))**2)
c			if((x-c(1)) .eq. 0.0 .and. (y-c(2)) .eq. 0.0 .and. 
c     .			(z-c(3)) .gt. 0.0 )then
c				th = 0.0
c			else
c				if((x-c(1)) .eq. 0.0 .and. (y-c(2)) .eq. 0.0 .and. 
c     .			(z-c(3)) .lt. 0.0 )then
c					th = Pi
c				else
c					th = acos((z-c(3))/tr)
c				end if
c			end if

			th = acos((z-c(3))/(tr + 1.0e-7))


			if(ph .ge. 2.0*Pi)ph = ph - 2.0*Pi

			if(ph .lt. 0.0)ph = ph + 2.0*Pi

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		real function intpsf(th,ph,S,NST,NSP)
			implicit none

			integer NST, NSP
		
			real	th ,ph, S(NST,NSP,3)

c		Local variables

			integer	nt, np, tt(2),pp(2)

			real	dtheta, dphi,Pi, rt(2),d1,d2,pv(2)

			Pi = 4.0 * atan(1.0)

			dtheta = S(2,1,2) - S(1,1,2)
			dphi = S(1,2,3) - S(1,1,3)

			nt = int(th/dtheta) + 1
			np = int(ph/dphi) + 1
c		The following commented out lines cannot happen since th <=Pi
c			if(nt .eq. NST)then
c				tt(1) = NST - 1	
c				tt(2) = NST
c			else
c				tt(1) = nt
c				tt(2) = nt + 1
c			endif

			if(np .eq. NSP)then
				pv(1) = 2.0*Pi - dphi
				pv(2) = 2.0*Pi
				pp(1) = NSP
				pp(2) = 1
			else
				pv(1) = float(np - 1)* dphi
				pv(2) = float(np) * dphi
				pp(1) = np
				pp(2) = np + 1
			endif

			d1 = (th - S(nt,pp(1),2))
			d2 = (th - S(nt+1,pp(1),2))

			rt(1) = (-d1*S(nt,pp(1),1) + d2*S(nt+1,pp(1),1))/dtheta

			d1 = (th - S(nt,pp(2),2))
			d2 = (th - S(nt+1,pp(2),2))

			rt(2) = (-d1*S(nt,pp(2),1) + d2*S(nt+1,pp(2),1))/dtheta
			

			d1 = (ph - pv(1))
			d2 = (ph - pv(2))
		
			intpsf = (-d1*rt(1) + d2*rt(2))/dphi

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	Routine to update the mask on the next slice

		subroutine updmsk4(S,Sc,grx,gry,grz,msk,mskcntr,center,
     .				itmpl,NST,NSP,NX,NY,NZ)
			implicit none

			integer	NST, NSP, NX, NY, NZ
			integer  itmpl(NST,NSP,4,4,2)

			real	S(NST,NSP,3), Sc(NST,NSP,3), center(3),
     .			grx(NX), gry(NY), grz(NZ),
     .			msk(NX,NY,NZ), mskcntr(3)

c		Local variables

			integer i, j, k, l, m

			real	rmax, rmin, rr, th, ph, srr, hx

			real	intpsf4

			hx = grx(2) - grx(1)


c		First thing is to find the rmax and rmin 

			rmin = 99999.0
			rmax = -99999.0
			do i = 1, NST
				do j = 1, NSP
						if(S(i,j,1) .gt. rmax)then
							rmax = S(i,j,1)
						end if
						if(S(i,j,1) .lt. rmin)then
							rmin = S(i,j,1)
						end if
				end do
			end do
	
			rmin = rmin - hx
			rmax = rmax + hx

			write(0,*)'min max -> ',rmin,rmax


c		Next locate points inside the surface and tag

			do i = 1, NX
				do j = 1, NY
					do k = 1, NZ

						rr = sqrt((grx(i) - center(1))**2 + 
     .		(gry(j) - center(2))**2 + (grz(k) - center(3))**2)

c			If outside of outer boundary of neighborhood of surface
						if( rr .ge. rmax)then
							msk(i,j,k) = 1.0
						end if

c			If inside of inner boundary of neighborhood of surface
c				also if rr = 0.0 relative to surface center
						if(rr .le. rmin .or. rr .eq. 0.0)then
							msk(i,j,k) = 0.0
						end if

c			If inside neighborhood of surface
						if( rr .le. rmax .and. rr .ge. rmin)then
c							Calculate theta and phi of point
							call calcsph(grx(i),gry(j), grz(k),th,ph,center)

c							Carry out linear interpolation to get surface r
							 srr = intpsf4(th,ph,S,itmpl,NST,NSP)

							if( rr .le. srr)then
								msk(i,j,k) = 0.0
							else
								msk(i,j,k) = 1.0
							end if

							
						end if

					end do
				end do
			end do

			mskcntr(1) = center(1)
			mskcntr(2) = center(2)
			mskcntr(3) = center(3)

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		real function intpsf4(th,ph,S,itmpl,NST,NSP)
			implicit none

			integer NST, NSP, itmpl(NST,NSP,4,4,2)
		
			real	th ,ph, S(NST,NSP,3)

c		Local variables

			integer	nt, np, tt(2),pp(2), ii, jj

			real	dtheta, dphi,Pi, rt(4),d1,d2,pv(2),trem,prem,
     .			newt,newp,rtt(4),rpp(4), rff(4),tff(4),reci

			reci = 1.0 / 6.0

			rtt(1) = 1.0
			rtt(2) = 2.0
			rtt(3) = 3.0
			rtt(4) = 4.0

			rpp(1) = 1.0
			rpp(2) = 2.0
			rpp(3) = 3.0
			rpp(4) = 4.0

			Pi = 4.0 * atan(1.0)

			dtheta = S(2,1,2) - S(1,1,2)
			dphi = S(1,2,3) - S(1,1,3)

			nt = int(th/dtheta) + 1
			np = int(ph/dphi) + 1

		   trem = th - S(nt,np,2)	
		   prem = ph - S(nt,np,3)	

			newt = trem/dtheta + 2.0
			newp = prem/dphi + 2.0

         jj=1
         do ii=1,4
            rff(ii) = S(itmpl(nt,np,ii,jj,1),
     &                     itmpl(nt,np,ii,jj,2),1)
         end do

         tff(jj) =
     #  -(newt-rtt(2))*(newt - rtt(3))
     #   *(newt -rtt(4))*rff(1)*reci+
     #   (newt-rtt(1))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(2)*0.5e0-
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(4))*rff(3)*0.5e0+
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(3))*rff(4)*reci

                jj=2
                do ii=1,4
                rff(ii) = S(itmpl(nt,np,ii,jj,1),
     &                  itmpl(nt,np,ii,jj,2),1)
                end do
         tff(jj) =
     #  -(newt-rtt(2))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(1)*reci+
     #   (newt-rtt(1))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(2)*0.5e0-
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(4))*rff(3)*0.5e0+
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(3))*rff(4)*reci

                jj=3
                do ii=1,4
                rff(ii) = S(itmpl(nt,np,ii,jj,1),
     &                  itmpl(nt,np,ii,jj,2),1)
                end do

         tff(jj) =
     #  -(newt-rtt(2))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(1)*reci+
     #   (newt-rtt(1))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(2)*0.5e0-
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(4))*rff(3)*0.5e0+
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(3))*rff(4)*reci


         jj=4
         do ii=1,4

           rff(ii) = S(itmpl(nt,np,ii,jj,1),
     &               itmpl(nt,np,ii,jj,2),1)
         end do  

         tff(jj) =
     #  -(newt-rtt(2))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(1)*reci+
     #   (newt-rtt(1))*(newt - rtt(3))*
     #   (newt -rtt(4))*rff(2)*0.5e0-
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(4))*rff(3)*0.5e0+
     #   (newt-rtt(1))*(newt - rtt(2))*
     #   (newt -rtt(3))*rff(4)*reci

         intpsf4 =
     #  -(newp-rpp(2))*(newp - rpp(3))*
     #   (newp -rpp(4))*tff(1)*reci+
     #   (newp-rpp(1))*(newp - rpp(3))*
     #   (newp -rpp(4))*tff(2)*0.5e0-
     #   (newp-rpp(1))*(newp - rpp(2))*
     #   (newp -rpp(4))*tff(3)*0.5e0+
     #   (newp-rpp(1))*(newp - rpp(2))*
     #   (newp -rpp(3))*tff(4)*reci

			return
		end
c-----------------------------------------------------------------------
c	Routine to update the mask on the next slice

		subroutine updmsk2(S,Sc,grx,gry,grz,msk,mskcntr,center,
     .				NST,NSP,NX,NY,NZ)
			implicit none

			integer	NST, NSP, NX, NY, NZ

			real	S(NST,NSP,3), Sc(NST,NSP,3), center(3),
     .			grx(NX), gry(NY), grz(NZ),
     .			msk(NX,NY,NZ), mskcntr(3)

c		Local variables

			integer i, j, k, l, m, imin1,imin2,jmin1,jmin2

			real	rmax, rmin, rr, th, ph, srr, hx, rs,rsmin1,rsmin2,
     .			xy, xymax, xymin

			real	intpsf

			hx = grx(2) - grx(1)


c		First thing is to find the rmax and rmin 

			xymin = 999999.0
			xymax = -999999.0
         rmin = 99999.0
         rmax = -99999.0
			do i = 1, NST
				do j = 1, NSP
						if(S(i,j,1) .gt. rmax)then
							rmax = S(i,j,1)
						end if
						if(S(i,j,1) .lt. rmin)then
							rmin = S(i,j,1)
						end if
						xy = sqrt((Sc(i,j,1)-center(1))**2 + 
     .					(Sc(i,j,2) -center(2))**2)
						if(xy .gt. xymax)xymax = xy
						if(xy .lt. xymax)xymin = xy
				end do
			end do
	
c			xymin = xymin - hx
c			xymax = xymax + hx

			write(0,*)'XY MIN/MAX ',xymin, xymax


c		Initialize msk
			do i = 1, NX
				do j = 1, NY
					do k = 1, NZ
						rs = sqrt((grx(i) - center(1))**2 + 
     .		(gry(j) - center(2))**2 + (grz(k) - center(3))**2)
						if(rs .ge. rmin)then
							msk(i,j,k) = 1.0
						else
							msk(i,j,k) = 0.0
						end if
					end do
				end do
			end do

c		Refine msk based on nearest surface points

c		First find nearest intersection of z-ray with x,y sets
			rsmin1 = 99999.0
			rsmin2 = 99999.0
			do i = 1, NX
				do j = 1, NY
					xy = sqrt((grx(i)-center(i))**2 + (gry(j)-center(2))**2)
					if(xy .le. xymax .and. xy .ge. xymin)then

					do l = 1, (NST-1)/2	
						do m =1, NSP	
							rs = sqrt((Sc(l,m,1)-grx(i))**2 + 
     .							(Sc(l,m,2)-gry(j))**2)
							if(rs .lt. rsmin1)then
								rsmin1 = rs
								imin1 = l
								jmin1 = m
							end if
						end do
					end do
					do l = 1 + (NST-1)/2	, NST
						do m =1, NSP	
							rs = sqrt((Sc(l,m,1)-grx(i))**2 + 
     .							(Sc(l,m,2)-gry(j))**2)
							if(rs .lt. rsmin2)then
								rsmin2 = rs
								imin2 = l
								jmin2 = m
							end if
						end do
					end do
					
					do k = 1, NZ
						if( grz(k) .le. sc(imin2,jmin2,3) .or. 
     $					grz(k) .ge. sc(imin2,jmin2,3))then
							msk(i,j,k) = 1.0
						else 
							if(grz(k) .gt. sc(imin2,jmin2,3) .and.
     $						grz(k) .lt. sc(imin1,jmin1,3))then
							msk(i,j,k) = 0.0
							end if
						end if
						
					end do

					end if
				end do
			end do

			mskcntr(1) = center(1)
			mskcntr(2) = center(2)
			mskcntr(3) = center(3)

			return
		end
c-----------------------------------------------------------------------
c	Routine to update metric and extrinsic curvature data on points that
c  previously might not be in the computational domain and now are.
c
c	o Assume at most NX points to pop up at most. 
c  o
c
c	Input

		subroutine chkcmpdom(xx,yy,zz,oldmsk,msk,
     .	g11, g12, g13, g22, g23, g33, k11, k12, k13, k22, k23, k33,
     .	NX, NY, NZ)
			implicit none

			integer	NX, NY, NZ

			real	xx(NX), yy(NY), zz(NZ), 
     .		oldmsk(NX,NY,NZ), msk(NX,NY,NZ),
     .		g11(NX,NY,NZ), g12(NX,NY,NZ), g13(NX,NY,NZ),
     .		g22(NX,NY,NZ), g23(NX,NY,NZ), g33(NX,NY,NZ),
     .		k11(NX,NY,NZ), k12(NX,NY,NZ), k13(NX,NY,NZ),
     .		k22(NX,NY,NZ), k23(NX,NY,NZ), k33(NX,NY,NZ)

c		Local variables

			integer	i, j, k, num,
     .		ix(NX*NY*NZ), iy(NY*NY*NZ), iz(NZ*NY*NZ)
			
			real	lx(NX*NY*NZ), ly(NX*NY*NZ), lz(NX*NY*NZ)


c			Count number of new points
		   num = 0
			do i = 1, NX
			  do j = 1, NY
			    do k = 1, NZ	
c				If defined on new mask and undefined on old mask then...
					if(msk(i,j,k) .eq. 1.0 .and. 
     .			   oldmsk(i,j,k) .eq. 0.0)then
					   num = num + 1 
					end if
			    end do
			  end do
			end do
			if(num .eq. 0) goto 999
			if(num .gt. NX)then
				write(0,*)'ERROR in chkcmpdom: Number exceeds NX'
				write(0,*)'     Increase list size'
				stop
			end if
			write(6,*)'Found ',num,' points undefined in domain'
			write(6,*)'  Updating via interpolation'
			

c			Collect list of undefined points 

			call getlst(lx,ly,lz,ix,iy,iz,
     .		  xx,yy,zz,oldmsk,msk,NX,NY,NZ,num)

c			Interpolate to given list of points using old mask

			call intpfncs(lx,ly,lz,xx,yy,zz,oldmsk,
     .   g11, g12, g13, g22, g23, g33, k11, k12, k13, k22, k23, k33,
     .	NX, NY, NZ, num)

999		continue

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	Subsidiary routine for chkcmpdom

c	Constructs list of points that have popped out

		subroutine getlst(lx,ly,lz,ix,iy,iz,
     .	xx,yy,zz,oldmsk,msk,NX,NY,NZ,NN)
			implicit none

			integer NX, NY, NZ, NN,
     .				ix(NN), iy(NN), iz(NN)
			
			real	lx(NN), ly(NN), lz(NN),
     .			xx(NX), yy(NY), zz(NZ),
     .			oldmsk(NX,NY,NZ), msk(NX,NY,NZ)

c		Local variables

			integer	i, j, k, p

			p = 0
			do i = 1, NX
			  do j = 1, NY
				do k = 1, NZ
               if(msk(i,j,k) .eq. 1.0 .and.
     .            oldmsk(i,j,k) .eq. 0.0)then
                  p = p + 1
c				Store coordinates of point
						lx(p) = xx(i)
						ly(p) = yy(j)
						lz(p) = zz(k)
c				Store grid indices for point
						ix(p) = i
						iy(p) = j
						iz(p) = k
               end if
			   end do
			  end do
			end do

			return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	Subsidiary routine for chkcmpdom
c	Interpolates using Extrapolation routines developed by Scott.

		subroutine intpfncs(lx,ly,lz,ix,iy,iz,xx,yy,zz,msk,
     .   g11, g12, g13, g22, g23, g33, k11, k12, k13, k22, k23, k33,
     .   NX, NY, NZ, NN)
			implicit none

			integer	NX, NY, NZ, NN,
     .            ix(NN), iy(NN), iz(NN)

			real	lx(NN), ly(NN), lz(NN), 
     .			xx(NX), yy(NY), zz(NZ), msk(NX,NY,NZ),
     .		g11(NX,NY,NZ), g12(NX,NY,NZ), g13(NX,NY,NZ),
     .		g22(NX,NY,NZ), g23(NX,NY,NZ), g33(NX,NY,NZ),
     .		k11(NX,NY,NZ), k12(NX,NY,NZ), k13(NX,NY,NZ),
     .		k22(NX,NY,NZ), k23(NX,NY,NZ), k33(NX,NY,NZ)

c		Local Variables	

			integer	i, j, k, ierr, d4lintxyzcd, d4lintxyzc,
     .		d2lintxyzcd, d2lintxyzc

			integer iwk(20*NN),swk(5*NN)

			integer diditf

			real	tfunc(NN)

			diditf = 0



c		Evaluate g11


			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .              lx,ly,lz,g11,tfunc,msk,iwk,swk,diditf,' ')
			
			do k = 1, NN
			  g11(ix(k), iy(k), iz(k)) = tfunc(k)
			end do

			
c		Evaluate g12
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,g12,tfunc,msk,iwk,swk,diditf,' ')
			
			do k = 1, NN
			  g12(ix(k), iy(k), iz(k)) = tfunc(k)
			end do


c		Evaluate g13
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,g13,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  g13(ix(k), iy(k), iz(k)) = tfunc(k)
			end do
			
c		Evaluate g22
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,g22,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  g22(ix(k), iy(k), iz(k)) = tfunc(k)
			end do


c		Evaluate g23
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,g23,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  g23(ix(k), iy(k), iz(k)) = tfunc(k)
			end do
			
c		Evaluate g33
			ierr= d4lintxyzcd( NX,NY,NZ,xx,xx,xx,NN,
     .         lx,ly,lz,g33,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  g33(ix(k), iy(k), iz(k)) = tfunc(k)
			end do

c		Next the Extrinsic curvature components

c		Evaluate k11

c		Evaluate k11


			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .              lx,ly,lz,k11,tfunc,msk,iwk,swk,diditf,' ')
			
			do k = 1, NN
			  k11(ix(k), iy(k), iz(k)) = tfunc(k)
			end do

			
c		Evaluate k12
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,k12,tfunc,msk,iwk,swk,diditf,' ')
			
			do k = 1, NN
			  k12(ix(k), iy(k), iz(k)) = tfunc(k)
			end do


c		Evaluate k13
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,k13,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  k13(ix(k), iy(k), iz(k)) = tfunc(k)
			end do
			
c		Evaluate k22
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,k22,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  k22(ix(k), iy(k), iz(k)) = tfunc(k)
			end do


c		Evaluate k23
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,k23,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  k23(ix(k), iy(k), iz(k)) = tfunc(k)
			end do
			
c		Evaluate k33
			ierr= d4lintxyzcd( NX,NY,NZ,xx,yy,zz,NN,
     .         lx,ly,lz,k33,tfunc,msk,iwk,swk,diditf,' ')

			do k = 1, NN
			  k33(ix(k), iy(k), iz(k)) = tfunc(k)
			end do

			return
		end
c-----------------------------------------------------------------------