!     program gray
!     use gray_params, only : ipass,igrad
!     implicit none
!! local variables
!     real(wp_) :: p0mw1
!! common/external functions/variables
!     integer :: ierr,index_rt
!     real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot,
!!
!     common/ierr/ierr
!     common/mode/sox
!     common/p0/p0mw
!     common/powrfl/powrfl
!     common/index_rt/index_rt
!     common/taumnx/taumn,taumx,pabstot,currtot
!!
!     if (ipass.gt.1) then
!!     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.0_wp_-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
!!
!     end program gray
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!         ins_pl=inside_plasma(rrm,zzm)
!         if (mod(iop(j,k),2).eq.0 .and. ins_pl) then
!           iop(j,k)=iop(j,k)+1
!           call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
!
!           if (ipass.gt.1 .and. index_rt.eq.1 .and. 
!    .          iowmax.gt.1 .and. istore(j,k).eq.0) then
!             istore(j,k)=istore(j,k)+1
!             yyrfl(j,k,1:3)=xv
!             yyrfl(j,k,4:6)=anv
!             ihcd(j,k)=0
!           end if
!         else if (mod(iop(j,k),2).eq.1.and.
!    .             .not.ins_pl) then
!           iop(j,k)=iop(j,k)+1
!           call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
!         end if
!         
!         if (ipass.gt.1) then
!           if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then
!             iow(j,k)=1
!           else if (iow(j,k).eq.1 .and. 
!    .               .not.inside(rlim,zlim,nlim,rrm,zzm)) then
!             iow(j,k)=2
!             if (ins_pl) then
!               iop(j,k)=iop(j,k)+1
!               call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
!             end if
!             call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)),
!    .                eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl)
!             istore(j,k)=istore(j,k)+1
!             yyrfl(j,k,1:3)=xvrfl
!             yyrfl(j,k,4:6)=anvrfl
!             tau1v(j,k)=tauv(j,k,iiv(j,k))
!             ext(j,k,iop(j,k))=extr
!             eyt(j,k,iop(j,k))=eytr
!             if (j.lt.jclosest) then
!               jclosest=j
!               anwcl=anw
!               xwcl=xvrfl
!             end if
!             xv=xvrfl
!             anv=anvrfl
!             rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2)
!             zzm=1.0e-2_wp_*xv(3)
!             ywrk(1:3,j,k)=xv
!             ywrk(4:6,j,k)=anv
!             igrad=0
!             call gwork(sox,xgcn,bres,j,k)
!             if (ins_pl) then
!               iop(j,k)=iop(j,k)+1
!               call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
!               if (index_rt.eq.1) ihcd(j,k)=0
!             end if
!           end if
!         end if
!
!         if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv
!         if(iop(j,k).lt.iopmin) iopmin=iop(j,k)
!         if(iow(j,k).lt.iowmin) iowmin=iow(j,k)
!         if(iow(j,k).gt.iowmax) iowmax=iow(j,k)
!
!         xvjk(:,j,k)=xv
!         anvjk(:,j,k)=anv
!
!       end do
!     end do
!     if(jclosest.le.nrayr) then
!       aknmin=1.0_wp_
!       do j=1,nrayr
!         kkk=nrayth
!         if(j.eq.1) kkk=1
!         do k=1,kkk
!           print*,i,j,k
!           print*,anwcl,xwcl,anvjk(1:2,j,k)
!           anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2))
!    .             /sqrt(xwcl(1)**2+xwcl(2)**2)
!           anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k))
!    .             /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2)
!           akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k)
!           if(akdotn.lt.aknmin) aknmin=akdotn
!         end do
!       end do
!     else
!       aknmin=-1.0_wp_
!     end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!     single pass is stopped when all the rays have crossed the plasma
!!         or complete absorption has occurred
!!     same for successive passes of multi-pass simulations (here exit
!!         from vessel is detected too
!!     first pass in multi-pass simulation is stopped when at least one
!!         ray has reflected and all rays are directed away from
!!         reflection point, or when no reflection has occurred and
!!         central ray re-enters the plasma
!     
!     if((ipass.eq.1 .and. ((iopmin.gt.1) .or.
!    . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr)))
!    . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or.
!    . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr)))) then
!       istop=1
!     else if(ipass.gt.1 .and. index_rt.eq.1 .and.
!    .   ((iowmin.gt.1 .and. aknmin.gt.0) .or.
!    .    (iowmax.le.1 .and. iop(1,1).gt.2))) then
!!     flag second pass mode coupling as unset
!       powrfl=-1.0_wp_
!       qqout=0.0_wp_
!       uuout=0.0_wp_
!       vvout=0.0_wp_
!       do j=1,nrayr
!         kkk=nrayth
!         if(j.eq.1) kkk=1
!         do k=1,kkk
!!     store missing initial conditions for the second pass
!           if (istore(j,k).eq.0) then
!             istore(j,k)=istore(j,k)+1
!             yyrfl(j,k,1:3)=xvjk(:,j,k)
!             yyrfl(j,k,4:6)=anvjk(:,j,k)
!             tau1v(j,k)=tauv(j,k,iiv(j,k))
!           end if
!!     determine mode coupling at the plasma boundary
!           if (powrfl.lt.0.0_wp_) then
!             call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac)
!!     look for first ray hitting the plasma, starting from the central
!!     and evaluate polarization
!             if (ivac.eq.1) then
!               y(1:3)=xvjk(:,j,k)
!               y(4:6)=anvjk(:,j,k)
!               call fwork(sox,xgcn,bres,y,dery)
!               call pol_limit(sox,exin2,eyin2)
!               call stokes(exin2,eyin2,qqin2,uuin2,vvin2)
!      powloop: do j1=1,nrayr
!                 kkkk=nrayth
!                 if(j1.eq.1) kkkk=1
!                 do k1=1,kkkk
!!     look for first ray which completed the first pass in the plasma
!                   if (iop(j1,k1).gt.1) then
!!     if found, use its polarization state to compute mode coupling
!                     call stokes(ext(j1,k1,2),eyt(j1,k1,2),
!    .                            qqout,uuout,vvout)
!                     exit powloop
!                   end if
!                 end do
!               end do powloop
!!     if no ray completed a first pass in the plasma, use central ray
!!     initial polarization (possibly reflected)
!               if (qqout.le.0.0_wp_) then
!                 call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout)
!               end if
!               powrfl=0.5_wp_*(1.0_wp_+vvout*vvin2+
!    .                              uuout*uuin2+qqout*qqin2)
!             end if
!           end if
!         end do
!       end do
!       strfl11=i*dst
!       write(6,*) '  '
!       write(6,*) 'Reflected power fraction =',powrfl
!       write(66,*) psipol,chipol,powrfl
!       istop=1
!     end if
!
!     return
!     end
!
!
!
!     subroutine ic_rt(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres,
!    . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0)
!! ray tracing initial conditions  igrad=0
!!
!     use const_and_precisions, only : wp_,izero,zero,one,pi,
!    . cvdr=>degree,ui=>im
!     use gray_params, only : ipol
!     use beamdata, only : nrayr,nrayth,rwmax,ywrk0=>ywrk,ypwrk0=>ypwrk,
!    . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt
!     implicit none
!! arguments
!     real(wp_), intent(in) :: x00,y00,z00,anx0c,any0c,anz0c
!     real(wp_), intent(in) :: ak0,xgcn,bres
!     real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir
!     real(wp_), intent(in) :: sox,psipol0,chipol0
!! local constants
!     integer, parameter :: ndim=6,ndimm=3
!! local variables
!     integer :: j,k,iv,jv,iproj,nfilp
!     real(wp_) :: csth,snth,csps,snps,phiwrad,csphiw,snphiw,dr,da,u,
!    . alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,
!    . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,vgradi,r0,
!    . x0m,y0m,r0m,z0m,ancsi,aneta,ppcsi,ppeta,deltapol,qq,uu,vv
!     real(wp_), dimension(ndim) :: ytmp,yptmp
!! common/external functions/variables
!     real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens,
!    . tekev,anpl,anpr,brr,bphi,bzz,ajphi,psipol,chipol,psinv11
!
!!
!     common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
!     common/nplr/anpl,anpr
!     common/psival/psinv
!     common/parpl/brr,bphi,bzz,ajphi
!     common/dens/dens,ddens
!     common/tete/tekev
!     common/polcof/psipol,chipol
!     common/psinv11/psinv11
!!
!     csth=anz0c
!     snth=sqrt(1.0_wp_-csth**2)
!     csps=1.0_wp_
!     snps=0.0_wp_
!     if(snth.gt.0.0_wp_) then
!       csps=any0c/snth
!       snps=anx0c/snth
!     end if
!!
!     phiwrad=phiw*cvdr
!     csphiw=cos(phiwrad)
!     snphiw=sin(phiwrad)
!!
!     dr=1.0_wp_
!     if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
!     da=2.0_wp_*pi/dble(nrayth)
!     z0t=0.0_wp_
!!
!     do j=1,nrayr
!       u=dble(j-1)
!       dffiu(j)=0.0_wp_
!       ddffiu(j)=0.0_wp_
!       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
!!
!!          csiw=u*dcsiw
!!          etaw=u*detaw
!!          csir=csiw
!!          etar=etaw
!!
!         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
!!
!         ppcsi=u*dr*cos(alfak)*rcicsi
!         ppeta=u*dr*sin(alfak)*rcieta
!!
!         anzt=1.0_wp_/sqrt(1.0_wp_+ppcsi**2+ppeta**2)
!         ancsi=ppcsi*anzt
!         aneta=ppeta*anzt
!!
!         anxt=ancsi*csphiw-aneta*snphiw
!         anyt=ancsi*snphiw+aneta*csphiw
!!
!         anx= anxt*csps+snps*(anyt*csth+anzt*snth)
!         any=-anxt*snps+csps*(anyt*csth+anzt*snth)
!         anz= anzt*csth-anyt*snth
!!
!         an20=1.0_wp_
!         an0=sqrt(an20)
!         anx0=anx
!         any0=any
!         anz0=anz
!!
!         xc0(1,j,k)=x0
!         xc0(2,j,k)=y0
!         xc0(3,j,k)=z0
!!
!         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
!!
!         ypwrk0(1,j,k) = anx0/an0
!         ypwrk0(2,j,k) = any0/an0
!         ypwrk0(3,j,k) = anz0/an0
!         ypwrk0(4,j,k) = 0.0_wp_
!         ypwrk0(5,j,k) = 0.0_wp_
!         ypwrk0(6,j,k) = 0.0_wp_
!!
!         ytmp=ywrk0(:,j,k)
!         yptmp=ypwrk0(:,j,k)
!         call fwork(sox,xgcn,bres,ytmp,yptmp)
!
!         if(ipol.eq.0) then
!           call pol_limit(sox,ext(j,k,0),eyt(j,k,0))
!           qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2
!           uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0)))
!           vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0)))
!           call polellipse(qq,uu,vv,psipol0,chipol0)
!         else
!           qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr)
!           uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr)
!           vv=sin(2.0_wp_*chipol0*cvdr)
!           if(qq**2.lt.1.0_wp_) then
!!             deltapol=phix-phiy, phix =0            
!             deltapol=atan2(vv,uu)
!             ext(j,k,0)= sqrt((1.0_wp_+qq)/2)
!             eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol)
!           else
!           if(qq.gt.0.0_wp_) then
!             ext(j,k,0)= 1.0_wp_
!             eyt(j,k,0)= 0.0_wp_
!             else
!             eyt(j,k,0)= 1.0_wp_
!             ext(j,k,0)= 0.0_wp_
!             end if
!           end if
!         endif
!         psipol=psipol0
!         chipol=chipol0
!!
!         do iv=1,3
!           gri(iv,j,k)=0.0_wp_
!           dgrad2v(iv,j,k)=0.0_wp_
!           du10(iv,j,k)=0.0_wp_
!           do jv=1,3
!             ggri(iv,jv,j,k)=0.0_wp_
!           end do
!         end do
!         grad2(j,k)=0.0_wp_
!!
!         dd=anx0**2+any0**2+anz0**2-an20
!         vgradi=0.0_wp_
!         ddi=2.0_wp_*vgradi
!!
!         r0=sqrt(x0**2+y0**2)
!         x0m=x0/1.0e2_wp_
!         y0m=y0/1.0e2_wp_
!         r0m=r0/1.0e2_wp_
!         z0m=z0/1.0e2_wp_
!         if(j.eq.nrayr) then
!            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
!    .          psinv,zero,anpl,zero,one
!         end if
!          if(j.eq.1.and.k.eq.1) then
!           psinv11=psinv
!           write(17,99) zero,zero,zero,zero
!           write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi,
!    .      psinv,one,dens,tekev,brr,bphi,bzz,
!    .      ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero,
!    .      zero,zero,zero,zero,zero,zero,zero,one
!         end if  
!       end do
!     end do
!
!     call pweigth
!!
!     if(nrayr.gt.1) then
!       iproj=0
!       nfilp=8
!       call projxyzt(iproj,nfilp)
!     end if
!!
!     return
!99   format(24(1x,e16.8e3))
!111   format(3i5,20(1x,e16.8e3))
!     end



  subroutine prfile
    implicit none
    write(4,*)' #sst R z phi psi rhot ne Te Btot '// &
        'Nperp Npl ki alpha tau Pt dIds nh iohkw index_rt ddr'
    write(8,*) ' #istep j k xt yt zt rt'
    write(9,*) ' #istep j k xt yt zt rt'
    write(17,*) ' #sst Dr_Nr1 Di_Nr1'
    write(33,*) ' #i jk sst x y R z psi tauv Npl alpha index_rt'
    write(12,*) ' #sst w1 w2'
    write(7,*)'#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav '// &
        'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt '// &
        'Jphimx dPdVmx drhotj drhotp'
    write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins'
