C     PROGRAM FOR COMPUTING THE BAYES FACTOR FOR GBM AND GSGJ
C     Q1=LIKELIHOOD RATIO OF GBM AND GSGJ
C     Q2=LIKELIHOOD RATIO OF GSGJ AND GBM

      PROGRAM BAYSEFACTOR
      INCLUDE 'parbf1.f'
      INTEGER K,L1,M1,W1,L2,M2,W2,LB2
      REAL*8 PX0(N),PM1(NM1),PS1(NS1),PR1(NR1),Q1(N,NM1,NS1,NR1)
      REAL*8 PM2(NM2),PS2(NS2),PR2(NR2),PL2(NL2),Q2(N,NM2,NS2,NL2,NR2),
     &BF1,BF2
C     -------------------------------------------------
C     COMMON FOR Q1 AND Q2
c     ASSIGN INITIAL DISTRIBUTION OF X_0 ---------------
      DO 10 K=1,N
         PX0(K)=0.0d0
 10   CONTINUE
      PX0(I0)=1.0d0
C     ###################################################
C     ASSIGN INITIAL DISTRIBUTION FOR Q1
C     ASSIGN PRIOR FOR MU1, PM(NM1) ----------------------
      DO 20 L1=1,NM1
         PM1(L1)=1.0d0/dble(NM1)
 20   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR SIGMA1, PS(NS1) -------------------
      DO 30 M1=1,NS1
         PS1(M1)=1.0d0/dble(NS1)
 30   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR RHO, PR(NR)
      DO 40 W1=1,NR1
         PR1(W1)=1.0D0/DBLE(NR1)
 40   CONTINUE
c     COMPUTE INITIAL Q1(X_0,MU1,SG1,R01) ----
      DO 600 K=1,N
         DO 400 L1=1,NM1
            DO 200 M1=1,NS1
            DO 190 W1=1,NR1
                  Q1(K,L1,M1,W1)=PX0(K)*PM1(L1)*PS1(M1)*PR1(W1)
c                  if (Q1(K,L1,M1,W1) .gt. 0d0) then
c                     write(*,*) K,L1,M1,W1, Q1(K,L1,M1,W1)
c                     pause 
c                     end if
 190           CONTINUE
 200           CONTINUE
 400           CONTINUE
 600           CONTINUE
C     ------------------------------------------------
c      WRITE(*,*) "Q1 INITIAL", Q1
c      pause "initial Q1"
C     ########################################################
C     ASSIGN INITIAL DISTRIBUTION FOR Q2
C     ASSIGN PRIOR FOR MU, PM(NM) ----------------------
      DO 22 L2=1,NM2
         PM2(L2)=1.0d0/dble(NM2)
 22   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR SIGMA, PS(NS) -------------------
      DO 32 M2=1,NS2
         PS2(M2)=1.0d0/dble(NS2)
 32   CONTINUE
c     --------------------------------------------------
C     ASSIGN PRIOR FOR RHO, PR(NR)----------------------
      DO 42 W2=1,NR2
         PR2(W2)=1.0D0/DBLE(NR2)
 42   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR Lambda, PL(NL)-------------------
      DO 52 LB2=1,NL2
         PL2(LB2)=1.0D0/DBLE(NL2)
 52   CONTINUE
c     COMPUTE INITIAL Q2(X_0,MU2,SG2,LB2,RO2) ----
      DO 602 K=1,N
         DO 402 L2=1,NM2
           DO 302 LB2=1,NL2
            DO 202 M2=1,NS2
            DO 192 W2=1,NR2
                  Q2(K,L2,M2,LB2,W2)=PX0(K)*PM2(L2)*PS2(M2)*
     &PL2(LB2)*PR2(W2)
c                  if (Q2(k,l2,m2,lB2,w2) .gt. 0d0) then
c                     write(*,*) k,l2,m2,lb2,w2,Q2(k,l2,m2,lb2,w2)
c                     pause
c                     end if
 192           CONTINUE
 202        CONTINUE
 302     CONTINUE
 402  CONTINUE
 602  CONTINUE
C     ###########################################################
C     ------------------------------------------------
      OPEN(UNIT=32,FILE='mid.rlt',STATUS="UNKNOWN")
      CALL RATIO(Q1,Q2)
