gray/src/graycore.f90

module graycore
  use const_and_precisions, only : wp_
  implicit none

contains

  subroutine gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, &
              eqp,psrad,terad,derad,zfc,prfp,  rlim,zlim,                      &
              p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0,           &
              psipol0,chipol0,  dpdv,jcd,pabs,icd,  outp,rtrp,hcdp,ierr, rhout)
    use const_and_precisions, only : zero, one
    use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
    use dispersion, only : expinit
    use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
       rtrparam_type, hcdparam_type, set_codepar, iequil, iprof, ieccd, &
       iwarm, ipec, istpr0, igrad
    use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
    use beamdata, only : pweight, rayi2jk
    use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl,    &
       zbinf, zbsup
    use errcodes, only : check_err, print_errn, print_errhcd
    use magsurf_data, only : flux_average
    use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
    use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
       rhop_tab, rhot_tab
    use limiter, only : set_lim
    use utils, only : vmaxmin
    implicit none
! arguments
    type(eqparam_type), intent(in) :: eqp
    type(prfparam_type), intent(in) :: prfp
    type(outparam_type), intent(in) :: outp
    type(rtrparam_type), intent(in) :: rtrp
    type(hcdparam_type), intent(in) :: hcdp

    real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc
    real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi
    real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim
    real(wp_), dimension(:,:), intent(in) :: psin
    real(wp_), intent(in) :: psia, rvac, rax, zax
    integer, intent(in) :: iox0
    real(wp_), intent(in) :: p0, fghz, psipol0, chipol0
    real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir
    real(wp_), dimension(3), intent(in) :: xv0

    real(wp_), intent(out) :: pabs,icd
    real(wp_), dimension(:), intent(out) :: dpdv,jcd
    real(wp_), dimension(:), intent(in), optional :: rhout

    integer, intent(out) :: ierr
! local variables
    real(wp_), parameter :: taucr = 12._wp_

    real(wp_), dimension(:), allocatable :: rhotn
  
    real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
    real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0
    real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
    real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip
    real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx

    real(wp_), dimension(3) :: xv,anv0,anv
    real(wp_), dimension(:,:), allocatable :: yw,ypw,gri
    real(wp_), dimension(:,:,:), allocatable :: xc,du1,ggri
    integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1
    logical :: ins_pl, somein, allout

    real(wp_),  dimension(:,:), allocatable :: psjki,ppabs,ccci
    real(wp_), dimension(:), allocatable :: tau0,alphaabs0,dids0,ccci0
    real(wp_),    dimension(:), allocatable :: p0jk
    complex(wp_), dimension(:), allocatable :: ext, eyt
    integer,      dimension(:), allocatable :: iiv

    real(wp_),    dimension(:), allocatable :: jphi,pins,currins

! ======= set environment BEGIN ======
    call set_codepar(eqp,prfp,outp,rtrp,hcdp)

    call set_lim(rlim,zlim)

    if(iequil<2) then
      call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3))
    else
      call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, &
                     rvac, rax,zax, rbnd,zbnd, eqp%ixp)
                     ! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol)
    end if
! compute flux surface averaged quantities
    call flux_average ! requires frhotor for dadrhot,dvdrhot

    if(iprof==0) then
      call set_prfan(terad,derad,zfc)
    else
      call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd)
    end if

    call xgygcoeff(fghz,ak0,bres,xgcn)
    call launchangles2n(alpha0,beta0,xv0,anv0)
    call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
                  tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)

    if(iwarm > 1) call expinit

! ======= set environment  END  ======

! ======= pre-proc prints BEGIN ======
    call print_headers
! print  psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
    call print_surfq((/1.5_wp_,2.0_wp_/))
! print  
    print*,' '
    print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
    print'(a,4f8.3)','x00, y00, z00 = ',xv0
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot
    call print_bres(bres)
    call print_prof
! ======= pre-proc prints  END  ======

! ======= main loop BEGIN ======
    iox=iox0
    sox=-1.0_wp_
    if(iox==2) sox=1.0_wp_
    call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
    call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri)
    
    psipol=psipol0
    chipol=chipol0
    call set_pol(yw,bres,sox,psipol,chipol,ext,eyt)
    call pweight(p0,p0jk)

    st0 = zero
    if(nray>1) call print_projxyzt(st0,yw,0)  ! iproj=0 ==> nfilp=8

    somein = .false. ! becomes true if at least part of the beam enters the plasma
! beam/ray propagation
    do i=1,nstep

! advance one step with "frozen" grad(S_I)
      st=i*dst+st0
      do jk=1,nray
        call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk))
      end do

! update position and grad
      if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)

      allout = .true. ! becomes false if at least part of the beam is inside the plsama
      ierr = 0
      do jk=1,nray
! compute derivatives with updated gradient and local plasma values
        xv = yw(1:3,jk)
        anv = yw(4:6,jk)
        call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), &
           psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn)
! update global error code and print message
        if (ierrn/=0) then
          ierr = ior(ierr,ierrn)
          call print_errn(ierrn,i,anpl)
        end if

        zzm = xv(3)*0.01_wp_
        ins_pl = (psinv>=zero .and. psinv<one .and. zzm>=zbinf .and. zzm<=zbsup)
! test if the beam is completely out of the plsama
        allout = allout .and. .not.ins_pl
! test if at least part of the beam has entered the plsama
        somein = somein .or. ins_pl

! compute ECRH&CD
        if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. tau0(jk)<=taucr) then
          tekev=temp(psinv)
          call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
             sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) 
          if (ierrhcd/=0) then
            ierr = ior(ierr,ierrhcd)
            call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
          end if
        else
          tekev=zero
          alpha=zero
          didp=zero
          anprim=zero
          anprre=anpr
          nharm=0
          nhf=0
          iokhawa=0
        end if
        if(nharm>0) iiv(jk)=i

        psjki(jk,i) = psinv
! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
        tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst
        pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk))
        ppabs(jk,i)=p0jk(jk)-pow

        dids=didp*pow*alpha
        ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst

        tau0(jk)=tau
        alphaabs0(jk)=alpha
        dids0(jk)=dids
        ccci0(jk)=ccci(jk,i)

        call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr,     &
           anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)

      end do

! print ray positions for j=nrayr in local reference system
      if (mod(i,istpr0) == 0) then
        if(nray > 1) call print_projxyzt(st,yw,0)
      end if

! check for any error code and stop if necessary
      call check_err(ierr,istop)
! test whether further trajectory integration is unnecessary
      call vmaxmin(tau0,nray,taumn,taumx)
      if ((taumn > taucr) .or. (somein .and. allout)) istop = 1

      if(istop == 1) exit
    end do

