      PROGRAM cro_ave

C  PROGRAMME POUR CALCULER LE SPECTRE D'UNE SERIE DE DONNEES
C  BRUITEE ET DETERMINER SI ELLE CONTIENT DES OSCILLATIONS.
C  LA METHODE CONSISTE:
C 1.  A CALCULER LE PERIODIGRAMME DE LA
C     SERIE EN UTILISANT UN PROGRAMME DE TRANSFOURMEE DE
C     FOURIER: RS
C 2.  LISSER LE PERIODIGRAMME POUR AUGMENTER LE RAPPORT
C     SIGNAL/BRUIT DANS LA TRANSFORMMEE DE FOURIER:
C     SMOOTHER
C 3.  FAIRE UN TEST DE MONTE-CARLO 
C     POUR EVALUER LA SIGNIFICATIVITE DES PICS DU SPECTRES
C     DE LA SERIE REELLE.


C1D   PARAMETER(NDATA=1512,NDEB=0,NWIN=1512,ntest=500)
      PARAMETER(NDATA=4383,NDEB=0,NWIN=4352,ntest=1000)
C     PARAMETER(NSMO=3,NTAP=1)
      PARAMETER(NSMO=18,NTAP=3)
c     PARAMETER(NSMO=18,NTAP=NSMO)
      PARAMETER(NJUMP=100)
C SERIES WITH AND WITHOUT ANNUAL CYCLE
      real time(nwin),force(nwin),forra(nwin),ranfor(nwin)
C RANDOM PROCESSES

      REAL FORCR(NWIN),FORCI(NWIN)
      INTEGER NNN(1),NDIM,IFRWD,ICPLX
      REAL WORKR(NWIN),WORKI(NWIN)
C 
      real FREQ(NWIN),CROSS(NWIN),CROST(NWIN),CROSR(NWIN,NTEST)
      real COPHA(NWIN),COHES(NWIN),COHET(NWIN),COHER(NWIN,NTEST)
      real SPECF(NWIN),SPECFR(NWIN,NTEST)
      real SPECFM(-NWIN:NWIN),SPECF95(-NWIN:NWIN),SPECF99(-NWIN:NWIN)
      real SPECP(NWIN),SPECPR(NWIN,NTEST)
      real SPECPM(-NWIN:NWIN),SPECP95(-NWIN:NWIN),SPECP99(-NWIN:NWIN)
      real TEMPF(NWIN),TEMPP(NWIN),SMOOT(-NSMO:NWIN+NSMO)
      real TEMPCR(NWIN),TEMPCI(NWIN)
      complex CI                                                
      real CROS99(-NWIN:NWIN),CROSM(-NWIN:NWIN),CROS95(-NWIN:NWIN)
      real COHE99(-NWIN:NWIN),COHEM(-NWIN:NWIN),COHE95(-NWIN:NWIN)
c
c  INPUT FILES        
c
C IN  FIRST: CAUSE

C   FORCING TERM
      open(11,file='../out_nhe/z700_atl1_3d.noa',form='formatted')
      open(12,file='../out_nhe/z700_atl1_3d.red',form='formatted')

c
c  OUTPUT FILES       
c
      open(17,file='../out_nhe/z700_atl1_3d.spe',form='formatted')

         NNN(1)=NWIN
         NDIM=1
         PI=ACOS(-1.)
         CI=(0.,1.)
         DT=3.       

C  LECTURE DES DONNEES, MOYENNE SUR 3 JOURS
      
       do 10 nd=1,nwin
       read(11,*)time(nd),force(nd)
 10    continue

C NORMALISATION AVEC STATISTIQUES POUR LE FENETRE D'INTERET

       formea=0.

       do 100 nd=1,ndata
       formea=formea+force(nd)/float(nwin)             
       forvar=forvar+force(nd)**2/float(nwin)             
 100   continue

       forvar=forvar-formea**2

       do 11 nd=1,ndata
          force(nd)=(force(nd)-formea)
 11    continue

       tt0=0.
       tt1=0.
       do 110 nd=ndeb+2,ndeb+nwin
       tt0=tt0+force(nd  )*force(nd)/float(nwin-1)
       tt1=tt1+force(nd-1)*force(nd)/float(nwin-1)
 110   continue

       print *,'tt0 tt1:',tt0,tt1

C TESTS AGAINST UNCORRELATED SERIES:

       za=tt1/tt0

       var1=0.
       do 111 nd=ndeb+1,nwin+ndeb
          temp1=force(nd)-za*force(nd-1)
          var1=var1+temp1*temp1/float(nwin)
          temp2=ALEAS(0.)-0.5
          var2=var2+temp2*temp2/float(nwin)
 111   continue

          alpha=sqrt(var1/var2)

          print *,'alpha:',alpha

C  FROM TIME TO TIME, I LIKE TO MAKE THE SPECTRA OF A RANDOM VARIABLE1

       nd=1
          forra(nd)=force(nd)
       write(12,*)time(nd),forra(nd)
          do nd=2,nwin
            forra(nd)=za*forra(nd-1)+alpha*(ALEAS(0.)-0.5)
       write(12,*)time(nd),forra(nd)
          enddo

