!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MODIFIED BY Lei Shi to have unified interface for 1D, 2D, and 3D quadratic spline
! Date: Nov. 19, 2016
! Modified by Lei Shi to have generic 1D smooth function
! Date: Jan 24, 2017
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!###### MODULE SMOOTH_FUNCTION ###########################
! Contains generic functions for smoothing

module smooth_function
use precision, only: lk
use global_parameters, only: gtcout
implicit none

private
public:: smooth_1d
contains

! Generic 1D smooth function
! Proxy to different types of 1D smoothing
! Arguments:
!     array: the data needs smoothing, will be smoothen IN PLACE
!     n: length of array
!     boundary_condition: optional, default to 0
!       chosen boundary condition:
!              0 for normal fix boundary, 
!              1 for periodic boundary.
!     smooth_type: optional default to 0
!       smoothing algorithm for choose:
!              0 for "121" smoothing,
!              1 for "14641" smoothing.
subroutine smooth_1d(array, n, boundary_condition, smooth_type)

  real(lk), intent(inout):: array(:)
  integer, intent(in):: n
  integer, intent(in), optional:: boundary_condition
  integer, intent(in), optional:: smooth_type
  integer :: bc, st ! local copies of the boundary_condition and 
                    ! smooth_type arguments 

  if (present(smooth_type)) then
    st = smooth_type
  else
    st = 0 ! default to use 1-2-1 smoothing function
  endif

  if (present(boundary_condition)) then
    bc = boundary_condition
  else
    bc = 0 ! default to use normal fix boundary
  endif

! Now switch through the smoothing types and boundary conditions

  select case (st)
    case (0)
      if (bc==0) then
        call smooth_fix_121(array, n)
      else if (bc==1) then
        call smooth_periodic_121(array, n)
      else
        print *, "Invalid BC option in smooth_1d: ", bc
      endif
    case (1)
      if (bc==0) then
        call smooth_fix_14641(array, n)
      else if (bc==1) then
        call smooth_periodic_14641(array, n)
      else
        print *, "Invalid BC option in smooth_1d: ", bc
      endif
    case default
      print *,"Invalid smoothing type in smooth_1d: ", st
  end select

end subroutine smooth_1d
! subroutine to smooth the fix boundary data using 121 formula
subroutine smooth_fix_121(array, n)
  real(lk), intent(inout) :: array(:)
  integer, intent (in) :: n
  integer :: i
  real(lk) :: temp(n)
  temp(:) = array(:)
  do i=2, n-1
     array(i) = 0.25_lk *(temp(i-1) + 2.0_lk*temp(i) + temp(i+1))
  enddo
end subroutine smooth_fix_121

! subroutine to smooth the periodic data using 121 formula
subroutine smooth_periodic_121(array, n)
  real(lk), intent(inout) :: array(:)
  integer, intent (in) :: n
  integer :: i
  real(lk) :: temp(n)
! if input array is not periodic, the array(n) will be assigned by array(1)
  if (ABS(array(1)-array(n))> 2*SPACING(array(1))) then ! test if array(1) equals array(n)
    write(gtcout,*) "SMOOTH_1D WARNING: periodic boundary condition &
        &applied to non-periodic data."//NEW_LINE('a')//"array(n)=array(1) is enforced."
    array(n)=array(1)
  endif
  temp(:) = array(:)
  do i=2, n-1
     array(i) = 0.25_lk *(temp(i-1) + 2.0_lk*temp(i) + temp(i+1))
  enddo
! smooth end points using periodic boundary condition
  array(1) = 0.25_lk *(temp(n-1) + 2.0_lk *temp(1) + temp(2))
  array(n) = array(1)
end subroutine smooth_periodic_121

! subroutine to smooth the fix boundary data using 14641 formula
subroutine smooth_fix_14641(array, n)
  real(lk), intent(inout) :: array(:)
  integer, intent (in) :: n
  integer :: i
  real(lk) :: temp(n), coeff
  temp(:) = array(:)
  coeff = 1.0_lk/12.0_lk
  do i=3, n-2
     array(i) = coeff * (-(temp(i-2)+temp(i+2)) + 4.0_lk*(temp(i-1)+temp(i+1)) + 6.0_lk*temp(i))
  enddo
end subroutine smooth_fix_14641

! subroutine to smooth the periodic data using 14641 formula
subroutine smooth_periodic_14641(array, n)
  real(lk), intent(inout) :: array(:)
  integer, intent (in) :: n
  integer :: i
  real(lk) :: temp(n), coeff
! if input array is not periodic, the array(n) will be assigned by array(1)
  if (ABS(array(1)-array(n))> 2*SPACING(array(1))) then ! test if array(1) equals array(n)
    write(gtcout,*) "SMOOTH_1D WARNING: periodic boundary condition &
        &applied to non-periodic data."//NEW_LINE('a')//"array(n)=array(1) is enforced."
    array(n)=array(1)
  endif
  temp(:) = array(:)
  coeff = 1.0_lk/12.0_lk
  do i=3, n-2
     array(i) = coeff * (-(temp(i-2)+temp(i+2)) + 4.0_lk*(temp(i-1)+temp(i+1)) + 6.0_lk*temp(i))
  enddo
! smooth end points using periodic boundary condition
  array(1) = coeff * (-(temp(n-2)+temp(3)) + 4.0_lk*(temp(n-1)+temp(2)) + 6.0_lk*temp(1))
  array(n) = array(1)
  array(2) = coeff * (-(temp(n-1)+temp(4)) + 4.0_lk*(temp(1)+temp(3)) + 6.0_lk*temp(2))
  array(n-1) = coeff * (-(temp(n-3)+temp(2)) + 4.0_lk*(temp(n-2)+temp(1)) + 6.0_lk*temp(n-1))
end subroutine smooth_periodic_14641

end module smooth_function


!##############################################################################
!####  MODULE SPLINE FUNCTIONS
!##############################################################################
module spline_function
use precision, only: lk
use global_parameters, only: gtcout
use smooth_function, only: smooth_1d
implicit none

public
private:: spline0, dspline0, spline1, dspline1, spline_periodic, dspline_periodic, &
      construct_spline0, construct_spline1, &
     construct_spline_periodic
contains


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic 1D spline evaluations
! Evaluate f(x) at given x 
! Boundary condition needs to be provided by programmer
! BC 0: use spline0
! BC 1: use spline1
! BC 2: use spline_periodic

