program GridGen2MHD
c      
c   2nd part of Grid generation for MHD calculations.
c   run GridGen1 first then GridGen2MHD.
c
c   grid generation of overlapping grids for a shpheroid
c   it consists of 2 programs: GridGen1 and GridGen2
c   GridGen1: generation of surface grids covering a
c             spheroidal surface, output: fort.7 
c   GridGen2MHD: generation of grids for the simulations
c             input file: fort.7 (from GridGen1)
c             ouput files: fort.71-fort.73, fort.81-fort.96.
c  


c 2 spheres + cube
      parameter (mx1=51,my1=51,mz1=51, mx2=71, my2=71, mz2=61+1)
      parameter (mx3=71, my3=71, mz3=61+1)
c for grid 5 and 5b and 5 ex
      parameter (nphi=64*2, nmu=nphi, nnu=nphi/2, nmuh=nmu/2)
      dimension amu(nmuh+2), anu(nnu+2), phi(nphi+2)
      dimension ig5(nphi+2, nmuh+2, nnu+2), ngp5(nphi+2, nmuh+2,nnu+2,3)
      dimension drngp5(nphi+2, nmuh+2, nnu+2, 3)
      dimension ig5b(nphi+2,nnu+2),ngp5b(nphi+2, nnu+2,3), 
     1 drngp5b(nphi+2,nnu+2,3)
c for 5 ex and Grid 1, 2, 3
      dimension ngp15(mx1,my1,mz1,3), drngp15(mx1,my1,mz1,3)      
      dimension ngp25(mx2,my2,mz2,3), drngp25(mx2,my2,mz2,3)      
      dimension ngp35(mx3,my3,mz3,3), drngp35(mx3,my3,mz3,3)    
      dimension ngp25b(mx2,my2,2), drngp25b(mx2,my2,2)      
      dimension ngp35b(mx3,my3,2), drngp35b(mx3,my3,2)    

c      
	dimension g1(mx1,my1,mz1,3), ig1(mx1,my1,mz1), 
     1	r1(mx1), s1(my1), t1(mz1) 
	dimension g2(mx2,my2,mz2,3), ig2(mx2,my2,mz2), 
     1 iig2(mx2,my2),
     1 r2(mx2), s2(my2), t2(mz2), trans(mx2,my2,mz2,3,3) 
	dimension g3(mx2,my2,mz2,3), ig3(mx2,my2,mz2), 
     1 iig3(mx3,my3), 
     1	r3(mx3), s3(my3), t3(mz3) 
	integer mth, mphi,mwsave
	parameter(mth=3*22,mphi=3*3*2**4,mwsave=2*mphi+15)
	parameter(imax=mphi, mmax=imax/3*2, jmax=mth,lmax=mmax,
     1 lmmax=lmax+mmax,lmmax1=lmmax+1)
      parameter (mx4=mth,my4=mphi)
	dimension g4(mx4,my4,2), ig4(mx4,my4), gausspt(mx4), gausswt(mx4) 

      common/parm/factor
c
c 1.
c grid 2 spheroidal shell coordinate 
c x=a sinphi costh, y= a sinphi sinth, z=b cosphi 
c npole=1; north pole
      radmax=1.
	factor=0.8
c     factor=1.
      radmin=0.4*factor

	nx2=mx2
	ny2=my2
	nz2=mz2
      npole=1
	call sphere(npole, g2, r2, s2, t2)

c grid 3 spheroidal shell coordinate 
c x=a sinphi costh, y= a sinphi sinth, z=b cosphi 
c npole=1; south pole

	nx3=mx3
	ny3=my3
	nz3=mz3
      npole=-1
	call sphere(npole, g3, r3, s3, t3)

c cubic Grid 1

	nx1=mx1
	ny1=my1
	nz1=mz1
	call cube(g1, r1, s1, t1)
	do i=1,nx1
	do j=1,ny1
	do k=1,nz1
	if(g1(i,j,k,3).ge.0.) then
	ig1(i,j,k)=2
	else
	ig1(i,j,k)=3
	end if
	end do
	end do
	end do
c
c 2. 
c 2.1 mark interior point and interpolating point on G2
      read(7) iig2, iig3
	do k=1,nz2
	do i=1,nx2
	do j=1,ny2
	ig2(i,j,k)=iig2(i,j)
c     ig2(i,j,k)=2
	if(k.eq.1 .and. ig2(i,j,k).eq.-3) ig2(i,j,k)=-1
	if(k.eq.1 .and. ig2(i,j,k).eq.2) ig2(i,j,k)=-1
	end do
	end do
	end do
c do the same for G3
c   
	do k=1,nz3
	do i=1,nx3
	do j=1,ny3
	ig3(i,j,k)=iig3(i,j)
c	ig3(i,j,k)=3
	if(k.eq.1 .and. ig3(i,j,k).eq.-2) ig3(i,j,k)=-1
	if(k.eq.1 .and. ig3(i,j,k).eq.3) ig3(i,j,k)=-1
	end do
	end do
	end do
	

c 3. mark interior point and interpolating point on G1
c    mark interpolating point
	do i=1,nx1
	do j=1,ny1
      do k=1,nz1
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
      r=sqrt(x*x+y*y+z*z)
      if (r.le.radmin) ig1(i,j,k)=1
      if (r.ge.radmax) ig1(i,j,k)=0
      end do
      end do
      end do
      
c now check if ig1(i,j,k)=2 points can be interpolated from G2
	do i=1,nx1
	do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.2) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
      npole=1
      call sphngp(npole, x, y, z, g2, r2, s2, t2, nnx, nny, nnz)
      if( nnz.lt.2 ) ig1(i,j,k)=1
	if( nnz.gt.nz2-1 ) ig1(i,j,k)=0
	end if
	end do
	end do
	end do

	do i=1,nx1
	do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.3) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
      npole=-1
      call sphngp(npole, x, y, z, g3, r3, s3, t3, nnx, nny, nnz)
      if( nnz.lt.2 ) ig1(i,j,k)=1
	if( nnz.gt.nz2-1 ) ig1(i,j,k)=0
	end if
	end do
	end do
	end do


c 4.1 mark interpolating points on G1 which are needed 
c for the boundary points on G2 and G3
      k=1
      do i=1,nx2
      do j=1,ny2
      x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c in G1 grid
     
      call cubengp(x, y, z, g1, r1, s1, t1, nnx, nny, nnz)
c mark the points on G1 that are required for this point on G2

  777 format(3i4,3f10.4,3i4)
      
      
      do k1=1,5
      kk1=k1-3
      do k2=1,5
      kk2=k2-3
      do k3=1,5
      kk3=k3-3

      ig1(nnx+kk1,nny+kk2,nnz+kk3)=
     1        -abs(ig1(nnx+kk1,nny+kk2,nnz+kk3))
      end do
	end do
      end do
      end do
	end do


      k=1
      do i=1,nx3
      do j=1,ny3
      x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c in G1 grid
      call cubengp(x, y, z, g1, r1, s1, t1, nnx, nny, nnz)
      
c mark the points on G1 that are required for this point on G3

      do k1=1,5
      kk1=k1-3
      do k2=1,5
      kk2=k2-3
      do k3=1,5
      kk3=k3-3

      ig1(nnx+kk1,nny+kk2,nnz+kk3)=
     1        -abs(ig1(nnx+kk1,nny+kk2,nnz+kk3))
      end do
	end do
      end do
	end do
      end do


c 5. delete interpolation points on G1 which are not needed
c 5.1 the interpolating points not needed by the interior points
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.2 .or. ig1(i,j,k).eq.3) then
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
	ii=i+kk1
	if(ii.eq.0) ii=i
	if(ii.eq.nx1+1) ii=nx1
	jj=j+kk2
	if(jj.eq.0) jj=j
	if(jj.eq.ny1+1) jj=ny1
	kk=k+kk3
	if(kk.eq.0) kk=k
	if(kk.eq.nz1+1) kk=nz1

c if ig1=1 it means the point is required by the interior point
      if(abs(ig1(ii,jj,kk)).eq.1) go to 100 
      end do
      end do
      end do
      ig1(i,j,k)=0
  100 continue
      end if
      end do
      end do
      end do

c 5.2 the interpolating points can be changed to interior points
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(abs(ig1(i,j,k)).eq.2 .or. abs(ig1(i,j,k)).eq.3) then
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
	ii=i+kk1
	if(ii.eq.0) ii=i
	if(ii.eq.nx1+1) ii=nx1
	jj=j+kk2
	if(jj.eq.0) jj=j
	if(jj.eq.ny1+1) jj=ny1
	kk=k+kk3
	if(kk.eq.0) kk=k
	if(kk.eq.nz1+1) kk=nz1
c if ig1=0 it means that the point cannot be an interior point
      if(abs(ig1(ii,jj,kk)).eq.0) go to 200 
      end do
      end do
      end do
      ig1(i,j,k)=1
  200 continue
      end if
      end do
      end do
      end do


c 5.3 on G2 mark the points reguired by G1 interpolating points
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(abs(ig1(i,j,k)).eq.2) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G2
      npole=1
      call sphngp(npole, x, y, z, g2, r2, s2, t2, nnx, nny, nnz)