C     ------------------------------------------------------
C     COMPUTE THE FINAL BAYES FACTORS
      BF1=0D0
         DO 180 L1=1,NM1
            DO 182 M1=1,NS1
               DO 184 W1=1,NR1
               DO 186 K=1,N
                 BF1=BF1+Q1(K,L1,M1,W1)
 186            CONTINUE
 184            CONTINUE
 182         CONTINUE
 180      CONTINUE
      WRITE(32,*) "BF1=", BF1
      BF2=0D0
         DO 170 L2=1,NM2
            DO 172 M2=1,NS2
               DO 174 W2=1,NR2
               DO 176 K=1,N
                DO 178 LB2=1,NL2
                   BF2=BF2+Q2(K,L2,M2,LB2,W2)
 178             CONTINUE
 176            CONTINUE
 174            CONTINUE
 172         CONTINUE
 170      CONTINUE
      WRITE(32,*) "BF2=", BF2
      WRITE(32,*) "BF1*BF2=", BF1*BF2
      CLOSE(32)
      END 

      SUBROUTINE RATIO(Q1,Q2)
      EXTERNAL LD11,LD12,LD21,LD22, LD31, LD32
      INCLUDE 'parbf1.f'
      INTEGER I,K,L,L1,L2,M1,M2,NTI,J,W1,W2,LB2
      REAL*8 Q1(N,NM1,NS1,NR1),Q2(N,NM2,NS2,NL2,NR2),
     &A1(NM1,NS1,N),B1(NM1,NS1,N),A2(NM2,NS2,N),B2(NM2,NS2,N)
      REAL*8 MU1,SG1,XB,XC,LD11,LD21,LD31,X2,DT,DDT
      REAL*8 MU2,SG2,LBD2,LD12,LD22,LD32
      REAL*8 QI(N),Q2SM(N,NM2,NL2,NR2)
c     asign A1(NM1,NS1,N),B1(NM1,NS1,N) accordingly------------------
      DO 20 J=1,NM1
        MU1=MU1MIN+(DBLE(J)-1.0000D0)*MU1L
       DO 40 K=1,NS1
         SG1=SG1MIN+(DBLE(K)-1.0000D0)*SG1L
        DO 60 L=1,N
          X=XMIN-LL+(DBLE(L)-1.0000D0)*XL
c          write(*,*) J,MU1,K,SG1,L,X,XL
            IF (MU1 .GE. 0D0) THEN
          A1(J,K,L)=(SG1*SG1*X*X)/(2.0d0*XL*XL)+MU1*X/XL
          B1(J,K,L)=(SG1*SG1*X*X)/(2.0d0*XL*XL)
            ELSE
          A1(J,K,L)=(SG1*SG1*X*X)/(2.0d0*XL*XL)
          B1(J,K,L)=(SG1*SG1*X*X)/(2.0d0*XL*XL)-MU1*X/XL
            end if
c          write(*,*) A1(J,K,L),B1(J,K,L)
c          pause
 60      CONTINUE
 40   CONTINUE
 20   CONTINUE
c     assign A2(NM2,NS2,N),B2(NM2,NS2,N) accordingly------------------
      DO 22 J=1,NM2
        MU2=MU2MIN+(DBLE(J)-1.0000D0)*MU2L
       DO 42 K=1,NS2
         SG2=SG2MIN+(DBLE(K)-1.0000D0)*SG2L
        DO 62 L=1,N
          X=XMIN-LL+(DBLE(L)-1.0000D0)*XL
c          write(*,*) J,MU2,K,SG2,L,X,XL
            IF (MU2 .GE. 0D0) THEN
          A2(J,K,L)=(SG2*SG2*X*X)/(2.0d0*XL*XL)+MU2*X/XL
          B2(J,K,L)=(SG2*SG2*X*X)/(2.0d0*XL*XL)
            ELSE
          A2(J,K,L)=(SG2*SG2*X*X)/(2.0d0*XL*XL)
          B2(J,K,L)=(SG2*SG2*X*X)/(2.0d0*XL*XL)-MU2*X/XL
            end if
c          write(*,*) A2(J,K,L),B2(J,K,L)
c          pause "A2,B2"
 62      CONTINUE
 42   CONTINUE
 22   CONTINUE
      OPEN(UNIT=25, FILE='bf1.dat', STATUS="UNKNOWN")
