program mtmsvdconf c c Copyright (C) 1995 c c Michael Mann c c determine signficance levels for LFV spectrum (and for c 2nd eigenvalue spectrum also) based on bootstrap c resampling calculations. c c v1.0 c c uses stripped version of original multitaper spectral code of c Jeffrey Park (c) c c mmax - maximum number of spatial dofs desired (ie, max # of time series) c ntapmax - maximum number of spectral dofs desired c ifmax - maximum padded DFT length (maximum=8192) c nscanmax - maximum length of time series c nitermax - maximum number of bootstrap resamples performed c c NOTE THE RANDOM NUMBER GENERATOR "RAND" USED CAN BE MACHINE c SPECIFIC, AND MAY NEED TO BE REPLACED BY AN ALTERNATIVE c RANDOM NUMBER GENERATOR c parameter (mmax=25,ntapmax=3,ifmax=2048,nscanmax=1200, $ nitermax=2000) c character*19 filename character*1 cfnum complex *8 AAA(mmax,ntapmax,ifmax) complex *16 AA(mmax,ntapmax),UU(mmax,mmax),VV(ntapmax,ntapmax) real*8 ell real *8 S(ntapmax),sum real *8 partsum(ntapmax),sum99(ntapmax),sum95(ntapmax), $ sum90(ntapmax),sum50(ntapmax) real partvar(nitermax,ntapmax,ifmax) real ar,ai real weight(mmax),a(nscanmax,mmax),anew(nscanmax,mmax) real datave(mmax),sd(mmax) integer iplace(nscanmax) common/npiprol/anpi common/nnaamm/iwflag common/work/b1(8194) common/taperz/ta(4100,16,2),tai(4100,16,2),dcf(4100,16,2), $ amu(4100,4) common/stap2/tas(8192,16) common/envel1/ell(20) c c iwflag=0 zero = 0.0 one = 1.0 c c c dt = time interval spacing in units of years c nscan = length of time series in units of "dt" c iabv = number of time series in analysis (ie, "spatial" array) c c presently set for 25 records of 1 month resolution c dt = 0.0833333 nscan = 1200 iabv = 25 c c f1 = minimum frequency for analysis = 0 c f2 = maximum frequency (in cyc/year -- this can be varied up to c the nyquist sampling frequency = 0.5/dt) c c ( presently set for the range 0 to 0.5 cycle year, or, periods greater c than tau=2 years. This is a smaller fraction of the full Nyquist c interval f=0 to f=6.0 cycle/year for monthly sampling) c f1 = 0.0 f2 = 0.5 c c default settings c npi = 2 anpi=float(npi) nwin=3 itype = 0 nbeg=1 nend = nscan npts = nscan niter = 1000 iseed = 999 isamp = 0 444 write (6,*) 'present settings' write (6,*) '1) use',nwin, ' x ',npi,' pi tapers' write (6,*) '2) data interval (',nbeg,' to ',nend,' )' write (6,*) '3) perform ',niter,' permutation realizations' write (6,*) '4) reinitialize random number seed' write (6,*) ' (presently, iseed=',iseed,')' if (isamp.eq.0) write (6,*) '5) resampling: uniform' if (isamp.eq.1) write (6,*) '5) resampling: block=',lblock write (6,*) 'item to change (0=continue)' read (5,*) icont if (icont.eq.1) then goto 1111 else if (icont.eq.2) then goto 2222 else if (icont.eq.3) then goto 3333 else if (icont.eq.4) then goto 4444 else if (icont.eq.5) then goto 5555 else goto 9999 endif endif endif endif endif 1111 write (6,*) 'time-frequency bandwidth product "p"' read (5,*) npinew if (ntapmax.lt.2*npinew-1) then write (6,*) 'sorry--exceeds maximum bandwidth' else write (6,*) 'number of tapers' read (5,*) nwinnew if (nwinnew.gt.2*npinew-1) then write (6,*) 'sorry--cannot exceed 2*p-1' else npi=npinew anpi = float(npi) nwin=nwinnew endif endif goto 444 2222 write (6,*) 'beginning time index' read (5,*) nbeg write (6,*) 'ending time index' read (5,*) nend nmove = (nend-nbeg+1)*dt npts = nend-nbeg+1 goto 444 3333 write (6,*) 'number of realizations' read (5,*) niter if (niter.gt.nitermax) then write (6,*) 'sorry-too many' goto 3333 endif goto 444 4444 write (6,*) 'seed (0-1000)?' read (5,*) iseed goto 444 5555 write (6,*) 'resampling:' write (6,*) '(0) uniform' write (6,*) '(1) block' read (5,*) isamp if (isamp.eq.1) then write (6,*) 'block (# samples)' read (5,*) lblock endif goto 444 c c npad = zero padding (power of two--must be > nscan) c c nfstep determines the frequency resolution at which the SVD c is performed c 9999 npad=2048 nfstep = 1 c c calculate frequency spacing, etc. c ddf=1./(npad*dt) eps=.1*ddf nf1=(f1+eps)/ddf nf2=(f2+eps)/ddf nf=nf2-nf1+1 nfmax = nf-1 c c set spatial, spectral dimensions for SVD matrix c M = iabv N = nwin NV = N NU = N IP = 0 c c read in data. c c * THIS SECTION WILL HAVE TO MODIFIED DEPENDING ON THE DATA c ANALYZED * c open (unit=9,file='synth1.dat', $ status='old') c open (unit=10,file='synth2.dat', $ status='old') c open (unit=11,file='synth3.dat', $ status='old') c open (unit=12,file='synth4.dat', $ status='old') c open (unit=13,file='synth5.dat', $ status='old') c c read in the iabv=25 gridpoints from the c five files of 5 gridpoint series/file c c do j=1,nscan read (9,*) adum,(a(j,jnode),jnode=1,5) read (10,*) adum,(a(j,jnode),jnode=6,10) read (11,*) adum,(a(j,jnode),jnode=11,15) read (12,*) adum,(a(j,jnode),jnode=16,20) read (13,*) adum,(a(j,jnode),jnode=21,25) end do close (unit=9) close (unit=10) close (unit=11) close (unit=12) close (unit=13) c c c write (6,*) 'read in: ' write (6,*) iabv, ' time series' write (6,*) nscan*dt, ' years of data...' c c for each time series, subtract mean, normalize by standard deviation, c also set weights for the different dataseries. c ( for the the present case, all dataseries are equally weighted for c analysis) c c open (unit=15,file='svd_permute-stats.out',status='unknown') do j=1,iabv weight(j)=1.0 sum=0.0 do i=nbeg,nend sum=sum+a(i,j) end do datave(j)=sum/float(npts) var=0.0 inew = 0 do i=nbeg,nend inew = inew+1 anew(inew,j)=a(i,j)-datave(j) var=var + anew(inew,j)**2 end do sd(j)=sqrt(var/float(npts)) do i=1,npts anew(i,j)=weight(j)*anew(i,j)/sd(j) end do write (15,*) j,sd(j),datave(j) end do close (unit=15) c c open output files c do i=1,nwin-1 cfnum = char(48+i) filename = 'evals-permute'//cfnum//'.out' open (unit=13+i-1,file=filename, $ status='unknown') end do c c construct tapers appropriate for window width c call taper(npts,nwin,ell) c write (6,*) 'nscan ',nscan write (6,*) 'npts ',npts write (6,*) 'nbeg ',nbeg write (6,*) 'nend ',nend write (6,*) 'dt ',dt c c irep = 2*iseed*rand(0)+1 c c reseed the random number generator c do i=1,irep adum = rand(0) end do n1 = 1 n2 = npts c c begin permutation resampling procedure c do it=1,niter c c permute the series randomly in time, keeping spatial structure c intact. Note that any time-domain structure is destroyed here, c some of which may be due to short-range serial correlation. c c generate randomized sequence of the numbers 1 to npts c do icount=1,npts 232 ipl = rand(0)*npts + 1 do j=1,icount-1 if (ipl.eq.iplace(j)) goto 232 end do iplace(icount)=ipl end do c c assign this new randomized sequence c do i=1,npts do j=1,iabv a(i,j)=anew(iplace(i),j) end do end do c c calculate DFTs for each gridpoint for new sequence c do ii=1,iabv do iwin=1,nwin do i=1,npad b1(i)=0. end do do i=n1,n2 b1(i-n1+1)=a(i,ii)*tas(i+1-n1,iwin) end do call realft(b1,npad/2,1) b1(npad-1)=b1(2) b1(npad)=0. b1(2)=0. do iif=1,nfmax iwhere = 2*(nf1+iif)-1 ar=b1(iwhere) ai=b1(iwhere+1) AAA(ii,iwin,iif)=cmplx(ar,ai) end do end do c c terminate DFT calculation c end do c c perform SVD for each frequency c do iif =1,nfmax,nfstep do ii=1,iabv do iwin=1,nwin AA(ii,iwin)=1.d0*AAA(ii,iwin,iif) end do end do c c perform SVD c call CSVD(AA,mmax,nwin,M,N,IP,NU,NV,S,UU,VV) c c determine local fractional variances, etc. c sum = 0.d0 do i=1,nwin-1 sum = sum + S(i)**2 partsum(i)=0.0 do j=i,nwin partsum(i)=partsum(i)+S(j)**2 end do end do do i=1,nwin-1 partvar(it,i,iif) = real(S(i)**2/partsum(i)) end do c c terminate frequency loop c end do c c provide progress update c open (unit=12,file='svdevals-permute.status', $ status='unknown') if ((niter.ge.10).and.(mod(it-1,niter/10).eq.0)) then write (12,*) it,' realizations of ',niter,' done...' write (6,*) it,' realizations of ',niter,' done...' endif close (unit=12) end do c c now sort -- calculate 99%,95%,90%, and 50% quantiles c for local fractional variance as function c of frequency and mode (1st mode corresponds to "local fractional c variance spectrum") c c we need to resolve the secular band (f< p/tau where p = time-frequency c bandwidth product (e.g., p=2 for 2pi tapers) and tau is length of c data series in years since the number of spectral degrees of freedom c varies from nwin-1 to nwin over that band, owing to vanishing of c complex part of spectrum at zero frequency. Outside the band, spectral c degrees of freedom (complex) are nwin and each independent c frequency of the DFT c provides an independent sample which we can use to calculate the c confidence levels with a factor of roughly nscan/2 greater precision. c c We average the estimates within and outside the secular band c separately c freq_sec = anpi/(npts*dt)+eps iavefreq = 0 iflag = 0 do j=1,nwin-1 sum99(j)=0.d0 sum95(j)=0.d0 sum90(j)=0.d0 sum50(j)=0.d0 end do iif0 = 1 do iif=1,nfmax,nfstep freq = (iif-1)*ddf+f1 do j=1,nwin-1 do i=1,niter do i2=i+1,niter if (partvar(i,j,iif).gt.partvar(i2,j,iif)) then dum = partvar(i,j,iif) partvar(i,j,iif)=partvar(i2,j,iif) partvar(i2,j,iif)=dum endif end do end do if ((freq.gt.freq_sec).and.(iflag.eq.0)) then iflag=1 inonsec=iif endif if (iflag.eq.0) then write (13+j-1,876) freq,partvar(niter-niter/100,j,iif0), $ partvar(niter-niter/20,j,iif0), $ partvar(niter-niter/10,j,iif0), $ partvar(niter-niter/2,j,iif0) else if (j.eq.1) iavefreq=iavefreq+1 sum99(j) = sum99(j)+partvar(niter-niter/100,j,iif) sum95(j) = sum95(j)+partvar(niter-niter/20,j,iif) sum90(j) = sum90(j)+partvar(niter-niter/10,j,iif) sum50(j) = sum50(j)+partvar(niter-niter/2,j,iif) endif end do end do do j=1,nwin-1 sum99(j)=sum99(j)/float(iavefreq) sum95(j)=sum95(j)/float(iavefreq) sum90(j)=sum90(j)/float(iavefreq) sum50(j)=sum50(j)/float(iavefreq) end do c c average over non-secular band c do iif=inonsec,nfmax freq = (iif-1)*ddf+f1 do j=1,nwin-1 write (13+j-1,876) freq,sum99(j),sum95(j),sum90(j), $ sum50(j) end do end do 876 format (5f8.4) c write (6,*) 'done.' c c done -- close files c do i=1,nwin-1 close (unit=13+i-1) end do 888 format (f9.4,i5,2f9.4,f12.4,i6) 999 format (f9.4,2f9.4,f12.4,i6) 6666 format (f9.4,i5,f9.4,f12.4,i6) 7777 format (f9.4,f9.4,f12.4,i6) stop end c c c subroutine taper(n,nwin,el) c c to generate slepian tapers c ta is a real*4 array c c j. park c real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb real*8 dfac,drat,gamma,bh,ell common/nnaamm/iwflag common/npiprol/anpi common/tapsum/tapsum(20),ell(20) common/work/ip(8194) common/taperzz/z(65536) ! we use this common block for scratch space common/stap2/ta(8192,16) dimension a(8192,8),el(10) data iwflag/0/,pi/3.14159265358979d0/ equivalence (a(1,1),ta(1,1)) an=dfloat(n) ww=dble(anpi)/an cs=dcos(2.d0*pi*ww) c initialize matrix for eispack subroutine c type *,'ww,cs,an',ww,cs,an do i=0,n-1 ai=dfloat(i) a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2 a(i+1,2)=-ai*(an-ai)/2.d0 a(i+1,3)=a(i+1,2)**2 ! required by eispack routine end do eps=1.e-13 m11=1 call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip, x ierr,a(1,4),a(1,5)) c type *,ierr,rlb,rlu c call tnou('eigenvalues for the eigentapers') c type *,(el(i),i=1,nwin) call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr, x a(1,4),a(1,5),a(1,6),a(1,7),a(1,8)) c we calculate the eigenvalues of the dirichlet-kernel problem c i.e. the bandwidth retention factors c from slepian 1978 asymptotic formula, gotten from thomson 1982 eq 2.5 c supplemented by the asymptotic formula for k near 2n from slepian 1978 eq 61 dfac=an*pi*ww drat=8.d0*dfac dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac) do k=1,nwin el(k)=1.d0-dfac dfac=dfac*drat/k ! is this correct formula? yes,but fails as k -> 2n end do c call tnou('eigenvalues for the eigentapers (small k)') c type *,(el(i),i=1,nwin) gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0 do k=1,nwin bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma ell(k)=1.d0/(1.d0+dexp(pi*bh)) end do c call tnou('eigenvalues for the eigentapers (k near 2n)') c type *,(ell(i),i=1,nwin) do i=1,nwin el(i)=dmax1(ell(i),el(i)) end do c call tnou('composite asymptotics for taper eigenvalues') c type *,(el(i),i=1,nwin) c normalize the eigentapers to preserve power for a white process c i.e. they have rms value unity do k=1,nwin kk=(k-1)*n tapsum(k)=0. tapsq=0. do i=1,n aa=z(kk+i) ta(i,k)=aa tapsum(k)=tapsum(k)+aa tapsq=tapsq+aa*aa end do aa=sqrt(tapsq/n) tapsum(k)=tapsum(k)/aa do i=1,n ta(i,k)=ta(i,k)/aa end do end do c type *,'tapsum',(tapsum(i),i=1,nwin) c refft preserves amplitudes with zeropadding c zum beispiel: a(i)=1.,i=1,100 will transform at zero freq to b(f=0.)=100 c no matter how much zero padding is done c therefore we need not doctor the taper normalization, c but wait until the fft to force the desired units iwflag=1 return end c c c SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V) C C Singular value decomposition of an M by N complex matrix A, C where M .GT. N . The singular values are stored in the vector C S. The first NU columns of the M by M unitary matrix U and the C first NV columns of the N by N unitary matrix V that minimize C Det(A-USV*) are also computed. C C C P.A. Businger and G.H. Golub, "Singular Value Decomposition C of a Complex Matrix," Communications of the ACM, vol. 12, C pp. 564-565, October 1969. C C This algorithm is reprinted by permission, Association for C Computing Machinery; copyright 1969. C parameter (nwork=10000) COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R REAL *8 S(NMAX),B(nwork),C(nwork),T(nwork),zero,one,ETA,TOL, $ EPS ETA = 1.2d-7 TOL = 2.4d-32 zero = 0.d0 one = 1.d0 NP=N+IP N1=N+1 C Householder reduction C(1)=zero K=1 10 K1=K+1 C Elimination of A(I,K) , I=K+1,...,M Z=zero DO 20 I=K,M 20 Z=Z+REAL(A(I,K))**2+dIMAG(A(I,K))**2 B(K)=zero IF (Z .LE. TOL) GO TO 70 Z=SQRT(Z) B(K)=Z W=CdABS(A(K,K)) Q=dcmplx(one,zero) IF (W .NE. zero) Q=A(K,K)/W A(K,K)=Q*(Z+W) IF (K .EQ. NP) GO TO 70 DO 50 J=K1,NP Q=dcmplx(zero,zero) DO 30 I=K,M 30 Q=Q+dCONJG(A(I,K))*A(I,J) Q=Q/(Z*(Z+W)) DO 40 I=K,M 40 A(I,J)=A(I,J)-Q*A(I,K) 50 CONTINUE C Phase transformation Q=-dCONJG(A(K,K))/CdABS(A(K,K)) DO 60 J=K1,NP 60 A(K,J)=Q*A(K,J) C Elimination of A(K,J) , J=K+2,...,N 70 IF (K .EQ. N) GO TO 140 Z=zero DO 80 J=K1,N 80 Z=Z+REAL(A(K,J))**2+dIMAG(A(K,J))**2 C(K1)=zero IF (Z .LE. TOL) GO TO 130 Z=SQRT(Z) C(K1)=Z W=CdABS(A(K,K1)) Q=dcmplx(one,zero) IF (W .NE. zero) Q=A(K,K1)/W A(K,K1)=Q*(Z+W) DO 110 I=K1,M Q=dcmplx(zero,zero) DO 90 J=K1,N 90 Q=Q+dCONJG(A(K,J))*A(I,J) Q=Q/(Z*(Z+W)) DO 100 J=K1,N 100 A(I,J)=A(I,J)-Q*A(K,J) 110 CONTINUE C Phase transformation Q=-dCONJG(A(K,K1))/CdABS(A(K,K1)) DO 120 I=K1,M 120 A(I,K1)=A(I,K1)*Q 130 K=K1 GO TO 10 C Tolerance for negligible elements 140 EPS=zero DO 150 K=1,N S(K)=B(K) T(K)=C(K) 150 EPS=DMAX1(EPS,S(K)+T(K)) EPS=EPS*ETA C Initialization of U and V IF (NU .EQ. 0) GO TO 180 DO 170 J=1,NU DO 160 I=1,M 160 U(I,J)=dcmplx(zero,zero) 170 U(J,J)=dcmplx(one,zero) 180 IF (NV .EQ. 0) GO TO 210 DO 200 J=1,NV DO 190 I=1,N 190 V(I,J)=dcmplx(zero,zero) 200 V(J,J)=dcmplx(one,zero) C QR diagonalization 210 DO 380 KK=1,N K=N1-KK C Test for split 220 DO 230 LL=1,K L=K+1-LL IF (ABS(T(L)) .LE. EPS) GO TO 290 IF (ABS(S(L-1)) .LE. EPS) GO TO 240 230 CONTINUE C Cancellation of B(L) 240 CS=zero SN=one L1=L-1 DO 280 I=L,K F=SN*T(I) T(I)=CS*T(I) IF (ABS(F) .LE. EPS) GO TO 290 H=S(I) W=SQRT(F*F+H*H) S(I)=W CS=H/W SN=-F/W IF (NU .EQ. 0) GO TO 260 DO 250 J=1,N X=REAL(U(J,L1)) Y=REAL(U(J,I)) U(J,L1)=dCMPLX(X*CS+Y*SN,0.) 250 U(J,I)=dCMPLX(Y*CS-X*SN,0.) 260 IF (NP .EQ. N) GO TO 280 DO 270 J=N1,NP Q=A(L1,J) R=A(I,J) A(L1,J)=Q*CS+R*SN 270 A(I,J)=R*CS-Q*SN 280 CONTINUE C Test for convergence 290 W=S(K) IF (L .EQ. K) GO TO 360 C Origin shift X=S(L) Y=S(K-1) G=T(K-1) H=T(K) F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y) G=SQRT(F*F+one) IF (F .LT. zero) G=-G F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X C QR step CS=one SN=one L1=L+1 DO 350 I=L1,K G=T(I) Y=S(I) H=SN*G G=CS*G W=SQRT(H*H+F*F) T(I-1)=W CS=F/W SN=H/W F=X*CS+G*SN G=G*CS-X*SN H=Y*SN Y=Y*CS IF (NV .EQ. 0) GO TO 310 DO 300 J=1,N X=REAL(V(J,I-1)) W=REAL(V(J,I)) V(J,I-1)=dCMPLX(X*CS+W*SN,0.) 300 V(J,I)=dCMPLX(W*CS-X*SN,0.) 310 W=SQRT(H*H+F*F) S(I-1)=W CS=F/W SN=H/W F=CS*G+SN*Y X=CS*Y-SN*G IF (NU .EQ. 0) GO TO 330 DO 320 J=1,N Y=REAL(U(J,I-1)) W=REAL(U(J,I)) U(J,I-1)=dCMPLX(Y*CS+W*SN,0.) 320 U(J,I)=dCMPLX(W*CS-Y*SN,0.) 330 IF (N .EQ. NP) GO TO 350 DO 340 J=N1,NP Q=A(I-1,J) R=A(I,J) A(I-1,J)=Q*CS+R*SN 340 A(I,J)=R*CS-Q*SN 350 CONTINUE T(L)=zero T(K)=F S(K)=X GO TO 220 C Convergence 360 IF (W .GE. zero) GO TO 380 S(K)=-W IF (NV .EQ. 0) GO TO 380 DO 370 J=1,N 370 V(J,K)=-V(J,K) 380 CONTINUE C Sort singular values DO 450 K=1,N G=-one J=K DO 390 I=K,N IF (S(I) .LE. G) GO TO 390 G=S(I) J=I 390 CONTINUE IF (J .EQ. K) GO TO 450 S(J)=S(K) S(K)=G IF (NV .EQ. 0) GO TO 410 DO 400 I=1,N Q=V(I,J) V(I,J)=V(I,K) 400 V(I,K)=Q 410 IF (NU .EQ. 0) GO TO 430 DO 420 I=1,N Q=U(I,J) U(I,J)=U(I,K) 420 U(I,K)=Q 430 IF (N .EQ. NP) GO TO 450 DO 440 I=N1,NP Q=A(J,I) A(J,I)=A(K,I) 440 A(K,I)=Q 450 CONTINUE C Back transformation IF (NU .EQ. 0) GO TO 510 DO 500 KK=1,N K=N1-KK IF (B(K) .EQ. zero) GO TO 500 Q=-A(K,K)/CdABS(A(K,K)) DO 460 J=1,NU 460 U(K,J)=Q*U(K,J) DO 490 J=1,NU Q=dcmplx(zero,zero) DO 470 I=K,M 470 Q=Q+dCONJG(A(I,K))*U(I,J) Q=Q/(CdABS(A(K,K))*B(K)) DO 480 I=K,M 480 U(I,J)=U(I,J)-Q*A(I,K) 490 CONTINUE 500 CONTINUE 510 IF (NV .EQ. 0) GO TO 570 IF (N .LT. 2) GO TO 570 DO 560 KK=2,N K=N1-KK K1=K+1 IF (C(K1) .EQ. zero) GO TO 560 Q=-dCONJG(A(K,K1))/CdABS(A(K,K1)) DO 520 J=1,NV 520 V(K1,J)=Q*V(K1,J) DO 550 J=1,NV Q=dcmplx(zero,zero) DO 530 I=K1,N 530 Q=Q+A(K,I)*V(I,J) Q=Q/(CdABS(A(K,K1))*C(K1)) DO 540 I=K1,N 540 V(I,J)=V(I,J)-Q*dCONJG(A(K,I)) 550 CONTINUE 560 CONTINUE 570 RETURN END