gray/src/gray.f

      implicit real*8 (a-h,o-z)
      common/istop/istop
      common/ierr/ierr
      common/igrad/igrad
      common/iovmin/iopmin,iowmin
      common/mode/sox
      common/p0/p0mw
      common/powrfl/powrfl
      common/index_rt/index_rt
      common/taumnx/taumn,taumx,pabstot,currtot
      common/ipass/ipass
      
c     read data plus initialization

      index_rt=1
      call prfile
      call paraminit
      call read_data
      call vectinit
      if(igrad.eq.0) call ic_rt
      if(igrad.gt.0) call ic_gb
      if(ierr.gt.0) then
        print*,'    IERR = ', ierr
        stop
      end if

c     beam/ray propagation
      call gray_integration

c     postprocessing
      
      call after_gray_integration
 
      pabstott=pabstot
      currtott=currtot
      powtr=p0mw-pabstot

      if (iowmin.eq.2.and.ipass.gt.1) then
c     second pass into plasma
        p0mw1=p0mw
        igrad=0

        index_rt=2
        p0mw=p0mw1*powrfl
        call prfile
        call vectinit2
        call paraminit
        call ic_rt2
        call gray_integration
        call after_gray_integration
        pabstott=pabstott+pabstot
        currtott=currtott+currtot
      
        index_rt=3
        sox=-sox
        p0mw=p0mw1*(1.0d0-powrfl)
        call prfile
        call vectinit2
        call paraminit
        call ic_rt2
        call gray_integration
        call after_gray_integration
        pabstott=pabstott+pabstot
        currtott=currtott+currtot
      end if
      print*,' '
      write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3

      end


      subroutine gray_integration
      implicit real*8 (a-h,o-z)

      common/ss/st
      common/dsds/dst
      common/istep/istep
      common/nstep/nstep
      common/istop/istop
      common/strfl11/strfl11
      common/index_rt/index_rt

c     ray integration: begin
      st0=0.0d0
      if(index_rt.gt.1) st0=strfl11
      do i=1,nstep
        istep=i
        st=i*dst+st0

c       advance one step 

        call rkint4

c       calculations after one step:

        call after_onestep(istep,istop)
        if(istop.eq.1) exit
c
      end do

c     ray integration: end

      return
      end


      subroutine after_gray_integration
      implicit real*8 (a-h,o-z)
      parameter(zero=0.0d0)
      character*24 filenmeqq,filenmprf,filenmbm
c
      common/ss/st
      common/ibeam/ibeam
      common/warm/iwarm,ilarm
      common/filesn/filenmeqq,filenmprf,filenmbm
      common/nray/nrayr,nrayth
      common/iieq/iequil
      common/iipr/iprof
      common/index_rt/index_rt
c 
      common/p0/p0mw
      common/factb/factb
      common/taumnx/taumn,taumx,pabstot,currtot
      common/scal/iscal
      common/facttn/factt,factn
c
c     print all ray positions in local reference system
c
      if(nrayr.gt.1) then
        iproj=1
        nfilp=9
        call projxyzt(iproj,nfilp)
      end if
c
c     print final results on screen
c
      print*,' '
      print'(a,f9.4)','final step (s, ct, Sr)  = ',st
      if(iwarm.gt.0) then
      print '(a,2e12.5)','taumn, taumx  = ', taumn,taumx
      else
      print '(a,2f9.4)','taumn, taumx  = ', zero,zero
      end if      
c
      print'(a,f9.4)','Pabs_tot (MW) = ',pabstot
      currtka =currtot*1.0d3
      print'(a,f9.4)','I_tot (kA)    = ',currtka
c
      if (index_rt.eq.1) then
      if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq
      if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM'
      if(iprof.eq.1) write(6,*)  'PROFILE file : ',filenmprf
      if(iprof.eq.0) write(6,*)  'ANALTYCAL PROFILES'
      if(ibeam.ge.1) write(6,*)  'LAUNCHER CASE : ',filenmbm

      end if
c
c       compute power and current density profiles for all rays
c
      pabs=pabstot
      currt=currtot
      call pec(pabs,currt)
c
      return
 99   format(20(1x,e16.8e3))
      end
c
c
c
      subroutine after_onestep(i,istop)
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36,nmx=8000)
      parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0)
      dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
      dimension iop(jmx,kmx),iow(jmx,kmx)
      dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
      dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
      
      common/pcjki/ppabs,ccci
      common/atjki/tauv,alphav
      common/tau1v/tau1v
c
      common/warm/iwarm,ilarm
      common/nray/nrayr,nrayth
      common/ist/istpr0,istpl0
      common/istgr/istpr,istpl
c
      common/iiv/iiv
      common/iov/iop,iow
      common/psjki/psjki
      common/psival/psinv
      common/psinv11/psinv11
      common/ierr/ierr
      common/taumnx/taumn,taumx,pabstot,currtot
      common/xv/xv
      common/anv/anv
      common/cent/btrcen,rcen
c
      common/p0/p0mw
      common/pol0/psipol0,chipol0
      common/ipol/ipolc
      common/iovmin/iopmin,iowmin
      common/densbnd/psdbnd
      common/yyrfl/yyrfl
      common/powrfl/powrfl
      common/dstvac/dstvac
      common/strfl11/strfl11
      common/dsds/dst
      common/index_rt/index_rt
      common/ipass/ipass
      common/rwallm/rwallm
      common/bound/zbmin,zbmax

      pabstot=0.0d0
      currtot=0.0d0
      taumn=1d+30
      taumx=-1d+30
      psinv11=1.0d0
      iopmin=100
      iowmin=100
c
      do j=1,nrayr
        kkk=nrayth
        if(j.eq.1) kkk=1
        do k=1,kkk
          call gwork(j,k)
c 
          if(ierr.gt.0) then
!            print*,'    IERR = ', ierr
            if(ierr.eq.97) then
c              igrad=0
c              ierr=0 
            else  
              istop=1
              exit
            end if
          end if

          psjki(j,k,i)=psinv
          rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0
          zzm=xv(3)/100.d0
          if(j.eq.1) rrm11=rrm

          if (iwarm.gt.0.and.i.gt.1) then
            if(psinv.ge.0.and.psinv.le.1.0d0.and.
     .         zzm.ge.zbmin.and.zzm.le.zbmax) then
              call pabs_curr(i,j,k)
              iiv(j,k)=i
            else
              if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1)
            end if
            if(tauv(j,k,i).le.taumn) taumn=tauv(j,k,i)
            if(tauv(j,k,i).ge.taumx) taumx=tauv(j,k,i)
            pabstot=pabstot+ppabs(j,k,iiv(j,k))
            currtot=currtot+ccci(j,k,iiv(j,k))
          end if 
          call print_output(i,j,k)

          if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) 
     .     iop(j,k)=1
          if(iop(j,k).eq.1.and.
     .       (psinv.ge.psdbnd.or.
     .        (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax))))
     .     iop(j,k)=2
c         iov=0 initially, iov=1 first entrance in plasma, 
c         iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
          
          if(index_rt.eq.1) then
          if(j.eq.1) then
            psinv11=psinv
            if(ipolc.eq.0.and.iop(j,k).eq.1) then
              call pol_limit(qqin,uuin,vvin)
              ipolc=1
               qqa=cos(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
               uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
               vva=sin(2.0d0*chipol0*cvdr)
               powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin)
c              p0mw=p0mw*powa
c              print*,' '
c              print*,'Coupled power fraction =',powa
c              print*,' '
c              print*,'Input coupled power (MW) =',p0mw
c              print*,' '
            end if
            if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1
     .          .and.ipolc.eq.1) then
              call pol_limit(qqout,uuout,vvout)
              ipolc=2
              call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
              strfl11=i*dst+dstvac
              call pol_limit(qqin2,uuin2,vvin2)
              if(irfl.gt.0) then
                powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2)
              else
                powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2)
              end if
              write(6,*) 'Reflected power fraction =',powrfl
              iop(j,k)=3
              yyrfl(j,k,1)=xvrfl(1)
              yyrfl(j,k,2)=xvrfl(2)
              yyrfl(j,k,3)=xvrfl(3)
              yyrfl(j,k,4)=anvrfl(1)
              yyrfl(j,k,5)=anvrfl(2)
              yyrfl(j,k,6)=anvrfl(3)
              tau1v(j,k)=tauv(j,k,iiv(j,k))
            end if
          else
            if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1)  then
              call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
              iop(j,k)=3
              yyrfl(j,k,1)=xvrfl(1)
              yyrfl(j,k,2)=xvrfl(2)
              yyrfl(j,k,3)=xvrfl(3)
              yyrfl(j,k,4)=anvrfl(1)
              yyrfl(j,k,5)=anvrfl(2)
              yyrfl(j,k,6)=anvrfl(3)
              tau1v(j,k)=tauv(j,k,iiv(j,k))
            end if  
          end if   
          end if

          if(iop(j,k).lt.iopmin) iopmin=iop(j,k)

        end do
      end do

      if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1

      psimin=psjki(1,1,i)
      if(nrayr.gt.1)
     .    psimin=min(psimin,minval(psjki(2:nrayr,1:nrayth,i)))
      if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
     .    istop=1

      if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1
      if(iopmin.eq.2.and.ipass.eq.1) istop=1
      
      if(iopmin.eq.3) istop=1
     
c     print ray positions for j=nrayr in local reference system

      istpr=istpr+1
      if (istpr.eq.istpr0) then
c        print*,'istep = ',i
        if(nrayr.gt.1) then
          iproj=0
          nfilp=8
          call projxyzt(iproj,nfilp)
        end if
        istpr=0
      end if
c
      if (istpl.eq.istpl0) istpl=0
      istpl=istpl+1
      return
      end
c
c
c
      subroutine print_output(i,j,k)
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,jmx=31,kmx=36,nmx=8000)
      parameter(taucr=12.0d0)
      parameter(pi=3.14159265358979d0)
      dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
      dimension psjki(jmx,kmx,nmx)
      dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
      dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension q(jmx),tau1v(jmx,kmx)
      complex*16 ex,ey,ez

c
      common/psjki/psjki
      common/atjki/tauv,alphav
      common/dpjjki/pdjki,currj
      common/pcjki/ppabs,ccci
      common/dcjki/didst
      common/nhn/nhn
      common/iokh/iohkawa
      common/p0/p0mw
      common/tau1v/tau1v
      common/qw/q
      common/index_rt/index_rt
      
      common/ss/st
      common/nray/nrayr,nrayth
      common/istgr/istpr,istpl
      common/ist/istpr0,istpl0
      common/iieq/iequil
c
      common/parpl/brr,bphi,bzz,ajphi
      common/btot/btot
      common/xgxg/xg
      common/ygyg/yg
      common/dens/dens,ddens
      common/tete/tekev
      common/absor/alpha,effjcd,akim,tau0
      common/densbnd/psdbnd
      common/epolar/ex,ey,ez
c
      common/nplr/anpl,anpr
      common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
      common/nprw/anprr,anpri
c
      common/wrk/ywrk,ypwrk
c
      x=ywrk(1,j,k)
      y=ywrk(2,j,k)
      z=ywrk(3,j,k)
      rr=sqrt(x*x+y*y)
c
      anx=ywrk(4,j,k)
      any=ywrk(5,j,k)
      anz=ywrk(6,j,k)
      anr=(anx*x+any*y)/rr
      anphi=(any*x-anx*y)/rr
c      cnphi=(any*x-anx*y)
c
      rrm=rr*1.0d-2
      zzm=z*1.0d-2
      stm=st*1.0d-2
      xxm=x*1.0d-2
      yym=y*1.0d-2
c
      if(index_rt.gt.1) then
        taujk=tauv(j,k,i)+tau1v(j,k)
      else
        taujk=tauv(j,k,i)
      end if

c     central ray only begin
      if(j.eq.1) then
        phi=acos(x/sqrt(y*y+x*x))
        if(y.lt.0.0d0) phi=-phi
        phideg=phi*180.0d0/pi
        psi=psjki(j,k,i)
        rhot=1.0d0
        bbr=0.0d0
        bbz=0.0d0
        bpol=0.0d0
        rhot=1.0d0
        dens11=0.0d0
        if(psi.ge.0.0d0) then
          if(iequil.eq.2) then
            if (psi.le.1.0d0) rhot=frhotor(psi)
            bbr=brr
            bbz=bzz
            bpol=sqrt(bbr**2+bbz**2)
          else
            rhot=sqrt(psi)
          end if
        else
          tekev=0.0d0
          akim=0.0d0
        end if
        ddr11=ddr
c        cutoff=xg-(1-yg)*(1-anpl**2)        

c     print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1

        if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens
        dids11=didst(j,k,i)*1.0d2/(p0mw*q(j))
        dpdv11=pdjki(j,k,i)/(p0mw*q(j))
        ajcd11=currj(j,k,i)/(p0mw*q(j))
        alpha11=alphav(j,k,i)
        tau11=taujk
        pt11=exp(-taujk)

        write(4,99) stm,rrm,zzm,phideg,
     .  psi,rhot,dens11,tekev,
     .  bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100,
     .  alpha11,tau11,pt11,ajcd11,dids11,
     .  dble(nhn),dble(iohkawa),dble(index_rt)

c        call polarcold(exf,eyif,ezf,elf,etf)
      end if
c     central ray only end

      if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi

c     print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
      if(istpl.eq.istpl0)  then
        if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
          write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
     .          psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
        end if
c      if(k.eq.nrayth) write(33,*) ' '
      end if
c
      return
 99   format(30(1x,e16.8e3))
111   format(3i5,16(1x,e16.8e3))
      end
c
c
c
      subroutine prfile
      implicit none
      integer*4 index_rt
      common/index_rt/index_rt

      If(index_rt.eq.1) then
      
      write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '//
     .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
      write(8,*) ' #istep j k xt yt zt rt psin'
      write(9,*) ' #istep j k xt yt zt rt psin'
      write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
      write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
      write(12,*) ' #i sst psi w1 w2'
      write(7,*)'#Icd Pa Jphimx dPdVmx '//
     .'rhotj rhotjava rhotp rhotpav '//
     .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
      write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'

      else

c      If(index_rt.eq.3) then
      write(4,*) ' '
      write(8,*) ' '
      write(9,*) ' '
      write(17,*) ' '
      write(12,*) ' '
      write(48,*) ' '

c      end if
      end if
      return
      end
c
c
c
      subroutine read_data
      use green_func_p, only:Setup_SpitzFunc
      implicit real*8 (a-h,o-z)
      real*8 me
      character*8 wdat
      character*10 wtim
      character*24 filenmeqq,filenmprf,filenmbm
      parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
      parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
      parameter(nmx=8000,nbb=1000)
      real*8 rlim(nbb),zlim(nbb)
c
      common/xgcn/xgcn

      common/ipec/ipec,nnd
      common/nstep/nstep
      common/ibeam/ibeam
      common/ist/istpr0,istpl0
      common/warm/iwarm,ilarm
      common/ieccd/ieccd
      common/idst/idst
c
      common/filesn/filenmeqq,filenmprf,filenmbm
c
      common/nray/nrayr,nrayth
      common/rwmax/rwmax
      common/dsds/dst
      common/igrad/igrad
      common/ipass/ipass
      common/rwallm/rwallm
      common/limiter/rlim,zlim,nlim
      common/iieq/iequil
      common/icocos/icocos
      common/ixp/ixp
      common/ipsn/ipsinorm
      common/sspl/sspl
      common/iipr/iprof
      common/factb/factb
      common/facttn/factt,factn
      common/cent/btrcen,rcen
c
      common/parwv/ak0,akinv,fhz
      common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
      common/anic/anx0c,any0c,anz0c
      common/mirr/x00,y00,z00
      common/pol0/psipol0,chipol0
c
      common/pardens/dens0,aln1,aln2
      common/parban/b0,rr0m,zr0m,rpam
      common/parqq/q0,qa,alq
      common/parqte/te0,dte0,alt1,alt2
      common/zz/Zeff
c
      common/parbres/bres
      common/densbnd/psdbnd
      common/nfile/neqdsk,nprof
      common/sgnib/sgnbphi,sgniphi
      common/p0/p0mw
c
      common/mode/sox
      common/angles/alpha0,beta0
      common/scal/iscal
c
      open(2,file='gray.data',status= 'unknown')
c
c     alpha0, beta0 (cartesian) launching angles
c     fghz      wave frequency (GHz)
c     p0mw      injected power (MW)
c
      read(2,*) alpha0,beta0
      read(2,*) fghz
      read(2,*) p0mw
c
c     nrayr     number of rays in radial direction
c     nrayth      number of rays in angular direction
c     rwmax      normalized maximum radius of beam power
c     rwmax=1 -> last ray at radius = waist
c
      read(2,*) nrayr,nrayth,rwmax
      if(nrayr.eq.1) nrayth=1
c
c     x00,y00,z00 coordinates of launching point
c
      read(2,*) x00,y00,z00
c
c     beams parameters in local reference system
c     w0 -> waist, d0 -> waist distance from launching point
c     phiw  angle of beam ellipse
c
      read(2,*) w0csi,w0eta,d0csi,d0eta,phiw
c
c     ibeam=0  :read data for beam as above
c     ibeam=1  :read data from file simple astigmatic beam
c     ibeam=2  :read data from file general astigmatic beam
      read(2,*) ibeam
      read(2,*) filenmbm
c
c     iequil=0  :vacuum
c     iequil=1  :analytical equilibrium
c     iequil=2  :read eqdsk
c     ixp=0,-1,+1  : no X point , bottom/up X point
c
      read(2,*) iequil,ixp
c
c     iprof=0  :analytical density and temp. profiles
c     iprof>0  :numerical density and temp. profiles
c
      read(2,*) iprof
c
c     iwarm=0  :no absorption and cd
c     iwarm=1  :weakly relativistic absorption
c     iwarm=2  :relativistic absorption, n<1 asymptotic expansion
c     iwarm=3  :relativistic absorption, numerical integration
c     ilarm    :order of larmor expansion
c
      read(2,*) iwarm,ilarm
c      if(iwarm.gt.2) iwarm=2
c
c     ieccd 0/1 NO/YES  ECCD calculation ieccd>0 different CD models
c
      read(2,*) ieccd
      if (ieccd.ge.10)     call Setup_SpitzFunc
c
c     ipec=0/1  :pec profiles grid in psi/rhop
c     nnd       :number of grid steps for pec profiles +1
c
      read(2,*) ipec,nnd
c
c     dst       integration step
c     nstep      maximum number of integration steps
c     istpr0     projection step = dsdt*istprj
c     istpl0      plot step = dsdt*istpl
c     igrad=0    optical ray-tracing, initial conditions as for beam
c     igrad=1    quasi-optical ray-tracing
c     igrad=-1   ray-tracing, init. condit.
c                from center of mirror and with angular spread
c     ipass=1/2 1 or 2 passes into plasma
c     iox=1/2    OM/XM
c     idst=0/1/2  0 integration in s, 1 integr. in ct, 2 integr. in Sr
c
      read(2,*) dst,nstep,istpr0,istpl0,idst
      read(2,*) igrad,ipass,rwallm
      read(2,*) iox, psipol0,chipol0
c
c     ipsinorm  0 standard EQDSK format, 1 format Portone summer 2004
c     sspl      spline parameter for psi interpolation
c
      read(2,*) filenmeqq
      read(2,*) ipsinorm,sspl
c
c     factb     factor for magnetic field (only for numerical equil)
c     scaling adopted: beta=const, qpsi=const, nustar=const
c     iscal    ne Te scaling 0: nustar=const,  1: n_greenw=const; 2 no rescaling
c     factT factn   factor for Te&ne scaling 
c
      read(2,*) factb, iscal,factt,factn
      if (factb.eq.0.0d0) factb=1.0d0
c
c     psbnd     value of psi ( > 1 ) of density boundary
c
      read(2,*) filenmprf
      read(2,*) psdbnd
      if(iprof.eq.0) psdbnd=1.0d0
c
c     signum of toroidal B and I
c     icocos index for equilibrium from COCOS - O. Sauter Feb 2012
      read(2,*) sgnbphi,sgniphi,icocos
c
c     read parameters for analytical density and temperature
c     profiles if iprof=0
c
      if(iprof.eq.0) then
        read(2,*) dens0,aln1,aln2
        read(2,*) te0,dte0,alt1,alt2
      else
        read(2,*) dummy,dummy,dummy
        read(2,*) dummy,dummy,dummy,dummy
      end if
      read(2,*) zeff
c
c     read parameters for analytical simple cylindical equilibrium
c     if iequil=0,1

      if(iequil.lt.2) then
        read(2,*) rr0,zr0,rpa
        read(2,*) b0
        read(2,*) q0,qa,alq
        rr0m=rr0/1.0d2
        zr0m=zr0/1.0d2
        rpam=rpa/1.0d2
      else
        read(2,*) dummy,dummy,dummy
        read(2,*) dummy
        read(2,*) dummy,dummy,dummy
      end if
c
      close(unit=2)
c
      if(nrayr.eq.1) igrad=0
      if (nrayr.lt.5) then
        igrad=0
        print*,'    nrayr < 5 ! => OPTICAL CASE ONLY'
        print*,' '
      end if
c
      fhz=fghz*1.0d9
      ak0=2.0d9*pi*fghz/vc
      akinv=1.0d0/ak0
c
      bresg=2.0d0*pi*fhz*me*vc/qe
      bres=bresg*1.0d-4
c
c     xg=xgcn*dens19
c
      xgcn=1.0d-5*qe**2/(pi*me*fghz**2)
c
      sox=-1.0d0
      if(iox.eq.2) sox=1.0d0
c
c     read data for beam from file if ibeam>0
c
      if(ibeam.gt.0) then
        call read_beams
      else
        zrcsi=0.5d0*ak0*w0csi**2
        zreta=0.5d0*ak0*w0eta**2
        if(igrad.gt.0) then
          wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2)
          weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2)
          rcicsi=-d0csi/(d0csi**2+zrcsi**2)
          rcieta=-d0eta/(d0eta**2+zreta**2)
          phir=phiw
        else
          d0eta=d0csi
          wcsi=w0csi*abs(d0csi/zrcsi)
          weta=w0eta*abs(d0eta/zreta)
          rcicsi=w0csi/zrcsi
          rcieta=w0eta/zreta
          if(d0csi.gt.0) then
            rcicsi=-rcicsi
            rcieta=-rcieta
          end if
          phir=phiw
        end if
      end if
c
      print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
c
      r00=sqrt(x00**2+y00**2)
      phi0=acos(x00/r00)
      if(y00.lt.0) phi0=-phi0
      print'(a,4f8.3)','x00, y00, R00, z00 = ',x00,y00,r00,z00
      print*,' '
c
c      anx0c=(anr0c*x00-anphi0c*y00)/r00
c      any0c=(anr0c*y00+anphi0c*x00)/r00
c
c      anr0c=(anx0c*x00+any0c*y00)/r00
c      anphi0c=(any0c*x00-anx0c*y00)/r00
c
c     angles alpha0, beta0 in a local reference system as proposed by Gribov et al
c
c      anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
c      anphi0c=sin(cvdr*beta0)
c      anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)

      anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0)
      anphi0c=sin(cvdr*beta0)
      anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
 
      anx0c=(anr0c*x00-anphi0c*y00)/r00
      any0c=(anr0c*y00+anphi0c*x00)/r00
c
c     read data for Te , ne , Zeff from file if iprof>0
c

      if (iprof.eq.1) then
        nprof=98
        lprf=length(filenmprf)
        open(file=filenmprf(1:lprf)//'.prf',
     .                    status= 'unknown',unit=nprof)
        call profdata
        close(nprof)
      end if
c
c     read equilibrium data  from file if iequil=2
c
      if (iequil.eq.2) then
        neqdsk=99
        leqq=length(filenmeqq)
         open(file=filenmeqq(1:leqq)//'.eqdsk',
     .        status= 'unknown', unit=neqdsk)
        call equidata
        close(neqdsk)
        
c     print density, temperature, safecty factor, toroidal current dens
c     versus psi, rhop, rhot
        call print_prof
      end if

      if (iequil.eq.1)   call surf_anal

      if (iequil.ne.2.or.ipass.lt.0) then
c     set simple limiter as two cylindrical walls at rwallm and r00
        nlim=5
        rlim(1)=rwallm
        rlim(2)=r00*1.d-2
        rlim(3)=rlim(2)
        rlim(4)=rlim(1)
        rlim(5)=rlim(1)
        zlim(1)=-dst*nmx*1.d-2
        zlim(2)=zlim(1)
        zlim(3)=dst*nmx*1.d-2
        zlim(4)=zlim(3)
        zlim(5)=zlim(1)
        ipass=abs(ipass)
      end if

      nfil=78
      
      open(nfil,file='headers.txt',status= 'unknown')
      call date_and_time(wdat,wtim)
      write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8),
     .                wtim(1:2),wtim(3:4),wtim(5:6)
      write(nfil,904) REVISION
      if (iequil.eq.2) then
        write(nfil,905) filenmeqq
      else
        if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq
        if (iequil.eq.0) write(nfil,'("# VACUUM CASE")') 
      end if    
      if (iprof.eq.1) then
        write(nfil,907) filenmprf
      else
        write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff
      end if    
      write(nfil,911) fghz,p0mw,iox
      write(nfil,902) x00,y00,z00      
      write(nfil,908) alpha0,beta0   
      if(ibeam.ge.1) write(nfil,909) filenmbm
      if(ibeam.eq.0) write(nfil,903) w0csi,w0eta,d0csi,d0eta,phiw
      write(nfil,900) nrayr, nrayth, rwmax
      write(nfil,901) igrad,iwarm,ilarm,ieccd,idst
      write(nfil,915) dst,nstep
      write(nfil,914) sgnbphi,sgniphi,icocos
      write(nfil,906) factb,factt,factn,iscal
      write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm
      write(nfil,*) ' '
      close(nfil) 
      
      return
      
