program incomp
      
c This program reads output data from fluid calculations and
c adds small B fields



        parameter( mtm=4, mn=4, mngp=20000, mx1=51, my1=51, mz1=51 )
        parameter( mx2=71, my2=71, mz2=61+1 )
        parameter( mx3=71, my3=71, mz3=61+1 )

c for grid 5 and 5b and 5 ex
      integer nphi, nmu, nnu, nmuh
      parameter (nphi=64*2, nmu=nphi, nnu=nphi/2, nmuh=nmu/2)
      real amu(nmu), aaanu(nnu), aphi(nphi), amuh(nmuh)
      integer ig5(nphi+2, nmuh+2, nnu+2), ngp5(nphi+2, nmuh+2, nnu+2, 3)
      real drngp5(nphi+2, nmuh+2, nnu+2, 3)
      integer ig5b(nphi+2,nnu+2), ngp5b(nphi+2, nnu+2,3)
      real drngp5b(nphi+2,nnu+2,3), g5b(nphi+2,nnu+2)
      real uc5(nphi+2, nmuh+2, nnu+2, 3), uc5b(nphi+2, nnu+2, 3)
      real uc5old(nphi+2, nmuh+2, nnu+2, 3)
      real cc, aa, aa2, amu0, philength, amulength, anulength
      real dphi, dmu, dnu, pi, pphi, aamu, aanu
c for double
      real uc5p(nphi, nmuh, nnu),uc5pdouble(nphi, nmu, nnu)
      real uc5pex(nphi+2, nmuh+2, nnu+2)
      real bndrynu(nphi,nmu,2), bndrymu(nphi,nnu)
      real bndrynumu(nphi,2), errmax, diff
      real uc5bndry(nphi,nnu,2)
      
c for 5 ex and Grid 1, 2, 3
c      real uc5ex(nphi+2, nmuh+2, nnu+2)
      integer ngp15(mx1,my1,mz1,3)
      real drngp15(mx1,my1,mz1,3)      
      integer ngp25(mx2,my2,mz2,3)
      real drngp25(mx2,my2,mz2,3)      
      integer ngp35(mx3,my3,mz3,3)
      real drngp35(mx3,my3,mz3,3)      
      integer ngp25b(mx2,my2,2)
      real drngp25b(mx2,my2,2)      
      integer ngp35b(mx3,my3,2)
      real drngp35b(mx3,my3,2)      

c
c
c for B bndry
	integer mth, mphi,mwsave,mx4,my4,mz4,imax,mmax,jmax,lmax,lmmax,lmmax1
	parameter(mth=3*22,mphi=3*3*2**4,mwsave=2*mphi+15,
     1 mx4=mth,my4=mphi)

	parameter(imax=mphi, mmax=imax/3*2, jmax=mth,lmax=mmax,
     1 lmmax=lmax+mmax,lmmax1=lmmax+1)
      real ub4(mth,mphi,2,6),rhsb4(mth,mphi,1,1,mtm)
     1 , g4(mth,mphi,2),  dr42(mngp,3), dr42b(mngp,3)
     1 , dr43(mngp,3), dr43b(mngp,3)
      real uc4(mth,mphi,3)
      integer ig4(mth,mphi), ngp42(mngp,3) ,ngp43(mngp,3)
     1 ,ngp42b(mngp,3) ,ngp43b(mngp,3)
      real wsave(mwsave)
 	real gausspt(jmax), snj(jmax)
c     1    ,rrr(0:mmax,jmax,jmax) ,hhh(0:mmax,jmax,jmax)
c     2      ,ttt(0:mmax,jmax,jmax),fff(0:mmax,jmax,jmax)
c      common/bbndryv/gausspt,snj, rrr,hhh,ttt,fff
      common/bbndryv/gausspt, snj