c nnx, nny, nnz on G2 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-3
      do k3=1,3
      kk3=k3-3
	ii=nnx+kk1
	jj=nny+kk2
	kk=nnz+kk3
	if(ii.eq.0 .or. ii.eq.nx2+1) then
       if(jj.eq.0 .or. jj.eq.ny2+1) then
	 if(kk.eq.0 .or. kk.eq.nz2+1) then
	 write(*,*) 'out of bounds in Sec 5.2 for G2'
	 end if
	 end if
	end if
      ig2(nnx+kk1,nny+kk2,nnz+kk3)=-ig2(nnx+kk1,nny+kk2,nnz+kk3) 
      end do
      end do
      end do
	end if
      end do
      end do
      end do

      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(abs(ig1(i,j,k)).eq.3) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G3
      npole=-1
      call sphngp(npole, x, y, z, g3, r3, s3, t3, nnx, nny, nnz)
c nnx, nny, nnz on G3 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-3
      do k3=1,3
      kk3=k3-3
	ii=nnx+kk1
	jj=nny+kk2
	kk=nnz+kk3
	if(ii.eq.0 .or. ii.eq.nx3+1) then
       if(jj.eq.0 .or. jj.eq.ny3+1) then
	 if(kk.eq.0 .or. kk.eq.nz3+1) then
	 write(*,*) 'out of bounds in Sec 5.2 for G3'
	 end if
	 end if
	end if

      ig3(nnx+kk1,nny+kk2,nnz+kk3)=-ig3(nnx+kk1,nny+kk2,nnz+kk3) 
      end do
      end do
      end do
	end if
      end do
      end do
      end do

c 6.0 
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(abs(ig1(i,j,k)).eq.1) ig1(i,j,k)=1
      if(abs(ig1(i,j,k)).ne.1) ig1(i,j,k)=-abs(ig1(i,j,k))
      end do
	end do
	end do

      do i=1,nx2
      do j=1,ny2
      do k=1,nz2
      if(abs(ig2(i,j,k)).eq.2) ig2(i,j,k)=2
      if(abs(ig2(i,j,k)).ne.2) ig2(i,j,k)=-abs(ig2(i,j,k))
      end do
	end do
	end do

      do i=1,nx3
      do j=1,ny3
      do k=1,nz3
      if(abs(ig3(i,j,k)).eq.3) ig3(i,j,k)=3
      if(abs(ig3(i,j,k)).ne.3) ig3(i,j,k)=-abs(ig3(i,j,k))
      end do
	end do
	end do            
       
c      write(*,*) 'ig1(1,17,21) before 6a', ig1(1,17,21)

c check the grid
c (a) check if G1 interpolating points depend on G2 and G3 
c     interpolating points
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.-2) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G2
      npole=1
      call sphngp(npole, x, y, z, g2, r2, s2, t2, nnx, nny, nnz)
      if(i.eq.1 .and. j.eq.17 .and. k.eq.21) then
	tmp=0.
	end if
	if(nnx.le.1 .or. nnx.ge.(nx2-1) .or.
     1  nny.le.1 .or. nny.ge.(ny2-1) .or.
	2  nnz.le.1 .or. nnz.ge.(nz2-1)) then
      write(*,*) 'out of bounds Sec. 6a', i, j, k, nnx, nny, nnz, 
     1	ig1(i,j,k)
	end if 

c nnx, nny, nnz on G2 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
cc
cc
      if(abs(ig2(nnx+kk1,nny+kk2,nnz+kk3)).eq.1) then
	 write(*,*) 'G2 overlap'
	end if 
      end do
      end do
      end do
	end if

      if(ig1(i,j,k).eq.-3) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G2
      npole=-1
      call sphngp(npole, x, y, z, g3, r3, s3, t3, nnx, nny, nnz)

c nnx, nny, nnz on G2 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      if(abs(ig3(nnx+kk1,nny+kk2,nnz+kk3)).eq.1) then
	write(*,*) 'G3 overlap'
	end if 
      end do
      end do
      end do
	end if
      end do
      end do
      end do
c (b) check if G2 (G3) interpolating points depend on G1
c     interpolating points
      k=1
      do i=1,nx2
      do j=1,ny2
      x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c in G1 grid   
      call cubengp(x, y, z, g1, r1, s1, t1, nnx, nny, nnz)
	if(nnx.le.1 .or. nnx.ge.(nx1-1) .or.
     1  nny.le.1 .or. nny.ge.(ny1-1) .or.
	2  nnz.le.1 .or. nnz.ge.(nz1-1)) write(*,*) 'out of bounds Sec. 6b'

      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      
      if(ig1(nnx+kk1,nny+kk2,nnz+kk3) .le. 0) then
      write(*,*) 'G1 (from G2) overlap'
	end if
      end do
	end do
      end do
      end do
	end do

      k=1
      do i=1,nx3
      do j=1,ny3
      x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c in G1 grid   
      call cubengp(x, y, z, g1, r1, s1, t1, nnx, nny, nnz)
c mark the points on G1 that are required for this point on G3     
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      if(ig1(nnx+kk1,nny+kk2,nnz+kk3) .le. 0) then
      write(*,*) 'G1 (from G3) overlap'
	end if
      end do
	end do
      end do
      end do
	end do

c (c)
c (c) check if G2 with ig2=-3 interpolating points depend on G3 
c     interpolating points
      do i=1,nx2
      do j=1,ny2
      do k=1,nz2
      if(ig2(i,j,k).eq.-3) then
	x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c find the nearest grid point on G3
      npole=-1
      call sphngp(npole, x, y, z, g3, r3, s3, t3, nnx, nny, nnz)
      if(nnx.le.1 .or. nnx.ge.(nx3-1) .or.
     1  nny.le.1 .or. nny.ge.(ny3-1) .or.
	2  nnz.le.1 .or. nnz.gt.(nz3)) then
      write(*,*) 'out of bounds Sec. 6c', nnx, nny, nnz
	end if

c nnx, nny, nnz on G3 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      kk3=0
cc
cc
      if(ig3(nnx+kk1,nny+kk2,nnz+kk3).le.0) then
	 write(*,*) 'G2 overlap in Section 6(c)'
	end if 

      end do
      end do

	end if

      end do
      end do
      end do


c (d) check if G3 with ig3=-2 interpolating points depend on G2 
c     interpolating points
      do i=1,nx3
      do j=1,ny3
      do k=1,nz3
      if(ig3(i,j,k).eq.-2) then
	x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c find the nearest grid point on G2
      npole=1
      call sphngp(npole, x, y, z, g2, r2, s2, t2, nnx, nny, nnz)
      if(nnx.le.1 .or. nnx.ge.(nx2-1) .or.
     1  nny.le.1 .or. nny.ge.(ny2-1) .or.
	2  nnz.le.1 .or. nnz.gt.(nz2)) then
      write(*,*) 'out of bounds Sec. 6d', nnx, nny, nnz
	end if

c nnx, nny, nnz on G2 is the nearest grid point
      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
c      do k3=1,3
c      kk3=k3-2
      kk3=0
cc
cc
      if(ig2(nnx+kk1,nny+kk2,nnz+kk3).le.0) then
	 write(*,*) 'G3 overlap in Section 6(d)'
	end if 
      end do
      end do

	end if

      end do
      end do
      end do

c 8. check consistency
c 8.1 Grid 1
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).ne.1 .and. ig1(i,j,k).ne.-2
     1 .and. ig1(i,j,k).ne.-3 .and. ig1(i,j,k).ne.0 ) then
	 write(*,*) 'G1 problem: wrong in ig1 Section 8.1'
      end if

      if(ig1(i,j,k).eq.1) then

      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      kk3=0
      if(ig1(i+kk1,j+kk2,k+kk3).eq.0) then
	 write(*,*) 'G1 problem: interior points Section 8.1'
	end if 
      end do
      end do
      end do

      end if
      end do
      end do
      end do

c 8.2 Grid 2
      do i=1,nx2
      do j=1,ny2
      do k=1,nz2
      if(ig2(i,j,k).ne.2 .and. ig2(i,j,k).ne.-1
     1 .and. ig2(i,j,k).ne.-3 .and. ig2(i,j,k).ne.0 ) then
	 write(*,*) 'G2 problem: wrong in ig2 Section 8.2'
      end if

      if(ig2(i,j,k).eq.2) then

      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      kk3=0
      if(ig2(i+kk1,j+kk2,k+kk3).eq.0) then
	 write(*,*) 'G2 problem: interior points Section 8.2'
	end if 
      end do
      end do
      end do

      end if
      end do
      end do
      end do


c 8.3 Grid 3
      do i=1,nx3
      do j=1,ny3
      do k=1,nz3
      if(ig3(i,j,k).ne.3 .and. ig3(i,j,k).ne.-1
     1 .and. ig3(i,j,k).ne.-2 .and. ig3(i,j,k).ne.0 ) then
	 write(*,*) 'G3 problem: wrong in ig3 Section 8.3'
      end if

      if(ig3(i,j,k).eq.3) then

      do k1=1,3
      kk1=k1-2
      do k2=1,3
      kk2=k2-2
      do k3=1,3
      kk3=k3-2
      kk3=0
      if(ig3(i+kk1,j+kk2,k+kk3).eq.0) then
	 write(*,*) 'G3 problem: interior points Section 8.3'
	end if 
      end do
      end do
      end do

      end if
      end do
      end do
      end do




