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_db
21 :
22 : use const_def, only: dp
23 :
24 : !implicit none
25 :
26 :
27 : contains
28 :
29 : ! from PSPLINE by Doug McCune (version as of February, 2004)
30 :
31 :
32 : ! PSPLINE Home Page:
33 : ! http://w3.pppl.gov/NTCC/PSPLINE
34 :
35 : ! Doug McCune, Princeton University
36 : ! dmccune@pppl.gov
37 :
38 :
39 : ! bcspeval -- eval bicubic spline function and/or derivatives
40 :
41 0 : subroutine bcspeval_db(xget,yget,iselect,fval,x,nx,y,ny,ilinx,iliny,f,inf3,ier)
42 :
43 : implicit none
44 : integer iselect(6)
45 : integer ilinx,iliny,nx,ny,inf3,ier
46 :
47 : real(dp) xget,yget
48 : real(dp) fval(6)
49 : real(dp) x(nx),y(ny),f(4,4,inf3,ny)
50 :
51 : ! modification -- dmc 11 Jan 1999 -- remove SAVE stmts;
52 : ! break routine into these parts:
53 :
54 : ! bcspevxy -- find grid cell of target pt.
55 : ! bcspevfn -- evaluate function using output of bcpsevxy
56 :
57 : ! in cases where multiple functions are defined on the same grid,
58 : ! time can be saved by using bcspevxy once and then bcspevfn
59 : ! multiple times.
60 :
61 : ! input:
62 : ! (xget,yget) location where interpolated value is desired
63 : ! x(1) <= xget <= x(nx) expected
64 : ! y(1) <= yget <= y(ny) expected
65 :
66 : ! iselect select desired output
67 :
68 : ! iselect(1)=1 -- want function value (f) itself
69 : ! iselect(2)=1 -- want df/dx
70 : ! iselect(3)=1 -- want df/dy
71 : ! iselect(4)=1 -- want d2f/dx2
72 : ! iselect(5)=1 -- want d2f/dy2
73 : ! iselect(6)=1 -- want d2f/dxdy
74 :
75 : ! example: iselect(1)=iselect(2)=iselect(3)=1
76 : ! f, df/dx, and df/dy all evaluated
77 : ! iselect(4)=iselect(5)=iselect(6)=0
78 : ! 2nd derivatives not evaluated.
79 :
80 : ! the number of non zero values iselect(1:6)
81 : ! determines the number of outputs...
82 : ! see fval (output) description.
83 :
84 : ! new dmc December 2005 -- access to higher derivatives (even if not
85 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
86 : ! if iselect(1)=3 -- want 3rd derivatives
87 : ! iselect(2)=1 for d3f/dx3
88 : ! iselect(3)=1 for d3f/dx2dy
89 : ! iselect(4)=1 for d3f/dxdy2
90 : ! iselect(5)=1 for d3f/dy3
91 : ! number of non-zero values iselect(2:5) gives no. of outputs
92 : ! if iselect(1)=4 -- want 4th derivatives
93 : ! iselect(2)=1 for d4f/dx3dy
94 : ! iselect(3)=1 for d4f/dx2dy2
95 : ! iselect(4)=1 for d4f/dxdy3
96 : ! number of non-zero values iselect(2:4) gives no. of outputs
97 : ! if iselect(1)=5 -- want 5th derivatives
98 : ! iselect(2)=1 for d5f/dx3dy2
99 : ! iselect(3)=1 for d5f/dx2dy3
100 : ! number of non-zero values iselect(2:3) gives no. of outputs
101 : ! if iselect(1)=6 -- want 6th derivatives
102 : ! d6f/dx3dy3 -- one value is returned.
103 :
104 : ! x(1...nx) independent coordinate x, strict ascending
105 : ! y(1...ny) independent coordinate y, strict ascending
106 :
107 : ! ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
108 : ! iliny -- =1: flag that y is linearly spaced (avoid search for speed)
109 :
110 : ! **CAUTION** actual even spacing of x, y is NOT CHECKED HERE!
111 :
112 : !
113 : ! f the function values (at grid points) and spline coefs
114 :
115 : ! evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
116 : ! and y btw y(j) and y(j+1), dy=y-y(j),
117 :
118 : ! spline value =
119 : ! 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)
120 : ! +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))
121 : ! +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))
122 : ! +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))
123 :
124 : ! where d2=dy**2 and d3=dy**3.
125 :
126 : ! output:
127 : ! up to 6 elements of fval, ordered as follows:
128 : ! fval(1)=function value or lowest order derivative requested
129 : ! fval(2)=next order derivative
130 : ! etc
131 : ! the ordering is a subset of the sequence given under the "iselect"
132 : ! description.
133 :
134 : ! ier = 0 -- successful completion; = 1 -- an error occurred.
135 :
136 : !-------------------------------------------------------------------
137 : ! local
138 :
139 : integer :: i=0
140 : integer :: j=0
141 :
142 0 : real(dp) :: dx, dy
143 :
144 : !--------------------------
145 :
146 : call bcspevxy_db(xget,yget,x,nx,y,ny,ilinx,iliny,
147 0 : & i,j,dx,dy,ier)
148 0 : if(ier /= 0) return
149 :
150 : call bcspevfn_db(iselect,1,1,fval,(/i/),(/j/),
151 0 : & (/dx/),(/dy/),f,inf3,ny)
152 :
153 0 : return
154 : end subroutine bcspeval_db
155 :
156 :
157 : !-------------------------------------------------------------------------
158 : ! bcspevxy -- look up x-y zone
159 :
160 : ! this is the "first part" of bcspeval, see comments, above.
161 :
162 0 : subroutine bcspevxy_db(xget,yget,x,nx,y,ny,ilinx,iliny,
163 : & i,j,dx,dy,ier)
164 :
165 : integer nx,ny ! array dimensions
166 :
167 : real(dp) xget,yget ! target point
168 : real(dp) x(nx),y(ny) ! indep. coords.
169 :
170 : integer ilinx ! =1: assume x evenly spaced
171 : integer iliny ! =1: assume y evenly spaced
172 :
173 : ! output of bcspevxy
174 :
175 : integer i,j ! index to cell containing target pt
176 : real(dp) dx,dy ! displacement of target pt w/in cell
177 : ! dx=x-x(i) dy=y-y(j)
178 :
179 : integer ier ! return ier /= 0 on error
180 :
181 : !------------------------------------
182 :
183 : real(dp) zxget, zyget
184 0 : ier=0
185 :
186 : ! range check
187 :
188 0 : zxget=xget
189 0 : zyget=yget
190 :
191 0 : if((xget < x(1)).or.(xget > x(nx))) then
192 0 : zxtol=4.0d-7*max(abs(x(1)),abs(x(nx)))
193 0 : if((xget < x(1)-zxtol).or.(xget > x(nx)+zxtol)) then
194 0 : ier=1
195 : ! write(6,1001) xget,x(1),x(nx)
196 : ! 1001 format(' ?bcspeval: xget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
197 : else
198 : ! if((xget < x(1)-0.5*zxtol).or. (xget > x(nx)+0.5*zxtol)) write(6,1011) xget,x(1),x(nx)
199 : ! 1011 format(' %bcspeval: xget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
200 0 : if(xget < x(1)) then
201 0 : zxget=x(1)
202 : else
203 0 : zxget=x(nx)
204 : end if
205 : end if
206 : end if
207 0 : if((yget < y(1)).or.(yget > y(ny))) then
208 0 : zytol=4.0d-7*max(abs(y(1)),abs(y(ny)))
209 0 : if((yget < y(1)-zytol).or.(yget > y(ny)+zytol)) then
210 0 : ier=1
211 : ! write(6,1002) yget,y(1),y(ny)
212 : ! 1002 format(' ?bcspeval: yget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
213 : else
214 : ! if((yget < y(1)-0.5*zytol).or.(yget > y(ny)+0.5*zytol)) write(6,1012) yget,y(1),y(ny)
215 : ! 1012 format(' %bcspeval: yget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
216 0 : if(yget < y(1)) then
217 0 : zyget=y(1)
218 : else
219 0 : zyget=y(ny)
220 : end if
221 : end if
222 : end if
223 0 : if(ier /= 0) return
224 :
225 : ! now find interval in which target point lies..
226 :
227 0 : nxm=nx-1
228 0 : nym=ny-1
229 :
230 0 : if(ilinx == 1) then
231 0 : ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
232 0 : i=min(nxm, ii)
233 0 : if(zxget < x(i)) then
234 0 : i=i-1
235 0 : else if(zxget > x(i+1)) then
236 0 : i=i+1
237 : end if
238 : else
239 0 : if((1 <= i).and.(i < nxm)) then
240 0 : if((x(i) <= zxget).and.(zxget <= x(i+1))) then
241 : continue ! already have the zone
242 : else
243 0 : call zonfind_db(x,nx,zxget,i)
244 : end if
245 : else
246 0 : call zonfind_db(x,nx,zxget,i)
247 : end if
248 : end if
249 :
250 0 : if(iliny == 1) then
251 0 : jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
252 0 : j=min(nym, jj)
253 0 : if(zyget < y(j)) then
254 0 : j=j-1
255 0 : else if(zyget > y(j+1)) then
256 0 : j=j+1
257 : end if
258 : else
259 0 : if((1 <= j).and.(j < nym)) then
260 0 : if((y(j) <= zyget).and.(zyget <= y(j+1))) then
261 : continue ! already have the zone
262 : else
263 0 : call zonfind_db(y,ny,zyget,j)
264 : end if
265 : else
266 0 : call zonfind_db(y,ny,zyget,j)
267 : end if
268 : end if
269 :
270 0 : dx=zxget-x(i)
271 0 : dy=zyget-y(j)
272 :
273 0 : return
274 : end subroutine bcspevxy_db
275 :
276 : !------------------------------------------------------------------------
277 : ! bcspevfn -- OK now evaluate the bicubic spline
278 :
279 0 : subroutine bcspevfn_db(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(dp) 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(dp) dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors
331 :
332 : integer inf3 ! 3rd dimension of f -- >= nx
333 : real(dp) f(4,4,inf3,ny) ! bicubic fcn spline coeffs array
334 :
335 : ! usage example:
336 :
337 : ! 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs,
338 : ! find the x and y zone indices and displacements with respect
339 : ! to the "lower left corner" of the zone; store these in vectors
340 : ! iv,jv and dxv,dyv.
341 :
342 : ! 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only
343 : ! the 1st derivatives.
344 :
345 : ! 3. ivec is the length of the vector; ivd is the 1st dimension
346 : ! of the array fval to receive the output
347 :
348 : ! real(dp) fval(ivd,6)
349 : ! real(dp) xv(ivd),yv(ivd)
350 : ! integer iv(ivd),jv(ivd)
351 : ! real(dp) dxv(ivd),dyv(ivd)
352 : ! integer ict(6)
353 :
354 : ! real(dp) fspline(4,4,nx,ny) ! spline coeffs
355 : ! data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives
356 : ! ! only ... these will be output to
357 : ! ! fval(*,1) fval(*,2)
358 : ! ...
359 : ! do iv=1,ivec
360 : ! ... ! find indices and displacements
361 : ! end do
362 : ! call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny)
363 :
364 : !-------------------------------------------------------------------
365 : ! local:
366 :
367 : integer v ! vector element index
368 :
369 : ! OK can now do evaluations
370 :
371 0 : iaval=0 ! fval addressing
372 :
373 0 : if(ict(1) <= 2) then
374 0 : if((ict(1) > 0).or.(ict(1) == -1)) then
375 : ! evaluate f
376 0 : iaval=iaval+1
377 0 : do v=1,ivec
378 0 : i=iv(v)
379 0 : j=jv(v)
380 0 : dx=dxv(v)
381 0 : dy=dyv(v)
382 : fval(v,iaval)=
383 : & f(1,1,i,j)+dy*(f(1,2,i,j)+dy*(f(1,3,i,j)+dy*f(1,4,i,j)))
384 : & +dx*(f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
385 : & +dx*(f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
386 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))))))
387 : end do
388 : end if
389 :
390 0 : if((ict(2) > 0).and.(ict(1) /= -1)) then
391 : ! evaluate df/dx
392 0 : iaval=iaval+1
393 0 : do v=1,ivec
394 0 : i=iv(v)
395 0 : j=jv(v)
396 0 : dx=dxv(v)
397 0 : dy=dyv(v)
398 : fval(v,iaval)=
399 : & f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
400 : & +2.d0*dx*(
401 : & f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
402 : & +1.5d0*dx*(
403 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j)))
404 0 : & ))
405 : end do
406 : end if
407 :
408 0 : if((ict(3) > 0).and.(ict(1) /= -1)) then
409 : ! evaluate df/dy
410 0 : iaval=iaval+1
411 0 : do v=1,ivec
412 0 : i=iv(v)
413 0 : j=jv(v)
414 0 : dx=dxv(v)
415 0 : dy=dyv(v)
416 : fval(v,iaval)=
417 : & f(1,2,i,j)+dy*(2.d0*f(1,3,i,j)+dy*3.d0*f(1,4,i,j))
418 : & +dx*(f(2,2,i,j)+dy*(2.d0*f(2,3,i,j)+dy*3.d0*f(2,4,i,j))
419 : & +dx*(f(3,2,i,j)+dy*(2.d0*f(3,3,i,j)+dy*3.d0*f(3,4,i,j))
420 : & +dx*(f(4,2,i,j)+dy*(2.d0*f(4,3,i,j)+dy*3.d0*f(4,4,i,j))
421 0 : & )))
422 : end do
423 : end if
424 :
425 0 : if((ict(4) > 0).or.(ict(1) == -1)) then
426 : ! evaluate d2f/dx2
427 0 : iaval=iaval+1
428 0 : do v=1,ivec
429 0 : i=iv(v)
430 0 : j=jv(v)
431 0 : dx=dxv(v)
432 0 : dy=dyv(v)
433 : fval(v,iaval)=
434 : & 2.d0*(
435 : & f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j))))
436 : & +6.d0*dx*(
437 0 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
438 : end do
439 : end if
440 :
441 0 : if((ict(5) > 0).or.(ict(1) == -1)) then
442 : ! evaluate d2f/dy2
443 0 : iaval=iaval+1
444 0 : do v=1,ivec
445 0 : i=iv(v)
446 0 : j=jv(v)
447 0 : dx=dxv(v)
448 0 : dy=dyv(v)
449 : fval(v,iaval)=
450 : & 2.d0*f(1,3,i,j)+6.d0*dy*f(1,4,i,j)
451 : & +dx*(2.d0*f(2,3,i,j)+6.d0*dy*f(2,4,i,j)
452 : & +dx*(2.d0*f(3,3,i,j)+6.d0*dy*f(3,4,i,j)
453 0 : & +dx*(2.d0*f(4,3,i,j)+6.d0*dy*f(4,4,i,j))))
454 : end do
455 : end if
456 :
457 0 : if((ict(6) > 0).and.(ict(1) /= -1)) then
458 : ! evaluate d2f/dxdy
459 0 : iaval=iaval+1
460 0 : do v=1,ivec
461 0 : i=iv(v)
462 0 : j=jv(v)
463 0 : dx=dxv(v)
464 0 : dy=dyv(v)
465 : fval(v,iaval)=
466 : & f(2,2,i,j)+dy*(2.d0*f(2,3,i,j)+dy*3.d0*f(2,4,i,j))
467 : & +2.d0*dx*(f(3,2,i,j)+dy*(2.d0*f(3,3,i,j)+dy*3.d0*f(3,4,i,j))
468 : &+1.5d0*dx*(f(4,2,i,j)+dy*(2.d0*f(4,3,i,j)+dy*3.d0*f(4,4,i,j))
469 0 : & ))
470 : end do
471 : end if
472 :
473 0 : if(ict(1) == -1) then
474 0 : iaval=iaval+1
475 0 : do v=1,ivec
476 0 : i=iv(v)
477 0 : j=jv(v)
478 0 : dx=dxv(v)
479 0 : dy=dyv(v)
480 : fval(v,iaval)=
481 : & 4.d0*f(3,3,i,j)+12.d0*dy*f(3,4,i,j)
482 0 : & +dx*(12.d0*f(4,3,i,j)+36.d0*dy*f(4,4,i,j))
483 : end do
484 : end if
485 :
486 : !-----------------------------------
487 : ! access to 3rd derivatives
488 :
489 0 : else if(ict(1) == 3) then
490 0 : if(ict(2) == 1) then
491 : ! evaluate d3f/dx3 (not continuous)
492 0 : iaval=iaval+1
493 0 : do v=1,ivec
494 0 : i=iv(v)
495 0 : j=jv(v)
496 0 : dy=dyv(v)
497 : fval(v,iaval)=
498 : & +6.d0*(
499 0 : & f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
500 : end do
501 : end if
502 :
503 0 : if(ict(3) == 1) then
504 : ! evaluate d3f/dx2dy
505 0 : iaval=iaval+1
506 0 : do v=1,ivec
507 0 : i=iv(v)
508 0 : j=jv(v)
509 0 : dx=dxv(v)
510 0 : dy=dyv(v)
511 : fval(v,iaval)=
512 : & 2.d0*(
513 : & f(3,2,i,j)+dy*(2.d0*f(3,3,i,j)+dy*3.d0*f(3,4,i,j)))
514 : & +6.d0*dx*(
515 0 : & f(4,2,i,j)+dy*(2.d0*f(4,3,i,j)+dy*3.d0*f(4,4,i,j)))
516 : end do
517 : end if
518 :
519 0 : if(ict(4) == 1) then
520 : ! evaluate d3f/dxdy2
521 0 : iaval=iaval+1
522 0 : do v=1,ivec
523 0 : i=iv(v)
524 0 : j=jv(v)
525 0 : dx=dxv(v)
526 0 : dy=dyv(v)
527 : fval(v,iaval)=
528 : & (2.d0*f(2,3,i,j)+6.d0*dy*f(2,4,i,j)
529 : & +2.d0*dx*(2.d0*f(3,3,i,j)+6.d0*dy*f(3,4,i,j)
530 : & +1.5d0*dx*(2.d0*f(4,3,i,j)+6.d0*dy*f(4,4,i,j))
531 0 : & ))
532 : end do
533 : end if
534 :
535 0 : if(ict(5) == 1) then
536 : ! evaluate d3f/dy3 (not continuous)
537 0 : iaval=iaval+1
538 0 : do v=1,ivec
539 0 : i=iv(v)
540 0 : j=jv(v)
541 0 : dx=dxv(v)
542 : fval(v,iaval)=6.d0*(f(1,4,i,j)+
543 0 : & dx*(f(2,4,i,j)+dx*(f(3,4,i,j)+dx*f(4,4,i,j))))
544 : end do
545 : end if
546 :
547 : !-----------------------------------
548 : ! access to 4th derivatives
549 :
550 0 : else if(ict(1) == 4) then
551 0 : if(ict(2) == 1) then
552 : ! evaluate d4f/dx3dy (not continuous)
553 0 : iaval=iaval+1
554 0 : do v=1,ivec
555 0 : i=iv(v)
556 0 : j=jv(v)
557 0 : dy=dyv(v)
558 : fval(v,iaval)=
559 : & +6.d0*(
560 0 : & f(4,2,i,j)+dy*2.d0*(f(4,3,i,j)+dy*1.5d0*f(4,4,i,j)))
561 : end do
562 : end if
563 :
564 0 : if(ict(3) == 1) then
565 : ! evaluate d4f/dx2dy2
566 0 : iaval=iaval+1
567 0 : do v=1,ivec
568 0 : i=iv(v)
569 0 : j=jv(v)
570 0 : dx=dxv(v)
571 0 : dy=dyv(v)
572 : fval(v,iaval)=
573 : & 4.d0*f(3,3,i,j)+12.d0*dy*f(3,4,i,j)
574 0 : & +dx*(12.d0*f(4,3,i,j)+36.d0*dy*f(4,4,i,j))
575 : end do
576 : end if
577 :
578 0 : if(ict(4) == 1) then
579 : ! evaluate d4f/dxdy3 (not continuous)
580 0 : iaval=iaval+1
581 0 : do v=1,ivec
582 0 : i=iv(v)
583 0 : j=jv(v)
584 0 : dx=dxv(v)
585 : fval(v,iaval)=
586 : & 6.d0*(f(2,4,i,j)
587 0 : & +2.d0*dx*(f(3,4,i,j)+1.5d0*dx*f(4,4,i,j)))
588 : end do
589 : end if
590 :
591 : !-----------------------------------
592 : ! access to 5th derivatives
593 :
594 0 : else if(ict(1) == 5) then
595 0 : if(ict(2) == 1) then
596 : ! evaluate d5f/dx3dy2 (not continuous)
597 0 : iaval=iaval+1
598 0 : do v=1,ivec
599 0 : i=iv(v)
600 0 : j=jv(v)
601 0 : dy=dyv(v)
602 : fval(v,iaval)=
603 0 : & +12.d0*(f(4,3,i,j)+dy*3.d0*f(4,4,i,j))
604 : end do
605 : end if
606 :
607 0 : if(ict(3) == 1) then
608 : ! evaluate d5f/dx3dy2 (not continuous)
609 0 : iaval=iaval+1
610 0 : do v=1,ivec
611 0 : i=iv(v)
612 0 : j=jv(v)
613 0 : dx=dxv(v)
614 : fval(v,iaval)=
615 0 : & 12.d0*(f(3,4,i,j)+dx*3.d0*f(4,4,i,j))
616 : end do
617 : end if
618 :
619 : !-----------------------------------
620 : ! access to 6th derivatives
621 :
622 0 : else if(ict(1) == 6) then
623 : ! evaluate d6f/dx3dy3 (not continuous)
624 0 : iaval=iaval+1
625 0 : do v=1,ivec
626 0 : i=iv(v)
627 0 : j=jv(v)
628 : fval(v,iaval)=
629 0 : & 36.d0*f(4,4,i,j)
630 : end do
631 : end if
632 :
633 0 : return
634 : end subroutine bcspevfn_db
635 :
636 : !----------------------
637 :
638 : ! bcspline -- dmc 30 May 1996
639 :
640 : ! set up coefficients for bicubic spline with following BC's:
641 : ! FULL BC CONTROL at all bdys
642 :
643 : ! inhomogeneous explicit BCs -- this means setting of 1st or 2nd
644 : ! derivative at boundary to a non-zero value.
645 :
646 : ! periodic, not-a-knot, zero derivative, and divided-difference based
647 : ! BCs are "homogeneous"-- i.e. if splines s & t satisfy the BC then
648 : ! the spline (c*s + t) formed as a linear combination of these two
649 : ! splines, also satisfies the BC.
650 :
651 : ! algorithm note -- handling of inhomogeneous explicit BC's while
652 : ! maintaining full C2 differentiability is delicate. Basic method: use
653 : ! a fully C2 method based on the "not-a-knot" BC, and then, correct to
654 : ! meet each user BC by calculating a C2 spline that is zero at all grid
655 : ! points but satisfies a BC which is the difference btw the user spec
656 : ! and the not-a-knot result; add the coeffs of this into the original.
657 :
658 : ! for this more workspace is needed: nwk >= 4*inx*inth +5*max(inx,inth)
659 :
660 0 : subroutine bcspline_db(x,inx,th,inth,fspl,inf3,
661 0 : & ibcxmin,bcxmin,ibcxmax,bcxmax,
662 0 : & ibcthmin,bcthmin,ibcthmax,bcthmax,
663 0 : & wk,nwk,ilinx,ilinth,ier)
664 :
665 : implicit none
666 : integer inx, inth, inf3, nwk, ibcxmin, ibcxmax, ibcthmin, ibcthmax, ilinx,ilinth,ier
667 : real(dp) x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
668 : real(dp) bcxmin(inth),bcxmax(inth)
669 : real(dp) bcthmin(inx),bcthmax(inx)
670 :
671 : ! input:
672 : ! x(1...inx) -- abscissae, first dimension of data
673 : ! th(1...inth) -- abscissae, second dimension of data f(x,th)
674 : ! fspl(1,1,1..inx,1..inth) -- function values
675 : ! inf3 -- fspl dimensioning, inf3 >= inx required.
676 :
677 : ! boundary conditions input:
678 : ! ibcxmin -- indicator for boundary condition at x(1):
679 : ! bcxmin(...) -- boundary condition data
680 : ! =-1 -- periodic boundary condition
681 : ! =0 -- use "not a knot", bcxmin(...) ignored
682 : ! =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
683 : ! =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
684 : ! =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
685 : ! =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
686 : ! =5 -- match 1st derivative df/dx to 1st divided difference
687 : ! =6 -- match 2nd derivative d2f/dx2 to 2nd divided difference
688 : ! =7 -- match 3rd derivative d3f/dx3 3rd divided difference
689 : ! (for more detailed definition of BCs 5-7, see the
690 : ! comments of subroutine mkspline)
691 : ! NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
692 :
693 : ! ibcxmax -- indicator for boundary condition at x(nx):
694 : ! bcxmax(...) -- boundary condition data
695 : ! (interpretation as with ibcxmin, bcxmin)
696 : ! NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
697 :
698 : ! ibcthmin -- indicator for boundary condition at th(1):
699 : ! bcthmin(...) -- boundary condition data
700 : ! (interpretation as with ibcxmin, bcxmin)
701 : ! ibcthmax -- indicator for boundary condition at th(inth):
702 : ! bcthmax(...) -- boundary condition data
703 : ! (interpretation as with ibcxmin, bcxmin)
704 : ! NOTE: if ibcthmin=-1, ibcthmax is ignored! ...and the BC is periodic.
705 :
706 : ! NOTE the bcxmin,bcxmax,bcthmin,bcthmax arrays are only used if the
707 : ! corresponding boundary condition flags are set to 1 or 2.
708 : ! Carefully note the dimensioning of these arrays!
709 :
710 : ! output:
711 : ! fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
712 : ! ...fspl(1,1,*,*) is not replaced.
713 :
714 : ! ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
715 : ! ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
716 :
717 : ! ier -- completion code, 0 for normal
718 :
719 : ! workspace:
720 : ! wk -- must be at least 5*max(inx,inth) large
721 : ! 5*max(inx,inth) + 4*inx*inth large
722 : ! if explicit non-zero d/dth or d2/dth2 BC's
723 : ! are supplied.
724 : ! nwk -- size of workspace of workspace provided
725 :
726 : !---------------------------------
727 : ! in what follows, "f" is an abbreviation for "fspl"...
728 :
729 : ! compute bicubic spline of 2d function, given values at the grid
730 : ! grid crossing points, f(1,1,i,j)=f(x(i),th(j)).
731 :
732 : ! on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
733 : ! and th btw th(j) and th(j+1), dt=th-th(j),
734 :
735 : ! spline =
736 : ! 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)
737 : ! +dt*(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))
738 : ! +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))
739 : ! +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))
740 :
741 : ! where d2=dt**2 and d3=dt**3.
742 :
743 : integer iselect1(10)
744 : integer iselect2(10)
745 :
746 :
747 : integer iflg2, ix, itest, ierx, inxo, ith, jth, ii, iadr, ia5w, iaspl, ierth, intho, ic
748 : integer ibcthmina, ibcthmaxa, iasc, iinc, iawk, jx
749 0 : real(dp) xo2, xo6, zcur, zdiff1, zhxn, zhth, zdiff2
750 0 : real(dp) fval(6)
751 :
752 : !---------------------------------
753 :
754 : ! see if 2nd pass is needed due to "non-linear" d/dth bdy cond.
755 :
756 0 : iflg2=0
757 0 : if(ibcthmin /= -1) then
758 0 : if((ibcthmin == 1).or.(ibcthmin == 2)) then
759 0 : do ix=1,inx
760 0 : if (bcthmin(ix) /= 0.d0) iflg2=1
761 : end do
762 : end if
763 0 : if((ibcthmax == 1).or.(ibcthmax == 2)) then
764 0 : do ix=1,inx
765 0 : if (bcthmax(ix) /= 0.d0) iflg2=1
766 : end do
767 : end if
768 : end if
769 :
770 0 : ier=0
771 0 : itest=5*max(inx,inth)
772 0 : if(iflg2 == 1) then
773 0 : itest=itest +4*inx*inth
774 : end if
775 :
776 0 : if(nwk < itest) then
777 0 : write(6,9901) nwk,itest
778 : 9901 format(' ?bcspline: workspace too small:'/
779 : & ' user supplied: nwk=',i6,'; need at least: ',i6/
780 : & ' nwk=4*nx*ny +5*max(nx,ny) will work for any user'/
781 : & ' choice of bdy conditions.')
782 0 : ier=1
783 : end if
784 0 : if(inx < 4) then
785 0 : write(6,'('' ?bcspline: at least 4 x points required.'')')
786 0 : ier=1
787 : end if
788 0 : if(inth < 4) then
789 0 : write(6,'('' ?bcspline: need at least 4 theta points.'')')
790 0 : ier=1
791 : end if
792 :
793 0 : call ibc_ck_db(ibcxmin,'bcspline','xmin',-1,7,ier)
794 0 : if(ibcxmin >= 0) call ibc_ck_db(ibcxmax,'bcspline','xmax',0,7,ier)
795 0 : call ibc_ck_db(ibcthmin,'bcspline','thmin',-1,7,ier)
796 0 : if(ibcthmin >= 0) call ibc_ck_db(ibcthmax,'bcspline','thmax',0,7,ier)
797 :
798 : ! check ilinx & x vector
799 :
800 0 : call splinck_db(x,inx,ilinx,1.0d-3,ierx)
801 0 : if(ierx /= 0) ier=2
802 :
803 0 : if(ier == 2) then
804 0 : write(6,'('' ?bcspline: x axis not strict ascending'')')
805 : end if
806 :
807 : ! check ilinth & th vector
808 :
809 0 : call splinck_db(th,inth,ilinth,1.0d-3,ierth)
810 0 : if(ierth /= 0) ier=3
811 :
812 : ! if(ier == 3) then
813 : ! write(6,'('' ?bcspline: th axis not strict ascending'')')
814 : ! end if
815 :
816 0 : if(ier /= 0) return
817 :
818 : !------------------------------------
819 :
820 0 : xo2=0.5d0
821 0 : xo6=1.d0/6.d0
822 :
823 : ! spline in x direction
824 :
825 0 : inxo=4*(inx-1)
826 0 : do ith=1,inth
827 :
828 : ! copy the function in
829 :
830 0 : do ix=1,inx
831 0 : wk(4*(ix-1)+1)=fspl(1,1,ix,ith)
832 : end do
833 :
834 0 : if(ibcxmin == 1) then
835 0 : wk(2)=bcxmin(ith)
836 0 : else if(ibcxmin == 2) then
837 0 : wk(3)=bcxmin(ith)
838 : end if
839 :
840 0 : if(ibcxmax == 1) then
841 0 : wk(inxo+2)=bcxmax(ith)
842 0 : else if(ibcxmax == 2) then
843 0 : wk(inxo+3)=bcxmax(ith)
844 : end if
845 :
846 : ! use Wayne's routine
847 :
848 0 : call v_spline_db(ibcxmin,ibcxmax,inx,x,wk,wk(4*inx+1))
849 :
850 : ! copy the coefficients out
851 :
852 0 : do ix=1,inx
853 0 : fspl(2,1,ix,ith)=wk(4*(ix-1)+2)
854 0 : fspl(3,1,ix,ith)=wk(4*(ix-1)+3)*xo2
855 0 : fspl(4,1,ix,ith)=wk(4*(ix-1)+4)*xo6
856 : end do
857 :
858 : end do
859 :
860 : !-----------------------------------
861 :
862 : ! spline in theta direction
863 :
864 0 : intho=4*(inth-1)
865 0 : do ix=1,inx
866 :
867 : ! spline each x coeff
868 :
869 0 : do ic=1,4
870 :
871 : ! copy ordinates in
872 :
873 0 : do ith=1,inth
874 0 : wk(4*(ith-1)+1)=fspl(ic,1,ix,ith)
875 : end do
876 :
877 : ! first pass: use a linear BC -- if flag indicates BC correction
878 : ! will be needed, it will be done later
879 :
880 0 : wk(2)=0.d0
881 0 : wk(3)=0.d0
882 0 : wk(intho+2)=0.d0
883 0 : wk(intho+3)=0.d0
884 :
885 0 : ibcthmina=ibcthmin
886 0 : ibcthmaxa=ibcthmax
887 0 : if(iflg2 == 1) then
888 0 : if((ibcthmin == 1).or.(ibcthmin == 2)) ibcthmina=0
889 0 : if((ibcthmax == 1).or.(ibcthmax == 2)) ibcthmaxa=0
890 : end if
891 :
892 0 : call v_spline_db(ibcthmina,ibcthmaxa,inth,th,wk,wk(4*inth+1))
893 :
894 : ! copy coeffs out
895 :
896 0 : do ith=1,inth
897 0 : fspl(ic,2,ix,ith)=wk(4*(ith-1)+2)
898 0 : fspl(ic,3,ix,ith)=wk(4*(ith-1)+3)*xo2
899 0 : fspl(ic,4,ix,ith)=wk(4*(ith-1)+4)*xo6
900 : end do
901 :
902 : end do
903 :
904 : end do
905 :
906 : ! now make correction for user BC's if needed
907 :
908 0 : if(iflg2 == 1) then
909 :
910 0 : iasc=1 ! wk addr for correction splines
911 0 : iinc=4*inth ! spacing btw correction splines
912 0 : iawk=iasc+4*inth*inx
913 :
914 : ! last grid zone widths
915 :
916 0 : zhxn=x(inx)-x(inx-1)
917 0 : jx=inx-1
918 0 : zhth=th(inth)-th(inth-1)
919 0 : jth=inth-1
920 :
921 0 : do ii=1,10
922 0 : iselect1(ii)=0
923 0 : iselect2(ii)=0
924 : end do
925 : if(ibcthmin == 1) iselect1(3)=1
926 : if(ibcthmin == 2) iselect1(5)=1
927 0 : if(ibcthmax == 1) iselect2(3)=1
928 0 : if(ibcthmax == 2) iselect2(5)=1
929 :
930 : ! loop over BC's
931 :
932 0 : do ix=1,inx
933 :
934 : ! (a) d/dth @ th(1) difference btw current BC and user request
935 :
936 0 : if(ibcthmin == 1) then
937 0 : if(ix < inx) then
938 0 : zcur=fspl(1,2,ix,1) ! 1st deriv.
939 : else
940 : zcur=fspl(1,2,jx,1)+zhxn*(fspl(2,2,jx,1)+zhxn*
941 0 : & (fspl(3,2,jx,1)+zhxn*fspl(4,2,jx,1)))
942 : end if
943 0 : zdiff1=bcthmin(ix)-zcur
944 0 : else if(ibcthmin == 2) then
945 0 : if(ix < inx) then
946 0 : zcur=2.d0*fspl(1,3,ix,1) ! 2nd deriv.
947 : else
948 : zcur=2.d0*(fspl(1,3,jx,1)+zhxn*(fspl(2,3,jx,1)+zhxn*
949 0 : & (fspl(3,3,jx,1)+zhxn*fspl(4,3,jx,1))))
950 : end if
951 0 : zdiff1=bcthmin(ix)-zcur
952 : else
953 : zdiff1=0.d0
954 : end if
955 :
956 : ! (b) d/dth @ th(inth) difference btw current BC and user request
957 :
958 0 : if(ibcthmax == 1) then
959 0 : if(ix < inx) then
960 : ! 1st deriv.
961 : zcur=fspl(1,2,ix,jth)+zhth*(2.d0*fspl(1,3,ix,jth)+
962 0 : & zhth*3.d0*fspl(1,4,ix,jth))
963 : else
964 : call bcspeval_db(x(inx),th(inth),iselect2,fval,
965 0 : & x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
966 0 : zcur=fval(1)
967 0 : if(ier /= 0) return
968 : end if
969 0 : zdiff2=bcthmax(ix)-zcur
970 0 : else if(ibcthmax == 2) then
971 0 : if(ix < inx) then
972 : ! 2nd deriv.
973 : zcur=2.d0*fspl(1,3,ix,jth)+
974 0 : & 6.d0*zhth*fspl(1,4,ix,jth)
975 : else
976 : call bcspeval_db(x(inx),th(inth),iselect2,fval,
977 0 : & x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
978 0 : zcur=fval(1)
979 0 : if(ier /= 0) return
980 : end if
981 0 : zdiff2=bcthmax(ix)-zcur
982 : else
983 : zdiff2=0.d0
984 : end if
985 :
986 : ! ok compute the theta spline with BC's to span the difference(s)
987 : ! these theta "correction splines" are zero at all the grid points
988 : ! but have at least one non-zero 1st or 2nd derivative BC
989 :
990 0 : iadr=iasc+(ix-1)*iinc
991 0 : do ith=1,inth
992 0 : wk(iadr+4*(ith-1))=0.d0
993 : end do
994 :
995 0 : wk(iadr+1)=0.d0
996 0 : wk(iadr+2)=0.d0
997 0 : wk(iadr+intho+1)=0.d0
998 0 : wk(iadr+intho+2)=0.d0
999 :
1000 0 : if(ibcthmin == 1) then
1001 0 : wk(iadr+1)=zdiff1
1002 0 : else if(ibcthmin == 2) then
1003 0 : wk(iadr+2)=zdiff1
1004 : end if
1005 :
1006 0 : if(ibcthmax == 1) then
1007 0 : wk(iadr+intho+1)=zdiff2
1008 0 : else if(ibcthmax == 2) then
1009 0 : wk(iadr+intho+2)=zdiff2
1010 : end if
1011 :
1012 0 : call v_spline_db(ibcthmin,ibcthmax,inth,th,wk(iadr),wk(iawk))
1013 : end do
1014 :
1015 : ! add in results to main array -- th spline coef corrections
1016 :
1017 0 : do ix=1,inx
1018 0 : iadr=iasc+(ix-1)*iinc
1019 0 : do ith=1,inth-1
1020 0 : wk(iadr+4*(ith-1)+2)=wk(iadr+4*(ith-1)+2)*xo2
1021 0 : wk(iadr+4*(ith-1)+3)=wk(iadr+4*(ith-1)+3)*xo6
1022 0 : if(ix < inx) then
1023 0 : fspl(1,2,ix,ith)=fspl(1,2,ix,ith)+wk(iadr+4*(ith-1)+1)
1024 0 : fspl(1,3,ix,ith)=fspl(1,3,ix,ith)+wk(iadr+4*(ith-1)+2)
1025 0 : fspl(1,4,ix,ith)=fspl(1,4,ix,ith)+wk(iadr+4*(ith-1)+3)
1026 : end if
1027 : end do
1028 : end do
1029 :
1030 : ! compute the x splines of the th spline correction coeffs
1031 :
1032 0 : ia5w=iawk+4*inx
1033 :
1034 0 : do ith=1,inth-1
1035 0 : do ic=2,4
1036 0 : do ix=1,inx
1037 0 : iaspl=iasc+iinc*(ix-1)
1038 0 : wk(iawk+4*(ix-1))=wk(iaspl+4*(ith-1)+(ic-1))
1039 : end do
1040 :
1041 : ! use zero BCs for this correction spline
1042 :
1043 0 : wk(iawk+1)=0.d0
1044 0 : wk(iawk+2)=0.d0
1045 0 : wk(iawk+inxo+1)=0.d0
1046 0 : wk(iawk+inxo+2)=0.d0
1047 :
1048 : ! periodic spline of correction spline higher coeffs (1st coeffs are
1049 : ! all zero by defn of the correction spline
1050 :
1051 0 : call v_spline_db(ibcxmin,ibcxmax,inx,x,wk(iawk),wk(ia5w))
1052 :
1053 0 : do ix=1,inx-1
1054 : fspl(2,ic,ix,ith)=fspl(2,ic,ix,ith)+
1055 0 : & wk(iawk+4*(ix-1)+1)
1056 : fspl(3,ic,ix,ith)=fspl(3,ic,ix,ith)+
1057 0 : & wk(iawk+4*(ix-1)+2)*xo2
1058 : fspl(4,ic,ix,ith)=fspl(4,ic,ix,ith)+
1059 0 : & wk(iawk+4*(ix-1)+3)*xo6
1060 : end do
1061 :
1062 : end do
1063 : end do ! ith
1064 :
1065 : end if ! BC correction needs test
1066 :
1067 : return
1068 : end subroutine bcspline_db
1069 :
1070 : ! cspline -- dmc 15 Feb 1999
1071 :
1072 : ! a standard interface to the 1d spline setup routine
1073 : ! modified dmc 3 Mar 2000 -- to use Wayne Houlberg's v_spline code.
1074 : ! new BC options added.
1075 :
1076 6806634 : subroutine cspline_db(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
1077 6806634 : & wk,iwk,ilinx,ier)
1078 :
1079 : implicit none
1080 : integer nx, iwk
1081 : real(dp) x(nx) ! x axis (in)
1082 : real(dp) fspl(4,nx) ! spline data (in/out)
1083 : integer ibcxmin ! x(1) BC flag (in, see comments)
1084 : real(dp) bcxmin ! x(1) BC data (in, see comments)
1085 : integer ibcxmax ! x(nx) BC flag (in, see comments)
1086 : real(dp) bcxmax ! x(nx) BC data (in, see comments)
1087 : real(dp) wk(iwk) ! workspace of size at least nx
1088 : integer ilinx ! even spacing flag (out)
1089 : integer ier ! output, =0 means OK
1090 :
1091 : ! ** note wk(...) array is not used unless ibcxmin=-1 (periodic spline
1092 : ! evaluation)
1093 :
1094 : ! this routine computes spline coefficients for a 1d spline --
1095 : ! evaluation of the spline can be done by cspeval.for subroutines
1096 : ! or directly by inline code.
1097 :
1098 : ! the input x axis x(1...nx) must be strictly ascending, i.e.
1099 : ! x(i+1) > x(i) is required for i=1 to nx-1. This is checked and
1100 : ! ier=1 is set and the routine exits if the test is not satisfied.
1101 :
1102 : ! on output, ilinx=1 is set if, to a reasonably close tolerance,
1103 : ! all grid spacings x(i+1)-x(i) are equal. This allows a speedier
1104 : ! grid lookup algorithm on evaluation of the spline. If on output
1105 : ! ilinx=2, this means the spline x axis is not evenly spaced.
1106 :
1107 : ! the input data for the spline are given in f[j] = fspl(1,j). The
1108 : ! output data are the spline coefficients fspl(2,j),fspl(3,j), and
1109 : ! fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
1110 : ! boundary conditions and with the properties
1111 :
1112 : ! s(x(j)) = fspl(1,j)
1113 : ! s'(x) is continuous even at the grid points x(j)
1114 : ! s''(x) is continuous even at the grid points x(j)
1115 :
1116 : ! the formula for evaluation of s(x) is:
1117 :
1118 : ! let dx = x-x(i), where x(i) <= x <= x(i+1). Then,
1119 : ! s(x)=fspl(1,i) + dx*(fspl(2,i) +dx*(fspl(3,i) + dx*fspl(4,i)))
1120 :
1121 : ! ==>boundary conditions. Complete specification of a 1d spline
1122 : ! requires specification of boundary conditions at x(1) and x(nx).
1123 :
1124 : ! this routine provides 4 options:
1125 :
1126 : ! -1 ***** PERIODIC BC
1127 : ! ibcxmin=-1 -- periodic boundary condition. This means the
1128 : ! boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
1129 : ! are imposed. Note that s(x(1))=s(x(nx)) (i.e. fspl(1,1)=fspl(1,nx))
1130 : ! is not required -- that is determined by the fspl array input data.
1131 : ! The periodic boundary condition is to be preferred for periodic
1132 : ! data. When splining periodic data f(x) with period P, the relation
1133 : ! x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
1134 : ! (ibcxmax, bcxmin, bcxmax are ignored).
1135 :
1136 : ! if a periodic boundary condition is set, this covers both boundaries.
1137 : ! for the other types of boundary conditions, the type of condition
1138 : ! chosen for the x(1) boundary need not be the same as the type chosen
1139 : ! for the x(nx) boundary.
1140 :
1141 : ! 0 ***** NOT A KNOT BC
1142 : ! ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
1143 : ! condition -- see cubsplb.for. This is a common way for inferring
1144 : ! a "good" spline boundary condition automatically from data in the
1145 : ! vicinity of the boundary. (bcxmin | bcxmax are ignored).
1146 :
1147 : ! 1 ***** BC: SPECIFIED SLOPE
1148 : ! ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
1149 : ! s'(x(nx)) match the passed value (bcxmin | bcxmax).
1150 :
1151 : ! 2 ***** BC: SPECIFIED 2nd DERIVATIVE
1152 : ! ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
1153 : ! s''(x(nx)) match the passed value (bcxmin | bcxmax).
1154 :
1155 : ! 3 ***** BC: SPECIFIED SLOPE = 0.0
1156 : ! ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
1157 : ! s'(x(nx)) equal to ZERO.
1158 :
1159 : ! 4 ***** BC: SPECIFIED 2nd DERIVATIVE = 0.0
1160 : ! ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
1161 : ! s''(x(nx)) equal to ZERO.
1162 :
1163 : ! 5 ***** BC: 1st DIVIDED DIFFERENCE
1164 : ! ibcxmin=5 | ibcxmax=5 -- boundary condition is to have s'(x(1)) |
1165 : ! s'(x(nx)) equal to the slope from the 1st|last 2 points
1166 :
1167 : ! 6 ***** BC: 2nd DIVIDED DIFFERENCE
1168 : ! ibcxmin=6 | ibcxmax=6 -- boundary condition is to have s''(x(1)) |
1169 : ! s''(x(nx)) equal to the 2nd derivative from the 1st|last 3 points
1170 :
1171 : ! 7 ***** BC: 3rd DIVIDED DIFFERENCE
1172 : ! ibcxmin=7 | ibcxmax=7 -- boundary condition is to have s'''(x(1)) |
1173 : ! s'''(x(nx)) equal to the 3rd derivative from the 1st|last 4 points
1174 :
1175 : !---------------------------------------------------------------------
1176 : real(dp), parameter :: half = 0.5d0
1177 : real(dp), parameter :: sixth = 1d0/6d0
1178 : integer inum,i,ierx
1179 :
1180 : ! error checks
1181 :
1182 6806634 : ier = 0
1183 6806634 : if(nx < 4) then
1184 : ! write(6,'('' ?cspline: at least 4 x points required.'')')
1185 0 : ier=1
1186 : end if
1187 6806634 : call ibc_ck_db(ibcxmin,'cspline','xmin',-1,7,ier)
1188 6806634 : if(ibcxmin >= 0) call ibc_ck_db(ibcxmax,'cspline','xmax',0,7,ier)
1189 :
1190 : ! x axis check
1191 :
1192 6806634 : call splinck_db(x,nx,ilinx,1.0d-3,ierx)
1193 6806634 : if(ierx /= 0) ier=2
1194 :
1195 : ! if(ier == 2) then
1196 : ! write(6,'('' ?cspline: x axis not strict ascending'')')
1197 : ! end if
1198 :
1199 6806634 : if(ibcxmin == -1) then
1200 0 : inum=nx
1201 0 : if(iwk < inum) then
1202 : ! write(6,1009) inum,iwk,nx
1203 : ! 1009 format(
1204 : ! & ' ?cspline: workspace too small. need: ',i6,' got: ',i6/
1205 : ! & ' (need = nx, nx=',i6)
1206 0 : ier=3
1207 : end if
1208 : end if
1209 :
1210 6806634 : if(ier /= 0) return
1211 :
1212 : ! OK -- evaluate spline
1213 :
1214 6806634 : if(ibcxmin == 1) then
1215 0 : fspl(2,1)=bcxmin
1216 6806634 : else if(ibcxmin == 2) then
1217 0 : fspl(3,1)=bcxmin
1218 : end if
1219 :
1220 6806634 : if(ibcxmax == 1) then
1221 0 : fspl(2,nx)=bcxmax
1222 6806634 : else if(ibcxmax == 2) then
1223 0 : fspl(3,nx)=bcxmax
1224 : end if
1225 :
1226 6806634 : call v_spline_db(ibcxmin,ibcxmax,nx,x,fspl,wk)
1227 :
1228 1583109912 : do i=1,nx
1229 1576303278 : fspl(3,i)=half*fspl(3,i)
1230 1583109912 : fspl(4,i)=sixth*fspl(4,i)
1231 : end do
1232 :
1233 : return
1234 : end subroutine cspline_db
1235 :
1236 :
1237 0 : subroutine evbicub_db(xget,yget,x,nx,y,ny,ilinx,iliny,f1,inf2,ict,fval,ier)
1238 :
1239 : ! use mkbicub to set up spline coefficients!
1240 :
1241 : ! evaluate a 2d cubic Spline interpolant on a rectilinear
1242 : ! grid -- this is C2 in both directions.
1243 :
1244 : ! this subroutine calls two subroutines:
1245 : ! herm2xy -- find cell containing (xget,yget)
1246 : ! fvbicub -- evaluate interpolant function and (optionally) derivatives
1247 :
1248 : ! input arguments:
1249 : ! ================
1250 :
1251 : implicit none
1252 : integer nx,ny ! grid sizes
1253 : real(dp) xget,yget ! target of this interpolation
1254 : real(dp) x(nx) ! ordered x grid
1255 : real(dp) y(ny) ! ordered y grid
1256 : integer ilinx ! ilinx=1 => assume x evenly spaced
1257 : integer iliny ! iliny=1 => assume y evenly spaced
1258 :
1259 : integer inf2
1260 : real(dp), pointer :: f1(:)
1261 :
1262 : ! f 2nd dimension inf2 must be >= nx
1263 : ! contents of f:
1264 :
1265 : ! f(0,i,j) = f @ x(i),y(j)
1266 : ! f(1,i,j) = d2f/dx2 @ x(i),y(j)
1267 : ! f(2,i,j) = d2f/dy2 @ x(i),y(j)
1268 : ! f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
1269 :
1270 : ! (these are spline coefficients selected for continuous 2-
1271 : ! diffentiability, see mkbicub[w].for)
1272 :
1273 : integer ict(6) ! code specifying output desired
1274 :
1275 : ! ict(1)=1 -- return f (0, don't)
1276 : ! ict(2)=1 -- return df/dx (0, don't)
1277 : ! ict(3)=1 -- return df/dy (0, don't)
1278 : ! ict(4)=1 -- return d2f/dx2 (0, don't)
1279 : ! ict(5)=1 -- return d2f/dy2 (0, don't)
1280 : ! ict(6)=1 -- return d2f/dxdy (0, don't)
1281 : ! the number of non zero values ict(1:6)
1282 : ! determines the number of outputs...
1283 :
1284 : ! new dmc December 2005 -- access to higher derivatives (even if not
1285 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
1286 : ! if ict(1)=3 -- want 3rd derivatives
1287 : ! ict(2)=1 for d3f/dx3
1288 : ! ict(3)=1 for d3f/dx2dy
1289 : ! ict(4)=1 for d3f/dxdy2
1290 : ! ict(5)=1 for d3f/dy3
1291 : ! number of non-zero values ict(2:5) gives no. of outputs
1292 : ! if ict(1)=4 -- want 4th derivatives
1293 : ! ict(2)=1 for d4f/dx3dy
1294 : ! ict(3)=1 for d4f/dx2dy2
1295 : ! ict(4)=1 for d4f/dxdy3
1296 : ! number of non-zero values ict(2:4) gives no. of outputs
1297 : ! if ict(1)=5 -- want 5th derivatives
1298 : ! ict(2)=1 for d5f/dx3dy2
1299 : ! ict(3)=1 for d5f/dx2dy3
1300 : ! number of non-zero values ict(2:3) gives no. of outputs
1301 : ! if ict(1)=6 -- want 6th derivatives
1302 : ! d6f/dx3dy3 -- one value is returned.
1303 :
1304 : ! output arguments:
1305 : ! =================
1306 :
1307 : real(dp) fval(6) ! output data
1308 : integer ier ! error code =0 ==> no error
1309 :
1310 : ! fval(1) receives the first output (depends on ict(...) spec)
1311 : ! fval(2) receives the second output (depends on ict(...) spec)
1312 : ! fval(3) receives the third output (depends on ict(...) spec)
1313 : ! fval(4) receives the fourth output (depends on ict(...) spec)
1314 : ! fval(5) receives the fourth output (depends on ict(...) spec)
1315 : ! fval(6) receives the fourth output (depends on ict(...) spec)
1316 :
1317 : ! examples:
1318 : ! on input ict = [1,1,1,0,0,1]
1319 : ! on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
1320 :
1321 : ! on input ict = [1,0,0,0,0,0]
1322 : ! on output fval= [f] ... elements 2 -- 6 never referenced.
1323 :
1324 : ! on input ict = [0,0,0,1,1,0]
1325 : ! on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
1326 :
1327 : ! on input ict = [0,0,1,0,0,0]
1328 : ! on output fval= [df/dy] ... elements 2 -- 6 never referenced.
1329 :
1330 : ! ier -- completion code: 0 means OK
1331 : !-------------------
1332 : ! local:
1333 :
1334 : integer i,j ! cell indices
1335 :
1336 : ! normalized displacement from (x(i),y(j)) corner of cell.
1337 : ! xparam=0 @x(i) xparam=1 @x(i+1)
1338 : ! yparam=0 @y(j) yparam=1 @y(j+1)
1339 :
1340 0 : real(dp) xparam,yparam
1341 :
1342 : ! cell dimensions and
1343 : ! inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
1344 :
1345 0 : real(dp) hx,hy
1346 0 : real(dp) hxi,hyi
1347 :
1348 : ! 0 <= xparam <= 1
1349 : ! 0 <= yparam <= 1
1350 :
1351 : ! ** the interface is very similar to herm2ev.for; can use herm2xy **
1352 : !---------------------------------------------------------------------
1353 :
1354 : call herm2xy_db(xget,yget,x,nx,y,ny,ilinx,iliny,
1355 0 : & i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
1356 0 : if(ier /= 0) return
1357 :
1358 : call fvbicub_db(ict,1,1,
1359 : & fval,(/i/),(/j/),(/xparam/),(/yparam/),
1360 0 : & (/hx/),(/hxi/),(/hy/),(/hyi/),f1,inf2,ny)
1361 :
1362 0 : return
1363 : end subroutine evbicub_db
1364 :
1365 : !---------------------------------------------------------------------
1366 : ! evaluate C1 cubic Hermite function interpolation -- 2d fcn
1367 : ! --vectorized-- dmc 10 Feb 1999
1368 :
1369 : ! use mkbicub to set up spline coefficients!
1370 :
1371 140 : subroutine fvbicub_db(ict,ivec,ivecd,
1372 140 : & fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
1373 : & f1,inf2,ny)
1374 :
1375 : integer ict(6) ! requested output control
1376 : integer ivec ! vector length
1377 : integer ivecd ! vector dimension (1st dim of fval)
1378 :
1379 : integer ii(ivec),jj(ivec) ! target cells (i,j)
1380 : real(dp) xparam(ivec),yparam(ivec)
1381 : ! normalized displacements from (i,j) corners
1382 :
1383 : real(dp) hx(ivec),hy(ivec) ! grid spacing, and
1384 : real(dp) hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
1385 : ! & 1/(y(j+1)-y(j))
1386 :
1387 : real(dp), pointer :: f1(:)
1388 :
1389 : real(dp) fval(ivecd,6) ! output returned
1390 :
1391 : ! for detailed description of fin, ict and fval see subroutine
1392 : ! evbicub comments. Note ict is not vectorized; the same output
1393 : ! is expected to be returned for all input vector data points.
1394 :
1395 : ! note that the index inputs ii,jj and parameter inputs
1396 : ! xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
1397 : ! output array fval has a vector ** 1st dimension ** whose
1398 : ! size must be given as a separate argument
1399 :
1400 : ! to use this routine in scalar mode, pass in ivec=ivecd=1
1401 :
1402 : C---------------
1403 : ! Spline evaluation consists of a "mixing" of the interpolant
1404 : ! data using the linear functionals xparam, xpi = 1-xparam,
1405 : ! yparam, ypi = 1-yparam, and the cubic functionals
1406 : ! xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
1407 : ! and their derivatives as needed.
1408 :
1409 : integer v
1410 140 : real(dp) sum
1411 :
1412 : real(dp), parameter :: sixth = 1.d0/6.d0
1413 :
1414 140 : real(dp), pointer :: fin(:,:,:)
1415 140 : fin(0:3,1:inf2,1:ny) => f1(1:4*inf2*ny)
1416 :
1417 : C---------------
1418 : ! ...in x direction
1419 :
1420 140 : z36th=sixth*sixth
1421 140 : iadr=0
1422 :
1423 140 : if(ict(1) <= 2) then
1424 :
1425 : ! get desired values:
1426 :
1427 140 : if(ict(1) == 1) then
1428 :
1429 : ! function value:
1430 :
1431 280 : iadr=iadr+1
1432 280 : do v=1,ivec
1433 140 : i=ii(v)
1434 140 : j=jj(v)
1435 :
1436 : ! in x direction...
1437 :
1438 280 : xp=xparam(v)
1439 280 : xpi=1.d0-xp
1440 280 : xp2=xp*xp
1441 280 : xpi2=xpi*xpi
1442 :
1443 280 : cx=xp*(xp2-1.d0)
1444 280 : cxi=xpi*(xpi2-1.d0)
1445 280 : hx2=hx(v)*hx(v)
1446 :
1447 : ! ...and in y direction
1448 :
1449 280 : yp=yparam(v)
1450 280 : ypi=1.d0-yp
1451 280 : yp2=yp*yp
1452 280 : ypi2=ypi*ypi
1453 :
1454 280 : cy=yp*(yp2-1.d0)
1455 280 : cyi=ypi*(ypi2-1.d0)
1456 280 : hy2=hy(v)*hy(v)
1457 :
1458 : sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
1459 140 : & xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
1460 :
1461 : sum=sum+sixth*hx2*(
1462 : & cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1463 140 : & cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1464 :
1465 : sum=sum+sixth*hy2*(
1466 : & xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
1467 140 : & xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
1468 :
1469 : sum=sum+z36th*hx2*hy2*(
1470 : & cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1471 140 : & cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1472 :
1473 280 : fval(v,iadr)=sum
1474 : end do
1475 : end if
1476 :
1477 140 : if(ict(2) == 1) then
1478 :
1479 : ! df/dx:
1480 :
1481 134 : iadr=iadr+1
1482 268 : do v=1,ivec
1483 134 : i=ii(v)
1484 134 : j=jj(v)
1485 :
1486 : ! in x direction...
1487 :
1488 134 : xp=xparam(v)
1489 134 : xpi=1.d0-xp
1490 134 : xp2=xp*xp
1491 134 : xpi2=xpi*xpi
1492 :
1493 274 : cxd=3.d0*xp2-1.d0
1494 274 : cxdi=-3.d0*xpi2+1.d0
1495 :
1496 : ! ...and in y direction
1497 :
1498 134 : yp=yparam(v)
1499 134 : ypi=1.d0-yp
1500 134 : yp2=yp*yp
1501 134 : ypi2=ypi*ypi
1502 :
1503 134 : cy=yp*(yp2-1.d0)
1504 134 : cyi=ypi*(ypi2-1.d0)
1505 134 : hy2=hy(v)*hy(v)
1506 :
1507 : sum=hxi(v)*(
1508 : & -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
1509 134 : & +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
1510 :
1511 : sum=sum+sixth*hx(v)*(
1512 : & cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1513 134 : & cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1514 :
1515 : sum=sum+sixth*hxi(v)*hy2*(
1516 : & -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
1517 134 : & +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
1518 :
1519 : sum=sum+z36th*hx(v)*hy2*(
1520 : & cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1521 134 : & cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1522 :
1523 268 : fval(v,iadr)=sum
1524 : end do
1525 : end if
1526 :
1527 140 : if(ict(3) == 1) then
1528 :
1529 : ! df/dy:
1530 :
1531 134 : iadr=iadr+1
1532 268 : do v=1,ivec
1533 134 : i=ii(v)
1534 134 : j=jj(v)
1535 :
1536 : ! in x direction...
1537 :
1538 134 : xp=xparam(v)
1539 134 : xpi=1.d0-xp
1540 134 : xp2=xp*xp
1541 134 : xpi2=xpi*xpi
1542 :
1543 134 : cx=xp*(xp2-1.d0)
1544 134 : cxi=xpi*(xpi2-1.d0)
1545 134 : hx2=hx(v)*hx(v)
1546 :
1547 : ! ...and in y direction
1548 :
1549 134 : yp=yparam(v)
1550 134 : ypi=1.d0-yp
1551 134 : yp2=yp*yp
1552 134 : ypi2=ypi*ypi
1553 :
1554 274 : cyd=3.d0*yp2-1.d0
1555 274 : cydi=-3.d0*ypi2+1.d0
1556 :
1557 : sum=hyi(v)*(
1558 : & xpi*(-fin(0,i,j) +fin(0,i,j+1))+
1559 134 : & xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
1560 :
1561 : sum=sum+sixth*hx2*hyi(v)*(
1562 : & cxi*(-fin(1,i,j) +fin(1,i,j+1))+
1563 134 : & cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1564 :
1565 : sum=sum+sixth*hy(v)*(
1566 : & xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
1567 134 : & xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1568 :
1569 : sum=sum+z36th*hx2*hy(v)*(
1570 : & cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1571 134 : & cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1572 :
1573 268 : fval(v,iadr)=sum
1574 : end do
1575 : end if
1576 :
1577 140 : if(ict(4) == 1) then
1578 :
1579 : ! d2f/dx2:
1580 :
1581 0 : iadr=iadr+1
1582 0 : do v=1,ivec
1583 0 : i=ii(v)
1584 0 : j=jj(v)
1585 :
1586 : ! in x direction...
1587 :
1588 0 : xp=xparam(v)
1589 0 : xpi=1.d0-xp
1590 :
1591 : ! ...and in y direction
1592 :
1593 0 : yp=yparam(v)
1594 0 : ypi=1.d0-yp
1595 0 : yp2=yp*yp
1596 0 : ypi2=ypi*ypi
1597 :
1598 0 : cy=yp*(yp2-1.d0)
1599 0 : cyi=ypi*(ypi2-1.d0)
1600 0 : hy2=hy(v)*hy(v)
1601 :
1602 : sum=(
1603 : & xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
1604 0 : & xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1605 :
1606 : sum=sum+sixth*hy2*(
1607 : & xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
1608 0 : & xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1609 :
1610 0 : fval(v,iadr)=sum
1611 : end do
1612 : end if
1613 :
1614 140 : if(ict(5) == 1) then
1615 :
1616 : ! d2f/dy2:
1617 :
1618 0 : iadr=iadr+1
1619 0 : do v=1,ivec
1620 0 : i=ii(v)
1621 0 : j=jj(v)
1622 :
1623 : ! in x direction...
1624 :
1625 0 : xp=xparam(v)
1626 0 : xpi=1.d0-xp
1627 0 : xp2=xp*xp
1628 0 : xpi2=xpi*xpi
1629 :
1630 0 : cx=xp*(xp2-1.d0)
1631 0 : cxi=xpi*(xpi2-1.d0)
1632 0 : hx2=hx(v)*hx(v)
1633 :
1634 : ! ...and in y direction
1635 :
1636 0 : yp=yparam(v)
1637 0 : ypi=1.d0-yp
1638 :
1639 : sum=(
1640 : & xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
1641 0 : & xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1642 :
1643 : sum=sum+sixth*hx2*(
1644 : & cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1645 0 : & cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1646 :
1647 0 : fval(v,iadr)=sum
1648 : end do
1649 : end if
1650 :
1651 140 : if(ict(6) == 1) then
1652 :
1653 : ! d2f/dxdy:
1654 :
1655 0 : iadr=iadr+1
1656 0 : do v=1,ivec
1657 0 : i=ii(v)
1658 0 : j=jj(v)
1659 :
1660 : ! in x direction...
1661 :
1662 0 : xp=xparam(v)
1663 0 : xpi=1.d0-xp
1664 0 : xp2=xp*xp
1665 0 : xpi2=xpi*xpi
1666 :
1667 0 : cxd=3.d0*xp2-1.d0
1668 0 : cxdi=-3.d0*xpi2+1.d0
1669 :
1670 : ! ...and in y direction
1671 :
1672 0 : yp=yparam(v)
1673 0 : ypi=1.d0-yp
1674 0 : yp2=yp*yp
1675 0 : ypi2=ypi*ypi
1676 :
1677 0 : cyd=3.d0*yp2-1.d0
1678 0 : cydi=-3.d0*ypi2+1.d0
1679 :
1680 : sum=hxi(v)*hyi(v)*(
1681 : & fin(0,i,j) -fin(0,i,j+1)
1682 0 : & -fin(0,i+1,j)+fin(0,i+1,j+1))
1683 :
1684 : sum=sum+sixth*hx(v)*hyi(v)*(
1685 : & cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
1686 0 : & cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1687 :
1688 : sum=sum+sixth*hxi(v)*hy(v)*(
1689 : & -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
1690 0 : & +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1691 :
1692 : sum=sum+z36th*hx(v)*hy(v)*(
1693 : & cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1694 0 : & cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1695 :
1696 0 : fval(v,iadr)=sum
1697 : end do
1698 : end if
1699 :
1700 : !-------------------------------------------------
1701 :
1702 0 : else if(ict(1) == 3) then
1703 0 : if(ict(2) == 1) then
1704 : ! evaluate d3f/dx3 (not continuous)
1705 0 : iadr=iadr+1
1706 0 : do v=1,ivec
1707 0 : i=ii(v)
1708 0 : j=jj(v)
1709 0 : yp=yparam(v)
1710 0 : ypi=1.d0-yp
1711 0 : yp2=yp*yp
1712 0 : ypi2=ypi*ypi
1713 0 : cy=yp*(yp2-1.d0)
1714 0 : cyi=ypi*(ypi2-1.d0)
1715 0 : hy2=hy(v)*hy(v)
1716 : sum=hxi(v)*(
1717 : & -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
1718 0 : & +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1719 :
1720 : sum=sum+sixth*hy2*hxi(v)*(
1721 : & -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
1722 0 : & +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1723 :
1724 0 : fval(v,iadr)=sum
1725 : end do
1726 : end if
1727 :
1728 0 : if(ict(3) == 1) then
1729 : ! evaluate d3f/dx2dy
1730 0 : iadr=iadr+1
1731 0 : do v=1,ivec
1732 0 : i=ii(v)
1733 0 : j=jj(v)
1734 0 : xp=xparam(v)
1735 0 : xpi=1.d0-xp
1736 0 : yp=yparam(v)
1737 0 : ypi=1.d0-yp
1738 0 : yp2=yp*yp
1739 0 : ypi2=ypi*ypi
1740 0 : cyd=3.d0*yp2-1.d0
1741 0 : cydi=-3.d0*ypi2+1.d0
1742 :
1743 : sum=hyi(v)*(
1744 : & xpi*(-fin(1,i,j) +fin(1,i,j+1))+
1745 0 : & xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
1746 :
1747 : sum=sum+sixth*hy(v)*(
1748 : & xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1749 0 : & xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1750 :
1751 0 : fval(v,iadr)=sum
1752 : end do
1753 : end if
1754 :
1755 0 : if(ict(4) == 1) then
1756 : ! evaluate d3f/dxdy2
1757 0 : iadr=iadr+1
1758 0 : do v=1,ivec
1759 0 : i=ii(v)
1760 0 : j=jj(v)
1761 0 : xp=xparam(v)
1762 0 : xpi=1.d0-xp
1763 0 : xp2=xp*xp
1764 0 : xpi2=xpi*xpi
1765 0 : cxd=3.d0*xp2-1.d0
1766 0 : cxdi=-3.d0*xpi2+1.d0
1767 0 : yp=yparam(v)
1768 0 : ypi=1.d0-yp
1769 :
1770 : sum=hxi(v)*(
1771 : & -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
1772 0 : & +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1773 :
1774 : sum=sum+sixth*hx(v)*(
1775 : & cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1776 0 : & cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1777 :
1778 0 : fval(v,iadr)=sum
1779 : end do
1780 : end if
1781 :
1782 0 : if(ict(5) == 1) then
1783 : ! evaluate d3f/dy3 (not continuous)
1784 0 : iadr=iadr+1
1785 0 : do v=1,ivec
1786 0 : i=ii(v)
1787 0 : j=jj(v)
1788 :
1789 0 : xp=xparam(v)
1790 0 : xpi=1.d0-xp
1791 0 : xp2=xp*xp
1792 0 : xpi2=xpi*xpi
1793 :
1794 0 : cx=xp*(xp2-1.d0)
1795 0 : cxi=xpi*(xpi2-1.d0)
1796 0 : hx2=hx(v)*hx(v)
1797 :
1798 : sum=hyi(v)*(
1799 : & xpi*(-fin(2,i,j) +fin(2,i,j+1))+
1800 0 : & xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
1801 :
1802 : sum=sum+sixth*hx2*hyi(v)*(
1803 : & cxi*(-fin(3,i,j) +fin(3,i,j+1))+
1804 0 : & cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1805 :
1806 0 : fval(v,iadr)=sum
1807 : end do
1808 : end if
1809 :
1810 : !-----------------------------------
1811 : ! access to 4th derivatives
1812 :
1813 0 : else if(ict(1) == 4) then
1814 0 : if(ict(2) == 1) then
1815 : ! evaluate d4f/dx3dy (not continuous)
1816 0 : iadr=iadr+1
1817 0 : do v=1,ivec
1818 0 : i=ii(v)
1819 0 : j=jj(v)
1820 0 : yp=yparam(v)
1821 0 : ypi=1.d0-yp
1822 0 : yp2=yp*yp
1823 0 : ypi2=ypi*ypi
1824 0 : cyd=3.d0*yp2-1.d0
1825 0 : cydi=-3.d0*ypi2+1.d0
1826 :
1827 : sum=hxi(v)*hyi(v)*(
1828 : & +( fin(1,i,j) -fin(1,i,j+1))
1829 0 : & +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1830 :
1831 : sum=sum+sixth*hy(v)*hxi(v)*(
1832 : & -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
1833 0 : & +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1834 :
1835 0 : fval(v,iadr)=sum
1836 : end do
1837 : end if
1838 :
1839 0 : if(ict(3) == 1) then
1840 : ! evaluate d4f/dx2dy2
1841 0 : iadr=iadr+1
1842 0 : do v=1,ivec
1843 0 : i=ii(v)
1844 0 : j=jj(v)
1845 :
1846 0 : xp=xparam(v)
1847 0 : xpi=1.d0-xp
1848 0 : yp=yparam(v)
1849 0 : ypi=1.d0-yp
1850 :
1851 : sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1852 0 : & xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
1853 :
1854 0 : fval(v,iadr)=sum
1855 : end do
1856 : end if
1857 :
1858 0 : if(ict(4) == 1) then
1859 : ! evaluate d4f/dxdy3 (not continuous)
1860 0 : iadr=iadr+1
1861 0 : do v=1,ivec
1862 0 : i=ii(v)
1863 0 : j=jj(v)
1864 :
1865 0 : xp=xparam(v)
1866 0 : xpi=1.d0-xp
1867 0 : xp2=xp*xp
1868 0 : xpi2=xpi*xpi
1869 :
1870 0 : cxd=3.d0*xp2-1.d0
1871 0 : cxdi=-3.d0*xpi2+1.d0
1872 :
1873 : sum=hyi(v)*hxi(v)*(
1874 : & +( fin(2,i,j) -fin(2,i,j+1))
1875 0 : & +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
1876 :
1877 : sum=sum+sixth*hx(v)*hyi(v)*(
1878 : & cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
1879 0 : & cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1880 :
1881 0 : fval(v,iadr)=sum
1882 : end do
1883 : end if
1884 :
1885 : !-----------------------------------
1886 : ! access to 5th derivatives
1887 :
1888 0 : else if(ict(1) == 5) then
1889 0 : if(ict(2) == 1) then
1890 : ! evaluate d5f/dx3dy2 (not continuous)
1891 0 : iadr=iadr+1
1892 0 : do v=1,ivec
1893 0 : i=ii(v)
1894 0 : j=jj(v)
1895 :
1896 0 : yp=yparam(v)
1897 0 : ypi=1.d0-yp
1898 :
1899 : sum=hxi(v)*(
1900 : & -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
1901 0 : & +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1902 :
1903 0 : fval(v,iadr)=sum
1904 : end do
1905 : end if
1906 :
1907 0 : if(ict(3) == 1) then
1908 : ! evaluate d5f/dx2dy3 (not continuous)
1909 0 : iadr=iadr+1
1910 0 : do v=1,ivec
1911 0 : i=ii(v)
1912 0 : j=jj(v)
1913 :
1914 0 : xp=xparam(v)
1915 0 : xpi=1.d0-xp
1916 :
1917 : sum=hyi(v)*(
1918 : & xpi*(-fin(3,i,j) +fin(3,i,j+1))+
1919 0 : & xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1920 :
1921 0 : fval(v,iadr)=sum
1922 : end do
1923 : end if
1924 :
1925 : !-----------------------------------
1926 : ! access to 6th derivatives
1927 :
1928 0 : else if(ict(1) == 6) then
1929 : ! evaluate d6f/dx3dy3 (not continuous)
1930 0 : iadr=iadr+1
1931 0 : do v=1,ivec
1932 0 : i=ii(v)
1933 0 : j=jj(v)
1934 : sum=hxi(v)*hyi(v)*(
1935 : & +( fin(3,i,j) -fin(3,i,j+1))
1936 0 : & +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1937 0 : fval(v,iadr)=sum
1938 : end do
1939 : end if
1940 :
1941 140 : return
1942 140 : end subroutine fvbicub_db
1943 :
1944 :
1945 0 : subroutine herm2ev_db(xget,yget,x,nx,y,ny,ilinx,iliny,f,inf2,ict,fval,ier)
1946 :
1947 : ! evaluate a 2d cubic Hermite interpolant on a rectilinear
1948 : ! grid -- this is C1 in both directions.
1949 :
1950 : ! this subroutine calls two subroutines:
1951 : ! herm2xy -- find cell containing (xget,yget)
1952 : ! herm2fcn -- evaluate interpolant function and (optionally) derivatives
1953 :
1954 : ! input arguments:
1955 : ! ================
1956 :
1957 : real(dp) xget,yget ! target of this interpolation
1958 : real(dp) x(nx) ! ordered x grid
1959 : real(dp) y(ny) ! ordered y grid
1960 : integer ilinx ! ilinx=1 => assume x evenly spaced
1961 : integer iliny ! iliny=1 => assume y evenly spaced
1962 :
1963 : real(dp) f(0:3,inf2,ny) ! function data
1964 :
1965 : ! f 2nd dimension inf2 must be >= nx
1966 : ! contents of f:
1967 :
1968 : ! f(0,i,j) = f @ x(i),y(j)
1969 : ! f(1,i,j) = df/dx @ x(i),y(j)
1970 : ! f(2,i,j) = df/dy @ x(i),y(j)
1971 : ! f(3,i,j) = d2f/dxdy @ x(i),y(j)
1972 :
1973 : integer ict(4) ! code specifying output desired
1974 :
1975 : ! ict(1)=1 -- return f (0, don't)
1976 : ! ict(2)=1 -- return df/dx (0, don't)
1977 : ! ict(3)=1 -- return df/dy (0, don't)
1978 : ! ict(4)=1 -- return d2f/dxdy (0, don't)
1979 :
1980 : ! output arguments:
1981 : ! =================
1982 :
1983 : real(dp) fval(4) ! output data
1984 : integer ier ! error code =0 ==> no error
1985 :
1986 : ! fval(1) receives the first output (depends on ict(...) spec)
1987 : ! fval(2) receives the second output (depends on ict(...) spec)
1988 : ! fval(3) receives the third output (depends on ict(...) spec)
1989 : ! fval(4) receives the fourth output (depends on ict(...) spec)
1990 :
1991 : ! examples:
1992 : ! on input ict = [1,1,1,1]
1993 : ! on output fval= [f,df/dx,df/dy,d2f/dxdy]
1994 :
1995 : ! on input ict = [1,0,0,0]
1996 : ! on output fval= [f] ... elements 2 & 3 & 4 never referenced
1997 :
1998 : ! on input ict = [0,1,1,0]
1999 : ! on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced
2000 :
2001 : ! on input ict = [0,0,1,0]
2002 : ! on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced.
2003 :
2004 : ! ier -- completion code: 0 means OK
2005 : !-------------------
2006 : ! local:
2007 :
2008 : integer i,j ! cell indices
2009 :
2010 : ! normalized displacement from (x(i),y(j)) corner of cell.
2011 : ! xparam=0 @x(i) xparam=1 @x(i+1)
2012 : ! yparam=0 @y(j) yparam=1 @y(j+1)
2013 :
2014 0 : real(dp) xparam,yparam
2015 :
2016 : ! cell dimensions and
2017 : ! inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
2018 :
2019 0 : real(dp) hx,hy
2020 0 : real(dp) hxi,hyi
2021 :
2022 : ! 0 <= xparam <= 1
2023 : ! 0 <= yparam <= 1
2024 :
2025 : !---------------------------------------------------------------------
2026 :
2027 : call herm2xy_db(xget,yget,x,nx,y,ny,ilinx,iliny,
2028 0 : & i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
2029 0 : if(ier /= 0) return
2030 :
2031 : call herm2fcn_db(ict,1,1,
2032 : & fval,(/i/),(/j/),(/xparam/),(/yparam/),
2033 0 : & (/hx/),(/hxi/),(/hy/),(/hyi/),f,inf2,ny)
2034 :
2035 0 : return
2036 : end subroutine herm2ev_db
2037 :
2038 : !---------------------------------------------------------------------
2039 : ! herm2xy -- look up x-y zone
2040 :
2041 : ! this is the "first part" of herm2ev, see comments, above.
2042 :
2043 140 : subroutine herm2xy_db(xget,yget,x,nx,y,ny,ilinx,iliny,
2044 : & i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
2045 :
2046 : ! input of herm2xy
2047 : ! ================
2048 : implicit none
2049 :
2050 : integer nx,ny ! array dimensions
2051 :
2052 : real(dp) xget,yget ! target point
2053 : real(dp) x(:) ! (nx) ! indep. coords., strict ascending
2054 : real(dp) y(:) ! (ny) ! indep. coords., strict ascending
2055 :
2056 : integer ilinx ! =1: x evenly spaced
2057 : integer iliny ! =1: y evenly spaced
2058 :
2059 : ! output of herm2xy
2060 : ! =================
2061 : integer i,j ! index to cell containing target pt
2062 : ! on exit: 1 <= i <= nx-1 1 <= j <= ny-1
2063 :
2064 : ! normalized position w/in (i,j) cell (btw 0 and 1):
2065 :
2066 : real(dp) xparam ! (xget-x(i))/(x(i+1)-x(i))
2067 : real(dp) yparam ! (yget-y(j))/(y(j+1)-y(j))
2068 :
2069 : ! grid spacing
2070 :
2071 : real(dp) hx ! hx = x(i+1)-x(i)
2072 : real(dp) hy ! hy = y(j+1)-y(j)
2073 :
2074 : ! inverse grid spacing:
2075 :
2076 : real(dp) hxi ! 1/hx = 1/(x(i+1)-x(i))
2077 : real(dp) hyi ! 1/hy = 1/(y(j+1)-y(j))
2078 :
2079 : integer ier ! return ier /= 0 on error
2080 :
2081 : !------------------------------------
2082 140 : real(dp) zxget,zyget,zxtol,zytol
2083 : integer nxm,nym,ii,jj
2084 :
2085 :
2086 140 : ier=0
2087 :
2088 : ! range check
2089 :
2090 140 : zxget=xget
2091 140 : zyget=yget
2092 140 : if((xget < x(1)).or.(xget > x(nx))) then
2093 0 : zxtol=4.0d-7*max(abs(x(1)),abs(x(nx)))
2094 0 : if((xget < x(1)-zxtol).or.(xget > x(nx)+zxtol)) then
2095 0 : ier=1
2096 : !write(6,1001) xget,x(1),x(nx)
2097 : ! 1001 format(' ?herm2ev: xget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
2098 : else
2099 : !if((xget < x(1)-0.5*zxtol).or.(xget > x(nx)+0.5*zxtol)) write(6,1011) xget,x(1),x(nx)
2100 : ! 1011 format(' %herm2ev: xget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
2101 0 : if(xget < x(1)) then
2102 0 : zxget=x(1)
2103 : else
2104 0 : zxget=x(nx)
2105 : end if
2106 : end if
2107 : end if
2108 140 : if((yget < y(1)).or.(yget > y(ny))) then
2109 0 : zytol=4.0d-7*max(abs(y(1)),abs(y(ny)))
2110 0 : if((yget < y(1)-zytol).or.(yget > y(ny)+zytol)) then
2111 0 : ier=1
2112 : ! write(6,1002) yget,y(1),y(ny)
2113 : ! 1002 format(' ?herm2ev: yget=',1pe11.4,' out of range ', 1pe11.4,' to ',1pe11.4)
2114 : else
2115 : ! if((yget < y(1)-0.5*zytol).or. (yget > y(ny)+0.5*zytol)) write(6,1012) yget,y(1),y(ny)
2116 : ! 1012 format(' %herm2ev: yget=',1pe15.8,' beyond range ', 1pe15.8,' to ',1pe15.8,' (fixup applied)')
2117 0 : if(yget < y(1)) then
2118 0 : zyget=y(1)
2119 : else
2120 0 : zyget=y(ny)
2121 : end if
2122 : end if
2123 : end if
2124 140 : if(ier /= 0) return
2125 :
2126 : ! now find interval in which target point lies..
2127 :
2128 140 : nxm=nx-1
2129 140 : nym=ny-1
2130 :
2131 140 : if(ilinx == 1) then
2132 134 : ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
2133 134 : i=min(nxm, ii)
2134 134 : if(zxget < x(i)) then
2135 0 : i=i-1
2136 134 : else if(zxget > x(i+1)) then
2137 0 : i=i+1
2138 : end if
2139 : else
2140 6 : if((1 <= i).and.(i < nxm)) then
2141 4 : if((x(i) <= zxget).and.(zxget <= x(i+1))) then
2142 : continue ! already have the zone
2143 : else
2144 4 : call zonfind_db(x,nx,zxget,i)
2145 : end if
2146 : else
2147 2 : call zonfind_db(x,nx,zxget,i)
2148 : end if
2149 : end if
2150 :
2151 140 : if(iliny == 1) then
2152 20 : jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
2153 20 : j=min(nym, jj)
2154 20 : if(zyget < y(j)) then
2155 0 : j=j-1
2156 20 : else if(zyget > y(j+1)) then
2157 0 : j=j+1
2158 : end if
2159 : else
2160 120 : if((1 <= j).and.(j < nym)) then
2161 102 : if((y(j) <= zyget).and.(zyget <= y(j+1))) then
2162 : continue ! already have the zone
2163 : else
2164 0 : call zonfind_db(y,ny,zyget,j)
2165 : end if
2166 : else
2167 18 : call zonfind_db(y,ny,zyget,j)
2168 : end if
2169 : end if
2170 :
2171 140 : hx=(x(i+1)-x(i))
2172 140 : hy=(y(j+1)-y(j))
2173 :
2174 140 : hxi=1.d0/hx
2175 140 : hyi=1.d0/hy
2176 :
2177 140 : xparam=(zxget-x(i))*hxi
2178 140 : yparam=(zyget-y(j))*hyi
2179 :
2180 140 : return
2181 : end subroutine herm2xy_db
2182 :
2183 :
2184 : !---------------------------------------------------------------------
2185 : ! evaluate C1 cubic Hermite function interpolation -- 2d fcn
2186 : ! --vectorized-- dmc 10 Feb 1999
2187 :
2188 0 : subroutine herm2fcn_db(ict,ivec,ivecd,
2189 0 : & fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
2190 0 : & fin,inf2,ny)
2191 :
2192 : integer ict(4) ! requested output control
2193 : integer ivec ! vector length
2194 : integer ivecd ! vector dimension (1st dim of fval)
2195 :
2196 : integer ii(ivec),jj(ivec) ! target cells (i,j)
2197 : real(dp) xparam(ivec),yparam(ivec)
2198 : ! normalized displacements from (i,j) corners
2199 :
2200 : real(dp) hx(ivec),hy(ivec) ! grid spacing, and
2201 : real(dp) hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
2202 : ! & 1/(y(j+1)-y(j))
2203 :
2204 : real(dp) fin(0:3,inf2,ny) ! interpolant data (cf "herm2ev")
2205 :
2206 : real(dp) fval(ivecd,4) ! output returned
2207 :
2208 : ! for detailed description of fin, ict and fval see subroutine
2209 : ! herm2ev comments. Note ict is not vectorized; the same output
2210 : ! is expected to be returned for all input vector data points.
2211 :
2212 : ! note that the index inputs ii,jj and parameter inputs
2213 : ! xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
2214 : ! output array fval has a vector ** 1st dimension ** whose
2215 : ! size must be given as a separate argument
2216 :
2217 : ! to use this routine in scalar mode, pass in ivec=ivecd=1
2218 :
2219 : !---------------
2220 : ! Hermite cubic basis functions
2221 : ! -->for function value matching
2222 : ! a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
2223 : ! abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
2224 :
2225 : ! a(x)=-2*x**3 + 3*x**2 = x*x*(-2.d0*x+3.0)
2226 : ! abar(x)=1-a(x)
2227 : ! a'(x)=-abar'(x) = 6.d0*x*(1.0-x)
2228 :
2229 : ! -->for derivative matching
2230 : ! b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
2231 : ! bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
2232 :
2233 : ! b(x)=x**3-x**2 b'(x)=3*x**2-2*x
2234 : ! bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
2235 :
2236 0 : real(dp) sum
2237 : integer v
2238 :
2239 0 : do v=1,ivec
2240 0 : i=ii(v)
2241 0 : j=jj(v)
2242 :
2243 : ! ...in x direction
2244 :
2245 0 : xp=xparam(v)
2246 0 : xpi=1.d0-xp
2247 0 : xp2=xp*xp
2248 0 : xpi2=xpi*xpi
2249 0 : ax=xp2*(3.d0-2.d0*xp)
2250 0 : axbar=1.d0-ax
2251 0 : bx=-xp2*xpi
2252 0 : bxbar=xpi2*xp
2253 :
2254 : ! ...in y direction
2255 :
2256 0 : yp=yparam(v)
2257 0 : ypi=1.d0-yp
2258 0 : yp2=yp*yp
2259 0 : ypi2=ypi*ypi
2260 0 : ay=yp2*(3.d0-2.d0*yp)
2261 0 : aybar=1.d0-ay
2262 0 : by=-yp2*ypi
2263 0 : bybar=ypi2*yp
2264 :
2265 : ! ...derivatives...
2266 :
2267 0 : axp=6.d0*xp*xpi
2268 0 : axbarp=-axp
2269 0 : bxp=xp*(3.d0*xp-2.d0)
2270 0 : bxbarp=xpi*(3.d0*xpi-2.d0)
2271 :
2272 0 : ayp=6.d0*yp*ypi
2273 0 : aybarp=-ayp
2274 0 : byp=yp*(3.d0*yp-2.d0)
2275 0 : bybarp=ypi*(3.d0*ypi-2.d0)
2276 :
2277 0 : iadr=0
2278 :
2279 : ! get desired values:
2280 :
2281 0 : if(ict(1) == 1) then
2282 :
2283 : ! function value:
2284 :
2285 0 : iadr=iadr+1
2286 : sum=axbar*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
2287 0 : & ax*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1))
2288 :
2289 : sum=sum+hx(v)*(
2290 : & bxbar*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
2291 0 : & bx*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1)))
2292 :
2293 : sum=sum+hy(v)*(
2294 : & axbar*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
2295 0 : & ax*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
2296 :
2297 : sum=sum+hx(v)*hy(v)*(
2298 : & bxbar*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
2299 0 : & bx*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
2300 :
2301 0 : fval(v,iadr)=sum
2302 : end if
2303 :
2304 0 : if(ict(2) == 1) then
2305 :
2306 : ! df/dx:
2307 :
2308 0 : iadr=iadr+1
2309 :
2310 : sum=hxi(v)*(
2311 : & axbarp*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
2312 0 : & axp*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1)))
2313 :
2314 : sum=sum+
2315 : & bxbarp*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
2316 0 : & bxp*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1))
2317 :
2318 : sum=sum+hxi(v)*hy(v)*(
2319 : & axbarp*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
2320 0 : & axp*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
2321 :
2322 : sum=sum+hy(v)*(
2323 : & bxbarp*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
2324 0 : & bxp*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
2325 :
2326 0 : fval(v,iadr)=sum
2327 : end if
2328 :
2329 0 : if(ict(3) == 1) then
2330 :
2331 : ! df/dy:
2332 :
2333 0 : iadr=iadr+1
2334 :
2335 : sum=hyi(v)*(
2336 : & axbar*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
2337 0 : & ax*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
2338 :
2339 : sum=sum+hx(v)*hyi(v)*(
2340 : & bxbar*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
2341 0 : & bx*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
2342 :
2343 : sum=sum+
2344 : & axbar*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
2345 0 : & ax*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1))
2346 :
2347 : sum=sum+hx(v)*(
2348 : & bxbar*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
2349 0 : & bx*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1)))
2350 :
2351 0 : fval(v,iadr)=sum
2352 : end if
2353 :
2354 0 : if(ict(4) == 1) then
2355 :
2356 : ! d2f/dxdy:
2357 :
2358 0 : iadr=iadr+1
2359 :
2360 : sum=hxi(v)*hyi(v)*(
2361 : & axbarp*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
2362 0 : & axp*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
2363 :
2364 : sum=sum+hyi(v)*(
2365 : & bxbarp*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
2366 0 : & bxp*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
2367 :
2368 : sum=sum+hxi(v)*(
2369 : & axbarp*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
2370 0 : & axp*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1)))
2371 :
2372 : sum=sum+
2373 : & bxbarp*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
2374 0 : & bxp*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1))
2375 :
2376 0 : fval(v,iadr)=sum
2377 : end if
2378 :
2379 : end do ! vector loop
2380 :
2381 0 : return
2382 : end subroutine herm2fcn_db
2383 :
2384 :
2385 : !---------------------------------------------------------------------
2386 : ! evaluate C1 cubic Hermite function interpolation -- 2d fcn
2387 : ! --vectorized-- dmc 10 Feb 1999
2388 :
2389 0 : subroutine herm2fcn_mesa_db(ict,ivec,ivecd,
2390 0 : & fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
2391 : & f1,inf2,ny)
2392 :
2393 : integer ny, inf2
2394 : integer ict(4) ! requested output control
2395 : integer ivec ! vector length
2396 : integer ivecd ! vector dimension (1st dim of fval)
2397 :
2398 : integer ii(ivec),jj(ivec) ! target cells (i,j)
2399 : real(dp) :: xparam(ivec),yparam(ivec)
2400 : ! normalized displacements from (i,j) corners
2401 :
2402 : real(dp) :: hx(ivec),hy(ivec) ! grid spacing, and
2403 : real(dp) :: hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
2404 : ! & 1/(y(j+1)-y(j))
2405 :
2406 : real(dp), pointer :: f1(:) ! =(0:3,inf2,ny) interpolant data (cf "evbicub")
2407 :
2408 : real(dp) :: fval(ivecd,4) ! output returned
2409 :
2410 : ! for detailed description of fin, ict and fval see subroutine
2411 : ! herm2ev comments. Note ict is not vectorized; the same output
2412 : ! is expected to be returned for all input vector data points.
2413 :
2414 : ! note that the index inputs ii,jj and parameter inputs
2415 : ! xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
2416 : ! output array fval has a vector ** 1st dimension ** whose
2417 : ! size must be given as a separate argument
2418 :
2419 : ! to use this routine in scalar mode, pass in ivec=ivecd=1
2420 :
2421 : !---------------
2422 : ! Hermite cubic basis functions
2423 : ! -->for function value matching
2424 : ! a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
2425 : ! abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
2426 :
2427 : ! a(x)=-2*x**3 + 3*x**2 = x*x*(-2.d0*x+3.0)
2428 : ! abar(x)=1-a(x)
2429 : ! a'(x)=-abar'(x) = 6.d0*x*(1.0-x)
2430 :
2431 : ! -->for derivative matching
2432 : ! b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
2433 : ! bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
2434 :
2435 : ! b(x)=x**3-x**2 b'(x)=3*x**2-2*x
2436 : ! bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
2437 :
2438 0 : real(dp) :: sum
2439 : integer v,iadr,i,j
2440 0 : real(dp) :: xp,yp,xpi,ypi,xp2,yp2,xpi2,ypi2
2441 0 : real(dp) :: ax,bx,axbar,bxbar,ay,by,aybar,bybar
2442 0 : real(dp) :: axp,axbarp,bxp,bxbarp,ayp,aybarp,bybarp,byp
2443 0 : real(dp), pointer :: fin(:,:,:)
2444 :
2445 0 : fin(0:3,1:inf2,1:ny) => f1(1:4*inf2*ny)
2446 :
2447 0 : do v=1,ivec
2448 0 : i=ii(v)
2449 0 : j=jj(v)
2450 :
2451 : ! ...in x direction
2452 :
2453 0 : xp=xparam(v)
2454 0 : xpi=1.d0-xp
2455 0 : xp2=xp*xp
2456 0 : xpi2=xpi*xpi
2457 0 : ax=xp2*(3.d0-2.d0*xp)
2458 0 : axbar=1.d0-ax
2459 0 : bx=-xp2*xpi
2460 0 : bxbar=xpi2*xp
2461 :
2462 : ! ...in y direction
2463 :
2464 0 : yp=yparam(v)
2465 0 : ypi=1.d0-yp
2466 0 : yp2=yp*yp
2467 0 : ypi2=ypi*ypi
2468 0 : ay=yp2*(3.d0-2.d0*yp)
2469 0 : aybar=1.d0-ay
2470 0 : by=-yp2*ypi
2471 0 : bybar=ypi2*yp
2472 :
2473 : ! ...derivatives...
2474 :
2475 0 : axp=6.d0*xp*xpi
2476 0 : axbarp=-axp
2477 0 : bxp=xp*(3.d0*xp-2.d0)
2478 0 : bxbarp=xpi*(3.d0*xpi-2.d0)
2479 :
2480 0 : ayp=6.d0*yp*ypi
2481 0 : aybarp=-ayp
2482 0 : byp=yp*(3.d0*yp-2.d0)
2483 0 : bybarp=ypi*(3.d0*ypi-2.d0)
2484 :
2485 0 : iadr=0
2486 :
2487 : ! get desired values:
2488 :
2489 0 : if(ict(1) == 1) then
2490 :
2491 : ! function value:
2492 :
2493 0 : iadr=iadr+1
2494 : sum=axbar*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
2495 0 : & ax*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1))
2496 :
2497 : sum=sum+hx(v)*(
2498 : & bxbar*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
2499 0 : & bx*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1)))
2500 :
2501 : sum=sum+hy(v)*(
2502 : & axbar*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
2503 0 : & ax*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
2504 :
2505 : sum=sum+hx(v)*hy(v)*(
2506 : & bxbar*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
2507 0 : & bx*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
2508 :
2509 0 : fval(v,iadr)=sum
2510 : end if
2511 :
2512 0 : if(ict(2) == 1) then
2513 :
2514 : ! df/dx:
2515 :
2516 0 : iadr=iadr+1
2517 :
2518 : sum=hxi(v)*(
2519 : & axbarp*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
2520 0 : & axp*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1)))
2521 :
2522 : sum=sum+
2523 : & bxbarp*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
2524 0 : & bxp*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1))
2525 :
2526 : sum=sum+hxi(v)*hy(v)*(
2527 : & axbarp*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
2528 0 : & axp*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
2529 :
2530 : sum=sum+hy(v)*(
2531 : & bxbarp*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
2532 0 : & bxp*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
2533 :
2534 0 : fval(v,iadr)=sum
2535 : end if
2536 :
2537 0 : if(ict(3) == 1) then
2538 :
2539 : ! df/dy:
2540 :
2541 0 : iadr=iadr+1
2542 :
2543 : sum=hyi(v)*(
2544 : & axbar*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
2545 0 : & ax*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
2546 :
2547 : sum=sum+hx(v)*hyi(v)*(
2548 : & bxbar*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
2549 0 : & bx*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
2550 :
2551 : sum=sum+
2552 : & axbar*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
2553 0 : & ax*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1))
2554 :
2555 : sum=sum+hx(v)*(
2556 : & bxbar*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
2557 0 : & bx*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1)))
2558 :
2559 0 : fval(v,iadr)=sum
2560 : end if
2561 :
2562 0 : if(ict(4) == 1) then
2563 :
2564 : ! d2f/dxdy:
2565 :
2566 0 : iadr=iadr+1
2567 :
2568 : sum=hxi(v)*hyi(v)*(
2569 : & axbarp*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
2570 0 : & axp*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
2571 :
2572 : sum=sum+hyi(v)*(
2573 : & bxbarp*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
2574 0 : & bxp*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
2575 :
2576 : sum=sum+hxi(v)*(
2577 : & axbarp*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
2578 0 : & axp*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1)))
2579 :
2580 : sum=sum+
2581 : & bxbarp*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
2582 0 : & bxp*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1))
2583 :
2584 0 : fval(v,iadr)=sum
2585 : end if
2586 :
2587 : end do ! vector loop
2588 :
2589 0 : return
2590 0 : end subroutine herm2fcn_mesa_db
2591 :
2592 :
2593 13745940 : subroutine ibc_ck_db(ibc,slbl,xlbl,imin,imax,ier)
2594 :
2595 : ! check that spline routine ibc flag is in range
2596 :
2597 : implicit none
2598 : integer ibc ! input -- flag value
2599 : character*(*) slbl ! input -- subroutine name
2600 : character*(*) xlbl ! input -- axis label
2601 :
2602 : integer imin ! input -- min allowed value
2603 : integer imax ! input -- max allowed value
2604 :
2605 : integer ier ! output -- set =1 if error detected
2606 :
2607 : !----------------------
2608 :
2609 13745940 : if((ibc < imin).or.(ibc > imax)) then
2610 0 : ier=1
2611 : ! write(6,1001) slbl,xlbl,ibc,imin,imax
2612 : ! 1001 format(' ?',a,' -- ibc',a,' = ',i9,' out of range ',
2613 : ! & i2,' to ',i2)
2614 : end if
2615 :
2616 0 : return
2617 : end subroutine ibc_ck_db
2618 :
2619 :
2620 33168 : subroutine do_mkbicub_db(x,nx,y,ny,f1,nf2,
2621 33168 : & ibcxmin,bcxmin,ibcxmax,bcxmax,
2622 33168 : & ibcymin,bcymin,ibcymax,bcymax,
2623 : & ilinx,iliny,ier)
2624 :
2625 : ! setup bicubic spline, dynamic allocation of workspace
2626 : ! fortran-90 fixed form source
2627 :
2628 : ! --NOTE-- dmc 22 Feb 2004 -- rewrite for direct calculation of
2629 : ! coefficients, to avoid large transient use of memory.
2630 :
2631 : !
2632 : implicit NONE
2633 :
2634 : ! input:
2635 : integer nx ! length of x vector
2636 : integer ny ! length of y vector
2637 : real(dp) x(:) ! (nx) ! x vector, strict ascending
2638 : real(dp) y(:) ! (ny) ! y vector, strict ascending
2639 :
2640 : integer nf2 ! 2nd dimension of f, nf2 >= nx
2641 : ! input/output:
2642 : real(dp), pointer :: f1(:) ! =(4,nf2,ny) ! data & spline coefficients
2643 :
2644 : ! on input: f(1,i,j) = f(x(i),y(j))
2645 : ! on output: f(1,i,j) unchanged
2646 : ! f(2,i,j) = d2f/dx2(x(i),y(j))
2647 : ! f(3,i,j) = d2f/dy2(x(i),y(j))
2648 : ! f(4,i,j) = d4f/dx2dy2(x(i),y(j))
2649 :
2650 : ! and the interpolation formula for (x,y) in (x(i),x(i+1))^(y(j),y(j+1))
2651 : ! is:
2652 : ! hx = x(i+1)-x(i) hy = y(j+1)-y(j)
2653 : ! dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
2654 : ! dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
2655 : ! dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
2656 :
2657 : ! finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
2658 : ! +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
2659 : ! +1/6*hx**2*
2660 : ! dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
2661 : ! +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
2662 : ! +1/6*hy**2*
2663 : ! dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
2664 : ! +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
2665 : ! +1/36*hx**2*hy**2*
2666 : ! dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
2667 : ! +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
2668 :
2669 : ! where the f(2:4,*,*) are cleverly evaluated to assure
2670 : ! (a) finterp is continuous and twice differentiable across all
2671 : ! grid cell boundaries, and
2672 : ! (b) all boundary conditions are satisfied.
2673 :
2674 : ! input bdy condition data:
2675 : integer ibcxmin ! bc flag for x=xmin
2676 : real(dp) bcxmin(:) ! (ny) ! bc data vs. y at x=xmin
2677 : integer ibcxmax ! bc flag for x=xmax
2678 : real(dp) bcxmax(:) ! (ny) ! bc data vs. y at x=xmax
2679 :
2680 : integer ibcymin ! bc flag for y=ymin
2681 : real(dp) bcymin(:) ! (nx) ! bc data vs. x at y=ymin
2682 : integer ibcymax ! bc flag for y=ymax
2683 : real(dp) bcymax(:) ! (nx) ! bc data vs. x at y=ymax
2684 :
2685 : ! with interpretation:
2686 : ! ibcxmin -- indicator for boundary condition at x(1):
2687 : ! bcxmin(...) -- boundary condition data
2688 : ! =-1 -- periodic boundary condition
2689 : ! =0 -- use "not a knot"
2690 : ! =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
2691 : ! =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
2692 : ! =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
2693 : ! =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
2694 : ! =5 -- match 1st derivative to 1st divided difference
2695 : ! =6 -- match 2nd derivative to 2nd divided difference
2696 : ! =7 -- match 3rd derivative to 3rd divided difference
2697 : ! (for more detailed definition of BCs 5-7, see the
2698 : ! comments of subroutine mkspline)
2699 : ! NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
2700 :
2701 : ! ibcxmax -- indicator for boundary condition at x(nx):
2702 : ! bcxmax(...) -- boundary condition data
2703 : ! (interpretation as with ibcxmin, bcxmin)
2704 : ! NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
2705 :
2706 : ! and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
2707 : ! (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
2708 :
2709 : ! output linear grid flags and completion code (ier=0 is normal):
2710 :
2711 : integer ilinx ! =1: x grid is "nearly" equally spaced
2712 : integer iliny ! =1: y grid is "nearly" equally spaced
2713 : ! ilinx and iliny are set to zero if corresponding grids are not equally
2714 : ! spaced
2715 :
2716 : integer ier ! =0 on exit if there is no error.
2717 :
2718 : ! if there is an error, ier is set and a message is output on unit 6.
2719 : ! these are considered programming errors in the calling routine.
2720 :
2721 : ! possible errors:
2722 : ! x(...) not strict ascending
2723 : ! y(...) not strict ascending
2724 : ! nx < 4
2725 : ! ny < 4
2726 : ! invalid boundary condition flag
2727 :
2728 : !-----------------------
2729 : integer ierx,iery
2730 :
2731 33168 : real(dp), dimension(:,:), allocatable :: fwk
2732 33168 : real(dp) :: zbcmin,zbcmax
2733 : integer ix,iy,ibcmin,ibcmax
2734 :
2735 33168 : real(dp), dimension(:,:,:), allocatable :: fcorr
2736 : integer iflg2
2737 99504 : real(dp) zdiff(2),hy
2738 :
2739 33168 : real(dp), pointer :: f(:,:,:) ! =(4,nf2,ny) ! data & spline coefficients
2740 33168 : f(1:4,1:nf2,1:ny) => f1(1:4*nf2*ny)
2741 :
2742 :
2743 : !-----------------------
2744 :
2745 : ! see if 2nd pass is needed due to inhomogeneous d/dy bdy cond.
2746 :
2747 33168 : iflg2=0
2748 33168 : if(ibcymin /= -1) then
2749 33168 : if((ibcymin == 1).or.(ibcymin == 2)) then
2750 0 : do ix=1,nx
2751 0 : if (bcymin(ix) /= 0.d0) iflg2=1
2752 : end do
2753 : end if
2754 33168 : if((ibcymax == 1).or.(ibcymax == 2)) then
2755 0 : do ix=1,nx
2756 0 : if (bcymax(ix) /= 0.d0) iflg2=1
2757 : end do
2758 : end if
2759 : end if
2760 :
2761 : ! check boundary condition specifications
2762 :
2763 33168 : ier=0
2764 :
2765 33168 : call ibc_ck_db(ibcxmin,'bcspline','xmin',-1,7,ier)
2766 66336 : if(ibcxmin >= 0) call ibc_ck_db(ibcxmax,'bcspline','xmax',0,7,ier)
2767 33168 : call ibc_ck_db(ibcymin,'bcspline','ymin',-1,7,ier)
2768 33168 : if(ibcymin >= 0) call ibc_ck_db(ibcymax,'bcspline','ymax',0,7,ier)
2769 :
2770 : ! check ilinx & x vector
2771 :
2772 33168 : call splinck_db(x,nx,ilinx,1.0d-3,ierx)
2773 33168 : if(ierx /= 0) ier=2
2774 :
2775 : ! if(ier == 2) then
2776 : ! write(6,'('' ?bcspline: x axis not strict ascending'')')
2777 : ! end if
2778 :
2779 : ! check iliny & y vector
2780 :
2781 33168 : call splinck_db(y,ny,iliny,1.0d-3,iery)
2782 33168 : if(iery /= 0) ier=3
2783 :
2784 : ! if(ier == 3) then
2785 : ! write(6,'('' ?bcspline: y axis not strict ascending'')')
2786 : ! end if
2787 :
2788 33168 : if(ier /= 0) return
2789 :
2790 : !------------------------------------
2791 33168 : allocate(fwk(2,max(nx,ny)))
2792 :
2793 : ! evaluate fxx (spline in x direction)
2794 :
2795 33168 : zbcmin=0
2796 33168 : zbcmax=0
2797 2287242 : do iy=1,ny
2798 527688500 : fwk(1,1:nx) = f(1,1:nx,iy)
2799 2254074 : if((ibcxmin == 1).or.(ibcxmin == 2)) zbcmin=bcxmin(iy)
2800 2254074 : if((ibcxmax == 1).or.(ibcxmax == 2)) zbcmax=bcxmax(iy)
2801 2254074 : call mkspline_db(x,nx,fwk,ibcxmin,zbcmin,ibcxmax,zbcmax,ilinx,ier)
2802 2254074 : if(ier /= 0) then
2803 0 : deallocate(fwk)
2804 0 : return
2805 : end if
2806 527721668 : f(2,1:nx,iy)=fwk(2,1:nx)
2807 : end do
2808 :
2809 : ! evaluate fyy (spline in y direction)
2810 : ! use homogeneous boundary condition; correction done later if necessary
2811 :
2812 33168 : zbcmin=0
2813 33168 : zbcmax=0
2814 33168 : ibcmin=ibcymin
2815 33168 : ibcmax=ibcymax
2816 2309448 : do ix=1,nx
2817 527710706 : fwk(1,1:ny) = f(1,ix,1:ny)
2818 2276280 : if(iflg2 == 1) then
2819 0 : if((ibcymin == 1).or.(ibcymin == 2)) ibcmin=0
2820 0 : if((ibcymax == 1).or.(ibcymax == 2)) ibcmax=0
2821 : end if
2822 2276280 : call mkspline_db(y,ny,fwk,ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
2823 2276280 : if(ier /= 0) then
2824 0 : deallocate(fwk)
2825 0 : return
2826 : end if
2827 527743874 : f(3,ix,1:ny)=fwk(2,1:ny)
2828 : end do
2829 :
2830 : ! evaluate fxxyy (spline fxx in y direction; BC simplified; avg
2831 : ! d2(d2f/dx2)/dy2 and d2(df2/dy2)/dx2
2832 :
2833 : zbcmin=0
2834 : zbcmax=0
2835 33168 : ibcmin=ibcymin
2836 33168 : ibcmax=ibcymax
2837 2309448 : do ix=1,nx
2838 527710706 : fwk(1,1:ny) = f(2,ix,1:ny)
2839 2276280 : if(iflg2 == 1) then
2840 0 : if((ibcymin == 1).or.(ibcymin == 2)) ibcmin=0
2841 0 : if((ibcymax == 1).or.(ibcymax == 2)) ibcmax=0
2842 : end if
2843 : call mkspline_db(y,ny,fwk,
2844 2276280 : & ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
2845 2276280 : if(ier /= 0) then
2846 0 : deallocate(fwk)
2847 0 : return
2848 : end if
2849 527743874 : f(4,ix,1:ny)= fwk(2,1:ny)
2850 : end do
2851 :
2852 33168 : if(iflg2 == 1) then
2853 0 : allocate(fcorr(2,nx,ny))
2854 :
2855 : ! correct for inhomogeneous y boundary condition
2856 :
2857 0 : do ix=1,nx
2858 : ! the desired inhomogeneouss BC is the difference btw the
2859 : ! requested derivative (1st or 2nd) and the current value
2860 :
2861 0 : zdiff(1)=0.d0
2862 0 : if(ibcymin == 1) then
2863 0 : hy=y(2)-y(1)
2864 : zdiff(1)=(f(1,ix,2)-f(1,ix,1))/hy +
2865 0 : & hy*(-2*f(3,ix,1)-f(3,ix,2))/6
2866 0 : zdiff(1)=bcymin(ix)-zdiff(1)
2867 0 : else if(ibcymin == 2) then
2868 0 : zdiff(1)=bcymin(ix)-f(3,ix,1)
2869 : end if
2870 :
2871 0 : zdiff(2)=0.d0
2872 0 : if(ibcymax == 1) then
2873 0 : hy=y(ny)-y(ny-1)
2874 : zdiff(2)=(f(1,ix,ny)-f(1,ix,ny-1))/hy +
2875 0 : & hy*(2*f(3,ix,ny)+f(3,ix,ny-1))/6
2876 0 : zdiff(2)=bcymax(ix)-zdiff(2)
2877 0 : else if(ibcymax == 2) then
2878 0 : zdiff(2)=bcymax(ix)-f(3,ix,ny)
2879 : end if
2880 :
2881 0 : fwk(1,1:ny)=0.d0 ! values are zero; only BC is not
2882 : call mkspline_db(y,ny,fwk,ibcymin,zdiff(1),ibcymax,zdiff(2),
2883 0 : & iliny,ier)
2884 0 : if(ier /= 0) then
2885 0 : deallocate(fwk,fcorr)
2886 0 : return
2887 : end if
2888 0 : fcorr(1,ix,1:ny)=fwk(2,1:ny) ! fyy-correction
2889 : end do
2890 :
2891 : zbcmin=0
2892 : zbcmax=0
2893 0 : do iy=1,ny
2894 0 : fwk(1,1:nx)=fcorr(1,1:nx,iy)
2895 0 : call mkspline_db(x,nx,fwk,ibcxmin,zbcmin,ibcxmax,zbcmax,ilinx,ier)
2896 0 : if(ier /= 0) then
2897 0 : deallocate(fwk,fcorr)
2898 0 : return
2899 : end if
2900 0 : fcorr(2,1:nx,iy)=fwk(2,1:nx) ! fxxyy-correction
2901 : end do
2902 :
2903 0 : f(3:4,1:nx,1:ny)=f(3:4,1:nx,1:ny)+fcorr(1:2,1:nx,1:ny)
2904 :
2905 0 : deallocate(fcorr)
2906 : end if
2907 :
2908 : ! correction spline -- f=fxx=zero; fyy & fxxyy are affected
2909 :
2910 33168 : deallocate(fwk)
2911 : !------------------------------------
2912 :
2913 : ! that's all
2914 :
2915 33168 : return
2916 33168 : end subroutine do_mkbicub_db
2917 :
2918 6806634 : subroutine mkspline_db(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
2919 : & ilinx,ier)
2920 : implicit none
2921 :
2922 : ! make a 2-coefficient 1d spline
2923 :
2924 : ! only 2 coefficients, the data and its 2nd derivative, are needed to
2925 : ! fully specify a spline. See e.g. Numerical Recipes in Fortran-77
2926 : ! (2nd edition) chapter 3, section on cubic splines.
2927 :
2928 : ! input:
2929 : integer nx ! no. of data points
2930 : real(dp) x(nx) ! x axis data, strict ascending order
2931 :
2932 : ! input/output:
2933 : real(dp) fspl(2,nx) ! f(1,*): data in; f(2,*): coeffs out
2934 :
2935 : ! f(1,j) = f(x(j)) on input (unchanged on output)
2936 : ! f(2,j) = f''(x(j)) (of interpolating spline) (on output).
2937 :
2938 : ! ...boundary conditions...
2939 :
2940 : ! input:
2941 :
2942 : integer ibcxmin ! b.c. flag @ x=xmin=x(1)
2943 : real(dp) bcxmin ! b.c. data @xmin
2944 :
2945 : integer ibcxmax ! b.c. flag @ x=xmax=x(nx)
2946 : real(dp) bcxmax ! b.c. data @xmax
2947 :
2948 : ! ibcxmin=-1 -- periodic boundary condition
2949 : ! (bcxmin,ibcxmax,bcxmax are ignored)
2950 :
2951 : ! the output spline s satisfies
2952 : ! s'(x(1))=s'(x(nx)) ..and.. s''(x(1))=s''(x(nx))
2953 :
2954 : ! if non-periodic boundary conditions are used, then the xmin and xmax
2955 : ! boundary conditions can be specified independently:
2956 :
2957 : ! ibcxmin (ibcxmax) = 0 -- this specifies a "not a knot" boundary
2958 : ! condition, see "cubsplb.for". This is a common way
2959 : ! for inferring a "good" spline boundary condition
2960 : ! automatically from data in the vicinity of the
2961 : ! boundary. ... bcxmin (bcxmax) are ignored.
2962 :
2963 : ! ibcxmin (ibcxmax) = 1 -- boundary condition is to have s'(x(1))
2964 : ! ( s'(x(nx)) ) match the passed value bcxmin (bcxmax).
2965 :
2966 : ! ibcxmin (ibcxmax) = 2 -- boundary condition is to have s''(x(1))
2967 : ! ( s''(x(nx)) ) match the passed value bcxmin (bcxmax).
2968 :
2969 : ! ibcxmin (ibcxmax) = 3 -- boundary condition is to have s'(x(1))=0.d0
2970 : ! ( s'(x(nx))=0.d0 )
2971 :
2972 : ! ibcxmin (ibcxmax) = 4 -- boundary condition is to have s''(x(1))=0.d0
2973 : ! ( s''(x(nx))=0.d0 )
2974 :
2975 : ! ibcxmin (ibcxmax) = 5 -- boundary condition is to have s'(x(1))
2976 : ! ( s'(x(nx)) ) match the 1st divided difference
2977 : ! e.g. at x(1): d(1)/h(1), where
2978 : ! d(j)=f(1,j+1)-f(1,j)
2979 : ! h(j)=x(j+1)-x(j)
2980 :
2981 : ! ibcxmin (ibcxmax) = 6 -- BC is to have s''(x(1)) ( s''(x(nx)) )
2982 : ! match the 2nd divided difference
2983 : ! e.g. at x(1):
2984 : ! e(1) = [d(2)/h(2) - d(1)/h(1)]/(0.5*(h(1)+h(2)))
2985 :
2986 : ! ibcxmin (ibcxmax) = 7 -- BC is to have s'''(x(1)) ( s'''(x(nx)) )
2987 : ! match the 3rd divided difference
2988 : ! e.g. at x(1): [e(2)-e(1)]/(0.33333*(h(1)+h(2)+h(3)))
2989 :
2990 : ! output:
2991 :
2992 : integer ilinx ! =1: hint, x axis is ~evenly spaced
2993 :
2994 : ! let dx[avg] = (x(nx)-x(1))/(nx-1)
2995 : ! let dx[j] = x(j+1)-x(j), for all j satisfying 1 <= j < nx
2996 :
2997 : ! if for all such j, abs(dx[j]-dx[avg]) <= (1.0d-3*dx[avg]) then
2998 : ! ilinx=1 is returned, indicating the data is (at least nearly)
2999 : ! evenly spaced. Even spacing is useful, for speed of zone lookup,
3000 : ! when evaluating a spline.
3001 :
3002 : ! if the even spacing condition is not satisfied, ilinx=2 is returned.
3003 :
3004 : integer ier ! exit code, 0=OK
3005 :
3006 : ! an error code is returned if the x axis is not strict ascending,
3007 : ! or if nx < 4, or if an invalid boundary condition specification was
3008 : ! input.
3009 :
3010 : !------------------------------------
3011 :
3012 : ! this routine calls traditional 4 coefficient spline software, and
3013 : ! translates the result to 2 coefficient form.
3014 :
3015 : ! this could be done more efficiently but we decided out of conservatism
3016 : ! to use the traditional software.
3017 :
3018 : !------------------------------------
3019 : integer i, inwk
3020 : ! workspaces -- f90 dynamically allocated
3021 :
3022 6806634 : real(dp), dimension(:,:), allocatable :: fspl4 ! traditional 4-spline
3023 6806634 : real(dp), dimension(:), allocatable :: wk ! cspline workspace
3024 :
3025 : !------------------------------------
3026 :
3027 6806634 : allocate(fspl4(4,nx),wk(nx))
3028 :
3029 : ! make the traditional call
3030 :
3031 1583109912 : do i=1,nx
3032 1576303278 : fspl4(1,i)=fspl(1,i)
3033 1583109912 : fspl(2,i)=0.d0 ! for now
3034 : end do
3035 :
3036 6806634 : inwk=nx
3037 :
3038 : ! boundary conditions imposed by cspline...
3039 :
3040 : call cspline_db(x,nx,fspl4,ibcxmin,bcxmin,ibcxmax,bcxmax,
3041 6806634 : & wk,inwk,ilinx,ier)
3042 :
3043 6806634 : if(ier == 0) then
3044 :
3045 : ! copy the output -- careful of end point.
3046 :
3047 1576303278 : do i=1,nx-1
3048 1576303278 : fspl(2,i) = 2.d0*fspl4(3,i)
3049 : end do
3050 6806634 : fspl(2,nx) = 2.d0*fspl4(3,nx-1) + (x(nx)-x(nx-1))*6.d0*fspl4(4,nx-1)
3051 : end if
3052 :
3053 6806634 : deallocate(fspl4,wk)
3054 :
3055 6806634 : return
3056 : end subroutine mkspline_db
3057 :
3058 :
3059 6872970 : subroutine splinck_db(x,inx,ilinx,ztol,ier)
3060 : implicit none
3061 :
3062 : ! check if a grid is strictly ascending and if it is evenly spaced
3063 : ! to w/in ztol
3064 :
3065 : integer inx
3066 : real(dp) x(inx) ! input -- grid to check
3067 :
3068 : integer ilinx ! output -- =1 if evenly spaced =2 O.W.
3069 :
3070 : real(dp) ztol ! input -- spacing check tolerance
3071 :
3072 : integer ier ! output -- =0 if OK
3073 :
3074 : ! ier=1 is returned if x(1...inx) is NOT STRICTLY ASCENDING...
3075 :
3076 : !-------------------------------
3077 :
3078 6872970 : real(dp) zeps, zdiffx, zdiff, dxavg
3079 : integer ix
3080 :
3081 6872970 : ier=0
3082 6872970 : ilinx=1
3083 6872970 : if(inx <= 1) return
3084 :
3085 6872970 : dxavg=(x(inx)-x(1))/(inx-1)
3086 6872970 : zeps=abs(ztol*dxavg)
3087 :
3088 1580833632 : do ix=2,inx
3089 1573960662 : zdiffx=(x(ix)-x(ix-1))
3090 1573960662 : if(zdiffx <= 0.0) ier=2
3091 1573960662 : zdiff=zdiffx-dxavg
3092 1580833632 : if(abs(zdiff) > zeps) then
3093 48131078 : ilinx=2
3094 : end if
3095 : end do
3096 : 10 continue
3097 :
3098 : return
3099 : end subroutine splinck_db
3100 :
3101 6806634 : subroutine V_SPLINE_db(k_bc1,k_bcn,n,x,f,wk)
3102 : !***********************************************************************
3103 : !V_SPLINE evaluates the coefficients for a 1d cubic interpolating spline
3104 : !References:
3105 : ! Forsythe, Malcolm, Moler, Computer Methods for Mathematical
3106 : ! Computations, Prentice-Hall, 1977, p.76
3107 : ! Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
3108 : ! 1996, p.251
3109 : ! W.A.Houlberg, D.McCune 3/2000
3110 : !Input:
3111 : ! k_bc1-option for BC at x(1)
3112 : ! =-1 periodic, ignore k_bcn
3113 : ! =0 not-a-knot
3114 : ! =1 s'(x1) = input value of f(2,1)
3115 : ! =2 s''(x1) = input value of f(3,1)
3116 : ! =3 s'(x1) = 0.0
3117 : ! =4 s''(x1) = 0.0
3118 : ! =5 match first derivative to first 2 points
3119 : ! =6 match second derivative to first 3 points
3120 : ! =7 match third derivative to first 4 points
3121 : ! =else use not-a-knot
3122 : ! k_bcn-option for boundary condition at x(n)
3123 : ! =0 not-a-knot
3124 : ! =1 s'(x1) = input value of f(2,1)
3125 : ! =2 s''(x1) = input value of f(3,1)
3126 : ! =3 s'(x1) = 0.0
3127 : ! =4 s''(x1) = 0.0
3128 : ! =5 match first derivative to first 2 points
3129 : ! =6 match second derivative to first 3 points
3130 : ! =7 match third derivative to first 4 points
3131 : ! =else use knot-a-knot
3132 : ! n-number of data points or knots-(n >= 2)
3133 : ! x(n)-abscissas of the knots in strictly increasing order
3134 : ! f(1,i)-ordinates of the knots
3135 : ! f(2,1)-input value of s'(x1) for k_bc1=1
3136 : ! f(2,n)-input value of s'(xn) for k_bcn=1
3137 : ! f(3,1)-input value of s''(x1) for k_bc1=2
3138 : ! f(3,n)-input value of s''(xn) for k_bcn=2
3139 : ! wk(n)-scratch work area for periodic BC
3140 : !Output:
3141 : ! f(2,i)=s'(x(i))
3142 : ! f(3,i)=s''(x(i))
3143 : ! f(4,i)=s'''(x(i))
3144 : !Comments:
3145 : ! s(x)=f(1,i)+f(2,i)*(x-x(i))+f(3,i)*(x-x(i))**2/2!
3146 : ! +f(4,i)*(x-x(i))**3/3! for x(i) <= x <= x(i+1)
3147 : ! W_SPLINE can be used to evaluate the spline and its derivatives
3148 : ! The cubic spline is twice differentiable (C2)
3149 :
3150 : ! bugfixes -- dmc 24 Feb 2004:
3151 : ! (a) fixed logic for not-a-knot:
3152 : ! ! Set f(3,1) for not-a-knot
3153 : ! if (k_bc1 <= 0.or.k_bc1 > 7) then ...
3154 : ! instead of
3155 : ! ! Set f(3,1) for not-a-knot
3156 : ! if (k_bc1 <= 0.or.k_bc1 > 5) then ...
3157 : ! and similarly for logic after cmt
3158 : ! ! Set f(3,n) for not-a-knot
3159 : ! as required since k_bc*=6 and k_bc*=7 are NOT not-a-knot BCs.
3160 :
3161 : ! (b) the BCs to fix 2nd derivative at end points did not work if that
3162 : ! 2nd derivative were non-zero. The reason is that in those cases
3163 : ! the off-diagonal matrix elements nearest the corners are not
3164 : ! symmetric; i.e. elem(1,2) /= elem(2,1) and
3165 : ! elem(n-1,n) /= elem(n,n-1) where I use "elem" to refer to
3166 : ! the tridiagonal matrix elements. The correct values for the
3167 : ! elements is: elem(1,2)=0, elem(2,1)=x(2)-x(1)
3168 : ! elem(n,n-1)=0, elem(n-1,n)=x(n)-x(n-1)
3169 : ! the old code in effect had these as all zeroes. Since this
3170 : ! meant the wrong set of linear equations was solved, the
3171 : ! resulting spline had a discontinuity in its 1st derivative
3172 : ! at x(2) and x(n-1). Fixed by introducing elem21 and elemnn1
3173 : ! to represent the non-symmetric lower-diagonal values. Since
3174 : ! elem21 & elemnn1 are both on the lower diagonals, logic to
3175 : ! use them occurs in the non-periodic forward elimination loop
3176 : ! only. DMC 24 Feb 2004.
3177 : !***********************************************************************
3178 : implicit none
3179 : !Declaration of input variables
3180 : integer :: k_bc1, k_bcn, n
3181 : real(dp) :: x(*), wk(*), f(4,*)
3182 : !Declaration in local variables
3183 : integer :: i, ib, imax, imin
3184 6806634 : real(dp) :: a1, an, b1, bn, q, t, hn
3185 6806634 : real(dp) :: elem21, elemnn1 ! (dmc)
3186 :
3187 : !Set default range
3188 6806634 : imin=1
3189 6806634 : imax=n
3190 : !Set first and second BC values
3191 6806634 : a1=0.d0
3192 6806634 : b1=0.d0
3193 6806634 : an=0.d0
3194 6806634 : bn=0.d0
3195 6806634 : if (k_bc1 == 1) then
3196 0 : a1=f(2,1)
3197 6806634 : else if(k_bc1 == 2) then
3198 0 : b1=f(3,1)
3199 6806634 : else if(k_bc1 == 5) then
3200 0 : a1=(f(1,2)-f(1,1))/(x(2)-x(1))
3201 6806634 : else if(k_bc1 == 6) then
3202 0 : b1=2.d0*((f(1,3)-f(1,2))/(x(3)-x(2)) -(f(1,2)-f(1,1))/(x(2)-x(1)))/(x(3)-x(1))
3203 : end if
3204 6806634 : if (k_bcn == 1) then
3205 0 : an=f(2,n)
3206 6806634 : else if(k_bcn == 2) then
3207 0 : bn=f(3,n)
3208 6806634 : else if(k_bcn == 5) then
3209 0 : an=(f(1,n)-f(1,n-1))/(x(n)-x(n-1))
3210 6806634 : else if(k_bcn == 6) then
3211 0 : bn=2.d0*((f(1,n)-f(1,n-1))/(x(n)-x(n-1)) -(f(1,n-1)-f(1,n-2))/(x(n-1)-x(n-2)))/(x(n)-x(n-2))
3212 : end if
3213 : !Clear f(2:4,n)
3214 6806634 : f(2,n)=0.d0
3215 6806634 : f(3,n)=0.d0
3216 6806634 : f(4,n)=0.d0
3217 6806634 : if (n == 2) then
3218 : !Coefficients for n=2
3219 0 : f(2,1)=(f(1,2)-f(1,1))/(x(2)-x(1))
3220 0 : f(3,1)=0.d0
3221 0 : f(4,1)=0.d0
3222 0 : f(2,2)=f(2,1)
3223 0 : f(3,2)=0.d0
3224 0 : f(4,2)=0.d0
3225 6806634 : else if(n > 2) then
3226 : !Set up tridiagonal system for A*y=B where y(i) are the second
3227 : ! derivatives at the knots
3228 : ! f(2,i) are the diagonal elements of A
3229 : ! f(4,i) are the off-diagonal elements of A
3230 : ! f(3,i) are the B elements/3, and will become c/3 upon solution
3231 6806634 : f(4,1)=x(2)-x(1)
3232 6806634 : f(3,2)=(f(1,2)-f(1,1))/f(4,1)
3233 1569496644 : do i=2,n-1
3234 1562690010 : f(4,i)=x(i+1)-x(i)
3235 1562690010 : f(2,i)=2.d0*(f(4,i-1)+f(4,i))
3236 1562690010 : f(3,i+1)=(f(1,i+1)-f(1,i))/f(4,i)
3237 1569496644 : f(3,i)=f(3,i+1)-f(3,i)
3238 : end do
3239 :
3240 : ! (dmc): save now:
3241 :
3242 6806634 : elem21=f(4,1)
3243 6806634 : elemnn1=f(4,n-1)
3244 :
3245 : ! BC's
3246 : ! Left
3247 : if (k_bc1 == -1) then
3248 0 : f(2,1)=2.d0*(f(4,1)+f(4,n-1))
3249 0 : f(3,1)=(f(1,2)-f(1,1))/f(4,1)-(f(1,n)-f(1,n-1))/f(4,n-1)
3250 0 : wk(1)=f(4,n-1)
3251 0 : do i=2,n-3
3252 0 : wk(i)=0.d0
3253 : end do
3254 0 : wk(n-2)=f(4,n-2)
3255 0 : wk(n-1)=f(4,n-1)
3256 : else if(k_bc1 == 1.or.k_bc1 == 3.or.k_bc1 == 5) then
3257 12825 : f(2,1)=2.d0*f(4,1)
3258 12825 : f(3,1)=(f(1,2)-f(1,1))/f(4,1)-a1
3259 : else if(k_bc1 == 2.or.k_bc1 == 4.or.k_bc1 == 6) then
3260 0 : f(2,1)=2.d0*f(4,1)
3261 0 : f(3,1)=f(4,1)*b1/3.d0
3262 0 : f(4,1)=0.d0 ! upper diagonal only (dmc: cf elem21)
3263 : else if(k_bc1 == 7) then
3264 0 : f(2,1)=-f(4,1)
3265 0 : f(3,1)=f(3,3)/(x(4)-x(2))-f(3,2)/(x(3)-x(1))
3266 0 : f(3,1)=f(3,1)*f(4,1)*f(4,1)/(x(4)-x(1))
3267 : else ! not a knot:
3268 6793809 : imin=2
3269 6793809 : f(2,2)=f(4,1)+2.d0*f(4,2)
3270 6793809 : f(3,2)=f(3,2)*f(4,2)/(f(4,1)+f(4,2))
3271 : end if
3272 : ! Right
3273 6806634 : if (k_bcn == 1.or.k_bcn == 3.or.k_bcn == 5) then
3274 12825 : f(2,n)=2.d0*f(4,n-1)
3275 12825 : f(3,n)=-(f(1,n)-f(1,n-1))/f(4,n-1)+an
3276 : else if(k_bcn == 2.or.k_bcn == 4.or.k_bcn == 6) then
3277 0 : f(2,n)=2.d0*f(4,n-1)
3278 0 : f(3,n)=f(4,n-1)*bn/3.d0
3279 : !xxx f(4,n-1)=0.d0 ! dmc: preserve f(4,n-1) for back subst.
3280 0 : elemnn1=0.d0 ! lower diagonal only (dmc)
3281 : else if(k_bcn == 7) then
3282 0 : f(2,n)=-f(4,n-1)
3283 0 : f(3,n)=f(3,n-1)/(x(n)-x(n-2))-f(3,n-2)/(x(n-1)-x(n-3))
3284 0 : f(3,n)=-f(3,n)*f(4,n-1)*f(4,n-1)/(x(n)-x(n-3))
3285 6793809 : else if(k_bc1 /= -1) then ! not a knot:
3286 6793809 : imax=n-1
3287 6793809 : f(2,n-1)=2.d0*f(4,n-2)+f(4,n-1)
3288 6793809 : f(3,n-1)=f(3,n-1)*f(4,n-2)/(f(4,n-1)+f(4,n-2))
3289 : end if
3290 : ! Limit solution for only three points in domain
3291 6806634 : if (n == 3) then
3292 0 : f(3,1)=0.d0
3293 0 : f(3,n)=0.d0
3294 : end if
3295 6806634 : if (k_bc1 == -1) then
3296 : !Solve system of equations for second derivatives at the knots
3297 : ! Periodic BC
3298 : ! Forward elimination
3299 0 : do i=2,n-2
3300 0 : t=f(4,i-1)/f(2,i-1)
3301 0 : f(2,i)=f(2,i)-t*f(4,i-1)
3302 0 : f(3,i)=f(3,i)-t*f(3,i-1)
3303 0 : wk(i)=wk(i)-t*wk(i-1)
3304 0 : q=wk(n-1)/f(2,i-1)
3305 0 : wk(n-1)=-q*f(4,i-1)
3306 0 : f(2,n-1)=f(2,n-1)-q*wk(i-1)
3307 0 : f(3,n-1)=f(3,n-1)-q*f(3,i-1)
3308 : end do
3309 : ! Correct the n-1 element
3310 0 : wk(n-1)=wk(n-1)+f(4,n-2)
3311 : ! Complete the forward elimination
3312 : ! wk(n-1) and wk(n-2) are the off-diag elements of the lower corner
3313 0 : t=wk(n-1)/f(2,n-2)
3314 0 : f(2,n-1)=f(2,n-1)-t*wk(n-2)
3315 0 : f(3,n-1)=f(3,n-1)-t*f(3,n-2)
3316 : ! Back substitution
3317 0 : f(3,n-1)=f(3,n-1)/f(2,n-1)
3318 0 : f(3,n-2)=(f(3,n-2)-wk(n-2)*f(3,n-1))/f(2,n-2)
3319 0 : do ib=3,n-1
3320 0 : i=n-ib
3321 0 : f(3,i)=(f(3,i)-f(4,i)*f(3,i+1)-wk(i)*f(3,n-1))/f(2,i)
3322 : end do
3323 0 : f(3,n)=f(3,1)
3324 : else
3325 : ! Non-periodic BC
3326 : ! Forward elimination
3327 : ! For Not-A-Knot BC the off-diagonal end elements are not equal
3328 1562715660 : do i=imin+1,imax
3329 1555909026 : if ((i == n-1).and.(imax == n-1)) then
3330 6793809 : t=(f(4,i-1)-f(4,i))/f(2,i-1)
3331 : else
3332 1549115217 : if(i == 2) then
3333 12825 : t=elem21/f(2,i-1)
3334 1549102392 : else if(i == n) then
3335 12825 : t=elemnn1/f(2,i-1)
3336 : else
3337 1549089567 : t=f(4,i-1)/f(2,i-1)
3338 : end if
3339 : end if
3340 1555909026 : if ((i == imin+1).and.(imin == 2)) then
3341 6793809 : f(2,i)=f(2,i)-t*(f(4,i-1)-f(4,i-2))
3342 : else
3343 1549115217 : f(2,i)=f(2,i)-t*f(4,i-1)
3344 : end if
3345 1562715660 : f(3,i)=f(3,i)-t*f(3,i-1)
3346 : end do
3347 : ! Back substitution
3348 6806634 : f(3,imax)=f(3,imax)/f(2,imax)
3349 1562715660 : do ib=1,imax-imin
3350 1555909026 : i=imax-ib
3351 1562715660 : if ((i == 2).and.(imin == 2)) then
3352 6793809 : f(3,i)=(f(3,i)-(f(4,i)-f(4,i-1))*f(3,i+1))/f(2,i)
3353 : else
3354 1549115217 : f(3,i)=(f(3,i)-f(4,i)*f(3,i+1))/f(2,i)
3355 : end if
3356 : end do
3357 : ! Reset d array to step size
3358 6806634 : f(4,1)=x(2)-x(1)
3359 6806634 : f(4,n-1)=x(n)-x(n-1)
3360 : ! Set f(3,1) for not-a-knot
3361 6806634 : if (k_bc1 <= 0.or.k_bc1 > 7) then
3362 6793809 : f(3,1)=(f(3,2)*(f(4,1)+f(4,2))-f(3,3)*f(4,1))/f(4,2)
3363 : end if
3364 : ! Set f(3,n) for not-a-knot
3365 6806634 : if (k_bcn <= 0.or.k_bcn > 7) then
3366 6793809 : f(3,n)=f(3,n-1)+(f(3,n-1)-f(3,n-2))*f(4,n-1)/f(4,n-2)
3367 : end if
3368 : end if
3369 : !f(3,i) is now the sigma(i) of the text and f(4,i) is the step size
3370 : !Compute polynomial coefficients
3371 1576303278 : do i=1,n-1
3372 1569496644 : f(2,i)= (f(1,i+1)-f(1,i))/f(4,i)-f(4,i)*(f(3,i+1)+2.d0*f(3,i))
3373 1569496644 : f(4,i)=(f(3,i+1)-f(3,i))/f(4,i)
3374 1569496644 : f(3,i)=6.d0*f(3,i)
3375 1576303278 : f(4,i)=6.d0*f(4,i)
3376 : end do
3377 6806634 : if (k_bc1 == -1) then
3378 0 : f(2,n)=f(2,1)
3379 0 : f(3,n)=f(3,1)
3380 0 : f(4,n)=f(4,1)
3381 : else
3382 6806634 : hn=x(n)-x(n-1)
3383 6806634 : f(2,n)=f(2,n-1)+hn*(f(3,n-1)+0.5d0*hn*f(4,n-1))
3384 6806634 : f(3,n)=f(3,n-1)+hn*f(4,n-1)
3385 6806634 : f(4,n)=f(4,n-1)
3386 : if (k_bcn == 1.or.k_bcn == 3.or.k_bcn == 5) then
3387 12825 : f(2,n)=an
3388 : else if(k_bcn == 2.or.k_bcn == 4.or.k_bcn == 6) then
3389 0 : f(3,n)=bn
3390 : end if
3391 : end if
3392 : end if
3393 6806634 : return
3394 : end subroutine V_SPLINE_db
3395 :
3396 :
3397 24 : subroutine zonfind_db(x,nx,zxget,i)
3398 : implicit none
3399 :
3400 : integer nx
3401 : real(dp) x(nx),zxget
3402 : integer i
3403 :
3404 : integer nxm, i1, i2, ii, ij
3405 24 : real(dp) dx
3406 :
3407 : ! find index i such that x(i) <= zxget <= x(i+1)
3408 :
3409 : ! x(1...nx) is strict increasing and x(1) <= zxget <= x(nx)
3410 : ! (this is assumed to already have been checked -- no check here!)
3411 :
3412 24 : nxm=nx-1
3413 24 : if((i < 1).or.(i > nxm)) then
3414 : i1=1
3415 : i2=nx-1
3416 : GOTO 10
3417 : end if
3418 :
3419 4 : if(x(i) > zxget) then
3420 : ! look down
3421 2 : dx=x(i+1)-x(i)
3422 2 : if((x(i)-zxget) > 4*dx) then
3423 2 : i1=1
3424 2 : i2=i-1
3425 2 : GOTO 10
3426 : else
3427 0 : i2=i-1
3428 0 : do ij=i2,1,-1
3429 0 : if((x(ij) <= zxget).and.(zxget <= x(ij+1))) then
3430 0 : i=ij
3431 0 : return
3432 : end if
3433 : end do
3434 0 : i=1
3435 0 : return
3436 : end if
3437 2 : else if(x(i+1) < zxget) then
3438 : ! look up
3439 2 : dx=x(i+1)-x(i)
3440 2 : if((zxget-x(i+1)) > 4*dx) then
3441 : i1=i+1
3442 : i2=nxm
3443 : GOTO 10
3444 : else
3445 0 : i2=i+1
3446 0 : do ij=i2,nxm
3447 0 : if((x(ij) <= zxget).and.(zxget <= x(ij+1))) then
3448 0 : i=ij
3449 0 : return
3450 : end if
3451 : end do
3452 24 : ij=nxm
3453 : return
3454 : end if
3455 : else
3456 : ! already there...
3457 : return
3458 : end if
3459 :
3460 : !---------------------------
3461 : ! binary search
3462 :
3463 : 10 continue
3464 :
3465 128 : if(i1 == i2) then
3466 : ! found by proc. of elimination
3467 10 : i=i1
3468 10 : return
3469 : end if
3470 :
3471 118 : ii=(i1+i2)/2
3472 :
3473 118 : if(zxget < x(ii)) then
3474 56 : i2=ii-1
3475 62 : else if(zxget > x(ii+1)) then
3476 : i1=ii+1
3477 : else
3478 : ! found
3479 14 : i=ii
3480 14 : return
3481 : end if
3482 :
3483 : GOTO 10
3484 :
3485 : end subroutine zonfind_db
3486 :
3487 : end module bicub_db
|
