module coreprofiles
  use const_and_precisions, only : wp_,zero,one
  implicit none

  INTEGER, SAVE :: npp,nsfd
  REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad
  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
  REAL(wp_), SAVE :: dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan

contains

  subroutine density(psin,dens,ddens)
    use gray_params, only : iprof
    use dierckx, only : splev,splder
    implicit none
! arguments
    real(wp_), intent(in) :: psin
    real(wp_), intent(out) :: dens,ddens
! local variables
    integer, parameter :: nn=3, nn1=nn+1, nn2=nn+2
    integer :: ier,nu
    real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
    real(wp_), dimension(1) :: xxs,ffs
    real(wp_), dimension(npp+4) :: wrkfd

!
!     computation of density [10^19 m^-3] and derivative wrt psi
!
    dens=zero
    ddens=zero
    if((psin >= psdbnd).or.(psin < zero)) return
!
    if(iprof == 0) then
      if(psin > one) return
      profd=(one-psin**aln1)**aln2
      dens=dens0*profd
      dprofd=-aln1*aln2*psin**(aln1-one)           &
              *(one-psin**aln1)**(aln2-one)
      ddens=dens0*dprofd
    else
      if(psin > psnpp) then

!     smooth interpolation for psnpp < psi < psdbnd
!     dens = fp * fh
!     fp: parabola matched at psi=psnpp with given profile density
!     fh=(1-t)^3(1+3t+6t^2) is a smoothing function: 
!     fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 
!
        dpsib=psin-psnpp
        fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
        dfp=ddenpp+dpsib*d2denpp
        tt=dpsib/(psdbnd-psnpp)
        fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2)
        dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp)
        dens=fp*fh
        ddens=dfp*fh+fp*dfh
      else
        xxs(1)=psin
        ier=0
        call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
        dens=ffs(1)
        nu=1
        ier=0
        call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
        ddens=ffs(1)
        if(abs(dens) < 1.0e-10_wp_) dens=zero
      end if
      if(dens < zero) print*,'psin = ',psin,':    DENSITY NEGATIVE ne=',dens
!      if(dens < zero) then
!        dens=zero
!        ddens=zero
!      end if  
    end if
  end subroutine density

  function temp(psin)
    use const_and_precisions, only : wp_,zero,one
    use gray_params, only : iprof
    use utils, only : locate
    use simplespline, only :spli
    implicit none
! arguments
    real(wp_), intent(in) :: psin
    real(wp_) :: temp
! local variables
    integer :: k
    real(wp_) :: proft,dps

    temp=zero
    if((psin >= one).or.(psin < zero)) return
    if(iprof == 0) then
      proft=(1.0_wp_-psin**alt1)**alt2
      temp=(te0-dte0)*proft+dte0
    else
      call locate(psrad,npp,psin,k)
      k=max(1,min(k,npp-1))
      dps=psin-psrad(k)
      temp=spli(ct,npp,k,dps)
    endif
  end function temp

  function fzeff(psin)
    use const_and_precisions, only : wp_,zero,one
    use gray_params, only : iprof
    use utils, only : locate
    use simplespline, only :spli
    implicit none
! arguments
    real(wp_), intent(in) :: psin
    real(wp_) :: fzeff
! local variables
    integer :: k
    real(wp_) :: dps

    fzeff=one
    if((psin >= one).or.(psin < zero)) return
    if(iprof == 0) then
      fzeff=zeffan
    else
      call locate(psrad,npp,psin,k)
      k=max(1,min(k,npp-1))
      dps=psin-psrad(k)
      fzeff=spli(cz,npp,k,dps)
    endif
  end function fzeff
  
  subroutine read_profiles(filenm,psin,te,ne,zeff,unit)
    use utils, only : get_free_unit
    implicit none
! arguments
    character(len=*), intent(in) :: filenm
    real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff
    integer, optional, intent(in) :: unit
! local variables
    integer :: u, i, n
    
    if (present(unit)) then
      u=unit
    else
      u=get_free_unit()
    end if
    open(file=trim(filenm),status='old',action='read',unit=u)
    read(u,*) n
    if(allocated(psin)) deallocate(psin)
    if(allocated(te)) deallocate(te)
    if(allocated(ne)) deallocate(ne)
    if(allocated(zeff)) deallocate(zeff)
    allocate(psin(n),te(n),ne(n),zeff(n))
    do i=1,n
      read(u,*) psin(i),te(i),ne(i),zeff(i)
    end do
    psin(1)=max(psin(1),zero)
    close(u)
  end subroutine read_profiles

  subroutine read_profiles_an(filenm,te,ne,zeff,unit)
    use utils, only : get_free_unit
     implicit none
! arguments
    character(len=*), intent(in) :: filenm
    real(wp_), dimension(:), allocatable, intent(out) :: te,ne,zeff
    integer, optional, intent(in) :: unit
! local variables
    integer :: u
        
    if (present(unit)) then
      u=unit
    else
      u=get_free_unit()
    end if

    if(allocated(te)) deallocate(te)
    if(allocated(ne)) deallocate(ne)
    if(allocated(zeff)) deallocate(zeff)
    allocate(te(4),ne(3),zeff(1))

    open(file=trim(filenm),status='old',action='read',unit=u)
    read(u,*) ne(1:3) ! dens0,aln1,aln2
    read(u,*) te(1:4) ! te0,dte0,alt1,alt2
    read(u,*) zeff(1) ! zeffan
    close(u)
  end subroutine read_profiles_an

  subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal,iprof)
    implicit none
