V3FIT
r8evintrp3d.f
1  subroutine r8evintrp3d(xget,yget,zget,x,nx,y,ny,z,nz,jspline,
2  > f,icoeff,ixdim,iydim,izdim,ict,fval,ier)
3 C
4 C
5 C use mkintrp3d to set up spline coefficients...
6 C
7 C evaluate a 3d Hybrid interpolant on a rectilinear grid
8 C
9 C this subroutine calls two subroutines:
10 C vecin3d_argchk -- error check
11 C herm3xyz -- find cell containing (xget,yget,zget)
12 C fvintrp3d -- evaluate the spline function (w/derivatives if req.)
13 C
14 C input arguments:
15 C ================
16 C
17  IMPLICIT NONE
18  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
19  real*8 xget,yget,zget ! target of this interpolation
20  integer nx,ny,nz ! grid sizes
21  real*8 x(nx) ! ordered x grid
22  real*8 y(ny) ! ordered y grid
23  real*8 z(nz) ! ordered z grid
24 C
25  integer :: jspline(3) ! interpolation method for each
26 C dimension: jspline(1) for x; jspline(2) for y; jspline(3) for z
27 C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline
28 C
29  integer :: icoeff ! no. of coefficients per data point
30 
31  integer :: ixdim,iydim,izdim ! dimensioning:
32  ! ixdim=nx-1 for zonal step in x; otherwise nx
33  ! iydim=ny-1 for zonal step in y; otherwise ny
34  ! izdim=nz-1 for zonal step in z; otherwise nz
35 C
36  real*8 f(icoeff,ixdim,iydim,izdim) ! function data
37 C
38  integer ict(10) ! code specifying output desired
39 C
40 C Note on derivatives: for dimensions along which zonal step function
41 C interpolation is done, ANY derivative request returns ZERO.
42 C For dimensions along which piecewise linear or Hermite interpolation
43 C are done, more than one differentiation returns ZERO!
44 C
45 C Derivative controls are the same as for the compact 3d spline evaluation
46 C routine (evtricub):
47 C
48 C ict(1)=1 -- return f (0, don't)
49 C ict(2)=1 -- return df/dx (0, don't)
50 C ict(3)=1 -- return df/dy (0, don't)
51 C ict(4)=1 -- return df/dz (0, don't)
52 C ict(5)=1 -- return d2f/dx2 (0, don't)
53 C ict(6)=1 -- return d2f/dy2 (0, don't)
54 C ict(7)=1 -- return d2f/dz2 (0, don't)
55 C ict(8)=1 -- return d2f/dxdy (0, don't)
56 C ict(9)=1 -- return d2f/dxdz (0, don't)
57 C ict(10)=1-- return d2f/dydz (0, don't)
58 c
59 c (new dmc Dec 2005 -- higher derivatives available)
60 c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
61 c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
62 c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
63 c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
64 c ict(2:4) select: fxxx fyyy fzzz
65 c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
66 c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
67 c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
68 c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
69 c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
70 c ict(2:4) select: fxxyyz fxxyzz fxyyzz
71 c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
72 c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
73 c fxyyyz fxyzzz fyyyzz fzzzyy
74 c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
75 c fxxyyzz
76 c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
77 c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
78 c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
79 c ict(1)=-7 --> 7th derivative
80 c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
81 c fxxyyyzz fxxyyzzz fxyyyzzz
82 c ict(1)=-8 --> 8th derivative
83 c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
84 c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
85 c
86 C
87 C output arguments:
88 C =================
89 C
90  real*8 fval(10) ! output data
91  integer ier ! error code =0 ==> no error
92 C
93 C fval(1) receives the first output (depends on ict(...) spec)
94 C fval(2) receives the second output (depends on ict(...) spec)
95 C fval(3) receives the third output (depends on ict(...) spec)
96 C fval(4) receives the 4th output (depends on ict(...) spec)
97 C fval(5-10) receive 5th thru 10th outputs (if required by ict(...) spec)
98 C
99 C examples:
100 C on input ict = [1,1,1,1,0,0,0,0,0,0,0]
101 C on output fval= [f,df/dx,df/dy,df/dz]
102 C
103 C on input ict = [1,0,0,0,0,0,0,0,0,0,0]
104 C on output fval= [f] ... elements 2-10 never referenced
105 C
106 C on input ict = [0,1,1,0,0,0,0,0,0,0,0]
107 C on output fval= [df/dx,df/dy] ... elements 3-10 never referenced
108 C
109 C on input ict = [0,0,0,0,1,0,0,0,0,0,0]
110 C on output fval= [d2f/dx2] ... elements 2-10 never referenced.
111 C
112 C ier -- completion code: 0 means OK
113 C-------------------
114 C local:
115 C
116  integer i,j,k ! cell indices
117 C
118 C normalized displacement from (x(i),y(j),z(k)) corner of cell.
119 C xparam=0 @x(i) xparam=1 @x(i+1)
120 C yparam=0 @y(j) yparam=1 @y(j+1)
121 C zparam=0 @z(k) zparam=1 @z(k+1)
122 C
123  real*8 xparam,yparam,zparam
124 C
125 C cell dimensions and
126 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
127 C
128  real*8 hx,hy,hz
129  real*8 hxi,hyi,hzi
130 C
131 C 0 .le. xparam .le. 1
132 C 0 .le. yparam .le. 1
133 C 0 .le. zparam .le. 1
134 C
135 C---------------------------------------------------------------------
136 C
137  call vecin3d_argchk('evintrp3d',jspline,
138  > icoeff,nx,ny,nz,ixdim,iydim,izdim,ier)
139  if(ier.ne.0) return
140 c
141 C use lookup routine as in Hermite interpolation
142 C
143  call r8herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,0,0,0,
144  > i,j,k,xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi,ier)
145  if(ier.ne.0) return
146 c
147  call r8fvintrp3d(ict,1,1,
148  > fval,i,j,k,xparam,yparam,zparam,
149  > hx,hxi,hy,hyi,hz,hzi,
150  > jspline,f,icoeff,ixdim,iydim,izdim)
151 
152 C
153  return
154  end
155 C---------------------------------------------------------------------
156 C evaluate Hybrid function interpolation -- 3d fcn
157 C NO ERROR CHECKING
158 C
159  subroutine r8fvintrp3d(ict,ivec,ivecd,
160  > fval,ii,jj,kk,xparam,yparam,zparam,
161  > hx,hxi,hy,hyi,hz,hzi,
162  > jspline,fin,icoeff,ixdim,iydim,izdim)
163 C
164 C
165 C use mkintrp3d to set up spline coefficients...
166 C
167  IMPLICIT NONE
168  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
169  integer ict(10) ! requested output control
170  integer ivec ! vector length
171  integer ivecd ! vector dimension (1st dim of fval)
172 C
173  integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
174  real*8 xparam(ivec),yparam(ivec),zparam(ivec)
175  ! normalized displacements from (i,j,k) corners
176 C
177  real*8 hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
178  real*8 hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
179  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
180 C
181  integer :: jspline(3) ! interpolation type control
182  integer :: icoeff ! coefficient dimension
183  integer :: ixdim,iydim,izdim ! x & y & z dimensions
184  REAL*8 fin(icoeff,ixdim,iydim,izdim) ! data on which to interpolate...
185 C
186  real*8 fval(ivecd,10) ! output returned
187 C
188 C for detailed description of fin, ict and fval see subroutine evintrp3d
189 C comments. Note ict is not vectorized; the same output
190 C is expected to be returned for all input vector data points.
191 C
192 C note that the index inputs ii,jj,kk and parameter inputs
193 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
194 C output array fval has a vector ** 1st dimension ** whose
195 C size must be given as a separate argument
196 C
197 C to use this routine in scalar mode, pass in ivec=ivecd=1
198 C
199 C---------------
200 C local:
201 C
202  integer :: i,j,k1,k2,imaxx,imaxy,imaxz,linrank,cubrank,zonrank
203  logical :: splin_flag
204 C
205  integer :: inum
206  integer :: iderivs(3,10),ict8(8),maxd_lin(8),maxd_cub(8)
207  integer :: idcub,idlin,idcub2(2),idlin2(2)
208 C
209  real*8 :: xp,yp,hhx,hhy,hhxi,hhyi,fac1,fac2(2)
210  real*8 :: h,hi,hhi,h2(2),h2i(2)
211  real*8 :: f22(2,2),f22b(2,2),f2222(2,2,2,2),ans22(2,2),ans1,ans2
212  real*8 :: f422(4,2,2),f422b(4,2,2)
213 C
214  real*8, parameter :: one = 1.0_r8
215  real*8, parameter :: sixth = 0.166666666666666667_r8
216  real*8, parameter :: z36th = sixth*sixth
217 C
218 C---------------
219 C
220  zonrank=0
221  linrank=0
222  cubrank=0
223 
224  do i=1,3
225  if(jspline(i).eq.-1) then
226  zonrank = zonrank + 1
227  else if(jspline(i).eq.0) then
228  linrank = linrank + 1
229  else
230  splin_flag = jspline(i).eq.2
231  cubrank = cubrank + 1
232  endif
233  enddo
234 C
235  if(cubrank.eq.3) then
236  if(splin_flag) then
237  ! bicubic spline
238  call r8fvtricub(ict,ivec,ivecd,
239  > fval,ii,jj,kk,xparam,yparam,zparam,
240  > hx,hxi,hy,hyi,hz,hzi,
241  > fin,ixdim,iydim,izdim)
242 
243  else
244  ! bicubic Hermite -- translate derivative code
245  call dtrans
246  do i=1,inum
247  if(maxval(iderivs(1:3,i)).le.1) then
248  call gen_ict8
249 
250  call r8herm3fcn(ict8,ivec,ivecd,
251  > fval(1,i),ii,jj,kk,xparam,yparam,zparam,
252  > hx,hxi,hy,hyi,hz,hzi,
253  > fin,ixdim,iydim,izdim)
254 
255  else
256  fval(1:ivec,i) = 0
257  endif
258  enddo
259  endif
260 
261  return
262  endif
263 
264  call dtrans
265 
266  if(cubrank.eq.0) then
267  if(jspline(1).eq.-1) then
268  imaxx=0
269  else
270  imaxx=1
271  endif
272  if(jspline(2).eq.-1) then
273  imaxy=0
274  else
275  imaxy=1
276  endif
277  if(jspline(3).eq.-1) then
278  imaxz=0
279  else
280  imaxz=1
281  endif
282 
283  do i=1,inum
284  if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
285  > .and.(iderivs(3,i).le.imaxz)) then
286  if(linrank.eq.3) then
287  ! trilinear interpolation
288  call gen_ict8
289 
290  call r8pc3fcn(ict8,ivec,ivecd,
291  > fval(1,i), ii,jj,kk, xparam,yparam,zparam,
292  > hx,hxi, hy,hyi, hz,hzi,
293  > fin, ixdim,iydim,izdim)
294  else if(zonrank.eq.3) then
295  ! 3d step fcn
296  do j=1,ivec
297  fval(j,i)=fin(1,ii(j),jj(j),kk(j))
298  enddo
299  else if(linrank.eq.1) then
300  if(jspline(1).eq.0) then
301  ! linear in x, zonal in y,z
302  do j=1,ivec
303  if(iderivs(1,i).eq.0) then
304  fval(j,i) =
305  > (one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
306  > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j))
307  else
308  fval(j,i) = hxi(j)*
309  > (fin(1,ii(j)+1,jj(j),kk(j))-
310  > fin(1,ii(j),jj(j),kk(j)))
311  endif
312  enddo
313  else if(jspline(2).eq.0) then
314  ! linear in y, zonal in x,z
315  do j=1,ivec
316  if(iderivs(2,i).eq.0) then
317  fval(j,i) =
318  > (one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
319  > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j))
320  else
321  fval(j,i) = hyi(j)*
322  > (fin(1,ii(j),jj(j)+1,kk(j))-
323  > fin(1,ii(j),jj(j),kk(j)))
324  endif
325  enddo
326  else
327  ! linear in z, zonal in x,y
328  do j=1,ivec
329  if(iderivs(3,i).eq.0) then
330  fval(j,i) =
331  > (one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
332  > +zparam(j)*fin(1,ii(j),jj(j),kk(j)+1)
333  else
334  fval(j,i) = hzi(j)*
335  > (fin(1,ii(j),jj(j),kk(j)+1)-
336  > fin(1,ii(j),jj(j),kk(j)))
337  endif
338  enddo
339  endif
340  else
341  if(jspline(1).eq.-1) then
342  ! zonal in x, linear in y,z
343  do j=1,ivec
344  if((iderivs(2,i).eq.0).and.
345  > (iderivs(3,i).eq.0)) then
346  fval(j,i) = (one-zparam(j))*
347  > ((one-yparam(j))*
348  > fin(1,ii(j),jj(j),kk(j))
349  > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)))
350  > +zparam(j)*
351  > ((one-yparam(j))*
352  > fin(1,ii(j),jj(j),kk(j)+1)
353  > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1))
354 
355  else if((iderivs(2,i).eq.1).and.
356  > (iderivs(3,i).eq.0)) then
357  fval(j,i) = hyi(j)*(
358  > (one-zparam(j))*fin(1,ii(j),jj(j)+1,kk(j))
359  > +zparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1)
360  > -(one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
361  > -zparam(j)*fin(1,ii(j),jj(j),kk(j)+1) )
362 
363  else if((iderivs(2,i).eq.0).and.
364  > (iderivs(3,i).eq.1)) then
365  fval(j,i) = hzi(j)*(
366  > (one-yparam(j))*fin(1,ii(j),jj(j),kk(j)+1)
367  > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1)
368  > -(one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
369  > -yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)) )
370 
371  else
372  fval(j,i) = hyi(j)*hzi(j)*(
373  > fin(1,ii(j),jj(j),kk(j))
374  > -fin(1,ii(j),jj(j)+1,kk(j))
375  > -fin(1,ii(j),jj(j),kk(j)+1)
376  > +fin(1,ii(j),jj(j)+1,kk(j)+1))
377 
378  endif
379  enddo
380  else if(jspline(2).eq.-1) then
381  ! zonal in y, linear in x,z
382  do j=1,ivec
383  if((iderivs(1,i).eq.0).and.
384  > (iderivs(3,i).eq.0)) then
385  fval(j,i) = (one-zparam(j))*
386  > ((one-xparam(j))*
387  > fin(1,ii(j),jj(j),kk(j))
388  > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)))
389  > +zparam(j)*
390  > ((one-xparam(j))*
391  > fin(1,ii(j),jj(j),kk(j)+1)
392  > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1))
393 
394  else if((iderivs(1,i).eq.1).and.
395  > (iderivs(3,i).eq.0)) then
396  fval(j,i) = hxi(j)*(
397  > (one-zparam(j))*fin(1,ii(j)+1,jj(j),kk(j))
398  > +zparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1)
399  > -(one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
400  > -zparam(j)*fin(1,ii(j),jj(j),kk(j)+1) )
401 
402  else if((iderivs(1,i).eq.0).and.
403  > (iderivs(3,i).eq.1)) then
404  fval(j,i) = hzi(j)*(
405  > (one-xparam(j))*fin(1,ii(j),jj(j),kk(j)+1)
406  > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1)
407  > -(one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
408  > -xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)) )
409 
410  else
411  fval(j,i) = hxi(j)*hzi(j)*(
412  > fin(1,ii(j),jj(j),kk(j))
413  > -fin(1,ii(j)+1,jj(j),kk(j))
414  > -fin(1,ii(j),jj(j),kk(j)+1)
415  > +fin(1,ii(j)+1,jj(j),kk(j)+1))
416 
417  endif
418  enddo
419  else
420  ! zonal in z, linear in x,y
421  do j=1,ivec
422  if((iderivs(1,i).eq.0).and.
423  > (iderivs(2,i).eq.0)) then
424  fval(j,i) = (one-yparam(j))*
425  > ((one-xparam(j))*
426  > fin(1,ii(j),jj(j),kk(j))
427  > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)))
428  > +yparam(j)*
429  > ((one-xparam(j))*
430  > fin(1,ii(j),jj(j)+1,kk(j))
431  > +xparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j)))
432 
433  else if((iderivs(1,i).eq.1).and.
434  > (iderivs(2,i).eq.0)) then
435  fval(j,i) = hxi(j)*(
436  > (one-yparam(j))*fin(1,ii(j)+1,jj(j),kk(j))
437  > +yparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j))
438  > -(one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
439  > -yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)) )
440 
441  else if((iderivs(1,i).eq.0).and.
442  > (iderivs(2,i).eq.1)) then
443  fval(j,i) = hyi(j)*(
444  > (one-xparam(j))*fin(1,ii(j),jj(j)+1,kk(j))
445  > +xparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j))
446  > -(one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
447  > -xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)) )
448 
449  else
450  fval(j,i) = hxi(j)*hyi(j)*(
451  > fin(1,ii(j),jj(j),kk(j))
452  > -fin(1,ii(j)+1,jj(j),kk(j))
453  > -fin(1,ii(j),jj(j)+1,kk(j))
454  > +fin(1,ii(j)+1,jj(j)+1,kk(j)))
455 
456  endif
457  enddo
458  endif
459  endif
460  else
461  fval(1:ivec,i) = 0
462  endif
463  enddo
464 
465  return
466  endif
467 
468  ! OK -- Hybrid: Hermite or Spline in one or two dimensions
469  ! pclin or step in the other dimension(s).
470 
471  if(jspline(1).eq.2) then
472  imaxx=3
473  else if(jspline(1).eq.-1) then
474  imaxx=0
475  else
476  imaxx=1
477  endif
478 
479  if(jspline(2).eq.2) then
480  imaxy=3
481  else if(jspline(2).eq.-1) then
482  imaxy=0
483  else
484  imaxy=1
485  endif
486 
487  if(jspline(3).eq.2) then
488  imaxz=3
489  else if(jspline(3).eq.-1) then
490  imaxz=0
491  else
492  imaxz=1
493  endif
494 
495  if(cubrank.eq.1) then
496  ! cubic in 1 dimension; bilinear or lookup in other dimensions
497 
498  do i=1,inum
499  if(jspline(1).ge.1) then
500  idcub=iderivs(1,i)
501  idlin2(1)=iderivs(2,i)
502  idlin2(2)=iderivs(3,i)
503  else if(jspline(2).ge.1) then
504  idcub=iderivs(2,i)
505  idlin2(1)=iderivs(1,i)
506  idlin2(2)=iderivs(3,i)
507  else if(jspline(3).ge.1) then
508  idcub=iderivs(3,i)
509  idlin2(1)=iderivs(1,i)
510  idlin2(2)=iderivs(2,i)
511  endif
512 
513  if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
514  > .and.(iderivs(3,i).le.imaxz)) then
515  if(linrank.eq.2) then
516  do j=1,ivec
517  if(jspline(1).ge.1) then
518  xp=xparam(j)
519  fac2(1)=yparam(j)
520  fac2(2)=zparam(j)
521  h=hx(j)
522  hi=hxi(j)
523  h2i(1)=hyi(j)
524  h2i(2)=hzi(j)
525  f2222(1:2,1,1:2,1:2) =
526  > fin(1:2,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
527  f2222(1:2,2,1:2,1:2) =
528  > fin(1:2,ii(j)+1,jj(j):jj(j)+1,kk(j):kk(j)+1)
529  else if(jspline(2).ge.1) then
530  xp=yparam(j)
531  fac2(1)=xparam(j)
532  fac2(2)=zparam(j)
533  h=hy(j)
534  hi=hyi(j)
535  h2i(1)=hxi(j)
536  h2i(2)=hzi(j)
537  f2222(1:2,1,1:2,1:2) =
538  > fin(1:2,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
539  f2222(1:2,2,1:2,1:2) =
540  > fin(1:2,ii(j):ii(j)+1,jj(j)+1,kk(j):kk(j)+1)
541  else if(jspline(3).ge.1) then
542  xp=zparam(j)
543  fac2(1)=xparam(j)
544  fac2(2)=yparam(j)
545  h=hz(j)
546  hi=hzi(j)
547  h2i(1)=hxi(j)
548  h2i(2)=hyi(j)
549  f2222(1:2,1,1:2,1:2) =
550  > fin(1:2,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
551  f2222(1:2,2,1:2,1:2) =
552  > fin(1:2,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j)+1)
553  endif
554 
555  ! evaluate 2x2 array of 1d splines
556  do k2=1,2
557  do k1=1,2
558  if(splin_flag) then
559  call sp1d(ans22(k1,k2),
560  > f2222(1:2,1:2,k1,k2))
561  else
562  call hm1d(ans22(k1,k2),
563  > f2222(1:2,1:2,k1,k2))
564  endif
565  enddo
566  enddo
567 
568  ! bilinear interpolation
569  if((idlin2(1).eq.0).and.(idlin2(2).eq.0)) then
570  ! value
571  ans1 = (one-fac2(2))*
572  > ((one-fac2(1))*ans22(1,1)+fac2(1)*ans22(2,1))
573  > +fac2(2)*
574  > ((one-fac2(1))*ans22(1,2)+fac2(1)*ans22(2,2))
575  else if((idlin2(1).eq.1).and.(idlin2(2).eq.0)) then
576  ! 1st deriv along 1st dimension
577  ans1 = (one-fac2(2))*(ans22(2,1)-ans22(1,1))
578  > +fac2(2)*(ans22(2,2)-ans22(1,2))
579  ans1 = ans1*h2i(1)
580  else if((idlin2(1).eq.0).and.(idlin2(2).eq.1)) then
581  ! 1st deriv along 2nd dimension
582  ans1 = (one-fac2(1))*(ans22(1,2)-ans22(1,1))
583  > +fac2(1)*(ans22(2,2)-ans22(2,1))
584  ans1 = ans1*h2i(2)
585  else
586  ! derivatives along both dimensions
587  ans1 = ans22(2,2)-ans22(1,2)
588  > -ans22(2,1)+ans22(1,1)
589  ans1 = ans1*h2i(1)*h2i(2)
590  endif
591  fval(j,i)=ans1
592  enddo
593 
594  else if(linrank.eq.0) then
595  do j=1,ivec
596  if(jspline(1).ge.1) then
597  xp=xparam(j)
598  h=hx(j)
599  hi=hxi(j)
600  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
601  f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
602  else if(jspline(2).ge.1) then
603  xp=yparam(j)
604  h=hy(j)
605  hi=hyi(j)
606  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
607  f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
608  else if(jspline(3).ge.1) then
609  xp=zparam(j)
610  h=hz(j)
611  hi=hzi(j)
612  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
613  f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
614  endif
615  if(splin_flag) then
616  call sp1d(ans1,f22)
617  else
618  call hm1d(ans1,f22)
619  endif
620  fval(j,i)=ans1
621  enddo
622 
623  else if(linrank.eq.1) then
624  do j=1,ivec
625  if(jspline(1).ge.1) then
626  xp=xparam(j)
627  h=hx(j)
628  hi=hxi(j)
629  if(jspline(2).eq.0) then
630  hhi=hyi(j)
631  fac1=yparam(j)
632  idlin=iderivs(2,i)
633  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
634  f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
635  f22b(1:2,1)=fin(1:2,ii(j),jj(j)+1,kk(j))
636  f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j)+1,kk(j))
637  else
638  hhi=hzi(j)
639  fac1=zparam(j)
640  idlin=iderivs(3,i)
641  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
642  f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
643  f22b(1:2,1)=fin(1:2,ii(j),jj(j),kk(j)+1)
644  f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j)+1)
645  endif
646  else if(jspline(2).ge.1) then
647  xp=yparam(j)
648  h=hy(j)
649  hi=hyi(j)
650  if(jspline(1).eq.0) then
651  hhi=hxi(j)
652  fac1=xparam(j)
653  idlin=iderivs(1,i)
654  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
655  f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
656  f22b(1:2,1)=fin(1:2,ii(j)+1,jj(j),kk(j))
657  f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j)+1,kk(j))
658  else
659  hhi=hzi(j)
660  fac1=zparam(j)
661  idlin=iderivs(3,i)
662  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
663  f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
664  f22b(1:2,1)=fin(1:2,ii(j),jj(j),kk(j)+1)
665  f22b(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j)+1)
666  endif
667  else if(jspline(3).ge.1) then
668  xp=zparam(j)
669  h=hz(j)
670  hi=hzi(j)
671  if(jspline(1).eq.0) then
672  hhi=hxi(j)
673  fac1=xparam(j)
674  idlin=iderivs(1,i)
675  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
676  f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
677  f22b(1:2,1)=fin(1:2,ii(j)+1,jj(j),kk(j))
678  f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j)+1)
679  else
680  hhi=hyi(j)
681  fac1=yparam(j)
682  idlin=iderivs(2,i)
683  f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
684  f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
685  f22b(1:2,1)=fin(1:2,ii(j),jj(j)+1,kk(j))
686  f22b(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j)+1)
687  endif
688  endif
689  if(splin_flag) then
690  call sp1d(ans1,f22)
691  call sp1d(ans2,f22b)
692  else
693  call hm1d(ans1,f22)
694  call hm1d(ans2,f22b)
695  endif
696  if(idlin.eq.0) then
697  fval(j,i)=(one-fac1)*ans1 + fac1*ans2
698  else
699  fval(j,i)=(ans2-ans1)*hhi
700  endif
701  enddo
702 
703  endif
704  else
705  fval(1:ivec,i) = 0
706  endif
707  enddo
708  return
709  endif
710 
711  ! OK: cubic in 2 dimensions; linear or lookup in other dimension.
712 
713  do i=1,inum
714  if(jspline(1).le.0) then
715  idlin=iderivs(1,i)
716  idcub2(1)=iderivs(2,i)
717  idcub2(2)=iderivs(3,i)
718  else if(jspline(2).le.0) then
719  idlin=iderivs(2,i)
720  idcub2(1)=iderivs(1,i)
721  idcub2(2)=iderivs(3,i)
722  else if(jspline(3).le.0) then
723  idlin=iderivs(3,i)
724  idcub2(1)=iderivs(1,i)
725  idcub2(2)=iderivs(2,i)
726  endif
727 
728  if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
729  > .and.(iderivs(3,i).le.imaxz)) then
730 
731  if(linrank.eq.1) then
732  ! linear interpolation in non-spline dimension
733  do j=1,ivec
734  if(jspline(1).le.0) then
735  fac1=xparam(j)
736  hi=hxi(j)
737  xp=yparam(j)
738  yp=zparam(j)
739  hhx=hy(j)
740  hhxi=hyi(j)
741  hhy=hz(j)
742  hhyi=hzi(j)
743  f422 = fin(1:4,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
744  f422b= fin(1:4,ii(j)+1,jj(j):jj(j)+1,kk(j):kk(j)+1)
745  else if(jspline(2).le.0) then
746  fac1=yparam(j)
747  hi=hyi(j)
748  xp=xparam(j)
749  yp=zparam(j)
750  hhx=hx(j)
751  hhxi=hxi(j)
752  hhy=hz(j)
753  hhyi=hzi(j)
754  f422 = fin(1:4,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
755  f422b= fin(1:4,ii(j):ii(j)+1,jj(j)+1,kk(j):kk(j)+1)
756  else if(jspline(3).le.0) then
757  fac1=zparam(j)
758  hi=hzi(j)
759  xp=xparam(j)
760  yp=yparam(j)
761  hhx=hx(j)
762  hhxi=hxi(j)
763  hhy=hy(j)
764  hhyi=hyi(j)
765  f422 = fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
766  f422b= fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j)+1)
767  endif
768  if(splin_flag) then
769  call sp2d(ans1,f422)
770  call sp2d(ans2,f422b)
771  else
772  call hm2d(ans1,f422)
773  call hm2d(ans2,f422b)
774  endif
775  if(idlin.eq.1) then
776  fval(j,i) = (ans2-ans1)*hi
777  else
778  fval(j,i) = (one-fac1)*ans1 + fac1*ans2
779  endif
780  enddo
781  else
782  ! lookup in the non-spline dimension
783  do j=1,ivec
784  if(jspline(1).le.0) then
785  xp=yparam(j)
786  yp=zparam(j)
787  hhx=hy(j)
788  hhxi=hyi(j)
789  hhy=hz(j)
790  hhyi=hzi(j)
791  f422 = fin(1:4,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
792  else if(jspline(2).le.0) then
793  xp=xparam(j)
794  yp=zparam(j)
795  hhx=hx(j)
796  hhxi=hxi(j)
797  hhy=hz(j)
798  hhyi=hzi(j)
799  f422 = fin(1:4,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
800  else if(jspline(3).le.0) then
801  xp=xparam(j)
802  yp=yparam(j)
803  hhx=hx(j)
804  hhxi=hxi(j)
805  hhy=hy(j)
806  hhyi=hyi(j)
807  f422 = fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
808  endif
809  if(splin_flag) then
810  call sp2d(ans1,f422)
811  else
812  call hm2d(ans1,f422)
813  endif
814  fval(j,i) = ans1
815  enddo
816  endif
817  else
818  fval(1:ivec,i) = 0
819  endif
820  enddo
821  contains
822  subroutine gen_ict8
823  ict8=0
824  if(iderivs(1,i).eq.0) then
825  if(iderivs(2,i).eq.0) then
826  if(iderivs(3,i).eq.0) then
827  ict8(1)=1 ! f
828  else
829  ict8(4)=1 ! df/dz
830  endif
831  else
832  if(iderivs(3,i).eq.0) then
833  ict8(3)=1 ! df/dy
834  else
835  ict8(7)=1 ! d2f/dydz
836  endif
837  endif
838  else
839  if(iderivs(2,i).eq.0) then
840  if(iderivs(3,i).eq.0) then
841  ict8(2)=2 ! df/dx
842  else
843  ict8(6)=6 ! d2f/dxdz
844  endif
845  else
846  if(iderivs(3,i).eq.0) then
847  ict8(5)=1 ! d2f/dxdy
848  else
849  ict8(8)=1 ! d3f/dxdydz
850  endif
851  endif
852  endif
853  end subroutine gen_ict8
854 
855  subroutine dtrans
856 
857  ! convert ict(...) codes into intermediate coding:
858  ! iderivs(i,j) = number of differentiations in coordinate i
859  ! of the j'th value evaluated
860 
861  ! maxd_lin(j) = max # of differentiations for a linear dimension
862  ! maxd_cub(j) = max # of differentiations for a spline dimension
863 
864  integer :: i
865 
866  inum=0 ! actual number of vectors to be evaluated (summed here)
867 
868  if(abs(ict(1)).le.(2)) then
869  if(ict(1).eq.1) then
870  call add1(0,0,0)
871  endif
872  if(ict(2).eq.1) then
873  call add1(1,0,0)
874  endif
875  if(ict(3).eq.1) then
876  call add1(0,1,0)
877  endif
878  if(ict(4).eq.1) then
879  call add1(0,0,1)
880  endif
881  if(ict(5).eq.1) then
882  call add1(2,0,0)
883  endif
884  if(ict(6).eq.1) then
885  call add1(0,2,0)
886  endif
887  if(ict(7).eq.1) then
888  call add1(0,0,2)
889  endif
890  if(ict(8).eq.1) then
891  call add1(1,1,0)
892  endif
893  if(ict(9).eq.1) then
894  call add1(1,0,1)
895  endif
896  if(ict(10).eq.1) then
897  call add1(0,1,1)
898  endif
899 
900  else if(ict(1).eq.3) then
901  if(ict(2).eq.1) then
902  call add1(2,1,0)
903  endif
904  if(ict(3).eq.1) then
905  call add1(2,0,1)
906  endif
907  if(ict(4).eq.1) then
908  call add1(1,2,0)
909  endif
910  if(ict(5).eq.1) then
911  call add1(1,1,1)
912  endif
913  if(ict(6).eq.1) then
914  call add1(1,0,2)
915  endif
916  if(ict(7).eq.1) then
917  call add1(0,2,1)
918  endif
919  if(ict(8).eq.1) then
920  call add1(0,1,2)
921  endif
922 
923  else if(ict(1).eq.-3) then
924  if(ict(2).eq.1) then
925  call add1(3,0,0)
926  endif
927  if(ict(3).eq.1) then
928  call add1(0,3,0)
929  endif
930  if(ict(4).eq.1) then
931  call add1(0,0,3)
932  endif
933 
934  else if(ict(1).eq.4) then
935  if(ict(2).eq.1) then
936  call add1(2,2,0)
937  endif
938  if(ict(3).eq.1) then
939  call add1(2,1,1)
940  endif
941  if(ict(4).eq.1) then
942  call add1(2,0,2)
943  endif
944  if(ict(5).eq.1) then
945  call add1(1,2,1)
946  endif
947  if(ict(6).eq.1) then
948  call add1(1,1,2)
949  endif
950  if(ict(7).eq.1) then
951  call add1(0,2,2)
952  endif
953 
954  else if(ict(1).eq.-4) then
955  if(ict(2).eq.1) then
956  call add1(3,1,0)
957  endif
958  if(ict(3).eq.1) then
959  call add1(3,0,1)
960  endif
961  if(ict(4).eq.1) then
962  call add1(1,3,0)
963  endif
964  if(ict(5).eq.1) then
965  call add1(1,0,3)
966  endif
967  if(ict(6).eq.1) then
968  call add1(0,3,1)
969  endif
970  if(ict(7).eq.1) then
971  call add1(0,1,3)
972  endif
973 
974  else if(ict(1).eq.5) then
975  if(ict(2).eq.1) then
976  call add1(2,2,1)
977  endif
978  if(ict(3).eq.1) then
979  call add1(2,1,2)
980  endif
981  if(ict(4).eq.1) then
982  call add1(1,2,2)
983  endif
984 
985  else if(ict(1).eq.-5) then
986  if(ict(2).eq.1) then
987  call add1(3,2,0)
988  endif
989  if(ict(3).eq.1) then
990  call add1(3,1,1)
991  endif
992  if(ict(4).eq.1) then
993  call add1(3,0,2)
994  endif
995  if(ict(5).eq.1) then
996  call add1(2,3,0)
997  endif
998  if(ict(6).eq.1) then
999  call add1(2,0,3)
1000  endif
1001  if(ict(7).eq.1) then
1002  call add1(1,3,1)
1003  endif
1004  if(ict(8).eq.1) then
1005  call add1(1,1,3)
1006  endif
1007  if(ict(9).eq.1) then
1008  call add1(0,3,2)
1009  endif
1010  if(ict(10).eq.1) then
1011  call add1(0,2,3)
1012  endif
1013 
1014  else if(ict(1).eq.6) then
1015  call add1(2,2,2)
1016 
1017  else if(ict(1).eq.-6) then
1018  if(ict(2).eq.1) then
1019  call add1(3,3,0)
1020  endif
1021  if(ict(3).eq.1) then
1022  call add1(3,2,1)
1023  endif
1024  if(ict(4).eq.1) then
1025  call add1(3,1,2)
1026  endif
1027  if(ict(5).eq.1) then
1028  call add1(3,0,3)
1029  endif
1030  if(ict(6).eq.1) then
1031  call add1(2,3,1)
1032  endif
1033  if(ict(7).eq.1) then
1034  call add1(2,1,3)
1035  endif
1036  if(ict(8).eq.1) then
1037  call add1(1,3,2)
1038  endif
1039  if(ict(9).eq.1) then
1040  call add1(1,2,3)
1041  endif
1042  if(ict(10).eq.1) then
1043  call add1(0,3,3)
1044  endif
1045 
1046  else if(abs(ict(1)).eq.7) then
1047  if(ict(2).eq.1) then
1048  call add1(3,3,1)
1049  endif
1050  if(ict(3).eq.1) then
1051  call add1(3,2,2)
1052  endif
1053  if(ict(4).eq.1) then
1054  call add1(3,1,3)
1055  endif
1056  if(ict(5).eq.1) then
1057  call add1(2,3,2)
1058  endif
1059  if(ict(6).eq.1) then
1060  call add1(2,2,3)
1061  endif
1062  if(ict(7).eq.1) then
1063  call add1(1,3,3)
1064  endif
1065 
1066  else if(abs(ict(1)).eq.8) then
1067  if(ict(2).eq.1) then
1068  call add1(3,3,2)
1069  endif
1070  if(ict(3).eq.1) then
1071  call add1(3,2,3)
1072  endif
1073  if(ict(4).eq.1) then
1074  call add1(2,3,3)
1075  endif
1076 
1077  else if(abs(ict(1)).eq.9) then
1078  call add1(3,3,3)
1079 
1080  endif
1081 
1082  end subroutine dtrans
1083 
1084  subroutine add1(idx,idy,idz)
1085  ! insert record of derivative d[idx+idy+idz]f/dx[idx]dy[idy]dz[idz]
1086 
1087  integer, intent(in) :: idx,idy,idz
1088 
1089  inum=inum+1
1090 
1091  iderivs(1,inum)=idx
1092  if(jspline(1).le.0) then
1093  maxd_lin(inum)=idx
1094  maxd_cub(inum)=0
1095  else
1096  maxd_lin(inum)=0
1097  maxd_cub(inum)=idx
1098  endif
1099 
1100  iderivs(2,inum)=idy
1101  if(jspline(2).le.0) then
1102  maxd_lin(inum)=max(maxd_lin(inum),idy)
1103  else
1104  maxd_cub(inum)=max(maxd_cub(inum),idy)
1105  endif
1106 
1107 
1108  iderivs(3,inum)=idz
1109  if(jspline(3).le.0) then
1110  maxd_lin(inum)=max(maxd_lin(inum),idz)
1111  else
1112  maxd_cub(inum)=max(maxd_cub(inum),idz)
1113  endif
1114 
1115  end subroutine add1
1116 
1117 C===========================================================================
1118  subroutine sp1d(ans,f22)
1119  real*8, intent(in) :: f22(2,2)
1120  real*8, intent(out) :: ans
1121 
1122 C----------------------
1123 C contained 1d spline evaluation
1124 C----------------------
1125 
1126  real*8 :: xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
1127 C
1128 C----------------------
1129 
1130  xpi=1.0_r8-xp
1131  xp2=xp*xp
1132  xpi2=xpi*xpi
1133 
1134  if(idcub.eq.0) then
1135  ! value
1136 
1137  cx=xp*(xp2-1.0_r8)
1138  cxi=xpi*(xpi2-1.0_r8)
1139  hx2=h*h
1140 
1141  ans=xpi*f22(1,1) +xp*f22(1,2)
1142  ans=ans+sixth*hx2*(cxi*f22(2,1) +cx*f22(2,2))
1143 
1144  else if(idcub.eq.1) then
1145  ! 1st derivative
1146 
1147  cxd=3.0_r8*xp2-1.0_r8
1148  cxdi=-3.0_r8*xpi2+1.0_r8
1149 
1150  ans=hi*(f22(1,2)-f22(1,1))
1151  ans=ans+sixth*h*(cxdi*f22(2,1) +cxd*f22(2,2))
1152 
1153  else if(idcub.eq.2) then
1154  ! 2nd derivative
1155 
1156  xpi=1.0_r8-xp
1157  ans=xpi*f22(2,1) + xp*f22(2,2)
1158 
1159  else
1160  ! 3rd derivative
1161 
1162  ans = hi*(f22(2,2)-f22(2,1))
1163 
1164  endif
1165 
1166  end subroutine sp1d
1167 
1168 C===========================================================================
1169  subroutine hm1d(ans,f22)
1170  real*8, intent(in) :: f22(2,2)
1171  real*8, intent(out) :: ans
1172 
1173 C----------------------
1174 C contained 1d hermite evaluation
1175 C----------------------
1176 
1177  real*8 :: xpi,xp2,xpi2,ax,axbar,bx,bxbar
1178  real*8 :: axp,axbarp,bxp,bxbarp
1179 
1180  xpi=1.0_r8-xp
1181  xp2=xp*xp
1182  xpi2=xpi*xpi
1183 
1184  if(idcub.eq.0) then
1185  ! value
1186 
1187  ax=xp2*(3.0_r8-2.0_r8*xp)
1188  axbar=1.0_r8-ax
1189  bx=-xp2*xpi
1190  bxbar=xpi2*xp
1191 
1192  ans=axbar*f22(1,1) + ax*f22(1,2)
1193  ans=ans + h*(bxbar*f22(2,1) + bx*f22(2,2))
1194 
1195  else
1196  ! 1st derivative
1197 
1198  axp=6.0_r8*xp*xpi
1199  axbarp=-axp
1200  bxp=xp*(3.0_r8*xp-2.0_r8)
1201  bxbarp=xpi*(3.0_r8*xpi-2.0_r8)
1202 
1203  ans=hi*(axbarp*f22(1,1) +axp*f22(1,2))
1204  ans=ans + bxbarp*f22(2,1) + bxp*f22(2,2)
1205 
1206  endif
1207 
1208  end subroutine hm1d
1209 
1210 C===========================================================================
1211  subroutine sp2d(ans,ff)
1212  real*8, intent(in) :: ff(0:3,2,2)
1213  real*8, intent(out) :: ans
1214 
1215 C----------------------
1216 C contained 2d spline evaluation
1217 C----------------------
1218 
1219  real*8 xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
1220  real*8 ypi,yp2,ypi2,cy,cyi,hy2,cyd,cydi
1221 
1222 C----------------------
1223 
1224  xpi=1.0_r8-xp
1225  xp2=xp*xp
1226  xpi2=xpi*xpi
1227 
1228  ypi=1.0_r8-yp
1229  yp2=yp*yp
1230  ypi2=ypi*ypi
1231 C
1232 C get desired values:
1233 C
1234  if((idcub2(1).eq.0).and.(idcub2(2).eq.0)) then
1235 
1236 C function value:
1237 
1238  cy=yp*(yp2-1.0_r8)
1239  cyi=ypi*(ypi2-1.0_r8)
1240  hy2=hhy*hhy
1241 
1242  cx=xp*(xp2-1.0_r8)
1243  cxi=xpi*(xpi2-1.0_r8)
1244  hx2=hhx*hhx
1245 
1246  ans=xpi*(ypi*ff(0,1,1) +yp*ff(0,1,2))+
1247  > xp*(ypi*ff(0,2,1)+yp*ff(0,2,2))
1248 
1249  ans=ans+sixth*hx2*(
1250  > cxi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1251  > cx*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1252 
1253  ans=ans+sixth*hy2*(
1254  > xpi*(cyi*ff(2,1,1) +cy*ff(2,1,2))+
1255  > xp*(cyi*ff(2,2,1)+cy*ff(2,2,2)))
1256 
1257  ans=ans+z36th*hx2*hy2*(
1258  > cxi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1259  > cx*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1260 
1261 
1262  else if((idcub2(1).eq.1).and.(idcub2(2).eq.0)) then
1263 
1264 C df/dx:
1265 
1266  cy=yp*(yp2-1.0_r8)
1267  cyi=ypi*(ypi2-1.0_r8)
1268  hy2=hhy*hhy
1269 
1270  cxd=3.0_r8*xp2-1.0_r8
1271  cxdi=-3.0_r8*xpi2+1.0_r8
1272 
1273  ans=hhxi*(
1274  > -(ypi*ff(0,1,1) +yp*ff(0,1,2))
1275  > +(ypi*ff(0,2,1)+yp*ff(0,2,2)))
1276 
1277  ans=ans+sixth*hhx*(
1278  > cxdi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1279  > cxd*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1280 
1281  ans=ans+sixth*hhxi*hy2*(
1282  > -(cyi*ff(2,1,1) +cy*ff(2,1,2))
1283  > +(cyi*ff(2,2,1)+cy*ff(2,2,2)))
1284 
1285  ans=ans+z36th*hhx*hy2*(
1286  > cxdi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1287  > cxd*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1288 
1289 
1290  else if((idcub2(1).eq.0).and.(idcub2(2).eq.1)) then
1291 
1292 C df/dy:
1293 
1294  cyd=3.0_r8*yp2-1.0_r8
1295  cydi=-3.0_r8*ypi2+1.0_r8
1296 
1297  cx=xp*(xp2-1.0_r8)
1298  cxi=xpi*(xpi2-1.0_r8)
1299  hx2=hhx*hhx
1300 
1301  ans=hhyi*(
1302  > xpi*(-ff(0,1,1) +ff(0,1,2))+
1303  > xp*(-ff(0,2,1)+ff(0,2,2)))
1304 
1305  ans=ans+sixth*hx2*hhyi*(
1306  > cxi*(-ff(1,1,1) +ff(1,1,2))+
1307  > cx*(-ff(1,2,1)+ff(1,2,2)))
1308 
1309  ans=ans+sixth*hhy*(
1310  > xpi*(cydi*ff(2,1,1) +cyd*ff(2,1,2))+
1311  > xp*(cydi*ff(2,2,1)+cyd*ff(2,2,2)))
1312 
1313  ans=ans+z36th*hx2*hhy*(
1314  > cxi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1315  > cx*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1316 
1317 
1318  else if((idcub2(1).eq.2).and.(idcub2(2).eq.0)) then
1319 
1320 C d2f/dx2:
1321 
1322  cy=yp*(yp2-1.0_r8)
1323  cyi=ypi*(ypi2-1.0_r8)
1324  hy2=hhy*hhy
1325 
1326  ans=(
1327  > xpi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1328  > xp*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1329 
1330  ans=ans+sixth*hy2*(
1331  > xpi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1332  > xp*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1333 
1334 
1335  else if((idcub2(1).eq.0).and.(idcub2(2).eq.2)) then
1336 
1337 C d2f/dy2:
1338 
1339  cx=xp*(xp2-1.0_r8)
1340  cxi=xpi*(xpi2-1.0_r8)
1341  hx2=hhx*hhx
1342 
1343  ans=(
1344  > xpi*(ypi*ff(2,1,1) +yp*ff(2,1,2))+
1345  > xp*(ypi*ff(2,2,1)+yp*ff(2,2,2)))
1346 
1347  ans=ans+sixth*hx2*(
1348  > cxi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1349  > cx*(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1350 
1351 
1352  else if((idcub2(1).eq.1).and.(idcub2(2).eq.1)) then
1353 
1354 C d2f/dxdy:
1355 
1356  cyd=3.0_r8*yp2-1.0_r8
1357  cydi=-3.0_r8*ypi2+1.0_r8
1358 
1359  cxd=3.0_r8*xp2-1.0_r8
1360  cxdi=-3.0_r8*xpi2+1.0_r8
1361 
1362  ans=hhxi*hhyi*(
1363  > ff(0,1,1) -ff(0,1,2)
1364  > -ff(0,2,1)+ff(0,2,2))
1365 
1366  ans=ans+sixth*hhx*hhyi*(
1367  > cxdi*(-ff(1,1,1) +ff(1,1,2))+
1368  > cxd*(-ff(1,2,1)+ff(1,2,2)))
1369 
1370  ans=ans+sixth*hhxi*hhy*(
1371  > -(cydi*ff(2,1,1) +cyd*ff(2,1,2))
1372  > +(cydi*ff(2,2,1)+cyd*ff(2,2,2)))
1373 
1374  ans=ans+z36th*hhx*hhy*(
1375  > cxdi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1376  > cxd*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1377 
1378  else if((idcub2(1).eq.3).and.(idcub2(2).eq.0)) then
1379 
1380 c evaluate d3f/dx3 (not continuous):
1381 
1382  cy=yp*(yp2-1.0_r8)
1383  cyi=ypi*(ypi2-1.0_r8)
1384  hy2=hhy*hhy
1385 
1386  ans=hhxi*(
1387  > -(ypi*ff(1,1,1) +yp*ff(1,1,2))
1388  > +(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1389 
1390  ans=ans+sixth*hy2*hhxi*(
1391  > -(cyi*ff(3,1,1) +cy*ff(3,1,2))
1392  > +(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1393 
1394 
1395  else if((idcub2(1).eq.2).and.(idcub2(2).eq.1)) then
1396 
1397 c evaluate d3f/dx2dy
1398 
1399  cyd=3.0_r8*yp2-1.0_r8
1400  cydi=-3.0_r8*ypi2+1.0_r8
1401 
1402  ans=hhyi*(
1403  > xpi*(-ff(1,1,1) +ff(1,1,2))+
1404  > xp*(-ff(1,2,1) +ff(1,2,2)))
1405 
1406  ans=ans+sixth*hhy*(
1407  > xpi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1408  > xp*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1409 
1410 
1411  else if((idcub2(1).eq.1).and.(idcub2(2).eq.2)) then
1412 
1413 c evaluate d3f/dxdy2
1414 
1415  cxd=3.0_r8*xp2-1.0_r8
1416  cxdi=-3.0_r8*xpi2+1.0_r8
1417 
1418  ans=hhxi*(
1419  > -(ypi*ff(2,1,1) +yp*ff(2,1,2))
1420  > +(ypi*ff(2,2,1)+yp*ff(2,2,2)))
1421 
1422  ans=ans+sixth*hhx*(
1423  > cxdi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1424  > cxd*(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1425 
1426  else if((idcub2(1).eq.0).and.(idcub2(2).eq.3)) then
1427 
1428 c evaluate d3f/dy3 (not continuous):
1429 
1430  cx=xp*(xp2-1.0_r8)
1431  cxi=xpi*(xpi2-1.0_r8)
1432  hx2=hhx*hhx
1433 
1434  ans=hhyi*(
1435  > xpi*(-ff(2,1,1) +ff(2,1,2))+
1436  > xp*(-ff(2,2,1) +ff(2,2,2)))
1437 
1438  ans=ans+sixth*hx2*hhyi*(
1439  > cxi*(-ff(3,1,1) +ff(3,1,2))+
1440  > cx*(-ff(3,2,1) +ff(3,2,2)))
1441 
1442 
1443  else if((idcub2(1).eq.3).and.(idcub2(2).eq.1)) then
1444 
1445 c evaluate d4f/dx3dy (not continuous):
1446 
1447  cyd=3.0_r8*yp2-1.0_r8
1448  cydi=-3.0_r8*ypi2+1.0_r8
1449 
1450  ans=hhxi*hhyi*(
1451  > +( ff(1,1,1) -ff(1,1,2))
1452  > +(-ff(1,2,1)+ff(1,2,2)))
1453 
1454  ans=ans+sixth*hhy*hhxi*(
1455  > -(cydi*ff(3,1,1) +cyd*ff(3,1,2))
1456  > +(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1457 
1458 
1459  else if((idcub2(1).eq.2).and.(idcub2(2).eq.2)) then
1460 
1461 c evaluate d4f/dx2dy2
1462 
1463  ans=xpi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1464  > xp*(ypi*ff(3,2,1) +yp*ff(3,2,2))
1465 
1466  else if((idcub2(1).eq.1).and.(idcub2(2).eq.3)) then
1467 
1468 c evaluate d4f/dxdy3 (not continuous):
1469 
1470  cxd=3.0_r8*xp2-1.0_r8
1471  cxdi=-3.0_r8*xpi2+1.0_r8
1472 
1473  ans=hhyi*hhxi*(
1474  > +( ff(2,1,1) -ff(2,1,2))
1475  > +(-ff(2,2,1)+ff(2,2,2)))
1476 
1477  ans=ans+sixth*hhx*hhyi*(
1478  > cxdi*(-ff(3,1,1) +ff(3,1,2))+
1479  > cxd*(-ff(3,2,1) +ff(3,2,2)))
1480 
1481 
1482  else if((idcub2(1).eq.3).and.(idcub2(2).eq.2)) then
1483 
1484 c evaluate d5f/dx3dy2 (not continuous):
1485 
1486  ans=hhxi*(
1487  > -(ypi*ff(3,1,1) +yp*ff(3,1,2))
1488  > +(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1489 
1490  else if((idcub2(1).eq.2).and.(idcub2(2).eq.3)) then
1491 
1492 c evaluate d5f/dx2dy3 (not continuous):
1493 
1494  ans=hhyi*(
1495  > xpi*(-ff(3,1,1) +ff(3,1,2))+
1496  > xp*(-ff(3,2,1)+ff(3,2,2)))
1497 
1498 
1499  else if((idcub2(1).eq.3).and.(idcub2(2).eq.3)) then
1500 
1501 c evaluate d6f/dx3dy3 (not continuous):
1502 
1503  ans=hhxi*hhyi*(
1504  > +( ff(3,1,1) -ff(3,1,2))
1505  > +(-ff(3,2,1) +ff(3,2,2)))
1506 
1507  endif
1508 
1509  end subroutine sp2d
1510 
1511 C===========================================================================
1512  subroutine hm2d(ans,ff)
1513  real*8, intent(in) :: ff(0:3,2,2)
1514  real*8, intent(out) :: ans
1515 
1516 C----------------------
1517 C contained 2d hermite evaluation
1518 C----------------------
1519 
1520  real*8 :: xpi,xp2,xpi2,ax,axbar,bx,bxbar
1521  real*8 :: axp,axbarp,bxp,bxbarp
1522 
1523  real*8 :: ypi,yp2,ypi2,ay,aybar,by,bybar
1524  real*8 :: ayp,aybarp,byp,bybarp
1525 
1526  xpi=1.0_r8-xp
1527  xp2=xp*xp
1528  xpi2=xpi*xpi
1529 
1530  ypi=1.0_r8-yp
1531  yp2=yp*yp
1532  ypi2=ypi*ypi
1533 
1534  if((idcub2(1).eq.0).and.(idcub2(2).eq.0)) then
1535  ! evaluation
1536 
1537  ax=xp2*(3.0_r8-2.0_r8*xp)
1538  axbar=1.0_r8-ax
1539  bx=-xp2*xpi
1540  bxbar=xpi2*xp
1541 
1542  ay=yp2*(3.0_r8-2.0_r8*yp)
1543  aybar=1.0_r8-ay
1544  by=-yp2*ypi
1545  bybar=ypi2*yp
1546 
1547  ans=axbar*(aybar*ff(0,1,1) +ay*ff(0,1,2))+
1548  > ax*(aybar*ff(0,2,1)+ay*ff(0,2,2))
1549 
1550  ans=ans+hhx*(
1551  > bxbar*(aybar*ff(1,1,1) +ay*ff(1,1,2))+
1552  > bx*(aybar*ff(1,2,1)+ay*ff(1,2,2)))
1553 
1554  ans=ans+hhy*(
1555  > axbar*(bybar*ff(2,1,1) +by*ff(2,1,2))+
1556  > ax*(bybar*ff(2,2,1)+by*ff(2,2,2)))
1557 
1558  ans=ans+hhx*hhy*(
1559  > bxbar*(bybar*ff(3,1,1) +by*ff(3,1,2))+
1560  > bx*(bybar*ff(3,2,1)+by*ff(3,2,2)))
1561 
1562  else if((idcub2(1).eq.1).and.(idcub2(2).eq.0)) then
1563  ! 1st deriv, 1st dim
1564 
1565  axp=6.0_r8*xp*xpi
1566  axbarp=-axp
1567  bxp=xp*(3.0_r8*xp-2.0_r8)
1568  bxbarp=xpi*(3.0_r8*xpi-2.0_r8)
1569 
1570  ay=yp2*(3.0_r8-2.0_r8*yp)
1571  aybar=1.0_r8-ay
1572  by=-yp2*ypi
1573  bybar=ypi2*yp
1574 
1575  ans=hhxi*(
1576  > axbarp*(aybar*ff(0,1,1) +ay*ff(0,1,2))+
1577  > axp*(aybar*ff(0,2,1)+ay*ff(0,2,2)))
1578 
1579  ans=ans+
1580  > bxbarp*(aybar*ff(1,1,1) +ay*ff(1,1,2))+
1581  > bxp*(aybar*ff(1,2,1)+ay*ff(1,2,2))
1582 
1583  ans=ans+hhxi*hhy*(
1584  > axbarp*(bybar*ff(2,1,1) +by*ff(2,1,2))+
1585  > axp*(bybar*ff(2,2,1)+by*ff(2,2,2)))
1586 
1587  ans=ans+hhy*(
1588  > bxbarp*(bybar*ff(3,1,1) +by*ff(3,1,2))+
1589  > bxp*(bybar*ff(3,2,1)+by*ff(3,2,2)))
1590 
1591  else if((idcub2(1).eq.0).and.(idcub2(2).eq.1)) then
1592  ! 1st deriv, 2nd dim
1593 
1594  ax=xp2*(3.0_r8-2.0_r8*xp)
1595  axbar=1.0_r8-ax
1596  bx=-xp2*xpi
1597  bxbar=xpi2*xp
1598 
1599  ayp=6.0_r8*yp*ypi
1600  aybarp=-ayp
1601  byp=yp*(3.0_r8*yp-2.0_r8)
1602  bybarp=ypi*(3.0_r8*ypi-2.0_r8)
1603 
1604  ans=hhyi*(
1605  > axbar*(aybarp*ff(0,1,1) +ayp*ff(0,1,2))+
1606  > ax*(aybarp*ff(0,2,1)+ayp*ff(0,2,2)))
1607 
1608  ans=ans+hhx*hhyi*(
1609  > bxbar*(aybarp*ff(1,1,1) +ayp*ff(1,1,2))+
1610  > bx*(aybarp*ff(1,2,1)+ayp*ff(1,2,2)))
1611 
1612  ans=ans+
1613  > axbar*(bybarp*ff(2,1,1) +byp*ff(2,1,2))+
1614  > ax*(bybarp*ff(2,2,1)+byp*ff(2,2,2))
1615 
1616  ans=ans+hhx*(
1617  > bxbar*(bybarp*ff(3,1,1) +byp*ff(3,1,2))+
1618  > bx*(bybarp*ff(3,2,1)+byp*ff(3,2,2)))
1619 
1620  else if((idcub2(1).eq.1).and.(idcub2(2).eq.1)) then
1621  ! 1st deriv, in both dimensions
1622 
1623  axp=6.0_r8*xp*xpi
1624  axbarp=-axp
1625  bxp=xp*(3.0_r8*xp-2.0_r8)
1626  bxbarp=xpi*(3.0_r8*xpi-2.0_r8)
1627 
1628  ayp=6.0_r8*yp*ypi
1629  aybarp=-ayp
1630  byp=yp*(3.0_r8*yp-2.0_r8)
1631  bybarp=ypi*(3.0_r8*ypi-2.0_r8)
1632 
1633  ans=hhxi*hhyi*(
1634  > axbarp*(aybarp*ff(0,1,1) +ayp*ff(0,1,2))+
1635  > axp*(aybarp*ff(0,2,1)+ayp*ff(0,2,2)))
1636 
1637  ans=ans+hhyi*(
1638  > bxbarp*(aybarp*ff(1,1,1) +ayp*ff(1,1,2))+
1639  > bxp*(aybarp*ff(1,2,1)+ayp*ff(1,2,2)))
1640 
1641  ans=ans+hhxi*(
1642  > axbarp*(bybarp*ff(2,1,1) +byp*ff(2,1,2))+
1643  > axp*(bybarp*ff(2,2,1)+byp*ff(2,2,2)))
1644 
1645  ans=ans+
1646  > bxbarp*(bybarp*ff(3,1,1) +byp*ff(3,1,2))+
1647  > bxp*(bybarp*ff(3,2,1)+byp*ff(3,2,2))
1648 
1649  endif
1650 
1651  end subroutine hm2d
1652 
1653  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17