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