! arguments
    real(wp_), dimension(:), intent(inout) :: te,ne
    real(wp_), intent(in) :: tfact,nfact,bfact
    integer, intent(in) :: iscal,iprof
! local variables
    real(wp_) :: aat,aan,ffact
    integer :: lastte,lastne

    if (iscal==0) then
      aat=2.0_wp_/3.0_wp_
      aan=4.0_wp_/3.0_wp_
    else
      aat=1.0_wp_
      aan=1.0_wp_
    end if
    if(iscal==2) then
      ffact=1.0_wp_
    else
      ffact=bfact
    end if
    if (iprof==0) then
      lastte=2
      lastne=1
    else
      lastte=size(te)
      lastne=size(ne)
    end if
    te(1:lastte)=te(1:lastte)*ffact**aat*tfact
    ne(1:lastne)=ne(1:lastne)*ffact**aan*nfact
  end subroutine tene_scal

  subroutine set_prfspl(psin,te,ne,zeff,ssplne,psdbndmx)
    use simplespline, only : difcs
    use dierckx, only : curfit, splev, splder
    implicit none
! arguments
    real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff
    real(wp_), intent(in) :: ssplne,psdbndmx
! local variables
    integer, parameter :: iopt=0, kspl=3
    integer :: n, npest, lwrkf, ier
    real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2,ssplne_loc
    real(wp_), dimension(:), allocatable :: wf, wrkf
    integer, dimension(:), allocatable :: iwrkf
    real(wp_), dimension(1) :: dedge,ddedge,d2dedge

    n=size(psin)
    npest=n+4
    lwrkf=n*4+npest*16
    allocate(wrkf(lwrkf),iwrkf(npest),wf(n))

    ssplne_loc=ssplne

!   if necessary, reallocate spline arrays
    if(.not.allocated(psrad)) then
      allocate(psrad(n),ct(n,4),cz(n,4))
    else
      if(size(psrad)<n) then
        deallocate(psrad,ct,cz)
        allocate(psrad(n),ct(n,4),cz(n,4))
      end if
    end if
    if(.not.allocated(cfn)) then
      allocate(tfn(npest),cfn(npest))
    else
      if(size(cfn)<npest) then
        deallocate(tfn,cfn)
        allocate(tfn(npest),cfn(npest))
      end if
    end if

!   spline approximation of  temperature and Zeff
    call difcs(psin,te,  n,iopt,ct,ier)
    call difcs(psin,zeff,n,iopt,cz,ier)
    psrad=psin
    npp=n

!   spline approximation of density
    xb=zero
    xe=psin(n)
    wf(:)=one
    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
       nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
! if ier=-1 data are re-fitted using sspl=0
    if(ier==-1) then
        write(*,*) 'density curfit: ier=-1. Re-fitting with interpolating spline'
        ssplne_loc=0.0_wp_
        call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
           nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
    end if

!   compute polinomial extrapolation matching the spline boundary up to the
!   2nd order derivative, extending the profile up to psi=psdbnd where
!   ne=ne'=ne''=0
!   spline value and derivatives at the edge
    call splev(tfn,nsfd,cfn,kspl,psin(n:n),dedge(1:1),1,ier)
    call splder(tfn,nsfd,cfn,kspl,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
    call splder(tfn,nsfd,cfn,kspl,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
!   determination of the boundary
    psdbnd=psdbndmx

    psnpp=psin(n)
    denpp=dedge(1)
    ddenpp=ddedge(1)
    d2denpp=d2dedge(1)

    delta2=(ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp
    xnv=psnpp-ddenpp/d2denpp
    if(delta2 < zero) then
!       if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
    else
      xxm=xnv-sqrt(delta2)
      xxp=xnv+sqrt(delta2)
      if(xxm > psnpp) then
        psdbnd=min(psdbnd,xxm)
      else if (xxp > psnpp) then
        psdbnd=min(psdbnd,xxp)
      end if    
      write(*,*) 'density psdbnd=',psdbnd
    end if

    deallocate(iwrkf,wrkf,wf)
  end subroutine set_prfspl

  subroutine unset_prfspl
    implicit none

    if(allocated(psrad)) deallocate(psrad)
    if(allocated(ct)) deallocate(ct)
    if(allocated(cz)) deallocate(cz)
    if(allocated(tfn)) deallocate(tfn)
    if(allocated(cfn)) deallocate(cfn)
  end subroutine unset_prfspl

  subroutine set_prfan(te,ne,zeff)
    implicit none
    REAL(wp_), dimension(:), intent(in) :: te,ne,zeff

    te0=te(1)
    dte0=te(2)
    alt1=te(3)
    alt2=te(4)
    dens0=ne(1)
    aln1=ne(2)
    aln2=ne(3)
    zeffan=zeff(1)

    psdbnd=1.0_wp_
  end subroutine set_prfan
  
end module coreprofiles