C   SPECTRE DE LA SERIE

       do 20 nw=1,nwin
          TEMPF(NW)=0.
          TEMPP(NW)=0.
 20    continue

         DO 31 nw=1,nwin
            FORCR(NW)=force(ndeb+nw)
            FORCI(NW)=0.
 31      CONTINUE

         call dtrend(forcr,ax,bx,nwin)

         CALL FOURT(FORCR,FORCI,NNN,NDIM,1,0,WORKR,WORKI)

         DO 32 nw=1,nwin/2
            TEMPF(NW)=FORCR(NW)**2+FORCI(NW)**2
            FREQ(NW)=(NW-1)/DT/FLOAT(NWIN)
 32      CONTINUE

            CALL SMOOTHER(NWIN,NSMO,NTAP,TEMPF,SMOOT)

         DO 35 nw=1,nwin/2
         SPECF(NW)=TEMPF(NW)
 35      CONTINUE

C SPECTRES DES SERIES ALEATOIRES

       FORCR(NWIN)=forra(100)
 
       do 4 NT=1,NTEST

       if(mod(nt,25).eq.0)print *,' spectra AR:',nt

       do 40 nw=1,nwin
          TEMPF(NW)=0.
 40    continue


            FORCR(1)=alpha*(ALEAS(0.)-0.5)
         DO 42 nw=2,nwin
            FORCR(NW)=za*FORCR(nw-1)+alpha*(ALEAS(0.)-0.5)
            FORCI(NW)=0.
 42      CONTINUE

       call dtrend(forcr,ax,bx,nwin)

         CALL FOURT(FORCR,FORCI,NNN,NDIM,1,0,WORKR,WORKI)

       DO 43 nw=1,nwin/2
          TEMPF(NW)=FORCR(NW)**2+FORCI(NW)**2
          FREQ(NW)=(NW-1)/DT/FLOAT(NWIN)
 43    CONTINUE

         CALL SMOOTHER(NWIN,NSMO,NTAP,TEMPF,SMOOT)

       DO 45 nw=1,nwin/2
          SPECFR(NW,NT)=TEMPF(NW)
 45    CONTINUE
       
 4     CONTINUE
       


C  SOME STATISTICS: MEAN
        
        DO 6 nw=1,nwin/2

        SPECFM(NW)=0.
        SPECF99(NW)=0.
        SPECF95(NW)=0.


        DO NT=1,NTEST
           SPECFM(NW)=SPECFM(NW)+SPECFR(NW,NT)/FLOAT(NTEST)
           SPECF99(NW)=SPECFM(NW)
           SPECF95(NW)=SPECFM(NW)
        ENDDO

 6      CONTINUE
      

C  TEST 99%

        DO 62 nw=njump,nwin/2,njump

        print *,' TESTS:',nw,nwin/2

        DO NT=1,NTEST
          SPECF99(NW)=AMAX1(SPECF99(NW),SPECFR(NW,NT))
        ENDDO

        IF(NTEST.GE.200)THEN

        DO MC=1,NTEST/100-1

           DO NT=1,NTEST
           IF(SPECFR(NW,NT).EQ.SPECF99(NW))THEN
           SPECFR(NW,NT)=0.
           SPECF99(NW)=0.
           ENDIF
           ENDDO
           DO NT=1,NTEST
           SPECF99(NW)=AMAX1(SPECF99(NW),SPECFR(NW,NT))
           ENDDO

        ENDDO

           SPECF95(NW)=SPECF99(NW)
         DO MC=NTEST/100+1,NTEST/20

           DO NT=1,NTEST
           IF(SPECFR(NW,NT).EQ.SPECF95(NW))THEN
           SPECFR(NW,NT)=0.
           SPECF95(NW)=0.
           ENDIF
           ENDDO
           DO NT=1,NTEST
           SPECF95(NW)=AMAX1(SPECF95(NW),SPECFR(NW,NT))
           ENDDO

        ENDDO


        ENDIF

62      CONTINUE      

        PRINT *,' END OF 99'


C  WRITE OUTPUT FILES.

