C23456789012345678901234567890123456789012345678901234567890123456789012
C
C S. D. Ellis and D. E. Soper, 23 August 1994
C Modified 9 May 1996 (PA & DES)
C Small modification 4 April 1997 (PA & DES) format G12.5 -> G13.5 
C
C    ROUTINE TO CALCULATE PARTON DATA TABLE FOR INTERPOLATION
C    F(NPARTON,Nx,Mmu) = DLOG( pdf + 1.0D-16 )
C
      PROGRAM MAKETABLE
C
      INTEGER VERSION
      PARAMETER(VERSION = 110496)
C
      INTEGER NIN,NOUT
C
      INTEGER NFL
      DOUBLE PRECISION X,LX,LXMIN,LXMAX,DELTALX
      DOUBLE PRECISION MU,LMU,LMUMIN,LMUMAX,DELTALMU
C
      DOUBLE PRECISION F(-2:6,256,256)
C
      CHARACTER*20 FILENAME
      CHARACTER*80 HEADER
      DOUBLE PRECISION LAMBDA
      INTEGER LAMBDAFLAVORS
      DOUBLE PRECISION XMIN,XMAX
      INTEGER NX
      DOUBLE PRECISION MUMIN,MUMAX
      INTEGER NMU
      INTEGER NPARTON
C
C Loop variables
C
      INTEGER N,M
C
      DOUBLE PRECISION PARTONIN,TRIAL
C
      NIN = 5
      NOUT = 6
C
      X = 0.1D0
      TRIAL = PARTONIN(0,X,100.0D0)
      WRITE(NOUT,*)'GLUON AT X = 0.1 AND MU = 100 IS',TRIAL
C
      TRIAL = 0.0D0
      DO N = 1,100
        X = N/100.0D0
        TRIAL = TRIAL + 0.01 * X * PARTONIN(0,X,100.0D0)
      ENDDO
      WRITE(NOUT,*)'GLUON MOMENTUM IS',TRIAL
C
C We use NFL for the number of quark flavors in the parton distribution set.
C We use LAMBDAFLAVORS for the number of quark flavors in the definition
C of the MSbar scale lambda in alpha_s.
C
      NFL = 6
      LAMBDAFLAVORS = 5
      WRITE(NOUT,501)LAMBDAFLAVORS
501   FORMAT('Please give MS-bar lambda for',I2,' flavors.')
      READ(NIN,*) LAMBDA
      WRITE(NOUT,502)LAMBDAFLAVORS,LAMBDA
502   FORMAT('LAMBDA for',I2,' flavors is',F10.4)
C
      XMIN =     0.1000000000D-04
      XMAX =     0.9500000000    
      NX   =         32
      MUMIN =      DSQRT(5.0D0)    
      MUMAX =      DSQRT(1310720.D0) 
      NMU =         32
C
      WRITE(NOUT,*)'Give file name for table'
      READ(NIN,98)FILENAME
98    FORMAT(A20)
      OPEN(UNIT = 99, FILE = FILENAME, STATUS = 'NEW')
      WRITE(NOUT,*)'Give header for table (up to 79 characters)'
      READ(NIN,99)HEADER
99    FORMAT(A80)
      HEADER = '!'//HEADER
      WRITE(99,99)HEADER
C                   12345678901234567890123456789012345678901234567890
      WRITE(99,100)'! This is a comment line, which could be modified.'
100   FORMAT(A50)
      WRITE(99,101)'PARTON DATA'
101   FORMAT(A11)
      WRITE(99,102)'   VERSION',' = ', VERSION
102   FORMAT(A10,A3,I12)
      WRITE(99,103)LAMBDA
103   FORMAT('    LAMBDA',' = ',G20.10)
      WRITE(99,104)LAMBDAFLAVORS
104   FORMAT('NFL_LAMBDA',' = ',I10)
      WRITE(99,105)XMIN
105   FORMAT('      XMIN',' = ',G20.10)
      WRITE(99,106)XMAX
106   FORMAT('      XMAX',' = ',G20.10)
      WRITE(99,107)NX
107   FORMAT(' N_XPOINTS',' = ',I10)
      WRITE(99,108)MUMIN
108   FORMAT('     MUMIN',' = ',G20.10)
      WRITE(99,109)MUMAX
109   FORMAT('     MUMAX',' = ',G20.10)
      WRITE(99,110)NMU
110   FORMAT('N_MUPOINTS',' = ',I10)
C
C
C Calculate lattice info.
C
      LXMIN  = DLOG(XMIN/(1.0D0 - XMIN))
      LXMAX  = DLOG(XMAX/(1.0D0 - XMAX))
      DELTALX = (LXMAX - LXMIN)/(NX-1)
      LMUMIN = DLOG(MUMIN)
      LMUMAX = DLOG(MUMAX)
      DELTALMU = (LMUMAX - LMUMIN)/(NMU-1)
C
      LX = LXMIN - DELTALX
      DO 20 N=1,NX 
        LX = LX + DELTALX
        X = DEXP(LX)/(1.0D0+DEXP(LX))
        LMU = LMUMIN - DELTALMU    
        DO 20 M=1,NMU
          LMU = LMU + DELTALMU
          MU = DEXP(LMU)
C
          DO NPARTON = -2,NFL
            TRIAL = PARTONIN(NPARTON,X,MU) 
            IF(TRIAL.LT.0.0D0) THEN
             WRITE(NOUT,*)'OOPS',NPARTON,X,MU,TRIAL
             TRIAL = 0.0D0
            ENDIF
            F(NPARTON,N,M) =  DLOG( TRIAL + 1.0D-16 )
          ENDDO
C
          DO NPARTON = NFL+1,6
            TRIAL = 0.0D0 
            F(NPARTON,N,M) =  DLOG( TRIAL + 1.0D-16 )
          ENDDO
C
 20   CONTINUE
C
      DO 30 N=1,NX
      DO 30 M=1,NMU
      WRITE(99,50)F(-2,N,M),F(-1,N,M),F(0,N,M),F(1,N,M),
     >            F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
 50   FORMAT(9G13.5)
C
 30   CONTINUE

C
      STOP
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      DOUBLE PRECISION FUNCTION PARTONIN(NPARTON,X,MU)
C
      INTEGER NPARTON
      DOUBLE PRECISION X,MU
C
C This is an example of a function to link to the CTEQ3M code provided
C by the CTEQ authors.
C
      INTEGER ISET,IRT
      DOUBLE PRECISION Ctq3Pd
C
      ISET = 1
      PARTONIN = Ctq3Pd(ISET, NPARTON, X, MU, IRT)
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012