c     read in data-------------
      OPEN(UNIT=1,
     &FILE='/usr/users/zeng/thesis/smlt/gsgj/gsgj.dat',      
     &STATUS='OLD')
      DO 6 I=1,NDATA
         READ(1,*) X2,DT
c         write(*,*) X2,dt
c     partition the waiting time into fine intervals
         NTI=IDNINT(DT/LT)+1
         DDT=DT/DBLE(NTI)
c         write(*,*) nti,ddt
c     -----------------------------------------------------
c     loop 101: updating Q1 for read-in datum from t_(i-1) to time t_i-
         DO 101 M1=1,NS1
            DO 111 L1=1,NM1
               DO 121 W1=1,NR1
         DO 125 J=1,NTI
          DO 51 K=1,N
            QI(K)=Q1(K,L1,M1,W1)
c            if (qi(k) .gt. 0) then
c               write(*,*) k,L1,M1,W1, qi(k)
c               pause
c             end if
 51      CONTINUE
c          pause "out loop 51"
                DO 131 K=1,N
            IF (K .EQ. 1) THEN
c               pause "k=1"
               Q1(K,L1,M1,W1)=QI(K)+(-(A1(L1,M1,K)+B1(L1,M1,K))*QI(K)
     &+B1(L1,M1,K+1)*QI(K+1))*DDT
               if (Q1(K,L1,M1,W1) .lt. 0) then
               write(*,*) "l1,m1,w1,j,k,nti,dt", l1,m1,w1,j,k,nti,dt
               write(*,*) "Q1(K,L1,M1,W1)",Q1(K,L1,M1,W1)
               write(*,*) "QI(k)",QI(k),"QI(k+1)",QI(k+1)
               write(*,*) "A1(L1,M1,K),B1(L1,M1,K),B1(L1,M1,K+1)",
     &A1(L1,M1,K),B1(L1,M1,K),B1(L1,M1,K+1)
               pause "Q1(K,L1,M1,W1) .lt. 0, k=1"
               end if
            ELSE IF (K .EQ. N) THEN
c               pause "k=n"
               Q1(K,L1,M1,W1)=QI(K)+(A1(L1,M1,K-1)*QI(K-1)-(A1(L1,M1,K)
     &+B1(L1,M1,K))*QI(K))*DDT
               if (Q1(K,L1,M1,W1) .lt. 0) then
              write(*,*) "l1,m1,w1,j,k,nti,dt", l1,m1,w1,j,k,nti,dt
              write(*,*) "Q1(K,L1,M1,W1)",Q1(K,L1,M1,W1),"QI(k)",QI(k)
              write(*,*) "QI(k-1),QI(k+1)", QI(k-1),QI(k+1)
      write(*,*) "A1(L1,M1,K-1),A1(L1,M1,K),B1(L1,M1,K),B1(L1,M1,K+1)",
     &A1(L1,M1,K-1),A1(L1,M1,K),B1(L1,M1,K)
     &,B1(L1,M1,K+1)
               pause "Q1(K,L1,M1,W1) .lt. 0, k=n"
               end if
            ELSE
               Q1(K,L1,M1,W1)=QI(K)+(A1(L1,M1,K-1)*QI(K-1)-(A1(L1,M1,K)
     &+B1(L1,M1,K))*QI(K)+B1(L1,M1,K+1)*QI(K+1))*DDT
               if (Q1(K,L1,M1,W1) .lt. 0) then
              write(*,*) "l1,m1,w1,j,k,nti,dt", l1,m1,w1,j,k,nti,dt
              write(*,*) "Q1(K,L1,M1,W1)",Q1(K,L1,M1,W1),"QI(k)",QI(k)
              write(*,*) "QI(k-1),QI(k+1)", QI(k-1),QI(k+1)
      write(*,*) "A1(L1,M1,K-1),A1(L1,M1,K),B1(L1,M1,K),B1(L1,M1,K+1)",
     &A1(L1,M1,K-1),A1(L1,M1,K),B1(L1,M1,K)
     &,B1(L1,M1,K+1)
               pause "Q1(K,L1,M1,W1) .lt. 0"
               end if
            END IF
 131     CONTINUE
 125  CONTINUE
 121  CONTINUE
 111  CONTINUE
 101  CONTINUE