C  WRITE SPECTRA OF SERIES 1 (FILE 17)


       DO nw=2,nwin/2
          WRITE(17,*) FREQ(NW),SPECF(NW)        
       ENDDO
       WRITE(17,*)' '
       DO nw=njump,nwin/2
         WRITE(17,*) FREQ(NW),SPECFM(NW)        
       ENDDO
        WRITE(17,*)' '
       DO nw=njump,nwin/2,njump
          WRITE(17,*) FREQ(NW),SPECF95(NW)        
       ENDDO
       WRITE(17,*)' '
       DO nw=njump,nwin/2,njump
        WRITE(17,*) FREQ(NW),SPECF99(NW)        
       ENDDO


        IF(NSMO.EQ.0)GOTO 710
        WRITE(17,*)' '
           WRITE(17,*)FREQ(NWin/4-1),2.*SPECFM(nwin/6)
        DO nw=nwin/4,nwin/4+2*nsmo
           WRITE(17,*) FREQ(NW),SMOOT(NW)*SPECFM(nwin/6)
        ENDDO
           WRITE(17,*)FREQ(nwin/4+2*nsmo+1),2.*SPECFM(nwin/6)
           WRITE(17,*)FREQ(NWin/4-1),2.*SPECFM(nwin/6)
 710   CONTINUE

      STOP
      END

       subroutine dtrend(x,ax,bx,n)
       real xxm,iim,xim,xm,im
       real x(n),det,ax,bx

       xxm=0.
       iim=0.
       xim=0.
       xm=0.
       im=0.
       do i=1,n
         xxm=xxm+x(i)*x(i)/float(n)
         iim=iim+float(i*i)/float(n)
         xim=xim+x(i)*float(i)/float(n)
         xm=xm+x(i)/float(n)
         im=im+float(i)/float(n)
       enddo

       det=iim-im**2
       ax=(xim-xm*im)/det
       bx=(iim*xm-xim*im)/det

       do i=1,n
          x(i)=x(i)-ax*float(i)-bx
       enddo

       return
       end

       subroutine taper(x,n)
       real pi

       real x(n)            
       pi=acos(-1.)
       do i=1,n
         x(i)=x(i)*0.25*(1.-cos(2.*pi*float(i-1)/float(n-1)))**2
       enddo

       return
       end

       subroutine smoother(nwin,nsmo,ntap,temp,smoo)
       integer nwin,nsmo,ntap
       real temp(nwin),smoo(-nsmo:nwin+nsmo)
       real weight(-100:100)

       xpi=acos(-1.)
       sum=0.

       if(nsmo.eq.0) return
       if(ntap.lt.nsmo)then
       do j=-nsmo,-nsmo+ntap
       weight(j)=0.5*(1.-cos(xpi*float(j+nsmo+1)/float(ntap+2)))
       sum=sum+weight(j)
       enddo
       do j=-nsmo+ntap+1,0
       weight(j)=1.
       sum=sum+weight(j)
       enddo
       do j=1,nsmo
       weight(j)=weight(-j)
       sum=sum+weight(j)
       enddo
       else
       do j=-nsmo,nsmo
       weight(j)=1.
       sum=sum+weight(j)
       enddo
       endif

       do j=-nsmo,nsmo
       weight(j)=weight(j)/sum
       enddo

       do i=-nsmo,0
       smoo(i)=0.
       enddo
       do i=nwin+1,nwin+nsmo
       smoo(i)=0.
       enddo
       do i=1,nwin
          smoo(i)=temp(i)
       enddo
       do i=1,nwin
          temp(i)=0.
          do j=-nsmo,nsmo
             temp(i)=temp(i)+smoo(i+j)*weight(j)
          enddo
       enddo
       do j=-nsmo,nsmo
       smoo(nwin/4+nsmo+j)=2.+weight(j)*sum
       enddo

       return
       end

       subroutine smoother_old(nwin,nsmo,temp,smoo)
       integer nwin,nsmo
       real temp(nwin),smoo(-nsmo:nwin+nsmo)
       real weight(-100:100)
 
       xpi=acos(-1.)

       sum=3.
       weight(-nsmo)=0.25
       weight( nsmo)=0.25
       weight(-nsmo+1)=0.5
       weight( nsmo-1)=0.5
       weight(-nsmo+2)=0.75
       weight( nsmo-2)=0.75
       do j=-nsmo+3,nsmo-3
c      weight(j)=0.5*(1.+cos(xpi*float(j)/float(nsmo+1)))
       weight(j)=1.
       sum=sum+weight(j)
       enddo
       do j=-nsmo,nsmo
       weight(j)=weight(j)/sum
       enddo
       do i=-nsmo,0
       smoo(i)=0.
       enddo
       do i=nwin+1,nwin+nsmo
       smoo(i)=0.
       enddo
       do i=1,nwin
          smoo(i)=temp(i)
       enddo
       do i=1,nwin
          temp(i)=0.
          do j=-nsmo,nsmo
             temp(i)=temp(i)+smoo(i+j)*weight(j)
          enddo
       enddo
       do j=-nsmo,nsmo
       smoo(nwin/4+nsmo+j)=2.+weight(j)*sum
       enddo

       return
       end



      SUBROUTINE FOURT(DATAR,DATAI,NN,NDIM,
     +                       IFRWD,ICPLX,WORKR,WORKI)
C
      DIMENSION DATAR(*),DATAI(*),NN(*),WORKR(*),WORKI(*),IFACT(20)
