
C     ------------------------------------------------------------------
      program Toolbox07
C     ------------------------------------------------------------------
C--   QCDNUM toolbox example #7
C--
C--   Full FFNS and VFNS evolution.
C--   Added an interpolation routine for gluon, q and qbar
C--
C--   QCDNUM  routines are in capitals
C--   Package routines are in mixed case (except the splitting functions)
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/0/, 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, usub
      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 = 231327 )
      dimension w(nw)
      character*12 key
      data key /'MyAddon v2.0'/
      dimension xxi(5), qqi(5), ffu(5), ffq(5)
      data xxi / 1.D-1, 2.D-1, 3.D-1, 4.D-1, 5.D-1 /
      data qqi / 1.D+2, 2.D+2, 3.D+2, 4.D+2, 5.D+2 /
      dimension qqb(-6:6)
      data qqb /13*0.D0/

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
      iqt  = IQFRMQ(q2t)
      call SETCBT(nfin,iqc,iqb,iqt)             !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--   ------------------------------------------------------------------

      iset    = 1              !1=unpol  2=pol

      call MyWeight(w, nw, iset, key)

      idebug  = 0
      call MyEvolve(w, iset, usub, def, iq0, nfmax, idebug)

      iq = nq                 !comparison grid point
      do id = 0,12
        call CompareF(w, iset, id, iq, 0, dif)
      enddo

C--   Interpolate and compare with QCDNUM routine FFLIST
      ichk = 1
      call GETINT('lunq',lun)
      write(lun,'(/)')

      do idf = -6,6
C--     My interpolation
        call Interpolate( w, iset, idf, xxi, qqi, ffu, 5, ichk )
C--     QCDNUM interpolation
        qqb(idf) = 1.D0   !select idf
        call FFLIST( iset, qqb, 0, xxi, qqi, ffq, 5, ichk )
        qqb(idf) = 0.D0   !deselect idf
C--     Compare
        dif = 0.D0
        do i = 1,5
          dif = max(dif,abs(ffu(i)-ffq(i)))
        enddo
        write(lun,'('' pdf'',I3,'' dif ='',(1PE10.3))') idf, dif
      enddo

      end

C--   ------------------------------------------------------------------
C--   User routines
C--   ------------------------------------------------------------------

C     ==============================================
      subroutine CompareF(ww,iset,idf,iq,nprint,dif)
C     ==============================================

C--   Compare pdfs
C--
C--   Remark: uses EVPDFIJ to access both local and QCDNUM pdfs
C--
C--   ww        (in)      workspace
C--   iset      (in)      1=unpolarised 2=polarised
C--   idf       (in)      pdf identifier [0-12]
C--   iq        (in)      q2 grid-point
C--   nprint    (in)      > 0 print first n comparisons
C--                       = 0 print only dif
C--                       < 0 run quiet
C--   dif       (out)     max absolute difference encountered

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

      dimension ww(*)

      call GRPARS(nx,xmi,xma,nq,qmi,qma,iosp)
      call GETINT('lunq',lun)
      call SetGetI('Get','isetpdf',isetpdf,iset)
      idw = 1000*isetpdf + 501 + idf     !address in local workspace
      idq = IPDFTAB(iset,idf)            !address in QCDNUM memory

      dif = 0.D0
      do ix = 1,nx
        uuu = EVPDFIJ(ww,idw,ix,iq,1)    !pdf from local workspace ww
        fff = EVPDFIJ(ww,idq,ix,iq,1)    !pdf from QCDNUM memory
        ddd = uuu-fff
        dif = max(dif,abs(ddd))
        if(ix.le.nprint)
     +  write(lun,'(I4,2E13.5,F10.6,1PE13.3)') ix,uuu,fff,ddd,ddd/fff
      enddo
      if(nprint.ge.0) then
        write(lun,'('' pdf'',I3,'' dif ='',(1PE10.3))') idf, dif
      endif

      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: User interface routines
C--   ------------------------------------------------------------------

C     ======================================
      subroutine MyWeight(w, nw, iset, key)
C     ======================================

