C***********************************************************************
C* This file contains the GRV parametrizations for the parton content  *
C* of the photon. Both leading and higher order versions are included. *
C* Notice that these parametrizations are valid down to Q^2 around 1/4 *
C* GeV^2. For the HO parametrization, separate functions for the total *
C* and leading contributions are given.                                *
C* As usual, double precision is used, and alpha_em has been factored  *
C* out.                                                                *
C***********************************************************************
      SUBROUTINE GRV
C  ***  The usual dummy subroutine.  ***
      WRITE(*,*) ' You are using the GRV parametrization!'
      RETURN
      END
 
      REAL*8 FUNCTION FL(X,S)
C******************************************************************
C* Functional form of the distribution function of light partons. *
C* S = log( log(Q**2/.232^2) / log(.25/.232^2) ) (LO),            *
C* S = log( log(Q**2/.248^2) / log(.3/.248^2) )  (HO).            *
C******************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      FL = ( X**A * (BA+BB*DSQRT(X)+C*X**B) +
     &       S**ALP * DEXP(-E+DSQRT(-EP*S**BET*DLOG(X))) ) *
     &     (1.D0-X)**D/X
      RETURN
      END
 
 
      REAL*8 FUNCTION FH(X,S)
C*************************************************
C* Same as FL, but for the heavy s, c, b quarks. *
C*************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      FH = ( (S-SP) * X**A * (BA+BB*DSQRT(X)+C*X**B) +
     &       (S-SP)**ALP * DEXP(-E+DSQRT(-EP*S**BET*DLOG(X))) ) *
     &     (1.D0-X)**D/X
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVUL(X,Q2)
C**********************************************************************
C* Leading order up-quark distributions. X is Bjorken-x, and Q2 is    *
C* the factorization scale. Recall that alpha_em has been factored    *
C* out. The parametrization is supposed to be valid for x > 10**-5,   *
C* Q**2 < 10**6 GeV**2.                                               *
C**********************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = 1.717D0
      BET = .641D0
      A = .5D0 - .176D0*S
      B = 15.D0 - 5.687D0*RTS - .552D0*S2
      BA = .235D0 + .046D0*RTS
      BB = .082D0 - .051D0*S + .168D0*S2
      C = .459D0*S
      D = .354D0 - .061D0*S
      E = 4.899D0 + 1.678D0*S
      EP = 2.046D0 + 1.389D0*S
      GRVUL = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVDL(X,Q2)
C************************************
C* Same as GRVUL, but for d-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = 1.549D0
      BET = .782D0
      A = .496D0 + .026D0*S
      B = .685D0 - .58D0*RTS + .608D0*S2
      BA = .233D0 + .302D0*S
      BB = -.818D0*S + .198D0*S2
      C = .114D0 + .154D0*S
      D = .405D0 - .195D0*S + .046D0*S2
      E = 4.807D0 + 1.226D0*S
      EP = 2.166D0 + .664D0*S
      GRVDL = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVGL(X,Q2)
C**********************************
C* Same as GRVUL, but for gluons. *
C**********************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = .676D0
      BET = 1.089D0
      A = .462D0 - .524D0*RTS
      B = 5.451D0 - .804D0*S2
      BA = .535D0 - .504D0*RTS + .288D0*S2
      BB = .364D0 - .52D0*S
      C = -.323D0 + .115D0*S2
      D = .233D0 + .79D0*S - .139D0*S2
      E = .893D0 + 1.968D0*S
      EP = 3.432D0 + .392D0*S
      GRVGL = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVSL(X,Q2)
C************************************
C* Same as GRVUL, but for s-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      RTS = DSQRT(S)
      S2 = S*S
      SP = 0.D0
      ALP = 1.609D0
      BET = .962D0
      A = .47D0 - .099D0*S2
      B = 3.246D0
      BA = .121D0 - .068D0*RTS
      BB = -.09D0 + .074D0*S
      C = .062D0 + .034D0*S
      D = .226D0*S - .06D0*S2
      E = 4.288D0 + 1.707D0*S
      EP = 2.122D0 + .656D0*S
      GRVSL = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVCL(X,Q2)
C************************************
C* Same as GRVUL, but for c-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      S2 = S*S
      SP = .888D0
      IF(S.LT.SP) THEN
	GRVCL = 0.D0
	RETURN
      ENDIF
      ALP = .97D0
      BET = .545D0
      A = 1.254D0 - .251D0*S
      B = 3.932D0 - .327D0*S2
      BA = .658D0 + .202D0*S
      BB = -.699D0
      C = .965D0
      D = .141D0*S - .027D0*S2
      E = 4.911D0 + .969D0*S
      EP = 2.796D0 + .952D0*S
      GRVCL = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVBL(X,Q2)