C
C     THE COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN
C     EVALUATES COMPLEX FOURIER SERIES FOR COMPLEX OR REAL FUNCTIONS.
C     THAT IS, IT COMPUTES
C     FTRAN(J1,J2,...)=SUM(DATA(I1,I2,...)*W1**(I1-1)*(J1-1)
C                                         *W2**(I2-1)*(J2-1)*...),
C     WHERE W1=EXP(-2*PI*SQRT(-1)/NN(1)), W2=EXP(-2*PI*SQRT(-1)/NN(2)),
C     ETC. AND I1 AND J1 RUN FROM 1 TO NN(1), I2 AND J2 RUN FROM 1 TO
C     NN(2), ETC.  THERE IS NO LIMIT ON THE DIMENSIONALITY (NUMBER OF
C     SUBSCRIPTS) OF THE ARRAY OF DATA.  THE PROGRAM WILL PERFORM
C     A THREE-DIMENSIONAL FOURIER TRANSFORM AS EASILY AS A ONE-DIMEN-
C     SIONAL ONE, THO IN A PROPORTIONATELY GREATER TIME.  AN INVERSE
C     TRANSFORM CAN BE PERFORMED, IN WHICH THE SIGN IN THE EXPONENTIALS
C     IS +, INSTEAD OF -.  IF AN INVERSE TRANSFORM IS PERFORMED UPON
C     AN ARRAY OF TRANSFORMED DATA, THE ORIGINAL DATA WILL REAPPEAR,
C     MULTIPLIED BY NN(1)*NN(2)*...  THE ARRAY OF INPUT DATA MAY BE
C     REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF
C     ABOUT THIRTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX.
C     (FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.)
C     THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE
C     ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA.  THE LENGTH
C     OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER.  THE
C     PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS
C     PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO.
C     TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA.  LET NTOT BE THE
C     TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT
C     IS, NTOT=NN(1)*NN(2)*...  DECOMPOSE NTOT INTO ITS PRIME FACTORS,
C     SUCH AS 2**K2 * 3**K3 * 5**K5 * ...  LET SUM2 BE THE SUM OF ALL
C     THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2.  LET SUMF BE
C     THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+..
C     THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA
C     IS T = T0 + T1*NTOT + T2*NTOT*SUM2 + T3*NTOT*SUMF.  FOR THE PAR-
C     TICULAR IMPLEMENTATION FORTRAN 32 ON THE CDC 3300 (FLOATING POINT
C     ADD TIME = SIX MICROSECONDS),
C     T = 3000 + 600*NTOT + 50*NTOT*SUM2 + 175*NTOT*SUMF MICROSECONDS
C     ON COMPLEX DATA.
C
C     IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME
C     PROPORTIONAL TO NTOT**2.  FOR HIGHLY COMPOSITE NTOT, THE SAVINGS
C     OFFERED BY COOLEY-TUKEY CAN BE DRAMATIC.  A MATRIX 100 BY 100 WILL
C     BE TRANSFORMED IN TIME PROPORTIONAL TO 10000*(2+2+2+2+5+5+5+5) =
C     280,000 (ASSUMING T2 AND T3 TO BE ROUGHLY COMPARABLE) VERSUS
C     10000**2 = 100,000,000 FOR THE STRAIGHTFORWARD TECHNIQUE.
C
C     THE COOLEY-TUKEY ALGORITHM PLACES TWO RESTRICTIONS UPON THE
C     NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT
C     THE DATA FROM ONE CYCLE OF A PERIODIC FUNCTION.  THEY ARE--
C     1.  THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES
C     MUST BE THE SAME.
C     2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN,
C     THEY MUST BE EQUI-SPACED AT INTERVALS OF DT.  FURTHER, THE TRANS-
C     FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI-
C     SPACED FROM 0 TO 2*PI*(NN(I)-1)/(NN(I)*DT) AT INTERVALS OF
C     2*PI/(NN(I)*DT) FOR EACH DIMENSION OF LENGTH NN(I).  OF COURSE,
C     DT NEED NOT BE THE SAME FOR EVERY DIMENSION.
C
C     THE CALLING SEQUENCE IS--
C     CALL FOURT(DATAR,DATAI,NN,NDIM,IFRWD,ICPLX,WORKR,WORKI)
C
C     DATAR AND DATAI ARE THE ARRAYS USED TO HOLD THE REAL AND IMAGINARY
C     PARTS OF THE INPUT DATA ON INPUT AND THE TRANSFORM VALUES ON
C     OUTPUT.  THEY ARE FLOATING POINT ARRAYS, MULTIDIMENSIONAL WITH
C     IDENTICAL DIMENSIONALITY AND EXTENT.  THE EXTENT OF EACH DIMENSION
C     IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM.  THAT IS,
C     NDIM IS THE DIMENSIONALITY OF THE ARRAYS DATAR AND DATAI.
C     IFRWD IS AN INTEGER USED TO INDICATE THE DIRECTION OF THE FOURIER
C     TRANSFORM.  IT IS NON-ZERO TO INDICATE A FORWARD TRANSFORM
C     (EXPONENTIAL SIGN IS -) AND ZERO TO INDICATE AN INVERSE TRANSFORM
C     (SIGN IS +).  ICPLX IS AN INTEGER TO INDICATE WHETHER THE DATA
C     ARE REAL OR COMPLEX.  IT IS NON-ZERO FOR COMPLEX, ZERO FOR REAL.
C     IF IT IS ZERO (REAL) THE CONTENTS OF ARRAY DATAI WILL BE ASSUMED
C     TO BE ZERO, AND NEED NOT BE EXPLICITLY SET TO ZERO.  AS EXPLAINED
C     ABOVE, THE TRANSFORM RESULTS ARE ALWAYS COMPLEX AND ARE STORED
C     IN DATAR AND DATAI ON RETURN.  WORKR AND WORKI ARE ARRAYS USED
C     FOR WORKING STORAGE.  THEY ARE NOT NECESSARY IF ALL THE DIMENSIONS
C     OF THE DATA ARE POWERS OF TWO.  IN THIS CASE, THE ARRAYS MAY BE
C     REPLACED BY THE NUMBER 0 IN THE CALLING SEQUENCE.  THUS, USE OF
C     POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE.  IF ANY DIMENSION
C     IS NOT A POWER OF TWO, THESE ARRAYS MUST BE SUPPLIED.  THEY ARE
C     FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO THE LARGEST
C     ARRAY DIMENSION, THAT IS, TO THE LARGEST VALUE OF NN(I).
C     WORKR AND WORKI, IF SUPPLIED, MUST NOT BE THE SAME ARRAYS AS DATAR
C     OR DATAI.  ALL SUBSCRIPTS OF ALL ARRAYS BEGIN AT 1.
C
C     EXAMPLE 1.  THREE-DIMENSIONAL FORWARD FOURIER TRANSFORM OF A
C     COMPLEX ARRAY DIMENSIONED 100 BY 16 BY 13.
C     DIMENSION DATAR(100,16,13),DATAI(100,16,13),WORKR(100),WORKI(100)
C     DIMENSION NN(3)
C     NN(1)=100
C     NN(2)=16
C     NN(3)=13
C     CALL FOURT(DATAR,DATAI,NN,3,1,1,WORKR,WORKI)
C
C     EXAMPLE 2.  ONE-DIMENSIONAL FORWARD TRANSFORM OF A REAL ARRAY OF
C     LENGTH 64.
C     DIMENSION DATAR(64),DATAI(64)
C     CALL FOURT(DATAR,DATAI,64,1,1,0,0,0)
C
C     THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM.  THE
C     PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE.
C
C     THE SINE AND COSINE VALUES REQUIRED FOR THE TRANSFORM ARE
C     GENERATED RECURSIVELY.  IF DOUBLE PRECISION IS AVAILABLE, IT IS
C     STRONGLY URGED THAT THE FOLLOWING VARIABLES BE SO DECLARED TO
C     REDUCE ACCUMULATION OF ROUNDOFF ERROR--
      DOUBLE PRECISION TWOPI,THETA,WSTPR,WSTPI,WMINR,WMINI,WR,WI,WTEMP
     *  ,THETM,WMSTR,WMSTI,TWOWR,SR,SI,OLDSR,OLDSI,STMPR,STMPI
