V3FIT
vectricub.f
1  subroutine vectricub(ict,ivec,xvec,yvec,zvec,ivd,fval,
2  > nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5,
3  > iwarn,ier)
4 c
5 c vectorized spline evaluation routine -- 3d *compact* spline
6 c 1. call vectorized zone lookup routine
7 c 2. call vectorized spline evaluation routine
8 c
9 c--------------------------
10 c input:
11  integer ict(10) ! selector:
12 c ict(1)=1 for f (don't evaluate if ict(1)=0)
13 c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
14 c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
15 c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
16 c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
17 c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
18 c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
19 c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
20 c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
21 c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
22 c
23  integer ivec ! vector dimensioning
24 c
25 c ivec-- number of vector pts (spline values to look up)
26 c
27 c list of (x,y,z) triples:
28 c
29  real xvec(ivec) ! x-locations at which to evaluate
30  real yvec(ivec) ! y-locations at which to evaluate
31  real zvec(ivec) ! z-locations at which to evaluate
32 c
33  integer ivd ! 1st dimension of output array
34 c
35 c ivd -- 1st dimension of fval, .ge.ivec
36 c
37 c output:
38  real fval(ivd,*) ! output array
39 c
40 c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
41 c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
42 c --etc--
43 c
44 c input:
45  integer nx,ny,nz ! dimension of spline grids
46  real xpkg(nx,4) ! x grid "package" (cf genxpkg)
47  real ypkg(ny,4) ! y grid "package" (cf genxpkg)
48  real zpkg(nz,4) ! z grid "package" (cf genxpkg)
49  integer inf4 ! fspl 4th array dimension, .ge.nx
50  integer inf5 ! fspl 5th array dimension, .ge.ny
51  real fspl(8,inf4,inf5,nz) ! (compact) spline coefficients
52 c
53 c output:
54 c condition codes, 0 = normal return
55  integer iwarn ! =1 if an x value was out of range
56  integer ier ! =1 if argument error detected
57 c
58 c---------------------------------------------------------------
59 c local arrays
60 c
61  integer, dimension(:), allocatable :: ix,iy,iz
62  real, dimension(:), allocatable :: dxn,dyn,dzn
63  real, dimension(:), allocatable :: hx,hxi,hy,hyi,hz,hzi
64 c
65 c---------------------------------------------------------------
66 c
67 c error checks
68 c
69  ier=0
70 c
71  if(nx.lt.2) then
72  write(6,*) .lt.' ?vectricub: nx2: nx = ',nx
73  ier=1
74  endif
75 c
76  if(ny.lt.2) then
77  write(6,*) .lt.' ?vectricub: ny2: ny = ',ny
78  ier=1
79  endif
80 c
81  if(nz.lt.2) then
82  write(6,*) .lt.' ?vectricub: nz2: nz = ',nz
83  ier=1
84  endif
85 c
86  if(ivec.le.0) then
87  write(6,*) .le.' ?vectricub: vector dimension 0: ivec = ',
88  > ivec
89  ier=1
90  endif
91 c
92  if(ivd.lt.ivec) then
93  write(6,*)
94  > ' ?vectricub: output vector dimension less than input ',
95  > 'vector dimension.'
96  write(6,*) ' ivec=',ivec,' ivd=',ivd
97  ier=1
98  endif
99 c
100  if(ier.ne.0) return
101 c
102  allocate(ix(ivec), iy(ivec), iz(ivec),
103  > dxn(ivec), dyn(ivec), dzn(ivec),
104  > hx(ivec), hy(ivec), hz(ivec),
105  > hxi(ivec), hyi(ivec), hzi(ivec), stat=ier)
106 c
107  if(ier.ne.0) then
108  write(6,*)
109  > ' ?vectricub: memory allocation failure.'
110  ier=99
111  endif
112 c
113  if(ier.ne.0) return
114 c
115 c vectorized lookup
116 c
117  ix=0; iy=0; iz=0
118  call xlookup(ivec,xvec,nx,xpkg,2,ix,dxn,hx,hxi,iwarn1)
119  call xlookup(ivec,yvec,ny,ypkg,2,iy,dyn,hy,hyi,iwarn2)
120  call xlookup(ivec,zvec,nz,zpkg,2,iz,dzn,hz,hzi,iwarn3)
121  iwarn=max(iwarn1,iwarn2,iwarn3)
122 c
123 c vectorized evaluation
124 c
125  call fvtricub(ict,ivec,ivd,fval,ix,iy,iz,dxn,dyn,dzn,
126  > hx,hxi,hy,hyi,hz,hzi,fspl,inf4,inf5,nz)
127 c
128  deallocate(ix,iy,iz,dxn,dyn,dzn,hx,hy,hz,hxi,hyi,hzi)
129 c
130  return
131  end