c
c
c    requires input files "eqb1" from equilibrium code
c    if solovev equilibrium is chosen the input file is not necessary
c    it generates output files "mapdsk" & "mpout1" for stability calculation
c    diagnostic output are 'output', 'f3dmap0x'
c
c      ... use parameters  ...   ns for number of surfaces - nsf0
c                                nt0 for number of thetas   - mth
c           to match input from equilibrium
c                                nq for number of thetas   > nthe
c                                ny for number of surfaces > npsi
c       nthe is the number of theta grid pts from the equilibrium code
c       npsi is the number of radial grid pts from the equilibrium code
c     from this mapping code, mthi is the number of theta grid pts
c     from this mapping code, nosuri is the number of radial grid pts
c       nosuri must be smaller than nsrf
ckg  Modified by N.N.Gorelenkov to make an output in portable cdf format. 
c    In addition to two files it now creates one mapdsk.cdf file. (2003)
ckg  Another modification is the grid accumulation points if igrid=4 
c    in one of the namelists. Works fine for NOVA and NOVA-KN. (2006)
c     program map(tape6=output)
      use ezcdf
ckg_cdf
      character*80 eqfile
c      character*40 work_cdf
      character*4 vartype
c      character*8 ntitle_cdf(20)
      character*80 Title
c
      integer*4 cdfid, dimlens(3), mth_cdf
c
      include 'clichpar.h'
      real*8 work_cdf(nt,ns),zmin,zmax,xmax,xmin, psibig_cdf(ns),
     &     p_cdf(ns),pp_cdf(ns),q_cdf(ns),qp_cdf(ns),x_cdf(nt,ns)
     &     ,z_cdf(nt,ns),xdth_cdf(nt,ns),zdth_cdf(nt,ns),xpsi_cdf(nt,ns)
     &     ,zpsi_cdf(nt,ns),grpssq_cdf(nt,ns),grthsq_cdf(nt,ns)
     &     ,grpsth_cdf(nt,ns),gptdth_cdf(nt,ns),gpsdth_cdf(nt,ns)
     &     ,xjprime_cdf(nt,ns),delta_cdf(nt,ns)
ckg
      integer*8 nadres,kshift,length,ierr,len
      integer*8 ntitle,nxx,nxy