c      ,dg0,
c     1  xxx,ttt,fff,Dxi,Dth,Dfi,dBth,dBfi
c       logical lrstrt,ldump


        dimension uc1(mx1,my1,mz1,mn), rhs1(mx1,my1,mz1,mn-1,mtm)
     1 , ig1(mx1,my1,mz1), g1(mx1,my1,mz1,3), rhsp1(mx1,my1,mz1)
     1 , ngp12(mngp,3), dr12(mngp,3)
     1 , ngp13(mngp,3), dr13(mngp,3)

        dimension uc2(mx2,my2,mz2,mn), rhs2(mx2,my2,mz2,mn-1,mtm)
     1 , ig2(mx2,my2,mz2), g2(mx2,my2,mz2,3), trans2(mx2,my2,mz2,3,3)
     1 , ngp21(mngp,3), dr21(mngp,3), rhsp2(mx2,my2,mz2)
     1 , ngp23(mngp,3), dr23(mngp,3)
     1 , uc2bndry(mx2,my2,3,9)

        dimension uc3(mx3,my3,mz3,mn), rhs3(mx3,my3,mz3,mn-1,mtm)
     1 , ig3(mx3,my3,mz3), g3(mx3,my3,mz3,3), trans3(mx3,my3,mz3,3,3)
     1 , ngp31(mngp,3), dr31(mngp,3), rhsp3(mx3,my3,mz3)
     1 , ngp32(mngp,3), dr32(mngp,3)
     1 , uc3bndry(mx3,my3,3,9)
     
        real  ub1(mx1,my1,mz1,mn), rhsb1(mx1,my1,mz1,mn-1,mtm)

        real  dr24(mngp,3),  uv2bndry(mx2,my2,3,3)
     1 , udv2bndry(mx2,my2,3)

        integer   ngp24(mngp,3)
        real  ub2(mx2,my2,mz2,mn), rhsb2(mx2,my2,mz2,mn-1,mtm)

        real dr34(mngp,3), uv3bndry(mx3,my3,3,3)
     1 , udv3bndry(mx3,my3,3)
        integer   ngp34(mngp,3)
        real  ub3(mx3,my3,mz3,mn), rhsb3(mx3,my3,mz3,mn-1,mtm)
        real ub5(nphi+2, nnu+2)

      character*12 real_clock(2)


      common /var1/ omega, omegap, omegapx,omegapz,ca,cb,anu,eta,factor 




c for grid 5
      cc=0.8
      amu0=0.5*(alog(1.+cc)-alog(1.-cc))
      aa2=1.-cc*cc
      aa=sqrt(aa2)
        pi=4.*atan(1.)
        philength=2.*pi
        amulength=2.*amu0
        anulength=pi

        dphi = philength / float(nphi)
        
        dmu = amulength / float(nmu)       

        dnu = anulength / float(nnu)       
c
        do i = 1, nphi
          aphi(i) = (i-1) * dphi
        enddo

        do j = 1, nmu
          amu(j) = -amu0 + (float(j)-0.5) * dmu
        enddo

        do j = 1, nmuh
          amuh(j) =  (float(j)-0.5) * dmu
        enddo


        do k = 1, nnu
          aaanu(k) = (float(k)-0.5) * dnu
        enddo
        
        ifirst=0
c



c parameters
      factor=cc
      ca=1.
      cb=ca*factor

	omega=1.
      omegap=0.25

      anu=0.0019*0.6
      omegapx=omegap

      dr1=2.*0.65/(mx1-1)
      ds1=dr1
      dt1=dr1

      dr2=2./float(mx2-1)
      ds2=2./float(my2-1)
      dt2=(1.-0.4)/float(mz2-2)

      dr3=2./float(mx3-1)
      ds3=2./float(my3-1)
      dt3=(1.-0.4)/float(mz3-2)
      
      ds4=2.*3.14159265359/(float(my4))
      dt4=dt3


c lrstrt=0, initialization; =1 restart
c ldump=0, no save; =1 save
c tf= final time; mt: order in time integration
      lrstrt=1
      ldump=0
      ntot=40002
      tf=600.*2.
      mt=2
      cfl=0.6