900   format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5)      
901   format('# igrad iwarm ilarm ieccd idst : ',6i5)      
902   format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5))
903   format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5))
904   format('# GRAY revision : ',a)
905   format('# EQUILIBRIUM file : ',a24)      
906   format('# fact_B  fact_T fact_n iscal : ',(3(1x,es12.5),i5))
907   format('# PROFILES file : ',a24)      
908   format('# alpha0 beta0 launch angles (deg) CYL : ',2(1x,es12.5))
909   format('# LAUNCHER file : ',a24)      
910   format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5)      
911   format('# fghz P0 IOX : ',2(1x,es12.5),i5) 
912   format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5))
913   format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : '
     . ,8(1x,es12.5))
914   format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5)      
915   format('# dst nstep : ',1x,es12.5,i5)
916   format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2))

 99   format(20(1x,e16.8e3))
      end
c
c
c
      subroutine surf_anal
      implicit real*8(a-h,o-z)
      parameter(pi=3.14159265358979d0)
      common/parban/b0,rr0m,zr0m,rpam
      common/parbres/bres
c
c     print circular magnetic surfaces iequil=1
c
      write(71,*) '#i psin R z'
      dal=2.0d-2*pi
      drn=0.1d0
      do i=1,10
        rni=i*drn
        do k=1,101
          drrm=rpam*rni*cos((k-1)*dal)
          dzzm=rpam*rni*sin((k-1)*dal)
          rrm=rr0m+drrm
          zzm=zr0m+dzzm
          write(71,111) i,rni,rrm,zzm
        end do
        write(71,*) ' '
        write(71,*) ' '
      end do
c
c     print resonant B iequil=1
c
      write(70,*)'#i Btot  R  z'
      rres=b0*rr0m/bres
      zmx=zr0m+rpam
      zmn=zr0m-rpam
      do i=1,21
          zzres=zmn+(i-1)*(zmx-zmn)/2.0d1
          write(70,111) i,bres,rres,zzres
      end do

      return
111   format(i4,20(1x,e16.8e3))
      end
c
c
      subroutine read_beams
      implicit real*8(a-h,o-z)
      character*24 filenmeqq,filenmprf,filenmbm
      parameter(nstrmx=50)
c
      dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4)
      dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx)
      dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4)
      dimension waist1v(nstrmx),waist2v(nstrmx)
      dimension rci1v(nstrmx),rci2v(nstrmx)
      dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4)
      dimension crci1(nstrmx,4),crci2(nstrmx,4)
      dimension phi1v(nstrmx),phi2v(nstrmx)
      dimension cphi1(nstrmx,4),cphi2(nstrmx,4)
c
      common/ibeam/ibeam
      common/filesn/filenmeqq,filenmprf,filenmbm
      common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
      common/mirr/x00,y00,z00
      common/angles/alpha0,beta0
      common/parwv/ak0,akinv,fhz
c
c     for given alpha0 -> beta0 + beam parameters
c
c     ibeam=1 simple astigmatic beam
c     ibeam=2 general astigmatic beam
c
c     initial beam data are measured in mm -> transformed to cm
c
      nfbeam=97
      lbm=length(filenmbm)
      open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam)
c
      read(nfbeam,*) nisteer
      do i=1,nisteer
        if(ibeam.eq.1) then
          read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm,
     .                   waist1,dvir1,waist2,dvir2,delta
          phi1=delta
          phi2=delta
          zr1=0.5d-1*ak0*waist1**2
          zr2=0.5d-1*ak0*waist2**2
          w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2)
          w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2)
          rci1=-dvir1/(dvir1**2+zr1**2)
          rci2=-dvir2/(dvir2**2+zr2**2)
        else
          read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm,
     .                   w1,w2,rci1,rci2,phi1,phi2
        end if
        x00v(i)=0.1d0*x00mm
        y00v(i)=0.1d0*y00mm
        z00v(i)=0.1d0*z00mm
        alphastv(i)=alphast
        betastv(i)=betast
        waist1v(i)=0.1d0*w1
        rci1v(i)=1.0d1*rci1
        waist2v(i)=0.1d0*w2
        rci2v(i)=1.0d1*rci2
        phi1v(i)=phi1
        phi2v(i)=phi2
      end do
c
      iopt=0
      call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier)
      call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier)
      call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier)
      call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier)
      call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier)
      call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier)
      call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier)
      call difcs(alphastv,x00v,nisteer,iopt,cx0,ier)
      call difcs(alphastv,y00v,nisteer,iopt,cy0,ier)
      call difcs(alphastv,z00v,nisteer,iopt,cz0,ier)
c
      if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
        call locate(alphastv,nisteer,alpha0,k)
        dal=alpha0-alphastv(k)
        betst=spli(cbeta,nisteer,k,dal)
        x00=spli(cx0,nisteer,k,dal)
        y00=spli(cy0,nisteer,k,dal)
        z00=spli(cz0,nisteer,k,dal)
        wcsi=spli(cwaist1,nisteer,k,dal)
        weta=spli(cwaist2,nisteer,k,dal)
        rcicsi=spli(crci1,nisteer,k,dal)
        rcieta=spli(crci2,nisteer,k,dal)
        phiw=spli(cphi1,nisteer,k,dal)
        phir=spli(cphi2,nisteer,k,dal)
      else
        print*,'    alpha0 outside table range !!!'
        if(alpha0.ge.alphastv(nisteer)) ii=nisteer
        if(alpha0.le.alphastv(1)) ii=1
        betst=betastv(ii)
        x00=x00v(ii)
        y00=y00v(ii)
        z00=z00v(ii)
        wcsi=waist1v(ii)
        weta=waist2v(ii)
        rcicsi=rci1v(ii)
        rcieta=rci2v(ii)
        phiw=phi1v(ii)
        phir=phi2v(ii)
      end if
      beta0=betst
c
      close(nfbeam)
c
      return
      end
c
c
c
      subroutine equidata
      implicit real*8 (a-h,o-z)
      parameter(nnw=501,nnh=501)
      parameter(pi=3.14159265358979d0)
      parameter(nbb=1000)
      character*48 stringa
      dimension fpol(nnw),pres(nnw),qpsi(nnw)
      dimension ffprim(nnw),pprim(nnw)
      dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw)
      dimension rbbbs(nbb),zbbbs(nbb)
      dimension rlim(nbb),zlim(nbb)
c
      parameter(nrest=nnw+4,nzest=nnh+4)
      parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54)
      parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3)
      dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh)
      dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk)
      parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh)
      dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
      parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh)
      parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh)
      parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh)
      dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11)
      dimension iwrkd(ldiwrk)
      parameter(lwrkf=nnw*4+nrest*16)
      dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw)
      dimension fpoli(nnw)
c
      common/pareq1/psia
      common/pareq1a/psiaxis0
      common/cent/btrcen,rcen
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      common/psinr/psinr
      common/qpsi/qpsi
      common/cpsin/rv,zv,psin
      common/cpsi/psi
      common/eqnn/nr,nz,npp,nintp
      common/ipsn/ipsinorm
      common/sspl/sspl
      common/nfile/neqdsk,nprof
      common/bound/zbmin,zbmax
      common/sgnib/sgnbphi,sgniphi
      common/factb/factb
      common/ixp/ixp
      common/icocos/icocos
      
      common/coffeqt/tr,tz
      common/coffeqtp/tfp
      common/coffeq/cc
      common/coffeqd/cc01,cc10,cc20,cc02,cc11
      common/coffeqn/nsrt,nszt,nsft
      common/cofffp/cfp
      common/fpas/fpolas
      common/rhotsx/rhotsx
      common/rrtor/rrtor
      common/cnt/rup,zup,rlw,zlw
      common/limiter/rlim,zlim,nlim
c
c     read from file eqdsk
c     (see http://fusion.gat.com/efit/g_eqdsk.html)
c
c     spline interpolation of psi and derivatives
c
      if(icocos.gt.0) then 
      read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz
      else
      read (neqdsk,*) nr,nz
      end if
      if(ipsinorm.eq.0) then
        read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid
        read (neqdsk,2020) rmaxis,zmaxis,psiax,psiedge,btrcen
        read (neqdsk,2020) current,simag,xdum,rmaxis,xdum
        read (neqdsk,2020) zmaxis,xdum,sibry,xdum,xdum
        read (neqdsk,2020) (fpol(i),i=1,nr)
        read (neqdsk,2020) (pres(i),i=1,nr)
        read (neqdsk,2020) (ffprim(i),i=1,nr)
        read (neqdsk,2020) (pprim(i),i=1,nr)
        read (neqdsk,2020) ((psi(i,j),i=1,nr),j=1,nz)
        read (neqdsk,2020) (qpsi(i),i=1,nr)
      else
        read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid
        read (neqdsk,*) rmaxis,zmaxis,psiax,psiedge,btrcen
        read (neqdsk,*) current,simag,xdum,rmaxis,xdum
        read (neqdsk,*) zmaxis,xdum,sibry,xdum,xdum
        read (neqdsk,*) (fpol(i),i=1,nr)
        read (neqdsk,*) (pres(i),i=1,nr)
        read (neqdsk,*) (ffprim(i),i=1,nr)
        read (neqdsk,*) (pprim(i),i=1,nr)
        read (neqdsk,*) ((psin(i,j),i=1,nr),j=1,nz)
        read (neqdsk,*) (qpsi(i),i=1,nr)
      end if
2020  format (5e16.9)

c
c  compensate for different reference systems
c
      icocosmod=mod(icocos,10)

      if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
c  icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
        btrcen=-btrcen
        current=-current
        do i=1,nr
          fpol(i)=-fpol(i)
        end do
      end if
c
      if (icocosmod.eq.1 .or. icocosmod.eq.4 .or. 
     &    icocosmod.eq.5 .or. icocosmod.eq.8) then
c  icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
c  icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
        psiedge=-psiedge
        psiax=-psiax
        if (ipsinorm.eq.0) then
          do j=1,nz
            do i=1,nr
              psi(i,j)=-psi(i,j)
            end do
          end do
        end if
      end if

c
c  add check for Ip/psi and B0/Fpol sign consistency?
c
      current=sign(current,psiax-psiedge)
      btrcen=sign(btrcen,fpol(nr))

c
c     length in m !!!
c
      dr=drnr1/dble(nr-1)
      dz=dznz1/dble(nz-1)
      rv(1)=rleft
      zv(1)=zmid-dznz1/2.0d0
      dpsinr=1.0d0/dble(nr-1)
c
      do i=1,nr
        psinr(i)=(i-1)*dpsinr
        rv(i)=rv(1)+(i-1)*dr
      end do
c
      do j=1,nz
        zv(j)=zv(1)+(j-1)*dz
      end do
c
      rmnm=rv(1)
      rmxm=rv(nr)
      zmnm=zv(1)
      zmxm=zv(nz)

c     psi function

      psia0=psiedge-psiax
c  icocos=0: adapt psi to force Ip sign, otherwise maintain psi
      if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0)
      current=sign(current,sgniphi)

      psia=-sgniphi*abs(psia0)*factb
c  icocos>10: input psi is in Wb
c  icocos<10: input psi is in Wb/rad (gray convention)
      if (icocos.ge.10) psia=psia/(2.0d0*pi)

      psiaxis0=0.0d0
      do j=1,nz
        do i=1,nr
          if(ipsinorm.eq.0) then
            psin(i,j)=(psi(i,j)-psiax)/(psia0)
            psi(i,j)=psin(i,j)*psia
          else
            psi(i,j)=psin(i,j)*psia
          end if
          ij=nz*(i-1)+j
          fvpsi(ij)=psin(i,j)
        enddo
      enddo
c
c     spline interpolation of psi(i,j) and derivatives
c
      iopt=0
      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
     .              kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
     .              wrk,lwrk,iwrk,liwrk,ier)
      if(ier.eq.-1) then
      sspl=0.0d0
      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,
     .              kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp,
     .              wrk,lwrk,iwrk,liwrk,ier)
      end if
      nsrt=nsr
      nszt=nsz
      if (sspl.gt.0.0d0) then
        call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
     .   wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
c
        do j=1,nz
          do i=1,nr
            ij=nz*(i-1)+j
            psin(i,j)=ffvpsi(ij)
            psi(i,j)=psin(i,j)*psia
c            write(79,2021) rv(i),zv(j),psin(i,j)
          enddo
c            write(79,*) ' '
        enddo
      end if
2021  format(5(1x,e16.9))
c
      nur=1
      nuz=0
      call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
     . ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier)
c
      nur=0
      nuz=1
      call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
     . ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier)
c
      nur=2
      nuz=0
      call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
     . ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier)
c
      nur=0
      nuz=2
      call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
     . ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier)
c
      nur=1
      nuz=1
      call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz,
     . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier)
c
      read (neqdsk,*) nbbbs,nlim
      if(nbbbs.gt.0) then
        if(ipsinorm.eq.1)
     .    read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
        if(ipsinorm.eq.0)
     .    read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
      end if
      if(nlim.gt.0) then
        if(ipsinorm.eq.1)
     .    read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim)
        if(ipsinorm.eq.0)
     .    read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim)
      end if
c
c     compute max and min z of last closed surface
c
      rbmin=rax
      rbmax=rax
      if (nbbbs.gt.1) then
        zbmin=1.0d+30
        zbmax=-1.0d+30
        do i=1,nbbbs
          if(zbbbs(i).le.zbmin) then
            zbmin=zbbbs(i)
            rbmin=rbbbs(i)
          end if
          if(zbbbs(i).ge.zbmax) then
            zbmax=zbbbs(i)
            rbmax=rbbbs(i)
          end if
        end do
      else
        zbmin=-1.0d+30
        zbmax=1.0d+30
      end if
      if(zbmin.le.zmnm) zbmin=zbmin+dz
      if(rbmin.le.rmnm) rbmin=rbmin+dr
      if(zbmax.ge.zmxm) zbmax=zbmax-dz
      if(rbmax.ge.rmxm) rbmax=rbmax-dr
c
c     scaling of f_poloidal
c
c  icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol
      if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr))
      btrcen=sign(btrcen,sgnbphi)

      do i=1,nr
        fpol(i)=sgnbphi*abs(fpol(i))*factb
        wf(i)=1.0d0
      end do
      wf(nr)=1.0d2

c
c     spline interpolation of fpol(i)
c
      iopt=0
      xb=0.0d0
      xe=1.0d0
      ssfp=0.01d0
      call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft,
     . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier)
c
      call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier)
      fpolas=fpoli(nr)
c
c     search for O-point
c
      call points_ox(rmaxis,zmaxis,rmop,zmop,psinop,info)
      rmaxis=rmop
      zmaxis=zmop
      print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinop
c
c     search for X-point   if ixp not = 0
c
      if(ixp.ne.0) then
        if(ixp.lt.0) then
          r10=rbmin
          z10=zbmin
          call points_ox(r10,z10,rxp,zxp,psinxp,info)
          if(psinxp.ne.-1.0d0) then
            print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxp
            rbmin=rxp
            zbmin=zxp
            delpsinox=(psinxp-psinop)
            psia=psia*delpsinox
            deltapsi=abs(psia)
            psiaxis0=psia*psinop
            psin1=1.0d0
            r10=rmaxis
            z10=(zbmax+zmaxis)/2.0d0
            call points_tgo(r10,z10,r1,z1,psin1,info)
            rbmax=r1
            zbmax=z1
          else
            ixp=0
c            print'(a)','no X-point'
           end if
        else
          r10=rmop
          z10=zbmax
          call points_ox(r10,z10,rxp,zxp,psinxp,info)
          if(psinxp.ne.-1.0d0) then
            print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxp
            zbmax=zxp
            delpsinox=(psinxp-psinop)
            psia=psia*delpsinox
            deltapsi=abs(psia)
            psiaxis0=psia*psinop
            psin1=1.0d0
            z10=(zbmin+zmaxis)/2.0d0
            call points_tgo(r10,z10,r1,z1,psin1,info)
            zbmin=z1
          else
            ixp=0
c           print'(a)','no X-point'
           end if
        end if
      end if
c
      if (ixp.eq.0) then
        psin1=1.0d0
        delpsinox=(psin1-psinop)
        psia=psia*delpsinox
        deltapsi=abs(psia)
        psiaxis0=psia*psinop
        r10=rmaxis
        z10=(zbmax+zmaxis)/2.0d0
        call points_tgo(r10,z10,r1,z1,psin1,info)
        zbmax=z1
        rbmax=r1
        z10=(zbmin+zmaxis)/2.0d0
        call points_tgo(r10,z10,r1,z1,psin1,info)
        zbmin=z1
        rbmin=r1
        print'(a,4f8.4)','no X-point  ',rbmin,zbmin,rbmax,zbmax
      end if
      print*,' '

c     compute B_toroidal on axis

      btaxis=fpol(1)/rmaxis
      btrcen=abs(btrcen)*factb
      print'(a,f8.4)','factb =   ',factb
      print'(a,f8.4)','|BT_centr|= ',btrcen
      print'(a,f8.4)','|BT_axis| = ',abs(btaxis)

c     compute normalized rho_tor from eqdsk q profile
      call rhotor(nr)
c      phitedge=deltapsi*rhotsx*2*pi
c      rrtor=sqrt(phitedge/abs(btrcen)/pi)

c     compute flux surface averaged quantities

      rup=rmaxis
      rlw=rmaxis
      zup=zmaxis+(zbmax-zmaxis)/10.0d0
      zlw=zmaxis-(zmaxis-zbmin)/10.0d0
      call flux_average

c     locate  psi surface for q=1.5 and q=2

      rup=rmaxis
      rlw=rmaxis
      zup=(zbmax+zmaxis)/2.0d0
      zlw=(zmaxis+zbmin)/2.0d0
      q2=2.0d0
      q15=1.5d0
      call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx)
      if (q15.gt.qmin.and.q15.lt.qmax) then
        call surfq(q15,psi15)
        rhot15=frhotor(psi15)
        print'(3(a,f8.5))','psi_1.5 = ',psi15,'  rhop_1.5 = '
     .  ,sqrt(psi15),'  rhot_1.5 = ',rhot15
      end if
      if (q2.gt.qmin.and.q2.lt.qmax) then
        call surfq(q2,psi2)
        rhot2=frhotor(psi2)
        print'(3(a,f8.5))','psi_2   = ',psi2,'  rhop_2   = '
     .   ,sqrt(psi2),'  rhot_2   = ',rhot2
      end if
c
c     locate  btot=bres
c
      call bfield_res
c
      return
      end
c
c
c
      subroutine points_ox(rz,zz,rf,zf,psinvf,info)
      implicit real*8 (a-h,o-z)
      parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
      dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
      external fcnox
      common/psival/psinv
      xvec(1)=rz
      xvec(2)=zz
      tol = sqrt(dpmpar(1))
      call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
      if(info.gt.1)  then
        print'(a,i2,a,2f8.4)','    info subr points_ox =',info,
     .   '  O/X coord.',xvec
      end if
      rf=xvec(1)
      zf=xvec(2)
      psinvf=psinv
      return
      end
c
c
c
      subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
      implicit real*8 (a-h,o-z)
      dimension x(n),fvec(n),fjac(ldfjac,n)
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      common/pareq1/psia
      call equinum(x(1),x(2))
      if (iflag.ne.2) then
          fvec(1) = dpsidr/psia
          fvec(2) = dpsidz/psia
      else
          fjac(1,1) = ddpsidrr/psia
          fjac(1,2) = ddpsidrz/psia
          fjac(2,1) = ddpsidrz/psia
          fjac(2,2) = ddpsidzz/psia
      end if
      return
      end
c
c
c
      subroutine points_tgo(rz,zz,rf,zf,psin,info)
      implicit real*8 (a-h,o-z)
      parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
      dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
      external fcntgo
      common/cnpsi/h
      h=psin
      xvec(1)=rz
      xvec(2)=zz
      tol = sqrt(dpmpar(1))
      call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
      if(info.gt.1)  then
      end if
      rf=xvec(1)
      zf=xvec(2)
      return
      end
c
c
c
      subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag)
      implicit real*8 (a-h,o-z)
      dimension x(n),fvec(n),fjac(ldfjac,n)
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      common/psival/psinv
      common/cnpsi/h
      common/pareq1/psia
      call equinum(x(1),x(2))
      if (iflag.ne.2) then
          fvec(1) = psinv-h
          fvec(2) = dpsidr/psia
      else
          fjac(1,1) = dpsidr/psia
          fjac(1,2) = dpsidz/psia
          fjac(2,1) = ddpsidrr/psia
          fjac(2,2) = ddpsidrz/psia
      end if
      return
      end
c
c
c
      subroutine print_prof
      implicit real*8 (a-h,o-z)
      parameter(nnw=501,eps=1.d-4)
      dimension psinr(nnw),qpsi(nnw)
c
      common/psinr/psinr
      common/qpsi/qpsi
      common/eqnn/nr,nz,npp,nintp
      common/dens/dens,ddens
c
      write(55,*) ' #psi rhot ne Te q Jphi'
      psin=0.0d0
      rhop=0.0d0
      rhot=0.0d0
      call density(psin)
      call tor_curr_psi(eps,ajphi)
      te=temp(psin)
      qq=qpsi(1)
c
      write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
c
      nst=nr
      do i=2,nst
        psin=dble(i-1)/dble(nst-1)
        rhop=sqrt(psin)
c
        call density(psin)
        te=temp(psin)
c
        ips=int((nr-1)*psin+1)
        if(i.lt.nst) then
          call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1),
     .            psin,qq)
        else
          qq=qpsi(nr)
        end if
        rhot=frhotor(psin)
        call tor_curr_psi(psin,ajphi)
        write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
      end do
c
      return
 111  format(12(1x,e12.5))
      end
c
c
c
      subroutine surfq(qval,psival)
      implicit real*8 (a-h,o-z)
      parameter(nnw=501)
      parameter(ncnt=100,ncntt=2*ncnt+1)
      dimension psinr(nnw),qpsi(nnw)
      dimension rcon(ncntt),zcon(ncntt)
c
      common/psinr/psinr
      common/qpsi/qpsi
      common/eqnn/nr,nz,npp,nintp
c
c     locate  psi surface for q=qval
c
      call locate(qpsi,nr,qval,i1)
      call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1),
     .  qval,psival)
      ipr=1
      call contours_psi(psival,ncnt,rcon,zcon,ipr)
      return
      end
c
c
c
      subroutine bfield_res
      implicit real*8 (a-h,o-z)
      parameter(nnw=501,nnh=501)
      dimension rv(nnw),zv(nnh),psin(nnw,nnh)
      dimension btotal(nnw,nnh)
      parameter(icmx=2002)
      dimension rrcb(icmx),zzcb(icmx),ncpts(10)
c
      common/cpsin/rv,zv,psin
      common/eqnn/nr,nz,npp,nintp
      common/parbres/bres
      common/btt/btotal
c
c     Btotal on psi grid
c
      btmx=-1.0d30
      btmn=1.0d30
      do j=1,nr
        rrj=rv(j)
        do k=1,nz
          zzk=zv(k)
          call bfield(rrj,zzk,bbphi,bbr,bbz)
          btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
          if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
          if(btotal(j,k).le.btmn) btmn=btotal(j,k)
c          write(90,113) j,rrj,zzk,btotal(j,k)
        enddo
c          write(90,*) ' '
       enddo
c
c     compute Btot=Bres and Btot=Bres/2
c
      write(70,*)'#i Btot  R  z'
      do n=1,5
        bbb=bres/dble(n)
        if (bbb.ge.btmn.and.bbb.le.btmx) then
          call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1)
          do inc=1,nctot
            write(70,113) inc,bbb,rrcb(inc),zzcb(inc)
          end do
        end if
        write(70,*) ' '
      end do
c
      return
 113  format(i6,12(1x,e12.5))
      end
c
c
      subroutine profdata
      implicit real*8 (a-h,o-z)
      parameter(npmx=250,npest=npmx+4)
      dimension psrad(npmx),terad(npmx),derad(npmx),zfc(npmx)
      dimension ct(npmx,4),cz(npmx,4)
      parameter(lwrkf=npmx*4+npest*16)
      dimension tfn(npest),cfn(npest),wrkf(lwrkf),iwrkf(npest),wf(npmx)
      dimension densi(npest),ddensi(npest),d2densi(npest),wrkfd(npest)