real(lk) function spline1d(x, deriv, nx, delx, y, bcx)
  integer, intent(in) :: nx, deriv, bcx
  real(lk), intent(in) :: x, delx, y(3,nx)
  integer i
  real(lk) dx
  
  if (bcx == 0) then
     if (deriv==0) spline1d = spline0(x, nx, delx, y)
     if (deriv==1) spline1d = dspline0(x, nx, delx, y)
  else if (bcx == 1) then
     if (deriv==0) spline1d = spline1(x, nx, delx, y)
     if (deriv==1) spline1d = dspline1(x, nx, delx, y)
  else if (bcx == 2) then
     if (deriv==0) spline1d = spline_periodic(x, nx, delx, y)
     if (deriv==1) spline1d = dspline_periodic(x, nx, delx, y)
  else
     write(gtcout,*) 'Spline1d Error: Invalid boundary condition bcx=', bcx
     stop
  endif

end function spline1d


! 1D spline for radial profiles
real(lk) function spline0(pdum,nsp,delx,y)
  integer nsp,i
  real(lk) pdum,y(3,nsp),delx,dpx

  i=max(1,min(nsp-1,ceiling(pdum/delx)))
  dpx=pdum-delx*real(i-1)
  spline0=y(1,i)+dpx*y(2,i)+dpx*dpx*y(3,i)

end function spline0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! derivative of 1D spline function
real(lk) function dspline0(pdum,nsp,delx,y)
  integer nsp,i
  real(lk) pdum,y(3,nsp),delx,dpx

  i=max(1,min(nsp-1,ceiling(pdum/delx)))
  dpx=pdum-delx*real(i-1)
  dspline0=y(2,i)+2.0*dpx*y(3,i)

end function dspline0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! 1D spline with first point being linear function y=sqrt(x)
real(lk) function spline1(pdum,nsp,delx,y)
  integer nsp,i
  real(lk) pdum,y(3,nsp),delx,dpx,dp2

  i=max(1,min(nsp-1,ceiling(pdum/delx)))
  dpx=pdum-delx*real(i-1)
! expand y(x) using sprt(x) near x=0
  if(i==1)dpx=sqrt(dpx)
  dp2=dpx*dpx
  spline1=y(1,i)+dpx*y(2,i)+dp2*y(3,i)

end function spline1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! derivative of 1D spline with first point being linear function y=sqrt(x)
real(lk) function dspline1(pdum,nsp,delx,y)
  integer nsp, i
  real(lk) pdum,y(3,nsp),delx,dpx

  i=max(1,min(nsp-1,ceiling(pdum/delx)))
  dpx=pdum-delx*real(i-1)
  if(i==1)dpx=sqrt(dpx)

  if(i==1)then
     dspline1=0.5*y(2,i)/dpx+y(3,i) ! y(x)=y1+y2*sprt(x)+y3*x near x=0
  else
     dspline1=y(2,i)+2.0*dpx*y(3,i) ! y(x)=y1+y2*x+y3*x*x otherwise
  endif

end function dspline1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! periodic spline evaluation 
! Same as spline0 but out of range values are wrapped back in to its periodic
! image location.
real(lk) function spline_periodic(x, nx, delx, y)
  integer, intent(in) :: nx
  real(lk), intent(in) :: x, delx, y(3,nx)
  integer i
  real(lk) dx, x_ev
  x_ev = modulo(x, (nx-1)*delx)
  i=max(1,min(nx-1,ceiling(x/delx)))
  dx=x-delx*real(i-1, lk)
  spline_periodic = y(1,i) + dx*y(2,i) + dx*dx*y(3,i)
end function spline_periodic

real(lk) function dspline_periodic(x, nx, delx, y)
  integer, intent(in) :: nx
  real(lk), intent(in) :: x, delx, y(3,nx)
  integer i
  real(lk) dx, x_ev
  x_ev = modulo(x, (nx-1)*delx)
  i=max(1,min(nx-1,ceiling(x/delx)))
  dx=x-delx*real(i-1, lk)
  dspline_periodic = y(2,i) + 2*dx*y(3,i)
end function dspline_periodic

!#############################################################################


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic 2D spline function
! Purpose: Evaluate f(x,y) or its partial derivatives using the created spline 
!          coefficients for f

! Written by Lei Shi

! Arguments:
!   x,y : the coordinates where the value of f is evaluated
!   deriv: derivative flag. 
!          0 function value, 
!          1 for partial derivative respect to x, 
!          2 for partial derivative respect to y
!          11 for partial derivative respect to x, and then x
!          12 for partial derivative respect to x, and then y
!          21 for partial derivative respect to y, and then x 
!          22 for partial derivative respect to y, and then y
!   nx, ny: number of spline nodes in x and y. 
!   delx, dely: stepsize for each spline cell
!               By default, spline starts from x=0,y=0, and get up to x=(nx-1)*delx,
!               and y=(ny-1)*dely. Out of range (x,y) are not checked
!   f: 2D spline coefficient array constructed by construct_spline2d subroutine
!   bcx, bcy: Boundary conditions used when f is created. Consistency is not checked. 

! WARNING: Boundary Condition consistency is not checked! It is the programer's 
!          responsibility to make sure bcx,bcy passed in is the same as those
!          used for construct_spline2d.

real(lk) function spline2d(x, y, deriv, nx, ny, delx, dely, f, bcx, bcy)
  integer, intent(in):: deriv, nx, ny, bcx, bcy
  real(lk), intent(in):: x, y, delx, dely
  real(lk), dimension(:,:,:), intent(in):: f
  integer i,j
  real(lk) dx, dx2, dy, dy2, &
           xtmp, ytmp
  
  xtmp=max(0.0_lk,x)
  ytmp=max(0.0_lk,y)
! locate the indexes for the spline cell, out of range values are not checked!
  i = max(1, min(nx-1, ceiling(xtmp/delx)))
  j = max(1, min(ny-1, ceiling(ytmp/dely)))
  dx = xtmp - delx*real(i-1, lk)
  dy = ytmp - dely*real(j-1, lk)
! Boundary Condition 1 requires first cell in form f = f1+f2*sqrt(x)+f3*x
  if (i==1 .and. bcx==1) dx=sqrt(dx)
  if (j==1 .and. bcy==1) dy=sqrt(dy)
  dx2 = dx*dx
  dy2 = dy*dy
  
  spline2d=0.0_lk
  if (deriv==0) then