c**************  begin reading in parameters  ***********************

      nx1=mx1
      ny1=my1
      nz1=mz1
      nxs1=1
      nxe1=nx1
      nys1=1
      nye1=ny1
      nzs1=1
      nze1=nz1

      nx2=mx2
      ny2=my2
      nz2=mz2
      nxs2=1
      nxe2=nx2
      nys2=1
      nye2=ny2
      nzs2=1
      nze2=nz2-1

      nx3=mx3
      ny3=my3
      nz3=mz3
      nxs3=1
      nxe3=nx3
      nys3=1
      nye3=ny3
      nzs3=1
      nze3=nz3-1
      
      nx4=mx4
      ny4=my4
      nz4=mz4
      nxs4=1
      nxe4=nx4
      nys4=1
      nye4=ny4



c read grid info
c G1
      read(81) ig1, g1
      icount=0

      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.-2) then
      icount=icount+1
      read(82, 821) (ngp12(icount,m),m=1,3), (dr12(icount,m),m=1,3)
  821 format(3i5,3e20.12)
	end if 
      end do
      end do
      end do
      numngp12=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount

      icount=0
      do i=1,nx1
      do j=1,ny1
      do k=1,nz1
      if(ig1(i,j,k).eq.-3) then
      icount=icount+1
      read(83, 821) (ngp13(icount,m),m=1,3), (dr13(icount,m),m=1,3)
	end if 
      end do
      end do
      end do
      numngp13=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount

c G2
      read(84) ig2, g2, trans2
       icount=0
      do i=1,nx2
      do j=1,ny2
      do k=1,nz2
      if(ig2(i,j,k).eq.-1) then
      icount=icount+1
      read(85, 821) (ngp21(icount,m),m=1,3), (dr21(icount,m),m=1,3)
	end if 
      end do
      end do
      end do
      numngp21=icount      

      icount=0
      do i=1,nx2
      do j=1,ny2
      do k=1,nz2
      if(ig2(i,j,k).eq.-3) then
      icount=icount+1
      read(86, 821) (ngp23(icount,m),m=1,3), (dr23(icount,m),m=1,3)
	end if 
      end do
      end do
      end do

      numngp23=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount



c G3

      read(87) ig3, g3, trans3

      icount=0
      do i=1,nx3
      do j=1,ny3
      do k=1,nz3
      if(ig3(i,j,k).eq.-1) then
      icount=icount+1
      read(88, 821) (ngp31(icount,m),m=1,3), (dr31(icount,m),m=1,3)
	end if 
      end do
      end do
      end do
      numngp31=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount

      icount=0
      do i=1,nx3
      do j=1,ny3
      do k=1,nz3
      if(ig3(i,j,k).eq.-2) then
      icount=icount+1
      read(89, 821) (ngp32(icount,m),m=1,3), (dr32(icount,m),m=1,3)
	end if 
      end do
      end do
      end do

      numngp32=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount


c for G5 and G5b
      read(71) ig5, ngp5, drngp5, ig5b, ngp5b, drngp5b
      read(72) ngp15, drngp15, ngp25, drngp25, ngp35, drngp35
      read(73) ngp25b, drngp25b, ngp35b, drngp35b



c set up coefficients in grids 2 and 3
      igrid=2
      call bndryvdd(uc2bndry,g2, trans2,ig2,mx2,my2,mz2,
     1 dr2,ds2,dt2, nxs2,nxe2,nys2,nye2,nzs2,nze2,igrid)
      igrid=3
      call bndryvdd(uc3bndry,g3, trans3,ig3,mx3,my3,mz3,
     1 dr3,ds3,dt3, nxs3,nxe3,nys3,nye3,nzs3,nze3,igrid)

c G4
      read(94) xi0,dxi,ig4, g4
c
c      write(*,*) 'xi0, coshxio', xi0, cosh(xi0), ca, cb
c      stop
      
