Version update.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / alphas.f
1       real*8 function alphasPDF(Q)
2       implicit none
3       integer nset
4       real*8 Q,a
5       call getnset(nset)
6       call evolveAs(nset,Q,a)
7       alphasPDF=a
8       return
9       end
10 *
11       real*8 function alphasPDFM(nset,Q)
12       implicit none
13       integer nset
14       real*8 Q,a
15       call evolveAs(nset,Q,a)
16       alphasPDFM=a
17       return
18       end
19 *
20       subroutine GetQmass(nnf,mass)
21       implicit none
22       real*8 mass
23       integer nnf,order,nset
24       nset = 1
25       call GetQmassM(nset,nnf,mass)
26       return
27       entry GetOrderAs(order)
28       nset = 1
29       call GetOrderAsM(nset,order)
30       return
31       end
32  
33  
34       subroutine evolveAs(nset,Q,alphas)
35       implicit none
36       include 'parmsetup.inc'
37       integer iset,imem
38       common/SET/iset,imem
39       integer nset
40       character*16 s1,s2,s3
41       real*8 Q,alphas,alfas0,scale0
42       real*8 Q0(nmxset)
43       real*8 AlfasQ(nmxset)
44       real*8 b0,b1,b2,L,As,mass
45       double precision CtLhALPI,CtLhAlphaNew
46       parameter (b0=1.2202,b1=0.4897,b2=0.1913)
47       integer n
48       integer order
49       integer EvlOrd(nmxset)
50       integer parm(nmxset)
51       integer Etype(nmxset)
52       integer Method(nmxset)
53       integer nnf,naf
54       real*8 cmass,bmass,tmass
55       common/masses_LHA/cmass,bmass,tmass
56       save EvlOrd,parm,Etype,alfasQ,Q0,Method
57
58       call setnset(nset)
59 c
60       if (method(nset).eq.0) then
61          if ((Etype(nset).eq.1).or.(Etype(nset).eq.2)) then
62             L=log(Q/Q0(nset))
63             As=alfasQ(nset)
64             if (Etype(nset).eq.2) call GetParmPDF(nset,parm,As)
65 c           print *,Etype,parm,As
66             if (EvlOrd(nset).eq.0) L=b0*L
67             if (EvlOrd(nset).eq.1) L=(b0+As*b1)*L
68             if (EvlOrd(nset).eq.2) L=(b0+As*b1+As**2*b2)*L
69      .                        -0.5*As**2*b0*b1/2d0*L**2
70             alphas=As/(1.0+As*L)
71          endif
72       endif
73       if (method(nset).eq.1) then
74          call alfasevolve(nset,alphas,Q)
75       elseif (method(nset).eq.2) then
76          call alphamrs(5,alphas,Q)
77       elseif (method(nset).eq.3) then
78          alphas = 3.1415926535898d0*CtLhALPI(Q)
79       elseif (method(nset).eq.4) then
80          alphas = CtLhAlphaNew(Q)
81       elseif (method(nset).eq.5) then
82          call alphamrs(3,alphas,Q)
83       elseif (method(nset).eq.6) then
84          call alphamrs(4,alphas,Q)
85       elseif (method(nset).eq.7) then
86          call alphacteq5f34(3,alphas,Q)
87       elseif (method(nset).eq.8) then
88          call alphacteq5f34(4,alphas,Q)
89       endif
90       return
91 *
92       entry GetQmassM(nset,nnf,mass)
93       n=abs(nnf)
94       mass=0d0
95       if (n.eq.4) mass=cmass
96       if (n.eq.5) mass=bmass
97       if (n.eq.6) mass=tmass
98       return
99     
100       entry GetNactive(naf,q)
101 c compute nnfn = number of active flavors at scale q.
102       n = 3
103       if(q .ge. cmass) n = 4
104       if(q .ge. bmass) n = n + 1
105       if(q .ge. tmass) n = n + 1
106       naf = n
107       return
108
109       entry GetAlfas(nset,alfas0,scale0)
110       scale0=Q0(nset)
111       alfas0=alfasQ(nset)
112       if (Etype(nset).eq.2) call GetParmPDF(nset,parm(nset),alfas0)
113       return
114 *
115       entry GetOrderAsM(nset,order)
116       order=EvlOrd(nset)
117       return
118 *
119       entry InitAlphasPDF(nset)
120       Etype(nset)=-1
121       EvlOrd(nset)=-1
122       read(1,*) s1,s2,s3
123       if (index(s2,'lo').eq.1) EvlOrd(nset)=0
124       if (index(s2,'nlo').eq.1) EvlOrd(nset)=1
125       if (index(s2,'nnlo').eq.1) EvlOrd(nset)=2
126       if (EvlOrd(nset).lt.0) then
127          write(*,*) 'File description error:'
128          write(*,*) 'Unknown alpha_s evolution order ',s2
129          stop
130       endif
131       if (index(s1,'Fixed').eq.1) then
132          Etype(nset)=1
133          parm(nset)=-1
134          read(1,*) alfasQ(nset),Q0(nset),cmass,bmass,tmass
135       endif
136       if (index(s1,'Variable').eq.1) then
137          Etype(nset)=2
138          alfasQ(nset)=0d0
139          read(1,*) parm(nset),Q0(nset),cmass,bmass,tmass
140       endif
141       if (Etype(nset).lt.0) then
142          write(*,*) 'File description error:'
143          write(*,*) 'Unknown alpha_s evolution method ',s1
144          stop
145       endif
146       Method(nset)=-1
147       if (index(s3,'Internal').eq.1) Method(nset)=0
148       if (index(s3,'EvolCode').eq.1) Method(nset)=1
149       if (index(s3,'MRSTalfa').eq.1) Method(nset)=2
150       if (index(s3,'CTEQalfa').eq.1) Method(nset)=3
151       if (index(s3,'CTEQABalfa').eq.1) Method(nset)=4
152       if (index(s3,'MRST3alfa').eq.1) Method(nset)=5
153       if (index(s3,'MRST4alfa').eq.1) Method(nset)=6
154       if (index(s3,'CTEQ5F3alfa').eq.1) Method(nset)=7
155       if (index(s3,'CTEQ5F4alfa').eq.1) Method(nset)=8
156       if (Method(nset).lt.0) then
157          write(*,*) 'File description error:'
158          write(*,*) 'Unknown alpha_s method ',s3
159          stop
160       endif
161 c      print *,'alpha_s evolution method = ',method
162       return
163       end
164  
165       subroutine alphamrs(nflav,alpha,Q)
166       IMPLICIT real*8 (a-h,o-z)
167       common/SET/iset,imem
168       common/masses_LHA/cmass,bmass,tmass
169       dimension parms(31)
170       DATA PI/3.14159/
171       DATA TOL/.0005/
172       data memold/-1/
173       integer nset
174       save alambda,memold 
175 c
176 c       print *,'calling alphamrs with ',nflav,Q
177        call getnset(nset)
178 c       
179 c
180 c      alambda=0.323
181 c      alambda=0.3342
182 c -- find alambda corresponding to alphas given in first parameter
183       call getnmem(nset,imem)
184       if(imem.ne.memold) then
185         call listPDF(nset,imem,parms)
186         alphas = parms(1)
187 c       print *,alphas
188         memold=imem
189         qz2=8315.
190         qz = dsqrt(qz2)
191         alambda = 0.3000
192         astep = 0.010
193         tol2 = 0.0000001
194         idir = +1
195   10    continue
196         alambda = alambda + idir*astep
197         call mrslambda(nflav,alpha,qz,alambda)
198
199           if(idir*(alphas-alpha).gt.0.0) then
200              go to 20
201           else
202              astep = 0.5*astep
203              idir = -1*idir
204              go to 20
205           endif  
206   20   continue
207         if(abs(alpha-alphas).gt.tol2) go to 10
208       endif
209 c - alambda found  -- save it !!!!
210 c - next call mrslambda to get alphas at q with the correct alambda
211       call mrslambda(nflav,alpha,q,alambda)
212       RETURN
213       END
214 *
215       subroutine rgras(alpha,Q2)
216       IMPLICIT real*8 (a-h,o-z)
217       real*8 mc,mb,mt
218       include 'parmsetup.inc'
219       character*16 name(nmxset)
220       integer nmem(nmxset),ndef(nmxset),mem
221       common/NAME/name,nmem,ndef,mem
222       q=dsqrt(q2)
223       call getnset(nset)
224       nflav=5
225       if(name(nset).eq.'QCDNUM_MRST3') nflav=3     
226       if(name(nset).eq.'QCDNUM_MRST4') nflav=4     
227       call alphamrs(nflav,alpha,Q)
228       return
229       end
230 *
231 * new vewrsion of mrslambda 13/5/2004 - includes lo, nlo, and nnlo with flags
232 *
233       subroutine mrslambda(nflav,alpha,Q,alambda)
234       IMPLICIT REAL*8(A-H,O-Z)
235       DATA PI/3.14159/
236       DATA TOL/.0005/
237 c  The value of Lambda required corresponds to nflav=4
238 c  iord=0 gives leading order, iord=1 gives NLO, iord=2 gives NNLO
239       qsdt=8.18    !!  This is the value of 4m_c^2
240       qsct=74.0    !!  This is the value of 4m_b^2
241 c      print *,'calling mrslambda with ',nflav,Q,alambda
242       al2=alambda*alambda
243       q2=q*q
244       t=dlog(q2/al2)
245 c
246 c next line to be replaced by a getiord in the final version
247 c      iord=1
248 c
249       call getnset(nset)
250       call GetOrderAsM(nset,iord)
251       ITH=0
252       TT=T
253       qsdtt=qsdt/4.
254       qsctt=qsct/4.
255       AL=ALAMBDA
256       AL2=AL*AL
257       FLAV=4.
258       if(nflav.eq.3) flav=3.
259       QS=AL2*dEXP(T)
260
261       if(qs.lt.0.5d0) then   !!  running stops below 0.5
262           qs=0.5d0
263           t=dlog(qs/al2)
264           tt=t
265       endif
266    
267 c      if(QS.lt.QSDTT.and.nflav.eq.4) then
268 c        flav=3.
269 c       go to 11
270 c      endif    
271         
272
273       IF(QS.gt.QSCTT.and.nflav.gt.4) GO TO 12  
274       IF(QS.lt.QSDTT.and.nflav.gt.3) GO TO 312  
275    11 CONTINUE
276       B0=11-2.*FLAV/3. 
277       X1=4.*PI/B0
278       IF(IORD.eq.0) then
279       ALPHA=X1/T
280       ELSE
281           if(iord.gt.1) then
282           alpha=qwikalf(t,iord,flav)
283           go to 51
284           endif 
285       B1=102.-38.*FLAV/3.
286       X2=B1/B0**2
287       AS2=X1/T*(1.-X2*dLOG(T)/T)
288     5 AS=AS2
289       F=-T+X1/AS-X2*dLOG(X1/AS+X2)
290       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
291       AS2=AS-F/FP
292       DEL=ABS(F/FP/AS)
293       IF((DEL-TOL).GT.0.) go to 5
294       ALPHA=AS2
295    51 continue
296       ENDIF
297       IF(ITH.EQ.0) RETURN
298       GO TO (13,14,15) ITH
299 c      GO TO 5
300    12 ITH=1
301       T=dLOG(QSCTT/AL2)
302       GO TO 11
303    13 ALFQC4=ALPHA
304       FLAV=5.   
305       ITH=2
306
307       GO TO 11
308    14 ALFQC5=ALPHA
309       ITH=3
310       T=TT
311       GO TO 11
312    15 ALFQS5=ALPHA
313       ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
314       ALPHA=1./ALFINV
315       RETURN
316
317   311 CONTINUE
318       B0=11-2.*FLAV/3. 
319       X1=4.*PI/B0
320       IF(IORD.eq.0) then
321       ALPHA=X1/T
322
323       ELSE
324           if(iord.gt.1) then
325           alpha=qwikalf(t,iord,flav)
326           go to 351
327           endif 
328       B1=102.-38.*FLAV/3.
329       X2=B1/B0**2
330       AS2=X1/T*(1.-X2*dLOG(T)/T)
331    35 AS=AS2
332       F=-T+X1/AS-X2*dLOG(X1/AS+X2)
333       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
334       AS2=AS-F/FP
335       DEL=ABS(F/FP/AS)
336       IF((DEL-TOL).GT.0.) go to 35
337       ALPHA=AS2
338
339   351 continue
340       endif
341       IF(ITH.EQ.0) RETURN
342       GO TO (313,314,315) ITH
343   312 ITH=1
344       T=dLOG(QSDTT/AL2)
345       GO TO 311
346   313 ALFQC4=ALPHA
347       FLAV=3.   
348       ITH=2
349       GO TO 311
350   314 ALFQC3=ALPHA
351       ITH=3
352       T=TT
353       GO TO 311
354   315 ALFQS3=ALPHA
355       ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
356       ALPHA=1./ALFINV
357       RETURN
358       END
359
360
361
362       double precision function qwikalf(t,iord,flav)
363       implicit real*8(a-h,o-z)
364       dimension z3(6),z4(6),z5(6),zz3(6),zz4(6),zz5(6)
365       data z3/ -.161667E+01,0.954244E+01,
366      .0.768623E+01,0.101523E+00,-.360127E-02,0.457867E-04/
367       data z4/ -.172239E+01,0.831185E+01,
368      .0.721463E+01,0.835531E-01,-.285436E-02,0.349129E-04/
369       data z5/ -.872190E+00,0.572816E+01,
370      .0.716119E+01,0.195884E-01,-.300199E-03,0.151741E-05/
371       data zz3/-.155611E+02,0.168406E+02,
372      .0.603014E+01,0.257682E+00,-.970217E-02,0.127628E-03/
373       data zz4/-.106762E+02,0.118497E+02,0.664964E+01,
374      .0.112996E+00,-.317551E-02,0.302434E-04/
375       data zz5/-.531860E+01,0.708503E+01,0.698352E+01,
376      .0.274170E-01,-.426894E-03,0.217591E-05/
377
378       data pi/3.14159/
379       nfm2=flav-2.
380       x=dsqrt(t)
381       x2=x*x
382       x3=x*x2
383       x4=x*x3
384       x5=x*x4
385       go to (1,2) iord
386     1 go to (3,4,5) nfm2
387     3 y=z3(1)+z3(2)*x+z3(3)*x2+z3(4)*x3+z3(5)*x4+z3(6)*x5
388       go to 10
389     4 y=z4(1)+z4(2)*x+z4(3)*x2+z4(4)*x3+z4(5)*x4+z4(6)*x5
390       go to 10
391     5 y=z5(1)+z5(2)*x+z5(3)*x2+z5(4)*x3+z5(5)*x4+z5(6)*x5
392       go to 10
393     2 go to (6,7,8) nfm2
394     6 y=zz3(1)+zz3(2)*x+zz3(3)*x2+zz3(4)*x3+zz3(5)*x4+zz3(6)*x5
395       go to 10
396     7 y=zz4(1)+zz4(2)*x+zz4(3)*x2+zz4(4)*x3+zz4(5)*x4+zz4(6)*x5
397       go to 10
398     8 y=zz5(1)+zz5(2)*x+zz5(3)*x2+zz5(4)*x3+zz5(5)*x4+zz5(6)*x5
399       go to 10
400    10 qwikalf=4.*pi/y
401       return
402       end
403 c=====================================================================
404 c next is alphas routine from PDFLIB
405 c
406 c      DOUBLE PRECISION FUNCTION ALPHAS2(SCALE,qcdl5)
407       subroutine aspdflib(alphas2,SCALE,iord,qcdl5)
408       implicit real*8 (a-h,o-z)
409       double precision 
410      +    NF
411 c      CHARACTER*20 PARM(NCHDIM)
412 c      double precision 
413 c     +       VAL(NCHDIM)
414 c      DATA XMC/1.5D0/,XMB/4.75D0/,XMT/100.D0/
415       DATA XMC/1.43D0/,XMB/4.30D0/,XMT/100.D0/
416       DATA ZEROD/0.D0/,PONED/0.001D0/,ONED/1.D0/,TWOD/2.D0/
417      
418 c      QCDL5=0.165d0
419       
420 c      print *,qcdl5,iord
421       LO = 1
422       if(iord.ne.0) LO = 2
423 c      if(LO.eq.1) then 
424 c         print *,'Calculating ALPHA_S at Leading Order',LO
425 c      else
426 c         print *,'Calculating ALPHA_S at Next to Leading Order',LO
427 c      endif
428
429       TMAS = 180.0d0
430     
431
432       ALPHAS2 = ZEROD
433 C.
434         PI=4.0D0*ATAN(ONED)
435         B6  = (33.D0-2.D0*6.D0)/PI/12.D0
436         BP6 = (153.D0 - 19.D0*6.D0) / PI / TWOD / (33.D0 - 2.D0*6.D0)
437         B5  = (33.D0-2.D0*5.D0)/PI/12.D0
438         BP5 = (153.D0 - 19.D0*5.D0) / PI / TWOD / (33.D0 - 2.D0*5.D0)
439         B4  = (33.D0-2.D0*4.D0)/PI/12.D0
440         BP4 = (153.D0 - 19.D0*4.D0) / PI / TWOD / (33.D0 - 2.D0*4.D0)
441         B3  = (33.D0-2.D0*3.D0)/PI/12.D0
442         BP3 = (153.D0 - 19.D0*3.D0) / PI / TWOD / (33.D0 - 2.D0*3.D0)
443         XLC = TWOD * LOG( XMC/QCDL5)
444         XLB = TWOD * LOG( XMB/QCDL5)
445         XLT = TWOD * LOG( XMT/QCDL5 * TMAS/XMT)
446         XLLC = LOG( XLC)
447         XLLB = LOG( XLB)
448         XLLT = LOG( XLT)
449         C65  =  ONED/( ONED/(B5 * XLT) - XLLT*BP5/(B5 * XLT)**2 )
450      +        - ONED/( ONED/(B6 * XLT) - XLLT*BP6/(B6 * XLT)**2 )
451         C45  =  ONED/( ONED/(B5 * XLB) - XLLB*BP5/(B5 * XLB)**2 )
452      +        - ONED/( ONED/(B4 * XLB) - XLLB*BP4/(B4 * XLB)**2 )
453         C35  =  ONED/( ONED/(B4 * XLC) - XLLC*BP4/(B4 * XLC)**2 )
454      +        - ONED/( ONED/(B3 * XLC) - XLLC*BP3/(B3 * XLC)**2 ) + C45
455 C
456       Q   = SCALE
457       XLQ = TWOD *  LOG( Q/QCDL5 )
458       XLLQ =  LOG( XLQ )
459      
460 c      IF  ( NF .LT. ZEROD) THEN
461         IF      ( Q .GT. XMT * TMAS/XMT) THEN
462           NF = 6.D0
463         ELSEIF  ( Q .GT. XMB ) THEN
464           NF = 5.D0
465         ELSEIF  ( Q .GT. XMC ) THEN
466           NF = 4.D0
467         ELSE
468           NF = 3.D0
469         ENDIF
470 c      ENDIF
471 c      print *,nf
472       IF(NF .GT. 6.D0) NF = 6.D0
473       IF      ( NF .EQ. 6.D0 ) THEN
474        ALF = ONED/(ONED/(ONED/(B6*XLQ)- BP6/(B6*XLQ)**2*XLLQ) + C65)
475        IF (LO.EQ.1) ALF = ONED/B6/XLQ
476       ELSEIF  ( NF .EQ. 5.D0 ) THEN
477        ALF = ONED/(B5 * XLQ) -  BP5/(B5 * XLQ)**2 * XLLQ
478        IF (LO.EQ.1) ALF = ONED/B5/XLQ
479       ELSEIF  ( NF .EQ. 4.D0 ) THEN
480        ALF = ONED/(ONED/(ONED/(B4*XLQ)- BP4/(B4*XLQ)**2*XLLQ) + C45)
481        IF (LO.EQ.1) ALF = ONED/B4/XLQ
482       ELSEIF  ( NF .EQ. 3.D0 ) THEN
483        ALF = ONED/(ONED/(ONED/(B3*XLQ)- BP3/(B3*XLQ)**2*XLLQ) + C35)
484        IF (LO.EQ.1) ALF = ONED/B3/XLQ
485       ELSE
486        WRITE(*,*)'Error in Alphas2'
487        STOP
488       ENDIF
489       ALPHAS2 = ALF
490       RETURN
491       END
492 c ========================================================================
493       SUBROUTINE CtLhAlphaNewSET(MC,MB,MT,Q0,ALPHA0,IORDER,IMODE)
494 c call to set quark masses for alpha_s, and choose lambda or its 
495 c equivalent to make alpha_s take the value alpha0 at scale Q0.
496
497       IMPLICIT NONE
498       DOUBLE PRECISION MC, MB, MT, Q0, ALPHA0
499       INTEGER IORDER, IMODE
500
501       DOUBLE PRECISION UDSCBT, QQ0, AALPHA0
502       INTEGER IIORDER, IIMODE
503       COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
504       COMMON /ALSCALE/ QQ0, AALPHA0
505
506       double precision Adummy
507       integer Nfl, Idummy
508       common / QCDtable /  Adummy, Nfl, Idummy
509
510         if((imode .lt. 1) .or. (imode .gt. 3)) then
511            print *,'CtLhAlphaNewSET: fatal imode=',imode
512            stop
513         endif
514
515         IIMODE = IMODE
516         QQ0 = Q0
517         AALPHA0 = ALPHA0
518         IIORDER = IORDER
519
520         UDSCBT(1) = .005D0
521         UDSCBT(2) = .010D0
522         UDSCBT(3) = .300D0
523         UDSCBT(4) = MC
524         UDSCBT(5) = MB
525         UDSCBT(6) = MT
526
527 c set artificial quark masses, if necessary, in alpha_s to enforce 
528 c the requested maximum number of flavors...
529         if(Nfl .le. 5) UDSCBT(6) = 6.d99
530         if(Nfl .le. 4) UDSCBT(5) = 5.d99
531         if(Nfl .le. 3) UDSCBT(4) = 4.d99
532         if(Nfl .le. 2) UDSCBT(3) = 3.d99
533         if(Nfl .le. 1) UDSCBT(2) = 2.d99
534         if(Nfl .le. 0) UDSCBT(1) = 1.d99
535
536 c       print *,'CtLhAlphaNewSET: imode=',imode,' Nfl=',Nfl             !*** temporary ***
537 c       print *,'CtLhAlphaNewSET: mc,mb,mt=', mc, mb, mt                !*** temporary ***
538 c       print *,'CtLhAlphaNewSET: Q0, alpha0=',Q0, alpha0               !*** temporary ***
539
540
541       RETURN
542       END
543
544       FUNCTION CtLhAlphaNew(Q)
545       IMPLICIT NONE
546       DOUBLE PRECISION CtLhAlphaNew, Q
547
548       DOUBLE PRECISION Q2, CtLhQALPHAS, UDSCBT
549       integer nf,ier
550
551       INTEGER IIMODE, IIORDER
552       COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
553
554         if((iimode .ge. 1) .and. (iimode .le. 3)) then
555            Q2=Q*Q
556            nf=6
557            if (Q2.lt. UDSCBT(6)**2) nf=5
558            if (Q2.lt. UDSCBT(5)**2) nf=4
559            if (Q2.lt. UDSCBT(4)**2) nf=3
560            if (Q2.lt. UDSCBT(3)**2) nf=2
561            if (Q2.lt. UDSCBT(2)**2) nf=1
562            if (Q2.lt. UDSCBT(1)**2) nf=0
563            
564 c external maximum number of flavors -- typically, nfmax=5 so 
565 c top quark is not a parton even at large Q...
566
567            CtLhAlphaNew=CtLhQALPHAS(Q2,nf,ier)
568            if(ier .ne. 0) then
569               print *,'warning in CtLhAlphaNew, Q=',Q,' nf=',nf,
570      &           ' ier=',ier,' CtLhAlphaNew=',CtLhAlphaNew
571            endif
572
573         else
574            print *,'CtLhAlphaNew: undefined mode=',iimode 
575            stop
576         endif
577
578
579       return
580       end
581
582       FUNCTION CtLhQALPHAS(QQ2,NF,IERR)
583       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
584  
585       COMMON /ALSCALE/ Q0, ALPHA0
586       COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
587
588       Q02 = Q0**2
589       ALP0 = ALPHA0
590       IOR = IIORDER
591
592       NFFF = NF
593       CtLhQALPHAS = CtLhA0TOA1(QQ2,Q02,ALP0,IOR,NFFF,IERR)
594
595       RETURN
596       END
597
598       FUNCTION CtLhA0TOA1(QSU,QS0,AS0,IORD,NFF,IERR)
599       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
600  
601       COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
602
603       QS1   = QSU
604         
605       QMU0  = SQRT(QS0)
606       QMU1  = SQRT(QS1)
607
608       DO I = 1, 6
609          IF(QMU0.GE.UDSCBT(I)) NF0 = I
610          IF(QMU1.GE.UDSCBT(I)) NF1 = I
611       ENDDO
612  
613       IF(NF1.LT.NF0) THEN
614         IST = -1
615         JST =  0
616       ELSE
617         IST = 1
618         JST = 1
619       ENDIF
620  
621       ALFA0 = AS0
622       Q00   = QS0
623  
624       IERR = 0
625       DO 50 NF = NF0,NF1,IST
626  
627       IF(NF.NE.NF1) THEN
628         Q21 = UDSCBT(NF+JST)**2
629       ELSE
630         Q21 = QS1
631       ENDIF
632       IMODE = IIMODE
633       ALFA1 = CtLhALPHAR(Q21,Q00,ALFA0,NF,IORD,IMODE,JERR)
634       IERR = IERR + JERR        !IERR is sum of all errors
635       ALFA0 = ALFA1
636       Q00   = Q21
637  
638   50  CONTINUE
639  
640       CtLhA0TOA1 = ALFA0
641       NFF    = NF1
642  
643       RETURN
644       END
645  
646       FUNCTION CtLhALPHAR(QSQ,QS0,AS0,NF,IORD,IMODE,IERR)
647 C calculate ALPHAS FROM RGE GIVEN AS0 AT QS0.
648 C
649 C IORD=1: LEADING ORDER DEFINED BY 
650 C                 Q*d(alpha)/d(Q) = c1*alpha**2 
651 C
652 C IORD=2,IMODE=1: QCD NUM CHOICE DEFINED BY 
653 C                 Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
654 C
655 C IORD=2,IMODE=2: AD HOC ALTERNATIVE DEFINED BY 
656 C                 Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
657 C
658 C IORD=2,IMODE=3: TRADITIONAL CTEQ CHOICE DEFINED BY 
659 C                 ALPHA = c3*(1 - c4*log(L)/L)/L, WHERE L=log((Q/lambda)**2)
660 C
661 C c1 = -beta0/(2*pi)     where beta0 =  11. -  (2./3.)*nf
662 C c2 = -beta1/(8*pi**2)  where beta1 = 102. - (38./3.)*nf
663 C
664 C c3 = -2/c1
665 C c4 = -2*c2/c1**2
666 C
667 C
668       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
669  
670  
671       DATA PI / 3.14159265358979d0 /
672
673       BET0 = 11.d0 -(2*NF)/3.d0
674       BET1 = 102.d0 - (38*NF)/3.d0
675       B0   = BET0/(4.d0*PI)
676       B1   = BET1/(4.d0*PI*BET0)
677       IERR = 0
678
679       TERM0 = 1.d0/AS0+B0*LOG(QSQ/QS0)
680       IF(TERM0.LE.0.) THEN
681         CtLhALPHAR = 100.
682         IERR   = 1
683         PRINT *,'CtLhALPHAR WARNING: RETURN 100.'
684         RETURN
685       ENDIF
686       ALFA0 = 1.d0/TERM0
687
688 C ORDER=1 IS LEADING ORDER, WHICH IS SAME FOR ALL IMODE.
689       IF(IORD.EQ.1) THEN
690         CtLhALPHAR = ALFA0
691         RETURN
692       ELSEIF(IORD.NE.2) THEN
693         PRINT *,'FATAL ERROR: UNDEFINED ORDER IN CtLhALPHAR'
694         STOP
695       ENDIF
696
697
698 C QCD NUM CHOICE: Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
699       IF(IMODE .EQ. 1) THEN
700
701 c use Newton's method to solve the equation, instead of the 
702 c simple iterative method used in qcdnum (jcp 9/01)
703          DO ITER = 1, 20
704             ARG   = (1.d0/ALFA0+B1)/(1.d0/AS0+B1)
705             IF(ARG.LE.0.) THEN
706               CtLhALPHAR = 10.
707               IERR   = 1
708               PRINT *,'CtLhALPHAR WARNING: RETURN 10.'
709               RETURN
710             ENDIF 
711
712             TERM  = TERM0 + B1*LOG(ARG) - 1.d0/ALFA0
713             ALFA1 = ALFA0/(1.d0 + ALFA0*(1.d0 + B1*ALFA0)*TERM)
714
715             IF(ABS(ALFA1-ALFA0).LT.1.E-12) GOTO 20
716             ALFA0 = ALFA1
717          ENDDO
718
719          CtLhALPHAR = 10.
720          IERR   = 1
721          RETURN
722  
723 20       CONTINUE
724          CtLhALPHAR = ALFA1
725          RETURN
726
727 C AD HOC ALTERNATIVE: Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
728       ELSEIF(IMODE .EQ. 2) THEN
729
730 c first get a good starting point, to be sure Newton's method doesn't go 
731 c to the wrong root.
732          BEST = 9.d99
733          DO ITRY = 0, 20
734             IF(ITRY.NE.0) ALFA0 = (1./B1)*(ITRY-0.5D0)/20.D0
735             F = -1.d0/ALFA0 + TERM0 - B1*LOG(ALFA0/AS0) 
736             IF(ABS(F) .LT. BEST) THEN
737                BEST = ABS(F)
738                ALBST = ALFA0
739             ENDIF
740          ENDDO
741
742          ALFA0 = ALBST
743          DO ITER=1, 20
744             F = -1.d0/ALFA0 + TERM0 - B1*LOG(ALFA0/AS0) 
745             ALFA1 = ALFA0/(1.d0 + ALFA0*F/(1.d0 - B1*ALFA0))
746             IF(ABS(ALFA1-ALFA0) .LT. 1.E-12) GOTO 30
747             ALFA0 = ALFA1
748          ENDDO
749          
750          CtLhALPHAR = 10.
751          IERR   = 1
752          RETURN
753  
754   30     CONTINUE
755          CtLhALPHAR = ALFA1
756          RETURN
757
758 C TRADITIONAL CTEQ CHOICE: ALPHA = c3*(1 - c4*log(L)/L)/L, WHERE L=log((Q/lambda)**2)
759       ELSEIF(IMODE .EQ. 3) THEN
760
761          Z = -LOG(B0*AS0)
762          TMP = BET1/BET0**2
763
764          DO ITER = 1, 20
765             F = EXP(Z) - (1.D0 - TMP*Z*EXP(-Z))/(B0*AS0)
766             FPRI = EXP(Z)  + TMP*(1.D0-Z)*EXP(-Z)/(B0*AS0)
767             ZNEW = Z - F/FPRI
768             IF(ABS(Z-ZNEW) .LT. 1.E-10) GOTO 40
769             Z = ZNEW
770          ENDDO
771
772          CtLhALPHAR = 10.
773          IERR   = 1
774          RETURN
775
776    40    CONTINUE
777          XLAMSQ = QS0 * EXP(-EXP(ZNEW))
778          XL = LOG(QSQ/XLAMSQ)
779
780 c return a fixed value if no solution...
781          IF(XL .LE. 0.D0) THEN
782             CtLhALPHAR = 10.D0
783             IERR = 1
784             RETURN
785          ENDIF
786
787          CtLhALPHAR = (1.d0 - TMP*LOG(XL)/XL)/(B0*XL)
788
789 c place a cutoff if comes out very large...
790            if(CtLhALPHAR .gt. 10.d0) then
791               CtLhALPHAR = 10.d0
792               IERR = 1
793            endif
794
795          RETURN
796
797       ELSE
798          PRINT *,'FATAL UNDEFINED IMODE=',IMODE
799          STOP
800       ENDIF
801
802       RETURN
803       END
804
805       FUNCTION CtLhALPInew (AMU)
806 C Returns effective g**2/(4pi**2) = alpha/pi.
807       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
808       data pi / 3.14159265358979d0 /
809
810       q = amu
811       alpha = CtLhAlphaNew(q)
812       CtLhalpinew = alpha/pi
813
814       RETURN
815       END
816
817       subroutine CTEQ6NewAlpha(nset,mem)
818       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
819       common / QCDtable /  Alambda, Nfl, Iorder
820       dimension parms(2)
821
822       q0 = 91.188d0
823 c      call GetLam4M(nset,mem,alpha0)
824 c      call GetLam5M(nset,mem,ximode)
825       call listPDF(nset,mem,parms)
826       alpha0 = parms(1)
827       ximode = parms(2)
828       imode = int(ximode)
829       iiorder = iorder
830
831 c **********************************************************
832 c the input data file should probably specify a very large 
833 c value for xmt, because we never allow top quark as a parton
834 c in the PDF fitting, so there are never more than 5 active 
835 c flavors.  Presume it is most consistent therefore to 
836 c also keep top quark loop out of running of alpha.
837 c **********************************************************
838       
839       call GetQmass(4,cmass)
840       call GetQmass(5,bmass)
841       call GetQmass(6,tmass)
842 c      write(6,*) 'CTEQ6NewAlpha: mc, mb, mt=',
843 c     &  cmass, bmass, tmass
844       call CtLhAlphaNewSET(XMC,XMB,XMT,Q0,ALPHA0,IIORDER,IMODE)
845       RETURN
846       END
847
848 c=================================================
849       subroutine getnset(nset)
850       integer iset,nset
851       save iset
852       nset = iset      
853 c      print *,'getting nset:',nset
854       return
855 c      
856       entry setnset(nset)
857       iset = nset
858 c      print *,'setting nset:',nset
859       return      
860       end
861 c=================================================
862       subroutine getnmem(nset,nmem)
863       include 'parmsetup.inc'
864       integer nmem,nset,member(nmxset)
865       save member
866       nmem = member(nset)      
867       return
868 c      
869       entry setnmem(nset,nmem)
870       member(nset) = nmem
871       return      
872       end
873 c=================================================
874       subroutine alphacteq5f34(nflav,alphas,q)
875       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
876       parameter(pi=3.14159265358979323846d0)
877       
878       alphas = pi*CtLhALPI34(nflav,q)
879       return
880       end
881 c===============================================      
882       FUNCTION CtLhALPI34 (nflav,AMU)
883       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
884       COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
885       COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
886       LOGICAL SET
887       PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
888       DATA IW1, IW2 / 2*0 /
889       IF(.NOT.SET) CALL CtLhLAMCWZ
890       NEFF = LhCtNFL(AMU)
891 cccccc      
892       if(neff.gt.nflav) neff=nflav
893 cccccc
894       ALM  = ALAM(NEFF)
895       if(neff.eq.3) alm = 0.395
896       if(neff.eq.4) alm = 0.309
897       CtLhALPI34 = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT)
898       IF (IRT .EQ. 1) THEN
899          CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI34', 'AMU', AMU,
900      >              ALM, BIG, 1)
901       ELSEIF (IRT .EQ. 2) THEN
902          CALL CtLhWARNR(IW2,'CtLhALPI34 > 3; Be aware!','CtLhALPI34',
903      >  CtLhALPI34, D0, D1, 0)
904       ENDIF
905       RETURN
906       END