c
      common/densbnd/psdbnd
      common/denspp/psnpp,aa,bb,cc
      common/eqnn/nr,nz,npp,nintp
      common/iipr/iprof
      common/nfile/neqdsk,nprof
      common/crad/psrad,derad,terad,zfc
      common/coffte/ct
      common/coffz/cz
      common/factb/factb
      common/coffdt/tfn
      common/coffdnst/nsfd
      common/cofffn/cfn
      common/scal/iscal
      common/facttn/factt,factn
c
c     read plasma profiles from file if iprof>0
c
      if (iscal.eq.0) then
        aat=2.0d0/3.0d0
        aan=4.0d0/3.0d0
      else
        aat=1.0d0
        aan=1.0d0
      end if
      ffact=factb
      if(iscal.eq.2) ffact=1.0d0
c
      if (iprof.gt.0) then
        read(nprof,*) npp
        read(nprof,*) psrad0,terad0,derad0,zfc0
        if(psrad0.ne.0.0d0) psrad0=0.0d0
        psrad(1)=psrad0
        terad(1)=terad0*ffact**aat*factt
        derad(1)=derad0*ffact**aan*factn
        zfc(1)=zfc0
        wf(1)=1.0d0
        do i=2,npp
          read(nprof,*) psradi,teradi,deradi,zfci
          psrad(i)=psradi
          terad(i)=teradi*ffact**aat*factt
          derad(i)=deradi*ffact**aan*factn
          zfc(i)=zfci
          wf(i)=1.0d0
        end do
c
c     spline approximation of  temperature and Zeff
c
        iopt=0
        call difcsn(psrad,terad,npmx,npp,iopt,ct,ier)
c
        iopt=0
        call difcsn(psrad,zfc,npmx,npp,iopt,cz,ier)
c
c     spline approximation of density
c
      iopt=0
      xb=0.0d0
      xe=psrad(npp)
      kspl=3
      sspl=.001d0
c
      call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
     . tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
c
      call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
      nu=1
      call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
      dnpp=densi(npp)
      ddnpp=ddensi(npp)
c
      nu=2
      call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
      d2dnpp=d2densi(npp)

        if(derad(npp).eq.0.0d0) then
          psdbnd=psrad(npp)
        else
          psnpp=psrad(npp)
          dpsb=-psnpp+psdbnd
          nn=3
          nn1=nn+1
          nn2=nn+2
          aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2)
          aa=aa/(-dpsb)**nn/2.0d0
          bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2)
          bb=bb/(-dpsb)**nn1
          cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2)
          cc=cc/(-dpsb)**nn2/2.0d0
        end if
c
      end if
c
      return
      end
c
c
c
      subroutine rhotor(nr)
      implicit real*8(a-h,o-z)
      parameter(nnw=501)
      dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4)
      common/psinr/psinr
      common/qpsi/qpsi
      common/rhotsx/rhotsx
      common/crhot/crhot
c
c     normalized toroidal rho : ~ Integral q(psi) dpsi
c
      iopt=0
      call difcsn(psinr,qpsi,nnw,nr,iopt,cq,ier)
c
      rhotnr(1)=0.0d0
      do k=1,nr-1
        dx=psinr(k+1)-psinr(k)
        drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0
     .        +dx*cq(k,4)/4.0d0)))
        rhotnr(k+1)=rhotnr(k)+drhot
      end do
      rhotsx=rhotnr(nr)
      do k=1,nr
        rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr))
      end do
c
c     spline interpolation of rhotor
c
      iopt=0
      call difcs(psinr,rhotnr,nr,iopt,crhot,ier)
      return
      end

      function frhotor_eq(psi)
      implicit real*8(a-h,o-z)
      parameter(nnw=501)
      dimension psinr(nnw),crhot(nnw,4)
      common/psinr/psinr
      common/eqnn/nr,nz,npp,nintp
      common/crhot/crhot
c
      irt=int((nr-1)*psi+1)
      if(irt.eq.0) irt=1
      if(irt.eq.nr) irt=nr-1
      dps=psi-psinr(irt)
      frhotor_eq=spli(crhot,nr,irt,dps)
      return
      end

      function frhotor(psi)
      implicit real*8(a-h,o-z)
      frhotor=frhotor_eq(psi)
      return
      end
      
      function frhotor_av(psi)
      implicit real*8(a-h,o-z)
      parameter(nnintp=101)
      dimension rpstab(nnintp),crhotq(nnintp,4)
      common/pstab/rpstab
      common/eqnn/nr,nz,npp,nintp
      common/crhotq/crhotq

      rpsi=sqrt(psi)
      ip=int((nintp-1)*rpsi+1)
      if(ip.eq.0) ip=1
      if(ip.eq.nintp) ip=nintp-1
      dps=rpsi-rpstab(ip)
      frhotor_av=spli(crhotq,nintp,ip,dps)

      return
      end

      subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi)
      implicit real*8 (a-h,o-z)
c
c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
c                   (based on an older code)
c
      parameter(nnw=501,nnh=501,icmx=2002,nna=nnw*nnh)
      dimension a(nna),ja(3,2),lx(1000),npts(10)
      dimension rcon(icmx),zcon(icmx)
      dimension rqgrid(nnw),zqgrid(nnh),psin(nnw,nnh),btotal(nnw,nnh)
c
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      common/cpsin/rqgrid,zqgrid,psin
      common/btt/btotal
      common/eqnn/nr,nz,npp,nintp
c
      data px/0.5d0/
c
      if(ichoi.eq.0) then
         do j=1,nz
          do i=1,nr
            a(nr*(j-1)+i)=psin(i,j)
          enddo
        enddo
      endif
c
      if(ichoi.eq.1) then
         do j=1,nz
          do i=1,nr
            a(nr*(j-1)+i)=btotal(i,j)
          enddo
        enddo
      endif
c
      do ico=1,icmx
        rcon(ico)=0.0d0
        zcon(ico)=0.0d0
      enddo
c
      nrqmax=nr
      nreq=nr
      nzeq=nz
      drgrd=dr
      dzgrd=dz
c
      ncon = 0
      do  i=1,10
        npts(i) = 0
      end do
      iclast = 0
      icount = 0
      mpl=0
      ix=0
      mxr = nrqmax * (nzeq - 1)
      n1 = nreq - 1
c
      do 3 jx=2,n1
      do 3 jm=jx,mxr,nrqmax
          j = jm + nrqmax
          ah=a(j)-h
          if(ah) 60,60,61
   60     if(a(jm)-h) 62,62,1
   61     if(a(jm)-h) 1,1,63
    1     ix=ix+1
          lx(ix)=-j
          if(ah) 62,62,63
   62     if(a(j-1)-h) 3,3,4
   63     if(a(j-1)-h) 4,4,3
    4     ix=ix+1
          lx(ix)=j
    3 continue
c
      do 79 jm=nreq,mxr,nrqmax
         j = jm + nrqmax
         ah=a(j)-h
         if(ah) 64,64,65
   64    if(a(j-1)-h) 66,66,76
   65    if(a(j-1)-h) 76,76,67
   76    ix=ix+1
         lx(ix)=j
         if(ah) 66,66,67
   66    if(a(jm)-h) 79,79,78
   67    if(a(jm)-h) 78,78,79
   78    ix=ix+1
         lx(ix)=-j
   79 continue
c
      do 75 jm=1,mxr,nrqmax
         j = jm + nrqmax
         if(a(j)-h) 68,68,69
   68    if(a(jm)-h)75,75,71
   69    if(a(jm)-h)71,71,75
   71    ix=ix+1
         lx(ix) =-j
   75 continue
c
      do 73 j=2,nreq
         if(a(j)-h) 58,58,59
   58    if(a(j-1)-h) 73,73,72
   59    if(a(j-1)-h) 72,72,73
   72    ix=ix+1
         lx(ix)=j
   73 continue
c
      if(ix) 50,50,8
  108 if(mpl.lt.4) go to 8
      ncon = ncon + 1
      npts(ncon) = icount - iclast
      iclast = icount
    8 in=ix
      jx=lx(in)
      jfor=0
      lda=1
      ldb=2
   30 if(jx) 21,22,22
   21 jabs=-jx
      jnb = jabs - nrqmax
      go to 23
   22 jabs=jx
      jnb=jabs-1
   23 adn=a(jabs)-a(jnb)
      if(adn) 24,9,24
   24 px=(a(jabs)-h)/adn
    9 kx = (jabs - 1) / nrqmax
      ikx = jabs - nrqmax * kx - 1
      if(jx) 25,26,26
   25 x = drgrd * ikx
      y = dzgrd * (kx - px)
      go to 27
   26 x = drgrd * (ikx - px)
      y = dzgrd * kx
   27 continue
      icount = icount + 1
      rcon(icount) = x + rqgrid(1)
      zcon(icount) = y + zqgrid(1)
      mpl= icount
      itm = 1
      ja(1,1) = jabs + nrqmax
      j=1
      if(jx) 10,10,11
   10 ja(1,1) = -jabs-1
      j=2
   11 ja(2,1)=-ja(1,1)
      ja(3,1) = -jx + 1 - nrqmax
      ja(3,2) = -jx
      ja(j,2) = jabs - nrqmax
      k= 3-j
      ja(k,2) = 1-jabs
      if (kx) 14,14,39
   39 if(ikx) 14,14,36
   36 if(ikx + 1 - nreq) 35, 37, 37
   37 if(jx) 38,38,35
   35 if(jfor) 28,29,28
   28 do 13 i=1,3
         if(jfor-ja(i,2)) 13,14,13
   13 continue
   38 lda=2
      go to 15
   14 lda=1
   15 ldb=lda
   29 do 18 k=1,3
         do 18 l=lda,ldb
            do 16 i=1,ix
            if(lx(i)-ja(k,l)) 16,17,16
   16    continue
         go to 18
   17    itm=itm+1
         inext= i
         if(jfor) 19,33,19
   33    if(itm .gt. 3) goto 20
   18 continue
   19 lx(in)=0
      if(itm .eq. 1) goto 6
   20 jfor=jx
      jx=lx(inext)
      in = inext
      go to 30
    6 if(lx(ix)) 108,7,108
    7 ix= ix-1
      if(ix) 51,51,6
   51 if(mpl.lt.4) go to  50
      ncon = ncon + 1
      npts(ncon) = icount - iclast
      iclast = icount
   50 continue
c
      return
      end
c
c
c
      subroutine contours_psi(h,np,rcn,zcn,ipr)
      implicit real*8 (a-h,o-z)
      parameter(pi=3.14159265358979d0)
      parameter(mest=4,kspl=3)
      parameter(nnw=501,nnh=501)
      parameter(nrest=nnw+4,nzest=nnh+4)
      dimension rcn(2*np+1),zcn(2*np+1)
      dimension cc(nnw*nnh),tr(nrest),tz(nzest)
      dimension czc(nrest),zeroc(mest)
c
      common/pareq1/psia
      common/pareq1a/psiaxis0
      common/coffeqn/nsr,nsz,nsft
      common/coffeq/cc
      common/coffeqt/tr,tz
      common/cnt/rup,zup,rlw,zlw
      common/rwallm/rwallm
c
      ra=rup
      rb=rlw
      za=zup
      zb=zlw
      call points_tgo(ra,za,rup,zup,h,info)
      call points_tgo(rb,zb,rlw,zlw,h,info)

      th=pi/dble(np)
      rcn(1)=rlw
      zcn(1)=zlw
      rcn(2*np+1)=rlw
      zcn(2*np+1)=zlw
      rcn(np+1)=rup
      zcn(np+1)=zup
      do ic=2,np
        zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0
        iopt=1
        call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
        if(ier.gt.0) print*,'    profil =',ier
        val=h+psiaxis0/psia
        call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)      
        if (zeroc(1).gt.rwallm) then
        rcn(ic)=zeroc(1)
        zcn(ic)=zc
        rcn(2*np+2-ic)=zeroc(2)
        zcn(2*np+2-ic)=zc
        else
        rcn(ic)=zeroc(2)
        zcn(ic)=zc
        rcn(2*np+2-ic)=zeroc(3)
        zcn(2*np+2-ic)=zc        
        end if
      end do
      if (ipr.gt.0) then
      do ii=1,2*np+1
        write(71,111) ii,h,rcn(ii),zcn(ii)
      end do
        write(71,*)
        write(71,*)
      end if
      return
111   format(i6,12(1x,e12.5))
99    format(12(1x,e12.5))
      end
c
c
c
      subroutine flux_average
      implicit real*8 (a-h,o-z)
      real*8 lam

      parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41)
      parameter(zero=0.0d0,one=1.0d0)
      parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
      parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
      parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
     .          njest*nnintp+nlest+54)
      parameter(kwrk=nnintp+nlam+njest+nlest+3)
      parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)

      dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp)
      dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp)
      dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp)
      dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp)
      dimension rcon(2*ncnt+1),zcon(2*ncnt+1)
      dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1)
      dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4)
      dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4)
      dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4)
      dimension cfc(nnintp,4),crhotq(nnintp,4)
      dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp)
      dimension alam(nlam),fhlam(nnintp,nlam)
      dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam)
      dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4))
      dimension iwrk(kwrk),wrk(lwrk)
      dimension ch01(lw01),weights(nlam)
c
      common/cent/btrcen,rcen
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      common/eqnn/nr,nz,npp,nintp
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
c
      common/pstab/rpstab
      common/flav/vvol,rri,rbav,bmxpsi,bmnpsi
      common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc
      common/cratj/cratja,cratjb,cratjpl
      common/crhotq/crhotq
      common/cnt/rup,zup,rlw,zlw
      common/bound/zbmin,zbmax
      common/rarea/rarea
c
      common/coffvl/ch
      common/coffdvl/ch01
      common/coffvlt/tjp,tlm
      common/coffvln/njpt,nlmt

c     computation of flux surface averaged quantities

      write(71,*)' #i psin R  z'

      nintp=nnintp
      ninpr=(nintp-1)/10
      
      dlam=1.0d0/dble(nlam-1)
      do l=1,nlam-1
        alam(l)=dble(l-1)*dlam
        fhlam(1,l)=sqrt(1.0d0-alam(l))
        ffhlam(l)=fhlam(1,l)
        dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l))
        weights(l)=1.0d0
      end do 
      weights(1)=0.5d0
      weights(nlam)=0.5d0
      alam(nlam)=1.0d0
      fhlam(1,nlam)=0.0d0
      ffhlam(nlam)=0.0d0
      dffhlam(nlam)=-99999.0d0

      jp=1
      anorm=2.0d0*pi*rmaxis/abs(btaxis)
      b2av=btaxis**2
      dvdpsi=2.0d0*pi*anorm
      dadpsi=2.0d0*pi/abs(btaxis)
      ratio_cdator=abs(btaxis/btrcen)
      ratio_cdbtor=1.0d0
      ratio_pltor=1.0d0
      qq=1.0d0
      fc=1.0d0
      
      pstab(1)=0.0d0
      rpstab(1)=0.0d0
      rhotqv(1)=0.0d0
      vcurrp(1)=0.0d0
      vajphiav(1)=0.0d0
      bmxpsi(1)=abs(btaxis)
      bmnpsi(1)=abs(btaxis)
      bav(1)=abs(btaxis)
      rbav(1)=1.0d0
      rri(1)=rmaxis
      varea(1)=0.0d0
      vvol(1)=0.0d0
      vratjpl(1)=ratio_pltor
      vratja(1)=ratio_cdator
      vratjb(1)=ratio_cdbtor
      ffc(1)=fc
      rhot2q=0.0d0
      
C      rup=rmaxis
C      rlw=rmaxis
C      zup=(zbmax+zmaxis)/2.0d0
C      zlw=(zmaxis+zbmin)/2.0d0

      do jp=2,nintp
        height=dble(jp-1)/dble(nintp-1)
        if(jp.eq.nintp) height=0.9999d0
        ipr=0
        jpr=mod(jp,ninpr)
        if(jpr.eq.1) ipr=1
        height2=height*height
        call contours_psi(height2,ncnt,rcon,zcon,ipr)
c
        r2iav=0.0d0
        anorm=0.0d0
        dadpsi=0.0d0
        currp=0.0d0
        b2av=0.0d0
        area=0.0d0
        volume=0.0d0
        ajphiav=0.0d0
        bbav=0.0d0
        bmmx=-1.0d+30
        bmmn=1.0d+30

        call bfield(rcon(1),zcon(1),bphi,brr,bzz)
        call tor_curr(rcon(1),zcon(1),ajphi0)
        btot0=sqrt(bphi**2+brr**2+bzz**2)
        bpoloid0=sqrt(brr**2+bzz**2)
        bv(1)=btot0
        bpv(1)=bpoloid0
        rpsim0=rcon(1)

        do inc=1,ncntt-1
          inc1=inc+1
          dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2)
          dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2)
          dlp=sqrt((rcon(inc1)-rcon(inc))**2+
     .             (zcon(inc1)-zcon(inc))**2)
          drc=(rcon(inc1)-rcon(inc))

c         compute length, area and volume defined by psi=height^2

          ph=0.5d0*(dla+dlb+dlp)
          area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp)
          area=area+sqrt(area2)
          rzp=rcon(inc1)*zcon(inc1)
          rz=rcon(inc)*zcon(inc)
          volume=pi*(rzp+rz)*drc+volume

c         compute line integral on the contour psi=height^2

          rpsim=rcon(inc1)
          zpsim=zcon(inc1)
          call bfield(rpsim,zpsim,bphi,brr,bzz)
          call tor_curr(rpsim,zpsim,ajphi)

          btot=sqrt(bphi**2+brr**2+bzz**2)
          bpoloid=sqrt(brr**2+bzz**2)
          dlpv(inc)=dlp
          bv(inc1)=btot
          bpv(inc1)=bpoloid
          
          dlph=0.5d0*dlp
          anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0)
          dadpsi=dadpsi+dlph*
     .          (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0))
          currp=currp+dlph*(bpoloid+bpoloid0)
          b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid)
          bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0)
          r2iav=r2iav+dlph*
     .          (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2))
          ajphiav=ajphiav+dlph*
     .            (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim))
          
          ajphi0=ajphi
          rpsim0=rpsim
          bpoloid0=bpoloid
          btot0=btot        

c       computation maximum/minimum B values on given flux surface

          if(btot.le.bmmn) bmmn=btot
          if(btot.ge.bmmx) bmmx=btot
        end do       

c       bav=<B> [T] ,  b2av=<B^2> [T^2] , rbav=<B>/b_min
c       anorm = int d l_p/B_p = dV/dpsi/(2pi)
c       r2iav=<1/R^2> [m^-2] ,
c       riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)),
c       rri = <B>/(|R B_tor|<1/R^2>) , used to compute I_tor  [m^-1]
c       currp = plasma current within psi=const

        bbav=bbav/anorm
        r2iav=r2iav/anorm
        dvdpsi=2.0d0*pi*anorm
        riav=dadpsi/anorm
        b2av=b2av/anorm
        vcurrp(jp)=ccj*currp
        vajphiav(jp)=ajphiav/dadpsi

c     area == varea, volume == vvol
c     flux surface minor radius == (area/pi)^1/2
c     ratio_cdator = Jcd_astra/J_phi      Jcd_astra = <Jcd.B>/B0
c     ratio_cdbtor = Jcd_jintrac/J_phi    Jcd_jintrac = <Jcd.B>/<B>
c     ratio_pltor  = Jcd_||/J_phi         Jcd_|| = <Jcd.b>

        pstab(jp)=height2
        rpstab(jp)=height
        vvol(jp)=abs(volume)
        varea(jp)=area
        bav(jp)=bbav
        rbav(jp)=bbav/bmmn
        bmxpsi(jp)=bmmx
        bmnpsi(jp)=bmmn
        rri(jp)=bav(jp)/abs(fpolv*r2iav)
        ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen))
        ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
        ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
        vratjpl(jp)=ratio_pltor
        vratja(jp)=ratio_cdator
        vratjb(jp)=ratio_cdbtor
        qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
        qqv(jp)=qq
c        
c      write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)

c       computation of rhot from calculated q profile 

        rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1))
     .         /dble(nintp-1)
        rhotqv(jp)=sqrt(rhot2q)
        
c  computation of fraction of circulating/trapped fraction fc, ft 
c  and of function H(lambda,rhop)
c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>

        fc=0.0d0
        shlam=0.0d0
        do l=nlam,1,-1
          lam=alam(l)
          srl=0.0d0
          rl2=1.0d0-lam*bv(1)/bmmx
          rl0=0.d0
          if(rl2.gt.0) rl0=sqrt(rl2)
          do inc=1,ncntt-1
            rl2=1.0d0-lam*bv(inc+1)/bmmx
            rl=0.0d0
            if(rl2.gt.0) rl=sqrt(rl2)
            srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc))
            rl0=rl
          end do
          srl=srl/anorm
          dhlam=0.5d0/srl
          fc=fc+lam/srl*weights(l)
          if(l.eq.nlam) then
            fhlam(jp,l)=0.0d0
            ffhlam(nlam*(jp-1)+l)=0.0d0
            dffhlam(nlam*(jp-1)+l)=-dhlam
            dhlam0=dhlam
          else
            shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam
            fhlam(jp,l)=shlam
            dffhlam(nlam*(jp-1)+l)=-dhlam
            dhlam0=dhlam
          end if
        end do
        fc=0.75d0*b2av/bmmx**2*fc*dlam
        ffc(jp)=fc
        
        ccfh=bmmn/bmmx/fc
        do l=1,nlam
          ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
          dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
        end do
      end do

      write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
     .'  Area Vol |I_pl| <J_phi> qq fc ratio_cdbtor'

      qqv(1)=qqv(2)
      vajphiav(1)=vajphiav(2)
      do jp=1,nintp
        rhotqv(jp)=rhotqv(jp)/rhotqv(nintp)
        if(jp.eq.nintp) then
              rhotqv(jp)=1.0d0
              rpstab(jp)=1.0d0
              pstab(jp)=1.0d0
        end if
        rhot_eq=frhotor_eq(pstab(jp))
        write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
     .    bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), 
     .    vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
      end do 

      rarea=sqrt(varea(nintp)/pi)

c     spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi
c     used for computations of dP/dV and J_cd
c     spline coefficients of rhot

        iopt=0
        call difcs(rpstab,vvol,nintp,iopt,cvol,ier)
        iopt=0
        call difcs(rpstab,rbav,nintp,iopt,crbav,ier)
        iopt=0
        call difcs(rpstab,rri,nintp,iopt,crri,ier)
        iopt=0
        call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier)
        iopt=0
        call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier)
        iopt=0
        call difcs(rpstab,vratja,nintp,iopt,cratja,ier)
        iopt=0
        call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier)
        iopt=0
        call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier)
        iopt=0
        call difcs(rpstab,varea,nintp,iopt,carea,ier)
        iopt=0
        call difcs(rpstab,ffc,nintp,iopt,cfc,ier)
        iopt=0
        call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier)

c   spline interpolation of H(lambda,rhop) and dH/dlambda

        iopt=0
        s=0.0d0
        call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam,
     .   zero,one,zero,one,ksp,ksp,s,
     .   njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier)       
        njpt=njp
        nlmt=nlm
        
        call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1,
     .                    ch01,lw01,ier)
 

      return
 99   format(20(1x,e12.5))
      end
c
c
c
      subroutine vectinit
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36,nmx=8000)
      dimension psjki(jmx,kmx,nmx)
      dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
      dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx)
      parameter(tmax=5,npts=500)
      real*8 ttv(npts+1),extv(npts+1)
      common/ttex/ttv,extv
c
      common/warm/iwarm,ilarm
      common/iiv/iiv
      common/iov/iop,iow
      common/psjki/psjki
      common/atjki/tauv,alphav
      common/dpjjki/pdjki,currj
      common/pcjki/ppabs,ccci
      common/dcjki/didst
      common/nray/nrayr,nrayth
      common/nstep/nstep
      common/tau1v/tau1v
c
      if(nstep.gt.nmx) nstep=nmx
      if(nrayr.gt.jmx) nrayr=jmx
      if(nrayth.gt.kmx) nrayth=kmx
      
c
      do i=1,nstep
        do k=1,nrayth
          do j=1,nrayr
            psjki(j,k,i)=0.0d0
            tauv(j,k,i)=0.0d0
            alphav(j,k,i)=0.0d0
            pdjki(j,k,i)=0.0d0
            ppabs(j,k,i)=0.0d0
            didst(j,k,i)=0.0d0
            ccci(j,k,i)=0.0d0
            currj(j,k,i)=0.0d0
            iiv(j,k)=1
            iop(j,k)=0
            iow(j,k)=0
            tau1v(j,k)=0.0d0
          end do
        end do
      end do
c
      if(iwarm.gt.1) then
            dt=2.0d0*tmax/dble(npts)
            do i = 1, npts+1
              ttv(i) = -tmax+dble(i-1)*dt
              extv(i)=exp(-ttv(i)*ttv(i))
            end do
      end if
