      program fbt
c
c               Free Boundary Tokamak Equilibrium Code (FBT)
c
c
c
c    Reference: F. Hofmann, Computer Physics Communications 48 (1988) 207
c
c
c
c  1) Input file "fbtd" (in namelist format, identifier = $data)
c
c
c                            control parameters
c
c  neqtcv              equilibrium number
c  nruns(1)            number of equilibria computed in one run
c  ievolv(0)           =0 no automatic parameter evolution
c                      =1 automatic parameter evolution
c  testa(1.0e-4)       convergence criterion
c
c                              computational grid
c
c  ri(0.58)            inner boundary of computational grid
c  ro(1.18)            outer boundary of computational grid
c  zu(0.794839)        upper boundary of computational grid
c  -zu                 lower boundary of computational grid
c  nr(30)              nr+1 = number of radial mesh points
c  nz(32)              (2*nz)+1 = number of axial mesh points
c
c                              boundary points
c
c  ilie                number of exact boundary points
c  rlim1(i),zlim1(i)   initial coordinates of exact boundary points
c  rlim2(i),zlim2(i)   final coordinates of exact boundary points
c  ilia                number of approximate boundary points
c  rlia1(i),zlia1(i)   initial coordinates of approximate boundary points
c  rlia2(i),zlia2(i)   final coordinates of approximate boundary points
c  eftlim              if((ilie.le.4).and.(ilia.le.4)),the initial estimate
c  ritlim                   for the plasma current distribution is set up
c  botlim                   in the rectangular area:
c  toplim                   eftlim < r < ritlim, botlim < z < toplim
c
c                                   X-points
c
c  ibro                number of points where B-r=0
c  rbro(i),zbro(i)     coordinates of B-r=0 points
c  ibzo                number of points where B-z=0
c  rbzo(i),zbzo(i)     coordinates of B-z=0 points
c                        
c                                 plasma current
c
c  placu1              initial plasma current (kA)
c  placu2              final plasma current (kA)
c  placex              if ievolv=1, plasma current is given by:
c                       Ip = placu1+(placu2-placu1)*frac**placex,
c                       where frac=(nrun-1)/(nruns-1)
c
c                              poloidal field coils
c
c  icur                number of coils
c  rvf(i),zvf(i)       coordinates of coil positions (and vessel filaments)
c  hvf(i),wvf(i)       height and width of coils (and vessel filaments)
c  tvf(i)              sign non integer number of turns for coils (and vessel filaments)
c  cvf(i)              predetermined coil currents (in kA)
c  icvx                number coil groups whose currents must be adjusted
c                       to obtain the desired plasma shape
c  cvx1(i)-cvx16(i)    cvxk(i) is the number of turns of the i-th coil
c                       in the k-th group
c
c                                  Volt-seconds
c
c  ivsec               if ivsec=1, the Volt secs due to all coil currents
c                        in the center of the plasma is given by vsec
c  vsec                prescribed Volt secs
c
c                                source functions
c
c  ttfac1,ttfac2       scaling coefficient for TT'
c  emm1,emm2           initial and final lowest exponent in TT' polinomial
c                       (emm1>1.,emm2>1.)
c  omega               parameter to modify li
c  ppfac               scaling coefficient for p' (ppfac > 1.)
c  ell1,ell2           initial and final lowest exponent in p' polinomial
c                       (ell1>1.,ell2>1.)
c  ppal                pressure profile parameter (0 < ppal < 1)
c
c          note: source functions are defined as follows:
c
c          ppri=ppal*phell+(1.-ppal)*(phell-phell1),
c               where phell=phi**ell, phell1=phell*phi
c          ttpr=phemm+ttbe*(phemm1-phemm+omega*(phemm2-phemm1))
c               where phemm=phi**emm, phemm1=phemm*phi,phemm2=phemm1*phi
c
c                          other physics parameters
c
c  bzeru               vacuum toroidal field half way between inner and
c                            outer edge of plasma
c  bzero	       vacuum toroidal field at R=0.88m  
c  qzero               q on magnetic axis
c  psirat              if the plasma is diverted, the plasma boundary is
c                       defined at psirat*100% of the separatrix flux
c
c              initial guess for plasma current distribution
c
c  capaj               parameter to vary the peaking factor of the initially
c                       assumed plasma current distribution
c
c                              least squares fit
c
c  weitam,weitex       parameters to adjust the weights of the approximate
c                       plasma boundary points automatically according to
c                       their vertical distance from the magnetic axis
c  dissi               parameter determining trade-off between shape
c                       accuracy and power dissipation in shaping coils
c  strki               parameter to discourage large differences between
c                       adjacent coil currents
c
c                                 measurements
c
c  nprob               number of magnetic field probes
c  rbee(n),zbee(n)     coordinates of magnetic field probes
c  thbee(n)            angle of probe with respect to horizontal axis
c  nloop               number of flux loops
c  rflu(n),zflu(n)     coordinates of flux loops
c
c                             analytic plasma boundary
c
c  iansha              =0 shape is given by coordinates of boundary points
c                        =1 shape is specified analytically
c  rmajo1,rmajo2       initial and final major radius (for iansha=1)
c  rmino1,rmino2       initial and final minor radius (for iansha=1)
c  cappa1,cappa2       initial and final elongation   (for iansha=1)
c  delta1,delta2       initial and final triangularity (for iansha=1)
c  hlamd1,hlamd2       initial and final shape parameter (for iansha=1)
c  zmajo1,zmajo2       initial and final z-coordinate of magnetic axis
c                            (iansha=1)
c
c   note: the analytic plasma boundary is defined by
c         r(i) = rmajo+rmino*cos(w(i)+delta*sin(w(i))-hlamd*sin(2.*w(i)))
c         z(i) = zmajo+rmino*cappa*sin(w(i))
c
c
c
c
c  2) Input file "fbtdcor" (in namelist format, identifier = $datacor)
c
c  npsi                number of flux surfaces in plasma
c  ilarg               if ilarg=1 output file "fbtequi" is written
c                        (standard output file)
c  ixdr                if ixdr=0 output file in FORTRAN (fbtequi)
c                         if ixdr=1 output file in XDR (fbtequi.xdr)
c  inova               if inova=1 output file "fbtnova" is written
c                        (input file for the NOVA-W code)
c  ierat               if ierat=1 output file "fbterat" is written
c                        (input file for ERATO)
c  ifour               number of Fourier terms used in the analytic
c                       representation of the plasma boundary. This is
c                       only used to create the ERATO input file
c  ikriz               if ikriz=1 output file "fbtkriz" is written
c                        (input file for A. Kritz's ray tracing code)
c  imeas               if imeas=1 output file "fbtmeas" is written
c                        (input file for the MGAMS code)
c  iwrida              if iwrida=1 output is displayed on the terminal
c  nvvel               number of vessel elements (must be 19 or 38)
c  nfilae              number of filaments in each E-coil (TCV: 1 or 2 or 4)
c  nfilaf              number of filaments in each F-coil (TCV: 1 or 2)
c  nfast               number of coils in fast group      (TCV: 2 or 6)
c
c
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
      namelist/data/ri,ro,zu,nr,nz,ilie,rlim1,rlim2,zlim1,zlim2,
     1 ilia,rlia1,rlia2,zlia1,zlia2,eftlim,ritlim,botlim,toplim,
     1 ibro,rbro,zbro,ibzo,rbzo,zbzo,placu1,placu2,placex,icur,rvf,
     2 zvf,cvf,ivsec,vsec,ttfac1,ttfac2,ell1,ell2,emm1,emm2,ppfac,ppal,
     2 itamax,testa,mei,ncosei,isaddl,ievolv,nruns,
     3 capaj,bzero,qzero,psirat,weitam,weitex,omega,
     3 icvx,cvx1,cvx2,cvx3,cvx4,cvx5,cvx6,cvx7,cvx8,cvx9,cvx10,cvx11,
     4 cvx12,cvx13,cvx14,cvx15,cvx16,icoilon,zeecorb,zeecort,
     4 dissi,strki,rbee,zbee,thbee,nprob,rflu,zflu,nloop,neqtcv,
     5 iansha,rmajo1,rmajo2,rmino1,rmino2,cappa1,cappa2,delta1,delta2,
     5 hlamd1,hlamd2,zmajo1,zmajo2,timeeq,inova,ierat,iwrida,ides,
     6 imeas,ilarg,ikriz,npsi,size,xor,yor,itcvpl,ifour,iplot,ixdr,
     6 nvvel,nfilae,nfilaf,nfast,rampt,nelz,xip,xop,zlp,zup,gain,
     7 ddd,eee,nmatrx,istop,iscramb,iprcinc,ohsame,ohsami,ipra,iprml,
     7 ipreps,iprvol,iprx,vscal,wscal,uc1,uc2,uc3,iprmea,iprmax,ipr16,
     8 ipripz,iscale,rshift1,rshift2,rshift3,zshift,rliaa,zliaa,ilimre,
     8 ohback,veback,gapin,gapout,emgain,aipgain,boost,amagic,zeecorr,
     9 ntoft,toft,deltat,hpla,wpla,f36same,psifac,numeq,relovo,ggain,
     9 ntmax,diohdt,mvloop,timefac,flattop,iohfb,if36fb,ohcorr,hgain,
     9 gainr,gainz,gainvz,gainze,gainvze,deltar,egain,fgain,ohgain,
     9 bzpzero,bzmzero,curfac,aipdel,dpszero,dpsfac,delipz,aipipz,
     9 brpzero,brmzero,bzpzerb,bzmzerb,brpzerb,brmzerb,
     9 gainext,fastm,alpha,nzaxel,nzaxre
c
      integer ishot
      character*8 strshot
      ishot = 1
      call get_command_argument(ishot,strshot)
      read (strshot,*) ishot
c
      call fbtmdsopen    (ishot)
      call fbtmdsgetconst(ishot)
      call fbtmdsgetvar  (ishot,1)
      call tcvdata
      if(iwrida.eq.1) write(6,data)
c
      do 2000 nrum = 1,nruns
      nrun = nrum
      if((nruns.eq.1).or.(ievolv.eq.1).or.(nrun.eq.1))   go to 30
      call fbtmdsgetvar(ishot,nrun)
      if(iwrida.eq.1) write(6,data)
   30 continue
c
      call start
      call check
      call psiext
c
      print 102
  102 format(71h  j1 k1 j2 k2     rax         zax         psax        cy
     1(2)       cy(3))
c
      do 500 ita = 1,itamax
      call psipla
      call lagrang
      call psinew
      call psiins
      call curden
c
      chmm = 0.
      do 415 j = 2,nr
      do 414 k = 2,km
      ch = abs(psi(j,k)-psn(j,k))
      if (ch-chmm)  413,413,412
  412 chmm = ch
  413 psn(j,k) = psi(j,k)
  414 continue
  415 continue
      ta = chmm/psi(jc,kc)
      if (ta-testa)  510,500,500
  500 continue
c
      print 502,nrun
  502 format (47h no convergence in outer loop for equilibrium #,i3)
      go to 2000
  510 print 530,ita
  530 format(5h ita=,i5)
c
1000  call flusurf
      call primat
      if((ilarg.eq.1).and.(ixdr.eq.0)) call wrieq
      if((ikriz.eq.1).and.(nruns.eq.1)) call wrikriz
c     [AM 16/06/2017]
c     Make computation and writing of measurements independent
      call measur
c     if(imeas.gt.0) call wrimeas
      call wrimeas
      if(inova.eq.1) call wrinova
      if(ierat.eq.1) call wrierat
      if((ierat.eq.3).and.(nrun.eq.numeq)) call wrichea
      call pseuliu
      call pseufbt
c
 2000 continue
c
      if ((ishot.eq.-1).or.(ishot.ge.100000)) call fbtmdsput(ishot)
      if (ishot.eq.0) close(unit=55)
c      
      stop
      end
c
      subroutine tcvdata
      include 'fbtmgams.inc'
      include 'fbt.inc'
      dimension rbeeh(38),zbeeh(38),thbeeh(38),rfluh(38),
     1 zfluh(38),rshap16(16),zshap16(16),tshap16(16),
     2 hshap16(16),wshap16(16),
     3 roh(20),zoh(20),toh(20),hoh(20),woh(20),
     4 rfast2(2),zfast2(2),tfast2(2),rfast6(6),zfast6(6),tfast6(6),
     5 hfast2(2),wfast2(2),hfast6(6),wfast6(6),
     6 rvvelb(39),zvvelb(39),wvvelh(38),confil(50),
     7 ifirst(38),ilast(38)
c
      ri     = 0.58
      ro     = 1.18
      zu     = 0.794839
      nr     = 30
      nz     = 32
c
      nfilae = 4
      nfilaf = 2
      filae  = float(nfilae)
      filaf  = float(nfilaf)
      nloop  = 38 
      nprob  = 38
      if(nfast.gt.1) nshafa=19
      if(nfast.le.1) nshafa=18
      ngroup = nshafa + nvvel
      icvx   = 16
      icur   = 8*nfilae + 8*nfilaf
c
      rhocu  = 2.0e-8
      filfac = 0.625
      rhovv  = 0.75e-6
c
      data roh/.4225,.4225,.4225,.4225,.4225,.4225,
     1 .4225,.4225,.4225,.4225,.4225,.4225,
     2 .4458,.4458,.4458,.4458,.6215,.6215,1.1765,1.1765/
      data zoh/-.726,-.594,-.462,-.330,-.198,-.066,
     1 .066,.198,.330,.462,.594,.726,
     2 -.983,-.889,.889,.983,-1.110,1.110,-1.170,1.170/
      data toh/11.917,11.917,11.917,11.917,11.917,11.917,
     1   11.917,11.917,11.917,11.917,11.917,11.917,
     2   14.5,14.5,14.5,14.5,12.0,12.0,8.0,8.0/
      data hoh/.132,.132,.132,.132,.132,.132,.132,.132,.132,.132,
     1   .132,.132,.094,.094,.094,.094,.120,.120,.090,.090/
      data woh/.063,.063,.063,.063,.063,.063,.063,.063,.063,.063,
     1   .063,.063,.1095,.1095,.1095,.1095,.066,.066,.066,.066/
c
c
      data rshap16/.505,.505,.505,.505,.505,.505,.505,.505,
     1 1.3095,1.3095,1.3095,1.3095,1.3095,1.3095,1.3095,1.3095/
      data zshap16/-.7,-.5,-.3,-.1,.1,.3,.5,.7,
     1 -.77,-.61,-.31,-.15,.15,.31,.61,.77/
      data tshap16/34.,34.,34.,34.,34.,34.,34.,34.,
     1   36.,36.,36.,36.,36.,36.,36.,36./
      data hshap16/.180,.180,.180,.180,.180,.180,.180,.180,
     1 .075,.075,.075,.075,.075,.075,.075,.075/
      data wshap16/.051,.051,.051,.051,.051,.051,.051,.051,
     1 .120,.120,.120,.120,.120,.120,.120,.120/
c
c
c
      data rfast6/1.099,1.114,1.129,1.129,1.114,1.129/
      data zfast6/-.648,-.630,-.612,0.612,0.630,0.648/
      data hfast6/.015,.015,.015,.015,.015,.015/
      data wfast6/.015,.015,.015,.015,.015,.015/
      data tfast6/-1.,-1.,-1.,1.,1.,1./
c
      data rfast2/1.114,1.114/
      data zfast2/-.630,0.630/
      data hfast2/.020,.020/
      data wfast2/.020,.020/
      data tfast2/-1.,1./
c
c
      data rbeeh/.606,.606,.606,.606,.606,.606,.606,.6294,.695,.790,
     1  .885,.980,1.0421,1.1188,1.154,1.154,1.154,1.154,1.154,
     2  1.154,1.154,1.154,1.154,1.154,1.154,1.1188,1.0421,.980,.885,
     3  .790,.695,.6294,.606,.606,.606,.606,.606,.606/
      data zbeeh/.000,.115,.230,.345,.460,.575,.690,.7456,.764,.764,
     1  .764,.764,.7348,.6705,.575,.460,.345,.230,.115,0.0,-.115,
     2  -.230,-.345,-.460,-.575,-.6705,-.7348,-.764,-.764,-.764,
     3  -.764,-.7456,-.690,-.575,-.460,-.345,-.230,-.115/
      data thbeeh/90.,90.,90.,90.,90.,90.,90.,45.,0.,0.,0.,0.,-40.,-40.,
     1  -90.,-90.,-90.,-90.,-90.,-90.,-90.,-90.,-90.,-90.,-90.,
     2  -140.,-140.,-180.,-180.,-180.,-180.,-225.,-270.,-270.,-270.,
     3  -270.,-270.,-270./
      data rfluh/.584,.584,.584,.584,.584,.584,.584,.600,.695,.790,.885,
     1  .980,1.0595,1.1361,1.176,1.176,1.176,1.176,1.176,1.176,1.176,
     2  1.176,1.176,1.176,1.176,1.1361,1.0595,.980,.885,.790,.695,
     3  .600,.584,.584,.584,.584,.584,.584/
      data zfluh/.000,.115,.230,.345,.460,.575,.690,.772,.791,.791,
     1  .791,.791,.7555,.6912,.575,.460,.345,.230,.115,0.0,
     2  -.115,-.230,-.345,-.460,-.575,-.6912,-.7555,-.791,-.791,
     3  -.791,-.791,-.772,-.690,-.575,-.460,-.345,-.230,-.115/
      do 5 i=1,nprob
       rbee(i) = rbeeh(i)
       zbee(i) = zbeeh(i)
       thbee(i)= thbeeh(i)
   5  continue
      do 6 i=1,nloop
       rflu(i) = rfluh(i)
       zflu(i) = zfluh(i)
   6  continue
c
c              vacuum vessel element boundaries
c
      data rvvelb/
     1 0.5925,0.5925,0.5925,0.5925,0.5925,0.5925,0.5925,0.5950,
     2 0.6280,0.7425,0.8375,0.9325,1.0130,1.0900,1.1675,1.1675,
     3 1.1675,1.1675,1.1675,1.1675,1.1675,1.1675,1.1675,1.1675,
     4 1.1675,1.1675,1.0900,1.0130,0.9325,0.8375,0.7425,0.6280,
     5 0.5950,0.5925,0.5925,0.5925,0.5925,0.5925,0.5925/
      data zvvelb/
     1 -.0575,0.0575,0.1725,0.2875,0.4025,0.5175,0.6325,0.7470,
     2 0.7800,0.7800,0.7800,0.7800,0.7800,0.7135,0.6470,0.5175,
     3 0.4025,0.2875,0.1725,0.0575,-.0575,-.1725,-.2875,-.4025,
     4 -.5175,-.6470,-.7135,-.7800,-.7800,-.7800,-.7800,-.7800,
     5 -.7470,-.6325,-.5175,-.4025,-.2875,-.1725,-.0575/
      data wvvelh/
     1 .015,.015,.015,.015,.015,.015,.016,.020,.020,.020,
     2 .020,.020,.020,.020,.016,.015,.015,.015,.015,.015,
     3 .015,.015,.015,.015,.016,.020,.020,.020,.020,.020,
     4 .020,.020,.016,.015,.015,.015,.015,.015/
c
c
c
c                       shaping coils
c
c
      do 9 i=1,mcoi
       do 8 k=1,mcog
        cvx(i,k) = 0.
    8  continue
    9 continue
c
      i=1
      do 12 k=1,8
       itype(k) = 1
       isvf(k) = i
       do 10 nf=1,nfilae
        hvf(i) = hshap16(k)/filae
        wvf(i) = wshap16(k)
        rvf(i) = rshap16(k)
        zvf(i) = zshap16(k)-0.5*hshap16(k)+(float(nf)-0.5)*hvf(i)
        tvf(i) = tshap16(k)/filae
        kgroup(i) = k
        cvx(i,k)  = 1.
        i=i+1
   10  continue
       isvl(k) = i-1
   12 continue
c
      do 16 k=9,16
       itype(k) = 1
       isvf(k) = i
       do 14 nf=1,nfilaf
        hvf(i) = hshap16(k)
        wvf(i) = wshap16(k)/filaf
        rvf(i) = rshap16(k)-0.5*wshap16(k)+(float(nf)-0.5)*wvf(i)
        zvf(i) = zshap16(k)
        tvf(i) = tshap16(k)/filaf
        kgroup(i) = k
        cvx(i,k)  = 1.
        i=i+1
   14  continue
       isvl(k) = i-1
   16 continue
c
c
c                     OH coils
c
c
      k=17
      itype(k) = 2
      isvf(k) = i
      do 18 nf=1,12
       rvf(i) = roh(nf)
       zvf(i) = zoh(nf)
       hvf(i) = hoh(nf)
       wvf(i) = woh(nf)
       tvf(i) = toh(nf)
       kgroup(i) = k
       i=i+1
   18 continue
      isvl(k) = i-1
c
      k=18
      itype(k) = 2
      isvf(k) = i
      do 20 nf=13,20
       rvf(i) = roh(nf)
       zvf(i) = zoh(nf)
       hvf(i) = hoh(nf)
       wvf(i) = woh(nf)
       tvf(i) = toh(nf)
       kgroup(i) = k
       i=i+1
   20 continue
      isvl(k) = i-1
c
c
c                       fast coils
c
c
      if(nfast.gt.1) then
       k=19
       itype(k) = 3
       isvf(k) = i
       if(nfast.eq.2) then
        do 22 nf=1,nfast     
         rvf(i) = rfast2(nf)  
         zvf(i) = zfast2(nf)  
         tvf(i) = tfast2(nf)  
         hvf(i) = hfast2(nf)  
         wvf(i) = wfast2(nf)  
         kgroup(i) = k        
         i=i+1  	      
   22   continue	     
       endif
c
       if(nfast.eq.6) then
        do 24 nf=1,nfast     
         rvf(i) = rfast6(nf)  
         zvf(i) = zfast6(nf)  
         tvf(i) = tfast6(nf)  
         hvf(i) = hfast6(nf)  
         wvf(i) = wfast6(nf)  
         kgroup(i) = k        
         i=i+1  	      
   24   continue	     
       endif		     
c
       isvl(k) = i-1
      endif
c
c
c                vacuum vessel filaments
c
c
      do 80 n=1,38
       k=n+nshafa
       itype(k) = 4
       isvf(k)  = i
       wvvel	= wvvelh(n)
       hvvel	= sqrt((rvvelb(n)-rvvelb(n+1))**2
     1  	      +(zvvelb(n)-zvvelb(n+1))**2)
       nfilam	= hvvel/wvvel
       afilam	= wvvel*(hvvel/float(nfilam))
       convv	= 0.
       do 70 nf = 1,nfilam
        hvf(i)   = hvvel/float(nfilam)
        wvf(i)   = wvvel
        filafr   = (float(nf)-0.5)/float(nfilam)
        rvf(i)   = rvvelb(n)+(rvvelb(n+1)-rvvelb(n))*filafr
        zvf(i)   = zvvelb(n)+(zvvelb(n+1)-zvvelb(n))*filafr
        kgroup(i)= k
        confil(nf) = afilam/rvf(i)
        convv = convv + confil(nf)
        i=i+1
   70  continue
       isvl(k) = i-1
       do 75 nf=1,nfilam
        ifil = isvf(k)-1+nf
        tvf(ifil) = confil(nf)/convv
   75  continue
   80 continue
      iallf = i-1
c
      if(nvvel.eq.19) then
       ibeg = isvf(nshafa+1)
       do 85 nh=1,19
        n=2*nh
        k=n+nshafa
        ilast(nh) = isvf(k)+(isvl(k)-isvf(k))/2
        ifirst(nh+1) = ilast(nh)+1
   85  continue
       ifirst(1) = ifirst(20)
       do 87 n=1,19
        k=n+nshafa
        isvf(k) = ifirst(n)
        isvl(k) = ilast(n)
   87  continue
       ishift = iallf-ifirst(20)+1
       do 89 n=1,19
        k=n+nshafa
        isvf(k) = isvf(k)+ishift
        isvl(k) = isvl(k)+ishift
   89  continue
       isvf(nshafa+1) = ibeg
       do 91 is=1,ishift
        rvf(iallf+is) = rvf(iallf-ishift+is)
        zvf(iallf+is) = zvf(iallf-ishift+is)
        hvf(iallf+is) = hvf(iallf-ishift+is)
        wvf(iallf+is) = wvf(iallf-ishift+is)
   91  continue
       do 93 i=ibeg,iallf
        im = iallf+ibeg-i
        if((im-ishift).ge.ibeg) ims=im-ishift
        if((im-ishift).lt.ibeg) ims=im+iallf-ibeg+1
        rvf(im) = rvf(ims)
        zvf(im) = zvf(ims)
        hvf(im) = hvf(ims)
        wvf(im) = wvf(ims)
   93  continue
       do 96 n=1,nvvel
        k=n+nshafa
        ifi = isvf(k)
        ila = isvl(k)
        convv = 0.
        do 94 i=ifi,ila
         nf = i-ifi+1
         confil(nf) = (hvf(i)*wvf(i))/rvf(i)
         convv = convv + confil(nf)
   94   continue
        do 95 i=ifi,ila
         nf = i-ifi+1
         tvf(i) = confil(nf)/convv
         kgroup(i) = k
   95  continue
   96  continue
      endif
c
      return
      end
c
c
c
c
c
      subroutine check
      include 'fbtmgams.inc'
      include 'fbt.inc'
      if(jmax.gt.mrmp)      go to 10
      if(kmax.gt.mamp)      go to 20
      if(icur.gt.mcoi)      go to 30
      if(icvx.gt.mcog)      go to 40
      if(ninsid.gt.minsi)   go to 50
      if(mei.gt.melli)      go to 60
      if(ncosei.gt.mcos)    go to 70
      if(ibro.gt.mbro)      go to 80
      if(ibzo.gt.mbzo)      go to 90
      if(ilia.gt.mabp)      go to 100
      if(icon.gt.icvy)      go to 110
c This check is only useful when ierat=1 (ifour is used for other purposes too!)
c   (SC, 13/04/05)
      if((ifour.gt.mfour).and.(ierat.eq.1))    go to 120
      if(nprob.gt.mbeep)    go to 130
      if(nloop.gt.mflup)    go to 140
      if((npsi+1).gt.mpsis) go to 150
      if(nruns.gt.mequi)    go to 160
      if(ell.le.0.)         go to 170
      if(emm.le.0.)         go to 180
      if(ppfac.le.0.)       go to 190
      if((ppal.le.0.).or.(ppal.gt.1.))     go to 200
      if((nvvel.ne.19).and.(nvvel.ne.38))  go to 210
      if(ievolv.eq.1)       go to 220
      go to 900
   10 print 11
   11 format(27h parameter mrmp too small  )
      go to 800
   20 print 21
   21 format(27h parameter mamp too small  )
      go to 800
   30 print 31
   31 format(27h parameter mcoi too small  )
      go to 800
   40 print 41
   41 format(27h parameter mcog too small  )
      go to 800
   50 print 51
   51 format(27h parameter minsi too small )
      go to 800
   60 print 61
   61 format(27h parameter melli too small )
      go to 800
   70 print 71
   71 format(27h parameter mcos too small  )
      go to 800
   80 print 81
   81 format(27h parameter mbro too small  )
      go to 800
   90 print 91
   91 format(27h parameter mbzo too small  )
      go to 800
  100 print 101
  101 format(27h parameter mabp too small  )
      go to 800
  110 print 111
  111 format(31h ilie or ibro or ibzo too large)
      go to 800
  120 print 121
  121 format(27h parameter mfour too small )
      go to 800
  130 print 131
  131 format(27h parameter mbeep too small )
      go to 800
  140 print 141
  141 format(27h parameter mflup too small )
      go to 800
  150 print 151
  151 format(27h parameter mpsis too small )
      go to 800
  160 print 161
  161 format(27h parameter mequi too small )
      go to 800
  170 print 171
  171 format(31h ell1 and ell2 must be positive)
      go to 800
  180 print 181
  181 format(31h emm1 and emm2 must be positive)
      go to 800
  190 print 191
  191 format(27h ppfac must be positive    )
      go to 800
  200 print 201
  201 format(31h ppal must be between 0. and 1.)
      go to 800
  210 print 211
  211 format(32h nvvel must be equal to 19 or 38)
      go to 800
  220 print *,'Option ievolv=1 removed'
      go to 800
c
  800 stop 8
  900 return
      end
c
c
c
      subroutine lagrang
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
      imat = icvy + icon
      if(icvy.eq.icon) imat = icon
      imau = imat + 1
      if(ilie)  228,228,222
  222 do 226 m1 = 1,ilie
      do 224 m2 = 1,icvy
      hx(m1,m2) = psea(lilia+m1,m2)
  224 continue
      hx(m1,icvz) = -psp(lilia+m1)-pse(lilia+m1)
  226 continue
  228 if(ibro)  236,236,230
  230 do 234 l = 1,ibro
      do 232 m2 = 1,icvy
      hx(l+ilie,m2) = berx(l,m2)
  232 continue
      hx(l+ilie,icvz) = -berp(l)-bere(l)
  234 continue
  236 if(ibzo)  244,244,238
  238 do 242 l = 1,ibzo
      do 240 m2 = 1,icvy
      hx(ilro+l,m2) = bezx(l,m2)
  240 continue
      hx(ilro+l,icvz) = -bezp(l)-beze(l)
  242 continue
  244 if(ivsec)  250,250,246
  246 do 248 m2 = 1,icvy
      hx(icon,m2) = psea(lee,m2)
  248 continue
      hx(icon,icvz) = vsec/(2.*pi)
c
c
c
  250 if(icon-icvy)  251,274,2000
  251 if(icon)  259,259,252
  252 do 258 m1 = 1,icon
      do 254 m2 = 1,icvy
      ax(m1,m2) = hx(m1,m2)
  254 continue
      do 256 m2 = icvz,imat
      ax(m1,m2) = 0.
  256 continue
      ax(m1,imau) = hx(m1,icvz)
  258 continue
  259 do 272 m1 = 1,icvy
      do 266 m2 = 1,icvy
      axp = 0.
      if(ilia)  262,262,260
  260 do 261 m = 1,ilia
      axp = axp + (weit(m)*psea(lmax+m,m2)*psea(lmax+m,m1))
  261 continue
  262 if((m1.eq.icvy).or.(m2.eq.icvy)) go to 265
c
      do 264 i = 1,icur
      axp = axp + (diss(i)*cvx(i,m2)*cvx(i,m1))
      if(i-1)  264,264,263
  263 axp = axp + (cstrk(i)*(cvx(i,m2)-cvx(i-1,m2))
     1                     *(cvx(i,m1)-cvx(i-1,m1)))
  264 continue
  265 continue
      ax(m1+icon,m2) = axp
  266 continue
      do 268 m2 = icvz,imat
      ax(m1+icon,m2) = hx(m2-icvy,m1)
  268 continue
      axp = 0.
      if(ilia)  271,271,269
  269 do 270 m = 1,ilia
      axp = axp-(weit(m)*(psp(lmax+m)+pse(lmax+m))*psea(lmax+m,m1))
  270 continue
  271 ax(m1+icon,imau) = axp
  272 continue
      go to 280
  274 do 278 m1=1,icvy
      do 276 m2=1,icvz
      ax(m1,m2) = hx(m1,m2)
  276 continue
  278 continue
  280 continue
c
c
c
      call gauss(imat)
c
c
c
      do 290 n = 1,icvy
      cy(n) = cx(n)
  290 continue
      do 300 l = 1,lmax
      pst(l) = psp(l) + pse(l)
      do 295 n = 1,icvy
      pst(l) = pst(l) + cy(n)*psea(l,n)
  295 continue
      j=jl(l)
      k=kl(l)
      psi(j,k) = pst(l)
  300 continue
c
c
      go to 3000
 2000 print 2020 
 2020 format(51h too many exact boundary points or to many x-points)
      stop
 3000 return
      end
c
c
c
c
c
      subroutine start
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
      pi = 3.141592653589793
      frac = 0.
      if(ievolv.eq.1)  frac = float(nrun-1)/float(nruns-1)
      placur = placu1 + (placu2-placu1)*frac**placex
      ttfac  = ttfac1 + (ttfac2-ttfac1)*frac
      ell    = ell1   + (ell2-ell1)*frac
      emm    = emm1   + (emm2-emm1)*frac
      rmajo  = rmajo1 + (rmajo2-rmajo1)*frac
      zmajo  = zmajo1 + (zmajo2-zmajo1)*frac
      rmino  = rmino1 + (rmino2-rmino1)*frac
      cappa  = cappa1 + (cappa2-cappa1)*frac
      delta  = delta1 + (delta2-delta1)*frac
      hlamd  = hlamd1 + (hlamd2-hlamd1)*frac
      do 40 i=1,ilie
      rlim(i) = rlim1(i) + (rlim2(i)-rlim1(i))*frac
      zlim(i) = zlim1(i) + (zlim2(i)-zlim1(i))*frac
   40 continue
      do 42 n=1,ilia
      rlia(n) = rlia1(n) + (rlia2(n)-rlia1(n))*frac
      zlia(n) = zlia1(n) + (zlia2(n)-zlia1(n))*frac
   42 continue
      if(iansha.eq.1) then
      do 43 n=1,ilia
      win     = (pi*2.*float(n-1))/float(ilia)
      rlia(n) = rmajo+rmino*cos(win+delta*sin(win)-hlamd*sin(win*2.))
      zlia(n) = zmajo + rmino*cappa*sin(win)
   43 continue
      endif
      if ((ilie.gt.4).or.(ilia.gt.4)) then
         eftlim =  100.
         ritlim = -100.
         botlim =  100.
         toplim = -100.
         if(isaddl.eq.1) then
         do 44 i=1,ilie
         if(rlim(i).lt.eftlim) eftlim=rlim(i)
         if(rlim(i).gt.ritlim) ritlim=rlim(i)
         if(zlim(i).lt.botlim) botlim=zlim(i)
         if(zlim(i).gt.toplim) toplim=zlim(i)
   44    continue
         do 46 i=1,ilia
         if(rlia(i).lt.eftlim) eftlim=rlia(i)
         if(rlia(i).gt.ritlim) ritlim=rlia(i)
         if(zlia(i).lt.botlim) botlim=zlia(i)
         if(zlia(i).gt.toplim) toplim=zlia(i)
   46    continue
         endif
         if(isaddl.eq.0) then
         do 47 i=1,ilia
         if(zlia(i).gt.0.) then
            if(rlia(i).lt.eftlim) eftlim=rlia(i)
            if(rlia(i).gt.ritlim) ritlim=rlia(i)
            if(zlia(i).lt.botlim) botlim=zlia(i)
            if(zlia(i).gt.toplim) toplim=zlia(i)
            endif
   47    continue
         endif
         endif
c
      if(nz.eq.2.or.nz.eq.4.or.nz.eq.8.or.nz.eq.16)       go to 5
      if(nz.eq.32.or.nz.eq.64.or.nz.eq.128.or.nz.eq.256)  go to 5
      nz = 32
    5 do 8 i =1,icur
      rvfsq(i) = rvf(i)*rvf(i)
    8 continue
c
c
      icvy = icvx + 1
      ilim = ilie + ilia
      icon = ilie + ibro + ibzo + ivsec
      icvz = icvy + 1
      lmax = 2*nr+4*nz
      lilia= lmax+ilia
      ilro = ilie+ibro
      lee  = lmax+ilim+1
c
c set up mesh
c
      um = (4.*pi)/10000.
      dr = (ro-ri)/float(nr)
      jmax = nr+1
      do 10 j = 1,jmax
      r(j) = ri+(float(j-1)*dr)
   10 rjsq(j) = r(j)*r(j)
      dz = zu/float(nz)
      kmax = (2*nz)+1
      km = kmax-1
      do 20 k = 1,kmax
   20 z(k) = float(k-1-nz)*dz
      drh = 0.5*dr
      dzsq = dz*dz
      drsq = dr*dr
      cns = 1./dzsq
      do 22 j = 1,jmax
      ce(j) =(r(j)/(drsq*(r(j)+drh)))/cns
      cw(j) =(r(j)/(drsq*(r(j)-drh)))/cns
   22 continue
      aaa = (ritlim-eftlim)*0.5
      rze = (ritlim+eftlim)*0.5
      bzeru = (bzero*0.88)/rze
      zze = (toplim+botlim)*0.5
      bbb = (toplim-botlim)*0.5
      do 32 j = 1,jmax
      if (r(j)-rze)  32,34,34
   32 continue
   34 jze = j
      do 160 l = 1,lmax
      if(l-kmax)  110,110,120
  110 j = 1
      k = l
      go to 150
  120 if (l-kmax-nr)  130,130,140
  130 j = l-kmax+1
      k = kmax
      go to 150
  140 if (l-kmax-nr-km)  142,142,144
  142 j = jmax
      k = nr+2*kmax-l
      go to 150
  144 j = 4*nz+2*nr+2-l
      k = 1
  150 rl(l) = r(j)
      zl(l) = z(k)
      jl(l)=j
      kl(l)=k
  160 continue
      if(ilia)  163,163,161
  161 do 162 l = 1,ilia
      rl(lmax+l) = rlia(l)
  162 zl(lmax+l) = zlia(l)
  163 if(ilie)  169,169,164
  164 do 165 l = 1,ilie
      rl(lmax+ilia+l) = rlim(l)
  165 zl(lmax+ilia+l) = zlim(l)
  169 rl(lee) = rze
      zl(lee) = zze
      do 170 l=1,lee
  170 rlsq(l) = rl(l)*rl(l)
c
      ii = 0
      do 172 i=1,icur
      inside(i) = 0
      if((rvf(i).gt.ri).and.(rvf(i).lt.ro).and.
     1   (zvf(i).gt.-zu).and.(zvf(i).lt.zu)) then
         inside(i) = 1
         ii = ii+1
         endif
  172 continue
      ninsid = ii
c
c
c
c
c
c  compute elliptic integrals
c
      if(nrun.gt.1)  go to 199
      dae = pi/float(ncosei)
      do 178 n = 1,ncosei
      coei(n) = cos((0.5+float(n-1))*dae)
  178 continue
      daks = 1./float(mei-1)
      do 180 m = 1,mei
      aks(m) = daks*float(m-1)
  180 continue
      do 184 m = 1,mei
      ei(m) = 0.
      do 182 n = 1,ncosei
      ei(m) = ei(m) + (coei(n)/sqrt(1.-(aks(m)*coei(n))))
  182 continue
      ei(m) = ei(m)*dae
  184 continue
c
c
      do 198 m = 1,mei
      eig(m) = 0.
      eih(m) = 0.
      do 197 n = 1,ncosei
      pren = 1.-(aks(m)*coei(n))
      deno = pren*sqrt(pren)
      eig(m) = eig(m) + coei(n)/deno
      eih(m) = eih(m) + 1./deno
  197 continue
      eig(m) = eig(m)*dae
      eih(m) = eih(m)*dae
  198 continue
  199 continue
c
c
c  guess plasma current density
c
      plcu = 0.
      do 208 j = 1,jmax
      do 206 k = 1,kmax
      psn(j,k) = 0.
      if(j.eq.1.or.k.eq.1.or.j.eq.jmax.or.k.eq.kmax) go to 202
      rsq = (r(j)-rze)**2 + (((z(k)-zze)*aaa)/bbb)**2
      xsq = rsq/(aaa*aaa)
      if (xsq-1.)  200,200,202
  200 dsh(j,k) = (1.-xsq)**capaj
      plcu = plcu + dsh(j,k)
      go to 206
  202 dsh(j,k) = 0.
  206 continue
  208 continue
      ratio = placur/(plcu*dr*dz)
      do 220 j = 2,nr
      do 218 k = 2,km
      if(isaddl.eq.1) dsh(j,k) = dsh(j,k)*ratio
      if(isaddl.eq.0) dsh(j,k) = dsh(j,k)*ratio*0.5
  218 continue
  220 continue
      if(isaddl.eq.0) then
      kzero=kmax/2
      do 230 j = 2,nr
      do 228 k = 2,kzero
      dsh(j,k) = dsh(j,km-k+2)
  228 continue
  230 continue
      endif
c
      ellp1 = ell+1.
      ellp2 = ell+2.
      emmp1 = emm+1.
      emmp2 = emm+2.
      emmp3 = emm+3.
      oell  =  1./ell
      oellp1 = oell+1.
      oellp2 = oell+2.
      oemm   = 1./emm
      oemmp1 = oemm+1.
      oemmp2 = oemm+2.
c
c
c
      if(ilia)  301,301,298
  298 do 300 n = 1,ilia
  300 weit(n) = 1.+weitam*
     1 ((((zlia(n)-zze)*(zlia(n)-zze))/(bbb*bbb))**weitex)
  301 do 303 k = 1,16
      ifi = isvf(k)
      ila = isvl(k)
      do 302 i=ifi,ila
      diss(i) = dissi
      if(icoilon(k).eq.0) diss(i) = diss(i)*1.0e6
  302 continue
  303 continue
c
c
c
      do 306 i = 1,icur
      strk(i) = strki
      if(i-1) 306,306,305
  305 cstrk(i) = strk(i)/(((rvf(i)-rvf(i-1))**2.)
     1                   +((zvf(i)-zvf(i-1))**2.))
  306 continue
c
c
c
      do 350 j = 1,jmax
      dzlog = alog((16.*r(j))/dz)
      dz3lg = alog((16.*r(j))/(3.*dz))
      drlog = alog((16.*r(j))/dr)
      dr3lg = alog((16.*r(j))/(3.*dr))
      dzo8r = (dz/(8.*r(j)))**2.
      dro8r = (dr/(8.*r(j)))**2.
      acprn = dzlog-1.+(dzo8r*(dzlog-(2./3.)))
      asprn = dz3lg-1.+(9.*dzo8r*(dz3lg-(2./3.)))
      bcprn = drlog-1.+(dro8r*(drlog-(2./3.)))
      bsprn = dr3lg-1.+(9.*dro8r*(dr3lg-(2./3.)))
      acpsi(j) = acprn
      aspsi(j) = (3.*asprn-acprn)*0.5
      bcpsi(j) = bcprn
      bspsi    = (3.*bsprn-bcprn)*0.5
      bsips(j) = bspsi*(1.-(dr/(2.*r(j))))
      bsops(j) = bspsi*(1.+(dr/(2.*r(j))))
  350 continue
c
c
      two(km/2) = 0.
      lo = km/2
  361 l  = lo/2
      two(l) = sqrt(2.+two(lo))
      lo = l
  362 two(km-l) = -two(l)
      l = l+2*lo
      if(((2*l)/km)*(2*lo-3))  364,363,361
  363 two(l) = (two(l+lo)+two(l-lo))/two(lo)
      go to 362
  364 continue
      return
      end
c
c
c
c
c
      subroutine gauss(mat)
      include 'fbtmgams.inc'
      include 'fbt.inc'
      ni=mat
      na=ni+1
      do 234 mm = 2,ni
      do 232 n = mm,na
  232 ax(mm-1,n) = ax(mm-1,n)/ax(mm-1,mm-1)
      do 233 m = mm,ni
      do 233 n = mm,na
  233 ax(m,n) = ax(m,n)-ax(mm-1,n)*ax(m,mm-1)
  234 continue
      ax(ni,na) = ax(ni,na)/ax(ni,ni)
      do 236 m = 1,ni
      mm = ni-m+1
      ml = mm-1
      cx(mm) = ax(mm,na)
      if(mm.eq.1)  go to 236
      do 235 mb = 1,ml
  235 ax(mb,na) = ax(mb,na) - cx(mm)*ax(mb,mm)
  236 continue
      return
      end
c
c
c
c
c
      subroutine psiext
      include 'fbtmgams.inc'
      include 'fbt.inc'
      par = float(mei-1)
      do 40 l = 1,lee
      pse(l) = 0.
      do 7 n = 1,icvx
    7 psea(l,n) = 0.
      b = rlsq(l)
      do 30 i = 1,icur
      if((inside(i).eq.1).and.(l.le.lmax)) go to 30
      a = zl(l)-zvf(i)
      c = rvfsq(i)
      d = rl(l)*rvf(i)*2.
c cvf no more used      e = cvf(i)*rvf(i)
      do 8 n = 1,icvx
    8 ux(n) = cvx(i,n)*rvf(i)
      h = (a*a)+b+c
      aksq = d/h
      m = (aksq*par) + 1.
      if(m.eq.mei)  m=m-1
      part = (ei(m+1)-ei(m))/daks
      eii = (ei(m)+((aksq-aks(m))*part))/sqrt(h)
c cvf no more used      pse(l) = pse(l) + (e*eii)
      do 9 n = 1,icvx
    9 psea(l,n) = psea(l,n) + (ux(n)*eii)
   30 continue
      pse(l) = pse(l)*rl(l)*0.0002
      do 31 n = 1,icvx
   31 psea(l,n) = psea(l,n)*rl(l)*0.0002
      psea(l,icvy) = 1.
   40 continue
      psea(lee,icvy) = 0.
c
      ii = 1
      do 50 i=1,icur
      if(inside(i).eq.0) go to 50
      c = rvfsq(i)
      do 48 j=1,jmax
      b = rjsq(j)
      d = r(j)*rvf(i)*2.
      do 46 k=1,kmax
      a = z(k)-zvf(i)
      h = (a*a)+b+c
      aksq = d/h
      m = (aksq*par)+1.
      if(m.eq.mei) m=m-1
      part = (ei(m+1)-ei(m))/daks
      eii  = (ei(m)+((aksq-aks(m))*part))/sqrt(h)
      psint(ii,j,k) = rvf(i)*eii*r(j)*0.0002
   46 continue
   48 continue
      ii = ii + 1
   50 continue
c
c
      if(ibro)  72,72,56
   56 do 70 l = 1,ibro
      bere(l) = 0.
      do 57 n = 1,icvx
   57 berx(l,n) = 0.
      b = rbro(l)*rbro(l)
      do 60 i = 1,icur
      a = zbro(l) - zvf(i)
      c = rvfsq(i)
      d = rbro(l)*rvf(i)*2.
c cvf no more used      e = cvf(i)*rvf(i)
      do 58 n = 1,icvx
   58 ux(n) = cvx(i,n)*rvf(i)
      h = (a*a) + b + c
      hh = h*sqrt(h)
      aksq = d/h
      m = (aksq*par)+1.
      if(m.eq.mei)  m=m-1
      partg = (eig(m+1)-eig(m))*par
      eiig = eig(m) + ((aksq-aks(m))*partg)
      eiii = (a*eiig*0.0002)/hh
c cvf no more used      bere(l) = bere(l) + (e*eiii)
      do 59 n = 1,icvx
   59 berx(l,n) = berx(l,n)+(ux(n)*eiii)
   60 continue
      berx(l,icvy) = 0.
   70 continue
   72 if(ibzo)  100,100,73
   73 do 90 l = 1,ibzo
      beze(l) = 0.
      do 77 n = 1,icvx
   77 bezx(l,n) = 0.
      b = rbzo(l)*rbzo(l)
      do 80 i = 1,icur
      a = zbzo(l)-zvf(i)
      c = rvfsq(i)
      d = rbzo(l)*rvf(i)*2.
c cvf no more used      e = cvf(i)*rvf(i)
      do 78 n = 1,icvx
   78 ux(n) = cvx(i,n)*rvf(i)
      h = (a*a)+b+c
      hh = h*sqrt(h)
      aksq = d/h
      m = (aksq*par)+1.
      if(m.eq.mei)  m=m-1
      partg = (eig(m+1)-eig(m))*par
      parth = (eih(m+1)-eih(m))*par
      eiig  = eig(m)+((aksq-aks(m))*partg)
      eiih  = eih(m)+((aksq-aks(m))*parth)
      eiii = ((rvf(i)*eiih-rbzo(l)*eiig)*0.0002)/hh
c cvf no more used      beze(l) = beze(l) + (e*eiii)
      do 79 n = 1,icvx
      bezx(l,n) = bezx(l,n)+(ux(n)*eiii)
   79 continue
   80 continue
      bezx(l,icvy) = 0.
   90 continue
  100 return
      end
c
c
c
c
c
      subroutine psipla
      include 'fbtmgams.inc'
      include 'fbt.inc'
      db = dr*dz
      par = float(mei-1)
      do 300 l = 1,lmax
      j = jl(l)
      k = kl(l)
      psi(j,k) = 0.
  300 continue
c
c
c
      call psinew
c
c
c
      um1 = 2.*um*dr*dr
      um2 = 2.*um*dz*dz
      do 312 l = 1,lmax
      j = jl(l)
      k = kl(l)
      if(l-kmax)  301,301,302
  301 rss(l) = (4.*psi(j+1,k)-psi(j+2,k))/(um1*r(j))
      go to 312
  302 if(l-kmax-nr)  304,304,306
  304 rss(l) = (4.*psi(j,k-1)-psi(j,k-2))/(um2*r(j))
      go to 312
  306 if(l-kmax-nr-km)  308,308,310
  308 rss(l) = (4.*psi(j-1,k)-psi(j-2,k))/(um1*r(j))
      go to 312
  310 rss(l) = (4.*psi(j,k+1)-psi(j,k+2))/(um2*r(j))
  312 continue
      do 404 l = 1,lmax
      psp(l) = 0.
      b = rlsq(l)
      j = jl(l)
      k = kl(l)
      if((j.eq.1).or.(j.eq.jmax))    iloc = 2
      if(k.eq.kmax)                  iloc = 1
      if(k.eq.1)                     iloc =-1
      if((l.eq.1).or.(l.eq.kmax).or.(l.eq.kmax+nr).
     1 or.(l.eq.kmax+nr+km))         iloc = 4
      do 402 ls = 1,lmax
      if(iloc-4)  324,320,324
  320 if(ls-l)    400,402,400
  324 if(iloc-2)  340,326,340
  326 if(ls-l)    330,328,330
  328 psp(l) = psp(l) + (acpsi(j)*rss(ls))
      go to 402
  330 if(ls-l+1)  332,334,332
  332 if(ls-l-1)  400,334,400
  334 psp(l) = psp(l) + (aspsi(j)*rss(ls))
      go to 402
  340 if(ls-l)    344,342,344
  342 psp(l) = psp(l) + (bcpsi(j)*rss(ls))
      go to 402
  344 if(ls-l+iloc)  348,346,348
  346 psp(l) = psp(l) + (bsips(j)*rss(ls))
      go to 402
  348 if(ls-l-iloc)  400,350,400
  350 psp(l) = psp(l) + (bsops(j)*rss(ls))
      go to 402
  400 c = rlsq(ls)
      d = rl(l)*rl(ls)*2.
      a = zl(l)-zl(ls)
      h = a*a+b+c
      e = rss(ls)*rl(ls)
      aksq = d/h
      m = (aksq*par) + 1.
      if(m.eq.mei) m=m-1
      part = (ei(m+1)-ei(m))/daks
      eii = ei(m) + ((aksq-aks(m))*part)
      psp(l) = psp(l) + ((e*eii)/sqrt(h))
  402 continue
      psp(l) = psp(l)*rl(l)*db*0.0002
  404 continue
c
c
c
      lst = lmax+1
      len=lee
      do 550 l = lst,len
      psp(l) = 0.
      b = rlsq(l)
      do 540 j = 2,nr
      c = rjsq(j)
      d = rl(l)*r(j)*2.
      do 530 k = 2,km
      if (dsh(j,k))  530,530,510
  510 a = zl(l)-z(k)
      h = (a*a)+b+c
      e = dsh(j,k)*r(j)
      aksq = d/h
      m = (aksq*par) + 1.
      if(m.eq.mei)  m=m-1
      part = (ei(m+1)-ei(m))/daks
      eii = ei(m) + ((aksq-aks(m))*part)
      psp(l) = psp(l) + ((e*eii)/sqrt(h))
  530 continue
  540 continue
      psp(l) = psp(l)*rl(l)*db*0.0002
  550 continue
c
c
c
      if(ibro)  570,570,560
  560 do 568 l = 1,ibro
      berp(l) = 0.
      b = rbro(l)*rbro(l)
      do 566 j = 2,nr
      c = rjsq(j)
      d = rbro(l)*r(j)*2.
      do 564 k = 2,km
      if(dsh(j,k))  564,564,562
  562 a = zbro(l)-z(k)
      e = dsh(j,k)*r(j)
      h = (a*a)+b+c
      hh = h*sqrt(h)
      aksq = d/h
      m = (aksq*par)+1.
      if(m.eq.mei)  m=m-1
      partg = (eig(m+1)-eig(m))*par
      eiig = eig(m)+((aksq-aks(m))*partg)
      berp(l) = berp(l) + (e*a*eiig)/hh
  564 continue
  566 continue
      berp(l) = berp(l)*0.0002*db
  568 continue
c
c
c
  570 if(ibzo)  700,700,572
  572 do 580 l = 1,ibzo
      bezp(l) = 0.
      b = rbzo(l)*rbzo(l)
      do 578 j = 2,nr
      c = rjsq(j)
      d = rbzo(l)*r(j)*2.
      do 576 k = 2,km
      if(dsh(j,k))  576,576,574
  574 a = zbzo(l)-z(k)
      e = dsh(j,k)*r(j)
      h = (a*a)+b+c
      hh = h*sqrt(h)
      aksq = d/h
      m = (aksq*par)+1.
      if(m.eq.mei)  m=m-1
      partg = (eig(m+1)-eig(m))*par
      parth = (eih(m+1)-eih(m))*par
      eiig  = eig(m)+((aksq-aks(m))*partg)
      eiih  = eih(m)+((aksq-aks(m))*parth)
      bezp(l) = bezp(l)+((e/hh)*(r(j)*eiih-rbzo(l)*eiig))
  576 continue
  578 continue
      bezp(l) = bezp(l)*0.0002*db
  580 continue
  700 return
      end
c
c
c
c
c
      subroutine psinew
      include 'fbtmgams.inc'
      include 'fbt.inc'
      dimension p(mamp),am(mamp),az(mamp),
     1 ap(mamp),cz(mamp),cp(mamp),d(mamp)
      do 3 k = 1,kmax
      do 2 j = 1,jmax
      i = k*jmax-j+1
      if(j.eq.1.or.k.eq.1.or.j.eq.jmax.or.k.eq.kmax) go to 1
      gg(i) = (um*r(j)*dsh(j,k))/(2.*cns)
      go to 2
    1 gg(i) = psi(j,k)
    2 continue
    3 continue
      n = km
      m = nr
c
c
c
      lo = n/2
      ju = (n-1)*(m+1)
      iu = m-1
      i = m+2
    5 gg(i+1) = gg(i+1) + (gg(i)*ce(nr))*0.5
      i = i+m+1
      gg(i-2) = gg(i-2) + (gg(i-1)*cw(2))*0.5
      if(i-1.le.ju)  go to 5
      mode = 2
   15 li = 2*lo
      iphase = (2*mode)-(li/n)
      jd = (m+1)*(n/li)
      jh = (m+1)*(n/(2*li))
      jt = jd+jh
      ji = 2*jd
      jo = jd*mode
      do 88 j = jo,ju,ji
      j1 = j+2
      jiu = j+iu+1
      go to (50,44,36,28),iphase
   28 do 29 i = j1,jiu
      pj = gg(i)-gg(i+jt)-gg(i-jt)
      gg(i) = gg(i)-gg(i+jh)-gg(i-jh)+gg(i+jd)+gg(i-jd)
   29 p(i-j-1) = pj+gg(i)
      go to 60
   36 do 37 i = j1,jiu
      p(i-j-1) = 2.*gg(i)
   37 gg(i) = gg(i+jd)+gg(i-jd)
      go to 60
   44 do 45 i = j1,jiu
      p(i-j-1) = 2.*gg(i)+gg(i+jd)+gg(i-jd)
   45 gg(i) = gg(i)-gg(i+jh)-gg(i-jh)
      go to 60
   50 do 53 i = j1,jiu
      p(i-j-1) = 2.*gg(i)+gg(i+jd)+gg(i-jd)
   53 gg(i) = 0.
   60 do 72 l = lo,n,li
c
      do 64 i = 1,iu
      am(i) = -ce(nr-i+1)
      ap(i) = -cw(nr-i+1)
      az(i) = 2.-am(i)-ap(i)-two(l)
   64 continue
      am(1) = 0.
      ap(iu) = 0.
      cz(1) = az(1)
      cp(1) = ap(1)
      d(1)  = p(1)
      do 65 i = 2,iu
      rat = am(i)/cz(i-1)
      cz(i) = az(i)-rat*cp(i-1)
      cp(i) = ap(i)
      d(i)  = p(i)-rat*d(i-1)
   65 continue
      p(iu) = d(iu)/cz(iu)
      do 66 i2 = 2,iu
      i = iu+1-i2
      p(i) = (d(i)-cp(i)*p(i+1))/cz(i)
   66 continue
c
   72 continue
      do 77 i = j1,jiu
   77 gg(i) = gg(i)+p(i-j-1)
   88 continue
      go to (94,93,92,92),iphase
   92 lo = lo/2
      if(lo.eq.1) mode=1
      go to 15
   93 lo = 2*lo
      if(lo.lt.n)  go to 15
   94 do 98 k = 2,km
      do 97 j = 2,nr
      i = k*jmax-j+1
      psi(j,k) = gg(i)
   97 continue
   98 continue
      return
      end
c
c
c
c
c
      subroutine psiins
      include 'fbtmgams.inc'
      include 'fbt.inc'
      ii = 1
      do 100 i=1,icur
      if(inside(i).eq.0) go to 100
c cvf no more used      curr(i) = cvf(i)
      curr(i) = 0.
      do 20 n=1,icvx
      curr(i) = curr(i) + cvx(i,n)*cy(n)
   20 continue
      do 90 j=1,jmax
      do 80 k=1,kmax
      psi(j,k) = psi(j,k) + psint(ii,j,k)*curr(i)
   80 continue 
   90 continue
      ii = ii + 1
  100 continue
      return
      end
c
c
c
c
c
      subroutine curden
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
c                  find the axis
c
      iax = 1
      do 10 j = 2,nr
      do 8 k = 2,km
      do 4 i = 1,icur
      if(inside(i).eq.0) go to 4
      if((abs(rvf(i)-r(j)).lt.(2.*dr)).and.
     1   (abs(zvf(i)-z(k)).lt.(2.*dz))) go to 8
    4 continue
      a = psi(j,k)
      p2 = psi(j+1,k)
      p3 = psi(j,k+1)
      p4 = psi(j-1,k)
      p5 = psi(j,k-1)
      if(p2.gt.a.or.p3.gt.a.or.p4.gt.a.or.p5.gt.a) go to 8
      b = (p2-p4)/(2.*dr)
      c = (p3-p5)/(2.*dz)
      d = (p2+p4-2.*a)/(2.*dr*dr)
      e = (p3+p5-2.*a)/(2.*dz*dz)
      if(iax.eq.1) then
         jc = j
         kc = k
         rax = r(jc)-(b/(2.*d))
         zax = z(kc)-(c/(2.*e))
         psax = a-((b*b)/(4.*d))-((c*c)/(4.*e))
         ellipt = sqrt(d/e)
         iax = iax + 1
         go to 8
      endif
      if(iax.eq.2) then
         jc2 = j
         kc2 = k
         rax2 = r(jc2)-(b/(2.*d))
         zax2 = z(kc2)-(c/(2.*e))
         psax2 = a-((b*b)/(4.*d))-((c*c)/(4.*e))
         ellipt2 = sqrt(d/e)
         iax = iax + 1
         go to 12
      endif
    8 continue
   10 continue
   12 if((iax.eq.3.and.psax2.gt.psax).and.(isaddl.eq.1)) then
         jc = jc2  
         kc = kc2
         rax = rax2
         zax = zax2
         psax = psax2
         ellipt = ellipt2
      endif
      if((zax.lt.zax2).and.(isaddl.eq.0)) then
      jcbid  = jc
      kcbid  = kc
      zaxbid = zax
      raxbid = rax
      psaxbid= psax
      elliptb= ellipt
      jc     = jc2
      kc     = kc2
      zax    = zax2
      rax    = rax2
      psax   = psax2
      ellipt = ellipt2
      jc2    = jcbid
      kc2    = kcbid
      zax2   = zaxbid
      rax2   = raxbid
      psax2  = psaxbid
      ellipt2= elliptb
      endif
c
c
c                  search for highest saddle point
c
      psione = -psax
      isone  = 0
      do 90 is=1,mcoi
      psisa(is) = 0.
      jsad(is)  = 0
      ksad(is)  = 0
      rxp(is)   = 0.
      zxp(is)   = 0.
   90 continue
      is        = 1
      do 330 k = 1,km
      do 320 j = 1,nr
      ps4 = psi(j,k)
      ps1 = psi(j+1,k)
      ps2 = psi(j+1,k+1)
      ps3 = psi(j,k+1)
      if((ps1.ge.ps4).and.(ps2.lt.ps1).and.
     1   (ps3.gt.ps2).and.(ps4.le.ps3)) go to 100
      if((ps1.le.ps4).and.(ps2.gt.ps1).and.
     1   (ps3.lt.ps2).and.(ps4.ge.ps3)) go to 100
      go to 200
  100 xxx = (dr*(ps4-ps3))/(ps2-ps3-ps1+ps4)
      yyy = (dz*(ps1-ps4))/(ps3-ps4-ps2+ps1)
      psix= (ps2*ps4-ps3*ps1)/(ps2-ps3-ps1+ps4)
      go to 310
  200 if((j.eq.1).or.(k.eq.1)) go to 320
      ps5 = psi(j-1,k)
      ps6 = psi(j,k-1)
      if((ps1.gt.ps4).and.(ps3.lt.ps4).and.
     1   (ps5.gt.ps4).and.(ps6.lt.ps4)) go to 210
      if((ps1.lt.ps4).and.(ps3.gt.ps4).and.
     1   (ps5.lt.ps4).and.(ps6.gt.ps4)) go to 210
      go to 320
  210 xxx = (dr*(ps5-ps1))/(2.*(ps1-2.*ps4+ps5))
      yyy = (dz*(ps6-ps3))/(2.*(ps3-2.*ps4+ps6))
      psix= ps4
     1 -((ps1-ps5)*(ps1-ps5))/(16.*(ps1-2.*ps4+ps5))
     2 -((ps3-ps6)*(ps3-ps6))/(16.*(ps3-2.*ps4+ps6))
  310 if(is.gt.mcoi) go to 900
      jsad(is)  = j
      ksad(is)  = k
      rxp(is)   = r(j)+xxx
      zxp(is)   = z(k)+yyy
      psisa(is) = psix
      if(psix.gt.psione) then
         psione = psix
         isone  = is
         endif
      is = is + 1
  320 continue
  330 continue
      ismax = is-1
      psitwo = -psax
      istwo  = 0
      do 332 is=1,ismax
      if(is.eq.isone) go to 332
      if(psisa(is).gt.psitwo) then
         psitwo = psisa(is)
         istwo  = is
         endif
  332 continue
      psisad = psione
      if(istwo.ne.0) then
         rxx = rxp(isone)
         zxx = zxp(isone)
         dshsum = 0.
         do 336 j=1,jmax
         do 334 k=1,kmax
         if((abs(z(k)-zxx).lt.(0.15*bbb)).and.
     1      (abs(r(j)-rxx).lt.(0.15*aaa)))
     2      dshsum = dshsum + abs(dsh(j,k))
  334    continue
  336    continue
         if(dshsum.lt.1.0e-10) then
            psisad = psitwo
            isonee = isone
            isone  = istwo
            istwo  = isonee
            endif
         endif
      psilim = psisad+(1.-psirat)*(psax-psisad)
      if((psilim.gt.0.).and.(isaddl.eq.1)) then
         do 342 j=1,jmax
         do 340 k=1,kmax
         psi(j,k) = psi(j,k)-psilim
  340    continue
  342    continue
         psax   = psax-psilim
         psione = psione-psilim
         psitwo = psitwo-psilim
         do 350 is=1,ismax
         psisa(is) = psisa(is)-psilim
  350    continue
         endif
      iacx = 0
      if((abs(psione).lt.(0.01*psax)).or.
     1   (abs(psitwo).lt.(0.01*psax))) iacx=1
      if((abs(psione).lt.(0.01*psax)).and.
     1   (abs(psitwo).lt.(0.01*psax))) iacx=2
c
c                  eliminate current outside separatrix
c
c
      do 405 j=1,jmax
      do 404 k=1,kmax
      dsh(j,k) = 0.
      if(psi(j,k).gt.0.) dsh(j,k)=1.
  404 continue
  405 continue
      if((isaddl.eq.1).and.(neqtcv.ne.888)) then
      do 425 is=1,ismax
      if(psisa(is).lt.(-0.1*psax))  go to 425
      rxx = rxp(is)
      zxx = zxp(is)
      deltur = rxx-rax
      deltuz = zxx-zax
      if(deltuz)  410,410,412
  410 isec = 3
      if(deltur.le.deltuz)  isec=4
      if(deltur.ge.-deltuz) isec=2
      go to 414
  412 isec = 1
      if(deltur.le.-deltuz) isec=4
      if(deltur.ge.deltuz)  isec=2
  414 if(isec.eq.1.or.isec.eq.3) then
         bsep = -(deltur/deltuz)
         asep = zxx - rxx*bsep
         do 420 j=1,jmax
         zcut = asep + r(j)*bsep
         do 418 k=1,kmax
         if((isec.eq.1).and.(z(k).gt.zcut)) dsh(j,k)=0.
         if((isec.eq.3).and.(z(k).lt.zcut)) dsh(j,k)=0.
  418    continue
  420    continue
      endif
      if(isec.eq.2.or.isec.eq.4) then
         dsep = -(deltuz/deltur)
         csep = rxx - zxx*dsep
         do 424 k=1,kmax
         rcut = csep + z(k)*dsep
         do 422 j=1,jmax
         if((isec.eq.2).and.(r(j).gt.rcut)) dsh(j,k)=0.
         if((isec.eq.4).and.(r(j).lt.rcut)) dsh(j,k)=0.
  422    continue
  424    continue
      endif
  425    continue
      endif
c
c
      do 520 k=1,kmax
      do 510 j=1,jmax
      if((psi(j,k).le.0.).or.(dsh(j,k).le.0.))  go to 512
  510 dsh(j,k) = 0.
  512 do 514 j1=1,jmax
      j = jmax+1-j1
      if((psi(j,k).le.0.).or.(dsh(j,k).le.0.))  go to 520
  514 dsh(j,k) = 0.
  520 continue
      do 540 j=1,jmax
      do 530 k=1,kmax
      if((psi(j,k).le.0.).or.(dsh(j,k).le.0.))  go to 532
  530 dsh(j,k) = 0.
  532 do 534 k1=1,kmax
      k = kmax+1-k1
      if((psi(j,k).le.0.).or.(dsh(j,k).le.0.))  go to 540
  534 dsh(j,k) = 0.
  540 continue
c
c
c                  compute plasma current
c
      ineg = 0
      plcu1 = 0.
      plcu2 = 0.
      do 620 j = 2,nr
      do 610 k = 2,km
      if(dsh(j,k))  610,610,600
  600 if (psi(j,k))  608,608,602
  602 phi = psi(j,k)/psax
      xxx = 1.-phi
      if(xxx.le.0.)  xxx=0.
      if(ell.gt.1.) then
         phell  = phi**ell
         phell1 = phell*phi
         ppri   = ppal*phell+(1.-ppal)*(phell-phell1)
         endif
      if(ell.le.1.) then
         xoell  = xxx**oell
         xoell1 = xoell*xxx
         ppri   = ppal*(1.-xoell)+(1.-ppal)*(xoell-xoell1)
         endif
      if(emm.ge.1.) then
         phemm  = phi**emm
         phemm1 = phemm*phi
         phemm2 = phemm1*phi
         ttpr1  = phemm
         ttpr2  = phemm1-phemm+omega*(phemm2-phemm1)
         endif
      if(emm.lt.1.) then
         xoemm  = xxx**oemm
         xoemm1 = xoemm*xxx
         ttpr1  = 1.-xoemm
         ttpr2  = xoemm1-xoemm
         endif
      dsha1 = r(j)*ppri*ppfac+((rze*rze)/r(j))*ttfac*ttpr1
      dsha2 = ((rze*rze)/r(j))*ttfac*ttpr2
      plcu1 = plcu1 + dsha1
      plcu2 = plcu2 + dsha2
      go to 610
  608 dsh(j,k) = 0.
  610 continue
  620 continue
      plcu1 = plcu1*dr*dz
      plcu2 = plcu2*dr*dz
      aa33  = -1./emmp1
      if(emm.ge.1.) aa44=1./emmp2-1./emmp1
     1              +omega*(1./emmp3-1./emmp2)
      if(emm.lt.1.) aa44=1./oemmp2-1./oemmp1
      aa55 = -((qzero*um*rax*rax*sqrt(d*e))/(rze*bzeru*(d+e)))
     1       *(rax*ppfac*ppal + ((rze*rze)/rax)*ttfac)
      aa66 = (um*psax*ttfac)/(bzeru*bzeru)
      ttbe = (placur*(aa55+aa66*aa33) - plcu1)/
     1       (plcu2+placur*aa66*aa44)
      ratio = placur/(plcu1+ttbe*plcu2)
      do 630 j = 2,nr
      do 628 k = 2,km
      if(dsh(j,k))  628,628,622
  622 if(psi(j,k))  626,626,624
  624 phi = psi(j,k)/psax
      radi = r(j)
      dsh(j,k) = radi*pprime(phi) + ttprim(phi)/radi
      if(dsh(j,k).lt.0.) ineg=1
      go to 628
  626 dsh(j,k) = 0.
  628 continue
  630 continue
c
      if(neqtcv.eq.888) then
      plcutop=0.
      plcubot=0.
      do 730 j = 2,nr
      do 728 k = 2,km
      if(z(k).gt.0.) plcutop=plcutop+dsh(j,k)
      if(z(k).lt.0.) plcubot=plcubot+dsh(j,k)
  728 continue
  730 continue
      factop=(plcutop+plcubot)/(2.*plcutop)
      facbot=(plcutop+plcubot)/(2.*plcubot)
      do 830 j = 2,nr
      do 828 k = 2,km
      if(z(k).gt.0.) dsh(j,k)=dsh(j,k)*factop
      if(z(k).lt.0.) dsh(j,k)=dsh(j,k)*facbot
  828 continue
  830 continue
      endif
c
      if(ineg.eq.1) print 903
  903 format(32h negative plasma current density)
      go to 910
c
  900 print 901
  901 format(18h too many x-points)
      stop
c
c
c
c
  910 print 1000,jsad(1),ksad(1),jsad(2),ksad(2),rax,zax,psax,
     1 cy(2),cy(3)
 1000 format(1x,4i3,5e12.5)
 2000 return
      end
c
c
c
c
c
      subroutine primat
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
c
      pi = 3.141592653589793
      pmax = 0.
      psum = 0.
      summu = 0.
      sumflu = 0.
      sumfla = 0.
      do 142 j = 1,jmax
      do 140 k = 1,kmax
      if (dsh(j,k))  130,130,132
  130 psn(j,k) = 0.
      go to 140
c
  132 phi = psi(j,k)/psax
      psn(j,k) = pres(phi)
      tsq = tsqa(phi)
      bsqa   = ((rze*rze*bzeru*bzeru)+tsq)/(r(j)*r(j))
      btor   = sqrt(bsqa)
      bvtor  = (rze*bzeru)/r(j)
      xxx = tsq/(bzeru*bzeru*rze*rze)
      summu  = summu + tsq/r(j)
      sumfla = sumfla + (tsq/r(j))*(1.-xxx/4.+(xxx*xxx)/8.)
      sumflu = sumflu + (btor-bvtor)
c
      psum = psum + psn(j,k)*r(j)
  137 if (psn(j,k)-pmax)  140,140,138
  138 pmax = psn(j,k)
  140 continue
  142 continue
      psum =(psum*dr*dz*1000.)/rze
      betai = (20.*psum)/(placur*placur)
      summu  = (summu*dr*dz)/rze
      diamag = (10000.*summu)/(placur*placur*um)
      torflu = sumflu*dr*dz
      torfla = (sumfla*dr*dz)/(2.*rze*bzeru)
      area = 0.
      bsqsum = 0.
      volum=0.
      do 261 k = 1,km
      do 260 j = 1,nr
      t = psi(j+1,k)
      u = psi(j+1,k+1)
      v = psi(j,k+1)
      w = psi(j,k)
      bz = (t+u-v-w)/(dr*(r(j)+r(j+1)))
      br =-(u+v-t-w)/(dz*(r(j)+r(j+1)))
      if((dsh(j,k)+dsh(j+1,k)+dsh(j,k+1)+dsh(j+1,k+1)).le.1.e-10)goto210
      if(t.le.0..and.u.le.0..and.v.le.0..and.w.le.0.)   go to 210
      if(t.gt.0..and.u.gt.0..and.v.gt.0..and.w.gt.0.)   go to 212
      if(t.gt.0..and.u.gt.0..and.v.le.0..and.w.le.0.)   go to 214
      if(u.gt.0..and.v.gt.0..and.w.le.0..and.t.le.0.)   go to 216
      if(v.gt.0..and.w.gt.0..and.t.le.0..and.u.le.0.)   go to 218
      if(w.gt.0..and.t.gt.0..and.u.le.0..and.v.le.0.)   go to 220
      if(t.gt.0..and.u.le.0..and.v.le.0..and.w.le.0.)   go to 224
      if(u.gt.0..and.v.le.0..and.w.le.0..and.t.le.0.)   go to 226
      if(v.gt.0..and.w.le.0..and.t.le.0..and.u.le.0.)   go to 228
      if(w.gt.0..and.t.le.0..and.u.le.0..and.v.le.0.)   go to 230
      if(t.gt.0..and.u.gt.0..and.v.gt.0..and.w.le.0.)   go to 234
      if(u.gt.0..and.v.gt.0..and.w.gt.0..and.t.le.0.)   go to 236
      if(v.gt.0..and.w.gt.0..and.t.gt.0..and.u.le.0.)   go to 238
      if(w.gt.0..and.t.gt.0..and.u.gt.0..and.v.le.0.)   go to 240
  210 db = 0.
      go to 260
  212 db = 1.
      go to 250
  214 db = 0.5*((t/(t-w))+(u/(u-v)))
      go to 250
  216 db = 0.5*((u/(u-t))+(v/(v-w)))
      go to 250
  218 db = 0.5*((v/(v-u))+(w/(w-t)))
      go to 250
  220 db = 0.5*((w/(w-v))+(t/(t-u)))
      go to 250
  224 db = 0.5*(t/(t-u))*(t/(t-w))
      go to 250
  226 db = 0.5*(u/(u-v))*(u/(u-t))
      go to 250
  228 db = 0.5*(v/(v-w))*(v/(v-u))
      go to 250
  230 db = 0.5*(w/(w-t))*(w/(w-v))
      go to 250
  234 db = 1.-(0.5*(w/(w-t))*(w/(w-v)))
      go to 250
  236 db = 1.-(0.5*(t/(t-u))*(t/(t-w)))
      go to 250
  238 db = 1.-(0.5*(u/(u-v))*(u/(u-t)))
      go to 250
  240 db = 1.-(0.5*(v/(v-w))*(v/(v-u)))
  250 area = area+db
      volum = volum + (db*0.5*(r(j)+r(j+1)))
      bsqsum= bsqsum+((bz*bz)+(br*br))*(db*0.5*(r(j)+r(j+1)))
  260 continue
  261 continue
      bsqsum = (bsqsum*dr*dz)/rze
      eli = (10000.*bsqsum)/(placur*placur*um)
      volum = volum*dr*dz
      area = area*dr*dz
      betat = (rze*psum*um)/(500.*volum*bzeru*bzeru)
      qcyl = (2.*area*bzeru)/(um*placur*rze)
c
      dsm = 0.
      do 270 j = 1,jmax
      do 268 k = 1,kmax
      if (dsh(j,k)-dsm)  268,268,266
  266 dsm = dsh(j,k)
  268 continue
  270 continue
      print 298,neqtcv
  298 format (33h equilibrium number              ,i12)
      print 360,iacx
  360 format (33h number of active X-points       ,i12)
      if(iacx.ge.1) print 362,rxp(isone),zxp(isone)
  362 format (33h r,z of first active X-point     ,2e12.4)
      if(iacx.eq.2) print 364,rxp(istwo),zxp(istwo)
  364 format (33h r,z of second active X-point    ,2e12.4)
      print 400,placur
  400 format (33h plasma current                  ,e12.4)
      print 402,betai
  402 format (33h beta poloidal (betai)           ,e12.4)
      print 403,betat
  403 format (33h beta toroidal                   ,e12.4)
      print 404,qq(npsi+2)
  404 format (33h q-95                            ,e12.4)
      print 406,rax
  406 format (33h radius of magnetic axis         ,e12.4)
      print 408,psax
  408 format (33h psi on axis                     ,e12.4)
      print 410,dsm
  410 format (33h max current density (ka/m**2)   ,e12.4)
      pmaxm= pmax*1000.
      print 412,pmaxm
  412 format (33h max pressure (mks)              ,e12.4)
      print 414,eli
  414 format (33h internal inductance             ,e12.4)
      print 416,diamag
  416 format (33h plasma diamagnetism             ,e12.4)
      print 418, torflu
  418 format (33h plasma toroidal flux            ,e12.4)
      apromu =(8.*pi*bzeru*torflu)/
     1        (um*um*placur*placur)
      print 420,apromu
  420 format (33h approximate diamagnetism        ,e12.4)
      print 422,qzero
  422 format (33h q-zero                          ,e12.4)
      print 424,bzero
  424 format (33h Toroidal Vacuum Field at R=0.88 ,e12.4)
      print 426,torfla
  426 format (33h approximate toroidal flux       ,e12.4)
      print 428,ellipt
  428 format (33h central ellipticity             ,e12.4)
      print 429,cappab
  429 format (33h edge elongation                 ,e12.4)
      print 430,cappa95
  430 format (33h elongation of 95% flux surface  ,e12.4)
      print 431,deltab
  431 format (33h edge triangularity              ,e12.4)
      print 432,delta95
  432 format (33h triangularity of 95% surface    ,e12.4)
      cursum = 0.
      do 502 i = 1,icur
c cvf no more used      curr(i)=cvf(i)
      curr(i) = 0.
      do 500 n = 1,icvx
  500 curr(i) = curr(i) + cvx(i,n)*cy(n)
      cursum = cursum + abs(curr(i))
  502 continue
      if (iwrida.ne.1) go to 515
      print 510
  510 format(14h coil currents/36h radius      height      current(ka))
      do 514 i = 1,icur
      print 512,rvf(i),zvf(i),curr(i)
  512 format(3e12.4)
  514 continue
  515 print 516,cursum
  516 format(33h sum of absolute coil currents   ,e12.4)
      shvsec = 0.
      do 520 n = 1,icvx
  520 shvsec = shvsec + (psea(lee,n)*cy(n))
      shvsec = shvsec*2.*pi
      ohvsec = pse(lee)*2.*pi
      print 522,ohvsec
  522 format(33h oh volt seconds                 ,e12.4)
      print 524,shvsec
  524 format(33h shaping volt seconds            ,e12.4)
      return
      end
c
c
c
c
c
      subroutine flusurf
      include 'fbtmgams.inc'
      include 'fbt.inc'
      print 1
    1 format (61h          psi         q           t**2        rsi        rso )
      ipmax  = 200
      delpsi = psax/float(npsi)
      npsi1  = npsi+1
      npsi2  = npsi+2
      kax    = ((zax+zu)/dz)+1.
      jax    = ((rax-ri)/dr)+1.
c
      print 910,psax,qzero,tsqa(1.),rax,rax
      do 1000 nn=2,npsi2
      pss = psax-float(nn-1)*delpsi
      if(nn.eq.npsi1) pss=0.
      if(nn.eq.npsi2) pss=psax*0.05
c
c look for starting point
c
      q = 0.
      qc= 0.
      qd= 0.
      do 300 k=kax,km
      psa = psi(jax,k)
      psb = psi(jax,k+1)
      if ((psb-pss)*(pss-psa))  300,180,200
  180 if (psb-pss)              200,300,200
  200 ksur = k
      go to 302
  300 continue
      go to 1000
  302 do 902 ipage = 1,2
      if(ipage.eq.1) then
         isur = 4
         jsur = jax
         jlo  = jax
         jup  = jmax
         endif
      if(ipage.eq.2) then
         isur = 2
         jsur = jax-1
         jlo  = 1
         jup  = jax
         endif
      j = jsur
      k = ksur
      is= isur
c
c    draw one field line
c
      do 800 ip = 1,ipmax
      if (ip-1)  400,400,410
  400 ie = is
      go to 650
  410 ip1 = 0
      ip2 = 0
      ip3 = 0
      ip4 = 0
      if (((psi(j+1,k)-pss)*(pss-psi(j,k))).gt.0.)  ip1 = 1
      if (((psi(j+1,k+1)-pss)*(pss-psi(j+1,k))).gt.0.)  ip2 = 1
      if (((psi(j+1,k+1)-pss)*(pss-psi(j,k+1))).gt.0.)  ip3 = 1
      if (((psi(j,k+1)-pss)*(pss-psi(j,k))).gt.0.)  ip4 = 1
      go to (610,620,630,640),is
  610 if (ip2)  612,612,611
  611 ie = 2
      go to 650
  612 if (ip3)  614,614,613
  613 ie = 3
      go to 650
  614 if (ip4) 990,990,615
  615 ie = 4
      go to 650
  620 if (ip3)  622,622,621
  621 ie = 3
      go to 650
  622 if (ip4)  624,624,623
  623 ie = 4
      go to 650
  624 if (ip1)  990,990,625
  625 ie = 1
      go to 650
  630 if (ip4)  632,632,631
  631 ie = 4
      go to 650
  632 if (ip1)  634,634,633
  633 ie = 1
      go to 650
  634 if (ip2)  990,990,635
  635 ie = 2
      go to 650
  640 if (ip1)  642,642,641
  641 ie = 1
      go to 650
  642 if (ip2)  644,644,643
  643 ie = 2
      go to 650
  644 if (ip3)  990,990,645
  645 ie = 3
  650 go to (710,720,730,740),ie
  710 rss(ip) = r(j)+(dr*((pss-psi(j,k))/(psi(j+1,k)-psi(j,k))))
      zss(ip) = z(k)
      go to 742
  720 rss(ip) = r(j+1)
      zss(ip) = z(k)+(dz*((pss-psi(j+1,k))/(psi(j+1,k+1)-psi(j+1,k))))
      go to 742
  730 rss(ip) = r(j)+(dr*((pss-psi(j,k+1))/(psi(j+1,k+1)-psi(j,k+1))))
      zss(ip) = z(k+1)
      go to 742
  740 rss(ip) = r(j)
      zss(ip) = z(k)+(dz*((pss-psi(j,k))/(psi(j,k+1)-psi(j,k))))
  742 if (ip-1)              800,800,746
  746 ps4 = psi(j,k)
      ps3 = psi(j,k+1)
      ps2 = psi(j+1,k+1)
      ps1 = psi(j+1,k)
c
c     compute q(psi)
c
      dpsdz = (ps2+ps3-ps1-ps4)/(2.*dz)
      dpsdr = (ps1+ps2-ps3-ps4)/(2.*dr)
      drq = rss(ip)-rss(ip-1)
      dzq = zss(ip)-zss(ip-1)
      dll = sqrt((drq*drq)+(dzq*dzq))
      graps = sqrt((dpsdr*dpsdr)+(dpsdz*dpsdz))
      arr = rss(ip)+rss(ip-1)
      dq  = dll/(graps*arr)
      dqc = (dll*arr)/graps
      dqd = (dll*graps)/arr
      q = q + dq
      qc= qc+ dqc
      qd= qd+ dqd
c
c
c
  750 go to (751,752,753,754),ie
  751 k = k-1
      is = 3
      go to 755
  752 j = j+1
      is = 4
      go to 755
  753 k = k+1
      is = 1
      go to 755
  754 j = j-1
      is = 2
  755 if(k-1)      810,760,760
  760 if(j-jlo)    810,770,770
  770 if(k-km)     780,780,810
  780 if(j-jup+1)  800,800,810
  800 continue
      go to 990
  810 if(abs(rss(1)-ri).lt.1.0e-6)  go to 902
      if(abs(rss(1)-ro).lt.1.0e-6)  go to 902
      if(abs(rss(ip)-ri).lt.1.0e-6) go to 902
      if(abs(rss(ip)-ro).lt.1.0e-6) go to 902
      if(abs(zss(1)+zu).lt.1.0e-6)  go to 902
      if(abs(zss(1)-zu).lt.1.0e-6)  go to 902
      if(abs(zss(ip)+zu).lt.1.0e-6) go to 902
      if(abs(zss(ip)-zu).lt.1.0e-6) go to 902
      if(nn.ne.npsi1) go to 860
      if(ipage.eq.1) then
         ipp1 = ip-1
         if((zss(1)-zss(ipp1)).gt.0.) then
            do 830 i=1,ipp1
            rbou(i) = rss(i)
            zbou(i) = zss(i)
  830       continue
            endif
         if((zss(1)-zss(ipp1)).le.0.) then
            do 835 i=1,ipp1
            rbou(i) = rss(ipp1+2-i)
            zbou(i) = zss(ipp1+2-i)
  835       continue
            endif
         endif
      if(ipage.eq.2) then
         ipp2 = ip-1
         if((zss(1)-zss(ipp2)).le.0.) then
            do 840 i=1,ipp2
            rbou(i+ipp1) = rss(i)
            zbou(i+ipp1) = zss(i)
  840       continue
            endif
         if((zss(1)-zss(ipp2)).gt.0.) then
            do 850 i=1,ipp2
            rbou(i+ipp1) = rss(ipp2+2-i)
            zbou(i+ipp1) = zss(ipp2+2-i)
  850       continue
            ippx = ipp1+ipp2
            ipps = ippx
            endif
         endif
c
  860 if(nn.ne.npsi2) go to 902
      if(ipage.eq.1) then
         ipp1 = ip-1
         if((zss(1)-zss(ipp1)).gt.0.) then
            do 862 i=1,ipp1
            rboug(i) = rss(i)
            zboug(i) = zss(i)
  862       continue
            endif
         if((zss(1)-zss(ipp1)).le.0.) then
            do 865 i=1,ipp1
            rboug(i) = rss(ipp1+2-i)
            zboug(i) = zss(ipp1+2-i)
  865       continue
            endif
         endif
      if(ipage.eq.2) then
         ipp2 = ip-1
         if((zss(1)-zss(ipp2)).le.0.) then
            do 870 i=1,ipp2
            rboug(i+ipp1) = rss(i)
            zboug(i+ipp1) = zss(i)
  870       continue
            endif
         if((zss(1)-zss(ipp2)).gt.0.) then
            do 872 i=1,ipp2
            rboug(i+ipp1) = rss(ipp2+2-i)
            zboug(i+ipp1) = zss(ipp2+2-i)
  872       continue
            ippx = ipp1+ipp2
            endif
         endif
c
  902 continue
c
c                 elongation and triangularity
c
      ippx = ipp1 + ipp2
      if(nn.eq.npsi1) then
         eftr =  100.
         ritr = -100.
         botz =  100.
         topz = -100.
         do 903 i=1,ippx
         if(rbou(i).lt.eftr) eftr=rbou(i)
         if(rbou(i).gt.ritr) ritr=rbou(i)
         if(zbou(i).lt.botz) then 
            botz = zbou(i)
            ibot = i
            endif
         if(zbou(i).gt.topz) then
            topz = zbou(i)
            itop = i
            endif
  903    continue
c
C AM: 20.06.2019 Initialise topr and bottomr
         topr = 0
         botr = 0
         do 904 iii=1,2
         if(iii.eq.1) i=ibot
         if(iii.eq.2) i=itop
            if((i.ne.1).and.(i.ne.ippx)) then
               rboum1 = rbou(i-1)
               zboum1 = zbou(i-1)
               rboup1 = rbou(i+1)
               zboup1 = zbou(i+1)
               endif
            if(i.eq.1) then
               rboum1 = rbou(ippx)
               zboum1 = zbou(ippx)
               rboup1 = rbou(i+1)
               zboup1 = zbou(i+1)
               endif
            if(i.eq.ippx) then
               rboum1 = rbou(i-1)
               zboum1 = zbou(i-1)
               rboup1 = rbou(1)
               zboup1 = zbou(1)
               endif
            rbouze = rbou(i)
            zbouze = zbou(i)
            rboums = rboum1*rboum1
            rbouzs = rbouze*rbouze
            rboups = rboup1*rboup1     
         if ((((rboums-rbouzs)*(rbouze-rboup1)-
     1      (rbouzs-rboups)*(rboum1-rbouze)).eq.0.).or.
     2      ((rboum1-rbouze).eq.0.)) then
C AM: 20.06.2019
C   In some cases, the test will fail and botr/topr will not be set.
C this would cause a problem later when computing cappa/delta
C   A quick and dirty fix was introduced by initialising those variables to 0 [see above]
C In this case, elongation and triangularity values are wrong.
         write (0,"(A,A,i3)") "WARNING: Failed to compute kappa/delta",
     1           " for equilibrium #", nrun
         igood = 0
         go to 904
         endif
            cccc =((zboum1-zbouze)*(rbouze-rboup1)-
     1             (zbouze-zboup1)*(rboum1-rbouze))/
     2            ((rboums-rbouzs)*(rbouze-rboup1)-
     3             (rbouzs-rboups)*(rboum1-rbouze))
            bbbb =(zboum1-zbouze-(cccc*(rboums-rbouzs)))/(rboum1-rbouze)
            if(iii.eq.1) botr = -bbbb/(2.*cccc)
            if(iii.eq.2) topr = -bbbb/(2.*cccc)
  904 continue 
c 
         cappab = (topz-botz)/(ritr-eftr)
         deltab = (eftr+ritr-botr-topr)/(ritr-eftr)
         deltbu = (eftr+ritr-topr-topr)/(ritr-eftr)
         deltbl = (eftr+ritr-botr-botr)/(ritr-eftr)
         endif
c
      if(nn.eq.npsi2) then
         eftr =  100.
         ritr = -100.
         botz =  100.
         topz = -100.
         do 905 i=1,ippx
         if(rboug(i).lt.eftr) eftr=rboug(i)
         if(rboug(i).gt.ritr) ritr=rboug(i)
         if(zboug(i).lt.botz) then 
            botz = zboug(i)
            ibot = i
            endif
         if(zboug(i).gt.topz) then
            topz = zboug(i)
            itop = i
            endif
  905    continue
c
         do 906 iii=1,2
         if(iii.eq.1) i=ibot
         if(iii.eq.2) i=itop
            if((i.ne.1).and.(i.ne.ippx)) then
               rboum1 = rboug(i-1)
               zboum1 = zboug(i-1)
               rboup1 = rboug(i+1)
               zboup1 = zboug(i+1)
               endif
            if(i.eq.1) then
               rboum1 = rboug(ippx)
               zboum1 = zboug(ippx)
               rboup1 = rboug(i+1)
               zboup1 = zboug(i+1)
               endif
            if(i.eq.ippx) then
               rboum1 = rboug(i-1)
               zboum1 = zboug(i-1)
               rboup1 = rboug(1)
               zboup1 = zboug(1)
               endif    
            rbouze = rboug(i)
            zbouze = zboug(i)
            rboums = rboum1*rboum1
            rbouzs = rbouze*rbouze
            rboups = rboup1*rboup1
            cccc =((zboum1-zbouze)*(rbouze-rboup1)-
     1             (zbouze-zboup1)*(rboum1-rbouze))/
     2            ((rboums-rbouzs)*(rbouze-rboup1)-
     3             (rbouzs-rboups)*(rboum1-rbouze))
            bbbb =(zboum1-zbouze-(cccc*(rboums-rbouzs)))/(rboum1-rbouze)
            if(iii.eq.1) botr = -bbbb/(2.*cccc)
            if(iii.eq.2) topr = -bbbb/(2.*cccc)
  906 continue   
         cappa95 = (topz-botz)/(ritr-eftr)
         delta95 = (eftr+ritr-botr-topr)/(ritr-eftr)
         delt95u = (eftr+ritr-topr-topr)/(ritr-eftr)
         delt95l = (eftr+ritr-botr-botr)/(ritr-eftr)
         endif
c
c     complete calculation of q(psi)
c
      phi = pss/psax
      tsq = tsqa(phi)
      pprim = pprime(phi)
      prpest0(nn) = pres(phi)*um
      pppest0(nn) = -pprim*um
      gsqa = (rze*rze*bzeru*bzeru)+tsq
      ajpest2(nn) = ((pprim*qc)/(4.*q)+
     1   (ttprim(phi)*(1.+(qd/(gsqa*q)))))*um
      qq(nn) = (q/pi)*sqrt(gsqa)
      psss(nn) = pss
      xsv2(nn) = -pss
      do 909 j=1,nr
      if(pss.gt.psi(j,kax).and.pss.le.psi(j+1,kax))
     1 rsi(nn)=r(j)+dr*((pss-psi(j,kax))/(psi(j+1,kax)-psi(j,kax)))
      if(pss.le.psi(j,kax).and.pss.gt.psi(j+1,kax))
     1 rso(nn)=r(j)+dr*((pss-psi(j,kax))/(psi(j+1,kax)-psi(j,kax)))
  909 continue
      if(nn.le.npsi1) print 910,pss,qq(nn),tsq,rsi(nn),rso(nn)
  910 format (5x,5e12.4)
      go to 1000
  990 print 992
  992 format (25h field line finds no exit)
 1000 continue
      return
      end
c
c
c
c
      subroutine measur
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
      do 835 n=1,nprob
      angl = (thbee(n)*pi)/180.
      coth = cos(angl)
      sith = sin(angl)
      do 830 j=2,jmax
      if(r(j)-rbee(n)) 830,831,831
  830 continue
      j=jmax
  831 jx = j-1
      rfr = (rbee(n)-r(jx))/dr
      do 832 k=2,kmax
      if(z(k)-zbee(n)) 832,833,833
  832 continue
      k=kmax
  833 kx = k-1
      zfr = (zbee(n)-z(kx))/dz
      h1 = psi(jx+1,kx)-psi(jx,kx)
      h2 = psi(jx+1,kx+1)-psi(jx,kx+1)
      h3 = psi(jx,kx+1)-psi(jx,kx)
      h4 = psi(jx+1,kx+1)-psi(jx+1,kx)
      bz = (h1+zfr*(h2-h1))/(dr*rbee(n))
      br =-(h3+rfr*(h4-h3))/(dz*rbee(n))
      bee(n) = br*coth + bz*sith
  835 continue
c
      do 845 n=1,nloop
      do 840 j=2,jmax
      if(r(j)-rflu(n)) 840,841,841
  840 continue
      j=jmax
  841 jx=j-1
      rfr = (rflu(n)-r(jx))/dr
      do 842 k=2,kmax
      if(z(k)-zflu(n)) 842,843,843
  842 continue
      k=kmax
  843 kx=k-1
      zfr = (zflu(n)-z(kx))/dz
      h1 = psi(jx+1,kx)-psi(jx,kx)
      h2 = psi(jx+1,kx+1)-psi(jx,kx+1)
      h3 = psi(jx,kx+1)-psi(jx,kx)
      h4 = psi(jx+1,kx+1)-psi(jx+1,kx)
      flux1=(psi(jx,kx)+rfr*h1)*(1.-zfr)+(psi(jx,kx+1)+rfr*h2)*zfr
      flux2=(psi(jx,kx)+zfr*h3)*(1.-rfr)+(psi(jx+1,kx)+zfr*h4)*rfr
      flux(n) = (flux1+flux2)*pi
  845 continue
c
      do 902 k=1,ngroup
      gcurr(k) = 0.
      ifi = isvf(k)
      if(ifi.le.icur) gcurr(k)=(curr(ifi)*1000.)/tvf(ifi)
  902 continue
c
      return
      end
c
c
c
c
c
      subroutine wrieq
      include 'fbtmgams.inc'
      include 'fbt.inc'
      if(nrun.eq.1) open(unit=46,file='fbtequi.dat',status='replace')
c
      call measur
c
      write(46,100)
  100 format(45h Free Boundary Tokamak Equilibrium Code (FBT)) 
      write(46,101)
  101 format(45h            Standard Output File             )
      write(46,102) neqtcv
  102 format (33h Equilibrium No                  ,i15)
      write(46,400) placur
  400 format (33h plasma current (kA)             ,e15.5)
      write(46,402) betai
  402 format (33h beta poloidal                   ,e15.5)
      write(46,403) betat
  403 format (33h beta toroidal                   ,e15.5)
      write(46,406) rax
  406 format (33h radius of magnetic axis         ,e15.5)
      write(46,407) zax
  407 format (33h z-position of magnetic axis     ,e15.5)
      write(46,408) psax
  408 format (33h psi on axis                     ,e15.5)
      write(46,410) dsm
  410 format (33h max current density (ka/m**2)   ,e15.5)
      pmaxm= pmax*1000.
      write(46,412) pmaxm
  412 format (33h max pressure (mks)              ,e15.5)
      write(46,414) eli
  414 format (33h internal inductance             ,e15.5)
      write(46,416) diamag
  416 format (33h plasma diamagnetism             ,e15.5)
      write(46,418) torflu
  418 format (33h plasma toroidal flux            ,e15.5)
      write(46,422) qzero
  422 format (33h q on axis                       ,e15.5)
      write(46,423) qq(npsi+2)
  423 format (33h q-95                            ,e15.5)
      write(46,424) bzero
  424 format (33h toroidal vacuum field at R=0.88 ,e15.5)
      write(46,425) rze
  425 format (33h rze                             ,e15.5)
      write(46,428) ellipt
  428 format (33h central ellipticity             ,e15.5)
      write(46,429) cappab
  429 format (33h edge elongation                 ,e15.5)
      write(46,430) cappa95
  430 format (33h elongation of 95% flux surface  ,e15.5)
      write(46,431) deltab
  431 format (33h edge triangularity              ,e15.5)
      write(46,432) delta95
  432 format (33h triangularity of 95% surface    ,e15.5)
      write(46,433) ri
  433 format (33h inner grid boundary (ri)        ,e15.5)
      write(46,434) ro
  434 format (33h outer grid boundary (ro)        ,e15.5)
      write(46,435) zu
  435 format (33h upper and lower grid boundaries ,e15.5)
      write(46,436) jmax
  436 format (33h number of radial mesh points    ,i15)
      write(46,438) kmax
  438 format (33h number of axial mesh points     ,i15)
      write(46,450)
  450 format (33h coordinates of mesh points      )
      write(46,452)
  452 format (33h r(j)                            )
      write(46,3000) (r(j),j=1,jmax)
      write(46,454)
  454 format (33h z(k)                            )
      write(46,3000) (z(k),k=1,kmax)
      write(46,460) iacx
  460 format (33h number of active X-points       ,i15)
      write(46,462) rxp(isone),zxp(isone)
  462 format (33h r,z of first active X-point     ,2e15.5)
      write(46,464) rxp(istwo),zxp(istwo)
  464 format (33h r,z of second active X-point    ,2e15.5)
      write(46,480)
  480 format (33h coil coordinates and currents   )
      write(46,482)
  482 format (41h      r              z              i(kA))
      do 500 i=1,icur
      write(46,502) rvf(i),zvf(i),curr(i)
  500 continue
  502 format (3e15.5)
      write(46,560)
  560 format (33h q-profile, q(psi)               )
      write(46,566)
  566 format (61h          psi            q              rsi               rso)
      npsip1  = npsi+1
      psss(1) = psax
      qq(1)   = qzero
      rsi(1)  = rax
      rso(1)  = rax
      do 570 i=1,npsip1
      write(46,910) psss(i),qq(i),rsi(i),rso(i)
  570 continue
      write(46,580)
  580 format (33h psi (poloidal flux per radian)  )
      do 600 k=1,kmax
      kq = kmax+1-k
      write(46,3000) (psi(j,kq),j=1,jmax)
  600 continue
      write(46,650)
  650 format (35h toroidal current density (kA/m**2))
      do 700 k=1,kmax
      kq = kmax+1-k
      write(46,3000) (dsh(j,kq),j=1,jmax)
  700 continue
      write(46,750)
  750 format (33h plasma pressure (kN/m**2)       )
      do 800 k=1,kmax
      kq = kmax+1-k
      write(46,3000) (psn(j,kq),j=1,jmax)
  800 continue
      write(46,820)
  820 format (33h toroidal field (T)              )
      do 824 k=1,kmax
      kq = kmax+1-k
      do 822 j=1,jmax
      phi = psi(j,k)/psax
      tsq = 0.
      if(psn(j,k).gt.0.) tsq=tsqa(phi)
      rss(j) = sqrt(rze*rze*bzeru*bzeru+tsq)/r(j)
  822 continue
      write(46,3000) (rss(j),j=1,jmax)
  824 continue
c
      write(46,846)    
  846 format (33h magnetic field measurements     )
      write(46,3000) (bee(n),n=1,nprob)
      write(46,847)
  847 format (33h flux loop measurements          )
      write(46,3000) (flux(n),n=1,nloop)
      write(46,860)
  860 format (46h r and z coordinates of plasma boundary points)
      do 862 i=1,ipps
      write(46,3000) rbou(i),zbou(i)
  862 continue
  910 format(4e15.5)
 3000 format(1x,6e12.5)
 3001 format(1x,12i6)
      if(nrun.eq.nruns) close(unit=46)
      return
      end
c
c
c
c
c
      subroutine wrikriz
      include 'fbtmgams.inc'
      include 'fbt.inc'
      open(unit=50,file='fbtkriz.dat',status='replace')
      write(50,426) jmax
  426 format (i10,3h nx)
      write(50,428) kmax
  428 format (i10,3h nz)
      write(50,430) ri
  430 format (f10.3,5h rmin)
      write(50,432) ro
  432 format (f10.3,5h rmax)
      write(50,434) -zu
  434 format (f10.3,5h zmin)
      write(50,436) zu
  436 format (f10.3,5h zmax)
      npsip1 = npsi+1
      write(50,438) npsip1
  438 format (i10,5h npsi)
      psax2p = psax*2.*pi
      write(50,448) psax2p
  448 format (f10.6,12h psi on axis)
      xxxxx = 0.
      write(50,450) xxxxx
  450 format (f10.6,23h psi on plasma boundary)
      write(50,452) rax
  452 format (f10.6,24h radius of magnetic axis)
      write(50,454) rze
  454 format (f10.6,6h rzero)
      write(50,456) bzero
  456 format (f10.6,6h bzero)
      do 500 i=1,npsip1
      phi = 1.-(float(i-1)/float(npsi))
      if(phi.lt.0.) phi=0.
      rss(i) = pres(phi)*1000.
      zss(i) = sqrt(tsqa(phi)+(rze*rze*bzeru*bzeru))
  500 continue
      write(50,502)
  502 format (11h p (N/m**2))
      write(50,504) (rss(i),i=1,npsip1)
      write(50,503)
  503 format (7h g (Tm))
      write(50,504) (zss(i),i=1,npsip1)
  504 format (5e15.6)
      do 600 k=1,kmax
      kq = kmax+1-k
      write(50,520) k
  520 format (18h psi (webers)   k=,i4)
      do 530 j=1,jmax
      rss(j) = psi(j,kq)*2.*pi
  530 continue
      write(50,504) (rss(j),j=1,jmax)
  600 continue
      close(unit=50)
      return
      end
c
c
c
c
      subroutine wrimeas
      include 'fbtmgams.inc'
      include 'fbt.inc'
      inquire(unit=47,opened=istat)
      if(.not.istat) then
      open(unit=47,file='fbtmeas.dat',status='replace')
      write(47,3001) nruns,icur,nprob,nloop,
     1               nr,nz,iallf,ngroup
      write(47,3000) ri,ro,zu,rhocu,filfac,rhovv
      write(47,3000) (rvf(i),i=1,iallf)
      write(47,3000) (zvf(i),i=1,iallf)
      write(47,3000) (hvf(i),i=1,iallf)
      write(47,3000) (wvf(i),i=1,iallf)
      write(47,3000) (tvf(i),i=1,iallf)
      write(47,3001) (kgroup(i),i=1,iallf)
      write(47,3001) (isvf(k),k=1,ngroup)
      write(47,3001) (isvl(k),k=1,ngroup)
      write(47,3001) (itype(k),k=1,ngroup)
      write(47,3000) (rbee(n),n=1,nprob)
      write(47,3000) (zbee(n),n=1,nprob)
      write(47,3000) (thbee(n),n=1,nprob)
      write(47,3000) (rflu(n),n=1,nloop)
      write(47,3000) (zflu(n),n=1,nloop)
      endif
c
c     [AM 16/06/2017]
c     Make computation and writing of measurements independent
c
c      call measur
      placamp = placur*1000.
c
      write(47,3001) nrun,ilie,ilia,ibro,ibzo
      write(47,3000) placamp,zax,bzero,rze,timeeq
      write(47,3000) rax,zax,rax2,zax2
      write(47,3000) zeecorb,zeecort
      write(47,3000) betai,eli,diamag,qq(npsi+1),qzero
      write(47,3000) area,cappab,shvsec
      write(47,3000) (gcurr(i),i=1,ngroup)
      write(47,3000) (bee(n),n=1,nprob)
      write(47,3000) (flux(n),n=1,nloop)
      if(ilie.gt.0) then
         write(47,3000) (rlim(i),i=1,ilie)
         write(47,3000) (zlim(i),i=1,ilie)
         endif
      if(ilia.gt.0) then
         write(47,3000) (rlia(n),n=1,ilia)
         write(47,3000) (zlia(n),n=1,ilia)
         endif
      if(ibro.gt.0)
     1   write(47,3000) (rbro(n),n=1,ibro),(zbro(n),n=1,ibro)
      if(ibzo.gt.0)
     1   write(47,3000) (rbzo(n),n=1,ibzo),(zbzo(n),n=1,ibzo)
      if(nrun.eq.nruns) close(unit=47)
 3000 format(1x,6e12.5)
 3001 format(1x,12i6)
      return
      end
c
c
c
c
c
      function tsqa(phi)
      include 'fbtmgams.inc'
      include 'fbt.inc'
      xxx = 1.-phi
      if(xxx.le.0.) xxx=0.
      if(emm.ge.1.) then
         phemm1 = phi**emmp1
         phemm2 = phemm1*phi
         phemm3 = phemm2*phi
         tsqa   = 2.*rze*rze*um*ratio*psax*ttfac*
     1            (((1.-ttbe)/emmp1)*phemm1
     2            +((ttbe*(1.-omega))/emmp2)*phemm2
     3            +((ttbe*omega)/emmp3)*phemm3)
         endif
      if(emm.lt.1.) then
         xoemm1 = xxx**oemmp1
         xoemm2 = xoemm1*xxx
         tsqa   = 2.*rze*rze*um*ratio*psax*ttfac*
     1            ((1.-xxx)
     2            +(((1.+ttbe)/oemmp1)*(xoemm1-1.))
     3            -((ttbe/oemmp2)*(xoemm2-1.)))
         endif
      return
      end
c
c
c
c
c
      function pres(phi)
      include 'fbtmgams.inc'
      include 'fbt.inc'
      xxx = 1.-phi
      if(xxx.le.0.) xxx=0.
      if(ell.gt.1.) then
         phell1 = phi**ellp1
         phell2 = phell1*phi
         pres   = ppfac*ratio*psax*
     1            ((phell1/ellp1)
     2            -((1.-ppal)*(phell2/ellp2)))
         endif
      if(ell.le.1.) then
         xoell1 = xxx**oellp1
         xoell2 = xoell1*xxx
         pres   = ppfac*ratio*psax*
     1            ((ppal*(1.-xxx)) + (((xoell1-1.)/oellp1)*
     2            (2.*ppal-1.))+(((1.-ppal)/oellp2)*(xoell2-1.)))
         endif
      return
      end
c
c
c
c
c
      function pprime(phi)
      include 'fbtmgams.inc'
      include 'fbt.inc'
      xxx = 1.-phi
      if(xxx.le.0.) xxx=0.
      if(ell.gt.1.) then
         phell  = phi**ell
         phell1 = phell*phi
         ppri   = ppal*phell+(1.-ppal)*(phell-phell1)
         endif
      if(ell.le.1.) then
         xoell  = xxx**oell
         xoell1 = xoell*xxx
         ppri   = ppal*(1.-xoell)+(1.-ppal)*(xoell-xoell1)
         endif
      pprime = ppri*ratio*ppfac
      return
      end
c
c
c
c
c
      function ttprim(phi)
      include 'fbtmgams.inc'
      include 'fbt.inc'
      xxx = 1.-phi
      if(xxx.le.0.) xxx=0.
      if(emm.ge.1.) then
         phemm  = phi**emm
         phemm1 = phemm*phi
         phemm2 = phemm1*phi
         ttpr   = phemm+ttbe*(phemm1-phemm+omega*(phemm2-phemm1))
         endif
      if(emm.lt.1.) then
         xoemm  = xxx**oemm
         xoemm1 = xoemm*xxx
         ttpr   = 1.-xoemm + ttbe*(xoemm1-xoemm)
         endif
      ttprim = ttpr*rze*rze*ratio*ttfac
      return
      end
c
c
c
c
c
      subroutine wrierat
      include 'fbtmgams.inc'
      include 'fbt.inc'
      dimension ikeep(mbou),theta(mbou)
      open(unit=49,file='fbterat.dat',status='replace')
c
c              look for max and min z
c
      ipp    =  ipp1+ipp2
      zbouma = -100.
      zboumi =  100.
      do 10 i=1,ipp
      if(zbou(i).gt.zbouma) then
         zbouma = zbou(i)
         ibouma = i
         endif
      if(zbou(i).lt.zboumi) then
         zboumi = zbou(i)
         iboumi = i
         endif
   10 continue
      zzero = (zbouma+zboumi)/2.
      zone  = (zbouma-zboumi)/2.
c
c               eliminate non-integer multiples of dz
c
      do 30 i=1,ipp
      ikeep(i) = 0
      if((i.eq.ibouma).or.(i.eq.iboumi)) then
         ikeep(i) = 1
         go to 30
         endif
      do 20 k=1,kmax
      if(abs(zbou(i)-z(k)).lt.1.0e-6) then
         ikeep(i) = 1
         go to 30
         endif
   20 continue
   30 continue
c
c                         reshuffling
c
      ii = 1
      do 50 i=1,ipp
      if(ikeep(i).eq.1) then
         if(i.ne.ii) rbou(ii)=rbou(i)
         if(i.ne.ii) zbou(ii)=zbou(i)
         if(i.eq.ibouma) ibouma=ii
         if(i.eq.iboumi) iboumi=ii
         ii = ii+1
         endif
   50 continue
      ipps = ii-1
c
c                 compute theta-values 
c
      do 60 i=1,ipps
      sinthet = (zbou(i)-zzero)/zone
      if(sinthet.gt.1.) sinthet=1.
      if(sinthet.lt.-1.) sinthet=-1.
      theta(i) = asin(sinthet)
      if(ibouma.lt.iboumi) then
      if((i.gt.iboumi).or.(i.le.ibouma)) theta(i)=-pi-theta(i)
      endif
      if(ibouma.ge.iboumi) then
      if((i.gt.iboumi).and.(i.le.ibouma)) theta(i)=-pi-theta(i)
      endif
   60 continue
c
c
c
      imat = 1+2*ifour
      imau = 1+imat
      ifous= 1+ifour
      do 80 i=1,ipps
      do 70 j=1,ifour
      sinma(j,i) = sin(theta(i)*float(j))
      cosma(j,i) = cos(theta(i)*float(j))
   70 continue
   80 continue
c
      ax(1,1) = float(ipps)
      do 90 l=1,ifour
      ax(1,1+l) = 0.
      do 88 i=1,ipps
      ax(1,1+l) = ax(1,1+l)+cosma(l,i)
   88 continue
   90 continue
      do 100 m=1,ifour
      ax(1,ifous+m) = 0.
      do 98 i=1,ipps
      ax(1,ifous+m) = ax(1,ifous+m)+sinma(m,i)
   98 continue
  100 continue
      ax(1,imau) = 0.
      do 110 i=1,ipps
      ax(1,imau) = ax(1,imau)+rbou(i)
  110 continue 
c
      do 200 k=1,ifour
      ax(1+k,1) = 0.
      do 120 i=1,ipps
      ax(1+k,1) = ax(1+k,1)+cosma(k,i)
  120 continue
      do 130 l=1,ifour
      ax(1+k,1+l) = 0.
      do 128 i=1,ipps
      ax(1+k,1+l) = ax(1+k,1+l)+cosma(l,i)*cosma(k,i)
  128 continue
  130 continue
      do 140 m=1,ifour
      ax(1+k,ifous+m) = 0.
      do 138 i=1,ipps
      ax(1+k,ifous+m) = ax(1+k,ifous+m)+sinma(m,i)*cosma(k,i)
  138 continue
  140 continue
      ax(1+k,imau) = 0.
      do 150 i=1,ipps
      ax(1+k,imau) = ax(1+k,imau)+rbou(i)*cosma(k,i)
  150 continue 
  200 continue
c
      do 300 n=1,ifour
      ax(ifous+n,1) = 0.
      do 220 i=1,ipps
      ax(ifous+n,1) = ax(ifous+n,1)+sinma(n,i)
  220 continue
      do 230 l=1,ifour
      ax(ifous+n,1+l) = 0.
      do 228 i=1,ipps
      ax(ifous+n,1+l) = ax(ifous+n,1+l)+cosma(l,i)*sinma(n,i)
  228 continue
  230 continue
      do 240 m=1,ifour
      ax(ifous+n,ifous+m) = 0.
      do 238 i=1,ipps
      ax(ifous+n,ifous+m)=ax(ifous+n,ifous+m)+sinma(m,i)*sinma(n,i)
  238 continue
  240 continue
      ax(ifous+n,imau) = 0.
      do 250 i=1,ipps
      ax(ifous+n,imau) = ax(ifous+n,imau)+rbou(i)*sinma(n,i)
  250 continue 
  300 continue
c
      call gauss(imat)
c
      alzero = cx(1)
      do 310 l=1,ifour
      albou(l) = cx(l+1)
      bebou(l) = cx(ifous+l)
  310 continue
c
      plac = placur*1000.
      effp = ppfac*ratio*1000.
      efft = ttfac*ratio*rze*rze*1000.
c
      write(49,400) neqtcv
      write(49,403) ppal,ell,effp
      write(49,403) ttbe,emm,efft
      write(49,403) rze,bzeru
      write(49,403) qzero,plac
      write(49,400) ifour
      write(49,403) zzero,zone,alzero
      write(49,401) (albou(i),i=1,ifour)
      write(49,401) (bebou(i),i=1,ifour)
c
  400 format(i5)
  401 format(e15.8)
  403 format(3e15.8)
c
      close(unit=49)
      return
      end
c
c
c
c
c
      subroutine wrinova
      include 'fbtmgams.inc'
      include 'fbt.inc'
      open(unit=48,file='fbtnova.dat',status='replace')
      ncycle = 0
      isyms  = 0
      ipest  = 1
      npsit  = npsi+1
      kmay   = ipp1+ipp2
      times  = 0.
      xaxes  = rax
      zaxes  = zax
      gzeros = rze*bzeru
      apls   = placur*um
      betas  = betat
      betaps = betai
      ali2s  = eli/2.
      qsaws  = 1.0
      psimins= -psax
      psilims= 0.
      prpest0(1) = pres(1.)*um
      pppest0(1) = -pprime(1.)*um
      ajpest2(1) = (rax*rax*pprime(1.)+ttprim(1.))*um
      xsv2(1)    = -psax
c
      kmayp1 = kmay+1
      kmay4  = ipp2/2
      kdiff  = kmay-kmay4
      kmay41 = kmay4+1
      do 10 k=1,kmay4
      xpass(k) = rbou(k+kdiff)
      zpass(k) = zbou(k+kdiff)
   10 continue
      do 20 k=kmay41,kmay
      xpass(k) = rbou(k-kmay4)
      zpass(k) = zbou(k-kmay4)
   20 continue
      xpass(kmayp1) = xpass(1)
      zpass(kmayp1) = zpass(1)
c
      write(48,100) neqtcv,isyms,ipest,npsit,kmay
      write(48,101) psirat,xaxes,zaxes
      write(48,101) gzeros,apls,betas
      write(48,101) betaps,ali2s,qsaws
      write(48,101) psimins,psilims
      write(48,102) (prpest0(j),j=1,npsit)
      write(48,102) (pppest0(j),j=1,npsit)
      write(48,102) (ajpest2(j),j=1,npsit)
      write(48,102) (xsv2(j),j=1,npsit)
      write(48,102) (xpass(k),k=1,kmayp1)
      write(48,102) (zpass(k),k=1,kmayp1)
c
  100 format(5i5)
  101 format(3e15.8)
  102 format(e15.8)
      close(unit=48)
      return
      end
c
       subroutine wrichea
       include 'fbtmgams.inc'
       include 'fbt.inc'
       dimension sss(200),pprims(200),ttprims(200)
       open(unit=49,file='fbtchea.dat',status='replace')
c
c              look for max and min r and z
c
       zbouma = -100.
       zboumi =  100.
       rbouma = -100.
       rboumi =  100.
       do 10 i=1,ipps
       if(zbou(i).gt.zbouma) zbouma = zbou(i)
       if(zbou(i).lt.zboumi) zboumi = zbou(i)
       if(rbou(i).gt.rbouma) rbouma = rbou(i)
       if(rbou(i).lt.rboumi) rboumi = rbou(i)
   10  continue
c
       zzero = (zbouma+zboumi)/2.
       zone  = (zbouma-zboumi)/2.
       rzero = (rbouma+rboumi)/2.
       rone  = (rbouma-rboumi)/2.
       aspi  = rone/rzero
       rnorm = 0.88
       znorm = zzero/rnorm
       pedge = 0.
       isource = 100
       isep  = 1
       bzerch = (bzeru*rze)/rnorm
c------------------------note that bzerch=bzero-------------------
       cpp   = -(4.*pi*1.0e-7*rnorm*rnorm*1000.)/bzerch
       ctt   = -(4.*pi*1.0e-4)/bzerch
       placuch = (placur*4.*pi*1.0e-7)/(rnorm*bzerch)
c
       do i=1,isource
       sss(i)= float(i-1)/float(isource-1)
       phi   = 1.-sss(i)*sss(i)
       xxx = 1.-phi
       if(xxx.le.0.) xxx=0.
       if(ell.gt.1.) then
          phell  = phi**ell
          phell1 = phell*phi
          ppri   = ppal*phell+(1.-ppal)*(phell-phell1)
          endif   
       if(ell.le.1.) then
          xoell  = xxx**oell
          xoell1 = xoell*xxx
          ppri   = ppal*(1.-xoell)+(1.-ppal)*(xoell-xoell1)
          endif
       pprims(i) = ppri*ratio*ppfac*cpp
       if(emm.ge.1.) then
       phemm  = phi**emm
          phemm1 = phemm*phi
          phemm2 = phemm1*phi
          ttpr   = phemm+ttbe*(phemm1-phemm+omega*(phemm2-phemm1))
          endif
       if(emm.lt.1.) then
          xoemm  = xxx**oemm
          xoemm1 = xoemm*xxx
          ttpr   = 1.-xoemm + ttbe*(xoemm1-xoemm)
          endif
       ttprims(i) = ttpr*rze*rze*ratio*ttfac*ctt
       enddo
c
       write(49,103) aspi
  103  format (e18.8,40h Inverse Aspect Ratio                   )
       write(49,104) znorm
  104  format (e18.8,40h Normalized Z-Position                  )
       write(49,105) pedge
  105  format (e18.8,40h Pressure at Edge                       )
       write(49,106) ipps
  106  format    (i5,40h Nb of Plasma Boundary Points           ) 
       do i=1,ipps
       write(49,107) rbou(i)/rnorm,zbou(i)/rnorm
       enddo
  107  format (2e18.8)
       write(49,108) isource
  108  format    (i5,40h Nb of Source Function Profile Points   )
       write(49,140) isep
  140  format(i5)
       do i=1,isource
       write(49,109) sss(i)
       enddo
       do i=1,isource
       write(49,109) pprims(i)
       enddo  
       do i=1,isource
       write(49,109) ttprims(i)
       enddo
  109  format (e18.8)
       write(49,120) rnorm
  120  format (e18.8,40h Radius for Normalization               )
       write(49,121) bzeru
  121  format (e18.8,40h Toroidal Vacuum Field at R=R_major     )
       write(49,132) bzerch
  132  format (e18.8,40h Toroidal Vacuum Field at R=0.88m       )
       write(49,122) placur*1000.
  122  format (e18.8,40h Plasma Current                         )
       write(49,131) placuch*1000.
  131  format (e18.8,40h CHEASE Plasma Current                  ) 
       write(49,123) betai
  123  format (e18.8,40h Beta Poloidal                          )
       write(49,124) betat
  124  format (e18.8,40h Beta Toroidal                          )
       write(49,125) eli
  125  format (e18.8,40h Internal Inductance                    )
       write(49,126) qzero
  126  format (e18.8,40h q_zero                                 )
       write(49,127) qq(npsi+1)
  127  format (e18.8,40h q_edge                                 )
       write(49,130) qq(npsi+2)
  130  format (e18.8,40h q_95                                   )
       write(49,128) rzero
  128  format (e18.8,40h Major Radius                           )
       write(49,129) rone
  129  format (e18.8,40h Minor Radius                           )
c
       write(49,100)
  100  format(45h Free Boundary Tokamak Equilibrium Code (FBT))
       write(49,101)
  101  format(45h             CHEASE Output File              )
       write(49,102) ishot
  102  format   (i15,40h TCV Shot No.                           )
c
       close(unit=49)
       return
       end        
c
c
c
      subroutine pseuliu
      include 'fbtmgams.inc'
      include 'fbt.inc'  

      pi = 3.141592653589793
c
      area_mg(nrun)      = area
      beta_pol(nrun)     = betai
      beta_tor(nrun)     = betat
      delta_95(nrun)     = delta95
      delta_95_bot(nrun) = delt95l
      delta_95_top(nrun) = delt95u
      delta_edge(nrun)   = deltab
      delta_ed_bot(nrun) = deltbl
      delta_ed_top(nrun) = deltbu
      diamag_mg(nrun)    = diamag
      indx_act_xp(nrun)  = isone
      indx_sec_xp(nrun)  = istwo
      kappa_95(nrun)     = cappa95
      kappa_edge(nrun)   = cappab
      kappa_zero(nrun)   = ellipt
      lambda(nrun)       = betai+eli/2.
      l_i(nrun)          = eli
      npts_contour(nrun) = ipps+1
      n_act_xpts(nrun)   = iacx
      n_xpts(nrun)       = ismax
      p_axis(nrun)       = pres(1.)*1000.
      q_95(nrun)         = qq(npsi+2)
      q_edge(nrun)       = qq(npsi+1)
      q_zero(nrun)       = qzero
      r_axis(nrun)       = rax
c  Found that volume was undercalculated by a factor of 2*pi until now
c    (SC, 01/09/10)
      volume(nrun)       = volum*2.*pi
      z_axis(nrun)       = zax
c
      do j = 1,nr+1
      r_psi(j)           = r(j)
      enddo
      do k = 1,km+1
      z_psi(k)           = z(k)
      enddo
c
      do k = 1,km
      do j = 1,nr
      pressure(j,k,nrun) = psn(j,k)*1000.
      enddo
      enddo
c
c  First point (magnetic axis) needs to be introduced here (before it was done
c    only when wrieq was called) (SC, 14/04/07)
c
      psss(1) = psax
      qq(1)   = qzero
      rsi(1)  = rax
      rso(1)  = rax
c
      do i = 1,npsi+1
      phi                = 1.-(float(i-1)/float(npsi))
      q_psi(i,nrun)      = qq(i)
      r_max_psi(i,nrun)  = rso(i)
      r_min_psi(i,nrun)  = rsi(i)
      enddo
c
      do is = 1,ismax
      r_xpts(is,nrun)     = rxp(is)
      z_xpts(is,nrun)     = zxp(is)
      enddo
c 
      ippt = ipps + 1
      rbou(ippt) = rbou(1)
      zbou(ippt) = zbou(1)
      do i = 1,ippt
      r_contour(i,nrun)   = rbou(i)
      z_contour(i,nrun)   = zbou(i)
      enddo
c
c  I_p-sign dependent quantities (SC, 13/04/07)
c
      if (iohfb.ge.0) then
        i_p(nrun)          = placur*1000.
        psi_axis(nrun)     = psax*2.*pi
c
c  Go to km+1 and nr+1, not km and nr (SC, 14/04/07)
c
        do k = 1,km+1
        do j = 1,nr+1
          psi_mg(j,k,nrun)   = psi(j,k)*2.*pi
          j_tor(j,k,nrun)    = dsh(j,k)*1000.
        enddo
        enddo
c
        do i = 1,npsi+1
          phi                = 1.-(float(i-1)/float(npsi))
c  The next two factors need to be multiplied by 1000 to go to MKSA 
c    (SC, 13/10/09)
          pprime_psi(i,nrun) = 1000.*pprime(phi)
          ttprime_psi(i,nrun)= 1000.*ttprim(phi)
          psi_values(i,nrun) = psss(i)*2.*pi
        enddo
c
        do is = 1,ismax
          psi_xpts(is,nrun)   = psisa(is)*2.*pi
        enddo
c
      else
        i_p(nrun)          = -1.*placur*1000.
        psi_axis(nrun)     = -1.*psax*2.*pi
c
c
c  Go to km+1 and nr+1, not km and nr (SC, 14/04/07)
c
        do k = 1,km+1
        do j = 1,nr+1
          psi_mg(j,k,nrun)   = -1.*psi(j,k)*2.*pi
          j_tor(j,k,nrun)    = -1.*dsh(j,k)*1000.
        enddo
        enddo
c
        do i = 1,npsi+1
          phi                = 1.-(float(i-1)/float(npsi))
c  The next two factors need to be multiplied by 1000 to go to MKSA 
c    (SC, 13/10/09)
          pprime_psi(i,nrun) = -1000.*pprime(phi)
          ttprime_psi(i,nrun)= -1000.*ttprim(phi)
          psi_values(i,nrun) = -1.*psss(i)*2.*pi
        enddo
c
        do is = 1,ismax
          psi_xpts(is,nrun)   = -1.*psisa(is)*2.*pi
        enddo
c
      endif
c
c  B_T-sign dependent quantities (SC, 13/04/07)
c
      if (if36fb.ge.0) then
        tor_flux(nrun)     = torflu
      else
        tor_flux(nrun)     = -1.*torflu
      endif
c
      return
      end
c
c
c
      subroutine pseufbt
      include 'fbtmgams.inc'
      include 'fbt.inc'
c
      placamp_t(nrun)  = placur*1000.
      rze_t    (nrun)  = rze
      zze_t    (nrun)  = zze
      rax_t    (nrun)  = rax
      zax_t    (nrun)  = zax
      rax2_t   (nrun)  = rax2
      zax2_t   (nrun)  = zax2
      shvsec_t (nrun)  = shvsec
      do i = 1,ngroup
       gcurr_t(i+(nrun-1)*ngroup) = gcurr(i)
      enddo
      print *,ngroup,nrun,gcurr
      do i = 1,nprob
       bee_t(i+(nrun-1)*nprob)   = bee(i)
      enddo
      do i = 1,nloop
       flux_t(i+(nrun-1)*nloop)  = flux(i)
      enddo

      return
      end

	SUBROUTINE FBTMDSPUT(SHOT)
	INTEGER SHOT
	INCLUDE 'fbtmgams.inc'
	INCLUDE 'fbt.inc'
	INCLUDE 'mdslib.inc'
	INTEGER MDSSTATUS MDSLENGTH
	
	REAL DUMMY(NPSI+1,NRUNS)
	
	MDSSTATUS = MDSSETDEFAULT('\MGAMS.FBTE'//CHAR(0))
	IF (MOD(MDSSTATUS,2).EQ.0) THEN
	 PRINT *,'FBTMDSPUT - error selecting MGAMS branch'
	 STOP
	ENDIF

	MDSSTATUS = MDSPUT('AREA'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,AREA_MG,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing AREA'
	MDSSTATUS = MDSPUT('BETA_POL'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,BETA_POL,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing BETA_POL'
	MDSSTATUS = MDSPUT('BETA_TOR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,BETA_TOR,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing BETA_TOR'
	MDSSTATUS = MDSPUT('DELTA_95'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_95,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_95'
	MDSSTATUS = MDSPUT('DELTA_95_BOT'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_95_BOT,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_95_BOT'
	MDSSTATUS = MDSPUT('DELTA_95_TOP'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_95_TOP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_95_TOP'
	MDSSTATUS = MDSPUT('DELTA_EDGE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_EDGE,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_EDGE'
	MDSSTATUS = MDSPUT('DELTA_ED_BOT'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_ED_BOT,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_ED_BOT'
	MDSSTATUS = MDSPUT('DELTA_ED_TOP'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DELTA_ED_TOP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DELTA_ED_TOP'
	MDSSTATUS = MDSPUT('DIAMAG'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DIAMAG_MG,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing DIAMAG'
	MDSSTATUS = MDSPUT('INDX_ACT_XP'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_LONG,INDX_ACT_XP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing INDX_ACT_XP'
	MDSSTATUS = MDSPUT('INDX_SEC_XP'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_LONG,INDX_SEC_XP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing INDX_SEC_XP'
	MDSSTATUS = MDSPUT('I_P'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,I_P,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing I_P'
	MDSSTATUS = MDSPUT('J_TOR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,J_TOR,MRMP,MAMP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing J_TOR'
	MDSSTATUS = MDSPUT('KAPPA_95'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,KAPPA_95,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing KAPPA_95'
	MDSSTATUS = MDSPUT('KAPPA_EDGE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,KAPPA_EDGE,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing KAPPA_EDGE'
	MDSSTATUS = MDSPUT('KAPPA_ZERO'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,KAPPA_ZERO,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing KAPPA_ZERO'
	MDSSTATUS = MDSPUT('LAMBDA'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,LAMBDA,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing LAMBDA'
	MDSSTATUS = MDSPUT('L_I'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,L_I,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing L_I'
	MDSSTATUS = MDSPUT('NPTS_CONTOUR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_LONG,NPTS_CONTOUR,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing NPTS_CONTOUR'
	MDSSTATUS = MDSPUT('N_ACT_XPTS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_LONG,N_ACT_XPTS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing N_ACT_XPTS'
	MDSSTATUS = MDSPUT('N_XPTS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_LONG,N_XPTS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing N_XPTS'
	DUMMY=PPRIME_PSI(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('PPRIME_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PPRIME_PSI'
	MDSSTATUS = MDSPUT('PRESSURE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,PRESSURE,MRMP,MAMP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PRESSURE'
	MDSSTATUS = MDSPUT('PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,PSI_MG,MRMP,MAMP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PSI'
	MDSSTATUS = MDSPUT('PSI_AXIS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,PSI_AXIS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PSI_AXIS'
	DUMMY=PSI_VALUES(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('PSI_VALUES'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PSI_VALUES'
	MDSSTATUS = MDSPUT('PSI_XPTS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,PSI_XPTS,MXPTS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PSI_XPTS'
	MDSSTATUS = MDSPUT('P_AXIS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,P_AXIS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing P_AXIS'
	MDSSTATUS = MDSPUT('Q_95'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Q_95,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Q_95'
	MDSSTATUS = MDSPUT('Q_EDGE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Q_EDGE,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Q_EDGE'
	DUMMY=Q_PSI(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('Q_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Q_PSI'
	MDSSTATUS = MDSPUT('Q_ZERO'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Q_ZERO,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Q_ZERO'
	MDSSTATUS = MDSPUT('R_AXIS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,R_AXIS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_AXIS'
	MDSSTATUS = MDSPUT('R_CONTOUR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,R_CONTOUR,MBOU,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_CONTOUR'
	DUMMY=R_MAX_PSI(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('R_MAX_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_MAX_PSI'
	DUMMY=R_MIN_PSI(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('R_MIN_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_MIN_PSI'
	MDSSTATUS = MDSPUT('R_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,R_PSI,MRMP,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_PSI'
	MDSSTATUS = MDSPUT('R_XPTS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,R_XPTS,MXPTS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing R_XPTS'
	MDSSTATUS = MDSPUT('TOR_FLUX'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,TOR_FLUX,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing TOR_FLUX'
	DUMMY=TTPRIME_PSI(1:NPSI+1,1:NRUNS)
	MDSSTATUS = MDSPUT('TTPRIME_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,DUMMY,NPSI+1,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing TTPRIME_PSI'
	MDSSTATUS = MDSPUT('VOLUME'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,VOLUME,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing VOLUME'
	MDSSTATUS = MDSPUT('Z_AXIS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Z_AXIS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Z_AXIS'
	MDSSTATUS = MDSPUT('Z_CONTOUR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Z_CONTOUR,MBOU,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Z_CONTOUR'
	MDSSTATUS = MDSPUT('Z_PSI'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Z_PSI,MAMP,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Z_PSI'
	MDSSTATUS = MDSPUT('Z_XPTS'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,Z_XPTS,MXPTS,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing Z_XPTS'

	MDSSTATUS = MDSSETDEFAULT('\MGAMS.FBT'//CHAR(0))
	IF (MOD(MDSSTATUS,2).EQ.0) THEN
	 PRINT *,'FBTMDSPUT - error selecting FBT branch'
	 STOP
	ENDIF

	MDSSTATUS = MDSPUT('PLACAMP'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,PLACAMP_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing PLACAMP'
	MDSSTATUS = MDSPUT('RZE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,RZE_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing RZE'
	MDSSTATUS = MDSPUT('ZZE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,ZZE_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing ZZE'
	MDSSTATUS = MDSPUT('RAX'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,RAX_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing RAX'
	MDSSTATUS = MDSPUT('ZAX'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,ZAX_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing ZAX'
	MDSSTATUS = MDSPUT('RAX2'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,RAX2_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing RAX2'
	MDSSTATUS = MDSPUT('ZAX2'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,ZAX2_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing ZAX2'
	MDSSTATUS = MDSPUT('SHVSEC'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,SHVSEC_T,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing SHVSEC'
	MDSSTATUS = MDSPUT('GCURR'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,GCURR_T,NGROUP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing GCURR'
	MDSSTATUS = MDSPUT('BEE'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,BEE_T,NPROB,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing BEE'
	MDSSTATUS = MDSPUT('FLUX'//CHAR(0),'$1'//CHAR(0),DESCR(IDTYPE_FLOAT,FLUX_T,NLOOP,NRUNS,0),0)
	IF (MOD(MDSSTATUS,2).EQ.0) PRINT *,'FBTMDSPUT - error writing FLUX'

	MDSSTATUS = MDSCLOSE('PCS'//CHAR(0),SHOT)
	IF (MOD(MDSSTATUS,2).EQ.0) THEN
	 PRINT *,'FBTMDSPUT - error closing PCS tree'
	ENDIF
 
	RETURN
	END