! compute total absorbed power and driven current
    if (i>nstep) i=nstep
    pabs = sum(ppabs(:,i))
    icd  = sum(ccci(:,i))
! ======= main loop  END  ======

! ======= post-proc BEGIN ======
! print final results on screen
    write(*,*)
    write(*,'(a,f9.4)') 'final step (s, ct, Sr)  = ',st
    write(*,'(a,2e12.5)') 'taumn, taumx  = ', taumn,taumx
    write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs
    write(*,'(a,f9.4)') 'I_tot (kA)    = ',icd*1.0e3_wp_
! print all ray positions in local reference system
    if(nray > 1) call print_projxyzt(st,yw,1)

! compute power and current density profiles for all rays
    call pec_init(ipec,rhout)
    nnd=size(rhop_tab)
    allocate(jphi(nnd),pins(nnd),currins(nnd))
    call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins)
! print power and current density profiles
    call print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
! compute profiles width
    call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi,           &
       rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,           &
       rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx)
! print 0D results
    call print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav,      &
        drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
        st,psipol,chipol,index_rt)


! ======= post-proc  END  ======

! ======= free memory BEGIN ======
!    call unset_eqspl
!    call unset_q
!    call unset_rhospl
!    call unset_prfspl
!    call unset_lim
!    call dealloc_surfvec
!    call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
!                      tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
!    call dealloc_pec
    deallocate(jphi,pins,currins)
! ======= free memory  END  ======
  end subroutine gray_main



  subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
    use const_and_precisions, only : wp_, zero
    implicit none
! arguments
    real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
    real(wp_), dimension(:), intent(out)   :: tau0,alphaabs0,dids0,ccci0
    integer, dimension(:), intent(out) :: iiv
!! common/external functions/variables
!      integer :: jclosest
!      real(wp_), dimension(3) :: anwcl,xwcl
!
!      common/refln/anwcl,xwcl,jclosest
!
!      jclosest=nrayr+1
!      anwcl(1:3)=0.0_wp_
!      xwcl(1:3)=0.0_wp_

      psjki     = zero
      ppabs     = zero
      ccci      = zero
      tau0      = zero
      alphaabs0 = zero
      dids0     = zero
      ccci0     = zero
      iiv       = 1

  end subroutine vectinit



  subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, &
      ywrk0,ypwrk0,xc0,du10,gri,ggri)
!   beam tracing initial conditions  igrad=1
!   !!!!!! check ray tracing initial conditions   igrad=0 !!!!!!
    use const_and_precisions, only : wp_,izero,zero,one,pi,half,two,degree,ui=>im
    use math, only : catand
    use gray_params, only : idst
    use beamdata, only : nray,nrayr,nrayth,rwmax
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: xv0c,anv0c
    real(wp_), intent(in) :: ak0
    real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir
    real(wp_), dimension(6,nray), intent(out) :: ywrk0,ypwrk0
    real(wp_), dimension(3,nray), intent(out) :: gri
    real(wp_), dimension(3,3,nray), intent(out) :: ggri
    real(wp_), dimension(3,nrayth,nrayr), intent(out) :: xc0,du10
      
! local variables
    integer :: j,k,jk
    real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,alfak
    real(wp_) :: wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy
    real(wp_) :: rcixx,rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy
    real(wp_) :: drcixx,drciyy,drcixy,dr,da,ddfu,dcsiw,detaw,dx0t,dy0t
    real(wp_) :: x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2
    real(wp_) :: gxxt,gyyt,gzzt,gxyt,gxzt,gyzt,dgr2xt,dgr2yt,dgr2zt
    real(wp_) :: dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy
    real(wp_) :: anzt,anxt,anyt,anx,any,anz,an20,an0
    real(wp_) :: du1tx,du1ty,du1tz,denom,ddr,ddi
    real(wp_), dimension(nrayr) :: uj
    real(wp_), dimension(nrayth) :: sna,csa
    complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2
    complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy

    csth=anv0c(3)
    snth=sqrt(one-csth**2)
    if(snth > zero) then
      csps=anv0c(2)/snth
      snps=anv0c(1)/snth
    else
      csps=one
      snps=zero
    end if

!  Gaussian beam: exp[-ik0 zt] exp[-i k0/2 S(xt,yt,zt)]
!  xt,yt,zt, cartesian coordinate system with zt along the beamline and xt in the z = 0 plane
!  S(xt,yt,zt) = S_real +i S_imag = Qxx(zt) xt^2 + Qyy(zt) yt^2 + 2 Qxy(zt) xt yt
!  (csiw, etaw) and (csiR, etaR) intensity and phase ellipse, rotated by angle phiw and phiR
!  S(xt,yt,zt) = csiR^2 / Rccsi +etaR^2 /Rceta - i (csiw^2 Wcsi +etaw^2 Weta)
!  Rccsi,eta curvature radius  at the launching point 
!  Wcsi,eta =2/(k0 wcsi,eta^2) with  wcsi,eta^2 beam size at the launching point

    phiwrad = phiw*degree
    phirrad = phir*degree
    csphiw = cos(phiwrad)
    snphiw = sin(phiwrad)
!    csphir = cos(phirrad)
!    snphir = sin(phirrad)

    wwcsi = two/(ak0*wcsi**2)
    wweta = two/(ak0*weta**2)

    if(phir/=phiw) then
      sk   =  rcicsi + rcieta
      sw   =  wwcsi  + wweta
      dk   =  rcicsi - rcieta
      dw   =  wwcsi  - wweta
      ts   = -(dk*sin(2*phirrad)       - ui*dw*sin(2*phiwrad))
      tc   =  (dk*cos(2*phirrad)       - ui*dw*cos(2*phiwrad))
      phic =  half*catand(ts/tc)
      ddd  =  dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic))
      sss  =  sk - ui*sw
      qi1  =  half*(sss + ddd)
      qi2  =  half*(sss - ddd)
      rci1 =  dble(qi1)
      rci2 =  dble(qi2)
      ww1  = -dimag(qi1)
      ww2  = -dimag(qi2)
    else
      rci1 = rcicsi
      rci2 = rcieta
      ww1 = wwcsi
      ww2 = wweta
      phic = -phiwrad
      qi1 = rci1 - ui*ww1
      qi2 = rci2 - ui*ww2
    end if
      