! evaluate f(x,y)
     spline2d = f(1,i,j)+f(2,i,j)*dx + f(3,i,j)*dx2 +&
          ( f(4,i,j)+f(5,i,j)*dx + f(6,i,j)*dx2)*dy +&
          (f(7,i,j) + f(8,i,j)*dx + f(9,i,j)*dx2)*dy2
  else if (deriv == 1) then
! evaluate df/dx
     if (i==1 .and. bcx==1) then
        spline2d = (f(2,i,j) + f(5,i,j)*dy + f(8,i,j)*dy2)/dx +&
             (f(3,i,j) + f(6,i,j)*dy + f(9,i,j)*dy2)
     else
        spline2d = (f(2,i,j) + f(5,i,j)*dy + f(8,i,j)*dy2) +&
             2.0_lk*(f(3,i,j) + f(6,i,j)*dy + f(9,i,j)*dy2)*dx
     endif
  else if (deriv == 2) then
     if (j==1 .and. bcy==1) then
        spline2d = (f(4,i,j) + f(5,i,j)*dx + f(6,i,j)*dx2)/dy +&
             (f(7,i,j) + f(8,i,j)*dy + f(9,i,j))
     else
        spline2d = (f(4,i,j) + f(5,i,j)*dx + f(6,i,j)*dx2) +&
             2.0_lk*(f(7,i,j) + f(8,i,j)*dx + f(9,i,j)*dx2)*dy
     endif
  else if (deriv == 11) then

        spline2d = 2.0_lk*(f(3,i,j) + f(6,i,j)*dy + f(9,i,j)*dy2)

  else if (deriv == 12) then

        spline2d = (f(5,i,j) + 2.0_lk*f(8,i,j)*dy) +2.0_lk*(f(6,i,j) + 2.0_lk*f(9,i,j)*dy)*dx

  else if (deriv == 21) then

        spline2d = (f(5,i,j) + 2.0_lk*f(6,i,j)*dx) +2.0_lk*(f(8,i,j) + 2.0_lk*f(9,i,j)*dx)*dy

  else if (deriv == 22) then
        spline2d =  2.0_lk*(f(7,i,j) + f(8,i,j)*dx + f(9,i,j)*dx2)
  else
     write(gtcout, *) "Spline2d Error: wrong derivative flag deriv=",deriv
     stop
  endif
end function spline2d

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic 3D spline function
! Purpose: Evaluate f(x,y,z) or its partial derivatives using the created spline 
!          coefficients for f

! Written by Lei Shi

! Arguments:
!   x,y,z : the coordinates where the value of f is evaluated
!   deriv: derivative flag. 0 function value, 1 for partial derivative respect
!          to x, 2 for partial derivative respect to y, 3 for partial derivative 
!          respect to z
!   nx, ny, nz: number of spline nodes in x, y and z. 
!   delx, dely, delz: stepsize for each spline cell
!               By default, spline starts from x=0,y=0,z=0 and get up to x=(nx-1)*delx,
!               and y=(ny-1)*dely, z=(nz-1)*delz. Out of range (x,y,z) are not checked
!   f: 3D spline coefficient array constructed by construct_spline3d subroutine
!   bcx, bcy, bcz: Boundary conditions used when f is created. 
!                  Consistency is not checked. 

! WARNING: Boundary Condition consistency is not checked! It is the programer's 
!          responsibility to make sure bcx,bcy,bcz passed in is the same as those
!          used for construct_spline3d.

real(lk) function spline3d(x, y, z, deriv, nx, ny, nz, delx, dely, delz,  f, bcx, bcy, bcz)
  integer, intent(in):: deriv, nx, ny, nz, bcx, bcy, bcz
  real(lk), intent(in):: x, y, z, delx, dely, delz, f(27, nx, ny, nz)
  integer i, j, k
  real(lk) dx, dxinv, dy, dyinv, dz, dzinv, dvec(27)
  
! locate the indexes for the spline cell, out of range values are not checked!
  i = max(1, min(nx-1, ceiling(x/delx)))
  j = max(1, min(ny-1, ceiling(y/dely)))
  k = max(1, min(nz-1, ceiling(z/delz)))
  dx = x - delx*real(i-1, lk)
  dy = y - dely*real(j-1, lk)
  dz = z - delz*real(k-1, lk)
  dxinv = 1.0_lk/dx
  dyinv = 1.0_lk/dy
  dzinv = 1.0_lk/dz

! Boundary Condition 1 requires first cell in form f = f1+f2*sqrt(x)+f3*x
  if (i==1 .and. bcx==1) dx=sqrt(dx)
  if (j==1 .and. bcy==1) dy=sqrt(dy)
  if (k==1 .and. bcz==1) dz=sqrt(dz)
  
! Construct the spline vector
! take care of the boundary condition
  if (deriv==0) then
     dvec(1) = 1
     dvec(2) = dx
     dvec(3) = dx*dx
     dvec(4:6) = dvec(1:3)*dy
     dvec(7:9) = dvec(4:6)*dy
     dvec(10:18) = dvec(1:9)*dz
     dvec(19:27) = dvec(10:18)*dz
  else if (deriv==1) then
! partial derivative in x
     if (bcx==1) then
        dvec(1) = 0
        dvec(2) = 0.5_lk*dxinv
        dvec(3) = 1
        dvec(4:6) = dvec(1:3)*dy
        dvec(7:9) = dvec(4:6)*dy
        dvec(10:18) = dvec(1:9)*dz
        dvec(19:27) = dvec(10:18)*dz      
     else
        dvec(1) = 0
        dvec(2) = dx
        dvec(3) = 2*dx
        dvec(4:6) = dvec(1:3)*dy
        dvec(7:9) = dvec(4:6)*dy
        dvec(10:18) = dvec(1:9)*dz
        dvec(19:27) = dvec(10:18)*dz        
     endif
  else if (deriv == 2) then
