module equilibrium
  use const_and_precisions, only : wp_
  implicit none
  
  real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi
  real(wp_), save :: btrcen ! used only for jcd_astra def.
  real(wp_), save :: rcen   ! computed as fpol(a)/btrcen
  real(wp_), save :: rmnm,rmxm,zmnm,zmxm
  real(wp_), save :: zbinf,zbsup
!  real(wp_), save :: rup,zup,rlw,zlw

  integer, parameter :: kspl=3,ksplp=kspl+1

! === 2d spline psi(r,z), normalization and derivatives ==========
  integer, save :: nsr, nsz
  real(wp_), save :: psia, psiant, psinop, psnbnd
  real(wp_), dimension(:), allocatable, save :: tr,tz
  real(wp_), dimension(:), allocatable, save :: cceq,   cceq01, cceq10, &
                                                cceq20, cceq02, cceq11

! === 1d spline fpol(psi) ========================================
! integer, save :: npsiest
  integer, save :: nsf
  real(wp_), dimension(:), allocatable, save :: tfp, cfp
  real(wp_), save :: fpolas

! === 1d spline rhot(rhop), rhop(rhot), q(psi) ===================
! computed on psinr,rhopnr [,rhotnr] arrays
  integer, save :: nq,nrho
  real(wp_), dimension(:), allocatable, save :: psinr,rhopr,rhotr
  real(wp_), dimension(:,:), allocatable, save :: cq,crhop,crhot
  real(wp_), save :: phitedge,aminor
  real(wp_), save :: q0,qa,alq

