      SUBROUTINE CHECKINPUT_SB
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/CONST5/ILO,INLO,IONLYLO
      COMMON/CUT1/PTMIN,PTMAX
      COMMON/CUT2/YMIN,YMAX
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

C***  CHECK MS,MG AND MT ( SHOULD BE ABOVE 1 GEV )
      IF (MS.LT.1D0) THEN
         MS = 1D0
         PRINT *,'MS CHANGED TO 1.'
      END IF
      IF (MG.LT.1D0) THEN
         MG = 1D0
         PRINT *,'MG CHANGED TO 1.'
      END IF
      IF (MT.LT.1D0) THEN
         MT = 1D0
         PRINT *,'MT CHANGED TO 1.'
      END IF

C***  CHECK ICOLL
      IF ((ICOLL.LT.0).OR.(ICOLL.GT.1)) THEN
         PRINT *,'COLLIDER TYPE SPECIFIED INCORRECTLY; ICOLL = ',ICOLL
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK ENERGY ( SHOULD BE ABOVE PRODUCTION THRESHOLD )
      IF (ENERGY.LT.2*MS) THEN
         PRINT *,'ENERGY TOO SMALL TO PRODUCE SQUARK-ANTISQUARK'
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK IFLAVOR
      IF ((IFLAVOR.LT.0).OR.(IFLAVOR.GT.5)) THEN
         PRINT *,'IFLAVOR IN THE WRONG RANGE; IFLAVOR = ',IFLAVOR
         PRINT *,'IFLAVOR CHANGED TO 0 (ALL)'
         IFLAVOR = 0
      END IF

C***  CHECK ITOTAL
      IF ((ITOTAL.LT.0).OR.(ITOTAL.GT.1)) THEN
         PRINT *,'ITOTAL IN THE WRONG RANGE; ITOTAL = ',ITOTAL
         PRINT *,'ITOTAL CHANGED TO 1'
         ITOTAL = 1
      END IF

C***  CHECK ISCAPT AND SCAFAC
      IF ((ISCAPT.LT.0).OR.(ISCAPT.GT.1)) THEN
         PRINT *,'ISCAPT IN THE WRONG RANGE; ISCAPT = ',ISCAPT
         PRINT *,'ISCAPT CHANGED TO 0'
         ISCAPT = 0
      END IF
      IF (SCAFAC.LT.1D-2) THEN
         PRINT *, 'SCAFAC TOO SMALL'
         PRINT *, 'SCAFAC CHANGED TO 1.'
         SCAFAC = 1.D0
      END IF
      IF ((ITOTAL.EQ.1).AND.(ISCAPT.EQ.1)) THEN
         PRINT *,'WARNING: TRANSVERSE MASS IS NOT USED AS A SCALE'
      END IF

C***  CHECK PTMIN AND PTMAX
      IF (PTMIN.LT.0D0) THEN
         PTMIN = 0D0
         PRINT *,'PTMIN NEGATIVE'
         PRINT *,'PTMIN CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.0D0) THEN
         PTMAX = 0D0
         PRINT *,'PTMAX NEGATIVE'
         PRINT *,'PTMAX CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.PTMIN) THEN
         PTMAX = PTMIN
         PRINT *,'PTMAX SMALLER THAN PTMIN'
         PRINT *,'PTMAX CHANGED TO PTMIN'
      END IF
      IF ((PTMAX.EQ.PTMIN).AND.(PTMAX.EQ.0D0)) THEN
         PTMAX = 0.1D0
         PTMIN = 0.1D0
         PRINT *,'PTMAX = PTMIN CHANGED TO 0.1'
      END IF

C***  CHECK YMIN AND YMAX
      IF (YMIN.LT.0D0) THEN
         YMIN = 0D0
         PRINT *,'YMIN NEGATIVE'
         PRINT *,'YMIN CHANGED TO 0.'
      END IF
      IF (YMAX.LT.0D0) THEN
         YMAX = 0D0
         PRINT *,'YMAX NEGATIVE'
         PRINT *,'YMAX CHANGED TO 0.'
      END IF
      IF (YMAX.LT.YMIN) THEN
         YMAX = YMIN
         PRINT *,'YMAX SMALLER THAN YMIN'
         PRINT *,'YMAX CHANGED TO YMIN'
      END IF

C***  CHECK ICHACONJ
      IF ((ICHARCONJ.LT.-1).OR.(ICHARCONJ.GT.1)) THEN
         PRINT *,'ICHARCONJ IN THE WRONG RANGE; ICHARCONJ = ',ICHARCONJ
         PRINT *,'ICHARCONJ CHANGED TO 0'
         ICHARCONJ = 0
      END IF

C***  CHECK IONLYLO
      IF ((IONLYLO.LT.0).OR.(IONLYLO.GT.1)) THEN
         PRINT *,'IONLYLO IN THE WRONG RANGE'
         PRINT *,'IONLYLO CHANGED TO 0'
         IONLYLO = 0
      END IF

C***  CHECK ILO, INLO
      IF (ILO.LT.10) THEN
         ILO = 10
         PRINT *,'ILO CHANGED TO 10'
      END IF
      IF (INLO.LT.10) THEN
         INLO = 10
         PRINT *,'INLO CHANGED TO 10'
      END IF

      RETURN
      END