! partial derivative in y
     if (bcy==1) then
        dvec(1:3) = 0
        dvec(7) = 1
        dvec(8) = dx
        dvec(9) = dx*dx
        dvec(4:6) = dvec(7:9)*0.5_lk*dyinv
        dvec(10:18) = dvec(1:9)*dz
        dvec(19:27) = dvec(10:18)*dz        
     else
        dvec(1:3) = 0
        dvec(4) = 1
        dvec(5) = dx
        dvec(6) = dx*dx
        dvec(7:9) = dvec(4:6)*2*dy
        dvec(10:18) = dvec(1:9)*dz
        dvec(19:27) = dvec(10:18)*dz  
     endif    
  else if (deriv == 3) then
     if (bcz==1) then
        dvec(1:9)=0
        dvec(19) = 1
        dvec(20) = dx
        dvec(21) = dx*dx
        dvec(22:24) = dvec(19:21)*dy
        dvec(25:27) = dvec(22:24)*dy
        dvec(10:18) = dvec(19:27)*dzinv*0.5_lk
     else
        dvec(1:9)=0
        dvec(10) = 1
        dvec(11) = dx
        dvec(12) = dx*dx
        dvec(13:15) = dvec(10:12)*dy
        dvec(16:18) = dvec(13:15)*dy
        dvec(19:27) = dvec(10:18)*dz*2
     endif
  else
     write(gtcout, *) "Spline2d Error: wrong derivative flag deriv=",deriv
     stop
  endif

  spline3d = sum(f(:,i,j,k)*dvec(:))

end function spline3d

!##############################################################################
! Construct spline subroutines
! Generic subroutines are provided for 1D, 2D, and 3D quadratic splines
!    construct_spline1d
!    construct_spline2d
!    construct_spline3d
! Three kinds of boundary conditions are available for each dimension
!    BC 0: Use 341 formula to obtain the estimated derivative at x=0
!    BC 1: Use 341 formula to obtain the estimated derivative at x=xmax, and use
!          f(x) = a + b*sqrt(x) + c*x formula for the first cell at x=0.
!    BC 2: Periodic condition, df/dx(0)=df/dx(xmax)
!###############################################################################

! Generic construct 1D spline 
subroutine construct_spline1d(nx, delx, f, bcx)
  integer, intent(in):: nx, bcx
  real(lk), intent(in):: delx
  real(lk), intent(inout):: f(3,nx)

  if (bcx == 0) then
     call construct_spline0(0,nx,delx,f)
  else if (bcx == 1) then
     call construct_spline1(nx,delx,f)
  else if (bcx == 2) then
     ! smooth the data before construct the spline
     ! periodic spline is very sensitive to the data smoothness 
     ! smooth_type using "14641" by default. Change it into "smooth_type=0" for "121" formula.
     call smooth_1d(f(1,:), nx, boundary_condition=1, smooth_type=1)
     call construct_spline_periodic(nx, delx, f)
  else
     print *,'Wrong spline boundary condition choice:', bcx
     stop 
  endif
end subroutine construct_spline1d

! Generic construct 2D spline
subroutine construct_spline2d(nx,ny,delx,dely,f, bcx, bcy)
  integer, intent(in):: nx, ny, bcx, bcy
  real(lk), intent(in):: delx, dely
  real(lk), intent(inout):: f(9,nx,ny)
  integer i,j,s
  real(lk) ddum1(3,nx),ddum2(3,ny)

  ddum1=0.0
  ddum2=0.0
  
! periodic condition needs to be enforced if bc is set to 2
  if (bcx==2) then
     f(1,1,:) = 0.5_lk*(f(1,1,:)+f(1,nx,:))
     f(1,nx,:) = f(1,1,:)
  endif
  if (bcy==2) then
     do j = 1, ny-1
        ddum1(1,:)=f(1,:,j)
        call construct_spline1d(nx, delx, ddum1, bcx)
        f(1,:,j)=ddum1(1,:) ! the function value may be altered due to smoothing
        f(2,:,j)=ddum1(2,:)
        f(3,:,j)=ddum1(3,:)
     enddo
     f(1,:,ny) = f(1,:,1)
     f(2,:,ny) = f(2,:,1)
     f(3,:,ny) = f(3,:,1)
  else
     do j = 1, ny
        ddum1(1,:)=f(1,:,j)
        call construct_spline1d(nx, delx, ddum1, bcx)
        f(1,:,j)=ddum1(1,:) ! the function value may be altered due to smoothing
        f(2,:,j)=ddum1(2,:)
        f(3,:,j)=ddum1(3,:)
     enddo
  endif

  do i = 1, nx
     do s = 1, 3
        ddum2(1,:)=f(s,i,:)
        call construct_spline1d(ny, dely, ddum2, bcy)
        f(s,i,:)=ddum2(1,:)
        f(s+3,i,:)=ddum2(2,:)
        f(s+6,i,:)=ddum2(3,:)
     enddo
  enddo
end subroutine construct_spline2d

! generic 3d quadratic spline
subroutine construct_spline3d(nx,ny,nz,delx,dely,delz,f,bcx,bcy,bcz)
  integer, intent(in):: nx, ny, nz, bcx, bcy, bcz
  real(lk), intent(in):: delx, dely, delz
  real(lk), intent(inout):: f(27,nx,ny,nz)

  integer i,j,k,s
  real(lk) temp1d(3,nz), temp2d(9,nx,ny)

! construct 2D spline on each x-y plane
! periodic in z will be enforced if bcz==2

  if (bcz==2) then
     do i=1,nz-1
        temp2d(1,:,:) = f(1,:,:,i)
        call construct_spline2d(nx,ny,delx,dely,temp2d,bcx,bcy)
        f(1:9,:,:,i) = temp2d(:,:,:)
     enddo
     f(1:9,:,:,nz)=f(1:9,:,:,1)
  else
     do i=1,nz
        temp2d(1,:,:) = f(1,:,:,i)
        call construct_spline2d(nx,ny,delx,dely,temp2d,bcx,bcy)
        f(1:9,:,:,i) = temp2d(:,:,:)
     enddo
  endif

! generate spline coefficient on z direction
  do i=1,nx
     do j=1,ny
        do s=1,9
           temp1d(1,:) = f(s,i,j,:)
           call construct_spline1d(nz, delz, temp1d, bcz)
           f(s,i,j,:) = temp1d(1,:)
           f(s+9,i,j,:) = temp1d(2,:)
           f(s+18,i,j,:) = temp1d(3,:)
        enddo
     enddo
  enddo
  
end subroutine construct_spline3d