contains

  subroutine read_eqdsk(params, data, unit)
    ! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm).
    ! If given, the file is opened in the `unit` number.
    ! For a description of the G-EQDSK, see the GRAY user manual.
    use const_and_precisions, only : one
    use gray_params,          only : equilibrium_parameters, equilibrium_data
    use utils,                only : get_free_unit
    use logger,               only : log_error

    implicit none

    ! subroutine arguments
    type(equilibrium_parameters), intent(in)  :: params
    type(equilibrium_data),       intent(out) :: data
    integer, optional,            intent(in)  :: unit

    ! local variables
    integer :: u, idum, i, j, nr, nz, nbnd, nlim
    character(len=48) :: string
    real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis
    real(wp_) :: xdum ! dummy variable, used to discard data
    integer :: err

    u = get_free_unit(unit)

    ! Open the G-EQDSK file
    open(file=params%filenm, status='old', action='read', unit=u, iostat=err)
    if (err /= 0) then
      call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', &
                     mod='equilibrium', proc='read_eqdsk')
      call exit(1)
    end if

    ! get size of main arrays and allocate them
    if (params%idesc == 1) then
      read (u,'(a48,3i4)') string,idum,nr,nz
    else
      read (u,*) nr, nz
    end if
    if (allocated(data%rv))    deallocate(data%rv)
    if (allocated(data%zv))    deallocate(data%zv)
    if (allocated(data%psin))  deallocate(data%psin)
    if (allocated(data%psinr)) deallocate(data%psinr)
    if (allocated(data%fpol))  deallocate(data%fpol)
    if (allocated(data%qpsi))  deallocate(data%qpsi)
    allocate(data%rv(nr), data%zv(nz), &
             data%psin(nr, nz),        &
             data%psinr(nr),           &
             data%fpol(nr),            &
             data%qpsi(nr))

    ! Store 0D data and main arrays
    if (params%ifreefmt==1) then
      read (u, *) dr, dz, data%rvac, rleft, zmid
      read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum
      read (u, *) xdum, xdum, xdum, xdum, xdum
      read (u, *) xdum, xdum, xdum, xdum, xdum
      read (u, *) (data%fpol(i), i=1,nr)
      read (u, *) (xdum,i=1, nr)
      read (u, *) (xdum,i=1, nr)
      read (u, *) (xdum,i=1, nr)
      read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz)
      read (u, *) (data%qpsi(i), i=1,nr)
    else
      read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid
      read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum
      read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
      read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
      read (u, '(5e16.9)') (data%fpol(i),i=1,nr)
      read (u, '(5e16.9)') (xdum,i=1,nr)
      read (u, '(5e16.9)') (xdum,i=1,nr)
      read (u, '(5e16.9)') (xdum,i=1,nr)
      read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz)
      read (u, '(5e16.9)') (data%qpsi(i),i=1,nr)
    end if

    ! Get size of boundary and limiter arrays and allocate them
    read (u,*) nbnd, nlim
    if (allocated(data%rbnd)) deallocate(data%rbnd)
    if (allocated(data%zbnd)) deallocate(data%zbnd)
    if (allocated(data%rlim)) deallocate(data%rlim)
    if (allocated(data%zlim)) deallocate(data%zlim)

    ! Load plasma boundary data
    if(nbnd > 0) then
      allocate(data%rbnd(nbnd), data%zbnd(nbnd))
      if (params%ifreefmt == 1) then
        read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd)
      else
        read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd)
      end if
    end if

    ! Load limiter data
    if(nlim > 0) then
      allocate(data%rlim(nlim), data%zlim(nlim))
      if (params%ifreefmt == 1) then
        read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim)
      else
        read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim)
      end if
    end if

    ! End of G-EQDSK file
    close(u)
    
    ! Build rv,zv,psinr arrays
    zleft = zmid-0.5_wp_*dz
    dr  = dr/(nr-1)
    dz  = dz/(nz-1)
    dps = one/(nr-1)
    do i=1,nr
      data%psinr(i) = (i-1)*dps
      data%rv(i)    = rleft + (i-1)*dr
    end do
    do i=1,nz
      data%zv(i) = zleft + (i-1)*dz
    end do

    ! Normalize psin
    data%psia = psiedge - psiaxis
    if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia

  end subroutine read_eqdsk


  subroutine read_equil_an(filenm, ipass, data, unit)
    ! Reads the MHD equilibrium `data` in the analytical format
    ! from params%filenm.
    ! If given, the file is opened in the `unit` number.
    !
    ! TODO: add format description
    use gray_params, only : equilibrium_data
    use utils,       only : get_free_unit
    use logger,      only : log_error

    implicit none

    ! subroutine arguments
    character(len=*),       intent(in)  :: filenm
    integer,                intent(in)  :: ipass
    type(equilibrium_data), intent(out) :: data
    integer, optional,      intent(in)  :: unit

    ! local variables
    integer :: i, u, nlim
    integer :: err
    real(wp_) :: rr0m, zr0m, rpam, b0, q0, qa, alq

    u = get_free_unit(unit)

    open(file=filenm, status='old', action='read', unit=u, iostat=err)
    if (err /= 0) then
      call log_error('opening equilibrium file ('//trim(filenm)//') failed!', &
                     mod='equilibrium', proc='read_equil_an')
      call exit(1)
    end if

    read(u, *) rr0m, zr0m, rpam
    read(u, *) b0
    read(u, *) q0, qa, alq

    if(allocated(data%rv))   deallocate(data%rv)
    if(allocated(data%zv))   deallocate(data%zv)
    if(allocated(data%fpol)) deallocate(data%fpol)
    if(allocated(data%qpsi)) deallocate(data%qpsi)
    allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3))

    data%rv   = [rr0m, rpam]
    data%zv   = [zr0m]
    data%fpol = [b0*rr0m]
    data%qpsi = [q0, qa, alq]

    if(ipass >= 2) then
      ! get size of boundary and limiter arrays and allocate them
      read (u,*) nlim
      if (allocated(data%rlim)) deallocate(data%rlim)
      if (allocated(data%zlim)) deallocate(data%zlim)

      ! store boundary and limiter data
      if(nlim > 0) then
        allocate(data%rlim(nlim), data%zlim(nlim))
        read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
      end if
    end if

    close(u)
  end subroutine read_equil_an


  subroutine change_cocos(data, cocosin, cocosout, error)
    ! Convert the MHD equilibrium data from one coordinate convention
    ! (COCOS) to another. These are specified by `cocosin` and
    ! `cocosout`, respectively.
    !
    ! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010
    use const_and_precisions, only : zero, one, pi
    use gray_params,          only : equilibrium_data

    implicit none

    ! subroutine arguments
    type(equilibrium_data), intent(inout) :: data
    integer, intent(in) :: cocosin, cocosout
    integer, intent(out), optional :: error

    ! local variables
    real(wp_) :: isign, bsign
    integer   :: exp2pi, exp2piout
    logical   :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout

    call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos)
    call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout)

    ! Check sign consistency
    isign = sign(one, data%psia)
    if (.not.psiincr) isign = -isign
    bsign = sign(one, data%fpol(size(data%fpol)))
    if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then
      ! Warning: sign inconsistency found among q, Ipla and Bref
      data%qpsi = -data%qpsi
      if (present(error)) error = 1
    else
      if (present(error)) error = 0
    end if

    ! Convert cocosin to cocosout

    ! Opposite direction of toroidal angle phi in cocosin and cocosout
    if (phiccw .neqv. phiccwout) data%fpol = -data%fpol

    ! q has opposite sign for given sign of Bphi*Ip
    if (qpos .neqv. qposout) data%qpsi = -data%qpsi

    ! psi and Ip signs don't change accordingly
    if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
      data%psia = -data%psia

    ! Convert Wb to Wb/rad or viceversa
    data%psia = data%psia * (2.0_wp_*pi)**(exp2piout - exp2pi)
  end subroutine change_cocos


  subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos)
    implicit none

    ! subroutine arguments
    integer, intent(in)  :: cocos
    integer, intent(out) :: exp2pi
    logical, intent(out) :: phiccw, psiincr, qpos

    ! local variables
    integer :: cmod10, cmod4

    cmod10 = mod(cocos, 10)
    cmod4  = mod(cmod10, 4)

    ! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad
    exp2pi = (cocos - cmod10)/10

    ! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW
    phiccw = (mod(cmod10, 2)== 1)

    ! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip
    psiincr = (cmod4==1 .or. cmod4==2)

    ! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip
    qpos = (cmod10<3 .or. cmod10>6)
  end subroutine decode_cocos


  subroutine eq_scal(params, data)
    ! Rescale the magnetic field (B) and the plasma current (I)
    ! and/or force their signs.
    !
    ! Notes:
    ! 1. signi and signb are ignored on input if equal to 0.
    !    They are used to assign the direction of Bphi and Ipla BEFORE scaling.
    ! 2. cocos=3 assumed: CCW direction is >0
    ! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged
    ! 4. factb<0 reverses the directions of Bphi and Ipla

    use const_and_precisions, only : one
    use gray_params,          only : equilibrium_parameters, equilibrium_data

    implicit none

    ! subroutine arguments
    type(equilibrium_parameters), intent(inout) :: params
    type(equilibrium_data),       intent(inout) :: data

    ! local variables
    real(wp_) :: last_fpol

    last_fpol = data%fpol(size(data%fpol))
    
    if (params%sgni /=0) &
      data%psia = sign(data%psia, real(-params%sgni, wp_))
    if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) &
      data%fpol = -data%fpol

    data%psia = data%psia * params%factb
    data%fpol = data%fpol * params%factb
    params%sgni = int(sign(one, -data%psia))
    params%sgnb = int(sign(one, last_fpol))
  end subroutine eq_scal


  subroutine set_equil_spline(params, data)
    ! Computes splines for the MHD equilibrium data and stores them
    ! in their respective global variables, see the top of this file.
    use const_and_precisions, only : zero, one
    use gray_params,          only : equilibrium_parameters, equilibrium_data
    use dierckx,              only : regrid, coeff_parder, curfit, splev
    use gray_params,          only : iequil
    use reflections,          only : inside
    use utils,                only : vmaxmin, vmaxmini
    use logger,               only : log_info

    implicit none

    ! subroutine arguments
    type(equilibrium_parameters), intent(in) :: params
    type(equilibrium_data),       intent(in) :: data

    ! local constants
    integer, parameter :: iopt=0

    ! local variables
    integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf
    integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup
    real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp
    real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1
    real(wp_), dimension(size(data%psinr)) :: rhotn
    real(wp_), dimension(1) :: fpoli
    real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk
    integer, dimension(:), allocatable :: iwrk
    integer :: ier,ixploc,info,i,j,ij
    character(256) :: msg ! for log messages formatting

    ! compute array sizes and prepare working space arrays
    nr=size(data%rv)
    nz=size(data%zv)
    nrz=nr*nz
    nrest=nr+ksplp
    nzest=nz+ksplp
    lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest)
    liwrk=nz+nr+nrest+nzest+kspl

    npsi=size(data%psinr)
    npsest=npsi+ksplp
    lwrkf=npsi*ksplp+npsest*(7+3*kspl)

    allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest)))

    ! spline fitting/interpolation of data%psin(i,j) and derivatives

    ! allocate knots and spline coefficients arrays
    if (allocated(tr)) deallocate(tr)
    if (allocated(tz)) deallocate(tz)
    if (allocated(cceq)) deallocate(cceq)
    allocate(tr(nrest),tz(nzest),cceq(nrz))

    ! length in m !!!
    rmnm=data%rv(1)
    rmxm=data%rv(nr)
    zmnm=data%zv(1)
    zmxm=data%zv(nz)

    if (iequil>2) then
      ! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO
      ! presence of boundary anticipated here to filter invalid data
      nbnd = min(size(data%rbnd), size(data%zbnd))

      ! determine number of valid grid points
      nrz=0
      do j=1,nz
        do i=1,nr
          if (nbnd.gt.0) then
            if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
          else
            if(data%psin(i,j).le.0.0d0) cycle
          end if
          nrz=nrz+1
        end do
      end do

      ! store valid data
      allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz))
      ij=0
      do j=1,nz
        do i=1,nr
          if (nbnd.gt.0) then
            if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
          else
            if(data%psin(i,j).le.0.0d0) cycle
          end if
          ij=ij+1
          rv1d(ij)=data%rv(i)
          zv1d(ij)=data%zv(j)
          fvpsi(ij)=data%psin(i,j)
          wf(ij)=1.0d0
        end do
      end do

      ! fit as a scattered set of points
      ! use reduced number of knots to limit memory comsumption ?
      nsr=nr/4+4
      nsz=nz/4+4
      sspln=params%ssplps
      call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln,               &
                      rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier)
      ! if ier=-1 data are fitted using params%ssplps=0    
      if(ier.eq.-1) then
        sspln=0.0_wp_
        nsr=nr/4+4
        nsz=nz/4+4
        call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln,             &
                        rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier)
      end if
      deallocate(rv1d,zv1d,wf,fvpsi)
      ! reset nrz to the total number of grid points for next allocations
      nrz=nr*nz
    else
      ! iequil==2: data are valid on the full R,z grid

      ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs
      allocate(fvpsi(nrz))
      fvpsi=reshape(transpose(data%psin),(/nrz/))
      sspln=params%ssplps
      call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm,         &
                  kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp,  &
                  wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
      ! if ier=-1 data are re-fitted using params%ssplps=0    
      if(ier==-1) then
        sspln=0.0_wp_
        call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm,        &
                    kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, &
                    wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
      end if
      deallocate(fvpsi)
    end if

    ! compute spline coefficients for psi partial derivatives
    lw10 = nr*(ksplp-1) + nz*ksplp     + nrz
    lw01 = nr*ksplp     + nz*(ksplp-1) + nrz
    lw20 = nr*(ksplp-2) + nz*ksplp     + nrz
    lw02 = nr*ksplp     + nz*(ksplp-2) + nrz
    lw11 = nr*(ksplp-1) + nz*(ksplp-1) + nrz
    if (allocated(cceq10)) deallocate(cceq10)
    if (allocated(cceq01)) deallocate(cceq01)
    if (allocated(cceq20)) deallocate(cceq20)
    if (allocated(cceq02)) deallocate(cceq02)
    if (allocated(cceq11)) deallocate(cceq11)
    allocate(cceq10(lw10),cceq01(lw01),cceq20(lw20),cceq02(lw02),cceq11(lw11))
    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,0,cceq10,lw10,ier)
    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,1,cceq01,lw01,ier)
    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier)
    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier)
    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier)

    ! Spline interpolation of data%fpol(i)

    ! Allocate knots and spline coefficients arrays
    if (allocated(tfp)) deallocate(tfp)
    if (allocated(cfp)) deallocate(cfp)
    allocate(tfp(npsest),cfp(npsest))
    allocate(wf(npsi))
    wf(1:npsi-1)=one
    wf(npsi)=1.0e2_wp_
    call curfit(iopt,npsi,data%psinr,data%fpol,wf,zero,one,kspl,params%ssplf,npsest,nsf, &
                tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier)
    call splev(tfp,nsf,cfp,kspl,data%psinr(npsi:npsi),fpoli,1,ier)
    ! Set vacuum value used outside 0<=data%psin<=1 range
    fpolas=fpoli(1)
    sgnbphi=sign(one,fpolas)
    ! Free temporary arrays
    deallocate(wrk,iwrk,wf)

    ! Re-normalize psi after spline computation

    ! Start with un-corrected psi
    psia=data%psia
    psinop=0.0_wp_
    psiant=1.0_wp_

    ! Use provided boundary to set an initial guess
    ! for the search of O/X points
    nbnd=min(size(data%rbnd), size(data%zbnd))
    if (nbnd>0) then
      call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup)
      rbinf=data%rbnd(ibinf)
      rbsup=data%rbnd(ibsup)
      call vmaxmin(data%rbnd,nbnd,rbmin,rbmax)
    else
      zbinf=data%zv(2)
      zbsup=data%zv(nz-1)
      rbinf=data%rv((nr+1)/2)
      rbsup=rbinf
      rbmin=data%rv(2)
      rbmax=data%rv(nr-1)
    end if

    ! Search for exact location of the magnetic axis
    rax0=data%rax
    zax0=data%zax
    call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)

    write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
      'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
    call log_info(msg, mod='equilibrium', proc='set_equil_spline')

    ! search for X-point if params%ixp /= 0

    ixploc = params%ixp
    if(ixploc/=0) then
      if(ixploc<0) then
        call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
        if(psinxptmp/=-1.0_wp_) then
          write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
            'r', r1, 'z', z1, 'ψ', psinxptmp
          call log_info(msg, mod='equilibrium', proc='set_equil_spline')

          zbinf=z1
          psinop=psinoptmp
          psiant=psinxptmp-psinop
          call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
          zbsup=z1
        else
          ixploc=0
        end if
      else
        call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
        if(psinxptmp.ne.-1.0_wp_) then
          write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
            'r', r1, 'z', z1, 'ψ', psinxptmp
          call log_info(msg, mod='equilibrium', proc='set_equil_spline')

          zbsup=z1
          psinop=psinoptmp
          psiant=psinxptmp-psinop
          call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
          zbinf=z1
        else
          ixploc=0
        end if
      end if
    end if

    if (ixploc==0) then
      psinop=psinoptmp
      psiant=one-psinop

      ! Find upper horizontal tangent point
      call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
      zbsup=z1
      rbsup=r1

      ! Find lower horizontal tangent point
      call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
      zbinf=z1
      rbinf=r1
      write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
        'r', rbinf, rbsup, 'z', zbinf, zbsup
      call log_info(msg, mod='equilibrium', proc='set_equil_spline')
    end if

    ! Save Bt value on axis (required in flux_average and used in Jcd def)
    ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def)

    call equinum_fpol(0.0_wp_,btaxis)
    btaxis = btaxis/rmaxis
    btrcen = fpolas/data%rvac
    rcen   = data%rvac

    write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis
    call log_info(msg, mod='equilibrium', proc='set_equil_spline')

    ! Compute rho_pol/rho_tor mapping based on input q profile
    call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn)
    call set_rhospl(sqrt(data%psinr),rhotn)

  end subroutine set_equil_spline


  subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
                        tx,nknt_x,ty,nknt_y,coeff,ierr)
    use const_and_precisions, only : wp_, comp_eps
    use dierckx, only : surfit
    implicit none