c
      return
      end



      subroutine vectinit2
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36,nmx=8000)
      dimension psjki(jmx,kmx,nmx)
      dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
      dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx)

      common/iiv/iiv
      common/iov/iop,iow
      common/psjki/psjki
      common/atjki/tauv,alphav
      common/dpjjki/pdjki,currj
      common/pcjki/ppabs,ccci
      common/dcjki/didst
      common/nray/nrayr,nrayth
      common/nstep/nstep
      
      do i=1,nstep
        do k=1,nrayth
          do j=1,nrayr
            psjki(j,k,i)=0.0d0
            tauv(j,k,i)=0.0d0
            alphav(j,k,i)=0.0d0
            pdjki(j,k,i)=0.0d0
            ppabs(j,k,i)=0.0d0
            didst(j,k,i)=0.0d0
            ccci(j,k,i)=0.0d0
            currj(j,k,i)=0.0d0
            iiv(j,k)=1
            iop(j,k)=0
            iow(j,k)=0
          end do
        end do
      end do

      return
      end
c
c
c
      subroutine paraminit
      implicit real*8(a-h,o-z)
c
      common/istep/istep
      common/istgr/istpr,istpl
      common/ierr/ierr
      common/istop/istop
      common/ipol/ipolc
c
      istpr=0
      istpl=1
      ierr=0
      istep=0
      istop=0
      ipolc=0
c
      return
      end
c
c
      subroutine updatepos
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36)
      dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
      dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
      dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c
      common/nray/nrayr,nrayth
      common/grco/xco,du1o
      common/grc/xc,du1
      common/wrk/ywrk,ypwrk
c
      do j=1,nrayr
        do k=1,nrayth
          if(j.eq.1.and.k.gt.1) then
            xco(1,j,k)=xco(1,j,1)
            xco(2,j,k)=xco(2,j,1)
            xco(3,j,k)=xco(3,j,1)
            xc(1,j,k)=xc(1,j,1)
            xc(2,j,k)=xc(2,j,1)
            xc(3,j,k)=xc(3,j,1)
          else
            xco(1,j,k)=xc(1,j,k)
            xco(2,j,k)=xc(2,j,k)
            xco(3,j,k)=xc(3,j,k)
            xc(1,j,k)=ywrk(1,j,k)
            xc(2,j,k)=ywrk(2,j,k)
            xc(3,j,k)=ywrk(3,j,k)
          end if
          du1o(1,j,k)=du1(1,j,k)
          du1o(2,j,k)=du1(2,j,k)
          du1o(3,j,k)=du1(3,j,k)
        end do
      end do
c
      return
      end
c
c
c
      subroutine gradi
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36)
      dimension dffiu(jmx),ddffiu(jmx)
      dimension xc(3,jmx,kmx),xco(3,jmx,kmx)
      dimension du1(3,jmx,kmx),du1o(3,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
      dimension dgrad2v(3,jmx,kmx)
      dimension grad2(jmx,kmx)
      dimension dxv1(3),dxv2(3),dxv3(3),dgu(3)
      dimension dgg1(3),dgg2(3),dgg3(3)
      dimension df1(3),df2(3),df3(3)
c
      common/nray/nrayr,nrayth
      common/fi/dffiu,ddffiu
      common/grco/xco,du1o
      common/grc/xc,du1
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
c
c     compute  grad u and grad(S_I)
c
      do k=1,nrayth
        do j=1,nrayr
          if(j.eq.1) then
            gri(1,j,k)=0.0d0
            gri(2,j,k)=0.0d0
            gri(3,j,k)=0.0d0
            jp=j+1
            km=k-1
            if(k.eq.1) km=nrayth
            kp=k+1
            if(k.eq.nrayth) kp=1
            do iv=1,3
              dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k)
              dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
              dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
            end do
            call solg0(dxv1,dxv2,dxv3,dgu)
          else
            jm=j-1
            km=k-1
            if(k.eq.1) km=nrayth
            kp=k+1
            if(k.eq.nrayth) kp=1
            do iv=1,3
              dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
              dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
              dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
            end do
            call solg0(dxv1,dxv2,dxv3,dgu)
          end if
          du1(1,j,k)=dgu(1)
          du1(2,j,k)=dgu(2)
          du1(3,j,k)=dgu(3)
          gri(1,j,k)=dgu(1)*dffiu(j)
          gri(2,j,k)=dgu(2)*dffiu(j)
          gri(3,j,k)=dgu(3)*dffiu(j)
          grad2(j,k)=gri(1,j,k)**2+gri(2,j,k)**2+gri(3,j,k)**2
        end do
      end do
c
c     compute derivatives of grad u and grad(S_I)
c
      do k=1,nrayth
        do j=1,nrayr
          if(j.eq.1) then
            jp=j+1
            km=k-1
            if(k.eq.1) km=nrayth
            kp=k+1
            if(k.eq.nrayth) kp=1
            do iv=1,3
              dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km)
              dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k)
              dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
            end do
            df1(1)=du1(1,jp,kp)-du1(1,jp,km)
            df1(2)=du1(1,jp,k)-du1(1,j,k)
            df1(3)=du1(1,j,k)-du1o(1,j,k)
            df2(1)=du1(2,jp,kp)-du1(2,jp,km)
            df2(2)=du1(2,jp,k)-du1(2,j,k)
            df2(3)=du1(2,j,k)-du1o(2,j,k)
            df3(1)=du1(3,jp,kp)-du1(3,jp,km)
            df3(2)=du1(3,jp,k)-du1(3,j,k)
            df3(3)=du1(3,j,k)-du1o(3,j,k)
            call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3)
          else
            jm=j-1
            km=k-1
            if(k.eq.1) km=nrayth
            kp=k+1
            if(k.eq.nrayth) kp=1
            do iv=1,3
              dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k)
              dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km)
              dxv3(iv)=xc(iv,j,k)-xco(iv,j,k)
            end do
            df1(1)=du1(1,j,k)-du1(1,jm,k)
            df1(2)=du1(1,j,kp)-du1(1,j,km)
            df1(3)=du1(1,j,k)-du1o(1,j,k)
            df2(1)=du1(2,j,k)-du1(2,jm,k)
            df2(2)=du1(2,j,kp)-du1(2,j,km)
            df2(3)=du1(2,j,k)-du1o(2,j,k)
            df3(1)=du1(3,j,k)-du1(3,jm,k)
            df3(2)=du1(3,j,kp)-du1(3,j,km)
            df3(3)=du1(3,j,k)-du1o(3,j,k)
            call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3)
          end if
c
c     derivatives of u
c
          ux=du1(1,j,k)
          uy=du1(2,j,k)
          uz=du1(3,j,k)
          uxx=dgg1(1)
          uyy=dgg2(2)
          uzz=dgg3(3)
          uxy=(dgg1(2)+dgg2(1))/2.0d0
          uxz=(dgg1(3)+dgg3(1))/2.0d0
          uyz=(dgg2(3)+dgg3(2))/2.0d0
c
c     derivatives of S_I and Grad(S_I)
c
          dfu=dffiu(j)
          dfuu=ddffiu(j)
          gx=ux*dfu
          gy=uy*dfu
          gz=uz*dfu
          gxx=dfuu*ux*ux+dfu*uxx
          gyy=dfuu*uy*uy+dfu*uyy
          gzz=dfuu*uz*uz+dfu*uzz
          gxy=dfuu*ux*uy+dfu*uxy
          gxz=dfuu*ux*uz+dfu*uxz
          gyz=dfuu*uy*uz+dfu*uyz
c
          ggri(1,1,j,k)=gxx
          ggri(2,2,j,k)=gyy
          ggri(3,3,j,k)=gzz
          ggri(1,2,j,k)=gxy
          ggri(2,1,j,k)=gxy
          ggri(1,3,j,k)=gxz
          ggri(3,1,j,k)=gxz
          ggri(2,3,j,k)=gyz
          ggri(3,2,j,k)=gyz
c
c     derivatives of |Grad(S_I)|^2
c
          dgrad2v(1,j,k)=2.0d0*(gx*gxx+gy*gxy+gz*gxz)
          dgrad2v(2,j,k)=2.0d0*(gx*gxy+gy*gyy+gz*gyz)
          dgrad2v(3,j,k)=2.0d0*(gx*gxz+gy*gyz+gz*gzz)
        end do
      end do
c
      return
      end
c
c     solution of the linear system of 3 eqs : dgg . dxv = dff
c     input vectors : dxv1, dxv2, dxv3, dff
c     output vector : dgg
c
c     dff=(1,0,0)
c
      subroutine solg0(dxv1,dxv2,dxv3,dgg)
      double precision denom,aa1,aa2,aa3
      double precision dxv1(3),dxv2(3),dxv3(3),dgg(3)
      aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3))
      aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2))
      aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2))
      denom = dxv1(1)*aa1-dxv2(1)*aa2+dxv3(1)*aa3
      dgg(1) = aa1/denom
      dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom
      dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom
      return
      end
c
c     three rhs vectors df1, df2, df3
c
      subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3)
      double precision denom,a11,a21,a31,a12,a22,a32,a13,a23,a33
      double precision dxv1(3),dxv2(3),dxv3(3)
      double precision df1(3),df2(3),df3(3)
      double precision dg1(3),dg2(3),dg3(3)
      a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3))
      a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2))
      a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2))
      a12=(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))
      a22=(dxv1(1)*dxv3(3)-dxv1(3)*dxv3(1))
      a32=(dxv1(1)*dxv2(3)-dxv1(3)*dxv2(1))
      a13=(dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))
      a23=(dxv1(1)*dxv3(2)-dxv1(2)*dxv3(1))
      a33=(dxv1(1)*dxv2(2)-dxv1(2)*dxv2(1))
      denom = dxv1(1)*a11-dxv2(1)*a21+dxv3(1)*a31
      dg1(1) = (df1(1)*a11-df1(2)*a21+df1(3)*a31)/denom
      dg1(2) = -(df1(1)*a12-df1(2)*a22+df1(3)*a32)/denom
      dg1(3) = (df1(1)*a13-df1(2)*a23+df1(3)*a33)/denom
      dg2(1) = (df2(1)*a11-df2(2)*a21+df2(3)*a31)/denom
      dg2(2) = -(df2(1)*a12-df2(2)*a22+df2(3)*a32)/denom
      dg2(3) = (df2(1)*a13-df2(2)*a23+df2(3)*a33)/denom
      dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom
      dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom
      dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom
      return
      end
c
c     Runge-Kutta integrator
c
      subroutine rkint4
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,jmx=31,kmx=36)
      dimension y(ndim),yy(ndim)
      dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim)
      dimension dgr2(3),dgr(3),ddgr(3,3)
      dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
      dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
      common/nray/nrayr,nrayth
      common/dsds/dst
c
      common/wrk/ywrk,ypwrk
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
      common/gr/gr2
      common/dgr/dgr2,dgr,ddgr
      common/igrad/igrad
c
      h=dst
      hh=h*0.5d0
      h6=h/6.0d0
c
      do j=1,nrayr
        kkk=nrayth
        if(j.eq.1) kkk=1
        do k=1,kkk
          do ieq=1,ndim
            y(ieq)=ywrk(ieq,j,k)
            fk1(ieq)=ypwrk(ieq,j,k)
            yy(ieq)=y(ieq)+fk1(ieq)*hh
          end do
          gr2=grad2(j,k)
          do iv=1,3
            dgr2(iv)=dgrad2v(iv,j,k)
            dgr(iv)=gri(iv,j,k)
            do jv=1,3
              ddgr(iv,jv)=ggri(iv,jv,j,k)
            end do
          end do
          call fwork(yy,fk2)
c
          do ieq=1,ndim
            yy(ieq)=y(ieq)+fk2(ieq)*hh
          end do
          call fwork(yy,fk3)
c
          do ieq=1,ndim
            yy(ieq)=y(ieq)+fk3(ieq)*h
          end do
          call fwork(yy,fk4)
c
          do ieq=1,ndim
            ywrk(ieq,j,k)=y(ieq)
     .        +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq))
          end do
        end do
      end do
c
      call updatepos
c
      if(igrad.eq.1) call gradi
c
      return
      end
c
c
c
      subroutine gwork(j,k)
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,jmx=31,kmx=36)
      dimension yy(ndim),yyp(ndim)
      dimension dgr2(3),dgr(3),ddgr(3,3)
      dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
      dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
      common/igrad/igrad
c
      common/wrk/ywrk,ypwrk
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
      common/gr/gr2
      common/dgr/dgr2,dgr,ddgr
c
c     begin --- update vector yy
c
      do ieq=1,ndim
        yy(ieq)=ywrk(ieq,j,k)
      end do
c
      if(igrad.eq.1) then
        gr2=grad2(j,k)
        do iv=1,3
          dgr2(iv)=dgrad2v(iv,j,k)
          dgr(iv)=gri(iv,j,k)
          do jv=1,3
            ddgr(iv,jv)=ggri(iv,jv,j,k)
          end do
        end do
      end if
c
      call fwork(yy,yyp)
c
      do ieq=1,ndim
         ypwrk(ieq,j,k)=yyp(ieq)
      end do
c
c     end --- update vector yy
c
      return
      end
c
c
c
      subroutine fwork(y,dery)
      implicit real*8 (a-h,o-z)
      parameter(ndim=6)
      dimension y(ndim),dery(ndim)
      dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
      dimension derdxv(3),danpldxv(3),derdnv(3)
      dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
c
      common/gr/gr2
      common/dgr/dgr2,dgr,ddgr
      common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
      common/mode/sox
      common/nplr/anpl,anpr
      common/bb/bv
      common/dbb/derbv
      common/igrad/igrad
      common/xgxg/xg
      common/ygyg/yg
      common/dxgyg/derxg,deryg
      common/vgv/vgm,derdnm
      common/dersdst/dersdst
      common/ierr/ierr
      common/anv/anv
      common/xv/xv
      common/idst/idst
c
      xx=y(1)
      yy=y(2)
      zz=y(3)
      
      xv(1)=y(1)
      xv(2)=y(2)
      xv(3)=y(3)
c
      anv(1) = y(4)
      anv(2) = y(5)
      anv(3) = y(6)
c
C      rr=sqrt(xx**2+yy**2)
C      phi=acos(xx/rr)
C      if (yy.lt.0.0d0) then
C        phi=-phi
C      end if
c
      call plas_deriv(xx,yy,zz)
c
      an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3)
      anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3)
c
      if(abs(anpl).gt.0.99d0) then
        if(abs(anpl).le.1.05d0) then
         ierr=97
        else 
         ierr=98
        end if 
      end if
c
      anpl2=anpl*anpl
      dnl=1.0d0-anpl2
      anpr2=an2-anpl2
      anpr=0.0d0
      if(anpr2.gt.0.0d0) anpr=sqrt(anpr2)
      yg2=yg**2
c
      an2s=1.0d0
      dan2sdxg=0.0d0
      dan2sdyg=0.0d0
      dan2sdnpl=0.0d0
      del=0.0d0
      fdia=0.0d0
      dfdiadnpl=0.0d0
      dfdiadxg=0.0d0
      dfdiadyg=0.0d0
c
      duh=1.0d0-xg-yg2
      if(xg.gt.0.0d0) then
        del=sqrt(dnl*dnl+4.0d0*anpl*anpl*(1.0d0-xg)/yg2)
        an2s=1.0d0-xg -xg*yg2*(1.0d0+anpl2+sox*del)/duh/2.0d0
c
        dan2sdxg=-yg2*(1.0d0-yg2)*(1.0d0+anpl2+sox*del)/duh**2/2.0d0
     .           +sox*xg*anpl2/(del*duh)-1.0d0
        dan2sdyg=-xg*yg*(1.0d0-xg)*(1.0d0+anpl2+sox*del)/duh**2
     .           +2.0d0*sox*xg*(1.0d0-xg)*anpl2/(yg*del*duh)
        dan2sdnpl=-xg*yg2*anpl/duh
     .            -sox*xg*anpl*(2.0d0*(1.0d0-xg)-yg2*dnl)/(del*duh)
c
        if(igrad.gt.0) then
        ddelnpl2=2.0d0*(2.0d0*(1.0d0-xg)*(1.0d0+3.0d0*anpl2**2)
     .           -yg2*dnl**3)/yg2/del**3
        fdia=-xg*yg2*(1.0d0+sox*ddelnpl2/2.0d0)/duh
        derdel=2.0d0*(1.0d0-xg)*anpl2*(1.0d0+3.0d0*anpl2**2)
     .         -dnl**2*(1.0d0+3.0d0*anpl2)*yg2
        derdel=4.0d0*derdel/(yg**5*del**5)
        ddelnpl2y=2.0d0*(1.0d0-xg)*derdel
        ddelnpl2x=yg*derdel
        dfdiadnpl=24.0d0*sox*xg*(1.0d0-xg)*anpl*(1.0d0-anpl**4)
     .            /(yg2*del**5)
        dfdiadxg=-yg2*(1.0d0-yg2)/duh**2-sox*yg2*((1.0d0-yg2)*ddelnpl2
     .           +xg*duh*ddelnpl2x)/(2.0d0*duh**2)
        dfdiadyg=-2.0d0*yg*xg*(1.0d0-xg)/duh**2
     .           -sox*xg*yg*(2.0d0*(1.0d0-xg)*ddelnpl2
     .           +yg*duh*ddelnpl2y)/(2.0d0*duh**2)
c
        end if
      end if
c
      bdotgr=0.0d0
      do iv=1,3
        bdotgr=bdotgr+bv(iv)*dgr(iv)
        dbgr(iv)=0.0d0
        do jv=1,3
          dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv)
        end do
      end do
c
      derdnm=0.0d0
c
      do iv=1,3
        danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv)
     .               +anv(3)*derbv(3,iv)
        derdxv(iv)=-(derxg(iv)*dan2sdxg
     .             +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl)
        derdxv(iv)=derdxv(iv)-igrad*dgr2(iv)
c
        derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5d0*bdotgr**2*
     .  (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl)
c
        derdnv(iv)=2.0d0*anv(iv)-bv(iv)*dan2sdnpl
     .             +0.5d0*bdotgr**2*bv(iv)*dfdiadnpl
c
        derdnm=derdnm+derdnv(iv)**2
      end do
c
      derdnm=sqrt(derdnm)
c
      derdom=-2.0d0*an2+2.0d0*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl
     .       +2.0d0*igrad*gr2
     .       -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0d0
     .                   +anpl*dfdiadnpl/2.0d0)
c
      if (idst.eq.0) then
c   integration variable: s
        denom=-derdnm
      else if (idst.eq.1) then
c   integration variable: c*t
        denom=derdom
      else
c   integration variable: Sr
        denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3))
      end if
c
c   coefficient for integration in s
c   ds/dst, where st is the integration variable
      dersdst=-derdnm/denom
c
      dery(1) = -derdnv(1)/denom
      dery(2) = -derdnv(2)/denom
      dery(3) = -derdnv(3)/denom
      dery(4) = derdxv(1)/denom
      dery(5) = derdxv(2)/denom
      dery(6) = derdxv(3)/denom
c
c     vgv :  ~ group velocity
c
      vgm=0
      do iv=1,3
        vgv(iv)=-derdnv(iv)/derdom
        vgm=vgm+vgv(iv)**2
      end do
      vgm=sqrt(vgm)
c
c     dd :  dispersion relation (real part)
c     ddi :  dispersion relation (imaginary part)
c
      dd=an2-an2s-igrad*(gr2-0.5d0*bdotgr**2*fdia)
      ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3)
c
      return
      end
c
c
c
      subroutine plas_deriv(xx,yy,zz)
      implicit real*8 (a-h,o-z)
      parameter(pi=3.14159265358979d0)
      dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3)
      dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3)
c
      common/parbres/bres
      common/parpl/brr,bphi,bzz,ajphi
      common/btot/btot
      common/bb/bv
      common/dbb/derbv
      common/xgxg/xg
      common/dxgdps/dxgdpsi
      common/ygyg/yg
      common/dxgyg/derxg,deryg
      common/iieq/iequil
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      common/psival/psinv
      common/sgnib/sgnbphi,sgniphi
c
      xg=0.0d0
      yg=9.9d1
c
      do iv=1,3
        derxg(iv)=0.0d0
        deryg(iv)=0.0d0
        bv(iv)=0.0d0
        dbtot(iv)=0.0d0
        do jv=1,3
          dbv(iv,jv)=0.0d0
          derbv(iv,jv)=0.0d0
          dbvcdc(iv,jv)=0.0d0
          dbvcdc(iv,jv)=0.0d0
          dbvdc(iv,jv)=0.0d0
        end do
      end do
c
      if(iequil.eq.0) return
c
c     cylindrical coordinates
c
      rr2=xx**2+yy**2
      rr=sqrt(rr2)
      csphi=xx/rr
      snphi=yy/rr
c
      bv(1)=-snphi*sgnbphi
      bv(2)=csphi*sgnbphi
c
c     convert from cm to meters
c
      zzm=1.0d-2*zz
      rrm=1.0d-2*rr
c
      if(iequil.eq.1) then
        call equian(rrm,zzm)
        yg=fpolv/(rrm*bres)
      end if
c
      if(iequil.eq.2) then
        call equinum(rrm,zzm)
        call sub_xg_derxg
        yg=fpolv/(rrm*bres)
        bphi=fpolv/rrm
        btot=abs(bphi)
        if(psinv.lt.0.0d0) return
      end if
c
c     B = f(psi)/R e_phi+  grad psi x e_phi/R
c
c     bvc(i) = B_i in cylindrical coordinates
c     bv(i)  = B_i in cartesian coordinates
c
      bphi=fpolv/rrm
      brr=-dpsidz/rrm
      bzz= dpsidr/rrm
c
      dfpolv=ffpv/fpolv
c
      bvc(1)=brr
      bvc(2)=bphi
      bvc(3)=bzz
c
      bv(1)=bvc(1)*csphi-bvc(2)*snphi
      bv(2)=bvc(1)*snphi+bvc(2)*csphi
      bv(3)=bvc(3)
c
      b2tot=bv(1)**2+bv(2)**2+bv(3)**2
      btot=sqrt(b2tot)
c
c     dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
c
      dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm
      dbvcdc(1,3)=-ddpsidzz/rrm
      dbvcdc(2,1)= dfpolv*dpsidr/rrm-bphi/rrm
      dbvcdc(2,3)= dfpolv*dpsidz/rrm
      dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm
      dbvcdc(3,3)= ddpsidrz/rrm
c
c     dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
c
      dbvdc(1,1) = dbvcdc(1,1)*csphi-dbvcdc(2,1)*snphi
      dbvdc(1,2) = -bv(2)
      dbvdc(1,3) = dbvcdc(1,3)*csphi-dbvcdc(2,3)*snphi
      dbvdc(2,1) = dbvcdc(1,1)*snphi+dbvcdc(2,1)*csphi
      dbvdc(2,2) = bv(1)
      dbvdc(2,3) = dbvcdc(1,3)*snphi+dbvcdc(2,3)*csphi
      dbvdc(3,1) = dbvcdc(3,1)
      dbvdc(3,2) = dbvcdc(3,2)
      dbvdc(3,3) = dbvcdc(3,3)
c
      drrdx=csphi
      drrdy=snphi
      dphidx=-snphi/rrm
      dphidy=csphi/rrm
c
c     dbv(iv,jv) = d Bcart(iv) / dxvcart(jv)
c
      dbv(1,1)=drrdx*dbvdc(1,1)+dphidx*dbvdc(1,2)
      dbv(1,2)=drrdy*dbvdc(1,1)+dphidy*dbvdc(1,2)
      dbv(1,3)=dbvdc(1,3)
      dbv(2,1)=drrdx*dbvdc(2,1)+dphidx*dbvdc(2,2)
      dbv(2,2)=drrdy*dbvdc(2,1)+dphidy*dbvdc(2,2)
      dbv(2,3)=dbvdc(2,3)
      dbv(3,1)=drrdx*dbvdc(3,1)+dphidx*dbvdc(3,2)
      dbv(3,2)=drrdy*dbvdc(3,1)+dphidy*dbvdc(3,2)
      dbv(3,3)=dbvdc(3,3)
c
      dbtot(1)=(bv(1)*dbv(1,1)+bv(2)*dbv(2,1)+bv(3)*dbv(3,1))/btot
      dbtot(2)=(bv(1)*dbv(1,2)+bv(2)*dbv(2,2)+bv(3)*dbv(3,2))/btot
      dbtot(3)=(bv(1)*dbv(1,3)+bv(2)*dbv(2,3)+bv(3)*dbv(3,3))/btot
c
      yg=btot/Bres
c
c     convert spatial derivatives from dummy/m -> dummy/cm
c     to be used in fwork
c
c     bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z
c
      do iv=1,3
        deryg(iv)=1.0d-2*dbtot(iv)/Bres
        bv(iv)=bv(iv)/btot
        do jv=1,3
          derbv(iv,jv)=1.0d-2*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot
        end do
      end do
c
      derxg(1)=1.0d-2*drrdx*dpsidr*dxgdpsi
      derxg(2)=1.0d-2*drrdy*dpsidr*dxgdpsi
      derxg(3)=1.0d-2*dpsidz*dxgdpsi
c
c     current density computation in Ampere/m^2
c     ccj==1/mu_0
c
      ccj=1.0d+7/(4.0d0*pi)
      ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1))
c      ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3))
c      ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2))
c
      return
      end
