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