!    write(66,*) "# psipol0 chipol0 powrfl"
  end subroutine prfile



  subroutine print_prof
    use const_and_precisions, only : wp_
    use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi
    use coreprofiles, only : density, temp
    implicit none
! local constants
    real(wp_), parameter :: eps=1.e-4_wp_
! local variables
    integer :: i
    real(wp_) :: psin,rhop,rhot,ajphi,te,qq
    real(wp_) :: dens,ddens

    write(55,*) ' #psi rhot ne Te q Jphi'
    do i=1,nq
      psin=psinr(i)
      rhop=sqrt(psin)

      call density(psin,dens,ddens)
      te=temp(psin)
      qq=fq(psin)
      rhot=frhotor(rhop)
      call tor_curr_psi(max(eps,psin),ajphi)
      write(55,"(12(1x,e12.5))") psin,rhot,dens,te,qq,ajphi*1.e-6_wp_
    end do
  end subroutine print_prof

  subroutine print_prof_an
    use const_and_precisions, only : wp_
    use coreprofiles, only : density, temp
    use equilibrium, only : frhotor
    implicit none
! local constants
    integer, parameter :: nst=51
! local variables
    integer :: i
    real(wp_) :: psin,rhop,rhot,te
    real(wp_) :: dens,ddens

    write(55,*) ' #psi rhot ne Te'
    do i=1,nst
      psin=dble(i-1)/dble(nst-1)
      rhop=sqrt(psin)
      rhot=frhotor(rhop)
      call density(psin,dens,ddens)
      te=temp(psin)
      write(55,"(12(1x,e12.5))") psin,rhot,dens,te
    end do
  end subroutine print_prof_an



  subroutine surfq(psinq,qpsi,nq,qval)
    use const_and_precisions, only : wp_
    use equilibrium, only : rmaxis,zmaxis,zbinf,zbsup,frhotor
    use magsurf_data, only : npoints,contours_psi
    use utils, only : locate, intlin
    implicit none