c from G2 for points at xi-dxi on Grid 4
      icount=0
      do i=1,nx4
      do j=1,ny4
      if(ig4(i,j).eq.-2) then
      icount=icount+1
      read(90, 821) (ngp42(icount,m),m=1,3), (dr42(icount,m),m=1,3)
	end if 
      end do
      end do
      
      numngp42=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount
c from G3 for points at xi-dxi on Grid 4
      icount=0
      do i=1,nx4
      do j=1,ny4
      if(ig4(i,j).eq.-3) then
      icount=icount+1
      read(91, 821) (ngp43(icount,m),m=1,3), (dr43(icount,m),m=1,3)
	end if 
      end do
      end do
      
      numngp43=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount
c
c from G2 for points at xi on Grid 4
      icount=0
      do i=1,nx4
      do j=1,ny4
      if(ig4(i,j).eq.-2) then
      icount=icount+1
      read(95, 821) (ngp42b(icount,m),m=1,3), (dr42b(icount,m),m=1,3)
	end if 
      end do
      end do
      
      numngp42b=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount
c from G3 for points at xi on Grid 4
      icount=0
      do i=1,nx4
      do j=1,ny4
      if(ig4(i,j).eq.-3) then
      icount=icount+1
      read(96, 821) (ngp43b(icount,m),m=1,3), (dr43b(icount,m),m=1,3)
	end if 
      end do
      end do
      
      numngp43b=icount      
      if(icount.gt.mngp) 
     1 write(*,*) 'trouble in interpolating points', icount
c


c G2 from G4
c from G2
      icount=0
      k=nz2-1
      do i=1,nx2
      do j=1,ny2
      if(ig2(i,j,k).eq.2) then
      icount=icount+1
      read(92, 821) (ngp24(icount,m),m=1,3), (dr24(icount,m),m=1,3)
	end if 
      end do
      end do
      numngp24=icount      

c
c G3 from G4
c from G3
      icount=0
      k=nz3-1
      do i=1,nx3
      do j=1,ny3
      if(ig3(i,j,k).eq.3) then
      icount=icount+1
      read(93, 821) (ngp34(icount,m),m=1,3), (dr34(icount,m),m=1,3)
	end if 
      end do
      end do
      numngp34=icount      
c


c Initial  state:
c read output data from fluid calculations and
c add B fields

c read output data from fluid calculations:
	  read(8) nts, tc, uc1, uc2, uc3
      write(*,*) 'nts= ', nts
	  close(8)


c add B
      igrid=1
      call setupb(ub1,ig1,g1,mx1,my1,mz1,
     1 nxs1,nxe1,nys1,nye1,nzs1,nze1,igrid)
      igrid=2
      call setupb(ub2,ig2,g2,mx2,my2,mz2,
     1 nxs2,nxe2,nys2,nye2,nzs2,nze2+1,igrid)
      igrid=3
      call setupb(ub3,ig3,g3,mx3,my3,mz3,
     1 nxs3,nxe3,nys3,nye3,nzs3,nze3+1,igrid)

c      igrid=4
      call setup4b(ub4,g4,mx4,my4,xi0,dxi)
      call setup5b(ub5,nphi+2,nnu+2,xi0,dphi,dnu)
      
       write(99) nts,tc,uc1,uc2,uc3, ub1, ub2, ub3, ub5

 1103 format(3f10.5,i3,3e16.6)
       
 1102 format(3f10.5,i3,4e13.4)

      stop
      end 





      subroutine bndryvdd(ucbndry,g,trans,ig,mx,my,mz,dr,ds,dt,
     1  nxs,nxe,nys,nye,nzs,nze, igrid)

