module pec
  use const_and_precisions, only : wp_

  implicit none

  real(wp_), dimension(:), allocatable, save :: rhop_tab,rhot_tab
  real(wp_), dimension(:), allocatable, save :: rtabpsi1
  real(wp_), dimension(:), allocatable, save :: dvol,darea
  real(wp_), dimension(:), allocatable, save :: ratjav,ratjbv,ratjplv

contains

  subroutine pec_init(params, equil, rt_in)
    use gray_params,  only : output_parameters, RHO_POL
    use gray_equil,   only : abstract_equil
    use magsurf_data, only : fluxval

    ! subroutine arguments
    type(output_parameters), intent(in)           :: params
    class(abstract_equil),   intent(in)           :: equil
    real(wp_),               intent(in), optional :: rt_in(params%nrho)

    ! local variables
    integer :: it
    real(wp_) :: drt, rt, rt1, rhop1
    real(wp_) :: ratjai, ratjbi, ratjpli
    real(wp_) :: voli0, voli1, areai0, areai1

    associate(n => params%nrho)

    ! rt_in present: read input grid
    ! else: build equidistant grid dimension n

    ! ipec=0  ρ_t grid
    ! ipec=1  ρ_p grid
    call dealloc_pec
    allocate(rhop_tab(n), rhot_tab(n), rtabpsi1(0:n))
    allocate(dvol(n), darea(n), ratjav(n), ratjbv(n), ratjplv(n))

    voli0  = 0
    areai0 = 0
    rtabpsi1(0) = 0

    do it = 1, n
      if(present(rt_in)) then
        ! read radial grid from input
        rt  = rt_in(it)
        if (it < n) then
          drt = rt_in(it+1)-rt
        end if
      else
        ! build equidistant radial grid
        drt = 1/dble(n - 1)
        rt  = (it - 1)*drt
      end if
      ! radial coordinate of i-(i+1) interval mid point
      if (it < n) then
        rt1 = rt + drt/2
      else
        rt1 = 1
      end if
      if (params%ipec == RHO_POL) then
        rhop_tab(it) = rt
        rhot_tab(it) = equil%pol2tor(rt)
        rhop1 = rt1
      else
        rhot_tab(it) = rt
        rhop_tab(it) = equil%tor2pol(rt)
        rhop1 = equil%tor2pol(rt1)
      end if

      ! psi grid at mid points, size n+1, for use in pec_tab
      rtabpsi1(it) = rhop1**2

      call fluxval(rhop1,area=areai1,vol=voli1)
      dvol(it)  = abs(voli1  - voli0)
      darea(it) = abs(areai1 - areai0)
      voli0  = voli1
      areai0 = areai1

      call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi,ratjpl=ratjpli)
      ratjav(it)  = ratjai
      ratjbv(it)  = ratjbi
      ratjplv(it) = ratjpli
    end do

    end associate
  end subroutine pec_init


  subroutine spec(psjki, ppabs, ccci, iiv, pabs, &
                  currt, dpdv, ajphiv, ajcd, pins, currins)
    ! subroutine arguments
    real(wp_), dimension(:, :), intent(in)  :: psjki, ppabs, ccci
    integer,                    intent(in)  :: iiv(:)
    real(wp_),                  intent(in)  :: pabs, currt
    real(wp_), dimension(:),    intent(out) :: dpdv, ajphiv, ajcd, pins, currins

    ! local variables
    integer   :: i, ii, jk
    real(wp_) :: spds, sccs, facpds, facjs
    real(wp_), dimension(size(psjki, dim=2)) :: xxi, ypt, yamp
    real(wp_), dimension(size(dpdv))         :: wdpdv, wajphiv

    ! calculation of dP and dI over radial grid
    dpdv   = 0
    ajphiv = 0
    do jk = 1, size(iiv)
      ii = iiv(jk)
      if (ii < size(psjki, dim=2)) then
        if (psjki(jk,ii+1) /= 0) ii = ii + 1  !!! CHECK
      end if
      xxi  = 0
      ypt  = 0
      yamp = 0
      do i = 1, ii
        xxi(i)=abs(psjki(jk,i))
        if(xxi(i) <= 1) then
          ypt(i) = ppabs(jk,i)
          yamp(i) = ccci(jk,i)
        end if
      end do
      call pec_tab(xxi, ypt, yamp, ii, rtabpsi1, wdpdv, wajphiv)
      dpdv   = dpdv   + wdpdv
      ajphiv = ajphiv + wajphiv
    end do

    ! here dpdv is still dP (not divided yet by dV)
    ! here ajphiv is still dI (not divided yet by dA)
    spds = 0
    sccs = 0
    do i = 1, size(dpdv)
      spds = spds + dpdv(i)
      sccs = sccs + ajphiv(i)
      pins(i) = spds
      currins(i) = sccs
    end do

    facpds = 1
    facjs = 1
    if(spds > 0) facpds = pabs / spds
    if(sccs /= 0) facjs = currt / sccs

    dpdv    = facpds * (dpdv / dvol)
    ajphiv  = facjs * (ajphiv / darea)
    ajcd    = ajphiv * ratjbv
    pins    = facpds * pins
    currins = facjs * currins

    ! now dpdv is  dP/dV [MW/m^3]
    ! now ajphiv is J_phi=dI/dA [MA/m^2]
  end subroutine spec


  subroutine pec_tab(xxi, ypt, yamp, ii, xtab1, wdpdv, wajphiv)
    ! Power and current projected on ψ grid - mid points
    use utils, only : locatex, intlin

    ! subroutine arguments
    real(wp_), dimension(:),                intent(in)  :: xxi, ypt, yamp
    integer,                                intent(in)  :: ii
    real(wp_), dimension(0:),               intent(in)  :: xtab1
    real(wp_), dimension(ubound(xtab1, 1)), intent(out) :: wdpdv, wajphiv

    ! local variables
    integer   :: isev(21)
    real(wp_) :: ppa1, ppa2, cci1, cci2, dppa, didst, rt1
    integer   :: i, is, ise0, idecr, iise0, iise, iis, iis1
    integer   :: ind1, ind2, iind, ind, indi, itb1

    isev  =  0
    ise0  =  0
    idecr = -1
    is    =  1
    wdpdv   = 0
    wajphiv = 0
    do i=1,ii
      if(ise0 == 0) then
        if(xxi(i) < 1) then
          ise0     = i
          isev(is) = max(i - 1, 1)
          is       = is + 1
        end if
      else
        if (idecr == -1) then
          if(xxi(i) > xxi(i-1)) then
            isev(is) = i  - 1
            is       = is + 1
            idecr    = 1
          end if
        else
          if(xxi(i) > 1) exit
          if(xxi(i) < xxi(i-1)) then
            isev(is) = i  - 1
            is       = is + 1
            idecr    = -1
          end if
        end if
      end if
    end do

    isev(is) = i-1
    ppa1 = 0
    cci1 = 0

    do iis = 1, is-1
      iis1  = iis + 1
      iise0 = isev(iis)
      iise  = isev(iis1)
      if (mod(iis,2) /= 0) then
        idecr = -1
        ind1  = size(wdpdv)
        ind2  =  2
        iind  = -1
      else
        idecr =  1
        ind1  =  1
        ind2  = size(wdpdv)
        iind  =  1
      end if
      do ind = ind1, ind2, iind
        indi = ind
        if (idecr == -1) indi = ind - 1
        rt1 = xtab1(indi)
        call locatex(xxi,iise,iise0,iise,rt1,itb1)
        if (itb1 >= iise0 .and. itb1 <= iise) then
          if(itb1 == iise) then
            ppa2=ypt(itb1)
            cci2=yamp(itb1)
          else
            call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2)
            call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2)
          end if
          dppa  = ppa2 - ppa1
          didst = cci2 - cci1
          wdpdv(ind)   = wdpdv(ind)   + dppa
          wajphiv(ind) = wajphiv(ind) + didst
          ppa1 = ppa2
          cci1 = cci2
        end if
      end do
    end do
  end subroutine pec_tab


  subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, &
        rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip,                 &
        rhotp,  drhotp,  rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx)

    ! Radial average values over power and current density profile
    use const_and_precisions, only : pi, comp_eps
    use gray_equil,           only : abstract_equil
    use magsurf_data,         only : fluxval

    class(abstract_equil), intent(in)  :: equil
    real(wp_),             intent(in)  :: pabs,currt
    real(wp_),             intent(in)  :: rhot_tab(:)
    real(wp_),             intent(in)  :: dpdv(:), ajphiv(:)
    real(wp_),             intent(out) :: rhotpav, drhotpav, dpdvp
    real(wp_),             intent(out) :: rhotjava, drhotjava, ajphip
    real(wp_),             intent(out) :: rhotp, drhotp, dpdvmx
    real(wp_),             intent(out) :: rhotjfi, drhotjfi, ajmxfi
    real(wp_),             intent(out) :: ratjamx, ratjbmx

    real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav
    real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava

    rhotpav   = 0
    rhot2pav  = 0
    rhotjav   = 0
    rhotjava  = 0
    rhot2java = 0

    if (pabs > 0) then
      rhotpav  = sum(rhot_tab   *dpdv*dvol)/pabs
      rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs
    end if

    if (abs(currt) > 0) then
      rhotjav = sum(rhot_tab*ajphiv*darea)/currt
    end if
    sccsa = sum(abs(ajphiv)*darea)
    if (sccsa > 0) then
      rhotjava  = sum(rhot_tab   *abs(ajphiv)*darea)/sccsa
      rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa
    end if

    ! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
    drhotpav  = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps)))
    drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps)))

    rhoppav  = equil%tor2pol(rhotpav)
    rhopjava = equil%tor2pol(rhotjava)

    if (pabs > 0) then
      call fluxval(rhoppav,dvdrhot=dvdrhotav)
      dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav)
      call profwidth(rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
    else
      dpdvp  = 0
      rhotp  = 0
      dpdvmx = 0
      drhotp = 0
    end if

    if (sccsa > 0) then
      call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx,ratjb=ratjbmx, &
                   ratjpl=ratjplmx)
      ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava)
      call profwidth(rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
    else
      ajphip   = 0
      rhotjfi  = 0
      ajmxfi   = 0
      drhotjfi = 0
      ratjamx  = 0
      ratjbmx  = 0
    end if
  end subroutine postproc_profiles


  subroutine profwidth(xx,yy,xpk,ypk,dxxe)
    use const_and_precisions, only : emn1
    use utils, only : locatex, locate, intlin, vmaxmini

    ! subroutine arguments
    real(wp_), intent(in)  :: xx(:), yy(:)
    real(wp_), intent(out) :: xpk, ypk, dxxe

    ! local variables
    integer :: imn,imx,ipk,ie
    real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
    real(wp_) :: ypkp,ypkm

    call vmaxmini(yy,size(yy),ymn,ymx,imn,imx)
    ypk = 0
    xmx = xx(imx)
    xmn = xx(imn)
    if (abs(ymx) > abs(ymn)) then
      ipk  = imx
      ypkp = ymx
      xpkp = xmx
      if(abs(ymn/ymx) < 1.0e-2_wp_) ymn = 0.0_wp_
      ypkm = ymn
      xpkm = xmn
    else
      ipk  = imn
      ypkp = ymn
      xpkp = xmn
      if(abs(ymx/ymn) < 1.0e-2_wp_) ymx = 0.0_wp_
      ypkm = ymx
      xpkm = xmx
    end if
    if(xpkp > 0) then
      xpk = xpkp
      ypk = ypkp
      yye = ypk*emn1
      call locatex(yy,size(yy),1,ipk,yye,ie)
      if(ie > 0 .and. ie < size(yy)) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1)
      else
        rte1 = 0
      end if
      call locatex(yy,size(yy),ipk,size(yy),yye,ie)
      if(ie > 0 .and. ie < size(yy)) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
      else
        rte2 = 0
      end if
    else
      ipk=2
      xpk=xx(2)
      ypk=yy(2)
      rte1=0.0_wp_
      yye=ypk*emn1
      call locate(yy,size(yy),yye,ie)
      if(ie > 0 .and. ie < size(yy)) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
      else
        rte2 = 0
      end if
    end if
    dxxe = rte2 - rte1
    if(ymx /= 0 .and. ymn /= 0) dxxe = -dxxe
  end subroutine profwidth


  subroutine dealloc_pec
    if (allocated(rhop_tab)) deallocate(rhop_tab)
    if (allocated(rhot_tab)) deallocate(rhot_tab)
    if (allocated(rtabpsi1)) deallocate(rtabpsi1)
    if (allocated(dvol)) deallocate(dvol)
    if (allocated(darea)) deallocate(darea)
    if (allocated(ratjav)) deallocate(ratjav)
    if (allocated(ratjbv)) deallocate(ratjbv)
    if (allocated(ratjplv)) deallocate(ratjplv)
  end subroutine dealloc_pec

end module pec