c      pause "out loop 101"
c     loop 8: updating Q2 for read-in datum from t_(i-1) to time t_i-
      DO 8 J=1,NTI
c         write(*,*)  "inside loop 8"
c     ------------------------------------------------------
c     Define the mean of prob of sigma (Q2SM) for K,L,LD,W of X,MU,
c     LD,RO at each step
c     ------------------------------------------------------
      DO 10 K=1,N
       DO 30 L2=1,NM2
        DO 50 LB2=1,NL2
         DO 70 W2=1,NR2
            Q2SM(K,L2,LB2,W2)=0D0
            DO 90 M2=1,NS2
             Q2SM(K,L2,LB2,W2)=Q2SM(K,L2,LB2,W2)+Q2(K,L2,M2,LB2,W2)
 90         CONTINUE
             Q2SM(K,L2,LB2,W2)=Q2SM(K,L2,LB2,W2)/DBLE(NS2)
c             if ( (j .eq. 2) .and. (Q2SM(K,L2,LB2,W2) .gt. 0)) then
c                write(*,*) K,L2,LB2,W2,Q2SM(K,L2,LB2,W2)
c                pause
c               end if
 70          CONTINUE
 50          CONTINUE
 30          CONTINUE
 10          CONTINUE
c             pause "out loop 10"
         DO 102 M2=1,NS2
            DO 112 L2=1,NM2
               DO 122 W2=1,NR2
                DO 123 LB2=1,NL2
                 LBD2=LBD2MIN+(LB2-1.0d0)*LBD2L
          DO 126 K=1,N
            QI(K)=Q2(K,L2,M2,LB2,W2)
 126     CONTINUE
                DO 132 K=1,N
            IF (K .EQ. 1) THEN
c               pause "k=1"
        Q2(K,L2,M2,LB2,W2)=QI(K)+(-(A2(L2,M2,K)+B2(L2,M2,K))*QI(K)
     &+B2(L2,M2,K+1)*QI(K+1)+LBD2*(Q2SM(K,L2,LB2,W2)-QI(K)))*DDT
               if (Q2(K,L2,M2,LB2,W2) .lt. 0) then
               write(*,*) "l2,m2,w2,j,k,nti,dt", l2,m2,w2,j,k,nti,dt
               write(*,*) "Q2(K,L2,M2,LB2,W2)",Q2(K,L2,M2,LB2,W2)
               write(*,*) "QI(k)",QI(k),"QI(k+1)",QI(k+1)
      write(*,*) "A2(L2,M2,K),B2(L2,M2,K),B2(L2,M2,K+1)",A2(L2,M2,K),
     &B2(L2,M2,K),B2(L2,M2,K+1)
               pause "Q2(K,L2,M2,LB2,W2) .lt. 0"
               end if
            ELSE IF (K .EQ. N) THEN
c               pause "k=n"
      Q2(K,L2,M2,LB2,W2)=QI(K)+(A2(L2,M2,K-1)*QI(K-1)-(A2(L2,M2,K)
     &+B2(L2,M2,K))*QI(K)+LBD2*(Q2SM(K,L2,LB2,W2)-QI(K)))*DDT
               if (Q2(K,L2,M2,LB2,W2) .lt. 0) then
               write(*,*) "l2,m2,w2,j,k,nti,dt", l2,m2,w2,j,k,nti,dt
      write(*,*) "Q2(K,L2,M2,LB2,W2)",Q2(K,L2,M2,LB2,W2),"QI(k)",QI(k)
               write(*,*) "QI(k-1)", QI(k-1)
               write(*,*) "A2(L2,M2,K-1),A2(L2,M2,K),B2(L2,M2,K),",
     &A2(L2,M2,K-1),A2(L2,M2,K),B2(L2,M2,K)
               pause "Q2(K,L2,M2,LB2,W2) .lt. 0"
               end if
            ELSE
      Q2(K,L2,M2,LB2,W2)=QI(K)+(A2(L2,M2,K-1)*QI(K-1)-(A2(L2,M2,K)
     &+B2(L2,M2,K))*QI(K)+B2(L2,M2,K+1)*QI(K+1)+
     &LBD2*(Q2SM(K,L2,LB2,W2)-QI(K)))*DDT
               if (Q2(K,L2,M2,LB2,W2) .lt. 0) then
               write(*,*) "l2,m2,w2,j,k,nti,dt", l2,m2,w2,j,k,nti,dt
      write(*,*) "Q2(K,L2,M2,LB2,W2)",Q2(K,L2,M2,LB2,W2),"QI(k)",QI(k)
               write(*,*) "QI(k-1),QI(k+1)", QI(k-1),QI(k+1)
      write(*,*) "A2(L2,M2,K-1),A2(L2,M2,K),B2(L2,M2,K),B2(L2,M2,K+1)",
     &A2(L2,M2,K-1),A2(L2,M2,K),B2(L2,M2,K),B2(L2,M2,K+1)
               pause "Q2(K,L2,M2,LB2,W2) .lt. 0"
               end if
            END IF
 132     CONTINUE
 123  CONTINUE
 122  CONTINUE
 112  CONTINUE
 102  CONTINUE