! construct 1D spline with on x=[0,xmax], grid x_i = (i-1)*delx, delx=xmax/(nsp-1)
! in domain i, y = y(1,i) + y(2,i)*delx + y(3,i)*delx**2
subroutine construct_spline0(iflag,nsp,delx,y)
  integer i,nsp,ipp,iflag
  real(lk) delx,y(3,nsp)

! first point
  if(iflag==0)then
! iflag=0: first point being y=y1+y2*x+y3*x*x
     y(2,1)=(4.0_lk*y(1,2)-y(1,3)-3.0_lk*y(1,1))/(2.0_lk*delx)
     y(3,1)=(y(1,2)-y(1,1)-y(2,1)*delx)/(delx*delx)

  elseif(iflag==1)then
! iflag=1: first point being linear function y=y_1+y2*x
     y(2,1)=(y(1,2)-y(1,1))/delx
     y(3,1)=0.0_lk

  elseif(iflag==2)then
! iflag=2: first point being quadratic function y=y1+y3*x*x
     y(2,1)=0.0_lk
     y(3,1)=(y(1,2)-y(1,1))/(delx*delx)
  endif

  do i=2,nsp-2
     ipp=min(i+2,nsp)
     y(2,i)=-y(2,i-1)+2.0_lk*(y(1,i)-y(1,i-1))/delx

! smooth f1
     y(1,i+1)=0.5_lk*delx*y(2,i)+0.25_lk*y(1,ipp)+0.75_lk*y(1,i)
  enddo

  y(2,nsp-1)=-y(2,nsp-2)+2.0_lk*(y(1,nsp-1)-y(1,nsp-2))/delx
  y(2,nsp)=-y(2,nsp-1)+2.0_lk*(y(1,nsp)-y(1,nsp-1))/delx

  do i=2,nsp-1
     y(3,i)=(y(2,i+1)-y(2,i))/(2.0_lk*delx)
  enddo

! last point is not used;
  y(3,nsp)=0.0_lk

end subroutine construct_spline0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Construct 1D Quadratic Spline using periodic boundary condition
! Written by Lei Shi, Oct 16, 2016

! Check the documentation for a detailed discussion
subroutine construct_spline_periodic(nsp,delx,y)
  integer, intent(in) :: nsp
  real(lk), intent(in) :: delx
  real(lk), intent(inout):: y(3, nsp)
! local variables
  integer :: i
  real(lk) :: dxinv

! Periodic boundary condition requires the first order derivative to be continuous at the end points
! Check the documentation for a detailed discussion about the spline algorithm

dxinv = 1.0_lk/delx

! First, let's check if the input data satisfies the requirements
! Check grid point number
if (mod(nsp, 2) /= 0) then
   write(gtcout,*) 'Periodic Spline Error: Total number of data must be even. nsp=',nsp,' doesn''t work. Use normal spline or request a new input data set.'
   stop
endif

! check periodicity in raw data
if (abs(y(1,1)-y(1,nsp)) > 2*SPACING(y(1,1)) ) then
   write(gtcout,*) 'Periodic Spline Error: function values at end points must be the same. y(1)=',y(1,1),', y(n)=',y(1,nsp),' doesn''t work. Use normal spline or request a new input data set.'
   stop
endif
! reset all the higher order spline coefficients
y(2:3,:) = 0.0_lk
! Now we let y = a + b*x + c*x**2 on each section, first step is to calculate b1

do i=1,nsp/2-1
   y(2,1) = y(2,1) + y(1,2*i) - y(1, 2*i+1)
enddo
y(2,1) = y(2,1)*2.0_lk*dxinv
! we use y(2,nsp) as a buffer to save the periodic b1
y(2, nsp)=y(2,1)
! Then, we can obtain bi one by one
do i=2, nsp-1
   y(2, i) = 2.0_lk*(y(1,i) - y(1,i-1))*dxinv - y(2,i-1)
enddo
! Then the ci's
do i=1, nsp-1
   y(3, i) = (y(2, i+1)-y(2, i))*0.5_lk*dxinv
enddo
end subroutine construct_spline_periodic


! inversion of spline0 y(x) to x(y)
subroutine invert_spline0(iflag,nsp,delx,dely,y,x)
  integer i,nsp,j,iflag
  real(lk) delx,dely,y(3,nsp),x(3,nsp),ydum,y0,y1,y2

! first point given by input
  x(1,1)=0.0_lk
! other points
  do i=2,nsp-1
! y grid
     ydum=dely*real(i-1)
! search x grid for ydum
     j=1
     do while (ydum>y(1,j+1))
        j=j+1
     enddo

! x(1,i)=j grid location + distance from j grid
     y0=y(1,j)
     y1=y(2,j)
     y2=y(3,j)
     if (abs(y2)>0.000001_lk) then
       x(1,i)=delx*real(j-1)+(sqrt(y1*y1+4.0_lk*y2*(ydum-y0))-y1)/(2.0_lk*y2)
     else
       x(1,i)=delx*real(j-1)+(ydum-y0)/y1
     endif

  enddo

! last point
  x(1,nsp)=x(1,1)+delx*real(nsp-1)

! spline fit x function
  if(iflag==0)call construct_spline0(0,nsp,dely,x)

end subroutine invert_spline0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! spline1 for cases with first point being function y=y1+y2*sqrt(x)+y3*x
subroutine construct_spline1(nsp,delx,f)
  integer i,nsp,ipp
  real(lk) delx,f(3,nsp)

! first point
  f(2,1)=(2.0_lk*f(1,2)-f(1,3)-f(1,1))/((2.0_lk-sqrt(2.0_lk))*sqrt(delx))
  f(3,1)=(f(1,2)-f(1,1)-f(2,1)*sqrt(delx))/delx

! second point
  f(2,2)=0.5_lk*f(2,1)/sqrt(delx)+f(3,1)
  f(3,2)=(f(1,3)-f(1,2)-delx*f(2,2))/(delx*delx)

  do i=3,nsp-2
     ipp=min(i+2,nsp)
     f(2,i)=-f(2,i-1)+2.0_lk*(f(1,i)-f(1,i-1))/delx

