bin/BCVEGPY/outerpdf.F

0001 c============================================================================

0002 c                             april 10, 2002, v6.01

0003 c                             february 23, 2003, v6.1

0004 c

0005 c   ref[1]: "new generation of parton distributions with uncertainties from global qcd analysis"

0006 c       by: j. pumplin, d.r. stump, j.huston, h.l. lai, p. nadolsky, w.k. tung

0007 c       jhep 0207:012(2002), hep-ph/0201195

0008 c

0009 c   ref[2]: "inclusive jet production, parton distributions, and the search for new physics"

0010 c       by : d. stump, j. huston, j. pumplin, w.k. tung, h.l. lai, s. kuhlmann, j. owens

0011 c       hep-ph/0303013

0012 c

0013 c   this package contains

0014 c   (1) 4 standard sets of cteq6 pdf's (cteq6m, cteq6d, cteq6l, cteq6l1) ;

0015 c   (2) 40 up/down sets (with respect to cteq6m) for uncertainty studies from ref[1];

0016 c   (3) updated version of the above: cteq6.1m and its 40 up/down eigenvector sets from ref[2].

0017 c

0018 c  the cteq6.1m set provides a global fit that is almost equivalent in every respect

0019 c  to the published cteq6m, ref[1], although some parton distributions (e.g., the gluon)

0020 c  may deviate from cteq6m in some kinematic ranges by amounts that are well within the

0021 c  specified uncertainties.

0022 c  the more significant improvements of the new version are associated with some of the

0023 c  40 eigenvector sets, which are made more symmetrical and reliable in (3), compared to (2).

0024
0025 c  details about calling convention are:

0026 c ---------------------------------------------------------------------------

0027 c  iset   pdf-set     description       alpha_s(mz)**lam4  lam5   table_file

0028 c ===========================================================================

0029 c standard, "best-fit", sets:

0030 c --------------------------

0031 c   1    cteq6m   standard msbar scheme   0.118     326   226    cteq6m.tbl

0032 c   2    cteq6d   standard dis scheme     0.118     326   226    cteq6d.tbl

0033 c   3    cteq6l   leading order           0.118**   326** 226    cteq6l.tbl

0034 c   4    cteq6l1  leading order           0.130**   215** 165    cteq6l1.tbl

0035 c ============================================================================

0036 c for uncertainty calculations using eigenvectors of the hessian:

0037 c ---------------------------------------------------------------

0038 c     central + 40 up/down sets along 20 eigenvector directions

0039 c                             -----------------------------

0040 c                original version, ref[1]:  central fit: cteq6m (=cteq6m.00)

0041 c                             -----------------------

0042 c  1xx  cteq6m.xx  +/- sets               0.118     326   226    cteq6m1xx.tbl

0043 c        where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.

0044 c        e.g. 100      is cteq6m.00 (=cteq6m),

0045 c             101/102 are cteq6m.01/02, +/- sets of 1st eigenvector, ... etc.

0046 c                              -----------------------

0047 c                updated version, ref[2]:  central fit: cteq6.1m (=cteq61.00)

0048 c                              -----------------------

0049 c  2xx  cteq61.xx  +/- sets               0.118     326   226    ctq61.xx.tbl

0050 c        where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.

0051 c        e.g. 200      is cteq61.00 (=cteq6.1m),

0052 c             201/202 are cteq61.01/02, +/- sets of 1st eigenvector, ... etc.

0053 c ===========================================================================

0054 c   ** all fits are obtained by using the same coupling strength

0055 c   \alpha_s(mz)=0.118 and the nlo running \alpha_s formula, except cteq6l1

0056 c   which uses the lo running \alpha_s and its value determined from the fit.

0057 c   for the lo fits, the evolution of the pdf and the hard cross sections are

0058 c   calculated at lo.  more detailed discussions are given in the references.

0059 c

0060 c   the table grids are generated for 10^-6 < x < 1 and 1.3 < q < 10,000 (gev).

0061 c   pdf values outside of the above range are returned using extrapolation.

0062 c   lam5 (lam4) represents lambda value (in mev) for 5 (4) flavors.

0063 c   the matching alpha_s between 4 and 5 flavors takes place at q=4.5 gev,

0064 c   which is defined as the bottom quark mass, whenever it can be applied.

0065 c

0066 c   the table_files are assumed to be in the working directory.

0067 c

0068 c   before using the pdf, it is necessary to do the initialization by

0069 c       call setctq6(iset)

0070 c   where iset is the desired pdf specified in the above table.

0071 c

0072 c   the function ctq6pdf (iparton, x, q)

0073 c   returns the parton distribution inside the proton for parton [iparton]

0074 c   at [x] bjorken_x and scale [q] (gev) in pdf set [iset].

0075 c   iparton  is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5)

0076 c                            for (b, c, s, d, u, g, u_bar, ..., b_bar),

0077 c

0078 c   for detailed information on the parameters used, e.q. quark masses,

0079 c   qcd lambda, ... etc.,  see info lines at the beginning of the

0080 c   table_files.

0081 c

0082 c   these programs, as provided, are in double precision.  by removing the

0083 c   "implicit double precision" lines, they can also be run in single

0084 c   precision.

0085 c

0086 c   if you have detailed questions concerning these cteq6 distributions,

0087 c   or if you find problems/bugs using this package, direct inquires to

0088 c   pumplin@pa.msu.edu or tung@pa.msu.edu.

0089 c

0090 c===========================================================================

0091
0092       double precision function ctq6pdf (iparton, x, q)
0093       implicit double precision (a-h,o-z)
0094       logical warn
0095       common
0096      > / ctqpar2 / nx, nt, nfmx
0097      > / qcdtable /  alambda, nfl, iorder
0098
0099       data warn /.true./
0100       save warn
0101
0102       if (x .lt. 0d0 .or. x .gt. 1d0) then
0103         print *, 'x out of range in ctq6pdf: ', x
0104         stop
0105       endif
0106       if (q .lt. alambda) then
0107         print *, 'q out of range in ctq6pdf: ', q
0108         stop
0109       endif
0110       if ((iparton .lt. -nfmx .or. iparton .gt. nfmx)) then
0111          if (warn) then
0112 c        put a warning for calling extra flavor.

0113              warn = .false.
0114              print *, 'warning: iparton out of range in ctq6pdf! '
0115              print *, 'iparton, mxflvn0: ', iparton, nfmx
0116          endif
0117          ctq6pdf = 0d0
0118          return
0119       endif
0120
0121       ctq6pdf = partonx6 (iparton, x, q)
0122
0123       if(ctq6pdf.lt.1.0d-16)  ctq6pdf = 0.0d0
0124
0125       return
0126
0127 c                             ********************

0128       end
0129
0130       subroutine setctq6 (iset)
0131       implicit double precision (a-h,o-z)
0132       parameter (isetmax0=5)
0133       character flnm(isetmax0)*6, nn*3, tablefile*140
0134       data (flnm(i), i=1,isetmax0)
0135      > / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l','ctq61.'/
0136       data isetold, isetmin0, isetmin1, isetmax1 /-987,1,100,140/
0137       data isetmin2,isetmax2 /200,240/
0138       save
0139 c             if data file not initialized, do so.

