V3FIT
r8pc2ev.f
1  subroutine r8pc2ev(xget,yget,x,nx,y,ny,ilinx,iliny,
2  > f,inf2,ict,fval,ier)
3 C
4 C evaluate a piecewise bilinear interpolant on a rectilinear
5 C grid -- this is C1 in both directions.
6 C
7 C this subroutine calls two subroutines:
8 C herm2xy -- find cell containing (xget,yget)
9 C pc2fcn -- evaluate interpolant function and (optionally) derivatives
10 C
11 C input arguments:
12 C ================
13 C
14 !============
15 ! idecl: explicitize implicit INTEGER declarations:
16  IMPLICIT NONE
17  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
18  INTEGER ny,inf2,nx
19 !============
20  real*8 xget,yget ! target of this interpolation
21  real*8 x(nx) ! ordered x grid
22  real*8 y(ny) ! ordered y grid
23  integer ilinx ! ilinx=1 => assume x evenly spaced
24  integer iliny ! iliny=1 => assume y evenly spaced
25 C
26  real*8 f(inf2,ny) ! function data
27 C
28 C f 2nd dimension inf2 must be .ge. nx
29 C contents of f:
30 C
31 C f(i,j) = f @ x(i),y(j)
32 C
33  integer ict(4) ! code specifying output desired
34 C
35 C ict(1)=1 -- return f (0, don't)
36 C ict(2)=1 -- return df/dx (0, don't)
37 C ict(3)=1 -- return df/dy (0, don't)
38 C ict(4)=1 -- return d2f/dxdy (0, don't)
39 C
40 C output arguments:
41 C =================
42 C
43  real*8 fval(*) ! output data
44  integer ier ! error code =0 ==> no error
45 C
46 C fval(1) receives the first output (depends on ict(...) spec)
47 C fval(2) receives the second output (depends on ict(...) spec)
48 C fval(3) receives the third output (depends on ict(...) spec)
49 C fval(4) receives the fourth output (depends on ict(...) spec)
50 C
51 C examples:
52 C on input ict = [1,1,1,1]
53 C on output fval= [f,df/dx,df/dy,d2f/dxdy]
54 C
55 C on input ict = [1,0,0,0]
56 C on output fval= [f] ... elements 2 & 3 & 4 never referenced
57 C
58 C on input ict = [0,1,1,0]
59 C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced
60 C
61 C on input ict = [0,0,1,0]
62 C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced.
63 C
64 C ier -- completion code: 0 means OK
65 C-------------------
66 C local:
67 C
68  integer i,j ! cell indices
69 C
70 C normalized displacement from (x(i),y(j)) corner of cell.
71 C xparam=0 @x(i) xparam=1 @x(i+1)
72 C yparam=0 @y(j) yparam=1 @y(j+1)
73 C
74  real*8 xparam,yparam
75 C
76 C cell dimensions and
77 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
78 C
79  real*8 hx,hy
80  real*8 hxi,hyi
81 C
82 C 0 .le. xparam .le. 1
83 C 0 .le. yparam .le. 1
84 C
85 C---------------------------------------------------------------------
86 C
87  call r8herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
88  > i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
89  if(ier.ne.0) return
90 c
91  call r8pc2fcn(ict,1,1,
92  > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,f,inf2,ny)
93 C
94  return
95  end
96 C---------------------------------------------------------------------
97 C evaluate piecewise bilinear function interpolation -- 2d fcn
98 C --vectorized-- dmc 10 Feb 1999
99 C
100  subroutine r8pc2fcn(ict,ivec,ivecd,
101  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
102  > fin,inf2,ny)
103 C
104 !============
105 ! idecl: explicitize implicit INTEGER declarations:
106  IMPLICIT NONE
107  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
108  INTEGER ny,inf2,i,j,iadr
109 !============
110 ! idecl: explicitize implicit REAL declarations:
111  real*8 xp,xpi,yp,ypi
112 !============
113  integer ict(4) ! requested output control
114  integer ivec ! vector length
115  integer ivecd ! vector dimension (1st dim of fval)
116 C
117  integer ii(ivec),jj(ivec) ! target cells (i,j)
118  real*8 xparam(ivec),yparam(ivec)
119  ! normalized displacements from (i,j) corners
120 C
121  real*8 hx(ivec),hy(ivec) ! grid spacing, and
122  real*8 hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
123  ! & 1/(y(j+1)-y(j))
124 C
125  real*8 fin(inf2,ny) ! interpolant data (cf "pc2ev")
126 C
127  real*8 fval(ivecd,*) ! output returned
128 C
129 C for detailed description of fin, ict and fval see subroutine
130 C pc2ev comments. Note ict is not vectorized; the same output
131 C is expected to be returned for all input vector data points.
132 C
133 C note that the index inputs ii,jj and parameter inputs
134 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
135 C output array fval has a vector ** 1st dimension ** whose
136 C size must be given as a separate argument
137 C
138 C to use this routine in scalar mode, pass in ivec=ivecd=1
139 C
140 C---------------
141 C
142  real*8 sum
143  integer v
144 C
145 C ...in x direction
146 C
147  do v=1,ivec
148  i=ii(v)
149  j=jj(v)
150 c
151  xp=xparam(v)
152  xpi=1.0_r8-xp
153 C
154 C ...in y direction
155 C
156  yp=yparam(v)
157  ypi=1.0_r8-yp
158 C
159  iadr=0
160 C
161 C get desired values:
162 C
163  if(ict(1).eq.1) then
164 C
165 C function value:
166 C
167  iadr=iadr+1
168  sum=ypi*(xpi*fin(i,j)+xp*fin(i+1,j))
169  > + yp*(xpi*fin(i,j+1)+xp*fin(i+1,j+1))
170 C
171  fval(v,iadr)=sum
172  endif
173 C
174  if(ict(2).eq.1) then
175 C
176 C df/dx:
177 C
178  iadr=iadr+1
179 C
180  sum=ypi*(fin(i+1,j)-fin(i,j))
181  > + yp*(fin(i+1,j+1)-fin(i,j+1))
182  fval(v,iadr)=sum*hxi(v)
183 C
184  endif
185 C
186  if(ict(3).eq.1) then
187 C
188 C df/dy:
189 C
190  iadr=iadr+1
191 C
192  sum=xpi*(fin(i,j+1)-fin(i,j))
193  > + xp*(fin(i+1,j+1)-fin(i+1,j))
194  fval(v,iadr)=sum*hyi(v)
195  endif
196 C
197  if(ict(4).eq.1) then
198 C
199 C d2f/dxdy:
200 C
201  iadr=iadr+1
202 C
203  sum=fin(i+1,j+1)-fin(i,j+1)-fin(i+1,j)+fin(i,j)
204  fval(v,iadr)=sum*hxi(v)*hyi(v)
205  endif
206 C
207  enddo ! vector loop
208 C
209  return
210  end
211 C---------------------------------------------------------------------
212 C evaluate piecewise bilinear function interpolation -- 2d fcn
213 C --vectorized-- dmc 10 Feb 1999
214 C --optimized for VARIATION along x axis ONLY--
215 C
216  subroutine r8pc2fcnx(ict,ivec,ivecd,
217  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
218  > fin,inf2,ny)
219 C
220 !============
221 ! idecl: explicitize implicit INTEGER declarations:
222  IMPLICIT NONE
223  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
224  INTEGER ny,inf2,j,i,iadr
225 !============
226 ! idecl: explicitize implicit REAL declarations:
227  real*8 yp,ypi,xp,xpi
228 !============
229  integer ict(4) ! requested output control
230  integer ivec ! vector length
231  integer ivecd ! vector dimension (1st dim of fval)
232 C
233  integer ii(ivec),jj ! target cells (i,j)
234  real*8 xparam(ivec),yparam
235  ! normalized displacements from (i,j) corners
236 C
237  real*8 hx(ivec),hy ! grid spacing, and
238  real*8 hxi(ivec),hyi ! inverse grid spacing 1/(x(i+1)-x(i))
239  ! & 1/(y(j+1)-y(j))
240 C
241  real*8 fin(inf2,ny) ! interpolant data (cf "pc2ev")
242 C
243  real*8 fval(ivecd,*) ! output returned
244 C
245 C for detailed description of fin, ict and fval see subroutine
246 C pc2ev comments. Note ict is not vectorized; the same output
247 C is expected to be returned for all input vector data points.
248 C
249 C note that the index inputs ii,jj and parameter inputs
250 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
251 C output array fval has a vector ** 1st dimension ** whose
252 C size must be given as a separate argument
253 C
254 C to use this routine in scalar mode, pass in ivec=ivecd=1
255 C
256 C---------------
257 C
258  real*8 sum
259  integer v
260 C
261 C ...in x direction
262 C
263  j=jj
264 C
265 C ...in y direction
266 C
267  yp=yparam
268  ypi=1.0_r8-yp
269 C
270  do v=1,ivec
271  i=ii(v)
272 c
273  xp=xparam(v)
274  xpi=1.0_r8-xp
275 C
276  iadr=0
277 C
278 C get desired values:
279 C
280  if(ict(1).eq.1) then
281 C
282 C function value:
283 C
284  iadr=iadr+1
285  sum=ypi*(xpi*fin(i,j)+xp*fin(i+1,j))
286  > + yp*(xpi*fin(i,j+1)+xp*fin(i+1,j+1))
287 C
288  fval(v,iadr)=sum
289  endif
290 C
291  if(ict(2).eq.1) then
292 C
293 C df/dx:
294 C
295  iadr=iadr+1
296 C
297  sum=ypi*(fin(i+1,j)-fin(i,j))
298  > + yp*(fin(i+1,j+1)-fin(i,j+1))
299  fval(v,iadr)=sum*hxi(v)
300 C
301  endif
302 C
303  if(ict(3).eq.1) then
304 C
305 C df/dy:
306 C
307  iadr=iadr+1
308 C
309  sum=xpi*(fin(i,j+1)-fin(i,j))
310  > + xp*(fin(i+1,j+1)-fin(i+1,j))
311  fval(v,iadr)=sum*hyi
312  endif
313 C
314  if(ict(4).eq.1) then
315 C
316 C d2f/dxdy:
317 C
318  iadr=iadr+1
319 C
320  sum=fin(i+1,j+1)-fin(i,j+1)-fin(i+1,j)+fin(i,j)
321  fval(v,iadr)=sum*hxi(v)*hyi
322  endif
323 C
324  enddo ! vector loop
325 C
326  return
327  end