! smooth f1
     f(1,i+1)=0.5_lk*delx*f(2,i)+0.25_lk*f(1,ipp)+0.75_lk*f(1,i)
  enddo

  f(2,nsp-1)=-f(2,nsp-2)+2.0_lk*(f(1,nsp-1)-f(1,nsp-2))/delx
  f(2,nsp)=-f(2,nsp-1)+2.0_lk*(f(1,nsp)-f(1,nsp-1))/delx

  do i=3,nsp-1  !!change start from 1 to 3
     f(3,i)=(f(2,i+1)-f(2,i))/(2.0_lk*delx)
  enddo

! last point is not used;
  f(3,nsp)=0.0_lk

end subroutine construct_spline1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! inversion of spline1 y(x) to x(y)
subroutine invert_spline1(nsp,delx,dely,y,x)
  integer i,nsp,j
  real(lk) delx,dely,y(3,nsp),x(3,nsp),ydum,y0,y1,y2

! first point is given by inputs
  x(1,1)=0.0_lk
! other points
  do i=2,nsp-1
! y grid
     ydum=dely*real(i-1,lk)
! search x grid for ydum
     j=1
     do while (ydum>y(1,j+1))
        j=j+1
     enddo

! x(1,i)=j grid location + distance from j grid
     y0=y(1,j)
     y1=y(2,j)
     y2=y(3,j)

     if (abs(y2)>0.000001_lk) then
       x(1,i)=delx*real(j-1,lk)+(sqrt(y1*y1+4.0_lk*y2*(ydum-y0))-y1)/(2.0_lk*y2)
     else
       x(1,i)=delx*real(j-1,lk)+(ydum-y0)/y1
     endif

     if (j==1) x(1,i)=x(1,i)*x(1,i)

  enddo
! last point
  x(1,nsp)=delx*real(nsp-1,lk)

! spline fit x function
! call spline with first point being ~ sqrt(r)
  call construct_spline1(nsp,dely,x)

end subroutine invert_spline1

! Deprecated, should be removed after other parts are changed.
subroutine construct_spline2d_dep(nx,ny,delx,dely,y)
  integer i,j,s,nx,ny,ipp
  real(lk) delx,dely,y(9,nx,ny),ddum1(3,nx),ddum2(3,ny)

  ddum1=0.0_lk
  ddum2=0.0_lk

  do j = 1, ny
     ddum1(1,:)=y(1,:,j)     
     call construct_spline1(nx,delx,ddum1)
     y(1,:,j)=ddum1(1,:)
     y(2,:,j)=ddum1(2,:)
     y(3,:,j)=ddum1(3,:)
  enddo
  do i = 2, nx
     do s = 1, 3
        ddum2(1,:)=y(s,i,:)
        call construct_spline0(0,ny,dely,ddum2)
        y(s,i,:)=ddum2(1,:)
        y(s+3,i,:)=ddum2(2,:)
        y(s+6,i,:)=ddum2(3,:)
     enddo
  enddo
end subroutine construct_spline2d_dep
!!!!!!!!!!!!!

! Deprecated, should be removed after other parts are changed.
! Construct 2D spline with the second dimension periodic
! Written by Lei Shi, Oct 16, 2016
! See subroutine construct_spline_periodic for more details
subroutine construct_spline2d_periodic_dep(nx,ny,delx,dely,z)
  integer i,j,s,nx,ny,ipp
  real(lk) delx,dely,z(9,nx,ny),ddum1(3,nx),ddum2(3,ny)

  ddum1=0.0_lk
  ddum2=0.0_lk

  do j = 1, ny-1
     ddum1(1,:)=z(1,:,j)
     call construct_spline1(nx,delx,ddum1)
     z(1,:,j)=ddum1(1,:)
     z(2,:,j)=ddum1(2,:)
     z(3,:,j)=ddum1(3,:)
  enddo
  ! Values at periodic boundary are forced to be equal
  z(1,:,ny)=z(1,:,1)
  z(2,:,ny)=z(2,:,1)
  z(3,:,ny)=z(3,:,1)

  do i = 2, nx
     do s = 1, 3
        ddum2(1,:)=z(s,i,:)
        call construct_spline_periodic(ny,dely,ddum2)
        z(s,i,:)=ddum2(1,:)
        z(s+3,i,:)=ddum2(2,:)
        z(s+6,i,:)=ddum2(3,:)
     enddo
  enddo
end subroutine construct_spline2d_periodic_dep

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine cross(vec1,vec2,vec3)
  real(lk):: vec1(3),vec2(3),vec3(3)
  
  vec3(1)=vec1(2)*vec2(3)-vec1(3)*vec2(2)
  vec3(2)=vec1(3)*vec2(1)-vec1(1)*vec2(3)
  vec3(3)=vec1(1)*vec2(2)-vec1(2)*vec2(1)
  
end subroutine cross
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

real(lk) function dots(vec1,vec2,n)
  integer i,n
  real(lk) ::temp,vec1(n),vec2(n)
  temp=0.0_lk
  do i=1,n
    temp=temp+vec1(i)*vec2(i)
  enddo
  dots=temp
end function dots
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Subroutine sppsiR1d
! Purpose: For FRC field line mesh creation, given Z and psi values, 
!          calculate the corresponding R coordinate based on the psi(R,Z) 
!          spline function.
! Arguments: x, y, z: spline z(x, y) mesh and coefficients
!            yi, zi: the y, z value required for calculate xi
!            xi: output, the calculated x satisfying z(xi, yi)=zi to certain 
!                accuracy
!            zidum: output, the actual z at xi, yi, zidum=z(xi,yi)
!            iflag: 0 is to find first root for Z function ; 1 is to find second root for Z function.

subroutine sppsiR1d(iflag,nx,ny,dx,dy,x,y,z,yi,zi,xi,zidum)
    use precision
    use global_parameters,only:gtcout
    implicit none

    integer,intent(in) :: nx,ny,iflag
    real(lk),intent(in) :: dx,dy,x(nx),y(ny),z(9,nx,ny),yi,zi
    real(lk),intent(out) :: xi,zidum

    integer :: i,j,ii,s
    real(lk) :: dpx,dpy,dx2,dy2,a,b,c,sign_r,zdum,phigh,plow,del_zdum,sperr,c1,c2
!Find the nearest j point (y direction)
    j=max(1,min(ny-1,ceiling(yi/dy)))
    dpy=yi-dy*real(j-1)
    dy2=dpy*dpy

    xi=0.0_lk
    zidum=0.0_lk