0140       character cmsdir*82
0141       call getenv('CMSSW_BASE',cmsdir)
0142       kspace=index(cmsdir,' ')
0143       if(iset.ne.isetold) then
0144          iu= nextun()
0145          if (iset.ge.isetmin0 .and. iset.le.3) then
0146             tablefile=cmsdir(1:kspace-1)
0147      $           //'/src/GeneratorInterface/GenExtensions/data/'
0148      $           //flnm(iset)//'.tbl'
0149          elseif (iset.eq.4) then
0150             tablefile=cmsdir(1:kspace-1)
0151      $           //'/src/GeneratorInterface/GenExtensions/data/'
0152      $           //flnm(iset)//'1.tbl'
0153          elseif (iset.ge.isetmin1 .and. iset.le.isetmax1) then
0154             write(nn,'(i3)') iset
0155             tablefile=cmsdir(1:kspace-1)
0156      $           //'/src/GeneratorInterface/GenExtensions/data/'
0157      $           //flnm(1)//nn//'.tbl'
0158          elseif (iset.ge.isetmin2 .and. iset.le.isetmax2) then
0159             write(nn,'(i3)') iset
0160             tablefile=cmsdir(1:kspace-1)
0161      $           //'/src/GeneratorInterface/GenExtensions/data/'
0162      $           //flnm(5)//nn(2:3)//'.tbl'
0163          else
0164             print *, 'invalid iset number in setctq6 :', iset
0165             stop
0166          endif
0167 c        print *, 'pdf file ', tablefile

0168          open(iu, file=tablefile, status='old', err=100)
0169          call readtbl (iu)
0170          close (iu)
0171          isetold=iset
0172       endif
0173       return
0174
0175  100  print *, ' data file ', tablefile, ' cannot be opened '
0176      >//'in setctq6!!'
0177       stop
0178 c                             ********************

0179       end
0180
0181       subroutine readtbl (nu)
0182       implicit double precision (a-h,o-z)
0183       character line*80
0184       parameter (mxx = 96, mxq = 20, mxf = 5)
0185       parameter (mxpqx = (mxf + 3) * mxq * mxx)
0186       common
0187      > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
0188      > / ctqpar2 / nx, nt, nfmx
0189      > / xqrange / qini, qmax, xmin
0190      > / qcdtable /  alambda, nfl, iorder
0191      > / masstbl / amass(6)
0192
0193       read  (nu, '(a)') line
0194       read  (nu, '(a)') line
0195       read  (nu, *) dr, fl, al, (amass(i),i=1,6)
0196       iorder = nint(dr)
0197       nfl = nint(fl)
0198       alambda = al
0199
0200       read  (nu, '(a)') line
0201       read  (nu, *) nx,  nt, nfmx
0202
0203       read  (nu, '(a)') line
0204       read  (nu, *) qini, qmax, (tv(i), i =0, nt)
0205
0206       read  (nu, '(a)') line
0207       read  (nu, *) xmin, (xv(i), i =0, nx)
0208
0209       do 11 iq = 0, nt
0210          tv(iq) = log(log (tv(iq) /al))
0211    11 continue
0212 c

0213 c                  since quark = anti-quark for nfl>2 at this stage,

0214 c                  we read  out only the non-redundent data points

0215 c     no of flavors = nfmx (sea) + 1 (gluon) + 2 (valence)

0216
0217       nblk = (nx+1) * (nt+1)
0218       npts =  nblk  * (nfmx+3)
0219       read  (nu, '(a)') line
0220       read  (nu, *, iostat=iret) (upd(i), i=1,npts)
0221
0222       return
0223 c                        ****************************

0224       end
0225
0226       function nextun()
0227 c                                 returns an unallocated fortran i/o unit.

0228       logical ex
0229 c

0230       do 10 n = 10, 300
0231          inquire (unit=n, opened=ex)
0232          if (.not. ex) then
0233             nextun = n
0234             return
0235          endif
0236  10   continue
0237       stop ' there is no available i/o unit. '
0238 c               *************************

0239       end
0240 c

0241
0242       subroutine polint (xa,ya,n,x,y,dy)
0243
0244       implicit double precision (a-h, o-z)
0245 c                                        adapted from "numerical recipes"

0246       parameter (nmax=10)
0247       dimension xa(n),ya(n),c(nmax),d(nmax)
0248       ns=1
0249       dif=abs(x-xa(1))
0250       do 11 i=1,n
0251         dift=abs(x-xa(i))
0252         if (dift.lt.dif) then
0253           ns=i
0254           dif=dift
0255         endif
0256         c(i)=ya(i)
0257         d(i)=ya(i)
0258 11    continue
0259       y=ya(ns)
0260       ns=ns-1
0261       do 13 m=1,n-1
0262         do 12 i=1,n-m
0263           ho=xa(i)-x
0264           hp=xa(i+m)-x
0265           w=c(i+1)-d(i)
0266           den=ho-hp
0267 c pause is not supported anymore

0268 c          if(den.eq.0.)pause

0269           if (den.eq.0.) then
0270             WRITE (*,*) 'den is 0'
0271             READ (*,'()')
0272           endif
0273           den=w/den
0274           d(i)=hp*den
0275           c(i)=ho*den
0276 12      continue
0277         if (2*ns.lt.n-m)then
0278           dy=c(ns+1)
0279         else
0280           dy=d(ns)
0281           ns=ns-1
0282         endif
0283         y=y+dy
0284 13    continue
0285       return
0286       end
0287
0288       double precision function partonx6 (iprtn, xx, qq)
0289
0290 c  given the parton distribution function in the array u in

0291 c  common / pevldt / , this routine interpolates to find

0292 c  the parton distribution at an arbitray point in x and q.

0293 c

0294       implicit double precision (a-h,o-z)
0295
0296       parameter (mxx = 96, mxq = 20, mxf = 5)
0297       parameter (mxqx= mxq * mxx,   mxpqx = mxqx * (mxf+3))
0298
0299       common
0300      > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
0301      > / ctqpar2 / nx, nt, nfmx
0302      > / xqrange / qini, qmax, xmin
0303
0304       dimension fvec(4), fij(4)
0305       dimension xvpow(0:mxx)
0306       data onep / 1.00001 /
0307       data xpow / 0.3d0 /       !**** choice of interpolation variable
0308       data nqvec / 4 /
0309       data ientry / 0 /
0310       save ientry,xvpow
0311
0312 c store the powers used for interpolation on first call...

0313       if(ientry .eq. 0) then
0314          ientry = 1
0315
0316          xvpow(0) = 0d0
0317          do i = 1, nx
0318             xvpow(i) = xv(i)**xpow
0319          enddo
0320       endif
0321
0322       x = xx
0323       q = qq
0324       tt = log(log(q/al))
0325
0326 c      -------------    find lower end of interval containing x, i.e.,

0327 c                       get jx such that xv(jx) .le. x .le. xv(jx+1)...

0328       jlx = -1
0329       ju = nx+1
0330  11   if (ju-jlx .gt. 1) then
0331          jm = (ju+jlx) / 2
0332          if (x .ge. xv(jm)) then
0333             jlx = jm
0334          else
0335             ju = jm
0336          endif
0337          goto 11
0338       endif
0339 c                     ix    0   1   2      jx  jlx         nx-2     nx

0340 c                           |---|---|---|...|---|-x-|---|...|---|---|

0341 c                     x     0  xmin               x                 1

0342 c

0343       if     (jlx .le. -1) then
0344         print '(a,1pe12.4)', 'severe error: x <= 0 in partonx6! x = ', x
0345         stop
0346       elseif (jlx .eq. 0) then
0347          jx = 0
0348       elseif (jlx .le. nx-2) then
0349
0350 c                for interrior points, keep x in the middle, as shown above

0351          jx = jlx - 1
0352       elseif (jlx.eq.nx-1 .or. x.lt.onep) then
0353
0354 c                  we tolerate a slight over-shoot of one (onep=1.00001),

0355 c              perhaps due to roundoff or whatever, but not more than that.

0356 c                                      keep at least 4 points >= jx

0357          jx = jlx - 2
0358       else
0359         print '(a,1pe12.4)', 'severe error: x > 1 in partonx6! x = ', x
0360         stop
0361       endif
0362 c          ---------- note: jlx uniquely identifies the x-bin; jx does not.

0363
0364 c                       this is the variable to be interpolated in

0365       ss = x**xpow
0366
0367       const1 = 0.
0368       const2 = 0.
0369       const3 = 0.
0370       const4 = 0.
0371       const5 = 0.
0372       const6 = 0.
0373       sy2 = 0.
0374       sy3 = 0.
0375       s23 = 0.
0376
0377       if (jlx.ge.2 .and. jlx.le.nx-2) then
0378
0379 c     initiation work for "interior bins": store the lattice points in s...