C--   Partition workspace and fill weight tables
C--
C--   w     (in) :    workspace
C--   nw    (in) :    number of words in w
C--   key   (in) :    weight-file tag

      implicit double precision (a-h,o-z)
      dimension w(*)
      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*(*) key
C--   Weight files
      character*24 fnamu(2)
C--                123456789012345678901234
      data fnamu /'../weights/QcdUsrUnp.wgt',
     +            '../weights/QcdUsrPol.wgt' /

C--   Create alfas table if it does not yet exist
      call SetGetI('Get','ialftab',ialftab,0)
      if(ialftab.eq.0) then
        call MAKETAB(w,nw,itypes_alf,0,0,isetalf,nwu)
        call SetGetI('Set','isetalf',isetalf,0)  !store table set id
        call SetGetI('Set','ialftab',   1   ,0)  !alfas table exists
      endif
C--   Weight tables
      lun = nxtlun(0)
      call READTAB(w,nw,lun,fnamu(iset),key,0,isetpij,nwu,ie)
      if(ie.ne.0) then
        call FillWtU(w,nw,isetpij)
        call DUMPTAB(w,isetpij,lun,fnamu(iset),key)
      endif
      call SetGetI('Set','isetpij',isetpij,iset) !store table set id
C--   Pdf tables
      call MAKETAB(w,nw,itypes_pdf,0,0,isetpdf, nwu)
      call SetGetI('Set','isetpdf',isetpdf,iset) !store table set id

      return
      end


C     ========================================================
      subroutine MyEvolve(w, iset, subr, def, iq0, nfmax, idb)
C     ========================================================

C--   Evolve all pdfs in the FFNS or VFNS
C--   Starting values are taken from user subroutine SUBR
C--
C--   w       (in) : workspace
C--   iset    (in) : set identifier 1=unpol 2=pol
C--   subr    (in) : user supplied subroutine returning starting pdfs
C--   def     (in) : def(i,j) contribution of flavour i to input pdf j
C--   iq0     (in) : starting q2 grid point
C--   nfmax  (out) : max number of flavours encountered
C--   idb     (in) : 1 = some debug printout

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

      dimension w(*)
      parameter (nxmax = 200) !change this when you use a larger xgrid
      dimension start(2,nxmax), stval(0:12,nxmax)
      dimension ith(4:6), def(-6:6,12)
      external  subr

      itab(i) = 1000*isetf + 501 + i             !statement function

C--   Get all start values
      call UsrStart(subr,def,iq0,stval)

      call GRPARS(nx,xmi,xma,nq,qmi,qma,iosp)
      if(nx.gt.nxmax) stop 'MyEvolve: too many x-grid points'
      call GETCBT(nf, qc, qb, qt)
      call SetGetI('Get','isetalf',iseta,0)
      call SetGetI('Get','isetpij',isetp,iset)
      call SetGetI('Get','isetpdf',isetf,iset)

C--   pdf   gl    si  e2+ e3+ e4+ e5+ e6+    va  e2- e3- e4- e5- e6-
C--    id    0     1   2   3   4   5   6      7   8   9  10  11  12
C--    nf          3   3   3   4   5   6      3   3   3   4   5   6

C--   Singlet/gluon
      call CpyStart(1, stval, start)       !singlet
      call CpyStart(0, stval, start)       !gluon
      call EvolSGNS(w, iseta, isetp, isetf, 1, iq0, start, 0, idb)

      if( nf.eq.0 ) then
        idmax1 = 3                             !VFNS light quarks only
        idmax2 = 9
      else
        idmax1 = nf                            !FFNS all quarks
        idmax2 = nf+6
      endif

      nfmax = nf

C--   Nonsinglets (only light nonsinglets in the VFNS)
      do idf = 2,idmax1
        call CpyStart(idf, stval, start)   !nonsinglet e^+
        call EvolSGNS(w, iseta, isetp, isetf, idf, iq0, start, 0, idb)
      enddo
      do idf = 7,idmax2
        call CpyStart(idf, stval, start)   !nonsinglet e^-
        call EvolSGNS(w, iseta, isetp, isetf, idf, iq0, start, 0, idb)
      enddo