C     IN ADDITION, TWOPI SHOULD BE ASSIGNED A SUFFICIENTLY PRECISE
C     VALUE AND THE VARIOUS CALLS TO THE FUNCTIONS COS AND SIN
C     SHOULD BE CHANGED TO DCOS AND DSIN.
C
C     PROGRAM BY NORMAN BRENNER FROM THE BASIC ALGORITHM BY CHARLES
C     RADER (BOTH OF MIT LINCOLN LABORATORY).  MAY 1967.  THE IDEA
C     FOR THE BIT REVERSAL WAS SUGGESTED BY RALPH ALTER (ALSO MIT LL).
C     ADAPTED FROM THE WORK OF JAMES W. COOLEY AND JOHN W. TUKEY,
C     AN ALGORITHM FOR THE MACHINE CALCULATION OF COMPLEX FOURIER
C     SERIES, MATH. COMPUT. 19, 90 (APRIL 1965), 297-301.
C
C     FOURT MODIFIED FROM A SUBROUTINE TO AN INTEGER FUNCTION BY
C     C. POTTER (NCSA). JAN 1988.
C
      IF(NDIM-1)920,1,1
1     NTOT=1
      DO 2 IDIM=1,NDIM
2     NTOT=NTOT*NN(IDIM)
      TWOPI=6.283185307179586476925286766559D0
C
C     MAIN LOOP FOR EACH DIMENSION
      NP1=1
      DO 910 IDIM=1,NDIM
      N=NN(IDIM)
      NP2=NP1*N
      IF(N-1)920,900,5
C
C     IS N A POWER OF TWO AND IF NOT, WHAT ARE ITS FACTORS
5     M=N
      NTWO=NP1
      IF=1
      IDIV=2
10    IQUOT=M/IDIV
      IREM=M-IDIV*IQUOT
      IF(IQUOT-IDIV)50,11,11
11    IF(IREM)20,12,20
12    NTWO=NTWO+NTWO
      IFACT(IF)=IDIV
      IF=IF+1
      M=IQUOT
      GO TO 10
20    IDIV=3
      INON2=IF
30    IQUOT=M/IDIV
      IREM=M-IDIV*IQUOT
      IF(IQUOT-IDIV)60,31,31
31    IF(IREM)40,32,40
32    IFACT(IF)=IDIV
      IF=IF+1
      M=IQUOT
      GO TO 30
40    IDIV=IDIV+2
      GO TO 30
50    INON2=IF
      IF(IREM)60,51,60
51    NTWO=NTWO+NTWO
      GO TO 70
60    IFACT(IF)=M
70    NON2P=NP2/NTWO
C
C     SEPARATE FOUR CASES--
C        1. COMPLEX TRANSFORM
C        2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION.  METHOD--
C           TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CON-
C           JUGATE SYMMETRY.
C        3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD.  METHOD--
C           SET THE IMAGINARY PARTS TO ZERO.
C        4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN.  METHOD--
C           TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS
C           ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS
C           ARE THE ODD-NUMBERED REAL VALUES.  UNSCRAMBLE AND SUPPLY
C           THE SECOND HALF BY CONJUGATE SYMMETRY.
      ICASE=1
      IFMIN=1
      IF(ICPLX)100,71,100
71    ICASE=2
      IF(IDIM-1)72,72,100
72    ICASE=3
      IF(NTWO-NP1)100,100,73
73    ICASE=4
      IFMIN=2
      NTWO=NTWO/2
      N=N/2
      NP2=NP2/2
      NTOT=NTOT/2
      I=1
      DO 80 J=1,NTOT
      DATAR(J)=DATAR(I)
      DATAI(J)=DATAR(I+1)
