Program evoldemo

C     A demo program to show how to generate CTEQ parton distributions from the
C    simple output file xxx.ini containing the basic QCD parameters and initial
C     parton shape parameters at Q = Qini.

C     The essential step is to call

C     Subroutine Evlini (Nini, Xmin,Qmax,Qmss, Nx,Nq, Lpr)

C   which (i) reads in the information in xxx.ini; (ii) performs QCD evolution
C according to the parameters in the xxx.ini file plus the input arguments, as
C     explained below; and (iii) gets ready for you to use the parton 
C     distributions.

C     Once Evlini is called, the function  

C                 ParDis (Iparton, X, Q)

C   will return the parton distribution function (probability distri.) with: 
C     Iparton = (6, 5, 4, 3, 2, 1, 0, -1, ......, -6)
C           for (t, b, c, s, d, u, g, u_bar, ..., t_bar)

C     X, Q are the usual x, Q.

C     ---------------------------------------------------------
C     This demo program does the following:

C  1. It reads in the arguments needed by Evlini from a file called evlini.inp;
C     (In your application, you can setup these parameters anyway
C      you like in the calling program.)

C  2  It opens the input xxx.ini file and associates it with Unit# = Nini;

C  3. It calls Evlini

C  4. It calls and prints out the values of the values of the parton 
C     distributions at some specific values of 
C     (X, Q) for d, u, g, u_bar, d_bar, strange.

C     In your application, step (1) can be replaced by simple assignment
C     statements in the calling program or be left the way it is; 
C     Step (2) & (3) are essential;
C     (The reason that (2) is not merged into (3) is because the
C     details of the "open" statement is occasionally machine dependent.)
C     Step (4), of course, is for demo only.  You replace it with your 
C     application.
C     ----------------------------------------------------------

C   Explanations of the arguments and functions of Evlini:

C * Read in xxx.ini file (Unit # = Nini) 
C   for input evolu parameters from fit xxx

C * Do evolution in the range (Xmin, 1) and (Qini, Qmax) 
C                                           (Qini is read in from xxx.ini)

C * Qmss(4:6) are heavy quark masses to be used for charm, bottom, and top.

C * Nx and Nq are # of grids in (x,Q) - tradeoff between accuracy/efficiency
C   Ranges           : 10 < Nx < 105 ;  3 < Nq < 21
C   Reasonable values: middle of the above ranges.

C * Set up Pdf(Iset, Ihadn, Iparton, x, Q, Irt)
C   Pdf's   to be used as Iset = 10

C * Lpr : print level -- for diagnostics and debugging


      IMPLICIT DOUBLE PRECISION (A-H, O-Z)

      CHARACTER FLIN*15, Flint*15

      Dimension Qmss(3), Df (-6:6, 20, 20), Ef (-6:6, 20, 20)

      DATA
     > IU1, IU3, Ihdn / 11, 12, 1 /
     > NNx, NNq /  3,  3 /

C      Dum = SetPdf()
C      Dum = SetQCD()
C                                        1  -----------------------
C                                           (necessary)
C     We need the following parameters to get started.  You can hardwire 
C     these parameters in your program or read them in from a file.
C     For simplicity, this demo program initializes these parameters inside 
C     the program.
C                       filename for QCD & initial parton shape parameters
      FlIn = 'cteq3m'
C                       Max value of Q needed. (Qmin = Q0 is fixed by input)
      Qmax = 1000.0
C                       Number of grid points in LnLnQ variable 
      Nq   = 10
C                       Min. value of X needed. (Xmax is, of course, 1.0)
      Xmin = 1.E-4
C                       Number of grid points in X
      Nx   = 80
C                       Print level control switch (for details to be printed)
      Lpr  = 2
C                       "heavy quark" mass thresholds
      Qmss(1) = 1.6d0
      Qmss(2) = 5.0d0
      Qmss(3) = 180.d0
C                                        2 ------------------------
C                                          (necessary)
      Print *,  'Opening xxx.ini file ......'
C                          Trim leading and trailing blanks and obtain length
      call trmstr(flin,len1)
C                                  (Flint, IU1) is the input Xxxx.ini file
      flint = flin(1:len1)//'.ini'
C                                       
      OPEN (IU1, FILE=FLINT, STATUS='OLD', ERR=901)
C                                        3 ------------------------
C                                          (necessary)
      Print *,                   'Calling Evlini ......'

      Call Evlini(Iu1, Xmin,Qmax,Qmss, Nx,Nq,Lpr)
      close(IU1)
C                                         4 -------------------------
C                                          demo applications section  
C                                     (to be replaced by user applications)
      Call ParPdf (2, 'Qini', Qini, Ir)

      Dx = (1.   - Xmin) / NNx
      Dq = (Qmax - Qini) / NNq

      Print '(A/2A)','Sample parton distributions:','   X       Q   ',
     > '     d         u         g       u_bar      d_bar      s'
      Do 11 Iq = 1, NNq
         Q = Qini + Dq * (Iq-1)
         Do 12 Ix = 1, NNx
            X = Xmin + Dx * (Ix-1)
            Do 13 Ip = 2, -3, -1
C CTEQ3 - evolved
               Df (Ip, Ix, Iq) = pdf (10, Ihdn, Ip, X, Q, Ir)
C CTEQ3 - parametrized
               Ef (Ip, Ix, Iq) = pdf ( 1, Ihdn, Ip, X, Q, Ir)
   13 Continue
      Print '(1X, F6.4, F8.3, 6(1pE10.3))',
     > X, Q, (Df(Ip, Ix, Iq), Ip=2,-3,-1)
      Print '(1X, F6.4, F8.3, 6(1pE10.3)/)',
     > X, Q, (Ef(Ip, Ix, Iq), Ip=2,-3,-1)
   12 Continue
   11 Continue
  99    stop
C
 901  Print *, 'Cannot open the file xxx.ini'
      Stop     'Stopped in Evlmak'

C                        ****************************
      END