V3FIT
r8vecintrp2d.f
1  subroutine r8vecintrp2d(ict,ivec,xvec,yvec,ivd,fval,
2  > nx,xpkg,ny,ypkg,jspline,fspl,icoeff,ixdim,iydim,
3  > iwarn,ier)
4 c
5 c
6 c vectorized hybrid spline evaluation routine -- 2d
7 c 1. call vectorized zone lookup routine
8 c 2. call vectorized hybrid spline evaluation routine
9 c
10 c--------------------------
11 c input:
12  IMPLICIT NONE
13  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
14  integer ict(6) ! selector:
15 c ict(1)=1 for f (don't evaluate if ict(1)=0)
16 c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
17 c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
18 c ict(4)=1 for d2f/dx2 (don't evaluate if ict(4)=0)
19 c ict(5)=1 for d2f/dy2 (don't evaluate if ict(5)=0)
20 c ict(6)=1 for d2f/dxdy (don't evaluate if ict(6)=0)
21 c
22 c ict(1)=3 followed by ict(2:5) = 1 or 0 allow 3rd derivatives to
23 c be selected: fxxx fxxy fxyy fyyy
24 c
25 c ict(1)=4 followed by ict(2:4) allow 4th derivatives to
26 c be selected: fxxxy fxxyy fxyyy; fxxxx=fyyyy=0
27 c
28 c ict(1)=5 followed by ict(2:4) allow 4th derivatives to
29 c be selected: fxxxyy fxxyyy
30 c
31 c ict(1)=6 specifies 6th derivative: fxxxyyy (step fcn)
32 c
33 c in hybrid spline evaluation, any derivative along a dimension
34 c with zonal interpolation (jspline(i)=-1) gives zero;
35 c
36 c piecewise linear and Hermite interpolation give zero if a 2nd or
37 c higher derivative is sought along any dimension.
38 c
39  integer ivec ! vector dimensioning
40 c
41 c ivec-- number of vector pts (spline values to look up)
42 c
43 c list of (x,y) pairs:
44 c
45  real*8 xvec(ivec) ! x-locations at which to evaluate
46  real*8 yvec(ivec) ! y-locations at which to evaluate
47 c
48  integer ivd ! 1st dimension of output array
49 c
50 c ivd -- 1st dimension of fval, .ge.ivec
51 c
52 c output:
53  real*8 fval(ivd,*) ! output array
54 c
55 c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
56 c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
57 c --etc--
58 c
59 c input:
60  integer nx,ny ! dimension of spline grids
61  REAL*8 xpkg(nx,4) ! x grid "package" (cf genxpkg)
62  real*8 ypkg(ny,4) ! y grid "package" (cf genxpkg)
63 c
64  integer :: jspline(2) ! interpolation method on each dim:
65 c jspline(1) for x; jspline(2) for y
66 c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline
67 c
68  integer icoeff ! #coeffs per data node
69  integer ixdim ! x dim for fspl
70  integer iydim ! y dim for fspl
71  real*8 fspl(icoeff,ixdim,iydim) ! hybrid spline coefficients
72 c
73 c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx
74 c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny
75 c
76 c output:
77 c condition codes, 0 = normal return
78  integer iwarn ! =1 if an x value was out of range
79  integer ier ! =1 if argument error detected
80 c
81 c---------------------------------------------------------------
82 c local variables and arrays
83 c
84  integer :: iwarn1,iwarn2
85 c
86  integer ix(ivec) ! zone indices {j}
87  REAL*8 dxn(ivec) ! normalized displacements w/in zones
88  real*8 hx(ivec) ! h(j) vector
89  real*8 hxi(ivec) ! 1/h(j) vector
90 c
91  integer iy(ivec) ! zone indices {j}
92  REAL*8 dyn(ivec) ! normalized displacements w/in zones
93  real*8 hy(ivec) ! h(j) vector
94  real*8 hyi(ivec) ! 1/h(j) vector
95 c
96 c---------------------------------------------------------------
97 c
98 c error checks
99 c
100  ier=0
101 c
102  if(nx.lt.2) then
103  write(6,*) .lt.' ?vecintrp2d: nx2: nx = ',nx
104  ier=1
105  endif
106 c
107  if(ny.lt.2) then
108  write(6,*) .lt.' ?vecintrp2d: ny2: ny = ',ny
109  ier=1
110  endif
111 c
112  call vecin2d_argchk('vecintrp2d',jspline,
113  > icoeff,nx,ny,ixdim,iydim,ier)
114 
115  if(ivec.le.0) then
116  write(6,*) .le.' ?vecintrp2d: vector dimension 0: ivec = ',
117  > ivec
118  ier=1
119  endif
120 c
121  if(ivd.lt.ivec) then
122  write(6,*)
123  > ' ?vecintrp2d: output vector dimension less than input ',
124  > 'vector dimension.'
125  write(6,*) ' ivec=',ivec,' ivd=',ivd
126  ier=1
127  endif
128 c
129  if(ier.ne.0) return
130 c
131 c vectorized lookup
132 c
133  ix=0
134  iy=0
135  call r8xlookup(ivec,xvec,nx,xpkg,2,ix,dxn,hx,hxi,iwarn1)
136  call r8xlookup(ivec,yvec,ny,ypkg,2,iy,dyn,hy,hyi,iwarn2)
137  iwarn=iwarn1+iwarn2
138 c
139 c vectorized evaluation
140 c
141  call r8fvintrp2d(ict,ivec,ivd,fval,ix,iy,dxn,dyn,
142  > hx,hxi,hy,hyi,jspline,fspl,icoeff,ixdim,iydim)
143 c
144  return
145  end