! This module handles the loading, interpolation and evaluation of the
! MHD equilibrium data (poloidal current function, safety factor,
! poloidal flux, magnetic field, etc.)
!
! Two kinds of data are supported: analytical (suffix `_an` in the
! subroutine names) or numerical (suffix `_spline`). For the latter, the
! the data is interpolated using splines.
module equilibrium
  use const_and_precisions, only : wp_
  use splines,              only : spline_simple, spline_1d, spline_2d, linear_1d

  implicit none

  ! Parameters of the analytical equilibrium model
  type analytic_model
    real(wp_) :: q0    ! Safety factor at the magnetic axis
    real(wp_) :: q1    ! Safety factor at the edge
    real(wp_) :: alpha ! Exponent for the q(ρ_p) power law
    real(wp_) :: R0    ! R of the magnetic axis (m)
    real(wp_) :: z0    ! z of the magnetic axis (m)
    real(wp_) ::  a    ! Minor radius (m)
    real(wp_) :: B0    ! Magnetic field at the magnetic axis (T)
  end type

  ! Order of the splines
  integer, parameter :: kspl=3, ksplp=kspl + 1

  ! Global variable storing the state of the module

  ! Splines
  type(spline_1d),     save :: fpol_spline              ! Poloidal current function F(ψ)
  type(spline_2d),     save :: psi_spline               ! Normalised poloidal flux ψ(R,z)
  type(spline_simple), save :: q_spline                 ! Safey factor q(ψ)
  type(linear_1d),     save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p)

  ! Analytic model
  type(analytic_model), save :: model

  ! More parameters
  real(wp_), save :: fpolas          ! Poloidal current at the edge (r=a), F_a
  real(wp_), save :: psia            ! Poloidal flux at the edge (r=a), ψ_a
  real(wp_), save :: phitedge        ! Toroidal flux at the edge
  real(wp_), save :: btaxis          ! B₀: B_φ = B₀R₀/R
  real(wp_), save :: rmaxis, zmaxis  ! R,z of the magnetic axis (O point)
  real(wp_), save :: sgnbphi         ! Sign of B_φ (>0 means CCW)
  real(wp_), save :: btrcen          ! used only for jcd_astra def.
  real(wp_), save :: rcen            ! Major radius R₀, computed as fpolas/btrcen
  real(wp_), save :: rmnm, rmxm      ! R interval of the equilibrium definition
  real(wp_), save :: zmnm, zmxm      ! z interval of the equilibrium definition
  real(wp_), save :: zbinf, zbsup

  private
  public read_eqdsk, read_equil_an           ! Reading data files
  public scale_equil, change_cocos           ! Transforming data
  public equian                              ! Analytical model
  public equinum_psi, equinum_fpol           ! Interpolated data
  public fq, bfield, tor_curr, tor_curr_psi  ! Accessing local quantities
  public frhotor, frhopol                    ! Converting between poloidal/toroidal flux
  public set_equil_spline, set_equil_an      ! Initialising internal state
  public unset_equil_spline                  ! Deinitialising internal state

  ! Members exposed for magsurf_data
  public kspl, psi_spline, q_spline, points_tgo, model

  ! Members exposed to gray_core and more
  public psia, phitedge
  public btaxis, rmaxis, zmaxis, sgnbphi
  public btrcen, rcen
  public rmnm, rmxm, zmnm, zmxm
  public zbinf, zbsup