c ucbndry(i,j,k,m)
c k=1,2,3: k=1 => k=nze-1 (inner point)
c          k=2  k=nze bndry point
c          k=3  k=nze+1 guard point
c m=1,2,3 vx,vy, vz
c m=4,5 coshksi, sinhksi
c m=6,7 costh, sinth
c m=8,9 cosphi sinphi


      parameter (mtm=4, mn=4)
      dimension ucbndry(mx,my,3,9), trans(mx,my,mz,3,3), ig(mx,my,mz)
      dimension g(mx,my,mz,3)
      dimension drdx(3,3), dudr(3,3), dudx(3,3), dpdr(3), dpdx(3)
      dimension d2uxdr(3,3,3), drxdr(3,3,3)
c ucbndry=vksi, vth, vphi (spheroidal) at 3 shells
c      dimension ucbndry(mx,my,3,3)
      common /var1/ omega, omegap, omegapx,omegapz,ca,cb,anu,eta,factor 

c ca: major radius (x, y); cb: minor radius (z)
c	
      cL=sqrt(ca**2-cb**2)
      cL2=cL*cL

      kk=0
      do 200 k = nze-1, nze+1
      kk=kk+1
      do 200 i = nxs, nxe
	do 200 j = nys, nye

      x = g(i,j,k,1)
      y = g(i,j,k,2)
      z = g(i,j,k,3)

c compute cosh(ksi), costh, and sinphi and cosphi
      rr=x*x+y*y
      r=sqrt(rr)
      zz=z*z
      bb=-(1+(rr+zz)/cL2)
      cc=rr/cL2
      Q=0.5* (-bb+sqrt(bb**2-4.*cc))
      coshksi=sqrt(Q)
      sinhksi=sqrt(Q-1.)
      sinth=r/(cL*coshksi)
      costh=z/(cL*sinhksi)
      if(r.ge.1.e-7) then
      cosphi=x/r
      sinphi=y/r
      else
      sinphi=0.
      cosphi=1.
      end if

c save them in ucbndry
      ucbndry(i,j,kk,4)=coshksi
      ucbndry(i,j,kk,5)=sinhksi
      ucbndry(i,j,kk,6)=costh
      ucbndry(i,j,kk,7)=sinth
      ucbndry(i,j,kk,8)=cosphi
      ucbndry(i,j,kk,9)=sinphi

200	continue


      return

      end



      subroutine setupb(uc,ig,g,mx,my,mz,
     1 nxs,nxe,nys,nye,nzs,nze,igrid)
      parameter (mtm=4, mn=4)
      dimension uc(mx,my,mz,mn), g(mx,my,mz,3), ig(mx,my,mz)
      common /var1/omega,omegap, omegapx,omegapz,ca,cb,anu,eta, factor 

c      factor=0.8
c      ca=1.
c      cb=ca*factor
      c2=cb**2
      pi=3.14159265359
      bbmag=1.e-10
      ak=pi

      do i=nxs, nxe 
      do j=nys, nye
      do k=nzs, nze
c      if(ig(i,j,k).eq.igrid) then
      x=g(i,j,k,1)
      y=g(i,j,k,2)
      z=g(i,j,k,3)
      rr=(x*x+y*y)
      rrr=sqrt(rr+z*z)
      if(rrr.eq.0.) then
      vx=2./3.*ak
      vy=0.
      vz=0.
      else

      rr=sqrt(rr)
      costh=z/rrr
      sinth=rr/rrr
      if(rr.ne.0.) then
      cosphi=x/rr
      sinphi=y/rr
      else
      cosphi=1.
      sinphi=0.
      end if
      akr=ak*rrr
      coskr=cos(akr)
      sinkr=sin(akr)
      vr=-2.*(akr*coskr-sinkr)*sinth*cosphi/(akr**2*rrr)
      vth=(akr*coskr-sinkr+akr**2*sinkr)*cosphi*costh/(akr**2*rrr)
      vphi=-(akr*coskr-sinkr+akr**2*sinkr)*sinphi/(akr**2*rrr)
      vx=sinth*cosphi*vr+costh*cosphi*vth-sinphi*vphi
      vy=sinth*sinphi*vr+costh*sinphi*vth+cosphi*vphi
      vz=costh*vr-sinth*vth

