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