 
C*********************************************************************
 
C...PYPDPR
C...Gives proton parton distributions according to a few different
C...parametrizations.
 
      SUBROUTINE PYPDPR(X,Q2,XPPR)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
C...Commonblocks.
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      SAVE /PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/
C...Arrays and data.
      DIMENSION XPPR(-6:6),Q2MIN(11)
      DATA Q2MIN/ 2.56D0, 2.56D0, 2.56D0, 0.4D0, 0.4D0, 0.4D0,
     &1.0D0, 1.0D0, 2*0D0,0.25D0/
 
C...Reset output array.
      DO 100 KFL=-6,6
        XPPR(KFL)=0D0
  100 CONTINUE
 
C...Common preliminaries.
      NSET=MAX(1,MIN(12,MSTP(51)))
      IF(NSET.EQ.9.OR.NSET.EQ.10) NSET=6
      VINT(231)=Q2MIN(NSET)
      IF(MSTP(57).EQ.0) THEN
        Q2L=Q2MIN(NSET)
      ELSE
        Q2L=MAX(Q2MIN(NSET),Q2)
      ENDIF
 
      IF(NSET.GE.1.AND.NSET.LE.3) THEN
C...Interface to the CTEQ 3 parton distributions.
        QRT=SQRT(MAX(1D0,Q2L))
 
C...Loop over flavours.
        DO 110 I=-6,6
          IF(I.LE.0) THEN
            XPPR(I)=PYCTEQ(NSET,I,X,QRT)
          ELSEIF(I.LE.2) THEN
            XPPR(I)=PYCTEQ(NSET,I,X,QRT)+XPPR(-I)
          ELSE
            XPPR(I)=XPPR(-I)
          ENDIF
  110   CONTINUE
 
      ELSEIF(NSET.GE.4.AND.NSET.LE.6) THEN
C...Interface to the GRV 94 distributions.
        IF(NSET.EQ.4) THEN
          CALL PYGRVL (X, Q2L, UV, DV, DEL, UDB, SB, CHM, BOT, GL)
        ELSEIF(NSET.EQ.5) THEN
          CALL PYGRVM (X, Q2L, UV, DV, DEL, UDB, SB, CHM, BOT, GL)
        ELSE
          CALL PYGRVD (X, Q2L, UV, DV, DEL, UDB, SB, CHM, BOT, GL)
        ENDIF
 
C...Put into output array.
        XPPR(0)=GL
        XPPR(-1)=0.5D0*(UDB+DEL)
        XPPR(-2)=0.5D0*(UDB-DEL)
        XPPR(-3)=SB
        XPPR(-4)=CHM
        XPPR(-5)=BOT
        XPPR(1)=DV+XPPR(-1)
        XPPR(2)=UV+XPPR(-2)
        XPPR(3)=SB
        XPPR(4)=CHM
        XPPR(5)=BOT
 
      ELSEIF(NSET.EQ.7) THEN
C...Interface to the CTEQ 5L parton distributions.
C...Range of validity 10^-6 < x < 1, 1 < Q < 10^4 extended by
C...freezing x*f(x,Q2) at borders. 
        QRT=SQRT(MAX(1D0,MIN(1D4,Q2L)))
        XIN=MAX(1D-6,MIN(1D0,X))
 
C...Loop over flavours (with u <-> d notation mismatch).
        SUMUDB=PYCT5L(-1,XIN,QRT)
        RATUDB=PYCT5L(-2,XIN,QRT)
        DO 120 I=-5,2
          IF(I.EQ.1) THEN
            XPPR(I)=XIN*PYCT5L(2,XIN,QRT)
          ELSEIF(I.EQ.2) THEN
            XPPR(I)=XIN*PYCT5L(1,XIN,QRT)
          ELSEIF(I.EQ.-1) THEN
            XPPR(I)=XIN*SUMUDB*RATUDB/(1D0+RATUDB)
          ELSEIF(I.EQ.-2) THEN
            XPPR(I)=XIN*SUMUDB/(1D0+RATUDB)
          ELSE
            XPPR(I)=XIN*PYCT5L(I,XIN,QRT)
            IF(I.LT.0) XPPR(-I)=XPPR(I)
          ENDIF
  120   CONTINUE
 
      ELSEIF(NSET.EQ.8) THEN
C...Interface to the CTEQ 5M1 parton distributions.
        QRT=SQRT(MAX(1D0,MIN(1D4,Q2L)))
        XIN=MAX(1D-6,MIN(1D0,X))
 
C...Loop over flavours (with u <-> d notation mismatch).
        SUMUDB=PYCT5M(-1,XIN,QRT)
        RATUDB=PYCT5M(-2,XIN,QRT)
        DO 130 I=-5,2
          IF(I.EQ.1) THEN
            XPPR(I)=XIN*PYCT5M(2,XIN,QRT)
          ELSEIF(I.EQ.2) THEN
            XPPR(I)=XIN*PYCT5M(1,XIN,QRT)
          ELSEIF(I.EQ.-1) THEN
            XPPR(I)=XIN*SUMUDB*RATUDB/(1D0+RATUDB)
          ELSEIF(I.EQ.-2) THEN
            XPPR(I)=XIN*SUMUDB/(1D0+RATUDB)
          ELSE
            XPPR(I)=XIN*PYCT5M(I,XIN,QRT)
            IF(I.LT.0) XPPR(-I)=XPPR(I)
          ENDIF
  130   CONTINUE
 
      ELSEIF(NSET.EQ.11) THEN