!Find the nearest ii point(x direction),only valid for the monotonos Z funtion
    do i=1,nx-1
        c1=z(1,i,j)+z(4,i,j)*dpy+z(7,i,j)*dy2
        c2=z(1,i+1,j)+z(4,i+1,j)*dpy+z(7,i+1,j)*dy2
        sign_r=(zi-c1)*(zi-c2)   ! find the equilibrium cell containing the solution
   
  
        if(sign_r .lt. 0.0_lk)then
            ii=i
  
!Find the x root of Z(x,y)=0 with known Z and y
  
            a=z(3,ii,j)+z(6,ii,j)*dpy+z(9,ii,j)*dy2
            b=z(2,ii,j)+z(5,ii,j)*dpy+z(8,ii,j)*dy2
            c=z(1,ii,j)+z(4,ii,j)*dpy+z(7,ii,j)*dy2
            phigh=dx
            plow=0.0_lk
            dpx=0.5_lk*(phigh+plow)
            dx2=dpx*dpx
            zdum=a*dx2+b*dpx+c
            sperr=0.000001_lk
            del_zdum=(zi-zdum)/abs(zi)  ! valid for both negative and positive zi
            s=1
            do while (abs(del_zdum)>sperr .and. s<1000)
                s=s+1
                if(zi>c2)then      ! case for monotonic decreasing
                    if(del_zdum<0)then
                        plow=dpx
                    else
                        phigh=dpx
                    endif
                elseif(zi<c2)then  ! case for monotonic increasing
                    if(del_zdum>0)then
                        plow=dpx
                    else
                        phigh=dpx
                    endif
                endif
  
                dpx=0.5_lk*(plow+phigh)
                dx2=dpx*dpx
                zdum=a*dx2+b*dpx+c
                del_zdum=(zi-zdum)/abs(zi)  ! valid for both negative and positive zi
            enddo

!Recalculate the Z value by using the x and y position
            xi=dpx+dx*real(ii-1)
            zidum= z(1,ii,j)    + z(2,ii,j)*dpx     + z(3,ii,j)*dx2     &
                   + z(4,ii,j)*dpy  + z(5,ii,j)*dpx*dpy + z(6,ii,j)*dx2*dpy &
                   + z(7,ii,j)*dy2 + z(8,ii,j)*dpx*dy2 + z(9,ii,j)*dx2*dy2
 
            if(iflag==0)then ! exit after find the first root
                exit
            endif
        endif
    enddo




end subroutine sppsiR1d



end module spline_function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module cylindrical_function
use precision, only: lk
use global_parameters, only: gtcout,pi,mype
implicit none

private
public :: geotheta,dgeotheta

contains



real(lk) function geotheta(x,y)
  real(lk), intent(in) :: x, y
  real(lk) tdum1,tdum2,rdum,zdum

  rdum=x
  zdum=y

  tdum1=zdum/sqrt((rdum-1.0_lk)**2+zdum**2)
  tdum2=(rdum-1.0_lk)/sqrt((rdum-1.0_lk)**2+zdum**2)
  if(zdum .ge. 0.0_lk .and. rdum .ge. 1.0_lk)then
     geotheta=asin(tdum1)
  elseif(zdum .ge. 0.0_lk .and. rdum .lt. 1.0_lk)then
     geotheta=pi-asin(tdum1)
  elseif(zdum .lt. 0.0_lk .and. rdum .le. 1.0_lk)then
     geotheta=abs(asin(tdum1))+pi
  elseif(zdum .lt. 0.0_lk .and. rdum .gt. 1.0_lk)then
     geotheta=2.0_lk*pi-acos(tdum2)
  endif

end function geotheta

real(lk) function dgeotheta(iflag,x,y)
  real(lk), intent(in) :: x, y
  real(lk) rdum,zdum,ldum2
  integer, intent(in) :: iflag

  rdum=x
  zdum=y
  ldum2=(rdum-1.0_lk)**2+zdum**2

  if(iflag==1)then
    dgeotheta=-zdum/ldum2
  elseif(iflag==2)then
    dgeotheta=(rdum-1.0_lk)/ldum2
  else
    if(mype==0)write(gtcout,*)"dthetadr/dthetadz option iflag is wrong"
    stop
  endif
end function dgeotheta

