V3FIT
r8vecintrp3d.f
1  subroutine r8vecintrp3d(ict,ivec,xvec,yvec,zvec,ivd,fval,
2  > nx,xpkg,ny,ypkg,nz,zpkg,jspline,fspl,icoeff,ixdim,iydim,izdim,
3  > iwarn,ier)
4 c
5 c
6 c vectorized hybrid spline evaluation routine -- 3d
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 !============
13 ! idecl: explicitize implicit REAL declarations:
14  IMPLICIT NONE
15  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
16  real*8 stat
17 !============
18  integer ict(10) ! selector:
19 c ict(1)=1 for f (don't evaluate if ict(1)=0)
20 c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
21 c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
22 c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
23 c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
24 c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
25 c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
26 c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
27 c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
28 c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
29 c
30 c higher derivatives can be selected by setting ict(1)=3,-3,4,-4,...
31 c see fvtricub.for comments.
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,z) triples:
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  real*8 zvec(ivec) ! z-locations at which to evaluate
48 c
49  integer ivd ! 1st dimension of output array
50 c
51 c ivd -- 1st dimension of fval, .ge.ivec
52 c
53 c output:
54  real*8 fval(ivd,*) ! output array
55 c
56 c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
57 c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
58 c --etc--
59 c
60 c input:
61  integer nx,ny,nz ! dimension of spline grids
62  REAL*8 xpkg(nx,4) ! x grid "package" (cf genxpkg)
63  real*8 ypkg(ny,4) ! y grid "package" (cf genxpkg)
64  real*8 zpkg(nz,4) ! z grid "package" (cf genxpkg)
65 c
66  integer :: jspline(3) ! interpolation method on each dim:
67 c jspline(1) for x; jspline(2) for y; jspline(3) for z
68 c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline
69 c
70  integer icoeff ! #coeffs per data node
71  integer ixdim ! x dim for fspl
72  integer iydim ! y dim for fspl
73  integer izdim ! z dim for fspl
74  real*8 fspl(icoeff,ixdim,iydim,izdim) ! hybrid spline coefficients
75 c
76 c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx
77 c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny
78 c izdim=nz-1 for zonal step function along z dimension; o.w. expect izdim=nz
79 c
80 c output:
81 c condition codes, 0 = normal return
82  integer iwarn ! =1 if an x value was out of range
83  integer ier ! =1 if argument error detected
84 c
85 c---------------------------------------------------------------
86 c local arrays
87 c
88  integer :: iwarn1,iwarn2,iwarn3
89 c
90  integer, dimension(:), allocatable :: ix,iy,iz
91  REAL*8, dimension(:), allocatable :: dxn,dyn,dzn
92  real*8, dimension(:), allocatable :: hx,hxi,hy,hyi,hz,hzi
93 c
94 c---------------------------------------------------------------
95 c
96 c error checks
97 c
98  ier=0
99 c
100  if(nx.lt.2) then
101  write(6,*) .lt.' ?vecintrp3d: nx2: nx = ',nx
102  ier=1
103  endif
104 c
105  if(ny.lt.2) then
106  write(6,*) .lt.' ?vecintrp3d: ny2: ny = ',ny
107  ier=1
108  endif
109 c
110  if(nz.lt.2) then
111  write(6,*) .lt.' ?vecintrp3d: nz2: nz = ',nz
112  ier=1
113  endif
114 c
115  call vecin3d_argchk('vecintrp3d',jspline,
116  > icoeff,nx,ny,nz,ixdim,iydim,izdim,ier)
117 c
118  if(ivec.le.0) then
119  write(6,*) .le.' ?vecintrp3d: vector dimension 0: ivec = ',
120  > ivec
121  ier=1
122  endif
123 c
124  if(ivd.lt.ivec) then
125  write(6,*)
126  > ' ?vecintrp3d: output vector dimension less than input ',
127  > 'vector dimension.'
128  write(6,*) ' ivec=',ivec,' ivd=',ivd
129  ier=1
130  endif
131 c
132  if(ier.ne.0) return
133 c
134  allocate(ix(ivec), iy(ivec), iz(ivec),
135  > dxn(ivec), dyn(ivec), dzn(ivec),
136  > hx(ivec), hy(ivec), hz(ivec),
137  > hxi(ivec), hyi(ivec), hzi(ivec), stat=ier)
138 c
139  if(ier.ne.0) then
140  write(6,*)
141  > ' ?vecintrp3d: memory allocation failure.'
142  ier=99
143  endif
144 c
145  if(ier.ne.0) return
146 c
147 c vectorized lookup
148 c
149  ix=0; iy=0; iz=0
150  call r8xlookup(ivec,xvec,nx,xpkg,2,ix,dxn,hx,hxi,iwarn1)
151  call r8xlookup(ivec,yvec,ny,ypkg,2,iy,dyn,hy,hyi,iwarn2)
152  call r8xlookup(ivec,zvec,nz,zpkg,2,iz,dzn,hz,hzi,iwarn3)
153  iwarn=max(iwarn1,iwarn2,iwarn3)
154 c
155 c vectorized evaluation
156 c
157  call r8fvintrp3d(ict,ivec,ivd,fval,ix,iy,iz,dxn,dyn,dzn,
158  > hx,hxi,hy,hyi,hz,hzi,jspline,fspl,icoeff,ixdim,iydim,izdim)
159 c
160  deallocate(ix,iy,iz,dxn,dyn,dzn,hx,hy,hz,hxi,hyi,hzi)
161 c
162  return
163  end