lagp.f
555 lines
| 17.8 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine lagp(plag,wl,r0,dr,nl,nrange) | |||
| c | ||||
| c predict lagged product matrix using state parameters stored | ||||
| c as cubic b-splines | ||||
| c read weights for numeric integration from a file | ||||
| c | ||||
| c faster version - assumes samples taken at half-integer lags, ranges | ||||
| c | ||||
| include 'fitter.h' | ||||
| include 'spline.h' | ||||
| include 'bfield.h' | ||||
| integer nl,nrange | ||||
| real plag(nl,nrange),taus(nl) | ||||
| parameter(nw=1000) | ||||
| integer ilag(nw) | ||||
| real weight(nw),dtau(nw),drange(nw) | ||||
| logical first | ||||
| real store(2*nl,2*nrange,2) | ||||
| data first/.true./ | ||||
| c | ||||
| c write(*,*) "dens_before: ",dens | ||||
| c call exit | ||||
| c zero out array | ||||
| c write(*,*) "Starting lagp" | ||||
| do i=1,nrange | ||||
| do j=1,nl | ||||
| plag(j,i)=0.0 | ||||
| end do | ||||
| end do | ||||
| c | ||||
| c compute all required acfs | ||||
| c | ||||
| do i=1,2*nrange | ||||
| ii=(i-1)/2+1 | ||||
| alt=r0+float(i-1)*dr*0.5 ! half integer ranges | ||||
| c write(*,*) "dens_before: ",dens | ||||
| c call exit | ||||
| c write(*,*) "alt: ",alt | ||||
| c write(*,*) "dens: ",dens | ||||
| c write(*,*) "te: ",te | ||||
| c write(*,*) "ti1: ",ti1 | ||||
| c write(*,*) "hf: ",hf | ||||
| c write(*,*) "hef: ",hef | ||||
| c call exit | ||||
| call get_spline(alt,dens,te,ti1,hf,hef) | ||||
| c write(*,*) "dens_after: ",dens | ||||
| c call exit | ||||
| do k=1,nion | ||||
| ti(k)=ti1 | ||||
| end do | ||||
| fi(2)=hf | ||||
| fi(3)=hef | ||||
| fi(1)=1.0-hf-hef | ||||
| do j=1,2*nl | ||||
| tau=float(j-1)*(dr/1.5e5)*0.5 ! half integer lags | ||||
| c just consider a single azimuth angle for now | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho) | ||||
| c store(j,i,1)=rho*dens*(100.0/alt)**2 | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho) | ||||
| c store(j,i,2)=rho*dens*(100.0/alt)**2 | ||||
| call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| & alpha_prof(ii),dens,bfld_prof(ii),rho) | ||||
| store(j,i,1)=rho*dens*(100.0/alt)**2 | ||||
| c write(*,*) "store", store(j,i,1) | ||||
| c call exit | ||||
| end do | ||||
| end do | ||||
| c construct lagged products | ||||
| c | ||||
| c write(*,*) "before weights" | ||||
| c write(*,*) first | ||||
| c write(*,*) "weight before: ",weight | ||||
| if(first) then | ||||
| open(unit=25,file='/usr/local/lib/faraday/weights.dat', | ||||
| & status='old') | ||||
| read (25,*) nwt | ||||
| do i=1,nwt | ||||
| read(25,*) weight(i),dtau(i),drange(i),ilag(i) | ||||
| c write(*,*) "ilag: ",ilag(i) | ||||
| end do | ||||
| c flush(25) | ||||
| close(25) | ||||
| c first=.false. | ||||
| end if | ||||
| c write(*,*) "after weights" | ||||
| c write(*,*) "weight(3): ",weight(3) | ||||
| c call exit | ||||
| c write(*,*) "nl: ",nl,"nrange: ",nrange,"nwt: ",nwt | ||||
| do i=1+nl,nrange ! don't worry about redundancy here | ||||
| do j=1,nwt | ||||
| c write(*,*) "j: ",j | ||||
| alt=r0+(float(i-1)+drange(j))*dr | ||||
| tau=abs(dtau(j)*dr/1.5e5) | ||||
| lag=ilag(j)+1 | ||||
| c write(*,*) tau | ||||
| k=1+2.0*(tau/(dr/1.5e5)) | ||||
| l=1+2.0*(alt-r0)/dr | ||||
| do m=1,1 ! 2 | ||||
| c call exit | ||||
| c write(*,*) i | ||||
| c call exit | ||||
| c if (i .eq. nrange) then | ||||
| c if (j .eq. 227) then | ||||
| c write(*,*) "lag: ",lag | ||||
| c write(*,*) "ilag: ",ilag(j) | ||||
| c call exit | ||||
| c write(*,*) "plag: ",plag(lag,i) | ||||
| c write(*,*) "weight: ",weight | ||||
| c | ||||
| c write(*,*) "weight: ",weight(j) | ||||
| c write(*,*) "store: ",store(k,l,m) | ||||
| c call exit | ||||
| c end if | ||||
| c end if | ||||
| plag(lag,i)=plag(lag,i)+weight(j)*store(k,l,m) | ||||
| c if (i .eq. nrange) then | ||||
| c write(*,*) "lag: ",lag | ||||
| c write(*,*) "plag: ",plag(lag,i) | ||||
| c write(*,*) "plag: ",plag(lag,i) | ||||
| c call exit | ||||
| c end if | ||||
| end do | ||||
| c write(*,*) "plag: ",plag(lag,i) | ||||
| c write(*,*) "plag: ",plag(12,72) | ||||
| c call exit | ||||
| end do | ||||
| c call exit | ||||
| c write(*,*) "plag: ",plag(12,72) | ||||
| end do | ||||
| c write(*,*) "plag: ",plag(12,72) | ||||
| c call exit | ||||
| c write(*,*) "End LAGP" | ||||
| return | ||||
| end | ||||
| subroutine lagp_old(plag,wl,r0,dr,nl,nrange) | ||||
| c | ||||
| c predict lagged product matrix using state parameters stored | ||||
| c as cubic b-splines | ||||
| c read weights for numeric integration from a file | ||||
| c | ||||
| c general version - samples can be anywhere | ||||
| c | ||||
| include 'fitter.h' | ||||
| include 'spline.h' | ||||
| include 'bfield.h' | ||||
| integer nl,nrange | ||||
| real plag(nl,nrange),taus(nl) | ||||
| parameter(nw=1000) | ||||
| integer ilag(nw) | ||||
| real weight(nw),dtau(nw),drange(nw) | ||||
| logical first | ||||
| data first/.true./ | ||||
| c | ||||
| c zero out array | ||||
| do i=1,nrange | ||||
| do j=1,nl | ||||
| plag(j,i)=0.0 | ||||
| end do | ||||
| end do | ||||
| c | ||||
| c construct lagged products | ||||
| c | ||||
| if(first) then | ||||
| open(unit=25,file='weights.dat',status='old') | ||||
| read (25,*) nwt | ||||
| do i=1,nwt | ||||
| read(25,*) weight(i),dtau(i),drange(i),ilag(i) | ||||
| end do | ||||
| close(25) | ||||
| first=.false. | ||||
| end if | ||||
| c can avoid recalculating splines and/or ACFs when parameters repeat | ||||
| altp=-1.0 | ||||
| taup=-1.0 | ||||
| do i=1+nl,nrange | ||||
| do j=1,nwt | ||||
| alt=r0+(float(i-1)+drange(j))*dr | ||||
| ii=(alt-r0)/dr+1 | ||||
| if(alt.ne.altp) then | ||||
| call get_spline(alt,dens,te,ti1,hf,hef) | ||||
| c write(*,*) alt,dens,te,ti1,hf,hef | ||||
| do k=1,nion | ||||
| ti(k)=ti1 | ||||
| end do | ||||
| fi(2)=hf | ||||
| fi(3)=hef | ||||
| fi(1)=1.0-hf-hef | ||||
| end if | ||||
| tau=abs(dtau(j)*dr/1.5e5) | ||||
| lag=ilag(j)+1 | ||||
| c write(*,*) i,alpha_prof(i),bfld_prof(i) | ||||
| if(tau.ne.taup.or.alt.ne.altp) then | ||||
| c two- or three- point Gauss Hermite quadrature rule | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho1) | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho2) | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii),dens,bfld_prof(ii),rho1) | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)+1.9e-2,dens,bfld_prof(ii),rho2) | ||||
| c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| c & alpha_prof(ii)-1.9e-2,dens,bfld_prof(ii),rho3) | ||||
| call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion, | ||||
| & alpha_prof(ii),dens,bfld_prof(ii),rho) | ||||
| c rho=(rho1+rho2)/2.0 | ||||
| c rho=(rho2+rho3)*0.29541+rho1*1.18164)/1.77246 | ||||
| end if | ||||
| plag(lag,i)=plag(lag,i)+weight(j)*rho*dens*(100.0/alt)**2 | ||||
| altp=alt | ||||
| taup=tau | ||||
| end do | ||||
| end do | ||||
| return | ||||
| end | ||||
| function atanh(x) | ||||
| c | ||||
| real atanh,x | ||||
| c | ||||
| atanh=log(sqrt((1.+x)/(1.-x))) | ||||
| return | ||||
| end | ||||
| subroutine get_spline(alt,dens,te,ti,hf,hef) | ||||
| c | ||||
| c routines for handling cubic b-spline interpolation | ||||
| c gets values consistent with stored coefficients | ||||
| c | ||||
| include 'spline.h' | ||||
| c | ||||
| c bspline values from 200-1500 km in 15-km intervals | ||||
| c must specify five knots above ceiling | ||||
| c note offset ... start accessing splines above bottom two knots | ||||
| c | ||||
| c five banks of splines, for Ne, Te, Ti, H+, and He+ versus altitude | ||||
| c need to initialize ta somewhere | ||||
| c | ||||
| c write(*,*) "ta: ",ta | ||||
| c write(*,*) "bcoef(1,1): ",bcoef(1,1) | ||||
| c write(*,*) "nspline: ",nspline | ||||
| c write(*,*) "norder: ",norder | ||||
| c write(*,*) "alt: ",alt | ||||
| dens=bvalue(ta,bcoef(1,1),nspline,norder,alt,0) | ||||
| c write(*,*) "dens_inside_get_spline: ",dens | ||||
| c write(*,*) dens | ||||
| c call exit | ||||
| dens=10.0**MAX(dens,2.0) | ||||
| c write(*,*) dens | ||||
| c call exit | ||||
| te=bvalue(ta,bcoef(1,2),nspline,norder,alt,0) | ||||
| te=t0+t1*(1.0+tanh(te))/2.0 ! DLH 10/14 was 3500 | ||||
| c tr=bvalue(ta,bcoef(1,3),nspline,norder,alt,0) | ||||
| c tr=exp(tr) | ||||
| c ti=te/tr | ||||
| ti=bvalue(ta,bcoef(1,3),nspline,norder,alt,0) | ||||
| ti=t0+t1*(1.0+tanh(ti))/2.0 ! DLH 10/14 was 3500 | ||||
| hf=bvalue(ta,bcoef(1,4),nspline,norder,alt,0) | ||||
| hf=(1.0+tanh(hf))/2.0 | ||||
| hef=bvalue(ta,bcoef(1,5),nspline,norder,alt,0) | ||||
| hef=(1.0+tanh(hef))/2.0 | ||||
| return | ||||
| end | ||||
| real function bvalue ( t, bcoef, n, k, x, jderiv ) | ||||
| c from * a practical guide to splines * by c. de boor | ||||
| calls interv | ||||
| c | ||||
| calculates value at x of jderiv-th derivative of spline from b-repr. | ||||
| c the spline is taken to be continuous from the right, EXCEPT at the | ||||
| c rightmost knot, where it is taken to be continuous from the left. | ||||
| c | ||||
| c****** i n p u t ****** | ||||
| c t, bcoef, n, k......forms the b-representation of the spline f to | ||||
| c be evaluated. specifically, | ||||
| c t.....knot sequence, of length n+k, assumed nondecreasing. | ||||
| c bcoef.....b-coefficient sequence, of length n . | ||||
| c n.....length of bcoef and dimension of spline(k,t), | ||||
| c a s s u m e d positive . | ||||
| c k.....order of the spline . | ||||
| c | ||||
| c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed | ||||
| c arbitrarily by the dimension statement for aj, dl, dr below, | ||||
| c but is n o w h e r e c h e c k e d for. | ||||
| c | ||||
| c x.....the point at which to evaluate . | ||||
| c jderiv.....integer giving the order of the derivative to be evaluated | ||||
| c a s s u m e d to be zero or positive. | ||||
| c | ||||
| c****** o u t p u t ****** | ||||
| c bvalue.....the value of the (jderiv)-th derivative of f at x . | ||||
| c | ||||
| c****** m e t h o d ****** | ||||
| c The nontrivial knot interval (t(i),t(i+1)) containing x is lo- | ||||
| c cated with the aid of interv . The k b-coeffs of f relevant for | ||||
| c this interval are then obtained from bcoef (or taken to be zero if | ||||
| c not explicitly available) and are then differenced jderiv times to | ||||
| c obtain the b-coeffs of (d**jderiv)f relevant for that interval. | ||||
| c Precisely, with j = jderiv, we have from x.(12) of the text that | ||||
| c | ||||
| c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) ) | ||||
| c | ||||
| c where | ||||
| c / bcoef(.), , j .eq. 0 | ||||
| c / | ||||
| c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1) | ||||
| c / ----------------------------- , j .gt. 0 | ||||
| c / (t(.+k-j) - t(.))/(k-j) | ||||
| c | ||||
| c Then, we use repeatedly the fact that | ||||
| c | ||||
| c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) ) | ||||
| c with | ||||
| c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1) | ||||
| c a(.,x) = --------------------------------------- | ||||
| c (x - t(.)) + (t(.+m-1) - x) | ||||
| c | ||||
| c to write (d**j)f(x) eventually as a linear combination of b-splines | ||||
| c of order 1 , and the coefficient for b(i,1,t)(x) must then be the | ||||
| c desired number (d**j)f(x). (see x.(17)-(19) of text). | ||||
| c | ||||
| parameter (kmax = 40) | ||||
| integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1 | ||||
| * ,mflag,nmi,jdrvp1 | ||||
| c integer kmax | ||||
| C real bcoef(n),t(1),x, aj(20),dl(20),dr(20),fkmj | ||||
| real bcoef(n),x, aj(kmax),dl(kmax),dr(kmax),fkmj | ||||
| dimension t(n+k) | ||||
| c former fortran standard made it impossible to specify the length of t | ||||
| c precisely without the introduction of otherwise superfluous addition- | ||||
| c al arguments. | ||||
| bvalue = 0. | ||||
| if (jderiv .ge. k) go to 99 | ||||
| c | ||||
| c *** Find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and | ||||
| c t(i) .le. x .lt. t(i+1) . If no such i can be found, x lies | ||||
| c outside the support of the spline f , hence bvalue = 0. | ||||
| c (The asymmetry in this choice of i makes f rightcontinuous, except | ||||
| c at t(n+k) where it is leftcontinuous.) | ||||
| call interv ( t, n+k, x, i, mflag ) | ||||
| if (mflag .ne. 0) go to 99 | ||||
| c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i). | ||||
| km1 = k - 1 | ||||
| if (km1 .gt. 0) go to 1 | ||||
| bvalue = bcoef(i) | ||||
| go to 99 | ||||
| c | ||||
| c *** store the k b-spline coefficients relevant for the knot interval | ||||
| c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j), | ||||
| c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable | ||||
| c from input to zero. set any t.s not obtainable equal to t(1) or | ||||
| c to t(n+k) appropriately. | ||||
| 1 jcmin = 1 | ||||
| imk = i - k | ||||
| if (imk .ge. 0) go to 8 | ||||
| jcmin = 1 - imk | ||||
| do 5 j=1,i | ||||
| 5 dl(j) = x - t(i+1-j) | ||||
| do 6 j=i,km1 | ||||
| aj(k-j) = 0. | ||||
| 6 dl(j) = dl(i) | ||||
| go to 10 | ||||
| 8 do 9 j=1,km1 | ||||
| 9 dl(j) = x - t(i+1-j) | ||||
| c | ||||
| 10 jcmax = k | ||||
| nmi = n - i | ||||
| if (nmi .ge. 0) go to 18 | ||||
| jcmax = k + nmi | ||||
| do 15 j=1,jcmax | ||||
| 15 dr(j) = t(i+j) - x | ||||
| do 16 j=jcmax,km1 | ||||
| aj(j+1) = 0. | ||||
| 16 dr(j) = dr(jcmax) | ||||
| go to 20 | ||||
| 18 do 19 j=1,km1 | ||||
| 19 dr(j) = t(i+j) - x | ||||
| c | ||||
| 20 do 21 jc=jcmin,jcmax | ||||
| 21 aj(jc) = bcoef(imk + jc) | ||||
| c | ||||
| c *** difference the coefficients jderiv times. | ||||
| if (jderiv .eq. 0) go to 30 | ||||
| do 23 j=1,jderiv | ||||
| kmj = k-j | ||||
| fkmj = float(kmj) | ||||
| ilo = kmj | ||||
| do 23 jj=1,kmj | ||||
| aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj | ||||
| 23 ilo = ilo - 1 | ||||
| c | ||||
| c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative, | ||||
| c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv). | ||||
| 30 if (jderiv .eq. km1) go to 39 | ||||
| jdrvp1 = jderiv + 1 | ||||
| do 33 j=jdrvp1,km1 | ||||
| kmj = k-j | ||||
| ilo = kmj | ||||
| do 33 jj=1,kmj | ||||
| aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj)) | ||||
| 33 ilo = ilo - 1 | ||||
| 39 bvalue = aj(1) | ||||
| c | ||||
| 99 return | ||||
| end | ||||
| subroutine interv ( xt, lxt, x, left, mflag ) | ||||
| c from * a practical guide to splines * by C. de Boor | ||||
| computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) . | ||||
| c | ||||
| c****** i n p u t ****** | ||||
| c xt.....a real sequence, of length lxt , assumed to be nondecreasing | ||||
| c lxt.....number of terms in the sequence xt . | ||||
| c x.....the point whose location with respect to the sequence xt is | ||||
| c to be determined. | ||||
| c | ||||
| c****** o u t p u t ****** | ||||
| c left, mflag.....both integers, whose value is | ||||
| c | ||||
| c 1 -1 if x .lt. xt(1) | ||||
| c i 0 if xt(i) .le. x .lt. xt(i+1) | ||||
| c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt) | ||||
| c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x | ||||
| c | ||||
| c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0 | ||||
| c indicates that x lies outside the CLOSED interval | ||||
| c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the | ||||
| c intervals is due to the decision to make all pp functions cont- | ||||
| c inuous from the right, but, by returning mflag = 0 even if | ||||
| C x = xt(lxt), there is the option of having the computed pp function | ||||
| c continuous from the left at xt(lxt) . | ||||
| c | ||||
| c****** m e t h o d ****** | ||||
| c The program is designed to be efficient in the common situation that | ||||
| c it is called repeatedly, with x taken from an increasing or decrea- | ||||
| c sing sequence. This will happen, e.g., when a pp function is to be | ||||
| c graphed. The first guess for left is therefore taken to be the val- | ||||
| c ue returned at the previous call and stored in the l o c a l varia- | ||||
| c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec- | ||||
| c essary since the present call may have nothing to do with the previ- | ||||
| c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left = | ||||
| c ilo and are done after just three comparisons. | ||||
| c Otherwise, we repeatedly double the difference istep = ihi - ilo | ||||
| c while also moving ilo and ihi in the direction of x , until | ||||
| c xt(ilo) .le. x .lt. xt(ihi) , | ||||
| c after which we use bisection to get, in addition, ilo+1 = ihi . | ||||
| c left = ilo is then returned. | ||||
| c | ||||
| integer left,lxt,mflag, ihi,ilo,istep,middle | ||||
| real x,xt(lxt) | ||||
| data ilo /1/ | ||||
| save ilo | ||||
| ihi = ilo + 1 | ||||
| if (ihi .lt. lxt) go to 20 | ||||
| if (x .ge. xt(lxt)) go to 110 | ||||
| if (lxt .le. 1) go to 90 | ||||
| ilo = lxt - 1 | ||||
| ihi = lxt | ||||
| c | ||||
| 20 if (x .ge. xt(ihi)) go to 40 | ||||
| if (x .ge. xt(ilo)) go to 100 | ||||
| c | ||||
| c **** now x .lt. xt(ilo) . decrease ilo to capture x . | ||||
| istep = 1 | ||||
| 31 ihi = ilo | ||||
| ilo = ihi - istep | ||||
| if (ilo .le. 1) go to 35 | ||||
| if (x .ge. xt(ilo)) go to 50 | ||||
| istep = istep*2 | ||||
| go to 31 | ||||
| 35 ilo = 1 | ||||
| if (x .lt. xt(1)) go to 90 | ||||
| go to 50 | ||||
| c **** now x .ge. xt(ihi) . increase ihi to capture x . | ||||
| 40 istep = 1 | ||||
| 41 ilo = ihi | ||||
| ihi = ilo + istep | ||||
| if (ihi .ge. lxt) go to 45 | ||||
| if (x .lt. xt(ihi)) go to 50 | ||||
| istep = istep*2 | ||||
| go to 41 | ||||
| 45 if (x .ge. xt(lxt)) go to 110 | ||||
| ihi = lxt | ||||
| c | ||||
| c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval. | ||||
| 50 middle = (ilo + ihi)/2 | ||||
| if (middle .eq. ilo) go to 100 | ||||
| c note. it is assumed that middle = ilo in case ihi = ilo+1 . | ||||
| if (x .lt. xt(middle)) go to 53 | ||||
| ilo = middle | ||||
| go to 50 | ||||
| 53 ihi = middle | ||||
| go to 50 | ||||
| c**** set output and return. | ||||
| 90 mflag = -1 | ||||
| left = 1 | ||||
| return | ||||
| 100 mflag = 0 | ||||
| left = ilo | ||||
| return | ||||
| 110 mflag = 1 | ||||
| if (x .eq. xt(lxt)) mflag = 0 | ||||
| left = lxt | ||||
| 111 if (left .eq. 1) return | ||||
| left = left - 1 | ||||
| if (xt(left) .lt. xt(lxt)) return | ||||
| go to 111 | ||||
| end | ||||