!    w01=sqrt(2.0_wp_/(ak0*ww1))
!    d01=-rci1/(rci1**2+ww1**2)
!    w02=sqrt(2.0_wp_/(ak0*ww2))
!    d02=-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(phic)*cos(phic)
    wwxx  = -dimag(qqxx)
    wwyy  = -dimag(qqyy)
    wwxy  = -dimag(qqxy)
    rcixx =  dble(qqxx)
    rciyy =  dble(qqyy)
    rcixy =  dble(qqxy)

    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(phic)*cos(phic)
    d2qqxx =  d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2
    d2qqyy =  d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2
    d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic)

    dwwxx  = -dimag(dqqxx)
    dwwyy  = -dimag(dqqyy)
    dwwxy  = -dimag(dqqxy)
    d2wwxx = -dimag(d2qqxx)
    d2wwyy = -dimag(d2qqyy)
    d2wwxy = -dimag(d2qqxy)
    drcixx =  dble(dqqxx)
    drciyy =  dble(dqqyy)
    drcixy =  dble(dqqxy)

    if(nrayr > 1) then
      dr = rwmax/dble(nrayr-1)
    else
      dr = one
    end if
    ddfu = two*dr**2/ak0   ! twodr2 = 2*dr**2 = 2*rwmax/dble(nrayr-1)
    do j = 1, nrayr
      uj(j)     = dble(j-1)
    end do  

    da=2*pi/dble(nrayth)
    do k=1,nrayth
      alfak  = (k-1)*da
      sna(k) = sin(alfak)
      csa(k) = cos(alfak)
    end do

!   central ray
    jk=1
    gri(:,1)    = zero
    ggri(:,:,1) = zero

    ywrk0(1:3,1)  = xv0c
    ywrk0(4:6,1)  = anv0c
    ypwrk0(1:3,1) = anv0c
    ypwrk0(4:6,1) = zero
      
    do k=1,nrayth
      dcsiw = dr*csa(k)*wcsi
      detaw = dr*sna(k)*weta
      dx0t  = dcsiw*csphiw - detaw*snphiw
      dy0t  = dcsiw*snphiw + detaw*csphiw
      du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu
      du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu

      xc0(:,k,1)  =  xv0c
      du10(1,k,1) =  du1tx*csps + snps*du1ty*csth
      du10(2,k,1) = -du1tx*snps + csps*du1ty*csth
      du10(3,k,1) = -du1ty*snth
    end do
    ddr = zero
    ddi = zero

!   loop on rays jk>1      
    j=2
    k=0
    do jk=2,nray
      k=k+1
      if(k > nrayth) then
        j=j+1
        k=1
      end if

!      csiw=u*dcsiw
!      etaw=u*detaw
!      csir=x0t*csphir+y0t*snphir
!      etar=-x0t*snphir+y0t*csphir
      dcsiw = dr*csa(k)*wcsi
      detaw = dr*sna(k)*weta
      dx0t  = dcsiw*csphiw - detaw*snphiw
      dy0t  = dcsiw*snphiw + detaw*csphiw
      x0t   = uj(j)*dx0t
      y0t   = uj(j)*dy0t
      z0t   = -(half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t)

      dx0 =  x0t*csps + snps*(y0t*csth + z0t*snth)
      dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth)
      dz0 =  z0t*csth -       y0t*snth
      x0 = xv0c(1) + dx0
      y0 = xv0c(2) + dy0
      z0 = xv0c(3) + dz0

      gxt  = x0t*wwxx + y0t*wwxy
      gyt  = x0t*wwxy + y0t*wwyy
      gzt  = half*(x0t**2*dwwxx  + y0t**2*dwwyy ) + x0t*y0t*dwwxy
      gr2  = gxt*gxt + gyt*gyt + gzt*gzt
      gxxt = wwxx
      gyyt = wwyy
      gzzt = half*(x0t**2*d2wwxx + y0t**2*d2wwyy) + x0t*y0t*d2wwxy
      gxyt = wwxy
      gxzt = x0t*dwwxx + y0t*dwwxy
      gyzt = x0t*dwwxy + y0t*dwwyy
      dgr2xt = 2*(gxt*gxxt + gyt*gxyt + gzt*gxzt)
      dgr2yt = 2*(gxt*gxyt + gyt*gyyt + gzt*gyzt)
      dgr2zt = 2*(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,jk) =  gxt*csps + snps*(gyt*csth + gzt*snth)
      gri(2,jk) = -gxt*snps + csps*(gyt*csth + gzt*snth)
      gri(3,jk) =  gzt*csth - gyt*snth
      ggri(1,1,jk) = gxxt*csps**2                                               &
                  +  snps**2  *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) &
                  +2*snps*csps*(gxyt*csth + gxzt*snth)
      ggri(2,1,jk) = csps*snps                                                  &
                  *(-gxxt+csth**2*gyyt + snth**2*gzzt + 2*csth*snth*gyzt)       &
                  +(csps**2 - snps**2)*(snth*gxzt + csth*gxyt)
      ggri(3,1,jk) = csth*snth*snps*(gzzt - gyyt) + (csth**2 - snth**2)         &
                  *snps*gyzt + csps*(csth*gxzt - snth*gxyt)
      ggri(1,2,jk) = ggri(2,1,jk)
      ggri(2,2,jk) = gxxt*snps**2                                               &
                  +  csps**2  *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) &
                  -2*snps*csps*(gxyt*csth    + gxzt*snth)
      ggri(3,2,jk) = csth*snth*csps*(gzzt - gyyt) + (csth**2-snth**2)           &
                  *csps*gyzt + snps*(snth*gxyt - csth*gxzt)
      ggri(1,3,jk) = ggri(3,1,jk)
      ggri(2,3,jk) = ggri(3,2,jk)
      ggri(3,3,jk) = gzzt*csth**2 + gyyt*snth**2 - 2*csth*snth*gyzt

      du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu
      du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu
      du1tz = half*uj(j)*(dx0t**2*dwwxx + dy0t**2*dwwyy + 2*dx0t*dy0t*dwwxy)/ddfu

      du10(1,k,j) =  du1tx*csps + snps*(du1ty*csth + du1tz*snth)
      du10(2,k,j) = -du1tx*snps + csps*(du1ty*csth + du1tz*snth)
      du10(3,k,j) =  du1tz*csth -       du1ty*snth

      pppx  = x0t*rcixx + y0t*rcixy
      pppy  = x0t*rcixy + y0t*rciyy
      denpp = pppx*gxt + pppy*gyt
      if (denpp/=zero) then
        ppx = -pppx*gzt/denpp
        ppy = -pppy*gzt/denpp
      else
        ppx =  zero
        ppy =  zero
      end if

      anzt = sqrt((one + gr2)/(one + ppx**2 + ppy**2))
      anxt = ppx*anzt
      anyt = ppy*anzt

      anx  = anxt*csps + snps*(anyt*csth + anzt*snth)
      any  =-anxt*snps + csps*(anyt*csth + anzt*snth)
      anz  = anzt*csth - anyt*snth

      an20 = one + gr2
      an0  = sqrt(an20)

      xc0(1,k,j) = x0
      xc0(2,k,j) = y0
      xc0(3,k,j) = z0

      ywrk0(1,jk) = x0
      ywrk0(2,jk) = y0
      ywrk0(3,jk) = z0
      ywrk0(4,jk) = anx
      ywrk0(5,jk) = any
      ywrk0(6,jk) = anz

      select case(idst) 
      case(1)