! arguments
    integer, intent(in) :: n
    real(wp_), dimension(n), intent(in) :: x, y, z
    real(wp_), dimension(n), intent(in) :: w
    integer, intent(in) :: kspl
    real(wp_), intent(in) :: sspl
    real(wp_), intent(in) :: xmin, xmax, ymin, ymax
    real(wp_), dimension(nknt_x), intent(inout) :: tx
    real(wp_), dimension(nknt_y), intent(inout) :: ty
    integer, intent(inout) :: nknt_x, nknt_y
    real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff
    integer, intent(out) :: ierr
! local variables
    integer :: iopt
    real(wp_) :: resid
    integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest
    real(wp_), dimension(:), allocatable :: wrk1, wrk2
    integer, dimension(:), allocatable :: iwrk

    nxest=nknt_x
    nyest=nknt_y
    ne = max(nxest,nyest)

    km = kspl+1
    u = nxest-km
    v = nyest-km
    b1 = kspl*min(u,v)+kspl+1
    b2 = (kspl+1)*min(u,v)+1
    lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1
    lwrk2 = u*v*(b2+1)+b2
    kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1)
    allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk))

    iopt=0
    call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, &
        sspl,nxest,nyest,ne,comp_eps,nknt_x,tx,nknt_y,ty, &
        coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr)

    deallocate(wrk1,wrk2,iwrk)

  end subroutine scatterspl



  subroutine setqphi_num(psinq,q,psia,rhotn)
    use const_and_precisions, only : pi
    use simplespline, only : difcs
    implicit none
