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