! integration variable: c*t
        denom = one
      case(2)
! integration variable: Sr
        denom = an20
      case default ! idst=0
! integration variable: s
        denom = an0
      end select
      ypwrk0(1,jk) = anx/denom
      ypwrk0(2,jk) = any/denom
      ypwrk0(3,jk) = anz/denom
      ypwrk0(4,jk) = dgr2x/(2*denom)
      ypwrk0(5,jk) = dgr2y/(2*denom)
      ypwrk0(6,jk) = dgr2z/(2*denom)

      ddr = anx**2 + any**2 + anz**2 - an20
      ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
      call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,ak0,zero,zero,zero, &
           zero,zero,zero,zero,zero,0,0,1,ddr,ddi) ! st=0, index_rt=1, Btot=0, psin=-1
    end do
  end subroutine ic_gb



  subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr)
!     Runge-Kutta integrator
    use const_and_precisions, only : wp_
!    use gray_params, only : igrad
    use beamdata, only : h,hh,h6
    implicit none
    real(wp_), intent(in) :: sox,bres,xgcn
    real(wp_), dimension(6), intent(inout) :: y
    real(wp_), dimension(6), intent(in) :: yp
    real(wp_), dimension(3), intent(in) :: dgr
    real(wp_), dimension(3,3), intent(in) :: ddgr

    real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
    real(wp_) :: gr2
    real(wp_), dimension(3) :: dgr2
    
!    if(igrad.eq.1) then
      gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
      dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
!    end if
    fk1 = yp

    yy = y + fk1*hh
    call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2)
    yy = y + fk2*hh
    call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3)
    yy = y + fk3*h
    call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4)

    y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
  end subroutine rkstep



  subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery)
! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator
    use const_and_precisions, only : wp_
    implicit none
! arguments
    real(wp_), dimension(6), intent(in) :: y
    real(wp_), intent(in) :: sox,bres,xgcn,gr2
    real(wp_), dimension(3), intent(in) :: dgr2,dgr
    real(wp_), dimension(3,3), intent(in) :: ddgr
    real(wp_), dimension(6), intent(out) :: dery
! local variables
    real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi
    real(wp_) :: ddr,ddi,dersdst,derdnm
    real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg
    real(wp_), dimension(3,3) :: derbv

    xv = y(1:3)
    call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg, &
       ajphi)

    anv = y(4:6)
    call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
       dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
  end subroutine rhs



  subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
     xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr)
! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update
    use errcodes, only : pnpl
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: xv,anv
    real(wp_), dimension(3), intent(in) :: dgr
    real(wp_), dimension(3,3), intent(in) :: ddgr
    real(wp_), intent(in) :: sox,bres,xgcn
    real(wp_), dimension(6), intent(out) :: dery
    real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr
    real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm
    integer, intent(out) :: ierr
! local variables
    real(wp_) :: gr2,ajphi
    real(wp_), dimension(3) :: dgr2,bv,derxg,deryg
    real(wp_), dimension(3,3) :: derbv
    real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_

    gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
    dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
    call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi)
    call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
       dery,anpl,anpr,ddr,ddi,dersdst,derdnm)

    ierr=0
    if(  abs(anpl) > anplth1) then
      if(abs(anpl) > anplth2) then
        ierr=ibset(ierr,pnpl+1)
      else 
        ierr=ibset(ierr,pnpl)
      end if 
    end if
  end subroutine ywppla_upd



  subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri)
    use const_and_precisions, only : wp_,zero,half
    use beamdata, only : nray,nrayr,nrayth,twodr2
    implicit none
    real(wp_), intent(in) :: ak0
    real(wp_), dimension(6,nray), intent(in) :: ywrk
    real(wp_), dimension(3,nrayth,nrayr), intent(inout) :: xc,du1
    real(wp_), dimension(3,nray), intent(out) :: gri
    real(wp_), dimension(3,3,nray), intent(out) :: ggri
! local variables
    real(wp_), dimension(3,nrayth,nrayr) :: xco,du1o
    integer :: jk,j,jm,jp,k,km,kp
    real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz
    real(wp_) :: dfuu,dffiu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz
    real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu
    real(wp_), dimension(3,3) :: dgg,dff

!     update position and du1 vectors
    xco = xc
    du1o = du1

    jk = 1
    do j=1,nrayr
      do k=1,nrayth
        if(j>1) jk=jk+1
        xc(1:3,k,j)=ywrk(1:3,jk)
      end do
    end do

!     compute  grad u1 for central ray
    j = 1
    jp = 2
    do k=1,nrayth
      if(k == 1) then
        km = nrayth
      else
        km = k-1
      end if
      if(k == nrayth) then
        kp = 1
      else
        kp = k+1
      end if
      dxv1 = xc(:,k ,jp) -  xc(:,k ,j)
      dxv2 = xc(:,kp,jp) -  xc(:,km,jp)
      dxv3 = xc(:,k ,j)  - xco(:,k ,j)
      call solg0(dxv1,dxv2,dxv3,dgu)
      du1(:,k,j) = dgu
    end do
    gri(:,1) = zero
      
!     compute  grad u1 and grad(S_I) for all the other rays
    dfuu=twodr2/ak0      ! twodr2 = 2*dr**2 = 2*(rwmax/(nrayr-1))**2
    jm=1
    j=2
    k=0
    dffiu = dfuu
    do jk=2,nray
      k=k+1
      if(k > nrayth) then
        jm = j
        j = j+1
        k = 1
        dffiu = dfuu*jm
      end if
      kp = k+1
      km = k-1
      if (k == 1) then
        km=nrayth
      else if (k == nrayth) then
        kp=1
      end if
      dxv1 = xc(:,k ,j) -  xc(:,k ,jm)
      dxv2 = xc(:,kp,j) -  xc(:,km,j)
      dxv3 = xc(:,k ,j) - xco(:,k ,j)
      call solg0(dxv1,dxv2,dxv3,dgu)
      du1(:,k,j) = dgu
      gri(:,jk) = dgu(:)*dffiu
    end do