c
c
c
      subroutine equian(rrm,zzm)
      implicit real*8 (a-h,o-z)
c
      common/parqq/q0,qa,alq
      common/parban/b0,rr0m,zr0m,rpam
      common/psival/psinv
      common/pareq1/psia
      common/densbnd/psdbnd
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      common/xgxg/xg
      common/dxgdps/dxgdpsi
      common/xgcn/xgcn
      common/dens/dens,ddens
      common/sgnib/sgnbphi,sgniphi
      common/bmxmn/bmxi,bmni
      common/fc/fci
      common/iipr/iprof
c
      if(iprof.eq.0) psdbnd=1.0d0
c
      dpsidrp=0.0d0
      d2psidrp=0.0d0
c
c     simple model for equilibrium
c
      rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2)
      rn=rpm/rpam
c
      snt=0.0d0
      cst=1.0d0
      if (rpm.gt.0.0d0) then
        snt=(zzm-zr0m)/rpm
        cst=(rrm-rr0m)/rpm
      end if
c
c     valore fittizio di psinv e di psia:
      psinv=rn**2
      psia=sgniphi
c
      sgn=sgniphi*sgnbphi
      if(rn.le.1.0d0) then
        qq=q0+(qa-q0)*rn**alq
        dpsidrp=B0*rpam*rn/qq*sgn
        dqq=alq*(qa-q0)*rn**(alq-1.0d0)
        d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn
      else
        dpsidrp=B0*rpam/qa/rn*sgn
        d2psidrp=-B0/qa/rn**2 *sgn
      end if
c
      fpolv=sgnbphi*b0*rr0m
      dfpolv=0.0d0
      ffpv=sgniphi*fpolv*dfpolv
c
      dpsidr=dpsidrp*cst
      dpsidz=dpsidrp*snt
      ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp
      ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm)
      ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp
c
      dens=0.0d0
      xg=0.0d0
      dxgdpsi=0.0d0
      if(psinv.lt.psdbnd) then
        call density(psinv)
      end if
      xg=xgcn*dens
      ddensdrp=2.0d0*rn*ddens
      if(dpsidrp.ne.0.0d0) dxgdpsi=xgcn*ddensdrp/(dpsidrp*rpam)
c
      bmax=sqrt(fpolv**2+dpsidrp**2)/(rr0m-rpam*rn)
      bmin=sqrt(fpolv**2+dpsidrp**2)/(rr0m+rpam*rn)
c
      return
      end
c
c
c
      subroutine equinum(rpsim,zpsim)
      implicit real*8 (a-h,o-z)
      parameter(nnw=501,nnh=501)
      parameter(nrest=nnw+4,nzest=nnh+4)
      parameter(lwrk=8,liwrk=2)
      dimension ccspl(nnw*nnh)
      dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk)
      dimension rrs(1),zzs(1),ffspl(1)
      parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh)
      parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh)
      parameter(lw11=nnw*3+nnh*3+nnw*nnh)
      parameter(nrs=1,nzs=1)
      dimension cc01(lw10),cc10(lw01),cc02(lw02),cc20(lw20),cc11(lw11)
      dimension tfp(nrest),cfp(nrest),wrkfd(nrest)
c
      common/eqnn/nr,nz,npp,nintp
      common/psival/psinv
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      common/pareq1/psia
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      common/pareq1a/psiaxis0
c
      common/coffeqt/tr,tz
      common/coffeqtp/tfp
      common/coffeq/ccspl
      common/coffeqd/cc01,cc10,cc20,cc02,cc11
      common/coffeqn/nsrt,nszt,nsft
      common/cofffp/cfp
      common/fpas/fpolas
c
      psinv=-1.0d0
c
      dpsidr=0.0d0
      dpsidz=0.0d0
      ddpsidrr=0.0d0
      ddpsidzz=0.0d0
      ddpsidrz=0.0d0
c
c     here lengths are measured in meters
c
      fpolv=fpolas
      ffpv=0.0d0
c
      if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return
      if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return
c
      rrs(1)=rpsim
      zzs(1)=zpsim
      nsr=nsrt
      nsz=nszt
      call fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
     .     rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
      psinv=ffspl(1)-psiaxis0/psia
c      if(psinv.lt.0.0d0)
c     .   print'(a,3e12.4)', '    psin < 0 , R, z  ',psinv,rpsim,zpsim
c
      nur=1
      nuz=0
      kkr=3-nur
      kkz=3-nuz
      iwr=1+(nr-nur-4)*(nz-nuz-4)
      iwz=iwr+4-nur
      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz,
     .     rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2))
      dpsidr= ffspl(1)*psia
c
      nur=0
      nuz=1
      kkr=3-nur
      kkz=3-nuz
      iwr=1+(nr-nur-4)*(nz-nuz-4)
      iwz=iwr+4-nur
      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz,
     .     rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2))
      dpsidz= ffspl(1)*psia
c
      nur=2
      nuz=0
      kkr=3-nur
      kkz=3-nuz
      iwr=1+(nr-nur-4)*(nz-nuz-4)
      iwz=iwr+4-nur
      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz,
     .     rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2))
      ddpsidrr= ffspl(1)*psia
c
      nur=0
      nuz=2
      kkr=3-nur
      kkz=3-nuz
      iwr=1+(nr-nur-4)*(nz-nuz-4)
      iwz=iwr+4-nur
      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz,
     .     rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2))
      ddpsidzz= ffspl(1)*psia
c
      nur=1
      nuz=1
      kkr=3-nur
      kkz=3-nuz
      iwr=1+(nr-nur-4)*(nz-nuz-4)
      iwz=iwr+4-nur
      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz,
     .     rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2))
      ddpsidrz= ffspl(1)*psia
c
      if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then
        rrs(1)=psinv
        call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
        fpolv=ffspl(1)
        nu=1
        call splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
        dfpolv=ffspl(1)
        ffpv=fpolv*dfpolv/psia
      end if
c
      return
      end
c
c
c
      subroutine bfield(rpsim,zpsim,bphi,brr,bzz)
      implicit real*8 (a-h,o-z)
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      call equinum(rpsim,zpsim)
      bphi=fpolv/rpsim
      brr=-dpsidz/rpsim
      bzz= dpsidr/rpsim
      return
      end
c
      subroutine tor_curr_psi(h,ajphi)
      implicit real*8 (a-h,o-z)
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      call psi_raxis(h,r1,r2)
      call tor_curr(r2,zmaxis,ajphi)
      return
      end

c
      subroutine tor_curr(rpsim,zpsim,ajphi)
      implicit real*8 (a-h,o-z)
      parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
      common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
      call equinum(rpsim,zpsim)
      bzz= dpsidr/rpsim
      dbvcdc13=-ddpsidzz/rpsim
      dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
      ajphi=ccj*(dbvcdc13-dbvcdc31)
      return
      end
c
      subroutine psi_raxis(h,r1,r2)
      implicit real*8 (a-h,o-z)
      parameter(mest=4,kspl=3)
      parameter(nnw=501,nnh=501)
      parameter(nrest=nnw+4,nzest=nnh+4)
      dimension cc(nnw*nnh),tr(nrest),tz(nzest)
      dimension czc(nrest),zeroc(mest)
c
      common/pareq1/psia
      common/pareq1a/psiaxis0
      common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
      common/coffeqn/nsr,nsz,nsft
      common/coffeq/cc
      common/coffeqt/tr,tz
c
      iopt=1
      zc=zmaxis
      call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
      if(ier.gt.0) print*,'    profil =',ier
      val=h+psiaxis0/psia
      call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
      r1=zeroc(1)
      r2=zeroc(2)
      return
      end
c
c
c
      subroutine sub_xg_derxg
      implicit real*8 (a-h,o-z)
      common/psival/psinv
      common/pareq1/psia
      common/xgxg/xg
      common/dxgdps/dxgdpsi
      common/xgcn/xgcn
      common/dens/dens,ddenspsin
      xg=0.0d0
      dxgdpsi=0.0d0
c      if(psinv.le.psdbnd.and.psinv.ge.0) then
        call density(psinv)
        xg=xgcn*dens
        dxgdpsi=xgcn*ddenspsin/psia
c      end if
      return
      end
c
c
c
      subroutine density(arg)
      implicit real*8 (a-h,o-z)
      parameter(npmx=250,npest=npmx+4)
      dimension xxs(1),ffs(1)
      dimension tfn(npest),cfn(npest),wrkfd(npest)
c
      common/densbnd/psdbnd
      common/denspp/psnpp,aad,bbd,ccd
      common/iipr/iprof
      common/pardens/dens0,aln1,aln2
      common/dens/dens,ddens
      common/coffdt/tfn
      common/coffdnst/nsfd
      common/cofffn/cfn
c
c     computation of density [10^19 m^-3] and derivative wrt psi
c
      dens=0.0d0
      ddens=0.0d0
      if(arg.gt.psdbnd.or.arg.lt.0.0d0) return
c
      if(iprof.eq.0) then
        if(arg.gt.1.0d0) return
        profd=(1.0d0-arg**aln1)**aln2
        dens=dens0*profd
        dprofd=-aln1*aln2*arg**(aln1-1.0d0)
     .         *(1.0d0-arg**aln1)**(aln2-1.0d0)
        ddens=dens0*dprofd
      else
        if(arg.le.psdbnd.and.arg.gt.psnpp) then
c
c     cubic interpolation for 1 < psi < psdbnd
c
          nn=3
          nn1=nn+1
          nn2=nn+2
          dpsib=arg-psdbnd
          dens=aad*dpsib**nn+bbd*dpsib**nn1+ccd*dpsib**nn2
          ddens=nn*aad*dpsib**(nn-1)+nn1*bbd*dpsib**nn
     .          +nn2*ccd*dpsib**nn1
        else
          xxs(1)=arg
          ier=0
          call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
          dens=ffs(1)
          nu=1
          ier=0
          call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
          ddens=ffs(1)
          if(ier.gt.0) print*,ier
          if(abs(dens).lt.1.0d-10) dens=0.0d0
        end if
        if(dens.lt.0.0d0) print*,'    DENSITY NEGATIVE',dens
c
      end if
      return
      end
c
c
c
      function temp(arg)
      implicit real*8 (a-h,o-z)
      parameter(npmx=250)
      dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4)
c
      common/parqte/te0,dte0,alt1,alt2
      common/iipr/iprof
      common/crad/psrad,derad,terad,zfc
      common/coffte/ct
      common/eqnn/nr,nz,npp,nintp
c
      temp=0.0d0
      if(arg.ge.1.0d0.or.arg.lt.0.0d0) return
      if(iprof.eq.0) then
        proft=(1.0d0-arg**alt1)**alt2
        temp=(te0-dte0)*proft+dte0
      else
        call locate(psrad,npp,arg,k)
        if(k.eq.0) k=1
        if(k.eq.npp) k=npp-1
        dps=arg-psrad(k)
        temp=spli(ct,npmx,k,dps)
      endif
      return
      end
c
c
c
      function fzeff(arg)
      implicit real*8 (a-h,o-z)
      parameter(npmx=250)
      dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4)
c
      common/iipr/iprof
      common/crad/psrad,derad,terad,zfc
      common/coffz/cz
      common/eqnn/nr,nz,npp,nintp
      common/zz/Zeff
c
      fzeff=1
      ps=arg
      if(ps.gt.1.0d0.and.ps.lt.0.0d0) return
      if(iprof.eq.0) then
        fzeff=zeff
      else
        call locate(psrad,npp,ps,k)
        if(k.eq.0) k=1
        if(k.eq.npp) k=npp-1
        dps=ps-psrad(k)
        fzeff=spli(cz,npmx,k,dps)
      endif
      return
      end
c
c     beam tracing initial conditions  igrad=1
c
      subroutine ic_gb
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,ndimm=3)
      parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
      parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
      dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
      dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
      dimension dffiu(jmx),ddffiu(jmx)
      dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
      complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy
      complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy
      complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
      complex*16 catand
      external catand
c      parameter(ui=(0.0d0,1.0d0))
c
      common/nray/nrayr,nrayth
      common/rwmax/rwmax
      common/parwv/ak0,akinv,fhz
      common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
      common/anic/anx0c,any0c,anz0c
      common/wrk/ywrk0,ypwrk0
      common/grc/xc0,du10
      common/fi/dffiu,ddffiu
      common/mirr/x00,y00,z00
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
      common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
c
      ui=(0.0d0,1.0d0)
      csth=anz0c
      snth=sqrt(1.0d0-csth**2)
      csps=1.0d0
      snps=0.0d0
      if(snth.gt.0.0d0) then
        csps=any0c/snth
        snps=anx0c/snth
      end if
c
      phiwrad=phiw*cvdr
      phirrad=phir*cvdr
      csphiw=cos(phiwrad)
      snphiw=sin(phiwrad)
c      csphir=cos(phirrad)
c      snphir=sin(phirrad)
c
      wwcsi=2.0d0*akinv/wcsi**2
      wweta=2.0d0*akinv/weta**2
c
      if(phir.ne.phiw) then
        sk=(rcicsi+rcieta)
        sw=(wwcsi+wweta)
        dk=(rcicsi-rcieta)
        dw=(wwcsi-wweta)
        ts=-(dk*sin(2.0d0*phirrad)-ui*dw*sin(2.0d0*phiwrad))
        tc=(dk*cos(2.0d0*phirrad)-ui*dw*cos(2.0d0*phiwrad))
        phic=0.5d0*catand(ts/tc)
        ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic))
        sss=sk-ui*sw
        qi1=0.5d0*(sss+ddd)
        qi2=0.5d0*(sss-ddd)
        rci1=dble(qi1)
        ww1=-dimag(qi1)
        rci2=dble(qi2)
        ww2=-dimag(qi2)
      else
        rci1=rcicsi
        rci2=rcieta
        ww1=wwcsi
        ww2=wweta
        phic=-phiwrad
        qi1=rci1-ui*ww1
        qi2=rci2-ui*ww2
      end if
      
c      w01=sqrt(2.0d0*akinv/ww1)
c      z01=-rci1/(rci1**2+ww1**2)
c      w02=sqrt(2.0d0*akinv/ww2)
c      z02=-rci2/(rci2**2+ww2**2)

      qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2
      qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2
      qqxy=-(qi1-qi2)*sin(2.0d0*phic)
      wwxx=-dimag(qqxx)
      wwyy=-dimag(qqyy)
      wwxy=-dimag(qqxy)/2.0d0
      rcixx=dble(qqxx)
      rciyy=dble(qqyy)
      rcixy=dble(qqxy)/2.0d0

      dqi1=-qi1**2
      dqi2=-qi2**2
      d2qi1=2*qi1**3
      d2qi2=2*qi2**3
      dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2
      dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2
      dqqxy=-(dqi1-dqi2)*sin(2.0d0*phic)
      d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2
      d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2
      d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic)
      dwwxx=-dimag(dqqxx)
      dwwyy=-dimag(dqqyy)
      dwwxy=-dimag(dqqxy)/2.0d0
      d2wwxx=-dimag(d2qqxx)
      d2wwyy=-dimag(d2qqyy)
      d2wwxy=-dimag(d2qqxy)/2.0d0
      drcixx=dble(dqqxx)
      drciyy=dble(dqqyy)
      drcixy=dble(dqqxy)/2.0d0

      dr=1.0d0
      if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
      da=2.0d0*pi/dble(nrayth)
c
      ddfu=2.0d0*dr**2*akinv
      do j=1,nrayr
        u=dble(j-1)
c        ffi=u**2*ddfu/2.0d0
        dffiu(j)=u*ddfu
        ddffiu(j)=ddfu
        do k=1,nrayth
          alfak=(k-1)*da
          dcsiw=dr*cos(alfak)*wcsi
          detaw=dr*sin(alfak)*weta
          dx0t=dcsiw*csphiw-detaw*snphiw
          dy0t=dcsiw*snphiw+detaw*csphiw
          x0t=u*dx0t
          y0t=u*dy0t
          z0t=-0.5d0*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t)
c
c          csiw=u*dcsiw
c          etaw=u*detaw
c          csir=x0t*csphir+y0t*snphir
c          etar=-x0t*snphir+y0t*csphir
c
          dx0= x0t*csps+snps*(y0t*csth+z0t*snth)
          dy0=-x0t*snps+csps*(y0t*csth+z0t*snth)
          dz0= z0t*csth-y0t*snth
          x0=x00+dx0
          y0=y00+dy0
          z0=z00+dz0

            gxt=x0t*wwxx+y0t*wwxy
            gyt=x0t*wwxy+y0t*wwyy
            gzt=0.5d0*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy
            gr2=gxt*gxt+gyt*gyt+gzt*gzt
            gxxt=wwxx
            gyyt=wwyy
            gzzt=0.5d0*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy
            gxyt=wwxy
            gxzt=x0t*dwwxx+y0t*dwwxy
            gyzt=x0t*dwwxy+y0t*dwwyy
            dgr2xt=2.0d0*(gxt*gxxt+gyt*gxyt+gzt*gxzt)
            dgr2yt=2.0d0*(gxt*gxyt+gyt*gyyt+gzt*gyzt)
            dgr2zt=2.0d0*(gxt*gxzt+gyt*gyzt+gzt*gzzt)
            dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth)
            dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth)
            dgr2z= dgr2zt*csth-dgr2yt*snth
            gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth)
            gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth)
            gri(3,j,k)=gzt*csth-gyt*snth
            ggri(1,1,j,k)=gxxt*csps**2+
     .        snps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)+
     .        2.0d0*snps*csps*(gxyt*csth+gxzt*snth)
            ggri(2,2,j,k)=gxxt*snps**2+
     .        csps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)-
     .        2.0d0*snps*csps*(gxyt*csth+gxzt*snth)
            ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2-2.0d0*csth*snth*gyzt
            ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt
     .        +2.0d0*csth*snth*gyzt)
     .        +(csps**2-snps**2)*(snth*gxzt+csth*gxyt)
            ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt)
     .        +(csth**2-snth**2)*snps*gyzt+csps*(csth*gxzt-snth*gxyt)
            ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt)
     .        +(csth**2-snth**2)*csps*gyzt+snps*(snth*gxyt-csth*gxzt)
            ggri(2,1,j,k)=ggri(1,2,j,k)
            ggri(3,1,j,k)=ggri(1,3,j,k)
            ggri(3,2,j,k)=ggri(2,3,j,k)
c
            du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu
            du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu
            du1tz=0.5d0*u*(dx0t**2*dwwxx+dy0t**2*dwwyy+
     .            2.0d0*dx0t*dy0t*dwwxy)/ddfu
c
          pppx=x0t*rcixx+y0t*rcixy
          pppy=x0t*rcixy+y0t*rciyy
          denpp=pppx*gxt+pppy*gyt
          if (denpp.ne.0.0d0) then
            ppx=-pppx*gzt/denpp
            ppy=-pppy*gzt/denpp
          else
            ppx=0.0d0
            ppy=0.0d0
          end if
c
          anzt=sqrt((1.0d0+gr2)/(1.0d0+ppx**2+ppy**2))
          anxt=ppx*anzt
          anyt=ppy*anzt
c
          anx= anxt*csps+snps*(anyt*csth+anzt*snth)
          any=-anxt*snps+csps*(anyt*csth+anzt*snth)
          anz= anzt*csth-anyt*snth
c
          an20=1.0d0+gr2
          an0=sqrt(an20)
          anx0=anx
          any0=any
          anz0=anz
c
          xc0(1,j,k)=x0
          xc0(2,j,k)=y0
          xc0(3,j,k)=z0
c
          ywrk0(1,j,k)=x0
          ywrk0(2,j,k)=y0
          ywrk0(3,j,k)=z0
          ywrk0(4,j,k)=anx0
          ywrk0(5,j,k)=any0
          ywrk0(6,j,k)=anz0
c
          ypwrk0(1,j,k) = anx0/an0
          ypwrk0(2,j,k) = any0/an0
          ypwrk0(3,j,k) = anz0/an0
          ypwrk0(4,j,k) = dgr2x/an0/2.0d0
          ypwrk0(5,j,k) = dgr2y/an0/2.0d0
          ypwrk0(6,j,k) = dgr2z/an0/2.0d0
c
          grad2(j,k)=gr2
          dgrad2v(1,j,k)=dgr2x
          dgrad2v(2,j,k)=dgr2y
          dgrad2v(3,j,k)=dgr2z
c
          du10(1,j,k)= du1tx*csps+snps*(du1ty*csth+du1tz*snth)
          du10(2,j,k)=-du1tx*snps+csps*(du1ty*csth+du1tz*snth)
          du10(3,j,k)= du1tz*csth-du1ty*snth
c
          dd=anx0**2+any0**2+anz0**2-an20
          vgradi=anxt*gxt+anyt*gyt+anzt*gzt
          ddi=2.0d0*vgradi
c
          if(j.eq.nrayr) then
            r0=sqrt(x0**2+y0**2)
            x0m=x0/1.0d2
            y0m=y0/1.0d2
            r0m=r0/1.0d2
            z0m=z0/1.0d2
            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
     .                    -1.0d0,zero,zero,zero,one
          end if
          if(j.eq.1.and.k.eq.1) then
            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
     .       zero,zero,zero,zero,
     .       zero,zero,zero,zero,one,zero,zero,
     .       zero,zero,one,zero,zero,zero,zero,one
            ddr110=dd
          end if  
          if(j.eq.nrayr.and.k.eq.1) then
            write(17,99) zero,ddr110,dd,ddi
          end if          
        end do
      end do
c
      call pweigth
c
      if(nrayr.gt.1) then
        iproj=0
        nfilp=8
        call projxyzt(iproj,nfilp)
      end if
c
      return
 99   format(24(1x,e16.8e3))
111   format(3i5,20(1x,e16.8e3))
      end
c
c     ray tracing initial conditions  igrad=0
c
      subroutine ic_rt
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,ndimm=3)
      parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
      parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
      dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
      dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
      dimension dffiu(jmx),ddffiu(jmx)
      dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
      common/nray/nrayr,nrayth
      common/rwmax/rwmax
      common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
      common/anic/anx0c,any0c,anz0c
      common/wrk/ywrk0,ypwrk0
      common/grc/xc0,du10
      common/fi/dffiu,ddffiu
      common/mirr/x00,y00,z00
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
      common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
c
      csth=anz0c
      snth=sqrt(1.0d0-csth**2)
      csps=1.0d0
      snps=0.0d0
      if(snth.gt.0.0d0) then
        csps=any0c/snth
        snps=anx0c/snth
      end if
c
      phiwrad=phiw*cvdr
      csphiw=cos(phiwrad)
      snphiw=sin(phiwrad)
c
      dr=1.0d0
      if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
      da=2.0d0*pi/dble(nrayth)
      z0t=0.0d0
c
      do j=1,nrayr
        u=dble(j-1)
        dffiu(j)=0.0d0
        ddffiu(j)=0.0d0
        do k=1,nrayth
          alfak=(k-1)*da
          dcsiw=dr*cos(alfak)*wcsi
          detaw=dr*sin(alfak)*weta
          dx0t=dcsiw*csphiw-detaw*snphiw
          dy0t=dcsiw*snphiw+detaw*csphiw
          x0t=u*dx0t
          y0t=u*dy0t
c
c          csiw=u*dcsiw
c          etaw=u*detaw
c          csir=csiw
c          etar=etaw
c
          dx0= x0t*csps+snps*(y0t*csth+z0t*snth)
          dy0=-x0t*snps+csps*(y0t*csth+z0t*snth)
          dz0= z0t*csth-y0t*snth
c
          x0=x00+dx0
          y0=y00+dy0
          z0=z00+dz0
c
          ppcsi=u*dr*cos(alfak)*rcicsi
          ppeta=u*dr*sin(alfak)*rcieta
c
          anzt=1.0d0/sqrt(1.0d0+ppcsi**2+ppeta**2)
          ancsi=ppcsi*anzt
          aneta=ppeta*anzt
c
          anxt=ancsi*csphiw-aneta*snphiw
          anyt=ancsi*snphiw+aneta*csphiw
c
          anx= anxt*csps+snps*(anyt*csth+anzt*snth)
          any=-anxt*snps+csps*(anyt*csth+anzt*snth)
          anz= anzt*csth-anyt*snth
c
          an20=1.0d0
          an0=sqrt(an20)
          anx0=anx
          any0=any
          anz0=anz
c
          xc0(1,j,k)=x0
          xc0(2,j,k)=y0
          xc0(3,j,k)=z0
c
          ywrk0(1,j,k)=x0
          ywrk0(2,j,k)=y0
          ywrk0(3,j,k)=z0
          ywrk0(4,j,k)=anx0
          ywrk0(5,j,k)=any0
          ywrk0(6,j,k)=anz0
c
          ypwrk0(1,j,k) = anx0/an0
          ypwrk0(2,j,k) = any0/an0
          ypwrk0(3,j,k) = anz0/an0
          ypwrk0(4,j,k) = 0.0d0
          ypwrk0(5,j,k) = 0.0d0
          ypwrk0(6,j,k) = 0.0d0
c
          do iv=1,3
            gri(iv,j,k)=0.0d0
            dgrad2v(iv,j,k)=0.0d0
            du10(iv,j,k)=0.0d0
            do jv=1,3
              ggri(iv,jv,j,k)=0.0d0
            end do
          end do
          grad2(j,k)=0.0d0