C--   Initialise the heavy quark pdfs
      do i = idmax1+1,6
        call copywgt(w, itab(1), itab(i), 0)   !copy singlet to ei^+
      enddo
      do i = idmax2+1,12
        call copywgt(w, itab(7), itab(i), 0)   !copy valence to ei^-
      enddo

C--   Heavy quark nonsinglets in the VFNS; evolve upward from threshold
      ith(4) = IQFRMQ(qc)            !zero if q2c above Q2 range
      ith(5) = IQFRMQ(qb)            !zero if q2b above Q2 range
      ith(6) = IQFRMQ(qt)            !zero id q2t above Q2 range
      if( nf.eq.0 )  then
        nfmax = 3
        do id = 4,6
          if(ith(id).ne.0) then
            i0    = ith(id)
            nfmax = id
            call SetStart(iset, 1, i0, start)
            call EvolSGNS(w, iseta, isetp, isetf, id, i0, start, 1, idb)
            jd = id+6
            call SetStart(iset, 7, i0, start)
            call EvolSGNS(w, iseta, isetp, isetf, jd, i0, start, 1, idb)
          endif
        enddo
      endif

      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: Interpolation
C--   ------------------------------------------------------------------

C     =========================================================
      subroutine Interpolate( w, iset, idf, x, q, f, n, ichk )
C     =========================================================

C--   Interpolate idf = u,d,s,..  to a list of interpolation points
C--
C--   w           (in) : local workspace
C--   iset        (in) : 1=unpolarised, 2=polarised
C--   idf         (in) : [-6,6] for (anti)quarks (0=gluon)
C--   x(n), q(n)  (in) : list of interpolation points
C--   f(n)       (out) : list of interpolated pdfs
C--   n           (in) : number of interpolation points
C--   ichk        (in) : yes/no check boundaries

      implicit double precision (a-h,o-z)
      dimension w(*), x(*), q(*), f(*)
      dimension coefg(3:6), coefd(3:6,12)
      data coefg /4*1.D0/


      id(i) = 1000*isetf + 501 + i                  !statement function

      call SetGetI('Get','isetpdf',isetf,iset)      !pdf set id
      call FASTINI(x, q, n, ichk)                   !interpolation mesh
      if( idf.eq.0 ) then 
        call FASTINP(w, id(0), coefg, -1, 0)        !put gluon in buf1
      else
        call getcoefs(idf, coefd)                   !get coefficients
        call FASTINP(w, id(1), coefd(3,1), -1, 0)   !put c1*id1 in buf1
        do i = 2,12
          call FASTINP(w, id(i), coefd(3,i), -1, 1) !add ci*idi to buf1
        enddo
      endif
      call FASTFXQ( 1, f, n )                       !interpolate buf1

      return
      end

C     =================================
      subroutine getcoefs( idf, coefd )
C     =================================

C--   Get coefficients of the expansion of u,d,s,.. to the basis pdfs
C--
C--   |idf> = Sumj Cj(nf) |ej>    j = 1,...,12
C--
C--   idf            (in) : (anti)quark identifier [-6,6]
C--   coefd(3:6,12) (out) : coefficients Cj(nf)

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

      dimension qvec(-6:6), temp(12), coefd(3:6,12)
      data qvec /13*0.D0/

      qvec(idf) = 1.D0                !select idf (and nothing else)
      do nf = 3,6
        call EFROMQQ(qvec, temp, nf)
        do i = 1,12
          coefd(nf,i) = temp(i)
        enddo
      enddo
      qvec(idf) = 0.D0                !deselect idf

      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: Evolution
C--   ------------------------------------------------------------------

C     ========================================
      subroutine CpyStart( idf, stval, start )
C     ========================================

C--   Copy stval(idf,ix) to start(1,ix) for quarks and to start(2,ix)
C--   for gluons. The array stval is taken from UsrStart and the
C--   output array start is fed into to EvolSGNS.

      implicit double precision (a-h,o-z)
      dimension stval(0:12,*), start(2,*)

      i = 1                     !copy quqarks
      if(idf.eq.0) i = 2        !copy gluon
      call GRPARS( nx, xmi, xma, nq, qmi, qma, iord )
      do ix = 1,nx
          start(i,ix) = stval(idf,ix)
      enddo

      return
      end

