C------------------------------------------------
C				                                       random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END
C--------------------------------------
C				normal random numbers
      FUNCTION GAUSS(M)
C--------------------------------------
      X=0.
      DO  I=1,5
         X=X+RANDd(M)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS=X2*(1.-0.01*(3.-X2*X2))
      RETURN
      END

C--------------------------------------   Fourier Transform
      SUBROUTINE SETF67(JBC1,JQ1,jbc,jp,jsl,ll1,k2,k3,k4,k7,
     &                  Qi,jndx,Uf,Vf)
c     -----------------------------------------------------
      INCLUDE 'PMPfour67.h'
      dimension Uf(marr) , Vf(marr)
      dimension Qi(mf67) , jndx(mf67)

      PI=DATAN(1.D0)*4.
       JBC=JBC1
       JQ=JQ1
       IF(JBC.LT.3) GO TO 101
       JQ=JQ-1
  101      K3=2**JQ
       K7=K3/2
       N5=K3/4
       I=1
       JNDX(I)=N5
       QI(I)=0.5*SQRT(2.)
       K=I
       I=I+1
  102     IL=I
       IF(I.EQ.K7) GO TO 104
  103     K1=JNDX(K)/2
       JNDX(I)=K1
       QI(I)=SIN(PI*FLOAT(K1)/FLOAT(K3))
       K1=K7-K1
       I=I+1
       JNDX(I)=K1
       QI(I)=SIN(PI*FLOAT(K1)/FLOAT(K3))
       K=K+1
       I=I+1
       IF(K.EQ.IL) GO TO 102
       GO TO 103
  104     RETURN
       END
C-----------------------------------------
       SUBROUTINE TFOLD(IS,L,ZZZ,k2)
      INCLUDE 'PMPfour67.h'

       DIMENSION ZZZ(MARR)
       IH2=K2/2-1
       DO 100 I=IS,IH2
       I1=I+L
       I2=K2-I+L
       A=ZZZ(I1)
       ZZZ(I1)=A-ZZZ(I2)
       ZZZ(I2)=A+ZZZ(I2)
  100     CONTINUE
       RETURN
       END
