V3FIT
r8evintrp2d.f
1  subroutine r8evintrp2d(xget,yget,x,nx,y,ny,jspline,
2  > f,icoeff,ixdim,iydim,ict,fval,ier)
3 C
4 C
5 C evaluate a 2d Hybrid interpolant on a rectilinear grid
6 C
7 C this subroutine calls three subroutines:
8 C vecin2d_argchk -- error check
9 C herm2xy -- find cell containing (xget,yget)
10 C fvintrp2d -- evaluate interpolant function and (optionally) derivatives
11 C
12 C input arguments:
13 C ================
14 C
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  real*8 xget,yget ! target of this interpolation
18  integer nx,ny ! grid sizes
19  real*8 x(nx) ! ordered x grid
20  real*8 y(ny) ! ordered y grid
21 C
22  integer :: jspline(2) ! interpolation method for each
23 C dimension: jspline(1) for x; jspline(2) for y
24 C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline
25 C
26  integer :: icoeff ! no. of coefficients per data point
27 
28  integer :: ixdim,iydim ! dimensioning:
29  ! ixdim=nx-1 for zonal step in x; otherwise nx
30  ! iydim=ny-1 for zonal step in y; otherwise ny
31 
32  real*8 f(icoeff,ixdim,iydim) ! function data w/coefficients
33 C
34  integer ict(6) ! code specifying output desired
35 C
36 C Note on derivatives: for dimensions along which zonal step function
37 C interpolation is done, ANY derivative request returns ZERO.
38 C For dimensions along which piecewise linear or Hermite interpolation
39 C are done, more than one differentiation returns ZERO!
40 C
41 C Derivative controls are the same as for the compact 2d spline evaluation
42 C routine (evbicub):
43 C
44 C ict(1)=1 -- return f (0, don't)
45 C ict(2)=1 -- return df/dx (0, don't)
46 C ict(3)=1 -- return df/dy (0, don't)
47 C ict(4)=1 -- return d2f/dx2 (0, don't)
48 C ict(5)=1 -- return d2f/dy2 (0, don't)
49 C ict(6)=1 -- return d2f/dxdy (0, don't)
50 c the number of non zero values ict(1:6)
51 c determines the number of outputs...
52 c
53 c new dmc December 2005 -- access to higher derivatives (even if not
54 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
55 c if ict(1)=3 -- want 3rd derivatives
56 c ict(2)=1 for d3f/dx3
57 c ict(3)=1 for d3f/dx2dy
58 c ict(4)=1 for d3f/dxdy2
59 c ict(5)=1 for d3f/dy3
60 c number of non-zero values ict(2:5) gives no. of outputs
61 c if ict(1)=4 -- want 4th derivatives
62 c ict(2)=1 for d4f/dx3dy
63 c ict(3)=1 for d4f/dx2dy2
64 c ict(4)=1 for d4f/dxdy3
65 c number of non-zero values ict(2:4) gives no. of outputs
66 c if ict(1)=5 -- want 5th derivatives
67 c ict(2)=1 for d5f/dx3dy2
68 c ict(3)=1 for d5f/dx2dy3
69 c number of non-zero values ict(2:3) gives no. of outputs
70 c if ict(1)=6 -- want 6th derivatives
71 c d6f/dx3dy3 -- one value is returned.
72 C
73 C output arguments:
74 C =================
75 C
76  real*8 fval(6) ! output data
77  integer ier ! error code =0 ==> no error
78 C
79 C fval(1) receives the first output (depends on ict(...) spec)
80 C fval(2) receives the second output (depends on ict(...) spec)
81 C fval(3) receives the third output (depends on ict(...) spec)
82 C fval(4) receives the fourth output (depends on ict(...) spec)
83 C fval(5) receives the fifth output (depends on ict(...) spec)
84 C fval(6) receives the sixth output (depends on ict(...) spec)
85 C
86 C examples:
87 C on input ict = [1,1,1,0,0,1]
88 C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
89 C
90 C on input ict = [1,0,0,0,0,0]
91 C on output fval= [f] ... elements 2 -- 6 never referenced.
92 C
93 C on input ict = [0,0,0,1,1,0]
94 C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
95 C
96 C on input ict = [0,0,1,0,0,0]
97 C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
98 C
99 C ier -- completion code: 0 means OK
100 C-------------------
101 C local:
102 C
103  integer i,j ! cell indices
104 C
105 C normalized displacement from (x(i),y(j)) corner of cell.
106 C xparam=0 @x(i) xparam=1 @x(i+1)
107 C yparam=0 @y(j) yparam=1 @y(j+1)
108 C
109  real*8 xparam,yparam
110 C
111 C cell dimensions and
112 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
113 C
114  real*8 hx,hy
115  real*8 hxi,hyi
116 C
117 C 0 .le. xparam .le. 1
118 C 0 .le. yparam .le. 1
119 C
120 C ** the interface is very similar to herm2ev.for; can use herm2xy **
121 C---------------------------------------------------------------------
122 C
123  call vecin2d_argchk('evintrp2d',jspline,
124  > icoeff,nx,ny,ixdim,iydim,ier)
125  if(ier.ne.0) return
126 c
127  call r8herm2xy(xget,yget,x,nx,y,ny,0,0,
128  > i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
129  if(ier.ne.0) return
130 c
131  call r8fvintrp2d(ict,1,1,
132  > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,
133  > jspline,f,icoeff,ixdim,iydim)
134 C
135  return
136  end
137 C---------------------------------------------------------------------
138 C evaluate Hybrid function interpolation -- 2d fcn
139 C NO ERROR CHECKING
140 C
141  subroutine r8fvintrp2d(ict,ivec,ivecd,
142  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
143  > jspline,fin,icoeff,ixdim,iydim)
144 C
145 C
146  IMPLICIT NONE
147  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
148  integer ict(6) ! requested output control
149  integer ivec ! vector length
150  integer ivecd ! vector dimension (1st dim of fval)
151 C
152  integer ii(ivec),jj(ivec) ! target cells (i,j)
153  real*8 xparam(ivec),yparam(ivec)
154  ! normalized displacements from (i,j) corners
155 C
156  real*8 hx(ivec),hy(ivec) ! grid spacing, and
157  real*8 hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
158  ! & 1/(y(j+1)-y(j))
159 C
160  integer :: jspline(2) ! interpolation type control
161  integer :: icoeff ! coefficient dimension
162  integer :: ixdim,iydim ! x & y dimensions
163  real*8 fin(icoeff,ixdim,iydim) ! data on which to interpolate...
164 C
165  real*8 fval(ivecd,6) ! output returned
166 C
167 C for detailed description of fin, ict and fval see subroutine
168 C evintrp2d comments. Note ict is not vectorized; the same output
169 C is expected to be returned for all input vector data points.
170 C
171 C note that the index inputs ii,jj and parameter inputs
172 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
173 C output array fval has a vector ** 1st dimension ** whose
174 C size must be given as a separate argument
175 C
176 C to use this routine in scalar mode, pass in ivec=ivecd=1
177 C
178 C---------------
179 C local:
180 C
181  integer :: i,j,i1,i2,zonrank,linrank,cubrank
182  integer :: imaxx,imaxy,imaxdlin,imaxdcub
183  logical :: splin_flag
184 C
185  integer :: inum,indx
186  integer :: iderivs(2,6),ict4(4),maxd_lin(4),maxd_cub(4)
187 C
188  integer :: idcub,idlin
189  real*8 :: xp,dxlin,dxlini,h,hi,ans1,ans2
190  real*8 :: f22(2,2)
191 C
192  real*8, parameter :: one = 1.0_r8
193 C---------------
194 C
195  linrank=0
196  cubrank=0
197  zonrank=0
198 
199  do i=1,2
200  if(jspline(i).eq.-1) then
201  zonrank = zonrank + 1
202 
203  else if(jspline(i).eq.0) then
204  linrank = linrank + 1
205 
206  else
207  splin_flag = jspline(i).eq.2
208  cubrank = cubrank + 1
209  endif
210  enddo
211 C
212  if(cubrank.eq.2) then
213  if(splin_flag) then
214  ! bicubic spline
215  call r8fvbicub(ict,ivec,ivecd,
216  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
217  > fin,ixdim,iydim)
218 
219  else
220  ! bicubic Hermite -- translate derivative code
221  call dtrans
222  do i=1,inum
223  if(maxval(iderivs(1:2,i)).le.1) then
224  ict4=0
225  indx=2*iderivs(2,i) + iderivs(1,i) + 1
226  ict4(indx)=1
227  call r8herm2fcn(ict4,ivec,ivecd,
228  > fval(1,i),ii,jj,xparam,yparam,hx,hxi,hy,hyi,
229  > fin,ixdim,iydim)
230 
231  else
232  fval(1:ivec,i) = 0
233  endif
234  enddo
235  endif
236 
237  return
238  endif
239 
240  call dtrans
241 
242  if(cubrank.eq.0) then
243  if(jspline(1).eq.0) then
244  imaxx=1
245  else
246  imaxx=0
247  endif
248 
249  if(jspline(2).eq.0) then
250  imaxy=1
251  else
252  imaxy=0
253  endif
254 
255  do i=1,inum
256  if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)) then
257  if(linrank.eq.2) then
258  ! bilinear interpolation
259  ict4=0
260  indx=2*iderivs(2,i) + iderivs(1,i) + 1
261  ict4(indx)=1
262  call r8pc2fcn(ict4,ivec,ivecd,
263  > fval(1,i),ii,jj,xparam,yparam,hx,hxi,hy,hyi,
264  > fin,ixdim,iydim)
265  else if(zonrank.eq.2) then
266  ! 2d step function
267  do j=1,ivec
268  fval(j,i)=fin(1,ii(j),jj(j))
269  enddo
270  else
271  ! linear-zone hybrid
272  if(jspline(1).eq.0) then
273  ! linear in x, zonal in y
274  do j=1,ivec
275  if(iderivs(1,i).eq.0) then
276  fval(j,i)=(one-xparam(j))*fin(1,ii(j),jj(j))
277  > +xparam(j)*fin(1,ii(j)+1,jj(j))
278  else
279  fval(j,i)=(fin(1,ii(j)+1,jj(j)) -
280  > fin(1,ii(j),jj(j)))*hxi(j)
281  endif
282  enddo
283  else
284  ! linear in y, zonal in x
285  do j=1,ivec
286  if(iderivs(2,i).eq.0) then
287  fval(j,i)=(one-yparam(j))*fin(1,ii(j),jj(j))
288  > +yparam(j)*fin(1,ii(j),jj(j)+1)
289  else
290  fval(j,i)=(fin(1,ii(j),jj(j)+1) -
291  > fin(1,ii(j),jj(j)))*hyi(j)
292  endif
293  enddo
294  endif
295  endif
296  else
297  fval(1:ivec,i) = 0
298  endif
299  enddo
300 
301  return
302  endif
303 
304  ! OK -- Hybrid: Hermite or Spline in one dimension
305  ! pclin or step in the other dimension
306 
307  if(splin_flag) then
308  imaxdcub=3
309  else
310  imaxdcub=1
311  endif
312 
313  if(linrank.eq.1) then
314  imaxdlin=1
315  else
316  imaxdlin=0
317  endif
318 
319  do i=1,inum
320  if(jspline(1).ge.1) then
321  idcub=iderivs(1,i)
322  idlin=iderivs(2,i)
323  else
324  idcub=iderivs(2,i)
325  idlin=iderivs(1,i)
326  endif
327 
328  if((maxd_lin(i).le.imaxdlin).and.
329  > (maxd_cub(i).le.imaxdcub)) then
330  if(linrank.eq.1) then
331  ! spline or Hermite in 1 dimension; pclin in other.
332  if(jspline(1).ge.1) then
333  do j=1,ivec
334  xp=xparam(j)
335  dxlin=yparam(j)
336  dxlini=1.0_r8-dxlin
337  h=hx(j)
338  hi=hxi(j)
339  i1=ii(j)
340  i2=jj(j)
341  f22(1:2,1)=fin(1:2,i1,i2)
342  f22(1:2,2)=fin(1:2,i1+1,i2)
343  if(splin_flag) then
344  call sp1d(ans1)
345  else
346  call hm1d(ans1)
347  endif
348  f22(1:2,1)=fin(1:2,i1,i2+1)
349  f22(1:2,2)=fin(1:2,i1+1,i2+1)
350  if(splin_flag) then
351  call sp1d(ans2)
352  else
353  call hm1d(ans2)
354  endif
355  if(idlin.eq.1) then
356  fval(j,i)=(ans2-ans1)*hyi(j)
357  else
358  fval(j,i)=ans2*dxlin + ans1*dxlini
359  endif
360  enddo
361 
362  else
363  do j=1,ivec
364  xp=yparam(j)
365  dxlin=xparam(j)
366  dxlini=1.0_r8-dxlin
367  h=hy(j)
368  hi=hyi(j)
369  i1=ii(j)
370  i2=jj(j)
371  f22(1:2,1)=fin(1:2,i1,i2)
372  f22(1:2,2)=fin(1:2,i1,i2+1)
373  if(splin_flag) then
374  call sp1d(ans1)
375  else
376  call hm1d(ans1)
377  endif
378  f22(1:2,1)=fin(1:2,i1+1,i2)
379  f22(1:2,2)=fin(1:2,i1+1,i2+1)
380  if(splin_flag) then
381  call sp1d(ans2)
382  else
383  call hm1d(ans2)
384  endif
385  if(idlin.eq.1) then
386  fval(j,i)=(ans2-ans1)*hxi(j)
387  else
388  fval(j,i)=ans2*dxlin + ans1*dxlini
389  endif
390  enddo
391  endif
392 
393  else
394  ! spline or Hermite in 1 dimension; step fcn in other.
395  if(jspline(1).ge.1) then
396  do j=1,ivec
397  xp=xparam(j)
398  h=hx(j)
399  hi=hxi(j)
400  i1=ii(j)
401  i2=jj(j)
402  f22(1:2,1)=fin(1:2,i1,i2)
403  f22(1:2,2)=fin(1:2,i1+1,i2)
404  if(splin_flag) then
405  call sp1d(ans1)
406  else
407  call hm1d(ans1)
408  endif
409  fval(j,i)=ans1
410  enddo
411  else
412  do j=1,ivec
413  xp=yparam(j)
414  h=hy(j)
415  hi=hyi(j)
416  i1=ii(j)
417  i2=jj(j)
418  f22(1:2,1)=fin(1:2,i1,i2)
419  f22(1:2,2)=fin(1:2,i1,i2+1)
420  if(splin_flag) then
421  call sp1d(ans1)
422  else
423  call hm1d(ans1)
424  endif
425  fval(j,i)=ans1
426  enddo
427  endif
428  endif
429  else
430  fval(1:ivec,i) = 0
431  endif
432  enddo
433 
434  return
435 
436  contains
437  subroutine dtrans
438 
439  ! convert ict(...) codes into intermediate coding:
440  ! iderivs(i,j) = number of differentiations in coordinate i
441  ! of the j'th value evaluated
442 
443  ! maxd_lin(j) = max # of differentiations for a linear dimension
444  ! maxd_cub(j) = max # of differentiations for a spline dimension
445 
446  integer :: i
447 
448  inum=0 ! actual number of vectors to be evaluated (summed here)
449 
450  if(ict(1).le.(2)) then
451  if(ict(1).eq.1) then
452  call add1(0,0)
453  endif
454  if(ict(2).eq.1) then
455  call add1(1,0)
456  endif
457  if(ict(3).eq.1) then
458  call add1(0,1)
459  endif
460  if(ict(4).eq.1) then
461  call add1(2,0)
462  endif
463  if(ict(5).eq.1) then
464  call add1(0,2)
465  endif
466  if(ict(6).eq.1) then
467  call add1(1,1)
468  endif
469 
470  else if(ict(1).eq.3) then
471  if(ict(2).eq.1) then
472  call add1(3,0)
473  endif
474  if(ict(3).eq.1) then
475  call add1(2,1)
476  endif
477  if(ict(4).eq.1) then
478  call add1(1,2)
479  endif
480  if(ict(5).eq.1) then
481  call add1(0,3)
482  endif
483 
484  else if(ict(1).eq.4) then
485  if(ict(2).eq.1) then
486  call add1(3,1)
487  endif
488  if(ict(3).eq.1) then
489  call add1(2,2)
490  endif
491  if(ict(4).eq.1) then
492  call add1(1,3)
493  endif
494 
495  else if(ict(1).eq.5) then
496  if(ict(2).eq.1) then
497  call add1(3,2)
498  endif
499  if(ict(3).eq.1) then
500  call add1(2,3)
501  endif
502 
503  else if(ict(1).eq.6) then
504  if(ict(2).eq.1) then
505  call add1(3,3)
506  endif
507  endif
508 
509  end subroutine dtrans
510 
511  subroutine add1(idx,idy)
512  ! insert record of derivative d[idx+idy]f/dx[idx]dy[idy]
513 
514  integer, intent(in) :: idx,idy
515 
516  inum=inum+1
517 
518  iderivs(1,inum)=idx
519  if(jspline(1).le.0) then
520  maxd_lin(inum)=idx
521  maxd_cub(inum)=0
522  else
523  maxd_lin(inum)=0
524  maxd_cub(inum)=idx
525  endif
526 
527  iderivs(2,inum)=idy
528  if(jspline(2).le.0) then
529  maxd_lin(inum)=max(maxd_lin(inum),idy)
530  else
531  maxd_cub(inum)=max(maxd_cub(inum),idy)
532  endif
533 
534  end subroutine add1
535 
536 C===========================================================================
537  subroutine sp1d(ans)
538  real*8, intent(out) :: ans
539 
540 C----------------------
541 C contained 1d spline evaluation
542 C----------------------
543 
544  real*8 :: xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
545 C
546  real*8 sixth
547  data sixth/0.166666666666666667_r8/
548 C----------------------
549 
550  if(idcub.eq.0) then
551  ! value
552 
553  xpi=1.0_r8-xp
554  xp2=xp*xp
555  xpi2=xpi*xpi
556 
557  cx=xp*(xp2-1.0_r8)
558  cxi=xpi*(xpi2-1.0_r8)
559  hx2=h*h
560 
561  ans=xpi*f22(1,1) +xp*f22(1,2)
562  ans=ans+sixth*hx2*(cxi*f22(2,1) +cx*f22(2,2))
563 
564  else if(idcub.eq.1) then
565  ! 1st derivative
566 
567  xpi=1.0_r8-xp
568  xp2=xp*xp
569  xpi2=xpi*xpi
570 
571  cxd=3.0_r8*xp2-1.0_r8
572  cxdi=-3.0_r8*xpi2+1.0_r8
573 
574  ans=hi*(f22(1,2)-f22(1,1))
575  ans=ans+sixth*h*(cxdi*f22(2,1) +cxd*f22(2,2))
576 
577  else if(idcub.eq.2) then
578  ! 2nd derivative
579 
580  xpi=1.0_r8-xp
581  ans=xpi*f22(2,1) + xp*f22(2,2)
582 
583  else
584  ! 3rd derivative
585 
586  ans = hi*(f22(2,2)-f22(2,1))
587 
588  endif
589 
590  end subroutine sp1d
591 
592 C===========================================================================
593  subroutine hm1d(ans)
594  real*8, intent(out) :: ans
595 
596 C----------------------
597 C contained 1d hermite evaluation
598 C----------------------
599 
600  real*8 :: xpi,xp2,xpi2,ax,axbar,bx,bxbar
601  real*8 :: axp,axbarp,bxp,bxbarp
602 
603  xpi=1.0_r8-xp
604  xp2=xp*xp
605  xpi2=xpi*xpi
606 
607  if(idcub.eq.0) then
608  ! value
609 
610  ax=xp2*(3.0_r8-2.0_r8*xp)
611  axbar=1.0_r8-ax
612  bx=-xp2*xpi
613  bxbar=xpi2*xp
614 
615  ans=axbar*f22(1,1) + ax*f22(1,2)
616  ans=ans + h*(bxbar*f22(2,1) + bx*f22(2,2))
617 
618  else
619  ! 1st derivative
620 
621  axp=6.0_r8*xp*xpi
622  axbarp=-axp
623  bxp=xp*(3.0_r8*xp-2.0_r8)
624  bxbarp=xpi*(3.0_r8*xpi-2.0_r8)
625 
626  ans=hi*(axbarp*f22(1,1) +axp*f22(1,2))
627  ans=ans + bxbarp*f22(2,1) + bxp*f22(2,2)
628 
629  endif
630 
631  end subroutine hm1d
632 
633  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17