C---5----0----5----0----5----0----5----0----5----0----5----0----5----0--
C       10        20        30        40        50        60        70
C___5____0____5____0____5____0____5____0____5____0____5____0____5____0__
C     
C     CDW111.F    :: (Nov. 15, 1993)  standard in DEC; put phi0.
C     CDW110.F    :: (Sep. 28, 1993)  delete unit=11; adjust size.
C     CDW109.F    :: (Sep. 23, 1993)  larger sample.
C     CDW108.F    :: (Aug. 26, 1993)  back to original average method.
C     CDW107.F    :: (Aug. 24, 1993)  put random from NumericalRecipes.
C     CDW107.F    :: (Aug. 16, 1993)  put fftpack for power spectrum.
C     CDW106.F    :: (Aug. 13, 1993)  adapt phasedia.f for step widths.
C                                     call readdata.
C     CDW105.F    :: (Aug. 11, 1993)  check fort.20,21.
C     CDW104.F    :: (Aug. 09, 1993)  to check diff btw 2d and 3d.
C     CDW103.F    :: (Aug. 06, 1993)  check strange diff btw 2d and 3d.
C     CDW102.F    :: (Aug. 03, 1993)  more steps for phase diagram.
C     CDW101.F    :: (Aug. 02, 1993)  correct phase diagram.
C     CDW99.F     :: (Jul. 29, 1993)  choose several ac or several reals.
C     CDW98.F     :: (Jul. 29, 1993)  activaticve run.
C     CDW97.F     :: (Jul. 01, 1993)  New average method.
C     CDW96.F     :: (Jun. 30, 1993)  Disable FFT, enable new average.
C     CDW95.F     :: (Jun. 28, 1993)  The first program for cdw.
C     CDW94.F     :: (Jun. 25, 1993)  The first program for cdw.
C     dPi/dt = Laplacian Pi - h*sin(Pi+Bi)+Edc+Eacsin(wt)
C     CTTERNPIRT4.F: (Jun. 10, 1993)  Giveaway program.
C     Use uniform dist. random #  <--- normal dist. ran. #
C     <X>g=<X>u=0,
C     <X*X>g=s*s,       when Pg(x)=exp(-x*x/2s*s)/s*sqrt(2*pi),
C     <X*X>u=a*a/3,     when Pu(x)=1/2a.
C     From <X*X>g=<X*X>u, we get a=s*sqrt(3) 
C
      PROGRAM CDW
      PARAMETER(ZERO=0.0,PI=3.141592653589793238462643 )
      PARAMETER(MTIM=30010,MSIZ=4096)
      INTEGER        BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /INTGR/ BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /DIMN1/ EE(MSIZ),EDCPRT(MSIZ)
      COMMON /DIMN2/ HH(MSIZ),GRA2(MSIZ),PHIDOT(MSIZ)
      COMMON /AVERG/ CURREN(MTIM),CMIDL(MTIM),CMID0(MTIM),CMIDR(MTIM)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /RAND1/ IHAND,IBETA
      COMMON /RAND2/ H0,DH,BETA(MSIZ)
      COMMON /COUNT/ KARRAY(MTIM)
      COMMON /NOIS1/ ISEED
      COMMON /NOIS2/ TNOISE(MSIZ)
      COMMON /FLAGS/ IPRT,ISAAD,MJCT,NFFT,IPICK,NPROB
      COMMON /EVALU/ ECS(-8192:8192),ESS(-8192:8192)
      COMMON /DATA1/ TEM00,TEMINC,TEM99,IPROB,ISCAN
      COMMON /DATA2/ ISEED0,LSEED0,MSEED0,NSEED0,NREALS
      COMMON /DATA3/ ISTEP,ISTAVG,NSTEP,ISTOFF,IBACK
      COMMON /DATA4/ EDC00,EDC33,EDC66,EDC99,EINC1,EINC2,EINC3
      COMMON /DATA5/ IFLAG,NCHRG,IEVA
      COMMON /DATA6/ EACINC,FREINC,ACPT,ROUND
      COMMON /INCND/ PHI0(MSIZ)
      DIMENSION PHI(MSIZ),AVGCUR(MSIZ),RNDD(MSIZ)
      DIMENSION Z00E(99),Z11E(99),Z21E(99),Z12E(99),Z32E(99),Z13E(99)
      DIMENSION Z23E(99),Z43E(99),Z14E(99),Z34E(99),Z54E(99),Z15E(99)
      DIMENSION Z25E(99),Z35E(99),Z45E(99)
C
C     read data file:
C
      CALL READDATA(IH,JH,KH,PI2,PNV,SQRT3,NNNEEE,NTEMP)
      CALL WRITEOUT1(IFLAG,IHAND,H0,DH,IBETA)
C
C     expected values of steps.
C
      CALL SHAPIRO(WE,EX00,EX11,EX21,EX12,EX32,EX13,EX23,
     *      EX43,EX14,EX34,EX54,EX15,EX25,EX35,EX45,EX65)
C
C
C     For many realizations:
C
      ISEED=ISEED0
      LSEED=LSEED0
      MSEED=MSEED0
      NSEED=NSEED0
      DO 9988 IRLS=1,NREALS  
      WRITE(6,*)  
      WRITE(6,*) ' #   ISEED      LSEED      MSEED      NSEED  '
      WRITE(6,'(I3,4I11)') IRLS,ISEED,LSEED,MSEED,NSEED
      EDC0=EDC00
      TEMP0=TEM00
C
C     disorder in beta's from impurities:
C
      IF(IBETA.EQ.0) THEN
         DO I=1,LXLYLZ
            BETA(I)=0.0
         ENDDO
      ELSE
         CALL RANINI(NSEED)
         CALL RANUNI(NSEED,LXLYLZ,RNDD)
         DO I=1,LXLYLZ
            BETA(I)=2.0*(RNDD(I)-0.50)*PI
         ENDDO
      ENDIF
C
C     disorder in hh's from defects:
C
      IF(IHAND.NE.0) CALL STRGTH(IRLS,IHAND,MSEED,H0,DH)
      CALL WRITEOUT2(ISCAN,ISAAD,IHAND,H0,DH,IBETA)
C
C     check histeresis:
C
      DO 9977 JTRACE=1,IBACK
      ISIGN=-(-1)**JTRACE
      SIGN =ISIGN
      IF(ISCAN.LT.0) THEN
         NNNNN=NNNEEE
      ELSEIF(ISCAN.GT.0) THEN
         NNNNN=NTEMP-1
      ENDIF
C
C     assign zero number of data for each step:
C
      CALL ZEROS(I00,I11,I21,I12,I32,I13,I23,I43,I14,I34,I54,
     *           I15,I25,I35,I45,I65)

C
C     scan horizontal values (dc fields or temperature):
C
      DO 9955 ICURVE=1,NNNNN
      IF(ISCAN.LT.0) THEN
         ICVSGN=(NNNNN+1)*(1-ISIGN)/2+ISIGN*ICURVE
         EDC0=EE(ICVSGN)
      ELSE
         EDC0=EDC00
      ENDIF
      WRITE(6,'('' Trace='',i2,i4,'': Edc='',f6.3,'', T='',f6.3)')
     *       JTRACE,ICURVE,Edc0,TEMP0
C
C     If TEMP = 0, then TNS = 0
C     If IPROB=+1, then use uniform   random noise.
C        IPROB=-1, then use Gaussian  random noise.
C
      NR=3*N*ISTEP
      IF(TEMP0.NE.ZERO) THEN
         CALL RANINI(ISEED)
         TEMP=SQRT(2.0*TEMP0/DELT)
      ELSE
         DO I=1,LXLYLZ
            TNOISE(I)=ZERO
         ENDDO
         TEMP=ZERO
      ENDIF
CC
      CALL INICOND(ICURVE,NCHRG,LSEED,PHI)
CC
      KKK=0
      IF(ISTAVG*ISTEP.EQ.1) GOTO 6900
      DO 690 K0=1,ISTAVG*ISTEP-1,ISTEP
      DO 660 KK=1,ISTEP
      KKK=KKK+1
      IF(TEMP.NE.0.0) CALL NOISDIS(TEMP,IPROB,SQRT3)
      TIMEKK=KKK*DELT
      IF(IEVA.GT.0) THEN 
         EIN=EDC0+EAC*SIN(WE*TIMEKK)
      ELSE
         ARGE=WE*TIMEKK*PNV
         IWE=2048.0*(ARGE-IFIX(ARGE))
         EIN=EDC0+EAC*ESS(IWE)
      ENDIF
      CALL GRAGRA(PHI,GRA2)
      IF(IEVA.GT.0) THEN
        DO I=1,N
           PHIDOT(I)=GRA2(I)-HH(I)*SIN(PHI(I)+BETA(I))+TNOISE(I)
        ENDDO
      ELSE
        DO I=1,N
           ARGU=(PHI(I)+BETA(I))*PNV
           IARG=2048.0*(ARGU-IFIX(ARGU))
           PHIDOT(I)=GRA2(I)-HH(I)*ESS(IARG)+TNOISE(I)
        ENDDO
      ENDIF
C
C     Use periodic boundary condition for all direction.
C
      IF(IPICK.GT.0) THEN			!reg vol-in and cur-pick.
         NPROB=N
         DO I=1,N
            PHIDOT(I)=PHIDOT(I)+EIN 
         ENDDO
      ELSEIF(IPICK.EQ.0.OR.IPICK.EQ.-1) THEN	!TopPlaneVol/BottomPlaneCur.
         NPROB=LXLY
         DO IJ=1,LXLY
            I1=(LZ-1)*LXLY+IJ
            PHIDOT(I1)=PHIDOT(I1)+EIN
         ENDDO
      ELSE					!reg vol-in and transcur-pick.
         DO I=1,N
            PHIDOT(I)=PHIDOT(I)+EIN 
         ENDDO
      ENDIF
C
C     Calculate new phases, PHI's.
C
      DO I=1,LXLYLZ
         PHI(I)=PHI(I)+DELT*PHIDOT(I)
C        PHI(I)=PHI(I)-PI2*IFIX(PHI(I)*PNV)
      ENDDO
  660 CONTINUE
  690 CONTINUE
 6900 CONTINUE
C
C     Average of PHIDOT over NDATA:
C
C     average = (1/N) \sum_{i=1}^{N} dPi/dt Dt 
C             = (1/N) (P_{N}-P_{0}) = (1/N) (P_{N+1}-P_{1})
C
      VEL$AV=0.0
      VEL$ML=0.0
      VEL$MD=0.0
      VEL$MR=0.0
      KCNT=0
      NLAST=NSTEP
      DO 695 K0=ISTAVG*ISTEP,NLAST,ISTEP
      DO 665 KK=1,ISTEP
      KKK=KKK+1
      IF(TEMP.NE.0.0) CALL NOISDIS(TEMP,IPROB,SQRT3)
      TIMEKK=KKK*DELT
      IF(IEVA.GT.0) THEN
        EIN=EDC0+EAC*SIN(WE*TIMEKK)
      ELSE 
        ARGE=WE*TIMEKK*PNV
        IWE=2048.0*(ARGE-IFIX(ARGE))
        EIN=EDC0+EAC*ESS(IWE)
      ENDIF