C------------------------------------
       SUBROUTINE NEG(I1,I3,I2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)
 
       DO 100 K=I1,I3,I2
       Vf(K)=-Vf(K)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE REVNEG(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       DO 100 I=1,K7
       J=K3+1+I
       K=K4+1-I
       A=Vf(J)
       Vf(J)=-Vf(K)
       Vf(K)=-A
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE ZEERO(L,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       DO 100 I=1,K2
       LI=L+I
       Uf(LI-1)=0.0
  100     CONTINUE
       RETURN
       END
C------------------------------------------
       SUBROUTINE TFOLD1(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       II=K2-1
       DO 100 I=1,II
       I1=JSL+I
       I2=LL1-I
       A=Uf(I1)
       Uf(I1)=A+Uf(I2)
       Uf(I2)=A-Uf(I2)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       JS1=K2
       I=1
       J5=JSL+K2
       IS1=JSL
       IC1=LL1
       JS1=JS1/2
       IF(JS1.NE.1) GO TO 200
       K1=JNDX(I)
       SN=QI(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=SN*(Uf(IC1)-Uf(IS1))
       Vf(K1+1)=Uf(IC0)+ODD1
       K3MK1=K3-K1
       Vf(K3MK1+1)=Uf(IC0)-ODD1
       IF(JBC.LT.3) GO TO 110
       ODD2=SN*(Uf(IC1)+Uf(IS1))
       K3PK1=K3+K1
       Vf(K3PK1+1)=Uf(IS0)+ODD2
       K4MK1=K4-K1
       Vf(K4MK1+1)=-Uf(IS0)+ODD2
  110     RETURN
  200     SN=QI(I)
       IS1=IS1+JS1
       IC1=IC1+JS1
       J3=IS1+JS1
  210     IS0=IS1-JS1
       IC0=IC1-JS1
       ODD1=SN*(Uf(IC1)-Uf(IS1))
       ODD2=SN*(Uf(IC1)+Uf(IS1))
       Uf(IC1)=Uf(IC0)-ODD1
       Uf(IS1)=-Uf(IS0)+ODD2
       Uf(IC0)=Uf(IC0)+ODD1
       Uf(IS0)=Uf(IS0)+ODD2
       IS1=IS1+1
       IC1=IC1+1
       IF(IS1.NE.J3) GO TO 210
       I=I+1
  300     IS1=JSL
       IC1=LL1
       JS1=JS1/2
       IF(JS1.EQ.1) GO TO 400
  310     SN=QI(I)
       I=I+1
       CS=QI(I)
       IS1=IS1+JS1
       IC1=IC1+JS1
       J3=IS1+JS1
  320     IS0=IS1-JS1
       IC0=IC1-JS1
       ODD1=CS*Uf(IC1)-SN*Uf(IS1)
       ODD2=SN*Uf(IC1)+CS*Uf(IS1)
       Uf(IC1)=Uf(IC0)-ODD1
       Uf(IC0)=Uf(IC0)+ODD1
       Uf(IS1)=-Uf(IS0)+ODD2
       Uf(IS0)=Uf(IS0)+ODD2
       IS1=IS1+1
       IC1=IC1+1
       IF(IS1.NE.J3) GO TO 320
       IS1=IS1+JS1
       IC1=IC1+JS1
       J3=IS1+JS1
  330     IS0=IS1-JS1
       IC0=IC1-JS1
       ODD1=SN*Uf(IC1)-CS*Uf(IS1)
       ODD2=CS*Uf(IC1)+SN*Uf(IS1)
       Uf(IC1)=Uf(IC0)-ODD1
       Uf(IC0)=Uf(IC0)+ODD1
       Uf(IS1)=-Uf(IS0)+ODD2
       Uf(IS0)=Uf(IS0)+ODD2
       IS1=IS1+1
       IC1=IC1+1
       IF(IS1.NE.J3) GO TO 330
       I=I+1
       IF(IS1.EQ.J5) GO TO 300
       GO TO 310
  400     K1=JNDX(I)
       SN=QI(I)
       I=I+1
       CS=QI(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=CS*Uf(IC1)-SN*Uf(IS1)
       Vf(K1+1)=Uf(IC0)+ODD1
       K3MK1=K3-K1
       Vf(K3MK1+1)=Uf(IC0)-ODD1
       IF(JBC.LT.3) GO TO 410
       ODD2=SN*Uf(IC1)+CS*Uf(IS1)
       K3PK1=K3+K1
       Vf(K3PK1+1)=Uf(IS0)+ODD2
       K4MK1=K4-K1
       Vf(K4MK1+1)=-Uf(IS0)+ODD2
  410     IS1=IS1+1
       IC1=IC1+1
       K1=JNDX(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=SN*Uf(IC1)-CS*Uf(IS1)
       Vf(K1+1)=Uf(IC0)+ODD1
       K3MK1=K3-K1
       Vf(K3MK1+1)=Uf(IC0)-ODD1
       IF(JBC.LT.3) GO TO 420
       ODD2=CS*Uf(IC1)+SN*Uf(IS1)
       K3PK1=K3+K1
       Vf(K3PK1+1)=Uf(IS0)+ODD2
       K4MK1=K4-K1
       Vf(K4MK1+1)=-Uf(IS0)+ODD2
  420     IS1=IS1+1
       IC1=IC1+1
       I=I+1
       IF(IS1.NE.J5) GO TO 400
       RETURN
       END

c      ------------------------------------------------------
       SUBROUTINE FOUR67(JBC1,JQ1,jp1,jsl1,ll11,k21,k71,
     &                   Qi,jndx,Uf,Vf)
c      ------------------------------------------------------
      INCLUDE 'PMPfour67.h'
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       JBC=JBC1
       JQ=JQ1
       jp = jp1
       jsl = jsl1
       ll1 = ll11
       k2 = k21
       k7 = k71
       A5=0.5*SQRT(2.0)
       K4=2**JQ
       K3=K4
       GO TO (103,103,101,102),JBC
  101     Uf(1)=Uf(1)/2.0
       Uf(K3+1)=Uf(1)
       K2=K3

       CALL TFOLD(0,1,Uf,k2)

       K3=K3/2
       JQ=JQ-1
       GO TO 103
  102     K3=K3/2
       JQ=JQ-1
  103     N5=K3/4
       K7=K3/2
       N11=3*K7
       K31=K3+1
       GO TO(300,400,500,600),JBC
  300     Uf(K31)=0.0
       Uf(1)=0.0
       K2=K3
       DO 301 I=2,JQ

       CALL TFOLD(1,1,Uf,k2)

  301     K2=K2/2
       Vf(K7+1)=Uf(2)
       JF=N5
       JSL=1
       DO 302 JP=2,JQ
       LL1=K2+1

       CALL ZEERO(1,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
       CALL KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       I1=3*JF+1
       I2=4*JF
       I3=I1+(K2/2-1)*I2

       CALL NEG(I1,I3,I2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       K2=K2+K2
  302     JF=JF/2
       RETURN
  400     Uf(1)=Uf(1)/2.0
       Uf(K31)=Uf(K31)/2.0
       K2=K3
       DO 401 I=2,JQ

       CALL TFOLD(0,K31-K2,Uf,k2)

  401     K2=K2/2
       LL1=K31-K2
       A=Uf(LL1)+Uf(LL1+2)
       Vf(1)=-A-Uf(LL1+1)
       Vf(K31)=-A+Uf(LL1+1)
       Vf(K7+1)=Uf(LL1)-Uf(LL1+2)
       DO 402 JP=2,JQ
       JSL=K31-K2
       LL1=JSL-K2
       idum = jsl

       CALL ZEERO(idum,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
       CALL KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

  402     K2=K2+K2

       CALL NEG(1,K31,2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       RETURN
  500     K2=K3
       L2=K4
       DO 501 JP=2,JQ

       CALL TFOLD(1,1,Uf,k2)
       CALL TFOLD(0,L2-K2+1,Uf,k2)

  501     K2=K2/2
       LL1=L2-K2+1
       A=Uf(LL1)+Uf(LL1+2)
       Vf(K7+1)=2.0*(-Uf(LL1)+Uf(LL1+2))
       Vf(1)=2.0*(A+Uf(LL1+1))
       Vf(K31)=2.0*(A-Uf(LL1+1))
       Vf(N11+1)=2.0*Uf(2)
       DO 502 JP=2,JQ
       Uf(K2+1)=2.0*Uf(K2+1)
       JSL=K2+1

       CALL TFOLD1(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       LL1=LL1-K2
       Uf(LL1)=-2.0*Uf(LL1)

       CALL KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

  502     K2=K2+K2
       Vf(1)=Vf(1)*A5
       Vf(K31)=Vf(K31)*A5
       RETURN
  600     Uf(1)=Uf(1)*A5
       Uf(K31)=Uf(K31)*A5
       K2=K3
       DO 601 JP=2,JQ

       CALL TFOLD(0,K31-K2,Uf,k2)
       CALL TFOLD(1,K31,Uf,k2)

  601     K2=K2/2
       LL1=K31-K2
       A=Uf(LL1)+Uf(LL1+2)
       Vf(1)=2.0*(A+Uf(LL1+1))
       Vf(K31)=2.0*(A-Uf(LL1+1))
       Vf(K7+1)=2.0*(-Uf(LL1)+Uf(LL1+2))
       Vf(N11+1)=2.0*Uf(K3+2)
       DO 602 JP=2,JQ
       JSL=K31+K2
       Uf(JSL)=2.0*Uf(JSL)

       CALL TFOLD1(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       LL1=LL1-K2
       Uf(LL1)=-2.0*Uf(LL1)

       CALL KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

  602     K2=K2+K2

       CALL REVNEG(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       idum = k3

       CALL NEG(2,idum,2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

       K2=K4

       CALL TFOLD(1,1,Vf,k2)

       Vf(K4+1)=0.0

       RETURN
       END