c done with the grid generation


c 10. output
c 10.1 G1
      write(81) ig1, g1
 
     
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.-2) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G2
      npole=1
      call sphngp31(npole, x,y,z, g2, r2, s2, t2, nnx, nny, nnz,
     1 dr, ds, dt)
      if(ig2(nnx,nny,nnz).ne.2) then
      write(*,*) 'in 10.1 not an interior point'
      stop
      end if

c nnx, nny, nnz is the nearest grid point on G2
c with displacement dr, ds, dt
      write(82, 821) nnx, nny, nnz, dr, ds, dt
  821 format(3i5,3e20.12)
	end if 
      end do
      end do
      end do
      nnx=9999
      write(82, 821) nnx, nny, nnz, dr, ds, dt
c close 82
c      
     
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.-3) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G3
      npole=-1
      call sphngp31(npole, x,y,z, g3, r3, s3, t3, nnx, nny, nnz,
     1 dr, ds, dt)
      if(ig3(nnx,nny,nnz).ne.3) then
      write(*,*) 'in 10.1 not an interior point'
      stop
      end if


c nnx, nny, nnz is the nearest grid point on G3
c with displacement dr, ds, dt
      write(83, 821) nnx, nny, nnz, dr, ds, dt
	end if 
      end do
      end do
      end do
      nnx=9999
      write(83, 821) nnx, nny, nnz, dr, ds, dt
c close 83
     
c 10.2 G2
      npole=1
      call transform(npole, trans, g2, r2, s2, t2)
c      stop 
      write(84) ig2, g2, trans

      k=1
      do i=1,nx2
      do j=1,ny2
      if(ig2(i,j,k).eq.-1) then
      x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
      
c in G1 grid   
      call cubengp3(x, y, z, g1, r1, s1, t1, nnx, nny, nnz, dr,ds,dt)
      write(85, 821) nnx, nny, nnz, dr, ds, dt

      if(ig1(nnx,nny,nnz).ne.1) then
      write(*,*) 'in 10.1 G2 in G1: G1 not an interior point'
      stop
      end if


	end if
      end do
	end do
      nnx=9999
      write(85, 821) nnx, nny, nnz, dr, ds, dt
     
c      dtmax=0.
       
      do i=1,nx2
      do j=1,ny2
c 
c      do k=1,nz2-1
      do k=1,nz2
      if(ig2(i,j,k).eq.-3) then
	x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c find the nearest grid point on G3
      npole=-1
      call sphngp32(npole, x,y,z, g3, r3, s3, t3, nnx, nny, nnz,
     1 dr, ds, dt)
c      dtmax=max(dtmax,abs(dt))
c nnx, nny, nnz is the nearest grid point on G3
c with displacement dr, ds, dt
      write(86, 821) nnx, nny, nnz, dr, ds, dt

      if(ig3(nnx,nny,nnz).ne.3) then
      write(*,*) 'in 10.1 G2 in G3: G3 not an interior point'
      stop
      end if

	end if 
      end do
      end do
      end do
      nnx=9999
      write(86, 821) nnx, nny, nnz, dr, ds, dt
c close 86

c 10.3 G3
      npole=-1
      call transform(npole, trans, g3, r3, s3, t3) 

      write(87) ig3, g3, trans

      k=1
      do i=1,nx3
      do j=1,ny3
      if(ig3(i,j,k).eq.-1) then
      x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
      
c in G1 grid   
      call cubengp3(x, y, z, g1, r1, s1, t1, nnx, nny, nnz, dr,ds,dt)
      write(88, 821) nnx, nny, nnz, dr, ds, dt

      if(ig1(nnx,nny,nnz).ne.1) then
      write(*,*) 'in 10.1 G3 in G1: G1 not an interior point'
      stop
      end if

	end if
      end do
	end do
      nnx=9999
      write(88, 821) nnx, nny, nnz, dr, ds, dt
c close 88
c      dtmax=0.
      do i=1,nx3
      do j=1,ny3
c
c      do k=1,nz3-1
      do k=1,nz3
      if(ig3(i,j,k).eq.-2) then
	x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c find the nearest grid point on G2
      npole= 1
      call sphngp32(npole, x,y,z, g2, r2, s2, t2, nnx, nny, nnz,
     1 dr, ds, dt)
c      dtmax=max(dtmax,abs(dt))
c nnx, nny, nnz is the nearest grid point on G2
c with displacement dr, ds, dt
      write(89, 821) nnx, nny, nnz, dr, ds, dt
      
      if(ig2(nnx,nny,nnz).ne.2) then
      write(*,*) 'in 10.1 G3 in G2: G2 not an interior point'
      stop
      end if
      
	end if 
      end do
      end do
      end do
      nnx=9999
      write(89, 821) nnx, nny, nnz, dr, ds, dt

c
c 11. G4 spherical shell
c find xi0 and dxi:
c      xi0=atanh(factor)
      xi0=1.0986122886

	xmax=cosh(xi0)
c 1-cosh(xi0-dxi)/xmax=da (da=grid spacing along x in the sphere routine)
c      dxi=xi0-acosh((1.-da)*xmax)
      dxi=0.0126	

c 11.1 define Gauss points
      call gquad(mx4,gausspt,gausswt)
      dphi=2.*3.14159265359/float(my4)
	nx4=mx4
	ny4=my4
      do i=1,nx4
      do j=1,ny4
      g4(i,j,1)=gausspt(i)
      g4(i,j,2)=dphi*float(j-1)
      end do
      end do


c 11.2 find ngp on G2 and G3 for points on G4 at R-dr
      
      xi=xi0-dxi
	pi=3.14159265359
	do i=1,nx4
      do j=1,ny4
      th=g4(i,j,1)
      pphi=g4(i,j,2)
	sinth=sin(th)
	x=cosh(xi)/xmax*sinth*cos(pphi)
	y=cosh(xi)/xmax*sinth*sin(pphi)
	z=sinh(xi)/xmax*cos(th)
	if(z.ge.0.) then
c from g2
      npole=1
      call sphngp31(npole,x,y,z,g2,r2,s2,t2,nnx,nny,nnz,
     1 dr, ds, dt)
c      call sphngp3(npole,th,pphi,r,g2,r2,s2,t2,nnx,nny,nnz,
c     1 dr, ds, dt)

      if(ig2(nnx,nny,nnz).eq.2) then
      ig4(i,j)=-2
c nnx, nny, nnz is the nearest grid point on G2
c with displacement dr, ds, dt
      write(90, 821) nnx, nny, nnz, dr, ds, dt
	end if 
	end if
      end do
      end do
      nnx=9999
      write(90, 821) nnx, nny, nnz, dr, ds, dt
c close 90
c
      do i=1,nx4
      do j=1,ny4
      th=g4(i,j,1)
      pphi=g4(i,j,2)
	sinth=sin(th)
	x=cosh(xi)/xmax*sinth*cos(pphi)
	y=cosh(xi)/xmax*sinth*sin(pphi)
	z=sinh(xi)/xmax*cos(th)
	if(z.lt.0.) then
c from g3
      npole=-1
      call sphngp31(npole,x,y,z,g3,r3,s3,t3,nnx,nny,nnz,
     1 dr, ds, dt)
c      call sphngp3(npole,th,pphi,r,g3,r3,s3,t3,nnx,nny,nnz,
c     1 dr, ds, dt)
      if(ig3(nnx,nny,nnz).eq.3) then
      ig4(i,j)=-3
c nnx, nny, nnz is the nearest grid point on G3
c with displacement dr, ds, dt
      write(91, 821) nnx, nny, nnz, dr, ds, dt
	end if 
	end if
      end do
      end do
      nnx=9999
      write(91, 821) nnx, nny, nnz, dr, ds, dt
c close 91
c
c 11.3 find ngp on G4 for points on G2 at R
      k=nz2-1
      do i=1,nx2
      do j=1,ny2
      if(ig2(i,j,k).eq.2) then
	x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c find the nearest grid point on G3
      call sphngp4(x,y,z,g4,gausspt,nnx,nny,
     1 dr, ds)

c nnx, nny, nnz is the nearest grid point on G4
c with displacement dr, ds, dt
      write(92, 821) nnx, nny, nnz, dr, ds, dt
	end if 
    
	
      end do
      end do
      nnx=9999
      write(92, 821) nnx, nny, nnz, dr, ds, dt
c close 92
c 11.4 find ngp on G4 for points on G3 at R
      k=nz3-1
      do i=1,nx3
      do j=1,ny3
      if(ig3(i,j,k).eq.3) then
	x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c find the nearest grid point on G3
      call sphngp4(x,y,z,g4,gausspt,nnx,nny,
     1 dr, ds)
c nnx, nny, nnz is the nearest grid point on G4
c with displacement dr, ds, dt
      write(93, 821) nnx, nny, nnz, dr, ds, dt
	end if 
      end do
      end do
      nnx=9999
      write(93, 821) nnx, nny, nnz, dr, ds, dt
c close 93
c
      write(94) xi0,dxi,ig4, g4


c