C
C     If ISTOFF<0: regular
C     If ISTOFF>0: off Edc ar ISTOFF*ISTEP
C
      CALL GRAGRA(PHI,GRA2)
      IF(IEVA.GT.0) THEN
        DO I=1,N
           PHIDOT(I)=GRA2(I)-HH(I)*SIN(PHI(I)+BETA(I))+TNOISE(I)
        ENDDO
      ELSE
        DO I=1,N
           ARGU=(PHI(I)+BETA(I))*PNV
           IARG=2048.0*(ARGU-IFIX(ARGU))
           PHIDOT(I)=GRA2(I)-HH(I)*ESS(IARG)+TNOISE(I)
        ENDDO
      ENDIF
C
C     Use periodic boundary condition for all direction.
C
      IF(IPICK.GT.0) THEN			!reg vol-in and cur-pick.
         DO I=1,N
            PHIDOT(I)=PHIDOT(I)+EIN 
         ENDDO
      ELSEIF(IPICK.EQ.0.OR.IPICK.EQ.-1) THEN	!TopPlaneVol/BottomPlaneCur.
         DO IJ=1,LXLY
            I1=(LZ-1)*LXLY+IJ
            PHIDOT(I1)=PHIDOT(I1)+EIN
         ENDDO
      ELSE					!reg vol-in and transcur-pick.
         DO I=1,N
            PHIDOT(I)=PHIDOT(I)+EIN 
         ENDDO
      ENDIF
C
C     Calculate new phases, PHI's.
C
      DO I=1,LXLYLZ
         PHI(I)=PHI(I)+DELT*PHIDOT(I)
C        PHI(I)=PHI(I)-PI2*IFIX(PHI(I)*PNV)
      ENDDO
  665 CONTINUE
C
      IF(IPICK.GT.0) THEN		!regular.
         KCNT=KCNT+1
         CMIDL(KCNT)=PHIDOT((KH-1)*LXLY+(JH-1)*LX+IH-1)
         CMID0(KCNT)=PHIDOT((KH-1)*LXLY+(JH-1)*LX+IH  )
         CMIDR(KCNT)=PHIDOT((KH-1)*LXLY+(JH-1)*LX+IH+1)
      ELSEIF(IPICK.EQ.0) THEN		!top voltage/bottom current.
         KCNT=KCNT+1
         CMIDL(KCNT)=PHIDOT((JH-1)*LX+IH-1)
         CMID0(KCNT)=PHIDOT((JH-1)*LX+IH  )
         CMIDR(KCNT)=PHIDOT((JH-1)*LX+IH+1)
      ELSEIF(IPICK.EQ.-1) THEN		!top voltage/trans-bottom current.
         WRITE(6,*) 'irrelevant in this business'
         STOP
      ELSE				!reg volt/transcur(Hall current?).
         WRITE(6,*) 'irrelevant in this business'
         STOP
      ENDIF
      KARRAY(KCNT)=KKK
      SUM=0.0
      DO I=1,NPROB
         SUM=SUM+PHIDOT(I)
      ENDDO
      CURREN(KCNT)=SUM/NPROB
      VEL$AV=VEL$AV+CURREN(KCNT)
      VEL$ML=VEL$ML+CMIDL(KCNT)
      VEL$MD=VEL$MD+CMID0(KCNT)
      VEL$MR=VEL$MR+CMIDR(KCNT)
  695 CONTINUE
C
C     Average the currents:all, midleft, mid, and midright.
C
      NDATA=KCNT
      AVG$A=VEL$AV/NDATA
      AVG$L=VEL$ML/NDATA
      AVG$M=VEL$MD/NDATA
      AVG$R=VEL$MR/NDATA
      IF(ISCAN.LT.0) THEN
C        WRITE(6,795) EDC0 ,AVG$A,AVG$L,AVG$M,AVG$R
         WRITE(9,795) EDC0 ,AVG$A,AVG$L,AVG$M,AVG$R
      ELSE
C        WRITE(6,795) TEMP0,AVG$A,AVG$L,AVG$M,AVG$R
         WRITE(9,795) TEMP0,AVG$A,AVG$L,AVG$M,AVG$R
      ENDIF
C
C     To get phase diagram for IFLAG=<0.
C
      AVGCUR(ICURVE)=AVG$A
      IF(IFLAG.LE.0) THEN
         IF(ABS(AVG$A).LT.ACPT) THEN
            I00=I00+1
            Z00E(I00)=EDC0
         ELSEIF(ABS(AVG$A-EX11).LT.ACPT) THEN
            I11=I11+1
            Z11E(I11)=EDC0
         ELSEIF(ABS(AVG$A-EX12).LT.ACPT) THEN
            I12=I12+1
            Z12E(I12)=EDC0
         ELSEIF(ABS(AVG$A-EX32).LT.ACPT) THEN
            I32=I32+1
            Z32E(I32)=EDC0
         ELSEIF(ABS(AVG$A-EX13).LT.ACPT) THEN
            I13=I13+1
            Z13E(I13)=EDC0
         ELSEIF(ABS(AVG$A-EX23).LT.ACPT) THEN
            I23=I23+1
            Z23E(I23)=EDC0
         ELSEIF(ABS(AVG$A-EX43).LT.ACPT) THEN
            I43=I43+1
            Z43E(I43)=EDC0
         ELSEIF(ABS(AVG$A-EX14).LT.ACPT) THEN
            I14=I14+1
            Z14E(I14)=EDC0
         ELSEIF(ABS(AVG$A-EX34).LT.ACPT) THEN
            I34=I34+1
            Z34E(I34)=EDC0
         ELSEIF(ABS(AVG$A-EX54).LT.ACPT) THEN
            I54=I54+1
            Z54E(I54)=EDC0
         ENDIF    
      ENDIF
C
C     Get last phases at this field to use them as an initial 
C     condition for the next calculation.
C
      DO I=1,N
         PHI(I)=PHI(I)-PI2*IFIX(PHI(I)*PNV)
         PHI0(I)=PHI(I)
      ENDDO
C11   WRITE(11,7530)  EDC0,TEMP0
C11   WRITE(11,7520) (PHI(I),I=1,N)
C11   CLOSE(UNIT=11)
C
      IF(IPRT.GT.99) THEN
        CALL CONVVTFR(NDATA)
      ELSEIF(IPRT.GT.0) THEN
        DO I=1,IPRT
           ABSPRT=ABS(EDC0-EDCPRT(I))
           IF(ABSPRT.LT.1.0E-05) THEN
              CALL CONVVTFR(NDATA)
              GOTO 994
           ENDIF
        ENDDO
      ENDIF
  994 CONTINUE
      IF(ISCAN.GT.0) THEN
        TEMP0=TEMP0+SIGN*TEMINC
      ENDIF
 9955 CONTINUE                ! end_of_field scan.
      EDCINC=ABS(EE(2)-EE(1))
      CALL INTERP(NNNNN,EE,AVGCUR,EN00,WD00,
     * EX11,BE11,EN11,WD11, EX12,BE12,EN12,WD12, EX32,BE32,EN32,WD32,
     * EX13,BE13,EN13,WD13, EX23,BE23,EN23,WD23, EX43,BE43,EN43,WD43,
     * EX14,BE14,EN14,WD14, EX34,BE34,EN34,WD34, EX54,BE54,EN54,WD54)
      IF(I00.GT.0) CALL MAXMIN(I00,Z00E,DUMY,EN00,WD00,EDCINC,ROUND)
      BE00=0.0    
      WD00=EN00-BE00
      IF(I11.GT.0) CALL MAXMIN(I11,Z11E,BE11,EN11,WD11,EDCINC,ROUND)
      IF(I21.GT.0) CALL MAXMIN(I21,Z21E,BE21,EN21,WD21,EDCINC,ROUND)
      IF(I12.GT.0) CALL MAXMIN(I12,Z12E,BE12,EN12,WD12,EDCINC,ROUND)
      IF(I32.GT.0) CALL MAXMIN(I32,Z32E,BE32,EN32,WD32,EDCINC,ROUND)
      IF(I13.GT.0) CALL MAXMIN(I13,Z13E,BE13,EN13,WD13,EDCINC,ROUND)
      IF(I23.GT.0) CALL MAXMIN(I23,Z23E,BE23,EN23,WD23,EDCINC,ROUND)
      IF(I43.GT.0) CALL MAXMIN(I43,Z43E,BE43,EN43,WD43,EDCINC,ROUND)
      IF(I14.GT.0) CALL MAXMIN(I14,Z14E,BE14,EN14,WD14,EDCINC,ROUND)
      IF(I34.GT.0) CALL MAXMIN(I34,Z34E,BE34,EN34,WD34,EDCINC,ROUND)
      IF(I54.GT.0) CALL MAXMIN(I54,Z54E,BE54,EN54,WD54,EDCINC,ROUND)
      IF(IFLAG.LT.0) THEN
         WRITE(20,8203) EAC,BE00,EN00,BE11,EN11,BE12,EN12,
     *                      BE13,EN13,BE23,EN23
         WRITE(21,8213) EAC,BE14,EN14,BE34,EN34,BE32,EN32,
     *                      BE43,EN43,BE54,EN54
         WRITE(22,8223) EAC,WD00,WD11,WD12,WD13,WD23,WD14,
     *                      WD34,WD32,WD43,WD54
      ELSE
         WRITE(20,8204) FREQE,BE00,EN00,BE11,EN11,BE12,EN12,
     *                        BE13,EN13,BE23,EN23
         WRITE(21,8214) FREQE,BE14,EN14,BE34,EN34,BE32,EN32,
     *                        BE43,EN43,BE54,EN54
         WRITE(22,8224) FREQE,WD00,WD11,WD12,WD13,WD23,WD14,
     *                        WD34,WD32,WD43,WD54
      ENDIF
 9977 CONTINUE			! end_of_trace scan.
      IF(IFLAG.GT.0) THEN
         ISEED=ISEED+10000
         LSEED=LSEED+10000
         MSEED=MSEED+10000
         NSEED=NSEED+10000
      ELSEIF(IFLAG.EQ.0) THEN	!scan different ac frequency.
         ISEED=ISEED0
         LSEED=LSEED0
         MSEED=MSEED0
         NSEED=NSEED0
         FREQF=FREQF+FREINC
      ELSEIF(IFLAG.LT.0) THEN	!scan different ac fields.
         ISEED=ISEED0
         LSEED=LSEED0
         MSEED=MSEED0
         NSEED=NSEED0
         EAC=EAC+EACINC
      ENDIF
 9988 CONTINUE			! end_of_realization scan.