! arguments
    integer, intent(in) :: nq
    real(wp_), dimension(nq), intent(in) :: psinq,qpsi
    real(wp_) :: qval
! local variables
    integer :: ncnt,i1,ipr
    real(wp_) :: rup,zup,rlw,zlw,rhot,psival
    real(wp_), dimension(npoints) :: rcn,zcn

    ncnt=(npoints-1)/2

! locate  psi surface for q=qval
    call locate(abs(qpsi),nq,qval,i1)
    if (i1>0.and.i1<nq) then
      call intlin(abs(qpsi(i1)),psinq(i1),abs(qpsi(i1+1)),psinq(i1+1),qval,psival)
      rup=rmaxis
      rlw=rmaxis
      zup=(zbsup+zmaxis)/2.0_wp_
      zlw=(zmaxis+zbinf)/2.0_wp_
      ipr=1
      call contours_psi(psival,rup,zup,rlw,zlw,rcn,zcn,ipr)
      rhot=frhotor(sqrt(psival))
      print'(4(a,f8.5))','q = ',qval,         '  psi = ',psival, &
                    '  rhop = ',sqrt(psival),'  rhot = ',rhot
    end if
  end



  subroutine bfield_res(rv,zv,nr,nz,bres)
    use const_and_precisions, only : wp_
    use equilibrium, only : bfield
    implicit none
