V3FIT
r8cspline.f
1 c cspline -- dmc 15 Feb 1999
2 c
3 c a standard interface to the 1d spline setup routine
4 c modified dmc 3 Mar 2000 -- to use Wayne Houlberg's v_spline code.
5 c new BC options added.
6 c
7  subroutine r8cspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
8  > wk,iwk,ilinx,ier)
9 c
10 !============
11 ! idecl: explicitize implicit INTEGER declarations:
12  IMPLICIT NONE
13  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
14  INTEGER iwk,nx,ierx,inum,i
15 !============
16 ! idecl: explicitize implicit REAL declarations:
17  real*8 half,sixth
18 !============
19  real*8 x(nx) ! x axis (in)
20  real*8 fspl(4,nx) ! spline data (in/out)
21  integer ibcxmin ! x(1) BC flag (in, see comments)
22  real*8 bcxmin ! x(1) BC data (in, see comments)
23  integer ibcxmax ! x(nx) BC flag (in, see comments)
24  real*8 bcxmax ! x(nx) BC data (in, see comments)
25  real*8 wk(iwk) ! workspace of size at least nx
26  integer ilinx ! even spacing flag (out)
27  integer ier ! output, =0 means OK
28 c
29 c ** note wk(...) array is not used unless ibcxmin=-1 (periodic spline
30 c evaluation)
31 c
32 c this routine computes spline coefficients for a 1d spline --
33 c evaluation of the spline can be done by cspeval.for subroutines
34 c or directly by inline code.
35 c
36 c the input x axis x(1...nx) must be strictly ascending, i.e.
37 c x(i+1).gt.x(i) is required for i=1 to nx-1. This is checked and
38 c ier=1 is set and the routine exits if the test is not satisfied.
39 c
40 c on output, ilinx=1 is set if, to a reasonably close tolerance,
41 c all grid spacings x(i+1)-x(i) are equal. This allows a speedier
42 c grid lookup algorithm on evaluation of the spline. If on output
43 c ilinx=2, this means the spline x axis is not evenly spaced.
44 c
45 c the input data for the spline are given in f[j] = fspl(1,j). The
46 c output data are the spline coefficients fspl(2,j),fspl(3,j), and
47 c fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
48 c boundary conditions and with the properties
49 c
50 c s(x(j)) = fspl(1,j)
51 c s'(x) is continuous even at the grid points x(j)
52 c s''(x) is continuous even at the grid points x(j)
53 c
54 c the formula for evaluation of s(x) is:
55 c
56 c let dx = x-x(i), where x(i).le.x.le.x(i+1). Then,
57 c s(x)=fspl(1,i) + dx*(fspl(2,i) +dx*(fspl(3,i) + dx*fspl(4,i)))
58 c
59 c ==>boundary conditions. Complete specification of a 1d spline
60 c requires specification of boundary conditions at x(1) and x(nx).
61 c
62 c this routine provides 4 options:
63 c
64 c -1 ***** PERIODIC BC
65 c ibcxmin=-1 -- periodic boundary condition. This means the
66 c boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
67 c are imposed. Note that s(x(1))=s(x(nx)) (i.e. fspl(1,1)=fspl(1,nx))
68 c is not required -- that is determined by the fspl array input data.
69 c The periodic boundary condition is to be preferred for periodic
70 c data. When splining periodic data f(x) with period P, the relation
71 c x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
72 c (ibcxmax, bcxmin, bcxmax are ignored).
73 c
74 c if a periodic boundary condition is set, this covers both boundaries.
75 c for the other types of boundary conditions, the type of condition
76 c chosen for the x(1) boundary need not be the same as the type chosen
77 c for the x(nx) boundary.
78 c
79 c 0 ***** NOT A KNOT BC
80 c ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
81 c condition -- see cubsplb.for. This is a common way for inferring
82 c a "good" spline boundary condition automatically from data in the
83 c vicinity of the boundary. (bcxmin | bcxmax are ignored).
84 c
85 c 1 ***** BC: SPECIFIED SLOPE
86 c ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
87 c s'(x(nx)) match the passed value (bcxmin | bcxmax).
88 c
89 c 2 ***** BC: SPECIFIED 2nd DERIVATIVE
90 c ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
91 c s''(x(nx)) match the passed value (bcxmin | bcxmax).
92 c
93 c 3 ***** BC: SPECIFIED SLOPE = 0.0
94 c ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
95 c s'(x(nx)) equal to ZERO.
96 c
97 c 4 ***** BC: SPECIFIED 2nd DERIVATIVE = 0.0
98 c ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
99 c s''(x(nx)) equal to ZERO.
100 c
101 c 5 ***** BC: 1st DIVIDED DIFFERENCE
102 c ibcxmin=5 | ibcxmax=5 -- boundary condition is to have s'(x(1)) |
103 c s'(x(nx)) equal to the slope from the 1st|last 2 points
104 c
105 c 6 ***** BC: 2nd DIVIDED DIFFERENCE
106 c ibcxmin=6 | ibcxmax=6 -- boundary condition is to have s''(x(1)) |
107 c s''(x(nx)) equal to the 2nd derivative from the 1st|last 3 points
108 c
109 c 7 ***** BC: 3rd DIVIDED DIFFERENCE
110 c ibcxmin=7 | ibcxmax=7 -- boundary condition is to have s'''(x(1)) |
111 c s'''(x(nx)) equal to the 3rd derivative from the 1st|last 4 points
112 c
113 c---------------------------------------------------------------------
114  data half/0.5_r8/
115  data sixth/0.166666666666666667_r8/
116 c
117 c error checks
118 c
119  ier = 0
120  if(nx.lt.2) then
121  write(6,'('' ?cspline: at least 2 x points required.'')')
122  ier=1
123  endif
124  call ibc_ck(ibcxmin,'cspline','xmin',-1,7,ier)
125  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'cspline','xmax',0,7,ier)
126 c
127 c x axis check
128 c
129  call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
130  if(ierx.ne.0) ier=2
131 c
132  if(ier.eq.2) then
133  write(6,'('' ?cspline: x axis not strict ascending'')')
134  endif
135 c
136  if(ibcxmin.eq.-1) then
137  inum=nx
138  if(iwk.lt.inum) then
139  write(6,1009) inum,iwk,nx
140  1009 format(
141  > ' ?cspline: workspace too small. need: ',i6,' got: ',i6/
142  > ' (need = nx, nx=',i6)
143  ier=3
144  endif
145  endif
146 c
147  if(ier.ne.0) return
148 c
149 c OK -- evaluate spline
150 c
151  if(ibcxmin.eq.1) then
152  fspl(2,1)=bcxmin
153  else if(ibcxmin.eq.2) then
154  fspl(3,1)=bcxmin
155  endif
156 c
157  if(ibcxmax.eq.1) then
158  fspl(2,nx)=bcxmax
159  else if(ibcxmax.eq.2) then
160  fspl(3,nx)=bcxmax
161  endif
162 c
163  call r8v_spline(ibcxmin,ibcxmax,nx,x,fspl,wk)
164 c
165  do i=1,nx
166  fspl(3,i)=half*fspl(3,i)
167  fspl(4,i)=sixth*fspl(4,i)
168  enddo
169 c
170  return
171  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17