80    I=I+2
C
C     SHUFFLE DATA BY BIT REVERSAL, SINCE N=2**K.  AS THE SHUFFLING
C     CAN BE DONE BY SIMPLE INTERCHANGE, NO WORKING ARRAY IS NEEDED
100   IF(NON2P-1)101,101,200
101   NP2HF=NP2/2
      J=1
      DO 150 I2=1,NP2,NP1
      IF(J-I2)121,130,130
121   I1MAX=I2+NP1-1
      DO 125 I1=I2,I1MAX
      DO 125 I3=I1,NTOT,NP2
      J3=J+I3-I2
      TEMPR=DATAR(I3)
      TEMPI=DATAI(I3)
      DATAR(I3)=DATAR(J3)
      DATAI(I3)=DATAI(J3)
      DATAR(J3)=TEMPR
125   DATAI(J3)=TEMPI
130   M=NP2HF
140   IF(J-M)150,150,141
141   J=J-M
      M=M/2
      IF(M-NP1)150,140,140
150   J=J+M
      GO TO 300
C
C     SHUFFLE DATA BY DIGIT REVERSAL FOR GENERAL N
200   DO 270 I1=1,NP1
      DO 270 I3=I1,NTOT,NP2
      J=I3
      DO 260 I=1,N
      IF(ICASE-3)210,220,210
210   WORKR(I)=DATAR(J)
      WORKI(I)=DATAI(J)
      GO TO 240
220   WORKR(I)=DATAR(J)
      WORKI(I)=0.
240   IFP2=NP2
      IF=IFMIN
250   IFP1=IFP2/IFACT(IF)
      J=J+IFP1
      IF(J-I3-IFP2)260,255,255
255   J=J-IFP2
      IFP2=IFP1
      IF=IF+1
      IF(IFP2-NP1)260,260,250
260   CONTINUE
      I2MAX=I3+NP2-NP1
      I=1
      DO 270 I2=I3,I2MAX,NP1
      DATAR(I2)=WORKR(I)
      DATAI(I2)=WORKI(I)
270   I=I+1
C
C     SPECIAL CASE--  W=1
300   I1RNG=NP1
      GO TO(302,301,302,302),ICASE
301   I1RNG=NP0*(1+NPREV/2)
302   IF(NTWO-NP1)600,600,303
303   DO 430 I1=1,I1RNG
      IMIN=NP1+I1
      ISTEP=2*NP1
      GO TO 330
310   J=I1
      DO 320 I=IMIN,NTOT,ISTEP
      TEMPR=DATAR(I)
      TEMPI=DATAI(I)
      DATAR(I)=DATAR(J)-TEMPR
      DATAI(I)=DATAI(J)-TEMPI
      DATAR(J)=DATAR(J)+TEMPR
      DATAI(J)=DATAI(J)+TEMPI
320   J=J+ISTEP
      IMIN=IMIN+IMIN-I1
      ISTEP=ISTEP+ISTEP
330   IF(ISTEP-NTWO)310,310,331
C
C     SPECIAL CASE--  W=-SQRT(-1)
331   IMIN=3*NP1+I1
      ISTEP=4*NP1
      GO TO 420
400   J=IMIN-ISTEP/2
      DO 410 I=IMIN,NTOT,ISTEP
      IF(IFRWD)401,402,401
401   TEMPR=DATAI(I)
      TEMPI=-DATAR(I)
      GO TO 403
402   TEMPR=-DATAI(I)
      TEMPI=DATAR(I)
403   DATAR(I)=DATAR(J)-TEMPR
      DATAI(I)=DATAI(J)-TEMPI
      DATAR(J)=DATAR(J)+TEMPR
      DATAI(J)=DATAI(J)+TEMPI
410   J=J+ISTEP
      IMIN=IMIN+IMIN-I1
      ISTEP=ISTEP+ISTEP
420   IF(ISTEP-NTWO)400,400,430
430   CONTINUE
C
C     MAIN LOOP FOR FACTORS OF TWO.  W=EXP(-2*PI*SQRT(-1)*M/MMAX)
      THETA=-TWOPI/8.
      WSTPR=0.
      WSTPI=-1.
      IF(IFRWD)502,501,502
501   THETA=-THETA
      WSTPI=1.
502   MMAX=8*NP1
      GO TO 540
500   WMINR=COS(THETA)
      WMINI=SIN(THETA)
      WR=WMINR
      WI=WMINI
      MMIN=MMAX/2+NP1
      MSTEP=NP1+NP1
      DO 530 M=MMIN,MMAX,MSTEP
      DO 525 I1=1,I1RNG
      ISTEP=MMAX
      IMIN=M+I1
510   J=IMIN-ISTEP/2
      DO 520 I=IMIN,NTOT,ISTEP
      TEMPR=DATAR(I)*WR-DATAI(I)*WI
      TEMPI=DATAR(I)*WI+DATAI(I)*WR
      DATAR(I)=DATAR(J)-TEMPR
      DATAI(I)=DATAI(J)-TEMPI
      DATAR(J)=DATAR(J)+TEMPR
      DATAI(J)=DATAI(J)+TEMPI
520   J=J+ISTEP
      IMIN=IMIN+IMIN-I1
      ISTEP=ISTEP+ISTEP
      IF(ISTEP-NTWO)510,510,525
