*    
*     bayesfactors-jsvgbm-gbm.f
*     (Bayes factors for a micromovement model with 
*      ^^^^^ ^^^^^^^  
*     Jumping Stochastic Volatility Geometric Brownian Motion as Value Process
*     ^       ^          ^          ^         ^        ^
*     and that with Geometric Brownian Motion as Value Process
*                   ^         ^        ^
*     Copyright (C) 2005
*
*     Yong Zeng
*     Department of Mathematics and Statistics
*     University of Missouri at Kansas City
*     5100 Rockhill Road
*     Kansas City, MO 64110-2499 USA.
*
*     Permission to use, copy, modify, and distribute this
*     software for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*----------------------------------------------------------------
C     the program computes the trade-by-trade Bayes factors 
c     for the JSV-GBM vs GBM models described
C     in the paper M. Kouritzin and Yong Zeng (2005), "Bayesian Model 
C     Selection via Filtering for a Class of Micro-movement Models of Asset 
C     Price" International Journal of Theoretical and Applied Finance (IJTAF) 
C     (2005) Vol.8 No.1, pp. 97-122. 
*----------------------------------------------------------------
c     "parbf.f" specifies the parameters of this program
c     simulated date set: smlt-jsvgbm.dat (it has two columes: price, and 
c     duration or waiting time between trades)
c     For a different data set such as MSFT940102.dat (MicroSoft,
c     Jan. and Feb., 1994), you need to reset the ranges for X, mu, and etc. 
c     All files should be in the same directory
c     > f77  bayesfactors-jsvgbm-gbm.f  (to generate "a.out")
c     > a.out  (to run)
c     Warning: it may take a long time to finish the program.
C----------------------------------------------------------------
C     PROGRAM FOR COMPUTING THE BAYES FACTOR FOR GBM AND JSVGBM
C     Q1=LIKELIHOOD RATIO OF GBM AND JSVGBM
C     Q2=LIKELIHOOD RATIO OF JSVGBM 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     ########################################################
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)
 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 'parbf.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
            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
 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
            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
 62      CONTINUE
 42   CONTINUE
 22   CONTINUE
c     read in data-------------
      OPEN(UNIT=1, FILE='smlt-jsvgbm.dat', STATUS='OLD')
      DO 6 I=1,NDATA
         READ(1,*) X2,DT
c     partition the waiting time into fine intervals
         NTI=IDNINT(DT/LT)+1
         DDT=DT/DBLE(NTI)
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)
 51      CONTINUE
                DO 131 K=1,N
            IF (K .EQ. 1) THEN
               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
            ELSE IF (K .EQ. N) THEN
               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
            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
            END IF
 131     CONTINUE
 125  CONTINUE
 121  CONTINUE
 111  CONTINUE
 101  CONTINUE
c     loop 8: updating Q2 for read-in datum from t_(i-1) to time t_i-
      DO 8 J=1,NTI
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)
 70          CONTINUE
 50          CONTINUE
 30          CONTINUE
 10          CONTINUE
         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
        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
            ELSE IF (K .EQ. N) THEN
      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
            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
            END IF
 132     CONTINUE
 123  CONTINUE
 122  CONTINUE
 112  CONTINUE
 102  CONTINUE
 8       continue
            CALL WOUTPUT(Q1,Q2,I,X2,DT)
c     ----------------------------------------------------------------
C     "IF STRUCTURE": ASSIGN FOR DIFFERENT INTENSITY FUNCTIONS TO OBTAIN
C     P(T_I)
       XB=X2*2D0
       XC=X2*4D0
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
            CALL Q12J(LD21,LD22,Q1,Q2,X2)
C     CASE 4: X2 IS ODD EIGHTH
         ELSE
            CALL Q12J(LD31,LD32,Q1,Q2,X2)
         END IF
c    For the points where the volatility is about to changes, 
c    middle-step output is written.
         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
 6           continue
             close(1)
      RETURN
      END

      SUBROUTINE Q12J(FUNC1,FUNC2,Q1,Q2,X)
      include 'parbf.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)
 186          CONTINUE
 184       CONTINUE
 182    CONTINUE
 180  CONTINUE
          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     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
 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
 298             CONTINUE
 296          CONTINUE
 294       CONTINUE
 292    CONTINUE
 290  CONTINUE
      RETURN
      END

      SUBROUTINE WOUTPUT(Q1,Q2,I,X2,DT)
      INCLUDE 'parbf.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
      
      FUNCTION LD11(X,Y,W)
      INCLUDE 'parbf.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
      RETURN
      END

      FUNCTION LD21(X,Y,W)
      INCLUDE 'parbf.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 'parbf.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 'parbf.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
      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
      RETURN
      END

      FUNCTION LD22(X,Y,W)
      INCLUDE 'parbf.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 'parbf.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