c 11.5 find ngp on G2 and G3 for points on G4 at R
      
      xi=xi0
	pi=3.14159265359
	do i=1,nx4
      do j=1,ny4
      th=g4(i,j,1)
      pphi=g4(i,j,2)
	sinth=sin(th)
	x=cosh(xi)/xmax*sinth*cos(pphi)
	y=cosh(xi)/xmax*sinth*sin(pphi)
	z=sinh(xi)/xmax*cos(th)
	if(z.ge.0.) then
c from g2
      npole=1
      call sphngp32(npole,x,y,z,g2,r2,s2,t2,nnx,nny,nnz,
     1 dr, ds, dt)
c      call sphngp3(npole,th,pphi,r,g2,r2,s2,t2,nnx,nny,nnz,
c     1 dr, ds, dt)

      if(ig2(nnx,nny,nnz).eq.2) then
      ig4(i,j)=-2
c nnx, nny, nnz is the nearest grid point on G2
c with displacement dr, ds, dt
      write(95, 821) nnx, nny, nnz, dr, ds, dt
	end if 
	end if
      end do
      end do
      nnx=9999
      write(95, 821) nnx, nny, nnz, dr, ds, dt
c close 95
c
      do i=1,nx4
      do j=1,ny4
      th=g4(i,j,1)
      pphi=g4(i,j,2)
	sinth=sin(th)
	x=cosh(xi)/xmax*sinth*cos(pphi)
	y=cosh(xi)/xmax*sinth*sin(pphi)
	z=sinh(xi)/xmax*cos(th)
	if(z.lt.0.) then
c from g3
      npole=-1
      call sphngp32(npole,x,y,z,g3,r3,s3,t3,nnx,nny,nnz,
     1 dr, ds, dt)
c      call sphngp3(npole,th,pphi,r,g3,r3,s3,t3,nnx,nny,nnz,
c     1 dr, ds, dt)
      if(ig3(nnx,nny,nnz).eq.3) then
      ig4(i,j)=-3
c nnx, nny, nnz is the nearest grid point on G3
c with displacement dr, ds, dt
      write(96, 821) nnx, nny, nnz, dr, ds, dt
	end if 
	end if
      end do
      end do
      nnx=9999
      write(96, 821) nnx, nny, nnz, dr, ds, dt
c close 96
c

c generating grid 5 and 5b
c Grid 5 and 5b
c  dimension phi(nphi+2), amu(nmu+2), anu(nnu+2)
c     define spheroid
      cc=0.8
      amu0=0.5*(alog(1.+cc)-alog(1.-cc))
      aa2=1.-cc*cc
      aa=sqrt(aa2)
      write(*,*) 'c, a, mu0=', cc, aa, amu0
c

c	nphi = 64*2
c	nmu = nphi
c      nnu = nphi/2 
c      nmuh = nmu/2

* Computational domain: [xleft,xright] x [yleft,yright]
c here x=th; y=r
        pi=4.*atan(1.)
        philength=2.*pi
        amulength=2.*amu0
        anulength=pi

        dphi = philength / float(nphi)
        
        do i = 1, nphi+2
          phi(i) = (i-2) * dphi
        enddo

c for Neumann bndry
        dmu = amulength / float(nmu)       
c        do j = 1, nmu
c          amu(j) = -amu0 + (float(j)-0.5) * dmu
c        enddo
c for mu > 0
        do j = 1, nmuh+2
          amu(j) =  (float(j)-1.5) * dmu
        enddo


        dnu = anulength / float(nnu)       
        do k = 1, nnu+2
          anu(k) = (float(k)-1.5) * dnu
        enddo

c for grid 5b (at amu=amu0)
c find ngp5
c (a) for points at g5 -> (x,y,z) -> is it in G1? yes find ngp5 and drngp5
c (b) is it in G2
c (c) is it in G3

c      icount=0
      do i = 2, nphi+1
      do j = 2, nmuh+1
      do k = 2, nnu+1
      
      
        pphi=phi(i)
        aamu=amu(j)
        aanu=anu(k)

      x=aa*cosh(aamu)* sin(aanu)* cos(pphi)
      y=aa*cosh(aamu)* sin(aanu)* sin(pphi)
      z=aa*sinh(aamu)* cos(aanu)

c is it in G1 ?
c nyes=1 : yes;  nyes=0 : no.
      nyes=1
      call cubengp5(x, y, z, g1, r1, s1, t1, nnx, nny, nnz, 
     1  dr,ds,dt,nyes)
      if (nyes.eq.1) then
c is ig1(nnx,nny,nnz) an interior point?  ig(,,)=1 -> interior point
      if(ig1(nnx,nny,nnz).ne.1) nyes=0
      end if
      if (nyes.eq.1) then
      ig5(i,j,k)=1
      ngp5(i,j,k,1)=nnx
      ngp5(i,j,k,2)=nny
      ngp5(i,j,k,3)=nnz
      drngp5(i,j,k,1)=dr
      drngp5(i,j,k,2)=ds
      drngp5(i,j,k,3)=dt
      go to 2001
      end if  
      if(z.ge.0.) then
c check if it is in G2
c find the nearest grid point on G2
      npole=1
      igrid=2
      call sphngp51(npole, x,y,z, g2, ig2, r2, s2, t2, nnx, nny, nnz,
     1 dr, ds, dt, igrid)
      ig5(i,j,k)=2
      ngp5(i,j,k,1)=nnx
      ngp5(i,j,k,2)=nny
      ngp5(i,j,k,3)=nnz
      drngp5(i,j,k,1)=dr
      drngp5(i,j,k,2)=ds
      drngp5(i,j,k,3)=dt
      if(ig2(nnx,nny,nnz).ne.2) then
      write(*,*) 'G5 points in G2 but not an interior points in G2',
     1 nnx, nny, nnz, ig2(nnx,nny,nnz) 
      write(*,*) 'i,j,k=',i,j,k 
      write(*,*) 'x,y,z=',x,y,z
      write(*,*) 'in G2=', (g2(nnx,nny,nnz,kk),kk=1,3)
      end if
      else
c in G3      
      npole=-1
      igrid=3
c      icount=icount+1
c      if(icount.ge.20) stop
      call sphngp51(npole, x,y,z, g3, ig3, r3, s3, t3, nnx, nny, nnz,
     1 dr, ds, dt, igrid)
      ig5(i,j,k)=3
      ngp5(i,j,k,1)=nnx
      ngp5(i,j,k,2)=nny
      ngp5(i,j,k,3)=nnz
      drngp5(i,j,k,1)=dr
      drngp5(i,j,k,2)=ds
      drngp5(i,j,k,3)=dt
      if(ig3(nnx,nny,nnz).ne.3) then
      write(*,*) 'G5 points in G3 but not an interior points in G3',
     1 nnx, nny, nnz, ig3(nnx,nny,nnz)
      write(*,*) 'x,y,z=',x,y,z
      end if

      end if      
 2001 continue
      end do      
      end do      
      end do      
      
      

c for bndry points on G5b
c
      aamu=amu0
      do i = 2, nphi+1
      do k = 2, nnu+1
          
        pphi=phi(i)
c        aamu=amu(j)
        aanu=anu(k)

      x=aa*cosh(aamu)* sin(aanu)* cos(pphi)
      y=aa*cosh(aamu)* sin(aanu)* sin(pphi)
      z=aa*sinh(aamu)* cos(aanu)

      if(z.ge.0.) then
c check if it is in G2
c find the nearest grid point on G2
      npole=1
      call sphngp31(npole, x,y,z, g2, r2, s2, t2, nnx, nny, nnz,
     1 dr, ds, dt)
      ig5b(i,k)=2
      ngp5b(i,k,1)=nnx
      ngp5b(i,k,2)=nny
      ngp5b(i,k,3)=nnz
      drngp5b(i,k,1)=dr
      drngp5b(i,k,2)=ds
      drngp5b(i,k,3)=dt
      if(ig2(nnx,nny,nnz).ne.2) then
      write(*,*) 'G5b points in G2 but not an interior points in G2'
      end if
      else
c in G3      
      npole=-1
      call sphngp31(npole, x,y,z, g2, r2, s2, t2, nnx, nny, nnz,
     1 dr, ds, dt)
      ig5b(i,k)=3
      ngp5b(i,k,1)=nnx
      ngp5b(i,k,2)=nny
      ngp5b(i,k,3)=nnz
      drngp5b(i,k,1)=dr
      drngp5b(i,k,2)=ds
      drngp5b(i,k,3)=dt
      if(ig3(nnx,nny,nnz).ne.3) then
      write(*,*) 'G5b points in G3 but not an interior points in G3'
      end if
      end if      
     
      end do      
      end do      
      
ccc

c output
      write(71) ig5, ngp5, drngp5, ig5b, ngp5b, drngp5b

c
c for points on G1, G2, G3 to find points on G5ex

c for points on G1
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).ne.0) then
	x=g1(i,j,k,1)
	y=g1(i,j,k,2)
	z=g1(i,j,k,3)
c find the nearest grid point on G5ex
      call spheroid5ex(x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)
      ngp15(i,j,k,1)=nnx
      ngp15(i,j,k,2)=nny
      ngp15(i,j,k,3)=nnz
      drngp15(i,j,k,1)=dr
      drngp15(i,j,k,2)=ds
      drngp15(i,j,k,3)=dt
	end if 
      end do
      end do
      end do
