C     ------------------------------------------------------------------
      program Toolbox01
C     ------------------------------------------------------------------
C--   QCDNUM toolbox example #1
C--   Shows how to partition the local workspace w(nw).
C--   The MAKETAB error messages will indicate how large nw should be.
C--   Disk dump and read of the (yet unfilled) weight tables.
C--
C--   QCDNUM  routines are in capitals
C--   Package routines are in mixed case
C     ------------------------------------------------------------------
      implicit double precision (a-h,o-z)

C--   ------------------------------------------------------------------
C--   QCDNUM declarations and parameter settings
C--   ------------------------------------------------------------------
C--   x-grid
      dimension xmin(5), iwt(5)
      data nxmi/5/, nxin/100/, iosp/3/
      data xmin/1.D-5,0.2D0,0.4D0,0.6D0,0.75D0/
      data iwt / 1, 2, 4, 8, 16 /
C--   mu2-grid
      dimension qq(2),wt(2)
      data qq/2.D0,1.D4/, wt/1.D0,1.D0/, nqin/60/
C--   FFNS or VFNS and flavour thresholds
      data nfin/3/, q2c/3.D0/, q2b/25.D0/, q2t/1.D11/
C--   Input alphas
      data as0/0.364/, r20/2.D0/
C--   LO/NLO/NNLO,  unpol/pol/frag and starting scale
      data iord/1/, q0/2.D0/
C--   Renormalisation scale and factorisation scale
      data aar/1.D0/, bbr/0.D0/, aaq/1.D0/, bbq/0.D0/
C--   Input parton distributions and flavour composition of these
      external func
      dimension def(-6:6,12)
      data def  /
C--   tb  bb  cb  sb  ub  db   g   d   u   s   c   b   t
C--   -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6  
     + 0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,   !dval
     + 0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,   !uval
     + 0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,   !sval
     + 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,   !dbar
     + 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,   !ubar
     + 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !sbar
     + 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,   !cval
     + 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !cbar
     + 0.,-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,   !bval
     + 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !bbar
     +-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,   !tval
     + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. /  !tbar
C--   Weight files
      character*24 fnamq(2)
C--               123456789012345678901234
      data fnamq /'../weights/QcdNumUnp.wgt',
     +            '../weights/QcdNumPol.wgt' /


C--   ------------------------------------------------------------------
C--   Put here your own declarations and data statements
C--   ------------------------------------------------------------------
      parameter( nw = 1 )   !The MAKETAB error messages will guide you
      dimension w(nw)
      dimension itypes_alf(6),itypes_pdf(6)
      data itypes_alf / 0, 0, 0, 0, 0, 1/
      data itypes_pdf / 0, 0, 0, 0,13, 0/
      character*12 key
      data key /'MyAddOn v1.0'/
      character*24 fnamu(2)
      data fnamu /'../weights/QcdUsrUnp.wgt',
     +            '../weights/QcdUsrPol.wgt' /

C--   ------------------------------------------------------------------
C--   QCDNUM set-up and pdf evolution
C--   ------------------------------------------------------------------

      lun    = 6            !lun = 6 stdout, -6 stdout w/out banner page
      lunout = abs(lun)      
      call QCINIT(lun,' ')
      call GXMAKE(xmin,iwt,nxmi,nxin,nx,iosp)
      call GQMAKE(qq,wt,2,nqin,nq)
      lun = NXTLUN(0)
      do iset = 1,2
        call READWT(lun,fnamq(iset),idmin,idmax,nwpdf,ierr)
        if(ierr.ne.0) then
          call FILLWT(iset,idmin,idmax,nwpdf)
          call DMPWGT(iset,lun,fnamq(iset))
        endif
      enddo
      call SETORD(iord)                         !LO, NLO, NNLO
      call SETALF(as0,r20)                      !input alphas
      iqc  = IQFRMQ(q2c)                        !charm threshold
      iqb  = IQFRMQ(q2b)                        !bottom threshold
      call SETCBT(nfin,iqc,iqb,999)             !thesholds in the vfns
      call SETABR(aar,bbr)                      !renormalisation scale
      iq0  = IQFRMQ(q0)                         !starting scale
      do iset = 1,2
        call EVOLFG(iset,func,def,iq0,eps)      !evolve all pdf's
      enddo

      write(lunout,'('' xmin ='',1PE11.2,3X,''subgrids  ='',I4)')
     &  xmin(1), nxmi
      write(lunout,'('' qmin ='',1PE11.2,3X,''qmax      ='',1PE11.2)')
     &  qq(1), qq(2)
      write(lunout,'('' nfin ='',I4,10X,''q2c, b, t ='',3(1PE11.2))')
     &  nfin, q2c, q2b, q2t
      write(lunout,'('' alfa ='',1PE11.2,3X,''r20       ='',1PE11.2)')
     &  as0, r20
      write(lunout,'('' aar  ='',1PE11.2,3X,''bbr       ='',1PE11.2)')
     &  aar, bbr
      write(lunout,'('' q20  ='',1PE11.2)') q0
      write(lunout,'(/)')

