Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module bicub_sg
21 :
22 : !implicit none
23 :
24 : contains
25 : ! from PSPLINE by Doug McCune (version as of February, 2004)
26 :
27 :
28 : ! PSPLINE Home Page:
29 : ! http://w3.pppl.gov/NTCC/PSPLINE
30 :
31 : ! Doug McCune, Princeton University
32 : ! dmccune@pppl.gov
33 :
34 :
35 : ! bcspeval -- eval bicubic spline function and/or derivatives
36 :
37 0 : subroutine bcspeval(xget,yget,iselect,fval,x,nx,y,ny,ilinx,iliny,f,inf3,ier)
38 :
39 : integer iselect(6)
40 : integer ilinx,iliny,nx,ny,inf3,ier
41 :
42 : real xget,yget
43 : real fval(6)
44 : !real x(nx),y(ny),f(4,4,inf3,ny)
45 : real x(:),y(:),f(:,:,:,:)
46 :
47 : ! modification -- dmc 11 Jan 1999 -- remove SAVE stmts; break routine
48 : ! into these parts:
49 :
50 : ! bcspevxy -- find grid cell of target pt.
51 : ! bcspevfn -- evaluate function using output of bcpsevxy
52 :
53 : ! in cases where multiple functions are defined on the same grid,
54 : ! time can be saved by using bcspevxy once and then bcspevfn
55 : ! multiple times.
56 :
57 : ! input:
58 : ! (xget,yget) location where interpolated value is desired
59 : ! x(1) <= xget <= x(nx) expected
60 : ! y(1) <= yget <= y(ny) expected
61 :
62 : ! iselect select desired output
63 :
64 : ! iselect(1)=1 -- want function value (f) itself
65 : ! iselect(2)=1 -- want df/dx
66 : ! iselect(3)=1 -- want df/dy
67 : ! iselect(4)=1 -- want d2f/dx2
68 : ! iselect(5)=1 -- want d2f/dy2
69 : ! iselect(6)=1 -- want d2f/dxdy
70 :
71 : ! example: iselect(1)=iselect(2)=iselect(3)=1
72 : ! f, df/dx, and df/dy all evaluated
73 : ! iselect(4)=iselect(5)=iselect(6)=0
74 : ! 2nd derivatives not evaluated.
75 :
76 : ! the number of non zero values iselect(1:6)
77 : ! determines the number of outputs...
78 : ! see fval (output) description.
79 :
80 : ! new dmc December 2005 -- access to higher derivatives (even if not
81 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
82 : ! if iselect(1)=3 -- want 3rd derivatives
83 : ! iselect(2)=1 for d3f/dx3
84 : ! iselect(3)=1 for d3f/dx2dy
85 : ! iselect(4)=1 for d3f/dxdy2
86 : ! iselect(5)=1 for d3f/dy3
87 : ! number of non-zero values iselect(2:5) gives no. of outputs
88 : ! if iselect(1)=4 -- want 4th derivatives
89 : ! iselect(2)=1 for d4f/dx3dy
90 : ! iselect(3)=1 for d4f/dx2dy2
91 : ! iselect(4)=1 for d4f/dxdy3
92 : ! number of non-zero values iselect(2:4) gives no. of outputs
93 : ! if iselect(1)=5 -- want 5th derivatives
94 : ! iselect(2)=1 for d5f/dx3dy2
95 : ! iselect(3)=1 for d5f/dx2dy3
96 : ! number of non-zero values iselect(2:3) gives no. of outputs
97 : ! if iselect(1)=6 -- want 6th derivatives
98 : ! d6f/dx3dy3 -- one value is returned.
99 :
100 : ! x(1...nx) independent coordinate x, strict ascending
101 : ! y(1...ny) independent coordinate y, strict ascending
102 :
103 : ! ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
104 : ! iliny -- =1: flag that y is linearly spaced (avoid search for speed)
105 :
106 : ! **CAUTION** actual even spacing of x, y is NOT CHECKED HERE!
107 :
108 : !
109 : ! f the function values (at grid points) and spline coefs
110 :
111 : ! evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
112 : ! and y btw y(j) and y(j+1), dy=y-y(j),
113 :
114 : ! spline value =
115 : ! f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
116 : ! +dy*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
117 : ! +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
118 : ! +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
119 :
120 : ! where d2=dy**2 and d3=dy**3.
121 :
122 : ! output:
123 : ! up to 6 elements of fval, ordered as follows:
124 : ! fval(1)=function value or lowest order derivative requested
125 : ! fval(2)=next order derivative
126 : ! etc
127 : ! the ordering is a subset of the sequence given under the "iselect"
128 : ! description.
129 :
130 : ! ier = 0 -- successful completion; = 1 -- an error occurred.
131 :
132 : !-------------------------------------------------------------------
133 : ! local
134 :
135 0 : real dx,dy
136 : integer ia(1), ja(1)
137 :
138 : !--------------------------
139 :
140 0 : ia(1) = 0
141 0 : ja(1) = 0
142 : call bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
143 0 : & 0,0,dx,dy,ier)
144 0 : if(ier /= 0) return
145 :
146 : call bcspevfn(iselect,1,1,fval,ia,ja,
147 0 : & (/dx/),(/dy/),f,inf3,ny)
148 :
149 0 : return
150 : end subroutine bcspeval
151 :
152 : !-------------------------------------------------------------------------
153 :
154 :
155 : ! bcspevxy -- look up x-y zone
156 :
157 : ! this is the "first part" of bcspeval, see comments, above.
158 :
159 0 : subroutine bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
160 : & i,j,dx,dy,ier)
161 :
162 : integer nx,ny ! array dimensions
163 :
164 : real xget,yget ! target point
165 : !real x(nx),y(ny) ! indep. coords.
166 : real x(:),y(:) ! indep. coords.
167 :
168 : integer ilinx ! =1: assume x evenly spaced
169 : integer iliny ! =1: assume y evenly spaced
170 :
171 : ! output of bcspevxy
172 :
173 : integer i,j ! index to cell containing target pt
174 : real dx,dy ! displacement of target pt w/in cell
175 : ! dx=x-x(i) dy=y-y(j)
176 :
177 : integer ier ! return ier /= 0 on error
178 :
179 : !------------------------------------
180 0 : real zxget, zyget, zxtol, zytol
181 : integer nxm, nym, ii, jj
182 :
183 0 : ier=0
184 :
185 : ! range check
186 :
187 0 : zxget=xget
188 0 : zyget=yget
189 :
190 0 : if((xget < x(1)).or.(xget > x(nx))) then
191 0 : zxtol=4.0e-7*max(abs(x(1)),abs(x(nx)))
192 0 : if((xget < x(1)-zxtol).or.(xget > x(nx)+zxtol)) then
193 0 : ier=1
194 : ! write(6,1001) xget,x(1),x(nx)
195 : ! 1001 format(' ?bcspeval: xget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
196 : else
197 : ! if((xget < x(1)-0.5*zxtol).or. (xget > x(nx)+0.5*zxtol)) write(6,1011) xget,x(1),x(nx)
198 : ! 1011 format(' %bcspeval: xget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
199 0 : if(xget < x(1)) then
200 0 : zxget=x(1)
201 : else
202 0 : zxget=x(nx)
203 : end if
204 : end if
205 : end if
206 0 : if((yget < y(1)).or.(yget > y(ny))) then
207 0 : zytol=4.0e-7*max(abs(y(1)),abs(y(ny)))
208 0 : if((yget < y(1)-zytol).or.(yget > y(ny)+zytol)) then
209 0 : ier=1
210 : ! write(6,1002) yget,y(1),y(ny)
211 : ! 1002 format(' ?bcspeval: yget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
212 : else
213 : ! if((yget < y(1)-0.5*zytol).or.(yget > y(ny)+0.5*zytol)) write(6,1012) yget,y(1),y(ny)
214 : ! 1012 format(' %bcspeval: yget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
215 0 : if(yget < y(1)) then
216 0 : zyget=y(1)
217 : else
218 0 : zyget=y(ny)
219 : end if
220 : end if
221 : end if
222 0 : if(ier /= 0) return
223 :
224 : ! now find interval in which target point lies..
225 :
226 0 : nxm=nx-1
227 0 : nym=ny-1
228 :
229 0 : if(ilinx == 1) then
230 0 : ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
231 0 : i=min(nxm, ii)
232 0 : if(zxget < x(i)) then
233 0 : i=i-1
234 0 : else if(zxget > x(i+1)) then
235 0 : i=i+1
236 : end if
237 : else
238 0 : if((1 <= i).and.(i < nxm)) then
239 0 : if((x(i) <= zxget).and.(zxget <= x(i+1))) then
240 : continue ! already have the zone
241 : else
242 0 : call zonfind(x,nx,zxget,i)
243 : end if
244 : else
245 0 : call zonfind(x,nx,zxget,i)
246 : end if
247 : end if
248 :
249 0 : if(iliny == 1) then
250 0 : jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
251 0 : j=min(nym, jj)
252 0 : if(zyget < y(j)) then
253 0 : j=j-1
254 0 : else if(zyget > y(j+1)) then
255 0 : j=j+1
256 : end if
257 : else
258 0 : if((1 <= j).and.(j < nym)) then
259 0 : if((y(j) <= zyget).and.(zyget <= y(j+1))) then
260 : continue ! already have the zone
261 : else
262 0 : call zonfind(y,ny,zyget,j)
263 : end if
264 : else
265 0 : call zonfind(y,ny,zyget,j)
266 : end if
267 : end if
268 :
269 0 : dx=zxget-x(i)
270 0 : dy=zyget-y(j)
271 :
272 0 : return
273 : end subroutine bcspevxy
274 :
275 :
276 : !------------------------------------------------------------------------
277 : ! bcspevfn -- OK now evaluate the bicubic spline
278 :
279 0 : subroutine bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,f,inf3,ny)
280 :
281 : ! input:
282 : integer ny
283 : integer ict(6) ! selector:
284 : ! ict(1)=1 for f (don't evaluate f if ict(1)=0)
285 : ! ict(2)=1 for df/dx ""
286 : ! ict(3)=1 for df/dy ""
287 : ! ict(4)=1 for d2f/dx2
288 : ! ict(5)=1 for d2f/dy2
289 : ! ict(6)=1 for d2f/dxdy
290 :
291 : ! note: if ict(1)=-1, evaluate f,d2f/dx2,d2f/dy2,d4f/dx2dy2
292 :
293 : ! the number of non zero values ict(1:6)
294 : ! determines the number of outputs...
295 :
296 : ! new dmc December 2005 -- access to higher derivatives (even if not
297 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
298 : ! if ict(1)=3 -- want 3rd derivatives
299 : ! ict(2)=1 for d3f/dx3
300 : ! ict(3)=1 for d3f/dx2dy
301 : ! ict(4)=1 for d3f/dxdy2
302 : ! ict(5)=1 for d3f/dy3
303 : ! number of non-zero values ict(2:5) gives no. of outputs
304 : ! if ict(1)=4 -- want 4th derivatives
305 : ! ict(2)=1 for d4f/dx3dy
306 : ! ict(3)=1 for d4f/dx2dy2
307 : ! ict(4)=1 for d4f/dxdy3
308 : ! number of non-zero values ict(2:4) gives no. of outputs
309 : ! if ict(1)=5 -- want 5th derivatives
310 : ! ict(2)=1 for d5f/dx3dy2
311 : ! ict(3)=1 for d5f/dx2dy3
312 : ! number of non-zero values ict(2:3) gives no. of outputs
313 : ! if ict(1)=6 -- want 6th derivatives
314 : ! d6f/dx3dy3 -- one value is returned.
315 :
316 : integer ivec,ivd ! vector dimensioning
317 :
318 : ! ivec-- number of vector pts (spline values to look up)
319 : ! ivd -- 1st dimension of fval, >= ivec
320 :
321 : ! output:
322 : real fval(ivd,6) ! output array
323 :
324 : ! v = index to element in vector;
325 : ! fval(v,1) = first item requested by ict(...),
326 : ! fval(v,2) = 2nd item requested, ...etc...
327 :
328 : ! input:
329 : !integer iv(ivec),jv(ivec) ! grid cell indices -- vectors
330 : !real dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors
331 : integer iv(:),jv(:) ! grid cell indices -- vectors
332 : real dxv(:),dyv(:) ! displacements w/in cell -- vectors
333 :
334 : integer inf3 ! 3rd dimension of f -- >= nx
335 : real f(:,:,:,:) ! (4,4,inf3,ny) ! bicubic fcn spline coeffs array
336 :
337 : ! usage example:
338 :
339 : ! 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs,
340 : ! find the x and y zone indices and displacements with respect
341 : ! to the "lower left corner" of the zone; store these in vectors
342 : ! iv,jv and dxv,dyv.
343 :
344 : ! 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only
345 : ! the 1st derivatives.
346 :
347 : ! 3. ivec is the length of the vector; ivd is the 1st dimension
348 : ! of the array fval to receive the output
349 :
350 : ! real fval(ivd,6)
351 : ! real xv(ivd),yv(ivd)
352 : ! integer iv(ivd),jv(ivd)
353 : ! real dxv(ivd),dyv(ivd)
354 : ! integer ict(6)
355 :
356 : ! real fspline(4,4,nx,ny) ! spline coeffs
357 : ! data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives
358 : ! ! only ... these will be output to
359 : ! ! fval(*,1) fval(*,2)
360 : ! ...
361 : ! do iv=1,ivec
362 : ! ... ! find indices and displacements
363 : ! end do
364 : ! call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny)
365 :
366 : !-------------------------------------------------------------------
367 : ! local:
368 :
369 : integer v ! vector element index
370 : integer iaval, i, j
371 0 : real dx, dy
372 :
373 : ! OK can now do evaluations
374 :
375 0 : iaval=0 ! fval addressing
376 :
377 0 : if(ict(1) <= 2) then
378 0 : if((ict(1) > 0).or.(ict(1) == -1)) then
379 : ! evaluate f
380 0 : iaval=iaval+1
381 0 : do v=1,ivec
382 0 : i=iv(v)
383 0 : j=jv(v)
384 0 : dx=dxv(v)
385 0 : dy=dyv(v)
386 : fval(v,iaval)=
387 : & f(1,1,i,j)+dy*(f(1,2,i,j)+dy*(f(1,3,i,j)+dy*f(1,4,i,j)))
388 : & +dx*(f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
389 : & +dx*(f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
390 0 : & +dx*(f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))))
391 : end do
392 : end if
393 :
394 0 : if((ict(2) > 0).and.(ict(1) /= -1)) then
395 : ! evaluate df/dx
396 0 : iaval=iaval+1
397 0 : do v=1,ivec
398 0 : i=iv(v)
399 0 : j=jv(v)
400 0 : dx=dxv(v)
401 0 : dy=dyv(v)
402 : fval(v,iaval)=
403 : & f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
404 : & +2.0*dx*(
405 : & f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
406 : & +1.5*dx*(
407 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j)))
408 0 : & ))
409 : end do
410 : end if
411 :
412 0 : if((ict(3) > 0).and.(ict(1) /= -1)) then
413 : ! evaluate df/dy
414 0 : iaval=iaval+1
415 0 : do v=1,ivec
416 0 : i=iv(v)
417 0 : j=jv(v)
418 0 : dx=dxv(v)
419 0 : dy=dyv(v)
420 : fval(v,iaval)=
421 : & f(1,2,i,j)+dy*(2.0*f(1,3,i,j)+dy*3.0*f(1,4,i,j))
422 : & +dx*(f(2,2,i,j)+dy*(2.0*f(2,3,i,j)+dy*3.0*f(2,4,i,j))
423 : & +dx*(f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j))
424 : & +dx*(f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j))
425 0 : & )))
426 : end do
427 : end if
428 :
429 0 : if((ict(4) > 0).or.(ict(1) == -1)) then
430 : ! evaluate d2f/dx2
431 0 : iaval=iaval+1
432 0 : do v=1,ivec
433 0 : i=iv(v)
434 0 : j=jv(v)
435 0 : dx=dxv(v)
436 0 : dy=dyv(v)
437 : fval(v,iaval)=
438 : & 2.0*(
439 : & f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j))))
440 : & +6.0*dx*(
441 0 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
442 : end do
443 : end if
444 :
445 0 : if((ict(5) > 0).or.(ict(1) == -1)) then
446 : ! evaluate d2f/dy2
447 0 : iaval=iaval+1
448 0 : do v=1,ivec
449 0 : i=iv(v)
450 0 : j=jv(v)
451 0 : dx=dxv(v)
452 0 : dy=dyv(v)
453 : fval(v,iaval)=
454 : & 2.0*f(1,3,i,j)+6.0*dy*f(1,4,i,j)
455 : & +dx*(2.0*f(2,3,i,j)+6.0*dy*f(2,4,i,j)
456 : & +dx*(2.0*f(3,3,i,j)+6.0*dy*f(3,4,i,j)
457 0 : & +dx*(2.0*f(4,3,i,j)+6.0*dy*f(4,4,i,j))))
458 : end do
459 : end if
460 :
461 0 : if((ict(6) > 0).and.(ict(1) /= -1)) then
462 : ! evaluate d2f/dxdy
463 0 : iaval=iaval+1
464 0 : do v=1,ivec
465 0 : i=iv(v)
466 0 : j=jv(v)
467 0 : dx=dxv(v)
468 0 : dy=dyv(v)
469 : fval(v,iaval)=
470 : & f(2,2,i,j)+dy*(2.0*f(2,3,i,j)+dy*3.0*f(2,4,i,j))
471 : & +2.*dx*(f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j))
472 : &+1.5*dx*(f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j))
473 0 : & ))
474 : end do
475 : end if
476 :
477 0 : if(ict(1) == -1) then
478 0 : iaval=iaval+1
479 0 : do v=1,ivec
480 0 : i=iv(v)
481 0 : j=jv(v)
482 0 : dx=dxv(v)
483 0 : dy=dyv(v)
484 : fval(v,iaval)=
485 : & 4.0*f(3,3,i,j)+12.0*dy*f(3,4,i,j)
486 0 : & +dx*(12.0*f(4,3,i,j)+36.0*dy*f(4,4,i,j))
487 : end do
488 : end if
489 :
490 : !-----------------------------------
491 : ! access to 3rd derivatives
492 :
493 0 : else if(ict(1) == 3) then
494 0 : if(ict(2) == 1) then
495 : ! evaluate d3f/dx3 (not continuous)
496 0 : iaval=iaval+1
497 0 : do v=1,ivec
498 0 : i=iv(v)
499 0 : j=jv(v)
500 0 : dy=dyv(v)
501 : fval(v,iaval)=
502 : & +6.0*(
503 0 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
504 : end do
505 : end if
506 :
507 0 : if(ict(3) == 1) then
508 : ! evaluate d3f/dx2dy
509 0 : iaval=iaval+1
510 0 : do v=1,ivec
511 0 : i=iv(v)
512 0 : j=jv(v)
513 0 : dx=dxv(v)
514 0 : dy=dyv(v)
515 : fval(v,iaval)=
516 : & 2.0*(
517 : & f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j)))
518 : & +6.0*dx*(
519 0 : & f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j)))
520 : end do
521 : end if
522 :
523 0 : if(ict(4) == 1) then
524 : ! evaluate d3f/dxdy2
525 0 : iaval=iaval+1
526 0 : do v=1,ivec
527 0 : i=iv(v)
528 0 : j=jv(v)
529 0 : dx=dxv(v)
530 0 : dy=dyv(v)
531 : fval(v,iaval)=
532 : & (2.0*f(2,3,i,j)+6.0*dy*f(2,4,i,j)
533 : & +2.0*dx*(2.0*f(3,3,i,j)+6.0*dy*f(3,4,i,j)
534 : & +1.5*dx*(2.0*f(4,3,i,j)+6.0*dy*f(4,4,i,j))
535 0 : & ))
536 : end do
537 : end if
538 :
539 0 : if(ict(5) == 1) then
540 : ! evaluate d3f/dy3 (not continuous)
541 0 : iaval=iaval+1
542 0 : do v=1,ivec
543 0 : i=iv(v)
544 0 : j=jv(v)
545 0 : dx=dxv(v)
546 : fval(v,iaval)=6.0*(f(1,4,i,j)+
547 0 : & dx*(f(2,4,i,j)+dx*(f(3,4,i,j)+dx*f(4,4,i,j))))
548 : end do
549 : end if
550 :
551 : !-----------------------------------
552 : ! access to 4th derivatives
553 :
554 0 : else if(ict(1) == 4) then
555 0 : if(ict(2) == 1) then
556 : ! evaluate d4f/dx3dy (not continuous)
557 0 : iaval=iaval+1
558 0 : do v=1,ivec
559 0 : i=iv(v)
560 0 : j=jv(v)
561 0 : dy=dyv(v)
562 : fval(v,iaval)=
563 : & +6.0*(
564 0 : & f(4,2,i,j)+dy*2.0*(f(4,3,i,j)+dy*1.5*f(4,4,i,j)))
565 : end do
566 : end if
567 :
568 0 : if(ict(3) == 1) then
569 : ! evaluate d4f/dx2dy2
570 0 : iaval=iaval+1
571 0 : do v=1,ivec
572 0 : i=iv(v)
573 0 : j=jv(v)
574 0 : dx=dxv(v)
575 0 : dy=dyv(v)
576 : fval(v,iaval)=
577 : & 4.0*f(3,3,i,j)+12.0*dy*f(3,4,i,j)
578 0 : & +dx*(12.0*f(4,3,i,j)+36.0*dy*f(4,4,i,j))
579 : end do
580 : end if
581 :
582 0 : if(ict(4) == 1) then
583 : ! evaluate d4f/dxdy3 (not continuous)
584 0 : iaval=iaval+1
585 0 : do v=1,ivec
586 0 : i=iv(v)
587 0 : j=jv(v)
588 0 : dx=dxv(v)
589 : fval(v,iaval)=
590 : & 6.0*(f(2,4,i,j)
591 0 : & +2.0*dx*(f(3,4,i,j)+1.5*dx*f(4,4,i,j)))
592 : end do
593 : end if
594 :
595 : !-----------------------------------
596 : ! access to 5th derivatives
597 :
598 0 : else if(ict(1) == 5) then
599 0 : if(ict(2) == 1) then
600 : ! evaluate d5f/dx3dy2 (not continuous)
601 0 : iaval=iaval+1
602 0 : do v=1,ivec
603 0 : i=iv(v)
604 0 : j=jv(v)
605 0 : dy=dyv(v)
606 : fval(v,iaval)=
607 0 : & +12.0*(f(4,3,i,j)+dy*3.0*f(4,4,i,j))
608 : end do
609 : end if
610 :
611 0 : if(ict(3) == 1) then
612 : ! evaluate d5f/dx3dy2 (not continuous)
613 0 : iaval=iaval+1
614 0 : do v=1,ivec
615 0 : i=iv(v)
616 0 : j=jv(v)
617 0 : dx=dxv(v)
618 : fval(v,iaval)=
619 0 : & 12.0*(f(3,4,i,j)+dx*3.0*f(4,4,i,j))
620 : end do
621 : end if
622 :
623 : !-----------------------------------
624 : ! access to 6th derivatives
625 :
626 0 : else if(ict(1) == 6) then
627 : ! evaluate d6f/dx3dy3 (not continuous)
628 0 : iaval=iaval+1
629 0 : do v=1,ivec
630 0 : i=iv(v)
631 0 : j=jv(v)
632 0 : fval(v,iaval) = 36.0*f(4,4,i,j)
633 : end do
634 : end if
635 :
636 0 : return
637 : end subroutine bcspevfn
638 : !----------------------
639 :
640 :
641 : ! cspline -- dmc 15 Feb 1999
642 :
643 : ! a standard interface to the 1d spline setup routine
644 : ! modified dmc 3 Mar 2000 -- to use Wayne Houlberg's v_spline code.
645 : ! new BC options added.
646 :
647 97 : subroutine cspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,wk,iwk,ilinx,ier)
648 :
649 : integer nx, iwk
650 : real x(nx) ! x axis (in)
651 : real fspl(4,nx) ! spline data (in/out)
652 : integer ibcxmin ! x(1) BC flag (in, see comments)
653 : real bcxmin ! x(1) BC data (in, see comments)
654 : integer ibcxmax ! x(nx) BC flag (in, see comments)
655 : real bcxmax ! x(nx) BC data (in, see comments)
656 : real wk(iwk) ! workspace of size at least nx
657 : integer ilinx ! even spacing flag (out)
658 : integer ier ! output, =0 means OK
659 :
660 : ! ** note wk(...) array is not used unless ibcxmin=-1 (periodic spline
661 : ! evaluation)
662 :
663 : ! this routine computes spline coefficients for a 1d spline --
664 : ! evaluation of the spline can be done by cspeval.for subroutines
665 : ! or directly by inline code.
666 :
667 : ! the input x axis x(1...nx) must be strictly ascending, i.e.
668 : ! x(i+1) > x(i) is required for i=1 to nx-1. This is checked and
669 : ! ier=1 is set and the routine exits if the test is not satisfied.
670 :
671 : ! on output, ilinx=1 is set if, to a reasonably close tolerance,
672 : ! all grid spacings x(i+1)-x(i) are equal. This allows a speedier
673 : ! grid lookup algorithm on evaluation of the spline. If on output
674 : ! ilinx=2, this means the spline x axis is not evenly spaced.
675 :
676 : ! the input data for the spline are given in f[j] = fspl(1,j). The
677 : ! output data are the spline coefficients fspl(2,j),fspl(3,j), and
678 : ! fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
679 : ! boundary conditions and with the properties
680 :
681 : ! s(x(j)) = fspl(1,j)
682 : ! s'(x) is continuous even at the grid points x(j)
683 : ! s''(x) is continuous even at the grid points x(j)
684 :
685 : ! the formula for evaluation of s(x) is:
686 :
687 : ! let dx = x-x(i), where x(i) <= x <= x(i+1). Then,
688 : ! s(x)=fspl(1,i) + dx*(fspl(2,i) +dx*(fspl(3,i) + dx*fspl(4,i)))
689 :
690 : ! ==>boundary conditions. Complete specification of a 1d spline
691 : ! requires specification of boundary conditions at x(1) and x(nx).
692 :
693 : ! this routine provides 4 options:
694 :
695 : ! -1 ***** PERIODIC BC
696 : ! ibcxmin=-1 -- periodic boundary condition. This means the
697 : ! boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
698 : ! are imposed. Note that s(x(1))=s(x(nx)) (i.e. fspl(1,1)=fspl(1,nx))
699 : ! is not required -- that is determined by the fspl array input data.
700 : ! The periodic boundary condition is to be preferred for periodic
701 : ! data. When splining periodic data f(x) with period P, the relation
702 : ! x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
703 : ! (ibcxmax, bcxmin, bcxmax are ignored).
704 :
705 : ! if a periodic boundary condition is set, this covers both boundaries.
706 : ! for the other types of boundary conditions, the type of condition
707 : ! chosen for the x(1) boundary need not be the same as the type chosen
708 : ! for the x(nx) boundary.
709 :
710 : ! 0 ***** NOT A KNOT BC
711 : ! ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
712 : ! condition -- see cubsplb.for. This is a common way for inferring
713 : ! a "good" spline boundary condition automatically from data in the
714 : ! vicinity of the boundary. (bcxmin | bcxmax are ignored).
715 :
716 : ! 1 ***** BC: SPECIFIED SLOPE
717 : ! ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
718 : ! s'(x(nx)) match the passed value (bcxmin | bcxmax).
719 :
720 : ! 2 ***** BC: SPECIFIED 2nd DERIVATIVE
721 : ! ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
722 : ! s''(x(nx)) match the passed value (bcxmin | bcxmax).
723 :
724 : ! 3 ***** BC: SPECIFIED SLOPE = 0.0
725 : ! ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
726 : ! s'(x(nx)) equal to ZERO.
727 :
728 : ! 4 ***** BC: SPECIFIED 2nd DERIVATIVE = 0.0
729 : ! ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
730 : ! s''(x(nx)) equal to ZERO.
731 :
732 : ! 5 ***** BC: 1st DIVIDED DIFFERENCE
733 : ! ibcxmin=5 | ibcxmax=5 -- boundary condition is to have s'(x(1)) |
734 : ! s'(x(nx)) equal to the slope from the 1st|last 2 points
735 :
736 : ! 6 ***** BC: 2nd DIVIDED DIFFERENCE
737 : ! ibcxmin=6 | ibcxmax=6 -- boundary condition is to have s''(x(1)) |
738 : ! s''(x(nx)) equal to the 2nd derivative from the 1st|last 3 points
739 :
740 : ! 7 ***** BC: 3rd DIVIDED DIFFERENCE
741 : ! ibcxmin=7 | ibcxmax=7 -- boundary condition is to have s'''(x(1)) |
742 : ! s'''(x(nx)) equal to the 3rd derivative from the 1st|last 4 points
743 :
744 : !---------------------------------------------------------------------
745 : real, parameter :: half = 0.5
746 : real, parameter :: sixth = 1.0/6.0
747 : integer::ierx,inum,i
748 :
749 : ! error checks
750 :
751 97 : ier = 0
752 97 : if(nx < 4) then
753 0 : write(6,'('' ?cspline: at least 4 x points required.'')')
754 0 : ier=1
755 : end if
756 97 : call ibc_ck(ibcxmin,'cspline','xmin',-1,7,ier)
757 97 : if(ibcxmin >= 0) call ibc_ck(ibcxmax,'cspline','xmax',0,7,ier)
758 :
759 : ! x axis check
760 :
761 97 : call splinck(x,nx,ilinx,1.0e-3,ierx)
762 97 : if(ierx /= 0) ier=2
763 :
764 97 : if(ier == 2) then
765 0 : write(6,'('' ?cspline: x axis not strict ascending'')')
766 : end if
767 :
768 97 : if(ibcxmin == -1) then
769 0 : inum=nx
770 0 : if(iwk < inum) then
771 0 : write(6,1009) inum,iwk,nx
772 : 1009 format(
773 : & ' ?cspline: workspace too small. need: ',i6,' got: ',i6/
774 : & ' (need = nx, nx=',i6)
775 0 : ier=3
776 : end if
777 : end if
778 :
779 97 : if(ier /= 0) return
780 :
781 : ! OK -- evaluate spline
782 :
783 97 : if(ibcxmin == 1) then
784 0 : fspl(2,1)=bcxmin
785 97 : else if(ibcxmin == 2) then
786 0 : fspl(3,1)=bcxmin
787 : end if
788 :
789 97 : if(ibcxmax == 1) then
790 0 : fspl(2,nx)=bcxmax
791 97 : else if(ibcxmax == 2) then
792 0 : fspl(3,nx)=bcxmax
793 : end if
794 :
795 97 : call v_spline(ibcxmin,ibcxmax,nx,x,fspl,wk)
796 :
797 3352 : do i=1,nx
798 3255 : fspl(3,i)=half*fspl(3,i)
799 3352 : fspl(4,i)=sixth*fspl(4,i)
800 : end do
801 :
802 : return
803 : end subroutine cspline
804 :
805 0 : subroutine evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,f1,inf2,ict,fval,ier)
806 :
807 : ! use mkbicub to set up spline coefficients!
808 :
809 : ! evaluate a 2d cubic Spline interpolant on a rectilinear
810 : ! grid -- this is C2 in both directions.
811 :
812 : ! this subroutine calls two subroutines:
813 : ! herm2xy -- find cell containing (xget,yget)
814 : ! fvbicub -- evaluate interpolant function and (optionally) derivatives
815 :
816 : ! input arguments:
817 : ! ================
818 :
819 : integer nx,ny ! grid sizes
820 : real xget,yget ! target of this interpolation
821 : real x(nx) ! ordered x grid
822 : real y(ny) ! ordered y grid
823 : integer ilinx ! ilinx=1 => assume x evenly spaced
824 : integer iliny ! iliny=1 => assume y evenly spaced
825 :
826 : integer inf2
827 : real, pointer :: f1(:) ! =(0:3,inf2,ny) interpolant data (cf "evbicub")
828 :
829 : ! f 2nd dimension inf2 must be >= nx
830 : ! contents of f:
831 :
832 : ! f(0,i,j) = f @ x(i),y(j)
833 : ! f(1,i,j) = d2f/dx2 @ x(i),y(j)
834 : ! f(2,i,j) = d2f/dy2 @ x(i),y(j)
835 : ! f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
836 :
837 : ! (these are spline coefficients selected for continuous 2-
838 : ! diffentiability, see mkbicub[w].for)
839 :
840 : integer ict(6) ! code specifying output desired
841 :
842 : ! ict(1)=1 -- return f (0, don't)
843 : ! ict(2)=1 -- return df/dx (0, don't)
844 : ! ict(3)=1 -- return df/dy (0, don't)
845 : ! ict(4)=1 -- return d2f/dx2 (0, don't)
846 : ! ict(5)=1 -- return d2f/dy2 (0, don't)
847 : ! ict(6)=1 -- return d2f/dxdy (0, don't)
848 : ! the number of non zero values ict(1:6)
849 : ! determines the number of outputs...
850 :
851 : ! new dmc December 2005 -- access to higher derivatives (even if not
852 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
853 : ! if ict(1)=3 -- want 3rd derivatives
854 : ! ict(2)=1 for d3f/dx3
855 : ! ict(3)=1 for d3f/dx2dy
856 : ! ict(4)=1 for d3f/dxdy2
857 : ! ict(5)=1 for d3f/dy3
858 : ! number of non-zero values ict(2:5) gives no. of outputs
859 : ! if ict(1)=4 -- want 4th derivatives
860 : ! ict(2)=1 for d4f/dx3dy
861 : ! ict(3)=1 for d4f/dx2dy2
862 : ! ict(4)=1 for d4f/dxdy3
863 : ! number of non-zero values ict(2:4) gives no. of outputs
864 : ! if ict(1)=5 -- want 5th derivatives
865 : ! ict(2)=1 for d5f/dx3dy2
866 : ! ict(3)=1 for d5f/dx2dy3
867 : ! number of non-zero values ict(2:3) gives no. of outputs
868 : ! if ict(1)=6 -- want 6th derivatives
869 : ! d6f/dx3dy3 -- one value is returned.
870 :
871 : ! output arguments:
872 : ! =================
873 :
874 : real fval(6) ! output data
875 : integer ier ! error code =0 ==> no error
876 :
877 : ! fval(1) receives the first output (depends on ict(...) spec)
878 : ! fval(2) receives the second output (depends on ict(...) spec)
879 : ! fval(3) receives the third output (depends on ict(...) spec)
880 : ! fval(4) receives the fourth output (depends on ict(...) spec)
881 : ! fval(5) receives the fourth output (depends on ict(...) spec)
882 : ! fval(6) receives the fourth output (depends on ict(...) spec)
883 :
884 : ! examples:
885 : ! on input ict = [1,1,1,0,0,1]
886 : ! on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
887 :
888 : ! on input ict = [1,0,0,0,0,0]
889 : ! on output fval= [f] ... elements 2 -- 6 never referenced.
890 :
891 : ! on input ict = [0,0,0,1,1,0]
892 : ! on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
893 :
894 : ! on input ict = [0,0,1,0,0,0]
895 : ! on output fval= [df/dy] ... elements 2 -- 6 never referenced.
896 :
897 : ! ier -- completion code: 0 means OK
898 : !-------------------
899 : ! local:
900 :
901 : integer i,j ! cell indices
902 :
903 : ! normalized displacement from (x(i),y(j)) corner of cell.
904 : ! xparam=0 @x(i) xparam=1 @x(i+1)
905 : ! yparam=0 @y(j) yparam=1 @y(j+1)
906 :
907 0 : real xparam,yparam
908 :
909 : ! cell dimensions and
910 : ! inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
911 :
912 0 : real hx,hy
913 0 : real hxi,hyi
914 :
915 : ! 0 <= xparam <= 1
916 : ! 0 <= yparam <= 1
917 :
918 : ! ** the interface is very similar to herm2ev.for; can use herm2xy **
919 : !---------------------------------------------------------------------
920 :
921 : call herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
922 0 : & i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
923 0 : if(ier /= 0) return
924 :
925 : call fvbicub(ict,1,1,
926 : & fval,(/i/),(/j/),(/xparam/),(/yparam/),
927 0 : & (/hx/),(/hxi/),(/hy/),(/hyi/),f1,inf2,ny)
928 :
929 0 : return
930 : end subroutine evbicub
931 :
932 : !---------------------------------------------------------------------
933 : ! evaluate C1 cubic Hermite function interpolation -- 2d fcn
934 : ! --vectorized-- dmc 10 Feb 1999
935 :
936 : ! use mkbicub to set up spline coefficients!
937 :
938 6 : subroutine fvbicub(ict,ivec,ivecd,
939 6 : & fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
940 : & f1,inf2,ny)
941 :
942 : integer ny
943 : integer ict(6) ! requested output control
944 : integer ivec ! vector length
945 : integer ivecd ! vector dimension (1st dim of fval)
946 :
947 : integer ii(ivec),jj(ivec) ! target cells (i,j)
948 : real xparam(ivec),yparam(ivec)
949 : ! normalized displacements from (i,j) corners
950 :
951 : real hx(ivec),hy(ivec) ! grid spacing, and
952 : real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
953 : ! & 1/(y(j+1)-y(j))
954 :
955 : integer inf2
956 : real, pointer :: f1(:) ! =(0:3,inf2,ny) interpolant data (cf "evbicub")
957 :
958 : real fval(ivecd,6) ! output returned
959 :
960 : ! for detailed description of fin, ict and fval see subroutine
961 : ! evbicub comments. Note ict is not vectorized; the same output
962 : ! is expected to be returned for all input vector data points.
963 :
964 : ! note that the index inputs ii,jj and parameter inputs
965 : ! xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
966 : ! output array fval has a vector ** 1st dimension ** whose
967 : ! size must be given as a separate argument
968 :
969 : ! to use this routine in scalar mode, pass in ivec=ivecd=1
970 :
971 : !---------------
972 : ! Spline evaluation consists of a "mixing" of the interpolant
973 : ! data using the linear functionals xparam, xpi = 1-xparam,
974 : ! yparam, ypi = 1-yparam, and the cubic functionals
975 : ! xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
976 : ! and their derivatives as needed.
977 :
978 : integer v,z36th,iadr,i,j
979 6 : real sum,xp,yp,xpi,ypi,xp2,yp2,xpi2,ypi2
980 6 : real cx,cy,cyd,cxi,cyi,cydi,hx2,hy2,cxd,cxdi
981 :
982 : real, parameter :: sixth = 0.166666666666666667
983 :
984 6 : real, pointer :: fin(:,:,:)
985 6 : fin(0:3,1:inf2,1:ny) => f1(1:4*inf2*ny)
986 :
987 :
988 : !---------------
989 : ! ...in x direction
990 :
991 6 : z36th=sixth*sixth
992 6 : iadr=0
993 :
994 6 : if(ict(1) <= 2) then
995 :
996 : ! get desired values:
997 :
998 6 : if(ict(1) == 1) then
999 :
1000 : ! function value:
1001 :
1002 12 : iadr=iadr+1
1003 12 : do v=1,ivec
1004 6 : i=ii(v)
1005 6 : j=jj(v)
1006 :
1007 : ! in x direction...
1008 :
1009 6 : xp=xparam(v)
1010 6 : xpi=1.0-xp
1011 6 : xp2=xp*xp
1012 6 : xpi2=xpi*xpi
1013 :
1014 6 : cx=xp*(xp2-1.0)
1015 6 : cxi=xpi*(xpi2-1.0)
1016 6 : hx2=hx(v)*hx(v)
1017 :
1018 : ! ...and in y direction
1019 :
1020 6 : yp=yparam(v)
1021 6 : ypi=1.0-yp
1022 6 : yp2=yp*yp
1023 6 : ypi2=ypi*ypi
1024 :
1025 6 : cy=yp*(yp2-1.0)
1026 6 : cyi=ypi*(ypi2-1.0)
1027 6 : hy2=hy(v)*hy(v)
1028 :
1029 : sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
1030 6 : & xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
1031 :
1032 : sum=sum+sixth*hx2*(
1033 : & cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1034 6 : & cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1035 :
1036 : sum=sum+sixth*hy2*(
1037 : & xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
1038 6 : & xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
1039 :
1040 : sum=sum+z36th*hx2*hy2*(
1041 : & cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1042 6 : & cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1043 :
1044 12 : fval(v,iadr)=sum
1045 : end do
1046 : end if
1047 :
1048 6 : if(ict(2) == 1) then
1049 :
1050 : ! df/dx:
1051 :
1052 6 : iadr=iadr+1
1053 12 : do v=1,ivec
1054 6 : i=ii(v)
1055 6 : j=jj(v)
1056 :
1057 : ! in x direction...
1058 :
1059 6 : xp=xparam(v)
1060 6 : xpi=1.0-xp
1061 6 : xp2=xp*xp
1062 6 : xpi2=xpi*xpi
1063 :
1064 6 : cxd=3.0*xp2-1.0
1065 6 : cxdi=-3.0*xpi2+1.0
1066 :
1067 : ! ...and in y direction
1068 :
1069 6 : yp=yparam(v)
1070 6 : ypi=1.0-yp
1071 6 : yp2=yp*yp
1072 6 : ypi2=ypi*ypi
1073 :
1074 6 : cy=yp*(yp2-1.0)
1075 6 : cyi=ypi*(ypi2-1.0)
1076 6 : hy2=hy(v)*hy(v)
1077 :
1078 : sum=hxi(v)*(
1079 : & -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
1080 6 : & +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
1081 :
1082 : sum=sum+sixth*hx(v)*(
1083 : & cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1084 6 : & cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1085 :
1086 : sum=sum+sixth*hxi(v)*hy2*(
1087 : & -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
1088 6 : & +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
1089 :
1090 : sum=sum+z36th*hx(v)*hy2*(
1091 : & cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1092 6 : & cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1093 :
1094 12 : fval(v,iadr)=sum
1095 : end do
1096 : end if
1097 :
1098 6 : if(ict(3) == 1) then
1099 :
1100 : ! df/dy:
1101 :
1102 6 : iadr=iadr+1
1103 12 : do v=1,ivec
1104 6 : i=ii(v)
1105 6 : j=jj(v)
1106 :
1107 : ! in x direction...
1108 :
1109 6 : xp=xparam(v)
1110 6 : xpi=1.0-xp
1111 6 : xp2=xp*xp
1112 6 : xpi2=xpi*xpi
1113 :
1114 6 : cx=xp*(xp2-1.0)
1115 6 : cxi=xpi*(xpi2-1.0)
1116 6 : hx2=hx(v)*hx(v)
1117 :
1118 : ! ...and in y direction
1119 :
1120 6 : yp=yparam(v)
1121 6 : ypi=1.0-yp
1122 6 : yp2=yp*yp
1123 6 : ypi2=ypi*ypi
1124 :
1125 6 : cyd=3.0*yp2-1.0
1126 6 : cydi=-3.0*ypi2+1.0
1127 :
1128 : sum=hyi(v)*(
1129 : & xpi*(-fin(0,i,j) +fin(0,i,j+1))+
1130 6 : & xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
1131 :
1132 : sum=sum+sixth*hx2*hyi(v)*(
1133 : & cxi*(-fin(1,i,j) +fin(1,i,j+1))+
1134 6 : & cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1135 :
1136 : sum=sum+sixth*hy(v)*(
1137 : & xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
1138 6 : & xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1139 :
1140 : sum=sum+z36th*hx2*hy(v)*(
1141 : & cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1142 6 : & cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1143 :
1144 12 : fval(v,iadr)=sum
1145 : end do
1146 : end if
1147 :
1148 6 : if(ict(4) == 1) then
1149 :
1150 : ! d2f/dx2:
1151 :
1152 0 : iadr=iadr+1
1153 0 : do v=1,ivec
1154 0 : i=ii(v)
1155 0 : j=jj(v)
1156 :
1157 : ! in x direction...
1158 :
1159 0 : xp=xparam(v)
1160 0 : xpi=1.0-xp
1161 :
1162 : ! ...and in y direction
1163 :
1164 0 : yp=yparam(v)
1165 0 : ypi=1.0-yp
1166 0 : yp2=yp*yp
1167 0 : ypi2=ypi*ypi
1168 :
1169 0 : cy=yp*(yp2-1.0)
1170 0 : cyi=ypi*(ypi2-1.0)
1171 0 : hy2=hy(v)*hy(v)
1172 :
1173 : sum=(
1174 : & xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1175 0 : & xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1176 :
1177 : sum=sum+sixth*hy2*(
1178 : & xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1179 0 : & xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1180 :
1181 0 : fval(v,iadr)=sum
1182 : end do
1183 : end if
1184 :
1185 6 : if(ict(5) == 1) then
1186 :
1187 : ! d2f/dy2:
1188 :
1189 0 : iadr=iadr+1
1190 0 : do v=1,ivec
1191 0 : i=ii(v)
1192 0 : j=jj(v)
1193 :
1194 : ! in x direction...
1195 :
1196 0 : xp=xparam(v)
1197 0 : xpi=1.0-xp
1198 0 : xp2=xp*xp
1199 0 : xpi2=xpi*xpi
1200 :
1201 0 : cx=xp*(xp2-1.0)
1202 0 : cxi=xpi*(xpi2-1.0)
1203 0 : hx2=hx(v)*hx(v)
1204 :
1205 : ! ...and in y direction
1206 :
1207 0 : yp=yparam(v)
1208 0 : ypi=1.0-yp
1209 :
1210 : sum=(
1211 : & xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
1212 0 : & xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1213 :
1214 : sum=sum+sixth*hx2*(
1215 : & cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1216 0 : & cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1217 :
1218 0 : fval(v,iadr)=sum
1219 : end do
1220 : end if
1221 :
1222 6 : if(ict(6) == 1) then
1223 :
1224 : ! d2f/dxdy:
1225 :
1226 0 : iadr=iadr+1
1227 0 : do v=1,ivec
1228 0 : i=ii(v)
1229 0 : j=jj(v)
1230 :
1231 : ! in x direction...
1232 :
1233 0 : xp=xparam(v)
1234 0 : xpi=1.0-xp
1235 0 : xp2=xp*xp
1236 0 : xpi2=xpi*xpi
1237 :
1238 0 : cxd=3.0*xp2-1.0
1239 0 : cxdi=-3.0*xpi2+1.0
1240 :
1241 : ! ...and in y direction
1242 :
1243 0 : yp=yparam(v)
1244 0 : ypi=1.0-yp
1245 0 : yp2=yp*yp
1246 0 : ypi2=ypi*ypi
1247 :
1248 0 : cyd=3.0*yp2-1.0
1249 0 : cydi=-3.0*ypi2+1.0
1250 :
1251 : sum=hxi(v)*hyi(v)*(
1252 : & fin(0,i,j) -fin(0,i,j+1)
1253 0 : & -fin(0,i+1,j)+fin(0,i+1,j+1))
1254 :
1255 : sum=sum+sixth*hx(v)*hyi(v)*(
1256 : & cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
1257 0 : & cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1258 :
1259 : sum=sum+sixth*hxi(v)*hy(v)*(
1260 : & -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
1261 0 : & +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1262 :
1263 : sum=sum+z36th*hx(v)*hy(v)*(
1264 : & cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1265 0 : & cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1266 :
1267 0 : fval(v,iadr)=sum
1268 : end do
1269 : end if
1270 :
1271 : !-------------------------------------------------
1272 :
1273 0 : else if(ict(1) == 3) then
1274 0 : if(ict(2) == 1) then
1275 : ! evaluate d3f/dx3 (not continuous)
1276 0 : iadr=iadr+1
1277 0 : do v=1,ivec
1278 0 : i=ii(v)
1279 0 : j=jj(v)
1280 0 : yp=yparam(v)
1281 0 : ypi=1.0-yp
1282 0 : yp2=yp*yp
1283 0 : ypi2=ypi*ypi
1284 0 : cy=yp*(yp2-1.0)
1285 0 : cyi=ypi*(ypi2-1.0)
1286 0 : hy2=hy(v)*hy(v)
1287 : sum=hxi(v)*(
1288 : & -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
1289 0 : & +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1290 :
1291 : sum=sum+sixth*hy2*hxi(v)*(
1292 : & -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
1293 0 : & +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1294 :
1295 0 : fval(v,iadr)=sum
1296 : end do
1297 : end if
1298 :
1299 0 : if(ict(3) == 1) then
1300 : ! evaluate d3f/dx2dy
1301 0 : iadr=iadr+1
1302 0 : do v=1,ivec
1303 0 : i=ii(v)
1304 0 : j=jj(v)
1305 0 : xp=xparam(v)
1306 0 : xpi=1.0-xp
1307 0 : yp=yparam(v)
1308 0 : ypi=1.0-yp
1309 0 : yp2=yp*yp
1310 0 : ypi2=ypi*ypi
1311 0 : cyd=3.0*yp2-1.0
1312 0 : cydi=-3.0*ypi2+1.0
1313 :
1314 : sum=hyi(v)*(
1315 : & xpi*(-fin(1,i,j) +fin(1,i,j+1))+
1316 0 : & xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
1317 :
1318 : sum=sum+sixth*hy(v)*(
1319 : & xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1320 0 : & xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1321 :
1322 0 : fval(v,iadr)=sum
1323 : end do
1324 : end if
1325 :
1326 0 : if(ict(4) == 1) then
1327 : ! evaluate d3f/dxdy2
1328 0 : iadr=iadr+1
1329 0 : do v=1,ivec
1330 0 : i=ii(v)
1331 0 : j=jj(v)
1332 0 : xp=xparam(v)
1333 0 : xpi=1.0-xp
1334 0 : xp2=xp*xp
1335 0 : xpi2=xpi*xpi
1336 0 : cxd=3.0*xp2-1.0
1337 0 : cxdi=-3.0*xpi2+1.0
1338 0 : yp=yparam(v)
1339 0 : ypi=1.0-yp
1340 :
1341 : sum=hxi(v)*(
1342 : & -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
1343 0 : & +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1344 :
1345 : sum=sum+sixth*hx(v)*(
1346 : & cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1347 0 : & cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1348 :
1349 0 : fval(v,iadr)=sum
1350 : end do
1351 : end if
1352 :
1353 0 : if(ict(5) == 1) then
1354 : ! evaluate d3f/dy3 (not continuous)
1355 0 : iadr=iadr+1
1356 0 : do v=1,ivec
1357 0 : i=ii(v)
1358 0 : j=jj(v)
1359 :
1360 0 : xp=xparam(v)
1361 0 : xpi=1.0-xp
1362 0 : xp2=xp*xp
1363 0 : xpi2=xpi*xpi
1364 :
1365 0 : cx=xp*(xp2-1.0)
1366 0 : cxi=xpi*(xpi2-1.0)
1367 0 : hx2=hx(v)*hx(v)
1368 :
1369 : sum=hyi(v)*(
1370 : & xpi*(-fin(2,i,j) +fin(2,i,j+1))+
1371 0 : & xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
1372 :
1373 : sum=sum+sixth*hx2*hyi(v)*(
1374 : & cxi*(-fin(3,i,j) +fin(3,i,j+1))+
1375 0 : & cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1376 :
1377 0 : fval(v,iadr)=sum
1378 : end do
1379 : end if
1380 :
1381 : !-----------------------------------
1382 : ! access to 4th derivatives
1383 :
1384 0 : else if(ict(1) == 4) then
1385 0 : if(ict(2) == 1) then
1386 : ! evaluate d4f/dx3dy (not continuous)
1387 0 : iadr=iadr+1
1388 0 : do v=1,ivec
1389 0 : i=ii(v)
1390 0 : j=jj(v)
1391 0 : yp=yparam(v)
1392 0 : ypi=1.0-yp
1393 0 : yp2=yp*yp
1394 0 : ypi2=ypi*ypi
1395 0 : cyd=3.0*yp2-1.0
1396 0 : cydi=-3.0*ypi2+1.0
1397 :
1398 : sum=hxi(v)*hyi(v)*(
1399 : & +( fin(1,i,j) -fin(1,i,j+1))
1400 0 : & +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1401 :
1402 : sum=sum+sixth*hy(v)*hxi(v)*(
1403 : & -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
1404 0 : & +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1405 :
1406 0 : fval(v,iadr)=sum
1407 : end do
1408 : end if
1409 :
1410 0 : if(ict(3) == 1) then
1411 : ! evaluate d4f/dx2dy2
1412 0 : iadr=iadr+1
1413 0 : do v=1,ivec
1414 0 : i=ii(v)
1415 0 : j=jj(v)
1416 :
1417 0 : xp=xparam(v)
1418 0 : xpi=1.0-xp
1419 0 : yp=yparam(v)
1420 0 : ypi=1.0-yp
1421 :
1422 : sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1423 0 : & xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
1424 :
1425 0 : fval(v,iadr)=sum
1426 : end do
1427 : end if
1428 :
1429 0 : if(ict(4) == 1) then
1430 : ! evaluate d4f/dxdy3 (not continuous)
1431 0 : iadr=iadr+1
1432 0 : do v=1,ivec
1433 0 : i=ii(v)
1434 0 : j=jj(v)
1435 :
1436 0 : xp=xparam(v)
1437 0 : xpi=1.0-xp
1438 0 : xp2=xp*xp
1439 0 : xpi2=xpi*xpi
1440 :
1441 0 : cxd=3.0*xp2-1.0
1442 0 : cxdi=-3.0*xpi2+1.0
1443 :
1444 : sum=hyi(v)*hxi(v)*(
1445 : & +( fin(2,i,j) -fin(2,i,j+1))
1446 0 : & +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
1447 :
1448 : sum=sum+sixth*hx(v)*hyi(v)*(
1449 : & cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
1450 0 : & cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1451 :
1452 0 : fval(v,iadr)=sum
1453 : end do
1454 : end if
1455 :
1456 : !-----------------------------------
1457 : ! access to 5th derivatives
1458 :
1459 0 : else if(ict(1) == 5) then
1460 0 : if(ict(2) == 1) then
1461 : ! evaluate d5f/dx3dy2 (not continuous)
1462 0 : iadr=iadr+1
1463 0 : do v=1,ivec
1464 0 : i=ii(v)
1465 0 : j=jj(v)
1466 :
1467 0 : yp=yparam(v)
1468 0 : ypi=1.0-yp
1469 :
1470 : sum=hxi(v)*(
1471 : & -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
1472 0 : & +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1473 :
1474 0 : fval(v,iadr)=sum
1475 : end do
1476 : end if
1477 :
1478 0 : if(ict(3) == 1) then
1479 : ! evaluate d5f/dx2dy3 (not continuous)
1480 0 : iadr=iadr+1
1481 0 : do v=1,ivec
1482 0 : i=ii(v)
1483 0 : j=jj(v)
1484 :
1485 0 : xp=xparam(v)
1486 0 : xpi=1.0-xp
1487 :
1488 : sum=hyi(v)*(
1489 : & xpi*(-fin(3,i,j) +fin(3,i,j+1))+
1490 0 : & xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1491 :
1492 0 : fval(v,iadr)=sum
1493 : end do
1494 : end if
1495 :
1496 : !-----------------------------------
1497 : ! access to 6th derivatives
1498 :
1499 0 : else if(ict(1) == 6) then
1500 : ! evaluate d6f/dx3dy3 (not continuous)
1501 0 : iadr=iadr+1
1502 0 : do v=1,ivec
1503 0 : i=ii(v)
1504 0 : j=jj(v)
1505 : sum=hxi(v)*hyi(v)*(
1506 : & +( fin(3,i,j) -fin(3,i,j+1))
1507 0 : & +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1508 0 : fval(v,iadr)=sum
1509 : end do
1510 : end if
1511 :
1512 6 : return
1513 6 : end subroutine fvbicub
1514 :
1515 :
1516 : !---------------------------------------------------------------------
1517 : ! herm2xy -- look up x-y zone
1518 :
1519 : ! this is the "first part" of herm2ev, see comments, above.
1520 :
1521 6 : subroutine herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
1522 : & i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
1523 :
1524 : ! input of herm2xy
1525 : ! ================
1526 :
1527 : integer nx,ny ! array dimensions
1528 :
1529 : real xget,yget ! target point
1530 : real x(:) ! (nx) ! indep. coords., strict ascending
1531 : real y(:) ! (ny) ! indep. coords., strict ascending
1532 :
1533 : integer ilinx ! =1: x evenly spaced
1534 : integer iliny ! =1: y evenly spaced
1535 :
1536 : ! output of herm2xy
1537 : ! =================
1538 : integer i,j ! index to cell containing target pt
1539 : ! on exit: 1 <= i <= nx-1 1 <= j <= ny-1
1540 :
1541 : ! normalized position w/in (i,j) cell (btw 0 and 1):
1542 :
1543 : real xparam ! (xget-x(i))/(x(i+1)-x(i))
1544 : real yparam ! (yget-y(j))/(y(j+1)-y(j))
1545 :
1546 : ! grid spacing
1547 :
1548 : real hx ! hx = x(i+1)-x(i)
1549 : real hy ! hy = y(j+1)-y(j)
1550 :
1551 : ! inverse grid spacing:
1552 :
1553 : real hxi ! 1/hx = 1/(x(i+1)-x(i))
1554 : real hyi ! 1/hy = 1/(y(j+1)-y(j))
1555 :
1556 : integer ier ! return ier /= 0 on error
1557 6 : real zxget,zyget,zxtol,zytol
1558 : integer nxm,nym,ii,jj
1559 :
1560 : !------------------------------------
1561 :
1562 6 : ier=0
1563 :
1564 : ! range check
1565 :
1566 6 : zxget=xget
1567 6 : zyget=yget
1568 6 : if((xget < x(1)).or.(xget > x(nx))) then
1569 0 : zxtol=4.0e-7*max(abs(x(1)),abs(x(nx)))
1570 0 : if((xget < x(1)-zxtol).or.(xget > x(nx)+zxtol)) then
1571 0 : ier=1
1572 : ! write(6,1001) xget,x(1),x(nx)
1573 : ! 1001 format(' ?herm2ev: xget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
1574 : else
1575 : ! if((xget < x(1)-0.5*zxtol).or. (xget > x(nx)+0.5*zxtol)) write(6,1011) xget,x(1),x(nx)
1576 : ! 1011 format(' %herm2ev: xget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
1577 0 : if(xget < x(1)) then
1578 0 : zxget=x(1)
1579 : else
1580 0 : zxget=x(nx)
1581 : end if
1582 : end if
1583 : end if
1584 6 : if((yget < y(1)).or.(yget > y(ny))) then
1585 0 : zytol=4.0e-7*max(abs(y(1)),abs(y(ny)))
1586 0 : if((yget < y(1)-zytol).or.(yget > y(ny)+zytol)) then
1587 0 : ier=1
1588 : ! write(6,1002) yget,y(1),y(ny)
1589 : ! 1002 format(' ?herm2ev: yget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
1590 : else
1591 : ! if((yget < y(1)-0.5*zytol).or. (yget > y(ny)+0.5*zytol)) write(6,1012) yget,y(1),y(ny)
1592 : ! 1012 format(' %herm2ev: yget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
1593 0 : if(yget < y(1)) then
1594 0 : zyget=y(1)
1595 : else
1596 0 : zyget=y(ny)
1597 : end if
1598 : end if
1599 : end if
1600 6 : if(ier /= 0) return
1601 :
1602 : ! now find interval in which target point lies..
1603 :
1604 6 : nxm=nx-1
1605 6 : nym=ny-1
1606 :
1607 6 : if(ilinx == 1) then
1608 6 : ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
1609 6 : i=min(nxm, ii)
1610 6 : if(zxget < x(i)) then
1611 0 : i=i-1
1612 6 : else if(zxget > x(i+1)) then
1613 0 : i=i+1
1614 : end if
1615 : else
1616 0 : if((1 <= i).and.(i < nxm)) then
1617 0 : if((x(i) <= zxget).and.(zxget <= x(i+1))) then
1618 : continue ! already have the zone
1619 : else
1620 0 : call zonfind(x,nx,zxget,i)
1621 : end if
1622 : else
1623 0 : call zonfind(x,nx,zxget,i)
1624 : end if
1625 : end if
1626 :
1627 6 : if(iliny == 1) then
1628 6 : jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
1629 6 : j=min(nym, jj)
1630 6 : if(zyget < y(j)) then
1631 0 : j=j-1
1632 6 : else if(zyget > y(j+1)) then
1633 0 : j=j+1
1634 : end if
1635 : else
1636 0 : if((1 <= j).and.(j < nym)) then
1637 0 : if((y(j) <= zyget).and.(zyget <= y(j+1))) then
1638 : continue ! already have the zone
1639 : else
1640 0 : call zonfind(y,ny,zyget,j)
1641 : end if
1642 : else
1643 0 : call zonfind(y,ny,zyget,j)
1644 : end if
1645 : end if
1646 :
1647 6 : hx=(x(i+1)-x(i))
1648 6 : hy=(y(j+1)-y(j))
1649 :
1650 6 : hxi=1.0/hx
1651 6 : hyi=1.0/hy
1652 :
1653 6 : xparam=(zxget-x(i))*hxi
1654 6 : yparam=(zyget-y(j))*hyi
1655 :
1656 6 : return
1657 : end subroutine herm2xy
1658 :
1659 0 : subroutine herm2fcn(ict,ivec,ivecd,
1660 0 : & fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
1661 0 : & fin,inf2,ny)
1662 :
1663 : integer ny, inf2
1664 : integer ict(4) ! requested output control
1665 : integer ivec ! vector length
1666 : integer ivecd ! vector dimension (1st dim of fval)
1667 :
1668 : integer ii(ivec),jj(ivec) ! target cells (i,j)
1669 : real xparam(ivec),yparam(ivec)
1670 : ! normalized displacements from (i,j) corners
1671 :
1672 : real hx(ivec),hy(ivec) ! grid spacing, and
1673 : real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
1674 : ! & 1/(y(j+1)-y(j))
1675 :
1676 : real fin(0:3,inf2,ny) ! interpolant data (cf "herm2ev")
1677 :
1678 : real fval(ivecd,4) ! output returned
1679 :
1680 : ! for detailed description of fin, ict and fval see subroutine
1681 : ! herm2ev comments. Note ict is not vectorized; the same output
1682 : ! is expected to be returned for all input vector data points.
1683 :
1684 : ! note that the index inputs ii,jj and parameter inputs
1685 : ! xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
1686 : ! output array fval has a vector ** 1st dimension ** whose
1687 : ! size must be given as a separate argument
1688 :
1689 : ! to use this routine in scalar mode, pass in ivec=ivecd=1
1690 :
1691 : !---------------
1692 : ! Hermite cubic basis functions
1693 : ! -->for function value matching
1694 : ! a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
1695 : ! abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
1696 :
1697 : ! a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
1698 : ! abar(x)=1-a(x)
1699 : ! a'(x)=-abar'(x) = 6.0*x*(1.0-x)
1700 :
1701 : ! -->for derivative matching
1702 : ! b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
1703 : ! bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
1704 :
1705 : ! b(x)=x**3-x**2 b'(x)=3*x**2-2*x
1706 : ! bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
1707 :
1708 0 : real sum
1709 : integer v,iadr,i,j
1710 0 : real xp,yp,xpi,ypi,xp2,yp2,xpi2,ypi2
1711 0 : real ax,bx,axbar,bxbar,ay,by,aybar,bybar
1712 0 : real axp,axbarp,bxp,bxbarp,ayp,aybarp,bybarp,byp
1713 :
1714 0 : do v=1,ivec
1715 0 : i=ii(v)
1716 0 : j=jj(v)
1717 :
1718 : ! ...in x direction
1719 :
1720 0 : xp=xparam(v)
1721 0 : xpi=1.0-xp
1722 0 : xp2=xp*xp
1723 0 : xpi2=xpi*xpi
1724 0 : ax=xp2*(3.0-2.0*xp)
1725 0 : axbar=1.0-ax
1726 0 : bx=-xp2*xpi
1727 0 : bxbar=xpi2*xp
1728 :
1729 : ! ...in y direction
1730 :
1731 0 : yp=yparam(v)
1732 0 : ypi=1.0-yp
1733 0 : yp2=yp*yp
1734 0 : ypi2=ypi*ypi
1735 0 : ay=yp2*(3.0-2.0*yp)
1736 0 : aybar=1.0-ay
1737 0 : by=-yp2*ypi
1738 0 : bybar=ypi2*yp
1739 :
1740 : ! ...derivatives...
1741 :
1742 0 : axp=6.0*xp*xpi
1743 0 : axbarp=-axp
1744 0 : bxp=xp*(3.0*xp-2.0)
1745 0 : bxbarp=xpi*(3.0*xpi-2.0)
1746 :
1747 0 : ayp=6.0*yp*ypi
1748 0 : aybarp=-ayp
1749 0 : byp=yp*(3.0*yp-2.0)
1750 0 : bybarp=ypi*(3.0*ypi-2.0)
1751 :
1752 0 : iadr=0
1753 :
1754 : ! get desired values:
1755 :
1756 0 : if(ict(1) == 1) then
1757 :
1758 : ! function value:
1759 :
1760 0 : iadr=iadr+1
1761 : sum=axbar*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
1762 0 : & ax*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1))
1763 :
1764 : sum=sum+hx(v)*(
1765 : & bxbar*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
1766 0 : & bx*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1)))
1767 :
1768 : sum=sum+hy(v)*(
1769 : & axbar*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
1770 0 : & ax*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
1771 :
1772 : sum=sum+hx(v)*hy(v)*(
1773 : & bxbar*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
1774 0 : & bx*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
1775 :
1776 0 : fval(v,iadr)=sum
1777 : end if
1778 :
1779 0 : if(ict(2) == 1) then
1780 :
1781 : ! df/dx:
1782 :
1783 0 : iadr=iadr+1
1784 :
1785 : sum=hxi(v)*(
1786 : & axbarp*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
1787 0 : & axp*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1)))
1788 :
1789 : sum=sum+
1790 : & bxbarp*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
1791 0 : & bxp*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1))
1792 :
1793 : sum=sum+hxi(v)*hy(v)*(
1794 : & axbarp*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
1795 0 : & axp*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
1796 :
1797 : sum=sum+hy(v)*(
1798 : & bxbarp*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
1799 0 : & bxp*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
1800 :
1801 0 : fval(v,iadr)=sum
1802 : end if
1803 :
1804 0 : if(ict(3) == 1) then
1805 :
1806 : ! df/dy:
1807 :
1808 0 : iadr=iadr+1
1809 :
1810 : sum=hyi(v)*(
1811 : & axbar*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
1812 0 : & ax*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
1813 :
1814 : sum=sum+hx(v)*hyi(v)*(
1815 : & bxbar*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
1816 0 : & bx*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
1817 :
1818 : sum=sum+
1819 : & axbar*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
1820 0 : & ax*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1))
1821 :
1822 : sum=sum+hx(v)*(
1823 : & bxbar*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
1824 0 : & bx*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1)))
1825 :
1826 0 : fval(v,iadr)=sum
1827 : end if
1828 :
1829 0 : if(ict(4) == 1) then
1830 :
1831 : ! d2f/dxdy:
1832 :
1833 0 : iadr=iadr+1
1834 :
1835 : sum=hxi(v)*hyi(v)*(
1836 : & axbarp*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
1837 0 : & axp*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
1838 :
1839 : sum=sum+hyi(v)*(
1840 : & bxbarp*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
1841 0 : & bxp*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
1842 :
1843 : sum=sum+hxi(v)*(
1844 : & axbarp*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
1845 0 : & axp*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1)))
1846 :
1847 : sum=sum+
1848 : & bxbarp*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
1849 0 : & bxp*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1))
1850 :
1851 0 : fval(v,iadr)=sum
1852 : end if
1853 :
1854 : end do ! vector loop
1855 :
1856 0 : return
1857 : end subroutine herm2fcn
1858 :
1859 :
1860 198 : subroutine ibc_ck(ibc,slbl,xlbl,imin,imax,ier)
1861 :
1862 : ! check that spline routine ibc flag is in range
1863 :
1864 : integer ibc ! input -- flag value
1865 : character*(*) slbl ! input -- subroutine name
1866 : character*(*) xlbl ! input -- axis label
1867 :
1868 : integer imin ! input -- min allowed value
1869 : integer imax ! input -- max allowed value
1870 :
1871 : integer ier ! output -- set =1 if error detected
1872 :
1873 : !----------------------
1874 :
1875 198 : if((ibc < imin).or.(ibc > imax)) then
1876 0 : ier=1
1877 : ! write(6,1001) slbl,xlbl,ibc,imin,imax
1878 : ! 1001 format(' ?',a,' -- ibc',a,' = ',i9,' out of range ', i2,' to ',i2)
1879 : end if
1880 :
1881 0 : return
1882 : end subroutine ibc_ck
1883 :
1884 :
1885 1 : subroutine do_mkbicub(x,nx,y,ny,f1,nf2,
1886 1 : & ibcxmin,bcxmin,ibcxmax,bcxmax,
1887 1 : & ibcymin,bcymin,ibcymax,bcymax,
1888 : & ilinx,iliny,ier)
1889 :
1890 : ! setup bicubic spline, dynamic allocation of workspace
1891 : ! fortran-90 fixed form source
1892 :
1893 : ! --NOTE-- dmc 22 Feb 2004 -- rewrite for direct calculation of
1894 : ! coefficients, to avoid large transient use of memory.
1895 :
1896 : !
1897 : implicit NONE
1898 :
1899 : ! input:
1900 : integer nx ! length of x vector
1901 : integer ny ! length of y vector
1902 : real x(:) ! (nx) ! x vector, strict ascending
1903 : real y(:) ! (ny) ! y vector, strict ascending
1904 :
1905 : integer nf2 ! 2nd dimension of f, nf2 >= nx
1906 : ! input/output:
1907 : real, pointer :: f1(:) ! =(4,nf2,ny) ! data & spline coefficients
1908 :
1909 : ! on input: f(1,i,j) = f(x(i),y(j))
1910 : ! on output: f(1,i,j) unchanged
1911 : ! f(2,i,j) = d2f/dx2(x(i),y(j))
1912 : ! f(3,i,j) = d2f/dy2(x(i),y(j))
1913 : ! f(4,i,j) = d4f/dx2dy2(x(i),y(j))
1914 :
1915 : ! and the interpolation formula for (x,y) in (x(i),x(i+1))^(y(j),y(j+1))
1916 : ! is:
1917 : ! hx = x(i+1)-x(i) hy = y(j+1)-y(j)
1918 : ! dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
1919 : ! dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
1920 : ! dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
1921 :
1922 : ! finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
1923 : ! +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
1924 : ! +1/6*hx**2*
1925 : ! dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
1926 : ! +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
1927 : ! +1/6*hy**2*
1928 : ! dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
1929 : ! +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
1930 : ! +1/36*hx**2*hy**2*
1931 : ! dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
1932 : ! +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
1933 :
1934 : ! where the f(2:4,*,*) are cleverly evaluated to assure
1935 : ! (a) finterp is continuous and twice differentiable across all
1936 : ! grid cell boundaries, and
1937 : ! (b) all boundary conditions are satisfied.
1938 :
1939 : ! input bdy condition data:
1940 : integer ibcxmin ! bc flag for x=xmin
1941 : real bcxmin(:) ! (ny) ! bc data vs. y at x=xmin
1942 : integer ibcxmax ! bc flag for x=xmax
1943 : real bcxmax(:) ! (ny) ! bc data vs. y at x=xmax
1944 :
1945 : integer ibcymin ! bc flag for y=ymin
1946 : real bcymin(:) ! (nx) ! bc data vs. x at y=ymin
1947 : integer ibcymax ! bc flag for y=ymax
1948 : real bcymax(:) ! (nx) ! bc data vs. x at y=ymax
1949 :
1950 : ! with interpretation:
1951 : ! ibcxmin -- indicator for boundary condition at x(1):
1952 : ! bcxmin(...) -- boundary condition data
1953 : ! =-1 -- periodic boundary condition
1954 : ! =0 -- use "not a knot"
1955 : ! =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
1956 : ! =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
1957 : ! =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
1958 : ! =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
1959 : ! =5 -- match 1st derivative to 1st divided difference
1960 : ! =6 -- match 2nd derivative to 2nd divided difference
1961 : ! =7 -- match 3rd derivative to 3rd divided difference
1962 : ! (for more detailed definition of BCs 5-7, see the
1963 : ! comments of subroutine mkspline)
1964 : ! NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
1965 :
1966 : ! ibcxmax -- indicator for boundary condition at x(nx):
1967 : ! bcxmax(...) -- boundary condition data
1968 : ! (interpretation as with ibcxmin, bcxmin)
1969 : ! NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
1970 :
1971 : ! and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
1972 : ! (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
1973 :
1974 : ! output linear grid flags and completion code (ier=0 is normal):
1975 :
1976 : integer ilinx ! =1: x grid is "nearly" equally spaced
1977 : integer iliny ! =1: y grid is "nearly" equally spaced
1978 : ! ilinx and iliny are set to zero if corresponding grids are not equally
1979 : ! spaced
1980 :
1981 : integer ier ! =0 on exit if there is no error.
1982 :
1983 : ! if there is an error, ier is set and a message is output on unit 6.
1984 : ! these are considered programming errors in the calling routine.
1985 :
1986 : ! possible errors:
1987 : ! x(...) not strict ascending
1988 : ! y(...) not strict ascending
1989 : ! nx < 4
1990 : ! ny < 4
1991 : ! invalid boundary condition flag
1992 :
1993 : !-----------------------
1994 : integer ierx,iery
1995 :
1996 1 : real, dimension(:,:), allocatable :: fwk
1997 1 : real :: zbcmin,zbcmax
1998 : integer ix,iy,ibcmin,ibcmax
1999 :
2000 1 : real, dimension(:,:,:), allocatable :: fcorr
2001 : integer iflg2
2002 3 : real zdiff(2),hy
2003 :
2004 1 : real, pointer :: f(:,:,:)
2005 1 : f(1:4,1:nf2,1:ny) => f1(1:4*nf2*ny)
2006 :
2007 : !-----------------------
2008 :
2009 : ! see if 2nd pass is needed due to inhomogeneous d/dy bdy cond.
2010 :
2011 1 : iflg2=0
2012 1 : if(ibcymin /= -1) then
2013 1 : if((ibcymin == 1).or.(ibcymin == 2)) then
2014 0 : do ix=1,nx
2015 0 : if (bcymin(ix) /= 0.0) iflg2=1
2016 : end do
2017 : end if
2018 1 : if((ibcymax == 1).or.(ibcymax == 2)) then
2019 0 : do ix=1,nx
2020 0 : if (bcymax(ix) /= 0.0) iflg2=1
2021 : end do
2022 : end if
2023 : end if
2024 :
2025 : ! check boundary condition specifications
2026 :
2027 1 : ier=0
2028 :
2029 1 : call ibc_ck(ibcxmin,'bcspline','xmin',-1,7,ier)
2030 2 : if(ibcxmin >= 0) call ibc_ck(ibcxmax,'bcspline','xmax',0,7,ier)
2031 1 : call ibc_ck(ibcymin,'bcspline','ymin',-1,7,ier)
2032 1 : if(ibcymin >= 0) call ibc_ck(ibcymax,'bcspline','ymax',0,7,ier)
2033 :
2034 : ! check ilinx & x vector
2035 :
2036 1 : call splinck(x,nx,ilinx,1.0e-3,ierx)
2037 1 : if(ierx /= 0) ier=2
2038 :
2039 1 : if(ier == 2) then
2040 0 : write(6,'('' ?bcspline: x axis not strict ascending'')')
2041 : end if
2042 :
2043 : ! check iliny & y vector
2044 :
2045 1 : call splinck(y,ny,iliny,1.0e-3,iery)
2046 1 : if(iery /= 0) ier=3
2047 :
2048 1 : if(ier == 3) then
2049 0 : write(6,'('' ?bcspline: y axis not strict ascending'')')
2050 : end if
2051 :
2052 1 : if(ier /= 0) return
2053 :
2054 : !------------------------------------
2055 1 : allocate(fwk(2,max(nx,ny)))
2056 :
2057 : ! evaluate fxx (spline in x direction)
2058 :
2059 1 : zbcmin=0
2060 1 : zbcmax=0
2061 36 : do iy=1,ny
2062 1120 : fwk(1,1:nx) = f(1,1:nx,iy)
2063 35 : if((ibcxmin == 1).or.(ibcxmin == 2)) zbcmin=bcxmin(iy)
2064 35 : if((ibcxmax == 1).or.(ibcxmax == 2)) zbcmax=bcxmax(iy)
2065 35 : call mkspline(x,nx,fwk,ibcxmin,zbcmin,ibcxmax,zbcmax,ilinx,ier)
2066 35 : if(ier /= 0) then
2067 0 : deallocate(fwk)
2068 0 : return
2069 : end if
2070 1121 : f(2,1:nx,iy)=fwk(2,1:nx)
2071 : end do
2072 :
2073 : ! evaluate fyy (spline in y direction)
2074 : ! use homogeneous boundary condition; correction done later if necessary
2075 :
2076 1 : zbcmin=0
2077 1 : zbcmax=0
2078 1 : ibcmin=ibcymin
2079 1 : ibcmax=ibcymax
2080 32 : do ix=1,nx
2081 1116 : fwk(1,1:ny) = f(1,ix,1:ny)
2082 31 : if(iflg2 == 1) then
2083 0 : if((ibcymin == 1).or.(ibcymin == 2)) ibcmin=0
2084 0 : if((ibcymax == 1).or.(ibcymax == 2)) ibcmax=0
2085 : end if
2086 31 : call mkspline(y,ny,fwk,ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
2087 31 : if(ier /= 0) then
2088 0 : deallocate(fwk)
2089 0 : return
2090 : end if
2091 1117 : f(3,ix,1:ny)=fwk(2,1:ny)
2092 : end do
2093 :
2094 : ! evaluate fxxyy (spline fxx in y direction; BC simplified; avg
2095 : ! d2(d2f/dx2)/dy2 and d2(df2/dy2)/dx2
2096 :
2097 : zbcmin=0
2098 : zbcmax=0
2099 1 : ibcmin=ibcymin
2100 1 : ibcmax=ibcymax
2101 32 : do ix=1,nx
2102 1116 : fwk(1,1:ny) = f(2,ix,1:ny)
2103 31 : if(iflg2 == 1) then
2104 0 : if((ibcymin == 1).or.(ibcymin == 2)) ibcmin=0
2105 0 : if((ibcymax == 1).or.(ibcymax == 2)) ibcmax=0
2106 : end if
2107 31 : call mkspline(y,ny,fwk,ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
2108 31 : if(ier /= 0) then
2109 0 : deallocate(fwk)
2110 0 : return
2111 : end if
2112 1117 : f(4,ix,1:ny)= fwk(2,1:ny)
2113 : end do
2114 :
2115 1 : if(iflg2 == 1) then
2116 0 : allocate(fcorr(2,nx,ny))
2117 :
2118 : ! correct for inhomogeneous y boundary condition
2119 :
2120 0 : do ix=1,nx
2121 : ! the desired inhomogeneous BC is the difference btw the
2122 : ! requested derivative (1st or 2nd) and the current value
2123 :
2124 0 : zdiff(1)=0.0
2125 0 : if(ibcymin == 1) then
2126 0 : hy=y(2)-y(1)
2127 : zdiff(1)=(f(1,ix,2)-f(1,ix,1))/hy +
2128 0 : & hy*(-2*f(3,ix,1)-f(3,ix,2))/6
2129 0 : zdiff(1)=bcymin(ix)-zdiff(1)
2130 0 : else if(ibcymin == 2) then
2131 0 : zdiff(1)=bcymin(ix)-f(3,ix,1)
2132 : end if
2133 :
2134 0 : zdiff(2)=0.0
2135 0 : if(ibcymax == 1) then
2136 0 : hy=y(ny)-y(ny-1)
2137 : zdiff(2)=(f(1,ix,ny)-f(1,ix,ny-1))/hy +
2138 0 : & hy*(2*f(3,ix,ny)+f(3,ix,ny-1))/6
2139 0 : zdiff(2)=bcymax(ix)-zdiff(2)
2140 0 : else if(ibcymax == 2) then
2141 0 : zdiff(2)=bcymax(ix)-f(3,ix,ny)
2142 : end if
2143 :
2144 0 : fwk(1,1:ny)=0.0 ! values are zero; only BC is not
2145 0 : call mkspline(y,ny,fwk,ibcymin,zdiff(1),ibcymax,zdiff(2),iliny,ier)
2146 0 : if(ier /= 0) then
2147 0 : deallocate(fwk,fcorr)
2148 0 : return
2149 : end if
2150 0 : fcorr(1,ix,1:ny)=fwk(2,1:ny) ! fyy-correction
2151 : end do
2152 :
2153 : zbcmin=0
2154 : zbcmax=0
2155 0 : do iy=1,ny
2156 0 : fwk(1,1:nx)=fcorr(1,1:nx,iy)
2157 0 : call mkspline(x,nx,fwk,ibcxmin,zbcmin,ibcxmax,zbcmax,ilinx,ier)
2158 0 : if(ier /= 0) then
2159 0 : deallocate(fwk,fcorr)
2160 0 : return
2161 : end if
2162 0 : fcorr(2,1:nx,iy)=fwk(2,1:nx) ! fxxyy-correction
2163 : end do
2164 :
2165 0 : f(3:4,1:nx,1:ny)=f(3:4,1:nx,1:ny)+fcorr(1:2,1:nx,1:ny)
2166 :
2167 0 : deallocate(fcorr)
2168 : end if
2169 :
2170 : ! correction spline -- f=fxx=zero; fyy & fxxyy are affected
2171 :
2172 1 : deallocate(fwk)
2173 : !------------------------------------
2174 :
2175 : ! that's all
2176 :
2177 1 : return
2178 2 : end subroutine do_mkbicub
2179 :
2180 97 : subroutine mkspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,ilinx,ier)
2181 :
2182 : ! make a 2-coefficient 1d spline
2183 :
2184 : ! only 2 coefficients, the data and its 2nd derivative, are needed to
2185 : ! fully specify a spline. See e.g. Numerical Recipes in Fortran-77
2186 : ! (2nd edition) chapter 3, section on cubic splines.
2187 :
2188 : ! input:
2189 : integer nx ! no. of data points
2190 : real x(:) ! (nx) ! x axis data, strict ascending order
2191 :
2192 : ! input/output:
2193 : real fspl(:,:) ! (2,nx) ! f(1,*): data in; f(2,*): coeffs out
2194 :
2195 : ! f(1,j) = f(x(j)) on input (unchanged on output)
2196 : ! f(2,j) = f''(x(j)) (of interpolating spline) (on output).
2197 :
2198 : ! ...boundary conditions...
2199 :
2200 : ! input:
2201 :
2202 : integer ibcxmin ! b.c. flag @ x=xmin=x(1)
2203 : real bcxmin ! b.c. data @xmin
2204 :
2205 : integer ibcxmax ! b.c. flag @ x=xmax=x(nx)
2206 : real bcxmax ! b.c. data @xmax
2207 :
2208 : ! ibcxmin=-1 -- periodic boundary condition
2209 : ! (bcxmin,ibcxmax,bcxmax are ignored)
2210 :
2211 : ! the output spline s satisfies
2212 : ! s'(x(1))=s'(x(nx)) ..and.. s''(x(1))=s''(x(nx))
2213 :
2214 : ! if non-periodic boundary conditions are used, then the xmin and xmax
2215 : ! boundary conditions can be specified independently:
2216 :
2217 : ! ibcxmin (ibcxmax) = 0 -- this specifies a "not a knot" boundary
2218 : ! condition, see "cubsplb.for". This is a common way
2219 : ! for inferring a "good" spline boundary condition
2220 : ! automatically from data in the vicinity of the
2221 : ! boundary. ... bcxmin (bcxmax) are ignored.
2222 :
2223 : ! ibcxmin (ibcxmax) = 1 -- boundary condition is to have s'(x(1))
2224 : ! ( s'(x(nx)) ) match the passed value bcxmin (bcxmax).
2225 :
2226 : ! ibcxmin (ibcxmax) = 2 -- boundary condition is to have s''(x(1))
2227 : ! ( s''(x(nx)) ) match the passed value bcxmin (bcxmax).
2228 :
2229 : ! ibcxmin (ibcxmax) = 3 -- boundary condition is to have s'(x(1))=0.0
2230 : ! ( s'(x(nx))=0.0 )
2231 :
2232 : ! ibcxmin (ibcxmax) = 4 -- boundary condition is to have s''(x(1))=0.0
2233 : ! ( s''(x(nx))=0.0 )
2234 :
2235 : ! ibcxmin (ibcxmax) = 5 -- boundary condition is to have s'(x(1))
2236 : ! ( s'(x(nx)) ) match the 1st divided difference
2237 : ! e.g. at x(1): d(1)/h(1), where
2238 : ! d(j)=f(1,j+1)-f(1,j)
2239 : ! h(j)=x(j+1)-x(j)
2240 :
2241 : ! ibcxmin (ibcxmax) = 6 -- BC is to have s''(x(1)) ( s''(x(nx)) )
2242 : ! match the 2nd divided difference
2243 : ! e.g. at x(1):
2244 : ! e(1) = [d(2)/h(2) - d(1)/h(1)]/(0.5*(h(1)+h(2)))
2245 :
2246 : ! ibcxmin (ibcxmax) = 7 -- BC is to have s'''(x(1)) ( s'''(x(nx)) )
2247 : ! match the 3rd divided difference
2248 : ! e.g. at x(1): [e(2)-e(1)]/(0.33333*(h(1)+h(2)+h(3)))
2249 :
2250 : ! output:
2251 :
2252 : integer ilinx ! =1: hint, x axis is ~evenly spaced
2253 :
2254 : ! let dx[avg] = (x(nx)-x(1))/(nx-1)
2255 : ! let dx[j] = x(j+1)-x(j), for all j satisfying 1 <= j < nx
2256 :
2257 : ! if for all such j, abs(dx[j]-dx[avg]) <= (1.0e-3*dx[avg]) then
2258 : ! ilinx=1 is returned, indicating the data is (at least nearly)
2259 : ! evenly spaced. Even spacing is useful, for speed of zone lookup,
2260 : ! when evaluating a spline.
2261 :
2262 : ! if the even spacing condition is not satisfied, ilinx=2 is returned.
2263 :
2264 : integer ier ! exit code, 0=OK
2265 :
2266 : ! an error code is returned if the x axis is not strict ascending,
2267 : ! or if nx < 4, or if an invalid boundary condition specification was
2268 : ! input.
2269 :
2270 : !------------------------------------
2271 :
2272 : ! this routine calls traditional 4 coefficient spline software, and
2273 : ! translates the result to 2 coefficient form.
2274 :
2275 : ! this could be done more efficiently but we decided out of conservatism
2276 : ! to use the traditional software.
2277 :
2278 : !------------------------------------
2279 : ! workspaces -- f90 dynamically allocated
2280 :
2281 97 : real, dimension(:,:), allocatable :: fspl4 ! traditional 4-spline
2282 97 : real, dimension(:), allocatable :: wk ! cspline workspace
2283 : integer i,inwk
2284 :
2285 : !------------------------------------
2286 :
2287 97 : allocate(fspl4(4,nx),wk(nx))
2288 :
2289 : ! make the traditional call
2290 :
2291 3352 : do i=1,nx
2292 3255 : fspl4(1,i)=fspl(1,i)
2293 3352 : fspl(2,i)=0.0 ! for now
2294 : end do
2295 :
2296 97 : inwk=nx
2297 :
2298 : ! boundary conditions imposed by cspline...
2299 :
2300 : call cspline(x,nx,fspl4,ibcxmin,bcxmin,ibcxmax,bcxmax,
2301 97 : & wk,inwk,ilinx,ier)
2302 :
2303 97 : if(ier == 0) then
2304 :
2305 : ! copy the output -- careful of end point.
2306 :
2307 3255 : do i=1,nx-1
2308 3255 : fspl(2,i) = 2.0*fspl4(3,i)
2309 : end do
2310 97 : fspl(2,nx) = 2.0*fspl4(3,nx-1) + (x(nx)-x(nx-1))*6.0*fspl4(4,nx-1)
2311 : end if
2312 :
2313 97 : deallocate(fspl4,wk)
2314 :
2315 97 : return
2316 : end subroutine mkspline
2317 :
2318 99 : subroutine splinck(x,inx,ilinx,ztol,ier)
2319 :
2320 : ! check if a grid is strictly ascending and if it is evenly spaced
2321 : ! to w/in ztol
2322 :
2323 : integer inx
2324 : real x(:) ! (inx) ! input -- grid to check
2325 :
2326 : integer ilinx ! output -- =1 if evenly spaced =2 O.W.
2327 :
2328 : real ztol ! input -- spacing check tolerance
2329 :
2330 : integer ier ! output -- =0 if OK
2331 :
2332 : ! ier=1 is returned if x(1...inx) is NOT STRICTLY ASCENDING...
2333 99 : real dxavg,zeps,zdiffx,zdiff
2334 : integer ix
2335 :
2336 : !-------------------------------
2337 :
2338 99 : ier=0
2339 99 : ilinx=1
2340 99 : if(inx <= 1) return
2341 :
2342 99 : dxavg=(x(inx)-x(1))/(inx-1)
2343 99 : zeps=abs(ztol*dxavg)
2344 :
2345 3321 : do ix=2,inx
2346 3222 : zdiffx=(x(ix)-x(ix-1))
2347 3222 : if(zdiffx <= 0.0) ier=2
2348 3222 : zdiff=zdiffx-dxavg
2349 3321 : if(abs(zdiff) > zeps) then
2350 0 : ilinx=2
2351 : end if
2352 : end do
2353 : 10 continue
2354 :
2355 : return
2356 : end subroutine splinck
2357 :
2358 97 : subroutine V_SPLINE(k_bc1,k_bcn,n,x,f,wk)
2359 : !***********************************************************************
2360 : !V_SPLINE evaluates the coefficients for a 1d cubic interpolating spline
2361 : !References:
2362 : ! Forsythe, Malcolm, Moler, Computer Methods for Mathematical
2363 : ! Computations, Prentice-Hall, 1977, p.76
2364 : ! Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
2365 : ! 1996, p.251
2366 : ! W.A.Houlberg, D.McCune 3/2000
2367 : !Input:
2368 : ! k_bc1-option for BC at x(1)
2369 : ! =-1 periodic, ignore k_bcn
2370 : ! =0 not-a-knot
2371 : ! =1 s'(x1) = input value of f(2,1)
2372 : ! =2 s''(x1) = input value of f(3,1)
2373 : ! =3 s'(x1) = 0.0
2374 : ! =4 s''(x1) = 0.0
2375 : ! =5 match first derivative to first 2 points
2376 : ! =6 match second derivative to first 3 points
2377 : ! =7 match third derivative to first 4 points
2378 : ! =else use not-a-knot
2379 : ! k_bcn-option for boundary condition at x(n)
2380 : ! =0 not-a-knot
2381 : ! =1 s'(x1) = input value of f(2,1)
2382 : ! =2 s''(x1) = input value of f(3,1)
2383 : ! =3 s'(x1) = 0.0
2384 : ! =4 s''(x1) = 0.0
2385 : ! =5 match first derivative to first 2 points
2386 : ! =6 match second derivative to first 3 points
2387 : ! =7 match third derivative to first 4 points
2388 : ! =else use knot-a-knot
2389 : ! n-number of data points or knots-(n >= 2)
2390 : ! x(n)-abscissas of the knots in strictly increasing order
2391 : ! f(1,i)-ordinates of the knots
2392 : ! f(2,1)-input value of s'(x1) for k_bc1=1
2393 : ! f(2,n)-input value of s'(xn) for k_bcn=1
2394 : ! f(3,1)-input value of s''(x1) for k_bc1=2
2395 : ! f(3,n)-input value of s''(xn) for k_bcn=2
2396 : ! wk(n)-scratch work area for periodic BC
2397 : !Output:
2398 : ! f(2,i)=s'(x(i))
2399 : ! f(3,i)=s''(x(i))
2400 : ! f(4,i)=s'''(x(i))
2401 : !Comments:
2402 : ! s(x)=f(1,i)+f(2,i)*(x-x(i))+f(3,i)*(x-x(i))**2/2!
2403 : ! +f(4,i)*(x-x(i))**3/3! for x(i) <= x <= x(i+1)
2404 : ! W_SPLINE can be used to evaluate the spline and its derivatives
2405 : ! The cubic spline is twice differentiable (C2)
2406 :
2407 : ! bugfixes -- dmc 24 Feb 2004:
2408 : ! (a) fixed logic for not-a-knot:
2409 : ! ! Set f(3,1) for not-a-knot
2410 : ! if (k_bc1 <= 0.or.k_bc1 > 7) then ...
2411 : ! instead of
2412 : ! ! Set f(3,1) for not-a-knot
2413 : ! if (k_bc1 <= 0.or.k_bc1 > 5) then ...
2414 : ! and similarly for logic after cmt
2415 : ! ! Set f(3,n) for not-a-knot
2416 : ! as required since k_bc*=6 and k_bc*=7 are NOT not-a-knot BCs.
2417 :
2418 : ! (b) the BCs to fix 2nd derivative at end points did not work if that
2419 : ! 2nd derivative were non-zero. The reason is that in those cases
2420 : ! the off-diagonal matrix elements nearest the corners are not
2421 : ! symmetric; i.e. elem(1,2) /= elem(2,1) and
2422 : ! elem(n-1,n) /= elem(n,n-1) where I use "elem" to refer to
2423 : ! the tridiagonal matrix elements. The correct values for the
2424 : ! elements is: elem(1,2)=0, elem(2,1)=x(2)-x(1)
2425 : ! elem(n,n-1)=0, elem(n-1,n)=x(n)-x(n-1)
2426 : ! the old code in effect had these as all zeroes. Since this
2427 : ! meant the wrong set of linear equations was solved, the
2428 : ! resulting spline had a discontinuity in its 1st derivative
2429 : ! at x(2) and x(n-1). Fixed by introducing elem21 and elemnn1
2430 : ! to represent the non-symmetric lower-diagonal values. Since
2431 : ! elem21 & elemnn1 are both on the lower diagonals, logic to
2432 : ! use them occurs in the non-periodic forward elimination loop
2433 : ! only. DMC 24 Feb 2004.
2434 : !***********************************************************************
2435 : implicit none
2436 : !Declaration of input variables
2437 : integer :: k_bc1, k_bcn, n
2438 : real :: x(*), wk(*), f(4,*)
2439 : !Declaration in local variables
2440 : integer :: i, ib, imax, imin
2441 97 : real :: a1, an, b1, bn, q, t, hn
2442 97 : real :: elem21, elemnn1 ! (dmc)
2443 :
2444 : !Set default range
2445 97 : imin=1
2446 97 : imax=n
2447 : !Set first and second BC values
2448 97 : a1=0.0
2449 97 : b1=0.0
2450 97 : an=0.0
2451 97 : bn=0.0
2452 97 : if (k_bc1 == 1) then
2453 0 : a1=f(2,1)
2454 97 : else if(k_bc1 == 2) then
2455 0 : b1=f(3,1)
2456 97 : else if(k_bc1 == 5) then
2457 0 : a1=(f(1,2)-f(1,1))/(x(2)-x(1))
2458 97 : else if(k_bc1 == 6) then
2459 : b1=2.0*((f(1,3)-f(1,2))/(x(3)-x(2))
2460 0 : & -(f(1,2)-f(1,1))/(x(2)-x(1)))/(x(3)-x(1))
2461 : end if
2462 97 : if (k_bcn == 1) then
2463 0 : an=f(2,n)
2464 97 : else if(k_bcn == 2) then
2465 0 : bn=f(3,n)
2466 97 : else if(k_bcn == 5) then
2467 0 : an=(f(1,n)-f(1,n-1))/(x(n)-x(n-1))
2468 97 : else if(k_bcn == 6) then
2469 : bn=2.0*((f(1,n)-f(1,n-1))/(x(n)-x(n-1))
2470 0 : & -(f(1,n-1)-f(1,n-2))/(x(n-1)-x(n-2)))/(x(n)-x(n-2))
2471 : end if
2472 : !Clear f(2:4,n)
2473 97 : f(2,n)=0.0
2474 97 : f(3,n)=0.0
2475 97 : f(4,n)=0.0
2476 97 : if (n == 2) then
2477 : !Coefficients for n=2
2478 0 : f(2,1)=(f(1,2)-f(1,1))/(x(2)-x(1))
2479 0 : f(3,1)=0.0
2480 0 : f(4,1)=0.0
2481 0 : f(2,2)=f(2,1)
2482 0 : f(3,2)=0.0
2483 0 : f(4,2)=0.0
2484 97 : else if(n > 2) then
2485 : !Set up tridiagonal system for A*y=B where y(i) are the second
2486 : ! derivatives at the knots
2487 : ! f(2,i) are the diagonal elements of A
2488 : ! f(4,i) are the off-diagonal elements of A
2489 : ! f(3,i) are the B elements/3, and will become c/3 upon solution
2490 97 : f(4,1)=x(2)-x(1)
2491 97 : f(3,2)=(f(1,2)-f(1,1))/f(4,1)
2492 3158 : do i=2,n-1
2493 3061 : f(4,i)=x(i+1)-x(i)
2494 3061 : f(2,i)=2.0*(f(4,i-1)+f(4,i))
2495 3061 : f(3,i+1)=(f(1,i+1)-f(1,i))/f(4,i)
2496 3158 : f(3,i)=f(3,i+1)-f(3,i)
2497 : end do
2498 :
2499 : ! (dmc): save now:
2500 :
2501 97 : elem21=f(4,1)
2502 97 : elemnn1=f(4,n-1)
2503 :
2504 : ! BC's
2505 : ! Left
2506 : if (k_bc1 == -1) then
2507 0 : f(2,1)=2.0*(f(4,1)+f(4,n-1))
2508 0 : f(3,1)=(f(1,2)-f(1,1))/f(4,1)-(f(1,n)-f(1,n-1))/f(4,n-1)
2509 0 : wk(1)=f(4,n-1)
2510 0 : do i=2,n-3
2511 0 : wk(i)=0.0
2512 : end do
2513 0 : wk(n-2)=f(4,n-2)
2514 0 : wk(n-1)=f(4,n-1)
2515 : else if(k_bc1 == 1.or.k_bc1 == 3.or.k_bc1 == 5) then
2516 0 : f(2,1)=2.0*f(4,1)
2517 0 : f(3,1)=(f(1,2)-f(1,1))/f(4,1)-a1
2518 : else if(k_bc1 == 2.or.k_bc1 == 4.or.k_bc1 == 6) then
2519 0 : f(2,1)=2.0*f(4,1)
2520 0 : f(3,1)=f(4,1)*b1/3.0
2521 0 : f(4,1)=0.0 ! upper diagonal only (dmc: cf elem21)
2522 : else if(k_bc1 == 7) then
2523 0 : f(2,1)=-f(4,1)
2524 0 : f(3,1)=f(3,3)/(x(4)-x(2))-f(3,2)/(x(3)-x(1))
2525 0 : f(3,1)=f(3,1)*f(4,1)**2/(x(4)-x(1))
2526 : else ! not a knot:
2527 97 : imin=2
2528 97 : f(2,2)=f(4,1)+2.0*f(4,2)
2529 97 : f(3,2)=f(3,2)*f(4,2)/(f(4,1)+f(4,2))
2530 : end if
2531 : ! Right
2532 97 : if (k_bcn == 1.or.k_bcn == 3.or.k_bcn == 5) then
2533 0 : f(2,n)=2.0*f(4,n-1)
2534 0 : f(3,n)=-(f(1,n)-f(1,n-1))/f(4,n-1)+an
2535 : else if(k_bcn == 2.or.k_bcn == 4.or.k_bcn == 6) then
2536 0 : f(2,n)=2.0*f(4,n-1)
2537 0 : f(3,n)=f(4,n-1)*bn/3.0
2538 : !xxx f(4,n-1)=0.0 ! dmc: preserve f(4,n-1) for back subst.
2539 0 : elemnn1=0.0 ! lower diagonal only (dmc)
2540 : else if(k_bcn == 7) then
2541 0 : f(2,n)=-f(4,n-1)
2542 0 : f(3,n)=f(3,n-1)/(x(n)-x(n-2))-f(3,n-2)/(x(n-1)-x(n-3))
2543 0 : f(3,n)=-f(3,n)*f(4,n-1)**2/(x(n)-x(n-3))
2544 97 : else if(k_bc1 /= -1) then ! not a knot:
2545 97 : imax=n-1
2546 97 : f(2,n-1)=2.0*f(4,n-2)+f(4,n-1)
2547 97 : f(3,n-1)=f(3,n-1)*f(4,n-2)/(f(4,n-1)+f(4,n-2))
2548 : end if
2549 : ! Limit solution for only three points in domain
2550 97 : if (n == 3) then
2551 0 : f(3,1)=0.0
2552 0 : f(3,n)=0.0
2553 : end if
2554 97 : if (k_bc1 == -1) then
2555 : !Solve system of equations for second derivatives at the knots
2556 : ! Periodic BC
2557 : ! Forward elimination
2558 0 : do i=2,n-2
2559 0 : t=f(4,i-1)/f(2,i-1)
2560 0 : f(2,i)=f(2,i)-t*f(4,i-1)
2561 0 : f(3,i)=f(3,i)-t*f(3,i-1)
2562 0 : wk(i)=wk(i)-t*wk(i-1)
2563 0 : q=wk(n-1)/f(2,i-1)
2564 0 : wk(n-1)=-q*f(4,i-1)
2565 0 : f(2,n-1)=f(2,n-1)-q*wk(i-1)
2566 0 : f(3,n-1)=f(3,n-1)-q*f(3,i-1)
2567 : end do
2568 : ! Correct the n-1 element
2569 0 : wk(n-1)=wk(n-1)+f(4,n-2)
2570 : ! Complete the forward elimination
2571 : ! wk(n-1) and wk(n-2) are the off-diag elements of the lower corner
2572 0 : t=wk(n-1)/f(2,n-2)
2573 0 : f(2,n-1)=f(2,n-1)-t*wk(n-2)
2574 0 : f(3,n-1)=f(3,n-1)-t*f(3,n-2)
2575 : ! Back substitution
2576 0 : f(3,n-1)=f(3,n-1)/f(2,n-1)
2577 0 : f(3,n-2)=(f(3,n-2)-wk(n-2)*f(3,n-1))/f(2,n-2)
2578 0 : do ib=3,n-1
2579 0 : i=n-ib
2580 0 : f(3,i)=(f(3,i)-f(4,i)*f(3,i+1)-wk(i)*f(3,n-1))/f(2,i)
2581 : end do
2582 0 : f(3,n)=f(3,1)
2583 : else
2584 : ! Non-periodic BC
2585 : ! Forward elimination
2586 : ! For Not-A-Knot BC the off-diagonal end elements are not equal
2587 3061 : do i=imin+1,imax
2588 2964 : if ((i == n-1).and.(imax == n-1)) then
2589 97 : t=(f(4,i-1)-f(4,i))/f(2,i-1)
2590 : else
2591 2867 : if(i == 2) then
2592 0 : t=elem21/f(2,i-1)
2593 2867 : else if(i == n) then
2594 0 : t=elemnn1/f(2,i-1)
2595 : else
2596 2867 : t=f(4,i-1)/f(2,i-1)
2597 : end if
2598 : end if
2599 2964 : if ((i == imin+1).and.(imin == 2)) then
2600 97 : f(2,i)=f(2,i)-t*(f(4,i-1)-f(4,i-2))
2601 : else
2602 2867 : f(2,i)=f(2,i)-t*f(4,i-1)
2603 : end if
2604 3061 : f(3,i)=f(3,i)-t*f(3,i-1)
2605 : end do
2606 : ! Back substitution
2607 97 : f(3,imax)=f(3,imax)/f(2,imax)
2608 3061 : do ib=1,imax-imin
2609 2964 : i=imax-ib
2610 3061 : if (abs(f(3,i)) > 1e20) then
2611 0 : f(3,i) = 0.0 ! protect against overflow
2612 2964 : else if((i == 2).and.(imin == 2)) then
2613 97 : f(3,i)=(f(3,i)-(f(4,i)-f(4,i-1))*f(3,i+1))/f(2,i)
2614 : else
2615 2867 : f(3,i)=(f(3,i)-f(4,i)*f(3,i+1))/f(2,i)
2616 : end if
2617 : end do
2618 : ! Reset d array to step size
2619 97 : f(4,1)=x(2)-x(1)
2620 97 : f(4,n-1)=x(n)-x(n-1)
2621 : ! Set f(3,1) for not-a-knot
2622 97 : if (k_bc1 <= 0.or.k_bc1 > 7) then
2623 97 : f(3,1)=(f(3,2)*(f(4,1)+f(4,2))-f(3,3)*f(4,1))/f(4,2)
2624 : end if
2625 : ! Set f(3,n) for not-a-knot
2626 97 : if (k_bcn <= 0.or.k_bcn > 7) then
2627 97 : f(3,n)=f(3,n-1)+(f(3,n-1)-f(3,n-2))*f(4,n-1)/f(4,n-2)
2628 : end if
2629 : end if
2630 : !f(3,i) is now the sigma(i) of the text and f(4,i) is the step size
2631 : !Compute polynomial coefficients
2632 3255 : do i=1,n-1
2633 : f(2,i)=
2634 3158 : & (f(1,i+1)-f(1,i))/f(4,i)-f(4,i)*(f(3,i+1)+2.0*f(3,i))
2635 3158 : f(4,i)=(f(3,i+1)-f(3,i))/f(4,i)
2636 3158 : f(3,i)=6.0*f(3,i)
2637 3255 : f(4,i)=6.0*f(4,i)
2638 : end do
2639 97 : if (k_bc1 == -1) then
2640 0 : f(2,n)=f(2,1)
2641 0 : f(3,n)=f(3,1)
2642 0 : f(4,n)=f(4,1)
2643 : else
2644 97 : hn=x(n)-x(n-1)
2645 97 : f(2,n)=f(2,n-1)+hn*(f(3,n-1)+0.5*hn*f(4,n-1))
2646 97 : f(3,n)=f(3,n-1)+hn*f(4,n-1)
2647 97 : f(4,n)=f(4,n-1)
2648 : if (k_bcn == 1.or.k_bcn == 3.or.k_bcn == 5) then
2649 0 : f(2,n)=an
2650 : else if(k_bcn == 2.or.k_bcn == 4.or.k_bcn == 6) then
2651 0 : f(3,n)=bn
2652 : end if
2653 : end if
2654 : end if
2655 97 : return
2656 : end subroutine V_SPLINE
2657 :
2658 :
2659 0 : subroutine zonfind(x,nx,zxget,i)
2660 :
2661 : integer nx
2662 : real zxget
2663 : real x(:) ! (nx)
2664 : integer i,nxm,i1,i2,ij,ii
2665 0 : real dx
2666 :
2667 : ! find index i such that x(i) <= zxget <= x(i+1)
2668 :
2669 : ! x(1...nx) is strict increasing and x(1) <= zxget <= x(nx)
2670 : ! (this is assumed to already have been checked -- no check here!)
2671 :
2672 0 : nxm=nx-1
2673 0 : if((i < 1).or.(i > nxm)) then
2674 : i1=1
2675 : i2=nx-1
2676 : GOTO 10
2677 : end if
2678 :
2679 0 : if(x(i) > zxget) then
2680 : ! look down
2681 0 : dx=x(i+1)-x(i)
2682 0 : if((x(i)-zxget) > 4*dx) then
2683 0 : i1=1
2684 0 : i2=i-1
2685 0 : GOTO 10
2686 : else
2687 0 : i2=i-1
2688 0 : do ij=i2,1,-1
2689 0 : if((x(ij) <= zxget).and.(zxget <= x(ij+1))) then
2690 0 : i=ij
2691 0 : return
2692 : end if
2693 : end do
2694 0 : i=1
2695 0 : return
2696 : end if
2697 0 : else if(x(i+1) < zxget) then
2698 : ! look up
2699 0 : dx=x(i+1)-x(i)
2700 0 : if((zxget-x(i+1)) > 4*dx) then
2701 : i1=i+1
2702 : i2=nxm
2703 : GOTO 10
2704 : else
2705 0 : i2=i+1
2706 0 : do ij=i2,nxm
2707 0 : if((x(ij) <= zxget).and.(zxget <= x(ij+1))) then
2708 0 : i=ij
2709 0 : return
2710 : end if
2711 : end do
2712 0 : ij=nxm
2713 : return
2714 : end if
2715 : else
2716 : ! already there...
2717 : return
2718 : end if
2719 :
2720 : !---------------------------
2721 : ! binary search
2722 :
2723 : 10 continue
2724 :
2725 0 : if(i1 == i2) then
2726 : ! found by proc. of elimination
2727 0 : i=i1
2728 0 : return
2729 : end if
2730 :
2731 0 : ii=(i1+i2)/2
2732 :
2733 0 : if(zxget < x(ii)) then
2734 0 : i2=ii-1
2735 0 : else if(zxget > x(ii+1)) then
2736 : i1=ii+1
2737 : else
2738 : ! found
2739 0 : i=ii
2740 0 : return
2741 : end if
2742 :
2743 : GOTO 10
2744 :
2745 : end subroutine zonfind
2746 :
2747 : end module bicub_sg
|