c      stop
c for points on G2
      do i=1,nx2
      do j=1,ny2
c excluding the bndry points at k=nz2+1
      do k=1,nz2
      if(ig2(i,j,k).ne.0) then
	x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c find the nearest grid point on G5
      call spheroid5ex(x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)
           

      ngp25(i,j,k,1)=nnx
      ngp25(i,j,k,2)=nny
      ngp25(i,j,k,3)=nnz
      drngp25(i,j,k,1)=dr
      drngp25(i,j,k,2)=ds
      drngp25(i,j,k,3)=dt
	end if 
      end do
      end do
      end do

c for points on G3
      do i=1,nx3
      do j=1,ny3
c excluding the bndry points at k=nz2+1
      do k=1,nz3
      if(ig3(i,j,k).ne.0) then
	x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c find the nearest grid point on G5
      call spheroid5ex(x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)


      ngp35(i,j,k,1)=nnx
      ngp35(i,j,k,2)=nny
      ngp35(i,j,k,3)=nnz
      drngp35(i,j,k,1)=dr
      drngp35(i,j,k,2)=ds
      drngp35(i,j,k,3)=dt
	end if 
      end do
      end do
      end do
c
c output
      write(72) ngp15, drngp15, ngp25, drngp25, ngp35, drngp35
c
c for points on the surfaces of G2, G3 to find points on G5bex

c for points on G2 surface
      k=nz2-1
      do i=1,nx2
      do j=1,ny2
      if(ig2(i,j,k).ne.0) then
	x=g2(i,j,k,1)
	y=g2(i,j,k,2)
	z=g2(i,j,k,3)
c find the nearest grid point on G5
      call spheroid5ex(x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)

c
c note nny should equal to nmuh+1 but not used 
      
      ngp25b(i,j,1)=nnx
      ngp25b(i,j,2)=nnz
      drngp25b(i,j,1)=dr
      drngp25b(i,j,2)=dt
	end if 
      end do
      end do
      

c for points on G3
      k=nz3-1
      do i=1,nx3
      do j=1,ny3
c excluding the bndry points at k=nz2+1
      
      if(ig3(i,j,k).ne.0) then
	x=g3(i,j,k,1)
	y=g3(i,j,k,2)
	z=g3(i,j,k,3)
c find the nearest grid point on G5
      call spheroid5ex(x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)


      ngp35b(i,j,1)=nnx
      ngp35b(i,j,2)=nnz
      drngp35b(i,j,1)=dr
      drngp35b(i,j,2)=dt
	end if 
      end do
      end do
c
c output
      write(73) ngp25b, drngp25b, ngp35b, drngp35b



      stop
	end	


	subroutine cube(g1, r1, s1, t1)
      parameter (mx1=51,my1=51,mz1=51)
	dimension g1(mx1,my1,mz1,3), r1(mx1), s1(my1), t1(mz1) 

c square Grid 1
	xs=-0.65
	xe=0.65
	ys=-0.65
	ye=0.65
	zs=-0.65
	ze=0.65

	nx1=mx1
	ny1=my1
        nz1=mz1
	nx11=nx1-1
	ny11=ny1-1
	nz11=nz1-1
	do i=1,nx1
	r1(i)=float(i-1)/float(nx11)
	end do
	do j=1,ny1
	s1(j)=float(j-1)/float(ny11)
	end do
	do k=1,nz1
	t1(k)=float(k-1)/float(nz11)
	end do
	do k=1,nz1
	do j=1,ny1
	do i=1,nx1
	g1(i,j,k,1)=xs+r1(i)*(xe-xs)
	g1(i,j,k,2)=ys+s1(j)*(ye-ys)
	g1(i,j,k,3)=zs+t1(k)*(ze-zs)
	end do
	end do
	end do
      return
	end

	subroutine cubengp(x, y,z,g1,r1,s1,t1,nnx,nny,nnz)
      parameter (mx1=51,my1=51,mz1=51)
	dimension g1(mx1,my1,mz1,3), r1(mx1), s1(my1), t1(mz1) 

	xs=-0.65
	xe=0.65
	ys=-0.65
	ye=0.65
	zs=-0.65
	ze=0.65

      nx1=mx1
      ny1=my1
      nz1=mz1
  
      if(x.lt.xs .or. x.gt.xe .or.	
     1 y.lt.ys .or. y.gt.ye .or.
     1 z.lt.zs .or. z.gt.ze ) then