C     =============================================
      subroutine UsrStart( subr, abmat, iq0, stval )
C     =============================================

C--   subr(ix,pdfusr) (in) : returns startfunction fj(ix) j=1,...2nf
C--   abmat(-6:6,12)  (in) : abmat(i,j) contribution of flavour i to fj
C--   iq0             (in) : start point
C--   stval(0:12,nx) (out) : stval(i,ix) value of basis pdf e_i(ix,iq0)

      implicit double precision (a-h,o-z)
      external  subr
      dimension abmat(-6:6,12), dinv(12,12), pdfusr(0:12)
      dimension stval(0:12,*)
      save nflast
      data nflast /0/

      nf = NFLAVS(iq0,ithrs)
      if(nf.ne.nflast) then
        call GetDinv( abmat, dinv, nf )
        nflast = nf
      endif
      call GRPARS( nx, xmi, xma, nq, qmi, qma, iord )
      do ix = 1,nx
        call subr( ix, pdfusr )
        stval(0,ix) = pdfusr(0)
        do i = 1,nf
          stval(i  ,ix) = 0.D0
          stval(i+6,ix) = 0.D0
          do j = 1,2*nf
            stval(i  ,ix) = stval(i  ,ix) + dinv(i  ,j)*pdfusr(j)
            stval(i+6,ix) = stval(i+6,ix) + dinv(i+6,j)*pdfusr(j)
          enddo
        enddo
      enddo
      
      return
      end

C     =====================================
      subroutine GetDinv( abmat, dinv, nf )
C     =====================================

C--   The basis pdfs |e> are given by a linear combination of the 2nf
C--   user input pdfs |f> :
C--
C--               |ei> = Sum_j dinv(i,j) |fj>  j = 1,...2nf
C--
C--   This routine calculates the coefficients dinv(i,j)
C--
C--   abmat(-6:6,12) (in) : abmat(i,j) contribution of flavour i to |fj>
C--   dinv(12,12)   (out) : dinv(i,j)  contribution of |fj> to |ei>
C--   nf             (in) : number of active flavours [3-6]

      implicit double precision (a-h,o-z)
      dimension abmat(-6:6,12), d(12), dmat(12,12), dinv(12,12)

C--   Build dmat
      do i = 1,2*nf
        call EFROMQQ( abmat(-6,i), d, nf )
        do j = 1,12
            dmat(i,j) = d(j)
        enddo
      enddo  
C--   Invert dmat
      call InvertD( dmat, dinv, nf, ierr )
      if( ierr.ne.0 ) stop 'GetDinv : singular matrix'

      return
      end

C     =========================================
      subroutine InvertD( dmat, dinv, nf ,ierr)
C     =========================================

C--   Copy dmat to a 2nf*2nf work matrix, invert, and copy back to dinv
C--
C--   Indexing      ep1 ep2 ep3 ep4 ep5 ep6 em1 em2 em3 em4 em5 em6
C--   dmat, dinv      1   2   3   4   5   6   7   8   9  10  11  12
C--   work  nf 3      1   2   3   7   8   9
C--   work  nf=4      1   2   3   4   7   8   9  10
C--   work  nf=5      1   2   3   4   5   7   8   9  10  11
C--   work  nf=6      1   2   3   4   5   6   7   8   9  10  11  12

      implicit double precision (a-h,o-z)
      dimension dmat(12,12), dinv(12,12), work(12,12), iw(12)

C--   Initialise dinv
      do i = 1,12
        do j = 1,12
          dinv(i,j) = 0.D0
        enddo
      enddo
C--   Copy dmat to work array
      do i = 1,2*nf
        do j = 1,nf
            work(i,j   ) = dmat(i,j  )
            work(i,j+nf) = dmat(i,j+6)
        enddo
      enddo  
C--   Invert work array        
      call SMB_DMINV( 2*nf, work, 12, iw, ierr )
      if( ierr.ne.0 ) return
C--   Copy work array to dinv 
      do i = 1,nf
        do j = 1,2*nf
          dinv(i  ,j) = work(i   ,j)
          dinv(i+6,j) = work(i+nf,j)
        enddo
      enddo
                   
      return
      end