0380       svec1 = xvpow(jx)
0381       svec2 = xvpow(jx+1)
0382       svec3 = xvpow(jx+2)
0383       svec4 = xvpow(jx+3)
0384
0385       s12 = svec1 - svec2
0386       s13 = svec1 - svec3
0387       s23 = svec2 - svec3
0388       s24 = svec2 - svec4
0389       s34 = svec3 - svec4
0390
0391       sy2 = ss - svec2
0392       sy3 = ss - svec3
0393
0394 c constants needed for interpolating in s at fixed t lattice points...

0395       const1 = s13/s23
0396       const2 = s12/s23
0397       const3 = s34/s23
0398       const4 = s24/s23
0399       s1213 = s12 + s13
0400       s2434 = s24 + s34
0401       sdet = s12*s34 - s1213*s2434
0402       tmp = sy2*sy3/sdet
0403       const5 = (s34*sy2-s2434*sy3)*tmp/s12
0404       const6 = (s1213*sy2-s12*sy3)*tmp/s34
0405
0406       endif
0407
0408 c         --------------now find lower end of interval containing q, i.e.,

0409 c                          get jq such that qv(jq) .le. q .le. qv(jq+1)...

0410       jlq = -1
0411       ju = nt+1
0412  12   if (ju-jlq .gt. 1) then
0413          jm = (ju+jlq) / 2
0414          if (tt .ge. tv(jm)) then
0415             jlq = jm
0416          else
0417             ju = jm
0418          endif
0419          goto 12
0420        endif
0421
0422       if     (jlq .le. 0) then
0423          jq = 0
0424       elseif (jlq .le. nt-2) then
0425 c                                  keep q in the middle, as shown above

0426          jq = jlq - 1
0427       else
0428 c                         jlq .ge. nt-1 case:  keep at least 4 points >= jq.

0429         jq = nt - 3
0430
0431       endif
0432 c                                   this is the interpolation variable in q

0433
0434       tdet = 0.
0435       t12 = 0.
0436       t13 = 0.
0437       t23 = 0.
0438       t24 = 0.
0439       t34 = 0.
0440       ty2 = 0.
0441       ty3 = 0.
0442       tmp1 = 0.
0443       tmp2 = 0.
0444
0445
0446       if (jlq.ge.1 .and. jlq.le.nt-2) then
0447 c                                        store the lattice points in t...

0448       tvec1 = tv(jq)
0449       tvec2 = tv(jq+1)
0450       tvec3 = tv(jq+2)
0451       tvec4 = tv(jq+3)
0452
0453       t12 = tvec1 - tvec2
0454       t13 = tvec1 - tvec3
0455       t23 = tvec2 - tvec3
0456       t24 = tvec2 - tvec4
0457       t34 = tvec3 - tvec4
0458
0459       ty2 = tt - tvec2
0460       ty3 = tt - tvec3
0461
0462       tmp1 = t12 + t13
0463       tmp2 = t24 + t34
0464
0465       tdet = t12*t34 - tmp1*tmp2
0466
0467       endif
0468
0469
0470 c get the pdf function values at the lattice points...

0471
0472       if (iprtn .ge. 3) then
0473          ip = - iprtn
0474       else
0475          ip = iprtn
0476       endif
0477       jtmp = ((ip + nfmx)*(nt+1)+(jq-1))*(nx+1)+jx+1
0478
0479       do it = 1, nqvec
0480
0481          j1  = jtmp + it*(nx+1)
0482
0483        if (jx .eq. 0) then
0484 c                          for the first 4 x points, interpolate x^2*f(x,q)

0485 c                           this applies to the two lowest bins jlx = 0, 1

0486 c            we can not put the jlx.eq.1 bin into the "interrior" section

0487 c                           (as we do for q), since upd(j1) is undefined.

0488          fij(1) = 0
0489          fij(2) = upd(j1+1) * xv(1)**2
0490          fij(3) = upd(j1+2) * xv(2)**2
0491          fij(4) = upd(j1+3) * xv(3)**2
0492 c

0493 c                 use polint which allows x to be anywhere w.r.t. the grid

0494
0495          call polint (xvpow(0), fij(1), 4, ss, fx, dfx)
0496
0497          if (x .gt. 0d0)  fvec(it) =  fx / x**2
0498 c                                              pdf is undefined for x.eq.0

0499        elseif  (jlx .eq. nx-1) then
0500 c                                                this is the highest x bin:

0501
0502         call polint (xvpow(nx-3), upd(j1), 4, ss, fx, dfx)
0503
0504         fvec(it) = fx
0505
0506        else
0507 c                       for all interior points, use jon's in-line function

0508 c                              this applied to (jlx.ge.2 .and. jlx.le.nx-2)

0509          sf2 = upd(j1+1)
0510          sf3 = upd(j1+2)
0511
0512          g1 =  sf2*const1 - sf3*const2
0513          g4 = -sf2*const3 + sf3*const4
0514
0515          fvec(it) = (const5*(upd(j1)-g1)
0516      &               + const6*(upd(j1+3)-g4)
0517      &               + sf2*sy3 - sf3*sy2) / s23
0518
0519        endif
0520
0521       enddo
0522 c                                   we now have the four values fvec(1:4)

0523 c     interpolate in t...

0524
0525       if (jlq .le. 0) then
0526 c                         1st q-bin, as well as extrapolation to lower q

0527         call polint (tv(0), fvec(1), 4, tt, ff, dfq)
0528
0529       elseif (jlq .ge. nt-1) then
0530 c                         last q-bin, as well as extrapolation to higher q

0531         call polint (tv(nt-3), fvec(1), 4, tt, ff, dfq)
0532       else
0533 c                         interrior bins : (jlq.ge.1 .and. jlq.le.nt-2)

0534 c       which include jlq.eq.1 and jlq.eq.nt-2, since upd is defined for

0535 c                         the full range qv(0:nt)  (in contrast to xv)

0536         tf2 = fvec(2)
0537         tf3 = fvec(3)
0538
0539         g1 = ( tf2*t13 - tf3*t12) / t23
0540         g4 = (-tf2*t34 + tf3*t24) / t23
0541
0542         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
0543      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
0544
0545         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
0546       endif
0547
0548       partonx6 = ff
0549
0550       return
0551       end
0552
0553       double precision function alpcteq(q2, naord)
0554 *********************************************************************

0555 *                                                                   *

0556 *   the alpha_s routine.                                            *

0557 *                                                                   *

0558 *   input :  q2    =  scale in gev**2  (not too low, of course);    *

0559 *            naord =  1 (lo),  2 (nlo).                             *

0560 *                                                                   *

0561 *   output:  alphas_s/(4 pi) for use with the grv(98) partons.      *

0562 *                                                                   *

0563 *******************************************************i*************

0564 *

0565       implicit double precision (a - z)
0566       integer nf, k, i, naord
0567       dimension lambdal (3:6),  lambdan (3:6), q2thr (3)
0568 *

0569 *...heavy quark thresholds and lambda values :

0570       data q2thr   /  1.690,  20.25,  32400. /
0571       data lambdal /  0.247,  0.215,  0.1650, 0.0847  /
0572       data lambdan /  0.350, 0.326, 0.2260, 0.151 /
0573 *

0574 *...determination of the appropriate number of flavours :

0575       nf = 3
0576       do 10 k = 1, 3
0577       if (q2 .gt. q2thr (k)) then
0578          nf = nf + 1
0579       else
0580           go to 20
0581        end if
0582   10   continue
0583 *

0584 *...lo alpha_s and beta function for nlo calculation :