c
          dd=anx0**2+any0**2+anz0**2-an20
          vgradi=0.0d0
          ddi=2.0d0*vgradi
c
          if(j.eq.nrayr) then
            r0=sqrt(x0**2+y0**2)
            x0m=x0/1.0d2
            y0m=y0/1.0d2
            r0m=r0/1.0d2
            z0m=z0/1.0d2
            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
     .                    -1.0d0,zero,zero,zero,one
          end if
           if(j.eq.1.and.k.eq.1) then
            write(17,99) zero,zero,zero,zero
            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
     .       zero,zero,zero,zero,
     .       zero,zero,zero,zero,one,zero,zero,
     .       zero,zero,one,zero,zero,zero,zero,one
          end if  
        end do
      end do
c
      call pweigth
c
      if(nrayr.gt.1) then
        iproj=0
        nfilp=8
        call projxyzt(iproj,nfilp)
      end if
c
      return
 99   format(24(1x,e16.8e3))
111   format(3i5,20(1x,e16.8e3))
      end



      subroutine ic_rt2
      implicit real*8 (a-h,o-z)
      parameter(ndim=6,ndimm=3)
      parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
      parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
      dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
      dimension yyrfl(jmx,kmx,ndim)
      dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
      dimension dffiu(jmx),ddffiu(jmx)
      dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
      dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
c
      common/nray/nrayr,nrayth
      common/wrk/ywrk0,ypwrk0
      common/grc/xc0,du10
      common/grad2jk/grad2
      common/dgrad2jk/dgrad2v
      common/gradjk/gri
      common/ggradjk/ggri
      common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
      common/yyrfl/yyrfl

      do j=1,nrayr
        do k=1,nrayth
          x0=yyrfl(j,k,1)
          y0=yyrfl(j,k,2)
          z0=yyrfl(j,k,3)
          anx0=yyrfl(j,k,4)
          any0=yyrfl(j,k,5)
          anz0=yyrfl(j,k,6)
          an20=anx0*anx0+any0*any0+anz0*anz0
          an0=sqrt(an20)
c
          xc0(1,j,k)=x0
          xc0(2,j,k)=y0
          xc0(3,j,k)=z0
c
          ywrk0(1,j,k)=x0
          ywrk0(2,j,k)=y0
          ywrk0(3,j,k)=z0
          ywrk0(4,j,k)=anx0/an0
          ywrk0(5,j,k)=any0/an0
          ywrk0(6,j,k)=anz0/an0
c
          ypwrk0(1,j,k) = anx0/an0
          ypwrk0(2,j,k) = any0/an0
          ypwrk0(3,j,k) = anz0/an0
          ypwrk0(4,j,k) = 0.0d0
          ypwrk0(5,j,k) = 0.0d0
          ypwrk0(6,j,k) = 0.0d0
c
          do iv=1,3
            gri(iv,j,k)=0.0d0
            dgrad2v(iv,j,k)=0.0d0
            du10(iv,j,k)=0.0d0
            do jv=1,3
              ggri(iv,jv,j,k)=0.0d0
            end do
          end do
          grad2(j,k)=0.0d0
c
          dd=anx0**2+any0**2+anz0**2-an20
          vgradi=0.0d0
          ddi=2.0d0*vgradi
c
          if(j.eq.nrayr) then
            r0=sqrt(x0**2+y0**2)
            x0m=x0/1.0d2
            y0m=y0/1.0d2
            r0m=r0/1.0d2
            z0m=z0/1.0d2
            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
     .                    -1.0d0,zero,zero,zero,one
          end if
          if(j.eq.1.and.k.eq.1) then
            write(17,99) zero,zero,zero,zero
            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
     .       zero,zero,zero,zero,
     .       zero,zero,zero,zero,one,zero,zero,
     .       zero,zero,one,zero,zero,zero,zero,one
          end if  
        end do
      end do
c
      call pweigth
c
      if(nrayr.gt.1) then
        iproj=0
        nfilp=8
        call projxyzt(iproj,nfilp)
      end if
c
      return
 99   format(24(1x,e16.8e3))
111   format(3i5,20(1x,e16.8e3))
      end
c
c     ray power weigth coefficient q(j)
c
      subroutine pweigth
      implicit real*8(a-h,o-z)
      parameter(jmx=31)
      dimension q(jmx)
      common/qw/q
      common/nray/nrayr,nrayth
      common/rwmax/rwmax
c
      dr=1.0d0
      if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
      r1=0.0d0
      summ=0.0d0
      q(1)=1.0d0
      if(nrayr.gt.1) then
        do j=1,nrayr
            r2=(dble(j)-0.5d0)*dr
          q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2))
          r1=r2
          summ=summ+q(j)
        end do
c
        q(1)=q(1)/summ
        sm=q(1)
        j=1
        k=1
        do j=2,nrayr
          q(j)=q(j)/nrayth/summ
          do k=1,nrayth
            sm=sm+q(j)
          end do
        end do
      end if
c
      return
      end
c
c
c
      subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi,
     .                     bmxi,bmni,fci,intp)
      implicit real*8 (a-h,o-z)
      parameter(nnintp=101)
      dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4)
      dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4)
      dimension carea(nnintp,4),cfc(nnintp,4)
c
      common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc
      common/pstab/rpstab
      common/eqnn/nr,nz,npp,nintp
c
      ip=int((nintp-1)*rpsi+1)
      if(ip.eq.0) ip=1
      if(ip.eq.nintp) ip=nintp-1
c
      dps=rpsi-rpstab(ip)
c
      areai=spli(carea,nintp,ip,dps)
      voli=spli(cvol,nintp,ip,dps)
      dervoli=splid(cvol,nintp,ip,dps)
      rrii=spli(crri,nintp,ip,dps)
c
      if(intp.eq.0) return
c
      rbavi=spli(crbav,nintp,ip,dps)
      bmxi=spli(cbmx,nintp,ip,dps)
      bmni=spli(cbmn,nintp,ip,dps)
      fci=spli(cfc,nintp,ip,dps)
c
      return
      end
c
c
c
      subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli)
      implicit real*8 (a-h,o-z)
      parameter(nnintp=101)
      dimension rpstab(nnintp)
      dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4)
      common/pstab/rpstab
      common/eqnn/nr,nz,npp,nintp
      common/cratj/cratja,cratjb,cratjpl
      ip=int((nintp-1)*rpsi+1)
      if(ip.eq.0) ip=1
      if(ip.eq.nintp) ip=nintp-1
      dps=rpsi-rpstab(ip)
      ratjai=spli(cratja,nintp,ip,dps)
      ratjbi=spli(cratjb,nintp,ip,dps)
      ratjpli=spli(cratjpl,nintp,ip,dps)
      return
      end
c
c
c
      subroutine pabs_curr(i,j,k)
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36,nmx=8000)
      parameter(pi=3.14159265358979d0)
      dimension psjki(jmx,kmx,nmx)
      dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
      dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension q(jmx),tau1v(jmx,kmx)
c
      common/psjki/psjki
      common/atjki/tauv,alphav
      common/dpjjki/pdjki,currj
      common/pcjki/ppabs,ccci
      common/dcjki/didst
      common/tau1v/tau1v
c
      common/qw/q
      common/p0/p0mw
c
      common/dsds/dst
      common/dersdst/dersdst
      common/iieq/iequil
c
      common/parban/b0,rr0m,zr0m,rpam
      common/absor/alpha,effjcd,akim,tau0
c
      common/psival/psinv
      common/sgnib/sgnbphi,sgniphi
      common/bmxmn/bmxi,bmni
      common/fc/fci
c
      dvvl=1.0d0
      rbavi=1.0d0
      rrii=rr0m
      rhop0=sqrt(psjki(j,k,i-1))
      if(iequil.eq.2) then
        intp=1
        psinv=psjki(j,k,i)
        rhop=sqrt(psinv)
        call valpsispl(rhop,voli,dervoli,area,rrii,
     .              rbavi,bmxi,bmni,fci,intp)
        dvvl=dervoli*abs(rhop-rhop0)
      else
        dvvl=4.0d0*pi*rr0m*pi*abs(rhop*rhop-rhop0*rhop0)*rpam**2
        rbavi=b0/bmni
        fci=1.0d0-1.469d0*sqrt(rpam/rr0m)
      end if
      if(dvvl.le.0.0d0) dvvl=1.0d0
c
      adnm=2.0d0*pi*rrii
c
c       calcolo della corrente cohen: currtot [MA]
c       calcolo della densita' corrente cohen: currj  [MA/m^2]
c       calcolo della densita' potenza:  pdjki       [MW/m^3]
c       calcolo della efficienza <j/p>:  effjcdav      [A m/W ]
c
      tau0=tauv(j,k,i-1)
c
      call ecrh_cd
c
      alphav(j,k,i)=alpha
      dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0
      tauv(j,k,i)=tauv(j,k,i-1)+dtau
      dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))*
     . alphav(j,k,i)*dersdst
      pdjki(j,k,i)=dpdst*dst/dvvl
      ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))*
     . (1.0d0-exp(-tauv(j,k,i)))
      effjcdav=rbavi*effjcd
      currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i)
      didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm
      dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0d0
      ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst
      return
      end
c
c
c
      subroutine ecrh_cd
      implicit real*8 (a-h,o-z)
      real*8 mc2,mesi
      parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8)
      parameter(qesi=1.602176487d-19,mesi=9.10938215d-31)
      parameter(vcsi=2.99792458d+8)
      parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3)
c
      common/ithn/ithn
      common/nharm/nharm,nhf
      common/warm/iwarm,ilarm
      common/ieccd/ieccd
c
      common/ygyg/yg
      common/nplr/anpl,anpr
      common/vgv/vgm,derdnm
      common/parwv/ak0,akinv,fhz
c
      common/nprw/anprre,anprim
c
      common/absor/alpha,effjcd,akim,tau
c
      common/psival/psinv
      common/tete/tekev
      common/amut/amu
      common/zz/Zeff
c
c     absorption computation
c
      alpha=0.0d0
      akim=0.0d0
      effjcd=0.0d0
c
      tekev=temp(psinv)
      amu=mc2/tekev
      zeff=fzeff(psinv)
c
      if(tekev.le.0.0d0.or.tau.gt.taucr) return
c
      dnl=1.0d0-anpl*anpl
      if(iwarm.eq.1) then
        ygc=1.0d0-anpl**2/2.0d0
      else
        ygc=sqrt(1.0d0-anpl**2)
      end if
c
      nharm=int(ygc/yg-eps)+1
c
      if(nharm.eq.0) return
c
      nhf=0
      do nn=nharm,nharm+4
        ygn=nn*yg
        if(iwarm.eq.1) then
          rdu2=2.0d0*(ygn-ygc)
          u1=anpl+sqrt(rdu2)
          u2=anpl-sqrt(rdu2)
          uu2=min(u1*u1,u2*u2)
          argex=amu*uu2/2.0d0
        else
          rdu2=ygn**2-ygc**2
          g1=(ygn+anpl*sqrt(rdu2))/dnl
          g2=(ygn-anpl*sqrt(rdu2))/dnl
          gg=min(g1,g2)
          argex=amu*(gg-1.0d0)
          u1=(ygn*anpl+sqrt(rdu2))/dnl
          u2=(ygn*anpl-sqrt(rdu2))/dnl
        end if
        if(argex.le.xxcr) nhf=nn
      end do
      if(nhf.eq.0) return
c
      lrm=ilarm
      if(lrm.lt.nhf) lrm=nhf
      call dispersion(lrm)
c
      akim=ak0*anprim
      ratiovgr=2.0d0*anpr/derdnm
c      ratiovgr=2.0d0*anpr/derdnm*vgm
      alpha=2.0d0*akim*ratiovgr
      if(alpha.lt.0.0d0) then
        ierr=94
        print*,'    IERR = ', ierr,'  alpha negative'
      end if
c
      ithn=1
      if(lrm.gt.nharm) ithn=2
      if(ieccd.gt.0) call eccd(effjcd)

      return
999   format(12(1x,e12.5))
      end
c
c
c
      subroutine dispersion(lrm)
      implicit real*8(a-h,o-z)
      parameter(imx=20)
      complex*16 cc0,cc2,cc4,rr
      complex*16 anpr2a,anpra
      complex*16 anpr,anpr2,ex,ey,ez,den
      complex*16 epsl(3,3,lrm),sepsl(3,3),e330
      complex*16 e11,e22,e12,e13,e23
c      complex*16 e33,e21,e31,e32
      complex*16 a13,a31,a23,a32,a33
c
      common/ygyg/yg
      common/nplr/anpl,anprf
c
      common/mode/sox
      common/warm/iwarm,ilarm
c
      common/nprw/anprr,anpri
      common/epolar/ex,ey,ez
      common/amut/amu
c
      errnpr=1.0d0
      anpr2a=anprf**2
      anpl2=anpl*anpl
c
      if (iwarm.eq.1) then
        call diel_tens_wr(e330,epsl,lrm)
      else
        call diel_tens_fr(e330,epsl,lrm)
      end if
c
      do i=1,imx
c
        do j=1,3
          do k=1,3
            sepsl(k,j)=dcmplx(0.0d0,0.0d0)
            do ilrm=1,lrm
              sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*anpr2a**(ilrm-1)
            end do
          end do
        end do
c
        anpra=sqrt(anpr2a)
c
        e11=sepsl(1,1)
        e22=sepsl(2,2)
        e12=sepsl(1,2)
        a33=sepsl(3,3)
        a13=sepsl(1,3)
        a23=sepsl(2,3)
        a31=a13
        a32=-a23
c        e33=e330+anpr2a*a33
        e13=anpra*a13
        e23=anpra*a23
c        e21=-e12
c        e31=e13
c        e32=-e23
c
        if(i.gt.2.and.errnpr.lt.1.0d-3) go to 999
c
        cc4=(e11-anpl2)*(1.0d0-a33)+(a13+anpl)*(a31+anpl)
        cc2=-e12*e12*(1.0d0-a33)
     .      -a32*e12*(a13+anpl)+a23*e12*(a31+anpl)
     .      -(a23*a32+e330+(e22-anpl2)*(1.0d0-a33))*(e11-anpl2)
     .      -(a13+anpl)*(a31+anpl)*(e22-anpl2)
        cc0=e330*((e11-anpl2)*(e22-anpl2)+e12*e12)
c
        rr=cc2*cc2-4.0d0*cc0*cc4
c
        if(yg.gt.1.0d0) then
          s=sox
          if(dimag(rr).LE.0.0d0) s=-s
        else
          s=-sox
          if(dble(rr).le.0.d0.and.dimag(rr).ge.0.d0) s=-s
        end if
c
        anpr2=(-cc2+s*sqrt(rr))/(2.0d0*cc4)
c
        if(dble(anpr2).lt.0.0d0.and.dimag(anpr2).lt.0.0d0) then
          anpr2=0.0d0
          print*,'    Y =',yg,'   nperp2 < 0'
c          ierr=99
          go to 999
        end if
c
        errnpr=abs(1.0d0-abs(anpr2)/abs(anpr2a))
        anpr2a=anpr2
      end do
c
  999 continue
       if(i.gt.imx) print*,'    i>imx ',yg,errnpr,i
c
      anpr=sqrt(anpr2)
      anprr=dble(anpr)
      anpri=dimag(anpr)
99    format(20(1x,e12.5))     
c
      ex=dcmplx(0.0d0,0.0d0)
      ey=dcmplx(0.0d0,0.0d0)
      ez=dcmplx(0.0d0,0.0d0)
c
      if (abs(anpl).gt.1.0d-6) then
        den=e12*e23-(e13+anpr*anpl)*(e22-anpr2-anpl2)
        ey=-(e12*(e13+anpr*anpl)+(e11-anpl2)*e23)/den
        ez=(e12*e12+(e22-anpr2-anpl2)*(e11-anpl2))/den
        ez2=abs(ez)**2
        ey2=abs(ey)**2
        enx2=1.0d0/(1.0d0+ey2+ez2)
        ex=dcmplx(sqrt(enx2),0.0d0)
        ez2=ez2*enx2
        ey2=ey2*enx2
        ez=ez*ex
        ey=ey*ex
      else
        if(sox.lt.0.0d0) then
          ez=dcmplx(1.0d0,0.0d0)
          ez2=abs(ez)**2
        else
          ex2=1.0d0/(1.0d0+abs(-e11/e12)**2)
          ex=sqrt(ex2)
          ey=-ex*e11/e12
          ey2=abs(ey)**2
          ez2=0.0d0
        end if
      end if
c
      return
      end
c
c     Fully relativistic case
c     computation of dielectric tensor elements
c     up to third order in Larmor radius for hermitian part
c
      subroutine diel_tens_fr(e330,epsl,lrm)
      implicit real*8(a-h,o-z)
      complex*16 ui
      complex*16 e330,epsl(3,3,lrm)
      complex*16 ca11,ca12,ca22,ca13,ca23,ca33
      complex*16 cq0p,cq0m,cq1p,cq1m,cq2p
      parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0)
      parameter(ui=(0.0d0,1.0d0))
      dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm)
c
      common/xgxg/xg
      common/ygyg/yg
      common/amut/amu
      common/nplr/anpl,anpr
      common/resah/anpl2,dnl
c
      common/cri/cr,ci
      common/warm/iwarm,ilarm
c
      anpl2=anpl**2
      dnl=1.0d0-anpl2
c
      cmxw=1.0d0+15.0d0/(8.0d0*amu)+105.0d0/(128.0d0*amu**2)
      cr=-amu*amu/(rpi*cmxw)
      ci=sqrt(2.0d0*pi*amu)*amu**2/cmxw
c
      do l=1,lrm
        do j=1,3
          do i=1,3
            epsl(i,j,l)=dcmplx(0.0d0,0.0d0)
          end do
        end do
      end do
c
      if(iwarm.eq.2) call hermitian(rr,lrm)
      if(iwarm.eq.4) call hermitian_2(rr,lrm)
c
      call antihermitian(ri,lrm)
c
      do l=1,lrm
        lm=l-1
        fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm))
        ca11=(0.0d0,0.0d0)
        ca12=(0.0d0,0.0d0)
        ca13=(0.0d0,0.0d0)
        ca22=(0.0d0,0.0d0)
        ca23=(0.0d0,0.0d0)
        ca33=(0.0d0,0.0d0)
        do is=0,l
          k=l-is
          w=(-1.0d0)**k
c
          asl=w/(fact(is+l)*fact(l-is))
          bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0))
c
          if(is.gt.0) then
          cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l)
          cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l)
          cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l)
          cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l)
          cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l)
          else
          cq0p=rr(is,0,l)
          cq0m=rr(is,0,l)
          cq1p=rr(is,1,l)
          cq1m=rr(is,1,l)
          cq2p=rr(is,2,l)
          end if
c
          ca11=ca11+is**2*asl*cq0p
          ca12=ca12+is*l*asl*cq0m
          ca22=ca22+bsl*cq0p
          ca13=ca13+is*asl*cq1m/yg
          ca23=ca23+l*asl*cq1p/yg
          ca33=ca33+asl*cq2p/yg**2
        end do
        epsl(1,1,l) =  - xg*ca11*fal
        epsl(1,2,l) =  + ui*xg*ca12*fal
        epsl(2,2,l) =  - xg*ca22*fal
        epsl(1,3,l) =  - xg*ca13*fal
        epsl(2,3,l) =  - ui*xg*ca23*fal
        epsl(3,3,l) =  - xg*ca33*fal
      end do
c
      cq2p=rr(0,2,0)
      e330=1.0d0+xg*cq2p
c
      epsl(1,1,1) = 1.d0 + epsl(1,1,1)
      epsl(2,2,1) = 1.d0 + epsl(2,2,1)
c
      do l=1,lrm
        epsl(2,1,l) = - epsl(1,2,l)
        epsl(3,1,l) =   epsl(1,3,l)
        epsl(3,2,l) = - epsl(2,3,l)
      end do
c
      return
      end
c
c
c
      subroutine hermitian(rr,lrm)
      implicit real*8(a-h,o-z)
      parameter(tmax=5,npts=500)
      dimension rr(-lrm:lrm,0:2,0:lrm)
      real*8 ttv(npts+1),extv(npts+1)
c
      common/ttex/ttv,extv
c
      common/ygyg/yg
      common/amut/amu
      common/nplr/anpl,anpr
      common/warm/iwarm,ilarm
      common/cri/cr,ci
c
      do n=-lrm,lrm
        do k=0,2
          do m=0,lrm
            rr(n,k,m)=0.0d0
          end do
        end do
      end do
c
      llm=min(3,lrm)
c
      dt=2.0d0*tmax/dble(npts)
      bth2=2.0d0/amu
      bth=sqrt(bth2)
      amu2=amu*amu
      amu4=amu2*amu2
      amu6=amu4*amu2
c
      do i = 1, npts+1
        t = ttv(i)
        rxt=sqrt(1.0d0+t*t/(2.0d0*amu))
        x = t*rxt
        upl2=bth2*x**2
        upl=bth*x
        gx=1.0d0+t*t/amu
        exdx=cr*extv(i)*gx/rxt*dt
c
        n1=1
        if(iwarm.gt.2) n1=-llm
c
        do n=n1,llm
          nn=abs(n)
          gr=anpl*upl+n*yg
          zm=-amu*(gx-gr)
          s=amu*(gx+gr)
          zm2=zm*zm
          zm3=zm2*zm
          call calcei3(zm,fe0m)
c
          do m=nn,llm
            if(n.eq.0.and.m.eq.0) then
            rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2
            else
            ffe=0.0d0
            if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))/amu2
            if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s
     .                      +s*s*(1.0d0+zm-zm2*fe0m))/amu4
            if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*
     .       (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/amu6
c
            rr(n,0,m) = rr(n,0,m) + exdx*ffe
            rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl
            rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2
            end if
c
          end do
        end do
      end do
c
      if(iwarm.gt.2) return
c
      sy1=1.0d0+yg
      sy2=1.0d0+yg*2.0d0
      sy3=1.0d0+yg*3.0d0
c
      bth4=bth2*bth2
      bth6=bth4*bth2
c
      anpl2=anpl*anpl
c
      rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2)
     .          +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2))
      rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2))
      rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2))
      rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2
     .             /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/
     .         sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3))
      rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+
     .             1.5d0*anpl2/sy1**2)) 
      rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*
     .              anpl2/sy1**2))
c
      if(llm.gt.1) then
c
      rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+
     .            bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2))
      rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2))
      rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2))
      rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2*
     .              (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4*
     .    (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2
     .    -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))     
      rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2*
     .              (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2*
     .              (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2*
     .              (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4*
     .    (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2
     .    -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))    
      rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2*
     .              (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2))
      rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2*
     .               (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2))
c
      if(llm.gt.2) then
c
      rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4*
     .            (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) 
      rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2))
      rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2))
      rr(-1,0,3) = -12.0d0*bth4/sy1*
     .             (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+
     .             bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2)
     .        /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
      rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2*
     .             (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-1,2,3) = -6.0d0*bth6/sy1*
     .             (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2))
c
      rr(-2,0,3) = -12.0d0*bth4/sy2*
     .              (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+
     .              bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2)
     .        /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
      rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2*
     .              (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2))
      rr(-2,2,3) = -6.0d0*bth6/sy2*
     .              (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2))
c
      rr(-3,0,3) = -12.0d0*bth4/sy3*
     .              (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+
     .              bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2)
     .        /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4))
      rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2*
     .              (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2))
      rr(-3,2,3) = -6.0d0*bth6/sy3*
     .              (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2))
c
      end if
      end if
c
      return
      end



      subroutine hermitian_2(rr,lrm)
      implicit real*8(a-h,o-z)
      parameter(tmax=5,npts=500)
      dimension rr(-lrm:lrm,0:2,0:lrm)
      parameter(epsa=0.0d0,epsr=1.0d-4)
      parameter (lw=5000,liw=lw/4)
      dimension w(lw),iw(liw)
      external fhermit
c
      common/ygyg/yg
      common/amut/amu
      common/nplr/anpl,anpr
      common/warm/iwarm,ilarm
      common/cri/cr,ci
      common/nmhermit/n,m,ih
c
      do n=-lrm,lrm
        do k=0,2
          do m=0,lrm
            rr(n,k,m)=0.0d0
          end do
        end do
      end do
c
      llm=min(3,lrm)
c
      bth2=2.0d0/amu
      bth=sqrt(bth2)
      amu2=amu*amu
      amu4=amu2*amu2
      amu6=amu4*amu2
c
        n1=1
        if(iwarm.gt.10) n1=-llm
c
        do n=n1,llm
          nn=abs(n)
          do m=nn,llm
            if(n.eq.0.and.m.eq.0) then
            ih=2
c              call dqagi(fhermit,bound,2,epsa,epsr,resfh,
              call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh,
     .             epp,neval,ier,liw,lw,last,iw,w)
              rr(n,2,m) = resfh
              else
              do ih=0,2
c              call dqagi(fhermit,bound,2,epsa,epsr,resfh,
              call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh,
     .             epp,neval,ier,liw,lw,last,iw,w)
                rr(n,ih,m) = resfh
              end do
            end if
          end do
        end do
