LHAPDF/lhapdf-5.9.1/src/wrapmstw-lite.f

   1 ! -*- F90 -*-
   2
   3 !-- LHAPDF version of MSTW interpolation code and alphaS routines.
   4 !-- 22/01/2009 by Graeme Watt <watt(at)hep.ucl.ac.uk>
   5
   6 !-- Modified to allow possibility of different heavy-quark masses.
   7 !-- 25/06/2010 by Graeme Watt <Graeme.Watt(at)cern.ch>
   8
   9 !-- Fix "NaN" bug for q <= m_c when m_c^2 < 1.25 GeV^2.
  10 !-- 25/01/2011 by Graeme Watt <Graeme.Watt(at)cern.ch>
  11
  12 subroutine MSTWevolve(x,Q,xpdf,xphoton)
  13   implicit none
  14   !
  15   include 'parmsetup.inc'
  16   character*16 name(nmxset)
  17   character*80 line
  18   character*512 setpath
  19   integer iset,mem,nset,nmem(nmxset),ndef(nmxset),mmem,imem
  20   COMMON/NAME/name,nmem,ndef,mmem
  21   double precision xpdf(-6:6),xvalence(6),xphoton, &
  22        &     x,Q,alfas,Qalfa,Eorder,Q2fit,MSTWALPHAS
  23   logical warn,fatal
  24   parameter(warn=.false.,fatal=.true.)
  25   !   Set warn=.true. to turn on warnings when extrapolating.
  26   !   Set fatal=.false. to return zero instead of terminating when
  27   !    invalid input values of x and q are used.
  28   integer ih,f,nhess,nx,nq,np,nqc0,nqb0,n,m,ip,io, &
  29        &     nExtraFlavours,i,j
  30   double precision xmin,xmax,qsqmin,qsqmax,mc2,mb2,eps, &
  31        &     qsq,xlog,qsqlog,res,res1,anom,ExtrapolateMSTWPDF, &
  32        &     InterpolateMSTWPDF
  33   parameter(nx=64,nq=48,np=12)
  34   parameter(xmin=1d-6,xmax=1d0,qsqmin=1d0,qsqmax=1d9,eps=1d-6)
  35   parameter(nhess=0)
  36   character dummyChar,dummyWord*50
  37   double precision gridx(nmxgridx),gridq(nmxgridq)
  38   integer ngridx,ngridq,jx,jq
  39   double precision ff(np,nx,nq)
  40   double precision qqorig(nq),qq(nq),xx(nx),cc(np,0:nhess,nx,nq,4,4,nmxset)
  41   double precision xxl(nx),qql(nq,nmxset,0:nhess)
  42   !   Store values of distance and tolerance along each eigenvector,
  43   !   heavy-quark masses and alphaS parameters in COMMON block.
  44   !   Allow the possibility of different alphaS parameters and
  45   !   heavy-quark masses for each member of the set.
  46   double precision distance(nmxset,0:nhess),tolerance(nmxset,0:nhess), &
  47        &     mCharm(nmxset,0:nhess),mBottom(nmxset,0:nhess), &
  48        &     alphaSQ0(nmxset,0:nhess),alphaSMZ(nmxset,0:nhess)
  49   integer alphaSorder(nmxset,0:nhess),alphaSnfmax(nmxset,0:nhess)
  50   common/mstwCommon/distance,tolerance, &
  51        &     mCharm,mBottom,alphaSQ0,alphaSMZ,alphaSorder,alphaSnfmax
  52   double precision mCharmSave,mBottomSave,mTopSave
  53   double precision cmass(nmxset),bmass(nmxset),tmass(nmxset)
  54   common/masses_LHA/cmass,bmass,tmass
  55   data xx/1d-6,2d-6,4d-6,6d-6,8d-6, &
  56        &     1d-5,2d-5,4d-5,6d-5,8d-5, &
  57        &     1d-4,2d-4,4d-4,6d-4,8d-4, &
  58        &     1d-3,2d-3,4d-3,6d-3,8d-3, &
  59        &     1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2, &
  60        &     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0, &
  61        &     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0, &
  62        &     .5d0,.525d0,.55d0,.575d0,.6d0,.625d0,.65d0,.675d0, &
  63        &     .7d0,.725d0,.75d0,.775d0,.8d0,.825d0,.85d0,.875d0, &
  64        &     .9d0,.925d0,.95d0,.975d0,1d0/
  65   data qqorig/1.d0, &
  66        &     1.25d0,1.5d0,0.d0,0.d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0, &
  67        &     1d1,1.2d1,0.d0,0.d0,2.6d1,4d1,6.4d1,1d2, &
  68        &     1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4, &
  69        &     1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6, &
  70        &     1.8d6,3.2d6,5.6d6,1d7,1.8d7,3.2d7,5.6d7,1d8, &
  71        &     1.8d8,3.2d8,5.6d8,1d9/
  72   save
  73
  74   call getnset(iset)
  75   call getnmem(iset,imem)
  76   ih = 0
  77
  78   qsq=q*q
  79   mc2=mCharm(iset,ih)**2
  80   mb2=mBottom(iset,ih)**2
  81   !   If mc2 < qsq < mc2+eps, then qsq = mc2+eps.
  82   if (qsq.gt.mc2.and.qsq.lt.mc2+eps) then
  83      qsq = mc2+eps
  84   end if
  85   !   If mb2 < qsq < mb2+eps, then qsq = mb2+eps.
  86   if (qsq.gt.mb2.and.qsq.lt.mb2+eps) then
  87      qsq = mb2+eps
  88   end if
  89
  90   xlog=log10(x)
  91   qsqlog=log10(qsq)
  92
  93   do f = 0, 13 ! loop over flavours
  94
  95      res = 0.d0
  96
  97      if (f.eq.0) then          ! gluon
  98         ip = 1
  99      else if (f.ge.1.and.f.le.5) then ! quarks
 100         ip = f+1
 101      else if (f.le.-1.and.f.ge.-5) then ! antiquarks
 102         ip = -f+1
 103      else if (f.ge.7.and.f.le.11) then ! valence quarks
 104         ip = f
 105      else if (f.eq.13) then    ! photon
 106         ip = 12
 107      else if (abs(f).ne.6.and.f.ne.12) then
 108         if (warn.or.fatal) print *,"Error in MSTWevolve: f = ",f
 109         if (fatal) stop
 110      end if
 111
 112      if (x.le.0.d0.or.x.gt.xmax.or.q.le.0.d0) then
 113
 114         if (warn.or.fatal) print *,"Error in MSTWevolve: x,qsq = ", &
 115              &        x,qsq
 116         if (fatal) stop
 117
 118      else if (abs(f).eq.6.or.f.eq.12) then ! set top quarks to zero
 119
 120         res = 0.d0
 121
 122      else if (qsq.lt.qsqmin) then ! extrapolate to low Q^2
 123
 124         if (warn) then
 125            print *, "Warning in MSTWevolve, extrapolating: f = ",f, &
 126                 &           ", x = ",x,", q = ",q
 127         end if
 128
 129         if (x.lt.xmin) then    ! extrapolate to low x
 130
 131            res = ExtrapolateMSTWPDF(ip,np,ih,nhess,xlog, &
 132                 & log10(qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 133            res1 = ExtrapolateMSTWPDF(ip,np,ih,nhess,xlog, &
 134                 & log10(1.01D0*qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 135            if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
 136               res = res - ExtrapolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 137                    & log10(qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 138               res1 = res1 - ExtrapolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 139                    & log10(1.01D0*qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 140            end if
 141
 142         else                   ! do usual interpolation
 143
 144            res = InterpolateMSTWPDF(ip,np,ih,nhess,xlog, &
 145                 & log10(qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 146            res1 = InterpolateMSTWPDF(ip,np,ih,nhess,xlog, &
 147                 & log10(1.01D0*qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 148            if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
 149               res = res - InterpolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 150                    & log10(qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 151               res1 = res1 - InterpolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 152                    & log10(1.01D0*qsqmin),nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 153            end if
 154
 155         end if
 156
 157         !   Calculate the anomalous dimension, dlog(xf)/dlog(qsq),
 158         !   evaluated at qsqmin.  Then extrapolate the PDFs to low
 159         !   qsq < qsqmin by interpolating the anomalous dimenion between
 160         !   the value at qsqmin and a value of 1 for qsq << qsqmin.
 161         !   If value of PDF at qsqmin is very small, just set
 162         !   anomalous dimension to 1 to prevent rounding errors.
 163         !   Impose minimum anomalous dimension of -2.5.
 164         if (abs(res).ge.1.D-5) then
 165            anom = max( -2.5D0, (res1-res)/res/0.01D0 )
 166         else
 167            anom = 1.D0
 168         end if
 169         res = res*(qsq/qsqmin)**(anom*qsq/qsqmin+1.D0-qsq/qsqmin)
 170
 171      else if (x.lt.xmin.or.qsq.gt.qsqmax) then ! extrapolate
 172
 173         if (warn) then
 174            print *, "Warning in MSTWevolve, extrapolating: f = ",f, &
 175                 &           ", x = ",x,", q = ",q
 176         end if
 177
 178         res = ExtrapolateMSTWPDF(ip,np,ih,nhess,xlog, &
 179              & qsqlog,nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 180
 181         if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
 182            res = res - ExtrapolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 183                 & qsqlog,nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 184         end if
 185
 186      else                      ! do usual interpolation
 187
 188         res = InterpolateMSTWPDF(ip,np,ih,nhess,xlog, &
 189              & qsqlog,nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 190
 191         if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
 192            res = res - InterpolateMSTWPDF(ip+5,np,ih,nhess,xlog, &
 193                 & qsqlog,nx,nq,xxl,qql(1,iset,ih),cc(1,0,1,1,1,1,iset))
 194         end if
 195
 196      end if
 197
 198      if (f.ge.7.and.f.le.12) then
 199         xvalence(f-6) = res
 200      else if (f.eq.13) then
 201         xphoton = res
 202      else
 203         xpdf(f) = res
 204      end if
 205
 206   end do
 207
 208   do f = 1, 6
 209      xpdf(-f) = xpdf(f) - xvalence(f) ! antiquarks
 210   end do
 211
 212   return
 213   !
 214   entry MSTWgetgrid(nset,ngridx,ngridq,gridx,gridq)
 215       do jx=1,nx
 216           gridx(jx)=xx(jx)
 217       enddo
 218       do jq=1,nq
 219           gridq(jq)=qq(jq)
 220       enddo
 221       ngridx=nx
 222       ngridq=nq
 223   return
 224
 225   entry MSTWread(nset)
 226
 227   ! Dummy read to get to the 'End:' (stream 1 is still open)
 228   read(1,*) nmem(nset),ndef(nset)
 229   do ih = 0, nmem(nset)
 230      !   Read header containing heavy-quark masses and alphaS values.
 231      read(1,*)
 232      read(1,*)
 233      read(1,*) dummyChar,dummyWord,dummyWord,dummyChar, &
 234           &        distance(nset,0),tolerance(nset,0)
 235      read(1,*) dummyChar,dummyWord,dummyChar,mCharm(nset,0)
 236      read(1,*) dummyChar,dummyWord,dummyChar,mBottom(nset,0)
 237      read(1,*) dummyChar,dummyWord,dummyChar,alphaSQ0(nset,0)
 238      read(1,*) dummyChar,dummyWord,dummyChar,alphaSMZ(nset,0)
 239      read(1,*) dummyChar,dummyWord,dummyWord,dummyChar, &
 240           &        alphaSorder(nset,0),alphaSnfmax(nset,0)
 241      read(1,*) dummyChar,dummyWord,dummyChar,nExtraFlavours
 242      read(1,*)
 243      read(1,*)
 244      !   Now read in the grids from the grid file.
 245      do n=1,nx-1
 246         do m=1,nq
 247            read(1,*,iostat=io)
 248            if (io.ne.0) then
 249               print *,"Error in MSTWevolve reading file"
 250               stop
 251            end if
 252         enddo
 253      enddo
 254   end do
 255
 256   return
 257   !
 258   entry MSTWalfa(alfas,Qalfa)
 259   call getnset(iset)
 260   call getnmem(iset,imem)
 261
 262   ! Check value of Qalfa is below appropriate thresholds.
 263   ! If not, redefine thresholds and reinitialise alphaS.
 264   IF (alphaSnfmax(iset,0).EQ.5) THEN ! maximum of five flavours
 265      IF (Qalfa.GT.mTopSave) THEN
 266         mTopSave = 2.D0*Qalfa
 267         CALL MSTWINITALPHAS(alphaSorder(iset,0), &
 268              &     1.D0,1.D0,alphaSQ0(iset,0), &
 269              &     mCharmSave,mBottomSave,mTopSave)
 270      END IF
 271   ELSE IF (alphaSnfmax(iset,0).EQ.4) THEN ! maximum of four flavours
 272      IF (Qalfa.GT.mBottomSave) THEN
 273         mBottomSave = 2.D0*Qalfa
 274         mTopSave = mBottomSave/0.9D0
 275         CALL MSTWINITALPHAS(alphaSorder(iset,0), &
 276              &     1.D0,1.D0,alphaSQ0(iset,0), &
 277              &     mCharmSave,mBottomSave,mTopSave)
 278      END IF
 279   ELSE IF (alphaSnfmax(iset,0).EQ.3) THEN ! maximum of three flavours
 280      IF (Qalfa.GT.mCharmSave) THEN
 281         mCharmSave = 2.D0*Qalfa
 282         mBottomSave = mCharmSave*0.9D0/0.8D0
 283         mTopSave = mBottomSave/0.9D0
 284         CALL MSTWINITALPHAS(alphaSorder(iset,0), &
 285              &     1.D0,1.D0,alphaSQ0(iset,0), &
 286              &     mCharmSave,mBottomSave,mTopSave)
 287      END IF
 288   END IF
 289   alfas = MSTWALPHAS(Qalfa)
 290   return
 291   !
 292   entry MSTWinit(nset,Eorder,Q2fit)
 293   return
 294   !
 295   entry MSTWpdf(mem)
 296   call getnset(iset)
 297   call setnmem(iset,mem)
 298
 299   ! here is the real read!!
 300
 301   ! have to reopen stream 1
 302   call getsetpath(setpath)
 303   open(1,file=setpath(1:len_trim(setpath)),action='READ')
 304
 305   line = ''
 306   do while (line(2:11).ne.'Evolution:')
 307      read(1,'(a)'),line
 308   enddo
 309   read(1,'(a)'),line
 310   read(1,'(a)'),line
 311
 312   read(1,*)nmem(iset),ndef(iset)
 313   ! Dummy read up to the member requested
 314   do i=0,mem-1
 315      do j=1,11
 316         read(1,*)
 317      enddo
 318      do n=1,nx-1
 319         do m=1,nq
 320            read(1,'(a)')line
 321         enddo
 322      enddo
 323   enddo
 324
 325   ih=0
 326   !   Read header containing heavy-quark masses and alphaS values.
 327   read(1,*)
 328   read(1,*)
 329   read(1,*) dummyChar,dummyWord,dummyWord,dummyChar, &
 330        &        distance(iset,ih),tolerance(iset,ih)
 331   read(1,*) dummyChar,dummyWord,dummyChar,mCharm(iset,ih)
 332   read(1,*) dummyChar,dummyWord,dummyChar,mBottom(iset,ih)
 333   read(1,*) dummyChar,dummyWord,dummyChar,alphaSQ0(iset,ih)
 334   read(1,*) dummyChar,dummyWord,dummyChar,alphaSMZ(iset,ih)
 335   read(1,*) dummyChar,dummyWord,dummyWord,dummyChar, &
 336        &        alphaSorder(iset,ih),alphaSnfmax(iset,ih)
 337   read(1,*) dummyChar,dummyWord,dummyChar,nExtraFlavours
 338   read(1,*)
 339   read(1,*)
 340
 341   !   Heavy-quark masses and alphaS values
 342   !   are stored in a COMMON block.
 343   !  WRITE(6,*) "mCharm = ",mCharm(iset,ih), &
 344   !       &     ", mBottom = ",mBottom(iset,ih)
 345   !  WRITE(6,*) "alphaS(Q0) = ",alphaSQ0(iset,ih),", alphaS(MZ) = ", &
 346   !       &     alphaSMZ(iset,ih),", alphaSorder = ",alphaSorder(iset,ih), &
 347   !       &     ", alphaSnfmax = ",alphaSnfmax(iset,ih)
 348
 349   mc2=mCharm(iset,ih)**2
 350   mb2=mBottom(iset,ih)**2
 351
 352   !   Check that the heavy-quark masses are sensible.
 353   !   Redistribute grid points if not in usual range.
 354   do m=1,nq
 355      qq(m) = qqorig(m)
 356   end do
 357   if (mc2.le.qq(1).or.mc2+eps.ge.qq(8)) then
 358      print *,"Error in MSTWevolve: invalid mCharm = ",mCharm(iset,ih)
 359      stop
 360   else if (mc2.lt.qq(2)) then
 361      nqc0=2
 362      qq(4)=qq(2)
 363      qq(5)=qq(3)
 364   else if (mc2.lt.qq(3)) then
 365      nqc0=3
 366      qq(5)=qq(3)
 367   else if (mc2.lt.qq(6)) then
 368      nqc0=4
 369   else if (mc2.lt.qq(7)) then
 370      nqc0=5
 371      qq(4)=qq(6)
 372   else
 373      nqc0=6
 374      qq(4)=qq(6)
 375      qq(5)=qq(7)
 376   end if
 377   if (mb2.le.qq(12).or.mb2+eps.ge.qq(17)) then
 378      print *,"Error in MSTWevolve: invalid mBottom = ",mBottom(iset,ih)
 379      stop
 380   else if (mb2.lt.qq(13)) then
 381      nqb0=13
 382      qq(15)=qq(13)
 383   else if (mb2.lt.qq(16)) then
 384      nqb0=14
 385   else
 386      nqb0=15
 387      qq(14)=qq(16)
 388   end if
 389   qq(nqc0)=mc2
 390   qq(nqc0+1)=mc2+eps
 391   qq(nqb0)=mb2
 392   qq(nqb0+1)=mb2+eps
 393
 394   !   The nExtraFlavours variable is provided to aid compatibility
 395   !   with future grids where, for example, a photon distribution
 396   !   might be provided (cf. the MRST2004QED PDFs).
 397   if (nExtraFlavours.lt.0.or.nExtraFlavours.gt.1) then
 398      print *,"Error in MSTWevolve: invalid nExtraFlavours = ", &
 399           &           nExtraFlavours
 400      stop
 401   end if
 402
 403   !   Now read in the grids from the grid file.
 404   do n=1,nx-1
 405      do m=1,nq
 406         if (nExtraFlavours.gt.0) then
 407            if (alphaSorder(iset,ih).eq.2) then ! NNLO
 408               read(1,'(12(1pe12.4))',iostat=io) &
 409                    & (ff(ip,n,m),ip=1,12)
 410            else          ! LO or NLO
 411               ff(10,n,m) = 0.d0 ! = chm-cbar
 412               ff(11,n,m) = 0.d0 ! = bot-bbar
 413               read(1,'(10(1pe12.4))',iostat=io) &
 414                    & (ff(ip,n,m),ip=1,9),ff(12,n,m)
 415            end if
 416         else             ! nExtraFlavours = 0
 417            if (alphaSorder(iset,ih).eq.2) then ! NNLO
 418               ff(12,n,m) = 0.d0 ! = photon
 419               read(1,'(11(1pe12.4))',iostat=io) &
 420                    & (ff(ip,n,m),ip=1,11)
 421            else          ! LO or NLO
 422               ff(10,n,m) = 0.d0 ! = chm-cbar
 423               ff(11,n,m) = 0.d0 ! = bot-bbar
 424               ff(12,n,m) = 0.d0 ! = photon
 425               read(1,'(9(1pe12.4))',iostat=io) &
 426                    & (ff(ip,n,m),ip=1,9)
 427            end if
 428         end if
 429         if (io.ne.0) then
 430            print *,"Error in MSTWevolve reading file"
 431            stop
 432         end if
 433      enddo
 434   enddo
 435
 436   !   PDFs are identically zero at x = 1.
 437   do m=1,nq
 438      do ip=1,np
 439         ff(ip,nx,m)=0d0
 440      enddo
 441   enddo
 442
 443   do n=1,nx
 444      xxl(n)=log10(xx(n))
 445   enddo
 446   do m=1,nq
 447      qql(m,iset,ih)=log10(qq(m))
 448   enddo
 449
 450   !   Initialise all parton flavours.
 451   do ip=1,np
 452      call InitialiseMSTWPDF(ip,np,ih,nhess,nx,nq,nqc0,nqb0, &
 453           & xxl,qql(1,iset,ih),ff,cc(1,0,1,1,1,1,iset))
 454   enddo
 455   ! end of section moved from the read routine
 456
 457   close(1)
 458
 459   ! Call the initialisation routine with alpha_S(Q_0).
 460   IF (alphaSnfmax(iset,ih).EQ.5) THEN ! maximum of five flavours
 461      mCharmSave = mCharm(iset,ih)
 462      mBottomSave = mBottom(iset,ih)
 463      mTopSave = 1.D10
 464   ELSE IF (alphaSnfmax(iset,ih).EQ.4) THEN ! maximum of four flavours
 465      mCharmSave = mCharm(iset,ih)
 466      mBottomSave = 0.9D10
 467      mTopSave = 1.D10
 468   ELSE IF (alphaSnfmax(iset,ih).EQ.3) THEN ! maximum of three flavours
 469      mCharmSave = 0.8D10
 470      mBottomSave = 0.9D10
 471      mTopSave = 1.D10
 472   END IF
 473
 474   ! Update masses stored in common/masses_LHA/.
 475   cmass(iset) = mCharm(iset,ih)
 476   bmass(iset) = mBottom(iset,ih)
 477   tmass(iset) = mTopSave
 478
 479   CALL MSTWINITALPHAS(alphaSorder(iset,ih),1.D0,1.D0,alphaSQ0(iset,ih), &
 480        &     mCharmSave,mBottomSave,mTopSave)
 481
 482   !  !   Check calculated value of alpha_S(M_Z) matches stored value.
 483   !  WRITE(6,'(" alphaS(MZ) = ",F7.5," = ",F7.5)') &
 484   !       &     MSTWALPHAS(91.1876D0),alphaSMZ(iset,ih)
 485
 486   return
 487   !
 488 END subroutine MSTWevolve
 489
 490
 491
 492 !----------------------------------------------------------------------
 493
 494 subroutine InitialiseMSTWPDF(ip,np,ih,nhess,nx,my,myc0,myb0, &
 495      &     xx,yy,ff,cc)
 496   implicit none
 497   integer nhess,ih,nx,my,myc0,myb0,j,k,l,m,n,ip,np
 498   double precision xx(nx),yy(my),ff(np,nx,my), &
 499        &     ff1(nx,my),ff2(nx,my),ff12(nx,my),ff21(nx,my), &
 500        &     yy0(4),yy1(4),yy2(4),yy12(4),z(16), &
 501        &     cl(16),cc(np,0:nhess,nx,my,4,4),iwt(16,16), &
 502        &     polderiv1,polderiv2,polderiv3,d1,d2,d1d2,xxd
 503
 504       data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, &
 505      &     0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0, &
 506      &     -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0, &
 507      &     2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0, &
 508      &     0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, &
 509      &     0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0, &
 510      &     0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1, &
 511      &     0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1, &
 512      &     -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0, &
 513      &     0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0, &
 514      &     9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2, &
 515      &     -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2, &
 516      &     2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0, &
 517      &     0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0, &
 518      &     -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1, &
 519      &     4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
 520
 521       do m=1,my
 522          ff1(1,m)=polderiv1(xx(1),xx(2),xx(3), &
 523      &        ff(ip,1,m),ff(ip,2,m),ff(ip,3,m))
 524          ff1(nx,m)=polderiv3(xx(nx-2),xx(nx-1),xx(nx), &
 525      &        ff(ip,nx-2,m),ff(ip,nx-1,m),ff(ip,nx,m))
 526          do n=2,nx-1
 527             ff1(n,m)=polderiv2(xx(n-1),xx(n),xx(n+1), &
 528      &           ff(ip,n-1,m),ff(ip,n,m),ff(ip,n+1,m))
 529          enddo
 530       enddo
 531
 532 !--   Calculate the derivatives at qsq=mc2,mc2+eps,mb2,mb2+eps
 533 !--   in a similar way as at the endpoints qsqmin and qsqmax.
 534       do n=1,nx
 535          do m=1,my
 536             if (myc0.eq.2.and.m.eq.1) then
 537                ff2(n,m)=(ff(ip,n,m+1)-ff(ip,n,m))/(yy(m+1)-yy(m))
 538             else if (myc0.eq.2.and.m.eq.2) then
 539                ff2(n,m)=(ff(ip,n,m)-ff(ip,n,m-1))/(yy(m)-yy(m-1))
 540             else if (m.eq.1.or.m.eq.myc0+1.or.m.eq.myb0+1) then
 541                ff2(n,m)=polderiv1(yy(m),yy(m+1),yy(m+2), &
 542      &              ff(ip,n,m),ff(ip,n,m+1),ff(ip,n,m+2))
 543             else if (m.eq.my.or.m.eq.myc0.or.m.eq.myb0) then
 544                ff2(n,m)=polderiv3(yy(m-2),yy(m-1),yy(m), &
 545      &              ff(ip,n,m-2),ff(ip,n,m-1),ff(ip,n,m))
 546             else
 547                ff2(n,m)=polderiv2(yy(m-1),yy(m),yy(m+1), &
 548      &              ff(ip,n,m-1),ff(ip,n,m),ff(ip,n,m+1))
 549             end if
 550          end do
 551       end do
 552
 553 !--   Calculate the cross derivatives (d/dx)(d/dy).
 554       do m=1,my
 555          ff12(1,m)=polderiv1(xx(1),xx(2),xx(3), &
 556      &        ff2(1,m),ff2(2,m),ff2(3,m))
 557          ff12(nx,m)=polderiv3(xx(nx-2),xx(nx-1),xx(nx), &
 558      &        ff2(nx-2,m),ff2(nx-1,m),ff2(nx,m))
 559          do n=2,nx-1
 560             ff12(n,m)=polderiv2(xx(n-1),xx(n),xx(n+1), &
 561      &           ff2(n-1,m),ff2(n,m),ff2(n+1,m))
 562          enddo
 563       enddo
 564
 565 !--   Calculate the cross derivatives (d/dy)(d/dx).
 566       do n=1,nx
 567          do m = 1, my
 568             if (myc0.eq.2.and.m.eq.1) then
 569                ff21(n,m)=(ff1(n,m+1)-ff1(n,m))/(yy(m+1)-yy(m))
 570             else if (myc0.eq.2.and.m.eq.2) then
 571                ff21(n,m)=(ff1(n,m)-ff1(n,m-1))/(yy(m)-yy(m-1))
 572             else if (m.eq.1.or.m.eq.myc0+1.or.m.eq.myb0+1) then
 573                ff21(n,m)=polderiv1(yy(m),yy(m+1),yy(m+2), &
 574      &              ff1(n,m),ff1(n,m+1),ff1(n,m+2))
 575             else if (m.eq.my.or.m.eq.myc0.or.m.eq.myb0) then
 576                ff21(n,m)=polderiv3(yy(m-2),yy(m-1),yy(m), &
 577      &              ff1(n,m-2),ff1(n,m-1),ff1(n,m))
 578             else
 579                ff21(n,m)=polderiv2(yy(m-1),yy(m),yy(m+1), &
 580      &              ff1(n,m-1),ff1(n,m),ff1(n,m+1))
 581             end if
 582          end do
 583       end do
 584
 585 !--   Take the average of (d/dx)(d/dy) and (d/dy)(d/dx).
 586       do n=1,nx
 587          do m = 1, my
 588             ff12(n,m)=0.5*(ff12(n,m)+ff21(n,m))
 589          end do
 590       end do
 591
 592       do n=1,nx-1
 593          do m=1,my-1
 594             d1=xx(n+1)-xx(n)
 595             d2=yy(m+1)-yy(m)
 596             d1d2=d1*d2
 597
 598             yy0(1)=ff(ip,n,m)
 599             yy0(2)=ff(ip,n+1,m)
 600             yy0(3)=ff(ip,n+1,m+1)
 601             yy0(4)=ff(ip,n,m+1)
 602
 603             yy1(1)=ff1(n,m)
 604             yy1(2)=ff1(n+1,m)
 605             yy1(3)=ff1(n+1,m+1)
 606             yy1(4)=ff1(n,m+1)
 607
 608             yy2(1)=ff2(n,m)
 609             yy2(2)=ff2(n+1,m)
 610             yy2(3)=ff2(n+1,m+1)
 611             yy2(4)=ff2(n,m+1)
 612
 613             yy12(1)=ff12(n,m)
 614             yy12(2)=ff12(n+1,m)
 615             yy12(3)=ff12(n+1,m+1)
 616             yy12(4)=ff12(n,m+1)
 617
 618             do k=1,4
 619                z(k)=yy0(k)
 620                z(k+4)=yy1(k)*d1
 621                z(k+8)=yy2(k)*d2
 622                z(k+12)=yy12(k)*d1d2
 623             enddo
 624
 625             do l=1,16
 626                xxd=0.d0
 627                do k=1,16
 628                   xxd=xxd+iwt(k,l)*z(k)
 629                enddo
 630                cl(l)=xxd
 631             enddo
 632             l=0
 633             do k=1,4
 634                do j=1,4
 635                   l=l+1
 636                   cc(ip,ih,n,m,k,j)=cl(l)
 637                enddo
 638             enddo
 639          enddo
 640       enddo
 641       return
 642       end
 643
 644 !----------------------------------------------------------------------
 645
 646       double precision function InterpolateMSTWPDF(ip,np,ih,nhess,x,y, &
 647      &     nx,my,xx,yy,cc)
 648       implicit none
 649       integer ih,nx,my,nhess,MSTWlocx,l,m,n,ip,np
 650       double precision xx(nx),yy(my),cc(np,0:nhess,nx,my,4,4), &
 651      &     x,y,z,t,u
 652
 653       n=MSTWlocx(xx,nx,x)
 654       m=MSTWlocx(yy,my,y)
 655
 656       t=(x-xx(n))/(xx(n+1)-xx(n))
 657       u=(y-yy(m))/(yy(m+1)-yy(m))
 658
 659       z=0.d0
 660       do l=4,1,-1
 661          z=t*z+((cc(ip,ih,n,m,l,4)*u+cc(ip,ih,n,m,l,3))*u &
 662      &        +cc(ip,ih,n,m,l,2))*u+cc(ip,ih,n,m,l,1)
 663       enddo
 664
 665       InterpolateMSTWPDF = z
 666
 667       return
 668       end
 669
 670 !----------------------------------------------------------------------
 671
 672       double precision function ExtrapolateMSTWPDF(ip,np,ih,nhess,x,y, &
 673      &     nx,my,xx,yy,cc)
 674       implicit none
 675       integer ih,nx,my,nhess,MSTWlocx,n,m,ip,np
 676       double precision xx(nx),yy(my),cc(np,0:nhess,nx,my,4,4), &
 677      &     x,y,z,f0,f1,z0,z1,InterpolateMSTWPDF
 678
 679       n=MSTWlocx(xx,nx,x)           ! 0: below xmin, nx: above xmax
 680       m=MSTWlocx(yy,my,y)           ! 0: below qsqmin, my: above qsqmax
 681
 682 !   If extrapolation in small x only:
 683       if (n.eq.0.and.m.gt.0.and.m.lt.my) then
 684          f0 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(1),y,nx,my,xx,yy,cc)
 685          f1 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(2),y,nx,my,xx,yy,cc)
 686          if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
 687             z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
 688          else
 689             z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
 690          end if
 691 !   If extrapolation into large q only:
 692       else if (n.gt.0.and.m.eq.my) then
 693          f0 = InterpolateMSTWPDF(ip,np,ih,nhess,x,yy(my),nx,my,xx,yy,cc)
 694          f1 = InterpolateMSTWPDF(ip,np,ih,nhess,x,yy(my-1),nx,my,xx,yy,cc)
 695          if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
 696             z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))* &
 697      &           (y-yy(my)))
 698          else
 699             z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 700          end if
 701 !   If extrapolation into large q AND small x:
 702       else if (n.eq.0.and.m.eq.my) then
 703          f0 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(1),yy(my),nx,my,xx,yy,cc)
 704          f1 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(1),yy(my-1),nx,my,xx,yy, &
 705      &        cc)
 706          if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
 707             z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))* &
 708      &           (y-yy(my)))
 709          else
 710             z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 711          end if
 712          f0 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(2),yy(my),nx,my,xx,yy,cc)
 713          f1 = InterpolateMSTWPDF(ip,np,ih,nhess,xx(2),yy(my-1),nx,my,xx,yy, &
 714      &        cc)
 715          if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
 716             z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))* &
 717      &           (y-yy(my)))
 718          else
 719             z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 720          end if
 721          if (z0.gt.1.d-3.and.z1.gt.1.d-3) then
 722             z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
 723          else
 724             z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
 725          end if
 726       else
 727          print *,"Error in ExtrapolateMSTWPDF"
 728          stop
 729       end if
 730
 731       ExtrapolateMSTWPDF = z
 732
 733       return
 734       end
 735
 736 !----------------------------------------------------------------------
 737
 738       integer function MSTWlocx(xx,nx,x)
 739 !   returns an integer j such that x lies inbetween xx(j) and xx(j+1).
 740 !   nx is the length of the array with xx(nx) the highest element.
 741       implicit none
 742       integer nx,jl,ju,jm
 743       double precision x,xx(nx)
 744       if(x.eq.xx(1)) then
 745          MSTWlocx=1
 746          return
 747       endif
 748       if(x.eq.xx(nx)) then
 749          MSTWlocx=nx-1
 750          return
 751       endif
 752       ju=nx+1
 753       jl=0
 754     1 if((ju-jl).le.1) goto 2
 755       jm=(ju+jl)/2
 756       if(x.ge.xx(jm)) then
 757          jl=jm
 758       else
 759          ju=jm
 760       endif
 761       goto 1
 762     2 MSTWlocx=jl
 763       return
 764       end
 765
 766 !----------------------------------------------------------------------
 767
 768       double precision function polderiv1(x1,x2,x3,y1,y2,y3)
 769 !--   returns the estimate of the derivative at x1 obtained by a
 770 !--   polynomial interpolation using the three points (x_i,y_i).
 771       implicit none
 772       double precision x1,x2,x3,y1,y2,y3
 773       polderiv1=(x3*x3*(y1-y2)+2.d0*x1*(x3*(-y1+y2)+x2*(y1-y3)) &
 774      &     +x2*x2*(-y1+y3)+x1*x1*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3))
 775       return
 776       end
 777
 778       double precision function polderiv2(x1,x2,x3,y1,y2,y3)
 779 !--   returns the estimate of the derivative at x2 obtained by a
 780 !--   polynomial interpolation using the three points (x_i,y_i).
 781       implicit none
 782       double precision x1,x2,x3,y1,y2,y3
 783       polderiv2=(x3*x3*(y1-y2)-2.d0*x2*(x3*(y1-y2)+x1*(y2-y3)) &
 784      &     +x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
 785       return
 786       end
 787
 788       double precision function polderiv3(x1,x2,x3,y1,y2,y3)
 789 !--   returns the estimate of the derivative at x3 obtained by a
 790 !--   polynomial interpolation using the three points (x_i,y_i).
 791       implicit none
 792       double precision x1,x2,x3,y1,y2,y3
 793       polderiv3=(x3*x3*(-y1+y2)+2.d0*x2*x3*(y1-y3)+x1*x1*(y2-y3) &
 794      &     +x2*x2*(-y1+y3)+2.d0*x1*x3*(-y2+y3))/ &
 795      &     ((x1-x2)*(x1-x3)*(x2-x3))
 796       return
 797       end
 798
 799 !----------------------------------------------------------------------
 800
 801 !----------------------------------------------------------------------
 802 !--   Stand-alone code for alpha_s cannibalised (with permission)
 803 !--   from Andreas Vogt's QCD-PEGASUS package (hep-ph/0408244).
 804 !--   The running coupling alpha_s is obtained at N^mLO (m = 0,1,2,3)
 805 !--   by solving the renormalisation group equation in the MSbar scheme
 806 !--   by a fourth-order Runge-Kutta integration.  Transitions from
 807 !--   n_f to n_f+1 flavours are made when the factorisation scale
 808 !--   mu_f equals the pole masses m_h (h = c,b,t).  At exactly
 809 !--   the thresholds m_{c,b,t}, the number of flavours n_f = {3,4,5}.
 810 !--   The top quark mass should be set to be very large to evolve with
 811 !--   a maximum of five flavours.  The factorisation scale mu_f may be
 812 !--   a constant multiple of the renormalisation scale mu_r.  The input
 813 !--   factorisation scale mu_(f,0) should be less than or equal to
 814 !--   the charm quark mass.  However, if it is greater than the
 815 !--   charm quark mass, the value of alpha_s at mu_(f,0) = 1 GeV will
 816 !--   be found using a root-finding algorithm.
 817 !--
 818 !--   Example of usage.
 819 !--   First call the initialisation routine (only needed once):
 820 !--
 821 !--    IORD = 2                  ! perturbative order (N^mLO,m=0,1,2,3)
 822 !--    FR2 = 1.D0                ! ratio of mu_f^2 to mu_r^2
 823 !--    MUR = 1.D0                ! input mu_r in GeV
 824 !--    ASMUR = 0.5D0             ! input value of alpha_s at mu_r
 825 !--    MC = 1.4D0                ! charm quark mass
 826 !--    MB = 4.75D0               ! bottom quark mass
 827 !--    MT = 1.D10                ! top quark mass
 828 !--    CALL MSTWINITALPHAS(IORD, FR2, MUR, ASMUR, MC, MB, MT)
 829 !--
 830 !--   Then get alpha_s at a renormalisation scale mu_r with:
 831 !--
 832 !--    MUR = 100.D0              ! renormalisation scale in GeV
 833 !--    ALFAS = MSTWALPHAS(MUR)
 834 !--
 835 !----------------------------------------------------------------------
 836 !--   Comments to Graeme Watt <watt(at)hep.ucl.ac.uk>
 837 !----------------------------------------------------------------------
 838
 839       SUBROUTINE MSTWINITALPHAS(IORD, FR2, MUR, ASMUR, MC, MB, MT)
 840 !--   IORD = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
 841 !--   FR2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
 842 !--   MUR = input renormalisation scale (in GeV) for alpha_s.
 843 !--   ASMUR = input value of alpha_s at the renormalisation scale MUR.
 844 !--   MC,MB,MT = heavy-quark masses in GeV.
 845       IMPLICIT NONE
 846       INTEGER IORD,IORDc,MAXF,MODE
 847       DOUBLE PRECISION FR2,MUR,ASMUR,MC,MB,MT,EPS,A,B,MSTWDZEROX, &
 848      &     R0c,FR2c,MURc,ASMURc,MCc,MBc,MTc,FINDALPHASR0,R0,ASI
 849       COMMON / MSTWDZEROXcommon / FR2c,MURc,ASMURc,MCc,MBc,MTc,R0c,IORDc
 850       PARAMETER(EPS=1.D-10,MAXF=10000,MODE=1)
 851       EXTERNAL FINDALPHASR0
 852
 853       IF (MUR*sqrt(FR2).LE.MC) THEN ! Check that MUF <= MC.
 854          R0 = MUR
 855          ASI = ASMUR
 856       ELSE                      ! Solve for alpha_s at R0 = 1 GeV.
 857 !--   Copy variables to common block.
 858          R0c = 1.D0/sqrt(FR2)
 859          IORDc = IORD
 860          FR2c = FR2
 861          MURc = MUR
 862          ASMURc = ASMUR
 863          MCc = MC
 864          MBc = MB
 865          MTc = MT
 866 !--   Now get alpha_s(R0) corresponding to alpha_s(MUR).
 867          A = 0.02D0              ! lower bound for alpha_s(R0)
 868          B = 2.00D0              ! upper bound for alpha_s(R0)
 869          R0 = R0c
 870          ASI = MSTWDZEROX(A,B,EPS,MAXF,FINDALPHASR0,MODE)
 871       END IF
 872
 873       CALL MSTWINITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
 874
 875       RETURN
 876       END
 877
 878 !----------------------------------------------------------------------
 879
 880 !--   Find the zero of this function using MSTWDZEROX.
 881       DOUBLE PRECISION FUNCTION FINDALPHASR0(ASI)
 882       IMPLICIT NONE
 883       INTEGER IORD
 884       DOUBLE PRECISION FR2, R0, ASI, MC, MB, MT, MUR, ASMUR, MSTWALPHAS
 885       COMMON / MSTWDZEROXcommon / FR2, MUR, ASMUR, MC, MB, MT, R0, IORD
 886
 887       CALL MSTWINITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
 888       FINDALPHASR0 = MSTWALPHAS(MUR) - ASMUR ! solve equal to zero
 889
 890       RETURN
 891       END
 892
 893 !----------------------------------------------------------------------
 894
 895       SUBROUTINE MSTWINITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
 896 !--   IORD = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
 897 !--   FR2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
 898 !--   R0 = input renormalisation scale (in GeV) for alphas_s.
 899 !--   ASI = input value of alpha_s at the renormalisation scale R0.
 900 !--   MC,MB,MT = heavy-quark masses in GeV.
 901 !--   Must have R0*sqrt(FR2) <= MC to call this subroutine.
 902       IMPLICIT NONE
 903       INTEGER IORD,NAORD,NASTPS,IVFNS,NFF
 904       DOUBLE PRECISION FR2,R0,ASI,MC,MB,MT,LOGFR,R20, &
 905      &     PI,ZETA,CF,CA,TR,AS0,M20,MC2,MB2,MT2
 906       PARAMETER(PI = 3.14159265358979D0)
 907
 908       COMMON / RZETA  / ZETA(6)
 909       COMMON / COLOUR / CF, CA, TR
 910       COMMON / ASINP  / AS0, M20
 911       COMMON / ASPAR  / NAORD, NASTPS
 912       COMMON / VARFLV / IVFNS
 913       COMMON / NFFIX  / NFF
 914       COMMON / FRRAT  / LOGFR
 915
 916 !
 917 ! ..QCD colour factors
 918 !
 919       CA = 3.D0
 920       CF = 4./3.D0
 921       TR = 0.5D0
 922 !
 923 ! ..The lowest integer values of the Zeta function
 924 !
 925       ZETA(1) = 0.57721566490153D0
 926       ZETA(2) = 1.644934066848226D0
 927       ZETA(3) = 1.202056903159594D0
 928       ZETA(4) = 1.082323233711138D0
 929       ZETA(5) = 1.036927755143370D0
 930       ZETA(6) = 1.017343061984449D0
 931
 932       IVFNS = 1                 ! variable flavour-number scheme (VFNS)
 933 !      IVFNS = 0                 ! fixed flavour-number scheme (FFNS)
 934       NFF = 4                   ! number of flavours for FFNS
 935       NAORD = IORD              ! perturbative order of alpha_s
 936       NASTPS = 20               ! num. steps in Runge-Kutta integration
 937       R20 = R0**2               ! input renormalisation scale
 938       MC2 = MC**2               ! mu_f^2 for charm threshold
 939       MB2 = MB**2               ! mu_f^2 for bottom threshold
 940       MT2 = MT**2               ! mu_f^2 for top threshold
 941       LOGFR = LOG(FR2)          ! log of ratio of mu_f^2 to mu_r^2
 942       M20 = R20 * FR2           ! input factorisation scale
 943
 944 !
 945 ! ..Stop some nonsense
 946 !
 947       IF ( (IVFNS .EQ. 0) .AND. (NFF .LT. 3) ) THEN
 948          WRITE (6,*) 'Wrong flavour number for FFNS evolution. STOP'
 949          STOP
 950       END IF
 951       IF ( (IVFNS .EQ. 0) .AND. (NFF .GT. 5) ) THEN
 952          WRITE (6,*) 'Wrong flavour number for FFNS evolution. STOP'
 953          STOP
 954       END IF
 955 !
 956       IF ( NAORD .GT. 3 ) THEN
 957          WRITE (6,*) 'Specified order in a_s too high. STOP'
 958          STOP
 959       END IF
 960 !
 961       IF ( (IVFNS .NE. 0) .AND. (FR2 .GT. 4.001D0) ) THEN
 962          WRITE (6,*) 'Too low mu_r for VFNS evolution. STOP'
 963          STOP
 964       END IF
 965 !
 966       IF ( (IVFNS .EQ. 1) .AND. (M20 .GT. MC2) ) THEN
 967          WRITE (6,*) 'Too high mu_0 for VFNS evolution. STOP'
 968          STOP
 969       END IF
 970 !
 971       IF ( (ASI .GT. 2.D0) .OR. (ASI .LT. 2.D-2) ) THEN
 972          WRITE (6,*) 'alpha_s out of range. STOP'
 973          STOP
 974       END IF
 975 !
 976       IF ( (IVFNS .EQ. 1) .AND. (MC2 .GT. MB2) ) THEN
 977          WRITE (6,*) 'Wrong charm-bottom mass hierarchy. STOP'
 978          STOP
 979       END IF
 980       IF ( (IVFNS .EQ. 1) .AND. (MB2 .GT. MT2) ) THEN
 981          WRITE (6,*) 'Wrong bottom-top mass hierarchy. STOP'
 982          STOP
 983       END IF
 984 !
 985
 986 !--   Store the beta function coefficients in a COMMON block.
 987       CALL BETAFCT
 988
 989 !--   Store a_s = alpha_s(mu_r^2)/(4 pi) at the input scale R0.
 990       AS0 = ASI / (4.D0* PI)
 991
 992 !--   Store alpha_s at the heavy flavour thresholds in a COMMON block.
 993        IF (IVFNS .NE. 0) THEN
 994           CALL EVNFTHR (MC2, MB2, MT2)
 995        END IF
 996
 997       RETURN
 998       END
 999
1000 !----------------------------------------------------------------------
1001
1002       DOUBLE PRECISION FUNCTION MSTWALPHAS(MUR)
1003       IMPLICIT NONE
1004       INTEGER NFF,IVFNS,NF
1005       DOUBLE PRECISION PI,LOGFR,AS0,M20,ASC,M2C,ASB,M2B,AST,M2T,M2,MUR, &
1006      &     R2,ASI,ASF,R20,R2T,R2B,R2C,AS
1007       PARAMETER ( PI = 3.14159265358979D0 )
1008 !
1009 ! ..Input common blocks
1010 !
1011        COMMON / NFFIX  / NFF
1012        COMMON / VARFLV / IVFNS
1013        COMMON / FRRAT  / LOGFR
1014        COMMON / ASINP  / AS0, M20
1015        COMMON / ASFTHR / ASC, M2C, ASB, M2B, AST, M2T
1016
1017        R2 = MUR**2
1018        M2 = R2 * EXP(+LOGFR)
1019        IF (IVFNS .EQ. 0) THEN
1020 !
1021 !   Fixed number of flavours
1022 !
1023           NF  = NFF
1024           R20 = M20 * R2/M2
1025           ASI = AS0
1026           ASF = AS (R2, R20, AS0, NF)
1027 !
1028        ELSE
1029 !
1030 ! ..Variable number of flavours
1031 !
1032           IF (M2 .GT. M2T) THEN
1033              NF = 6
1034              R2T = M2T * R2/M2
1035              ASI = AST
1036              ASF = AS (R2, R2T, AST, NF)
1037 !
1038           ELSE IF (M2 .GT. M2B) THEN
1039              NF = 5
1040              R2B = M2B * R2/M2
1041              ASI = ASB
1042              ASF = AS (R2, R2B, ASB, NF)
1043 !
1044           ELSE IF (M2 .GT. M2C) THEN
1045              NF = 4
1046              R2C = M2C * R2/M2
1047              ASI = ASC
1048              ASF = AS (R2, R2C, ASC, NF)
1049 !
1050           ELSE
1051              NF = 3
1052              R20 = M20 * R2/M2
1053              ASI = AS0
1054              ASF = AS (R2, R20, AS0, NF)
1055 !
1056           END IF
1057 !
1058        END IF
1059 !
1060 ! ..Final value of alpha_s
1061 !
1062        MSTWALPHAS = 4.D0*PI*ASF
1063
1064        RETURN
1065        END
1066 !
1067 ! =================================================================av==
1068
1069
1070 ! =====================================================================
1071 !
1072 ! ..The threshold matching of the QCD coupling in the MS(bar) scheme,
1073 !    a_s = alpha_s(mu_r^2)/(4 pi),  for  NF -> NF + 1  active flavours
1074 !    up to order a_s^4 (NNNLO).
1075 !
1076 ! ..The value  ASNF  of a_s for NF flavours at the matching scale, the
1077 !    logarithm  LOGRH = ln (mu_r^2/m_H^2) -- where m_H is the pole mass
1078 !    of the heavy quark -- and  NF  are passed as arguments to the
1079 !    function  ASNF1.  The order of the expansion  NAORD  (defined as
1080 !    the 'n' in N^nLO) is provided by the common-block  ASPAR.
1081 !
1082 ! ..The matching coefficients are inverted from Chetyrkin, Kniehl and
1083 !    Steinhauser, Phys. Rev. Lett. 79 (1997) 2184. The QCD colour
1084 !    factors have been hard-wired in these results. The lowest integer
1085 !    values of the Zeta function are given by the common-block  RZETA.
1086 !
1087 ! =====================================================================
1088 !
1089 !
1090       DOUBLE PRECISION FUNCTION ASNF1 (ASNF, LOGRH, NF)
1091 !
1092       IMPLICIT NONE
1093       INTEGER NF, NAORD, NASTPS, PRVCLL, K1, K2
1094       DOUBLE PRECISION ASNF,LOGRH,ZETA,CMC,CMCI30,CMCF30,CMCF31, &
1095      &     CMCI31,ASP,LRHP
1096
1097       DIMENSION CMC(3,0:3)
1098 !
1099 ! ---------------------------------------------------------------------
1100 !
1101 ! ..Input common-blocks
1102 !
1103       COMMON / ASPAR  / NAORD, NASTPS
1104       COMMON / RZETA  / ZETA(6)
1105 !
1106 ! ..Variables to be saved for the next call
1107 !
1108       SAVE CMC, CMCI30, CMCF30, CMCF31, CMCI31, PRVCLL
1109 !
1110 ! ---------------------------------------------------------------------
1111 !
1112 ! ..The coupling-constant matching coefficients (CMC's) up to NNNLO
1113 !   (calculated and saved in the first call of this routine)
1114 !
1115        IF (PRVCLL .NE. 1) THEN
1116 !
1117          CMC(1,0) =  0.D0
1118          CMC(1,1) =  2./3.D0
1119 !
1120          CMC(2,0) = 14./3.D0
1121          CMC(2,1) = 38./3.D0
1122          CMC(2,2) =  4./9.D0
1123 !
1124          CMCI30 = + 80507./432.D0 * ZETA(3) + 58933./1944.D0 &
1125      &            + 128./3.D0 * ZETA(2) * (1.+ DLOG(2.D0)/3.D0)
1126          CMCF30 = - 64./9.D0 * (ZETA(2) + 2479./3456.D0)
1127          CMCI31 =   8941./27.D0
1128          CMCF31 = - 409./27.D0
1129          CMC(3,2) = 511./9.D0
1130          CMC(3,3) = 8./27.D0
1131 !
1132          PRVCLL = 1
1133 !
1134        END IF
1135 !
1136 ! ---------------------------------------------------------------------
1137 !
1138 ! ..The N_f dependent CMC's, and the alpha_s matching at order NAORD
1139 !
1140        CMC(3,0) = CMCI30 + NF * CMCF30
1141        CMC(3,1) = CMCI31 + NF * CMCF31
1142 !
1143        ASNF1 = ASNF
1144        IF (NAORD .EQ. 0) GOTO 1
1145        ASP   = ASNF
1146 !
1147        DO 11 K1 = 1, NAORD
1148          ASP = ASP * ASNF
1149          LRHP = 1.D0
1150 !
1151        DO 12 K2 = 0, K1
1152          ASNF1 = ASNF1 + ASP * CMC(K1,K2) * LRHP
1153          LRHP = LRHP * LOGRH
1154 !
1155   12   CONTINUE
1156   11   CONTINUE
1157 !
1158 ! ---------------------------------------------------------------------
1159 !
1160    1   RETURN
1161        END
1162
1163 !
1164 ! =================================================================av==
1165 !
1166 ! ..The subroutine  EVNFTHR  for the evolution of  a_s = alpha_s/(4 pi)
1167 !    from a three-flavour initial scale to the four- to six-flavour
1168 !    thresholds (identified with the squares of the corresponding quark
1169 !    masses).  The results are written to the common-block  ASFTHR.
1170 !
1171 ! ..The input scale  M20 = mu_(f,0)^2  and the corresponding value
1172 !    AS0  of a_s  are provided by  ASINP.  The fixed scale logarithm
1173 !    LOGFR = ln (mu_f^2/mu_r^2) is specified in  FRRAT.  The alpha_s
1174 !    matching is done by the function ASNF1.
1175 !
1176 ! =====================================================================
1177 !
1178 !
1179        SUBROUTINE EVNFTHR (MC2, MB2, MT2)
1180 !
1181        IMPLICIT NONE
1182        DOUBLE PRECISION MC2, MB2, MT2, M20, M2C, M2B, M2T, R20, R2C, &
1183      &                  R2B, R2T, AS, ASNF1, AS0, ASC, ASB, AST, &
1184      &                  ASC3, ASB4, AST5, LOGFR, SC, SB, ST
1185 !
1186 ! ---------------------------------------------------------------------
1187 !
1188 ! ..Input common blocks
1189 !
1190        COMMON / ASINP  / AS0, M20
1191        COMMON / FRRAT  / LOGFR
1192 !
1193 ! ..Output common blocks
1194 !
1195        COMMON / ASFTHR / ASC, M2C, ASB, M2B, AST, M2T
1196
1197 ! ---------------------------------------------------------------------
1198 !
1199 ! ..Coupling constants at and evolution distances to/between thresholds
1200 !
1201        R20 = M20 * EXP(-LOGFR)
1202 !
1203 ! ..Charm
1204 !
1205        M2C  = MC2
1206        R2C  = M2C * R20/M20
1207        ASC3 = AS (R2C, R20, AS0, 3)
1208        SC   = LOG (AS0 / ASC3)
1209        ASC  = ASNF1 (ASC3, -LOGFR, 3)
1210 !
1211 ! ..Bottom
1212 !
1213        M2B  = MB2
1214        R2B  = M2B * R20/M20
1215        ASB4 = AS (R2B, R2C, ASC, 4)
1216        SB   = LOG (ASC / ASB4)
1217        ASB  = ASNF1 (ASB4, -LOGFR, 4)
1218 !
1219 ! ..Top
1220 !
1221        M2T  = MT2
1222        R2T  = M2T * R20/M20
1223        AST5 = AS (R2T, R2B, ASB, 5)
1224        ST   = LOG (ASB / AST5)
1225        AST  = ASNF1 (AST5, -LOGFR, 5)
1226
1227        RETURN
1228        END
1229
1230 !
1231 ! =================================================================av==
1232 !
1233 ! ..The running coupling of QCD,
1234 !
1235 !         AS  =  a_s  =  alpha_s(mu_r^2)/(4 pi),
1236 !
1237 !    obtained by integrating the evolution equation for a fixed number
1238 !    of massless flavours  NF.  Except at leading order (LO),  AS  is
1239 !    obtained using a fourth-order Runge-Kutta integration.
1240 !
1241 ! ..The initial and final scales  R20  and  R2,  the value  AS0  at
1242 !    R20, and  NF  are passed as function arguments.  The coefficients
1243 !    of the beta function up to  a_s^5 (N^3LO)  are provided by the
1244 !    common-block  BETACOM.  The order of the expansion  NAORD (defined
1245 !    as the 'n' in N^nLO) and the number of steps  NASTPS  for the
1246 !    integration beyond LO are given by the common-block  ASPAR.
1247 !
1248 ! =====================================================================
1249 !
1250 !
1251       DOUBLE PRECISION FUNCTION AS (R2, R20, AS0, NF)
1252 !
1253       IMPLICIT NONE
1254       INTEGER NFMIN, NFMAX, NF, NAORD, NASTPS, K1
1255       DOUBLE PRECISION R2, R20, AS0, SXTH, BETA0, BETA1, BETA2, BETA3, &
1256      &     FBETA1,FBETA2,FBETA3,A,LRRAT,DLR,XK0,XK1,XK2,XK3
1257       PARAMETER (NFMIN = 3, NFMAX = 6)
1258       PARAMETER ( SXTH = 0.166666666666666D0 )
1259 !
1260 ! ---------------------------------------------------------------------
1261 !
1262 ! ..Input common-blocks
1263 !
1264        COMMON / ASPAR  / NAORD, NASTPS
1265        COMMON / BETACOM   / BETA0 (NFMIN:NFMAX), BETA1 (NFMIN:NFMAX), &
1266      &                   BETA2 (NFMIN:NFMAX), BETA3 (NFMIN:NFMAX)
1267 !
1268 ! ..The beta functions FBETAn at N^nLO for n = 1, 2, and 3
1269 !
1270        FBETA1(A) = - A**2 * ( BETA0(NF) + A *   BETA1(NF) )
1271        FBETA2(A) = - A**2 * ( BETA0(NF) + A * ( BETA1(NF) &
1272      &                        + A * BETA2(NF) ) )
1273        FBETA3(A) = - A**2 * ( BETA0(NF) + A * ( BETA1(NF) &
1274      &                        + A * (BETA2(NF) + A * BETA3(NF)) ) )
1275 !
1276 ! ---------------------------------------------------------------------
1277 !
1278 ! ..Initial value, evolution distance and step size
1279 !
1280        AS = AS0
1281        LRRAT = LOG (R2/R20)
1282        DLR = LRRAT / NASTPS
1283 !
1284 ! ..Solution of the evolution equation depending on  NAORD
1285 !   (fourth-order Runge-Kutta beyond the leading order)
1286 !
1287        IF (NAORD .EQ. 0) THEN
1288 !
1289          AS = AS0 / (1.+ BETA0(NF) * AS0 * LRRAT)
1290 !
1291        ELSE IF (NAORD .EQ. 1) THEN
1292 !
1293        DO 2 K1 = 1, NASTPS
1294          XK0 = DLR * FBETA1 (AS)
1295          XK1 = DLR * FBETA1 (AS + 0.5 * XK0)
1296          XK2 = DLR * FBETA1 (AS + 0.5 * XK1)
1297          XK3 = DLR * FBETA1 (AS + XK2)
1298          AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
1299   2    CONTINUE
1300 !
1301        ELSE IF (NAORD .EQ. 2) THEN
1302 !
1303        DO 3 K1 = 1, NASTPS
1304          XK0 = DLR * FBETA2 (AS)
1305          XK1 = DLR * FBETA2 (AS + 0.5 * XK0)
1306          XK2 = DLR * FBETA2 (AS + 0.5 * XK1)
1307          XK3 = DLR * FBETA2 (AS + XK2)
1308          AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
1309   3    CONTINUE
1310 !
1311        ELSE IF (NAORD .EQ. 3) THEN
1312 !
1313        DO 4 K1 = 1, NASTPS
1314          XK0 = DLR * FBETA3 (AS)
1315          XK1 = DLR * FBETA3 (AS + 0.5 * XK0)
1316          XK2 = DLR * FBETA3 (AS + 0.5 * XK1)
1317          XK3 = DLR * FBETA3 (AS + XK2)
1318          AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
1319   4    CONTINUE
1320        END IF
1321 !
1322 ! ---------------------------------------------------------------------
1323 !
1324        RETURN
1325        END
1326
1327 !
1328 ! =================================================================av==
1329 !
1330 ! ..The subroutine BETAFCT for the coefficients  BETA0...BETA3  of the
1331 !    beta function of QCD up to order alpha_s^5 (N^3LO), normalized by
1332 !
1333 !        d a_s / d ln mu_r^2  =  - BETA0 a_s^2 - BETA1 a_s^3 - ...
1334 !
1335 !    with  a_s = alpha_s/(4*pi).
1336 !
1337 ! ..The MSbar coefficients are written to the common-block BETACOM for
1338 !   NF = 3...6  (parameters NFMIN, NFMAX) quark flavours.
1339 !
1340 ! ..The factors CF, CA and TF  are taken from the common-block  COLOUR.
1341 !    Beyond NLO the QCD colour factors are hard-wired in this routine,
1342 !    and the numerical coefficients are truncated to six digits.
1343 !
1344 ! =====================================================================
1345 !
1346 !
1347        SUBROUTINE BETAFCT
1348 !
1349        IMPLICIT DOUBLE PRECISION (A - Z)
1350        INTEGER NFMIN, NFMAX, NF
1351        PARAMETER (NFMIN = 3, NFMAX = 6)
1352 !
1353 ! ---------------------------------------------------------------------
1354 !
1355 ! ..Input common-block
1356 !
1357        COMMON / COLOUR / CF, CA, TR
1358 !
1359 ! ..Output common-block
1360 !
1361        COMMON / BETACOM   / BETA0 (NFMIN:NFMAX), BETA1 (NFMIN:NFMAX), &
1362      &                   BETA2 (NFMIN:NFMAX), BETA3 (NFMIN:NFMAX)
1363
1364 !
1365 ! ---------------------------------------------------------------------
1366 !
1367 ! ..The full LO and NLO coefficients
1368 !
1369        B00 =  11./3.D0 * CA
1370        B01 =  -4./3.D0 * TR
1371        B10 =  34./3.D0 * CA**2
1372        B11 = -20./3.D0 * CA*TR - 4.* CF*TR
1373 !
1374 ! ..Flavour-number loop and output to the array
1375 !
1376        DO 1 NF = NFMIN, NFMAX
1377 !
1378        BETA0(NF) = B00 + B01 * NF
1379        BETA1(NF) = B10 + B11 * NF
1380 !
1381        BETA2(NF) = 1428.50 - 279.611 * NF + 6.01852 * NF**2
1382        BETA3(NF) = 29243.0 - 6946.30 * NF + 405.089 * NF**2 &
1383      &             + 1.49931 * NF**3
1384 !
1385 ! ---------------------------------------------------------------------
1386 !
1387   1    CONTINUE
1388 !
1389        RETURN
1390        END
1391 !
1392 ! =================================================================av==
1393
1394
1395 !--   G.W. MSTWDZEROX taken from CERNLIB to find the zero of a function.
1396       DOUBLE PRECISION FUNCTION MSTWDZEROX(A0,B0,EPS,MAXF,F,MODE)
1397       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1398 !     Based on
1399 !
1400 !        J.C.P. Bus and T.J. Dekker, Two Efficient Algorithms with
1401 !        Guaranteed Convergence for Finding a Zero of a Function,
1402 !        ACM Trans. Math. Software 1 (1975) 330-345.
1403 !
1404 !        (MODE = 1: Algorithm M;    MODE = 2: Algorithm R)
1405       CHARACTER*80 ERRTXT
1406       LOGICAL LMT
1407       DIMENSION IM1(2),IM2(2),LMT(2)
1408       PARAMETER (Z1 = 1, HALF = Z1/2)
1409       DATA IM1 /2,3/, IM2 /-1,3/
1410       MSTWDZEROX = 0.D0             ! G.W. to prevent compiler warning
1411       IF(MODE .NE. 1 .AND. MODE .NE. 2) THEN
1412        C=0
1413        WRITE(ERRTXT,101) MODE
1414        WRITE(6,*) ERRTXT
1415        GOTO 99
1416       ENDIF
1417       FA=F(B0)
1418       FB=F(A0)
1419       IF(FA*FB .GT. 0) THEN
1420        C=0
1421        WRITE(ERRTXT,102) A0,B0
1422        WRITE(6,*) ERRTXT
1423        GOTO 99
1424       ENDIF
1425       ATL=ABS(EPS)
1426       B=A0
1427       A=B0
1428       LMT(2)=.TRUE.
1429       MF=2
1430     1 C=A
1431       FC=FA
1432     2 IE=0
1433     3 IF(ABS(FC) .LT. ABS(FB)) THEN
1434        IF(C .NE. A) THEN
1435         D=A
1436         FD=FA
1437        END IF
1438        A=B
1439        B=C
1440        C=A
1441        FA=FB
1442        FB=FC
1443        FC=FA
1444       END IF
1445       TOL=ATL*(1+ABS(C))
1446       H=HALF*(C+B)
1447       HB=H-B
1448       IF(ABS(HB) .GT. TOL) THEN
1449        IF(IE .GT. IM1(MODE)) THEN
1450         W=HB
1451        ELSE
1452         TOL=TOL*SIGN(Z1,HB)
1453         P=(B-A)*FB
1454         LMT(1)=IE .LE. 1
1455         IF(LMT(MODE)) THEN
1456          Q=FA-FB
1457          LMT(2)=.FALSE.
1458         ELSE
1459          FDB=(FD-FB)/(D-B)
1460          FDA=(FD-FA)/(D-A)
1461          P=FDA*P
1462          Q=FDB*FA-FDA*FB
1463         END IF
1464         IF(P .LT. 0) THEN
1465          P=-P
1466          Q=-Q
1467         END IF
1468         IF(IE .EQ. IM2(MODE)) P=P+P
1469         IF(P .EQ. 0 .OR. P .LE. Q*TOL) THEN
1470          W=TOL
1471         ELSEIF(P .LT. HB*Q) THEN
1472          W=P/Q
1473         ELSE
1474          W=HB
1475         END IF
1476        END IF
1477        D=A
1478        A=B
1479        FD=FA
1480        FA=FB
1481        B=B+W
1482        MF=MF+1
1483        IF(MF .GT. MAXF) THEN
1484         WRITE(6,*) "Error in MSTWDZEROX: TOO MANY FUNCTION CALLS"
1485         GOTO 99
1486        ENDIF
1487        FB=F(B)
1488        IF(FB .EQ. 0 .OR. SIGN(Z1,FC) .EQ. SIGN(Z1,FB)) GOTO 1
1489        IF(W .EQ. HB) GOTO 2
1490        IE=IE+1
1491        GOTO 3
1492       END IF
1493       MSTWDZEROX=C
1494    99 CONTINUE
1495       RETURN
1496   101 FORMAT('Error in MSTWDZEROX: MODE = ',I3,' ILLEGAL')
1497   102 FORMAT('Error in MSTWDZEROX: F(A) AND F(B) HAVE THE SAME SIGN, A = ', &
1498      &     1P,D15.8,', B = ',D15.8)
1499       END
1500
1501 ! ---------------------------------------------------------------------