!******************** start file spl1dr_nb.for ; group spl1dr *****************
! jed 8/6/01
! this file was made as a copy of spl1dr.for in the trulib dir
! the copy was made on 8/6/01 and reflects the state of spl1dr at that time
! the function was renamed to spl1dr_nb to remove nubeam dependence on trulib
! all nubeam calls to spl1dr will now be replaced by calls to spl1dr_nb
! this was done to facilitate the moving of msplin.for from trulib to nubeam

subroutine spl1dr_nb(nmax,npts,dx,ca, &
     nbc1,alpha, &
     nbcn,beta, &
     lunerr,ierr)
  use iso_c_binding, only: fp => c_double
  !
  !
  !       bob mccann,  4-jun-84
  !
  !
  !       --------
  !       comments
  !       --------
  !
  !	**************************************************
  !	***this is the same as spl1d.for except for the***
  !	***order of the subscripts of ca	       ***
  !	**************************************************
  !
  !	this is from conte&deboor,p.238.  this solves for
  !   cubic spline coefficients assuming f''(j) = f''(j+1)
  !   and continuous first derivatives of the
  !   intepolating function:
  !
  !	p(j)(x) = (((ca(4,j)*h + ca(3,j))*h + ca(2,j))*h + ca(1,j)
  !
  !		h = (x-j*dx) .le. dx
  !
  !	p'(j)( (j+1)*dx ) = p'(j+1)( (j+1)*dx )
  !
  !	p''(j)( (j+1)*dx ) = p''(j)( (j+1)*dx )
  !
  !	p', p'', or p''' are given at 1 and n
  !		one can also give p'(1) and p'''(n), etc.
  !
  !   where the ca(1,j) are the specified functional values f(j*dx).
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  !***
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !============
  implicit none
  integer :: nm1,nm2,j
  !============
  real(fp) :: rdx,a,b,rhs1,c,d,rhsn,g
  !============
  integer ::	nmax		!the second dim. of c
  integer ::	npts		!the number of functional values
  !
  real(fp) ::	dx		!the constant step size
  real(fp) ::	ca(4,nmax)	!ca(1,j) are specified
  !				!ca(2-4,j) are returned
  integer ::	nbc1, nbcn	!left and right boundary
  !					!condition specifiers
  !					!=(first,second,or third)
  real(fp) ::	alpha, beta	!specified left or right side
  !				!derivative
  !
  integer ::	lunerr		!unit number for error messages
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !	ca(2-4,j) are returned
  !
  integer ::	ierr		!=0 => ok, error otherwise
  !
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !
  !	ca(3,j) is used as the diagonal element in section 1.0
  !
  !	ca(4,j) is used as the right hand side in 1.0
  !
  !
  !       -------------------
  !       statement functions
  !       -------------------
  !
  !
  !***
  !
  !
  !=========================================
  !
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  !
  ierr = 0
  nm1 = npts - 1
  nm2 = npts - 2
  rdx = 1.0_fp/dx
  !
  if(npts.lt.3)	then
     write(lunerr,80301) npts
80301 format(' spl1dr_nb error: the number of points ',&
          'specified =',i2,&
          ' is less than 3.')
     ierr = 3
     return
  endif
  !
  !
  !	---------------------------
  !	0.2	boundary conditions
  !	---------------------------
  !
  !
  if(nbc1.eq.1)	then	!set the left first deriv.
     !
     a = 1.0_fp		!(1,1) element of matrix
     b = 0.0_fp		!(1,2) element of matrix
     rhs1 = alpha
     !
  elseif(nbc1.eq.2)	then	!set left second deriv.
     !
     a = 2.0_fp
     b = 1.0_fp
     rhs1       = 3.0_fp*rdx*(ca(1,2)-ca(1,1)) &
          - dx*alpha/2.0_fp
     !
  elseif(nbc1.eq.3)  then	!set left third deriv.
     !
     a = 1.0_fp
     b = 1.0_fp
     rhs1       = 2.0_fp*rdx*(ca(1,2)-ca(1,1)) &
          + dx*dx*alpha/6.0_fp
     !
  else				!error
     !
     write(lunerr,80101) nbc1
80101 format(' spl1dr_nb error: the left boundary condition ',&
           'specifier=',i5,/,&
           ' this is not one of (1, 2, 3).')
     ierr = 1
     return
     !
  end if
  !
  !   now set the right boundary condition
  !
  if(nbcn.eq.1)	then	!set right first deriv.
     !
     c = 0.0_fp	!(n,n-1) element of the matrix
     d = 1.0_fp	!(n,n) element of the matrix
     rhsn = beta
     !
  elseif(nbcn.eq.2)	then	!set right second deriv.
     !
     c = 1.0_fp
     d = 2.0_fp
     rhsn   = 3.0_fp*rdx*(ca(1,npts)-ca(1,npts-1)) &
          + dx*beta/2.0_fp
     !
  elseif(nbcn.eq.3)	then	!set right third deriv.
     !
     c = 1.0_fp
     d = 1.0_fp
     rhsn   = 2.0_fp*rdx*(ca(1,npts)-ca(1,npts-1)) &
          + dx*dx*beta/6.0_fp
     !
  else				!error
     !
     write(lunerr,80201) nbcn
80201 format(' spl1dr_nb error: the right boundary condition ',&
           'specifier=',i5,/,&
           ' this is not one of (1, 2, 3).')
     ierr = 2
     return
  end if
  !
  !   now define the diagonal elements of the tridiagonal system
  !
   ca(3,2:npts-1) = 4.0_fp
  !
  ca(3,1) = a
  ca(3,npts) = d
  !
  !   all the super-diagonal and sub-diagonal elements are 1.0
  !   except the first super-diagonal element (1,2) and the
  !   last subdiagonal element (n,n-1).
  !
  !
  !   now define the right hand side
  !
  ca(4,2:npts-1) = 3.0_fp*rdx*(ca(1,3:npts) - ca(1,1:npts-2))
  !
  ca(4,1) = rhs1
  ca(4,npts) = rhsn
  !
  !
  !       ----------------------------------------------
  !       1.0     tridiagonal gauss elimination solution
  !       ----------------------------------------------
  !
  !	here we solve the system for ca(2,j)
  !
  !	ca(2,j-1) + ca(3,j)*ca(2,j) + ca(2,j+1) = ca(4,j)
  !
  !   now do the gaussian elimination for the left boundary condition
  !
  ca(3,2) = 4.0_fp- b/ca(3,1)
  ca(4,2) = ca(4,2) - ca(4,1)/ca(3,1)
  !
  !
  !   now do the internal gaussian elimination
  !
  if(nm1.ge.3) then
     do j = 3, nm1
        g = -1.0_fp/ca(3,j-1)
        ca(3,j) = ca(3,j) + g
        ca(4,j) = ca(4,j) + g*ca(4,j-1)
     end do
  endif
  !
  !   now do the right side boundary condition
  !
  ca(3,npts) = ca(3,npts) - c/ca(3,nm1)
  ca(4,npts) = ca(4,npts) - (c/ca(3,nm1))*ca(4,nm1)
  !
  !   now do the back substitution:
  !
  !   first to the right hand endpoint
  !
  ca(2,npts) = ca(4,npts)/ca(3,npts)
  !
  !   then do the interior
  !
  do j = nm1, 2, -1
     ca(2,j) = (ca(4,j) - ca(2,j+1))/ca(3,j)
  end do
  !
  !   finally do the left hand point
  !
  ca(2,1) = (ca(4,1) - b*ca(2,2))/ca(3,1)
  !
  !
  !	-----------------------------------
  !	3.0	calculate ca(3,j) and ca(4,j)
  !	-----------------------------------
  !
  !
  do j = 1, nm1
     !
     ca(3,j) = rdx*( 3.0_fp*rdx*(ca(1,j+1)-ca(1,j)) &
          - 2.0_fp*ca(2,j) - ca(2,j+1) )
     !
     ca(4,j) = rdx*rdx*( ca(2,j) + ca(2,j+1) &
          - 2.0_fp*rdx*(ca(1,j+1)-ca(1,j))  )
     !
  end do
  !
  !   now extrapolate to get the third and fourth coefficients
  !   at the last node (assume continuity of all derivatives).
  !
  ca(4,npts) = ca(4,nm1)		!dx doesn't enter here
  ca(3,npts) = ca(3,nm1) + 3.0_fp*dx*ca(4,nm1) !extrapolated
  
  return
end subroutine spl1dr_nb
!******************** end file spl1dr_nb.for ; group spl1dr ******************