c      write(83,'(12(1x,e12.5))') 
c     .       yg,rr(1,0,1),rr(1,1,1),rr(1,2,1)
      
      if(iwarm.gt.10) return
c
      sy1=1.0d0+yg
      sy2=1.0d0+yg*2.0d0
      sy3=1.0d0+yg*3.0d0
c
      bth4=bth2*bth2
      bth6=bth4*bth2
c
      anpl2=anpl*anpl
c
      rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2)
     .          +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2))
      rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2))
      rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2))
      rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2
     .             /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/
     .         sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3))
      rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+
     .             1.5d0*anpl2/sy1**2)) 
      rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*
     .              anpl2/sy1**2))
c
      if(llm.gt.1) then
c
      rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+
     .            bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2))
      rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2))
      rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2))
      rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2*
     .              (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4*
     .    (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2
     .    -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))     
      rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2*
     .              (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2*
     .              (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2*
     .              (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4*
     .    (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2
     .    -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))    
      rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2*
     .              (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2))
      rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2*
     .               (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2))
c
      if(llm.gt.2) then
c
      rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4*
     .            (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) 
      rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2))
      rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2))
      rr(-1,0,3) = -12.0d0*bth4/sy1*
     .             (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+
     .             bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2)
     .        /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4))
      rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2*
     .             (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2))
      rr(-1,2,3) = -6.0d0*bth6/sy1*
     .             (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2))
c
      rr(-2,0,3) = -12.0d0*bth4/sy2*
     .              (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+
     .              bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2)
     .        /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4))
      rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2*
     .              (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2))
      rr(-2,2,3) = -6.0d0*bth6/sy2*
     .              (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2))
c
      rr(-3,0,3) = -12.0d0*bth4/sy3*
     .              (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+
     .              bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2)
     .        /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4))
      rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2*
     .              (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2))
      rr(-3,2,3) = -6.0d0*bth6/sy3*
     .              (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2))
c
      end if
      end if
c
      return
      end



      function fhermit(t)
      implicit real*8 (a-h,o-z)
c
      common/ygyg/yg
      common/amut/amu
      common/nplr/anpl,anpr
      common/cri/cr,ci
      common/nmhermit/n,m,ih
      
      bth2=2.0d0/amu
      bth=sqrt(bth2)
      amu2=amu*amu
      amu4=amu2*amu2
      amu6=amu4*amu2
     
      rxt=sqrt(1.0d0+t*t/(2.0d0*amu))
      x = t*rxt
      upl2=bth2*x**2
      upl=bth*x
      gx=1.0d0+t*t/amu
      exdxdt=cr*exp(-t*t)*gx/rxt
      nn=abs(n)
      gr=anpl*upl+n*yg
      zm=-amu*(gx-gr)
      s=amu*(gx+gr)
      zm2=zm*zm
      zm3=zm2*zm
      call calcei3(zm,fe0m)
      ffe=0.0d0
      uplh=upl**ih
      if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2
      if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2
      if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s
     .                 +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4
      if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2)
     .                +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6
      fhermit= exdxdt*ffe
      return
      end
c
c
c
      subroutine antihermitian(ri,lrm)
      implicit none
      integer lmx,nmx,lrm,n,k,m,mm
      real*8 rpi
      parameter(rpi=1.77245385090552d0)
      parameter(lmx=20,nmx=lmx+2)
      real*8 fsbi(nmx)
      real*8 ri(lrm,0:2,lrm)
      real*8 yg,amu,anpl,cmu,anpl2,dnl,cr,ci,ygn,rdu2,rdu,anpr
      real*8 du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim
      real*8 fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0
      real*8 fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m
      real*8 fact
c
      common/ygyg/yg
      common/amut/amu
      common/nplr/anpl,anpr
      common/uu/ub,cmu
      common/resah/anpl2,dnl
c
      common/cri/cr,ci
c
      do n=1,lrm
        do k=0,2
          do m=1,lrm
            ri(n,k,m)=0.0d0
          end do
        end do
      end do
c
      dnl=1.0d0-anpl2
      cmu=anpl*amu
c
      do n=1,lrm
        ygn=n*yg
        rdu2=ygn**2-dnl
        if(rdu2.gt.0.0d0) then
          rdu=sqrt(rdu2)
            du=rdu/dnl
            ub=anpl*ygn/dnl
          aa=amu*anpl*du
          if (abs(aa).gt.5.0d0) then
            up=ub+du
            um=ub-du
            gp=anpl*up+ygn
            gm=anpl*um+ygn
            xp=up+1.0d0/cmu
            xm=um+1.0d0/cmu
            eem=exp(-amu*(gm-1.0d0))
            eep=exp(-amu*(gp-1.0d0))
            fi0p0=-1.0d0/cmu
            fi1p0=-xp/cmu
            fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu
            fi0m0=-1.0d0/cmu
            fi1m0=-xm/cmu
            fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu
            do m=1,lrm
            fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu
            fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu
            fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0
     .               +up*um*fi0p0)/cmu
            fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0
     .               +up*um*fi0m0)/cmu
              fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m*
     .             (ub*fi2p0-up*um*fi1p0))/cmu
              fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m*
     .             (ub*fi2m0-up*um*fi1m0))/cmu
            if(m.ge.n) then
              ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem)
              ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem)
              ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem)
            end if
            fi0p0=fi0p1
            fi1p0=fi1p1
            fi2p0=fi2p1
            fi0m0=fi0m1
            fi1m0=fi1m1
            fi2m0=fi2m1
            end do
          else
            ee=exp(-amu*(ygn-1.0d0+anpl*ub))
            call ssbi(aa,n,lrm,fsbi)
            do m=n,lrm
              cm=rpi*fact(m)*du**(2*m+1)
              cim=0.5d0*ci*dnl**m
              mm=m-n+1
              fi0m=cm*fsbi(mm)
              fi1m=-0.5d0*aa*cm*fsbi(mm+1)
              fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2))
              ri(n,0,m)=cim*ee*fi0m
              ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m)
              ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m)
            end do
          end if
        end if
      end do
c
      return
      end
c
      subroutine ssbi(zz,n,l,fsbi)
      implicit none
      integer n,l,nmx,lmx,k,m,mm
      real*8 eps,c0,c1,sbi,zz
      real*8 gamm
      parameter(eps=1.0d-10,lmx=20,nmx=lmx+2)
      real*8 fsbi(nmx)
      do m=n,l+2
        c0=1.0d0/gamm(dble(m)+1.5d0)
        sbi=c0
        do k=1,50
          c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k)
          sbi=sbi+c1
          if(c1/sbi.lt.eps) exit
          c0=c1
        end do
        mm=m-n+1
        fsbi(mm)=sbi
      end do
      return
      end
c
c     Weakly relativistic dielectric tensor
c     computation of dielectric tensor elements
c     Krivenki and Orefice, JPP 30,125 (1983)
c
      subroutine diel_tens_wr(ce330,cepsl,lrm)
      implicit real*8 (a-b,d-h,o-z)
      implicit complex*16 (c)
      dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm)
      parameter(cui=(0.0d0,1.0d0))
c
      common/xgxg/xg
      common/ygyg/yg
      common/nplr/anpl,anprf
      common/amut/amu
c
      anpl2=anpl*anpl
c
      call fsup(cefp,cefm,lrm)
c
      do l=1,lrm
        lm=l-1
        fcl=0.5d0**l*((1.0d0/yg)**2/amu)**lm*fact(2*l)/fact(l)
        ca11=(0.d0,0.d0)
        ca12=(0.d0,0.d0)
        ca13=(0.d0,0.d0)
        ca22=(0.d0,0.d0)
        ca23=(0.d0,0.d0)
        ca33=(0.d0,0.d0)
        do is=0,l
          k=l-is
          w=(-1.0d0)**k
c
          asl=w/(fact(is+l)*fact(l-is))
          bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0))
c
          cq0p=amu*cefp(is,0)
          cq0m=amu*cefm(is,0)
          cq1p=amu*anpl*(cefp(is,0)-cefp(is,1))
          cq1m=amu*anpl*(cefm(is,0)-cefm(is,1))
          cq2p=cefp(is,1)+amu*anpl2*(cefp(is,2)
     .        +cefp(is,0)-2.0d0*cefp(is,1))
c
          ca11=ca11+is**2*asl*cq0p
          ca12=ca12+is*l*asl*cq0m
          ca22=ca22+bsl*cq0p
          ca13=ca13+is*asl*cq1m/yg
          ca23=ca23+l*asl*cq1p/yg
          ca33=ca33+asl*cq2p/yg**2
        end do
        cepsl(1,1,l) =  - xg*ca11*fcl
        cepsl(1,2,l) =  + cui*xg*ca12*fcl
        cepsl(2,2,l) =  - xg*ca22*fcl
        cepsl(1,3,l) =  - xg*ca13*fcl
        cepsl(2,3,l) =  - cui*xg*ca23*fcl
        cepsl(3,3,l) =  - xg*ca33*fcl
      end do
c
      cq2p=cefp(0,1)+amu*anpl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1))
      ce330=1.0d0-xg*amu*cq2p
c
      cepsl(1,1,1) = 1.d0 + cepsl(1,1,1)
      cepsl(2,2,1) = 1.d0 + cepsl(2,2,1)
c
      do l=1,lrm
        cepsl(2,1,l) = - cepsl(1,2,l)
        cepsl(3,1,l) =   cepsl(1,3,l)
        cepsl(3,2,l) = - cepsl(2,3,l)
      end do
c
      return
      end
c
c
c
      subroutine fsup(cefp,cefm,lrm)
      implicit real*8 (a-b,d-h,o-z)
      implicit complex*16 (c)
      parameter(apsicr=0.7d0)
      parameter(cui=(0.0d0,1.0d0))
      dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2)
c
      common/ygyg/yg
      common/nplr/anpl,anprf
      common/amut/amu
c
      psi=sqrt(0.5d0*amu)*anpl
      apsi=abs(psi)
c
      do is=0,lrm
        alpha=anpl*anpl/2.0d0+is*yg-1.0d0
        phi2=amu*alpha
        phim=sqrt(abs(phi2))
        if (alpha.ge.0) then
          xp=psi-phim
          yp=0.0d0
          xm=-psi-phim
          ym=0.0d0
          x0=-phim
          y0=0.0d0
        else
          xp=psi
          yp=phim
          xm=-psi
          ym=phim
          x0=0.0d0
          y0=phim
        end if
        call zetac (xp,yp,zrp,zip,iflag)
        call zetac (xm,ym,zrm,zim,iflag)
c
        czp=dcmplx(zrp,zip)
        czm=dcmplx(zrm,zim)
        cf12=(0.0d0,0.0d0)
        if (alpha.ge.0) then
          if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim)
        else
          cf12=-cui*(czp+czm)/(2.0d0*phim)
        end if
c
        if(apsi.gt.apsicr) then
          cf32=-(czp-czm)/(2.0d0*psi)
        else
          cphi=phim
          if(alpha.lt.0) cphi=-cui*phim
          call zetac (x0,y0,zr0,zi0,iflag)
          cz0=dcmplx(zr0,zi0)
          cdz0=2.0d0*(1.0d0-cphi*cz0)
          cf32=cdz0
        end if
c
        cf0=cf12
        cf1=cf32
        cefp(is,0)=cf32
        cefm(is,0)=cf32
        do l=1,is+2
          iq=l-1
          if(apsi.gt.apsicr) then
            cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2
          else
            cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0)
          end if
          ir=l-is
          if(ir.ge.0) then
            cefp(is,ir)=cf2
            cefm(is,ir)=cf2
          end if
          cf0=cf1
          cf1=cf2
        end do
c
        if(is.ne.0) then
c
        alpha=anpl*anpl/2.0d0-is*yg-1.0d0
        phi2=amu*alpha
        phim=sqrt(abs(phi2))
        if (alpha.ge.0.0d0) then
          xp=psi-phim
          yp=0.0d0
          xm=-psi-phim
          ym=0.0d0
          x0=-phim
          y0=0.0d0
        else
          xp=psi
          yp=phim
          xm=-psi
          ym=phim
          x0=0.0d0
          y0=phim
        end if
        call zetac (xp,yp,zrp,zip,iflag)
        call zetac (xm,ym,zrm,zim,iflag)
c
        czp=dcmplx(zrp,zip)
        czm=dcmplx(zrm,zim)
c
        cf12=(0.0d0,0.0d0)
        if (alpha.ge.0) then
          if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim)
        else
          cf12=-cui*(czp+czm)/(2.0d0*phim)
        end if
        if(apsi.gt.apsicr) then
          cf32=-(czp-czm)/(2.0d0*psi)
        else
          cphi=phim
          if(alpha.lt.0) cphi=-cui*phim
          call zetac (x0,y0,zr0,zi0,iflag)
          cz0=dcmplx(zr0,zi0)
          cdz0=2.0d0*(1.0d0-cphi*cz0)
          cf32=cdz0
        end if
c
        cf0=cf12
        cf1=cf32
        do l=1,is+2
          iq=l-1
          if(apsi.gt.apsicr) then
            cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2
          else
            cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0)
          end if
          ir=l-is
          if(ir.ge.0) then
            cefp(is,ir)=cefp(is,ir)+cf2
            cefm(is,ir)=cefm(is,ir)-cf2
          end if
          cf0=cf1
          cf1=cf2
        end do
c
        end if
c
      end do
c
      return
      end
      
      
      
      subroutine eccd(effjcd)
      use green_func_p    
      implicit none
      real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev
      real*8 anucc,canucc,ddens
      real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff
      real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc
      real*8 fjch,fjncl,fjch0,fconic
      real*8 cst,cst2
      integer ieccd,nharm,nhf,nhn
      external fjch,fjncl,fjch0

      parameter(qesi=1.602176487d-19,mesi=9.10938215d-31)
      parameter(vcsi=2.99792458d+8)
      parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2)
      parameter(canucc=4.0d13*pi*qe**4/(me**2*vc**3))
      parameter(ceff=qesi/(mesi*vcsi))
      
      common/nharm/nharm,nhf
      common/nhn/nhn
      
      common/ieccd/ieccd
      common/tete/tekev
      common/dens/dens,ddens
      common/zz/Zeff
      common/btot/btot
      common/bmxmn/bmax,bmin
      common/fc/fc
      common/ncl/rbx
      common/cohen/rbn,alams,pa,fp0s
      common/cst/cst,cst2
      
      anum=0.0d0
      denom=0.0d0
      effjcd=0.0d0

      coullog=24.0d0-log(1.0d4*sqrt(0.1d0*dens)/tekev)
      anucc=canucc*dens*coullog
      
c      nhf=nharm
      
      select case(ieccd) 
      
      case(1)
c     cohen model      
c     rbn=B/B_min
c     rbx=B/B_max
c     cst2=1.0d0-B/B_max
c     alams=sqrt(1-B_min/B_max)
c     Zeff < 31 !!!
c     fp0s= P_a (alams)   
        rbn=btot/bmin
        rbx=btot/bmax
        cst2=1.0d0-rbx
        if(cst2.lt.1d-6) cst2=0.0d0
        cst=sqrt(cst2)
        alams=sqrt(1.0d0-bmin/bmax)
        pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0
        fp0s=fconic(alams,pa,0)
        do nhn=nharm,nhf
          call curr_int(fjch,resj,resp)
          anum=anum+resj
          denom=denom+resp
        end do
       
      case(2:9)
        cst=0.0d0
        cst2=0.0d0
        do nhn=nharm,nhf
          call curr_int(fjch0,resj,resp)
          anum=anum+resj
          denom=denom+resp
        end do
     
      case(10:11)
c     neoclassical model:
c     ft=1-fc trapped particle fraction
c     rzfc=(1+Zeff)/fc     
        rbn=btot/bmin
        rbx=btot/bmax
        cst2=1.0d0-rbx
        if(cst2.lt.1d-6) cst2=0.0d0
        cst=sqrt(cst2)
        call SpitzFuncCoeff(Tekev,Zeff,fc)
        do nhn=nharm,nhf
          call curr_int(fjncl,resj,resp)
          anum=anum+resj
          denom=denom+resp
      end do
      nhn=nhn-1
       
       CASE DEFAULT
       print*,'ieccd undefined'
       
      end select
c
c     effjpl = <J_parallel>/<p_d> /(B_min/<B>)  [A m /W]
c
      if(denom.gt.0.0d0) effjcd=-ceff*anum/(anucc*denom)
      return
      end
c
c
c
      subroutine curr_int(fcur,resj,resp)
      implicit real*8(a-h,o-z)
      parameter(epsa=0.0d0,epsr=1.0d-2)
      parameter(xxcr=16.0d0)
      parameter (lw=5000,liw=lw/4)
      dimension w(lw),iw(liw)
      external fcur,fpp

      common/nhn/nhn
      common/ygyg/yg
      common/nplr/anpl,anpr
      common/amut/amu
      common/gg/uplp,uplm,ygn
      common/ierr/ierr
      common/iokh/iokhawa
      common/cst/cst,cst2

c     EC power and current densities

      anpl2=anpl*anpl
      dnl=1.0d0-anpl2
      ygn=nhn*yg
      ygn2=ygn*ygn

      resj=0.0d0
      resj1=0.0d0
      resj2=0.0d0
      resp=0.0d0
c
      rdu2=anpl2+ygn2-1.0d0
      uplp=0.0d0
      uplm=0.0d0
      upltp=0.0d0
      upltm=0.0d0
c
      if (rdu2.ge.0.0d0) then
        rdu=sqrt(rdu2)
        uplp=(anpl*ygn+rdu)/dnl
        uplm=(anpl*ygn-rdu)/dnl
c
        uu1=uplm
        uu2=uplp
        xx1=amu*(anpl*uu1+ygn-1.0d0)
        xx2=amu*(anpl*uu2+ygn-1.0d0)
c
        if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
        if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
        duu=abs(uu1-uu2)
c
        if(duu.gt.1.d-6) then
          call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier,
     .             liw,lw,last,iw,w)
          if (ier.gt.0) ierr=90
        end if
c
        rdu2t=cst2*anpl2+ygn2-1.0d0
c
        if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then
c
c       resonance curve does not cross the trapping region
c
          iokhawa=0
          if(duu.gt.1.d-4) then
            call dqags(fcur,uu1,uu2,epsa,epsr,
     .               resj,ej,neval,ier,liw,lw,last,iw,w)
            if (ier.gt.0) ierr=91
          end if
        else
c
c       resonance curve crosses the trapping region
c
          iokhawa=1
          rdut=sqrt(rdu2t)
          upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2)
          upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2)
c
          uu1=uplm
          uu2=upltm
          xx1=amu*(anpl*uu1+ygn-1.0d0)
          xx2=amu*(anpl*uu2+ygn-1.0d0)
          if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
            if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
            if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
            duu=abs(uu1-uu2)
            if(duu.gt.1.d-6) then
              call dqags(fcur,uu1,uu2,epsa,epsr,
     .        resj1,ej1,neval,ier,liw,lw,last,iw,w)
              if (ier.gt.0)   then
              if (abs(resj1).lt.1.0d-10) then
                resj1=0.0d0
              else
                ierr=92
              end if
              end if
            end if
          end if
c
          uu1=upltp
          uu2=uplp
          xx1=amu*(anpl*uu1+ygn-1.0d0)
          xx2=amu*(anpl*uu2+ygn-1.0d0)
          if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
            if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl
            if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl
            duu=abs(uu1-uu2)
            if(duu.gt.1.d-6) then
              call dqags(fcur,uu1,uu2,epsa,epsr,
     .        resj2,ej2,neval,ier,liw,lw,last,iw,w)
              if (ier.gt.0) then
                if(ier.ne.2) ierr=93
            end if
            end if
          end if
          resj=resj1+resj2
        end if
      end if

      return
      end
c
c     computation of integral for power density, integrand function fpp
c
c     ith=0 : polarization term = const
c     ith=1 : polarization term  Larmor radius expansion to lowest order
c     ith=2 : full polarization term (J Bessel)
c
      function fpp(upl)
      implicit real*8 (a-h,o-z)
      complex*16 ex,ey,ez,ui,emxy,epxy
      parameter(ui=(0.0d0,1.0d0))
      
      common/ygyg/yg
      common/nplr/anpl,anpr
      common/amut/amu
      common/gg/uplp,uplm,ygn
      common/epolar/ex,ey,ez
      common/nprw/anprre,anprim
      common/ithn/ithn
      common/nhn/nhn

      upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
      gam=anpl*upl+ygn
      ee=exp(-amu*(gam-1))

      thn2=1.0d0
      thn2u=upr2*thn2
      if(ithn.gt.0) then
        emxy=ex-ui*ey
        epxy=ex+ui*ey
        if(upr2.gt.0.0d0) then
          upr=sqrt(upr2)
          bb=anprre*upr/yg
          if(ithn.eq.1)  then
c     Larmor radius expansion polarization term at lowest order
            cth=1.0d0
            if(nhn.gt.1) cth=(0.5d0*bb)**(nhn-1)*nhn/fact(nhn)
            thn2=(0.5d0*cth*abs(emxy+ez*anprre*upl/ygn))**2
            thn2u=upr2*thn2
          else
c     Full polarization term
            nm=nhn-1
            np=nhn+1
            ajbnm=dbesjn(nm, bb)
            ajbnp=dbesjn(np, bb)
            ajbn=dbesjn(nhn, bb)
           thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0d0))**2
          end if
        end if
      end if
      
      fpp=ee*thn2u
      return
      end

c     computation of integral for current density
c     fjch integrand for Cohen model with trapping
c     fjch0 integrand for Cohen model without trapping    
c     fjncl integrand for momentum conserv. model K(u) from Maruschenko
c     gg=F(u)/u with F(u) as in Cohen paper

      function fjch(upl)
      implicit real*8 (a-h,o-z)

      common/gg/uplp,uplm,ygn
      common/nplr/anpl,anpr
      common/zz/Zeff
      common/cohen/rb,alams,pa,fp0s
      
      upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm)
      gam=anpl*upl+ygn
      u2=gam*gam-1.0d0
      u=sqrt(u2)
      z5=Zeff+5.0d0
      xi=1.0d0/z5**2
      xib=1.0d0-xi
      xibi=1.0d0/xib
      fu2b=1.0d0+xib*u2
      fu2=1.0d0+xi*u2
      gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2)
      gg=u*gu/z5
      dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5

      alam=sqrt(1.0d0-upr2/u2/rb)
      fp0=fconic(alam,pa,0)
      dfp0=-(pa*pa/2.0d0+0.125d0)
      if (alam.lt.1.0d0) then
        dfp0=-fconic(alam,pa,1)/sqrt(1.0d0-alam**2)
      end if
      fh=alam*(1.0d0-alams*fp0/(alam*fp0s))
      dfhl=1.0d0-alams*dfp0/fp0s

      eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam)

      if(upl.lt.0.0d0) eta=-eta
      fjch=eta*fpp(upl)
      return
      end



      function fjch0(upl)
      implicit real*8 (a-h,o-z)
      common/gg/uplp,uplm,ygn
      common/nplr/anpl,anpr
      common/zz/Zeff
      gam=anpl*upl+ygn
      u2=gam*gam-1.0d0
      u=sqrt(u2)
      z5=Zeff+5.0d0
      xi=1.0d0/z5**2
      xib=1.0d0-xi
      xibi=1.0d0/xib
      fu2b=1.0d0+xib*u2
      fu2=1.0d0+xi*u2
      gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2)
      gg=u*gu/z5
      dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
      eta=anpl*gg+gam*upl*dgg/u
      fjch0=eta*fpp(upl)
      return
      end



      function fjncl(upl)
      use  green_func_p
      implicit real*8 (a-h,o-z)    

      common/gg/uplp,uplm,ygn
      common/nplr/anpl,anpr
      common/fc/fc
      common/ncl/hb
      common/psival/psinv
      common/amut/amu
      common/tete/tekev
      common/zz/Zeff

      gam=anpl*upl+ygn
      u2=gam*gam-1.0d0
      u=sqrt(u2)
      upr2=u2-upl*upl
      bth=sqrt(2.0d0/amu)
      uth=u/bth
      call GenSpitzFunc(Tekev,Zeff,fc,uth,u,gam,fk,dfk)
      fk=fk*(4.0d0/amu**2)
      dfk=dfk*(2.0d0/amu)*bth
      
      alam=upr2/u2/hb
      psi=psinv
      call vlambda(alam,psi,fu,dfu)
      
      eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb     
      if(upl.lt.0) eta=-eta
      fjncl=eta*fpp(upl)
      return
      end

      subroutine vlambda(alam,psi,fv,dfv)
      implicit real*8 (a-h,o-z)
      parameter (nnintp=101,nlam=41)
      parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
      parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
     .          njest*nnintp+nlest+54)
      parameter(kwrk=nnintp+nlam+njest+nlest+3)
      parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)
      
      external fpbisp

      dimension xxs(1),yys(1),ffs(1)
      dimension ch01(lw01),ch((njest-4)*(nlest-4))
      dimension tjp(njest),tlm(nlest)
      dimension iwrk(kwrk),wrk(lwrk)

      common/coffvl/ch
      common/coffdvl/ch01
      common/coffvlt/tjp,tlm
      common/coffvln/njpt,nlmt

      njp=njpt
      nlm=nlmt
      xxs(1)=sqrt(psi)
      yys(1)=alam
      
      call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs,
     .            wrk(1),wrk(5),iwrk(1),iwrk(2))
      fv=ffs(1)

      iwp=1+(njp-4)*(nlm-5)
      iwl=iwp+4
      call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2,
     .     xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2))
      dfv=ffs(1)

      return
      end  


      subroutine projxyzt(iproj,nfile)
      implicit real*8 (a-h,o-z)
      parameter(jmx=31,kmx=36)
      dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c
      common/nray/nrayr,nrayth
      common/wrk/ywrk,ypwrk
      common/psinv11/psinv11
      common/istep/istep
      common/ss/st