0585   20   b0 = 11.- 2./3.* nf
0586        b1 = 102.- 38./3.* nf
0587        b10 = b1 / (b0*b0)
0588        if (naord .eq. 1) then
0589          lam2 = lambdal (nf) * lambdal (nf)
0590          alp  = 1./(b0 * dlog (q2/lam2))
0591          go to 1
0592        else if (naord .eq. 2) then
0593          lam2 = lambdan (nf) * lambdan (nf)
0594          b1 = 102.- 38./3.* nf
0595          b10 = b1 / (b0*b0)
0596        else
0597          write (6,91)
0598   91     format ('invalid choice for order in alpha_s')
0599          stop
0600        end if
0601 *

0602 *...start value for nlo iteration :

0603        lq2 = dlog (q2 / lam2)
0604        alp = 1./(b0*lq2) * (1.- b10*dlog(lq2)/lq2)
0605 *

0606 *...exact nlo value, found via newton procedure :

0607        do 2 i = 1, 6
0608        xl  = dlog (1./(b0*alp) + b10)
0609        xlp = dlog (1./(b0*alp*1.01) + b10)
0610        xlm = dlog (1./(b0*alp*0.99) + b10)
0611        y  = lq2 - 1./ (b0*alp) + b10 * xl
0612        y1 = (- 1./ (b0*alp*1.01) + b10 * xlp
0613      1       + 1./ (b0*alp*0.99) - b10 * xlp) / (0.02d0*alp)
0614        alp = alp - y/y1
0615   2    continue
0616 *

0617 *...output :

0618   1    alpcteq = alp
0619        return
0620        end
0621 c*******************************************************************

0622 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

0623 *                                                                   *

0624 *     g r v  -  p r o t o n  - p a r a m e t r i z a t i o n s      *

0625 *                                                                   *

0626 *                          1998 update                              *

0627 *                                                                   *

0628 *                  for a detailed explanation see                   *

0629 *                   m. glueck, e. reya, a. vogt :                   *

0630 *        hep-ph/9806404  =  do-th 98/07  =  wue-itp-98-019          *

0631 *                  (to appear in eur. phys. j. c)                   *

0632 *                                                                   *

0633 *   this package contains subroutines returning the light-parton    *

0634 *   distributions in nlo (for the msbar and dis schemes) and lo;    *

0635 *   the respective light-parton, charm, and bottom contributions    *

0636 *   to f2(electromagnetic); and the scale dependence of alpha_s.    *

0637 *                                                                   *

0638 *   the parton densities and f2 values are calculated from inter-   *

0639 *   polation grids covering the regions                             *

0640 *         q^2/gev^2  between   0.8   and  1.e6 ( 1.e4 for f2 )      *

0641 *            x       between  1.e-9  and   1.                       *

0642 *   any call outside these regions stops the program execution.     *

0643 *                                                                   *

0644 *   at q^2 = mz^2, alpha_s reads  0.114 (0.125) in nlo (lo); the    *

0645 *   heavy quark thresholds, qh^2 = mh^2, in the beta function are   *

0646 *            mc = 1.4 gev,  mb = 4.5 gev,  mt = 175 gev.            *

0647 *   note that the nlo alpha_s running is different from grv(94).    *

0648 *                                                                   *

0649 *    questions, comments etc to:  avogt@physik.uni-wuerzburg.de     *

0650 *                                                                   *

0651 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

0652 *

0653 *

0654 *

0655 *

0656       subroutine grv98pa (iset, x, q2, uv, dv, us, ds, ss, gl)
0657 *********************************************************************

0658 *                                                                   *

0659 *   the parton routine.                                             *

0660 *                                     __                            *

0661 *   input:   iset =  1 (lo),  2 (nlo, ms), or  3 (nlo, dis)         *

0662 *            x  =  bjorken-x        (between  1.e-9 and 1.)         *

0663 *            q2 =  scale in gev**2  (between  0.8 and 1.e6)         *

0664 *                                                                   *

0665 *   output:  uv = u - u(bar),  dv = d - d(bar),  us = u(bar),       *

0666 *            ds = d(bar),  ss = s = s(bar),  gl = gluon.            *

0667 *            always x times the distribution is returned.           *

0668 *                                                                   *

0669 *   common:  the main program or the calling routine has to have    *

0670 *            a common block  common / intinip / iinip , and the     *

0671 *            integer variable  iinip  has always to be zero when    *

0672 *            grv98pa is called for the first time or when  iset     *

0673 *            has been changed.                                      *

0674 *                                                                   *

0675 *   grids:   1. grv98lo.grid, 2. grv98nlm.grid, 3. grv98nld.grid,   *

0676 *            (1+1809 lines with 6 columns, 4 significant figures)   *

0677 *                                                                   *

0678 *******************************************************i*************

0679 *

0680       implicit double precision (a-h, o-z)
0681       parameter (npart=6, nx=68, nq=27, narg=2)
0682       dimension xuvf(nx,nq), xdvf(nx,nq), xdef(nx,nq), xudf(nx,nq),
0683      1          xsf(nx,nq), xgf(nx,nq), parton (npart,nq,nx-1),
0684      2          qs(nq), xb(nx), xt(narg), na(narg), arrf(nx+nq)
0685       character*80 line
0686       common /intinip/iinip
0687         save xuvf, xdvf, xdef, xudf, xsf, xgf, na, arrf
0688 *

0689 *...bjorken-x and q**2 values of the grid :

0690        data qs / 0.8e0,
0691      1           1.0e0, 1.3e0, 1.8e0, 2.7e0, 4.0e0, 6.4e0,
0692      2           1.0e1, 1.6e1, 2.5e1, 4.0e1, 6.4e1,
0693      3           1.0e2, 1.8e2, 3.2e2, 5.7e2,
0694      4           1.0e3, 1.8e3, 3.2e3, 5.7e3,
0695      5           1.0e4, 2.2e4, 4.6e4,
0696      6           1.0e5, 2.2e5, 4.6e5,
0697      7           1.e6 /
0698        data xb / 1.0e-9, 1.8e-9, 3.2e-9, 5.7e-9,
0699      a           1.0e-8, 1.8e-8, 3.2e-8, 5.7e-8,
0700      b           1.0e-7, 1.8e-7, 3.2e-7, 5.7e-7,
0701      c           1.0e-6, 1.4e-6, 2.0e-6, 3.0e-6, 4.5e-6, 6.7e-6,
0702      1           1.0e-5, 1.4e-5, 2.0e-5, 3.0e-5, 4.5e-5, 6.7e-5,
0703      2           1.0e-4, 1.4e-4, 2.0e-4, 3.0e-4, 4.5e-4, 6.7e-4,
0704      3           1.0e-3, 1.4e-3, 2.0e-3, 3.0e-3, 4.5e-3, 6.7e-3,
0705      4           1.0e-2, 1.4e-2, 2.0e-2, 3.0e-2, 4.5e-2, 0.06, 0.08,
0706      5           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
0707      6           0.3, 0.325, 0.35, 0.375, 0.4,  0.45, 0.5, 0.55,
0708      7           0.6, 0.65,  0.7,  0.75,  0.8,  0.85, 0.9, 0.95, 1. /
0709 *

0710 *...check of x and q2 values :

0711       if ( (x.lt.0.99d-9) .or. (x.gt.1.d0) ) then
0712          write(6,91)
0713   91     format (2x,'parton interpolation: x out of range')
0714          stop
0715       endif
0716       if ( (q2.lt.0.799) .or. (q2.gt.1.01e6) ) then
0717 c         write(6,92)

0718 c  92     format (2x,'parton interpolation: q2 out of range')

0719          q2=1.0d6
0720 c          stop

0721       endif
0722       if (iinip .ne. 0) goto 16
0723 *

0724 *...initialization, if required :

0725 *

0726 *    selection and reading of the grid :

0727 *    (comment: first number in the first line of the grid)