c      pause "out loop 102"
 8       continue
c         pause "out loop 8"
c     ----------------------------------------------------------------
C     "IF STRUCTURE": ASSIGN FOR DIFFERENT INTENSITY FUNCTIONS TO OBTAIN
C     P(T_I)
       XB=X2*2D0
       XC=X2*4D0
c       write (*,*) 'Xb',xb,XC
C     CASE 1: X2 IS AN INTEGER OR HALF
       if (DMOD(XB,1d0) .EQ. 0) THEN
            CALL Q12J(LD11,LD12,Q1,Q2,X2)
C     CASE 3: X2 IS ODD QUARTERS
         ELSE IF (DMOD(XC,1d0) .EQ. 0) THEN
c            write(*,*) "INSIDE quarters"
            CALL Q12J(LD21,LD22,Q1,Q2,X2)
C     CASE 4: X2 IS ODD EIGHTH
         ELSE
C            write (*,*) "INSIDE 1/8"
            CALL Q12J(LD31,LD32,Q1,Q2,X2)
         END IF
         if ((i .eq. 1) .or. (i .eq. 10) .or. (i .eq. 1000) .or. 
     & (i .eq. 2166) .or. (i .eq. 2676)
     & .or. (i .eq. 7790) .or. (i .eq. 8113) .or. (i .eq. 8870)
     & .or. (i .eq. 11793) .or. (i .eq. 11972) .or. (i .eq. 13517)
     & .or. (i .eq. 15267) .or. (i .eq. 17161) .or. (i .eq. 19235)
     & .or. (i .eq. 20012) .or. (i .eq. 20299) .or. (i .eq. 22355)
     & .or. (i .eq. 24964) .or. (i .eq. 29856) .or. (i .eq. 31732)
     & .or. (i .eq. 31830) .or. (i .eq. 32053) .or. (i .eq. 34738)
     & .or. (i .eq. 35231) .or. (i .eq. 35767) .or. (i .eq. 38157)
     & .or. (i .eq. 39158) .or. (i .eq. 39607) .or. (i .eq. 41360)
     & .or. (i .eq. 43226) .or. (i .eq. 44478) .or. (i .eq. 51559)
     & .or. (i .eq. 53894) .or. (i .eq. 57307) .or. (i .eq. 58145)
     & .or. (i .eq. 61812) .or. (i .eq. 62334) .or. (i .eq. 62627)
     & .or. (i .eq. 64298) .or. (i .eq. 64900) .or. (i .eq. 71690)
     & .or. (i .eq. 73699) .or. (i .eq. 74971) .or. (i .eq. 75661)
     & .or. (i .eq. 76336) .or. (i .eq. 81416) .or. (i .eq. 84814)
     & .or. (i .eq. 90000))  then
            CALL WOUTPUT(Q1,Q2,I,X2,DT)
         end if
            CALL WOUTPUTbf(Q2,X2,DT)
 6           continue
             close(1)
             close(25)
      RETURN
      END

      SUBROUTINE Q12J(FUNC1,FUNC2,Q1,Q2,X)
      include 'parbf1.f'
      INTEGER K,L1,M1,W1,L2,M2,LB2,W2
      REAL*8 FUNC1,Q1(N,NM1,NS1,NR1),X,XA,TLK1, TL1,TLK2,TL2,
     &Q2(N,NM2,NS2,NL2,NR2),FUNC2