cc dipole
       r5=(sqrt(x**2+y**2+z**2))**5
c      if(r5.le.0.0001) r5=0.0001
c pointing to z direction
c      vx=3.*x*z/r5
c      vy=3.*y*z/r5
c      vz=(2.*z*z-x*x-y*y)/r5	
c pointing to x direction
c      vy=3.*y*x/r5
c      vz=3.*z*x/r5
c      vx=(2.*x*x-z*z-y*y)/r5	      
      end if
c here v->B pointing in x direction
c      uc(i,j,k,1)=vx*bbmag
c      uc(i,j,k,2)=vy*bbmag
c      uc(i,j,k,3)=vz*bbmag

c here v->B pointing in z direction
c      uc(i,j,k,1)=vy*bbmag
c      uc(i,j,k,2)=vz*bbmag
c      uc(i,j,k,3)=vx*bbmag
      
c here v->B pointing in z direction
      uc(i,j,k,1)=0.
      uc(i,j,k,2)=0.
      uc(i,j,k,3)=bbmag
      


c dipole field pole along z
c      r5=(sqrt(x*x+y*y+z*z))**5
c      if(r5.eq.0.) r5=0.01
c      vx=3.*x*z/r5
c      vy=3.*y*z/r5
c      vz=(2.*z*z-x*x-y*y)/r5	

c      uc(i,j,k,1)=vx*bbmag 
c      uc(i,j,k,2)=vy*bbmag
c      uc(i,j,k,3)=vz*bbmag
      
c      end if
      end do
      end do
      end do
      return
      end

      subroutine setup5b(uc,mx,my,xi0,dphi,dth)
      
      dimension uc(mx,my)
      common /var1/omega,omegap, omegapx,omegapz,ca,cb,anu,eta, factor 
ccc
c only for the normal component
ccc


c dipole field
c     B=(3 (m.rr) rr -m)/r^3;  m=zz if pointing in z direction
c                              rr unit vector; rr=(xx,yy,zz)
c                              r=Radius
c      Bx=3.*x*z/r5
c      By=3.*y*z/r5
c      Bz=(2.*z*z-x*x-y*y)/r5	

c
      bbmag=1.e-10

      pi=3.14159265359
      ak=pi
      coshxi=cosh(xi0)
      xmax=cosh(xi0)
      
      sinhxi=sinh(xi0)
      sinhxi2=sinh(xi0)**2


      do j=1, my
      do i=1, mx
      
c      if(ig(i,j,k).eq.igrid) then
      th=(float(j)-1.5)*dth     
      phi=float(i-2)*dphi
      sinth=sin(th)
      costh=cos(th)
      sinphi=sin(phi)
      cosphi=cos(phi)
	x=cosh(xi0)/xmax*sinth*cosphi
	y=cosh(xi0)/xmax*sinth*sinphi
	z=sinh(xi0)/xmax*costh
	
      rr=(x*x+y*y)
      rrr=sqrt(rr+z*z)
      
      if(rrr.eq.0.) then
      vx=2./3.*ak
      vy=0.
      vz=0.
      else
      rr=sqrt(rr)
      costh=z/rrr
      sinth=rr/rrr
      if(rr.ne.0.) then
      cosphi=x/rr
      sinphi=y/rr
      else
      cosphi=1.
      sinphi=0.
      end if
cc
cc      rrr=1.
cc            
      akr=ak*rrr
      coskr=cos(akr)
      sinkr=sin(akr)
      vr=-2.*(akr*coskr-sinkr)*sinth*cosphi/(akr**2*rrr)
      vth=(akr*coskr-sinkr+akr**2*sinkr)*cosphi*costh/(akr**2*rrr)
      vphi=-(akr*coskr-sinkr+akr**2*sinkr)*sinphi/(akr**2*rrr)
      vx=sinth*cosphi*vr+costh*cosphi*vth-sinphi*vphi
      vy=sinth*sinphi*vr+costh*sinphi*vth+cosphi*vphi
      vz=costh*vr-sinth*vth