C     ============================================
      subroutine SetStart( iset, idf, iq0, start )
C     ============================================

C--   Pickup qcdnum pdf at iq0 and put into start array
C--   start(1,ix)  : singlet or nonsinglet quark pdf
C--   start(2,ix)  : gluon

      implicit double precision (a-h,o-z)
      dimension start(2,*)

      i = 1                    !quarks
      if( idf.eq.0 ) i = 2     !gluon
      call GRPARS(nx, xmi, xma, nq, qmi, qma, iord)
      do ix = 1,nx
        start(i,ix) = BVALIJ(iset, idf, ix, iq0, 1)
      enddo

      return
      end

C     =================================================================
      subroutine EvolSGNS(w,iseta,isetp,isetf,idfin,iq0,start,iup,idbg)
C     =================================================================

C--   Singlet/gluon or non-singlet evolution

C--   w         (in) : local workspace
C--   iseta,p,f (in) : table set identifiers for alfas, pij and pdf
C--   idfin     (in) : 1 = singlet/gluon 2-12 = nonsinglet
C--   iq0       (in) : start point
C--   start     (in) : start array
C--   iup       (in) : = 1 upward evolution only, otherwise up and down
C--   idbg      (in) : = 1 some debug printout

      implicit double precision (a-h,o-z)
      logical   AtabInvalid
      external  AlfasFun
      dimension w(*), start(2,*)
      parameter (nxmax = 200) !change this when you use a larger xgrid
      dimension sbuf(2,nxmax)
      dimension idw(2,2,1), ida(2,2,1), idf(2), iqlim(2)

      call GRPARS( nx, xmi, xma, nq, qmi, qma, iosp )
      if( nx.gt.nxmax) stop 'EvolSGNS: too many x-grid points'
      call GETORD(io)
      if(io.ne.1) stop 'EvolSGNS: can only handle LO'
      call GETINT('lunq',lun)  !for debug printout
      
      idalf      = 1000*iseta + 601
      idw(1,1,1) = 1000*isetp + 101           !PQQ
      idw(1,2,1) = 1000*isetp + 201           !PQG
      idw(2,1,1) = 1000*isetp + 102           !PGQ
      idw(2,2,1) = 1000*isetp + 202           !PGG
      ida(1,1,1) = idalf
      ida(1,2,1) = idalf
      ida(2,1,1) = idalf
      ida(2,2,1) = idalf
      idf(1)     = 1000*isetf + 501 + idfin   !quark table
      idf(2)     = 1000*isetf + 501           !gluon table 

      if(AtabInvalid()) call EVFILLA( w, idalf, AlfasFun )      

      n = 1                          !non-singlet evolution      
      if( idfin.eq.1 ) n = 2         !singlet/gluon evolution

      do i = 1,n
        do j = 1,nx
           sbuf(i,j) = start(i,j)    !copy start array
        enddo   
      enddo
      iqlim(2) = iq0
      nf       = 1
      do while( nf.gt.0 )            !upward evolution
        iqlim(1) = iqlim(2)
        iqlim(2) = 99999
        call EVDGLAP( w, idw, ida, idf, sbuf, 2, n, iqlim, nf, eps )
        if(idbg.eq.1) then
          write(lun,'(''Up evolve'',I3,'' from'',I3,'' to'',I3)')
     +    idfin,iqlim(1),iqlim(2)
        endif
      enddo  

      if( iup.ne.1 ) then
        do i = 1,n
          do j = 1,nx
             sbuf(i,j) = start(i,j)    !copy start array
          enddo
        enddo
        iqlim(2) = iq0
        nf       = 1
        do while( nf.gt.0 )            !downward evolution
          iqlim(1) = iqlim(2)
          iqlim(2) = -99999
          call EVDGLAP( w, idw, ida, idf, sbuf, 2, n, iqlim, nf, eps )
          if(idbg.eq.1) then
             write(lun,'(''Dn evolve'',I3,'' from'',I3,'' to'',I3)')
     +       idfin,iqlim(1),iqlim(2)
          endif
        enddo
      endif

      if(idbg.eq.1) write(lun,'(1X)')
          
      return
      end

