module interp_eqprof
  use const_and_precisions, only : wp_
  implicit none
  
! equidata
  ! === 2D array psi(R,z) ==========================================
  INTEGER, SAVE :: nr,nz
  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rv,zv
  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: psin
  ! === 1D array Fpol(psi), q(psi), rhotor(psi) ====================
  !INTEGER, SAVE :: npsieq
  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopnr,fpol,qpsi
  
  ! === 1D array plasma boundary Rbnd_i, Zbnd_i ====================
  INTEGER, SAVE :: nbbbs
  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rbbbs,zbbbs
  ! === 1D array limiter Rlim_i, Zlim_i   ==> move in wall/reflections
  INTEGER, SAVE :: nlim
  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim

  REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis
  REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm,dr,dz
  REAL(wp_), SAVE :: zbmin,zbmax
  REAL(wp_), SAVE :: phitedge
  REAL(wp_), SAVE :: rup,zup,rlw,zlw
  REAL(wp_), SAVE :: rcen,btrcen ! rcen unused, btrcen used only for Jcd_ASTRA def.

  ! === 2D spline psi(R,z), normalization and derivatives ==========
  INTEGER, SAVE :: nrest, nzest, lw10, lw01, lw20, lw02, lw11
  INTEGER :: nsr,nsz
  REAL(wp_), SAVE :: psia, psiant, psinop
  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
  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhot,cq

  ! !!! 2D B(R,z) array. Computed in bfield_res,
  ! !!! used by cniteq to plot resonant field contour lines.
  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: btotal

contains

  subroutine alloc_equilvec(ier)
    implicit none
    integer, intent(out) :: ier

    if(nr.le.0.or.nz.le.0) then
      ier = -1
      return
    end if

    nrest=nr+4
    nzest=nz+4
    lw10=nr*3+nz*4+nr*nz
    lw01=nr*4+nz*3+nr*nz
    lw20=nr*2+nz*4+nr*nz
    lw02=nr*4+nz*2+nr*nz
    lw11=nr*3+nz*3+nr*nz
    
    call dealloc_equilvec
    allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), &
     btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), &
     cceq20(lw20),cceq11(lw11),psin(nr,nz),psinr(nr),rhopnr(nr),fpol(nr), &
     qpsi(nr), stat=ier)
    if (ier/=0) call dealloc_equilvec
  end subroutine alloc_equilvec

  subroutine dealloc_equilvec
    implicit none
    if(allocated(rv)) deallocate(rv)
    if(allocated(zv)) deallocate(zv)
    if(allocated(tr)) deallocate(tr)
    if(allocated(tz)) deallocate(tz)
    if(allocated(tfp)) deallocate(tfp)
    if(allocated(cfp)) deallocate(cfp)
    if(allocated(btotal)) deallocate(btotal)
    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)
    if(allocated(psin)) deallocate(psin)
    if(allocated(psinr)) deallocate(psinr)
    if(allocated(fpol)) deallocate(fpol)
    if(allocated(rhopnr)) deallocate(rhopnr)
    if(allocated(qpsi)) deallocate(qpsi)

  end subroutine dealloc_equilvec


  subroutine alloc_bnd(ier)
    implicit none
    integer, intent(out) :: ier
    
    if(nlim.lt.0.or.nbbbs.lt.0) then
      ier = -1
      return
    end if
    
    call dealloc_bnd
    allocate(rlim(nlim),zlim(nlim),rbbbs(nbbbs),zbbbs(nbbbs), &
     stat=ier)
    if (ier/=0) call dealloc_bnd
  end subroutine alloc_bnd

  subroutine dealloc_bnd
    implicit none
    if(allocated(rlim))   deallocate(rlim)
    if(allocated(zlim))   deallocate(zlim)
    if(allocated(rbbbs))   deallocate(rbbbs)
    if(allocated(zbbbs))   deallocate(zbbbs)
  end subroutine dealloc_bnd


  subroutine alloc_lim(ier)
    implicit none
    integer, intent(out) :: ier
    
    if(nlim.lt.0) then
      ier = -1
      return
    end if
    
    call dealloc_lim
    allocate(rlim(nlim),zlim(nlim), &
     stat=ier)
    if (ier/=0) call dealloc_lim
  end subroutine alloc_lim

  subroutine dealloc_lim
    implicit none
    if(allocated(rlim))   deallocate(rlim)
    if(allocated(zlim))   deallocate(zlim)
  end subroutine dealloc_lim

  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,lw10)
      end if
      if (present(dpsidz)) then
        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01,lw01)
      end if
      if (present(ddpsidrr)) then
        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20,lw20)
      end if
      if (present(ddpsidzz)) then
        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02,lw02)
      end if
      if (present(ddpsidrz)) then
        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11,lw11)
      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,lw)
    use dierckx, only : fpbisp
    implicit none
! local constants
    integer, parameter :: liwrk=2,nrs=1,nzs=1
! arguments
    integer :: nur,nuz,lw
    real(wp_) :: rpsim,zpsim,derpsi
    real(wp_), dimension(lw) :: cc
! local variables
    integer :: iwr,iwz
    integer, dimension(liwrk) :: iwrk
    real(wp_), dimension(1) :: rrs,zzs,ffspl

    rrs(1)=rpsim
    zzs(1)=zpsim
    iwr=1+(nr-nur-4)*(nz-nuz-4)
    iwz=iwr+4-nur
    call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,3-nur,3-nuz, &
         rrs,nrs,zzs,nzs,ffspl,cc(iwr),cc(iwz),iwrk(1),iwrk(2))
    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.gt.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
  
end module interp_eqprof