c      Bx=3.*x*z/r5
c      By=3.*y*z/r5
c      Bz=(2.*z*z-x*x-y*y)/r5	
cc
c     r5=(sqrt(x**2+y**2+z**2))**5
c pointing to z direction
c      vx=3.*x*z/r5
c      vy=3.*y*z/r5
c      vz=(2.*z*z-x*x-y*y)/r5	
c pointing to x direction
c      vy=3.*y*x/r5
c      vz=3.*z*x/r5
c      vx=(2.*x*x-z*z-y*y)/r5	

      
      end if

c      write(21,11211) i,j,vr, vth,vphi
c11211 format(2i,3e13.4)
c here v->B in x direction
      bx=vx*bbmag
      by=vy*bbmag
      bz=vz*bbmag

c here v->B pointing in z direction
c      bx=vy*bbmag
c      by=vz*bbmag
c      bz=vx*bbmag

c here v->B pointing in z direction
      bx=0.
      by=0.
      bz=bbmag

c dipole field pole along z
      r5=(sqrt(x*x+y*y+z*z))**5
      if(r5.eq.0.) r5=0.01
c pointing to z direction
c      vx=3.*x*z/r5
c      vy=3.*y*z/r5
c      vz=(2.*z*z-x*x-y*y)/r5	
c pointing to x direction
c      vy=3.*y*x/r5
c      vz=3.*z*x/r5
c      vx=(2.*x*x-z*z-y*y)/r5 	      
c      bx=vx*bbmag
c      by=vy*bbmag
c      bz=vz*bbmag





c now into spheroidal coordinate

	costh=cos(th)
      sinth=sin(th)
      cosphi=cos(phi)
      sinphi=sin(phi)
      gg=1./sqrt(costh**2+sinhxi2)

      a1=gg*sinhxi*sinth*cosphi
 	b1=gg*sinhxi*sinth*sinphi
	c1=gg*coshxi*costh
      a2=gg*coshxi*costh*cosphi
	b2=gg*coshxi*costh*sinphi
	c2=-gg*sinhxi*sinth
      a3=-sinphi
      b3=cosphi
      c3=0.
c here br=bxi
	br=a1*bx+b1*by+c1*bz
	bth=a2*bx+b2*by+c2*bz
      bphi=a3*bx+b3*by
      
      uc(i,j)=br
c      uc(i,j,1,2)=bth
c      uc(i,j,1,3)=bphi

c      uc(i,j,1,1)=vr
c      uc(i,j,1,2)=vth
c      uc(i,j,1,3)=vphi
      


c      uc(i,j,1,1)=bx
c      uc(i,j,1,2)=by
c      uc(i,j,1,3)=bz
c      if(i.le.4.and.j.le.4) then
c      write(*,'(2i4,5f10.5,3e13.4)') i,j,th,phi,x,y,z,bx,by,bz
c      end if

      end do
      end do
      
      

      return
      end
      


      subroutine setup4b(uc,g,mx,my,xi0,dxi)
      
      dimension uc(mx,my,2,6), g(mx,my,2)
      common /var1/omega,omegap, omegapx,omegapz,ca,cb,anu,eta, factor 

c
c
c dipole field
c     B=(3 (m.rr) rr -m)/r^3;  m=zz if pointing in z direction
c                              rr unit vector; rr=(xx,yy,zz)
c                              r=Radius
c      Bx=3.*x*z/r5
c      By=3.*y*z/r5
c      Bz=(2.*z*z-x*x-y*y)/r5	

c
      bbmag=1.e-10

      pi=3.14159265359
      ak=pi
      coshxi=cosh(xi0)
      xmax=cosh(xi0)
      
      sinhxi=sinh(xi0)
      sinhxi2=sinh(xi0)**2

c      do i=1, mx
      do j=1, my
      do i=1, mx
      