!     compute derivatives of grad u and grad(S_I) for rays jk>1
    ggri(:,:,1) = zero
    jm=1
    j=2
    k=0
    dffiu = dfuu
    do jk=2,nray
      k=k+1
      if(k > nrayth) then
        jm=j
        j=j+1
        k=1
        dffiu = dfuu*jm
      end if
      kp=k+1
      km=k-1
      if (k == 1) then
        km=nrayth
      else if (k == nrayth) then
        kp=1
      end if
      dxv1 = xc(:,k ,j) -  xc(:,k ,jm)
      dxv2 = xc(:,kp,j) -  xc(:,km,j)
      dxv3 = xc(:,k ,j) - xco(:,k ,j)
      dff(:,1) = du1(:,k ,j) -  du1(:,k ,jm)
      dff(:,2) = du1(:,kp,j) -  du1(:,km,j)
      dff(:,3) = du1(:,k ,j) - du1o(:,k ,j)
      call solg3(dxv1,dxv2,dxv3,dff,dgg)

!     derivatives of u
      ux = du1(1,k,j)
      uy = du1(2,k,j)
      uz = du1(3,k,j)
      uxx = dgg(1,1)
      uyy = dgg(2,2)
      uzz = dgg(3,3)
      uxy = (dgg(1,2) + dgg(2,1))*half
      uxz = (dgg(1,3) + dgg(3,1))*half
      uyz = (dgg(2,3) + dgg(3,2))*half

!     derivatives of S_I and Grad(S_I)
      gx = ux*dffiu
      gy = uy*dffiu
      gz = uz*dffiu
      gxx = dfuu*ux*ux + dffiu*uxx
      gyy = dfuu*uy*uy + dffiu*uyy
      gzz = dfuu*uz*uz + dffiu*uzz
      gxy = dfuu*ux*uy + dffiu*uxy
      gxz = dfuu*ux*uz + dffiu*uxz
      gyz = dfuu*uy*uz + dffiu*uyz

      ggri(1,1,jk)=gxx
      ggri(2,1,jk)=gxy
      ggri(3,1,jk)=gxz
      ggri(1,2,jk)=gxy
      ggri(2,2,jk)=gyy
      ggri(3,2,jk)=gyz
      ggri(1,3,jk)=gxz
      ggri(2,3,jk)=gyz
      ggri(3,3,jk)=gzz
    end do
      
  end subroutine gradi_upd



  subroutine solg0(dxv1,dxv2,dxv3,dgg)
!     solution of the linear system of 3 eqs : dgg . dxv = dff
!     input vectors : dxv1, dxv2, dxv3, dff
!     output vector : dgg
!     dff=(1,0,0)
    use const_and_precisions, only : wp_
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
    real(wp_), dimension(3), intent(out) :: dgg
! local variables
    real(wp_) :: denom,aa1,aa2,aa3

    aa1    = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3))
    aa2    = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3))
    aa3    = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3))

    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
  end subroutine solg0

  subroutine solg3(dxv1,dxv2,dxv3,dff,dgg)
!      rhs "matrix" dff, result in dgg
    use const_and_precisions, only : wp_
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
    real(wp_), dimension(3,3), intent(in) :: dff
    real(wp_), dimension(3,3), intent(out) :: dgg
! local variables
    real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33

    a11 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3))
    a21 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3))
    a31 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3))

    a12 = (dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))
    a22 = (dxv1(1)*dxv3(3) - dxv3(1)*dxv1(3))
    a32 = (dxv1(1)*dxv2(3) - dxv2(1)*dxv1(3))

    a13 = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))
    a23 = (dxv1(1)*dxv3(2) - dxv3(1)*dxv1(2))
    a33 = (dxv1(1)*dxv2(2) - dxv2(1)*dxv1(2))

    denom = dxv1(1)*a11 - dxv2(1)*a21 + dxv3(1)*a31

    dgg(:,1) = ( dff(:,1)*a11 - dff(:,2)*a21 + dff(:,3)*a31)/denom
    dgg(:,2) = (-dff(:,1)*a12 + dff(:,2)*a22 - dff(:,3)*a32)/denom
    dgg(:,3) = ( dff(:,1)*a13 - dff(:,2)*a23 + dff(:,3)*a33)/denom
  end subroutine solg3



  subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, &
     xg,yg,derxg,deryg,ajphi)
    use const_and_precisions, only : wp_,zero,pi,ccj=>mu0inv
    use gray_params, only : iequil
    use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi
    use coreprofiles, only : density
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: xv
    real(wp_), intent(in) :: xgcn,bres
    real(wp_), intent(out) :: psinv,dens,btot,xg,yg
    real(wp_), dimension(3), intent(out) :: bv,derxg,deryg
    real(wp_), dimension(3,3), intent(out) :: derbv
! local variables
    integer :: jv
    real(wp_) ::  xx,yy,zz
    real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm
    real(wp_), dimension(3) :: dbtot,bvc
    real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
    real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi
    real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin

    xg = zero
    yg = 99._wp_
    psinv = -1._wp_
    dens = zero
    btot = zero
    ajphi = zero
    derxg = zero
    deryg = zero
    bv = zero
    derbv = zero

    if(iequil==0) return

    dbtot = zero
    dbv = zero
    dbvcdc = zero
    dbvcdc = zero
    dbvdc = zero

    xx = xv(1)
    yy = xv(2)
    zz = xv(3)

!     cylindrical coordinates
    rr2 = xx**2 + yy**2
    rr = sqrt(rr2)
    csphi = xx/rr
    snphi = yy/rr

    bv(1) = -snphi*sgnbphi
    bv(2) =  csphi*sgnbphi

!     convert from cm to meters
    zzm = 1.0e-2_wp_*zz
    rrm = 1.0e-2_wp_*rr

    if(iequil==1) then
      call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, &
                  ddpsidrr,ddpsidzz,ddpsidrz)
    else
      call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz)
      call equinum_fpol(psinv,fpolv,dfpv)
    end if

!     compute yg and derivative
    if(psinv < zero) then
      bphi = fpolv/rrm
      btot = abs(bphi)
      yg   = btot/bres
      return
    end if

!     compute xg and derivative
    call density(psinv,dens,ddenspsin)
    xg      = xgcn*dens
    dxgdpsi = xgcn*ddenspsin/psia
      
!     B = f(psi)/R e_phi+  grad psi x e_phi/R
    bphi =  fpolv/rrm
    brr  =-dpsidz/rrm
    bzz  = dpsidr/rrm

!     bvc(i) = B_i in cylindrical coordinates
    bvc(1) = brr
    bvc(2) = bphi
    bvc(3) = bzz

