V3FIT
r8cspeval.f
1 c cspeval -- eval cubic spline function and/or derivatives
2 c
3  subroutine r8cspeval(xget,iselect,fval,x,nx,ilinx,f,ier)
4 c
5 !============
6 ! idecl: explicitize implicit INTEGER declarations:
7  IMPLICIT NONE
8  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
9  INTEGER ier,nx
10 !============
11  real*8 xget ! interpolation target
12  real*8 fval(3) ! output values
13  real*8 x(nx),f(4,nx) ! spline data
14 c
15  integer iselect(3) ! output selector
16  integer ilinx ! =1 if x(...) is evenly spaced
17 c
18 c modification -- dmc 11 Jan 1999 -- remove SAVE stmts; break routine
19 C into these parts:
20 C
21 C cspevx -- find grid cell of target pt.
22 C cspevfn -- evaluate function using output of bcpsevxy
23 C
24 C in cases where multiple functions are defined on the same grid,
25 C time can be saved by using cspevx once and then cspevfn
26 C multiple times.
27 c
28 c input:
29 c (xget) location where interpolated value is desired
30 c x(1).le.xget.le.x(nx) expected
31 c
32 c iselect select desired output
33 c
34 c iselect(1)=1 -- want function value (f) itself
35 c iselect(2)=1 -- want df/dx
36 c iselect(3)=1 -- want d2f/dx2
37 c
38 c example: iselect(1)=iselect(2)=iselect(3)=1
39 c f, df/dx, and d2f/dx2 are all evaluated
40 c and returned in order in fval(1), fval(2),
41 c and fval(3)
42 c iselect(1)=0, iselect(2)=1, iselect(3)=0
43 c only the 1st derivative is evaluated
44 c and returned in fval(1).
45 c
46 c set iselect(1)=3 to get d3f/dx3, 1 value only.
47 c
48 c see fval (output) description.
49 c
50 c x(1...nx) independent coordinate x, strict ascending
51 c
52 c ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
53 c
54 c **CAUTION** actual even spacing of x, is NOT CHECKED HERE!
55 c
56 c f the function values (at grid points) and spline coefs
57 c
58 c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
59 c
60 c spline value =
61 c f(1,i) + dx*f(2,i) + dx**2*f(3,i) + dx**3*f(4,i)
62 c
63 c output:
64 c up to 3 elements of fval, ordered as follows:
65 c fval(1)=function value or lowest order derivative requested
66 c fval(2)=next order derivative
67 c etc
68 c the ordering is a subset of the sequence given under the "iselect"
69 c description.
70 c
71 c ier = 0 -- successful completion; = 1 -- an error occurred.
72 c
73 c-------------------------------------------------------------------
74 c local
75 c
76  integer :: ia(1) = (/ 0 /)
77  real*8 :: dxa(1)
78 c
79 c--------------------------
80 c
81  call r8cspevx(xget,x,nx,ilinx,ia(1),dxa(1),ier)
82  if(ier.ne.0) return
83 c
84  call r8cspevfn(iselect,1,1,fval,ia,dxa,f,nx)
85 c
86  return
87  end
88 c
89 c-------------------------------------------------------------------------
90 c cspevx -- look up x zone
91 c
92 c this is the "first part" of cspeval, see comments, above.
93 c
94  subroutine r8cspevx(xget,x,nx,ilinx,i,dx,ier)
95 c
96 !============
97 ! idecl: explicitize implicit INTEGER declarations:
98  IMPLICIT NONE
99  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
100  INTEGER nxm,ii
101 !============
102 ! idecl: explicitize implicit REAL declarations:
103  real*8 zxget,zxtol
104 !============
105  integer nx ! x array dimension
106 c
107  real*8 xget ! target point
108  real*8 x(nx) ! independent coord. array
109 c
110  integer ilinx ! =1: assume x evenly spaced
111 c
112 c output of cspevx
113 c
114  integer i ! index to cell containing target pt
115  real*8 dx ! displacement of target pt w/in cell
116  ! dx = x-x(i)
117 c
118 c the input argument range is checked...
119 c
120  integer ier ! return ier.ne.0 on error
121 c
122 c------------------------------------
123 c
124  ier=0
125 c
126 c range check
127 c
128  zxget=xget
129  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
130  zxtol=4.0e-7_r8*max(abs(x(1)),abs(x(nx)))
131  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
132  ier=1
133  write(6,1001) xget,x(1),x(nx)
134  1001 format(' ?cspeval: xget=',1pe11.4,' out of range ',
135  > 1pe11.4,' to ',1pe11.4)
136  else
137  if((xget.lt.x(1)-0.5_r8*zxtol).or.
138  > (xget.gt.x(nx)+0.5_r8*zxtol))
139  > write(6,1011) xget,x(1),x(nx)
140  1011 format(' %cspeval: xget=',1pe15.8,' beyond range ',
141  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
142  if(xget.lt.x(1)) then
143  zxget=x(1)
144  else
145  zxget=x(nx)
146  endif
147  endif
148  endif
149  if(ier.ne.0) return
150 c
151 c now find interval in which target point lies..
152 c
153  nxm=nx-1
154 c
155  if(ilinx.eq.1) then
156  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
157  i=min(nxm, ii)
158  if(zxget.lt.x(i)) then
159  i=i-1
160  else if(zxget.gt.x(i+1)) then
161  i=i+1
162  endif
163  else
164  if((1.le.i).and.(i.lt.nxm)) then
165  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
166  continue ! already have the zone
167  else
168  call r8zonfind(x,nx,zxget,i)
169  endif
170  else
171  i=nx/2
172  call r8zonfind(x,nx,zxget,i)
173  endif
174  endif
175 c
176  dx=zxget-x(i)
177 c
178  return
179  end
180 c------------------------------------------------------------------------
181 c cspevfn -- OK now evaluate the cubic spline
182 c
183  subroutine r8cspevfn(ict,ivec,ivd,fval,iv,dxv,f,nx)
184 c
185 c input:
186 !============
187 ! idecl: explicitize implicit INTEGER declarations:
188  IMPLICIT NONE
189  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
190  INTEGER nx,iaval,i
191 !============
192 ! idecl: explicitize implicit REAL declarations:
193  real*8 dx
194 !============
195  integer ict(3) ! selector:
196 c ict(1)=1 for f (don't evaluate f if ict(1)=0)
197 c ict(2)=1 for df/dx ""
198 c ict(3)=1 for d2f/dx2
199 c
200 c set ict(1)=3 to get d3f/dx3 (only)
201 c
202  integer ivec,ivd ! vector dimensioning
203 c
204 c ivec-- number of vector pts (spline values to look up)
205 c ivd -- 1st dimension of fval, .ge.ivec
206 c
207 c output:
208  real*8 fval(ivd,*) ! output array
209 c
210 c v = index to element in vector;
211 c fval(v,1) = first item requested by ict(...),
212 c fval(v,2) = 2nd item requested, ...etc...
213 c
214 c input:
215  integer iv(ivec) ! grid cell indices -- vectors
216  real*8 dxv(ivec) ! displacements w/in cell -- vectors
217 c
218  real*8 f(4,nx) ! cubic fcn spline coeffs array
219 c
220 c usage example:
221 c
222 c 1. for each element xx(v) in a vector of x values:
223 c find the x zone index and displacement with respect to the
224 c lower end point of the zone; store thes in vectors iv and dxv.
225 c
226 c 2. set ict(1)=0, ict(2)=1, ict(3)=0 -- get only the 1st derivative
227 c
228 c 3. ivec is the length of the vector; ivd is the 1st dimension of the
229 c array fval to receive the output
230 c
231 c real fval(ivd,3)
232 c real xv(ivd)
233 c real dxv(ivd)
234 c integer iv(ivd)
235 c integer ict(3)
236 c
237 c real fspline(4,nx) ! spline coeffs
238 c data ict/0,1,0/ ! this call: want 1st derivative only
239 c ! this will be output to fval(*,1).
240 c ...
241 c do iv=1,ivec
242 c ... ! comput indices & displacements
243 c enddo
244 c call cspevfn(ict,ivec,ivd,fval,iv,dxv,fspline,nx)
245 c
246 c--------------------------------------------------------------------
247 c local:
248 c
249  integer v ! vector element index
250 c
251 c OK can now do evaluations
252 c
253  iaval=0 ! fval addressing
254 c
255  if(ict(1).eq.3) then
256 c
257 c fxxx = d3f/dx3
258 c
259  iaval=iaval+1
260  do v=1,ivec
261  i=iv(v)
262  fval(v,iaval)=6.0_r8*f(4,i)
263  enddo
264  else
265 c
266 c normal call...
267 c
268  if(ict(1).gt.0) then
269 c evaluate f
270  iaval=iaval+1
271  do v=1,ivec
272  i=iv(v)
273  dx=dxv(v)
274  fval(v,iaval)=f(1,i)+dx*(f(2,i)+dx*(f(3,i)+dx*f(4,i)))
275  enddo
276  endif
277 c
278  if(ict(2).gt.0) then
279 c evaluate df/dx
280  iaval=iaval+1
281  do v=1,ivec
282  i=iv(v)
283  dx=dxv(v)
284  fval(v,iaval)=f(2,i)+dx*(2.0_r8*f(3,i)+dx*3.0_r8*f(4,i))
285  enddo
286  endif
287 c
288  if(ict(3).gt.0) then
289 c evaluate d2f/dx2
290  iaval=iaval+1
291  do v=1,ivec
292  i=iv(v)
293  dx=dxv(v)
294  fval(v,iaval)=2.0_r8*f(3,i)+dx*6.0_r8*f(4,i)
295  enddo
296  endif
297  endif
298 c
299  return
300  end
301 c----------------------