! arguments
    real(wp_), dimension(:), intent(in) :: psinq,q
    real(wp_), intent(in) :: psia
    real(wp_), dimension(:), intent(out), optional :: rhotn
! local variables
    real(wp_), dimension(size(q)) :: phit
    real(wp_) :: dx
    integer, parameter :: iopt=0
    integer :: k,ier

    nq=size(q)
    if(allocated(psinr)) deallocate(psinr)
    if(allocated(cq)) deallocate(cq)
    allocate(psinr(nq),cq(nq,4))

    psinr=psinq
    call difcs(psinr,q,nq,iopt,cq,ier)

! Toroidal flux phi = 2*pi*Integral q dpsi
    phit(1)=0.0_wp_
    do k=1,nq-1
      dx=psinr(k+1)-psinr(k)
      phit(k+1)=phit(k) + dx*(cq(k,1)         + dx*(cq(k,2)/2.0_wp_ +  &
                          dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) )
    end do
    phitedge=phit(nq)
    if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge)
    phitedge=2*pi*psia*phitedge
  end subroutine setqphi_num


  subroutine set_equil_an(data, n)
    ! Computes the analytical equilibrium data and stores them
    ! in their respective global variables, see the top of this file.
    use const_and_precisions, only : pi, zero, one
    use gray_params,          only : equilibrium_data

    implicit none

    ! subroutine arguments
    type(equilibrium_data), intent(in) :: data
    integer, optional,      intent(in) :: n

    ! local variables
    integer, parameter :: nqdef=101
    integer :: i
    real(wp_) :: dr, fq0, fq1, qq, res, rn
    real(wp_) :: rax, zax, a, bax, qax, q1, qexp
    real(wp_), dimension(:), allocatable :: rhotn,rhopn

    rax  = data%rv(1)
    zax  = data%zv(1)
    a    = data%rv(2)
    bax  = data%fpol(1) / rax
    qax  = data%qpsi(1)
    q1   = data%qpsi(2)
    qexp = data%qpsi(3)

    btaxis=bax
    rmaxis=rax
    zmaxis=zax
    btrcen=bax
    rcen=rax
    aminor=a
    zbinf=zmaxis-a
    zbsup=zmaxis+a
    q0=qax
    qa=q1
    alq=qexp
    sgnbphi=sign(one,bax)

    rmxm=rmaxis+aminor
    rmnm=rmaxis-aminor
    zmxm=zbsup
    zmnm=zbinf

    if (present(n)) then
      nq=n
    else
      nq=nqdef
    end if
    if (allocated(psinr)) deallocate(psinr)
    allocate(psinr(nq),rhotn(nq),rhopn(nq))
    
    dr=one/(nq-1)
    rhotn(1)=zero
    psinr(1)=zero
    res=zero
    fq0=zero
    do i=2,nq
      rn=(i-1)*dr
      qq=q0+(q1-q0)*rn**qexp
      fq1=rn/qq
      res=res+0.5_wp_*(fq1+fq0)*dr
      fq0=fq1
      rhotn(i)=rn
      psinr(i)=res
    end do
    
    phitedge=btaxis*aminor**2 ! temporary
    psia=res*phitedge
    phitedge=pi*phitedge    ! final
    psinr=psinr/res
    rhopn=sqrt(psinr)
    
    call set_rhospl(rhopn,rhotn)
  end subroutine set_equil_an


  subroutine set_rhospl(rhop,rhot)
    use simplespline, only : difcs
    implicit none