0728       if (iset .eq. 1) then
0729         open (11,file='grv98lo.grid',status='old')   ! 7.332e-05
0730       else if (iset .eq. 2) then
0731         open (11,file='grv98nlm.grid',status='old')  ! 1.015e-04
0732       else if (iset .eq. 3) then
0733         open (11,file='grv98nld.grid',status='old')  ! 1.238e-04
0734       else
0735         write(6,93)
0736   93    format (2x,'no or invalid parton set choice')
0737         stop
0738       end if
0739       iinip = 1
0740       read(11,89) line
0741   89  format(a80)
0742       do 15 m = 1, nx-1
0743       do 15 n = 1, nq
0744       read(11,90) parton(1,n,m), parton(2,n,m), parton(3,n,m),
0745      1            parton(4,n,m), parton(5,n,m), parton(6,n,m)
0746   90  format (6(1pe10.3))
0747   15  continue
0748       close(11)
0749 *

0750 *....arrays for the interpolation subroutine :

0751       do 10 iq = 1, nq
0752       do 20 ix = 1, nx-1
0753         xb0v = xb(ix)**0.5
0754         xb0s = xb(ix)**(-0.2)
0755         xb1 = 1.-xb(ix)
0756         xuvf(ix,iq) = parton(1,iq,ix) / (xb1**3 * xb0v)
0757         xdvf(ix,iq) = parton(2,iq,ix) / (xb1**4 * xb0v)
0758         xdef(ix,iq) = parton(3,iq,ix) / (xb1**7 * xb0v)
0759         xudf(ix,iq) = parton(4,iq,ix) / (xb1**7 * xb0s)
0760         xsf(ix,iq)  = parton(5,iq,ix) / (xb1**7 * xb0s)
0761         xgf(ix,iq)  = parton(6,iq,ix) / (xb1**5 * xb0s)
0762   20  continue
0763         xuvf(nx,iq) = 0.e0
0764         xdvf(nx,iq) = 0.e0
0765         xdef(nx,iq) = 0.e0
0766         xudf(nx,iq) = 0.e0
0767         xsf(nx,iq)  = 0.e0
0768         xgf(nx,iq)  = 0.e0
0769   10  continue
0770       na(1) = nx
0771       na(2) = nq
0772       do 30 ix = 1, nx
0773         arrf(ix) = dlog(xb(ix))
0774   30  continue
0775       do 40 iq = 1, nq
0776         arrf(nx+iq) = dlog(qs(iq))
0777   40  continue
0778 *

0779 *...continuation, if initialization was done previously.

0780 *

0781   16  continue
0782 *

0783 *...interpolation :

0784       xt(1) = dlog(x)
0785       xt(2) = dlog(q2)
0786       x1 = 1.- x
0787       xv = x**0.5
0788       xs = x**(-0.2)
0789       uv = fint(narg,xt,na,arrf,xuvf) * x1**3 * xv
0790       dv = fint(narg,xt,na,arrf,xdvf) * x1**4 * xv
0791       de = fint(narg,xt,na,arrf,xdef) * x1**7 * xv
0792       ud = fint(narg,xt,na,arrf,xudf) * x1**7 * xs
0793       us = 0.5 * (ud - de)
0794       ds = 0.5 * (ud + de)
0795       ss = fint(narg,xt,na,arrf,xsf)  * x1**7 * xs
0796       gl = fint(narg,xt,na,arrf,xgf)  * x1**5 * xs
0797 *

0798       return
0799       end
0800 *

0801 *

0802 *

0803 *

0804       subroutine grv98f2 (iset, x, q2, f2l, f2c, f2b, f2)
0805 *********************************************************************

0806 *                                                                   *

0807 *   the f2 routine.                                                 *

0808 *                                                                   *

0809 *   input :  iset = 1 (lo),  2 (nlo).                               *

0810 *            x  =  bjorken-x        (between  1.e-9 and 1.)         *

0811 *            q2 =  scale in gev**2  (between  0.8 and 1.e4)         *

0812 *                                                                   *

0813 *   output:  f2l = f2(light), f2c = f2(charm), f2b = f2(bottom,)    *

0814 *            f2  = sum, all given for electromagnetic proton dis.   *

0815 *                                                                   *

0816 *   common:  the main program or the calling routine has to have    *

0817 *            a common block  common / intinif / iinif , and the     *

0818 *            integer variable  iinif  has always to be zero when    *

0819 *            grv98f2 is called for the first time or when  iset     *

0820 *            has been changed.                                      *

0821 *                                                                   *

0822 *   grids:   1. grv98lof.grid, 2. grv98nlf.grid.                    *

0823 *            (1+1407 lines with 3 columns, 4 significant figures)   *

0824 *                                                                   *

0825 *******************************************************i*************

0826 *

0827       implicit double precision (a-h, o-z)
0828       parameter (nstrf=3, nx=68, nq=21, narg=2)
0829       dimension xf2lf(nx,nq), xf2cf(nx,nq), xf2bf(nx,nq),
0830      1          strfct (nstrf,nq,nx-1), qs(nq), xb(nx),
0831      3          xt(narg), na(narg), arrf(nx+nq)
0832       character*80 line
0833       common / intinif / iinif
0834       save xf2lf, xf2cf, xf2bf, na, arrf
0835 *

0836 *...bjorken-x and q**2 values of the grid :

0837        data qs / 0.8e0,
0838      1           1.0e0, 1.3e0, 1.8e0, 2.7e0, 4.0e0, 6.4e0,
0839      2           1.0e1, 1.6e1, 2.5e1, 4.0e1, 6.4e1,
0840      3           1.0e2, 1.8e2, 3.2e2, 5.7e2,
0841      4           1.0e3, 1.8e3, 3.2e3, 5.7e3,
0842      5           1.0e4/
0843        data xb / 1.0e-9, 1.8e-9, 3.2e-9, 5.7e-9,
0844      a           1.0e-8, 1.8e-8, 3.2e-8, 5.7e-8,
0845      b           1.0e-7, 1.8e-7, 3.2e-7, 5.7e-7,
0846      c           1.0e-6, 1.4e-6, 2.0e-6, 3.0e-6, 4.5e-6, 6.7e-6,
0847      1           1.0e-5, 1.4e-5, 2.0e-5, 3.0e-5, 4.5e-5, 6.7e-5,
0848      2           1.0e-4, 1.4e-4, 2.0e-4, 3.0e-4, 4.5e-4, 6.7e-4,
0849      3           1.0e-3, 1.4e-3, 2.0e-3, 3.0e-3, 4.5e-3, 6.7e-3,
0850      4           1.0e-2, 1.4e-2, 2.0e-2, 3.0e-2, 4.5e-2, 0.06, 0.08,
0851      5           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
0852      6           0.3, 0.325, 0.35, 0.375, 0.4,  0.45, 0.5, 0.55,
0853      7           0.6, 0.65,  0.7,  0.75,  0.8,  0.85, 0.9, 0.95, 1. /
0854 *

0855 *...check of x and q2 values :

0856       if ( (x.lt.0.99d-9) .or. (x.gt.1.d0) ) then
0857          write(6,91)
0858   91     format (2x,'str.fct. interpolation: x out of range')
0859          stop
0860       endif
0861       if ( (q2.lt.0.799) .or. (q2.gt.1.01e4) ) then
0862          write(6,92)
0863   92     format (2x,'str.fct. interpolation: q2 out of range')
0864          stop
0865       endif
0866       if (iinif .ne. 0) goto 16
0867 *

0868 *...initialization, if required :

0869 *

0870 *    selection and reading of the grid :

0871 *    (comment: first number in the first line of the grid)

0872       if (iset .eq. 1) then
0873         open (11,file='grv98lof.grid',status='old')  !  7.907e-01
0874       else if (iset.eq.2) then
0875         open (11,file='grv98nlf.grid',status='old')  !  9.368e-01
0876       else
0877         write(6,93)
0878   93    format (2x,'no or invalid str.fct. set choice')
0879         stop
0880       end if
0881       iinif = 1
0882       read(11,89) line
0883   89  format(a80)
0884       do 15 m = 1, nx-1
0885       do 15 n = 1, nq
0886       read(11,90) strfct(1,n,m), strfct(2,n,m), strfct(3,n,m)
0887   90  format (3(1pe10.3))
0888   15  continue
0889       close(11)
0890 *

