1 subroutine r8xlookup(ivec,xvec,nx,xpkg,imode,iv,dxn,hv,hiv,iwarn)
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
29 real*8 stat,ztola,period,hav,havi,hloci,hloc,zdelta,xfac
30 real*8 zindx0,zdindx,zindex
65 real*8,
dimension(:),
allocatable :: xuse
66 logical,
dimension(:),
allocatable :: iok
67 integer,
dimension(:),
allocatable :: imina,imaxa
74 write(6,*) .lt.
' ?? xlookup: nx2, nx=',nx
79 allocate(xuse(inum),stat=istat)
82 write(6,*)
' ?? xlookup "xuse" vector allocation failure!'
91 if(xpkg(3,4).eq.0.0_r8)
then
95 if(xpkg(3,4).gt.2.5_r8)
then
97 else if(xpkg(3,4).gt.1.5_r8)
then
105 if(xpkg(2,4).ne.0.0_r8)
then
113 if(xpkg(1,4).ge.0.0_r8)
then
121 if(xpkg(4,4).gt.0.0_r8)
then
123 iprev=min((nx-1),max(1,iv(ivec)))
139 if(xvec(i).lt.xpkg(1,1))
then
141 if((xpkg(1,1)-xvec(i)).gt.ztola) iwarn=iwarn+1
142 else if(xvec(i).gt.xpkg(nx,1))
then
144 if((xvec(i)-xpkg(nx,1)).gt.ztola) iwarn=iwarn+1
154 period=xpkg(nx,1)-xpkg(1,1)
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)))
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)
184 if(init_guess.eq.0)
then
189 iv(i)=1+havi*(xuse(i)-xpkg(1,1))
190 iv(i)=max(1,min((nx-1),iv(i)))
192 dxn(i)=(xuse(i)-xpkg(iv(i),1))
194 dxn(i)=(xuse(i)-xpkg(iv(i),1))*havi
205 if((xpkg(iprev,1).le.xuse(i)).and.
206 > (xuse(i).le.xpkg(iprev+1,1)))
then
209 iv(i)=1+havi*(xuse(i)-xpkg(1,1))
210 iv(i)=max(1,min((nx-1),iv(i)))
214 dxn(i)=(xuse(i)-xpkg(iv(i),1))
216 dxn(i)=(xuse(i)-xpkg(iv(i),1))*havi
231 allocate(iok(inum),stat=istat)
234 write(6,*)
' ?? xlookup "iok" vector allocation failure!'
239 allocate(imina(inum),stat=istat)
242 write(6,*)
' ?? xlookup "imina" vector allocation failure!'
246 allocate(imaxa(inum),stat=istat)
249 write(6,*)
' ?? xlookup "imaxa" vector allocation failure!'
259 if(init_guess.eq.0)
then
260 go to (100,200,300) ialg
262 go to (150,250,350) ialg
279 iv(i)=1+havi*(xuse(i)-xpkg(1,1))
282 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
284 if(xuse(i).lt.xpkg(iv(i),1))
then
287 imaxa(i)=max(1,(iv(i)-1))
288 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
291 imina(i)=min((nx-1),(iv(i)+1))
303 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
304 if(i_sign*zdelta.le.hloc)
then
314 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
316 if(xuse(i).lt.xpkg(iv(i),1))
then
319 imaxa(i)=max(1,(iv(i)-1))
320 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
323 imina(i)=min((nx-1),(iv(i)+1))
333 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
334 if(i_sign*zdelta.le.hloc)
then
344 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
346 if(xuse(i).lt.xpkg(iv(i),1))
then
349 imaxa(i)=max(1,(iv(i)-1))
350 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
353 imina(i)=min((nx-1),(iv(i)+1))
363 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
364 if(i_sign*zdelta.le.hloc)
then
374 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
376 if(xuse(i).lt.xpkg(iv(i),1))
then
379 imaxa(i)=max(1,(iv(i)-1))
380 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
383 imina(i)=min((nx-1),(iv(i)+1))
393 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
394 if(i_sign*zdelta.le.hloc)
then
404 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
406 if(xuse(i).lt.xpkg(iv(i),1))
then
409 imaxa(i)=max(1,(iv(i)-1))
410 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
413 imina(i)=min((nx-1),(iv(i)+1))
418 if(.not.iok(i)) iprob=iprob+1
443 if((xpkg(iprev,1).le.xuse(i)).and.
444 > (xuse(i).le.xpkg(iprev+1,1)))
then
452 iv(i)=1+havi*(xuse(i)-xpkg(1,1))
453 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
455 if(xuse(i).lt.xpkg(iv(i),1))
then
458 imaxa(i)=max(1,(iv(i)-1))
459 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
462 imina(i)=min((nx-1),(iv(i)+1))
475 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
476 if(i_sign*zdelta.le.hloc)
then
484 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
486 if(xuse(i).lt.xpkg(iv(i),1))
then
489 imaxa(i)=max(1,(iv(i)-1))
490 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
493 imina(i)=min((nx-1),(iv(i)+1))
503 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
504 if(i_sign*zdelta.le.hloc)
then
512 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
514 if(xuse(i).lt.xpkg(iv(i),1))
then
517 imaxa(i)=max(1,(iv(i)-1))
518 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
521 imina(i)=min((nx-1),(iv(i)+1))
531 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
532 if(i_sign*zdelta.le.hloc)
then
540 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
542 if(xuse(i).lt.xpkg(iv(i),1))
then
545 imaxa(i)=max(1,(iv(i)-1))
546 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
549 imina(i)=min((nx-1),(iv(i)+1))
559 zdelta=(xuse(i)-xpkg(iv(i)+isrch,1))
560 if(i_sign*zdelta.le.hloc)
then
568 iv(i)=max(imina(i),min(imaxa(i),iv(i)))
570 if(xuse(i).lt.xpkg(iv(i),1))
then
573 imaxa(i)=max(1,(iv(i)-1))
574 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
577 imina(i)=min((nx-1),(iv(i)+1))
582 if(.not.iok(i)) iprob=iprob+1
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))
626 iv(i)=(imina(i)+imaxa(i))/2
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))
640 iv(i)=(imina(i)+imaxa(i))/2
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))
654 iv(i)=(imina(i)+imaxa(i))/2
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))
668 iv(i)=(imina(i)+imaxa(i))/2
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))
678 if(.not.iok(i)) iprob=iprob+1
702 if((xpkg(iprev,1).le.xuse(i)).and.
703 > (xuse(i).le.xpkg(iprev+1,1)))
then
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))
728 iv(i)=(imina(i)+imaxa(i))/2
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))
742 iv(i)=(imina(i)+imaxa(i))/2
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))
756 iv(i)=(imina(i)+imaxa(i))/2
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))
770 iv(i)=(imina(i)+imaxa(i))/2
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))
780 if(.not.iok(i)) iprob=iprob+1
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
813 if(iv(i).lt.nx-1)
then
814 zdindx=xpkg(iv(i)+1,2)-zindx0
818 zindex=zindx0+xfac*zdindx
821 iv(i)=max(1,min(nx-1,iv(i)))
823 if(xuse(i).lt.xpkg(iv(i),1))
then
825 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
837 if(xuse(i).lt.xpkg(iv(i),1))
then
839 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
849 if(xuse(i).lt.xpkg(iv(i),1))
then
851 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
861 if(xuse(i).lt.xpkg(iv(i),1))
then
863 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
873 if(xuse(i).lt.xpkg(iv(i),1))
then
875 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
881 if(.not.iok(i)) iprob=iprob+1
905 if((xpkg(iprev,1).le.xuse(i)).and.
906 > (xuse(i).le.xpkg(iprev+1,1)))
then
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
919 if(iv(i).lt.nx-1)
then
920 zdindx=xpkg(iv(i)+1,2)-zindx0
924 zindex=zindx0+xfac*zdindx
927 iv(i)=max(1,min(nx-1,iv(i)))
929 if(xuse(i).lt.xpkg(iv(i),1))
then
931 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
944 if(xuse(i).lt.xpkg(iv(i),1))
then
946 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
956 if(xuse(i).lt.xpkg(iv(i),1))
then
958 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
968 if(xuse(i).lt.xpkg(iv(i),1))
then
970 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
980 if(xuse(i).lt.xpkg(iv(i),1))
then
982 else if(xuse(i).gt.xpkg(iv(i)+1,1))
then
988 if(.not.iok(i)) iprob=iprob+1
1008 if(iprob.gt.0)
go to 5
1013 dxn=(xuse-xpkg(iv,1))
1014 else if(ialg.ne.3)
then
1015 dxn=(xuse-xpkg(iv,1))*xpkg(iv,3)
1019 dxn=(xuse-xpkg(iv,1))*xpkg(iv,3)
1020 hv=xpkg(iv+1,1)-xpkg(iv,1)