| 2 |
michaesp |
1 |
c *************************************************************
|
|
|
2 |
c * This package provides a general interpolaton routine *
|
|
|
3 |
c *************************************************************
|
|
|
4 |
|
|
|
5 |
c The main interface routines are:
|
|
|
6 |
c get_index3,4 : get the grid indices for interpolation
|
|
|
7 |
c int_index3,4 : interpolate to the grid position
|
|
|
8 |
|
|
|
9 |
c -------------------------------------------------------------
|
|
|
10 |
c Get index in grid space for interpolation
|
|
|
11 |
c -------------------------------------------------------------
|
|
|
12 |
|
|
|
13 |
subroutine get_index4 (rid,rjd,rkd,xpo,ypo,ppo,rtp,
|
|
|
14 |
> vert0,vert1,surf0,surf1,mode,
|
|
|
15 |
> nx,ny,nz,lonw,lats,dlon,dlat,misdat)
|
|
|
16 |
|
|
|
17 |
c Purpose:
|
|
|
18 |
c This subroutine determines the indices (rid,rjd,rkd) in grid
|
|
|
19 |
c space for a point in physical space (xpo,ypo,ppo). The
|
|
|
20 |
c horizontal grid is specified by the south-west point (lats,lonw)
|
|
|
21 |
c and the grid spacing (dlat,dlon). The vertical grid is given
|
|
|
22 |
c by <vert(n1,n2,n3)>. The lower boundary (typicall surface
|
|
|
23 |
c pressure) is given by <surf(n1,n2)>.
|
|
|
24 |
c Arguments:
|
|
|
25 |
c rid,rjd,rkd real output grid location to be interpolated to
|
|
|
26 |
c xpo,ypo,ppo real input physical coordinates
|
|
|
27 |
c rtp real input relative time position (0=beginning, 1=end)
|
|
|
28 |
c n1,n2,n3 int input grid dimensions in x-, y- and p-direction
|
|
|
29 |
c lats,lonw real input south and west boundary of grid space
|
|
|
30 |
c vert real input vertical coordinate grid
|
|
|
31 |
c surf real input lower boundary (surface pressure)
|
|
|
32 |
c mode int input direction of vertical axis (p=1,th=-1)
|
|
|
33 |
c 1: linear, 1 -> nz (th)
|
|
|
34 |
c 2: linear, nz -> 1 (pv)
|
|
|
35 |
c 3: binary (z)
|
|
|
36 |
c 4: binary (p)
|
|
|
37 |
|
|
|
38 |
implicit none
|
|
|
39 |
|
|
|
40 |
c Declartion of function parameters
|
|
|
41 |
integer nx,ny,nz
|
|
|
42 |
real xpo,ypo,ppo,rtp
|
|
|
43 |
real vert0(nx*ny*nz),vert1(nx*ny*nz)
|
|
|
44 |
real surf0(nx*ny) ,surf1(nx*ny*nz)
|
|
|
45 |
real rid,rjd,rkd
|
|
|
46 |
real dlat,dlon,lats,lonw
|
|
|
47 |
real misdat
|
|
|
48 |
integer mode
|
|
|
49 |
|
|
|
50 |
c Set numerical parameters
|
|
|
51 |
real eps
|
|
|
52 |
parameter (eps=1.e-8)
|
|
|
53 |
|
|
|
54 |
c Auxiliary variables
|
|
|
55 |
real rid0,rjd0,rkd0,rid1,rjd1,rkd1
|
|
|
56 |
|
|
|
57 |
c Externals
|
|
|
58 |
real int_time
|
|
|
59 |
external int_time
|
|
|
60 |
|
|
|
61 |
c Get the inidices
|
|
|
62 |
if (abs(rtp).lt.eps) then
|
|
|
63 |
call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
|
|
|
64 |
> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)
|
|
|
65 |
elseif (abs(rtp-1.).lt.eps) then
|
|
|
66 |
call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
|
|
|
67 |
> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)
|
|
|
68 |
else
|
|
|
69 |
call get_index3 (rid0,rjd0,rkd0,xpo,ypo,ppo,mode,
|
|
|
70 |
> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)
|
|
|
71 |
call get_index3 (rid1,rjd1,rkd1,xpo,ypo,ppo,mode,
|
|
|
72 |
> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)
|
|
|
73 |
rid = int_time (rid0,rid1,rtp,misdat)
|
|
|
74 |
rjd = int_time (rjd0,rjd1,rtp,misdat)
|
|
|
75 |
rkd = int_time (rkd0,rkd1,rtp,misdat)
|
|
|
76 |
|
|
|
77 |
endif
|
|
|
78 |
|
|
|
79 |
end
|
|
|
80 |
|
|
|
81 |
c -------------------------------------------------------------
|
|
|
82 |
c Interpolate to an arbitrary position in grid space and time
|
|
|
83 |
c -------------------------------------------------------------
|
|
|
84 |
|
|
|
85 |
real function int_index4 (ar0,ar1,n1,n2,n3,rid,rjd,rkd,rtp,misdat)
|
|
|
86 |
|
|
|
87 |
c Purpose:
|
|
|
88 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
89 |
c location within the grid including a linear time-interpolation.
|
|
|
90 |
c Arguments:
|
|
|
91 |
c rid,rjd,rkd real output grid location to be interpolated to
|
|
|
92 |
c xpo,ypo,ppo real input physical coordinates
|
|
|
93 |
c n1,n2,n3 int input grid dimensions in x-, y- and p-direction
|
|
|
94 |
c lats,lonw real input south and west boundary of grid space
|
|
|
95 |
c vert real input vertical coordinate grid
|
|
|
96 |
c surf real input lower boundary (surface pressure)
|
|
|
97 |
|
|
|
98 |
implicit none
|
|
|
99 |
|
|
|
100 |
c Declartion of function parameters
|
|
|
101 |
integer n1,n2,n3
|
|
|
102 |
real ar0(n1*n2*n3),ar1(n1*n2*n3)
|
|
|
103 |
real rid,rjd,rkd
|
|
|
104 |
real rtp
|
|
|
105 |
real misdat
|
|
|
106 |
|
|
|
107 |
c Set numerical parameters
|
|
|
108 |
real eps
|
|
|
109 |
parameter (eps=1.e-8)
|
|
|
110 |
|
|
|
111 |
c Externals
|
|
|
112 |
real int_index3,int_time
|
|
|
113 |
external int_index3,int_time
|
|
|
114 |
|
|
|
115 |
c Auxiliary variables
|
|
|
116 |
real val0,val1,val
|
|
|
117 |
|
|
|
118 |
c Do the 3d-interpolation
|
|
|
119 |
if (abs(rtp).lt.eps) then
|
|
|
120 |
val = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
121 |
elseif (abs(rtp-1.).lt.eps) then
|
|
|
122 |
val = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
123 |
else
|
|
|
124 |
val0 = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
125 |
val1 = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
126 |
val = int_time (val0,val1,rtp,misdat)
|
|
|
127 |
endif
|
|
|
128 |
|
|
|
129 |
c Return value
|
|
|
130 |
int_index4 = val
|
|
|
131 |
|
|
|
132 |
return
|
|
|
133 |
end
|
|
|
134 |
|
|
|
135 |
|
|
|
136 |
c -------------------------------------------------------------
|
|
|
137 |
c Interpolate to an arbitrary position in grid space
|
|
|
138 |
c -------------------------------------------------------------
|
|
|
139 |
|
|
|
140 |
real function int_index3 (ar,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
141 |
|
|
|
142 |
c Purpose:
|
|
|
143 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
144 |
c location within the grid. The interpolation includes the
|
|
|
145 |
c testing of the missing data flag 'misdat'. If one dimension
|
|
|
146 |
c is 1, a 2d-interpolation is performed; if two dimensions
|
|
|
147 |
c are 1, it is a 1d-interpolation; if all three dimensions are
|
|
|
148 |
c 1, no interpolation is performed and the input value is
|
|
|
149 |
c returned.
|
|
|
150 |
c Arguments:
|
|
|
151 |
c ar real input input data array
|
|
|
152 |
c n1,n2,n3 int input dimensions of ar
|
|
|
153 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
154 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
155 |
|
|
|
156 |
implicit none
|
|
|
157 |
|
|
|
158 |
c Declartion of function parameters
|
|
|
159 |
integer n1,n2,n3
|
|
|
160 |
real ar(n1*n2*n3)
|
|
|
161 |
real rid,rjd,rkd
|
|
|
162 |
real misdat
|
|
|
163 |
|
|
|
164 |
c Set numerical parameters
|
|
|
165 |
real eps
|
|
|
166 |
parameter (eps=1.e-8)
|
|
|
167 |
|
|
|
168 |
c Local variables
|
|
|
169 |
integer i,j,k,ip1,jp1,kp1
|
|
|
170 |
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k
|
|
|
171 |
real ri,rj,rk
|
|
|
172 |
real val000,val001,val010,val011,val100,val101,val110,val111
|
|
|
173 |
real frc000,frc001,frc010,frc011,frc100,frc101,frc110,frc111
|
|
|
174 |
real frc
|
|
|
175 |
real mdv
|
|
|
176 |
real val
|
|
|
177 |
|
|
|
178 |
c Elementary test for dimensions
|
|
|
179 |
if ( (n1.lt.1).or.(n2.lt.1).or.(n3.lt.1) ) then
|
|
|
180 |
print*,'Invalid grid dimensions ',n1,n2,n3
|
|
|
181 |
stop
|
|
|
182 |
endif
|
|
|
183 |
|
|
|
184 |
c Activate or inactive the missing data check (quick and dirty)
|
|
|
185 |
if (misdat.ne.0.) then
|
|
|
186 |
mdv = misdat
|
|
|
187 |
else
|
|
|
188 |
mdv = 257.22725394015
|
|
|
189 |
endif
|
|
|
190 |
|
|
|
191 |
c Bring the indices into the grid space
|
|
|
192 |
ri = amax1(1.,amin1(float(n1),rid))
|
|
|
193 |
rj = amax1(1.,amin1(float(n2),rjd))
|
|
|
194 |
rk = amax1(1.,amin1(float(n3),rkd))
|
|
|
195 |
|
|
|
196 |
c Get the index of the west-south-bottom corner of the box
|
|
|
197 |
i = min0(int(ri),n1-1)
|
|
|
198 |
ip1 = i+1
|
|
|
199 |
j = min0(int(rj),n2-1)
|
|
|
200 |
jp1 = j+1
|
|
|
201 |
k = min0(int(rk),n3-1)
|
|
|
202 |
kp1 = k+1
|
|
|
203 |
|
|
|
204 |
c Special handling for 2d arrays
|
|
|
205 |
if (n3.eq.1) then
|
|
|
206 |
k=1
|
|
|
207 |
kp1=1
|
|
|
208 |
endif
|
|
|
209 |
|
|
|
210 |
c Get location relative to grid box
|
|
|
211 |
if ( i.ne.ip1 ) then
|
|
|
212 |
frac0i = ri-float(i)
|
|
|
213 |
frac1i = 1.-frac0i
|
|
|
214 |
else
|
|
|
215 |
frac0i = 0.
|
|
|
216 |
frac1i = 1.
|
|
|
217 |
endif
|
|
|
218 |
if ( j.ne.jp1 ) then
|
|
|
219 |
frac0j = rj-float(j)
|
|
|
220 |
frac1j = 1.-frac0j
|
|
|
221 |
else
|
|
|
222 |
frac0j = 0.
|
|
|
223 |
frac1j = 1.
|
|
|
224 |
endif
|
|
|
225 |
if ( k.ne.kp1 ) then
|
|
|
226 |
frac0k = rk-float(k)
|
|
|
227 |
frac1k = 1.-frac0k
|
|
|
228 |
else
|
|
|
229 |
frac0k = 0.
|
|
|
230 |
frac1k = 1.
|
|
|
231 |
endif
|
|
|
232 |
|
|
|
233 |
c On a grid point - take the grid point value
|
|
|
234 |
if ( ( abs(frac0i).lt.eps ).and.
|
|
|
235 |
> ( abs(frac0j).lt.eps ).and.
|
|
|
236 |
> ( abs(frac0k).lt.eps ) ) then
|
|
|
237 |
|
|
|
238 |
val = ar( i + n1*(j -1) + n1*n2*(k -1) )
|
|
|
239 |
goto 100
|
|
|
240 |
|
|
|
241 |
endif
|
|
|
242 |
|
|
|
243 |
c Init the fractions
|
|
|
244 |
frc000 = frac1i * frac1j * frac1k
|
|
|
245 |
frc001 = frac0i * frac1j * frac1k
|
|
|
246 |
frc010 = frac1i * frac0j * frac1k
|
|
|
247 |
frc011 = frac0i * frac0j * frac1k
|
|
|
248 |
frc100 = frac1i * frac1j * frac0k
|
|
|
249 |
frc101 = frac0i * frac1j * frac0k
|
|
|
250 |
frc110 = frac1i * frac0j * frac0k
|
|
|
251 |
frc111 = frac0i * frac0j * frac0k
|
|
|
252 |
|
|
|
253 |
c Init the values
|
|
|
254 |
val000 = ar( i + n1*(j -1) + n1*n2*(k -1) )
|
|
|
255 |
val001 = ar( ip1 + n1*(j -1) + n1*n2*(k -1) )
|
|
|
256 |
val010 = ar( i + n1*(jp1-1) + n1*n2*(k -1) )
|
|
|
257 |
val011 = ar( ip1 + n1*(jp1-1) + n1*n2*(k -1) )
|
|
|
258 |
val100 = ar( i + n1*(j -1) + n1*n2*(kp1-1) )
|
|
|
259 |
val101 = ar( ip1 + n1*(j -1) + n1*n2*(kp1-1) )
|
|
|
260 |
val110 = ar( i + n1*(jp1-1) + n1*n2*(kp1-1) )
|
|
|
261 |
val111 = ar( ip1 + n1*(jp1-1) + n1*n2*(kp1-1) )
|
|
|
262 |
|
|
|
263 |
c Handle missing data
|
|
|
264 |
if ( abs(val000-mdv).lt.eps ) frc000 = 0.
|
|
|
265 |
if ( abs(val001-mdv).lt.eps ) frc001 = 0.
|
|
|
266 |
if ( abs(val010-mdv).lt.eps ) frc010 = 0.
|
|
|
267 |
if ( abs(val011-mdv).lt.eps ) frc011 = 0.
|
|
|
268 |
if ( abs(val100-mdv).lt.eps ) frc100 = 0.
|
|
|
269 |
if ( abs(val101-mdv).lt.eps ) frc101 = 0.
|
|
|
270 |
if ( abs(val110-mdv).lt.eps ) frc110 = 0.
|
|
|
271 |
if ( abs(val111-mdv).lt.eps ) frc111 = 0.
|
|
|
272 |
|
|
|
273 |
c Build the final value
|
|
|
274 |
frc = frc000 + frc001 + frc010 + frc011 +
|
|
|
275 |
> frc100 + frc101 + frc110 + frc111
|
|
|
276 |
if ( frc.gt.0. ) then
|
|
|
277 |
val = 1./frc * ( frc000 * val000 + frc001 * val001 +
|
|
|
278 |
> frc010 * val010 + frc011 * val011 +
|
|
|
279 |
> frc100 * val100 + frc101 * val101 +
|
|
|
280 |
> frc110 * val110 + frc111 * val111 )
|
|
|
281 |
else
|
|
|
282 |
val = misdat
|
|
|
283 |
endif
|
|
|
284 |
|
|
|
285 |
c Return the value
|
|
|
286 |
100 continue
|
|
|
287 |
|
|
|
288 |
int_index3 = val
|
|
|
289 |
|
|
|
290 |
end
|
|
|
291 |
|
|
|
292 |
|
|
|
293 |
c -------------------------------------------------------------
|
|
|
294 |
c Time interpolation
|
|
|
295 |
c -------------------------------------------------------------
|
|
|
296 |
|
|
|
297 |
real function int_time (val0,val1,reltpos,misdat)
|
|
|
298 |
|
|
|
299 |
c Purpose:
|
|
|
300 |
c This subroutine interpolates linearly in time between two
|
|
|
301 |
c values.
|
|
|
302 |
c Arguments:
|
|
|
303 |
c val0 real input value at time 0
|
|
|
304 |
c val1 real input value at time 1
|
|
|
305 |
c reltpos real input relative time (between 0 and 1)
|
|
|
306 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
307 |
|
|
|
308 |
implicit none
|
|
|
309 |
|
|
|
310 |
c Declaration of parameters
|
|
|
311 |
real val0
|
|
|
312 |
real val1
|
|
|
313 |
real reltpos
|
|
|
314 |
real misdat
|
|
|
315 |
|
|
|
316 |
c Numerical epsilon
|
|
|
317 |
real eps
|
|
|
318 |
parameter (eps=1.e-8)
|
|
|
319 |
|
|
|
320 |
c Local variables
|
|
|
321 |
real val
|
|
|
322 |
real mdv
|
|
|
323 |
|
|
|
324 |
c Activate or inactive the missing data check (quick and dirty)
|
|
|
325 |
if (misdat.ne.0.) then
|
|
|
326 |
mdv = misdat
|
|
|
327 |
else
|
|
|
328 |
mdv = 257.22725394015
|
|
|
329 |
endif
|
|
|
330 |
|
|
|
331 |
c Do the linear interpolation
|
|
|
332 |
if ( abs(reltpos).lt.eps ) then
|
|
|
333 |
val = val0
|
|
|
334 |
elseif ( abs(reltpos-1.).lt.eps ) then
|
|
|
335 |
val = val1
|
|
|
336 |
elseif ( (abs(val0-mdv).gt.eps).and.
|
|
|
337 |
> (abs(val1-mdv).gt.eps) ) then
|
|
|
338 |
val = (1.-reltpos)*val0+reltpos*val1
|
|
|
339 |
else
|
|
|
340 |
val = mdv
|
|
|
341 |
endif
|
|
|
342 |
|
|
|
343 |
c Return value
|
|
|
344 |
int_time = val
|
|
|
345 |
|
|
|
346 |
end
|
|
|
347 |
|
|
|
348 |
|
|
|
349 |
c -------------------------------------------------------------
|
|
|
350 |
c Get the position of a physical point in grid space
|
|
|
351 |
c -------------------------------------------------------------
|
|
|
352 |
|
|
|
353 |
subroutine get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
|
|
|
354 |
> vert,surf,nx,ny,nz,lonw,lats,dlon,dlat)
|
|
|
355 |
|
|
|
356 |
c Purpose:
|
|
|
357 |
c This subroutine determines the indices (rid,rjd,rkd) in grid
|
|
|
358 |
c space for a point in physical space (xpo,ypo,ppo). The
|
|
|
359 |
c horizontal grid is specified by the south-west point (lats,lonw)
|
|
|
360 |
c and the grid spacing (dlat,dlon). The vertical grid is given
|
|
|
361 |
c by <vert(n1,n2,n3)>. The lower boundary (typicall surface
|
|
|
362 |
c pressure) is given by <surf(n1,n2)>.
|
|
|
363 |
c Arguments:
|
|
|
364 |
c rid,rjd,rkd real output grid location to be interpolated to
|
|
|
365 |
c xpo,ypo,ppo real input physical coordinates
|
|
|
366 |
c n1,n2,n3 int input grid dimensions in x-, y- and p-direction
|
|
|
367 |
c lats,lonw real input south and west boundary of grid space
|
|
|
368 |
c vert real input vertical coordinate grid
|
|
|
369 |
c surf real input lower boundary (surface pressure)
|
|
|
370 |
c mode int input direction of vertical axis
|
|
|
371 |
c 1: linear, 1 -> nz (th)
|
|
|
372 |
c 2: linear, nz -> 1 (pv)
|
|
|
373 |
c 3: binary (z)
|
|
|
374 |
c 4: binary (p)
|
|
|
375 |
|
|
|
376 |
implicit none
|
|
|
377 |
|
|
|
378 |
c Declartion of function parameters
|
|
|
379 |
integer nx,ny,nz
|
|
|
380 |
real vert(nx*ny*nz)
|
|
|
381 |
real surf(nx*ny)
|
|
|
382 |
real rid,rjd,rkd
|
|
|
383 |
real xpo,ypo,ppo
|
|
|
384 |
real dlat,dlon,lats,lonw
|
|
|
385 |
integer mode
|
|
|
386 |
|
|
|
387 |
c Numerical epsilon
|
|
|
388 |
real eps
|
|
|
389 |
parameter (eps=1.e-8)
|
|
|
390 |
|
|
|
391 |
c Local variables
|
|
|
392 |
integer i,j,k
|
|
|
393 |
real ppo0,ppo1,ppom,psur
|
|
|
394 |
integer i0,im,i1
|
|
|
395 |
|
|
|
396 |
c Externals
|
|
|
397 |
real int_index3
|
|
|
398 |
external int_index3
|
|
|
399 |
|
|
|
400 |
c Get the horizontal grid indices
|
|
|
401 |
rid=(xpo-lonw)/dlon+1.
|
|
|
402 |
rjd=(ypo-lats)/dlat+1.
|
|
|
403 |
|
|
|
404 |
c Two-dimensional interpolation on horizontal plane: return
|
|
|
405 |
if ( nz.eq.1 ) then
|
|
|
406 |
rkd = 1.
|
|
|
407 |
goto 100
|
|
|
408 |
endif
|
|
|
409 |
|
|
|
410 |
c Lowest-level interpolation: return
|
|
|
411 |
if ( abs(ppo).lt.eps ) then
|
|
|
412 |
rkd = 1.
|
|
|
413 |
goto 100
|
|
|
414 |
endif
|
|
|
415 |
|
|
|
416 |
c Get the pressure at the lowest level and at the surface
|
|
|
417 |
ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)
|
|
|
418 |
psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)
|
|
|
419 |
|
|
|
420 |
c The point is between the surface and the lowest level: return
|
|
|
421 |
if ( (ppo.ge.ppo0).and.(ppo.le.psur).or.
|
|
|
422 |
> (ppo.le.ppo0).and.(ppo.ge.psur) )
|
|
|
423 |
>then
|
|
|
424 |
psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)
|
|
|
425 |
rkd = (psur-ppo)/(psur-ppo0)
|
|
|
426 |
goto 100
|
|
|
427 |
endif
|
|
|
428 |
|
|
|
429 |
c Full-level search (TH): linear ascending scanning through all levels
|
|
|
430 |
if ( mode.eq.1 ) then
|
|
|
431 |
|
|
|
432 |
ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)
|
|
|
433 |
rkd=0
|
|
|
434 |
do i=1,nz-1
|
|
|
435 |
ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(i+1),0.)
|
|
|
436 |
if ( (ppo0.lt.ppo).and.(ppo1.ge.ppo) ) then
|
|
|
437 |
rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)
|
|
|
438 |
goto 100
|
|
|
439 |
endif
|
|
|
440 |
ppo0 = ppo1
|
|
|
441 |
enddo
|
|
|
442 |
|
|
|
443 |
c Full-level search (PV): linear descending scanning through all levels
|
|
|
444 |
elseif ( mode.eq.2 ) then
|
|
|
445 |
|
|
|
446 |
ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
|
|
|
447 |
rkd=0
|
|
|
448 |
do i=nz-1,1,-1
|
|
|
449 |
ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(i),0.)
|
|
|
450 |
if ( (ppo1.gt.ppo).and.(ppo0.le.ppo) ) then
|
|
|
451 |
rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)
|
|
|
452 |
goto 100
|
|
|
453 |
endif
|
|
|
454 |
ppo1 = ppo0
|
|
|
455 |
enddo
|
|
|
456 |
|
|
|
457 |
c Full-level search (z): binary search
|
|
|
458 |
elseif ( mode.eq.3 ) then
|
|
|
459 |
|
|
|
460 |
rkd = 0
|
|
|
461 |
i0 = 1
|
|
|
462 |
i1 = nz
|
|
|
463 |
ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real( 1),0.)
|
|
|
464 |
ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
|
|
|
465 |
|
|
|
466 |
do while ( i1.gt.(i0+1) )
|
|
|
467 |
im = (i0+i1)/2
|
|
|
468 |
ppom = int_index3(vert,nx,ny,nz,rid,rjd,real(im),0.)
|
|
|
469 |
if (ppom.gt.ppo) then
|
|
|
470 |
i1 = im
|
|
|
471 |
ppo1 = ppom
|
|
|
472 |
else
|
|
|
473 |
i0 = im
|
|
|
474 |
ppo0 = ppom
|
|
|
475 |
endif
|
|
|
476 |
|
|
|
477 |
enddo
|
|
|
478 |
|
|
|
479 |
rkd=real(i0)+(ppo0-ppo)/(ppo0-ppo1)
|
|
|
480 |
|
|
|
481 |
c Full-level search (P): binary search
|
|
|
482 |
elseif ( mode.eq.4 ) then
|
|
|
483 |
|
|
|
484 |
rkd = 0
|
|
|
485 |
i0 = 1
|
|
|
486 |
i1 = nz
|
|
|
487 |
ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real( 1),0.)
|
|
|
488 |
ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
|
|
|
489 |
|
|
|
490 |
do while ( i1.gt.(i0+1) )
|
|
|
491 |
im = (i0+i1)/2
|
|
|
492 |
ppom = int_index3(vert,nx,ny,nz,rid,rjd,real(im),0.)
|
|
|
493 |
if (ppom.lt.ppo) then
|
|
|
494 |
i1 = im
|
|
|
495 |
ppo1 = ppom
|
|
|
496 |
else
|
|
|
497 |
i0 = im
|
|
|
498 |
ppo0 = ppom
|
|
|
499 |
endif
|
|
|
500 |
enddo
|
|
|
501 |
|
|
|
502 |
rkd=real(i0)+(ppo0-ppo)/(ppo0-ppo1)
|
|
|
503 |
|
|
|
504 |
endif
|
|
|
505 |
|
|
|
506 |
c Exit point for subroutine
|
|
|
507 |
100 continue
|
|
|
508 |
|
|
|
509 |
end
|
|
|
510 |
|