! arguments
    integer, intent(in) :: nr, nz
    real(wp_), intent(in) :: rv(nr), zv(nz), bres
! local constants
    integer, parameter :: icmx=2002
! local variables
    integer :: j,k,n,nconts,inc,nctot
    integer, dimension(10) :: ncpts
    real(wp_) :: btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
    real(wp_), dimension(icmx) :: rrcb,zzcb
    real(wp_), dimension(nr,nz) :: btotal

!   Btotal on psi grid
    btmx=-1.0e30_wp_
    btmn=1.0e30_wp_
    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)
      enddo
    enddo

!   compute Btot=Bres/n with n=1,5
    write(70,*)'#i Btot  R  z'
    do n=1,5
      bbb=bres/dble(n)
      if (bbb.ge.btmn.and.bbb.le.btmx) then
        nconts=size(ncpts)
        nctot=size(rrcb)
        call cniteq(rv,zv,btotal,nr,nz,bbb,nconts,ncpts,nctot,rrcb,zzcb)
        do inc=1,nctot
          write(70,'(i6,12(1x,e12.5))') inc,bbb,rrcb(inc),zzcb(inc)
        end do
      end if
      write(70,*)
    end do
  end subroutine bfield_res

  subroutine bres_anal(bres)
    use const_and_precisions, only : wp_,pi
    use equilibrium, only : aminor,rmaxis,zmaxis
    implicit none