C
  760 FORMAT(' LX LY LZ   *CCX  *CCY  *CCZ   BX  BY  BZ   *DELT')
  770 FORMAT(1X,I2,I3,I3,1X,3F6.2,1X,3I4,1X,F7.3)
  780 FORMAT('  *Edc*  Temp  *Ho*  *dH* IHAND IBETA  *Eac*  Frequ'/
     *          F7.3, F6.2, F6.2, F6.2,  I4,   I6,   F9.3,  F7.3   )
  790 FORMAT('  ***Edc*** **Current** **Midleft***Middle***Midright**')
  791 FORMAT('  ***Temp** **Current** **Midleft***Middle***Midright**')
  795 FORMAT(F10.5,F12.6,3F10.4)
 7029 FORMAT(2X,'Idc',9X,'m=1/1',7X,'m=2/1',7X,'m=3/1',7X,'m=4/1')
 7030 FORMAT(2X,'Idc',9X,'m=1/2',7X,'m=3/2',7X,'m=5/2' )
 7031 FORMAT(2X,'Idc',9X,'m=1/3',7X,'m=2/3',7X,'m=4/3' )
 7032 FORMAT(2X,'Idc',9X,'m=1/4',7X,'m=3/4',7X,'m=5/4' )
 7530 FORMAT(' Current =',F10.5,',  Temperature =',F10.5)
 7520 FORMAT(7F11.6)
 8200 FORMAT(' **Eac*  **Step 0/0****  **Step 1/1****  **Step 1/2****',
     *                              '  **Step 1/3****  **Step 2/3****')
 8203 FORMAT(F7.4,10F8.4) 
 8210 FORMAT(' **Eac*  **Step 1/4****  **Step 3/4****  **Step 3/2****',
     *                              '  **Step 4/3****  **Step 5/4****')
 8213 FORMAT(F7.4,10F8.4)
 8220 FORMAT('  **Eac*  **0/0**  **1/1**  **1/2**  **1/3**  **2/3**',
     *               '  **1/4**  **3/4**  **3/2**  **4/3**  **5/4**')
 8223 FORMAT(F8.4,10F9.5)
 8201 FORMAT(' *Freq*  **Step 0/0****  **Step 1/1****  **Step 1/2****',
     *       '  **Step 1/3****  **Step 2/3****')
 8204 FORMAT(F7.4,10F8.4)
 8211 FORMAT(' *Freq*  **Step 1/4****  **Step 3/4****  **Step 3/2****',
     *                              '  **Step 4/3****  **Step 5/4****')
 8214 FORMAT(F7.4,10F8.4) 
 8221 FORMAT('  *Freq*  **0/0**  **1/1**  **1/2**  **1/3**  **2/3**',
     *               '  **1/4**  **3/4**  **3/2**  **4/3**  **5/4**')
 8224 FORMAT(F8.4,10F9.5)
      STOP
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     calculate Laplacian:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE GRAGRA(PHI,GRA2)
      PARAMETER(MSIZ=4096)
      INTEGER        BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /INTGR/ BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      DIMENSION GRA2(MSIZ),PHI(MSIZ)
      IF(LZ.EQ.1) THEN            !Two dimensional sample.
        DO J=2,LY-1
        DO I=2,LX-1
           IJ=(J-1)*LX+I
           GRA2(IJ)=PHI(IJ+1)+PHI(IJ-1)+PHI(IJ+LX)+PHI(IJ-LX)-4*PHI(IJ)
        ENDDO
        ENDDO
        DO I=2,LX-1   !  bottom-front and bottom-back row:
           I1=I
           I2=(LY-1)*LX+I
           GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I2)-4*PHI(I1)
           GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I1)+PHI(I2-LX)-4*PHI(I2)
        ENDDO
        DO J=2,LY-1   !  bottom-left and bottom-right column:
           I1=(J-1)*LX+1
           I2=J*LX
           GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I1-LX)-4*PHI(I1)
           GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)-4*PHI(I2)
        ENDDO
        GRA2(BFL)=PHI(BFL+1)+PHI(BFR)+PHI(BFL+LX)+PHI(BBL)-4*PHI(BFL)
        GRA2(BFR)=PHI(BFL)+PHI(BFR-1)+PHI(BFR+LX)+PHI(BBR)-4*PHI(BFR)
        GRA2(BBL)=PHI(BBL+1)+PHI(BBR)+PHI(BFL)+PHI(BBL-LX)-4*PHI(BBL)
        GRA2(BBR)=PHI(BBL)+PHI(BBR-1)+PHI(BFR)+PHI(BBR-LX)-4*PHI(BBR)
        IF(LY.EQ.1) THEN   ! one dimension with e-field on x-axis.
           DO I=2,LX-1
              GRA2(I)=PHI(I+1)+PHI(I-1)-2.0*PHI(I)
           ENDDO
           GRA2(1)=PHI(2)+PHI(LX)-2*PHI(1)
           GRA2(LX)=PHI(1)+PHI(LX-1)-2*PHI(LX)
        ENDIF
      ELSE                 ! Three dimensional sample.
        IF(IBCZ.LT.0) THEN
          DO K=2,LZ-1
          DO J=2,LY-1
          DO I=2,LX-1
             IJK=(K-1)*LXLY+(J-1)*LX+I
             GRA2(IJK)=PHI(IJK+1)+PHI(IJK-1)+PHI(IJK+LX)+PHI(IJK-LX)
     *                +PHI(IJK+LXLY)+PHI(IJK-LXLY)-6*PHI(IJK)
          ENDDO
          ENDDO
          ENDDO
          DO K=2,LZ-1                ! front and back face:
          DO I=2,LX-1
             I1=(K-1)*LXLY+I
             I2=(K-1)*LXLY+(LY-1)*LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I2)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I1)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
          ENDDO
          ENDDO
          DO K=2,LZ-1                ! left and right face.
          DO J=2,LY-1
             I1=(K-1)*LXLY+(J-1)*LX+1
             I2=(K-1)*LXLY+J*LX
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
          ENDDO
          ENDDO
          DO J=2,LY-1                ! bottom and top face:
          DO I=2,LX-1
             I1=(J-1)*LX+I
             I2=(LZ-1)*LXLY+(J-1)*LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)+PHI(I2)-6*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I1)+PHI(I2-LXLY)-6*PHI(I2)
          ENDDO
          ENDDO
          DO I=2,LX-1                 !  lines || x-axis:
             I1=I
             I2=LXLY-LX+I
             I3=LXLYLZ-LXLY+I
             I4=LXLYLZ-LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I2)
     *               +PHI(I1+LXLY)+PHI(I3)-6*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I1)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I4)-6*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I3-1)+PHI(I3+LX)+PHI(I4)
     *               +PHI(I1)+PHI(I3-LXLY)-6*PHI(I3)
             GRA2(I4)=PHI(I4+1)+PHI(I4-1)+PHI(I3)+PHI(I4-LX)
     *               +PHI(I2)+PHI(I4-LXLY)-6*PHI(I4) 
          ENDDO
          DO J=2,LY-1                 !  lines || y-axis:
             I1=(J-1)*LX+1
             I2=J*LX
             I3=(LZ-1)*LXLY+(J-1)*LX+1
             I4=(LZ-1)*LXLY+J*LX
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)+PHI(I3)-6*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I4)-6*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I4)+PHI(I3+LX)+PHI(I3-LX)
     *               +PHI(I1)+PHI(I3-LXLY)-6*PHI(I3)
             GRA2(I4)=PHI(I3)+PHI(I4-1)+PHI(I4+LX)+PHI(I4-LX)
     *               +PHI(I2)+PHI(I4-LXLY)-6*PHI(I4)
          ENDDO
          DO K=2,LZ-1                 !  lines || z-axis:
             I1=(K-1)*LXLY+1
             I2=(K-1)*LXLY+LX
             I3=(K-1)*LXLY+(LY-1)*LX+1
             I4=(K-1)*LXLY+LXLY
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I3)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I4)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I4)+PHI(I1)+PHI(I3-LX)
     *               +PHI(I3+LXLY)+PHI(I3-LXLY)-6*PHI(I3)
             GRA2(I4)=PHI(I3)+PHI(I4-1)+PHI(I2)+PHI(I4-LX)
     *               +PHI(I4+LXLY)+PHI(I4-LXLY)-6*PHI(I4)
          ENDDO
          GRA2(BFL)=PHI(BFL+1)+PHI(BFR)+PHI(BFL+LX)+PHI(BBL)
     *             +PHI(BFL+LXLY)+PHI(TFL)-6*PHI(BFL)
          GRA2(BFR)=PHI(BFL)+PHI(BFR-1)+PHI(BFR+LX)+PHI(BBR)
     *             +PHI(BFR+LXLY)+PHI(TFR)-6*PHI(BFR)
          GRA2(BBL)=PHI(BBL+1)+PHI(BBR)+PHI(BFL)+PHI(BBL-LX)
     *             +PHI(BBL+LXLY)+PHI(TBL)-6*PHI(BBL)
          GRA2(BBR)=PHI(BBL)+PHI(BBR-1)+PHI(BFR)+PHI(BBR-LX)
     *             +PHI(BBR+LXLY)+PHI(TBR)-6*PHI(BBR)
          GRA2(TFL)=PHI(TFL+1)+PHI(TFR)+PHI(TFL+LX)+PHI(TBL)
     *             +PHI(BFL)+PHI(TFL-LXLY)-6*PHI(TFL)
          GRA2(TFR)=PHI(TFL)+PHI(TFR-1)+PHI(TFR+LX)+PHI(TBR)
     *             +PHI(BFR)+PHI(TFR-LXLY)-6*PHI(TFR)
          GRA2(TBL)=PHI(TBL+1)+PHI(TBR)+PHI(TFL)+PHI(TBL-LX)
     *             +PHI(BBL)+PHI(TBL-LXLY)-6*PHI(TBL)
          GRA2(TBR)=PHI(TBL)+PHI(TBR-1)+PHI(TFR)+PHI(TBR-LX)
     *             +PHI(BBR)+PHI(TBR-LXLY)-6*PHI(TBR)
        ELSE  ! 3-dimenaion and free boundary on z-axis
          IF(LZ.EQ.2) GOTO 2530
          DO K=2,LZ-1
          DO J=2,LY-1
          DO I=2,LX-1
             IJK=(K-1)*LXLY+(J-1)*LX+I
             GRA2(IJK)=PHI(IJK+1)+PHI(IJK-1)+PHI(IJK+LX)+PHI(IJK-LX)
     *                +PHI(IJK+LXLY)+PHI(IJK-LXLY)-6*PHI(IJK)
          ENDDO
          ENDDO
          ENDDO
          DO K=2,LZ-1                ! front and back face:
          DO I=2,LX-1
             I1=(K-1)*LXLY+I
             I2=(K-1)*LXLY+(LY-1)*LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I2)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I1)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
          ENDDO
          ENDDO
          DO K=2,LZ-1                ! left and right face.
          DO J=2,LY-1
             I1=(K-1)*LXLY+(J-1)*LX+1
             I2=(K-1)*LXLY+J*LX
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
          ENDDO
          ENDDO
 2530     CONTINUE
          DO J=2,LY-1                ! bottom and top face:
          DO I=2,LX-1
             I1=(J-1)*LX+I
             I2=(LZ-1)*LXLY+(J-1)*LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)-5*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I2-LXLY)-5*PHI(I2)
          ENDDO
          ENDDO
          DO I=2,LX-1                 !  lines || x-axis:
             I1=I
             I2=LXLY-LX+I
             I3=LXLYLZ-LXLY+I
             I4=LXLYLZ-LX+I
             GRA2(I1)=PHI(I1+1)+PHI(I1-1)+PHI(I1+LX)+PHI(I2)
     *               +PHI(I1+LXLY)-5*PHI(I1)
             GRA2(I2)=PHI(I2+1)+PHI(I2-1)+PHI(I1)+PHI(I2-LX)
     *               +PHI(I2+LXLY)-5*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I3-1)+PHI(I3+LX)+PHI(I4)
     *               +PHI(I3-LXLY)-5*PHI(I3)
             GRA2(I4)=PHI(I4+1)+PHI(I4-1)+PHI(I3)+PHI(I4-LX)
     *               +PHI(I4-LXLY)-5*PHI(I4) 
          ENDDO
          DO J=2,LY-1                 !  lines || y-axis:
             I1=(J-1)*LX+1
             I2=J*LX
             I3=(LZ-1)*LXLY+(J-1)*LX+1
             I4=(LZ-1)*LXLY+J*LX
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I1-LX)
     *               +PHI(I1+LXLY)-5*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I2-LX)
     *               +PHI(I2+LXLY)-5*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I4)+PHI(I3+LX)+PHI(I3-LX)
     *               +PHI(I3-LXLY)-5*PHI(I3)
             GRA2(I4)=PHI(I3)+PHI(I4-1)+PHI(I4+LX)+PHI(I4-LX)
     *               +PHI(I4-LXLY)-5*PHI(I4)
          ENDDO
          IF(LZ.EQ.2) GOTO 2550
          DO K=2,LZ-1                 !  lines || z-axis:
             I1=(K-1)*LXLY+1
             I2=(K-1)*LXLY+LX
             I3=(K-1)*LXLY+(LY-1)*LX+1
             I4=(K-1)*LXLY+LXLY
             GRA2(I1)=PHI(I1+1)+PHI(I2)+PHI(I1+LX)+PHI(I3)
     *               +PHI(I1+LXLY)+PHI(I1-LXLY)-6*PHI(I1)
             GRA2(I2)=PHI(I1)+PHI(I2-1)+PHI(I2+LX)+PHI(I4)
     *               +PHI(I2+LXLY)+PHI(I2-LXLY)-6*PHI(I2)
             GRA2(I3)=PHI(I3+1)+PHI(I4)+PHI(I1)+PHI(I3-LX)
     *               +PHI(I3+LXLY)+PHI(I3-LXLY)-6*PHI(I3)
             GRA2(I4)=PHI(I3)+PHI(I4-1)+PHI(I2)+PHI(I4-LX)
     *               +PHI(I4+LXLY)+PHI(I4-LXLY)-6*PHI(I4)
          ENDDO
 2550     CONTINUE
          GRA2(BFL)=PHI(BFL+1)+PHI(BFR)+PHI(BFL+LX)+PHI(BBL)
     *             +PHI(BFL+LXLY)-5*PHI(BFL)
          GRA2(BFR)=PHI(BFL)+PHI(BFR-1)+PHI(BFR+LX)+PHI(BBR)
     *             +PHI(BFR+LXLY)-5*PHI(BFR)
          GRA2(BBL)=PHI(BBL+1)+PHI(BBR)+PHI(BFL)+PHI(BBL-LX)
     *             +PHI(BBL+LXLY)-5*PHI(BBL)
          GRA2(BBR)=PHI(BBL)+PHI(BBR-1)+PHI(BFR)+PHI(BBR-LX)
     *             +PHI(BBR+LXLY)-5*PHI(BBR)
          GRA2(TFL)=PHI(TFL+1)+PHI(TFR)+PHI(TFL+LX)+PHI(TBL)
     *             +PHI(TFL-LXLY)-5*PHI(TFL)
          GRA2(TFR)=PHI(TFL)+PHI(TFR-1)+PHI(TFR+LX)+PHI(TBR)
     *             +PHI(TFR-LXLY)-5*PHI(TFR)
          GRA2(TBL)=PHI(TBL+1)+PHI(TBR)+PHI(TFL)+PHI(TBL-LX)
     *             +PHI(TBL-LXLY)-5*PHI(TBL)
          GRA2(TBR)=PHI(TBL)+PHI(TBR-1)+PHI(TFR)+PHI(TBR-LX)
     *             +PHI(TBR-LXLY)-5*PHI(TBR)
        ENDIF
      ENDIF
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     disorder in h's from defects (not from impurities).
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE STRGTH(IRLS,IHAND,MSEED,H0,DH)
      PARAMETER(MSIZ=4096)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /DIMN2/ HH(MSIZ),GRA2(MSIZ),PHIDOT(MSIZ)
      DIMENSION RNDD(MSIZ)
      IF(IHAND.EQ.1) THEN
        WRITE(6,*) 'uniformly distributed random strength h'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LXLYLZ,RNDD)
        DO I=1,LXLYLZ
           HH(I)=2.0*(RNDD(I)-0.50)*DH+H0
        ENDDO
      ELSEIF(IHAND.EQ.21) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on x'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LYLZ,RNDD)
        DO K=1,LZ
        DO J=1,LY
           JK=(K-1)*LY+J
           RANH=2.0*(RNDD(JK)-0.50)*DH+H0
           DO I=1,LX
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH
           ENDDO
        ENDDO
        ENDDO
      ELSEIF(IHAND.EQ.22) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on y'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LXLZ,RNDD)
        DO K=1,LZ
        DO I=1,LX
           IK=(K-1)*LX+I
           RANH=2.0*(RNDD(IK)-0.50)*DH+H0
           DO J=1,LY
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH
           ENDDO
        ENDDO
        ENDDO
      ELSEIF(IHAND.EQ.23) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on z'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LXLY,RNDD)
        DO J=1,LY
        DO I=1,LX
           IJ=(J-1)*LX+I
           RANH=2.0*(RNDD(IJ)-0.50)*DH+H0
           DO K=1,LZ
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH 
           ENDDO
        ENDDO
        ENDDO
      ELSEIF(IHAND.EQ.31) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on xy'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LZ,RNDD)
        DO K=1,LZ
           RANH=2.0*(RNDD(K)-0.50)*DH+H0
           DO J=1,LY
           DO I=1,LX
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH 
           ENDDO
        ENDDO
        ENDDO
      ELSEIF(IHAND.EQ.32) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on yz'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LX,RNDD)
        DO I=1,LX
           RANH=2.0*(RNDD(I)-0.50)*DH+H0
           DO K=1,LZ
           DO J=1,LY
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH 
           ENDDO
        ENDDO
        ENDDO
      ELSEIF(IHAND.EQ.33) THEN
        WRITE(6,*) 'uniformly dis ran strength h, but correlated on zx'
        CALL RANINI(MSEED)
        CALL RANUNI(MSEED,LY,RNDD)
        DO J=1,LY
           RANH=2.0*(RNDD(J)-0.50)*DH+H0
           DO K=1,LZ
           DO I=1,LX
              IJK=(K-1)*LXLY+(J-1)*LX+I
              HH(IJK)=RANH 
           ENDDO
        ENDDO
        ENDDO
      ELSE   !(IHAND=other than 0,1,21,22,23,31,32,33)
        IF(IRLS.GT.1) GOTO 45
        WRITE(6,*)
        WRITE(6,*) 'Check IHAND.'
        WRITE(6,*)
        STOP
      ENDIF
   45 RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     expected values of Shapiro steps: (2*pi*freq*p/q)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SHAPIRO(WE,EX00,EX11,EX21,EX12,EX32,EX13,EX23,
     *            EX43,EX14,EX34,EX54,EX15,EX25,EX35,EX45,EX65)
      EX00=0.0
      EX11=WE*1.000
      EX21=WE*2.000
      EX12=WE*0.500
      EX32=WE*1.500
      EX13=WE*0.333333333333
      EX23=WE*0.666666666667
      EX43=WE*1.333333333333
      EX14=WE*0.250
      EX34=WE*0.750
      EX54=WE*1.250
      EX15=WE*0.200
      EX25=WE*0.400
      EX35=WE*0.600
      EX45=WE*0.800
      EX65=WE*1.200
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     initial interpolated values of Edc's corresponding to starting
C     and ending of each step.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INTERP(INET,ARREDC,ARRCUR,EN00,WD00,
     *  EX11,BE11,EN11,WD11, EX12,BE12,EN12,WD12, EX32,BE32,EN32,WD32,
     *  EX13,BE13,EN13,WD13, EX23,BE23,EN23,WD23, EX43,BE43,EN43,WD43,
     *  EX14,BE14,EN14,WD14, EX34,BE34,EN34,WD34, EX54,BE54,EN54,WD54)
      DIMENSION ARREDC(10000), ARRCUR(10000)
      EN00=0.0
      WD00=0.0
      ESIGN=ARREDC(INET)-ARREDC(INET-1)
      DO I=1,INET-1
         EDC=ARREDC(I)
         CUR=ARRCUR(I)
         EDC2=ARREDC(I+1)
         CUR2=ARRCUR(I+1)
         IF(ESIGN.LT.0) THEN
           IF(CUR2.LE.EX11.AND.EX11.LT.CUR) THEN
             EN11=EDC2+(EDC-EDC2)*(EX11-CUR2)/(CUR-CUR2)
             BE11=EN11
             WD11=0.0
           ELSEIF(CUR2.LE.EX12.AND.EX12.LT.CUR) THEN
             EN12=EDC2+(EDC-EDC2)*(EX12-CUR2)/(CUR-CUR2)
             BE12=EN12
             WD12=0.0
           ELSEIF(CUR2.LE.EX32.AND.EX32.LT.CUR) THEN
             EN32=EDC2+(EDC-EDC2)*(EX32-CUR2)/(CUR-CUR2)
             BE32=EN32
             WD32=0.0
           ELSEIF(CUR2.LE.EX13.AND.EX13.LT.CUR) THEN
             EN13=EDC2+(EDC-EDC2)*(EX13-CUR2)/(CUR-CUR2)
             BE13=EN13
             WD13=0.0
           ELSEIF(CUR2.LE.EX23.AND.EX23.LT.CUR) THEN
             EN23=EDC2+(EDC-EDC2)*(EX23-CUR2)/(CUR-CUR2)
             BE23=EN23
             WD23=0.0
           ELSEIF(CUR2.LE.EX43.AND.EX43.LT.CUR) THEN
             EN43=EDC2+(EDC-EDC2)*(EX43-CUR2)/(CUR-CUR2)
             BE43=EN43
             WD43=0.0
           ELSEIF(CUR2.LE.EX14.AND.EX14.LT.CUR) THEN
             EN14=EDC2+(EDC-EDC2)*(EX14-CUR2)/(CUR-CUR2)
             BE14=EN14
             WD14=0.0
           ELSEIF(CUR2.LE.EX34.AND.EX34.LT.CUR) THEN
             EN34=EDC2+(EDC-EDC2)*(EX34-CUR2)/(CUR-CUR2)
             BE34=EN34
             WD34=0.0
           ELSEIF(CUR2.LE.EX54.AND.EX54.LT.CUR) THEN
             EN54=EDC2+(EDC-EDC2)*(EX54-CUR2)/(CUR-CUR2)
             BE54=EN54
             WD54=0.0
           ENDIF
         ELSE
           IF(CUR.LE.EX11.AND.EX11.LT.CUR2) THEN
             EN11=EDC+(EDC-EDC2)*(EX11-CUR)/(CUR-CUR2)
             BE11=EN11
             WD11=0.0
           ELSEIF(CUR.LE.EX12.AND.EX12.LT.CUR2) THEN
             EN12=EDC+(EDC-EDC2)*(EX12-CUR)/(CUR-CUR2)
             BE12=EN12
             WD12=0.0
           ELSEIF(CUR.LE.EX32.AND.EX32.LT.CUR2) THEN
             EN32=EDC+(EDC-EDC2)*(EX32-CUR)/(CUR-CUR2)
             BE32=EN32
             WD32=0.0
           ELSEIF(CUR.LE.EX13.AND.EX13.LT.CUR2) THEN
             EN13=EDC+(EDC-EDC2)*(EX13-CUR)/(CUR-CUR2)
             BE13=EN13
             WD13=0.0
           ELSEIF(CUR.LE.EX23.AND.EX23.LT.CUR2) THEN
             EN23=EDC+(EDC-EDC2)*(EX23-CUR)/(CUR-CUR2)
             BE23=EN23
             WD23=0.0
           ELSEIF(CUR.LE.EX43.AND.EX43.LT.CUR2) THEN
             EN43=EDC+(EDC-EDC2)*(EX43-CUR)/(CUR-CUR2)
             BE43=EN43
             WD43=0.0
           ELSEIF(CUR.LE.EX14.AND.EX14.LT.CUR2) THEN
             EN14=EDC+(EDC-EDC2)*(EX14-CUR)/(CUR-CUR2)
             BE14=EN14
             WD14=0.0
           ELSEIF(CUR.LE.EX34.AND.EX34.LT.CUR2) THEN
             EN34=EDC+(EDC-EDC2)*(EX34-CUR)/(CUR-CUR2)
             BE34=EN34
             WD34=0.0
           ELSEIF(CUR.LE.EX54.AND.EX54.LT.CUR2) THEN
             EN54=EDC+(EDC-EDC2)*(EX54-CUR)/(CUR-CUR2)
             BE54=EN54
             WD54=0.0
           ENDIF
        ENDIF
      ENDDO
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     initial assigned zeros for the number of data for each step.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ZEROS(I00,I11,I21,I12,I32,I13,I23,I43,I14,I34,I54,
     *                                     I15,I25,I35,I45,I65)
      I00=0
      I11=0
      I21=0
      I12=0
      I32=0
      I13=0
      I23=0
      I43=0
      I14=0
      I34=0
      I54=0
      I15=0
      I25=0
      I35=0
      I45=0
      I65=0
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     find maximum, minimum, and its width for each step.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE MAXMIN(III,ZDC,EDCBEG,EDCEND,EDCWID,EDCINC,ROUND)
      DIMENSION ZDC(99)
      ZDCMIN=+1.0E+10
      ZDCMAX=-1.0E+10
      DO I=1,III
         IF(ZDC(I).LE.ZDCMIN) ZDCMIN=ZDC(I)
         IF(ZDC(I).GE.ZDCMAX) ZDCMAX=ZDC(I)
      ENDDO
      EDCBEG=ZDCMIN-EDCINC*0.05
      EDCEND=ZDCMAX+EDCINC*0.05
      EDCWID=EDCEND-EDCBEG
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     calculate evaluated sine and cosine functions.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ECSESS(PI)
      COMMON /EVALU/ ECS(-8192:8192),ESS(-8192:8192)
      DO I=0,2048
        ESS(I)=SIN(2.0*PI*I/2048.)
        ECS(I)=COS(2.0*PI*I/2048.)
      ENDDO
      DO I=0,2048
        ESS(I+2048)=ESS(I)
        ESS(I+4096)=ESS(I)
        ESS(I+6144)=ESS(I)
        ESS(I-2048)=ESS(I)
        ESS(I-4096)=ESS(I)
        ESS(I-6144)=ESS(I)
        ESS(I-8192)=ESS(I)
        ECS(I+2048)=ECS(I)
        ECS(I+4096)=ECS(I)
        ECS(I+6144)=ECS(I)
        ECS(I-2048)=ECS(I)
        ECS(I-4096)=ECS(I)
        ECS(I-6144)=ECS(I)
        ECS(I-8192)=ECS(I)
      ENDDO
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Initial phase configuration.
C     If NCHRG=1, use random/previous
C     If NCHRG=2, use random all
C     If NCHRG=3, use parallel/previous
C     If NCHRG=4, use parallel/all.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INICOND(ICURVE,NCHRG,LSEED,PHI)
      PARAMETER(MSIZ=4096,PI=3.141592653589793238462643 )
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /INCND/ PHI0(MSIZ)
      DIMENSION PHI(MSIZ),RNDD(MSIZ)
      PI2=2.0*PI
      PNV=1.0/PI2
      IF(ICURVE.EQ.1) THEN
         IF(NCHRG.EQ.3 .OR. NCHRG.EQ.4) THEN
            DO I=1,N
               PHI(I)=0.0
            ENDDO
         ELSE              ! IF(NCHRG.EQ.1 .OR. NCHRG.EQ.2) THEN
            CALL RANINI(LSEED)
            CALL RANUNI(LSEED,N,RNDD)
            DO I=1,N
               PHI(I)=2.0*(RNDD(I)-0.50)*PI
            ENDDO
         ENDIF
      ELSE
         IF(NCHRG.EQ.4) THEN
            DO I=1,N
               PHI(I)=0.0
            ENDDO
         ELSEIF(NCHRG.GE.2) THEN
            CALL RANINI(LSEED)
            CALL RANUNI(LSEED,N,RNDD)
            DO I=1,N
               PHI(I)=2.0*(RNDD(I)-0.50)*PI
            ENDDO
         ELSE
            DO I=1,N
               PHI(I)=PHI0(I)
            ENDDO
