! rfxqlo_toric4_util.f90 ... collection of utility routines to support 
! access to the toric4 wave code  (dmc Feb. 2005)

subroutine rfxqlo_profnt_write(ilun,ierr)

  ! write the ascii temperature and density profiles used as input by TORIC4
  ! MOD DMC March 2006 -- write on normalized Psi grid not normalized Phi grid.

  use rfxqlo_private

  implicit NONE

  integer, intent(in) :: ilun   ! i/o unit number
  integer, intent(out) :: ierr  ! exit status, 0=OK

  ! separate concentration and temperature profiles for each ion species...
  integer, parameter :: kdiff_conc = 1, kdiff_tempi = 1

  integer :: i,id_rho=0,ingrid,ingrid2,idum,ia,iz,iopen,imainsp
  integer :: id_xpsi=0
  real*8, dimension(:), allocatable :: zxgrid,zdata,zne,zxgrid2,zconc,zxpsi
  character*20 zlbl
  real*8 :: zdelta,zzeff

  real*8, parameter :: zmin_density = 1.e10    ! minimum density m^3
  real*8, parameter :: zmin_temp    = 0.005    ! minimum temperature keV
  type (rfxqlo_rfgrid) :: zgrid
  !---------------------------

  ierr=0
  if(profnt_file.eq.' ') return

  ! fetch from rfxqlo: main species info
  allocate(zconc(rfx_ions%n_ions))
  call rfxqlo_iconc(rfx_ions%n_ions,zconc,imainsp,zzeff,ierr)
  if(ierr.ne.0) return

  ! fetch from xplasma: grid for data
  call eq_ganum('__RHO',id_rho)
  call eq_ngrid(id_rho,ingrid)
  ingrid2=2*ingrid-1            ! double density grid

  allocate(zxgrid(ingrid),zxgrid2(ingrid2),zdata(ingrid2),zne(ingrid2))
  allocate(zxpsi(ingrid2))

  call rfxqlo_wavefield_id_xpsi(id_xpsi,ierr)
  if(ierr.ne.0) return

  iopen = 0

  do
     call eq_grid(id_rho,zxgrid,ingrid,idum,ierr)  ! fetch grid
     if(ierr.ne.0) exit
     zxgrid2(1)=zxgrid(1)                          ! construct doubled grid
     zdelta=0
     do i=2,ingrid
        zxgrid2(2*i-1)=zxgrid(i)
        zxgrid2(2*i-2)=(zxgrid(i)+zxgrid(i-1))/2
        zdelta=max(zdelta,(zxgrid(i)-zxgrid(i-1))/2)
     enddo

     !  get Sqrt(Psi/Psilim) x...

     call eq_rgetf(ingrid2,zxgrid2,id_xpsi,0,zxpsi,ierr)
     if(ierr.ne.0) exit
     
     open(unit=ilun,file=trim(rfx_sys%filepath)//trim(profnt_file), &
          status='unknown',iostat=ierr)
     if(ierr.ne.0) exit

     iopen = 1

     ! write header -- grid size, no. of species, A & Z of species...

     write(ilun,'(A10,4i4)') 'Rfxqlo_Pro', ingrid2, rfx_ions%n_ions, &
          imainsp, kdiff_tempi

     do i=1,rfx_ions%n_ions

        ia=rfx_ions%a_ions(i) + 0.5d0   ! convert to integer
        iz=rfx_ions%z_ions(i) + 0.5d0   ! convert to integer

        write(ilun,'(2i4)') ia,iz
     enddo

     ! write grid

     zlbl = 'x: sqrt(Psi/Psilim)'; write(ilun,'(A20)') zlbl
     write(ilun,'(5E16.9)') zxpsi

     ! fetch & write ne & Te

     call eq_rgetf(ingrid2,zxgrid2,rfx_prof%id_ne,0,zne,ierr)
     if(ierr.ne.0) exit
     call msmoo(zne,zmin_density,"Ne",1.0d-6)

     zne = zne*1.0d-6  ! convert to cm^-3

     zlbl = 'ne [cm^-3]'; write(ilun,'(A20)') zlbl
     write(ilun,'(5E16.9)') zne

     call eq_rgetf(ingrid2,zxgrid2,rfx_prof%id_te,0,zdata,ierr)
     if(ierr.ne.0) exit
     call msmoo(zdata,zmin_temp,"Te",1.0d0)

     zlbl = 'Te [keV]'; write(ilun,'(A20)') zlbl
     write(ilun,'(5E16.9)') zdata

     ! fetch & write ion densities; skip main species...

     do i=1,rfx_ions%n_ions
        if(i.eq.imainsp) cycle

        zlbl = "Ni_"//trim(rfx_ions%name_ions(i))//' [cm^-3]'
        write(ilun,'(A20)') zlbl

        if(zconc(i).le.0.0d0) then
           zdata = zmin_density

        else
           call eq_rgetf(ingrid2,zxgrid2,rfx_ions%id_ni(i),0,zdata,ierr)
           if(ierr.ne.0) exit

           call msmoo(zdata,zmin_density,zlbl,1.0d-6)
        endif

        zdata = zdata*1.0d-6
        write(ilun,'(5E16.9)') zdata
     enddo
     if(ierr.ne.0) exit

     ! fetch & write ion temperatures

     do i=1,rfx_ions%n_ions
        if(zconc(i).le.0.0d0) then
           zdata = zmin_temp
        else
           call eq_rgetf(ingrid2,zxgrid2,rfx_ions%id_ti(i),0,zdata,ierr)
           if(ierr.ne.0) exit
           call msmoo(zdata,zmin_temp,"Ti "//rfx_ions%name_ions(i),1.0d0)
        endif
           zlbl=' '
           write(zlbl,'("TEMPI(",i1,") ")') i; write(ilun,'(A20)') zlbl
        write(ilun,'(5E16.9)') zdata
     enddo
     if(ierr.ne.0) exit

     exit

  enddo

  if(iopen.eq.1) then
     close(unit=ilun)
     call rfxqlo_rfgrid_get(zgrid,ierr)
     if(ierr.ne.0) then
!        write(lunmsg,*) ' ?rfxqlo_profnt_write: rfxqlo_rfgrid_get failure.'
        return
     endif
     ! 03/19/2013 CLF: Mark if serial
     if (zgrid%ncpu .gt. 1) then
        call trmpi_requestor('BROADCAST_FILE', &
          trim(rfx_sys%filepath)//trim(profnt_file))
!     else
!        call trmpi_requestor('!BROADCAST_FILE', &
!          trim(rfx_sys%filepath)//trim(profnt_file))
     endif
  endif

  deallocate(zxgrid,zxgrid2,zdata,zne,zxpsi,zconc)

  contains

    subroutine msmoo(zdata, zmin_data, zname_data, znorm_data)
      real*8, intent(inout)     :: zdata(ingrid2)
      real*8, intent(in)        :: zmin_data      ! force data above this minimum
      character*(*), intent(in) :: zname_data     ! data name for warning message
      real*8, intent(in)        :: znorm_data     ! normalization of data for messages

      real*8, parameter :: zbc0=1.0d0  ! max smooth -> d(zdata)/dx near axis
      real*8, parameter :: zbc1=0.0d0  ! min smooth at edge -- fixed point.

      integer :: ier
      real*8  :: zmin_got   ! minimum of the data

      zmin_got = minval(zdata)
      if (zmin_got<zmin_data) then
         print '(a,es12.3,a,es12.3)', '!rfxqlo_profnt_write: minimum '//trim(zname_data) &
              //' toric input data is ',znorm_data*zmin_got,' -- reset to ',znorm_data*zmin_data
         zdata = max(zdata,zmin_data)
      end if
      call r8_qksmooth(ingrid2,zxgrid2,zdata,zdelta,zbc0,zbc1,ier)
      if(ier.ne.0) then
         call errmsg_exit(' ? fatal error in msmoo, rfxqlo_profnt_write')
      endif
    end subroutine msmoo

end subroutine rfxqlo_profnt_write

subroutine rfxqlo_xx_write(ilun,ierr)
  ! DMC March 2006 -- write file containing mapping btw xpsi and xphi
  !   sqrt(Psi/Psilim), sqrt(Phi/Philim)
  ! This will define the wavefield output grid for t4_mod_tofpp in TORIC5

  use rfxqlo_private
  use rfxqlo_interface, only: rfxqlo_wavefield_nx,rfxqlo_wavefield_xb, &
       rfxqlo_wavefield_xpsi,rfxqlo_rfgrid_get

  implicit NONE

  integer, intent(in) :: ilun   ! i/o unit number
  integer, intent(out) :: ierr  ! exit status, 0=OK

  !------------------

  integer :: inx,inxb
  real*8, dimension(:), allocatable :: xb,xpsi
  character*20 :: zlbl
  type (rfxqlo_rfgrid) :: zgrid

  !------------------

  call rfxqlo_wavefield_nx(inx,ierr)  ! #zones
  if(ierr.ne.0) return

  open(unit=ilun,file=trim(rfx_sys%filepath)//trim(xx_file), &
          status='unknown',iostat=ierr)
  if(ierr.ne.0) return

  ! OK have grid size & open file...

  inxb=inx+1                          ! #bdys
  allocate(xb(inxb),xpsi(inxb))     

  do
     call rfxqlo_wavefield_xb(xb,ierr)
     if(ierr.ne.0) exit

     call rfxqlo_wavefield_xpsi(xpsi,ierr)
     if(ierr.ne.0) exit

     write(ilun,'(A10,i4)') 'Rfxqlo_XX ',inxb

     zlbl = 'x: sqrt(Psi/Psilim)'; write(ilun,'(A20)') zlbl
     write(ilun,'(5E16.9)') xpsi

     zlbl = 'x: sqrt(Phi/Philim)'; write(ilun,'(A20)') zlbl
     write(ilun,'(5E16.9)') xb

     deallocate(xb,xpsi)

     exit
  enddo

  close(unit=ilun)
  ! 03/19/2013 CLF: Mark if serial
  call rfxqlo_rfgrid_get(zgrid,ierr)
  if(ierr.ne.0) then
!     write(lunmsg,*) ' ? rfxqlo_xx_write: rfxqlo_rfgrid_get failure.'
     return
  endif
  if (zgrid%ncpu .gt. 1) then
     call trmpi_requestor('BROADCAST_FILE', &
       trim(rfx_sys%filepath)//trim(xx_file))
!  else
!     call trmpi_requestor('!BROADCAST_FILE', &
!       trim(rfx_sys%filepath)//trim(xx_file))
  endif
end subroutine rfxqlo_xx_write

subroutine rfxqlo_equigs_write(ilun,ierr)

  ! write the ascii equilibrium data file used as input by TORIC4
  ! mod dmc Sept 2005 -- require at least 80 radial grid points
  ! MOD DMC March 2006 -- write on normalized Psi grid not normalized Phi grid.

  use rfxqlo_private
  implicit NONE

  integer, intent(in) :: ilun   ! i/o unit number
  integer, intent(out) :: ierr  ! exit status, 0=OK

  !---------------------------
  !  DMC Apr 2007: actual write moved to old_xplasma/equigs_write.f90
  !   where it can be called as a "geqxpl" program option as well
  !---------------------------
  integer :: id_q,id_xpsi
  type (rfxqlo_rfgrid) :: zgrid
  !---------------------------

  ! also get "Q" fcn id (for q profile)

  ierr = 1
  call eq_gfnum('Q',id_q)
  if(id_q.eq.0) return

  ! and get xpsi id
  call rfxqlo_wavefield_id_xpsi(id_xpsi,ierr)
  if(id_xpsi.eq.0) return

  call equigs_write(id_xpsi, id_q, ilun, &
       trim(rfx_sys%filepath)//trim(equigs_file), &
       ierr)
  ! 03/19/2013 CLF: Mark if serial
  call rfxqlo_rfgrid_get(zgrid,ierr)
  if(ierr.ne.0) then
!     write(lunmsg,*) ' ?rfxqlo_profnt_write: rfxqlo_rfgrid_get failure.'
     return
  endif
  if (zgrid%ncpu .gt. 1) then
     call trmpi_requestor('BROADCAST_FILE', &
       trim(rfx_sys%filepath)//trim(equigs_file))
!  else
!     call trmpi_requestor('!BROADCAST_FILE', &
!       trim(rfx_sys%filepath)//trim(equigs_file))
  endif
end subroutine rfxqlo_equigs_write

subroutine rfxqlo_genppj
  use rfxqlo_private
  implicit NONE

  integer, parameter :: inpp = 51
  real*8 :: zppj(inpp),zfppj(inpp)
  real*8, parameter :: zminppj=0.2d0, zmaxppj=5.0d0
  integer :: i,k
  real*8 :: zsum25,zsumtot,ztarg25,zj,zi,zip
  real*8, parameter :: zxtarg25=0.25d0

  !-------------------------
  !  find target ratio:  Iphi(0.25)/Iphi(1.0)

  zsumtot = Iphi_qlo(nx_qlo)
  do k=1,nx_qlo
     if((xb_qlo(k-1).le.zxtarg25).and.(zxtarg25.le.xb_qlo(k))) then
        ztarg25=Iphi_qlo(k-1)+(Iphi_qlo(k)-Iphi_qlo(k-1))* &
             (zxtarg25-xb_qlo(k-1))/(xb_qlo(k)-xb_qlo(k-1))
        ztarg25 = ztarg25/zsumtot
        exit
     endif
  enddo

  !-------------------------
  !  form grid of j exponent values for j = j0*(1-x**ppjti)**ppjte expression

  ppjti = 2.0d0

  do i=1,inpp
     zppj(i)=(zminppj*(inpp-i)+zmaxppj*(i-1))/(inpp-1)

  !  j = j0*(1-x**ppjti)**ppjte
  !      ppjti=2.0 forced; ppjte set to get right averaged current in
  !      central region [0,zxtarg25]

     zi=0.0;
     do k=1,nx_qlo
        zip=zi
        zj=(1-x_qlo(k)**ppjti)**zppj(i)
        zi = zi + zj*da_qlo(k)
        if((xb_qlo(k-1).le.zxtarg25).and.(zxtarg25.le.xb_qlo(k))) then
           zsum25=zip+(zi-zip)* &
             (zxtarg25-xb_qlo(k-1))/(xb_qlo(k)-xb_qlo(k-1))
        endif
     enddo
     zsumtot=zi
     zfppj(i)=zsum25/zsumtot
  enddo

  !  find match

  if(ztarg25.le.zfppj(1)) then
     ppjte=zfppj(1)

  else if(ztarg25.ge.zfppj(inpp)) then
     ppjte=zfppj(inpp)

  else
     do i=2,inpp
        if((zfppj(i-1).le.ztarg25).and.(ztarg25.le.zfppj(i))) then
           zi=(ztarg25-zfppj(i-1))/(zfppj(i)-zfppj(i-1))
           ppjte = zppj(i-1)*(1-zi) + zppj(i)*zi
           exit
        endif
     enddo
  endif

end subroutine rfxqlo_genppj

subroutine rfxqlo_toric4_files(z_profnt_file,z_equigs_file,z_xx_file)
  ! return filenames of TORIC4 input files

  ! see rfxqlo_input_verify (rfxqlo_verify.f90)
  ! and rfxqlo_profnt_write (above)
  ! and rfxqlo_equigs_write (above)

  use rfxqlo_private
  implicit NONE
  character*(*), intent(out) :: z_profnt_file
  character*(*), intent(out) :: z_equigs_file
  character*(*), intent(out) :: z_xx_file

  z_profnt_file = profnt_file
  z_equigs_file = equigs_file
  z_xx_file = xx_file

end subroutine rfxqlo_toric4_files

subroutine rfxqlo_toric4_ppjs(zicurr,z_ppjti,z_ppjte)
  ! return fit information
  !   j = J0*(1-x**ppjti)**ppjte
  !   with J0 chosen to yield the right total current
  !
  ! this information calculated by rfxqlo_genppj (above)

  use rfxqlo_private
  implicit NONE
  real*8, intent(out) :: zicurr  ! total plasma current (Amps)
  real*8, intent(out) :: z_ppjti ! internal exponent
  real*8, intent(out) :: z_ppjte ! external exponent

  zicurr = Iphi_qlo(nx_qlo)
  z_ppjte = ppjte
  z_ppjti = ppjti

end subroutine rfxqlo_toric4_ppjs

subroutine rfxqlo_toric4_fpread(ilun,zfile,ksym,ierr)
  !
  !  read TORIC4 wave field results & power and current profile information.
  !
  use rfxqlo_interface
  implicit NONE

  integer, intent(in) :: ilun   ! channel for file i/o
  character*(*), intent(in) :: zfile  ! filename

  integer, intent(in) :: ksym   ! n_phi symmetry flag =2 => symmetric
  ! if ksym = 1, spectrum is asymmetric and current calculation is for a
  ! single signed value of n_phi; if ksym = 2 the current calculation is 
  ! for +/- the current n_phi value, and the power has to be doubled.
  ! and also the driven current set to zero

  ! NB dmc Aug 2005 -- per advice of TORIC group, nphi symmetrization 
  ! is applied for the RF QLO regardless of whether the initial spectrum
  ! is symmetric or not.  This is known to be necessary for the "numerical
  ! survival" of fppmod.  In anticipation of this +/- symmetrization of
  ! the QLO contribution for *every* nphi value, the power seen by TORIC
  ! was halved; so, the net power and current drive have to be DOUBLED here.
  ! hence the mysterious factor of TWO below...

  integer, intent(out) :: ierr  ! status code on exit (0=OK)

  type (rfxqlo_ion) :: zions
  type (rfxqlo_outpro) :: zpro

  integer :: inxi,inti,inspeci, idum
  integer :: inttp1,i,j,k,id_xpx=0,id_rho_orig=0
  integer :: inx, intt, inspec
  integer :: istat,ifound,ieq,ic1,ic2,ic
  character*20 zfdcod
  real*8 :: zpwc

  real*8, parameter :: ZERO = 0.0d0
  real*8, parameter :: TWO = 2.0d0

  character*80 ztext
  real*8 :: ztmp,zdum,zdum1,zdum2,ztol

  !-------------------------------------------------
  real*8, dimension(:), allocatable :: zbufr,zjbuf,zbufn,zxbuf,zxtr
  real*8, dimension(:,:), allocatable :: zkperp,zeplus,zpbuf
  complex*16, dimension(:), allocatable :: zbufc1,zbufc2
  complex*16, dimension(:,:), allocatable :: zcrfpol
  !-------------------------------------------------

  call rfxqlo_ion_get(zions,ierr)
  if(ierr.ne.0) return

  call rfxqlo_outpro_init(zpro,ierr)
  if(ierr.ne.0) return

  inspec = zions%n_ions
  call rfxqlo_wavefield_dims(inx,inttp1,ierr)
  if(ierr.ne.0) return

  intt = inttp1 - 1

  open(unit=ilun,file=trim(zfile),status='old',iostat=ierr)
  if(ierr.ne.0) return

  allocate(zbufr(intt),zbufc1(intt),zbufc2(intt),zbufn(inspec))
  allocate(zkperp(inx,inttp1),zeplus(inx,inttp1),zcrfpol(inx,inttp1))
  allocate(zxbuf(0:inx),zjbuf(0:inx),zpbuf(0:inx,inspec+2),zxtr(0:inx))
  call rfxqlo_wavefield_xb(zxtr,ierr)

  do

     read(ilun,'(3I5)') inxi,inti,inspeci
     if(inxi.ne.inx) then
        ierr=101
     else if(inti.ne.intt) then
        ierr=102
     else if(inspeci.ne.inspec) then
        ierr=103
     endif
     if(ierr.ne.0) exit

     read(ilun,'(A)') ztext
     read(ilun,'(6E15.6)') zbufr

     !  read data for each zone

     zjbuf(0)=0; zpbuf(0,:)=0; zxbuf(0)=0

     do k=1,inx

        zjbuf(k)=zjbuf(k-1)
        zpbuf(k,:)=zpbuf(k-1,:)

        read(ilun,'(A)',iostat=istat) ztext
        if(istat.ne.0) then
           ierr=999  ! short fppdata file -- indicates TORIC error
           exit
        endif

        read(ilun,'(3E15.6)') zdum1,zdum2,ztmp
        zxbuf(k)=ztmp  ! psi of outer bdy of zone

        read(ilun,'(A)') ztext
        read(ilun,'(6E15.6)') zbufr   ! k_perp^2  this zone,(see t4_mod_tofppmod)
        read(ilun,'(A)') ztext
        read(ilun,'(6E15.6)') zbufc1  ! E_+ this zone (complex)
        read(ilun,'(A)') ztext
        read(ilun,'(6E15.6)') zbufc2  ! E_- this zone (complex)
        read(ilun,'(A)') ztext
        read(ilun,'(A)') ztext  ! skip volume
        read(ilun,'(A)') ztext
        read(ilun,'(E15.6)') ztmp
        zpbuf(k,inspec+1)=zpbuf(k,inspec+1) + ztmp  ! electron heating

        read(ilun,'(A)') ztext
        read(ilun,'(E15.6)') zbufn
        zpbuf(k,1:inspec)=zpbuf(k,1:inspec)+zbufn   ! fundamental resonance htg

        read(ilun,'(A)') ztext
        read(ilun,'(E15.6)') zbufn
        zpbuf(k,1:inspec)=zpbuf(k,1:inspec)+zbufn   ! higher harmonic heating

        read(ilun,'(A)') ztext
        read(ilun,'(A)') ztext  ! skip area
        read(ilun,'(A)') ztext
        read(ilun,'(E15.6)') ztmp
        zjbuf(k)=zjbuf(k)+ztmp  ! driven current

        zkperp(k,1:intt)=sqrt(max(ZERO,zbufr)); zkperp(k,inttp1)=zkperp(k,1)! k_perp, m-1
        zeplus(k,1:intt)=abs(zbufc1); zeplus(k,inttp1)=zeplus(k,1)
        zcrfpol(k,1:intt)=zbufc2/zbufc1; zcrfpol(k,inttp1)=zcrfpol(k,1)
     enddo

     call rfxqlo_wavefield_set1(zkperp,zeplus,zcrfpol,ierr)
     if(ierr.ne.0) return

     ! verify that x axis is available; if not, create it...
     call eq_ganum('x_toric4',id_xpx)
     if(id_xpx.eq.0) then
        call eq_ganum('__RHO',id_rho_orig)
        ztol=zxtr(inx)/100000
        call eqm_uaxis('x_toric4',id_rho_orig,0,zxtr,inx+1,ztol,id_xpx,ierr)
        if(ierr.ne.0) exit
     endif

     idum=0
     zdum=0
     do i=1,inspec+2
        !  store the powers in xplasma -- apply ksym now --
        !  dmc -- no, apply TWO, regardless of ksym value -- Aug 2005
        !  so:
        !  call eqm_rhofun(100,id_xpx,zpro%xn_pwrs(i),ksym*zpbuf(0:inx,i), &
        !  --> replaced:

        call eqm_rhofun(100,id_xpx,zpro%xn_pwrs(i),TWO*zpbuf(0:inx,i), &
             idum,zdum,idum,zdum, zpro%id_pwrs(i), ierr)
        if(ierr.ne.0) exit
     enddo
     if(ierr.ne.0) exit

     !  ... and the driven current

     !  dmc -- apply TWO here also

     call eqm_rhofun(100,id_xpx,zpro%xn_jrfe,TWO*zjbuf(0:inx), &
          zdum,zdum,idum,zdum, zpro%id_jrfe, ierr)
     if(ierr.ne.0) exit

     call rfxqlo_outpro_store(zpro,ierr)
     if(ierr.ne.0) exit

     exit
  enddo

  !-----------------------------
  ! DMC Jan 2009 -- look for mode coupling data

  ifound = 0
  do
     ztext = ' '
     read(unit=ilun,fmt='(A)',iostat=istat) ztext
     if(istat.ne.0) exit

     if(index(ztext,'__Supplemental_Scalars__').gt.0) then
        ifound = 1
        exit
     endif
  enddo

  if(ifound.eq.1) then
     do
        ztext = ' '
        read(unit=ilun,fmt='(A)') ztext
        if(index(ztext,'tpw').gt.0) exit
     enddo

     ieq=index(ztext,'=')
     if(ieq.le.0) then
        call errmsg_exit('rfxqlo_toric4_util parse error 1: '// &
             trim(ztext))
     endif

     ic1=0
     ic2=0
     do ic=ieq+1,len(ztext)
        if(ic1.eq.0) then
           if(ztext(ic:ic).ne.' ') then
              ic1=ic
           endif
        else
           if((ztext(ic:ic).eq.' ').or.(ztext(ic:ic).eq.char(9))) then
              ic2=ic-1
              exit
           endif
        endif
     enddo

     if(ic2.le.0) then
        call errmsg_exit('rfxqlo_toric4_util parse error 2: '// &
             trim(ztext))
     endif

     ic=ic2-ic1+1  ! non-blank length of numeric field
     if(ic.gt.20) then
        call errmsg_exit('rfxqlo_toric4_util parse error 3: '// &
             trim(ztext))
     endif

     zfdcod=' '
     zfdcod(20-ic+1:20) = ztext(ic1:ic2)

     read(zfdcod,'(G20.0)') zpwc

     ! OK found the coupling now set it in rfxqlo

     call rfxqlo_set_coupling(zpwc,ierr)

  endif

  !-----------------------------
  close(unit=ilun)

  deallocate(zbufr,zbufc1,zbufc2,zbufn)
  deallocate(zkperp,zeplus,zcrfpol)
  deallocate(zjbuf,zpbuf)

end subroutine rfxqlo_toric4_fpread