! arguments
    real(wp_) :: bres
! local variables
    integer :: i
    integer, parameter :: ngrid=51
    real(wp_) :: dxgrid
    real(wp_), dimension(ngrid) :: rv,zv
      
    dxgrid=2.0_wp_*aminor/dble(ngrid-1)
    do i=1,ngrid
      rv(i) = rmaxis - aminor + dxgrid*(i-1)
      zv(i) = zmaxis - aminor + dxgrid*(i-1)
    end do  
    call bfield_res(rv,zv,ngrid,ngrid,bres)
  end  subroutine bres_anal



  subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
    use const_and_precisions, only : wp_
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
!                   (based on an older code)
    use const_and_precisions, only : wp_
    implicit none
! arguments
    integer, intent(in) :: nr,nz
    real(wp_), dimension(nr), intent(in) :: rqgrid
    real(wp_), dimension(nz), intent(in) :: zqgrid
    real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid
    real(wp_), intent(in) :: h
    integer, intent(inout) :: ncon, icount
    integer, dimension(ncon), intent(out) :: npts
    real(wp_), dimension(icount), intent(out) :: rcon,zcon
! local variables
    integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
    integer :: jabs,jnb,kx,ikx,itm,inext,in
    integer, dimension(3,2) :: ja
    integer, dimension(icount/2-1) :: lx
    real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
    real(wp_), dimension(nr*nz) :: a
    logical :: flag1

    px = 0.5_wp_

    a = reshape(matr2dgrid,(/nr*nz/))

    rcon = 0.0_wp_
    zcon = 0.0_wp_

    nrqmax = nr
    drgrd = rqgrid(2) - rqgrid(1)
    dzgrd = zqgrid(2) - zqgrid(1)

    ncon = 0

    npts = 0

    iclast = 0
    icount = 0
    mpl = 0
    ix = 0
    mxr = nrqmax * (nz - 1)
    n1 = nr - 1

    do jx=2,n1
      do jm=jx,mxr,nrqmax
        j = jm + nrqmax
        ah=a(j)-h
        if (ah <= 0.0_wp_ .and. a(jm) >  h .or. &
            ah >  0.0_wp_ .and. a(jm) <= h) then
          ix=ix+1
          lx(ix)=-j
        end if
        if (ah <= 0.0_wp_ .and. a(j-1) >  h .or. &
            ah >  0.0_wp_ .and. a(j-1) <= h) then
          ix=ix+1
          lx(ix)=j
        end if
      end do
    end do

    do jm=nr,mxr,nrqmax
      j = jm + nrqmax
      ah=a(j)-h
      if (ah <= 0.0_wp_ .and. a(j-1) >  h .or. &
          ah >  0.0_wp_ .and. a(j-1) <= h) then
        ix=ix+1
        lx(ix)=j
      end if
      if (ah <= 0.0_wp_ .and. a(jm) >  h .or. &
          ah >  0.0_wp_ .and. a(jm) <= h) then
        ix=ix+1
        lx(ix)=-j
      end if
    end do

    do jm=1,mxr,nrqmax
      j = jm + nrqmax
      if (a(j) <= h .and. a(jm) >  h .or. &
          a(j) >  h .and. a(jm) <= h) then
        ix=ix+1
        lx(ix) =-j
      end if
    end do

    do j=2,nr
      if (a(j) <= h .and. a(j-1) >  h .or. &
          a(j) >  h .and. a(j-1) <= h) then
        ix=ix+1
        lx(ix)=j
      end if
    end do

    if(ix<=0) return
  