C--   ------------------------------------------------------------------
C--   User code
C--   ------------------------------------------------------------------

      lun = nxtlun(0)
      call READTAB(w,nw,lun,fnamu(1),key,0,iset_piju,nwu,ie)
      if(ie.ne.0) then
        call FillWtU(w,nw,iset_piju)
        call DUMPTAB(w,iset_piju,lun,fnamu(1),key)
      endif  

      call READTAB(w,nw,lun,fnamu(2),key,0,iset_pijp,nwu,ie)
      if(ie.ne.0) then
        call FillWtP(w,nw,iset_pijp)
        call DUMPTAB(w,iset_pijp,lun,fnamu(2),key)
      endif
        
      call MAKETAB(w,nw,itypes_alf,0,0,iset_alfa, nwu)
      call MAKETAB(w,nw,itypes_pdf,0,0,iset_pdfu, nwu)
      call MAKETAB(w,nw,itypes_pdf,0,0,iset_pdfp, nwu)

      write(lunout,'('' iset_piju ='',I3)') iset_piju
      write(lunout,'('' iset_pijp ='',I3)') iset_pijp
      write(lunout,'('' iset_alfa ='',I3)') iset_alfa
      write(lunout,'('' iset_pdfu ='',I3)') iset_pdfu
      write(lunout,'('' iset_pdfp ='',I3)') iset_pdfp

      end

C--   ------------------------------------------------------------------
C--   Package routines
C--   ------------------------------------------------------------------

C     ======================================
      subroutine FillWtU( w, nw, iset_piju )
C     ======================================
      implicit double precision (a-h,o-z)
      dimension w(nw)
      dimension itypes(6)
      data itypes / 2, 2, 0, 0, 0, 0 / 
      
      call MAKETAB(w, nw, itypes, 0, 0, iset_piju, nwused)

C--   Put here the code to fill unpolarised weight tables

      return
      end

C     ======================================
      subroutine FillWtP( w, nw, iset_pijp )
C     ======================================
      implicit double precision (a-h,o-z)
      dimension w(nw)
      dimension itypes(6)
      data itypes / 2, 2, 0, 0, 0, 0 / 
      
      call MAKETAB(w, nw, itypes, 0, 0, iset_pijp, nwused)

C--   Put here the code to fill polarised weight tables

      return
      end
      
C--   ------------------------------------------------------------------
C--   Pdf start functions
C--   ------------------------------------------------------------------

C     ======================================
      double precision function func(ipdf,x)
C     ======================================

      implicit double precision (a-h,o-z)

                     func = 0.D0
      if(ipdf.eq. 0) func = xglu(x)
      if(ipdf.eq. 1) func = xdnv(x)
      if(ipdf.eq. 2) func = xupv(x)
      if(ipdf.eq. 3) func = 0.D0
      if(ipdf.eq. 4) func = xdbar(x)
      if(ipdf.eq. 5) func = xubar(x)
      if(ipdf.eq. 6) func = xsbar(x)
      if(ipdf.eq. 7) func = 0.D0
      if(ipdf.eq. 8) func = 0.D0
      if(ipdf.eq. 9) func = 0.D0
      if(ipdf.eq.10) func = 0.D0
      if(ipdf.eq.11) func = 0.D0
      if(ipdf.eq.12) func = 0.D0

      return
      end
 
C     =================================
      double precision function xupv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data au /5.107200D0/

      xupv = au * x**0.8D0 * (1.D0-x)**3.D0

      return
      end

C     =================================
      double precision function xdnv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ad /3.064320D0/

      xdnv = ad * x**0.8D0 * (1.D0-x)**4.D0

      return
      end

C     =================================
      double precision function xglu(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ag    /1.7D0 /

      common /msum/ glsum, uvsum, dvsum

      xglu = ag * x**(-0.1D0) * (1.D0-x)**5.D0

      return
      end

C     ==================================
      double precision function xdbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      data adbar /0.1939875D0/

      xdbar = adbar * x**(-0.1D0) * (1.D0-x)**6.D0

      return
      end

C     ==================================
      double precision function xubar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xubar = xdbar(x) * (1.D0-x)

      return
      end

C     ==================================
      double precision function xsbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xsbar = 0.2D0 * (xdbar(x)+xubar(x))

      return
      end







