[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM.f

      subroutine QCDNUMevolve(x,Q,pdf)
      implicit none
      character*64 gridname
      character*16 s1,dummy
      integer Eorder,index
      real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
      real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
      real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
      real*8 QALFAS,QPDFXQ,XFROMIX
      real*8 f(-6:6),pdf(-6:6)
      integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
      integer nx,nq
      integer nset,iset,isetlast
      data isetlast/-1/
      real*8 xmin,xmax,qmin,qmax,S
      save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
      save nx,xmin,xmax,nq,qmin,qmax,gridname,isetlast
*
      call getnset(iset)
c fix 14.7.2006 to speed up sue of QCDNUM
      if(iset.ne.isetlast) then
        call get_pdfqcd(iset)
        isetlast = iset
      endif
c      
      Q2=Q*Q
      nf=5
      if (Q2.lt.tb2) nf=4
      if (Q2.lt.tc2) nf=3
      pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
      singlet=QPDFXQ('SINGLET',x,Q2,iflag)
      dm= QPDFXQ('DM',x,Q2,IFLAG)
      um= QPDFXQ('UM',x,Q2,IFLAG)
      dp= QPDFXQ('DP',x,Q2,IFLAG)
      up= QPDFXQ('UP',x,Q2,IFLAG)
      sp= QPDFXQ('SP',x,Q2,IFLAG)
      ub=0.5d0*(up-um+singlet/dble(nf))
      db=0.5d0*(dp-dm+singlet/dble(nf))
      sb=0.5d0*(sp+singlet/dble(nf))
      cb=0.d0
      if (nf.ge.4) then
         cp= QPDFXQ('CP',X,Q2,IFLAG)
         cb=0.5d0*(cp+singlet/dble(nf))
      end if
      bb=0.d0
      if (nf.ge.5) then
         bp= QPDFXQ('BP',X,Q2,IFLAG)            
         bb=0.5d0*(bp+singlet/dble(nf))
      end if

      pdf(1)=dm+db
      pdf(2)=um+ub
      pdf(3)=sb
      pdf(4)=cb
      pdf(5)=bb
      pdf(6)=0d0
      pdf(-1)=db
      pdf(-2)=ub
      pdf(-3)=sb
      pdf(-4)=cb
      pdf(-5)=bb
      pdf(-6)=0d0
      return
*
      entry QCDNUMalfa(alfas,Q)
      Q2=Q*Q
      nf=6
      if (Q2.lt.mt2) nf=5
      if (Q2.lt.mb2) nf=4
      if (Q2.lt.mc2) nf=3
      alfas=QALFAS(Q2,Qlam,nf,iflag)
      return
*
      entry QCDNUMread(nset)
c      read(1,*) s1
c      print *,s1
c         gridname='large.grid'
c         nx=400
c         xmin=1d-6
c         xmax=1d0
c         nq=112
c         qmin=1d0
c         qmax=1d10
c         read(1,*) dummy
*      else
         read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
c*      endif
c      iset = nset
      return
*
      entry QCDNUMinit(nset,Eorder,Q2fit)
      call qninit
      call QNISET('ORDER',Eorder+1)
      call grxdef(nx,xmin)
      call grqdef(nq,qmin,qmax)
      call grqinp(Q2fit,1)
      call getQMassM(nset,4,mc)
      mc2=mc*mc
      CALL QNRSET('CMASS',mc)
      CALL QNRSET('MCALF',mc)
      call getThresholdM(nset,4,tc)
      tc2=tc*tc
      call GRQINP(tc2,1)
      CALL GRQINP(tc2-1.0d-4,1)
      call getQMassM(nset,5,mb)
      mb2=mb*mb
      CALL QNRSET('BMASS',mb)
      CALL QNRSET('MBALF',mb)
      call getThresholdM(nset,5,tb)
      tb2=tb*tb
      call GRQINP(tb2,1)
      CALL GRQINP(tb2-1.0d-4,1)
      call getQMassM(nset,6,mt)
      mt2=mt*mt
      CALL QNRSET('TMASS',mt)
      CALL QNRSET('MTALF',mt)
      call qthres(tc2,tb2)
      iq0=IQFROMQ(Q2fit)
      iqc=IQFROMQ(tc2)
      iqb=IQFROMQ(tb2)
      nf0= NFLGET(iq0)
c      print *,nf0
      CALL QNBOOK(2,'dm')
      CALL QNBOOK(3,'um')
      CALL QNBOOK(4,'dp')
      CALL QNBOOK(5,'up')
      CALL QNBOOK(6,'sp')
      CALL QNBOOK(7,'cp')
      CALL QNBOOK(8,'bp')
      CALL QNLSET('W1ANA',.TRUE.)
      CALL QNLSET('W2NUM',.TRUE.)
      CALL QNLSET('W2STF',.FALSE.)
      if (index(gridname,'none').eq.1) then
         call qnfilw(0,0)
      else
         qnerr=-1
         open(unit=2,status='old',file=gridname,
     .        form='unformatted',err=1)
         call QNREAD(2,1,qnerr)
 1       close(2)
         if (qnerr.ne.0) then
            write(*,*) 'Grid file problem: ',gridname
            if (qnerr.lt.0) then 
               write(*,*) 'Grid file does not exist'
               write(*,*) 'Calculating and creating grid file'
               call qnfilw(0,0)
               open(unit=2,status='unknown',file=gridname,
     .              form='unformatted')
               call QNDUMP(2)
               close(2)
            else
               write(*,*) 'Existing grid file is inconsistent'
               if (qnerr.eq.1)
     .              write(*,*) 'Defined grid different'
               if (qnerr.eq.2)
     .              write(*,*) 'Heavy quark weight table different'
               if (qnerr.eq.3)
     .              write(*,*) 'Charm mass different'
               if (qnerr.eq.4)
     .              write(*,*) 'Bottom mass different'
               stop
            endif
         endif
      endif
      return
*
      entry QCDNUMpdf(nset)
c      print *,'entering QCDNUMpdf',nset
      call GetAlfas(nset,alfas0,scale0)
c      print *,alfas0,scale0
      Q2=scale0*scale0
      CALL QNRSET('ALFAS',alfas0)      
      CALL QNRSET('ALFQ0',Q2)
      DO ix = 1,nx
         xx = XFROMIX(ix)
c	 print *,'calling parmPDF',ix
         call parmPDF(nset,xx,f)
c	 if(ix.lt.6) print *,nset,xx,f
         singlet=0.d0
         do i=1,nf0
            singlet=singlet+f(i)+f(-i)
c	    print  *,i,singlet
         end do
         gl=f(0)
         dm=f(1)-f(-1)
         um=f(2)-f(-2)
         dp=f(1)+f(-1)-singlet/dble(nf0)
         up=f(2)+f(-2)-singlet/dble(nf0)
         sp=f(3)+f(-3)-singlet/dble(nf0)
c      print *,'calling qpnset sg...',ix,iq0,singlet
         CALL QNPSET('SINGLET',ix,iq0,singlet)
c      print *,'calling qpnset gl...',ix,iq0,gl
         CALL QNPSET('GLUON',ix,iq0,gl)
c      print *,'calling qpnset dm...',ix,iq0,dm
         CALL QNPSET('DM',ix,iq0,DM)
c      print *,'calling qpnset um...',ix,iq0,um
         CALL QNPSET('UM',ix,iq0,UM)
c      print *,'calling qpnset dp...',ix,iq0,dp
         CALL QNPSET('DP',ix,iq0,DP)
c      print *,'calling qpnset up...',ix,iq0,up
         CALL QNPSET('UP',ix,iq0,UP)
c      print *,'calling qpnset sp...',ix,iq0,sp
         CALL QNPSET('SP',ix,iq0,SP)
      ENDDO
c      print *,'calling evols...'
      CALL EVOLSG(iq0,1,nq)
c      print *,'calling evols dm...'
      CALL EVOLNM('DM',iq0,1,nq)
c      print *,'calling evols um...'
      CALL EVOLNM('UM',iq0,1,nq)
c      print *,'calling evols dp...',iq0,iqc,iqb,nq
      CALL EVLSEA('dp',iq0,iqc,iqb,nq)
c      print *,'calling evols up...'
      CALL EVLSEA('up',iq0,iqc,iqb,nq)
c      print *,'calling evols sp...'
      CALL EVLSEA('sp',iq0,iqc,iqb,nq)

c      print *,'calling evols - heavy...'

C--   Heavy quark evolution

      CALL EVOLCP('cp',iqc,iqb,nq)
      CALL EVOLBP('bp',iqb,nq)
c
      call getnset(iset)
      call save_pdfqcd(iset)
      return
*
      end
      subroutine EVLSEA(name,IQ0,IQC,IQB,NQGRI)
      implicit none
      
      CHARACTER*(*) name
      integer iq0,iqc,iqb,nqgri
      real*8 f34,f45,f43,f54,factor
      parameter(f34=1.d0/12.d0,f45=1.d0/20.d0,f43=-1.d0/12.d0,
     >          f54=-1.d0/20.d0)


      If(IQ0.le.IQC)then
         CAll EVPLUS(name,IQ0,1,IQC)
         CALL QADDSI(name,IQC,f34)
         CAll EVPLUS(name,IQC,IQC,IQB)
         CALL QADDSI(name,IQB,f45)
         CAll EVPLUS(name,IQB,IQB,NQGRI)
      else if(IQ0.gt.IQC.and.IQ0.le.IQB)then
         CAll EVPLUS(name,IQ0,IQC,IQB)
         CALL QADDSI(name,IQC,f43)
         CAll EVPLUS(name,IQC,1,IQC)
         CALL QADDSI(name,IQC,f34)
         CALL QADDSI(name,IQB,f45)
         CAll EVPLUS(name,IQB,IQB,NQGRI)
      else if(IQ0.gt.IQB)then
         CAll EVPLUS(name,IQ0,IQB,NQGRI)
         CALL QADDSI(name,IQB,f54)
         CAll EVPLUS(name,IQB,IQC,IQB)
         CALL QADDSI(name,IQB,f45)
         CALL QADDSI(name,IQC,f43)
         CAll EVPLUS(name,IQC,1,IQC)
         CALL QADDSI(name,IQC,f34)
      end if         

      end


      subroutine EVOLCP(name,IQC,IQB,NQGRI)
      implicit none
      
      CHARACTER*(*) name
      integer iq0,iqc,iqb,nqgri
      real*8 f4,f45
      parameter(f4=-1.d0/4.d0,f45=1.d0/20.d0)

c
c     First set to zero to avoid adding -1/4Singl at each iteration
c

      CAll QNPNUL(name)
      CALL QADDSI(name,IQC,f4)      
      CAll EVPLUS(name,IQC,IQC,IQB)
      CALL QADDSI(name,IQB,f45)
      CAll EVPLUS(name,IQB,IQB,NQGRI)
      
      end

      subroutine EVOLBP(name,IQB,NQGRI)
      implicit none
      
      CHARACTER*(*) name
      integer iqb,nqgri
      real*8 f5
      parameter(f5=-1.d0/5.d0)

c
c     First set to zero to avoid adding -1/5Singl at each iteration
c

      CAll QNPNUL(name)
      CALL QADDSI(name,IQB,f5)      
      CAll EVPLUS(name,IQB,IQB,NQGRI)
      
      end
c
c routine to swap in/out the array PDFQCD(MXX,MQ2,0:10) in QCDNUM as required 
c by nset
c
      subroutine get_pdfqcd(nset)
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c      
      include 'parmsetup.inc'     
c
      PARAMETER ( MXX = 410 )
      PARAMETER ( MQ2 =  120 )
C--   Do not set the following parameter to zero!
      PARAMETER ( NDFMAX = 20)
c
      COMMON/QCPASS/
     +ALPHA0, Q0ALFA, ASLAST, QALAST,
     +ALFASQ(MQ2), ALFAPQ(MQ2), ALFA2Q(MQ2),
     +DELUP(MQ2), DELDN(MQ2), PDFQCD(MXX,MQ2,0:10),
     +FNSQCD(MXX,MQ2),DNSQCD(MXX,MQ2),
     +FSIQCD(MXX,MQ2),DSIQCD(MXX,MQ2),
     +FGLQCD(MXX,MQ2),DGGQCD(MXX,MQ2),
     +FSTORE(MXX,MQ2,31:30+NDFMAX),IDFAST(7,30),NDFAST,
     +MARKFF(MXX,MQ2),MARKFH(MXX,MQ2),MARKQQ(MQ2),
     +ISTFID(31:30+NDFMAX),IPDFID(31:30+NDFMAX),IEALFA(MQ2),
     +IQL_LAST(10),IQ0_LAST(10),IQH_LAST(10)
c
      real*8  PDFSAV(MXX,MQ2,0:10,nmxset)
      save PDFSAV
c      
      do k=0,10
        do j=1,mq2
	  do i=1,mxx
	     PDFQCD(i,j,k)=PDFSAV(i,j,k,nset) 
	  enddo     
	enddo     
      enddo 
c
      return
c
      entry save_pdfqcd(nset)
c
      do k=0,10
        do j=1,mq2
	  do i=1,mxx
	     PDFSAV(i,j,k,nset)=PDFQCD(i,j,k) 
	  enddo     
	enddo     
      enddo 
c
      return
      end
Commit	Line	Data
4e9e3152	1	subroutine QCDNUMevolve(x,Q,pdf)
	2	implicit none
	3	character*64 gridname
	4	character*16 s1,dummy
	5	integer Eorder,index
	6	real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
	7	real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
	8	real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
	9	real*8 QALFAS,QPDFXQ,XFROMIX
	10	real*8 f(-6:6),pdf(-6:6)
	11	integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
	12	integer nx,nq
	13	integer nset,iset,isetlast
	14	data isetlast/-1/
	15	real*8 xmin,xmax,qmin,qmax,S
	16	save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
	17	save nx,xmin,xmax,nq,qmin,qmax,gridname,isetlast
	18	*
	19	call getnset(iset)
	20	c fix 14.7.2006 to speed up sue of QCDNUM
	21	if(iset.ne.isetlast) then
	22	call get_pdfqcd(iset)
	23	isetlast = iset
	24	endif
	25	c
	26	Q2=Q*Q
	27	nf=5
	28	if (Q2.lt.tb2) nf=4
	29	if (Q2.lt.tc2) nf=3
	30	pdf(0)=QPDFXQ('GLUON' ,x,Q2,iflag)
	31	singlet=QPDFXQ('SINGLET',x,Q2,iflag)
	32	dm= QPDFXQ('DM',x,Q2,IFLAG)
	33	um= QPDFXQ('UM',x,Q2,IFLAG)
	34	dp= QPDFXQ('DP',x,Q2,IFLAG)
	35	up= QPDFXQ('UP',x,Q2,IFLAG)
	36	sp= QPDFXQ('SP',x,Q2,IFLAG)
	37	ub=0.5d0*(up-um+singlet/dble(nf))
	38	db=0.5d0*(dp-dm+singlet/dble(nf))
	39	sb=0.5d0*(sp+singlet/dble(nf))
	40	cb=0.d0
	41	if (nf.ge.4) then
	42	cp= QPDFXQ('CP',X,Q2,IFLAG)
	43	cb=0.5d0*(cp+singlet/dble(nf))
	44	end if
	45	bb=0.d0
	46	if (nf.ge.5) then
	47	bp= QPDFXQ('BP',X,Q2,IFLAG)
	48	bb=0.5d0*(bp+singlet/dble(nf))
	49	end if
	50
	51	pdf(1)=dm+db
	52	pdf(2)=um+ub
	53	pdf(3)=sb
	54	pdf(4)=cb
	55	pdf(5)=bb
	56	pdf(6)=0d0
	57	pdf(-1)=db
	58	pdf(-2)=ub
	59	pdf(-3)=sb
	60	pdf(-4)=cb
	61	pdf(-5)=bb
	62	pdf(-6)=0d0
	63	return
	64	*
65	entry QCDNUMalfa(alfas,Q)
66	Q2=Q*Q
67	nf=6
68	if (Q2.lt.mt2) nf=5
69	if (Q2.lt.mb2) nf=4
70	if (Q2.lt.mc2) nf=3
71	alfas=QALFAS(Q2,Qlam,nf,iflag)
72	return
73	*
74	entry QCDNUMread(nset)
75	c read(1,*) s1
76	c print *,s1
77	c gridname='large.grid'
78	c nx=400
79	c xmin=1d-6
80	c xmax=1d0
81	c nq=112
82	c qmin=1d0
83	c qmax=1d10
84	c read(1,*) dummy
85	* else
86	read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
87	c* endif
88	c iset = nset
89	return
90	*
91	entry QCDNUMinit(nset,Eorder,Q2fit)
92	call qninit
93	call QNISET('ORDER',Eorder+1)
94	call grxdef(nx,xmin)
95	call grqdef(nq,qmin,qmax)
96	call grqinp(Q2fit,1)
97	call getQMassM(nset,4,mc)
98	mc2=mc*mc
99	CALL QNRSET('CMASS',mc)
100	CALL QNRSET('MCALF',mc)
101	call getThresholdM(nset,4,tc)
102	tc2=tc*tc
103	call GRQINP(tc2,1)
104	CALL GRQINP(tc2-1.0d-4,1)
105	call getQMassM(nset,5,mb)
106	mb2=mb*mb
107	CALL QNRSET('BMASS',mb)
108	CALL QNRSET('MBALF',mb)
109	call getThresholdM(nset,5,tb)
110	tb2=tb*tb
111	call GRQINP(tb2,1)
112	CALL GRQINP(tb2-1.0d-4,1)
113	call getQMassM(nset,6,mt)
114	mt2=mt*mt
115	CALL QNRSET('TMASS',mt)
116	CALL QNRSET('MTALF',mt)
117	call qthres(tc2,tb2)
118	iq0=IQFROMQ(Q2fit)
119	iqc=IQFROMQ(tc2)
120	iqb=IQFROMQ(tb2)
121	nf0= NFLGET(iq0)
122	c print *,nf0
123	CALL QNBOOK(2,'dm')
124	CALL QNBOOK(3,'um')
125	CALL QNBOOK(4,'dp')
126	CALL QNBOOK(5,'up')
127	CALL QNBOOK(6,'sp')
128	CALL QNBOOK(7,'cp')
129	CALL QNBOOK(8,'bp')
130	CALL QNLSET('W1ANA',.TRUE.)
131	CALL QNLSET('W2NUM',.TRUE.)
132	CALL QNLSET('W2STF',.FALSE.)
133	if (index(gridname,'none').eq.1) then
134	call qnfilw(0,0)
135	else
136	qnerr=-1
137	open(unit=2,status='old',file=gridname,
138	. form='unformatted',err=1)
139	call QNREAD(2,1,qnerr)
140	1 close(2)
141	if (qnerr.ne.0) then
142	write(,) 'Grid file problem: ',gridname
143	if (qnerr.lt.0) then
144	write(,) 'Grid file does not exist'
145	write(,) 'Calculating and creating grid file'
146	call qnfilw(0,0)
147	open(unit=2,status='unknown',file=gridname,
148	. form='unformatted')
149	call QNDUMP(2)
150	close(2)
151	else
152	write(,) 'Existing grid file is inconsistent'
153	if (qnerr.eq.1)
154	. write(,) 'Defined grid different'
155	if (qnerr.eq.2)
156	. write(,) 'Heavy quark weight table different'
157	if (qnerr.eq.3)
158	. write(,) 'Charm mass different'
159	if (qnerr.eq.4)
160	. write(,) 'Bottom mass different'
161	stop
162	endif
163	endif
164	endif
165	return
166	*
167	entry QCDNUMpdf(nset)
168	c print *,'entering QCDNUMpdf',nset
169	call GetAlfas(nset,alfas0,scale0)
170	c print *,alfas0,scale0
171	Q2=scale0*scale0
172	CALL QNRSET('ALFAS',alfas0)
173	CALL QNRSET('ALFQ0',Q2)
174	DO ix = 1,nx
175	xx = XFROMIX(ix)
176	c print *,'calling parmPDF',ix
177	call parmPDF(nset,xx,f)
178	c if(ix.lt.6) print *,nset,xx,f
179	singlet=0.d0
180	do i=1,nf0
181	singlet=singlet+f(i)+f(-i)
182	c print *,i,singlet
183	end do
184	gl=f(0)
185	dm=f(1)-f(-1)
186	um=f(2)-f(-2)
187	dp=f(1)+f(-1)-singlet/dble(nf0)
188	up=f(2)+f(-2)-singlet/dble(nf0)
189	sp=f(3)+f(-3)-singlet/dble(nf0)
190	c print *,'calling qpnset sg...',ix,iq0,singlet
191	CALL QNPSET('SINGLET',ix,iq0,singlet)
192	c print *,'calling qpnset gl...',ix,iq0,gl
193	CALL QNPSET('GLUON',ix,iq0,gl)
194	c print *,'calling qpnset dm...',ix,iq0,dm
195	CALL QNPSET('DM',ix,iq0,DM)
196	c print *,'calling qpnset um...',ix,iq0,um
197	CALL QNPSET('UM',ix,iq0,UM)
198	c print *,'calling qpnset dp...',ix,iq0,dp
199	CALL QNPSET('DP',ix,iq0,DP)
200	c print *,'calling qpnset up...',ix,iq0,up
201	CALL QNPSET('UP',ix,iq0,UP)
202	c print *,'calling qpnset sp...',ix,iq0,sp
203	CALL QNPSET('SP',ix,iq0,SP)
204	ENDDO
205	c print *,'calling evols...'
206	CALL EVOLSG(iq0,1,nq)
207	c print *,'calling evols dm...'
208	CALL EVOLNM('DM',iq0,1,nq)
209	c print *,'calling evols um...'
210	CALL EVOLNM('UM',iq0,1,nq)
211	c print *,'calling evols dp...',iq0,iqc,iqb,nq
212	CALL EVLSEA('dp',iq0,iqc,iqb,nq)
213	c print *,'calling evols up...'
214	CALL EVLSEA('up',iq0,iqc,iqb,nq)
215	c print *,'calling evols sp...'
216	CALL EVLSEA('sp',iq0,iqc,iqb,nq)
217
218	c print *,'calling evols - heavy...'
219
220	C-- Heavy quark evolution
221
222	CALL EVOLCP('cp',iqc,iqb,nq)
223	CALL EVOLBP('bp',iqb,nq)
224	c
225	call getnset(iset)
226	call save_pdfqcd(iset)
227	return
228	*
229	end
230	subroutine EVLSEA(name,IQ0,IQC,IQB,NQGRI)
231	implicit none
232
233	CHARACTER() name
234	integer iq0,iqc,iqb,nqgri
235	real*8 f34,f45,f43,f54,factor
236	parameter(f34=1.d0/12.d0,f45=1.d0/20.d0,f43=-1.d0/12.d0,
237	> f54=-1.d0/20.d0)
238
239
240	If(IQ0.le.IQC)then
241	CAll EVPLUS(name,IQ0,1,IQC)
242	CALL QADDSI(name,IQC,f34)
243	CAll EVPLUS(name,IQC,IQC,IQB)
244	CALL QADDSI(name,IQB,f45)
245	CAll EVPLUS(name,IQB,IQB,NQGRI)
246	else if(IQ0.gt.IQC.and.IQ0.le.IQB)then
247	CAll EVPLUS(name,IQ0,IQC,IQB)
248	CALL QADDSI(name,IQC,f43)
249	CAll EVPLUS(name,IQC,1,IQC)
250	CALL QADDSI(name,IQC,f34)
251	CALL QADDSI(name,IQB,f45)
252	CAll EVPLUS(name,IQB,IQB,NQGRI)
253	else if(IQ0.gt.IQB)then
254	CAll EVPLUS(name,IQ0,IQB,NQGRI)
255	CALL QADDSI(name,IQB,f54)
256	CAll EVPLUS(name,IQB,IQC,IQB)
257	CALL QADDSI(name,IQB,f45)
258	CALL QADDSI(name,IQC,f43)
259	CAll EVPLUS(name,IQC,1,IQC)
260	CALL QADDSI(name,IQC,f34)
261	end if
262
263	end
264
265
266	subroutine EVOLCP(name,IQC,IQB,NQGRI)
267	implicit none
268
269	CHARACTER() name
270	integer iq0,iqc,iqb,nqgri
271	real*8 f4,f45
272	parameter(f4=-1.d0/4.d0,f45=1.d0/20.d0)
273
274	c
275	c First set to zero to avoid adding -1/4Singl at each iteration
276	c
277
278	CAll QNPNUL(name)
279	CALL QADDSI(name,IQC,f4)
280	CAll EVPLUS(name,IQC,IQC,IQB)
281	CALL QADDSI(name,IQB,f45)
282	CAll EVPLUS(name,IQB,IQB,NQGRI)
283
284	end
285
286	subroutine EVOLBP(name,IQB,NQGRI)
287	implicit none
288
289	CHARACTER() name
290	integer iqb,nqgri
291	real*8 f5
292	parameter(f5=-1.d0/5.d0)
293
294	c
295	c First set to zero to avoid adding -1/5Singl at each iteration
296	c
297
298	CAll QNPNUL(name)
299	CALL QADDSI(name,IQB,f5)
300	CAll EVPLUS(name,IQB,IQB,NQGRI)
301
302	end
303	c
304	c routine to swap in/out the array PDFQCD(MXX,MQ2,0:10) in QCDNUM as required
305	c by nset
306	c
307	subroutine get_pdfqcd(nset)
308	c
309	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
310	c
311	include 'parmsetup.inc'
312	c
313	PARAMETER ( MXX = 410 )
314	PARAMETER ( MQ2 = 120 )
315	C-- Do not set the following parameter to zero!
316	PARAMETER ( NDFMAX = 20)
317	c
318	COMMON/QCPASS/
319	+ALPHA0, Q0ALFA, ASLAST, QALAST,
320	+ALFASQ(MQ2), ALFAPQ(MQ2), ALFA2Q(MQ2),
321	+DELUP(MQ2), DELDN(MQ2), PDFQCD(MXX,MQ2,0:10),
322	+FNSQCD(MXX,MQ2),DNSQCD(MXX,MQ2),
323	+FSIQCD(MXX,MQ2),DSIQCD(MXX,MQ2),
324	+FGLQCD(MXX,MQ2),DGGQCD(MXX,MQ2),
325	+FSTORE(MXX,MQ2,31:30+NDFMAX),IDFAST(7,30),NDFAST,
326	+MARKFF(MXX,MQ2),MARKFH(MXX,MQ2),MARKQQ(MQ2),
327	+ISTFID(31:30+NDFMAX),IPDFID(31:30+NDFMAX),IEALFA(MQ2),
328	+IQL_LAST(10),IQ0_LAST(10),IQH_LAST(10)
329	c
330	real*8 PDFSAV(MXX,MQ2,0:10,nmxset)
331	save PDFSAV
332	c
333	do k=0,10
334	do j=1,mq2
335	do i=1,mxx
336	PDFQCD(i,j,k)=PDFSAV(i,j,k,nset)
337	enddo
338	enddo
339	enddo
340	c
341	return
342	c
343	entry save_pdfqcd(nset)
344	c
345	do k=0,10
346	do j=1,mq2
347	do i=1,mxx
348	PDFSAV(i,j,k,nset)=PDFQCD(i,j,k)
349	enddo
350	enddo
351	enddo
352	c
353	return
354	end
355
356