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