bb: do
      in=ix
      jx=lx(in)
      jfor=0
      lda=1
      ldb=2

      do
        if(jx<0) then
          jabs=-jx
          jnb = jabs - nrqmax
        else
          jabs=jx
          jnb=jabs-1
        end if

        adn=a(jabs)-a(jnb)
        if(adn/=0) px=(a(jabs)-h)/adn
        kx = (jabs - 1) / nrqmax
        ikx = jabs - nrqmax * kx - 1

        if(jx<0) then
          x = drgrd * ikx
          y = dzgrd * (kx - px)
        else
          x = drgrd * (ikx - px)
          y = dzgrd * kx
        end if

        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<=0) then
          ja(1,1) = -jabs-1
          j=2
        end if

        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<=0 .or. ikx<=0) then
          lda=1
          ldb=lda
        else if (ikx + 1 - nr >= 0 .and. jx <= 0) then
          lda=2
          ldb=lda
        else if(jfor/=0) then
          lda=2
          do i=1,3
            if(jfor==ja(i,2)) then
              lda=1
              exit
            end if
          end do
          ldb=lda
        end if

        flag1=.false.
   aa:  do k=1,3
          do l=lda,ldb
            do i=1,ix
              if(lx(i)==ja(k,l)) then
                itm=itm+1
                inext= i
                if(jfor/=0)  exit aa
                if(itm .gt. 3) then
                  flag1=.true.
                  exit aa
                end if
              end if
            end do
          end do
        end do aa

        if(.not.flag1) then
          lx(in)=0
          if(itm .eq. 1) exit
        end if  

        jfor=jx
        jx=lx(inext)
        in = inext
      end do

      do
        if(lx(ix)/=0) then
          if(mpl>=4) then
            ncon = ncon + 1
            npts(ncon) = icount - iclast
            iclast = icount
          end if
          exit
        end if
        ix= ix-1
        if(ix<=0) exit bb
      end do

    end do bb

    if(mpl >= 4) then
      ncon = ncon + 1
      npts(ncon) = icount - iclast
      iclast = icount
    end if
  end subroutine cniteq



  logical function inside_plasma(rrm,zzm)
    use const_and_precisions, only : wp_, zero, one
    use gray_params, only : iequil
    use coreprofiles, only : psdbnd
    use equilibrium, only : zbinf,zbsup,equinum_psi,equian
    implicit none