C11         READ(11,7510) PRECUR,PRETEM
C11         READ(11,7520) (PHI(I),I=1,N)
C11         CLOSE(UNIT=11)
         ENDIF
      ENDIF
C     DO IJK=1,LXLYM
C        PHI(IJK)=PHI(IJK)-PI2*IFIX(PHI(IJK)*PNV)
C     ENDDO
 7510 FORMAT(10X,F10.5,16X,F10.5)
 7520 FORMAT(7F11.6)
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     random noise at each node:
C     IPROB=+1, then use uniform  random noise.
C     IPROB=-1, then use Gaussian  random noise.
C     IPROB=otherwise, then use IPROB=+1
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE NOISDIS(TEMP,IPROB,SQRT3)
      PARAMETER(MSIZ=4096)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /NOIS1/ ISEED
      COMMON /NOIS2/ TNOISE(MSIZ)
      DIMENSION RNDD(MSIZ)
      IF(IPROB.EQ.-1) THEN	! Gaussian dist.
         CALL RANGAS(ISEED,N,RNDD)
         DO I=1,N
            TNOISE(I)=RNDD(I)*TEMP
         ENDDO
      ELSE			! regular uniform dist.
         CALL RANUNI(ISEED,N,RNDD)
         DO I=1,N
            TNOISE(I)=2.0*(RNDD(I)-0.50)*TEMP*SQRT3
        ENDDO
      ENDIF
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     write the time-dependent quantity and its Fourier spectrum:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CONVVTFR(NDATA)
      PARAMETER(MTIM=30010,MSIZ=4096)
      COMMON /AVERG/ CURREN(MTIM),CMIDL(MTIM),CMID0(MTIM),CMIDR(MTIM)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /RAND1/ IHAND,IBETA
      COMMON /RAND2/ H0,DH,BETA(MSIZ)
      COMMON /COUNT/ KARRAY(MTIM)
      COMMON /FLAGS/ IPRT,ISAAD,MJCT,NFFT,IPICK,NPROB
      DIMENSION  SPT(MTIM),SPAVG(MTIM)
      IF(IPRT.GT.99) GOTO 500
      WRITE(15,760)
      WRITE(15,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
      WRITE(15,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
      WRITE(15,785)
      DO I=1,NDATA
         SPT(I)=DELT*KARRAY(I)
         WRITE(15,860) SPT(I),CURREN(I),CMIDL(I),CMID0(I),CMIDR(I)
         IF(MJCT.GT.0) THEN
            SPAVG(I)=CURREN(I)
         ELSEIF(MJCT.EQ.-1) THEN
            SPAVG(I)=CMIDL(I)
         ELSEIF(MJCT.EQ.-2) THEN
            SPAVG(I)=CMID0(I)
         ELSEIF(MJCT.EQ.-3) THEN
            SPAVG(I)=CMIDR(I)
         ENDIF
      ENDDO
  500 CONTINUE
      CALL FFTPACK(NDATA,SPT,SPAVG)
  760 FORMAT(' LX LY LZ   *CCX  *CCY  *CCZ   BX  BY  BZ   *DELT')
  770 FORMAT(1X,I2,I3,I3,1X,3F6.2,1X,3I4,1X,F7.3)
  780 FORMAT('  *Edc*  Temp  *Ho*  *dH* IHAND IBETA  *Eac*  Frequ'/
     *          F7.3, F6.2, F6.2, F6.2,  I4,   I6,   F9.3,  F7.3   )
  785 FORMAT('  ****Time****   ***Current*** ',
     *         '**Midleft***Middle***Midright')
  860 FORMAT(1X,1PE13.6,1X,1PE15.6,1P3E10.2)
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Fourier transform and power spectrum:
C     (using fftpack from research.att.com)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE FFTPACK(NDATA,TIME,V)
C     FFTCI,F2TCF,F2TCB (IMSL);  REAL WORK(4*N+15), REAL CPY(2*N)  
      PARAMETER(MTIM=30010,MSIZ=4096,PI=3.141592653589793238462643)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /RAND1/ IHAND,IBETA
      COMMON /RAND2/ H0,DH,BETA(MSIZ)
      COMMON /FLAGS/ IPRT,ISAAD,MJCT,NFFT,IPICK,NPROB
      COMPLEX   CO(2*MTIM),CI(2*MTIM)
      DIMENSION TIME(MTIM),V(MTIM),FREQU(MTIM),SW(MTIM),WORK(4*MTIM+15)
      NTAU=NDATA
      NTAU1=NTAU
      NTAU2=ALOG(FLOAT(NDATA))/ALOG(2.)+0.001
      IF(NTAU2.GT.17) NTAU2=17
      NTAU2=2**NTAU2
      IF(NFFT.GT.0) THEN	! using FFT.
         NTAU=NTAU2
      ELSE			! using regular FT.
         NTAU=NTAU1
      ENDIF
      DT=(TIME(NTAU)-TIME(1))/(NTAU-1.0)
      DELF=1.0/(DT*NTAU)
      WRITE(6,*) 'delta_time =', DT
      WRITE(6,*) '# of data, delta_frequency =', NTAU,DELF
      DO ITAU=1,NTAU
         CI(ITAU)=V(ITAU)
         CO(ITAU)=V(ITAU)
      ENDDO
CIMSL CALL FFTCI(NTAU,WORK); CALL F2TCF(NTAU,CI,CO,WORK,CPY)
      CALL CFFTI(NTAU,WORK)
      CALL CFFTF(NTAU,CO,WORK)
      FACT0=1.0/(NTAU*NTAU)
      NSP=NTAU
      IF(NTAU.GT.4000) NSP=3000
C
C     power spectrum:
C
      DO I=1,NSP
         FREQU(I)=DELF*(I-1)
         SW(I)=ABS(CO(I))
         SW(I)=SW(I)*SW(I)*FACT0
         IF(SW(I).EQ.0.0) SW(I)=1.0E-20
      ENDDO
      SMAL=DELF*0.5
      IF(ISAAD.LE.0) THEN
        WRITE(28,760)
        WRITE(28,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
        WRITE(28,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
        WRITE(28,793)  NTAU
        WRITE(28,795)  (FREQU(I),SW(I),I=1,NSP)
      ENDIF
      IF(ISAAD.GE.0) THEN
        DOU=FREQE/2.0
        D1=0.0
        D2=0.0
        D3=0.0
        TRI=FREQE/3.0
        T1=0.0
        T2=0.0
        T3=0.0
        QUA=FREQE/4.0
        Q1=0.0
        Q2=0.0
        Q3=0.0
        STD=FREQE
        S1=0.0
        S2=0.0
        S3=0.0
        S4=0.0
        DO I=1,NSP
           FQ=FREQU(I)
           IF(FQ.GT.  STD-SMAL.AND.FQ.LE.  STD+SMAL) S1=S1+SW(I)
           IF(FQ.GT.2*STD-SMAL.AND.FQ.LE.2*STD+SMAL) S2=S2+SW(I)
           IF(FQ.GT.3*STD-SMAL.AND.FQ.LE.3*STD+SMAL) S3=S3+SW(I)
           IF(FQ.GT.4*STD-SMAL.AND.FQ.LE.4*STD+SMAL) S4=S4+SW(I)
           IF(FQ.GT.  DOU-SMAL.AND.FQ.LE.  DOU+SMAL) D1=D1+SW(I)
           IF(FQ.GT.3*DOU-SMAL.AND.FQ.LE.3*DOU+SMAL) D2=D2+SW(I)
           IF(FQ.GT.5*DOU-SMAL.AND.FQ.LE.5*DOU+SMAL) D3=D3+SW(I)
           IF(FQ.GT.  TRI-SMAL.AND.FQ.LE.  TRI+SMAL) T1=T1+SW(I)
           IF(FQ.GT.2*TRI-SMAL.AND.FQ.LE.2*TRI+SMAL) T2=T2+SW(I)
           IF(FQ.GT.4*TRI-SMAL.AND.FQ.LE.4*TRI+SMAL) T3=T3+SW(I)
           IF(FQ.GT.  QUA-SMAL.AND.FQ.LE.  QUA+SMAL) Q1=Q1+SW(I)
           IF(FQ.GT.3*QUA-SMAL.AND.FQ.LE.3*QUA+SMAL) Q2=Q2+SW(I)
           IF(FQ.GT.5*QUA-SMAL.AND.FQ.LE.5*QUA+SMAL) Q3=Q3+SW(I)
        ENDDO
        WRITE(29,710) EDC0,S1,S2,S3,S4
        WRITE(30,710) EDC0,D1,D2,D3
        WRITE(31,710) EDC0,T1,T2,T3
        WRITE(32,710) EDC0,Q1,Q2,Q3
      ENDIF
  710 FORMAT(1P10E12.4)
  760 FORMAT(' LX LY LZ   *CCX  *CCY  *CCZ   BX  BY  BZ   *DELT')
  770 FORMAT(1X,I2,I3,I3,1X,3F6.2,1X,3I4,1X,F7.3)
  780 FORMAT('  *Edc*  Temp  *Ho*  *dH* IHAND IBETA  *Eac*  Frequ'/
     *          F7.3, F6.2, F6.2, F6.2,  I4,   I6,   F9.3,  F7.3   )
  793 FORMAT(' **Frequency**  ****S(w)*****  ( N  =',I6,' )') 
  795 FORMAT(1P2E14.6)
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     read data file:

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE READDATA(IH,JH,KH,PI2,PNV,SQRT3,NNNEEE,NTEMP)
      PARAMETER(ZERO=0.0,PI=3.141592653589793238462643 )
      PARAMETER(MTIM=30010,MSIZ=4096)
      INTEGER        BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /INTGR/ BFL,BFR,BBL,BBR,TFL,TFR,TBL,TBR
      COMMON /DIMN1/ EE(MSIZ),EDCPRT(MSIZ)
      COMMON /DIMN2/ HH(MSIZ),GRA2(MSIZ),PHIDOT(MSIZ)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      COMMON /RAND1/ IHAND,IBETA
      COMMON /RAND2/ H0,DH,BETA(MSIZ)
      COMMON /NOIS1/ ISEED
      COMMON /NOIS2/ TNOISE(MSIZ)
      COMMON /FLAGS/ IPRT,ISAAD,MJCT,NFFT,IPICK,NPROB
      COMMON /DATA1/ TEM00,TEMINC,TEM99,IPROB,ISCAN
      COMMON /DATA2/ ISEED0,LSEED0,MSEED0,NSEED0,NREALS
      COMMON /DATA3/ ISTEP,ISTAVG,NSTEP,ISTOFF,IBACK
      COMMON /DATA4/ EDC00,EDC33,EDC66,EDC99,EINC1,EINC2,EINC3
      COMMON /DATA5/ IFLAG,NCHRG,IEVA
      COMMON /DATA6/ EACINC,FREINC,ACPT,ROUND
      COMMON /EVALU/ ECS(-8192:8192),ESS(-8192:8192)
      PI2=2.0*PI
      PNV=1.0/PI2
      SQRT3=SQRT(3.0)
      READ(51,*)
      READ(51,*)
      READ(51,*)  LX,LY,LZ
      READ(51,*)
      READ(51,*)  IBCX,IBCY,IBCZ
      IF(LZ.EQ.2.AND.IBCZ.LT.0) THEN
         WRITE(6,*) 'Lz=2 and pbc on z=axis: meaningless!'
         STOP
      ENDIF
      READ(51,*)
      READ(51,*)  CCX,CCY,CCZ    ! anisotropic coupling constant.
      READ(51,*)
      READ(51,*)  H0,DH
      READ(51,*)
      READ(51,*)  DELT
      READ(51,*)
C
C     H0     :  absolute value of strength h.
C     DH     :  standard deviation of strength H0. 
C     Temperature = TEM00, ....., TEM99
C     ISCAN < 0 : scan currents,
C           > 0 : scan temperature.
C     NREALS    : the number of realizations.
C
      READ(51,*)  TEM00,TEMINC,TEM99
      READ(51,*)
      READ(51,*)  ISCAN,NREALS
      READ(51,*)
      READ(51,*)  IPROB
      READ(51,*)
      READ(51,*)  ISEED0,LSEED0
      READ(51,*)
C
C     IHAND: 0=ordered strength (fixed Ho for all sites).
C            1=uniformly distributed ranom strength from Ho-dH to Ho+dH.
C            21,22,23=same as 1 but correlated on x,y,z (columnar defects).
C            31,32,33=same as 1 but correlated on xy,yz,zx (planar defects).
C
      READ(51,*)  MSEED0,IHAND
      READ(51,*)
C
C     IBETA=0:  no positional disorder of impurities
C          >0:  positional disorder of impuritires
C          <0:  the sameas above.
C
C     average from ISTEP*ISTAVG to NSTEP-ISTEP with increment ISTEP.
C
C
      READ(51,*)  NSEED0,IBETA         !random seed for beta
      READ(51,*)
      READ(51,*)  ISTEP,ISTAVG,NSTEP
      READ(51,*)
      READ(51,*)  ISTOFF
      READ(51,*)
      READ(51,*)  Edc00,EINC1
      READ(51,*)
      READ(51,*)  Edc33,EINC2
      READ(51,*)
      READ(51,*)  Edc66,EINC3
      READ(51,*)
      READ(51,*)  Edc99,IBACK
      READ(51,*)
      READ(51,*)  EAC,FREQE,IPRT
      IF(IPRT.GT.0 .AND. IPRT.LT.100) THEN
        READ(51,*)
        READ(51,*)  (EdcPRT(I),I=1,IPRT)
        ELSE
        READ(51,*)
        READ(51,*)
      ENDIF
C
C     IFLAG>0:  scan many realizations.
C          <0:  scan many ac electric fields.
C          =0:  scan many ac frequencies.
C     EACINC: ac electric field increment for IFLAG<0.
C     FREINC: ac frequency increment for IFLAG<0.
C     (CAUTION: use the same quenched disorder. disable changing seed)
C
C
      READ(51,*)
      READ(51,*)  IFLAG,EACINC,FREINC
C
C     IPICK > 0 : regular Edc and current.
C     IPICK = 0 : top-plane Edc / bottom_plane current.
C     IPICK =-1 : top-plane Edc / transverse bottom_plane current.
C     IPICK <-1 : regular Edc   / transverse current (Hall geometry).
C
      READ(51,*)
      READ(51,*)  IPICK
C
      IBREAK = ISTOFF*ISTEP
      INTERV = ISTAVG*ISTEP
C
C     ISTOFF < 0 : normal procedure 
C            > 0 : turn off Edc abruptly at ISTEP*ISToff
C                  to get remanent voltage decay.
C
C     MJCT > 0  choose an average voltage  
C          =-1  pick a current at a mid-left junction,
C          =-2  pick a current at a mid junction,
C          =-3  pick a current at a mid-right junction,
C     ISAAD< 0  get S(w),
C          > 0  get S(Idc),
C          = 0  get both S(w) and S(Idc). 
C     NFFT < 0  regular FT,
C          > 0  effective FFT. 
C
      READ(51,*)
      READ(51,*)  MJCT
      READ(51,*)
      READ(51,*)  ISAAD, NFFT
      READ(51,*)
      READ(51,*)  NCHRG 
      READ(51,*)
      READ(51,*)  IEVA
C
C     ACPT : accept for a step by |Observed-Expected|<ACPT. (def=0.003).
C     ROUND: extend step widths at edges. (def=0.02 times Einc1).
C
      READ(51,*)
      READ(51,*)  ACPT,ROUND
C
C     NCHRG = 1 random/previous initial condition
C           = 2 random/all
C           = 3 parallel/previous
C           = 4 parallel/all 
C 
      N  = LX*LY*LZ
      LXLYLZ= N
      LXLY  = LX*LY
      LXLZ  = LX*LZ
      LYLZ  = LY*LZ
      IH    = (LX+1)/2
      JH    = (LY+1)/2
      KH    = (LZ+1)/2
      BFL=1
      BFR=LX
      BBL=LXLY-LX+1
      BBR=LXLY
      TFL=LXLYLZ-LXLY+1
      TFR=LXLYLZ-LXLY+LX
      TBL=LXLYLZ-LX+1
      TBR=LXLYLZ
      DO I=1,N
         HH(I)=H0
      ENDDO
C
C     make a table of sine and cosine function for IEVA<0.
C
      IF(IEVA.LT.0) CALL ECSESS(PI)  !Make tables for evaluted functions:
      TDEL=1.0/DELT
      WE=2.0*PI*FREQE
      NALL=(NSTEP-ISTAVG*ISTEP)/ISTEP
      IF(MTIM.LT.NALL) THEN
         WRITE(6,*) 'Small space for time evolution.'
         WRITE(6,*) 'Adjust dimensions!'
         STOP
      ENDIF
C
      EDC9=EDC99
      ESIGN1=EINC1*1.0001/ABS(EINC1)
      EE(1)=EDC00
      IF(ISCAN.GE.0) THEN
         NNNEEE=1
         GOTO 1851
      ENDIF
      DO I1=2,1000
         EE(I1)=EE(I1-1)+EINC1
         IF(ESIGN1.GE.0.0) THEN
           IF(EE(I1).GT.EDC33) GOTO 165
         ELSE
           IF(EE(I1).LT.EDC33) GOTO 165
         ENDIF
      ENDDO
  165 NNEE=I1-1
      ESIGN2=EINC2/ABS(EINC2)
      EE(NNEE)=EDC33
      DO I2=NNEE+1,1000
         EE(I2)=EE(I2-1)+EINC2
         IF(ESIGN2.GE.0.0) THEN
           IF(EE(I2).GT.EDC66) GOTO 175
         ELSE
           IF(EE(I2).LT.EDC66) GOTO 175
         ENDIF
      ENDDO
  175 NNEE=I2-1
      ESIGN3=EINC3/ABS(EINC3)
      EE(NNEE)=EDC66
      DO I3=NNEE+1,1000
         EE(I3)=EE(I3-1)+EINC3
         IF(ESIGN3.GE.0.0) THEN
           IF(EE(I3).GT.EDC99) GOTO 185
         ELSE
           IF(EE(I3).LT.EDC99) GOTO 185
         ENDIF
      ENDDO
  185 NNNEEE=I3-1-1
 1851 WRITE(6,*) 'The number of Edc =',NNNEEE
      WRITE(6,'(1X,15F5.2)') (EE(I3),I3=1,NNNEEE)
  195 TEMP9=TEM99
      NTEMP=( ABS(TEMP9-TEMP0)+0.0000001 )/ABS(TEMINC)
      NTEMP=NTEMP+1
      IF(NTEMP.EQ.1) NTEMP=2
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     write output-1:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE WRITEOUT1(IFLAG,IHAND,H0,DH,IDETA)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      IF(IFLAG.LE.0) THEN   !assuming ISCAN<0 (electric field scan).
         WRITE(20,760)
         WRITE(20,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
         WRITE(20,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
         WRITE(21,760)
         WRITE(21,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
         WRITE(21,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
         WRITE(22,760)
         WRITE(22,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
         WRITE(22,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
         IF(IFLAG.LT.0) THEN
            WRITE(20,8200)
            WRITE(21,8210)
            WRITE(22,8220)
         ELSE
            WRITE(20,8201)
            WRITE(21,8211)
            WRITE(22,8221)
         ENDIF
      ENDIF
  760 FORMAT(' LX LY LZ   *CCX  *CCY  *CCZ   BX  BY  BZ   *DELT')
  770 FORMAT(1X,I2,I3,I3,1X,3F6.2,1X,3I4,1X,F7.3)
  780 FORMAT('  *Edc*  Temp  *Ho*  *dH* IHAND IBETA  *Eac*  Frequ'/
     *          F7.3, F6.2, F6.2, F6.2,  I4,   I6,   F9.3,  F7.3   )
 8200 FORMAT(' **Eac*  **Step 0/0****  **Step 1/1****  **Step 1/2****',
     *                              '  **Step 1/3****  **Step 2/3****')
 8210 FORMAT(' **Eac*  **Step 1/4****  **Step 3/4****  **Step 3/2****',
     *                              '  **Step 4/3****  **Step 5/4****')
 8220 FORMAT('  **Eac*  **0/0**  **1/1**  **1/2**  **1/3**  **2/3**',
     *               '  **1/4**  **3/4**  **3/2**  **4/3**  **5/4**')
 8201 FORMAT(' *Freq*  **Step 0/0****  **Step 1/1****  **Step 1/2****',
     *       '  **Step 1/3****  **Step 2/3****')
 8211 FORMAT(' *Freq*  **Step 1/4****  **Step 3/4****  **Step 3/2****',
     *                              '  **Step 4/3****  **Step 5/4****')
 8221 FORMAT('  *Freq*  **0/0**  **1/1**  **1/2**  **1/3**  **2/3**',
     *               '  **1/4**  **3/4**  **3/2**  **4/3**  **5/4**')
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     write output-2:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE WRITEOUT2(ISCAN,ISAAD,IHAND,H0,DH,IBETA)
      COMMON /PARA1/ LX,LY,LZ,IBCX,IBCY,IBCZ,CCX,CCY,CCZ
      COMMON /PARA2/ LXLYLZ,N,LXLY,LXLZ,LYLZ
      COMMON /PARA3/ DELT,EIN,TEMP0,EDC0,EAC,FREQE,TDEL,WE
      WRITE(9,760)
      WRITE(9,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
      WRITE(9,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
      IF(ISCAN.LT.0) THEN 
         WRITE(9,790)
      ELSE
         WRITE(9,791)
      ENDIF
      IF(ISAAD.GE.0) THEN
        WRITE(29,760)
        WRITE(29,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
        WRITE(29,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
        WRITE(29,7029)
        WRITE(30,760)
        WRITE(30,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
        WRITE(30,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
        WRITE(30,7030)
        WRITE(31,760)
        WRITE(31,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
        WRITE(31,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
        WRITE(31,7031)
        WRITE(32,760)
        WRITE(32,770)  LX,LY,LZ,CCX,CCY,CCZ,IBCX,IBCY,IBCZ,DELT
        WRITE(32,780)  EDC0,TEMP0,H0,DH,IHAND,IBETA,EAC,FREQE
        WRITE(32,7032) 
      ENDIF
  760 FORMAT(' LX LY LZ   *CCX  *CCY  *CCZ   BX  BY  BZ   *DELT')
  770 FORMAT(1X,I2,I3,I3,1X,3F6.2,1X,3I4,1X,F7.3)
  780 FORMAT('  *Edc*  Temp  *Ho*  *dH* IHAND IBETA  *Eac*  Frequ'/
     *          F7.3, F6.2, F6.2, F6.2,  I4,   I6,   F9.3,  F7.3   )
  790 FORMAT('  ***Edc*** **Current** **Midleft***Middle***Midright**')
  791 FORMAT('  ***Temp** **Current** **Midleft***Middle***Midright**')
  795 FORMAT(F10.5,F12.6,3F10.4)
 7029 FORMAT(2X,'Idc',9X,'m=1/1',7X,'m=2/1',7X,'m=3/1',7X,'m=4/1')
 7030 FORMAT(2X,'Idc',9X,'m=1/2',7X,'m=3/2',7X,'m=5/2' )
 7031 FORMAT(2X,'Idc',9X,'m=1/3',7X,'m=2/3',7X,'m=4/3' )
 7032 FORMAT(2X,'Idc',9X,'m=1/4',7X,'m=3/4',7X,'m=5/4' )
 7530 FORMAT(' Current =',F10.5,',  Temperature =',F10.5)
      RETURN
      END
CC
CC    This is the end of the regular part.
CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     ranini.f for initialize random seed:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE RANINI(ISEED)
      IF(ISEED.GT.0) ISEED=-ISEED
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     ranuni.f  for uniform random numbers:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE RANUNI(ISEED,N,R)
      DIMENSION R(N)
      DO I=1,N
         R(I)=RAN2(ISEED)
      ENDDO
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     rangas.f for Gaussian random numbers:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE RANGAS(ISEED,N,R)
      DIMENSION R(N)
      DO I=1,N
         R(I)=GASDEV(ISEED)
      ENDDO
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     gasdev.f from numerical recipes:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION GASDEV(IDUM)
      DATA ISET/0/
      IF (ISET.EQ.0) THEN
1       V1=2.*RAN1(IDUM)-1.
        V2=2.*RAN1(IDUM)-1.
        R=V1**2+V2**2
        IF(R.GE.1.)GO TO 1
        FAC=SQRT(-2.*LOG(R)/R)
        GSET=V1*FAC
        GASDEV=V2*FAC
        ISET=1
      ELSE
        GASDEV=GSET
        ISET=0
      ENDIF
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     ran1.f from numerical recipes:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION RAN1(IDUM)
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1./M1)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
      DATA IFF /0/
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO J=1,97
                IX1=MOD(IA1*IX1+IC1,M1)
                IX2=MOD(IA2*IX2+IC2,M2)
                R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
        ENDDO
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
      IF(J.GT.97.OR.J.LT.1)PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     ran2.f from numerical recipes:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION RAN2(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1./M)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO J=1,97
                IDUM=MOD(IA*IDUM+IC,M)
                IR(J)=IDUM
        ENDDO
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     cffti.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CFFTI (N,WSAVE)
      DIMENSION       WSAVE(1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     cfftf.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CFFTF (N,C,WSAVE)
      DIMENSION       C(1)       ,WSAVE(1)
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     cffti1.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CFFTI1 (N,WA,IFAC)
      DIMENSION       WA(1)      ,IFAC(1)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      IFAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         IFAC(IB+2) = IFAC(IB+1)
  106 CONTINUE
      IFAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      IFAC(1) = N
      IFAC(2) = NF
      TPI = 6.28318530717959
      ARGH = TPI/FLOAT(N)
      I = 2
      L1 = 1
      DO 110 K1=1,NF
         IP = IFAC(K1+2)
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IDOT = IDO+IDO+2
         IPM = IP-1
         DO 109 J=1,IPM
            I1 = I
            WA(I-1) = 1.
            WA(I) = 0.
            LD = LD+L1
            FI = 0.
            ARGLD = FLOAT(LD)*ARGH
            DO 108 II=4,IDOT,2
               I = I+2
               FI = FI+1.
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IF (IP .LE. 5) GO TO 109
            WA(I1-1) = WA(I-1)
            WA(I1) = WA(I)
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     cfftf1.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CFFTF1 (N,C,CH,WA,IFAC)
      DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(1)
      NF = IFAC(2)
      NA = 0
      L1 = 1
      IW = 1
      DO 116 K1=1,NF
         IP = IFAC(K1+2)
         L2 = IP*L1
         IDO = N/L2
         IDOT = IDO+IDO
         IDL1 = IDOT*L1
         IF (IP .NE. 4) GO TO 103
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IF (NA .NE. 0) GO TO 101
         CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
         GO TO 102
  101    CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
         GO TO 115
  103    IF (IP .NE. 2) GO TO 106
         IF (NA .NE. 0) GO TO 104
         CALL PASSF2 (IDOT,L1,C,CH,WA(IW))
         GO TO 105
  104    CALL PASSF2 (IDOT,L1,CH,C,WA(IW))
  105    NA = 1-NA
         GO TO 115
  106    IF (IP .NE. 3) GO TO 109
         IX2 = IW+IDOT
         IF (NA .NE. 0) GO TO 107
         CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
         GO TO 108
  107    CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
         GO TO 115
  109    IF (IP .NE. 5) GO TO 112
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IX4 = IX3+IDOT
         IF (NA .NE. 0) GO TO 110
         CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 111
  110    CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
         GO TO 115
  112    IF (NA .NE. 0) GO TO 113
         CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
         GO TO 114
  113    CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    IF (NAC .NE. 0) NA = 1-NA
  115    L1 = L2
         IW = IW+(IP-1)*IDOT
  116 CONTINUE
      IF (NA .EQ. 0) RETURN
      N2 = N+N
      DO 117 I=1,N2
         C(I) = CH(I)
  117 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     passf.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PASSF (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
     2                CH2(IDL1,IP)
      IDOT = IDO/2
      NT = IP*IDL1
      IPP2 = IP+2
      IPPH = (IP+1)/2
      IDP = IP*IDO
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 J=2,IPPH
         JC = IPP2-J
         DO 102 K=1,L1
            DO 101 I=1,IDO
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       CONTINUE
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 112
  106 DO 109 J=2,IPPH
         JC = IPP2-J
         DO 108 I=1,IDO
            DO 107 K=1,L1
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       CONTINUE
  108    CONTINUE
  109 CONTINUE
      DO 111 I=1,IDO
         DO 110 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  110    CONTINUE
  111 CONTINUE
  112 IDL = 2-IDO
      INC = 0
      DO 116 L=2,IPPH
         LC = IPP2-L
         IDL = IDL+IDO
         DO 113 IK=1,IDL1
            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
            C2(IK,LC) = -WA(IDL)*CH2(IK,IP)
  113    CONTINUE
         IDLJ = IDL
         INC = INC+IDO
         DO 115 J=3,IPPH
            JC = IPP2-J
            IDLJ = IDLJ+INC
            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
            WAR = WA(IDLJ-1)
            WAI = WA(IDLJ)
            DO 114 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC)
  114       CONTINUE
  115    CONTINUE
  116 CONTINUE
      DO 118 J=2,IPPH
         DO 117 IK=1,IDL1
            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  117    CONTINUE
  118 CONTINUE
      DO 120 J=2,IPPH
         JC = IPP2-J
         DO 119 IK=2,IDL1,2
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  119    CONTINUE
  120 CONTINUE
      NAC = 1
      IF (IDO .EQ. 2) RETURN
      NAC = 0
      DO 121 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  121 CONTINUE
      DO 123 J=2,IP
         DO 122 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  122    CONTINUE
  123 CONTINUE
      IF (IDOT .GT. L1) GO TO 127
      IDIJ = 0
      DO 126 J=2,IP
         IDIJ = IDIJ+2
         DO 125 I=4,IDO,2
            IDIJ = IDIJ+2
            DO 124 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  124       CONTINUE
  125    CONTINUE
  126 CONTINUE
      RETURN
  127 IDJ = 2-IDO
      DO 130 J=2,IP
         IDJ = IDJ+IDO
         DO 129 K=1,L1
            IDIJ = IDJ
            DO 128 I=4,IDO,2
               IDIJ = IDIJ+2
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  128       CONTINUE
  129    CONTINUE
  130 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     passf2.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PASSF2 (IDO,L1,CC,CH,WA1)
      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)           ,
     1                WA1(1)
      IF (IDO .GT. 2) GO TO 102
      DO 101 K=1,L1
         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            TI2 = CC(I,1,K)-CC(I,2,K)
            CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2
            CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2
  103    CONTINUE
  104 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     passf3.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PASSF3 (IDO,L1,CC,CH,WA1,WA2)
      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)           ,
     1                WA1(1)     ,WA2(1)
      DATA TAUR,TAUI /-.5,-.866025403784439/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TR2 = CC(1,2,K)+CC(1,3,K)
         CR2 = CC(1,1,K)+TAUR*TR2
         CH(1,K,1) = CC(1,1,K)+TR2
         TI2 = CC(2,2,K)+CC(2,3,K)
         CI2 = CC(2,1,K)+TAUR*TI2
         CH(2,K,1) = CC(2,1,K)+TI2
         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
         CH(1,K,2) = CR2-CI3
         CH(1,K,3) = CR2+CI3
         CH(2,K,2) = CI2+CR3
         CH(2,K,3) = CI2-CR3
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,2,K)+CC(I,3,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
  103    CONTINUE
  104 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     passf4.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PASSF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)           ,
     1                WA1(1)     ,WA2(1)     ,WA3(1)
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI1 = CC(2,1,K)-CC(2,3,K)
         TI2 = CC(2,1,K)+CC(2,3,K)
         TR4 = CC(2,2,K)-CC(2,4,K)
         TI3 = CC(2,2,K)+CC(2,4,K)
         TR1 = CC(1,1,K)-CC(1,3,K)
         TR2 = CC(1,1,K)+CC(1,3,K)
         TI4 = CC(1,4,K)-CC(1,2,K)
         TR3 = CC(1,2,K)+CC(1,4,K)
         CH(1,K,1) = TR2+TR3
         CH(1,K,3) = TR2-TR3
         CH(2,K,1) = TI2+TI3
         CH(2,K,3) = TI2-TI3
         CH(1,K,2) = TR1+TR4
         CH(1,K,4) = TR1-TR4
         CH(2,K,2) = TI1+TI4
         CH(2,K,4) = TI1-TI4
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI1 = CC(I,1,K)-CC(I,3,K)
            TI2 = CC(I,1,K)+CC(I,3,K)
            TI3 = CC(I,2,K)+CC(I,4,K)
            TR4 = CC(I,2,K)-CC(I,4,K)
            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
            TI4 = CC(I-1,4,K)-CC(I-1,2,K)
            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
            CH(I-1,K,1) = TR2+TR3
            CR3 = TR2-TR3
            CH(I,K,1) = TI2+TI3
            CI3 = TI2-TI3
            CR2 = TR1+TR4
            CR4 = TR1-TR4
            CI2 = TI1+TI4
            CI4 = TI1-TI4
            CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2
            CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2
            CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3
            CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3
            CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4
            CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4
  103    CONTINUE
  104 CONTINUE
      RETURN
      END
CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     passf5.f from fftpack in research.netlib.com:
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PASSF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)           ,
     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
      DATA TR11,TI11,TR12,TI12 /.309016994374947,-.951056516295154,
     1-.809016994374947,-.587785252292473/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI5 = CC(2,2,K)-CC(2,5,K)
         TI2 = CC(2,2,K)+CC(2,5,K)
         TI4 = CC(2,3,K)-CC(2,4,K)
         TI3 = CC(2,3,K)+CC(2,4,K)
         TR5 = CC(1,2,K)-CC(1,5,K)
         TR2 = CC(1,2,K)+CC(1,5,K)
         TR4 = CC(1,3,K)-CC(1,4,K)
         TR3 = CC(1,3,K)+CC(1,4,K)
         CH(1,K,1) = CC(1,1,K)+TR2+TR3
         CH(2,K,1) = CC(2,1,K)+TI2+TI3
         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
         CR5 = TI11*TR5+TI12*TR4
         CI5 = TI11*TI5+TI12*TI4
         CR4 = TI12*TR5-TI11*TR4
         CI4 = TI12*TI5-TI11*TI4
         CH(1,K,2) = CR2-CI5
         CH(1,K,5) = CR2+CI5
         CH(2,K,2) = CI2+CR5
         CH(2,K,3) = CI3+CR4
         CH(1,K,3) = CR3-CI4
         CH(1,K,4) = CR3+CI4
         CH(2,K,4) = CI3-CR4
         CH(2,K,5) = CI2-CR5
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI5 = CC(I,2,K)-CC(I,5,K)
            TI2 = CC(I,2,K)+CC(I,5,K)
            TI4 = CC(I,3,K)-CC(I,4,K)
            TI3 = CC(I,3,K)+CC(I,4,K)
            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
            CH(I,K,1) = CC(I,1,K)+TI2+TI3
            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
            CR5 = TI11*TR5+TI12*TR4
            CI5 = TI11*TI5+TI12*TI4
            CR4 = TI12*TR5-TI11*TR4
            CI4 = TI12*TI5-TI11*TI4
            DR3 = CR3-CI4
            DR4 = CR3+CI4
            DI3 = CI3+CR4
            DI4 = CI3-CR4
            DR5 = CR2+CI5
            DR2 = CR2-CI5
            DI5 = CI2-CR5
            DI2 = CI2+CR5
            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
            CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4
            CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4
            CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5
            CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5
  103    CONTINUE
  104 CONTINUE
      RETURN
      END
CC
CCC