ckg_cdf
      real*8 axx,axy
      common ntitle(20),dat,nxx(5),axx(13),nxy(10),axy(10)
      common/var/ p(ny),pp(ny),q(ny),qp(ny),
     1 g(ny),gp(ny),fb(ny),fbp(ny),f(ny),fp(ny),
     1 psival(ny),x(nq,ny),z(nq,ny),aj3(nq,ny),aj(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var2/dtst(nq,ny),dpts(nq,ny)
      common/var3/fg(nsgrd),psig(nsgrd),fd(ny),psivad(ny),fpd(ny)
      common/var4/tb(nq,ns)
      common/arry1/work1(ns),work2(nt,ns),work3(nt,ns)
     1,gg(ns),ggp(ns),workx(ns),workz(ns),xinf(nt),zinf(nt)
      common/arry2/w1(ny),w2(nq,ny),w3(nq,ny)
      common/arry3/dtsx(nq,ny),dtsz(nq,ny),gps(nq,ny),t(nq,ny),
     1dpst(nq,ny),tp(nq,ny),dttp(nq,ny),dtx(nq,ny),dtz(nq,ny),
     1dpx(nq,ny),dpz(nq,ny),xjacob(nq,ny),delta(nq,ny),
     1xsq(nq,ny),xmap(nt,ns),zmap(nt,ns),bmod(nq,ny)
      common/arry4/dpxsq(nq,ny),dpxj(nq,ny),dpqdel(nq,ny)
      character*8 hedr(2),hdr,ekdsk
      common / head / hedr
      character*1 pqkey
c
ckg      namelist /hdr/ hdrn
      namelist/wdat/nsf,nths,nosuri,mthi,njacg,mjacx,nptsi,jndexi
     &,igrid,ljacb,ekdsk,n00,nnos,nos,rs,xm1,gdxm1,d1f1,xm2,gdxm2
     &     ,d1f2,xm3,gdxm3,d1f3
c     Replaced
c      call link ( 'unit6=(output,create,text),
c     .             unit59 = terminal//' )
c
ckg   check the consistency of the grid dimension in particular the
c     one used for psig. If it is smaller than ns or ny stop:
      if(nsgrd.lt.ns.or.nsgrd.lt.ny) stop
     &     'WARNING check the dimension of nsgrd vs ns and ny'
ckg      istat=ishell('rm -f output')
ckg
      open(16,file='mapoutput',status='unknown',form='formatted')
ckg      open(101,file='mapprint',status='unknown',form='formatted')
c     Replace input from 59 with *(100-stdin) and output with 101(stdout)
c
      go to 911
ckg      write ( 101, 800 )
c      write ( *, 800 )
c  800 format ( 1x, 'type header for mapping - hdrn = 8h...' )
c      write(*,*) "HDR NAMELIST"
c      read ( *, hdr )
c      hedr(1) = 'box p45'
c      hedr(2) = hdrn
ckg      call setbox ( hedr, 16 )
  911 continue
c
c     Replaced:
c      call gfsize(3,500000b)
c      call keep80('dmap',3)
c      call fr80id('dm map',1,1)
      call ncarcgm(1,'plotmap')

c     call ddi(wdat,12,6,1)
c
c.......  use namelist from tty.
c
c........ jacobian proportional to x**m  / ( grdpsi**n * b**l )
c
      call anfun(lanal)
      if(lanal.eq.1) go to 910
c
      ekdsk='eqb1'
c     Replaced call zop(inp1,ekdsk,nsz2,id2,st,999)
      write(*,*)
     &     'What equilibrium P(EST) or Q(-solver) will be maped'
      read(*,'(a)')pqkey
      kshift=0
      if(pqkey.eq.'p') kshift=1 
      write(*,*)'1',ny,nq
      length=(ny*11+nq*ny*4+2)/512
      nadres=0
      call wopen("eqb1",length,nadres,ierr)
      lgivup=1
      nadres=1+kshift
      length=49
c     Replaced call zrd(inp1,ntitle(1),length,nadres,lgivup,999)
      call getwa("eqb1",ntitle(1),nadres,length,ierr)
      write(*,'(17a,20A8)') "Kg check ntitle:",ntitle
c      write(*,*) nxx,dat,axx
      nadres=50+kshift
      length=10
      call getwa("eqb1",axy(1),nadres,length,ierr)
cntitle(20),dat,nxx(5),axx(13),nxy(10),axy(10)
      write(*,*)'dat',dat
      write(*,*)'nxx',nxx
      write(*,*)'axx',axx
      write(*,*)'nxy',nxy
      write(*,*)'axy',axy
      nadres=60+kshift
      length=ny*11+nq*ny*4
c     Replaced call zrd(inp1,p(1),length,nadres,lgivup,999)
      call getwa("eqb1",p(1),nadres,length,ierr)
c      nadres=60+kshift+kshift+length
c      length=nq*ny*4
c      call getwa("eqb1",x(1,1),nadres,length,ierr)
c these are read
c p(ny),pp(ny),q(ny),qp(ny),
c g(ny),gp(ny),fb(ny),fbp(ny),f(ny),fp(ny),
c psival(ny),x(nq,ny),z(nq,ny),aj3(nq,ny),aj(nq,ny)
      nadres=nadres+nq*ny
      length=1
      call getwa("eqb1",bp2sq,nadres,length,ierr)
      nadres=nadres+length
c      write(*,*) 'bp2sq',bp2sq
c      stop '1 stop in dmapb2.f'
      call getwa("eqb1",volume,nadres,length,ierr)
      nadres=nadres+length
ccc
      nthe=nxx(1)
      npsi=nxx(2)
      npsim=npsi-1
      nthe1=nthe+1
      nthe2=nthe1+1
      nthe3=nthe2+1
      nthe4=nthe3+1
      npsip=npsi+1
ccc
ccc  renormalized AJ & AJ3 so that AJ's are the actual Jacobian
      do 40 is=2,npsi
      do 40 it=1,nthe4
      aj3(it,is)=aj3(it,is)/((f(is-1)+f(is))/2.0)
      aj(it,is)=aj(it,is)/((f(is-1)+f(is))/2.0)
   40 continue
      do 41 it=2,nthe3
      aj(it,1)=aj(it,1)/f(1)*2.0
   41 continue
      do 42 it=3,nthe2
      aj3(it,1)=(aj(it,1)+aj(it-1,1))*0.5
   42 continue
c
cc  extroplate aj, aj3 to very near magnetic axis
      nptsi=5
      jndexi=1
      mth=nthe2
c     call exi2(aj,3,0)
c     call exi2(aj3,3,0)
ckg
      do 600 j=1,npsip
      psivad(j)=psival(j)
      fpd(j)=fp(j)
      fd(j)=f(j)
600   continue
c
cc  define grid points for computing quantities on the original equilibrium grid points
      npsip2=npsi+2
      call grid(1,npsi,1,2,2,0.,xm1,gdxm1,d1f1,xm2,gdxm2
     &     ,d1f2,xm3,gdxm3,d1f3)
      call splprp(npsi)
      call splprp1(npsi,npsip2,psig)
c
  910 continue
c
c
ccc   input parameters for mapping
ccc     use cubic spline for interpolation
c
      nsf=ns
      nths=nt
      nosuri=nmbrgrd
      mthi=mnths0
      igrid=2
      nptsi=5
      jndexi=1
      mjacx=1
      njacg=1
      ljacb=0
      n00=4
      nos=2
      nnos=5
      rs=.5
ckg This input is for mesh accumulation if igrid=4 at up to three points 
ckg in radius they are followed as xm{1,2,3} with default width gdxm{1,2,3} and 
ckg density of accumulation d1f{1,2,3}, which is between 0.1 and 1. By default
ckg we have only one accumulation point. In the namelist that can
ckg be changed. (N.N. Gorelenkov 2006)
      xm1=0.6
      gdxm1=0.1
      d1f1=0.4
      xm2=0.15
      gdxm2=0.
      d1f2=1.
      xm3=0.8
      gdxm3=0.
      d1f3=1.
c
      do 211 k0=1,npsi
      if(q(k0).le.1.0.and.q(k0+1).gt.1.0) go to 212
      go to 211
  212 continue
      rs=(k0+(1.0-q(k0+1))/(q(k0+1)-q(k0)))/(npsi-1.)
  211 continue
c
ccc... options for choosing grids
ccc    igrid=1; uniform sqrt(psi) grid
ccc     igrid=2 uniform psi grid
ccc     igrid=3 tailored grid with dense mesh at q=1 surface
ckg     igrid=4 see right above
      write ( *, 801 )nmbrgrd,mthi
ckg      write ( 101, 801 )
 801  format (1x,'nsf=ns nths=nt, nosuri=',i4,', mthi=',i3,/,
     .  1x, ', igrid=2 (Psi), (choose igrid=1 for sqrt(Psi))',/,
     .  1x, 'nptsi=5 jndexi=1 mjacx=1 njacg=1 ljacb=0 ekdsk=eqb1',/,
     .  1x, 'nnos=5 n00=4 nos=2 rs=.5 and type in desired value',//)
c
      write ( *, wdat )
      write(*,*) "WDAT NAMELIST"
      read ( *, wdat )
      write ( 16, wdat )
c
      npi = njacg
      mjac = mjacx
      ljac = ljacb
c
      np=npi
      nxy(1)=np
      nxy(2)=mjac
      nxy(3) = ljac
      dth=axy(1)
      dr=axy(2)
      pi=axy(3)
      xzero = axx(5)
      mth=nxx(4)
      nosurf=nosuri
      mth=mthi
      nxx(3)=nosurf
      nxx(4)=mth
      nxx(5)=igrid
      mth1 = mth + 1
      mth2=mth+2
ckg      do 600 j=1,npsip
c      psivad(j)=psival(j)
c      fpd(j)=fp(j)
c      fd(j)=f(j)
c600   continue
c
      do 601 j=1,npsi
      do 601 i=1,nthe4
      z(i,j)=-z(i,j)
601   continue
      call dts(x,dtsx,-1)
      call dts(z,dtsz,+1)
      do 1 j=2,npsi
      do 1 i=3,nthe2
      ggrd = ( g(j) + g(j-1) ) / 2.0
      gps(i,j)=(dtsx(i,j)**2
     1+dtsz(i,j)**2)*x(i,j)**2/aj(i,j)**2
      bmod(i,j) = sqrt ( gps(i,j) + (xzero*ggrd)**2 ) / x(i,j)
1     continue
c
      do 111 i=1,nthe4
      gps(i,1)=0.0
      bmod(i,1)=g(1)
  111 continue
c
c
      call sym(gps,+1)
      call sym ( bmod, +1 )
c
      do 32 j=2,npsi
      const=0.
      do 33 i=4,nthe2
c     gpsh=(gps(i,j)+gps(i-1,j))/2.
      gpsh=(sqrt(gps(i,j))+sqrt(gps(i-1,j)))/2.
      bmodh = ( bmod(i,j) + bmod(i-1,j) ) / 2.0
      const=const + bmodh**ljac * aj3(i,j)*dth/
     1(((x(i,j)+x(i-1,j))**mjac)/
     1(gpsh**np))
33    continue
      t(3,j)=0.
      do 34 i=4,nthe2
c     gpsh=(gps(i,j)+gps(i-1,j))/2.
      gpsh=(sqrt(gps(i,j))+sqrt(gps(i-1,j)))/2.
      bmodh = ( bmod(i,j) + bmod(i-1,j) ) / 2.0
      t(i,j)=t(i-1,j) + bmodh**ljac * (aj3(i,j)*dth/
     1(((x(i,j)+x(i-1,j))**mjac)/
     1(gpsh**np)))/(const/pi)
34    continue
32    continue
c     pause 1
c
      do 2 j=1,npsi
      constp=0.
      do 3 i=4,nthe2
c     gpsh=(gps(i,j)+gps(i-1,j))/2.
      gpsh=(sqrt(gps(i,j))+sqrt(gps(i-1,j)))/2.
      bmodh = ( bmod(i,j) + bmod(i-1,j) ) / 2.0
      if(j.eq.1.and.kshift.eq.1) then 
         constp=constp+aj3(i,2)*dth/
     1        ((x(i,j)+x(i-1,j))**2)
      else
         constp=constp+aj3(i,j)*dth/
     1        ((x(i,j)+x(i-1,j))**2)
      endif
3     continue
      t(3,j)=0.
      tp(3,j)=0.
      do 4 i=4,nthe2
c     gpsh=(gps(i,j)+gps(i-1,j))/2.
      gpsh=(sqrt(gps(i,j))+sqrt(gps(i-1,j)))/2.
      bmodh = ( bmod(i,j) + bmod(i-1,j) ) / 2.0
      if(j.eq.1.and.kshift.eq.1) then 
         tp(i,1)=tp(i-1,2)+(aj3(i,2)*dth/
     1        ((x(i,j)+x(i-1,j))**2))/(constp/pi)
      else
         tp(i,j)=tp(i-1,j)+(aj3(i,j)*dth/
     1        ((x(i,j)+x(i-1,j))**2))/(constp/pi)
      endif
      if(j.eq.1) t(i,1)=t(i,2)
      if(ljac.eq.0.and.mjac.eq.1.and.np.eq.1) t(i,j)=(i-3.)*dth
4     continue
ccc   note t and tp do not have parity
      t(2,j)=-t(4,j)
      t(nthe3,j)=2.*t(nthe2,j)-t(nthe1,j)
      tp(2,j)=-tp(4,j)
      tp(nthe3,j)=2.*tp(nthe2,j)-tp(nthe1,j)
2     continue
c     pause 2
c
      call dts(t,dtst,+1)
      call dps(t,dpst,-1)
      do 7 j=1,npsi
      do 7 i=2,nthe3
      dpts(i,j)=-dpst(i,j)/dtst(i,j)
7     continue
      call dt(tp,dttp,+1)
c
      call dp(x,dpx,+1)
      call dp(z,dpz,-1)
      do 504 j=1,npsi
      do 504 i=3,nthe2
      xsq(i,j)=x(i,j)**2
504   continue
      call sym(xsq,+1)
      do 505 j=1,npsi
      do 505 i=3,nthe2
      delta(i,j)=tp(i,j)-t(i,j)
      xjacob(i,j)=(dttp(i,j)*xsq(i,j))/fb(j)
      if(ljac.eq.0.and.mjac.eq.1.and.np.eq.1) xjacob(i,j)=aj(i,j)
505   continue
      call sym(delta,-1)
      call sym(xjacob,+1)
      call dp(xsq,dpxsq,+1)
      call dp(xjacob,dpxj,+1)
      call sym(dpxj,+1)
      do 509 j=1,npsi
      do 509 i=3,nthe2
      w2(i,j)=q(j)*delta(i,j)
509   continue
      call dp(w2,dpqdel,-1)
      call sym(dpqdel,-1)
ccc....   define the new mapping grids
      call grid(igrid,nosurf,n00,nnos,nos,rs,xm1,gdxm1,d1f1,xm2,gdxm2
     &     ,d1f2,xm3,gdxm3,d1f3)
      npsip2=npsi+2
c
      call splprp(npsi)
      call splprp1(nosurf,npsip2,psig)
c
      write(16,1001)
 1001 format('  xjacob near magnetic axis ')
      write(16,1000)(aj(i,1),i=1,nthe4)
      write(16,1000)(aj(i,2),i=1,nthe4)
1000  format(1x,10e12.4)
c
ccccccc........
ccc  in intp2s tb is defined for j=1,nosurf
      call intp2s(t,tb,0)
      do 8 j=1,nosurf
      tb(2,j)=-tb(4,j)
      tb(nthe3,j)=2.*tb(nthe2,j)-tb(nthe1,j)
8     continue
c     pause 3
      len=50+9*nosurf+6*nsf*nths+5
      idisk=1
c     Replaced call zowc(inp2,6rmpout1,len,idisk,ff,999)
ckg      istat=ishell('rm -f mpout1')
      length=len/512
      nadres=0
      call wopen("mpout1",length,nadres,ierr)
      length=49
      lgivup=1
      nadres=1
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,ntitle(1),length,nadres,lgivup,999)
      call putwa("mpout1",ntitle(1),nadres,length,ierr)
      length=nosurf
      nadres=50
      call intp1(p,p_cdf,1)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
ccc since p(nosurf) was not calculated correctly in equilibrium code, set
      p_cdf(nosurf)=0.01*p_cdf(nosurf-1)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",p_cdf,nadres,length,ierr)
      write(16,4200)
      call plotx(psig,p_cdf,work2,0,1,1,1)
      write(string,4200)
      call gtext(string,-1,-1)
4200  format('  p')
c
      nadres=nadres+length
      call intp1(pp,pp_cdf,0)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",pp_cdf,nadres,length,ierr)
      write(16,4201)
      call plotx(psig,pp_cdf,work2,0,1,1,2)
      write(string,4201)
      call gtext(string,-1,-1)
4201  format('  pp')
c
      nadres=nadres+length
      call intp1(q,q_cdf,0)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",q_cdf,nadres,length,ierr)
      write(16,4202)
      call plotx(psig,q_cdf,work2,0,1,1,3)
      write(string,4202)
      call gtext(string,-1,-1)
4202  format('  q')
c
      nadres=nadres+length
      call intp1(qp,qp_cdf,0)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",qp_cdf,nadres,length,ierr)
      write(16,4203)
      call plotx(psig,qp_cdf,work2,0,1,1,4)
      write(string,4203)
      call gtext(string,-1,-1)
      call framep
4203  format('  qp')
c
      nadres=nadres+length
      call intp1(g,gg,1)
ccc   call exi1(gg,1,4)
ccc   call exo1(gg,nosurf)
c     Replaced call zwr(inp2,gg(1),length,nadres,lgivup,999)
      call putwa("mpout1",gg(1),nadres,length,ierr)
      write(16,4204)
      call plotx(psig,gg,work2,0,1,1,1)
      write(string,4204)
      call gtext(string,-1,-1)
4204  format('  g')
c
c
      nadres=nadres+length
      call intp1(gp,ggp,0)
ccc   call exi1(ggp,1,4)
ccc   call exo1(ggp,nosurf)
c     Replaced call zwr(inp2,ggp(1),length,nadres,lgivup,999)
      call putwa("mpout1",ggp(1),nadres,length,ierr)
      write(16,4205)
      call plotx(psig,ggp,work2,0,1,1,2)
      write(string,4205)
      call gtext(string,-1,-1)
4205  format('  gp')
c
      nadres=nadres+length
      call intp1(fb,work1,0)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",work1,nadres,length,ierr)
      write(16,4206)
      call plotx(psig,work1,work2,0,1,1,3)
      write(string,4206)
      call gtext(string,-1,-1)
4206  format('  fb')
c
      nadres=nadres+length
      call intp1(fbp,work1,0)
ccc   call exi1(work1,1,4)
ccc   call exo1(work1,nosurf)
c     Replaced call zwr(inp2,work1,length,nadres,lgivup,999)
      call putwa("mpout1",work1,nadres,length,ierr)
      write(16,4207)
      call plotx(psig,work1,work2,0,1,1,4)
      write(string,4207)
      call gtext(string,-1,-1)
      call framep
4207  format('  fbp')
c
cccxxxxxxx  This is different from standard PEST "mpout1" file
c  add "psig" grid points to the "mpout1" to take into account variable grid information
      nadres=nadres+length
c     Replaced call zwr(inp2,psig,nosurf,nadres,lgivup,999)
      length=nosurf
      call putwa("mpout1",psig,nadres,length,ierr)
ccccxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
      nadres=nadres+length
      length=nsf*nths
      call int2spl(x,work2,0,+1,0)
c
      do 5100 j = 1, nosurf
      do 5100 i = 1, mth2
      xmap(i,j) = work2(i,j)
      x_cdf(i,j) = work2(i,j)
 5100 continue
ccccc
c
c     Replaced call zwr(inp2,work2,length,nadres,lgivup,999)
      call putwa("mpout1",work2,nadres,length,ierr)
      write(16,4208)
      call plotx(psig,work1,work2,0,1,2,1)
      write(string,4208)
      call gtext(string,-1,-1)
4208  format('  x')
      nadres=nadres+length
      do 501 i=1,mth2
      xinf(i)=work2(i,nosurf)
501   continue
      call int2spl(z,work2,0,-1,0)
      do 5200 i = 1, mth2
      do 5200 j = 1, nosurf
      zmap(i,j) = work2(i,j)
      z_cdf(i,j) = work2(i,j)
 5200 continue
c
c     Replaced call zwr(inp2,work2,length,nadres,lgivup,999)
      call putwa("mpout1",work2,nadres,length,ierr)
      write(16,4209)
      call plotx(psig,work1,work2,0,2,2,2)
      write(string,4209)
      call gtext(string,-1,-1)
4209  format('  z')
      nadres=nadres+length
      do 502 i=1,mth2
      zinf(i)=work2(i,nosurf)
502   continue
c     pause 4
c
      call dt(x,dtx,-1)
      call int2spl(dtx,xdth_cdf,0,-1,0)
c     Replaced call zwr(inp2,xdth_cdf,length,nadres,lgivup,999)
      call putwa("mpout1",xdth_cdf,nadres,length,ierr)
      write(16,4210)
      call plotx(psig,work1,xdth_cdf,0,2,2,3)
      write(string,4210)
      call gtext(string,-1,-1)
4210  format('  xdth')
c
      nadres=nadres+length
      call dt(z,dtz,+1)
      call int2spl(dtz,zdth_cdf,0,+1,0)
c     Replaced call zwr(inp2,work2,length,nadres,lgivup,999)
      call putwa("mpout1",zdth_cdf,nadres,length,ierr)
      write(16,4211)
      call plotx(psig,work1,zdth_cdf,0,1,2,4)
      write(string,4211)
      call gtext(string,-1,-1)
      call framep
4211  format('  zdth')
ccc  plot flux surfaces
      call plotxz ( xmap,zmap, nt,ns, mth1,nosurf, 4,5 )
      call framep
c
      nadres=nadres+length
      call int2spl(dpx,xpsi_cdf,0,+1,0)
      call dpsinor(xpsi_cdf)
c     Replaced call zwr(inp2,work2,length,nadres,lgivup,999)
      call putwa("mpout1",xpsi_cdf,nadres,length,ierr)
      write(16,4212)
      call plotx(psig,work1,xpsi_cdf,1,1,2,1)
      write(string,4212)
      call gtext(string,-1,-1)
4212  format('  xpsi note ln scale and +-')
c
      nadres=nadres+length
      call int2spl(dpz,zpsi_cdf,0,-1,0)
      call dpsinor(zpsi_cdf)
c     Replaced call zwr(inp2,work2,length,nadres,lgivup,999)
      call putwa("mpout1",zpsi_cdf,nadres,length,ierr)
      write(16,4213)
      call plotx(psig,work1,zpsi_cdf,1,2,2,2)
      write(string,4213)
      call gtext(string,-1,-1)
4213  format('  zpsi note ln scale and +-')
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      len=(12*nsf*nths+nosurf+2*mth2+50+5)/512
c     Replaced call zowc(inp3,6rmapdsk,len,idisk,ff,999)
ckg      istat=ishell('rm -f mapdsk')
      nadres=0
      call wopen("mapdsk",len,nadres,ierr)
      length=49
      nadres=1
c     Replaced call zwr(inp3,ntitle(1),length,nadres,lgivup,999)
      call putwa("mapdsk",ntitle(1),nadres,length,ierr)
      length=mth2
      nadres=50
c     Replaced call zwr(inp3,xinf(1),length,nadres,lgivup,999)
      call putwa("mapdsk",xinf(1),nadres,length,ierr)
      nadres=nadres+length
c     Replaced call zwr(inp3,zinf(1),length,nadres,lgivup,999)
      call putwa("mapdsk",zinf(1),nadres,length,ierr)
c
      nadres=nadres+length
      do 503 j=1,npsi
      w1(j)=psival(j)*2.*pi
503   continue
      call intp1(w1,work1,0)
ckg_cdf
      xmin=1000.
      xmax=0.
      zmin=1000.
      zmax=-1000.
      do i=1,mth2
         work_cdf(i,1)=xinf(i)
         work_cdf(i,2)=zinf(i)
         zmin=min(zmin,zinf(i))
         zmax=max(zmax,zinf(i))
         xmin=min(xmin,xinf(i))
         xmax=max(xmax,xinf(i))
      enddo
      psibig_cdf(nosurf)=psival(1)*2.*pi
      do j=1,nosurf
         psibig_cdf(j)=work1(j)-psibig_cdf(nosurf)
      enddo
ckg_cdf
      length=nosurf
c     Replaced call zwr(inp3,work1,length,nadres,lgivup,999)
      call putwa("mapdsk",work1,nadres,length,ierr)
      write(16,4214)
      call plotx(psig,work1,work2,0,1,1,3)
      write(string,4214)
      call gtext(string,-1,-1)
4214  format('  bigpsi')
c
      nadres=nadres+length
      length=nths*nsf
      call int2spl(gps,grpssq_cdf,0,+1,0)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",grpssq_cdf,nadres,length,ierr)
      write(16,4215)
      call plotx(psig,work1,grpssq_cdf,0,1,2,4)
      write(string,4215)
      call gtext(string,-1,-1)
      call framep
4215  format('  grpssq')
c
c...... plot bp/bt at various surfaces.
c
      mthh1 = mth/2 + 1
c
      do 5115 i = 1, mthh1
      xinf(i) = i
      zinf(i) = sqrt(abs(grpssq_cdf(i,2))) / (xzero*gg(2))
 5115 continue
      call posplt ( xinf,zinf, 1,mthh1, 1, 'bp/bt at 2')
c
      jsurf = (nosurf-1)/3
      do 5215 i = 1, mthh1
 5215 zinf(i) = sqrt(abs(grpssq_cdf(i,jsurf))) / (xzero*gg(jsurf))
      call posplt ( xinf,zinf, 1,mthh1, 2, 'bp/bt at 1/3' )
c
      jsurf = 2*(nosurf-1)/3
      do 5315 i = 1, mthh1
 5315 zinf(i) = sqrt(abs(grpssq_cdf(i,jsurf))) / (xzero*gg(jsurf))
      call posplt ( xinf,zinf, 1,mthh1, 3, 'bp/bt at 2/3' )
c
      jsurf = nosurf
      do 5415 i = 1, mthh1
 5415 zinf(i) = sqrt(abs(grpssq_cdf(i,jsurf))) / (xzero*gg(jsurf))
      call posplt ( xinf,zinf, 1,mthh1, 4, 'bp/bt at edge' )
      call framep
c
      nadres=nadres+length
      call int2spl(xsq,xmap,0,+1,0)
c     Replaced call zwr(inp3,xmap,length,nadres,lgivup,999)
      call putwa("mapdsk",xmap,nadres,length,ierr)
      write(16,4216)
      call plotx(psig,work1,xmap,0,1,2,1)
      write(string,4216)
      call gtext(string,-1,-1)
4216  format('  xsq')
c
      nadres=nadres+length
      do 506 j=1,npsi
      do 506 i=3,nthe2
         if(j.eq.1.and.kshift.eq.1)xjacob(i,1)=xjacob(i,2)
         w2(i,j)=(dpx(i,j)**2+dpz(i,j)**2)
     1        *xsq(i,j)/xjacob(i,j)**2
506   continue
      call sym(w2,+1)
      call int2spl(w2,grthsq_cdf,0,+1,0)
      call dpsinor2(grthsq_cdf)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",grthsq_cdf,nadres,length,ierr)
      write(16,4217)
      call plotx(psig,work1,grthsq_cdf,1,2,2,2)
      write(string,4217)
      call gtext(string,-1,-1)
4217  format('  grthsq note ln scale and+-')
c
      nadres=nadres+length
      do 507 j=1,npsi
      do 507 i=3,nthe2
         if(j.eq.1.and.kshift.eq.1)xjacob(i,1)=xjacob(i,2)
         w2(i,j)=-(dtx(i,j)*dpx(i,j)+dtz(i,j)*dpz(i,j))
     1        *xsq(i,j)/xjacob(i,j)**2
507   continue
      call sym(w2,-1)
      call int2spl(w2,grpsth_cdf,0,-1,0)
      call dpsinor(grpsth_cdf)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",grpsth_cdf,nadres,length,ierr)
      write(16,4218)
      call plotx(psig,work1,grpsth_cdf,0,2,2,3)
      write(string,4218)
      call gtext(string,-1,-1)
4218  format('  grpsth')
c
      nadres=nadres+length
      call dt(w2,w3,+1)
      call int2spl(w3,gptdth_cdf,0,+1,0)
      call dpsinor(gptdth_cdf)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",gptdth_cdf,nadres,length,ierr)
      write(16,4219)
      call plotx(psig,work1,gptdth_cdf,0,1,2,4)
      write(string,4219)
      call gtext(string,-1,-1)
      call framep
4219  format('  gptdth')
c
      nadres=nadres+length
      call int2spl(dpxsq,work3,0,+1,0)
      call dpsinor(work3)
c     Replaced call zwr(inp3,work3,length,nadres,lgivup,999)
      call putwa("mapdsk",work3,nadres,length,ierr)
      write(16,4220)
      call plotx(psig,work1,work3,1,1,2,1)
      write(string,4220)
      call gtext(string,-1,-1)
4220  format('  xsqdps note ln scale and +-')
c
      nadres=nadres+length
      call dt(gps,w2,-1)
      call int2spl(w2,gpsdth_cdf,0,-1,0)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",gpsdth_cdf,nadres,length,ierr)
      write(16,4221)
      call plotx(psig,work1,gpsdth_cdf,0,2,2,2)
      write(string,4221)
      call gtext(string,-1,-1)
4221  format('  gpsdth')
c
      nadres=nadres+length
      call dt(xsq,w2,-1)
      call int2spl(w2,work2,0,-1,0)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",work2,nadres,length,ierr)
      write(16,4222)
      call plotx(psig,work1,work2,0,2,2,3)
      write(string,4222)
      call gtext(string,-1,-1)
4222  format('  xsqdth')
c
      nadres=nadres+length
      call int2spl(xjacob,zmap,0,+1,0)
c     Replaced call zwr(inp3,zmap,length,nadres,lgivup,999)
      call putwa("mapdsk",zmap,nadres,length,ierr)
      write(16,4223)
      call plotx(psig,work1,zmap,0,1,2,4)
      write(string,4223)
      call gtext(string,-1,-1)
      call framep
4223  format('  xjacob')
c
c...... plot xjacob at various angles.
c
      mthh1 = mth/8
      jang=mthh1 + 1
c
      do 5135 i = 1, nosurf
      workz(i) = zmap(jang,i)
 5135 continue
      call posplt ( psig,workz, 1,nosurf, 1, 'xjacob at 1/8')
c
      jang = jang + mthh1
      do 5235 i = 1, nosurf
 5235 workz(i) = zmap(jang,i)
      call posplt ( psig,workz, 1,nosurf, 2, 'xjacob at 1/4' )
c
      jang = jang + mthh1
      do 5335 i = 1, nosurf
 5335 workz(i) = zmap(jang,i)
      call posplt ( psig,workz, 1,nosurf, 3, 'xjacob at 3/8' )
c
      jang = jang + mthh1
      do 5435 i = 1, nosurf
 5435 workz(i) = zmap(jang,i)
      call posplt ( psig,workz, 1,nosurf, 4, 'xjacob at 1/2' )
      call framep
c
c
      nadres=nadres+length
      call int2spl(dpxj,xjprime_cdf,0,+1,0)
      call dpsinor(xjprime_cdf)
ccc correction of xjprym near magnetic axis
c
      if(np.eq.1.and.mjac.eq.1.and.ljac.eq.0) go to 123
      write(16,112) (xjprime_cdf(17,k),k=1,7)
c     call exi2(work2,2,0)
      write(16,112) (xjprime_cdf(17,k),k=1,7)
  123 continue
      do 91 i=1,mth2
         yc=(xjprime_cdf(i,2)-xjprime_cdf(i,3))/(xjprime_cdf(i,2)
     &        +xjprime_cdf(i,3))
ckg      if(abs(2.*yc).gt.(1./nosurf)) go to 91
ckg two points interpolation
ckg      work2(i,1)=2.*work2(i,2)-work2(i,3)
         xjprime_cdf(i,1)=3.*xjprime_cdf(i,2)-3.*xjprime_cdf(i,3)
     &        +xjprime_cdf(i,4)
      write(16,92) i,(xjprime_cdf(i,j),j=1,6)
   91 continue
   92 format('xjprym for i=',i3,1x,6e13.5)
c
  112 format(8e12.5)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",xjprime_cdf,nadres,length,ierr)
c
c...... plot xjprym at various angles.
c
      mthh1 = mth/8
      jang=mthh1 + 1
c
      do 5125 i = 1, nosurf
      workz(i) = xjprime_cdf(jang,i)
 5125 continue
      call posplt ( psig,workz, 1,nosurf, 1, 'xjprym at 1/8')
c
      jang = jang + mthh1
      do 5225 i = 1, nosurf
 5225 workz(i) = xjprime_cdf(jang,i)
      call posplt ( psig,workz, 1,nosurf, 2, 'xjprym at 1/4' )
c
      jang = jang + mthh1
      do 5325 i = 1, nosurf
 5325 workz(i) = xjprime_cdf(jang,i)
      call posplt ( psig,workz, 1,nosurf, 3, 'xjprym at 3/8' )
c
      jang = jang + mthh1
      do 5425 i = 1, nosurf
 5425 workz(i) = xjprime_cdf(jang,i)
      call posplt ( psig,workz, 1,nosurf, 4, 'xjprym at 1/2' )
      call framep
c
      write(16,4224)
      call plotx(psig,work1,xjprime_cdf,1,1,2,1)
      write(string,4224)
      call gtext(string,-1,-1)
4224  format('  xjprym')
c
c    check qp by integrating (d/dpsi)(rg*xjacob/xsq) over theta
c
      do 969 j=1,nosurf
      workx(j)=0.0
      workz(j)=0.0
      do 968 i=1,mth
         work1(j)=xjprime_cdf(i,j)/zmap(i,j)-work3(i,j)/xmap(i,j)
         workx(j)=workx(j)+zmap(i,j)/xmap(i,j)
         workz(j)=workz(j)+work1(j)*zmap(i,j)/xmap(i,j)
  968 continue
      workz(j)=xzero*(workx(j)*ggp(j)+workz(j)*gg(j))
      workz(j)=workz(j)/mth
      work1(j)=xzero*workx(j)*gg(j)/mth
  969 continue
      call posplt ( psig,work1, 1,nosurf, 2, 'average q')
      call posplt ( psig,workz, 1,nosurf, 3, 'average qp')
c     write(string,4227)
c     call gtext(string,-1,-1)
4227  format('  qp')
c
      nadres=nadres+length
      call int2spl(delta,delta_cdf,0,-1,0)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",delta_cdf,nadres,length,ierr)
      write(16,4225)
      call plotx(psig,work1,delta_cdf,0,2,2,4)
      write(string,4225)
      call gtext(string,-1,-1)
4225  format('  delta')
      call framep
c
      nadres=nadres+length
      call int2spl(dpqdel,work2,0,-1,0)
      call dpsinor(work2)
c     Replaced call zwr(inp3,work2,length,nadres,lgivup,999)
      call putwa("mapdsk",work2,nadres,length,ierr)
      write(16,4226)
      call plotx(psig,work1,work2,1,2,2,1)
      write(string,4226)
      call gtext(string,-1,-1)
4226  format('  qdelp')
      
      do i=1,nosurf
         work1(i)=i
      enddo
      call posplt (work1, psig, 1,nosurf, 2, 'Psi(i)' )
c      call plotx(psig,work1,work2,1,2,2,1)
c      write(string,4227)
c      call gtext(string,-1,-1)
c 4227 format('  psig(i)')
      call framep
      nadres=nadres+length
      iss=0
c     Added
      call wclose("eqb1",ierr)
c     call close(6)
c     Replaced call zcl(inp2,iss,997)
      call wclose("mpout1",ierr)
c     Replaced call zcl(inp3,iss,998)
      call wclose("mapdsk",ierr)
ckg_cdf this is interface for the mapdsk.cdf file for i2mex
      eqfile='mapdsk.cdf'
      CALL cdfopn(cdfid,eqfile,'w')
c      do i=1,20
c         write(ntitle_cdf(i),'(A)') ntitle(i)
c      enddo
      write(Title,'(10A)') (ntitle(i),i=1,10)
      vartype = 'CHAR'
      dimlens(1) = 80
      dimlens(2) = 1
      dimlens(3) = 1
      CALL cdfDefVar(cdfid,'Title',dimlens,vartype)
      CALL cdfDefVar(cdfid,'date',dimlens,vartype)
      CALL cdfDefVar(cdfid,'comments',dimlens,vartype)
      dimlens(1) = 1
      vartype = 'INT'
      CALL cdfDefVar(cdfid,'mth',dimlens,vartype)
      CALL cdfDefVar(cdfid,'nosurf',dimlens,vartype)
      CALL cdfDefVar(cdfid,'remap',dimlens,vartype)
      CALL cdfDefVar(cdfid,'mx',dimlens,vartype)
      CALL cdfDefVar(cdfid,'npsi',dimlens,vartype)
      CALL cdfDefVar(cdfid,'kb',dimlens,vartype)
      vartype = 'R8'
      CALL cdfDefVar(cdfid,'Xmin',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Xmax',dimlens,vartype)
      CALL cdfDefVar(cdfid,'X0',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Xmag',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Zmin',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Zmax',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Z0',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Zmag',dimlens,vartype)
      CALL cdfDefVar(cdfid,'B0',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Ip',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Beta',dimlens,vartype)
      CALL cdfDefVar(cdfid,'BetaStar',dimlens,vartype)
      CALL cdfDefVar(cdfid,'BetaN',dimlens,vartype)
      CALL cdfDefVar(cdfid,'li',dimlens,vartype)
      CALL cdfDefVar(cdfid,'PPF',dimlens,vartype)
      CALL cdfDefVar(cdfid,'Upsiln',dimlens,vartype)
      dimlens(1) = nosurf
ckg_cdf some of the valiables below are in the mapdsk file but we do not
ckg_cdf put them there because they are not used in nova of i2mex 
ckg_cdf you can generate them by runing i2all(i2mex) with our mapdsk file
call cdfDefVar(iu, 'JdotBOverBSquare', dims, 'R8')
call cdfDefVar(iu, 'ballooning_alpha', dims, 'R8')
call cdfDefVar(iu, 'local_shear_s', dims, 'R8')
call cdfDefVar(iu, 'Di', dims, 'R8')
call cdfDefVar(iu, 'Dr', dims, 'R8')
call cdfDefVar(iu, 'surface_averaged_radius', dims, 'R8')
      call cdfDefVar(cdfid, 'PsiBig', dimlens, vartype)
call cdfDefVar(iu, 'V_of_psi', dims, 'R8')
call cdfDefVar(iu, 'Vp_of_psi', dims, 'R8')
      call cdfDefVar(cdfid, 'pa', dimlens, vartype)
      call cdfDefVar(cdfid, 'ppa',dimlens, vartype)
      call cdfDefVar(cdfid, 'qa', dimlens, vartype)
      call cdfDefVar(cdfid, 'qpa',dimlens, vartype)
call cdfDefVar(iu, 'qppa', dims, 'R8')
      call cdfDefVar(cdfid, 'ga', dimlens, vartype)
      call cdfDefVar(cdfid, 'gpa', dimlens, vartype)
call cdfDefVar(iu, 'spest1', dims, 'R8')
      dimlens(1) = mth1
      call cdfDefVar(cdfid, 'xinf', dimlens, vartype)
      call cdfDefVar(cdfid, 'zinf', dimlens, vartype)
      dimlens(2) = nosurf
      call cdfDefVar(cdfid, 'xa', dimlens, vartype)
      call cdfDefVar(cdfid, 'za', dimlens, vartype)
      call cdfDefVar(cdfid, 'xpth', dimlens, vartype)
      call cdfDefVar(cdfid, 'zpth', dimlens, vartype)
      call cdfDefVar(cdfid, 'xpsi', dimlens, vartype)
      call cdfDefVar(cdfid, 'zpsi', dimlens, vartype)
      call cdfDefVar(cdfid, 'grpssq', dimlens, vartype)
      call cdfDefVar(cdfid, 'grthsq', dimlens, vartype)
c
      call cdfDefVar(cdfid, 'grpsth', dimlens, vartype)
      call cdfDefVar(cdfid, 'gptdth', dimlens, vartype)
      call cdfDefVar(cdfid, 'gpsdth', dimlens, vartype)
      call cdfDefVar(cdfid, 'xjacob', dimlens, vartype)
      call cdfDefVar(cdfid, 'xjprym', dimlens, vartype)
      call cdfDefVar(cdfid, 'delta', dimlens, vartype)
      call cdfDefVar(cdfid, 'qdelp', dimlens, vartype)
call cdfDefVar(cdfid, 'qdelt', dimlens, vartype)

c      write(*,*) work_cdf(1:mth1,1)
c      write(*,*) work_cdf(1:mth1,2)
c write the output 
      CALL cdfPutVar(cdfid,'Title',Title)
      Title='today'
      CALL cdfPutVar(cdfid,'date',Title)
      Title='none'
      CALL cdfPutVar(cdfid,'comments',Title)
      mth_cdf=mth
      CALL cdfPutVar(cdfid,'mth',mth_cdf)
      mth_cdf=nosurf
      CALL cdfPutVar(cdfid,'nosurf',mth_cdf)
      mth_cdf=0
      CALL cdfPutVar(cdfid,'remap',mth_cdf)
      mth_cdf=nxy(2)  !=mjac
      CALL cdfPutVar(cdfid,'mx',mth_cdf)
      mth_cdf=nxy(1)  !=np
      CALL cdfPutVar(cdfid,'npsi',mth_cdf)
      mth_cdf=nxy(3)  !=ljac
      CALL cdfPutVar(cdfid,'kb',mth_cdf)
      CALL cdfPutVar(cdfid,'Xmin',xmin)
      CALL cdfPutVar(cdfid,'Xmax',xmax)
      xinf(2)=(work_cdf(1,1)+min(work_cdf(mth1/2,1),work_cdf(mth1/2+1,1)
     &     ))*0.5
      xinf(2)=xzero
      CALL cdfPutVar(cdfid,'X0',xinf(2))
      CALL cdfPutVar(cdfid,'Xmag',x(1,1))
      CALL cdfPutVar(cdfid,'Zmin',zmin)
      CALL cdfPutVar(cdfid,'Zmax',zmax)
      xinf(1)=(work_cdf(1,2)+work_cdf(mth1/2+1,2))*0.5
      CALL cdfPutVar(cdfid,'Z0',xinf(1))
      CALL cdfPutVar(cdfid,'Zmag',z(1,1))
      xinf(4)=1.
      CALL cdfPutVar(cdfid,'B0',xinf(4))
      xinf(3)=axx(11)*xinf(4)*xinf(2)
      CALL cdfPutVar(cdfid,'Ip',xinf(3))
      CALL cdfPutVar(cdfid,'Beta',axx(12))
ckg_debug
      CALL cdfPutVar(cdfid,'BetaStar',axx(12)*2) !is not correct
      xinf(1)=100.*axx(12)*(xmax-xmin)*0.5*0.4*pi/xinf(3)
      CALL cdfPutVar(cdfid,'BetaN',xinf(1))
      xinf(1)=2.*bp2sq*volume/xinf(3)**2/xinf(2)
      CALL cdfPutVar(cdfid,'li',xinf(1))
      xinf(1)=p(1)/(axx(12)*xinf(4)**2/2.)
      CALL cdfPutVar(cdfid,'PPF',xinf(1))
      CALL cdfPutVar(cdfid,'Upsiln',axx(13))
ckg_cdf write arrays
      call cdfPutVar(cdfid, 'PsiBig', psibig_cdf(1:nosurf))
      call cdfPutVar(cdfid, 'pa',  p_cdf(1:nosurf))
      call cdfPutVar(cdfid, 'ppa',pp_cdf(1:nosurf))
      call cdfPutVar(cdfid, 'qa',  q_cdf(1:nosurf))
      call cdfPutVar(cdfid, 'qpa', qp_cdf(1:nosurf))
      call cdfPutVar(cdfid, 'ga', gg(1:nosurf))
      call cdfPutVar(cdfid, 'gpa', ggp(1:nosurf))
      call cdfPutVar(cdfid, 'xinf',work_cdf(1:mth1,1))
      call cdfPutVar(cdfid, 'zinf',work_cdf(1:mth1,2))
      call cdfPutVar(cdfid, 'xa',x_cdf(1:mth1,1:nosurf)) 
      call cdfPutVar(cdfid, 'za',z_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'xpth',xdth_cdf(1:mth1,1:nosurf)) 
      call cdfPutVar(cdfid, 'zpth',zdth_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'xpsi',xpsi_cdf(1:mth1,1:nosurf)) 
      call cdfPutVar(cdfid, 'zpsi',zpsi_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'grpssq',grpssq_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'grthsq',grthsq_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'grpsth',grpsth_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'gptdth',gptdth_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'gpsdth',gpsdth_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'xjacob',zmap(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'xjprym',xjprime_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'delta',delta_cdf(1:mth1,1:nosurf))
      call cdfPutVar(cdfid, 'qdelp',work2(1:mth1,1:nosurf))
c
      write(*,*)'first values', psig(1),psibig_cdf(1)/(psibig_cdf(nosurf
     &     )-psibig_cdf(1)),psibig_cdf(nosurf)/2/pi,psival(npsi)
     &     -psival(1),psival(1)
      CALL cdfcls(cdfid)
ckg_cdf
c     Close plot file
      call plote
c
      stop
  997 call errmes(101,4hmain,997)
  998 call errmes(101,4hmain,998)
      end
c
c****************************************
c
ckg Looks like output of this program is the grid for the equilibrium 
ckg mapping. Two arrays are used on the output psig and fg: psig is the 
ckg square root of normalized poloidal flux, and fg is 2.*psitot*psig
      subroutine grid(igrid,no,n00,nnos,nos,rs,xm1,gdxm1,d1f1,xm2,gdxm2
     &     ,d1f2,xm3,gdxm3,d1f3)
      include 'clichpar.h'
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
ccc    igrid=1; uniform sqrt(psi) grid
ccc     igrid=2 uniform psi grid
ccc     igrid=3 tailored grid with dense mesh at q=1 surface
ckg     igrid=4 is nonuniform grid (see comment above).
      psitot=psival(npsi)-psival(1)
      if(igrid.le.2)then
ckg igrid = 1 or 2
         drg=1./(no-1.)
         do 11 jj=1,no
            rgrd=(jj-1.)*drg
            if(igrid.eq.2) rgrd=sqrt(rgrd)
cc  psig is the square root of normalized poloidal flux
            psig(jj)=rgrd
            fg(jj)=2.*psitot*psig(jj)
 11      continue
         psig(no)=1.
         fg(no)=2.*psitot
         psig(1)=.1*psig(2)
         fg(1)=2.*psitot*psig(1)
         return
      elseif(igrid.ge.4)then
         drg=1./(no-3.)
Ckg place three points near origin, NOVAKN may have a bug in the
C   boundary condition, hard to find. With this fix it converges fine.
C
         psig(1)=.1*drg
         psig(2)=psig(1)+2.*psig(1)
         psig(3)=psig(2)+4.*psig(1)
         drg=(1.-psig(3))*drg
         do jj=4,no
Ckg this changes from 0 to 1 between 3rd and no's points
            rgrd=(jj-3.)*drg
Ckg compute the derivative first
            d1fx=1.
            if(rgrd.gt.(xm1-gdxm1).and.rgrd.lt.(xm1+gdxm1))then
               hhh=max(0.1,d1f1)
               d1fx=min(d1fx,hhh+(1.-hhh)*((rgrd-xm1)/gdxm1)**2)
            endif
            if(rgrd.gt.(xm2-gdxm2).and.rgrd.lt.(xm2+gdxm2))then
               hhh=max(0.1,d1f2)
               d1fx=min(d1fx,hhh+(1.-hhh)*((rgrd-xm2)/gdxm2)**2)
            endif
            if(rgrd.gt.(xm3-gdxm3).and.rgrd.lt.(xm3+gdxm3))then
               hhh=max(0.1,d1f3)
               d1fx=min(d1fx,hhh+(1.-hhh)*((rgrd-xm3)/gdxm3)**2)
            endif
Ckg make the next grid point according to the derivative
            psig(jj)=psig(jj-1)+d1fx*drg
c            psig(jj)=psig(3)+(rgrd**(0.8)
         enddo
         hhh=(1.-psig(3))/(psig(no)-psig(3))
         do jj=1,no
            if(jj.gt.3)psig(jj)=(psig(jj)-psig(3))*hhh+psig(3)
            if(jj.eq.no)psig(jj)=1.
            fg(jj)=2.*psitot*psig(jj)
         enddo
         return
      endif
ccc   tailored grids
ccc  n00 is the # of subdivision in the first uniform grid spacing
ccc  nos is one half # of uniform grid near q=1 surface for dense grids
ccc  nnos is one half # of dense grid points near q=1 surface
ccc  no must be greater than (n00+2*nnos+2)
ccccccccccccccccccccccccccccc
c     n00=1
c     nos=2
c     nnos=3
c     rs=.6
      drs=1./nnos
      nno=no-2*nnos+2*nos-n00
      dr=1./(nno-1)
      als=nos*dr
ccc   get radial q=1 surface rs and ns1, corresp uniform surface #
c
      ns1=rs/dr+1
      nns=ns1-nos+nnos+n00
      nns1=nns-nnos-1
      nns2=nns+nnos+1
c
      psig(nns)=rs
      fg(nns)=2.*psitot*psig(nns)
      do 102 j=1,nnos
      rrs=(j*drs)
      rrs=rrs*als
      psig(nns-j)=rs-rrs
      psig(nns+j)=rs+rrs
      fg(nns-j)=2.*psitot*psig(nns-j)
      fg(nns+j)=2.*psitot*psig(nns+j)
  102 continue
c
      do 105 m=1,n00
      psig(m)=(m-1)*dr/n00
      fg(m)=2.*psitot*psig(m)
  105 continue
      psig(1)=.001
      fg(1)=2.*psitot*psig(1)
      n00p1=n00+1
      do 103 k=n00p1,nns1
      rrs=((k-n00)*dr)
      psig(k)=rrs
      fg(k)=2.*psitot*psig(k)
  103 continue
c
      do 104 l=nns2,no
      rrs=((l-2*nnos+2*nos-n00-1)*dr)
      psig(l)=rrs
      fg(l)=2.*psitot*psig(l)
  104 continue
      psig(no)=1.0
      fg(no)=2.*psitot
      return
      end
c****************************************
      subroutine dts(dum1,dum2,nsignx)
      include 'clichpar.h'
      dimension dum1(nq,ny),dum2(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      do 1 j=1,npsi
      do 1 i=3,nthe2
      dum2(i,j)=(dum1(i+1,j)-dum1(i-1,j))/
     1(2.*dth)
1     continue
      call sym(dum2,nsignx)
      return
      end
      subroutine dt(dum1,dum2,nsignx)
      include 'clichpar.h'
      dimension dum1(nq,ny),dum2(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var2/dtst(nq,ny),dpts(nq,ny)
      do 1 j=1,npsi
      do 1 i=3,nthe2
      dum2(i,j)=(dum1(i+1,j)-dum1(i-1,j))/
     1(2.*dth*dtst(i,j))
1     continue
      call sym(dum2,nsignx)
      return
      end
      subroutine dp(dum1,dum2,nsignx)
      include 'clichpar.h'
      dimension dum1(nq,ny),dum2(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var2/dtst(nq,ny),dpts(nq,ny)
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      call dps(dum1,dum2,nsignx)
      do 1 j=1,npsi
      do 1 i=3,nthe2
      dum2(i,j)=(dum1(i+1,j)-dum1(i-1,j))*
     1dpts(i,j)/(2.*dth)+dum2(i,j)
1     continue
      call sym(dum2,nsignx)
      return
      end
      subroutine dps(dum1,dum2,nsignx)
      include 'clichpar.h'
      dimension dum1(nq,ny),dum2(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      common/rcl/rpsi(ns)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      npsip2=npsi+2
      do 6 i=2,nthe3
      do 11 kk=1,npsi
      xtemp(kk)=dum1(i,kk)
   11 continue
c
      call intsplp(npsi,npsip2,xtemp,ytemp,rpsi,ztemp)
c
      do 4 jj=1,npsi
      dum2(i,jj)=ytemp(jj)
4     continue
c
6     continue
      call sym(dum2,nsignx)
      return
      end
      subroutine dpsinor2(work)
      include 'clichpar.h'
      dimension work(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      mth1=mth+1
      mth2=mth+2
      do 1 j=1,nosurf
       do 1 i=1,mth2
      work(i,j)=work(i,j)/fg(j)**2
    1 continue
       return
      end
      subroutine dpsinor(work)
      include 'clichpar.h'
      dimension work(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      mth1=mth+1
      mth2=mth+2
      do 1 j=1,nosurf
       do 1 i=1,mth2
      work(i,j)=work(i,j)/fg(j)
    1 continue
       return
      end
      subroutine intp1(px,pxa,nhf)
      include 'clichpar.h'
      dimension px(ny),pxa(ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      do 2 j=1,npsi
      ytemp(j)=(j-1.)/npsim
      if(nhf.eq.1)
     &ytemp(j)=(j-0.5)/npsim
2     continue
      j=4
      do 3 jj=1,nosurf
      pval=psig(jj)
4     continue
      if(j.eq.npsim) go to 5
      if(pval.lt.ytemp(j)) go to 5
      j=j+1
      if(j.ge.npsim) go to 5
      go to 4
5     continue
      call cubic
     1 (px(j-2),px(j-1),px(j),px(j+1)
     2,ytemp(j-2),ytemp(j-1),ytemp(j),ytemp(j+1)
     3,pval,pxa(jj),dum,0,ierr)
3     continue
      return
      end
      subroutine intp2(p2x,p2xa,nhf,nsignx,inv)
      include 'clichpar.h'
      dimension p2x(nq,ny),p2xb(nq,ns),p2xa(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      do 2 j=1,npsi
      ytemp(j)=(j-1.)/npsim
      if(nhf.eq.1)
     &ytemp(j)=(j-0.5)/npsim
2     continue
      imin=3
      imax=nthe2
      if(nsignx.eq.-1) imin=4
      if(nsignx.eq.-1) imax=nthe2-1
      do 1 i=imin,imax
      do 11 kk=2,npsi
      xtemp(kk)=p2x(i,kk)
      if(inv.eq.1) xtemp(kk)=1./xtemp(kk)
   11 continue
      j=4
      do 3 jj=1,nosurf
      pval=psig(jj)
4     continue
      if(j.eq.npsim) go to 5
      if(pval.lt.ytemp(j)) go to 5
      j=j+1
      if(j.ge.npsim) go to 5
      go to 4
5     continue
6     call cubic
     1(xtemp(j-2),xtemp(j-1),xtemp(j),xtemp(j+1)
     2,ytemp(j-2),ytemp(j-1),ytemp(j),ytemp(j+1)
     3,pval,xnew,dum,0,ierr)
      if(inv.eq.1) xnew=1./xnew
      p2xb(i,jj)=xnew
3     continue
1     continue
      do 8 j=1,nosurf
      p2xb(2,j)=p2xb(4,j)*nsignx
      p2xb(1,j)=p2xb(5,j)*nsignx
      p2xb(nthe3,j)=p2xb(nthe2-1,j)*nsignx
      p2xb(nthe4,j)=p2xb(nthe2-2,j)*nsignx
      if(nsignx.eq.-1) p2xb(3,j)=0.0
      if(nsignx.eq.-1) p2xb(nthe2,j)=0.0
8     continue
7     continue
      call intt(p2xb,p2xa,nsignx)
      return
      end
      subroutine int2spl(p2x,p2xa,nhf,nsignx,inv)
      include 'clichpar.h'
      dimension p2x(nq,ny),p2xb(nq,ns),p2xa(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      imin=3
      imax=nthe2
      if(nsignx.eq.-1) imin=4
      if(nsignx.eq.-1) imax=nthe2-1
      npsip2=npsi+2
      do 1 i=imin,imax
      do 11 kk=1,npsi
      xtemp(kk)=p2x(i,kk)
      if(inv.eq.1) xtemp(kk)=1./xtemp(kk)
   11 continue
c
      call intspl(nosurf,npsip2,xtemp,ytemp,psig,ztemp)
c
      if(inv.ne.1) go to 2
      do 3 jj=1,nosurf
      p2xb(i,jj)=1./ytemp(jj)
3     continue
      go to 1
    2 continue
      do 4 jj=1,nosurf
      p2xb(i,jj)=ytemp(jj)
4     continue
1     continue
c
      do 8 j=1,nosurf
      p2xb(2,j)=p2xb(4,j)*nsignx
      p2xb(1,j)=p2xb(5,j)*nsignx
      p2xb(nthe3,j)=p2xb(nthe2-1,j)*nsignx
      p2xb(nthe4,j)=p2xb(nthe2-2,j)*nsignx
      if(nsignx.eq.-1) p2xb(3,j)=0.0
      if(nsignx.eq.-1) p2xb(nthe2,j)=0.0
8     continue
c
      call intt(p2xb,p2xa,nsignx)
      return
      end
      subroutine intt(p2xb,p2xa,nsignx)
      include 'clichpar.h'
      dimension p2xb(nq,ns) ,p2xa(nt,ns)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var4/tb(nq,ns)
      mthx=mth/2+1
      rmth=mth
      mth1=mth+1
      mth2=mth+2
      do 1 j=1,nosurf
      do 2 i=2,nthe3
      ytemp(i)=tb(i,j)
2     continue
      i=4
      do 3 ii=1,mthx
      tval=(ii-1.)*2*pi/rmth
4     continue
      if(i.eq.nthe2) go to 5
      if(tval.lt.ytemp(i)) go to 5
      i=i+1
      if(i.eq.nthe2) go to 5
      go to 4
5     continue
      call cubic
     1(p2xb(i-2,j),p2xb(i-1,j),p2xb(i,j),p2xb(i+1,j)
     2,ytemp(i-2),ytemp(i-1),ytemp(i),ytemp(i+1)
     3,tval,p2xa(ii,j),dub,0,ierr)
3     continue
1     continue
      mthxm=mthx-1
      do 7 j=1,nosurf
      do 6 i=2,mthxm
      ii=mth-i+2
      p2xa(ii,j)=p2xa(i,j)*nsignx
6     continue
      if(nsignx.eq.-1) p2xa(1,j)=0.0
      if(nsignx.eq.-1) p2xa(mthx,j)=0.0
      p2xa(mth1,j)=p2xa(1,j)
      p2xa(mth2,j)=p2xa(2,j)
7     continue
      return
      end
      subroutine extap(x1,x2,x3,x4)
      include 'clichpar.h'
      x4=3.*x3-3.*x2+x1
      ddx1=x3-x2
      ddx2=x2-x1
      ddx=x4-x3
      pm=ddx*ddx1
      if(pm.gt.0.) return
      if(ddx2.eq.0.) x4=2.*x3-x2
      if(ddx2.eq.0.) return
      ddx=(ddx1*ddx1)/ddx2
      x4=x3+ddx
      return
      end
      subroutine ex(d,n)
      include 'clichpar.h'
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      dimension d(nq,ny)
      do 1 i=2,nthe3
      if(n.eq.1)
     1call extap(d(i,4),d(i,3),d(i,2),d(i,1))
      if(n.eq.2)
     1call extap(d(i,5),d(i,4),d(i,3),d(i,2))
      if(n.eq.npsi)
     1call extap(d(i,npsi-3),d(i,npsi-2),d(i,npsim),d(i,npsi))
1     continue
      return
      end
      subroutine sym(d,nsignx)
      include 'clichpar.h'
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      dimension d(nq,ny)
      do 1 j=1,npsi
      d(2,j)=d(4,j)*nsignx
      d(nthe3,j)=d(nthe2-1,j)*nsignx
      if(nsignx.ne.-1) go to 1
      d(3,j)=0.0
      d(nthe2,j)=0.0
1     continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     * * * number 3j * * *                                c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine cubic(x1,x2,x3,x4,y1,y2,y3,y4,y,ans,dans,is,ierr)
      include 'clichpar.h'
cray  lcm(cubic)
c      lagrange 4-points interpolation formula
c      this routine returns interpolated value and first derivative
c
      d1 = (y1-y2)*(y1-y3)*(y1-y4)
      d2 = (y2-y1)*(y2-y3)*(y2-y4)
      d3 = (y3-y1)*(y3-y2)*(y3-y4)
      d4 = (y4-y1)*(y4-y2)*(y4-y3)
      ans = x1*(y-y2)*(y-y3)*(y-y4)/d1
     1    + x2*(y-y1)*(y-y3)*(y-y4)/d2
     2    + x3*(y-y1)*(y-y2)*(y-y4)/d3
     3    + x4*(y-y1)*(y-y2)*(y-y3)/d4
c
      if ( is .eq. 0 )  return
      if(is.eq.2) go to 5
c
      dans = x1/d1*((y-y3)*(y-y4)+(y-y2)*(y-y4)+(y-y2)*(y-y3))
     1     + x2/d2*((y-y3)*(y-y4)+(y-y1)*(y-y4)+(y-y1)*(y-y3))
     2     + x3/d3*((y-y2)*(y-y4)+(y-y1)*(y-y4)+(y-y1)*(y-y2))
     3     + x4/d4*((y-y2)*(y-y3)+(y-y1)*(y-y3)+(y-y1)*(y-y2))
      return
    5 continue
      dans = x1/d1*((y-y3)+(y-y4)+(y-y2)+(y-y4)+(y-y2)+(y-y3))
     1     + x2/d2*((y-y3)+(y-y4)+(y-y1)+(y-y4)+(y-y1)+(y-y3))
     2     + x3/d3*((y-y2)+(y-y4)+(y-y1)+(y-y4)+(y-y1)+(y-y2))
     3     + x4/d4*((y-y2)+(y-y3)+(y-y1)+(y-y3)+(y-y1)+(y-y2))
      return
   11 continue
      if(y.gt.y3) go to 10
c
      d1 = (x3-x1)/(y3-y1)
      d2 = (x4-x2)/(y4-y2)
      fac2 = (y-y2)/(y3-y2)
      fac1 = 1.-fac2
      dans = fac1*d1 + fac2*d2
c
      return
   10 continue
      anum = (y3-y2)*(x4-x3) + (y4-y3)*(x2-x3)
      adem = (y4-y2)*(y4-y3)**2
      a = anum/adem
      b = (x4-x2)/(y4-y2) - 2.*y3*a
      dans = 2.*a*y + b
      return
   20 continue
      anum = (y2-y3)*(x1-x2)+(y1-y2)*(x3-x2)
      adem = (y1-y3)*(y1-y2)**2
      a = anum/adem
      b = (x3-x1)/(y3-y1) - 2.*y2*a
      dans = 2.*a*y + b
      return
      end
      subroutine plotx(xtemp,work1,work2,ip4,ip3,ip2,ip)
      include 'clichpar.h'
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      dimension work1(ns),work2(nt,ns),xtemp(ns)
ccc   return
      iii=ip3
      if(ip3.gt.1) iii=(mth+4.)/4.+.001
      if(ip2.eq.1) go to 3
      do 1 j=1,nosurf
      work1(j)=work2(iii,j)
1     continue
3     continue
      ymax2=-1.e32
      ymin2=1.e32
      do 951 j=1,nosurf
      if(work1(j).gt.ymax2) ymax2 = work1(j)
      if(work1(j).lt.ymin2) ymin2 = work1(j)
951   continue
      xmin2 = xtemp(1)
      xmax2 = xtemp(nosurf)
c     if(ymax2.ne.ymin2)go to 310
      ay1=abs(ymax2-ymin2)
      if(ay1.gt.1.e-10) go to 310
      ymax2=2.*ymax2
      ymin2=0.0
      if(ymax2.eq.0.0) ymax2=1.0
      if(ymax2.lt.0.0) ymin2=2.*ymax2
      if(ymax2.lt.0.0) ymax2=0.0
310   continue
      iiix=ip4
      if(ymax2*ymin2.lt.0.) iiix=0
      if(iiix.eq.1) go to 40
      if(ip.eq.1)
     1call mapg(xmin2,xmax2,ymin2,ymax2,.12,.42,.65,.95)
      if(ip.eq.2)
     1call mapg(xmin2,xmax2,ymin2,ymax2,.60,.90,.65,.95)
      if(ip.eq.3)
     1call mapg(xmin2,xmax2,ymin2,ymax2,.12,.42,.15,.45)
      if(ip.eq.4)
     1call mapg(xmin2,xmax2,ymin2,ymax2,.60,.90,.15,.45)
      call trace(xtemp,work1,nosurf,-1,-1,0,0)
      write(16,1000)(work1(j),j=1,nosurf)
      call setlch(xmin2+(xmax2-xmin2)/3.,ymin2-(ymax2-ymin2)/20.,
     $           0,2,0,-1)
      go to 41
40    continue
      if(ymax2.gt.0) go to 45
      ymin3=-ymax2
      ymax3=-ymin2
      ymin2=ymin3
      ymax2=ymax3
      do 46 j=1,nosurf
      work1(j)=-work1(j)
46    continue
45    continue
      if (ymin2.eq.ymax2) ymin2=ymax2-1.
      if(ip.eq.1)
     1call mapgsl(xmin2,xmax2,ymin2,ymax2,.12,.42,.65,.95)
      if(ip.eq.2)
     1call mapgsl(xmin2,xmax2,ymin2,ymax2,.60,.90,.65,.95)
      if(ip.eq.3)
     1call mapgsl(xmin2,xmax2,ymin2,ymax2,.12,.42,.15,.45)
      if(ip.eq.4)
     1call mapgsl(xmin2,xmax2,ymin2,ymax2,.60,.90,.15,.45)
      call trace(xtemp,work1,nosurf,-1,-1,0,0)
      write(16,1000)(work1(j),j=1,nosurf)
1000  format(1x,10e12.4)
      if(ip.eq.1) rix=11
      if(ip.eq.1) riy=30
      if(ip.eq.2) rix=53
      if(ip.eq.2) riy=30
      if(ip.eq.3) rix=11
      if(ip.eq.3) riy=10
      if(ip.eq.4) rix=53
      if(ip.eq.4) riy=10
      call setch(rix,riy,0,1,0,-1)
41    continue
      return
      end
      subroutine intp2s(p2x,p2xb,nhf)
      include 'clichpar.h'
      dimension p2x(nq,ny),p2xb(nq,ns),p2xa(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      common/work/ytemp(ns),xtemp(ns),ztemp(ns)
      do 2 j=1,npsi
      ytemp(j)=(j-1.)/npsim
      if(nhf.eq.1)
     &ytemp(j)=(j-0.5)/npsim
2     continue
      do 1 i=3,nthe2
      j=4
      do 3 jj=1,nosurf
      pval=psig(jj)
4     continue
      if(j.eq.npsim) go to 5
      if(pval.lt.ytemp(j)) go to 5
      j=j+1
      if(j.ge.npsim) go to 5
      go to 4
5     continue
6     call cubic
     1(p2x(i,j-2),p2x(i,j-1),p2x(i,j),p2x(i,j+1)
     2,ytemp(j-2),ytemp(j-1),ytemp(j),ytemp(j+1)
     3,pval,p2xb(i,jj),dum,0,ierr)
3     continue
1     continue
      return
      end
      subroutine exi1(w1,jj,nst)
      include 'clichpar.h'
      dimension w1(ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      do 1 jj=1,nst
      j=nst-jj
      j4=j+4
      j3=j+3
      j2=j+2
      j1=j+1
      call extap(w1(j4),w1(j3),w1(j2),w1(j1))
1     continue
      return
      end
      subroutine exi2(ww,ip,inv)
      include 'clichpar.h'
c......................................................................
c      interpolates array arrf from psia to psinew. uses a polynomial
c      fit with leading power in sqrt(psi)=e0 and npts points/powers
c      near the axis. this determines an jndex beyond which a cubic
c      spline interpolation takes over.
c      march 1881 jm.
c.......................................................................
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),f(ny),psival(ny),fp(ny)
      dimension iop(2),coef(10),arrf(ns),psia(ns),temp3(ns)
      dimension ww(nt,ns)
      data iop/5,5/
      npts=nptsi
      ipn=ip
      if(inv.eq.1.and.ip.eq.1) ipn=5
      if(inv.eq.1.and.ip.eq.2) ipn=4
      e0=ipn-3.
      psitot=psival(npsi)-psival(1)
      do 3 jj=1,nosurf
      psia(jj)=psig(jj)**2*psitot
    3 continue
      imin=1
      imax=mth+2
      do 90 i=imin,imax
      do 120 jj=1,nosurf
      arrf(jj)=ww(i,jj)
      if(inv.eq.1) arrf(jj)=1./arrf(jj)
120   continue
c      set up jndex
c
      jndex =jndexi
   10 jndex = jndex + 1
      if(jndex+npts .lt. nosurf) go to 20
c      cannot fit to npts polynomial. increase npts and start over
      npts = npts + 1
      jndex =jndexi
c      check on npts . if equal to 10 then error.
      if(npts .gt. 9) go to 100
      go to 10
   20 call fitter(npts,e0,psia(jndex),arrf(jndex),coef)
c
c      check the next point to verify goodness of polynomial
c
      jtest = jndex + npts
      arrnxt = 0.0
      do 30 j=1,npts
      eval = (e0 + j -1.) / 2.
   30 arrnxt = arrnxt + coef(j) * psia(jtest)**eval
c
c      now compare with known value
c
      if(abs(arrf(jtest)) .lt. 1.0e-10) go to 40
      tstval = (arrnxt-arrf(jtest))/arrf(jtest)
      if(abs(tstval).gt.1.e-2) write(16,1012) i,tstval
      if(abs(tstval).gt.1.e-2) write(16,1013) (arrf(k),k=1,jtest)
 1013 format(10e12.4)
1012  format(1x,'itheta=',i4,1x,'eii=',1e12.4)
ccc40 jndex = jtest
40    continue
cc    note manny takes jndex=jtest at this point
c      now interpolate to values of psinew. for psinew less than
c      psia(jndex) use polynomial approximation. beyond that use spline.
c     note: jndex=jndexi+1
c
      do 70 jj = 1, jndex-1
      psit = psia(jj)
      temp3(jj) = 0.0
      do 50 j = 1, npts
      eval = (e0 + j -1.) / 2.
   50 temp3(jj) = temp3(jj) + coef(j) * psit**eval
   70 continue
c
      do 80 jj=1, jndex-1
      if(inv.eq.1) temp3(jj)=1./temp3(jj)
   80 ww(i,jj) = temp3(jj)
90    continue
c
      return
  100 call errmes(101,6hinterp , 111)
      end
      subroutine fitter(npts,e0,psif,arrf,coef)
      include 'clichpar.h'
c......................................................................
c      fits npts points to a polynomial in sqrt(psi). the leading term
c      of the polynomial is (psi)**(e0/2). the coefficients coef are
c      determined by matching the values of the array arrf at the points
c      psif.      jm. march 1881
c......................................................................
      dimension a(10,10),b(10),ip(10),psif(1),arrf(1),
     .          coef(1), as(10,10)
c      load the work matrix a
      do 20 i=1, npts
c      set the exponent
      pexp = (e0 + i -1.) / 2.
      do 10 j=1, npts
      a(i,j) = 1.0
      if(pexp .eq. 0.) go to 10
      a(i,j) = (psif(j))**pexp
   10 continue
   20 continue
c
c      now decompose a
c
      call decomp(npts,10,a,ip)
c      check for singularity
c
      det = 0.0
      do 30 i=1, npts
   30 det = det + alog10(abs(a(i,i)))
      idet = ifix(det)
      det  = 10.**(det-idet)
      s    = 1.0
      do 40 i=1, npts
   40 s    = s * a(i,i)/abs(a(i,i))
      s    = s * ip(npts)
      if(s .gt. 0.) go to 50
      write(101,8000)s,det,idet
c     call errmes(101,6hfitter)
      return
   50 continue
c
c      compute coefficients
c
      do 70 i=1, npts
      do 60 j=1, npts
      b(j) = 0.0
      if(i .eq. j)b(j) = 1.0
   60 continue
c
      call solve(npts,10,a,b,ip)
c
      do 65 j=1, npts
   65 as(j,i) = b(j)
   70 continue
c
      do 80 j=1, npts
      coef(j) = 0.0
      do 80 i=1, npts
   80 coef(j) = coef(j) + as(i,j)*arrf(i)
c
      return
 8000 format(' determinant = ',f4.1,f5.2,i5,' in invert le zero',/)
      end
c......................................................................
      subroutine errmes ( nout,mesage,is )
      include 'clichpar.h'
c......................................................................
      dimension mesage(2)
      write(16,100)mesage,is
      write(nout,100)mesage,is
  100 format(' error in subroutine',1x,2a5/,1x,' error # ',i3 )
c      call exit()
      end
      subroutine solve( n, ndim, a, b, ip )
      include 'clichpar.h'
c     lcm(solve)
      dimension a(ndim,ndim), b(ndim)
      integer ip(ndim)
      if( n .eq. 1 ) go to 9
      nm1 = n - 1
      do 7   k = 1, nm1
      kp1 = k + 1
      m = ip(k)
      t = b(m)
      b(m) = b(k)
      b(k) = t
      do 7   i = kp1, n
    7 b(i) = b(i) + a(i,k) * t
      do 8   kb = 1, nm1
      km1 = n - kb
      k = km1 + 1
      b(k) = b(k) / a(k,k)
      t = -b(k)
      do 8   i = 1, km1
    8 b(i) = b(i) + a(i,k) * t
    9 b(1) = b(1) / a(1,1)
      return
      end
      subroutine decomp( n, ndim, a, ip )
      include 'clichpar.h'
c     lcm(decomp)
c     communications of the acm, volume 15, #4, april, 1972, page 274
c     algorithm 423 by cleve b. moler
      dimension a(ndim,ndim)
      integer ip(ndim)
c
      ip(n) = 1
      do 6   k = 1, n
      if( k .eq. n ) go to 5
      kp1 = k + 1
      m = k
      do 1   i = kp1, n
      if( abs(a(i,k)) .gt. abs(a(m,k))) m = i
    1 continue
      ip(k) = m
      if(m . ne. k) ip(n) = - ip(n)
      t = a(m,k)
      a(m,k) = a(k,k)
      a(k,k) = t
      if( t .eq. 0.0 ) go to 5
      do 2   i = kp1, n
    2 a(i,k) = -a(i,k)/t
      do 4   j = kp1, n
      t = a(m,j)
      a(m,j) = a(k,j)
      a(k,j) = t
      if( t .eq. 0.0 ) go to 4
      do 3   i = kp1, n
    3 a(i,j) = a(i,j) + a(i,k) * t
    4 continue
    5 if( a(k,k) .eq. 0.0 ) ip(n) = 0
    6 continue
      return
      end
      subroutine exi2s(w2,nst)
      include 'clichpar.h'
      dimension w2(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      mth2=mth+2
      do 2 i=1,mth2
      do 1 jj=1,nst
      j=nst-jj
      j4=j+4
      j3=j+3
      j2=j+2
      j1=j+1
      call extap(w2(i,j4),w2(i,j3),w2(i,j2),w2(i,j1))
1     continue
2     continue
      return
      end
      subroutine exo1(w1,nst)
      include 'clichpar.h'
      dimension w1(ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      do 1 j=nst,nosurf
      j4=j-3
      j3=j-2
      j2=j-1
      j1=j
      call extap(w1(j4),w1(j3),w1(j2),w1(j1))
1     continue
      return
      end
      subroutine exo2(w2,nst)
      include 'clichpar.h'
      dimension w2(nt,ns)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      mth2=mth+2
      do 2 i=1,mth2
      do 1 j=nst,nosurf
      j4=j-3
      j3=j-2
      j2=j-1
      j1=j
      call extap(w2(i,j4),w2(i,j3),w2(i,j2),w2(i,j1))
1     continue
2     continue
      return
      end
c
      subroutine posplt ( x,y, jmin,jmax, ipos, wlabel )
c
      dimension x(1), y(1), wlabel(2)
c     temp var for write 100-gtext
      character*100 string
c
      ymax = -1.0e32
      ymin =  1.0e32
c
      jd = jmax - jmin + 1
c
      do 50 j = jmin, jmax
c
      if ( y(j) .gt. ymax ) ymax = y(j)
      if ( y(j) .lt. ymin ) ymin = y(j)
c
   50 continue
c
      xmin = x(jmin)
      xmax = x(jmax)
      if ( ymax .eq. ymin ) ymax = ymin + 1.0
      dx = xmax - xmin
      dy = ymax - ymin
      xw = xmin + dx/3.0
      yw = ymax + dy/20.0
c
      go to (110,120,130,140), ipos
c
  110 call mapg(xmin,xmax, ymin,ymax, 0.12,0.42, 0.65,0.95)
      go to 200
  120 call mapg(xmin,xmax, ymin,ymax, 0.60,0.90, 0.65,0.95)
      go to 200
  130 call mapg(xmin,xmax, ymin,ymax, 0.12,0.42, 0.15,0.45)
      go to 200
  140 call mapg(xmin,xmax, ymin,ymax, 0.60,0.90, 0.15,0.45)
c
  200 continue
c
      call trace ( x, y, jd,-1,-1,0,0 )
      call setlch ( xw,yw,0,1,0,-1 )
      write (16, 300 ) wlabel
      write(16,1000)(y(j),j=jmin,jmax)
1000  format(1x,10e12.4)
c
      write ( string, 300 ) wlabel
      call gtext(string,-1,-1)
  300 format ( 2a8 )
c
      return
      end
c
c
      subroutine plotxz ( x,z, ndt,ndp, nt,np, ndth,ndps )
c
      dimension x(ndt,ndp), z(ndt,ndp)
c
      xmax = - 1.0e10
      zmax = - 1.0e10
      xmin =   1.0e10
      zmin =   1.0e10
c
      do 20 i = 1, nt
      if ( x(i,np) .gt. xmax ) xmax = x(i,np)
      if ( x(i,np) .lt. xmin ) xmin = x(i,np)
      if ( z(i,np) .gt. zmax ) zmax = z(i,np)
      if ( z(i,np) .lt. zmin ) zmin = z(i,np)
   20 continue
c
      zmin = - zmax
c
      dx = xmax - xmin
      dz = zmax - zmin
      dxz = amax1(dx,dz)
c
      x2 = (xmax+xmin)/2.0 + 0.55*dxz
      x1 = (xmax+xmin)/2.0 - 0.55*dxz
      z2 = (zmax+zmin)/2.0 + 0.55*dxz
      z1 = (zmax+zmin)/2.0 - 0.55*dxz
c
      call maps ( x1,x2,z1,z2, .175, .825, .35, 1. )
c
      npsid = np + ndps - 1
c
      do 100 j = 2, npsid,ndps
c
      jj = j
      if ( j .gt. np ) jj = np
c
      i = 1
      xx = x(i,jj)
      zz = z(i,jj)
      call setcrt(xx,zz)
c
      do 50 i = 2, nt
c
      xx = x(i,jj)
      zz = z(i,jj)
      call vector(xx,zz)
c
   50 continue
  100 continue
c
      nthid = nt + ndth - 1
c
      do 200 i = 1, nthid, ndth
c
      ii = i
      if ( i .gt. nt ) ii = nt
c
      j = 1
c
      xx = x(ii,j)
      zz = z(ii,j)
      call setcrt(xx,zz)
c
      do 150 j = 2, np
c
      xx = x(ii,j)
      zz = z(ii,j)
      call vector(xx,zz)
c
  150 continue
  200 continue
c
      return
      end
c*******************************************************************
      subroutine framep
c
      character*8 hedr(2),logo
      common / head / hedr
      dimension logo(2)
c
      data logo/8heqsjdmmc,5hblank/
      data icall /0/
c
c      placate compiler
c
      data ntim,ndate,mach,nsfx/0,0,0,0/
c
      call map(0.,1.,0.,1.,0.88,1.,0.85,1.)
cc    now plot the box
      call setcrt(0.,0.)
      call vector(0.,1.)
      call vector(1.,1.)
      call vector(1.,0.)
      call vector(0.,0.)
c     now enter the logo
      call setlch(0.1,0.75,0,1,0,-1)
c      call timedate(ntim,ndate,mach,nsfx)
c      call date(ndate)
c      logo(1) = ntim
c      logo(2) = mach
c      logo(3) = ndate
c      logo(4) = nsfx
c       replaced crtbcd ( logo, 2 )
      call gtext(logo,2*8,-1)
c      call setlch(0.1,0.55,0,1,0,-1)
c       replaced crtbcd(ndate)
c      call gtext(ndate,-1,-1)
      call setlch(0.1,0.35,0,1,0,-1)
c       replaced crtbcd(7hp p p l )
      call gtext("p p p l",-1,-1)
      call setlch(0.1,0.15,0,1,-1,-1)
c       replaced crtbcd(hedr(2),1)
      call gtext(hedr(2),1*8,-1)
      call dders(1)
      call frame(0)
c
      return
      end
c*****************************
c
      subroutine framepp
c
      character*8 hedr(2), logo(2)
      common / head / hedr
c
      data logo/8hmpsjdmmc,5hblank/
      data icall /0/
c
c      placate compiler
c
      data ntim,ndate,mach,nsfx/0,0,0,0/
c
c
c     go to 100
c
      call map(0.,1.,0.,1.,0.88,1.,0.85,1.)
cc    now plot the box
      call setcrt(0.,0.)
      call vector(0.,1.)
      call vector(1.,1.)
      call vector(1.,0.)
      call vector(0.,0.)
c     now enter the logo
      call setlch(0.1,0.75,0,1,0,-1)
c     Removed call timedate(ntim,ndate,mach,nsfx)
c      call date(ndate)
c      logo(1) = ntim
c      logo(2) = mach
c     Replaced call crtbcd ( logo, 2 )
c      call gtext(logo,2*8,-1)
      call setlch(0.1,0.55,0,1,-1,-1)
c     Replaced call crtbcd(ndate)
      call gtext(ndate,-1,-1)
      call setlch(0.1,0.35,0,1,-1,-1)
c     Replaced call crtbcd(7hp p p l )
      call gtext("p p p l",-1,-1)
      call setlch(0.1,0.15,0,1,-1,-1)
c     Replaced call crtbcd(hedr(2),1)
      call gtext(hedr(2),1*8,-1)
      call dders(1)
c
  100 continue
c
      call frame(0)
c
      return
      end
c
c
cccccc****************************************
ccccccc****************************************
c
      subroutine anfun(llanal)
ccc
ccc   solovev equilibrium
ccc
      include 'clichpar.h'
      integer*8 ntitle,nxx,nxy
      common ntitle(20),dat,nxx(5),axx(13),nxy(10),axy(10)
      common/var/ p(ny),pp(ny),q(ny),qp(ny),
     1 g(ny),gp(ny),fb(ny),fbp(ny),f(ny),fp(ny),
     1 psival(ny),x(nq,ny),z(nq,ny),aj3(nq,ny),aj(nq,ny)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/var3/fg(nsgrd),psig(nsgrd),fd(ny),psivad(ny),fpd(ny)
      common/panal/elip, epsl, q00
      namelist/dfun/elip,aspect,q00,nthe,npsi,lanal
      elip=2.
      aspect=3.
      q00=1.0
      npsi=51
      nthe=65
      lanal=0
      write(*,111)
  111 format(1x,'elip=2, aspect=3, q00=1., npsi=51 , nthe=65'/,
     &  1x,'lanal=0; choose lanal=1  for SOLOVEV EQUIL',//)
c
      write(16,dfun)
      write(*,*) "DFUN NAMELIST"
      read(*,dfun)
      llanal=lanal
c
      if(lanal.ne.1) return
c
      npsip2=npsi+2
      call grid(1,npsi,1,2,2,.0,xm1,gdxm1,d1f1,xm2,gdxm2
     &     ,d1f2,xm3,gdxm3,d1f3)
      call splprp(npsi)
      call splprp1(npsi,npsip2,psig)
c
      epsl=1./aspect
      psib=(epsl)**2*elip/(2.*q00)
      npsim=npsi-1
      npsip=npsi+1
      drg=1./npsim
      p0=(1.+elip**2)/(elip*q00)
      pi=acos(-1.)
      nthe1=nthe+1
      nthe2=nthe+2
      nthe3=nthe+3
      nthe4=nthe+4
c
cc   define variables in 'ntitle' common block
      xzero=1.0
      dth=pi/(nthe-1.)
      dr=sqrt(2.*pi)/npsim
      nxx(1)=nthe
      nxx(2)=npsi
      axx(1)=1.0
      axx(2)=1.0
      axx(3)=xzero
      axx(4)=xzero
      axx(5)=xzero
      axx(8)=0.0
      axx(9)=psib*2.*pi
      nxy(1)=1
      nxy(2)=1
      axy(1)=dth
      axy(2)=dr
      axy(3)=pi
c
cc  calculate equilibrium on a uniform rt(psi)-grid
cc  equilibrium input 'eqb1' is also calculated on a uniform rt(psi)-grid
cc
      do 10 is=1,npsip
      ism=is-1
      rgrd=ism*drg
      grd=rgrd
      if(is.eq.1) grd=0.01*drg
      psig(is)=grd
      fg(is)=2.*grd*psib
      psi=rgrd**2*psib
      psival(is)=psi
      psivad(is)=psi
      if(is.eq.1) psi=1.e-6*grd**2*psib
      call surf(psi,q(is),x(1,is),z(1,is),aj(1,is),aj3(1,is))
cccc...  calculate p, pp, g, gp, etc.
ccc...  p, g, and f are defined on the off-center grids
      p(is)=p0*(psib-psi)
      g(is)=1.0
      pp(is)=-p0
      gp(is)=0.0
      if(ism.gt.0) f(ism)=(psival(ism+1)-psival(ism))/dr
   10 continue
cccc  calculate at the magnetic axis
      q(1)=q00
      do 12 it=1,nthe4
      x(it,1)=xzero
      z(it,1)=0.0
   12 continue
cccc
ccc  calculate qp numerically
      call intsplp(npsi,npsip2,q,qp,psig,fpd)
ccc
ccc  calculate fb, fbp, etc.
ccc
      do 30 is=1,npsi
      qp(is)=qp(is)/fg(is)
      fb(is)=xzero*g(is)/q(is)
      fbp(is)=xzero*(gp(is)*q(is)-g(is)*qp(is))/q(is)**2
   30 continue
      return
      end
c
      subroutine surf(psi,q,x,z,aj,aj3)
      common/var1/ npsi,nthe2,dth,dr,nptsi,jndexi
     1,nosurf,npsim,npsip,nthe,nthe1,nthe3,nthe4,mth,pi
      common/panal/elip, epsl, q00
      dimension xb(501),zb(501),as(501),asb(501),bet(501)
      dimension x(1),z(1),aj(1),aj3(1)
      mt=500
      pi2=pi*2.0
      dbth=pi2/mt
      xr=sqrt(8.*q00*psi/elip)
      zr=sqrt(2.*elip*q00*psi)
      do 10 it=1,mt
      bet(it)=(it-1)*dbth
      sint=sin(bet(it))
      cost=cos(bet(it))
      xb(it)=sqrt(1.+xr*cost)
      zb(it)=zr*sint/xb(it)
      xb1=-xr*sint/(2.*xb(it))
      zb1=zr*cost/xb(it)
      zb1=zb1-xb1*zb(it)/xb(it)
      asb(it)=sqrt(xb1**2+zb1**2)
   10 continue
ccc  calculate arc length
      st=0.0
      as(1)=0.0
      do 20 it=2,mt
      st=st+(asb(it)+asb(it-1))*dbth/2.
      as(it)=st
   20 continue
      as(mt+1)=as(mt)+(asb(1)+asb(mt))*dbth/2.
      alfpsi=pi2/as(mt+1)
ccc
ccc  calculate q
ccc
      ce=(elip**2)/2.0
      q=0.0
      do 21 it=1,mt
      xbs=xb(it)**2
      zbs=zb(it)**2
      grps=xbs*zbs+(zbs+ce*(xbs-1.))**2
      q=q+asb(it)/(xbs*sqrt(grps))
   21 continue
      q=elip*q00*q/mt
c
cccc   calculate x, z, aj, aj3 on the equal arc grids
      darc=as(mt+1)/(2.*(nthe-1.))
      do 60 i=2,nthe
      sarc=(i-1-0.5)*darc
      do 70 it=1,mt
      if(sarc.ge.as(it).and.sarc.lt.as(it+1)) go to 71
      go to 70
   71 continue
      theta=bet(it)+dbth*(sarc-as(it))/(as(it+1)-as(it))
      go to 61
   70 continue
   61 continue
      ii=i+2
      xx=sqrt(1.+xr*cos(theta))
      zz=zr*sin(theta)/xx
      xsq=xx**2
      zsq=zz**2
      ajs=(zsq+(xsq-1.)*elip**2/2.)**2+xsq*zsq
      aj3(ii)=elip*q00/(alfpsi*sqrt(ajs))
   60 continue
ccc
      do 30 i=1,nthe
      sarc=(i-1)*darc
      do 40 it=1,mt
      if(sarc.ge.as(it).and.sarc.lt.as(it+1)) go to 41
      go to 40
   41 continue
      theta=bet(it)+dbth*(sarc-as(it))/(as(it+1)-as(it))
      go to 31
   40 continue
   31 continue
      ii=i+2
      x(ii)=sqrt(1.+xr*cos(theta))
      z(ii)=zr*sin(theta)/x(ii)
      xsq=x(ii)**2
      zsq=z(ii)**2
      ajs=(zsq+(xsq-1.)*elip**2/2.)**2+xsq*zsq
      aj(ii)=elip*q00/(alfpsi*sqrt(ajs))
   30 continue
ccc
ccc  symmetry conditions
cccc
      x(2)=x(4)
      x(1)=x(5)
      x(nthe3)=x(nthe1)
      x(nthe4)=x(nthe)
      z(2)=-z(4)
      z(1)=-z(5)
      z(nthe3)=-z(nthe1)
      z(nthe4)=-z(nthe)
      aj(2)=aj(4)
      aj(1)=aj(5)
      aj(nthe3)=aj(nthe1)
      aj(nthe4)=aj(nthe)
      aj3(3)=aj3(4)
      aj3(2)=aj3(5)
      aj3(1)=aj3(6)
      aj3(nthe3)=aj3(nthe2)
      aj3(nthe4)=aj3(nthe1)
      return
      end
c
cccccc****************************************
ccccccc****************************************
c
      subroutine intspl(n,npp2,f,ff,xg,cf)
      include 'gridparam'
      parameter(nn=nmbrgrd,nnp2=nn+2)
      common/splinc/gspl(nnp2,nn),gsplp(nnp2,nn)
      dimension f(n),ff(n),xg(n),cf(npp2)
c
      call depose(f,cf)
c
      do 1 i=1,n
      ff(i)=0.0
      do 2 isp=1,npp2
      ff(i)=ff(i)+cf(isp)*gspl(isp,i)
    2 continue
    1 continue
c
      return
      end
      subroutine intsplp(n,npp2,f,ff,xg,cf)
      include 'gridparam'
      parameter(nn=nmbrgrd,nnp2=nn+2)
      common/splinc/gspl(nnp2,nn),gsplp(nnp2,nn)
      dimension f(n),ff(n),xg(n),cf(npp2)
c
      call depose(f,cf)
c
      do 1 i=1,n
      ff(i)=0.0
      do 2 isp=1,npp2
      ff(i)=ff(i)+cf(isp)*gsplp(isp,i)
    2 continue
    1 continue
c
      return
      end
      subroutine splprp(npsi)
      include 'gridparam'
      parameter(n=nmbrgrd,np2=n+2)
      common/rcl/r(n)
      common/sacl/sa(4,4,np2)
      common/phi01cl/p01(n)/phi02cl/p02(n)/phi03cl/p03(n)
      common/phi11cl/p11(n)/phi12cl/p12(n)/phi13cl/p13(n)
      common/phi21cl/p21(n)/phi22cl/p22(n)/phi23cl/p23(n)
      common/phii1cl/pi1(n)/phii2cl/pi2(n)
      common/phii3cl/pi3(n)/phii4cl/pi4(n)
      common/ncl/nm2,nm1,nnn,np1,nnnp2
      common/phi4cl/dum4(np2)
      common/phi5cl/dum5(np2)
ccc
ccc  set up spline grids
c
      nm1=npsi-1
      nm2=npsi-2
      nnn=npsi
      np1=npsi+1
      nnnp2=npsi+2
      dr=1./nm1
      do 1 i=1,nnn
      r(i)=(i-1.)*dr
    1 continue
      r(nnn)=1.0
c
      call bsplcof
      call deset(0)
c
      return
      end
      subroutine splprp1(n,npp2,xg)
      include 'gridparam'
      parameter(nn=nmbrgrd,nnp2=nn+2)
      common/splinc/gspl(nnp2,nn),gsplp(nnp2,nn)
      dimension xg(n)
c
      do 1 i=1,n
      do 2 isp=1,npp2
      call spval(isp,xg(i),spl,1)
      call spvalp(isp,xg(i),splp,1)
      gspl(isp,i)=spl
      gsplp(isp,i)=splp
    2 continue
    1 continue
      return
      end
      subroutine spval(i,x,y,m)
c---computes m values of the i-th b-spline at the m x-values and
c   returns them in y. note that the x-values are assumed to be
c---ordered.
c---
      include 'gridparam'
      parameter(npp=nmbrgrd+2)
      dimension x(m),y(m)
      common/rcl/r(nmbrgrd)
      common/sacl/sa(4,4,npp)
      common/ncl/nm2,nm1,n,np1,np2
      kk=1
      do 10 j=1,m
      y(j)=0.0
      do 11 k=kk,n
      if(x(j).le.r(k)) go to 12
   11 continue
   12 l=max0(2,k)-i+3
      kk=k
      if(l.lt.1.or.l.gt.4) go to 10
      y(j)=(sa(1,l,i)+x(j)*(sa(2,l,i)+x(j)*(sa(3,l,i)+x(j)*sa(4,l,i))))
   10 continue
      return
      entry spvalp(i,x,y,m)
      kk=1
      do 20 j=1,m
      y(j)=0.0
      do 21 k=kk,n
      if(x(j).le.r(k)) go to 22
   21 continue
   22 l=max0(2,k)-i+3
      kk=k
      if(l.lt.1.or.l.gt.4) go to 20
      y(j)=sa(2,l,i)+x(j)*(2.0*sa(3,l,i)+3.0*x(j)*sa(4,l,i))
   20 continue
      return
      entry spvalpp(i,x,y,m)
      kk=1
      do 30 j=1,m
      y(j)=0.0
      do 31 k=kk,n
      if(x(j).le.r(k)) go to 32
   31 continue
   32 l=max0(2,k)-i+3
      kk=k
      if(l.lt.1.or.l.gt.4) go to 30
      y(j)=2.0*sa(3,l,i)+6.0*sa(4,l,i)*x(j)
   30 continue
      return
      end
      subroutine depose(fn,c)
      include 'gridparam'
      dimension fn(1),c(1)
      common/rcl/r(nmbrgrd)
      common/ncl/nm2,nm1,n,np1,np2
      common/phi01cl/p01(nmbrgrd)/phi02cl/p02(nmbrgrd)/phi03cl
     &     /p03(nmbrgrd)
      common/phi11cl/p11(nmbrgrd)/phi12cl/p12(nmbrgrd)/phi13cl
     &     /p13(nmbrgrd)
      common/phi4cl/e(nmbrgrd+2)/phi5cl/f(nmbrgrd+2)
      nm3=n-3
      if(abs(fn(1)).gt.abs(10.*fn(2))) go to 20
c---find derivative of fcn at end pts.
       d11=(fn(2)-fn(1))/(r(2)-r(1))
       d12=(fn(3)-fn(2))/(r(3)-r(2))
       d13=(fn(4)-fn(3))/(r(4)-r(3))
       d21=(d12-d11)/(r(3)-r(1))
       d22=(d13-d12)/(r(4)-r(2))
       d31=(d22-d21)/(r(4)-r(1))
      fp0=d11-d21*r(2)+d31*r(2)*r(3)
      c(1)=fn(1)/p01(1)
   22 continue
c---
       d11=(fn(nm2)-fn(nm3))/(r(nm2)-r(nm3))
       d12=(fn(nm1)-fn(nm2))/(r(nm1)-r(nm2))
       d13=(fn(n)-fn(nm1))/(r(n)-r(nm1))
       d21=(d12-d11)/(r(nm1)-r(nm3))
       d22=(d13-d12)/(r(n)-r(nm2))
       d31=(d22-d21)/(r(n)-r(nm3))
       fp1=+r(n)*(d13+(r(n)-r(nm1))*(d22+d31*(r(n)-r(nm2))))
c---compute e and f.
      f(2)=(fp0-p11(1)*c(1))/p12(1)
      e(2)=0.0
      e(3)=-p03(2)/p02(2)
      f(3)=(fn(2)-p01(2)*f(2))/p02(2)
      do 10 i=4,n
       e(i)=1.0/(p01(i-1)*e(i-1)+p02(i-1))
       f(i)=(fn(i-1)-p01(i-1)*f(i-1))*e(i)
       e(i)=-p03(i-1)*e(i)
   10 continue
c---compute c.
      c(np2)=fn(n)/p03(n)
      c(np1)=(fp1-p13(n)*c(np2))/p12(n)
       do 12 i=2,n
       j=np2-i
       c(j)=e(j)*c(j+1)+f(j)
   12 continue
      return
   20 continue
       d11=(fn(3)-fn(2))/(r(3)-r(2))
       d12=(fn(4)-fn(3))/(r(4)-r(3))
       d13=(fn(5)-fn(4))/(r(5)-r(4))
       d21=(d12-d11)/(r(4)-r(2))
       d22=(d13-d12)/(r(5)-r(3))
       d31=(d22-d21)/(r(5)-r(2))
      fn0=fn(2)+(r(1)-r(2))*(d11+(r(1)-r(3))*(d21
     2                          +(r(1)-r(4))* d31))
      fp0=d11+d21*((r(1)-r(3))+(r(1)-r(2)))+d31*((r(1)-r(3))*
     2    (r(1)-r(4))+(r(1)-r(2))*(r(1)-r(4))+
     3 (r(1)-r(2))*(r(1)-r(3)))
      c(1)=fn0/p01(1)
      go to 22
      end
      subroutine repose(c,v)
      include 'gridparam'
      common/rcl/r(nmbrgrd)
      dimension c(1),v(1)
      common/ncl/nm2,nm1,n,np1,np2
      common/phi01cl/p01(nmbrgrd)/phi02cl/p02(nmbrgrd)/phi03cl
     &     /p03(nmbrgrd)
      common/phi11cl/p11(nmbrgrd)/phi12cl/p12(nmbrgrd)/phi13cl
     &     /p13(nmbrgrd)
      common/phi21cl/p21(nmbrgrd)/phi22cl/p22(nmbrgrd)/phi23cl
     &     /p23(nmbrgrd)
      common/phii1cl/pi1(nmbrgrd)/phii2cl/pi2(nmbrgrd)/phii3cl
     &     /pi3(nmbrgrd)
      common/phii4cl/pi4(nmbrgrd)
c---computes the function values at the test points from the spline coef
      do 10 i=1,n
      v(i)=c(i)*p01(i)+c(i+1)*p02(i)+c(i+2)*p03(i)
   10 continue
      return
      entry reposp(c,v)
      do 11 i=1,n
      v(i)=c(i)*p11(i)+c(i+1)*p12(i)+c(i+2)*p13(i)
   11 continue
      return
      entry repospp(c,v)
      do 12 i=1,n
      v(i)=c(i)*p21(i)+c(i+1)*p22(i)+c(i+2)*p23(i)
   12 continue
      return
      entry reposi(c,v)
      v(1)=0.0
      do 13 i=2,n
      v(i)=v(i-1)+c(i-1)*pi1(i)+c(i)*pi2(i)+c(i+1)*pi3(i)+
     2            c(i+2)*pi4(i)
   13 continue
      return
      end
      subroutine splerr(c,er)
      dimension er(1),c(1)
c---relative error estimate of spline fit.
c---note# the routine only returns a value in er(i) if
c---      the error is greater than the initial value of er(i).
      include 'gridparam'
      common/rcl/r(nmbrgrd)
      parameter(npp=nmbrgrd+2)
      common/ncl/nm2,nm1,n,np1,np2
      common/sacl/sa(4,4,npp)
      fcn(i,l,x)=sa(1,l,i)+x*(sa(2,l,i)+x*(sa(3,l,i)+x*sa(4,l,i)))
      fval=0.0
      do 11 i=1,n
      fval=fval+abs(c(i)*fcn(i,4,r(i))+c(i+1)*fcn(i+1,3,r(i))+
     2 c(i+2)*fcn(i+2,2,r(i)))
   11 continue
      fval=fval/n
      if(fval.eq.0.0) fval=1.
      do 10 i=2,nm1
      err=0.02625*(c(i-1)*(           -sa(4,4,i-1))
     2              +c(i)  *(sa(4,4,i)  -sa(4,3,i)  )
     3              +c(i+1)*(sa(4,3,i+1)-sa(4,2,i+1))
     4              +c(i+2)*(sa(4,2,i+2)-sa(4,1,i+2))
     5              +c(i+3)*(sa(4,1,i+3)))
     6         *(amax1(r(i+1)-r(i),r(i)-r(i-1))**3)
      err=abs(err/fval)
      if(err.gt.er(i))er(i)=err
   10 continue
      return
      end
      subroutine bsplcof
c---routine computes the coefficients of the b-splines which form
c   a basis over the set of knots, r, with repeated knots
c   at the end points.  s(j,l,i) j-power, l-segment, i-spline.
c---
      include 'gridparam'
      common/rcl/r(nmbrgrd)
      parameter(npp=(nmbrgrd+2)*16)
      common/ncl/nm2,nm1,n,np1,np2
      common/sacl/sa(npp)
c---i=1
      c4=4.00/((r(2)-r(1))**4)
      idx=16
      sa(idx)   =              -c4
      sa(idx-1) =           3.0*c4*r(2)
      sa(idx-2) =          -3.0*c4*r(2)*r(2)
      sa(idx-3) =               c4*r(2)*r(2)*r(2)
c---i=2
      c4=4.0/((r(3)-r(2))*(r(3)-r(1))**3)
      c3=4.0/((r(2)-r(3))*(r(2)-r(1))**3)
      idx=32
      sa(idx)   =              -c4
      sa(idx-1) =           3.0*c4*r(3)
      sa(idx-2) =          -3.0*c4*r(3)*r(3)
      sa(idx-3) =               c4*r(3)*r(3)*r(3)
      sa(idx-4) =sa(idx)       -c3
      sa(idx-5) =sa(idx-1) +3.0*c3*r(2)
      sa(idx-6) =sa(idx-2) -3.0*c3*r(2)*r(2)
      sa(idx-7) =sa(idx-3)     +c3*r(2)*r(2)*r(2)
c---i=3
      c4=4.0/((r(4)-r(2))*(r(4)-r(3))*(r(4)-r(1))**2)
      c3=4.0/((r(3)-r(2))*(r(3)-r(4))*(r(3)-r(1))**2)
      c2=4.0/((r(2)-r(3))*(r(2)-r(4))*(r(2)-r(1))**2)
      idx=48
      sa(idx)   =              -c4
      sa(idx-1) =          +3.0*c4*r(4)
      sa(idx-2) =          -3.0*c4*r(4)*r(4)
      sa(idx-3) =               c4*r(4)*r(4)*r(4)
      sa(idx-4) =sa(idx)       -c3
      sa(idx-5) =sa(idx-1) +3.0*c3*r(3)
      sa(idx-6) =sa(idx-2) -3.0*c3*r(3)*r(3)
      sa(idx-7) =sa(idx-3)     +c3*r(3)*r(3)*r(3)
      sa(idx-8) =sa(idx-4)     -c2
      sa(idx-9) =sa(idx-5) +3.0*c2*r(2)
      sa(idx-10)=sa(idx-6) -3.0*c2*r(2)*r(2)
      sa(idx-11)=sa(idx-7)     +c2*r(2)*r(2)*r(2)
c---i=4,n-1
      do 10 i=4,nm1
      m1=i-3
      m2=i-2
      m3=i-1
      m4=i
      m5=i+1
      c4=4.0/((r(m5)-r(m1))*(r(m5)-r(m2))*(r(m5)-r(m3))*(r(m5)-r(m4)))
      c3=4.0/((r(m4)-r(m1))*(r(m4)-r(m2))*(r(m4)-r(m3))*(r(m4)-r(m5)))
      c2=4.0/((r(m3)-r(m1))*(r(m3)-r(m2))*(r(m3)-r(m4))*(r(m3)-r(m5)))
      c1=4.0/((r(m2)-r(m1))*(r(m2)-r(m3))*(r(m2)-r(m4))*(r(m2)-r(m5)))
      idx=16*i
      sa(idx)   =              -c4
      sa(idx-1) =          +3.0*c4*r(m5)
      sa(idx-2) =          -3.0*c4*r(m5)*r(m5)
      sa(idx-3) =               c4*r(m5)*r(m5)*r(m5)
      sa(idx-4) =sa(idx)       -c3
      sa(idx-5) =sa(idx-1) +3.0*c3*r(m4)
      sa(idx-6) =sa(idx-2) -3.0*c3*r(m4)*r(m4)
      sa(idx-7) =sa(idx-3)     +c3*r(m4)*r(m4)*r(m4)
      sa(idx-8) =sa(idx-4)     -c2
      sa(idx-9) =sa(idx-5) +3.0*c2*r(m3)
      sa(idx-10)=sa(idx-6) -3.0*c2*r(m3)*r(m3)
      sa(idx-11)=sa(idx-7)     +c2*r(m3)*r(m3)*r(m3)
      sa(idx-12)=sa(idx-8)     -c1
      sa(idx-13)=sa(idx-9) +3.0*c1*r(m2)
      sa(idx-14)=sa(idx-10)-3.0*c1*r(m2)*r(m2)
      sa(idx-15)=sa(idx-11)    +c1*r(m2)*r(m2)*r(m2)
   10 continue
c---i=n
      m1=n-3
      m2=n-2
      m3=n-1
      m4=n
      c4=12.0/((r(m4)-r(m1))*(r(m4)-r(m2))*(r(m4)-r(m3)))
      c3=-4.0/((r(m4)-r(m2))*(r(m4)-r(m3))*(r(m4)-r(m1))**2)
     2   -4.0/((r(m4)-r(m1))*(r(m4)-r(m3))*(r(m4)-r(m2))**2)
     3   -4.0/((r(m4)-r(m1))*(r(m4)-r(m2))*(r(m4)-r(m3))**2)
      c2= 4.0/((r(m3)-r(m1))*(r(m3)-r(m2))*(r(m3)-r(m4))**2)
      c1= 4.0/((r(m2)-r(m1))*(r(m2)-r(m3))*(r(m2)-r(m4))**2)
      idx=16*n
      sa(idx-4) =              -c3
      sa(idx-5) =c4        +3.0*c3*r(m4)
      sa(idx-6) = -2.0*c4  -3.0*c3*r(m4)*r(m4)
      sa(idx-7) =c4*r(m4)*r(m4)+c3*r(m4)*r(m4)*r(m4)
      sa(idx-8) =sa(idx-4)     -c2
      sa(idx-9) =sa(idx-5) +3.0*c2*r(m3)
      sa(idx-10)=sa(idx-6) -3.0*c2*r(m3)*r(m3)
      sa(idx-11)=sa(idx-7)     +c2*r(m3)*r(m3)*r(m3)
      sa(idx-12)=sa(idx-8)     -c1
      sa(idx-13)=sa(idx-9) +3.0*c1*r(m2)
      sa(idx-14)=sa(idx-10)-3.0*c1*r(m2)*r(m2)
      sa(idx-15)=sa(idx-11)    +c1*r(m2)*r(m2)*r(m2)
c---i=n+1
      c4= 12.0/((r(m4)-r(m2))    *(r(m4)-r(m3)))
      c3=-12.0/((r(m4)-r(m3))    *(r(m4)-r(m2))**2)
     2   -12.0/((r(m4)-r(m2))    *(r(m4)-r(m3))**2)
      c2=  4.0/((r(m4)-r(m3))    *(r(m4)-r(m2))**3)
     2    +4.0/((r(m4)-r(m2))**2 *(r(m4)-r(m3))**2)
     3    +4.0/((r(m4)-r(m2))    *(r(m4)-r(m3))**3)
      c1=  4.0/((r(m3)-r(m2))    *(r(m3)-r(m4))**3)
      idx=16*(n+1)
      sa(idx-8) =                        -c2
      sa(idx-9) =          c3        +3.0*c2*r(m4)
      sa(idx-10)=-c4  -2.0*c3*r(m4)  -3.0*c2*r(m4)*r(m4)
      sa(idx-11)= c4*r(m4)+c3*r(m4)*r(m4)+c2*r(m4)*r(m4)*r(m4)
      sa(idx-12)=sa(idx-8)     -c1
      sa(idx-13)=sa(idx-9) +3.0*c1*r(m3)
      sa(idx-14)=sa(idx-10)-3.0*c1*r(m3)*r(m3)
      sa(idx-15)=sa(idx-11)    +c1*r(m3)*r(m3)*r(m3)
c---i=n+2
      c4= +4.0/(r(m4)-r(m3))
      c3=-12.0/(r(m4)-r(m3))**2
      c2=+12.0/(r(m4)-r(m3))**3
      c1= -4.0/(r(m4)-r(m3))**4
      idx=16*(n+2)
      sa(idx-12)=                          -c1
      sa(idx-13)=            c2        +3.0*c1*r(m4)
      sa(idx-14)=  -c3  -2.0*c2*r(m4)  -3.0*c1*r(m4)*r(m4)
      sa(idx-15)=c4+c3*r(m4)+c2*r(m4)*r(m4)+c1*r(m4)*r(m4)*r(m4)
      return
      end
      subroutine deset(jh)
      include 'gridparam'
      common/rcl/r(nmbrgrd)
      common/ncl/nm2,nm1,n,np1,np2
      common/phi01cl/p01(nmbrgrd)/phi02cl/p02(nmbrgrd)/phi03cl
     &     /p03(nmbrgrd)
      common/phi11cl/p11(nmbrgrd)/phi12cl/p12(nmbrgrd)/phi13cl
     &     /p13(nmbrgrd)
      common/phi21cl/p21(nmbrgrd)/phi22cl/p22(nmbrgrd)/phi23cl
     &     /p23(nmbrgrd)
      common/phii1cl/pi1(nmbrgrd)/phii2cl/pi2(nmbrgrd)/phii3cl
     &     /pi3(nmbrgrd)
      common/phii4cl/pi4(nmbrgrd)
      common/bvpcl/d1,d2,d3,d4,d5,d6
c      data pi1(1),pi2(1),pi3(1),pi4(1)/4*0.0/
      pi1(1)=0.
      pi2(1)=0.
      pi3(1)=0.
      pi4(1)=0.
c---subroutine sets up values of splines at the knots.
      do 10 i=1,n
      call spval(i,r(i),p01(i),1)
      call spval(i+1,r(i),p02(i),1)
      call spval(i+2,r(i),p03(i),1)
      call spvalp(i,r(i),p11(i),1)
      call spvalp(i+1,r(i),p12(i),1)
      call spvalp(i+2,r(i),p13(i),1)
      call spvalpp(i,r(i),p21(i),1)
      call spvalpp(i+1,r(i),p22(i),1)
      call spvalpp(i+2,r(i),p23(i),1)
   10 continue
      p02(1)=0.
      p03(1)=0.
      p01(n)=0.
      p02(n)=0.
      p13(1)=0.
      p11(n)=0.
      do 12 i=2,n
      pi1(i)=gauss6(i-1,i,jh)
      pi2(i)=gauss6(i,i,jh)
      pi3(i)=gauss6(i+1,i,jh)
      pi4(i)=gauss6(i+2,i,jh)
   12 continue
      d1=p11(1)
      d2=p12(1)
      d3=p12(n)
      d4=p13(n)
      d5=p03(n)
      d6=p01(1)
      return
      end
      function gauss6(k,i,jh)
      dimension f(6),x(6)
      include 'gridparam'
      common/rcl/r(nmbrgrd)
      data w1,w2,w3/.467913934572691,.360761573048139,.171324492379170/
      data x1,x2,x3/.238619186083197,.661209386466265,.932469514203152/
      rd=0.5*(r(i)-r(i-1))
      rp=0.5*(r(i)+r(i-1))
      x(1)=rp-rd*x3
      x(2)=rp-rd*x2
      x(3)=rp-rd*x1
      x(4)=rp+rd*x1
      x(5)=rp+rd*x2
      x(6)=rp+rd*x3
      call spval(k,x(1),f(1),6)
      gauss6=rd*(w3*(x(1)**jh*f(1)+x(6)**jh*f(6))
     2          +w2*(x(2)**jh*f(2)+x(5)**jh*f(5))
     4          +w1*(x(3)**jh*f(3)+x(4)**jh*f(4)))
      return
      end
c