contains

  subroutine read_eqdsk(params, data, err, 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,                      intent(out) :: err
    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

    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')
      err = 1
      return
    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, err, 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,                intent(out) :: err
    integer, optional,      intent(in)  :: unit

    ! local variables
    integer :: i, u, nlim

    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')
      err = 1
      return
    end if

    read(u, *) model%R0, model%z0, model%a
    read(u, *) model%B0
    read(u, *) model%q0, model%q1, model%alpha

    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)
    ! Extracts the sign and units conventions from a COCOS index
    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 scale_equil(params, data)
    ! Rescale the magnetic field (B) and the plasma current (I_p)
    ! 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 B_φ and I_p BEFORE scaling.
    ! 2. cocos=3 assumed: positive toroidal direction is CCW from above
    ! 3. B_φ and I_p 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
    use gray_params,          only : iequil

    implicit none

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

    ! Notes on cocos=3
    ! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0.
    ! 2. ψ_a = ψ_edge - ψ_axis < 0.
    ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
    ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).

    if (iequil < 2) then
      ! Analytical model

      ! Apply signs
      if (params%sgni /= 0) then
        model%q0 = sign(model%q0, -params%sgni*one)
        model%q1 = sign(model%q1, -params%sgni*one)
      end if
      if (params%sgnb /= 0) then
        model%B0 = sign(model%B0, +params%sgnb*one)
      end if
      ! Rescale
      model%B0 = model%B0 * params%factb
    else
      ! Numeric data

      ! Apply signs
      if (params%sgni /= 0) &
        data%psia = sign(data%psia, -params%sgni*one)
      if (params%sgnb /= 0) &
        data%fpol = sign(data%fpol, +params%sgnb*one)
      ! Rescale
      data%psia = data%psia * params%factb
      data%fpol = data%fpol * params%factb

      ! Compute the signs to be shown in the outputs header when cocos≠0,10.
      ! Note: In these cases the values sgni,sgnb from gray.ini are unused.
      params%sgni = int(sign(one, -data%psia))
      params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
    end if
  end subroutine scale_equil


  subroutine set_equil_spline(params, data, err)
    ! 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 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
    integer,                      intent(out) :: err

    ! local variables
    integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup
    real(wp_) :: tension, rax0, zax0, psinoptmp, psinxptmp
    real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1
    real(wp_) :: psiant, psinop
    real(wp_),  dimension(size(data%psinr)) :: rhotn
    real(wp_),  dimension(:),  allocatable :: rv1d, zv1d, fvpsi, wf
    integer :: ixploc, info, i, j, ij
    character(256) :: msg ! for log messages formatting

    ! compute array sizes
    nr = size(data%rv)
    nz = size(data%zv)
    nrz  = nr*nz
    npsi = size(data%psinr)
    nrest  = nr + ksplp
    nzest  = nz + ksplp
    npsest = npsi + ksplp

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

    ! Spline interpolation of ψ(R, z)

    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))

      ! allocate knots and spline coefficients arrays
      if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x)
      if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y)
      if (allocated(psi_spline%coeffs))  deallocate(psi_spline%coeffs)
      allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest))
      allocate(psi_spline%coeffs(nrz))

      ! 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 ?
      psi_spline%nknots_x=nr/4+4
      psi_spline%nknots_y=nz/4+4
      tension = params%ssplps
      call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
                      rmnm, rmxm, zmnm, zmxm,                    &
                      psi_spline%knots_x, psi_spline%nknots_x,   &
                      psi_spline%knots_y, psi_spline%nknots_y,   &
                      psi_spline%coeffs, err)

      ! if failed, re-fit with an interpolating spline (zero tension)
      if(err == -1) then
        err = 0
        tension = 0
        psi_spline%nknots_x=nr/4+4
        psi_spline%nknots_y=nz/4+4
        call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
                        rmnm, rmxm, zmnm, zmxm,                    &
                        psi_spline%knots_x, psi_spline%nknots_x,   &
                        psi_spline%knots_y, psi_spline%nknots_y,   &
                        psi_spline%coeffs, err)
      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 ψ array to 1D (transposed)
      allocate(fvpsi(nrz))
      fvpsi = reshape(transpose(data%psin), [nrz])

      ! compute spline coefficients
      call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
                           range=[rmnm, rmxm, zmnm, zmxm],  &
                           tension=params%ssplps, err=err)

      ! if failed, re-fit with an interpolating spline (zero tension)
      if(err == -1) then
        call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
                             range=[rmnm, rmxm, zmnm, zmxm],  &
                             tension=zero)
        err = 0
      end if
      deallocate(fvpsi)
    end if

    if (err /= 0) then
      err = 2
      return
    end if

    ! compute spline coefficients for ψ(R,z) partial derivatives
    call psi_spline%init_deriv(nr, nz, 1, 0)  ! ∂ψ/∂R
    call psi_spline%init_deriv(nr, nz, 0, 1)  ! ∂ψ/∂z
    call psi_spline%init_deriv(nr, nz, 1, 1)  ! ∂²ψ/∂R∂z
    call psi_spline%init_deriv(nr, nz, 2, 0)  ! ∂²ψ/∂R²
    call psi_spline%init_deriv(nr, nz, 0, 2)  ! ∂²ψ/∂z²

    ! Spline interpolation of F(ψ)

    ! give a small weight to the last point
    allocate(wf(npsi))
    wf(1:npsi-1) = 1
    wf(npsi)     = 1.0e2_wp_
    call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], &
                          weights=wf, tension=params%ssplf)
    deallocate(wf)

    ! set vacuum value used outside 0 ≤ ψ ≤ 1 range
    fpolas  = fpol_spline%eval(data%psinr(npsi))
    sgnbphi = sign(one ,fpolas)

    ! Re-normalize ψ_n after the spline computation
    ! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma

    ! Start with un-corrected ψ_n
    psia   = data%psia
    psinop = 0  ! ψ_n(O point)
    psiant = 1  ! ψ_n(X point) - ψ_n(O point)

    ! 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

    ! Adjust all the B-spline coefficients
    ! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n
    ! it's enough to correct (shift and scale) the c_ij
    psi_spline%coeffs  = (psi_spline%coeffs - psinop) / psiant

    ! Same for the derivatives
    psi_spline%partial(1,0)%ptr = psi_spline%partial(1,0)%ptr / psiant
    psi_spline%partial(0,1)%ptr = psi_spline%partial(0,1)%ptr / psiant
    psi_spline%partial(2,0)%ptr = psi_spline%partial(2,0)%ptr / psiant
    psi_spline%partial(0,2)%ptr = psi_spline%partial(0,2)%ptr / psiant
    psi_spline%partial(1,1)%ptr = psi_spline%partial(1,1)%ptr / psiant

    ! Do the opposite scaling to preserve un-normalised values
    ! Note: this is only used for the poloidal magnetic field
    psia = psia * psiant

    ! 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(zero, 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 ρ_p/ρ_t mapping based on the input q profile
    call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn)
    call set_rho_spline(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)
    ! Computes the spline interpolation of a surface when
    ! the data points are irregular, i.e. not on a uniform grid
    use const_and_precisions, only : comp_eps
    use dierckx,              only : surfit

    implicit none

    ! subroutine 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)
    ! Computes the spline of the safety factor q(ψ)
    use const_and_precisions, only : pi

    implicit none

    ! subroutine 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 :: k

    call q_spline%init(psinq, q, size(q))

    ! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ
    phit(1)=0
    do k=1,q_spline%ndata-1
      dx=q_spline%data(k+1)-q_spline%data(k)
      phit(k+1)=phit(k) + dx*(q_spline%coeffs(k,1)   + dx*(q_spline%coeffs(k,2)/2 +  &
                          dx*(q_spline%coeffs(k,3)/3 + dx* q_spline%coeffs(k,4)/4) ) )
    end do
    phitedge=phit(q_spline%ndata)
    if(present(rhotn)) rhotn(1:q_spline%ndata)=sqrt(phit/phitedge)
    phitedge=2*pi*psia*phitedge
  end subroutine setqphi_num


  subroutine set_equil_an
    ! 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, one

    implicit none

    real(wp_) :: dq

    btaxis  = model%B0
    rmaxis  = model%R0
    zmaxis  = model%z0
    btrcen  = model%B0
    rcen    = model%R0
    zbinf   = model%z0 - model%a
    zbsup   = model%z0 + model%a
    sgnbphi = sign(one, model%B0)

    rmxm = model%R0 + model%a
    rmnm = model%R0 - model%a
    zmxm = zbsup
    zmnm = zbinf

    phitedge = model%B0 * pi * model%a**2
    dq   = (model%q1 - model%q0) / (model%alpha/2 + 1)
    psia = 1/(2*pi) * phitedge / (model%q0 + dq)
  end subroutine set_equil_an


  subroutine set_rho_spline(rhop, rhot)
    ! Computes the splines for converting between the poloidal (ρ_p)
    ! and toroidal (ρ_t) normalised radii

    implicit none

    ! subroutine arguments
    real(wp_), dimension(:), intent(in) :: rhop, rhot

    call rhop_spline%init(rhot, rhop, size(rhop))
    call rhot_spline%init(rhop, rhot, size(rhot))
  end subroutine set_rho_spline


  subroutine equinum_psi(R, z, psi, dpsidr, dpsidz,  &
                         ddpsidrr, ddpsidzz, ddpsidrz)
    ! Computes the normalised poloidal flux ψ (and its derivatives)
    ! as a function of (R, z) using the numerical equilibrium

    implicit none

    ! subroutine arguments
    real(wp_), intent(in)            :: R, z
    real(wp_), intent(out), optional :: psi, dpsidr, dpsidz
    real(wp_), intent(out), optional :: ddpsidrr, ddpsidzz, ddpsidrz

    ! Note: here lengths are measured in meters
    if (R <= rmxm .and. R >= rmnm .and. &
        z <= zmxm .and. z >= zmnm) then

      if (present(psi))      psi      = psi_spline%eval(R, z)
      if (present(dpsidr))   dpsidr   = psi_spline%deriv(R, z, 1, 0)
      if (present(dpsidz))   dpsidz   = psi_spline%deriv(R, z, 0, 1)
      if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0)
      if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2)
      if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1)
    else
      if (present(psi))      psi      = -1
      if (present(dpsidr))   dpsidr   = 0
      if (present(dpsidz))   dpsidz   = 0
      if (present(ddpsidrr)) ddpsidrr = 0
      if (present(ddpsidzz)) ddpsidzz = 0
      if (present(ddpsidrz)) ddpsidrz = 0
    end if
  end subroutine equinum_psi


  subroutine equinum_fpol(psi, fpol, dfpol)
    ! Computes the poloidal current function F(ψ)
    ! (and its derivative) using the numerical equilibrium

    implicit none

    ! subroutine arguments
    real(wp_), intent(in)            :: psi
    real(wp_), intent(out)           :: fpol
    real(wp_), intent(out), optional :: dfpol

    if(psi <= 1 .and. psi >= 0) then
      fpol = fpol_spline%eval(psi)
      if (present(dfpol)) dfpol = fpol_spline%deriv(psi) / psia
    else
      fpol = fpolas
      if (present(dfpol)) dfpol = 0
    end if
  end subroutine equinum_fpol


  subroutine equian(R, z, psin, fpolv, dfpv, dpsidr, dpsidz,  &
                    ddpsidrr, ddpsidzz, ddpsidrz)
    ! Computes all MHD equilibrium parameters from a simple analytical model
    implicit none

    ! subroutine arguments
    real(wp_), intent(in) :: R, z
    real(wp_), intent(out), optional :: psin, fpolv, dfpv, dpsidr, dpsidz,  &
                                        ddpsidrr, ddpsidzz, ddpsidrz

    ! variables
    real(wp_) :: r_p   ! poloidal radius
    real(wp_) :: rho_p ! poloidal radius normalised

    ! ρ_p is just the geometrical poloidal radius divided by a
    r_p = hypot(R - model%R0, z - model%z0)
    rho_p = r_p / model%a

    ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
    if (present(psin)) psin = rho_p**2

    ! ∂ψ/∂R, ∂ψ/∂z
    if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2
    if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2

    ! ∂²ψ/∂R², ∂²ψ/∂z², ∂²ψ/∂R∂z
    if (present(ddpsidrr)) ddpsidrr = 2/model%a**2
    if (present(ddpsidzz)) ddpsidzz = 2/model%a**2
    if (present(ddpsidrz)) ddpsidrz = 0

    ! F(ψ) = B₀⋅R₀, a constant
    if (present(fpolv)) fpolv = model%B0 * model%R0
    if (present(dfpv))   dfpv = 0
  end subroutine equian


  function frhotor(rho_p)
    ! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius
    use gray_params, only : iequil

    implicit none

    ! function arguments
    real(wp_), intent(in) :: rho_p
    real(wp_)             :: frhotor

    if (iequil < 2) then
      ! Analytical model
      block
        ! The change of variable is obtained by integrating
        !
        !   q(ψ) = 1/2π ∂Φ/∂ψ
        !
        ! and defining ψ = ψ_a ρ_p²,  Φ = Φ_a ρ_t².
        ! The result is:
        !
        !   - ψ_a = 1/2π Φ_a / [q₀ + Δq]
        !
        !   - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)]
        !
        ! where Δq = (q₁ - q₀)/(α/2 + 1)
        real(wp_) :: dq
        dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
        frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) &
                               / (model%q0 + dq))
      end block
    else
      ! Numerical data
      frhotor = rhot_spline%eval(rho_p)
    end if
  end function frhotor


  function frhopol(rho_t)
    ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
    use gray_params,          only : iequil
    use const_and_precisions, only : comp_eps

    implicit none

    ! function arguments
    real(wp_), intent(in) :: rho_t
    real(wp_)             :: frhopol

    if (iequil < 2) then
      ! Analytical model
      block
        ! In general there is no closed form for ρ_p(ρ_t) in the
        ! analytical model, we thus solve numerically the equation
        ! ρ_t(ρ_p) = ρ_t₀ for ρ_p.
        use minpack, only : hybrj1

        real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
        integer   :: info

        rho_p = [rho_t]  ! first guess, ρ_p ≈ ρ_t
        call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, &
                    ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7)
        frhopol = rho_p(1)
      end block
    else
      ! Numerical data
      frhopol = rhop_spline%eval(rho_t)
    end if

    contains

      subroutine equation(n, x, f, df, ldf, flag)
        ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0
        implicit none

        ! optimal step size
        real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)

        ! subroutine arguments
        integer,   intent(in)    :: n, ldf, flag
        real(wp_), intent(in)    :: x(n)
        real(wp_), intent(inout) :: f(n), df(ldf,n)

        if (flag == 1) then
          ! return f(x)
          f(1) = frhotor(x(1)) - rho_t
        else
          ! return f'(x), computed numerically
          df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e)
        end if
      end subroutine

  end function frhopol


  function fq(psin)
    ! Computes the safety factor q as a function of the normalised
    ! poloidal flux ψ.
    !
    ! Note: this returns the absolute value of q.
    use gray_params, only : iequil

    implicit none

    ! function arguments
    real(wp_), intent(in) :: psin
    real(wp_) :: fq

    if (iequil < 2) then
      ! Analytical model
      ! The safety factor is a power law in ρ_p:
      !   q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
      block
        real(wp_) :: rho_p
        rho_p = sqrt(psin)
        fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha)
      end block
    else
      ! Numerical data
      fq = q_spline%eval(psin)
    end if
  end function fq


  subroutine bfield(rpsim, zpsim, bphi, br, bz)
    ! Computes the magnetic field as a function of
    ! (R, z) in cylindrical coordinates
    use gray_params, only : iequil

    implicit none

    ! subroutine 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
      ! Analytical model
      call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
      if (present(bphi)) bphi=bphi/rpsim
    else
      ! Numerical data
      call equinum_psi(rpsim,zpsim,psi=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


  function tor_curr(rpsim,zpsim) result(jphi)
    ! Computes the toroidal current Jφ as a function of (R, z)
    use const_and_precisions, only : ccj=>mu0inv
    use gray_params,          only : iequil

    implicit none

    ! function arguments
    real(wp_), intent(in) :: rpsim, zpsim
    real(wp_)             :: jphi

    ! local variables
    real(wp_) :: bzz,dbvcdc13,dbvcdc31
    real(wp_) :: dpsidr,ddpsidrr,ddpsidzz

    if(iequil < 2) then
      ! Analytical model
      call equian(rpsim,zpsim,dpsidr=dpsidr, &
                  ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
    else
      ! Numerical data
      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, &
                       ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
    end if
    bzz= dpsidr/rpsim
    dbvcdc13=-ddpsidzz/rpsim
    dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
    jphi=ccj*(dbvcdc13-dbvcdc31)
  end function tor_curr


  function tor_curr_psi(psin) result(jphi)
    ! Computes the toroidal current Jφ as a function of ψ

    implicit none

    ! function arguments
    real(wp_), intent(in) :: psin
    real(wp_)             :: jphi

    ! local variables
    real(wp_) :: r1, r2
    call psi_raxis(psin, r1, r2)
    jphi = tor_curr(r2, zmaxis)
  end function tor_curr_psi


  subroutine psi_raxis(psin,r1,r2)
    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(psi_spline%nknots_x) :: czc
    character(64) :: msg

    if (iequil < 2) then
      ! Analytical model
      val = sqrt(psin)
      r1 = model%R0 - val*model%a
      r2 = model%R0 + val*model%a
    else
      ! Numerical data
      iopt=1
      zc=zmaxis
      call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &
                  psi_spline%knots_y, psi_spline%nknots_y,       &
                  psi_spline%coeffs, kspl, kspl, zc,             &
                  psi_spline%nknots_x, 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

      call sproota(psin, psi_spline%knots_x, psi_spline%nknots_x, &
                   czc, zeroc, mest, m, ier)
      r1=zeroc(1)
      r2=zeroc(2)
    end if
  end subroutine psi_raxis


  subroutine points_ox(rz,zz,rf,zf,psinvf,info)
    ! Finds the location of the O,X points
    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
      fvec(2) = dpsidz
    case(2)
      call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, &
           ddpsidrz=ddpsidrz)
      fjac(1,1) = ddpsidrr
      fjac(1,2) = ddpsidrz
      fjac(2,1) = ddpsidrz
      fjac(2,2) = ddpsidzz
    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 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-f0(2)
    case(2)
      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, &
           ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
      fjac(1,1) = dpsidr
      fjac(1,2) = dpsidz
      fjac(2,1) = ddpsidrr
      fjac(2,2) = ddpsidrz
    case default
      write (msg, '("invalid iflag: ",g0)')
      call log_error(msg, mod='equilibrium', proc='fcntgo')
    end select
  end subroutine fcntgo


  subroutine unset_equil_spline
    ! Unsets the splines global variables, see the top of this file.
    implicit none

    call fpol_spline%deinit
    call psi_spline%deinit
    call q_spline%deinit
    call rhop_spline%deinit
    call rhot_spline%deinit
  end subroutine unset_equil_spline

end module equilibrium