c out of bounds
      write(*,*) 'in cubengp, out of bounds, x,y,z=',x,y,z
      return 
      end if

      rr1=(x-xs)/(xe-xs)
      ss1=(y-ys)/(ye-ys)
      tt1=(z-zs)/(ze-zs)
      nnx1=int(rr1*float(nx1-1)*1.0)+1
      nny1=int(ss1*float(ny1-1)*1.0)+1
	nnz1=int(tt1*float(nz1-1)*1.0)+1

      dd=(rr1-r1(nnx1))**2+(ss1-s1(nny1))**2
     1	+(tt1-t1(nnz1))**2
       nnx=nnx1
       nny=nny1
       nnz=nnz1
	do k1=1,2
	do k2=1,2
	do k3=1,2
      ddd=(rr1-r1(nnx1+k1-1))**2+(ss1-s1(nny1+k2-1))**2
     1	+(tt1-t1(nnz1+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx1+k1-1
       nny=nny1+k2-1
       nnz=nnz1+k3-1
      end if
      end do
	end do
	end do

	return
	end

	subroutine cubengp3(x,y,z,g1,r1,s1,t1,nnx,nny,nnz,dr,ds,dt)
      parameter (mx1=51,my1=51,mz1=51)
	dimension g1(mx1,my1,mz1,3), r1(mx1), s1(my1), t1(mz1) 

	xs=-0.65
	xe=0.65
	ys=-0.65
	ye=0.65
	zs=-0.65
	ze=0.65

      nx1=mx1
      ny1=my1
      nz1=mz1
  
      if(x.lt.xs .or. x.gt.xe .or.	
     1 y.lt.ys .or. y.gt.ye .or.
     1 z.lt.zs .or. z.gt.ze ) then
c out of bounds
      write(*,*) 'in cubengp, out of bounds, x,y,z=',x,y,z
      return 
      end if

      rr1=(x-xs)/(xe-xs)
      ss1=(y-ys)/(ye-ys)
      tt1=(z-zs)/(ze-zs)
      nnx1=int(rr1*float(nx1-1)*1.0)+1
      nny1=int(ss1*float(ny1-1)*1.0)+1
	nnz1=int(tt1*float(nz1-1)*1.0)+1

      dd=(rr1-r1(nnx1))**2+(ss1-s1(nny1))**2
     1	+(tt1-t1(nnz1))**2
       nnx=nnx1
       nny=nny1
       nnz=nnz1
	do k1=1,2
	do k2=1,2
	do k3=1,2
      ddd=(rr1-r1(nnx1+k1-1))**2+(ss1-s1(nny1+k2-1))**2
     1	+(tt1-t1(nnz1+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx1+k1-1
       nny=nny1+k2-1
       nnz=nnz1+k3-1
      end if
      end do
	end do
	end do

       dr=(rr1-r1(nnx))*(xe-xs)
       ds=(ss1-s1(nny))*(ye-ys)
       dt=(tt1-t1(nnz))*(ze-zs)

	return
	end



      subroutine sphere(npole, g, r2, s2, t2)
      parameter (mx=71,my=71,mz=61+1)
      dimension g(mx,my,mz,3),  r2(mx), s2(my), t2(mz) 
	common/parm/factor
c orthographic transformation
  
	sa=1.8
	sb=1.8
	nx=mx
	ny=my
      nz=mz
	nx1=nx-1
	ny1=ny-1
      nz1=nz-2
	rs=-1.
	re=1.
      as=0.4
      ae=1.
c      factor=299./300.

	pi=3.14159265359
      do k=1,nz
      t2(k)=float(k-1)/float(nz1)
      end do
	do i=1,nx
	r2(i)=float(i-1)/float(nx1)
	end do
	do j=1,ny
	s2(j)=float(j-1)/float(ny1)
	end do

      do k=1,nz
      a=as+t2(k)*(ae-as)
      b=a*factor
	do i=1,nx
	rr1=(rs+r2(i)*(re-rs))*sa
	do j=1,ny
	rr2=(rs+s2(j)*(re-rs))*sb
      r=sqrt(rr1**2+rr2**2)
	sig=r
	sig2=sig**2
      cosphi=(1.-sig2)/(1.+sig2)*float(npole)
	sinphi=2.*sig/(1.+sig2)
	if(sig.ne.0.) then
	costh=rr1/sig
      sinth=rr2/sig*float(npole)
	else
c the following 2 statements do not matter.
	costh=0.
	sinth=0.
	end if
	g(i,j,k,1)=a*sinphi*costh
	g(i,j,k,2)=a*sinphi*sinth
      g(i,j,k,3)=b*cosphi
	end do
	end do
	end do
	return
	end

      subroutine sphngp(npole,x,y,z,g,r2,s2,t2,nnx,nny,nnz)
      parameter (mx=71,my=71,mz=61+1)
      dimension g(mx,my,mz,3),  r2(mx), s2(my), t2(mz) 
      common/parm/factor
	nx=mx
	ny=my
      nz=mz
      nx1=nx-1
      ny1=ny-1
      nz1=nz-2
c	a=1.
c	b=1.
	sa=1.8
	sb=1.8
	rs=-1.
	re=1.
      as=0.4
      ae=1.
C factor=b/a
c      factor=299./300.
c
c convert it to G3
	rr=sqrt(x**2+y**2)
c tan phi= rr factor/z
c
      phi=atan2(rr*factor, z)     
	sinphi=sin(phi)
	cosphi=cos(phi)
      a=sqrt(rr*rr+(z/factor)**2)
      b=a*factor
	if(rr.ne.0.) then
	costh=x/rr
	sinth=y/rr
	else
	costh=0.
	sinth=0.
	end if
	fact=1.+float(npole)*cosphi
	rr1=sinphi*costh/fact
	rr2=float(npole)*sinphi*sinth/fact
	rr1=(rr1/sa-rs)/(re-rs)
	rr2=(rr2/sb-rs)/(re-rs)
	rr3=(a-as)/(ae-as)

c find the nearest grid point on G3
      nnx2=int(rr1*float(nx1)*1.0)+1
      nny2=int(rr2*float(ny1)*1.0)+1
      nnz2=int(rr3*float(nz1)*1.0)+1
	if(nnx2.le.1 .or. nnx2.ge.nx-1 .or.
     1nny2.le.1 .or. nny2.ge.ny-1) then

	nnx=nnx2
	nny=nny2
	nnz=nnz2
      write(*,*) 'stop in sphngp'
      stop
      else
      if(nnz2.ge.nz) nnz2=nz
      if(nnz2.le.1) nnz2=1
	dd=(rr1-r2(nnx2))**2+(rr2-s2(nny2))**2+(rr3-t2(nnz2))**2
      nnx=nnx2
      nny=nny2
      nnz=nnz2
c
      k3up=2
      if(nnz2.ge.nz) k3up=1      
	do k1=1,2
	do k2=1,2
	do k3=1,k3up
	ddd=(rr1-r2(nnx2+k1-1))**2+(rr2-s2(nny2+k2-1))**2
     1   +(rr3-t2(nnz2+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx2+k1-1
       nny=nny2+k2-1
       nnz=nnz2+k3-1
      end if
      end do
	end do
	end do
	end if

      return
	end

      subroutine sphngp31(npole,x,y,z,g,r2,s2,t2,nnx,nny,nnz,
     1 dr, ds, dt)
      parameter (mx=71,my=71,mz=61+1)
      dimension g(mx,my,mz,3),  r2(mx), s2(my), t2(mz) 
      common/parm/factor
	nx=mx
	ny=my
      nz=mz
      nx1=nx-1
      ny1=ny-1
      nz1=nz-2
c	a=1.
c	b=1.
	sa=1.8
	sb=1.8
	rs=-1.
	re=1.
      as=0.4
      ae=1.
C factor=b/a
c      factor=299./300.
c
c convert it to G3
	rr=sqrt(x**2+y**2)
c tan phi= rr factor/z
c
      phi=atan2(rr*factor, z)     
	sinphi=sin(phi)
	cosphi=cos(phi)
      a=sqrt(rr*rr+(z/factor)**2)
      b=a*factor
	if(rr.ne.0.) then
	costh=x/rr
	sinth=y/rr
	else
	costh=0.
	sinth=0.
	end if
	fact=1.+float(npole)*cosphi
	rr1=sinphi*costh/fact
	rr2=float(npole)*sinphi*sinth/fact
	rr1=(rr1/sa-rs)/(re-rs)
	rr2=(rr2/sb-rs)/(re-rs)
	rr3=(a-as)/(ae-as)

c find the nearest grid point on G3
      nnx2=int(rr1*float(nx1)*1.0)+1
      nny2=int(rr2*float(ny1)*1.0)+1
      nnz2=int(rr3*float(nz1)*1.0)+1
	if(nnx2.le.1 .or. nnx2.ge.nx-1 .or.
     1nny2.le.1 .or. nny2.ge.ny-1) then

	nnx=nnx2
	nny=nny2
	nnz=nnz2
      write(*,*) 'stop in sphngp'
      stop
      else
      if(nnz2.ge.nz) nnz2=nz
      if(nnz2.le.1) nnz2=1
	dd=(rr1-r2(nnx2))**2+(rr2-s2(nny2))**2+(rr3-t2(nnz2))**2
      nnx=nnx2
      nny=nny2
      nnz=nnz2
c
      k3up=2
      if(nnz2.ge.nz) k3up=1      
	do k1=1,2
	do k2=1,2
	do k3=1,k3up
	ddd=(rr1-r2(nnx2+k1-1))**2+(rr2-s2(nny2+k2-1))**2
     1   +(rr3-t2(nnz2+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx2+k1-1
       nny=nny2+k2-1
       nnz=nnz2+k3-1
      end if
      end do
	end do
	end do
	end if
       dr=(rr1-r2(nnx))*(re-rs)
       ds=(rr2-s2(nny))*(re-rs)
       dt=(rr3-t2(nnz))*(ae-as)
      return
	end
      subroutine sphngp32(npole,x,y,z,g,r2,s2,t2,nnx,nny,nnz,
     1 dr, ds, dt)
      parameter (mx=71,my=71,mz=61+1)
      dimension g(mx,my,mz,3),  r2(mx), s2(my), t2(mz) 
      common/parm/factor
	nx=mx
	ny=my
      nz=mz
      nx1=nx-1
      ny1=ny-1
      nz1=nz-2
c	a=1.
c	b=1.
	sa=1.8
	sb=1.8
	rs=-1.
	re=1.
      as=0.4
      ae=1.
C factor=b/a
c      factor=299./300.
c
c convert it to G3
	rr=sqrt(x**2+y**2)
c tan phi= rr factor/z
c
      
      phi=atan2(rr*factor, z)     
	sinphi=sin(phi)
	cosphi=cos(phi)
      a=sqrt(rr*rr+(z/factor)**2)
      b=a*factor
	if(rr.ne.0.) then
	costh=x/rr
	sinth=y/rr
	else
	costh=0.
	sinth=0.
	end if
	fact=1.+float(npole)*cosphi
	rr1=sinphi*costh/fact
	rr2=float(npole)*sinphi*sinth/fact
	rr1=(rr1/sa-rs)/(re-rs)
	rr2=(rr2/sb-rs)/(re-rs)
	rr3=(a-as)/(ae-as)

c find the nearest grid point on G3
      nnx2=int(rr1*float(nx1)*1.0)+1
      nny2=int(rr2*float(ny1)*1.0)+1
      nnz2=int(rr3*float(nz1)*1.0)+1
	if(nnx2.le.1 .or. nnx2.ge.nx-1 .or.
     1nny2.le.1 .or. nny2.ge.ny-1) then

	nnx=nnx2
	nny=nny2
	nnz=nnz2
      write(*,*) 'stop in sphngp'
      stop
      else
      if(nnz2.ge.nz) nnz2=nz
      if(nnz2.le.1) nnz2=1
	dd=(rr1-r2(nnx2))**2+(rr2-s2(nny2))**2+(rr3-t2(nnz2))**2
      nnx=nnx2
      nny=nny2
      nnz=nnz2
c
      k3up=2
      if(nnz2.ge.nz) k3up=1      
	do k1=1,2
	do k2=1,2
	do k3=1,k3up
	ddd=(rr1-r2(nnx2+k1-1))**2+(rr2-s2(nny2+k2-1))**2
     1   +(rr3-t2(nnz2+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx2+k1-1
       nny=nny2+k2-1
       nnz=nnz2+k3-1
      end if
      end do
	end do
	end do
	end if
       dr=(rr1-r2(nnx))*(re-rs)
       ds=(rr2-s2(nny))*(re-rs)
c       dt=(rr3-t2(nnz))*(ae-as)
c dt=0 because on the same shell
       dt=0      
      return
	end


      subroutine transform(npole, trans, g2, r2, s2, t2)
      parameter (mx=71,my=71,mz=61+1)
      dimension  trans(mx,my,mz,3,3), r2(mx), s2(my), t2(mz) 
      dimension  g2(mx,my,mz,3) 
      common/parm/factor 

c orthographic transformation
  
	sa=1.8
	sb=1.8
	nx=mx
	ny=my
      nz=mz
	nx1=nx-1
	ny1=ny-1
      nz1=nz-2
	rs=-1.
	re=1.
      as=0.4
      ae=1.
c      factor=299./300.

	pi=3.14159265359
      do k=1,nz
      t2(k)=float(k-1)/float(nz1)
      end do
	do i=1,nx
	r2(i)=float(i-1)/float(nx1)
	end do
	do j=1,ny
	s2(j)=float(j-1)/float(ny1)
	end do

      do k=1,nz
      a=as+t2(k)*(ae-as)
      b=a*factor
	do i=1,nx
	rr1=(rs+r2(i)*(re-rs))*sa
	do j=1,ny
	rr2=(rs+s2(j)*(re-rs))*sb

c change notation r=rr1, s=rr2, t=a
      r=rr1
      s=rr2
      t=a
c r:[rs, re]; s:[rs, re]; t:[as, ae]
c trans(1,1)= dr/dx= -(-1 + r^2 - s^2)/(2*t*sa)
c trans(1,2)= dr/dy= -((r*s)/(npole*t*sa)
c trans(1,3)= dr/dz= -(r/(fac*npole*t*sa))
c trans(2,1)= ds/dx= -((r*s)/(t*sb))
c trans(2,2)= ds/dy= (1 + r^2 - s^2)/(2*npole*t*sb)
c trans(2,3)= ds/dz= -(s/(fac*npole*t*sb))
c trans(3,1)= dt/dx= (2*r)/(1 + r^2 + s^2)
c trans(3,2)= dt/dy= (2*s)/(npole*(1 + r^2 + s^2))
c trans(3,3)= dt/dz= -((-1 + r^2 + s^2)/(fac*npole*(1 + r^2 + s^2)))

      fpole=npole
      sig= 1. + r*r + s*s      
      trans(i,j,k,1,1)=  -(-1. + r*r - s*s)/(2.*t*sa)
      trans(i,j,k,1,2)=  -(r*s)/(fpole*t*sa)
      trans(i,j,k,1,3)=  -r/(factor*fpole*t*sa)
      trans(i,j,k,2,1)=  -(r*s)/(t*sb)
      trans(i,j,k,2,2)=  (1. + r*r - s*s)/(2.*fpole*t*sb)
      trans(i,j,k,2,3)=  -s/(factor*fpole*t*sb)
      trans(i,j,k,3,1)=  (2.*r)/sig
      trans(i,j,k,3,2)=  (2.*s)/(fpole*sig)
      trans(i,j,k,3,3)=  -(-1. + r*r + s*s)/(factor*fpole*sig)


	end do
	end do
	end do
c
c check
c a. t2:[as ae]; r2:[rs re]; s2:[rs re]
c      do k=1,nz
c      a=as+t2(k)*(ae-as)
c      b=a*factor
c      t2(k)=a
c	end do
c	do i=1,nx
c	rr1=(rs+r2(i)*(re-rs))
c      r2(i)=rr1
c	end do
c	do j=1,ny
c	rr2=(rs+s2(j)*(re-rs))
c     s2(j)=rr2
c      end do


	return
	end



	
      subroutine spheroid5ex (x,y,z,nnx,nny,nnz,dr,ds,dt, nphi,nmuh,nnu
     1 ,phi,amu,anu, aa, dphi, dmu, dnu)
      dimension amu(nmuh+2), anu(nnu+2), phi(nphi+2)
c convert it to G5ex
      pi=4.*atan(1. )
      aa2=aa*aa
      a2=aa2
      r2=x**2+y**2
	rr=sqrt(r2)
c tan phi= rr factor/z
c
c (a) find mu such that (x^2+y^2)/cosh(mu)^2+ z^2/sinh(mu)^2 close to a^2
c (b) given mu find nu from z=a sinh(mu) cos nu. z>0 for 0<nu<pi/2; z<0 for pi/2<nu<pi.
c (c) from ratio of x/y find phi if nu ne 0 or pi

c (a) 
c the following trick works
c      if(z.eq.0.) z=z+0.0000001 
c
      z2=z**2
   
      if(z.eq.0.) then
c      write(*,*) 'z=0'
      
      if(r2.gt.a2) then
      aanu=pi/2.
      coshmu=sqrt(r2/a2)
      aamu=acosh(coshmu)
      end if
      if(r2.lt.a2) then
      aamu=0.
      sinnu=sqrt(r2/a2)
      aanu=asin(sinnu)
c here anu is in [0 pi/2] degenerate      
      end if
      if(r2.eq.a2) then
      aamu=0.
      aanu=pi/2.
      end if


      
c      write(*,*) 'z=0, r2 a2, amu, anu='
c      write(*,*) r2, a2, amu, anu
      
      else
c q=coshhmu^2
      qq=(a2+z2+r2)**2-4.*a2*r2
      if(qq .le. 0.) then
      write(*,*) 'qq < 0 abort.'
      stop
      else
      qq=sqrt(qq)
      q1=0.5/a2*(a2+z2+r2+qq)
      q2=0.5/a2*(a2+z2+r2-qq)
c      write(*,*) 'q1,q2=', q1, q2
c      if(q1.ge.qmax) then
c      write(*,*) 'q1 too large'
c      end if
       coshmu=sqrt(q1)
       aamu=acosh(coshmu)
       sinnu=sqrt(q2)
       aanu=asin(sinnu)
       if(z.lt.0.) aanu=pi-aanu
      end if
      
      end if

      
c find nu; q=sinnu^2
c      q= r2/a2/coshmu**2
c      anu=asin(sqrt(q))
      
c find phi
              
      pphi=atan2(y, x)
c      write(*,*) 'phi2=', phi     

      if(pphi.lt.0.) pphi=2.*pi+pphi

c      write(*,*) aamu, aanu, pphi
c      stop
c??      
c given aamu, aanu and pphi find the nearest grid
c  grid in mu = (j-1.5)*dmu:  j=1 -> -1/2 dmu
c                               j=2     1/2 dmu
c                               j=nmu/2+1 -> amu0-1/2 dmu
c                               j=nmu/2+2 -> amu0+1/2 dmu
c      jj=int(1.5+amu/dmu)

      jj=0
      do j=2, nmuh+2
      if(aamu.ge.(amu(j)-0.5*dmu) .and. aamu.lt.(amu(j)+0.5*dmu)) jj=j
      if(aamu.eq.amu0+0.5*dmu) jj=nmuh+1
      end do
      if(jj.eq.0) then
      write(*,*) 'something is wrong, out of bound in mu'
      write(*,*) aamu, dmu
      stop
      end if
      if(jj.eq.nmuh+2) jj=nmuh+1
c nu= (k-1.5) dnu  k=1,...,nnu+2
c      kk=int(1.5+anu/dnu)
      anu0=pi
      kk=0
      do k=2, nnu+1
      if(aanu.ge.(anu(k)-0.5*dnu) .and. aanu.lt.(anu(k)+0.5*dnu)) kk=k
      if(aanu.eq.anu0) kk=nnu+1
      end do
      if(kk.eq.0) then
      write(*,*) 'something is wrong, out of bound in nu'
      write(*,*) 'nu, dnu=', aanu, dnu
      write(*,*) 'x,y,z=',x,y,z
      stop
      end if


c phi=(i-2)*dphi; i=1,...,nphi+2; dphi=2*pi/nphi
c     i=1, phi=-dphi; i=2, phi=0; i=nphi+1 phi=2pi-dphi; i=nphi+2 phi=2pi            
c      ii=2+int(phi/dphi)
      pphi0=2.*pi+dphi
      ii=0
      do i=2, nphi+2
      if(pphi.ge.(phi(i)-0.5*dphi) .and. pphi.lt.(phi(i)+0.5*dphi)) ii=i
c      if(pphi.eq.pphi0) ii=nphi+1
      end do
      
      if(ii.eq.nphi+2) then
      ii=2
      pphi=2.*pi-pphi
      end if
            
             
      if(ii.eq.0) then
      write(*,*) 'something is wrong, out of bound in phi'
      write(*,*) 'phi=', pphi, 'pphi0, dphi', pphi0, dphi
      stop
c      write(*,*) (phi(i),i=1,nphi+2)
      end if

      nnx=ii
      nny=jj
      nnz=kk


c find the nearest grid point on G3
       dr=pphi-phi(ii)
       ds=aamu-amu(jj)
       dt=aanu-anu(kk)
      return
	end
	
	



	subroutine cubengp5(x,y,z,g1,r1,s1,t1,nnx,nny,nnz,dr,ds,dt,nyes)
      parameter (mx1=51,my1=51,mz1=51)
	dimension g1(mx1,my1,mz1,3), r1(mx1), s1(my1), t1(mz1) 

	xs=-0.65
	xe=0.65
	ys=-0.65
	ye=0.65
	zs=-0.65
	ze=0.65

      nx1=mx1
      ny1=my1
      nz1=mz1
  
      if(x.lt.xs .or. x.gt.xe .or.	
     1 y.lt.ys .or. y.gt.ye .or.
     1 z.lt.zs .or. z.gt.ze ) then
c out of bounds
      nyes=0
c      write(*,*) 'in cubengp, out of bounds, x,y,z=',x,y,z
      return 
      end if

      rr1=(x-xs)/(xe-xs)
      ss1=(y-ys)/(ye-ys)
      tt1=(z-zs)/(ze-zs)
      nnx1=int(rr1*float(nx1-1)*1.0)+1
      nny1=int(ss1*float(ny1-1)*1.0)+1
	nnz1=int(tt1*float(nz1-1)*1.0)+1

      dd=(rr1-r1(nnx1))**2+(ss1-s1(nny1))**2
     1	+(tt1-t1(nnz1))**2
       nnx=nnx1
       nny=nny1
       nnz=nnz1
	do k1=1,2
	do k2=1,2
	do k3=1,2
      ddd=(rr1-r1(nnx1+k1-1))**2+(ss1-s1(nny1+k2-1))**2
     1	+(tt1-t1(nnz1+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx1+k1-1
       nny=nny1+k2-1
       nnz=nnz1+k3-1
      end if
      end do
	end do
	end do

       dr=(rr1-r1(nnx))*(xe-xs)
       ds=(ss1-s1(nny))*(ye-ys)
       dt=(tt1-t1(nnz))*(ze-zs)

	return
	end


      subroutine sphngp51(npole,x,y,z,g,ig,r2,s2,t2,nnx,nny,nnz,
     1 dr, ds, dt, igrid)
c
c find nearest grid point on Grid 2 (3) for the points on Grid 1     
      parameter (mx=71,my=71,mz=61+1)
      dimension g(mx,my,mz,3), ig(mx,my,mz), r2(mx), s2(my), t2(mz) 
      common/parm/factor
	nx=mx
	ny=my
      nz=mz
      nx1=nx-1
      ny1=ny-1
      nz1=nz-2
c	a=1.
c	b=1.
	sa=1.8
	sb=sa
	rs=-1.
	re=1.
      as=0.4
      ae=1.
C factor=b/a
c      factor=299./300.
c
c convert it to G3
	rr=sqrt(x**2+y**2)
c tan phi= rr factor/z
c
      phi=atan2(rr*factor, z)     
	sinphi=sin(phi)
	cosphi=cos(phi)
      a=sqrt(rr*rr+(z/factor)**2)
      b=a*factor
	if(rr.ne.0.) then
	costh=x/rr
	sinth=y/rr
	else
	costh=0.
	sinth=0.
	end if
	fact=1.+float(npole)*cosphi
	rr1=sinphi*costh/fact
	rr2=float(npole)*sinphi*sinth/fact
	rr1=(rr1/sa-rs)/(re-rs)
	rr2=(rr2/sb-rs)/(re-rs)
	rr3=(a-as)/(ae-as)

c find the nearest grid point on G3
      nnx2=int(rr1*float(nx1)*1.0)+1
      nny2=int(rr2*float(ny1)*1.0)+1
      nnz2=int(rr3*float(nz1)*1.0)+1
	if(nnx2.le.1 .or. nnx2.ge.nx-1 .or.
     1nny2.le.1 .or. nny2.ge.ny-1) then

	nnx=nnx2
	nny=nny2
	nnz=nnz2
      write(*,*) 'stop in sphngp'
      stop
      else
      if(nnz2.gt.nz-2) then
       write(*,*) 'sphngp51', nnz2, nz-1, nz
       nnz2=nz-2
      end if 
      if(nnz2.le.1) nnz2=1
	dd=(rr1-r2(nnx2))**2+(rr2-s2(nny2))**2+(rr3-t2(nnz2))**2
      nnx=nnx2
      nny=nny2
      nnz=nnz2
c
      k3up=2
      if(nnz2.eq.nz-2) k3up=1      
c      if(nnz2.eq.nz-3) k3up=2      
c      if(nnz2.eq.nz-4) k3up=3      
	do k1=1,2
	do k2=1,2
	do k3=1,k3up
	
	ddd=(rr1-r2(nnx2+k1-1))**2+(rr2-s2(nny2+k2-1))**2
     1   +(rr3-t2(nnz2+k3-1))**2
      if(ddd.lt.dd) then
       dd=ddd
       nnx=nnx2+k1-1
       nny=nny2+k2-1
       nnz=nnz2+k3-1
      end if
      end do
	end do
	end do
	end if
c check if ig(nnx,nny,nnz).eq.igrid
      if (ig(nnx,nny,nnz).ne.igrid) then
      write(*,*) 'ig ne igrid'
c      k3up=3
c      if(nnz.eq.nz-2) k3up=1      
c      if(nnz.eq.nz-3) k3up=2      
c	do k1=1,3
c	do k2=1,3
c	do k3=1,k3up
c       nnx2=nnx+k1-2
c       nny2=nny+k2-2
c       nnz2=nnz+k3-2
c      if(ig(nnx2,nny2,nnz2).eq.igrid) go to 300
c      end do
c	end do
c	end do
c	write(*,*) 'no luck'
c  300 continue	
c      nnx=nnx2
c      nny=nny2
c      nnz=nnz2
      end if
      
      

       dr=(rr1-r2(nnx))*(re-rs)
       ds=(rr2-s2(nny))*(re-rs)
       dt=(rr3-t2(nnz))*(ae-as)
      return
	end

      subroutine gquad(l,root,w)
c
c     Finds the l roots (in theta) and associated gaussian weights 
c     for the legendre polynomial of degree l > 1.
c     Note: usually Gauss integration is done in mu=cos(theta)
c     with -1.le.mu.le.1. In converting from mu to theta
c     an extra factor of sin(theta) is required and this is
c     incorporated by the subroutine below.
c     Subroutine required: pbar 
c
      implicit none
	integer l,l1,l2,l3,l22,k,i
      real root(l),w(l),pi,del,co,p,p1,p2,s,t1,t2,theta
c
      pi=4.*atan(1. )
      del=pi/float(4*l)
      l1=l+1
      co=float(2*l+3)/float(l1**2)
      p2=1.
      t2=-del
      l2=l/2
      k=1
c
      do 10 i=1,l2
   20 t1=t2
      t2=t1+del
      theta=t2
      call pbar(theta,l,0,p)
      p1=p2
      p2=p 
      if ((k*p2).gt.0. ) go to 20
      k=-k
   40 s=(t2-t1)/(p2-p1)
      t1=t2
      t2=t2-s*p2
      theta=t2
      call pbar(theta,l,0,p)
      p1=p2
      p2=p 
      if (abs(p).le.1.e-10) go to 30
      if (p2.eq.p1) then
         write(6,*) 'sub gquad: zero = ',p,' at i = ',i
         go to 30
      endif
      go to 40
   30 root(i)=theta
      call pbar(theta,l1,0,p)
      w(i)=co*(sin(theta)/p)**2
   10 continue
c
      l22=2*l2
      if (l22.eq.l) go to 70
      l2=l2+1
      theta=pi/2.
      root(l2)=theta
      call pbar(theta,l1,0,p)
      w(l2)=co/p**2
   70 continue
c
      l3=l2+1
      do 50 i=l3,l
      root(i)=pi-root(l-i+1)
      w(i)=w(l-i+1)
   50 continue
c
      return
      end
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      subroutine pbar(the,l,m,p)
c
c     pbar calculates the value of the normalized associated
c     legendre function of the first kind, of degree l,
c     of order m, for the real argument cos(the),
c     0 .le. m .le. l
c
      implicit none
      integer l,m,m1,i,j
      real the,p,p1,p2,s,c
c
      s=sin(the)
      c=cos(the)
      p=1. /sqrt(2. )
      if (m.eq.0) go to 22
      do 20 i=1,m
        p=sqrt(float(2*i+1)/float(2*i))*s*p
   20 continue
   22 continue
      if (l.eq.m) return
      p1=1. 
      m1=m+1
      do 30 j=m1,l
        p2=p1
        p1=p
        p=2.0*sqrt((float(j**2)-0.25 )/float(j**2-m**2))*c*p1
     $   -sqrt(float((2*j+1)*(j-m-1)*(j+m-1))/
     $          float((2*j-3)*(j-m)*(j+m)))*p2
   30 continue
c
      return
      end
c
      subroutine sphngp4(x,y,z,g,thgauss,nnx2,nny2,
     1 dr, ds)
c find on G4 the corresponding ngp of a grid point (x,y,z) on G2 or G3
c mx=#of gaussian points; my=#of modes in phi
	parameter(mth=3*22,mphi=3*3*2**4)
      parameter (mx4=mth,my4=mphi)
      dimension g(mx4,my4,2),  thgauss(mx4)
      common/parm/factor 
	nx=mx4
	ny=my4
      pi=3.14159265359
	rr=sqrt(x**2+y**2)
c tan phi= rr factor/z
c th from 0 to pi; phi from 0 to 2p
      th=atan2(rr*factor, z)     
      phi=atan2(y,x)
	if(phi.lt.0.) phi=phi+2.*pi
      if(phi.ge.2.*pi) phi=0.
c
c phi=m*dphi
      dphi=2.*pi/float(ny)
c nny2
      nny2=int(phi/dphi)+1
      ds=phi-dphi*float(nny2-1)
c
c find nnx2
      if(th.lt.thgauss(1)) then
      nnx2=0
      dr=th
      end if

      if(th.ge.thgauss(nx)) then
      nnx2=nx
      dr=th-thgauss(nx)
      end if

      do i=1,nx-1
       if(th.ge.thgauss(i).and.th.lt.thgauss(i+1)) then
        nnx2=i
        dr=th-thgauss(i)
       end if
      end do
            
c      thn=thgauss(nnx2)
c	phin=dphi*float(nny2-1)
c	costh=cos(thn)
c	sinth=sin(thn)
c     cosphi=cos(phin)
c	sinphi=sin(phin)
c	xn=sinth*cosphi
c	yn=sinth*sinphi
c	zn=costh
c	write(*,*) x,y,z,xn,yn,zn,th,phi,thn,phin
c
c     
      return
	end