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