! arguments
    real(wp_), dimension(:), intent(in) :: rhop, rhot
! local variables
    integer, parameter :: iopt=0
    integer :: ier
    
    nrho=size(rhop)

    if(allocated(rhopr)) deallocate(rhopr)
    if(allocated(rhotr)) deallocate(rhotr)
    if(allocated(crhop)) deallocate(crhop)
    if(allocated(crhot)) deallocate(crhot)
    allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4))
    
    rhopr=rhop
    rhotr=rhot
    call difcs(rhotr,rhopr,nrho,iopt,crhop,ier)
    call difcs(rhopr,rhotr,nrho,iopt,crhot,ier)
  end subroutine set_rhospl



  subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, &
        ddpsidrr,ddpsidzz,ddpsidrz)
    use dierckx, only : fpbisp
    implicit none
! local constants
    integer, parameter :: lwrk=8,liwrk=2
! arguments
    real(wp_), intent(in) :: rpsim,zpsim
    real(wp_), intent(out), optional :: psinv,dpsidr,dpsidz, &
        ddpsidrr,ddpsidzz,ddpsidrz

! local variables
    integer, dimension(liwrk) :: iwrk
    real(wp_), dimension(1) :: rrs,zzs,ffspl
    real(wp_), dimension(lwrk) :: wrk