525   CONTINUE
      WTEMP=WR*WSTPI
      WR=WR*WSTPR-WI*WSTPI
530   WI=WI*WSTPR+WTEMP
      WSTPR=WMINR
      WSTPI=WMINI
      THETA=THETA/2.
      MMAX=MMAX+MMAX
540   IF(MMAX-NTWO)500,500,600
C
C     MAIN LOOP FOR FACTORS NOT EQUAL TO TWO.
C     W=EXP(-2*PI*SQRT(-1)*(J2-I3)/IFP2)
600   IF(NON2P-1)700,700,601
601   IFP1=NTWO
      IF=INON2
610   IFP2=IFACT(IF)*IFP1
      THETA=-TWOPI/FLOAT(IFACT(IF))
      IF(IFRWD)612,611,612
611   THETA=-THETA
612   THETM=THETA/FLOAT(IFP1/NP1)
      WSTPR=COS(THETA)
      WSTPI=SIN(THETA)
      WMSTR=COS(THETM)
      WMSTI=SIN(THETM)
      WMINR=1.
      WMINI=0.
      DO 660 J1=1,IFP1,NP1
      I1MAX=J1+I1RNG-1
      DO 650 I1=J1,I1MAX
      DO 650 I3=I1,NTOT,NP2
      I=1
      WR=WMINR
      WI=WMINI
      J2MAX=I3+IFP2-IFP1
      DO 640 J2=I3,J2MAX,IFP1
      TWOWR=WR+WR
      JMIN=I3
      J3MAX=J2+NP2-IFP2
      DO 630 J3=J2,J3MAX,IFP2
      J=JMIN+IFP2-IFP1
      SR=DATAR(J)
      SI=DATAI(J)
      OLDSR=0.
      OLDSI=0.
      J=J-IFP1
620   STMPR=SR
      STMPI=SI
      SR=TWOWR*SR-OLDSR+DATAR(J)
      SI=TWOWR*SI-OLDSI+DATAI(J)
      OLDSR=STMPR
      OLDSI=STMPI
      J=J-IFP1
      IF(J-JMIN)621,621,620
621   WORKR(I)=WR*SR-WI*SI-OLDSR+DATAR(J)
      WORKI(I)=WI*SR+WR*SI-OLDSI+DATAI(J)
      JMIN=JMIN+IFP2
630   I=I+1
      WTEMP=WR*WSTPI
      WR=WR*WSTPR-WI*WSTPI
640   WI=WI*WSTPR+WTEMP
      I=1
      DO 650 J2=I3,J2MAX,IFP1
      J3MAX=J2+NP2-IFP2
      DO 650 J3=J2,J3MAX,IFP2
      DATAR(J3)=WORKR(I)
      DATAI(J3)=WORKI(I)
650   I=I+1
      WTEMP=WMINR*WMSTI
      WMINR=WMINR*WMSTR-WMINI*WMSTI
660   WMINI=WMINI*WMSTR+WTEMP
      IF=IF+1
      IFP1=IFP2
      IF(IFP1-NP2)610,700,700
C
C     COMPLETE A REAL TRANSFORM IN THE 1ST DIMENSION, N EVEN, BY CON-
C     JUGATE SYMMETRIES.
700   GO TO (900,800,900,701),ICASE
701   NHALF=N
      N=N+N
      THETA=-TWOPI/FLOAT(N)
      IF(IFRWD)703,702,703
702   THETA=-THETA
703   WSTPR=COS(THETA)
      WSTPI=SIN(THETA)
      WR=WSTPR
      WI=WSTPI
      IMIN=2
      JMIN=NHALF
      GO TO 725
710   J=JMIN
      DO 720 I=IMIN,NTOT,NP2
      SUMR=(DATAR(I)+DATAR(J))/2.
      SUMI=(DATAI(I)+DATAI(J))/2.
      DIFR=(DATAR(I)-DATAR(J))/2.
      DIFI=(DATAI(I)-DATAI(J))/2.
      TEMPR=WR*SUMI+WI*DIFR
      TEMPI=WI*SUMI-WR*DIFR
      DATAR(I)=SUMR+TEMPR
      DATAI(I)=DIFI+TEMPI
      DATAR(J)=SUMR-TEMPR
      DATAI(J)=-DIFI+TEMPI
720   J=J+NP2
      IMIN=IMIN+1
      JMIN=JMIN-1
      WTEMP=WR*WSTPI
      WR=WR*WSTPR-WI*WSTPI
      WI=WI*WSTPR+WTEMP
725   IF(IMIN-JMIN)710,730,740
730   IF(IFRWD)731,740,731
731   DO 735 I=IMIN,NTOT,NP2
735   DATAI(I)=-DATAI(I)
740   NP2=NP2+NP2
      NTOT=NTOT+NTOT
      J=NTOT+1
      IMAX=NTOT/2+1
745   IMIN=IMAX-NHALF
      I=IMIN
      GO TO 755
750   DATAR(J)=DATAR(I)
      DATAI(J)=-DATAI(I)
755   I=I+1
      J=J-1
      IF(I-IMAX)750,760,760
760   DATAR(J)=DATAR(IMIN)-DATAI(IMIN)
      DATAI(J)=0.
      IF(I-J)770,780,780