0891 *....arrays for the interpolation subroutine :

0892       do 10 iq = 1, nq
0893       do 20 ix = 1, nx-1
0894         xbs = xb(ix)**0.2
0895         xb1 = 1.-xb(ix)
0896         xf2lf(ix,iq) = strfct(1,iq,ix) / xb1**3 * xbs
0897         xf2cf(ix,iq) = strfct(2,iq,ix) / xb1**7 * xbs
0898         xf2bf(ix,iq) = strfct(3,iq,ix) / xb1**7 * xbs
0899   20  continue
0900         xf2lf(nx,iq) = 0.e0
0901         xf2cf(nx,iq) = 0.e0
0902         xf2bf(nx,iq) = 0.e0
0903   10  continue
0904       na(1) = nx
0905       na(2) = nq
0906       do 30 ix = 1, nx
0907         arrf(ix) = dlog(xb(ix))
0908   30  continue
0909       do 40 iq = 1, nq
0910         arrf(nx+iq) = dlog(qs(iq))
0911   40  continue
0912 *

0913 *...continuation, if initialization was done previously.

0914 *

0915   16  continue
0916 *

0917 *...interpolation :

0918       xt(1) = dlog(x)
0919       xt(2) = dlog(q2)
0920       x1 = 1.- x
0921       xs = x**(-0.2)
0922       f2l = fint(narg,xt,na,arrf,xf2lf) * x1**3 * xs
0923       f2c = fint(narg,xt,na,arrf,xf2cf) * x1**7 * xs
0924       f2b = fint(narg,xt,na,arrf,xf2bf) * x1**7 * xs
0925       f2  = f2l + f2c + f2b
0926 *

0927       return
0928       end
0929 *

0930 *

0931       double precision function fint(narg,arg,nent,ent,table)
0932 *********************************************************************

0933 *                                                                   *

0934 *   the interpolation routine (cern library routine e104)           *

0935 *                                                                   *

0936 *********************************************************************

0937       implicit double precision (a-h, o-z)
0938       dimension arg(narg),nent(narg),ent(*),table(*)
0939       dimension d(5),ncomb(5),ient(5)
0940
0941       kd=1
0942       m=1
0943       ja=1
0944
0945       do 5 ii=1,narg
0946         ncomb(ii)=1
0947         jb=ja-1+nent(ii)
0948         do 2 j=ja,jb
0949           if (arg(ii).le.ent(j)) go to 3
0950     2   continue
0951         j=jb
0952     3   if (j.ne.ja) go to 4
0953         j=j+1
0954     4   jr=j-1
0955         d(ii)=(ent(j)-arg(ii))/(ent(j)-ent(jr))
0956         ient(ii)=j-ja
0957         kd=kd+ient(ii)*m
0958         m=m*nent(ii)
0959         ja=jb+1
0960     5 continue
0961
0962       fint=0.
0963    10 fac=1.
0964       iadr=kd
0965       ifadr=1
0966
0967       do 15 i=1,narg
0968         if (ncomb(i).eq.0) go to 12
0969         fac=fac*(1.-d(i))
0970         go to 15
0971    12   fac=fac*d(i)
0972         iadr=iadr-ifadr
0973         ifadr=ifadr*nent(i)
0974    15 continue
0975
0976       fint=fint+fac*table(iadr)
0977       il=narg
0978    40 if (ncomb(il).eq.0) go to 80
0979       ncomb(il)=0
0980       if (il.eq.narg) go to 10
0981       il=il+1
0982       do 50  k=il,narg
0983    50    ncomb(k)=1
0984       go to 10
0985    80 il=il-1
0986       if(il.ne.0) go to 40
0987
0988         return
0989       end
0990 *

0991 *

0992       double precision function alpgrv(q2, naord)
0993 *********************************************************************

0994 *                                                                   *

0995 *   the alpha_s routine.                                            *

0996 *                                                                   *

0997 *   input :  q2    =  scale in gev**2  (not too low, of course);    *

0998 *            naord =  1 (lo),  2 (nlo).                             *

0999 *                                                                   *

1000 *   output:  alphas_s/(4 pi) for use with the grv(98) partons.      *

1001 *                                                                   *

1002 *******************************************************i*************

1003 *

1004       implicit double precision (a - z)
1005       integer nf, k, i, naord
1006       dimension lambdal (3:6),  lambdan (3:6), q2thr (3)
1007 *

1008 *...heavy quark thresholds and lambda values :

1009       data q2thr   /  1.960,  20.25,  30625. /
1010       data lambdal / 0.2041, 0.1750, 0.1320, 0.0665 /
1011       data lambdan / 0.2994, 0.2460, 0.1677, 0.0678 /
1012 *

1013 *...determination of the appropriate number of flavours :

1014       nf = 3
1015       do 10 k = 1, 3
1016       if (q2 .gt. q2thr (k)) then
1017          nf = nf + 1
1018       else
1019           go to 20
1020        end if
1021   10   continue
1022 *

1023 *...lo alpha_s and beta function for nlo calculation :

1024   20   b0 = 11.- 2./3.* nf
1025        b1 = 102.- 38./3.* nf
1026        b10 = b1 / (b0*b0)
1027        if (naord .eq. 1) then
1028          lam2 = lambdal (nf) * lambdal (nf)
1029          alp  = 1./(b0 * dlog (q2/lam2))
1030          go to 1
1031        else if (naord .eq. 2) then
1032          lam2 = lambdan (nf) * lambdan (nf)
1033          b1 = 102.- 38./3.* nf
1034          b10 = b1 / (b0*b0)
1035        else
1036          write (6,91)
1037   91     format ('invalid choice for order in alpha_s')
1038          stop
1039        end if
1040 *

1041 *...start value for nlo iteration :

1042        lq2 = dlog (q2 / lam2)
1043        alp = 1./(b0*lq2) * (1.- b10*dlog(lq2)/lq2)
1044 *

1045 *...exact nlo value, found via newton procedure :

1046        do 2 i = 1, 6
1047        xl  = dlog (1./(b0*alp) + b10)
1048        xlp = dlog (1./(b0*alp*1.01) + b10)
1049        xlm = dlog (1./(b0*alp*0.99) + b10)
1050        y  = lq2 - 1./ (b0*alp) + b10 * xl
1051        y1 = (- 1./ (b0*alp*1.01) + b10 * xlp
1052      1       + 1./ (b0*alp*0.99) - b10 * xlp) / (0.02d0*alp)
1053        alp = alp - y/y1
1054   2    continue
1055 *

1056 *...output :

1057   1    alpgrv = alp
1058        return
1059        end
1060 c*******************************************************************