!   here lengths are measured in meters

    if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. &
        zpsim.le.zmxm .and. zpsim.ge.zmnm) then

      if (present(psinv)) then
        rrs(1)=rpsim
        zzs(1)=zpsim
        call fpbisp(tr,nsr,tz,nsz,cceq,3,3,rrs,1,zzs,1,ffspl, &
                   wrk(1),wrk(5),iwrk(1),iwrk(2))
        psinv=(ffspl(1)-psinop)/psiant
      end if
      if (present(dpsidr)) then
        call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10)
      end if
      if (present(dpsidz)) then
        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01)
      end if
      if (present(ddpsidrr)) then
        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20)
      end if
      if (present(ddpsidzz)) then
        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02)
      end if
      if (present(ddpsidrz)) then
        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11)
      end if
    else
      if(present(psinv)) psinv=-1.0_wp_
      if(present(dpsidr)) dpsidr=0.0_wp_
      if(present(dpsidz)) dpsidz=0.0_wp_
      if(present(ddpsidrr)) ddpsidrr=0.0_wp_
      if(present(ddpsidzz)) ddpsidzz=0.0_wp_
      if(present(ddpsidrz)) ddpsidrz=0.0_wp_
    end if
  end subroutine equinum_psi



  subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc)
    use dierckx, only : fpbisp
    implicit none
! arguments
    real(wp_), intent(in) :: rpsim,zpsim
    integer, intent(in) :: nur,nuz
    real(wp_), intent(out) :: derpsi
    real(wp_), dimension(:) :: cc
! local variables
    integer, dimension(1) :: iwrkr,iwrkz
    real(wp_), dimension(1) :: rrs,zzs,ffspl
    real(wp_), dimension(1,ksplp) :: wrkr
    real(wp_), dimension(1,ksplp) :: wrkz

    rrs(1)=rpsim
    zzs(1)=zpsim
    call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kspl-nur,kspl-nuz, &
         rrs,1,zzs,1,ffspl,wrkr,wrkz,iwrkr,iwrkz)
    derpsi=ffspl(1)*psia
  end subroutine sub_derpsi



  subroutine equinum_fpol(psinv,fpolv,dfpv)
    use dierckx, only : splev,splder
    implicit none
! arguments
    real(wp_), intent(in) :: psinv
    real(wp_), intent(out) :: fpolv
    real(wp_), intent(out), optional :: dfpv
! local variables
    integer :: ier
    real(wp_), dimension(1) :: rrs,ffspl
    real(wp_), dimension(nsf) :: wrkfd
!    
    if(psinv.le.1.0_wp_.and.psinv.ge.0.0_wp_) then
      rrs(1)=psinv
      call splev(tfp,nsf,cfp,3,rrs,ffspl,1,ier)
      fpolv=ffspl(1)
      if(present(dfpv)) then       
        call splder(tfp,nsf,cfp,3,1,rrs,ffspl,1,wrkfd,ier)
        dfpv=ffspl(1)/psia
      end if
    else
      fpolv=fpolas
      if (present(dfpv)) dfpv=0._wp_
    end if
  end subroutine equinum_fpol



  subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, &
      ddpsidrr,ddpsidzz,ddpsidrz)
    use const_and_precisions, only : wp_
    implicit none
! arguments
    real(wp_), intent(in) :: rrm,zzm
    real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, &
      ddpsidrr,ddpsidzz,ddpsidrz
! local variables
    real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot,btaxqq

!  simple model for equilibrium: large aspect ratio
!  outside plasma: analytical continuation, not solution Maxwell eqs

    rpm=sqrt((rrm-rmaxis)**2+(zzm-zmaxis)**2)  !!! rpm==rho_tor[m], rn=rho_tor_norm
    rn=rpm/aminor

    snt=0.0_wp_
    cst=1.0_wp_
    if (rpm > 0.0_wp_) then
      snt=(zzm-zmaxis)/rpm
      cst=(rrm-rmaxis)/rpm
    end if

    if (present(psinv)) then
      rhot=rn
      if(rn <= 1.0_wp_) then
        rhop=frhopol(rhot)
        psinv=rhop**2
      else
        psinv=1.0_wp_+btaxis/(2.0_wp_*psia*qa)*(rpm**2-aminor**2)
        rhop=sqrt(psinv)
      end if
    end if

    if(rn <= 1.0_wp_) then
      qq=q0+(qa-q0)*rn**alq
      btaxqq=btaxis/qq
      dpsidrp=btaxqq*rpm
      dqq=alq*(qa-q0)*rn**(alq-1.0_wp_)
      d2psidrp=btaxqq*(1.0_wp_-rn*dqq/qq)
    else
      btaxqq=btaxis/qa
      dpsidrp=btaxqq*rpm
      d2psidrp=btaxqq
    end if

    if(present(fpolv)) fpolv=btaxis*rmaxis
    if(present(dfpv))  dfpv=0.0_wp_      

    if(present(dpsidr)) dpsidr=dpsidrp*cst
    if(present(dpsidz)) dpsidz=dpsidrp*snt
    if(present(ddpsidrr)) ddpsidrr=btaxqq*snt**2+cst**2*d2psidrp
    if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-btaxqq)
    if(present(ddpsidzz)) ddpsidzz=btaxqq*cst**2+snt**2*d2psidrp
  end subroutine equian



  function frhopol(rhot)
    use utils, only : locate
    use simplespline, only : spli
    implicit none
! arguments
    real(wp_), intent(in) :: rhot
    real(wp_) :: frhopol
