module pec
  use const_and_precisions, only : wp_,zero,one

  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(ipec,rt_in)
    use equilibrium, only : frhotor,frhopol
    use gray_params, only : nnd
    use magsurf_data, only : fluxval
! arguments
    integer, intent(in) :: ipec
    real(wp_), dimension(nnd), intent(in), optional ::  rt_in
! local variables
    integer :: it
    real(wp_) :: drt,rt,rt1,rhop1
    real(wp_) :: ratjai,ratjbi,ratjpli
    real(wp_) :: voli0,voli1,areai0,areai1

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

!     ipec=0  rho_tor  grid
!     ipec=1  rho_pol  grid
    call dealloc_pec
    allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), &
      ratjav(nnd),ratjbv(nnd),ratjplv(nnd))
    
    voli0  = zero
    areai0 = zero
    rtabpsi1(0) = zero

    do it=1,nnd
      if(present(rt_in)) then 
!     read radial grid from input
        rt  = rt_in(it)
        if(it<nnd) then
          drt = rt_in(it+1)-rt
        end if
      else
!     build equidistant radial grid       
        drt = one/dble(nnd-1)
        rt  = dble(it-1)*drt
      end if
!     radial coordinate of i-(i+1) interval mid point
      if(it < nnd) then
        rt1 = rt + drt/2.0_wp_
      else
        rt1 = one
      end if
      if (ipec == 1) then
        rhop_tab(it) = rt
        rhot_tab(it) = frhotor(rt)
        rhop1 = rt1
      else
        rhot_tab(it) = rt
        rhop_tab(it) = frhopol(rt)
        rhop1 = frhopol(rt1)
      end if  
!     psi grid at mid points, dimension nnd+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 subroutine pec_init


  subroutine spec(psjki,ppabs,ccci,iiv,pabs,currt,dpdv,ajphiv,ajcd,pins,currins)
    use gray_params, only : nnd
    use beamdata, only : nray,nstep
! arguments
    real(wp_), dimension(nray,nstep), intent(in) :: psjki,ppabs,ccci
    integer, dimension(nray), intent(in) :: iiv
    real(wp_), intent(in) :: pabs,currt
    real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv,ajcd,pins,currins
! local variables
    integer :: i,ii,jk
    real(wp_) :: spds,sccs,facpds,facjs
    real(wp_), dimension(nstep):: xxi,ypt,yamp
    real(wp_), dimension(nnd) :: wdpdv,wajphiv
    
! calculation of dP and dI over radial grid 
    dpdv=zero
    ajphiv=zero
    do jk=1,nray
      ii=iiv(jk)
      if (ii < nstep ) then
        if(psjki(jk,ii+1) /= zero) ii=ii+1  !!! CHECK
      end if
      xxi=zero
      ypt=zero
      yamp=zero
      do i=1,ii
        xxi(i)=abs(psjki(jk,i))
        if(xxi(i) <= one) 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=zero
    sccs=zero
    do i=1,nnd
      spds=spds+dpdv(i)
      sccs=sccs+ajphiv(i)
      pins(i)=spds
      currins(i)=sccs
    end do
   
    facpds=one
    facjs=one
    if(spds > zero) facpds=pabs/spds
    if(sccs /= zero) 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 psi grid - mid points
    use const_and_precisions, only : wp_,one,zero
    use gray_params, only : nnd
    use utils, only : locatex,intlin
! arguments
    integer, intent(in) :: ii
    real(wp_), dimension(ii), intent(in) :: xxi,ypt,yamp
    real(wp_), dimension(0:nnd), intent(in) :: xtab1
    real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv 
! local variables
    integer, parameter :: llmx = 21
    integer, dimension(llmx) ::isev
    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   = zero
    wajphiv = zero
    do i=1,ii
      if(ise0 == 0) then
        if(xxi(i) < one) 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) > one) 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 = zero
    cci1 = zero
    
    do iis=1,is-1
      iis1  = iis + 1
      iise0 = isev(iis)
      iise  = isev(iis1)
      if (mod(iis,2) /= 0) then
        idecr = -1
        ind1  = nnd
        ind2  =  2
        iind  = -1
      else
        idecr =  1
        ind1  =  1
        ind2  = nnd
        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(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_params, only : nnd
    use equilibrium, only : frhopol
    use magsurf_data, only : fluxval
    real(wp_), intent(in) ::  pabs,currt
    real(wp_), dimension(nnd), intent(in) :: rhot_tab
    real(wp_), dimension(nnd), 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=zero
    rhot2pav=zero
    rhotjav=zero
    rhotjava=zero
    rhot2java=zero
         
    if (pabs > zero) then
      rhotpav  = sum(rhot_tab   *dpdv*dvol)/pabs
      rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs
    end if
    
    if (abs(currt) > zero) then
      rhotjav = sum(rhot_tab*ajphiv*darea)/currt
    end if
    sccsa = sum(abs(ajphiv)*darea)
    if (sccsa > zero) 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  = frhopol(rhotpav)
    rhopjava = frhopol(rhotjava)
    
    if (pabs > zero) then
      call fluxval(rhoppav,dvdrhot=dvdrhotav)
      dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav)
      call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
    else
      dpdvp  = zero
      rhotp  = zero
      dpdvmx = zero
      drhotp = zero
    end if

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



  subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe)
    use const_and_precisions, only : wp_,emn1
    use utils, only : locatex, locate, intlin, vmaxmini
! arguments
    integer :: nd
    real(wp_), dimension(nd) :: 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,nd,ymn,ymx,imn,imx)
    ypk = zero
    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 > zero) then
      xpk = xpkp
      ypk = ypkp
      yye = ypk*emn1
      call locatex(yy,nd,1,ipk,yye,ie)
      if(ie > 0 .and. ie < nd) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1)
      else
        rte1 = zero
      end if
      call locatex(yy,nd,ipk,nd,yye,ie)
      if(ie > 0 .and. ie < nd) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
      else
        rte2 = zero
      end if
    else
      ipk=2
      xpk=xx(2)
      ypk=yy(2)
      rte1=0.0_wp_
      yye=ypk*emn1
      call locate(yy,nd,yye,ie)
      if(ie > 0 .and. ie < nd) then
        call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
      else
        rte2 = zero
      end if
    end if
    dxxe = rte2 - rte1
    if(ymx /= zero .and. ymn /= zero) 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