1061       subroutine mrstlo(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
1062 c***************************************************************c

1063 c                                                               c

1064 c  this is a package for the new mrst 2001 lo parton            c

1065 c  distributions.                                               c

1066 c  reference: a.d. martin, r.g. roberts, w.j. stirling and      c

1067 c  r.s. thorne, hep-ph/0201xxx                                  c

1068 c                                                               c

1069 c  there is 1 pdf set corresponding to mode = 1                 c

1070 c                                                               c

1071 c  mode=1 gives the default set with lambda(msbar,nf=4) = 0.220 c

1072 c  corresponding to alpha_s(m_z) of 0.130                       c

1073 c  this set reads a grid whose first number is 0.02868          c

1074 c                                                               c

1075 c   this subroutine uses an improved interpolation procedure    c

1076 c   for extracting values of the pdf's from the grid            c

1077 c                                                               c

1078 c         comments to : w.j.stirling@durham.ac.uk               c

1079 c                                                               c

1080 c***************************************************************c

1081       implicit real*8(a-h,o-z)
1082       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
1083       q2=q*q
1084       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print *, 99,q2
1085       if(x.lt.xmin.or.x.gt.xmax)       print *, 98,x
1086           if(mode.eq.1) then
1087         call mrst1(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu)
1088       endif
1089 c  99  format('  warning:  q^2 value is out of range   ','q2= ',e10.5)

1090 c  98  format('  warning:   x  value is out of range   ','x= ',e10.5)

1091       return
1092       end
1093 c********************************************************

1094       subroutine mrst1(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
1095       implicit real*8(a-h,o-z)
1096       parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
1097       real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
1098      .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
1099       real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
1100      .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
1101      .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
1102       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
1103       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
1104      .        1d-4,2d-4,4d-4,6d-4,8d-4,
1105      .        1d-3,2d-3,4d-3,6d-3,8d-3,
1106      .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
1107      .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
1108      .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
1109      .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
1110      .     .8d0,.9d0,1d0/
1111       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
1112      .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
1113      .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
1114      .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
1115      .        1.8d6,3.2d6,5.6d6,1d7/
1116       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
1117       data init/0/
1118       save
1119       xsave=x
1120       q2save=qsq
1121       if(init.ne.0) goto 10
1122         open(unit=33,file='lo2002.dat',status='old')
1123         do 20 n=1,nx-1
1124         do 20 m=1,nq
1125         read(33,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
1126      .            f5(n,m),f7(n,m),f6(n,m),f8(n,m)
1127 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea

1128   20  continue
1129       do 40 m=1,nq
1130       f1(nx,m)=0.d0
1131       f2(nx,m)=0.d0
1132       f3(nx,m)=0.d0
1133       f4(nx,m)=0.d0
1134       f5(nx,m)=0.d0
1135       f6(nx,m)=0.d0
1136       f7(nx,m)=0.d0
1137       f8(nx,m)=0.d0
1138   40  continue
1139       do n=1,nx
1140       xxl(n)=dlog(xx(n))
1141       enddo
1142       do m=1,nq
1143       qql(m)=dlog(qq(m))
1144       enddo
1145
1146       call jeppe1(nx,nq,xxl,qql,f1,cc1)
1147       call jeppe1(nx,nq,xxl,qql,f2,cc2)
1148       call jeppe1(nx,nq,xxl,qql,f3,cc3)
1149       call jeppe1(nx,nq,xxl,qql,f4,cc4)
1150       call jeppe1(nx,nq,xxl,qql,f6,cc6)
1151       call jeppe1(nx,nq,xxl,qql,f8,cc8)
1152
1153       emc2=2.045
1154       emb2=18.5
1155
1156       do 44 m=1,nqc
1157       qqlc(m)=qql(m+nqc0)
1158       do 44 n=1,nx
1159       fc(n,m)=f5(n,m+nqc0)
1160    44 continue
1161       qqlc(1)=dlog(emc2)
1162       call jeppe1(nx,nqc,xxl,qqlc,fc,ccc)
1163
1164       do 45 m=1,nqb
1165       qqlb(m)=qql(m+nqb0)
1166       do 45 n=1,nx
1167       fb(n,m)=f7(n,m+nqb0)
1168    45 continue
1169       qqlb(1)=dlog(emb2)
1170       call jeppe1(nx,nqb,xxl,qqlb,fb,ccb)
1171
1172
1173       init=1
1174    10 continue
1175
1176       xlog=dlog(x)
1177       qsqlog=dlog(qsq)
1178
1179       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
1180       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
1181       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
1182       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
1183       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
1184       call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
1185
1186       chm=0.d0
1187       if(qsq.gt.emc2) then
1188       call jeppe2(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
1189       endif
1190
1191       bot=0.d0
1192       if(qsq.gt.emb2) then
1193       call jeppe2(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
1194       endif
1195
1196       x=xsave
1197       qsq=q2save
1198
1199         close(33)
1200
1201       return
1202    50 format(8f10.5)
1203       end
1204 c************************************************

1205       subroutine jeppe1(nx,my,xx,yy,ff,cc)
1206       implicit real*8(a-h,o-z)
1207       parameter(nnx=49,mmy=37)
1208       dimension xx(nx),yy(my),ff(nx,my),ff1(nnx,mmy),ff2(nnx,mmy),
1209      xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),
1210      xcl(16),cc(nx,my,4,4),iwt(16,16)
1211
1212       data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1213      x            0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
1214      x            -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
1215      x            2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
1216      x            0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
1217      x            0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
1218      x            0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
1219      x            0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
1220      x            -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
1221      x            0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
1222      x            9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
1223      x            -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
1224      x            2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
1225      x            0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
1226      x            -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
1227      x            4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
1228
1229
1230       do 42 m=1,my
1231       dx=xx(2)-xx(1)
1232       ff1(1,m)=(ff(2,m)-ff(1,m))/dx
1233       dx=xx(nx)-xx(nx-1)
1234       ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
1235       do 41 n=2,nx-1
1236       ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
1237      xff(n+1,m))
1238    41 continue
1239    42 continue
1240
1241       do 44 n=1,nx
1242       dy=yy(2)-yy(1)
1243       ff2(n,1)=(ff(n,2)-ff(n,1))/dy
1244       dy=yy(my)-yy(my-1)
1245       ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
1246       do 43 m=2,my-1
1247       ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
1248      xff(n,m+1))
1249    43 continue
1250    44 continue
1251
1252       do 46 m=1,my
1253       dx=xx(2)-xx(1)
1254       ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
1255       dx=xx(nx)-xx(nx-1)
1256       ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
1257       do 45 n=2,nx-1
1258       ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
1259      xff2(n+1,m))
1260    45 continue
1261    46 continue
1262
1263       do 53 n=1,nx-1
1264       do 52 m=1,my-1
1265       d1=xx(n+1)-xx(n)
1266       d2=yy(m+1)-yy(m)
1267       d1d2=d1*d2
1268
1269       yy0(1)=ff(n,m)
1270       yy0(2)=ff(n+1,m)
1271       yy0(3)=ff(n+1,m+1)
1272       yy0(4)=ff(n,m+1)
1273
1274       yy1(1)=ff1(n,m)
1275       yy1(2)=ff1(n+1,m)
1276       yy1(3)=ff1(n+1,m+1)
1277       yy1(4)=ff1(n,m+1)
1278
1279       yy2(1)=ff2(n,m)
1280       yy2(2)=ff2(n+1,m)
1281       yy2(3)=ff2(n+1,m+1)
1282       yy2(4)=ff2(n,m+1)
1283
1284       yy12(1)=ff12(n,m)
1285       yy12(2)=ff12(n+1,m)
1286       yy12(3)=ff12(n+1,m+1)
1287       yy12(4)=ff12(n,m+1)
1288
1289       do 47 k=1,4
1290       z(k)=yy0(k)
1291       z(k+4)=yy1(k)*d1
1292       z(k+8)=yy2(k)*d2
1293       z(k+12)=yy12(k)*d1d2
1294    47 continue
1295
1296       do 49 l=1,16
1297       xxd=0.
1298       do 48 k=1,16
1299       xxd=xxd+iwt(k,l)*z(k)
1300    48 continue
1301       cl(l)=xxd
1302    49 continue
1303       l=0
1304       do 51 k=1,4
1305       do 50 j=1,4
1306       l=l+1
1307       cc(n,m,k,j)=cl(l)
1308    50 continue
1309    51 continue
1310    52 continue
1311    53 continue
1312       return
1313       end
1314 c***************************************

1315       subroutine jeppe2(x,y,nx,my,xx,yy,cc,z)
1316       implicit real*8(a-h,o-z)
1317       dimension xx(nx),yy(my),cc(nx,my,4,4)
1318
1319       n=locx(xx,nx,x)
1320       m=locx(yy,my,y)
1321
1322       t=(x-xx(n))/(xx(n+1)-xx(n))
1323       u=(y-yy(m))/(yy(m+1)-yy(m))
1324
1325       z=0.
1326       do 1 l=4,1,-1
1327       z=t*z+((cc(n,m,l,4)*u+cc(n,m,l,3))*u
1328      .       +cc(n,m,l,2))*u+cc(n,m,l,1)
1329     1 continue
1330       return
1331       end
1332 c**************************************

1333       integer function locx(xx,nx,x)
1334       implicit real*8(a-h,o-z)
1335       dimension xx(nx)
1336       if(x.le.xx(1)) then
1337       locx=1
1338       return
1339       endif
1340       if(x.ge.xx(nx)) then
1341       locx=nx-1
1342       return
1343       endif
1344       ju=nx+1
1345       jl=0
1346     1 if((ju-jl).le.1) go to 2
1347       jm=(ju+jl)/2
1348       if(x.ge.xx(jm)) then
1349       jl=jm
1350       else
1351       ju=jm
1352       endif
1353       go to 1
1354     2 locx=jl
1355       return
1356       end
1357 c*******************************************

1358       real*8 function  polderiv(x1,x2,x3,y1,y2,y3)
1359       implicit real*8(a-h,o-z)
1360       polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
1361      .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
1362       return
1363       end
1364 c**************************************************

1365       double precision function alpmsrt(scale,lambda,iorder)
1366 c.

1367 c.   the qcd coupling used in the mrs analysis.

1368 c.   automatically corrects for nf=3,4,5 with thresholds at m_q.

1369 c.   the following parameters must be supplied:

1370 c.       scale   =  the qcd scale in gev (real)

1371 c.       lambda  =  the 4 flavour msbar lambda parameter in gev (real)

1372 c.       iorder  =  0 for lo, 1 for nlo and 2 for nnlo

1373 c.

1374       implicit real*8 (a-h,o-z)
1375       real*8 lambda
1376       data tol,qsct,qsdt/5d-4,74d0,8.18d0/
1377
1378 c-ap initialize to avoid compiler warnings

1379       alpmsrt =0.
1380
1381       pi=4d0*datan(1d0)
1382       iord=iorder
1383       ith=0
1384       flav=4d0
1385       al=lambda
1386       al2=lambda*lambda
1387       qs=scale*scale
1388       t=dlog(qs/al2)
1389       tt=t
1390 c.

1391       qsctt=qsct/4.
1392       qsdtt=qsdt/4.
1393       if(qs.lt.0.5d0) then   !!  running stops below 0.5
1394           qs=0.5d0
1395           t=dlog(qs/al2)
1396           tt=t
1397       endif
1398
1399       alfqc5=0.
1400       alfqc4=0.
1401       alfqc3=0.
1402
1403       if(qs.gt.qsctt) go        to 12
1404       if(qs.lt.qsdtt) go        to 312
1405    11 continue
1406       b0=11-2.*flav/3.
1407       x1=4.*pi/b0
1408   5   continue
1409       if(iord.eq.0) then
1410       alpmsrt=x1/t
1411       elseif(iord.eq.1) then
1412       b1=102.-38.*flav/3.
1413       x2=b1/b0**2
1414       as2=x1/t*(1.-x2*dlog(t)/t)
1415   95    as=as2
1416         f=-t+x1/as-x2*dlog(x1/as+x2)
1417         fp=-x1/as**2*(1.-x2/(x1/as+x2))
1418         as2=as-f/fp
1419         del=abs(f/fp/as)
1420         if((del-tol).gt.0.) go to 95
1421       alpmsrt=as2
1422       elseif(iord.eq.2) then
1423       alpmsrt=qwikalf(t,2,flav)
1424       endif
1425       if(ith.eq.0) return
1426       go to (13,14,15) ith
1427       print *, 'here'
1428       go to 5
1429    12 ith=1
1430       t=dlog(qsctt/al2)
1431       go to 11
1432    13 alfqc4=alpmsrt
1433       flav=5.
1434       ith=2
1435       go to 11
1436    14 alfqc5=alpmsrt
1437       ith=3
1438       t=tt
1439       go to 11
1440    15 alfqs5=alpmsrt
1441       alfinv=1./alfqs5+1./alfqc4-1./alfqc5
1442       alpmsrt=1./alfinv
1443       return
1444
1445   311 continue
1446       b0=11-2.*flav/3.
1447       x1=4.*pi/b0
1448       if(iord.eq.0) then
1449       alpmsrt=x1/t
1450       elseif(iord.eq.1) then
1451       b1=102.-38.*flav/3.
1452       x2=b1/b0**2
1453       as2=x1/t*(1.-x2*dlog(t)/t)
1454    35 as=as2
1455       f=-t+x1/as-x2*dlog(x1/as+x2)
1456       fp=-x1/as**2*(1.-x2/(x1/as+x2))
1457       as2=as-f/fp
1458       del=abs(f/fp/as)
1459       if((del-tol).gt.0.) go to 35
1460       elseif(iord.eq.2) then
1461       alpmsrt=qwikalf(t,2,flav)
1462       endif
1463       if(ith.eq.0) return
1464
1465       go to (313,314,315) ith
1466   312 ith=1
1467       t=dlog(qsdtt/al2)
1468       go to 311
1469   313 alfqc4=alpmsrt
1470       flav=3.
1471       ith=2
1472       go to 311
1473   314 alfqc3=alpmsrt
1474       ith=3
1475       t=tt
1476       go to 311
1477   315 alfqs3=alpmsrt
1478       alfinv=1./alfqs3+1./alfqc4-1./alfqc3
1479       alpmsrt=1./alfinv
1480       return
1481       end
1482 c***********************************

1483       double precision function qwikalf(t,iord,flav)
1484       implicit real*8(a-h,o-z)
1485       dimension z3(6),z4(6),z5(6),zz3(6),zz4(6),zz5(6)
1486       data z3/ -.161667e+01,0.954244e+01,
1487      .0.768623e+01,0.101523e+00,-.360127e-02,0.457867e-04/
1488       data z4/ -.172239e+01,0.831185e+01,
1489      .0.721463e+01,0.835531e-01,-.285436e-02,0.349129e-04/
1490       data z5/ -.872190e+00,0.572816e+01,
1491      .0.716119e+01,0.195884e-01,-.300199e-03,0.151741e-05/
1492       data zz3/-.155611e+02,0.168406e+02,
1493      .0.603014e+01,0.257682e+00,-.970217e-02,0.127628e-03/
1494       data zz4/-.106762e+02,0.118497e+02,0.664964e+01,
1495      .0.112996e+00,-.317551e-02,0.302434e-04/
1496       data zz5/-.531860e+01,0.708503e+01,0.698352e+01,
1497      .0.274170e-01,-.426894e-03,0.217591e-05/
1498       data pi/3.14159/
1499       nfm2=flav-2.
1500       x=dsqrt(t)
1501       x2=x*x
1502       x3=x*x2
1503       x4=x*x3
1504       x5=x*x4
1505       go to (1,2) iord  !!iord change!!
1506     1 go to (3,4,5) nfm2
1507     3 y=z3(1)+z3(2)*x+z3(3)*x2+z3(4)*x3+z3(5)*x4+z3(6)*x5
1508       go to 10
1509     4 y=z4(1)+z4(2)*x+z4(3)*x2+z4(4)*x3+z4(5)*x4+z4(6)*x5
1510       go to 10
1511     5 y=z5(1)+z5(2)*x+z5(3)*x2+z5(4)*x3+z5(5)*x4+z5(6)*x5
1512       go to 10
1513     2 go to (6,7,8) nfm2
1514     6 y=zz3(1)+zz3(2)*x+zz3(3)*x2+zz3(4)*x3+zz3(5)*x4+zz3(6)*x5
1515       go to 10
1516     7 y=zz4(1)+zz4(2)*x+zz4(3)*x2+zz4(4)*x3+zz4(5)*x4+zz4(6)*x5
1517       go to 10
1518     8 y=zz5(1)+zz5(2)*x+zz5(3)*x2+zz5(4)*x3+zz5(5)*x4+zz5(6)*x5
1519       go to 10
1520    10 qwikalf=4.*pi/y
1521       return
1522       end