C************************************
C* Same as GRVUL, but for b-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG( DLOG(18.58D0*Q2)/1.536D0 )
      S2 = S*S
      SP = 1.351D0
      IF(S.LT.SP) THEN
	GRVBL = 0.D0
	RETURN
      ENDIF
      ALP = 1.016D0
      BET = .338D0
      A = 1.961D0 - .37D0*S
      B = .923D0 + .119D0*S
      BA = .815D0 + .207D0*S
      BB = -2.275D0
      C = 1.48D0
      D = -.223D0 + .173D0*S
      E = 5.426D0 + .623D0*S
      EP = 3.819D0 + .901D0*S
      GRVBL = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVUH(X,Q2)
C**********************************************************************
C* Higher order up-quark distributions. X is Bjorken-x, and Q2 is the *
C* factorization scale. Recall that alpha_em has been factored out.   *
C* The parametrization is supposed to be  valid for x > 10**-5,       *
C* Q**2 < 10**6 GeV**2.                                               *
C* WARNING: These distributions are in the DIS_gamma scheme; to trans-*
C* late them into MSbar, a term proportional to B_gamma must be added.*
C**********************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = .583D0
      BET = .688D0
      A = .449D0 - .025D0*S - .071D0*S2
      B = 5.06D0 - 1.116D0*RTS
      BA = .103D0
      BB = .319D0 + .422D0*S
      C = 1.508D0 + 4.792D0*S - 1.963D0*S2
      D = 1.075D0 +.222D0*RTS - .193D0*S2
      E = 4.147D0 + 1.131D0*S
      EP = 1.661D0 + .874D0*S
      GRVUH = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVDH(X,Q2)
C************************************
C* Same as GRVUH, but for d-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = .591D0
      BET = .698D0
      A = .442D0 - .132D0*S - .058D0*S2
      B = 5.437D0 - 1.916D0*RTS
      BA = .099D0
      BB = .311D0 - .059D0*S
      C = .8D0 + .078D0*S - .1D0*S2
      D = .862D0 +.294D0*RTS - .184D0*S2
      E = 4.202D0 + 1.352D0*S
      EP = 1.841D0 + .99D0*S
      GRVDH = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVGH(X,Q2)
C****************************************************************
C* Same as GRVUH, but for gluons. Note that DIS_gamma and MSbar *
C* scheme are identical for gluons (to the given order).        *
C****************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = 1.161D0
      BET = 1.591D0
      A = .53D0 - .742D0*RTS + .025D0*S2
      B = 5.662D0
      BA = .533D0 - .281D0*RTS + .218D0*S2
      BB = .025D0 - .518D0*S + .156D0*S2
      C = -.282D0 + .209D0*S2
      D = .107D0 +1.058D0*S - .218D0*S2
      E = 2.704D0*S
      EP = 3.071D0 - .378D0*S
      GRVGH = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVSH(X,Q2)
C************************************
C* Same as GRVUH, but for s-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      SP = 0.D0
      ALP = .635D0
      BET = .456D0
      A = 1.77D0 - .735D0*RTS - .079D0*S2
      B = 3.832D0
      BA = .084D0 -.023D0*S
      BB = .136D0
      C = 2.119D0 - .942D0*S + .063D0*S2
      D = 1.271D0*S + .076D0*S - .19D0*S2
      E = 4.604D0 + .737D0*S
      EP = 1.641D0 + .976D0*S
      GRVSH = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVCH(X,Q2)
C************************************
C* Same as GRVUH, but for c-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      S2 = S*S
      SP = .82D0
      ALP = .926D0
      BET = .152D0
      A = 1.142D0 - .175D0*S
      B = 3.276D0
      BA = .504D0 +.317D0*S
      BB = -.433D0
      C = 3.334D0
      D = .398D0 + .326D0*S - .107D0*S2
      E = 5.493D0 + .408D0*S
      EP = 2.426D0 + 1.277D0*S
      GRVCH = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVBH(X,Q2)
