C     ----------------------------------------------------------------
      program example
C     ----------------------------------------------------------------
C--   12-03-07  Basic QCDNUM example job
C     ----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      data as0/0.364/, r20/2.D0/, iord/3/, nfin/0/ !alphas, NNLO, VFNS
      external func                                !input parton dists
      dimension def(-6:6,12)                       !flavor composition
      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
     + 52*0.    /
      data xmin/1.D-4/, nxin/100/, iosp/3/             !x grid, splord
      dimension qq(2),wt(2)                                  !mu2 grid
      data qq/2.D0,1.D4/, wt/1.D0,1.D0/, nqin/60/            !mu2 grid
      data q2c/3.D0/, q2b/25.D0/, q0/2.0/             !thresh and mu20
      data x/1.D-3/, q/1.D3/, qmz2/8315.25D0/           !output scales
      dimension pdf(-6:6)                                      !pdfout

C--   lun = 6 stdout, -6 stdout w/out banner page
      lun    = 6
      lunout = abs(lun)      
      call qcinit(lun,' ')                                 !initialize
      call gxmake(xmin,1,1,nxin,nx,iosp)                       !x-grid
      call gqmake(qq,wt,2,nqin,nq)                           !mu2-grid
      call fillwt(1,id1,id2,nw)                     !calculate weights
      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
      iq0  = iqfrmq(q0)                                !starting scale
      call evolfg(1,func,def,iq0,eps)                !evolve all pdf's
      call allfxq(1,x,q,pdf,0,1)                !interpolate all pdf's

      csea = 2.D0*pdf(-4)                           !charm sea at x,Q2
      asmz = asfunc(qmz2,nfout,ierr)                      !alphas(mz2)

      write(lunout,'('' x, q, CharmSea ='',3E13.5)'), x,q,csea
      write(lunout,'('' as(mz2)        ='', E13.5)'), asmz

      end
      
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