!     bv(i)  = B_i in cartesian coordinates
    bv(1)=bvc(1)*csphi - bvc(2)*snphi
    bv(2)=bvc(1)*snphi + bvc(2)*csphi
    bv(3)=bvc(3)

!     dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv)
    dbvcdc(1,1) =   -ddpsidrz/rrm - brr/rrm
    dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm
    dbvcdc(3,1) =    ddpsidrr/rrm - bzz/rrm
    dbvcdc(1,3) =   -ddpsidzz/rrm
    dbvcdc(2,3) = dfpv*dpsidz/rrm
    dbvcdc(3,3) =    ddpsidrz/rrm

!     dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv)
    dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi
    dbvdc(2,1) = dbvcdc(1,1)*snphi + dbvcdc(2,1)*csphi
    dbvdc(3,1) = dbvcdc(3,1)
    dbvdc(1,2) = -bv(2)
    dbvdc(2,2) =  bv(1)
    dbvdc(3,2) = dbvcdc(3,2)
    dbvdc(1,3) = dbvcdc(1,3)*csphi - dbvcdc(2,3)*snphi
    dbvdc(2,3) = dbvcdc(1,3)*snphi + dbvcdc(2,3)*csphi
    dbvdc(3,3) = dbvcdc(3,3)

    drrdx  =  csphi
    drrdy  =  snphi
    dphidx = -snphi/rrm
    dphidy =  csphi/rrm

!     dbv(iv,jv) = d Bcart(iv) / dxvcart(jv)
    dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2)
    dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2)
    dbv(:,3) = dbvdc(:,3)

!     B magnitude and derivatives 
    b2tot = bv(1)**2 + bv(2)**2 + bv(3)**2
    btot  = sqrt(b2tot)

    dbtot = (bv(1)*dbv(1,:) + bv(2)*dbv(2,:) + bv(3)*dbv(3,:))/btot

    yg = btot/Bres

!     convert spatial derivatives from dummy/m -> dummy/cm
!     to be used in rhs

!     bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z
    deryg = 1.0e-2_wp_*dbtot/Bres
    bv    = bv/btot
    do jv=1,3
      derbv(:,jv) = 1.0e-2_wp_*(dbv(:,jv) - bv(:)*dbtot(jv))/btot
    end do

    derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi
    derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi
    derxg(3) = 1.0e-2_wp_*dpsidz      *dxgdpsi

!     current density computation in Ampere/m^2, ccj==1/mu_0
    ajphi = ccj*(dbvcdc(1,3) - dbvcdc(3,1))
!      ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3))
!      ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2))

  end subroutine plas_deriv



  subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
     dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
    use const_and_precisions, only : wp_,zero,one,half,two
    use gray_params, only : idst,igrad
    implicit none
! arguments
    real(wp_), intent(in) :: xg,yg,gr2,sox
    real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst
    real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg
    real(wp_), dimension(3), intent(in) :: dgr2,dgr
    real(wp_), dimension(3,3), intent(in)  :: ddgr,derbv
    real(wp_), dimension(6), intent(out) :: dery
! local variables
    integer :: iv
    real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s
    real(wp_) :: dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel
    real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm
    real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv

    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)

    anpl2 = anpl**2
    dnl   = one - anpl2
    anpr2 = max(an2-anpl2,zero)
    anpr  = sqrt(anpr2)
    yg2   = yg**2

    an2s      = one
    dan2sdxg  = zero
    dan2sdyg  = zero
    dan2sdnpl = zero
    del       = zero
    fdia      = zero
    dfdiadnpl = zero
    dfdiadxg  = zero
    dfdiadyg  = zero

    duh = one - xg - yg2
    if(xg > zero) then
      del  = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2)
      an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh

      dan2sdxg  = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2    &
                  + sox*xg*anpl2/(del*duh) - one
      dan2sdyg  = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2        &
                  + two*sox*xg*(one - xg)*anpl2/(yg*del*duh)
      dan2sdnpl = - xg*yg2*anpl/duh                                        &
                  - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh)

      if(igrad > 0) then
        ddelnpl2  = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2)           &
                         - yg2*dnl**3)/yg2/del**3
        fdia      = - xg*yg2*(one + half*sox*ddelnpl2)/duh
        derdel    = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2)          &
                    - dnl**2*(one + 3.0_wp_*anpl2)*yg2
        derdel    = 4.0_wp_*derdel/(yg*del)**5
        ddelnpl2y = two*(one - xg)*derdel
        ddelnpl2x = yg*derdel
        dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2)       &
                    /(yg2*del**5)
        dfdiadxg  = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2)        &
                    *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2)
        dfdiadyg  = - two*yg*xg*(one - xg)/duh**2                          &
                    - sox*xg*yg*(two*(one - xg)*ddelnpl2                   &
                                 + yg*duh*ddelnpl2y)/(two*duh**2)
      end if
    end if

    bdotgr = bv(1)*dgr(1) + bv(2)*dgr(2) + bv(3)*dgr(3)
    do iv=1,3
      dbgr(iv)     = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv)                 &
                   + dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv)                 &
                   + dgr(3)*derbv(3,iv) + bv(3)*ddgr(3,iv)
      danpldxv(iv) = anv(1)*derbv(1,iv) + anv(2)*derbv(2,iv) + anv(3)*derbv(3,iv)
    end do

    derdxv   = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl +    &
                 igrad*dgr2)                                               &
               + fdia*bdotgr*dbgr + half*bdotgr**2                         &
               *(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl)
    derdnv   = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv

    derdnm   = sqrt(derdnv(1)**2 + derdnv(2)**2 + derdnv(3)**2)

    derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl     &
             + two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg               &
                                          + half*yg*dfdiadyg               &
                                          + half*anpl*dfdiadnpl)

    if (idst == 0) then
!   integration variable: s
      denom = derdnm
    else if (idst == 1) then
!   integration variable: c*t
      denom = -derdom
    else
!   integration variable: Sr
      denom = anv(1)*derdnv(1) + anv(2)*derdnv(2) + anv(3)*derdnv(3)
    end if

!   coefficient for integration in s
!   ds/dst, where st is the integration variable
    dersdst = derdnm/denom

!   rhs vector
    dery(1:3) =  derdnv(:)/denom
    dery(4:6) = -derdxv(:)/denom

!     vgv :  ~ group velocity
!      vgm=0
!      do iv=1,3
!        vgv(iv)=-derdnv(iv)/derdom
!        vgm=vgm+vgv(iv)**2
!      end do
!      vgm=sqrt(vgm)