! local variables
    integer :: i
    real(wp_) :: dr

    call locate(rhotr,nrho,rhot,i)
    i=min(max(1,i),nrho-1)
    dr=rhot-rhotr(i)
    frhopol=spli(crhop,nrho,i,dr)
  end function frhopol



  function frhopolv(rhot)
    use utils, only : locate
    use simplespline, only : spli
    implicit none
! arguments
    real(wp_), dimension(:), intent(in) :: rhot
    real(wp_), dimension(size(rhot)) :: frhopolv
! local variables
    integer :: i,i0,j
    real(wp_) :: dr

    i0=1
    do j=1,size(rhot)
      call locate(rhotr(i0:),nrho-i0+1,rhot(j),i)
      i=min(max(1,i),nrho-i0)+i0-1
      dr=rhot(j)-rhotr(i)
      frhopolv(j)=spli(crhop,nrho,i,dr)
      i0=i
    end do
  end function frhopolv



  function frhotor(rhop)
    use utils, only : locate
    use simplespline, only : spli
    implicit none
! arguments
    real(wp_), intent(in) :: rhop
    real(wp_) :: frhotor
! local variables
    integer :: i
    real(wp_) :: dr

    call locate(rhopr,nrho,rhop,i)
    i=min(max(1,i),nrho-1)
    dr=rhop-rhopr(i)
    frhotor=spli(crhot,nrho,i,dr)
  end function frhotor



  function fq(psin)
    use const_and_precisions, only : wp_
    use gray_params, only : iequil
    use simplespline, only :spli
    use utils, only : locate
    implicit none
! arguments
    real(wp_), intent(in) :: psin
    real(wp_) :: fq
! local variables
    integer :: i
    real(wp_) :: dps,rn

    if (iequil<2) then
      rn=frhotor(sqrt(psin))
      fq=q0+(qa-q0)*rn**alq
    else
      call locate(psinr,nq,psin,i)
      i=min(max(1,i),nq-1)
      dps=psin-psinr(i)
      fq=spli(cq,nq,i,dps)
    end if
  end function fq



  subroutine bfield(rpsim,zpsim,bphi,br,bz)
    use gray_params, only : iequil
    implicit none
! arguments
    real(wp_), intent(in) :: rpsim,zpsim
    real(wp_), intent(out), optional :: bphi,br,bz
! local variables
    real(wp_) :: psin,fpol

    if (iequil < 2) then
      call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
      if (present(bphi)) bphi=bphi/rpsim
    else
      call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br)
      if (present(bphi)) then
        psin=bphi
        call equinum_fpol(psin,fpol)
        bphi=fpol/rpsim
      end if
    end if
    if (present(br)) br=-br/rpsim
    if (present(bz)) bz= bz/rpsim
  end subroutine bfield



  subroutine tor_curr(rpsim,zpsim,ajphi)
    use const_and_precisions, only : wp_,ccj=>mu0inv
    use gray_params, only : iequil
    implicit none
! arguments
    real(wp_) :: rpsim,zpsim,ajphi
! local variables
    real(wp_) :: bzz,dbvcdc13,dbvcdc31
    real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
    
    if(iequil < 2) then
      call equian(rpsim,zpsim,dpsidr=dpsidr,                     &
                  ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
    else 
      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,                &
                       ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
    end if 
    bzz= dpsidr/rpsim
    dbvcdc13=-ddpsidzz/rpsim
    dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
    ajphi=ccj*(dbvcdc13-dbvcdc31)
  end subroutine tor_curr



  subroutine psi_raxis(psin,r1,r2)
    use const_and_precisions, only : wp_
    use gray_params,          only : iequil
    use dierckx,              only : profil, sproota
    use logger,               only : log_error

    implicit none

    ! subroutine arguments
    real(wp_) :: psin,r1,r2

    ! local constants
    integer, parameter :: mest=4

    ! local variables
    integer :: iopt,ier,m
    real(wp_) :: zc,val
    real(wp_), dimension(mest) :: zeroc
    real(wp_), dimension(nsr) :: czc
    character(64) :: msg

    if (iequil < 2) then
      val=frhotor(sqrt(psin))
      r1=rmaxis-val*aminor
      r2=rmaxis+val*aminor
    else
      iopt=1
      zc=zmaxis
      call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier)
      if (ier > 0) then
        write (msg, '("profil failed with error ",g0)') ier
        call log_error(msg, mod='equilibrium', proc='psi_raxis')
      end if

      val=psin*psiant+psinop
      call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
      r1=zeroc(1)
      r2=zeroc(2)
    end if  
  end subroutine psi_raxis



  subroutine tor_curr_psi(psin,ajphi)
    use const_and_precisions, only : wp_
    implicit none
! arguments
    real(wp_) :: psin,ajphi