C     =====================================================
      double precision function AlfasFun( iq, nf, ithresh )
C     =====================================================
      implicit double precision (a-h,o-z)

      idum = nf                    !avoid compiler warning
      if( ithresh .eq. -1 ) then
        AlfasFun = ALTABN( 0, -iq, 1, ierr )  !alfas/2pi
      else
        AlfasFun = ALTABN( 0,  iq, 1, ierr )  !alfas/2pi
      endif    

      return
      end

C     ==============================
      logical function AtabInvalid()
C     ==============================
      implicit double precision (a-h,o-z)
      save ioL, asL, r2L, nfL, qcL, qbL, qtL, arL, brL
      data ioL/0/, asL/-1.D0/, r2L/-1.D0/, nfL/9/, qcL/-1.D0/
      data qbL/-1.D0/, qtL/-1.D0/, arL/-1.D0/, brL/-1.D0/

      call GETORD( io )
      call GETALF( as, r2 )
      call GETCBT( nf, qc, qb, qt )
      call GETABR( ar, br)

      if(io.eq.ioL .and. as.eq.asL .and. r2.eq.r2L .and.
     &   nf.eq.nfL .and. qc.eq.qcL .and. qb.eq.qbL .and. qt.eq.qtL .and.
     &   ar.eq.arL .and. br.eq.brL) then
        AtabInvalid = .false.
      else
        AtabInvalid = .true.
        ioL = io
        asL = as
        r2L = r2
        nfL = nf
        qcL = qc
        qbL = qb
        qtL = qt
        arL = ar
        brL = br
      endif

      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: Weights
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 /

      external ACHI,PQQR,PQQS,PQQD,PQGA,PGQA,PGGR,PGGS,PGGA,PGGD
      
      call MAKETAB( w, nw, itypes, 0, 0, iset_piju, nwused ) 

C--   Fill unpolarised weight tables

      call getint('lunq',lun)
      write(lun,'(/'' FILLWTU: start unpolarised weight calculation'')')

      idPQQ = 1000*iset_piju + 101
      idPQG = 1000*iset_piju + 201
      idPGQ = 1000*iset_piju + 102
      idPGG = 1000*iset_piju + 202
      
      call MAKEWRS( w, idPQQ, PQQR, PQQS, ACHI, 0 )
      call MAKEWTD( w, idPQQ, PQQD, ACHI )
      
      call MAKEWTA( w, idPQG, PQGA, ACHI )
      
      call MAKEWTA( w, idPGQ, PGQA, ACHI )
      
      call MAKEWRS( w, idPGG, PGGR, PGGS, ACHI, 0 )
      call MAKEWTA( w, idPGG, PGGA, ACHI )
      call MAKEWTD( w, idPGG, PGGD, ACHI )

      write(lun,'( '' FILLWTU: end   unpolarised weight calculation'')')
      
      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
      call getint('lunq',lun)
      write(lun,'(/'' FILLWTP: empty routine'')')

      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: Unpolarised splitting functions
C--   ------------------------------------------------------------------

C     ====================================
      double precision function ACHI(qmu2)
C     ====================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      ACHI = 1.D0
      return
      end

