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

    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

    u = get_free_unit(unit)

    ! Open the G-EQDSK file
    open(file=trim(params%filenm), status='old', action='read', unit=u)

    ! 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,rv,zv,fpol,q,rlim,zlim,unit)
    use utils, only : get_free_unit
    implicit none
! arguments
    character(len=*), intent(in) :: filenm
    integer, intent(in) :: ipass
    integer, optional, intent(in) :: unit
    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim
! local variables
    integer :: i, u, nlim
    real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen

    u = get_free_unit(unit)

    open(file=trim(filenm),status='old',action='read',unit=u)
    read(u,*) rr0m,zr0m,rpam
    read(u,*) b0
    read(u,*) q0,qa,alq
!    rcen=rr0m
!    btrcen=b0
!    b0=isgnbphi*b0*factb
!    rvac=rr0m
!    rax=rr0m
!    zax=zr0m
!    zbmin=(zr0-rpa)/1.0e2_wp_
!    zbmax=(zr0+rpa)/1.0e2_wp_
    if(allocated(rv)) deallocate(rv)
    if(allocated(zv)) deallocate(zv)
    if(allocated(fpol)) deallocate(fpol)
    if(allocated(q)) deallocate(q)
    allocate(rv(2),zv(1),fpol(1),q(3))
    rv(1)=rr0m
    rv(2)=rpam
    zv(1)=zr0m
    fpol(1)=b0*rr0m
    q(1)=q0
    q(2)=qa
    q(3)=alq

  if(ipass.ge.2) then

  ! get size of boundary and limiter arrays and allocate them
      read (u,*) nlim
      if (allocated(rlim)) deallocate(rlim)
      if (allocated(zlim)) deallocate(zlim)
  
  ! store boundary and limiter data
      if(nlim>0) then
        allocate(rlim(nlim),zlim(nlim))
        read(u,*) (rlim(i),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_eqspl(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

    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

    ! 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)
    print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp

    ! 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
          print'(a,2f8.4,es12.5)','X-point',r1,z1,psinxptmp
          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
          print'(a,2f8.4,e16.8)','X-point',r1,z1,psinxptmp
          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
      print'(a,4f8.4)','no X-point  ',rbinf,zbinf,rbsup,zbsup
    end if
    print*,' '

    ! 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
    print '(a,f8.4)', 'BT_centr=', btrcen
    print '(a,f8.4)', 'BT_axis =', btaxis

    ! 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_eqspl


  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_equian(rax,zax,a,bax,qax,q1,qexp,n)
    use const_and_precisions, only : pi,zero,one
    implicit none
! arguments
    real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp
    integer, intent(in), optional :: n
! local variables
    integer, parameter :: nqdef=101
    integer :: i
    real(wp_) :: dr,fq0,fq1,qq,res,rn
    real(wp_), dimension(:), allocatable :: rhotn,rhopn

    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_equian



  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
    implicit none
! local constants
    integer, parameter :: mest=4
! arguments
    real(wp_) :: psin,r1,r2
! local variables
    integer :: iopt,ier,m
    real(wp_) :: zc,val
    real(wp_), dimension(mest) :: zeroc
    real(wp_), dimension(nsr) :: czc

    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.gt.0) print*,'    profil =',ier
      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
    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

    xvec(1)=rz
    xvec(2)=zz
    tol = sqrt(comp_eps)
    call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
      if(info.gt.1)  then
        print'(a,i2,a,2f8.4)','    info subr points_ox =',info, &
          '  O/X coord.',xvec
      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)
    implicit none
! 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

    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
      print*,'iflag undefined'
    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
    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
! 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.gt.1)  then
      print'(a,i2,a,5f8.4)','    info subr points_tgo =',info, &
        '  R,z coord.',xvec,rz,zz,psin0
    end if
    rf=xvec(1)
    zf=xvec(2)
  end



  subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
    use const_and_precisions, only : wp_
    implicit none
! 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
! internal variables
    real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz

    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
      print*,'iflag undefined'
    end select
  end subroutine fcntgo



  subroutine unset_eqspl
    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_eqspl



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



  subroutine unset_rhospl
    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_rhospl

end module equilibrium