!     ddr :  dispersion relation (real part)
!     ddi :  dispersion relation (imaginary part)
    ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia)
    ddi = derdnv(1)*dgr(1) + derdnv(2)*dgr(2) + derdnv(3)*dgr(3)

  end subroutine disp_deriv



  subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
     sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr)
    use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_
    use gray_params, only : iwarm,ilarm,ieccd,imx
    use coreprofiles, only : fzeff
    use equilibrium, only : sgnbphi
    use dispersion, only : harmnumber, warmdisp
    use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl
    use errcodes, only : palph
    use magsurf_data, only : fluxval
    implicit none
! arguments
    real(wp_),intent(in) ::psinv,ak0,bres
    real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox 
    real(wp_),intent(out) :: anprre,anprim,alpha,didp
    integer, intent(out) :: nhmin,nhmax,iokhawa
    integer, intent(out) :: ierr
! local constants
    real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_
! local variables
    real(wp_) :: rbavi,rrii,rhop
    integer :: lrm,ithn,ierrcd
    real(wp_) :: amu,ratiovgr,rbn,rbx
    real(wp_) :: zeff,cst2,bmxi,bmni,fci
    real(wp_), dimension(:), allocatable :: eccdpar
    real(wp_) :: effjcd,effjcdav,akim,btot
    complex(wp_) :: ex,ey,ez

    alpha=zero
    anprim=zero
    anprre=zero
    didp=zero
    nhmin=0
    nhmax=0
    iokhawa=0
    ierr=0

    if(tekev>zero) then
!     absorption computation
      amu=mc2/tekev
      call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm)
      if(nhmin.gt.0) then
        lrm=max(ilarm,nhmax)
        call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim,  &
                      iwarm,imx,ex,ey,ez)
        akim=ak0*anprim
        ratiovgr=2.0_wp_*anpr/derdnm!*vgm
        alpha=2.0_wp_*akim*ratiovgr
        if(alpha<zero) then
          ierr=ibset(ierr,palph)
          return
        end if
          
!       calcolo della efficienza <j/p>:  effjcdav      [A m/W ]
        if(ieccd>0) then
!     current drive computation
          zeff=fzeff(psinv)
          ithn=1
          if(lrm>nhmin) ithn=2
          rhop=sqrt(psinv)
          call fluxval(rhop,rri=rrii,rbav=rbavi,bmn=bmni,bmx=bmxi,fc=fci)
          btot=yg*bres
          rbn=btot/bmni
          rbx=btot/bmxi

          select case(ieccd) 
          case(1)
!             cohen model      
            call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar)
            call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
                         ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd)
          case(2)
!             no trapping
            call setcdcoeff(zeff,cst2,eccdpar)
            call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
                         ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd)
          case default
!             neoclassical model
            call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar)
            call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
                         ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd)
          end select
          ierr=ierr+ierrcd
          !deallocate(eccdpar)
          effjcdav=rbavi*effjcd
          didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii)
        end if
      end if
    end if
  end subroutine alpha_effj



  subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0)
    use const_and_precisions, only : wp_,degree,zero,one,half,im
    use beamdata, only : nray,nrayth
    use equilibrium, only : bfield
    use gray_params, only : ipol
    use polarization, only : pol_limit, polellipse, stokes_ce, stokes_ell
    implicit none
! arguments
    real(wp_), dimension(6,nray), intent(in) :: ywrk0
    real(wp_), intent(in) :: sox,bres
    real(wp_), intent(inout) :: psipol0, chipol0
    complex(wp_), dimension(nray), intent(out) :: ext0, eyt0
! local variables
    integer :: j,k,jk
    real(wp_), dimension(3) :: xmv, anv, bv
    real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol

    j=1
    k=0
    do jk=1,nray
      k=k+1
      if(jk == 2 .or. k > nrayth) then
        j=j+1
        k=1
      end if

      if(ipol == 0) then
        xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m
        anv=ywrk0(4:6,jk)
        rm=sqrt(xmv(1)**2+xmv(2)**2)
        csphi=xmv(1)/rm
        snphi=xmv(2)/rm
        call bfield(rm,xmv(3),bphi,br,bz)
! bv(i)  = B_i in cartesian coordinates
        bv(1)=br*csphi-bphi*snphi
        bv(2)=br*snphi+bphi*csphi
        bv(3)=bz
        call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk))

        if (jk == 1) then
          call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv)
          call polellipse(qq,uu,vv,psipol0,chipol0)
          psipol0=psipol0/degree  ! convert from rad to degree
          chipol0=chipol0/degree
        end if
      else
        call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv)
        if(qq**2 < one) then
          deltapol=asin(vv/sqrt(one - qq**2))
          ext0(jk)= sqrt(half*(one + qq))
          eyt0(jk)= sqrt(half*(one - qq))*exp(-im*deltapol)
        else
          ext0(jk)= one
          eyt0(jk)= zero
        end if
      endif
    end do
  end subroutine set_pol



!   logical function inside_plasma(rrm,zzm)
!     use const_and_precisions, only : wp_, zero, one
!     use gray_params, only : iequil
!     use equilibrium, only : equian,equinum_psi,zbinf,zbsup
!     use coreprofiles, only : psdbnd
!     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 beamdata, only : dst
!     use limiter, only : rlim,zlim,nlim
!     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
! 
! ! 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



  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



  subroutine print_headers
    use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm
    implicit none
    write(uprj0,*)   ' #sst j k xt yt zt rt'
    write(uprj0+1,*) ' #sst j k xt yt zt rt'
    write(uwbm,*)    ' #sst w1 w2'
    write(udisp,*)   ' #sst Dr_Nr Di_Nr'
    write(ucenr,*)   ' #sst R z phi psin rhot ne Te Btot Nperp Npl ki '// &
         'alpha tau Pt dIds nhmax iohkw index_rt ddr'
    write(uoutr,*)   ' #i k sst x y R z psin tau Npl alpha index_rt'
    write(upec,*)    ' #rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
    write(usumm,*)   ' #Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
        'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
        'Jphimx dPdVmx drhotj drhotp'
  end subroutine print_headers



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

    write(uprfin,*)  ' #psi rhot ne Te q Jphi'
    do i=1,nq
      psin=psinr(i)
      rhot=frhotor(sqrt(psin))
      call density(psin,dens,ddens)
      call tor_curr_psi(max(eps,psin),ajphi)

      write(uprfin,"(12(1x,e12.5))") psin,rhot,dens,temp(psin),fq(psin),ajphi*1.e-6_wp_
    end do
  end subroutine print_prof



  subroutine print_bres(bres)
    use const_and_precisions, only : wp_
    use gray_params, only : iequil
    use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq
    use units, only : ubres
    implicit none
! arguments
    real(wp_) :: bres
! local constants
    integer, parameter :: icmx=2002