C************************************
C* Same as GRVUH, but for b-quarks. *
C************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      S2 = S*S
      SP = 1.297D0
      ALP = .969D0
      BET = .266D0
      A = 1.953D0 - .391D0*S
      B = 1.657D0 - .161D0*S
      BA = 1.076D0 + .034D0*S
      BB = -2.015D0
      C = 1.662D0
      D = .353D0 + .016D0*S
      E = 5.713D0 + .249D0*S
      EP = 3.456D0 + .673D0*S
      GRVBH = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVUH0(X,Q2)
C**********************************************************************
C* The leading order part of the total HO up-quark distribution. This *
C* part is needed since, strictly speaking, only this part should be  *
C* multiplied with the HO part of a NLO cross section or Wilson co-   *
C* efficient. For the meaning of the variables, see GRVUH.            *
C**********************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = 1.447D0
      BET = .848D0
      A = .527D0 + .2D0*S - .107D0*S2
      B = 7.106D0 - .31D0*RTS - .786D0*S2
      BA = .197D0 + .533D0*S
      BB = .062D0 - .398D0*S + .109D0*S2
      C = .755D0*S - .112D0*S2
      D = .318D0 - .059D0*S
      E = 4.225D0 + 1.708D0*S
      EP = 1.752D0 + .866D0*S
      GRVUH0 = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVDH0(X,Q2)
C*************************************
C* Same as GRVUH0, but for d-quarks. *
C*************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = 1.424D0
      BET = .77D0
      A = .5D0 + .067D0*RTS - .055D0*S2
      B = .376D0 - .453D0*RTS + .405D0*S2
      BA = .156D0 + .184D0*S
      BB = -.528D0*S + .146D0*S2
      C = .121D0 + .092D0*S
      D = .379D0 - .301D0*S + .081D0*S2
      E = 4.346D0 + 1.638D0*S
      EP = 1.645D0 + 1.016D0*S
      GRVDH0 = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVGH0(X,Q2)
C*****************************************************************
C* Same as GRVUH0, but for gluons. Note that DIS_gamma and MSbar *
C* scheme are identical for gluons (to the given order).         *
C*****************************************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      ALP = .661D0
      BET = .793D0
      A = .537D0 - .6D0*RTS
      B = 6.389D0 - .953D0*S2
      BA = .558D0 - .383D0*RTS + .261D0*S2
      BB = -.305D0*S
      C = -.222D0 + .078D0*S2
      D = .153D0 +.978D0*S - .209D0*S2
      E = 1.429D0 + 1.772D0*S
      EP = 3.331D0 + .806D0*S
      GRVGH0 = FL(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVSH0(X,Q2)
C*************************************
C* Same as GRVUH0, but for s-quarks. *
C*************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      RTS = DSQRT(S)
      S2 = S*S
      SP = 0.D0
      ALP = 1.578D0
      BET = .863D0
      A = .622D0 + .332D0*S - .3D0*S2
      B = 2.469D0
      BA = .211D0 -.064D0*RTS - .018D0*S2
      BB = -.215D0 + .122D0*S
      C = .153D0
      D = .253D0*S - .081D0*S2
      E = 3.99D0 + 2.014D0*S
      EP = 1.72D0 + .986D0*S
      GRVSH0 = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVCH0(X,Q2)
C*************************************
C* Same as GRVUH0, but for c-quarks. *
C*************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      S2 = S*S
      SP = .82D0
      ALP = .929D0
      BET = .381D0
      A = 1.228D0 - .231D0*S
      B = 3.806D0 - .337D0*S2
      BA = .932D0 +.15D0*S
      BB = -.906D0
      C = 1.133D0
      D = .138D0*S - .028D0*S2
      E = 5.588D0 + .628D0*S
      EP = 2.665D0 + 1.054D0*S
      GRVCH0 = FH(X,S)
      RETURN
      END
 
 
      REAL*8 FUNCTION GRVBH0(X,Q2)
C*************************************
C* Same as GRVUH0, but for b-quarks. *
C*************************************
      IMPLICIT REAL*8 (A-H,M,O-Z)
      COMMON /GRVPAR/ ALP,BET,A,B,BA,BB,C,D,E,EP,SP
      S = DLOG ( DLOG(Q2*16.26D0)/1.585D0 )
      S2 = S*S
      SP = 1.297D0
      ALP = .97D0
      BET = .207D0
      A = 1.719D0 - .292D0*S
      B = .928D0 + .096D0*S
      BA = .845D0 +.178D0*S
      BB = -2.31D0
      C = 1.558D0
      D = -.191D0 + .151D0*S
      E = 6.089D0 + .282D0*S
      EP = 3.379D0 + 1.062D0*S
      GRVBH0 = FH(X,S)
      RETURN
      END