! local variables
    real(wp_) :: r1,r2
    call psi_raxis(psin,r1,r2)
    call tor_curr(r2,zmaxis,ajphi)
  end subroutine tor_curr_psi



  subroutine points_ox(rz,zz,rf,zf,psinvf,info)
    use const_and_precisions, only : comp_eps
    use minpack,              only : hybrj1
    use logger,               only : log_error, log_debug

    implicit none

    ! local constants
    integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2

    ! arguments
    real(wp_), intent(in) :: rz,zz
    real(wp_), intent(out) :: rf,zf,psinvf
    integer, intent(out) :: info

    ! local variables
    real(wp_) :: tol
    real(wp_), dimension(n) :: xvec,fvec
    real(wp_), dimension(lwa) :: wa
    real(wp_), dimension(ldfjac,n) :: fjac
    character(256) :: msg

    xvec(1)=rz
    xvec(2)=zz
    tol = sqrt(comp_eps)
    call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
    if(info /= 1) then
      write (msg, '("O,X coordinates:",2(x,", ",g0.3))') xvec
      call log_debug(msg, mod='equilibrium', proc='points_ox')
      write (msg, '("hybrj1 failed with error ",g0)') info
      call log_error(msg, mod='equilibrium', proc='points_ox')
    end if
    rf=xvec(1)
    zf=xvec(2)
    call equinum_psi(rf,zf,psinvf)
  end subroutine points_ox



  subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
    use logger, only : log_error

    implicit none

    ! subroutine arguments
    integer, intent(in) :: n,iflag,ldfjac
    real(wp_), dimension(n), intent(in) :: x
    real(wp_), dimension(n), intent(inout) :: fvec
    real(wp_), dimension(ldfjac,n), intent(inout) :: fjac

    ! local variables
    real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
    character(64) :: msg

    select case(iflag)
    case(1)
      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
      fvec(1) = dpsidr/psia
      fvec(2) = dpsidz/psia    
    case(2)
      call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, &
           ddpsidrz=ddpsidrz)
      fjac(1,1) = ddpsidrr/psia
      fjac(1,2) = ddpsidrz/psia
      fjac(2,1) = ddpsidrz/psia
      fjac(2,2) = ddpsidzz/psia
    case default
      write (msg, '("invalid iflag: ",g0)')
      call log_error(msg, mod='equilibrium', proc='fcnox')
    end select
  end subroutine fcnox



  subroutine points_tgo(rz,zz,rf,zf,psin0,info)
    use const_and_precisions, only : comp_eps
    use minpack,              only : hybrj1mv
    use logger,               only : log_error, log_debug

    implicit none

    ! local constants
    integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2

    ! arguments
    real(wp_), intent(in) :: rz,zz,psin0
    real(wp_), intent(out) :: rf,zf
    integer, intent(out) :: info
    character(256) :: msg

    ! local variables
    real(wp_) :: tol
    real(wp_), dimension(n) :: xvec,fvec,f0
    real(wp_), dimension(lwa) :: wa
    real(wp_), dimension(ldfjac,n) :: fjac

    xvec(1)=rz
    xvec(2)=zz
    f0(1)=psin0
    f0(2)=0.0_wp_
    tol = sqrt(comp_eps)
    call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
    if(info /= 1) then
      write (msg, '("R,z coordinates:",5(x,g0.3))') xvec, rz, zz, psin0
      call log_debug(msg, mod='equilibrium', proc='points_tgo')
      write (msg, '("hybrj1mv failed with error ",g0)') info
      call log_error(msg, mod='equilibrium', proc='points_tgo')
    end if
    rf=xvec(1)
    zf=xvec(2)
  end



  subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
    use const_and_precisions, only : wp_
    use logger,               only : log_error

    implicit none

    ! subroutine arguments
    integer, intent(in) :: n,ldfjac,iflag
    real(wp_), dimension(n), intent(in) :: x,f0
    real(wp_), dimension(n), intent(inout) :: fvec
    real(wp_), dimension(ldfjac,n), intent(inout) :: fjac

    ! local variables
    real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
    character(64) :: msg

    select case(iflag)
    case(1)
      call equinum_psi(x(1),x(2),psinv,dpsidr)
      fvec(1) = psinv-f0(1)
      fvec(2) = dpsidr/psia-f0(2)
    case(2)
      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, &
           ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
      fjac(1,1) = dpsidr/psia
      fjac(1,2) = dpsidz/psia
      fjac(2,1) = ddpsidrr/psia
      fjac(2,2) = ddpsidrz/psia
    case default
      write (msg, '("invalid iflag: ",g0)')
      call log_error(msg, mod='equilibrium', proc='fcntgo')
    end select
  end subroutine fcntgo



  subroutine unset_equil_spline
    implicit none

    if (allocated(tr))     deallocate(tr)
    if (allocated(tz))     deallocate(tz)
    if (allocated(tfp))    deallocate(tfp)
    if (allocated(cfp))    deallocate(cfp)
    if (allocated(cceq))   deallocate(cceq)
    if (allocated(cceq01)) deallocate(cceq01)
    if (allocated(cceq10)) deallocate(cceq10)
    if (allocated(cceq02)) deallocate(cceq02)
    if (allocated(cceq20)) deallocate(cceq20)
    if (allocated(cceq11)) deallocate(cceq11)
    nsr = 0
    nsz = 0
    nsf = 0
  end subroutine unset_equil_spline



  subroutine unset_q
    implicit none
    
    if (allocated(psinr)) deallocate(psinr)
    if (allocated(cq))    deallocate(cq)
    nq = 0
  end subroutine unset_q



  subroutine unset_rho_spline
    implicit none
    
    if(allocated(rhopr)) deallocate(rhopr)
    if(allocated(rhotr)) deallocate(rhotr)
    if(allocated(crhop)) deallocate(crhop)
    if(allocated(crhot)) deallocate(crhot)
    nrho=0
  end subroutine unset_rho_spline

end module equilibrium