C     =========================================
      double precision function PQQR(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PQQR = (4.D0/3.D0)*(1.D0+x*x)
      return
      end

C     =========================================
      double precision function PQQS(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PQQS = 1.D0/(1.D0-x)
      return
      end

C     =========================================
      double precision function PQQD(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = x             !avoid compiler warning
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PQQD = 2.D0
      return
      end

C     =========================================
      double precision function PQGA(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      PQGA = nf * ( x**2 + (1.D0-x)**2 )
      return
      end

C     =========================================
      double precision function PGQA(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PGQA = (4.D0/3.D0) * ( 1.D0 + (1.D0-x)**2 ) / x
      return
      end

C     =========================================
      double precision function PGGR(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PGGR = 6.D0 * x
      return
      end

C     =========================================
      double precision function PGGS(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PGGS = 1.D0/(1.D0-x)
      return
      end

C     =========================================
      double precision function PGGA(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = qmu2          !avoid compiler warning
      idum = nf            !avoid compiler warning
      PGGA = 6.D0 * ( (1.D0-x)/x + x*(1.D0-x) )
      return
      end

C     =========================================
      double precision function PGGD(x,qmu2,nf)
C     =========================================
      implicit double precision (a-h,o-z)
      ddum = x             !avoid compiler warning
      ddum = qmu2          !avoid compiler warning
      PGGD = 11.D0/2.D0 - nf/3.D0
      return
      end

C--   ------------------------------------------------------------------
C--   Package routines: Store and retreive
C--   ------------------------------------------------------------------

C     ======================================
      subroutine SetGetI(action,ivar,ival,i)
C     ======================================

C--   Store or retreive value of variable ivar or array element ivar(i)
C--   Variable names and array dimensions are hardwired in this code
C--
C--   action   (in) : 'Set' or 'Get' first character in upper case
C--   ivar     (in) :  Character string with variable name in lower case
C--   ival  (inout) :  Integer value passed
C--   i        (in) :  i = irrelevant when ivar is not an array
C--                    i = in range     --> transfer ivar(i) <--> ival
C--                    i = out of range --> do nothing

      character*(*) action, ivar

      save isetalf, isetpij, isetpdf
      dimension isetpij(2), isetpdf(2)
      data isetalf/0/, isetpij/2*0/ isetpdf/2*0/

      save ialftab, ipijtab, ipdftab
      dimension ipijtab(2), ipdftab(2)
      data ialftab/0/, ipijtab/2*0/ ipdftab/2*0/

C---  Set
      if(action(1:1).eq.'S') then
        if    (ivar(1:7).eq.'isetalf') then
          isetalf = ival
        elseif(ivar(1:7).eq.'isetpij') then
          if(i.ge.1 .and. i.le.2) isetpij(i) = ival
        elseif(ivar(1:7).eq.'isetpdf') then
          if(i.ge.1 .and. i.le.2) isetpdf(i) = ival

        elseif(ivar(1:7).eq.'ialftab') then
          ialftab = ival
        elseif(ivar(1:7).eq.'ipijtab') then
          if(i.ge.1 .and. i.le.2) ipijtab(i) = ival
        elseif(ivar(1:7).eq.'ipdftab') then
          if(i.ge.1 .and. i.le.2) ipdftab(i) = ival

        else
          stop 'SetGetI : unknown variable'
        endif

C--   Get
      elseif(action(1:1).eq.'G') then
        if    (ivar(1:7).eq.'isetalf') then
          ival = isetalf
        elseif(ivar(1:7).eq.'isetpij') then
          if(i.ge.1 .and. i.le.2) ival = isetpij(i)
        elseif(ivar(1:7).eq.'isetpdf') then
          if(i.ge.1 .and. i.le.2) ival = isetpdf(i)

        elseif(ivar(1:7).eq.'ialftab') then
          ival = ialftab
        elseif(ivar(1:7).eq.'ipijtab') then
          if(i.ge.1 .and. i.le.2) ival = ipijtab(i)
        elseif(ivar(1:7).eq.'ipdftab') then
          if(i.ge.1 .and. i.le.2) ival = ipdftab(i)

        else
          stop 'SetGetI : unknown variable'
        endif

      else
        stop 'SetGetI : action not S(et) or G(et)'
      endif

      return
      end

C--   ------------------------------------------------------------------
C--   User routines: Pdf start functions
C--   ------------------------------------------------------------------

C     =======================
      subroutine usub(ix,pdf)
C     =======================

      implicit double precision (a-h,o-z)
      dimension pdf(0:12)

      x       = xfrmix(ix)
      pdf( 0) = xglu(x)
      pdf( 1) = xdnv(x)
      pdf( 2) = xupv(x)
      pdf( 3) = 0.D0
      pdf( 4) = xdbar(x)
      pdf( 5) = xubar(x)
      pdf( 6) = xsbar(x)
      pdf( 7) = 0.D0
      pdf( 8) = 0.D0
      pdf( 9) = 0.D0
      pdf(10) = 0.D0
      pdf(11) = 0.D0
      pdf(12) = 0.D0

      return
      end

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