c      if(ig(i,j,k).eq.igrid) then
      th=g(i,j,1)     
      phi=g(i,j,2)
      sinth=sin(th)
      costh=cos(th)
      sinphi=sin(phi)
      cosphi=cos(phi)
	x=cosh(xi0)/xmax*sinth*cosphi
	y=cosh(xi0)/xmax*sinth*sinphi
	z=sinh(xi0)/xmax*costh
	
      rr=(x*x+y*y)
      rrr=sqrt(rr+z*z)
      
      if(rrr.eq.0.) then
      vx=2./3.*ak
      vy=0.
      vz=0.
      else
      rr=sqrt(rr)
      costh=z/rrr
      sinth=rr/rrr
      if(rr.ne.0.) then
      cosphi=x/rr
      sinphi=y/rr
      else
      cosphi=1.
      sinphi=0.
      end if
cc
cc      rrr=1.
cc            
      akr=ak*rrr
      coskr=cos(akr)
      sinkr=sin(akr)
      vr=-2.*(akr*coskr-sinkr)*sinth*cosphi/(akr**2*rrr)
      vth=(akr*coskr-sinkr+akr**2*sinkr)*cosphi*costh/(akr**2*rrr)
      vphi=-(akr*coskr-sinkr+akr**2*sinkr)*sinphi/(akr**2*rrr)
      vx=sinth*cosphi*vr+costh*cosphi*vth-sinphi*vphi
      vy=sinth*sinphi*vr+costh*sinphi*vth+cosphi*vphi
      vz=costh*vr-sinth*vth

c      Bx=3.*x*z/r5
c      By=3.*y*z/r5
c      Bz=(2.*z*z-x*x-y*y)/r5	
cc
c     r5=(sqrt(x**2+y**2+z**2))**5
c pointing to z direction
c     vx=3.*x*z/r5
c     vy=3.*y*z/r5
c     vz=(2.*z*z-x*x-y*y)/r5	
c pointing to x direction
c      vy=3.*y*x/r5
c      vz=3.*z*x/r5
c      vx=(2.*x*x-z*z-y*y)/r5	

      
      end if

c      write(21,11211) i,j,vr, vth,vphi
c11211 format(2i,3e13.4)
c here v->B in x direction
c      bx=vx*bbmag
c      by=vy*bbmag
c      bz=vz*bbmag

c here v->B pointing in z direction
c      bx=vy*bbmag
c      by=vz*bbmag
c      bz=vx*bbmag

c here v->B pointing in z direction
      bx=0.
      by=0.
      bz=bbmag


c now into spheroidal coordinate

      th=g(i,j,1)     
      phi=g(i,j,2)

	costh=cos(th)
      sinth=sin(th)
      cosphi=cos(phi)
      sinphi=sin(phi)
      gg=1./sqrt(costh**2+sinhxi2)

      a1=gg*sinhxi*sinth*cosphi
 	b1=gg*sinhxi*sinth*sinphi
	c1=gg*coshxi*costh
      a2=gg*coshxi*costh*cosphi
	b2=gg*coshxi*costh*sinphi
	c2=-gg*sinhxi*sinth
      a3=-sinphi
      b3=cosphi
      c3=0.
c here br=bxi
	br=a1*bx+b1*by+c1*bz
	bth=a2*bx+b2*by+c2*bz
      bphi=a3*bx+b3*by
      
      uc(i,j,1,1)=br
      uc(i,j,1,2)=bth
      uc(i,j,1,3)=bphi

c      uc(i,j,1,1)=vr
c      uc(i,j,1,2)=vth
c      uc(i,j,1,3)=vphi
      


c      uc(i,j,1,1)=bx
c      uc(i,j,1,2)=by
c      uc(i,j,1,3)=bz
c      if(i.le.4.and.j.le.4) then
c      write(*,'(2i4,5f10.5,3e13.4)') i,j,th,phi,x,y,z,bx,by,bz
c      end if

      end do
      end do

      return
      end