C     SUM UP THE TOTEL TLK1,TL1
      TLK1=0.0d0
      TL1=0.0D0
         DO 180 L1=1,NM1
            DO 182 M1=1,NS1
               DO 184 W1=1,NR1
               DO 186 K=1,N
                  XA=(DBLE(K)-1D0)*XL+XMIN-LL
                  TL1=TL1+Q1(K,L1,M1,W1)
                  Q1(K,L1,M1,W1)=Q1(K,L1,M1,W1)*FUNC1(X,XA,W1)
                 TLK1=TLK1+Q1(K,L1,M1,W1)
c                     IF (P(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,x,xa,P1(K,L,M,W)
c                      write(32,*) func(x,xa,w),func(x,xa,w)*p1(k,l,m,w)
c                      END IF
 186          CONTINUE
 184       CONTINUE
 182    CONTINUE
 180  CONTINUE
c      write(32,*) TL1,TLK1
          IF ((TL1 .EQ. 0D0) .or. (TLK1 .eq. 0)) THEN
             write(32,*) "TL1,TLK1",TL1,TLK1
             PAUSE "TL1=0"
          END IF
          TLK2=0.0d0
          TL2=0.0d0
         DO 190 L2=1,NM2
            DO 192 M2=1,NS2
               DO 194 W2=1,NR2
               DO 196 K=1,N
                  XA=(DBLE(K)-1D0)*XL+(XMIN-LL)
                DO 198 LB2=1,NL2
                  TL2=TL2+Q2(K,L2,M2,LB2,W2)
                  Q2(K,L2,M2,LB2,W2)=Q2(K,L2,M2,LB2,W2)*FUNC2(X,XA,W2)
                 TLK2=TLK2+Q2(K,L2,M2,LB2,W2)
 198          CONTINUE
 196       CONTINUE
 194    CONTINUE
 192  CONTINUE
 190  CONTINUE
c      write(32,*) TL2,TLK2
c     Obtain Q1 at jump time t_i
         DO 280 L1=1,NM1
            DO 282 M1=1,NS1
               DO 284 W1=1,NR1
               DO 286 K=1,N
                 Q1(K,L1,M1,W1)=Q1(K,L1,M1,W1)*TL2/TLK2
c                     IF (P(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,xa,P(K,L,M,W)
c                     END IF
 286          CONTINUE
 284       CONTINUE
 282    CONTINUE
 280  CONTINUE
c     Obtain Q2 at jump time t_i
         DO 290 L2=1,NM2
            DO 292 M2=1,NS2
               DO 294 W2=1,NR2
               DO 296 K=1,N
                 DO 298 LB2=1,NL2
                 Q2(K,L2,M2,LB2,W2)=Q2(K,L2,M2,LB2,W2)*TL1/TLK1
c                     IF (PJ(K,L,M,LD,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,xa,PJ(K,L,M,LD,W)
c                     END IF
 298             CONTINUE
 296          CONTINUE
 294       CONTINUE
 292    CONTINUE
 290  CONTINUE
      RETURN
      END

      SUBROUTINE WOUTPUT(Q1,Q2,I,X2,DT)
      INCLUDE 'parbf1.f'
      INTEGER I,K,L1,M1,W1,L2,M2,LB2,W2
      REAL*8 X2,DT,BF1,BF2,Q1(N,NM1,NS1,NR1),Q2(N,NM2,NS2,NL2,NR2)
      BF1=0D0
         DO 380 L1=1,NM1
            DO 382 M1=1,NS1
               DO 384 W1=1,NR1
               DO 386 K=1,N
                 BF1=BF1+Q1(K,L1,M1,W1)
 386          CONTINUE
 384       CONTINUE
 382    CONTINUE
 380  CONTINUE
      BF2=0D0
         DO 370 L2=1,NM2
            DO 372 M2=1,NS2
               DO 374 W2=1,NR2
               DO 376 K=1,N
                DO 378 LB2=1,NL2
                   BF2=BF2+Q2(K,L2,M2,LB2,W2)
 378            CONTINUE
 376         CONTINUE
 374      CONTINUE
 372   CONTINUE
 370  CONTINUE
      WRITE(32,*) I,X2,DT
      WRITE(32,*) BF1,BF2, BF1*BF2
      RETURN
      END
      
      SUBROUTINE WOUTPUTbf(Q2,X2,DT)
      INCLUDE 'parbf1.f'
      INTEGER K,L2,M2,LB2,W2
      REAL*8 X2,DT,BF2,Q2(N,NM2,NS2,NL2,NR2)
      BF2=0D0
         DO 360 L2=1,NM2
            DO 362 M2=1,NS2
               DO 364 W2=1,NR2
               DO 366 K=1,N
                DO 368 LB2=1,NL2
                   BF2=BF2+Q2(K,L2,M2,LB2,W2)
 368            CONTINUE
 366         CONTINUE
 364      CONTINUE
 362   CONTINUE
 360  CONTINUE
        write(25,*) BF2,X2,DT
      RETURN
      END

      FUNCTION LD11(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,LD11,Y,RHO
      INTEGER XA,W
c     add 0.000001 for protection of roundoff error
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO1MIN+(DBLE(W)-1.0D0)*RO1L
C      WRITE(*,*) "INSIDE LD1"
      IF (XA .EQ. 0) THEN
         LD11=(1.0D0-RHO)*(1.0D0+BETA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD11=0.50d0*(1.0d0-rho)*rho*rho*(1.0d0+rho*beta)
     &+(1.0d0-rho)*beta
      ELSE
         LD11=0.50d0*(1.0D0-RHO)*(rho**(IABS(XA)))*
     & (1.0D0+BETA*(RHO+1.0d0/RHO))
      END IF
C      write(*,*) X,Y,W,XA,RHO, ROMIN,ROL,BETA,LD1
C      PAUSE
      RETURN
      END

      FUNCTION LD21(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,LD21,Y,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO1MIN+(DBLE(W)-1.0D0)*RO1L
      IF (XA .EQ. 0) THEN
         LD21=(1.0D0-RHO)*(1.0D0+ALPHA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD21=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*alpha)
     &+(1.0d0-rho)*alpha
      ELSE
         LD21=0.50D0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+ALPHA*(RHO+1/RHO))
      END IF
      RETURN
      END

      FUNCTION LD31(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,Y,LD31,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO1MIN+(DBLE(W)-1.0D0)*RO1L
      IF (XA .EQ. 0) THEN
        LD31=(1.0D0-RHO)*(1.0D0-ALPHA-BETA)
      ELSE
        LD31=0.50D0*(1.0D0-ALPHA-BETA)*(1.0D0-RHO)
     &*(RHO**(IABS(XA)))
      END IF
      RETURN
      END

      FUNCTION LD12(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,LD12,Y,RHO
      INTEGER XA,W
c     add 0.000001 for protection of roundoff error
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO2MIN+(DBLE(W)-1.0D0)*RO2L
C      WRITE(*,*) "INSIDE LD1"
      IF (XA .EQ. 0) THEN
         LD12=(1.0D0-RHO)*(1.0D0+BETA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD12=0.50d0*(1.0d0-rho)*rho*rho*(1.0d0+rho*beta)
     &+(1.0d0-rho)*beta
      ELSE
         LD12=0.50d0*(1.0D0-RHO)*(rho**(IABS(XA)))*
     & (1.0D0+BETA*(RHO+1.0d0/RHO))
      END IF
C      write(*,*) X,Y,W,XA,RHO, ROMIN,ROL,BETA,LD1
C      PAUSE
      RETURN
      END

      FUNCTION LD22(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,LD22,Y,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO2MIN+(DBLE(W)-1.0D0)*RO2L
      IF (XA .EQ. 0) THEN
         LD22=(1.0D0-RHO)*(1.0D0+ALPHA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD22=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*alpha)
     &+(1.0d0-rho)*alpha
      ELSE
         LD22=0.50D0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+ALPHA*(RHO+1/RHO))
      END IF
      RETURN
      END

      FUNCTION LD32(X,Y,W)
      INCLUDE 'parbf1.f'
      REAL*8 X,Y,LD32,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=RO2MIN+(DBLE(W)-1.0D0)*RO2L
      IF (XA .EQ. 0) THEN
        LD32=(1.0D0-RHO)*(1.0D0-ALPHA-BETA)
      ELSE
        LD32=0.50D0*(1.0D0-ALPHA-BETA)*(1.0D0-RHO)
     &*(RHO**(IABS(XA)))
      END IF
      RETURN
      END