C ======================================================================

      SUBROUTINE INTEGSB(RESLO,ERRLO,RESNLO,ERRNLO)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      integer  ivegas(1:4)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST5/ILO,INLO,IONLYLO

      EXTERNAL SIGLOSB, SIGNLOSB

      CALL CHECKINPUT_SB

      S = ENERGY**2

      CALL INILO(ivegas)
      idim = 3
      acc = 1.d-4
      call INTEG(SIGLOSB,idim,ivegas,acc,reslo,errlo)

      IF (IONLYLO.EQ.0) THEN
         CALL ININLO(ivegas)
         idim = 5 
         call INTEG(SIGNLOSB,idim,ivegas,acc,resnlo,errnlo)
      ELSE 
         RESNLO = 0D0
         ERRNLO = 0D0
      END IF

      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN LO

      REAL*8 FUNCTION SIGLOSB(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(3)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ
      CALL DVEGLO_SB(VAR)
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYSB
         CALL DEFSCASB
         CALL DDILOSB
         SIGLOSB1 = ESBGGLO()
         SIGLOSB2 = ESBQBLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            YVAL = -YVAL
            CALL DDILOSB
            SIGLOSB1 = 2 * SIGLOSB1
            SIGLOSB2 = SIGLOSB2 + ESBQBLO()
         ELSE 
            SIGLOSB1 = 2 * SIGLOSB1
            SIGLOSB2 = 2 * SIGLOSB2
         END IF
      ELSE
         CALL DEFSCASB
         CALL DTOLOSB
         SIGLOSB1 = FSBGGLO()
         SIGLOSB2 = FSBQBLO()
      END IF
      SIGLOSB =  SUMLOSB(SIGLOSB1,SIGLOSB2,IFLAVOR)
      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN NLO

      REAL*8 FUNCTION SIGNLOSB(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(5)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

      CALL DVEGNLO_SB(VAR)
      CALL DMINLSB
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYSB
         CALL DEFSCASB
         CALL DDINLSB
         CALL DD4NLSB
         CALL DD3NLSB
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4))      
     +        SIGNLOSB1 = ESBGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOSB2 = ESBQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +        SIGNLOSB3 = ESBGBNLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            SIGNLOSB1 = 2 * SIGNLOSB1
            YVAL = -YVAL
            CALL DDINLSB
            CALL DD4NLSB
            CALL DD3NLSB
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +           SIGNLOSB2 = SIGNLOSB2 + ESBQBNLO()
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +           SIGNLOSB3 = SIGNLOSB3 + ESBGBNLO()
         ELSE 
            SIGNLOSB1 = 2 * SIGNLOSB1
            SIGNLOSB2 = 2 * SIGNLOSB2
            SIGNLOSB3 = 2 * SIGNLOSB3
         END IF
      ELSE
         CALL DEFSCASB
         CALL DTONLSB
         CALL DTSNLSB
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOSB1 = FSBGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOSB2 = FSBQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5)) 
     +        SIGNLOSB3 = FSBGBNLO()
      END IF
      SIGNLOSB =  SUMNLOSB(SIGNLOSB1,SIGNLOSB2,SIGNLOSB3,IFLAVOR)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGLO_SB(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 3 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(3)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGNLO_SB(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 5 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(5)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      Y4 = VAR(4)
      Y5 = VAR(5)
      RETURN
      END

C ======================================================================

      SUBROUTINE SWAP_SB(X,Y)
C***  INTERCHANGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      Z = X
      X = Y
      Y = Z
      RETURN
      END

C ======================================================================

      SUBROUTINE AVER_SB(X,Y)
C***  AVERAGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      X = ( X + Y )/2D0
      Y = X
      RETURN
      END

C ======================================================================

      SUBROUTINE DEFSCASB()
C***  DEFINES THE RENORMALIZATION AND FACTORIZATION SCALE
C***  FILLS UP PART OF THE COMMON BLOCKS CONST1 AND CONST2
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL

      IF ((ITOTAL.EQ.0).AND.(ISCAPT.EQ.1)) THEN
         SCALE = MTRANS * SCAFAC
      ELSE 
         SCALE = MS * SCAFAC
      END IF

C***  THE SCALE SHOULD BE BELOW 1 TEV, BECAUSE SOME PARTON DENSITIES
C***  ARE RESTRICTED TO THIS REGION

      IF (SCALE.GE.1000D0) SCALE = 1000D0

      IF (SCALE.LE.5D0) SCALE = 5D0

      ALPHAS = ALPS(SCALE)

      RETURN
      END

C ======================================================================

      SUBROUTINE DPTYSB()
C***  DEFINES PT AND Y (INCLUDING CUTS) 
C***  FILLS UP THE COMMON BLOCK PTYSB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CUT1/PTMINC,PTMAXC
      COMMON/CUT2/YMINC,YMAXC
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      PT2MAX = S/4.D0 -MS**2
      PTMAX = SQRT( DMIN1(PT2MAX,PTMAXC**2) )
      PTMIN = PTMINC
      IF (PTMAX.GT.PTMIN) THEN
         PT = (PTMAX-PTMIN) * Y1 +PTMIN
         PREPT = 2.D0 * PT * (PTMAX -PTMIN)
      ELSE IF (PTMAX.EQ.PTMIN) THEN
         PT = PTMIN
         PREPT = 2.D0 * PT
      ELSE IF (PTMAX.LT.PTMIN) THEN 
         PT = 0D0
         PREPT = 0D0
      END IF
      MTRANS = SQRT(PT**2 +MS**2)
      YMAXK = +DACOSH(SQRT(S)/2.D0/MTRANS)
      YMINK = -DACOSH(SQRT(S)/2.D0/MTRANS)
      YMAX = DMIN1( YMAXK, YMAXC)
      YMIN = DMAX1( YMINK, YMINC, 0D0)

      IF (YMAX.GT.YMIN) THEN
         YVAL = (YMAX -YMIN) * Y2 + YMIN
         PREY = YMAX - YMIN
      ELSE IF (YMAX.EQ.YMIN) THEN         
         YVAL = YMIN
         PREY = 1.D0
      ELSE IF (YMAX.LT.YMIN) THEN
         YVAL = 0D0
         PREY = 0D0
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE DDILOSB()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK DILOSB AND PDLOSB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/DILOSB/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      X2B= -X1*U1/(X1*S +T1)
      SHATB  = X1 *X2B*S
      T1HATB = X2B*T1
      U1HATB = X1 *U1
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFLOSB(X1,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTOLOSB()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK TOLOSB AND PDLOSB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/TOLOSB/SHAT,T1HAT,PREF,X1,X2
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      X1 = (4*MS**2/S) ** Y1
      X2 = (4*MS**2/S/X1) ** Y2
      SHAT = X1*X2*S
      BETA = SQRT(1D0 -4*MS**2/SHAT)
      T1HAT = -SHAT/2D0 +(Y3 -0.5D0)*SHAT*BETA
      PREF = SHAT*BETA *LOG(4*MS**2/S)*LOG(4*MS**2/S/X1)
      CALL PDFLOSB(X1,X2)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBGGLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DILOSB/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      ESBGGLO = 
     +     PDGGB *PREFB * DSBGGB(ALPHAS,SHATB,T1HATB,MS,MG) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBQBLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DILOSB/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/CHARCONJ/ICHARCONJ
      
C***  DISTINGUISH BETWEEN SQUARKS AND ANTISQUARKS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- SQUARKS
C***  ICHARCONJ = -1  <-- ANTISQUARKS
C***  ICHARCONJ =  0  <-- AVERAGE OF SQUARKS AND ANTISQUARKS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_SB(PDQBB, PDBQB)
         CALL SWAP_SB(PDQBPB, PDBQPB)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_SB(PDQBB, PDBQB)
         CALL AVER_SB(PDQBPB, PDBQPB)
      END IF

      ESBQBLO = 
     +     +  PDQBB * PREFB*DSBQBB(ALPHAS,SHATB,T1HATB,MS,MG,1) 
     +     +  PDQBPB* PREFB*DSBQBB(ALPHAS,SHATB,T1HATB,MS,MG,0) 
     +     +  PDBQB * PREFB*DSBQBB(ALPHAS,SHATB,U1HATB,MS,MG,1) 
     +     +  PDBQPB* PREFB*DSBQBB(ALPHAS,SHATB,U1HATB,MS,MG,0) 

      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSBGGLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/TOLOSB/SHAT,T1HAT,PREF,X1,X2
      FSBGGLO = 
     +     PDGGB * PREF*DSBGGB(ALPHAS,SHAT,T1HAT,MS,MG) 
      RETURN
      END

C ====================================================================== 
      REAL*8 FUNCTION FSBQBLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/TOLOSB/SHAT,T1HAT,PREF,X1,X2
      FSBQBLO = 
     +     +  (PDQBB +PDBQB ) * PREF*DSBQBB(ALPHAS,SHAT,T1HAT,MS,MG,1) 
     +     +  (PDQBPB+PDBQPB) * PREF*DSBQBB(ALPHAS,SHAT,T1HAT,MS,MG,0) 
      RETURN
      END

C ======================================================================
 
      REAL*8 FUNCTION SUMLOSB(SIGLOSB1,SIGLOSB2,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN LO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.4))
     +                  SUMLOSB = SIGLOSB1 + SIGLOSB2
      IF (IFLAVOR.EQ.1) SUMLOSB = SIGLOSB1
      IF (IFLAVOR.EQ.2) SUMLOSB = SIGLOSB2
      IF ((IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5)) SUMLOSB = 0D0
      RETURN
      END

C ======================================================================

      SUBROUTINE DMINLSB()
C***  DEFINES SOME MISCELLANEOUS VARIABLES FOR CROSS-SECTION IN NLO
C***  FILLS UP THE COMMON BLOCK MINLSB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
C***  EPS = ( GAMMA(GLUINO)/MASS(GLUINO) )**2
      EPS = 1.D-5
      DEL = 1.D-4 * MS**2 
      SCA = SCALE**2/MS**2
      M2 = MG**2 -MS**2
      GAM = SQRT(EPS)*MG**2
      RETURN
      END

C ======================================================================

      SUBROUTINE DDINLSB()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCK DINLSB1-2 AND PDLOSB, PDNLSB1
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/PDNLSB1/PDGG,PDQB,PDQBP ,PDBQ ,PDBQP,PDQG, PDGQ,PDBG, PDGB

      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      S4MAX = X1 * (S+U1) +T1
      S4 = (S4MAX -DEL)* Y4 +DEL
      X2 = (S4 -X1*U1)/(X1*S +T1)
      X2B= -X1*U1/(X1*S +T1)
      SHAT   = X1 *X2 *S
      SHATB  = X1 *X2B*S
      T1HAT  = X2 *T1
      T1HATB = X2B*T1
      U1HAT  = X1 *U1
      U1HATB = X1 *U1
      PREF  = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2 )
     +     *(S4MAX -DEL)/(S4-X1*U1)
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFNLOSB1(X1,X2,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DD4NLSB()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN NLO
C***  INCLUDING SUBTRACTION IN S4
C***  FILLS UP THE COMMON BLOCKS DINLSB3-4 AND PDNLSB2A-B
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB3/SHAT0,T1HAT0,U1HAT0,S40,X10,X20,PREF0,S4G02
      COMMON/DINLSB4/SHAT1,T1HAT1,U1HAT1,S41,X11,X21,PREF1,S4UP
      COMMON/PDNLSB2A/PDGB0,PDGB1, PDBG0,PDBG1
      COMMON/PDNLSB2B/PDGQ0,PDGQ1, PDQG0,PDQG1

      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      S4MAX0 = S + T1 +U1
      TYMIN = ATAN(-M2/GAM)
      TYMAX = ATAN((S4MAX0-M2)/GAM)
      TY = (TYMAX -TYMIN) * Y3 + TYMIN 
      PREFTY = (TYMAX -TYMIN)/GAM 
      S40 = GAM*TAN(TY) + M2
      S41 = M2

      X1MIN0 = (S40 -T1)/(S +U1)
      X1MIN1 = (S41 -T1)/(S +U1)
      X10 = X1MIN0**Y4
      X11 = X1MIN1**Y4
      X20 = (S40 -X10*U1)/(X10*S +T1)
      X21 = (S41 -X11*U1)/(X11*S +T1)
      SHAT0  = X10 * X20 * S
      SHAT1  = X11 * X21 * S
      T1HAT0 = X20 * T1
      T1HAT1 = X21 * T1
      U1HAT0 = X10 * U1
      U1HAT1 = X11 * U1
      S4G02  = (S40 -MG**2 +MS**2)**2 + EPS * MG**4
      S4UP   = SHAT1 - 2D0*SQRT(SHAT1*MS**2)
      PREF0  = S * PREPT * PREY * LOG(X1MIN0)*(-X10*X20)/(S40-X10*U1)
     +     *PREFTY
      PREF1  = S * PREPT * PREY * LOG(X1MIN1)*(-X11*X21)/(S41-X11*U1)
     +     *PREFTY
      CALL PDFNLOSB2(X10,X20,X11,X21)
      RETURN
      END

C ======================================================================

      SUBROUTINE DD3NLSB()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN NLO
C***  INCLUDING SUBTRACTION IN S3
C***  FILLS UP THE COMMON BLOCKS DINLSB5-6 AND PDNLSB3A-B
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PTYSB/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/DINLSB5/SHAT2,T1HAT2,U1HAT2,S42,X12,X22,PREF2,S32,S3G2
      COMMON/DINLSB6/SHAT3,T1HAT3,U1HAT3,S43,X13,X23,PREF3,S3G3,S3UP
      COMMON/PDNLSB3A/PDQG2,PDQG3,PDGQ2,PDGQ3 
      COMMON/PDNLSB3B/PDBG2,PDBG3,PDGB2,PDGB3 

      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      S3MAX2 = (S+T1+U1)/2D0/(S+T1+U1+MS**2) * 
     +     (-T1 -U1 -2*MS**2 +SQRT((T1+U1)**2 -4*MS**2 *S))

      TZMIN = ATAN(-M2/GAM)
      TZMAX = ATAN((S3MAX2-M2)/GAM)
      TZ = (TZMAX -TZMIN) * Y3 + TZMIN 
      PREFTZ = (TZMAX -TZMIN)/GAM 
      S32 = GAM*TAN(TZ) + M2
      S33 = M2

      X1MIN2 = -(2*MS**2*S*S32 +S*S32**2 +2*MS**2*S*T1 +S*S32*T1 
     +     +2*MS**2*S32*U1 +S32**2*U1 +2*MS**2*T1*U1 +2*S32*T1*U1
     +     +S32*SQRT(4*MS**2*S**2*S32 +S**2*S32**2 +4*MS**2*S**2*T1 
     +     +2*S**2*S32*T1 +S**2*T1**2 +4*MS**2*S*S32*U1 +2*S*S32**2*U1 
     +     +4*MS**2*S*T1*U1 +2*S*S32*T1*U1 +S32**2*U1**2))/
     +     2D0/(S+U1)/(MS**2*S+MS**2*U1 +S32*U1)

      X1MIN3 = -(2*MS**2*S*S33 +S*S33**2 +2*MS**2*S*T1 +S*S33*T1 
     +     +2*MS**2*S33*U1 +S33**2*U1 +2*MS**2*T1*U1 +2*S33*T1*U1
     +     +S33*SQRT(4*MS**2*S**2*S33 +S**2*S33**2 +4*MS**2*S**2*T1 
     +     +2*S**2*S33*T1 +S**2*T1**2 +4*MS**2*S*S33*U1 +2*S*S33**2*U1 
     +     +4*MS**2*S*T1*U1 +2*S*S33*T1*U1 +S33**2*U1**2))/
     +     2D0/(S+U1)/(MS**2*S+MS**2*U1 +S33*U1)

      X12 = X1MIN2**Y4
      X13 = X1MIN3**Y4

      SQRS42 = SQRT(
     +  (-2*MS**2*T1 -S32*T1 -2*MS**2*X12*S -X12*S*S32 -X12**2*S*U1)**2 
     +  -4*MS**2*(T1+S*X12)*(MS**2*T1 +S32*T1 +MS**2*X12*S) )
      SQRS43 = SQRT(
     +  (-2*MS**2*T1 -S33*T1 -2*MS**2*X13*S -X13*S*S33 -X13**2*S*U1)**2 
     +  -4*MS**2*(T1+S*X13)*(MS**2*T1 +S33*T1 +MS**2*X13*S) )

      S42MIN = S32/2D0/(MS**2*(T1+X12*S) +S32*T1) * 
     +     (-2*MS**2*T1 -S32*T1 -2*MS**2*X12*S -X12*S*S32 -X12**2*S*U1 
     +     -SQRS42 )

      S43MIN = S33/2D0/(MS**2*(T1+X13*S) +S33*T1) * 
     +     (-2*MS**2*T1 -S33*T1 -2*MS**2*X13*S -X13*S*S33 -X13**2*S*U1 
     +     -SQRS43 )

      S32DEC = (T1 +S*X12 +U1*X12)/2D0/(T1 +S*X12 +U1*X12 +MS**2) * (
     +     -2*MS**2 -T1 -X12*U1 -SQRT((T1+U1*X12)**2 -4*MS**2*S*X12))
      S33DEC = (T1 +S*X13 +U1*X13)/2D0/(T1 +S*X13 +U1*X13 +MS**2) * (
     +     -2*MS**2 -T1 -X13*U1 -SQRT((T1+U1*X13)**2 -4*MS**2*S*X13))
      
      IF (S32.GE.S32DEC) THEN
         S42MAX = X12*(S+U1) +T1
      ELSE 
         S42MAX = S32/2D0/(MS**2*(T1+X12*S) +S32*T1) * 
     +     (-2*MS**2*T1 -S32*T1 -2*MS**2*X12*S -X12*S*S32 -X12**2*S*U1 
     +     +SQRS42 )
      END IF
      IF (S33.GE.S33DEC) THEN
         S43MAX = X13*(S+U1) +T1
      ELSE 
         S43MAX = S33/2D0/(MS**2*(T1+X13*S) +S33*T1) * 
     +     (-2*MS**2*T1 -S33*T1 -2*MS**2*X13*S -X13*S*S33 -X13**2*S*U1 
     +     +SQRS43 )
      END IF

      PREFS42 = S42MAX -S42MIN
      PREFS43 = S43MAX -S43MIN
      S42 = PREFS42 * Y5 +S42MIN
      S43 = PREFS43 * Y5 +S43MIN

      S3G2 = S32 - M2
      S3G3 = 0D0
      S3G22 = S3G2**2 + EPS*MG**4

      X22 = (S42 -X12*U1)/(X12*S +T1)
      X23 = (S43 -X13*U1)/(X13*S +T1)
      SHAT2 = X12 * X22 * S
      SHAT3 = X13 * X23 * S
      T1HAT2 = X22 * T1
      T1HAT3 = X23 * T1
      U1HAT2 = X12 * U1
      U1HAT3 = X13 * U1
      PREX2 = S42/2D0/(S42+MS**2)*SQRT((SHAT2-S42)**2-4*MS**2*SHAT2)
      PREX3 = S43/2D0/(S43+MS**2)*SQRT((SHAT3-S43)**2-4*MS**2*SHAT3)
      S3UP   = SHAT3 - 2D0*SQRT(SHAT3*MS**2)

      PREF2  = S * PREPT * PREY * PREFTZ/PREX2 * LOG(X1MIN2)*(-X12*X22)
     +    * PREFS42/(S42-X12*U1)
      PREF3  = S * PREPT * PREY * PREFTZ/PREX3 * LOG(X1MIN3)*(-X13*X23)
     +    * PREFS43/(S43-X13*U1)
      CALL PDFNLOSB3(X12,X22,X13,X23)
      RETURN
      END

C ======================================================================
      
      SUBROUTINE DTONLSB()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCKS DINLSB1-2 AND PDNLSBT
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDNLSBT/PDGGT, PDQBT, PDQBPT, PDGQT
      X1 = (4*MS**2/S) ** Y1
      X2 = (4*MS**2/S/X1) ** Y2
      SHAT = X1*X2*S
      DEL = 1.D-4 * MS**2 * (1D0 -4*MS**2/SHAT)
      BETA = SQRT(1D0 -4*MS**2/SHAT)
      T1HAT = 
     +   -(SHAT -DEL)/2D0 +(Y3 -0.5D0)*SQRT((SHAT-DEL)**2-4*MS**2*SHAT)
      S4MAX = SHAT +T1HAT +SHAT*MS**2/T1HAT
      S4 = (S4MAX -DEL)*Y4 +DEL
      PREF  = SHAT*BETA *S4MAX*LOG(4*MS**2/S)*LOG(4*MS**2/S/X1)
      PREFB = SHAT*BETA *      LOG(4*MS**2/S)*LOG(4*MS**2/S/X1)
      CALL PDFNLOSBT(X1,X2)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTSNLSB()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN NLO
C***  INCLUDING SUBTRACTION
C***  FILLS UP THE COMMON BLOCK DINLSB3-4 
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4,Y5
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/DINLSB3/SHAT0,T1HAT0,U1HAT0,S40,X10,X20,PREF0,S4G02
      COMMON/DINLSB4/SHAT1,T1HAT1,U1HAT1,S41,X11,X21,PREF1,S4UP

      S4MAX0 = SHAT - 2.D0*SQRT(SHAT*MS**2)
      TYMIN = ATAN(-M2/GAM)
      TYMAX = ATAN((S4MAX0-M2)/GAM)
      TY = (TYMAX -TYMIN) * Y3 + TYMIN 
      PREFTY = (TYMAX -TYMIN)/GAM 
      S40 = GAM*TAN(TY) + M2
      S41 = M2
      PRET10 = SQRT( (SHAT-S40)**2 -4*MS**2*SHAT )
      PRET11 = 1
      IF((SHAT-S41)**2 -4*MS**2*SHAT.GT.0.D0)THEN
       PRET11 = SQRT( (SHAT-S41)**2 -4*MS**2*SHAT )
      ENDIF
      T1HAT0 = -(SHAT-S40)/2.D0 +(Y4 -0.5D0)*PRET10
      T1HAT1 = -(SHAT-S41)/2.D0 +(Y4 -0.5D0)*PRET11
      PREF0 = PRET10 * PREFTY*LOG(4*MS**2/S)*LOG(4*MS**2/S/X1)
      PREF1 = PRET11 * PREFTY*LOG(4*MS**2/S)*LOG(4*MS**2/S/X1)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBGGNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/PDNLSB1/PDGG,PDQB,PDQBP ,PDBQ ,PDBQP,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      ESBGGNLO = 
     +        +PDGGB * PREFB * DSBGGB(ALPHAS,SHATB,T1HATB,MS,MG) 
     +        +PDGGB * PREFB * DSBGGV(ALPHAS,SHATB,T1HATB,MS,MG,MT) 
     +        +PDGGB * PREFB * DSBGG1(ALPHAS,SHATB,T1HATB,MS,MG,SCA) 
     +        +PDGGB * PREFB * 
     +        DSBGGD(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX)
     +        +PDGGB * PREFB * 
     +        DSBGG2(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX,SCA)
     +        +PDGG  * PREF  * DSBGGH(ALPHAS,SHAT,T1HAT,S4,MS,MG) 
     +        +PDGG  * PREF  * DSBGG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBQBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/PDNLSB1/PDGG,PDQB,PDQBP ,PDBQ ,PDBQP,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN SQUARKS AND ANTISQUARKS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- SQUARKS
C***  ICHARCONJ = -1  <-- ANTISQUARKS
C***  ICHARCONJ =  0  <-- AVERAGE OF SQUARKS AND ANTISQUARKS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_SB(PDQBB, PDBQB)
         CALL SWAP_SB(PDQB , PDBQ )
         CALL SWAP_SB(PDQBPB, PDBQPB)
         CALL SWAP_SB(PDQBP , PDBQP )
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_SB(PDQBB, PDBQB)
         CALL AVER_SB(PDQB , PDBQ )
         CALL AVER_SB(PDQBPB, PDBQPB)
         CALL AVER_SB(PDQBP , PDBQP )
      END IF

      ESBQBNLO = 
     +        +PDQBB * PREFB * DSBQBB(ALPHAS,SHATB,T1HATB,MS,MG,1) 
     +        +PDQBB * PREFB * DSBQBV(ALPHAS,SHATB,T1HATB,MS,MG,MT,1) 
     +        +PDQBB * PREFB * DSBQB1(ALPHAS,SHATB,T1HATB,MS,MG,SCA,1) 
     +        +PDQBB * PREFB * 
     +        DSBQBD(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX,1)
     +        +PDQBB * PREFB * 
     +        DSBQB2(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX,SCA,1)
     +        +PDQB  * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,T1HAT,S4,MS,MG,1,EPS) 
     +        +PDQB  * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA,1) 

      ESBQBNLO = ESBQBNLO
     +        +PDQBPB* PREFB * DSBQBB(ALPHAS,SHATB,T1HATB,MS,MG,0) 
     +        +PDQBPB* PREFB * DSBQBV(ALPHAS,SHATB,T1HATB,MS,MG,MT,0) 
     +        +PDQBPB* PREFB * DSBQB1(ALPHAS,SHATB,T1HATB,MS,MG,SCA,0) 
     +        +PDQBPB* PREFB * 
     +        DSBQBD(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX,0)
     +        +PDQBPB* PREFB * 
     +        DSBQB2(ALPHAS,SHATB,T1HATB,S4,MS,MG,DEL,S4MAX,SCA,0)
     +        +PDQBP  * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,T1HAT,S4,MS,MG,0,EPS) 
     +        +PDQBP  * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA,0) 

      ESBQBNLO = ESBQBNLO
     +        +PDBQB * PREFB * DSBQBB(ALPHAS,SHATB,U1HATB,MS,MG,1) 
     +        +PDBQB * PREFB * DSBQBV(ALPHAS,SHATB,U1HATB,MS,MG,MT,1) 
     +        +PDBQB * PREFB * DSBQB1(ALPHAS,SHATB,U1HATB,MS,MG,SCA,1) 
     +        +PDBQB * PREFB * 
     +        DSBQBD(ALPHAS,SHATB,U1HATB,S4,MS,MG,DEL,S4MAX,1)
     +        +PDBQB * PREFB * 
     +        DSBQB2(ALPHAS,SHATB,U1HATB,S4,MS,MG,DEL,S4MAX,SCA,1)
     +        +PDBQ  * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,U1HAT,S4,MS,MG,1,EPS) 
     +        +PDBQ  * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,U1HAT,S4,MS,MG,SCA,1) 

      ESBQBNLO = ESBQBNLO
     +        +PDBQPB* PREFB * DSBQBB(ALPHAS,SHATB,U1HATB,MS,MG,0) 
     +        +PDBQPB* PREFB * DSBQBV(ALPHAS,SHATB,U1HATB,MS,MG,MT,0) 
     +        +PDBQPB* PREFB * DSBQB1(ALPHAS,SHATB,U1HATB,MS,MG,SCA,0) 
     +        +PDBQPB* PREFB * 
     +        DSBQBD(ALPHAS,SHATB,U1HATB,S4,MS,MG,DEL,S4MAX,0)
     +        +PDBQPB* PREFB * 
     +        DSBQB2(ALPHAS,SHATB,U1HATB,S4,MS,MG,DEL,S4MAX,SCA,0)
     +        +PDBQP  * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,U1HAT,S4,MS,MG,0,EPS) 
     +        +PDBQP  * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,U1HAT,S4,MS,MG,SCA,0) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBGBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G QB AND G Q
      IMPLICIT REAL*8 (A-H,M-Z)

      ESBGBNLO = ESBGBNLOA() + ESBGBNLOB()
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBGBNLOA()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/PDNLSB2A/PDGB0,PDGB1, PDBG0,PDBG1
      COMMON/PDNLSB2B/PDGQ0,PDGQ1, PDQG0,PDQG1
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DINLSB3/SHAT0,T1HAT0,U1HAT0,S40,X10,X20,PREF0,S4G02
      COMMON/DINLSB4/SHAT1,T1HAT1,U1HAT1,S41,X11,X21,PREF1,S4UP
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN SQUARKS AND ANTISQUARKS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- SQUARKS
C***  ICHARCONJ = -1  <-- ANTISQUARKS
C***  ICHARCONJ =  0  <-- AVERAGE OF SQUARKS AND ANTISQUARKS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_SB(PDGB0, PDGQ0)
         CALL SWAP_SB(PDGB1, PDGQ1)
         CALL SWAP_SB(PDBG0, PDQG0)
         CALL SWAP_SB(PDBG1, PDQG1)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_SB(PDGB0, PDGQ0)
         CALL AVER_SB(PDGB1, PDGQ1)
         CALL AVER_SB(PDBG0, PDQG0)
         CALL AVER_SB(PDBG1, PDQG1)
      END IF

      IF ((SHAT1.GT.(MS+MG)**2).AND.(MG.GT.MS).AND.(S40.LE.S4UP)) 
     +     THEN
         ESBGBNLOA = 
     +        + PDGB0 * PREF0 * S4G02 *
     +        DSBGBH(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,EPS)
     +        + PDGB0 * PREF0 * S4G02 * 
     +        DSBGB3(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,SCA)
     +        + PDGB0 * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,EPS)
     +        - PDGB1 * PREF1 * 
     +        DSBGBS(ALPHAS,SHAT1,T1HAT1,S41,MS,MG,EPS)
     +        
     +        + PDBG0 * PREF0 * S4G02 *
     +        DSBGBH(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,EPS)
     +        + PDBG0 * PREF0 * S4G02 * 
     +        DSBGB3(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,SCA)
     +        + PDBG0 * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,EPS)
     +        - PDBG1 * PREF1 * 
     +        DSBGBS(ALPHAS,SHAT1,U1HAT1,S41,MS,MG,EPS)
      ELSE
         ESBGBNLOA = 
     +        + PDGB0 * PREF0 * S4G02 *
     +        DSBGBH(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,EPS)
     +        + PDGB0 * PREF0 * S4G02 * 
     +        DSBGB3(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,SCA)
     +        + PDGB0 * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT0,T1HAT0,S40,MS,MG,EPS)
     +        
     +        + PDBG0 * PREF0 * S4G02 *
     +        DSBGBH(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,EPS)
     +        + PDBG0 * PREF0 * S4G02 * 
     +        DSBGB3(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,SCA)
     +        + PDBG0 * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT0,U1HAT0,S40,MS,MG,EPS)
         
      END IF
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESBGBNLOB()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G Q
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/PDNLSB1/PDGG,PDQB,PDQBP ,PDBQ ,PDBQP,PDQG, PDGQ,PDBG, PDGB
      COMMON/PDNLSB3A/PDQG2,PDQG3,PDGQ2,PDGQ3 
      COMMON/PDNLSB3B/PDBG2,PDBG3,PDGB2,PDGB3 

      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/DINLSB5/SHAT2,T1HAT2,U1HAT2,S42,X12,X22,PREF2,S32,S3G2
      COMMON/DINLSB6/SHAT3,T1HAT3,U1HAT3,S43,X13,X23,PREF3,S3G3,S3UP
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN SQUARKS AND ANTISQUARKS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- SQUARKS
C***  ICHARCONJ = -1  <-- ANTISQUARKS
C***  ICHARCONJ =  0  <-- AVERAGE OF SQUARKS AND ANTISQUARKS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_SB(PDQG,   PDBG)
         CALL SWAP_SB(PDQG2,  PDBG2)
         CALL SWAP_SB(PDQG3,  PDBG3)
         CALL SWAP_SB(PDGQ,   PDGB)
         CALL SWAP_SB(PDGQ2,  PDGB2)
         CALL SWAP_SB(PDGQ3,  PDGB3)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_SB(PDQG,   PDBG)
         CALL AVER_SB(PDQG2,  PDBG2)
         CALL AVER_SB(PDQG3,  PDBG3)
         CALL AVER_SB(PDGQ,   PDGB)
         CALL AVER_SB(PDGQ2,  PDGB2)
         CALL AVER_SB(PDGQ3,  PDGB3)
      END IF

      IF ((SHAT3.GT.(MS+MG)**2).AND.(MG.GT.MS).AND.(S32.LE.S3UP)) 
     +     THEN
         ESBGBNLOB = 
     +        +PDQG * PREF *
     +        DSBQGH(ALPHAS,SHAT,T1HAT,S4,MS,MG,EPS)
     +        +PDQG * PREF *
     +        DSBQG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA)
     +        + PDQG2 * PREF2 *
     +        DSBQGT(ALPHAS,SHAT2,T1HAT2,S42,S3G2,MS,MG,EPS) 
     +        - PDQG3 * PREF3 *
     +        DSBQGT(ALPHAS,SHAT3,T1HAT3,S43,S3G3,MS,MG,EPS) 
     +        
     +        +PDGQ * PREF *
     +        DSBQGH(ALPHAS,SHAT,U1HAT,S4,MS,MG,EPS)
     +        +PDGQ * PREF *
     +        DSBQG3(ALPHAS,SHAT,U1HAT,S4,MS,MG,SCA)
     +        + PDGQ2 * PREF2 *
     +        DSBQGT(ALPHAS,SHAT2,U1HAT2,S42,S3G2,MS,MG,EPS) 
     +        - PDGQ3 * PREF3 *
     +        DSBQGT(ALPHAS,SHAT3,U1HAT3,S43,S3G3,MS,MG,EPS) 
      ELSE
         ESBGBNLOB = 
     +        +PDQG * PREF * 
     +        DSBQGH(ALPHAS,SHAT,T1HAT,S4,MS,MG,EPS)
     +        +PDQG * PREF * 
     +        DSBQG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA)
     +        + PDQG2 * PREF2 *
     +        DSBQGT(ALPHAS,SHAT2,T1HAT2,S42,S3G2,MS,MG,EPS)      
     +        
     +        +PDGQ * PREF *
     +        DSBQGH(ALPHAS,SHAT,U1HAT,S4,MS,MG,EPS)
     +        +PDGQ * PREF *
     +        DSBQG3(ALPHAS,SHAT,U1HAT,S4,MS,MG,SCA)
     +        + PDGQ2 * PREF2 *
     +        DSBQGT(ALPHAS,SHAT2,U1HAT2,S42,S3G2,MS,MG,EPS) 
      END IF
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSBGGNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PDNLSBT/PDGGT, PDQBT, PDQBPT, PDGQT
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSBGGNLO = 
     +        +PDGGT * PREFB * DSBGGB(ALPHAS,SHAT,T1HAT,MS,MG) 
     +        +PDGGT * PREFB * DSBGGV(ALPHAS,SHAT,T1HAT,MS,MG,MT) 
     +        +PDGGT * PREFB * DSBGG1(ALPHAS,SHAT,T1HAT,MS,MG,SCA) 
     +        +PDGGT * PREFB * 
     +        DSBGGD(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX)
     +        +PDGGT * PREFB * 
     +        DSBGG2(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX,SCA)
     +        +PDGGT * PREF  * DSBGGH(ALPHAS,SHAT,T1HAT,S4,MS,MG) 
     +        +PDGGT * PREF  * DSBGG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSBQBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDNLSBT/PDGGT, PDQBT, PDQBPT, PDGQT
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSBQBNLO = 
     +        +PDQBT * PREFB * DSBQBB(ALPHAS,SHAT,T1HAT,MS,MG,1) 
     +        +PDQBT * PREFB * DSBQBV(ALPHAS,SHAT,T1HAT,MS,MG,MT,1) 
     +        +PDQBT * PREFB * DSBQB1(ALPHAS,SHAT,T1HAT,MS,MG,SCA,1) 
     +        +PDQBT * PREFB * 
     +        DSBQBD(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX,1)
     +        +PDQBT * PREFB * 
     +        DSBQB2(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX,SCA,1)
     +        +PDQBT * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,T1HAT,S4,MS,MG,1,EPS) 
     +        +PDQBT * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA,1) 

      FSBQBNLO = FSBQBNLO
     +        +PDQBPT* PREFB * DSBQBB(ALPHAS,SHAT,T1HAT,MS,MG,0) 
     +        +PDQBPT* PREFB * DSBQBV(ALPHAS,SHAT,T1HAT,MS,MG,MT,0) 
     +        +PDQBPT* PREFB * DSBQB1(ALPHAS,SHAT,T1HAT,MS,MG,SCA,0) 
     +        +PDQBPT* PREFB * 
     +        DSBQBD(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX,0)
     +        +PDQBPT* PREFB * 
     +        DSBQB2(ALPHAS,SHAT,T1HAT,S4,MS,MG,DEL,S4MAX,SCA,0)
     +        +PDQBPT * PREF  * 
     +        DSBQBH(ALPHAS,SHAT,T1HAT,S4,MS,MG,0,EPS) 
     +        +PDQBPT * PREF  * 
     +        DSBQB3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA,0) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSBGBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G QB AND G Q 
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/PDNLSBT/PDGGT, PDQBT, PDQBPT, PDGQT
      COMMON/MINLSB/EPS,DEL,SCA,M2,GAM
      COMMON/DINLSB1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLSB2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/DINLSB3/SHAT0,T1HAT0,U1HAT0,S40,X10,X20,PREF0,S4G02
      COMMON/DINLSB4/SHAT1,T1HAT1,U1HAT1,S41,X11,X21,PREF1,S4UP

      IF ((MG.GT.MS).AND.(SHAT.GT.(MS+MG)**2)) THEN
         FSBGBNLO = 
     +        +PDGQT * PREF *
     +        DSBQGH(ALPHAS,SHAT,T1HAT,S4,MS,MG,EPS)
     +        +PDGQT * PREF *
     +        DSBQG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA)
     +        +PDGQT * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT,T1HAT0,S40,MS,MG,EPS)
     +        -PDGQT * PREF1 * 
     +        DSBGBS(ALPHAS,SHAT,T1HAT1,S41,MS,MG,EPS)
      ELSE
         FSBGBNLO = 
     +        +PDGQT * PREF *
     +        DSBQGH(ALPHAS,SHAT,T1HAT,S4,MS,MG,EPS)
     +        +PDGQT * PREF *
     +        DSBQG3(ALPHAS,SHAT,T1HAT,S4,MS,MG,SCA)
     +        +PDGQT * PREF0 * 
     +        DSBGBS(ALPHAS,SHAT,T1HAT0,S40,MS,MG,EPS)
      END IF
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION SUMNLOSB(SIGNLOSB1,SIGNLOSB2,SIGNLOSB3,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN NLO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF (IFLAVOR.EQ.0) SUMNLOSB = SIGNLOSB1 +SIGNLOSB2 +SIGNLOSB3
      IF (IFLAVOR.EQ.1) SUMNLOSB = SIGNLOSB1
      IF (IFLAVOR.EQ.2) SUMNLOSB = SIGNLOSB2
      IF (IFLAVOR.EQ.3) SUMNLOSB = SIGNLOSB3
      IF (IFLAVOR.EQ.4) SUMNLOSB = SIGNLOSB1 +SIGNLOSB2
      IF (IFLAVOR.EQ.5) SUMNLOSB = SIGNLOSB3
      IF (.NOT.(ABS(SUMNLOSB).LT.1.D35)) THEN
         SUMNLOSB = 0.D0
         PRINT *,'NAN'
      END IF
      RETURN
      END
      
C ======================================================================

      SUBROUTINE PDFLOSB(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDQBB = 0D0
      PDQBPB= 0D0
      PDBQB = 0D0
      PDBQPB= 0D0
      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I) 
            DO 1000 J =1,5
               IF (I.NE.J) PDQBPB  = PDQBPB  + XPDF( I) * YPDF( J) 
               IF (I.NE.J) PDBQPB  = PDBQPB  + XPDF(-I) * YPDF(-J) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF( I) 
            DO 1010 J =1,5
               IF (I.NE.J) PDQBPB  = PDQBPB  + XPDF( I) * YPDF(-J) 
               IF (I.NE.J) PDBQPB  = PDBQPB  + XPDF(-I) * YPDF( J) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSBT(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDNLSBT/PDGGT, PDQBT, PDQBPT, PDGQT

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGT = XPDF(0) * YPDF(0)
      PDQBT = 0D0
      PDQBPT= 0D0
      PDGQT = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF( I) + XPDF(-I) * YPDF(-I) 
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF( 0) 
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF( 0) 
            DO 1000 J =1,5
               IF (I.NE.J) PDQBPT  = PDQBPT  
     +              + XPDF( I) * YPDF( J) + XPDF(-I) * YPDF(-J) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF(-I) + XPDF(-I) * YPDF(I) 
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF(0) 
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF(0) 
            DO 1010 J =1,5
               IF (I.NE.J) PDQBPT  = PDQBPT  
     +              + XPDF( I) * YPDF(-J) + XPDF(-I) * YPDF( J) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSB1(X1,X2,X2B)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      DIMENSION ZPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOSB/PDGGB,PDQBB,PDQBPB,PDBQB,PDBQPB
      COMMON/PDNLSB1/PDGG,PDQB,PDQBP ,PDBQ ,PDBQP,PDQG, PDGQ,PDBG, PDGB

      CALL PARTONDF(X1 ,SCALE,XPDF)
      CALL PARTONDF(X2B,SCALE,YPDF)
      CALL PARTONDF(X2 ,SCALE,ZPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDGG  = XPDF(0) * ZPDF(0)
      PDQBB = 0D0
      PDQB  = 0D0
      PDQBPB= 0D0
      PDQBP = 0D0
      PDBQB = 0D0
      PDBQ  = 0D0
      PDBQPB= 0D0
      PDBQP = 0D0
      PDQG  = 0D0
      PDGQ  = 0D0
      PDBG  = 0D0
      PDGB  = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I) 
            PDQB  = PDQB  + XPDF( I) * ZPDF( I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I) 
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF(-I) 
            PDQG  = PDQG  + XPDF(+I) * ZPDF( 0) 
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF(-I) 
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0) 
            PDGB  = PDGB  + XPDF( 0) * ZPDF( I) 
            
            DO 1000 J =1,5
               IF (I.NE.J) PDQBPB  = PDQBPB  + XPDF( I) * YPDF( J) 
               IF (I.NE.J) PDQBP   = PDQBP   + XPDF( I) * ZPDF( J) 
               IF (I.NE.J) PDBQPB  = PDBQPB  + XPDF(-I) * YPDF(-J) 
               IF (I.NE.J) PDBQP   = PDBQP   + XPDF(-I) * ZPDF(-J) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I) 
            PDQB  = PDQB  + XPDF( I) * ZPDF(-I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF( I) 
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF( I) 
            PDQG  = PDQG  + XPDF( I) * ZPDF( 0) 
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF( I) 
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0) 
            PDGB  = PDGB  + XPDF( 0) * ZPDF(-I) 

            DO 1010 J =1,5
               IF (I.NE.J) PDQBPB  = PDQBPB  + XPDF( I) * YPDF(-J) 
               IF (I.NE.J) PDQBP   = PDQBP   + XPDF( I) * ZPDF(-J) 
               IF (I.NE.J) PDBQPB  = PDBQPB  + XPDF(-I) * YPDF( J) 
               IF (I.NE.J) PDBQP   = PDBQP   + XPDF(-I) * ZPDF( J) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSB2(X10,X20,X11,X21)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION WPDF(-6:6)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      DIMENSION ZPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDNLSB2A/PDGB0,PDGB1, PDBG0,PDBG1
      COMMON/PDNLSB2B/PDGQ0,PDGQ1, PDQG0,PDQG1

      CALL PARTONDF(X10,SCALE,WPDF)
      CALL PARTONDF(X20,SCALE,XPDF)
      CALL PARTONDF(X11,SCALE,YPDF)
      CALL PARTONDF(X21,SCALE,ZPDF)

      PDGB0 = 0D0
      PDGB1 = 0D0
      PDBG0 = 0D0
      PDBG1 = 0D0
      PDGQ0 = 0D0
      PDGQ1 = 0D0
      PDQG0 = 0D0
      PDQG1 = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDGB0 = PDGB0 + WPDF( 0) * XPDF( I) 
            PDGB1 = PDGB1 + YPDF( 0) * ZPDF( I) 
            PDBG0 = PDBG0 + WPDF(-I) * XPDF( 0) 
            PDBG1 = PDBG1 + YPDF(-I) * ZPDF( 0) 
            PDGQ0 = PDGQ0 + WPDF( 0) * XPDF(-I) 
            PDGQ1 = PDGQ1 + YPDF( 0) * ZPDF(-I) 
            PDQG0 = PDQG0 + WPDF( I) * XPDF( 0) 
            PDQG1 = PDQG1 + YPDF( I) * ZPDF( 0) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDGB0 = PDGB0 + WPDF( 0) * XPDF(-I) 
            PDGB1 = PDGB1 + YPDF( 0) * ZPDF(-I) 
            PDBG0 = PDBG0 + WPDF(-I) * XPDF( 0) 
            PDBG1 = PDBG1 + YPDF(-I) * ZPDF( 0) 
            PDGQ0 = PDGQ0 + WPDF( 0) * XPDF( I) 
            PDGQ1 = PDGQ1 + YPDF( 0) * ZPDF( I) 
            PDQG0 = PDQG0 + WPDF( I) * XPDF( 0) 
            PDQG1 = PDQG1 + YPDF( I) * ZPDF( 0) 
 1010    CONTINUE
      END IF

      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSB3(X10,X20,X11,X21)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION WPDF(-6:6)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      DIMENSION ZPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDNLSB3A/PDQG2,PDQG3,PDGQ2,PDGQ3 
      COMMON/PDNLSB3B/PDBG2,PDBG3,PDGB2,PDGB3 

      CALL PARTONDF(X10,SCALE,WPDF)
      CALL PARTONDF(X20,SCALE,XPDF)
      CALL PARTONDF(X11,SCALE,YPDF)
      CALL PARTONDF(X21,SCALE,ZPDF)

      PDQG2 = 0D0
      PDQG3 = 0D0
      PDGQ2 = 0D0
      PDGQ3 = 0D0
      PDBG2 = 0D0
      PDBG3 = 0D0
      PDGB2 = 0D0
      PDGB3 = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQG2 = PDQG2 + WPDF(+I) * XPDF( 0) 
            PDQG3 = PDQG3 + YPDF(+I) * ZPDF( 0) 
            PDGQ2 = PDGQ2 + WPDF( 0) * XPDF(-I) 
            PDGQ3 = PDGQ3 + YPDF( 0) * ZPDF(-I) 
            PDBG2 = PDBG2 + WPDF(-I) * XPDF( 0) 
            PDBG3 = PDBG3 + YPDF(-I) * ZPDF( 0) 
            PDGB2 = PDGB2 + WPDF( 0) * XPDF( I) 
            PDGB3 = PDGB3 + YPDF( 0) * ZPDF( I) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQG2 = PDQG2 + WPDF(+I) * XPDF( 0) 
            PDQG3 = PDQG3 + YPDF(+I) * ZPDF( 0) 
            PDGQ2 = PDGQ2 + WPDF( 0) * XPDF(+I) 
            PDGQ3 = PDGQ3 + YPDF( 0) * ZPDF(+I) 
            PDBG2 = PDBG2 + WPDF(-I) * XPDF( 0) 
            PDBG3 = PDBG3 + YPDF(-I) * ZPDF( 0) 
            PDGB2 = PDGB2 + WPDF( 0) * XPDF(-I) 
            PDGB3 = PDGB3 + YPDF( 0) * ZPDF(-I) 
 1010    CONTINUE
      END IF

      RETURN
      END

C ======================================================================
