V3FIT
r8xlookup.f
1  subroutine r8xlookup(ivec,xvec,nx,xpkg,imode,iv,dxn,hv,hiv,iwarn)
2 c
3 c vector lookup routine
4 c
5 c given a set of x points xvec(...) and an x grid (xpkg, nx x pts)
6 c return the vector of indices iv(...) and displacements dxv(...)
7 c within the indexed zones, corresponding to each x point.
8 c
9 c if any of the x points in the vector are out of range the warning
10 c flag iwarn is set.
11 c
12 c MOD DMC Feb 2010: changes related to supporting nx=2 and nx=3 small grids.
13 c Changes are consistent with Feb 2010 changes in genxpkg:
14 c meanings of xpkg(1,4) and xpkg(3,4) interchanged.
15 c
16 c if nx.eq.2: xpkg(3,4) and xpkg(4,4) never referenced;
17 c if nx.eq.3: xpkg(4,4) never referenced.
18 c
19 c----------------------------
20 c input:
21 !============
22 ! idecl: explicitize implicit INTEGER declarations:
23  IMPLICIT NONE
24  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
25  INTEGER inum,istat,ilin,ialg,iper,imsg,init_guess,iprev,i
26  INTEGER init,iprob,isrch,i_sign,inc
27 !============
28 ! idecl: explicitize implicit REAL declarations:
29  real*8 stat,ztola,period,hav,havi,hloci,hloc,zdelta,xfac
30  real*8 zindx0,zdindx,zindex
31 !============
32  integer ivec ! size of vector
33  real*8 xvec(ivec) ! x points to lookup on xpkg grid
34 c
35  integer nx ! size of grid
36  real*8 xpkg(nx,4) ! grid data
37 c
38  integer imode ! output control flag
39 c
40 c imode=1: return indices iv(...) and un-normalized displacements dxn(...)
41 c ignore hv and hiv
42 c dxn(j)=xvec(j)-xpkg(iv(j),1)
43 c
44 c imode=2 return indices iv(...) and *normalized* displacements dxn(...)
45 c and hv(...) and hiv(...)
46 c dxn(j)=(xvec(j)-xpkg(iv(j),1))*hiv(j)
47 c
48 c output:
49  integer iv(ivec) ! index into grid for each xvec(j)
50  ! note: old values of iv(...) may be used as start point for grid
51  ! searches, depending on xpkg controls
52 
53  real*8 dxn(ivec) ! displacement w/in zone (see imode)
54  real*8 hv(ivec) ! zone width (if imode=2)
55  real*8 hiv(ivec) ! inverse zone width (if imode=2)
56 c
57  integer iwarn ! =0: OK; =n: n points out of range
58 c
59 c----------------------------
60 c
61 c xpkg is a "structure" constructed by subroutine genxpkg, which
62 c contains the x grid, xpkg(1:nx,1), spacing and error handling
63 c information -- see genxpkg.for (& r8genxpkg.for)
64 c
65  real*8, dimension(:), allocatable :: xuse
66  logical, dimension(:), allocatable :: iok
67  integer, dimension(:), allocatable :: imina,imaxa
68 c
69 cdbg data idbg/0/
70 c
71 c----------------------------
72  if(nx.lt.2) then
73  iwarn=1
74  write(6,*) .lt.' ?? xlookup: nx2, nx=',nx
75  go to 1100
76  endif
77 
78  inum=ivec
79  allocate(xuse(inum),stat=istat)
80  if(istat.ne.0) then
81  iwarn=1
82  write(6,*) ' ?? xlookup "xuse" vector allocation failure!'
83  go to 1000
84  endif
85 c
86  if(nx.eq.2) then
87  ilin=1
88  endif
89 c
90  if(nx.gt.2) then
91  if(xpkg(3,4).eq.0.0_r8) then
92  ilin=1 ! evenly spaced grid
93  else
94  ilin=0
95  if(xpkg(3,4).gt.2.5_r8) then
96  ialg=3
97  else if(xpkg(3,4).gt.1.5_r8) then
98  ialg=2
99  else
100  ialg=1
101  endif
102  endif
103  endif
104 c
105  if(xpkg(2,4).ne.0.0_r8) then
106  iper=1 ! periodic grid
107  else
108  iper=0
109  endif
110 c
111  ztola=abs(xpkg(1,4)) ! tolerance for range checking
112 c
113  if(xpkg(1,4).ge.0.0_r8) then
114  imsg=1 ! write message on range error
115  else
116  imsg=0
117  endif
118 c
119  init_guess=0
120  if(nx.gt.3) then
121  if(xpkg(4,4).gt.0.0_r8) then
122  init_guess=1
123  iprev=min((nx-1),max(1,iv(ivec)))
124  else
125  init_guess=0
126  endif
127  endif
128 c
129  iwarn=0
130 c
131 c---------------------
132 c range check
133 c
134  if(iper.eq.0) then
135 c
136 c check min/max with tolerance
137 c
138  do i=1,ivec
139  if(xvec(i).lt.xpkg(1,1)) then
140  xuse(i)=xpkg(1,1)
141  if((xpkg(1,1)-xvec(i)).gt.ztola) iwarn=iwarn+1
142  else if(xvec(i).gt.xpkg(nx,1)) then
143  xuse(i)=xpkg(nx,1)
144  if((xvec(i)-xpkg(nx,1)).gt.ztola) iwarn=iwarn+1
145  else
146  xuse(i)=xvec(i)
147  endif
148  enddo
149 c
150  else
151 c
152 c normalize to interval
153 c
154  period=xpkg(nx,1)-xpkg(1,1)
155  do i=1,ivec
156  if((xvec(i).lt.xpkg(1,1)).or.(xvec(i).gt.xpkg(nx,1))) then
157  xuse(i)=mod(xvec(i)-xpkg(1,1),period)
158  if(xuse(i).lt.0.0_r8) xuse(i)=xuse(i)+period
159  xuse(i)=xuse(i)+xpkg(1,1)
160  xuse(i)=max(xpkg(1,1),min(xpkg(nx,1),xuse(i)))
161  else
162  xuse(i)=xvec(i)
163  endif
164  enddo
165 c
166  endif
167 c
168  if((imsg.eq.1).and.(iwarn.gt.0)) then
169  write(6,*) ' %xlookup: ',iwarn,' points not in range: ',
170  > xpkg(1,1),' to ',xpkg(nx,1)
171  endif
172 c
173 c---------------------
174 c actual lookup -- initially, assume even spacing
175 c
176 c initial index guess: 1 + <1/h>*(x-x0) ... then refine
177 c
178  hav=xpkg(nx,2)
179  havi=xpkg(nx,3)
180  if(ilin.eq.1) then
181 c
182 c faster lookup OK: even spacing
183 c
184  if(init_guess.eq.0) then
185 c
186 c even spacing lookup, no initial guess from previous iteration
187 c
188  do i=1,ivec
189  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
190  iv(i)=max(1,min((nx-1),iv(i)))
191  if(imode.eq.1) then
192  dxn(i)=(xuse(i)-xpkg(iv(i),1))
193  else
194  dxn(i)=(xuse(i)-xpkg(iv(i),1))*havi
195  hiv(i)=havi
196  hv(i)=hav
197  endif
198  enddo
199 c
200  else
201 c
202 c even spacing lookup, do use initial guess from previous iteration
203 c
204  do i=1,ivec
205  if((xpkg(iprev,1).le.xuse(i)).and.
206  > (xuse(i).le.xpkg(iprev+1,1))) then
207  iv(i)=iprev
208  else
209  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
210  iv(i)=max(1,min((nx-1),iv(i)))
211  iprev=iv(i)
212  endif
213  if(imode.eq.1) then
214  dxn(i)=(xuse(i)-xpkg(iv(i),1))
215  else
216  dxn(i)=(xuse(i)-xpkg(iv(i),1))*havi
217  hiv(i)=havi
218  hv(i)=hav
219  endif
220  enddo
221 c
222  endif
223 c
224  go to 1000
225 c
226  endif
227 c
228  4 continue
229 c
230  init=-1
231  allocate(iok(inum),stat=istat)
232  if(istat.ne.0) then
233  iwarn=1
234  write(6,*) ' ?? xlookup "iok" vector allocation failure!'
235  go to 1000
236  endif
237 c
238  if(ialg.lt.3) then
239  allocate(imina(inum),stat=istat)
240  if(istat.ne.0) then
241  iwarn=1
242  write(6,*) ' ?? xlookup "imina" vector allocation failure!'
243  go to 1000
244  endif
245 c
246  allocate(imaxa(inum),stat=istat)
247  if(istat.ne.0) then
248  iwarn=1
249  write(6,*) ' ?? xlookup "imaxa" vector allocation failure!'
250  go to 1000
251  endif
252  endif
253 c
254  5 continue ! re-entry -- hit problem cases...
255 c
256  iprob=0 ! count "problem" cases...
257  init=init+1
258 c
259  if(init_guess.eq.0) then
260  go to (100,200,300) ialg
261  else
262  go to (150,250,350) ialg
263  endif
264 c-----------------------------------------------------------
265 c Newton like algorithm: use local spacing to estimate
266 c step to next zone guess; don't use prior guess
267 c
268  100 continue
269  do i=1,ivec
270 cdbg jdbg=0
271 c
272 c first iteration
273 c
274  if(init.eq.0) then
275  iok(i)=.false.
276  imina(i)=1
277  imaxa(i)=nx-1
278 c
279  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
280 cdbg jdbg=jdbg+1
281 cdbg if(jdbg.le.idbg) iv(i)=1
282  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
283 c
284  if(xuse(i).lt.xpkg(iv(i),1)) then
285  isrch=0
286  i_sign=-1
287  imaxa(i)=max(1,(iv(i)-1))
288  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
289  isrch=1
290  i_sign=+1
291  imina(i)=min((nx-1),(iv(i)+1))
292  else
293  iok(i)=.true.
294  endif
295 c
296  endif
297 c
298 c second iteration
299 c
300  if(.not.iok(i)) then
301  hloci=xpkg(iv(i),3)
302  hloc=xpkg(iv(i),2)
303  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
304  if(i_sign*zdelta.le.hloc) then
305  inc=i_sign
306  else
307  inc=zdelta*hloci
308  inc=inc+i_sign
309  endif
310 c
311  iv(i)=iv(i)+inc
312 cdbg jdbg=jdbg+1
313 cdbg if(jdbg.le.idbg) iv(i)=1
314  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
315 c
316  if(xuse(i).lt.xpkg(iv(i),1)) then
317  isrch=0
318  i_sign=-1
319  imaxa(i)=max(1,(iv(i)-1))
320  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
321  isrch=1
322  i_sign=+1
323  imina(i)=min((nx-1),(iv(i)+1))
324  else
325  iok(i)=.true.
326  endif
327 c
328 c third iteration
329 c
330  if(.not.iok(i)) then
331  hloci=xpkg(iv(i),3)
332  hloc=xpkg(iv(i),2)
333  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
334  if(i_sign*zdelta.le.hloc) then
335  inc=i_sign
336  else
337  inc=zdelta*hloci
338  inc=inc+i_sign
339  endif
340 c
341  iv(i)=iv(i)+inc
342 cdbg jdbg=jdbg+1
343 cdbg if(jdbg.le.idbg) iv(i)=1
344  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
345 c
346  if(xuse(i).lt.xpkg(iv(i),1)) then
347  isrch=0
348  i_sign=-1
349  imaxa(i)=max(1,(iv(i)-1))
350  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
351  isrch=1
352  i_sign=+1
353  imina(i)=min((nx-1),(iv(i)+1))
354  else
355  iok(i)=.true.
356  endif
357 c
358 c fourth iteration
359 c
360  if(.not.iok(i)) then
361  hloci=xpkg(iv(i),3)
362  hloc=xpkg(iv(i),2)
363  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
364  if(i_sign*zdelta.le.hloc) then
365  inc=i_sign
366  else
367  inc=zdelta*hloci
368  inc=inc+i_sign
369  endif
370 c
371  iv(i)=iv(i)+inc
372 cdbg jdbg=jdbg+1
373 cdbg if(jdbg.le.idbg) iv(i)=1
374  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
375 c
376  if(xuse(i).lt.xpkg(iv(i),1)) then
377  isrch=0
378  i_sign=-1
379  imaxa(i)=max(1,(iv(i)-1))
380  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
381  isrch=1
382  i_sign=+1
383  imina(i)=min((nx-1),(iv(i)+1))
384  else
385  iok(i)=.true.
386  endif
387 c
388 c fifth iteration
389 c
390  if(.not.iok(i)) then
391  hloci=xpkg(iv(i),3)
392  hloc=xpkg(iv(i),2)
393  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
394  if(i_sign*zdelta.le.hloc) then
395  inc=i_sign
396  else
397  inc=zdelta*hloci
398  inc=inc+i_sign
399  endif
400 c
401  iv(i)=iv(i)+inc
402 cdbg jdbg=jdbg+1
403 cdbg if(jdbg.le.idbg) iv(i)=1
404  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
405 c
406  if(xuse(i).lt.xpkg(iv(i),1)) then
407  isrch=0
408  i_sign=-1
409  imaxa(i)=max(1,(iv(i)-1))
410  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
411  isrch=1
412  i_sign=+1
413  imina(i)=min((nx-1),(iv(i)+1))
414  else
415  iok(i)=.true.
416  endif
417 c
418  if(.not.iok(i)) iprob=iprob+1
419 c
420 c end chain of iteration if-then-else blocks
421 c
422  endif
423  endif
424  endif
425  endif
426 c
427 c end of loop
428 c
429  enddo
430 c
431  go to 500
432 c-----------------------------------------------------------
433 c Newton like algorithm: use local spacing to estimate
434 c step to next zone guess; DO use prior guess
435 c
436  150 continue
437 c
438  do i=1,ivec
439 c
440 c first iteration
441 c
442  if(init.eq.0) then
443  if((xpkg(iprev,1).le.xuse(i)).and.
444  > (xuse(i).le.xpkg(iprev+1,1))) then
445  iok(i)=.true.
446  iv(i)=iprev
447  else
448  iok(i)=.false.
449  imina(i)=1
450  imaxa(i)=nx-1
451 c
452  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
453  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
454 c
455  if(xuse(i).lt.xpkg(iv(i),1)) then
456  isrch=0
457  i_sign=-1
458  imaxa(i)=max(1,(iv(i)-1))
459  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
460  isrch=1
461  i_sign=+1
462  imina(i)=min((nx-1),(iv(i)+1))
463  else
464  iok(i)=.true.
465  endif
466  endif
467 c
468  endif
469 c
470 c second iteration
471 c
472  if(.not.iok(i)) then
473  hloci=xpkg(iv(i),3)
474  hloc=xpkg(iv(i),2)
475  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
476  if(i_sign*zdelta.le.hloc) then
477  inc=i_sign
478  else
479  inc=zdelta*hloci
480  inc=inc+i_sign
481  endif
482 c
483  iv(i)=iv(i)+inc
484  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
485 c
486  if(xuse(i).lt.xpkg(iv(i),1)) then
487  isrch=0
488  i_sign=-1
489  imaxa(i)=max(1,(iv(i)-1))
490  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
491  isrch=1
492  i_sign=+1
493  imina(i)=min((nx-1),(iv(i)+1))
494  else
495  iok(i)=.true.
496  endif
497 c
498 c third iteration
499 c
500  if(.not.iok(i)) then
501  hloci=xpkg(iv(i),3)
502  hloc=xpkg(iv(i),2)
503  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
504  if(i_sign*zdelta.le.hloc) then
505  inc=i_sign
506  else
507  inc=zdelta*hloci
508  inc=inc+i_sign
509  endif
510 c
511  iv(i)=iv(i)+inc
512  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
513 c
514  if(xuse(i).lt.xpkg(iv(i),1)) then
515  isrch=0
516  i_sign=-1
517  imaxa(i)=max(1,(iv(i)-1))
518  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
519  isrch=1
520  i_sign=+1
521  imina(i)=min((nx-1),(iv(i)+1))
522  else
523  iok(i)=.true.
524  endif
525 c
526 c fourth iteration
527 c
528  if(.not.iok(i)) then
529  hloci=xpkg(iv(i),3)
530  hloc=xpkg(iv(i),2)
531  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
532  if(i_sign*zdelta.le.hloc) then
533  inc=i_sign
534  else
535  inc=zdelta*hloci
536  inc=inc+i_sign
537  endif
538 c
539  iv(i)=iv(i)+inc
540  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
541 c
542  if(xuse(i).lt.xpkg(iv(i),1)) then
543  isrch=0
544  i_sign=-1
545  imaxa(i)=max(1,(iv(i)-1))
546  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
547  isrch=1
548  i_sign=+1
549  imina(i)=min((nx-1),(iv(i)+1))
550  else
551  iok(i)=.true.
552  endif
553 c
554 c fifth iteration
555 c
556  if(.not.iok(i)) then
557  hloci=xpkg(iv(i),3)
558  hloc=xpkg(iv(i),2)
559  zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
560  if(i_sign*zdelta.le.hloc) then
561  inc=i_sign
562  else
563  inc=zdelta*hloci
564  inc=inc+i_sign
565  endif
566 c
567  iv(i)=iv(i)+inc
568  iv(i)=max(imina(i),min(imaxa(i),iv(i)))
569 c
570  if(xuse(i).lt.xpkg(iv(i),1)) then
571  isrch=0
572  i_sign=-1
573  imaxa(i)=max(1,(iv(i)-1))
574  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
575  isrch=1
576  i_sign=+1
577  imina(i)=min((nx-1),(iv(i)+1))
578  else
579  iok(i)=.true.
580  endif
581 c
582  if(.not.iok(i)) iprob=iprob+1
583 c
584 c end chain of iteration if-then-else blocks
585 c
586  endif
587  endif
588  endif
589  endif
590 c
591 c end of loop
592 c
593  iprev=iv(i)
594  enddo
595 c
596  go to 500
597 c-----------------------------------------------------------
598 c Binary search algorithm
599 c
600  200 continue
601  do i=1,ivec
602 c
603 c first iteration
604 c
605  if(init.eq.0) then
606  iok(i)=.false.
607  imina(i)=1
608  imaxa(i)=nx-1
609 c
610  iv(i)=nx/2
611 c
612  if(xuse(i).lt.xpkg(iv(i),1)) then
613  imaxa(i)=max(1,(iv(i)-1))
614  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
615  imina(i)=min((nx-1),(iv(i)+1))
616  else
617  iok(i)=.true.
618  endif
619 c
620  endif
621 c
622 c second iteration
623 c
624  if(.not.iok(i)) then
625 c
626  iv(i)=(imina(i)+imaxa(i))/2
627 c
628  if(xuse(i).lt.xpkg(iv(i),1)) then
629  imaxa(i)=max(1,(iv(i)-1))
630  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
631  imina(i)=min((nx-1),(iv(i)+1))
632  else
633  iok(i)=.true.
634  endif
635 c
636 c third iteration
637 c
638  if(.not.iok(i)) then
639 c
640  iv(i)=(imina(i)+imaxa(i))/2
641 c
642  if(xuse(i).lt.xpkg(iv(i),1)) then
643  imaxa(i)=max(1,(iv(i)-1))
644  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
645  imina(i)=min((nx-1),(iv(i)+1))
646  else
647  iok(i)=.true.
648  endif
649 c
650 c fourth iteration
651 c
652  if(.not.iok(i)) then
653 c
654  iv(i)=(imina(i)+imaxa(i))/2
655 c
656  if(xuse(i).lt.xpkg(iv(i),1)) then
657  imaxa(i)=max(1,(iv(i)-1))
658  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
659  imina(i)=min((nx-1),(iv(i)+1))
660  else
661  iok(i)=.true.
662  endif
663 c
664 c fifth iteration
665 c
666  if(.not.iok(i)) then
667 c
668  iv(i)=(imina(i)+imaxa(i))/2
669 c
670  if(xuse(i).lt.xpkg(iv(i),1)) then
671  imaxa(i)=max(1,(iv(i)-1))
672  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
673  imina(i)=min((nx-1),(iv(i)+1))
674  else
675  iok(i)=.true.
676  endif
677 c
678  if(.not.iok(i)) iprob=iprob+1
679 c
680 c end chain of iteration if-then-else blocks
681 c
682  endif
683  endif
684  endif
685  endif
686 c
687 c end of loop
688 c
689  enddo
690 c
691  go to 500
692 c-----------------------------------------------------------
693 c Binary search algorithm
694 c
695  250 continue
696 c
697  do i=1,ivec
698 c
699 c first iteration
700 c
701  if(init.eq.0) then
702  if((xpkg(iprev,1).le.xuse(i)).and.
703  > (xuse(i).le.xpkg(iprev+1,1))) then
704  iok(i)=.true.
705  iv(i)=iprev
706  else
707  iok(i)=.false.
708  imina(i)=1
709  imaxa(i)=nx-1
710 c
711  iv(i)=nx/2
712 c
713  if(xuse(i).lt.xpkg(iv(i),1)) then
714  imaxa(i)=max(1,(iv(i)-1))
715  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
716  imina(i)=min((nx-1),(iv(i)+1))
717  else
718  iok(i)=.true.
719  endif
720  endif
721 c
722  endif
723 c
724 c second iteration
725 c
726  if(.not.iok(i)) then
727 c
728  iv(i)=(imina(i)+imaxa(i))/2
729 c
730  if(xuse(i).lt.xpkg(iv(i),1)) then
731  imaxa(i)=max(1,(iv(i)-1))
732  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
733  imina(i)=min((nx-1),(iv(i)+1))
734  else
735  iok(i)=.true.
736  endif
737 c
738 c third iteration
739 c
740  if(.not.iok(i)) then
741 c
742  iv(i)=(imina(i)+imaxa(i))/2
743 c
744  if(xuse(i).lt.xpkg(iv(i),1)) then
745  imaxa(i)=max(1,(iv(i)-1))
746  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
747  imina(i)=min((nx-1),(iv(i)+1))
748  else
749  iok(i)=.true.
750  endif
751 c
752 c fourth iteration
753 c
754  if(.not.iok(i)) then
755 c
756  iv(i)=(imina(i)+imaxa(i))/2
757 c
758  if(xuse(i).lt.xpkg(iv(i),1)) then
759  imaxa(i)=max(1,(iv(i)-1))
760  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
761  imina(i)=min((nx-1),(iv(i)+1))
762  else
763  iok(i)=.true.
764  endif
765 c
766 c fifth iteration
767 c
768  if(.not.iok(i)) then
769 c
770  iv(i)=(imina(i)+imaxa(i))/2
771 c
772  if(xuse(i).lt.xpkg(iv(i),1)) then
773  imaxa(i)=max(1,(iv(i)-1))
774  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
775  imina(i)=min((nx-1),(iv(i)+1))
776  else
777  iok(i)=.true.
778  endif
779 c
780  if(.not.iok(i)) iprob=iprob+1
781 c
782 c end chain of iteration if-then-else blocks
783 c
784  endif
785  endif
786  endif
787  endif
788 c
789 c end of loop
790 c
791  iprev=iv(i)
792  enddo
793 c
794  go to 500
795 c-----------------------------------------------------------
796 c algorithm: piecewise linear indexing function lookup & correction
797 c
798  300 continue
799  do i=1,ivec
800 c
801 c first iteration
802 c
803  if(init.eq.0) then
804  iok(i)=.false.
805 c
806 c piecewise linear indexing function on even spaced grid
807 c (not same grid as x axis itself)
808 c
809  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
810  iv(i)=max(1,min(nx-1,iv(i)))
811  xfac=(xuse(i)-(xpkg(1,1)+(iv(i)-1)*hav))*havi
812  zindx0=xpkg(iv(i),2)
813  if(iv(i).lt.nx-1) then
814  zdindx=xpkg(iv(i)+1,2)-zindx0
815  else
816  zdindx=nx-zindx0
817  endif
818  zindex=zindx0+xfac*zdindx
819 c
820  iv(i)=zindex
821  iv(i)=max(1,min(nx-1,iv(i)))
822 c
823  if(xuse(i).lt.xpkg(iv(i),1)) then
824  i_sign=-1
825  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
826  i_sign=+1
827  else
828  iok(i)=.true.
829  endif
830 c
831  endif
832 c
833 c second iteration
834 c
835  if(.not.iok(i)) then
836  iv(i)=iv(i)+i_sign
837  if(xuse(i).lt.xpkg(iv(i),1)) then
838  i_sign=-1
839  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
840  i_sign=+1
841  else
842  iok(i)=.true.
843  endif
844 c
845 c third iteration
846 c
847  if(.not.iok(i)) then
848  iv(i)=iv(i)+i_sign
849  if(xuse(i).lt.xpkg(iv(i),1)) then
850  i_sign=-1
851  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
852  i_sign=+1
853  else
854  iok(i)=.true.
855  endif
856 c
857 c fourth iteration
858 c
859  if(.not.iok(i)) then
860  iv(i)=iv(i)+i_sign
861  if(xuse(i).lt.xpkg(iv(i),1)) then
862  i_sign=-1
863  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
864  i_sign=+1
865  else
866  iok(i)=.true.
867  endif
868 c
869 c fifth iteration
870 c
871  if(.not.iok(i)) then
872  iv(i)=iv(i)+i_sign
873  if(xuse(i).lt.xpkg(iv(i),1)) then
874  i_sign=-1
875  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
876  i_sign=+1
877  else
878  iok(i)=.true.
879  endif
880 c
881  if(.not.iok(i)) iprob=iprob+1
882 c
883 c end chain of iteration if-then-else blocks
884 c
885  endif
886  endif
887  endif
888  endif
889 c
890 c end of loop
891 c
892  enddo
893 c
894  go to 500
895 c
896 c-----------------------------------------------------------
897 c algorithm: piecewise linear indexing function lookup & correction
898 c
899  350 continue
900  do i=1,ivec
901 c
902 c first iteration
903 c
904  if(init.eq.0) then
905  if((xpkg(iprev,1).le.xuse(i)).and.
906  > (xuse(i).le.xpkg(iprev+1,1))) then
907  iok(i)=.true.
908  iv(i)=iprev
909  else
910  iok(i)=.false.
911 c
912 c piecewise linear indexing function on even spaced grid
913 c (not same grid as x axis itself)
914 c
915  iv(i)=1+havi*(xuse(i)-xpkg(1,1))
916  iv(i)=max(1,min(nx-1,iv(i)))
917  xfac=(xuse(i)-(xpkg(1,1)+(iv(i)-1)*hav))*havi
918  zindx0=xpkg(iv(i),2)
919  if(iv(i).lt.nx-1) then
920  zdindx=xpkg(iv(i)+1,2)-zindx0
921  else
922  zdindx=nx-zindx0
923  endif
924  zindex=zindx0+xfac*zdindx
925 c
926  iv(i)=zindex
927  iv(i)=max(1,min(nx-1,iv(i)))
928 c
929  if(xuse(i).lt.xpkg(iv(i),1)) then
930  i_sign=-1
931  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
932  i_sign=+1
933  else
934  iok(i)=.true.
935  endif
936  endif
937 c
938  endif
939 c
940 c second iteration
941 c
942  if(.not.iok(i)) then
943  iv(i)=iv(i)+i_sign
944  if(xuse(i).lt.xpkg(iv(i),1)) then
945  i_sign=-1
946  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
947  i_sign=+1
948  else
949  iok(i)=.true.
950  endif
951 c
952 c third iteration
953 c
954  if(.not.iok(i)) then
955  iv(i)=iv(i)+i_sign
956  if(xuse(i).lt.xpkg(iv(i),1)) then
957  i_sign=-1
958  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
959  i_sign=+1
960  else
961  iok(i)=.true.
962  endif
963 c
964 c fourth iteration
965 c
966  if(.not.iok(i)) then
967  iv(i)=iv(i)+i_sign
968  if(xuse(i).lt.xpkg(iv(i),1)) then
969  i_sign=-1
970  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
971  i_sign=+1
972  else
973  iok(i)=.true.
974  endif
975 c
976 c fifth iteration
977 c
978  if(.not.iok(i)) then
979  iv(i)=iv(i)+i_sign
980  if(xuse(i).lt.xpkg(iv(i),1)) then
981  i_sign=-1
982  else if(xuse(i).gt.xpkg(iv(i)+1,1)) then
983  i_sign=+1
984  else
985  iok(i)=.true.
986  endif
987 c
988  if(.not.iok(i)) iprob=iprob+1
989 c
990 c end chain of iteration if-then-else blocks
991 c
992  endif
993  endif
994  endif
995  endif
996 c
997 c end of loop
998 c
999  iprev=iv(i)
1000  enddo
1001 c
1002  go to 500
1003 c
1004 c--------------------------------------------------------------------
1005 c any "problems" left? if so, re-enter the loop...
1006 c
1007  500 continue
1008  if(iprob.gt.0) go to 5
1009 c
1010 c OK -- all zones found; complete output stats
1011 c
1012  if(imode.eq.1) then
1013  dxn=(xuse-xpkg(iv,1)) ! un-normalized
1014  else if(ialg.ne.3) then
1015  dxn=(xuse-xpkg(iv,1))*xpkg(iv,3) ! normalized to 1/h in each zone
1016  hv=xpkg(iv,2)
1017  hiv=xpkg(iv,3)
1018  else
1019  dxn=(xuse-xpkg(iv,1))*xpkg(iv,3) ! normalized to 1/h in each zone
1020  hv=xpkg(iv+1,1)-xpkg(iv,1)
1021  hiv=xpkg(iv,3)
1022  endif
1023 c
1024  deallocate(iok)
1025  if(ialg.lt.3) then
1026  deallocate(imina)
1027  deallocate(imaxa)
1028  endif
1029 c
1030 c all done -- generalized lookup
1031 c
1032  1000 continue
1033  deallocate(xuse)
1034 c
1035  1100 continue
1036  return
1037  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17