C...GRV 92 LO distributions (for backwards checks).
C...Allowed variable range: 0.25 GeV^2 < Q^2 < 10^8 GeV^2 and
C...10^-5 < x < 1.
 
C...Determine s expansion variable and some x expressions.
        IF(MSTP(57).LE.0) THEN
          SD=0D0
        ELSE
          Q2IN=MIN(1D8,MAX(0.25D0,Q2))
          SD=LOG(LOG(Q2IN/0.232D0**2)/LOG(0.25D0/0.232D0**2))
        ENDIF
        SD2=SD**2
        XL=-LOG(X)
        XS=SQRT(X)
 
C...Evaluate valence, gluon and sea distributions.
        XFVUD=(0.663D0+0.191D0*SD-0.041D0*SD2+0.031D0*SD**3)*
     &  X**0.326D0*(1D0+(-1.97D0+6.74D0*SD-1.96D0*SD2)*XS+
     &  (24.4D0-20.7D0*SD+4.08D0*SD2)*X)*
     &  (1D0-X)**(2.86D0+0.70D0*SD-0.02D0*SD2)
        XFVDD=(0.579D0+0.283D0*SD+0.047D0*SD2)*X**(0.523D0-0.015D0*SD)*
     &  (1D0+(2.22D0-0.59D0*SD-0.27D0*SD2)*XS+(5.95D0-6.19D0*SD+
     &  1.55D0*SD2)*X)*(1D0-X)**(3.57D0+0.94D0*SD-0.16D0*SD2)
        XFGLU=(X**(1.00D0-0.17D0*SD)*((4.879D0*SD-1.383D0*SD2)+
     &  (25.92D0-28.97D0*SD+5.596D0*SD2)*X+(-25.69D0+23.68D0*SD-
     &  1.975D0*SD2)*X**2)+SD**0.558D0*EXP(-(0.595D0+2.138D0*SD)+
     &  SQRT(4.066D0*SD**1.218D0*XL)))*
     &  (1D0-X)**(2.537D0+1.718D0*SD+0.353D0*SD2)
        XFSEA=(X**(0.412D0-0.171D0*SD)*(0.363D0-1.196D0*X+(1.029D0+
     &  1.785D0*SD-0.459D0*SD2)*X**2)*XL**(0.566D0-0.496D0*SD)+
     &  SD**1.396D0*EXP(-(3.838D0+1.944D0*SD)+SQRT(2.845D0*SD**1.331D0*
     &  XL)))*(1D0-X)**(4.696D0+2.109D0*SD)
        XFSTR=SD**0.803D0*(1D0+(-3.055D0+1.024D0*SD**0.67D0)*XS+
     &  (27.4D0-20.0D0*SD**0.154D0)*X)*(1D0-X)**6.22D0*
     &  EXP(-(4.33D0+1.408D0*SD)+SQRT((8.27D0-0.437D0*SD)*
     &  SD**0.563D0*XL))/XL**(2.082D0-0.577D0*SD)
        IF(SD.LE.0.888D0) THEN
          XFCHM=0D0
        ELSE
          XFCHM=(SD-0.888D0)**1.01D0*(1.+(4.24D0-0.804D0*SD)*X)*
     &    (1D0-X)**(3.46D0+1.076D0*SD)*EXP(-(4.61D0+1.49D0*SD)+
     &    SQRT((2.555D0+1.961D0*SD)*SD**0.37D0*XL))
        ENDIF
        IF(SD.LE.1.351D0) THEN
          XFBOT=0D0
        ELSE
          XFBOT=(SD-1.351D0)*(1D0+1.848D0*X)*(1D0-X)**(2.929D0+
     &    1.396D0*SD)*EXP(-(4.71D0+1.514D0*SD)+
     &    SQRT((4.02D0+1.239D0*SD)*SD**0.51D0*XL))
        ENDIF
 
C...Put into output array.
        XPPR(0)=XFGLU
        XPPR(1)=XFVDD+XFSEA
        XPPR(2)=XFVUD-XFVDD+XFSEA
        XPPR(3)=XFSTR
        XPPR(4)=XFCHM
        XPPR(5)=XFBOT
        XPPR(-1)=XFSEA
        XPPR(-2)=XFSEA
        XPPR(-3)=XFSTR
        XPPR(-4)=XFCHM
        XPPR(-5)=XFBOT
 
C...Symmetric choice for debugging only
      ELSEIF(NSET.EQ.12) THEN
        XPPR(0)=.5D0/X
        XPPR(1)=.05D0/X
        XPPR(2)=.05D0/X
        XPPR(3)=.05D0/X
        XPPR(4)=.05D0/X
        XPPR(5)=.05D0/X
        XPPR(-1)=.05D0/X
        XPPR(-2)=.05D0/X
        XPPR(-3)=.05D0/X
        XPPR(-4)=.05D0/X
        XPPR(-5)=.05D0/X
 
      ENDIF
 
      RETURN
      END