c
      rtimn=1.d+30
      rtimx=-1.d-30
c
      jd=1
      if(iproj.eq.0) jd=nrayr-1
      do j=1,nrayr,jd
        kkk=nrayth
        if(j.eq.1) kkk=1
        do k=1,kkk
c
          dx=ywrk(1,j,k)-ywrk(1,1,1)
          dy=ywrk(2,j,k)-ywrk(2,1,1)
          dz=ywrk(3,j,k)-ywrk(3,1,1)
c
          dirx=ywrk(4,j,k)
          diry=ywrk(5,j,k)
          dirz=ywrk(6,j,k)
          dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
c
          if(j.eq.1.and.k.eq.1) then
            csth1=dirz/dir
            snth1=sqrt(1.0d0-csth1**2)
            csps1=1.0d0
            snps1=0.0d0
            if(snth1.gt.0.0d0) then
              csps1=diry/(dir*snth1)
              snps1=dirx/(dir*snth1)
            end if
          end if
          xti= dx*csps1-dy*snps1
          yti=(dx*snps1+dy*csps1)*csth1-dz*snth1
          zti=(dx*snps1+dy*csps1)*snth1+dz*csth1
          rti=sqrt(xti**2+yti**2)
c
          dirxt= (dirx*csps1-diry*snps1)/dir
          diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir
          dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir
c
          if(k.eq.1) then
            xti1=xti
            yti1=yti
            zti1=zti
            rti1=rti
          end if

          if(.not.(iproj.eq.0.and.j.eq.1))
     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
c
          if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
          if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
c
        end do
c
c        if(.not.(iproj.eq.0.and.j.eq.1))
c     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
        if(iproj.eq.1) write(nfile,*) ' '
      end do
c
      write(nfile,*) ' '
c
      write(12,99) istep,st,psinv11,rtimn,rtimx
      return
 99   format(i5,12(1x,e16.8e3))
111   format(3i5,12(1x,e16.8e3))
      end
c
c
c
      subroutine pec(pabs,currt)
      implicit real*8(a-h,o-z)
      parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000)
      parameter(rtbc=1.0d0)
      parameter(pi=3.14159265358979d0)
      dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx)
      dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
      dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx)
      dimension xxi(nmx),ypt(nmx),yamp(nmx)
      dimension rtab(nndmx),rhotv(nndmx)
      dimension rtab1(0:nndmx)
      dimension ajphiv(nndmx),dpdv(nndmx)
      dimension dvolt(nndmx),darea(nndmx)
      dimension ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx)
      dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx)
      dimension pins(nndmx),currins(nndmx),fi(nndmx)
      parameter(llmx=21)
      dimension isev(llmx)
c
      common/nray/nrayr,nrayth
      common/istep/istep
      common/dsds/dst
      common/ipec/ipec,nnd
      common/ieccd/ieccd
      common/index_rt/index_rt
c
      common/iiv/iiv
      common/psjki/psjki
      common/pcjki/ppabs,ccci
      common/dpjjki/pdjki,currj
c
      common/cent/btrcen,rcen
      common/angles/alpha0,beta0
      common/iieq/iequil
      common/parban/b0,rr0m,zr0m,rpam
      common/taumnx/taumn,taumx,pabstot,currtot
      common/jmxmn/rhot1,rhot2,aj1,aj2
c
      common/polcof/psipol,chipol
c
      stf=istep*dst
      nd=nnd

      if(pabstot.gt.0.0d0) then
      do ll=1,llmx
        isev(ll)=0
      end do

      intp=0
      voli0=0.0d0
      areai0=0.0d0
      rtab1(0)=0.0d0
      do it=1,nd
        drt=1.0d0/dble(nd-1)
        rt=dble(it-1)*drt
        if(it.lt.nd) then
          rt1=rt+drt/2.0d0
        else
          rt1=rt
        end if
        rtab(it)=rt
        rtab1(it)=rt1
        dpdv(it)=0.0d0
        ajphiv(it)=0.0d0
        if (ipec.eq.0) then
          psit=rt
          psit1=rt1
        else
          psit=rt**2
          psit1=rt1**2
        end if
        if(iequil.eq.2) then
          rhotv(it)=frhotor(psit)
          call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii,
     .                  rbavi,bmxi,bmni,fci,intp)
          dvolt(it)=abs(voli1-voli0)
          darea(it)=abs(areai1-areai0)
          voli0=voli1
          areai0=areai1
          call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli)
          ratjav(it)=ratjai
          ratjbv(it)=ratjbi
          ratjplv(it)=ratjpli
        else
          rhotv(it)=sqrt(psit)
          area1=pi*psit1*rpam**2
          voli1=2.0d0*pi*rr0m*area1
          dvolt(it)=abs(voli1-voli0)
          darea(it)=abs(areai1-areai0)
          voli0=voli1
          areai0=areai1
        end if
      end do

      kkk=1
      do j=1,nrayr
        if(j.gt.1) kkk=nrayth
        do k=1,kkk
          ise0=0
          ii=iiv(j,k)
          if (ii.lt.nmx.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1
          idecr=-1
          is=1
          do i=1,ii
            if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i))
            if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i)))
            if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then
              ypt(i)=ppabs(j,k,i)
              yamp(i)=ccci(j,k,i)
            else
              ypt(i)=0.0d0
              yamp(i)=0.0d0
            end if
            if(ise0.eq.0) then
              if(xxi(i).lt.rtbc) then
              ise0=i
              isev(is)=i-1
              is=is+1
              end if
            else
              if (idecr.eq.-1) then
                if(xxi(i).gt.xxi(i-1)) then
                  isev(is)=i-1
                  is=is+1
                  idecr=1
                 end if
              else
                if(xxi(i).gt.rtbc) exit
                if(xxi(i).lt.xxi(i-1)) then
                  isev(is)=i
                  is=is+1
                  idecr=-1
                end if
              end if
            end if
          end do
c
          isev(is)=i-1
          ppa1=0.0d0
          cci1=0.0d0
          do iis=1,is-1
            iis1=iis+1
            iise0=isev(iis)
            iise=isev(iis1)
            if (mod(iis,2).ne.0) then
              idecr=-1
              ind1=nd
              ind2=2
              iind=-1
            else
              idecr=1
              ind1=1
              ind2=nd
              iind=1
            end if
            do ind=ind1,ind2,iind
              iind=ind
              if (idecr.eq.-1) iind=ind-1
              rt1=rtab1(iind)
              call locatex(xxi,iise,iise0,iise,rt1,itb1)
              if(itb1.gt.iise0.and.itb1.lt.iise) then
                call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),
     .                      ypt(itb1+1),rt1,ppa2)
                call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),
     .                      yamp(itb1+1),rt1,cci2)
                dppa=ppa2-ppa1
                dpdv(ind)=dpdv(ind)+dppa
                didst=(cci2-cci1)
                ajphiv(ind)=ajphiv(ind)+didst
                ppa1=ppa2
                cci1=cci2
              end if
            end do
          end do
        end do
      end do
 
      h=1.0d0/dble(nd-1)
      rhotpav=0.0d0
      drhotpav=0.0d0
      rhotjav=0.0d0
      rhotjava=0.0d0
      rhot2java=0.0d0
     
      fi=dpdv/h
      call simpson (nd,h,fi,spds)
      fi=rhotv*dpdv/h
      call simpson (nd,h,fi,rhotpav)
      fi=rhotv*rhotv*dpdv/h
      call simpson (nd,h,fi,rhot2pav)
      rhotpav=rhotpav/spds
      rhot2pav=rhot2pav/spds

      if (ieccd.ne.0) then      
      fi=ajphiv/h
      call simpson (nd,h,fi,sccs)
      fi=rhotv*ajphiv/h
      call simpson (nd,h,fi,rhotjav)    
      fi=abs(ajphiv)/h
      call simpson (nd,h,fi,sccsa)
      fi=rhotv*abs(ajphiv)/h
      call simpson (nd,h,fi,rhotjava)
      fi=rhotv*rhotv*abs(ajphiv)/h
      call simpson (nd,h,fi,rhot2java)     
      rhotjav=rhotjav/sccs
      rhotjava=rhotjava/sccsa
      rhot2java=rhot2java/sccsa
      end if

c     factor sqrt(8)=2 sqrt(2) to match with full width
c     of gaussian profile
      drhot2pav=rhot2pav-rhotpav**2
      drhotpav=sqrt(8.d0*drhot2pav)
      drhot2java=rhot2java-rhotjava**2
      drhotjava=sqrt(8.d0*drhot2java)
      
      spds=0.0d0
      sccs=0.0d0
      do i=1,nd
        spds=spds+dpdv(i)
        sccs=sccs+ajphiv(i)
        pins(i)=spds
        currins(i)=sccs
        dpdv(i)=dpdv(i)/dvolt(i)
        ajphiv(i)=ajphiv(i)/darea(i)
      end do
      
      facpds=1.0d0
      facjs=1.0d0
      if(spds.gt.0.0d0) facpds=pabs/spds
      if(sccs.ne.0.0d0) facjs=currt/sccs

      do i=1,nd
        dpdv(i)=facpds*dpdv(i)
        ajphiv(i)=facjs*ajphiv(i)
        ajcdav(i)=ajphiv(i)*ratjav(i)
        ajcdbv(i)=ajphiv(i)*ratjbv(i)
        ajplv(i)=ajphiv(i)*ratjplv(i)
      end do

      call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
     .                 drhotp,drhopp)
      if(ieccd.ne.0) then
        call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
     .                 drhotjfi,drhopfi)
        xps=rhopfi
        call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
      else
        rhotjfi=0.0d0
        rhopfi=0.0d0
        ajmxfi=0.0d0
        drhotjfi=0.0d0
        drhopfi=0.0d0
        xps=rhopp
        ratjamx=1.0d0
        ratjbmx=1.0d0
        ratjplmx=1.0d0
      end if

       iif1=iiv(1,1)
       istmx=1
       do i=2,iif1
        if(psjki(1,1,i).ge.0.0d0) then
         if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i
        end if
      end do
      stmx=istmx*dst
      
      pins_02=0.0d0
      pins_05=0.0d0
      pins_085=0.0d0

      xrhot=0.2d0
      call locate(rhotv,nd,xrhot,i1)
      call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
     .  xrhot,pins_02)
      xrhot=0.5d0
      call locate(rhotv,nd,xrhot,i1)
      call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
     .  xrhot,pins_05)
      xrhot=0.85d0
      call locate(rhotv,nd,xrhot,i1)
      call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
     .  xrhot,pins_085)
     
      else
      ajmxfi=0.0d0
      dpdvmx=0.0d0
      rhotjfi=1.0d0
      rhotjav=1.0d0
      rhotjava=1.0d0
      rhotp=1.0d0
      rhotpav=1.0d0
      drhotjfi=0.0d0
      drhotjava=0.0d0
      drhotp=0.0d0
      drhotpav=0.0d0
      taumn=0
      taumx=0
      stmx=stf
      pins_02=0.0d0
      pins_05=0.0d0
      pins_085=0.0d0     
c  end of pabstot > 0   
      end if     

c     dPdV [MW/m^3], Jcd [MA/m^2]

      if(ieccd.eq.0) currt=0.0d0
      currtka=currt*1.0d3

      write(6,*)' '
      write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '//
     .'rhotj rhotjava rhotp rhotpav '//
     .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
      write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx,
     .  rhotjfi,rhotjava,rhotp,rhotpav,
     .  drhotjava,drhotpav,
     .  stmx,psipol,chipol,real(index_rt)
     
      write(7,99) currtka,pabstot,ajmxfi,dpdvmx,
     .  rhotjfi,rhotjava,rhotp,rhotpav,
     .  drhotjava,drhotpav,ratjbmx,
     .  stmx,psipol,chipol,real(index_rt)

      do i=1,nd
        if (ipec.eq.0) then
          psin=rtab(i)
          rhop=sqrt(rtab(i))
        else
          psin=rtab(i)**2
          rhop=rtab(i)
        end if
        pinsr=0.0d0
        if(pabstot.gt.0) pinsr=pins(i)/pabstot
        write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
     .   currins(i),pins(i),pinsr,real(index_rt)
      end do

      return
 49   format(i5,20(1x,e12.5))
 99   format(30(1x,e12.5))
      end
c
c
c
      subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop)
      implicit real*8(a-h,o-z)
      parameter(emn1=0.367879441171442d0)
      dimension xx(nd),yy(nd)
      common/jmxmn/rhotp,rhotm,ypkp,ypkm
      common/ipec/ipec,nnd
      common/iieq/iequil
c
      call vmaxmini(yy,nd,ymn,ymx,imn,imx)
      ypk=0.0d0
      xmx=xx(imx)
      xmn=xx(imn)
      if (abs(ymx).gt.abs(ymn)) then
        ipk=imx
        ypkp=ymx
        xpkp=xmx
        if(abs(ymn/ymx).lt.1d-2) ymn=0.0d0
        ypkm=ymn
        xpkm=xmn
      else
        ipk=imn
        ypkp=ymn
        xpkp=xmn
        if(abs(ymx/ymn).lt.1d-2) ymx=0.0d0
        ypkm=ymx
        xpkm=xmx
      end if
      if(xpkp.gt.0.0d0) then
        xpk=xpkp
        ypk=ypkp
        yye=ypk*emn1
        call locatex(yy,nd,1,ipk,yye,ie1)
        if(ie1.gt.0.and.ie1.lt.nd) then
          call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1)
        else
          rte1=0.0d0
        end if
        call locatex(yy,nd,ipk,nd,yye,ie2)
        call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2)
      else
        ipk=2
        xpk=xx(2)
        ypk=yy(2)
        rte1=0.0d0
        yye=ypk*emn1
        call locate(yy,nd,yye,ie2)
        call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
      end if
c
      ipk=1
      if(ymx.ne.0.and.ymn.ne.0) ipk=2
c
      drhop=0.0d0
      drhot=0.0d0
      psie1=0.0d0
      psie2=1.0d0
      rhopmx=1.0d0
      rhopmn=0.0d0
      if (ie1.gt.0.or.ie2.gt.0) then
        if(ipec.eq.0) then
          rhopmx=sqrt(xpkp)
          rhopmn=sqrt(xpkm)
          psie2=rte2
          psie1=rte1
          drhop=sqrt(rte2)-sqrt(rte1)
        else
          rhopmx=xpkp
          rhopmn=xpkm
          drhop=rte2-rte1
          psie2=rte2**2
          psie1=rte1**2
        end if
      end if
c
      if(iequil.eq.2) then
        rhotmx=frhotor(rhopmx**2)
        rhotmn=frhotor(rhopmn**2)
        rhote2=frhotor(psie2)
        rhote1=frhotor(psie1)
        drhot=rhote2-rhote1
      else
        rhotmx=rhopmx
        rhotmn=rhopmn
        drhot=drhop
      end if
c
      if(ipk.eq.2) then
        drhop=-drhop
        drhot=-drhot
      end if

      ypk=ypkp
      rhotp=rhotmx
      rhotm=rhotmn
c
      return
      end
c
      subroutine polarcold(exf,eyif,ezf,elf,etf)
      implicit real*8(a-h,o-z)
      common/nplr/anpl,anpr
      common/xgxg/xg
      common/ygyg/yg
      common/mode/sox
c
c     dcold dispersion
c     dielectric tensor (transposed)
c
c      exf=0.0d0
c      eyif=0.0d0
c      ezf=0.0d0
c      if(xg.le.0) return
c
      anpl2=anpl*anpl
      anpr2=anpr*anpr
      an2=anpl2+anpr2
      yg2=yg**2
      dy2=1.0d0-yg2
      aa=1.0d0-xg-yg2
      e3=1.0d0-xg
c
      if(xg.gt.0.0d0) then
        if (anpl.ne.0.0d0) then
          qq=xg*yg/(an2*dy2-aa)
          p=(anpr2-e3)/(anpl*anpr)
          ezf=1.0d0/sqrt(1.0d0+p*p*(1.0d0+qq*qq))
          exf=p*ezf
          eyif=qq*exf
        else
          if(sox.lt.0.d0) then
            ezf=1
            exf=0
            eyif=0
          else
            ezf=0
            qq=-aa/(xg*yg)
            exf=1.0d0/sqrt(1.0d0+qq*qq)
            eyif=qq*exf
          end if
        end if
        elf=(anpl*ezf+anpr*exf)/sqrt(an2)
        etf=sqrt(1.0d0-elf*elf)
      else
        if(sox.lt.0.0d0) then
          ezf=1
          exf=0
          eyif=0
        else
          ezf=0
          exf=0.0d0
          eyif=1.0d0
        end if
        elf=0
        etf=1
      end if
c
      return
      end

      subroutine pol_limit(qq,uu,vv)
      implicit none
      integer*4 ipolc
      real*8 bv(3),anv(3)
      real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm
      real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2
      real*8 deno,denx,anpl2,dnl,del0
      real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv
      real*8 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm
      real*8 pi,beta0,alpha0,gam
      real*8 sox,psipol,chipol
      complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt
      parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
c
      common/anv/anv
      common/nplr/anpl,anpr
      common/ygyg/yg
      common/bb/bv
      common/angles/alpha0,beta0
      common/mode/sox
      common/polcof/psipol,chipol
      common/ipol/ipolc
      common/evt/ext,eyt
c
      anx=anv(1)
      any=anv(2)
      anz=anv(3)
      anpl2=anpl*anpl
      dnl=1.0d0-anpl2
      del0=sqrt(dnl**2+4.0d0*anpl2/yg**2)
      ffo=0.5d0*yg*(dnl+del0)
      ffx=0.5d0*yg*(dnl-del0)
      an2=anx*anx+any*any+anz*anz
      an=sqrt(an2)
      anxy=sqrt(anx*anx+any*any)
      sngam=(anz*anpl-an2*bv(3))/(an*anxy*anpr)
      csgam=-(any*bv(1)-anx*bv(2))/anxy/anpr
      csg2=csgam**2
      sng2=sngam**2
      ffo2=ffo*ffo
      ffx2=ffx*ffx
      deno=ffo2+anpl2
      denx=ffx2+anpl2
      xe2om=(ffo2*csg2+anpl2*sng2)/deno
      ye2om=(ffo2*sng2+anpl2*csg2)/deno
      xe2xm=(ffx2*csg2+anpl2*sng2)/denx
      ye2xm=(ffx2*sng2+anpl2*csg2)/denx
c
      exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno)
      eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno)
      qqom=abs(exom)**2-abs(eyom)**2
      uuom=2.0d0*dble(exom*dconjg(eyom))
      vvom=2.0d0*dimag(exom*dconjg(eyom))
      llmom=sqrt(qqom**2+uuom**2)
      aaom=sqrt((1+llmom)/2.0d0)
      bbom=sqrt((1-llmom)/2.0d0)
      ellom=bbom/aaom
      psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi
      chiom=0.5d0*asin(vvom)*180.d0/pi
c
      exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx)
      eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx)
      qqxm=abs(exxm)**2-abs(eyxm)**2
      uuxm=2.0d0*dble(exxm*dconjg(eyxm))
      vvxm=2.0d0*dimag(exxm*dconjg(eyxm))
      llmxm=sqrt(qqxm**2+uuxm**2)
      aaxm=sqrt((1+llmxm)/2.0d0)
      bbxm=sqrt((1-llmxm)/2.0d0)
      ellxm=bbxm/aaxm
      psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi
      chixm=0.5d0*asin(vvxm)*180.d0/pi
c
      if (sox.lt.0.0d0) then
        psipol=psiom
        chipol=chiom
        ext=exom
        eyt=eyom
        qq=qqom
        vv=vvom
        uu=uuom
      else
        psipol=psixm
        chipol=chixm
        ext=exxm
        eyt=eyxm
        qq=qqxm
        vv=vvxm
        uu=uuxm
      endif
      gam=atan(sngam/csgam)*180.d0/pi

      return
 111  format(20(1x,e12.5))      
      end



      subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl)
      implicit none
      integer*4 ivac,irfl
      real*8 anv(3),xv(3),xvrfl(3)
      real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3)
      real*8 uutr,vvtr,qqtr,qq,uu,vv
      real*8 pi
      real*8 psipol,chipol,psitr,chitr
      complex*16 ui,extr,eytr,eztr,ext,eyt
      complex*16 evin(3),evrfl(3) 
      parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)

      common/xv/xv
      common/anv/anv
      common/polcof/psipol,chipol
      common/evt/ext,eyt
      common/wrefl/walln
      
      anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
      
c     computation of reflection coordinates and normal to the wall
      irfl=1
      ivac=1
      call vacuum_rt(xv,xvout,ivac)
      
      if(ivac.lt.0) then
        irfl=0
        xvrfl=xvout
        xv=xvout
        anvrfl=anv
        return
      end if  
      
c     rotation matrix from local to lab frame
      vv1(1)=anv(2)
      vv1(2)=-anv(1)
      vv1(3)=0.0d0
      vv2(1)=anv(1)*anv(3)
      vv2(2)=anv(2)*anv(3)
      vv2(3)=-anv(1)*anv(1)-anv(2)*anv(2)
      vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
      vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
      vv3=anv
      
      evin=ext*vv1+eyt*vv2
c     wave vector and electric field after reflection in lab frame
      anvrfl=anv-2.0d0*
     .       (anv(1)*walln(1)+anv(2)*walln(2)+anv(3)*walln(3))*walln
      evrfl=-evin+2.0d0*
     .     (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln

      vv1(1)=anvrfl(2)
      vv1(2)=-anvrfl(1)
      vv1(3)=0.0d0
      vv2(1)=anvrfl(1)*anvrfl(3)
      vv2(2)=anvrfl(2)*anvrfl(3)
      vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
      vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
      vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
      vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)

      extr=dot_product(vv1,evrfl)
      eytr=dot_product(vv2,evrfl)
      eztr=dot_product(vv3,evrfl)

      qqtr=abs(extr)**2-abs(eytr)**2
      uutr=2.0d0*dble(extr*dconjg(eytr))
      vvtr=2.0d0*dimag(extr*dconjg(eytr))
      psitr=0.5d0*atan2(uutr,qqtr)*180.d0/pi
      chitr=0.5d0*asin(vvtr)*180.d0/pi

      ivac=2
      anv=anvrfl
      xvrfl=xvout
      xv=xvout
      
      call vacuum_rt(xv,xvout,ivac)
      xv=xvout
      
      return
 111  format(20(1x,e12.5))      
      end
      
      subroutine vacuum_rt(xvstart,xvend,ivac)
      use reflections, only : inters_linewall,inside
      implicit none
      integer*4 ivac
      integer nbb,nlim,i,imax
      parameter(nbb=1000)
      real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax
      real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
      real*8 xv0(3)
      real*8 rlim(nbb),zlim(nbb)
      
      common/wrefl/walln
      common/limiter/rlim,zlim,nlim
      common/anv/anv
      common/dsds/dst      
      common/psival/psinv
      common/densbnd/psdbnd
      common/dstvac/dstvac
c     ivac=1 first interface plasma-vacuum
c     ivac=2 second interface vacuum-plasma after wall reflection
c     ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
      
!     the real argument associated to xvstart is in a common block
!     used by fwork!!!
      xv0=xvstart
      call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
     .                     nlim,smax,walln)
      smax=smax*1.d2
      rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2)
      zzm=1.d-2*xv0(3)
      if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then
        ! first wall interface is outside-inside
        if (dot_product(walln,walln)<tiny(walln)) then
          ! wall never hit
          dstvac=0.0d0
          xvend=xv0
          ivac=-1
          return
        end if
        ! search second wall interface (inside-outside)
        st=smax
        xvend=xv0+st*anv
        call inters_linewall(xvend/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
     .                       nlim,smax,walln)
        smax=smax*1.d2+st
      end if
      i=0
      do
        st=i*dst
        xvend=xv0+st*anv
        y(1)=xvend(1)
        y(2)=xvend(2)
        y(3)=xvend(3)
        y(4)=anv(1)
        y(5)=anv(2)
        y(6)=anv(3)
        call fwork(y,dery)
        if (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit
        i=i+1
      end do
      
      if (st.lt.smax) then
        ivac=-1
        dstvac=st
      else
        dstvac=smax
        xvend=xv0+smax*anv
      end if

!     the real argument associated to xvstart is in a common block
!     used by fwork!!!
      xvstart=xv0
      
      return
 111  format(20(1x,e12.5))      
      end