V3FIT
spline.f
1  SUBROUTINE spline(x,y,n,yp1,ypn,y2)
2  USE stel_kinds, ONLY: dp
3  IMPLICIT NONE
4 !-----------------------------------------------
5 ! D u m m y A r g u m e n t s
6 !-----------------------------------------------
7  INTEGER n
8  REAL(dp) :: yp1,ypn,x(n),y(n),y2(n)
9 !Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., yi = f(xi), with
10 !x1 < x2 < .. . < xN, and given values yp1 and ypn for the first derivative of the interpolating
11 !function at points 1 and n, respectively, this routine returns an array y2(1:n) of
12 !length n which contains the second derivatives of the interpolating function at the tabulated
13 !points xi. If yp1 and/or ypn are equal to 1 × 10**30 or larger, the routine is signaled to set
14 !the corresponding boundary condition for a natural spline, with zero second derivative on
15 !that boundary. If yp1/ypn are -1 X 10**30 or smaller, set the bdy conditions so y2(1) = y2(2) or
16 !y2(n) = y2(n-1) (added by SPH)
17 !-----------------------------------------------
18 ! L o c a l V a r i a b l e s
19 !-----------------------------------------------
20  INTEGER i,k
21  REAL(dp) :: p,qn,sig,un,u(n)
22 !-----------------------------------------------
23  IF (yp1.gt..99e30_dp) THEN !The lower boundary condition is set either to be "natural"
24  y2(1)=0
25  u(1)=0
26  ELSE IF (yp1.lt.-0.99e30_dp) THEN !continuous 2nd derivative at boundary, y2(1) = y2(2)
27  y2(1) = 1; u(1) = 0
28  ELSE !or else to have a specified first derivative.
29  y2(1)=-0.5_dp
30  u(1)=(3.0_dp/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
31  ENDIF
32 
33  DO i=2,n-1 !This is the decomposition loop of the tridiagonal
34  !algorithm. y2 and u are used for temporary
35  !storage of the decomposed factors.
36  sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
37  p=sig*y2(i-1)+2
38  y2(i)=(sig-1._dp)/p
39  u(i)=(6*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
40  1 /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
41  ENDDO
42  IF (ypn.gt..99e30_dp) THEN !The upper boundary condition is set either to be "natural"
43  qn=0
44  un=0
45  ELSE IF (ypn.lt.-0.99e30_dp) THEN !continuous 2nd derivative at bdy
46  qn=-1; un=0
47  ELSE !or else to have a specified first derivative.
48  qn=0.5_dp
49  un=(3._dp/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
50  ENDIF
51  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1)
52  DO k=n-1,1,-1 !This is the backsubstitution loop of the tridiagonal algorithm.
53  y2(k)=y2(k)*y2(k+1)+u(k)
54  ENDDO
55 
56  END SUBROUTINE spline