end module cylindrical_function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Bessel I_0(x) function in double precision
! http://www.kurims.kyoto-u.ac.jp/~ooura/bessel.html
function besi0d(x)
    implicit real*8 (a - h, o - z)
    dimension a(0 : 64), b(0 : 69), c(0 : 44)
    data (a(i), i = 0, 12) / &
    8.5246820682016865877d-11, 2.5966600546497407288d-9, &
    7.9689994568640180274d-8, 1.9906710409667748239d-6, &
    4.0312469446528002532d-5, 6.4499871606224265421d-4, &
    7.9012345761930579108d-3, 7.1111111109207045212d-2, &
    4.4444444444472490900d-1, 1.7777777777777532045d0, &
    4.0000000000000011182d0, 3.9999999999999999800d0, &
    1.0000000000000000001d0 / 
    data (a(i), i = 13, 25) / &
    1.1520919130377195927d-10, 2.2287613013610985225d-9, &
    8.1903951930694585113d-8, 1.9821560631611544984d-6, &
    4.0335461940910133184d-5, 6.4495330974432203401d-4, &
    7.9013012611467520626d-3, 7.1111038160875566622d-2, &
    4.4444450319062699316d-1, 1.7777777439146450067d0, &
    4.0000000132337935071d0, 3.9999999968569015366d0, &
    1.0000000003426703174d0 / 
    data (a(i), i = 26, 38) / &
    1.5476870780515238488d-10, 1.2685004214732975355d-9, &
    9.2776861851114223267d-8, 1.9063070109379044378d-6, &
    4.0698004389917945832d-5, 6.4370447244298070713d-4, &
    7.9044749458444976958d-3, 7.1105052411749363882d-2, &
    4.4445280640924755082d-1, 1.7777694934432109713d0, &
    4.0000055808824003386d0, 3.9999977081165740932d0, &
    1.0000004333949319118d0 / 
    data (a(i), i = 39, 51) / &
    2.0675200625006793075d-10, -6.1689554705125681442d-10, &
    1.2436765915401571654d-7, 1.5830429403520613423d-6, &
    4.2947227560776583326d-5, 6.3249861665073441312d-4, &
    7.9454472840953930811d-3, 7.0994327785661860575d-2, &
    4.4467219586283000332d-1, 1.7774588182255374745d0, &
    4.0003038986252717972d0, 3.9998233869142057195d0, &
    1.0000472932961288324d0 / 
    data (a(i), i = 52, 64) / &
    2.7475684794982708655d-10, -3.8991472076521332023d-9, &
    1.9730170483976049388d-7, 5.9651531561967674521d-7, &
    5.1992971474748995357d-5, 5.7327338675433770752d-4, &
    8.2293143836530412024d-3, 6.9990934858728039037d-2, &
    4.4726764292723985087d-1, 1.7726685170014087784d0, &
    4.0062907863712704432d0, 3.9952750700487845355d0, &
    1.0016354346654179322d0 / 
    data (b(i), i = 0, 13) / &
    6.7852367144945531383d-8, 4.6266061382821826854d-7, &
    6.9703135812354071774d-6, 7.6637663462953234134d-5, &
    7.9113515222612691636d-4, 7.3401204731103808981d-3, &
    6.0677114958668837046d-2, 4.3994941411651569622d-1, &
    2.7420017097661750609d0, 14.289661921740860534d0, &
    59.820609640320710779d0, 188.78998681199150629d0, &
    399.87313678256011180d0, 427.56411572180478514d0 / 
    data (b(i), i = 14, 27) / &
    1.8042097874891098754d-7, 1.2277164312044637357d-6, &
    1.8484393221474274861d-5, 2.0293995900091309208d-4, &
    2.0918539850246207459d-3, 1.9375315654033949297d-2, &
    1.5985869016767185908d-1, 1.1565260527420641724d0, &
    7.1896341224206072113d0, 37.354773811947484532d0, &
    155.80993164266268457d0, 489.52113711585409180d0, &
    1030.9147225169564806d0, 1093.5883545113746958d0 / 
    data (b(i), i = 28, 41) / &
    4.8017305613187493564d-7, 3.2613178439123800740d-6, &
    4.9073137508166159639d-5, 5.3806506676487583755d-4, &
    5.5387918291051866561d-3, 5.1223717488786549025d-2, &
    4.2190298621367914765d-1, 3.0463625987357355872d0, &
    18.895299447327733204d0, 97.915189029455461554d0, &
    407.13940115493494659d0, 1274.3088990480582632d0, &
    2670.9883037012547506d0, 2815.7166284662544712d0 / 
    data (b(i), i = 42, 55) / &
    1.2789926338424623394d-6, 8.6718263067604918916d-6, &
    1.3041508821299929489d-4, 1.4282247373727478920d-3, &
    1.4684070635768789378d-2, 1.3561403190404185755d-1, &
    1.1152592585977393953d0, 8.0387088559465389038d0, &
    49.761318895895479206d0, 257.26842323135291380d0, &
    1066.8543146269566231d0, 3328.3874581009636362d0, &
    6948.8586598121634874d0, 7288.4893398212481055d0 / 
    data (b(i), i = 56, 69) / &
    3.4093503681970328930d-6, 2.3079025203103376076d-5, &
    3.4691373283901830239d-4, 3.7949949772229085450d-3, &
    3.8974209677945602145d-2, 3.5949483804148783710d-1, &
    2.9522878893539528226d0, 21.246564609514287056d0, &
    131.28727387146173141d0, 677.38107093296675421d0, &
    2802.3724744545046518d0, 8718.5731420798254081d0, &
    18141.348781638832286d0, 18948.925349296308859d0 / 
    data (c(i), i = 0, 8) / &
    2.5568678676452702768d-15, 3.0393953792305924324d-14, &
    6.3343751991094840009d-13, 1.5041298011833009649d-11, &
    4.4569436918556541414d-10, 1.7463930514271679510d-8, &
    1.0059224011079852317d-6, 1.0729838945088577089d-4, &
    5.1503226936425277380d-2 / 
    data (c(i), i = 9, 17) / &
    5.2527963991711562216d-15, 7.2021184814210056410d-15, &
    7.2561421229904797156d-13, 1.4823121466731042510d-11, &
    4.4602670450376245434d-10, 1.7463600061788679671d-8, &
    1.0059226091322347560d-6, 1.0729838937545111487d-4, &
    5.1503226936437300716d-2 / 
    data (c(i), i = 18, 26) / &
    1.3365917359358069908d-14, -1.2932643065888544835d-13, &
    1.7450199447905602915d-12, 1.0419051209056979788d-11, &
    4.5804788198059832600d-10, 1.7442405450073548966d-8, &
    1.0059461453281292278d-6, 1.0729837434500161228d-4, &
    5.1503226940658446941d-2 / 
    data (c(i), i = 27, 35) / &
    5.3771611477352308649d-14, -1.1396193006413731702d-12, &
    1.2858641335221653409d-11, -5.9802086004570057703d-11, &
    7.3666894305929510222d-10, 1.6731837150730356448d-8, &
    1.0070831435812128922d-6, 1.0729733111203704813d-4, &
    5.1503227360726294675d-2 / 
    data (c(i), i = 36, 44) / &
    3.7819492084858931093d-14, -4.8600496888588034879d-13, &
    1.6898350504817224909d-12, 4.5884624327524255865d-11, &
    1.2521615963377513729d-10, 1.8959658437754727957d-8, &
    1.0020716710561353622d-6, 1.0730371198569275590d-4, &
    5.1503223833002307750d-2 / 
    w = abs(x)
    if (w .lt. 8.5d0) then
    t = w * w * 0.0625d0
    k = 13 * (int(t))
    y = (((((((((((a(k) * t + a(k + 1)) * t + &
          a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + &
          a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + &
          a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + &
          a(k + 11)) * t + a(k + 12)
    else if (w .lt. 12.5d0) then
    k = int(w)
    t = w - k
    k = 14 * (k - 8)
    y = ((((((((((((b(k) * t + b(k + 1)) * t + &
          b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + &
          b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + &
          b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + &
          b(k + 11)) * t + b(k + 12)) * t + b(k + 13)
    else
    t = 60 / w
    k = 9 * (int(t))
    y = ((((((((c(k) * t + c(k + 1)) * t + &
          c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + &
          c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + &
          c(k + 8)) * sqrt(t) * exp(w)
    end if
    besi0d = y
end function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