! local variables
    integer :: j,k,n,nconts,nctot
    integer, dimension(10) :: ncpts
    real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
    real(wp_), dimension(icmx) :: rrcb,zzcb
    real(wp_) :: rv(nq), zv(nq)
    real(wp_), dimension(nq,nq) :: btotal

    
    dr = (rmxm-rmnm)/(nq-1)
    dz = (zmxm-zmnm)/(nq-1)
    do j=1,nq
      rv(j) = rmnm + dr*(j-1)
      zv(j) = zmnm + dz*(j-1)
    end do  

!   Btotal on psi grid
    btmx=-1.0e30_wp_
    btmn=1.0e30_wp_
    do k=1,nq
      zzk=zv(k)
      do j=1,nq
        rrj=rv(j)
        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(ubres,*)'#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,nq,nq,bbb,nconts,ncpts,nctot,rrcb,zzcb)
        do j=1,nctot
          write(ubres,'(i6,12(1x,e12.5))') j,bbb,rrcb(j),zzcb(j)
        end do
      end if
      write(ubres,*)
    end do
  end subroutine print_bres



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

! build q profile on psin grid
    do i=1,nq
      qpsi(i)  = fq(psinr(i))
    end do
! locate  psi surface for q=qval
    print*
    do i=1,size(qval)
      call locate(abs(qpsi),nq,qval(i),i1)  !!!! check for non monotonous q profile
      if (i1>0.and.i1<nq) then
        call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1), &
          qval(i),psival)
        rup=rmaxis
        rlw=rmaxis
        zup=(zbsup+zmaxis)/2.0_wp_
        zlw=(zmaxis+zbinf)/2.0_wp_
        call contours_psi(psival,rcn,zcn,rup,zup,rlw,zlw)
        call print_contour(psival,rcn,zcn)
        rhot=frhotor(sqrt(psival))
        print'(4(a,f8.5))','q = ',qval(i),     '  psi = ',psival, &
                      '  rhop = ',sqrt(psival),'  rhot = ',rhot
      end if
    end do
  end subroutine print_surfq



  subroutine print_projxyzt(st,ywrk,iproj)
    use const_and_precisions, only : wp_, comp_huge, zero, one
    use beamdata, only : nray, nrayr, nrayth, rayi2jk
    use units, only : uprj0,uwbm
    implicit none
! arguments
    real(wp_), intent(in) :: st
    real(wp_), dimension(:,:), intent(in) :: ywrk
    integer, intent(in) :: iproj
! local variables
    integer :: jk,jkz,uprj
    integer, dimension(2) ::jkv
    real(wp_), dimension(3) :: xv1,dir,dxv
    real(wp_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti
! common/external functions/variables

    uprj = uprj0 + iproj

    xv1  = ywrk(1:3,1)
    dir  = ywrk(4:6,1)
    dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2)
    dir  = dir/dirm
    csth1 = dir(3)
    snth1 = sqrt(one - csth1**2)
    if(snth1 > zero) then
      csps1=dir(2)/snth1
      snps1=dir(1)/snth1
    else
      csps1=one
      snps1=zero
    end if

    if(iproj==0) then
      jkz = nray - nrayth + 1
    else
      jkz = 1
    end if

    rtimn = comp_huge
    rtimx = zero
    do jk = jkz, nray
      dxv = ywrk(1:3,jk) - xv1
      xti = dxv(1)*csps1 - dxv(2)*snps1
      yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1
      zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1
      rti = sqrt(xti**2 + yti**2)

      jkv=rayi2jk(jk)
      if(.not.(iproj==0 .and. jk==1))                  &
         write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') st,jkv,xti,yti,zti,rti
      if(iproj==1 .and. jkv(2)==nrayth) write(uprj,*)

      if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti
      if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti
    end do
    write(uprj,*)
    write(uwbm,'(3(1x,e16.8e3))') st,rtimn,rtimx
  end subroutine print_projxyzt



  subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, &
           dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
    use const_and_precisions, only : degree,zero,one
    use equilibrium, only : frhotor
    use gray_params, only : istpl0
    use beamdata, only : nray,nrayth,jkray1
    use units, only : ucenr,uoutr,udisp
    implicit none
! arguments
    integer, intent(in) :: i,jk,nhf,iokhawa,index_rt
    real(wp_), dimension(3), intent(in) :: xv
    real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim
    real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi
! local variables
    real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn
    integer :: k

    stm=st*1.0e-2_wp_
    xxm=xv(1)*1.0e-2_wp_
    yym=xv(2)*1.0e-2_wp_
    zzm=xv(3)*1.0e-2_wp_
    rrm=sqrt(xxm**2 + yym**2)

!     print central ray trajectory. dIds in A/m/W, ki in m^-1
    if(jk.eq.1) then
      phideg=atan2(yym,xxm)/degree
      if(psinv>=zero .and. psinv<=one) then
        rhot=frhotor(psinv)
      else
        rhot=1.0_wp_
      end if
      akim=anprim*ak0*1.0e2_wp_
      pt=exp(-tau)
      didsn=dids*1.0e2_wp_/qj

      write(ucenr,'(16(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, &
         psinv,rhot,dens,tekev,btot,anpr,anpl,akim,alpha,tau,pt,didsn,        &
         nhf,iokhawa,index_rt,ddr
    end if

!     print conservation of dispersion relation
    if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi

!     print outer trajectories
    if(mod(i,istpl0)==0)  then
      k = jk - jkray1 + 1
      if(k>0 .and. k<=nrayth) then
        write(uoutr,'(2i5,9(1x,e16.8e3),i5)') i,k,stm,xxm,yym,rrm,zzm, &
           psinv,tau,anpl,alpha,index_rt
      end if
    end if
  end subroutine print_output



  subroutine print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
    use const_and_precisions, only : wp_
    use units, only : upec
    implicit none
! arguments
    real(wp_), dimension(:), intent(in) :: rhop_tab,rhot_tab,jphi,jcd,dpdv, &
                                           currins,pins
    integer, intent(in) :: index_rt
! local variables
    integer :: i
    
    do i=1,size(rhop_tab)
      write(upec,'(7(1x,e16.8e3),i5)') rhop_tab(i),rhot_tab(i), &
          jphi(i),jcd(i),dpdv(i),currins(i),pins(i),index_rt
    end do  
  end subroutine print_pec



  subroutine print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav,  &
        drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
        stmx,psipol,chipol,index_rt)
    use const_and_precisions, only : wp_
    use units, only : usumm
    implicit none
    real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav,  &
        drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx,    &
        stmx,psipol,chipol
    integer, intent(in) :: index_rt
    
    write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp, &
         rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
         stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp
  end subroutine print_finals

end module graycore