program makepomdel
      implicit none

C Initial data for evolution
C    Q0 = 1.6 GeV is starting point for evolution.
C    alam = Lambda_MSbar for CTEQ 4M
C    dknl = 1 to use 1-loop unpolarized kernel
C    dfl is number of flavors
C    dnxevl is number of points in x-grid.
C    ordn is the order in the strong coupling constant to be used.
C N.B. Use real*8, even for integer data, since that is the type needed by
C      the routine PARPDF used for setting parameters
      real*8 q0, alam, dfl, Xmin, qqmax, dknl, dnxevl, Xmax,ordn 
      integer nnq
      data q0, alam, dknl, dfl, dnxevl, Xmax,
     > Xmin, nnq, qqmax,ordn/1.60, 0.202, 1.0, 5.0, 80.00, 0.95
     > , 0.0001, 20.0, 10.00, 1.0/

C Parameters for initial pdf at Q0. Initial values set in the block data 
C statement below.
 

C  Called functions
      external PomZEUS2
      real*8 PomZEUS2

C  Local variables:
      integer iret
     
C      real*8 sigmapom, sigmaglue, sigmaquark, cg, cq

C Input for the number of points in x and Q

      print *, 'Enter: nx, nq: '
      read (*,*) dnxevl, nnq


C Set parameters for the initial pdfs (nfl etc).

      call parpdf (1, 'iknl', DKNL, iret)
      call parpdf (1, 'nx',   dnxevl, iret)
      call parpdf (1, 'QINI', Q0, iret)
      call parpdf (1, 'NFL',  dfl, iret)
      call parpdf (1, 'LAM',  alam, iret)
      call parpdf (1, 'Xmin', Xmin, iret)
      call parpdf (1, 'Qmax', qqmax, iret)
      call parqcd (1, 'ORDR', ordn, iret)
      call parqcd (1, 'NFL', dfl, iret)
  
C Do the evolution:
C PomZEUS2 is the initial distribution.
      call evolve (PomZEUS2, iret)
      if (iret .ne. 0) then
         print *, 'Error in evolution, IRET = ', iret, 
     >            '. No file written.'
         goto 90
      endif

C Save into file tmppom.dat
       call testpdfpom(10,Xmin,Xmax,dnxevl,Q0,qqmax,nnq)
C     All purpose exit:
 90   continue

      end

C=======================================================
C=======================================================
C writes the result of the evolution in a data file
C In the format first line Q and then the values for x, gluon and the sum
C of light quarks.
 
      subroutine testpdfpom (iset,
     >     xmin, xmax, nx, qmin, qmax, nq)
      implicit none     
      integer iset, nq, iret, ix, iq, nnx
      double precision xmin, xmax, qmin, qmax, u, d, s,
     >                 nx,glue, qall, xglue, ubar, dbar,
     >                 sbar, delx
      double precision x, q, pdf
      character*78 filename
      external pdf
      nnx = nx/1.
      delx = 0.0
C  Name of output file
      filename = 'tmppom.dat'
     
      Open(2,File = filename, Form = 'Formatted',status='New'
     >              ,iostat=iret)
      do 30 iq = 0, nq-1
            q = qmin + (qmax - qmin) * real(iq) / (nq-1)
            write(2,'(2F11.3)') q
      do 20 ix = 0, nnx-1
            x = xmin + (xmax - xmin) * real(ix) / (nx-1)
     >      - delx/2. 
C
C The 1 after iset specifies that one is interested in the distributions within
C a proton. 2 is for neutrons. See the pdf subroutine in pdfcpac3.f for a 
C complete list of this flag.

            u    = pdf(iset, 1, 1, x, q, iret)
            d    = pdf(iset, 1, 2, x, q, iret)
            s    = pdf(iset, 1, 3, x, q, iret)
            glue = pdf(iset, 1, 0, x, q, iret)
            ubar = pdf(iset, 1,-1, x, q, iret)
            dbar = pdf(iset, 1,-2, x, q, iret)
            sbar = pdf(iset, 1,-3, x, q, iret)
            if (iret .ne. 0) then
               write (*,*) 'Error code in pdf is ', iret
            endif
            qall = (u + d + s + ubar + dbar + sbar)
            xglue = glue       
            write(2,'(1F12.10,7E13.6)') x, xglue, qall
 20   continue
 30   continue

      close(2)

      end



C================================================================
C================================================================
C Specifies the initial parton distribution in the normalization point Q_0!
C Here we use CTEQ 4M.
C
      real*8 FUNCTION PomZEUS2 (LP, X)
C Initial values for parton densities
      Implicit none
      
      real*8 x
      integer lp
      real*8 result,del
      Data del/0.0/

      if ( (x .le. 0) .or. (x .ge. 1) .or. (abs(lp) .gt. 5) ) then
C        x out of range, or quark is top:
         result = 0
      else if ((abs(lp).eq.4) .or. (abs(lp).eq.5)) then
C        Charm or bottom quark is zero
         result = 0
      else if (lp .eq. 0) then
C        Gluon: 
         result = 1.123*(x-del/2.)**(-.206)
     >*(1-(x-del/2.))**(4.673)
     >*(1+4.269*(x-del/2.)**(1.508))
      
C     u-quark

       else if (lp. eq. 1) then

      result=1.344*(x-del/2.)**(.501)*(1-x+del/2.)**
     >3.689*(1+6.402*(x-del/2.)**(0.873))+
     >.5*(.255*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))-
     >.071*(x-del/2.)**(.501)*(1-x+del/2.)**8.041)

C     d-quark

      else if (lp. eq. 2) then 
         
      result=0.64*(x-del/2.)**(.501)*(1-x+del/2.)**
     >4.247*(1+2.69*(x-del/2.)**(0.333))+
     >.5*(.255*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))+
     >.071*(x-del/2.)**(.501)*(1-x+del/2.)**8.041)
     
      else if (lp. eq. -1) then
      
C  u_sea = anti-u in the normalization point

      result=.5*(.255*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))-
     >.071*(x-del/2.)**(.501)*(1-x+del/2.)**8.041)

      else if (lp. eq. -2) then
      
C  d_sea = anti-d in the normalization point
  
      result=.5*(.255*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))+
     >.071*(x-del/2.)**(.501)*(1-x+del/2.)**8.041)
      
      else if (lp. eq. 3) then
      
C  s-quark = s_sea in the normalization point

      result= .064*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))
    
      
      else if (lp. eq. -3) then

C  anti-s quark = anti-s_sea in the normalization point

      result= .064*(x-del/2.)**(-.143)*(1-x+del/2.)**8.041*
     >(1+6.112*(x-del/2.))
      endif
      
       PomZEUS2 = result
       End