! arguments
    real(wp_), intent(in) :: rrm,zzm
! local variables
    real(wp_) :: psinv

    if(iequil.eq.1) then
      call equian(rrm,zzm,psinv)
    else
      call equinum_psi(rrm,zzm,psinv)
    end if

    inside_plasma = (psinv >= zero .and. psinv < psdbnd) .and. &
                    (psinv >= one  .or. (zzm >= zbinf .and. zzm <= zbsup))
  end function inside_plasma



  subroutine vacuum_rt(xv0,anv0,xvend,dstvac,ivac)
    use const_and_precisions, only : wp_
    use reflections, only : inters_linewall,inside,rlim,zlim,nlim
    use beamdata, only : dst
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: xv0,anv0
    real(wp_), dimension(3), intent(out) :: xvend
    real(wp_), intent(out) :: dstvac
    integer, intent(out) :: ivac
! local variables
    integer :: i
    real(wp_) :: st,rrm,zzm,smax
    real(wp_), dimension(3) :: walln
    logical :: plfound
! common/external functions/variables
    logical, external :: inside_plasma

! ivac=1 plasma hit before wall reflection
! ivac=2 wall hit before plasma
! ivac=-1 vessel (and thus plasma) never crossed
      
    call inters_linewall(xv0/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim),  &
                         nlim,smax,walln)
    smax=smax*1.0e2_wp_
    rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2)
    zzm=1.0e-2_wp_*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.0_wp_
        xvend=xv0
        ivac=-1
        return
      end if
      ! search second wall interface (inside-outside)
      st=smax
      xvend=xv0+st*anv0
      call inters_linewall(xvend/1.0e2_wp_,anv0,rlim(1:nlim),  &
                           zlim(1:nlim),nlim,smax,walln)
      smax=smax*1.0e2_wp_+st
    end if
    i=0
    do
      st=i*dst
      xvend=xv0+st*anv0
      rrm=1.0e-2_wp_*sqrt(xvend(1)**2+xvend(2)**2)
      zzm=1.0e-2_wp_*xvend(3)
      plfound=inside_plasma(rrm,zzm)
      if (st.ge.smax.or.plfound) exit
      i=i+1
    end do
      
    if (plfound) then
      ivac=1
      dstvac=st
    else
      ivac=2
      dstvac=smax
      xvend=xv0+smax*anv0
    end if
  end subroutine vacuum_rt