765   DATAR(J)=DATAR(I)
      DATAI(J)=DATAI(I)
770   I=I-1
      J=J-1
      IF(I-IMIN)775,775,765
775   DATAR(J)=DATAR(IMIN)+DATAI(IMIN)
      DATAI(J)=0.
      IMAX=IMIN
      GO TO 745
780   DATAR(1)=DATAR(1)+DATAI(1)
      DATAI(1)=0.
      GO TO 900
C
C     COMPLETE A REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION BY
C     CONJUGATE SYMMETRIES.
800   IF(NPREV-2)900,900,805
805   DO 860 I3=1,NTOT,NP2
      I2MAX=I3+NP2-NP1
      DO 860 I2=I3,I2MAX,NP1
      IMAX=I2+NP1-1
      IMIN=I2+I1RNG
      JMAX=I3+I3+NP1-IMIN
      IF(I2-I3)820,820,810
810   JMAX=JMAX+NP2
820   IF(IDIM-2)850,850,830
830   J=JMAX+NP0
      DO 840 I=IMIN,IMAX
      DATAR(I)=DATAR(J)
      DATAI(I)=-DATAI(J)
840   J=J-1
850   J=JMAX
      DO 860 I=IMIN,IMAX,NP0
      DATAR(I)=DATAR(J)
      DATAI(I)=-DATAI(J)
860   J=J-NP0
C
C     END OF LOOP ON EACH DIMENSION
900   NP0=NP1
      NP1=NP2
910   NPREV=N
920   RETURN
      END

      FUNCTION ALEAS (R)
C***BEGIN PROLOGUE  RAND
C***PURPOSE  Generate a uniformly distributed random number.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  L6A21
C***TYPE      SINGLE PRECISION (RAND-S)
C***KEYWORDS  FNLIB, RANDOM NUMBER, SPECIAL FUNCTIONS, UNIFORM
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C      This pseudo-random number generator is portable among a wide
C variety of computers.  RAND(R) undoubtedly is not as good as many
C readily available installation dependent versions, and so this
C routine is not recommended for widespread usage.  Its redeeming
C feature is that the exact same random numbers (to within final round-
C off error) can be generated from machine to machine.  Thus, programs
C that make use of random numbers can be easily transported to and
C checked in a new environment.
C
C      The random numbers are generated by the linear congruential
C method described, e.g., by Knuth in Seminumerical Methods (p.9),
C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
C sequence, the I+1 -st number is generated from
C             X(I+1) = (A*X(I) + C) MOD M,
C where here M = 2**22 = 4194304, C = 1731 and several suitable values
C of the multiplier A are discussed below.  Both the multiplier A and
C random number X are represented in double precision as two 11-bit
C words.  The constants are chosen so that the period is the maximum
C possible, 4194304.
C
C      In order that the same numbers be generated from machine to
C machine, it is necessary that 23-bit integers be reducible modulo
C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
C integers be multiplied exactly.  Furthermore, if the restart option
C is used (where R is between 0 and 1), then the product R*2**22 =
C R*4194304 must be correct to the nearest integer.
C
C      The first four random numbers should be .0004127026,
C .6750836372, .1614754200, and .9086198807.  The tenth random number
C is .5527787209, and the hundredth is .3600893021 .  The thousandth
C number should be .2176990509 .
C
C      In order to generate several effectively independent sequences
C with the same generator, it is necessary to know the random number
C for several widely spaced calls.  The I-th random number times 2**22,
C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
C still of the form L*P/8.  In particular we find the I-th random
C number multiplied by 2**22 is given by
C I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
C RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
C
C      Several multipliers have been subjected to the spectral test
C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
C goodness according to the spectral test are
C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
C
C      In the table below LOG10(NU(I)) gives roughly the number of
C random decimal digits in the random numbers considered I at a time.
C C is the primary measure of goodness.  In both cases bigger is better.
C
C                   LOG10 NU(I)              C(I)
C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
C
C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
C   Best
C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
C
C             Input Argument --
C R      If R=0., the next random number of the sequence is generated.
C        If R .LT. 0., the last generated number will be returned for
C          possible use in a restart procedure.
C        If R .GT. 0., the sequence of random numbers will start with
C          the seed R mod 1.  This seed is also returned as the value of
C          RAND provided the arithmetic is done exactly.
C
C             Output Value --
C RAND   a pseudo-random number between 0. and 1.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  RAND
      SAVE IA1, IA0, IA1MA0, IC, IX1, IX0
      DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
      DATA IC /1731/
      DATA IX1, IX0 /0, 0/
C***FIRST EXECUTABLE STATEMENT  RAND
      IF (R.LT.0.) GO TO 10
      IF (R.GT.0.) GO TO 20
C
C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
C                   + IA0*IX0) + IA0*IX0
C
      IY0 = IA0*IX0
      IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
      IY0 = IY0 + IC
      IX0 = MOD (IY0, 2048)
      IY1 = IY1 + (IY0-IX0)/2048
      IX1 = MOD (IY1, 2048)
C
 10   ALEAS = IX1*2048 + IX0
      ALEAS = ALEAS / 4194304.
      RETURN
C
 20   IX1 = MOD(R,1.)*4194304. + 0.5
      IX0 = MOD (IX1, 2048)
      IX1 = (IX1-IX0)/2048
      GO TO 10
C
      END
