LHAPDF/lhapdf5.2.2/wrapwhitg.f

   1       subroutine WHITGevolvep(xin,qin,p2in,ip2in,pdf)
   2       include 'parmsetup.inc'
   3       real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5
   4       real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu
   5       character*16 name(nmxset)
   6       integer nmem(nmxset),ndef(nmxset),mmem
   7       common/NAME/name,nmem,ndef,mmem
   8       integer nset
   9
  10       save
  11       call getnset(iset)
  12       call getnmem(iset,imem)
  13
  14       if(imem.eq.1.or.imem.eq.0) then
  15         call SFWHI1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  16
  17       elseif(imem.eq.2) then
  18         call SFWHI2(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  19
  20       elseif(imem.eq.3) then
  21         call SFWHI3(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  22
  23       elseif(imem.eq.4) then
  24         call SFWHI4(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  25
  26       elseif(imem.eq.5) then
  27         call SFWHI5(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  28
  29       elseif(imem.eq.6) then
  30         call SFWHI6(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  31
  32       else
  33         CONTINUE
  34       endif
  35
  36       pdf(-6)= 0.0d0
  37       pdf(6)= 0.0d0
  38       pdf(-5)= 0.0d0
  39       pdf(5 )= 0.0d0
  40       pdf(-4)= chm
  41       pdf(4 )= chm
  42       pdf(-3)= str
  43       pdf(3 )= str
  44       pdf(-2)= usea
  45       pdf(2 )= upv
  46       pdf(-1)= dsea
  47       pdf(1 )= dnv
  48       pdf(0 )= glu
  49
  50       return
  51 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  52       entry WHITGread(nset)
  53       read(1,*)nmem(nset),ndef(nset)
  54       return
  55 c
  56 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  57       entry WHITGalfa(alfas,qalfa)
  58         call getnset(iset)
  59         call getnmem(iset,imem)
  60         call GetOrderAsM(iset,iord)
  61         call Getlam4M(iset,imem,qcdl4)
  62         call Getlam5M(iset,imem,qcdl5)
  63         call aspdflib(alfas,Qalfa,iord,qcdl5)
  64
  65       return
  66 c
  67 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  68       entry WHITGinit(Eorder,Q2fit)
  69       return
  70 c
  71 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  72       entry WHITGpdf(mem)
  73       call getnset(iset)
  74       call setnmem(iset,mem)
  75 c      imem = mem
  76       return
  77 c
  78  1000 format(5e13.5)
  79       end
  80 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  81 c-------------------------------------------------------
  82       subroutine SFWHI1(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
  83 c-------------------------------------------------------
  84 c     WHIT1 parton distribution in the photon
  85 c
  86 c     INPUT:  integer ic  : if ic=0 then qc=0
  87 c                           else qc is calculated
  88 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
  89 c             DOUBLE PRECISION  x   : energy fraction
  90 c
  91 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
  92 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
  93 c             DOUBLE PRECISION  qc  : charm-quark dist.
  94 c             DOUBLE PRECISION  g   : gluon dist.
  95 c-------------------------------------------------------
  96 c     Modified by M.Tanaka on July 22, 1994.
  97 c     The bug pointed out by M.Drees is fixed.
  98 c-------------------------------------------------------
  99 c     Modified by I.Watanabe on July 22, 1994.
 100 c-------------------------------------------------------
 101       implicit none
 102       external WHIT1G
 103       double precision
 104      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
 105 c arg
 106       integer ic
 107       DOUBLE PRECISION Q2,x
 108       DOUBLE PRECISION qu,qd,qc,g
 109 c const
 110       DOUBLE PRECISION q42it,q52it,lam42,lam52
 111       DOUBLE PRECISION alinv,mc,PI
 112 c local
 113       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
 114       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
 115      $       A0sea,B0sea,BB0sea,C0sea
 116       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
 117       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
 118       DOUBLE PRECISION x1,x2,mc2q2
 119       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
 120       DOUBLE PRECISION WHIT1G
 121 c parameters
 122       parameter(lam42=0.16d0, lam52=0.091411319d0)
 123       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
 124       parameter(alinv=137.036d0, mc=1.5d0)
 125       parameter(pi=3.14159265358979323846d0)
 126       common /scale/ s,s2,s3,s4,prsccf
 127 c
 128 c begin
 129       x=ZX
 130       Q2=ZQ*ZQ
 131       ic=1
 132 c
 133       x1=1.0d0-x
 134       x2=x**2
 135       mc2q2=mc**2/Q2
 136 c
 137       if(Q2.lt.100.0d0) then
 138 c  under 100 GeV^2
 139 c
 140 c  set  scale s
 141          if(Q2.lt.4.0d0) then
 142 cccc  for under 4GeV^2 prescription
 143             s=  0.0d0
 144             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
 145             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
 146          else
 147             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
 148             prsccf = 1.0d0
 149             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
 150          endif
 151             s2=s**2
 152             s3=s2*s
 153             s4=s2**2
 154 c
 155 cccccc   WHIT1 quark (U100)
 156 c
 157       A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
 158       A1val=              s*(-2.361000d+00)+s2*(-1.136000d+00)
 159       A2val=              s*( 5.280000d-01)+s2*( 2.406000d+00)
 160       Bval = 5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
 161       Cval = 2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
 162      $           +s3*(-5.040000d-02)
 163       A0sea= 6.510000d-01+s*( 1.291000d+00)+s2*(-4.470000d+00)
 164      $           +s3*( 5.140000d+00)+s4*(-2.091000d+00)
 165       B0sea=-3.820000d-02+s*( 9.010000d-02)+s2*(-1.356000d+00)
 166      $           +s3*( 1.582000d+00)+s4*(-6.440000d-01)
 167       BB0sea=2.084000d+00+s*( 7.740000d+00)+s2*(-2.970000d+01)
 168      $           +s3*( 3.860000d+01)+s4*(-1.705000d+01)
 169       C0sea= 7.000000d+00+s*(-1.608000d+01)+s2*( 4.670000d+01)
 170      $           +s3*(-5.710000d+01)+s4*( 2.386000d+01)
 171 c
 172          qv  = prsccf/alinv/x*
 173      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 174          qsea= prsccf/alinv/x*
 175      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 176 c
 177          qu  = qv/3.0d0  + qsea/6.0d0
 178          qu  = qu*x
 179          ZUV=qu
 180          ZUB=qu
 181          qd  = qv/12.0d0 + qsea/6.0d0
 182          qd  = qd*x
 183          ZDV=qd
 184          ZDB=qd
 185          ZSB=qd
 186 c
 187          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 188             call WHIT1Q(x,mc2q2,cv,cs)
 189             qc = cv/alinv/2.0d0/PI + cs*alstpi
 190             qc  = qc*x
 191             ZCB=qc
 192          else
 193             qc = 0.0d0
 194             ZCB=qc
 195          endif
 196 c
 197          g   = WHIT1G(x,Q2)
 198          g   = g*x
 199          ZGL=g
 200 c
 201       else
 202 c over 100 GeV^2
 203 c
 204 c set scale s
 205          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
 206          prsccf = 1.0d0
 207          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
 208             s2=s**2
 209             s3=s2*s
 210             s4=s2**2
 211 c
 212 cccccc   WHIT1 quark (O100)
 213 c
 214       A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
 215       A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.251000d-01)
 216       A2val= 1.522000d+00+s*( 4.310000d+00)+s2*( 1.314000d+00)
 217       Bval = 5.170000d-01+s*( 4.040000d-02)+s2*(-2.100000d-02)
 218       Cval = 1.655000d-01+s*(-2.062000d-02)+s2*( 5.360000d-02)
 219       A0sea= 6.250000d-01+s*(-5.890000d-01)+s2*( 4.180000d+00)
 220      $           +s3*(-1.206000d+01)+s4*( 1.257000d+01)
 221       B0sea=-2.492000d-01+s*(-4.110000d-01)+s2*( 9.660000d-01)
 222      $           +s3*(-2.584000d+00)+s4*( 2.670000d+00)
 223       BB0sea=2.100000d+00+s*(-5.750000d+00)+s2*( 4.780000d+01)
 224      $           +s3*(-1.407000d+02)+s4*( 1.476000d+02)
 225       C0sea= 4.780000d+00+s*( 4.860000d+00)+s2*(-4.890000d+01)
 226      $           +s3*( 1.477000d+02)+s4*(-1.602000d+02)
 227 c
 228          qv  = 1.0d0/alinv/x*
 229      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 230          qsea= 1.0d0/alinv/x*
 231      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 232 c
 233          qu  = qv/3.0d0  + qsea/6.0d0
 234          qu  = qu*x
 235          ZUV=qu
 236          ZUB=qu
 237          qd  = qv/12.0d0 + qsea/6.0d0
 238          qd  = qd*x
 239          ZDV=qd
 240          ZDB=qd
 241          ZSB=qd
 242          g   = WHIT1G(x,Q2)
 243          g   = g*x
 244          ZGL=g
 245 c
 246          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 247       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
 248      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
 249       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
 250      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
 251       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
 252      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
 253       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
 254      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
 255       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
 256      $           +s3*( 1.967000d+01)
 257       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
 258      $           +s3*( 4.630000d+01)
 259       Adcs =              s*(-1.815000d-02)+s2*( 2.043000d-03)
 260      $           +s3*(-4.130000d-03)
 261       B0dcs=-3.086000d-01+s*(-2.565000d-01)+s2*( 9.840000d-02)
 262       B1dcs= 1.376000d+00+s*(-4.630000d-01)+s2*( 1.232000d+00)
 263       Cdcs = 3.650000d+00+s*( 7.290000d-01)+s2*(-7.570000d+00)
 264      $           +s3*( 7.790000d+00)
 265 c
 266          dcv = 1.0d0/alinv/x*
 267      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
 268          dcs = 1.0d0/alinv/x*
 269      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
 270 c
 271            call WHIT1Q(x,mc*mc/Q2,cv,cs)
 272            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
 273            qc  = qc*x
 274            ZCB=qc
 275          else
 276            qc = 0.0d0
 277            ZCB=qc
 278          endif
 279       endif
 280 c
 281       return
 282       end
 283 c-------------------------------------------------------
 284       subroutine SFWHI2(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
 285 c-------------------------------------------------------
 286 c     WHIT2 parton distribution in the photon
 287 c
 288 c     INPUT:  integer ic  : if ic=0 then qc=0
 289 c                           else qc is calculated
 290 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
 291 c             DOUBLE PRECISION  x   : energy fraction
 292 c
 293 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
 294 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
 295 c             DOUBLE PRECISION  qc  : charm-quark dist.
 296 c             DOUBLE PRECISION  g   : gluon dist.
 297 c-------------------------------------------------------
 298 c     Modified by M.Tanaka on July 22, 1994.
 299 c     The bug pointed out by M.Drees is fixed.
 300 c-------------------------------------------------------
 301 c     Modified by I.Watanabe on July 22, 1994.
 302 c-------------------------------------------------------
 303       implicit none
 304       external WHIT2G
 305       double precision
 306      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
 307 c arg
 308       integer ic
 309       DOUBLE PRECISION Q2,x
 310       DOUBLE PRECISION qu,qd,qc,g
 311 c const
 312       DOUBLE PRECISION q42it,q52it,lam42,lam52
 313       DOUBLE PRECISION alinv,mc,PI
 314 c local
 315       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
 316       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
 317      $       A0sea,B0sea,BB0sea,C0sea
 318       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
 319       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
 320       DOUBLE PRECISION x1,x2,mc2q2
 321       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
 322       DOUBLE PRECISION WHIT2G
 323 c parameters
 324       parameter(lam42=0.16d0, lam52=0.091411319d0)
 325       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
 326       parameter(alinv=137.036d0, mc=1.5d0)
 327       parameter(pi=3.14159265358979323846d0)
 328       common /scale/ s,s2,s3,s4,prsccf
 329 c
 330 c begin
 331       x=ZX
 332       Q2=ZQ*ZQ
 333       ic=1
 334 c
 335       x1=1.0d0-x
 336       x2=x**2
 337       mc2q2=mc**2/Q2
 338 c
 339       if(Q2.lt.100.0d0) then
 340 c  under 100 GeV^2
 341 c
 342 c  set  scale s
 343          if(Q2.lt.4.0d0) then
 344 cccc  for under 4GeV^2 prescription
 345             s=  0.0d0
 346             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
 347             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
 348          else
 349             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
 350             prsccf = 1.0d0
 351             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
 352          endif
 353             s2=s**2
 354             s3=s2*s
 355             s4=s2**2
 356 c
 357 cccccc   WHIT2 quark (U100)
 358 c
 359       A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
 360       A1val=              s*(-2.361000d+00)+s2*(-1.136000d+00)
 361       A2val=              s*( 5.280000d-01)+s2*( 2.406000d+00)
 362       Bval=  5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
 363       Cval=  2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
 364      $           +s3*(-5.040000d-02)
 365       A0sea= 1.237000d+00+s*( 3.390000d+00)+s2*(-1.075000d+01)
 366      $           +s3*( 1.246000d+01)+s4*(-5.580000d+00)
 367       B0sea=-7.270000d-02+s*( 1.748000d-01)+s2*(-1.392000d+00)
 368      $           +s3*( 1.711000d+00)+s4*(-7.960000d-01)
 369       BB0sea=4.290000d+00+s*( 1.787000d+01)+s2*(-5.810000d+01)
 370      $           +s3*( 8.190000d+01)+s4*(-4.140000d+01)
 371       C0sea= 1.434000d+01+s*(-4.490000d+01)+s2*( 1.197000d+02)
 372      $           +s3*(-1.585000d+02)+s4*( 7.530000d+01)
 373 c
 374          qv  = prsccf/alinv/x*
 375      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 376          qsea= prsccf/alinv/x*
 377      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 378 c
 379          qu  = qv/3.0d0  + qsea/6.0d0
 380          qu  = qu*x
 381          ZUV=qu
 382          ZUB=qu
 383          qd  = qv/12.0d0 + qsea/6.0d0
 384          qd  = qd*x
 385          ZDV=qd
 386          ZDB=qd
 387          ZSB=qd
 388 c
 389          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 390             call WHIT2Q(x,mc2q2,cv,cs)
 391             qc = cv/alinv/2.0d0/PI + cs*alstpi
 392             qc  = qc*x
 393             ZCB=qc
 394          else
 395             qc = 0.0d0
 396             ZCB=qc
 397          endif
 398 c
 399          g   = WHIT2G(x,Q2)
 400          g   = g*x
 401          ZGL=g
 402 c
 403       else
 404 c over 100 GeV^2
 405 c
 406 c set scale s
 407          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
 408          prsccf = 1.0d0
 409          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
 410             s2=s**2
 411             s3=s2*s
 412             s4=s2**2
 413 c
 414 cccccc   WHIT2 quark (O100)
 415 c
 416       A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
 417       A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.259000d-01)
 418       A2val= 1.522000d+00+s*( 4.300000d+00)+s2*( 1.315000d+00)
 419       Bval = 5.170000d-01+s*( 4.030000d-02)+s2*(-2.098000d-02)
 420       Cval = 1.655000d-01+s*(-2.063000d-02)+s2*( 5.370000d-02)
 421       A0sea= 1.287000d+00+s*(-2.069000d+00)+s2*( 1.157000d+01)
 422      $           +s3*(-3.570000d+01)+s4*( 3.740000d+01)
 423       B0sea=-2.340000d-01+s*(-4.430000d-01)+s2*( 1.235000d+00)
 424      $           +s3*(-3.720000d+00)+s4*( 3.840000d+00)
 425       BB0sea=6.460000d+00+s*(-1.048000d+01)+s2*( 8.980000d+01)
 426      $           +s3*(-2.847000d+02)+s4*( 2.998000d+02)
 427       C0sea= 5.350000d+00+s*( 1.011000d+01)+s2*(-1.337000d+02)
 428      $           +s3*( 4.270000d+02)+s4*(-4.570000d+02)
 429 c
 430          qv  = 1.0d0/alinv/x*
 431      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 432          qsea= 1.0d0/alinv/x*
 433      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 434 c
 435          qu  = qv/3.0d0  + qsea/6.0d0
 436          qu  = qu*x
 437          ZUV=qu
 438          ZUB=qu
 439          qd  = qv/12.0d0 + qsea/6.0d0
 440          qd  = qd*x
 441          ZDV=qd
 442          ZDB=qd
 443          ZSB=qd
 444          g   = WHIT2G(x,Q2)
 445          g   = g*x
 446          ZGL=g
 447 c
 448          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 449       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
 450      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
 451       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
 452      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
 453       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
 454      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
 455       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
 456      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
 457       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
 458      $           +s3*( 1.967000d+01)
 459       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
 460      $           +s3*( 4.630000d+01)
 461       Adcs =              s*(-2.786000d-02)+s2*( 3.490000d-02)
 462      $           +s3*(-2.223000d-02)
 463       B0dcs=-3.141000d-01+s*(-4.250000d-01)+s2*( 1.564000d-01)
 464       B1dcs= 4.720000d+00+s*(-5.480000d+00)+s2*( 2.686000d+00)
 465       Cdcs = 2.961000d+00+s*( 7.760000d-01)+s2*(-8.280000d+00)
 466      $           +s3*( 9.780000d+00)
 467 c
 468          dcv = 1.0d0/alinv/x*
 469      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
 470          dcs = 1.0d0/alinv/x*
 471      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
 472 c
 473            call WHIT2Q(x,mc*mc/Q2,cv,cs)
 474            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
 475            qc  = qc*x
 476            ZCB=qc
 477          else
 478            qc = 0.0d0
 479            ZCB=qc
 480          endif
 481       endif
 482 c
 483       return
 484       end
 485 c-------------------------------------------------------
 486       subroutine SFWHI3(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
 487 c-------------------------------------------------------
 488 c     WHIT3 parton distribution in the photon
 489 c
 490 c     INPUT:  integer ic  : if ic=0 then qc=0
 491 c                           else qc is calculated
 492 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
 493 c             DOUBLE PRECISION  x   : energy fraction
 494 c
 495 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
 496 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
 497 c             DOUBLE PRECISION  qc  : charm-quark dist.
 498 c             DOUBLE PRECISION  g   : gluon dist.
 499 c-------------------------------------------------------
 500 c     Modified by M.Tanaka on July 22, 1994.
 501 c     The bug pointed out by M.Drees is fixed.
 502 c-------------------------------------------------------
 503 c     Modified by I.Watanabe on July 22, 1994.
 504 c-------------------------------------------------------
 505       implicit none
 506       external whit3g
 507       double precision
 508      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
 509 c arg
 510       integer ic
 511       DOUBLE PRECISION Q2,x
 512       DOUBLE PRECISION qu,qd,qc,g
 513 c const
 514       DOUBLE PRECISION q42it,q52it,lam42,lam52
 515       DOUBLE PRECISION alinv,mc,PI
 516 c local
 517       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
 518       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
 519      $       A0sea,B0sea,BB0sea,C0sea
 520       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
 521       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
 522       DOUBLE PRECISION x1,x2,mc2q2
 523       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
 524       DOUBLE PRECISION WHIT3g
 525 c parameters
 526       parameter(lam42=0.16d0, lam52=0.091411319d0)
 527       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
 528       parameter(alinv=137.036d0, mc=1.5d0)
 529       parameter(pi=3.14159265358979323846d0)
 530       common /scale/ s,s2,s3,s4,prsccf
 531 c
 532 c begin
 533       x=ZX
 534       Q2=ZQ*ZQ
 535       ic=1
 536 c
 537       x1=1.0d0-x
 538       x2=x**2
 539       mc2q2=mc**2/Q2
 540 c
 541       if(Q2.lt.100.0d0) then
 542 c  under 100 GeV^2
 543 c
 544 c  set  scale s
 545          if(Q2.lt.4.0d0) then
 546 cccc  for under 4GeV^2 prescription
 547             s=  0.0d0
 548             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
 549             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
 550          else
 551             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
 552             prsccf = 1.0d0
 553             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
 554          endif
 555             s2=s**2
 556             s3=s2*s
 557             s4=s2**2
 558 c
 559 cccccc   WHIT3 quark (U100)
 560 c
 561       A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
 562       A1val=              s*(-2.361000d+00)+s2*(-1.136000d+00)
 563       A2val=              s*( 5.280000d-01)+s2*( 2.406000d+00)
 564       Bval = 5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
 565       Cval = 2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
 566      $           +s3*(-5.040000d-02)
 567       A0sea= 1.587000d+00+s*( 5.050000d+00)+s2*(-1.126000d+01)
 568      $           +s3*( 7.560000d+00)+s4*(-1.471000d+00)
 569       B0sea=-1.006000d-01+s*( 2.259000d-01)+s2*(-1.195000d+00)
 570      $           +s3*( 1.175000d+00)+s4*(-4.460000d-01)
 571       BB0sea=5.730000d+00+s*( 2.564000d+01)+s2*(-5.870000d+01)
 572      $           +s3*( 6.320000d+01)+s4*(-2.577000d+01)
 573       C0sea= 2.136000d+01+s*(-7.290000d+01)+s2*( 1.532000d+02)
 574      $           +s3*(-1.679000d+02)+s4*( 6.740000d+01)
 575 c
 576          qv  = prsccf/alinv/x*
 577      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 578          qsea= prsccf/alinv/x*
 579      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 580 c
 581          qu  = qv/3.0d0  + qsea/6.0d0
 582          qu  = qu*x
 583          ZUV=qu
 584          ZUB=qu
 585          qd  = qv/12.0d0 + qsea/6.0d0
 586          qd  = qd*x
 587          ZDV=qd
 588          ZDB=qd
 589          ZSB=qd
 590 c
 591          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 592             call WHIT3Q(x,mc2q2,cv,cs)
 593             qc = cv/alinv/2.0d0/PI + cs*alstpi
 594             qc  = qc*x
 595             ZCB=qc
 596          else
 597             qc = 0.0d0
 598             ZCB=qc
 599          endif
 600 c
 601          g   = WHIT3G(x,Q2)
 602          g   = g*x
 603          ZGL=g
 604 c
 605       else
 606 c over 100 GeV^2
 607 c
 608 c set scale s
 609          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
 610          prsccf = 1.0d0
 611          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
 612             s2=s**2
 613             s3=s2*s
 614             s4=s2**2
 615 c
 616 cccccc   WHIT3 quark (O100)
 617 c
 618       A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
 619       A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.264000d-01)
 620       A2val= 1.522000d+00+s*( 4.300000d+00)+s2*( 1.315000d+00)
 621       Bval = 5.170000d-01+s*( 4.030000d-02)+s2*(-2.097000d-02)
 622       Cval = 1.655000d-01+s*(-2.064000d-02)+s2*( 5.370000d-02)
 623       A0sea= 1.850000d+00+s*(-3.670000d+00)+s2*( 2.714000d+01)
 624      $           +s3*(-1.066000d+02)+s4*( 1.309000d+02)
 625       B0sea=-2.299000d-01+s*(-4.970000d-01)+s2*( 2.464000d+00)
 626      $           +s3*(-9.950000d+00)+s4*( 1.232000d+01)
 627       BB0sea=1.042000d+01+s*(-1.074000d+01)+s2*( 1.327000d+02)
 628      $           +s3*(-5.390000d+02)+s4*( 6.560000d+02)
 629       C0sea= 4.070000d+00+s*( 4.110000d+00)+s2*(-1.719000d+02)
 630      $           +s3*( 7.070000d+02)+s4*(-8.590000d+02)
 631 c
 632          qv  = 1.0d0/alinv/x*
 633      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 634          qsea= 1.0d0/alinv/x*
 635      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 636 c
 637          qu  = qv/3.0d0  + qsea/6.0d0
 638          qu  = qu*x
 639          ZUV=qu
 640          ZUB=qu
 641          qd  = qv/12.0d0 + qsea/6.0d0
 642          qd  = qd*x
 643          ZDV=qd
 644          ZDB=qd
 645          ZSB=qd
 646          g   = WHIT3G(x,Q2)
 647          g   = g*x
 648          ZGL=g
 649 c
 650          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 651       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
 652      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
 653       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
 654      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
 655       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
 656      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
 657       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
 658      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
 659       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
 660      $           +s3*( 1.967000d+01)
 661       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
 662      $           +s3*( 4.630000d+01)
 663       Adcs =              s*(-1.948000d-02)+s2*( 2.861000d-02)
 664      $           +s3*(-2.036000d-02)
 665       B0dcs=-4.130000d-01+s*(-4.390000d-01)+s2*( 1.810000d-01)
 666       B1dcs= 5.190000d+00+s*(-7.400000d+00)+s2*( 3.400000d+00)
 667       Cdcs = 2.359000d+00+s*( 9.770000d-01)+s2*(-7.730000d+00)
 668      $           +s3*( 9.480000d+00)
 669 c
 670          dcv = 1.0d0/alinv/x*
 671      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
 672          dcs = 1.0d0/alinv/x*
 673      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
 674 c
 675            call WHIT3Q(x,mc*mc/Q2,cv,cs)
 676            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
 677            qc  = qc*x
 678            ZCB=qc
 679          else
 680            qc = 0.0d0
 681            ZCB=qc
 682          endif
 683       endif
 684 c
 685       return
 686       end
 687 c-------------------------------------------------------
 688       subroutine SFWHI4(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
 689 c-------------------------------------------------------
 690 c     WHIT4 parton distribution in the photon
 691 c
 692 c     INPUT:  integer ic  : if ic=0 then qc=0
 693 c                           else qc is calculated
 694 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
 695 c             DOUBLE PRECISION  x   : energy fraction
 696 c
 697 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
 698 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
 699 c             DOUBLE PRECISION  qc  : charm-quark dist.
 700 c             DOUBLE PRECISION  g   : gluon dist.
 701 c-------------------------------------------------------
 702 c     Modified by M.Tanaka on July 22, 1994.
 703 c     The bug pointed out by M.Drees is fixed.
 704 c-------------------------------------------------------
 705 c     Modified by I.Watanabe on July 22, 1994.
 706 c-------------------------------------------------------
 707       implicit none
 708       external WHIT4G
 709       double precision
 710      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
 711 c arg
 712       integer ic
 713       DOUBLE PRECISION Q2,x
 714       DOUBLE PRECISION qu,qd,qc,g
 715 c const
 716       DOUBLE PRECISION q42it,q52it,lam42,lam52
 717       DOUBLE PRECISION alinv,mc,PI
 718 c local
 719       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
 720       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
 721      $       A0sea,B0sea,BB0sea,C0sea
 722       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
 723       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
 724       DOUBLE PRECISION x1,x2,mc2q2
 725       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
 726       DOUBLE PRECISION WHIT4G
 727 c parameters
 728       parameter(lam42=0.16d0, lam52=0.091411319d0)
 729       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
 730       parameter(alinv=137.036d0, mc=1.5d0)
 731       parameter(pi=3.14159265358979323846d0)
 732       common /scale/ s,s2,s3,s4,prsccf
 733 c
 734 c begin
 735       x=ZX
 736       Q2=ZQ*ZQ
 737       ic=1
 738 c
 739       x1=1.0d0-x
 740       x2=x**2
 741       mc2q2=mc**2/Q2
 742 c
 743       if(Q2.lt.100.0d0) then
 744 c  under 100 GeV^2
 745 c
 746 c  set  scale s
 747          if(Q2.lt.4.0d0) then
 748 cccc  for under 4GeV^2 prescription
 749             s=  0.0d0
 750             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
 751             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
 752          else
 753             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
 754             prsccf = 1.0d0
 755             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
 756          endif
 757             s2=s**2
 758             s3=s2*s
 759             s4=s2**2
 760 c
 761 cccccc   WHIT4 quark (U100)
 762 c
 763       A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
 764       A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
 765       A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
 766      $           +s3*(-5.200000d+00)
 767       Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
 768       Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
 769      $           +s3*(-9.950000d-01)
 770       A0sea= 1.308000d+00+s*( 2.315000d+00)+s2*(-7.880000d+00)
 771      $           +s3*( 8.260000d+00)+s4*(-3.004000d+00)
 772       B0sea=-3.730000d-02+s*( 5.630000d-02)+s2*(-1.133000d+00)
 773      $           +s3*( 1.185000d+00)+s4*(-4.180000d-01)
 774       BB0sea=2.103000d+00+s*( 4.850000d+00)+s2*(-1.781000d+01)
 775      $           +s3*( 2.062000d+01)+s4*(-7.940000d+00)
 776       C0sea= 7.000000d+00+s*(-1.017000d+01)+s2*( 2.600000d+01)
 777      $           +s3*(-2.960000d+01)+s4*( 1.227000d+01)
 778 c
 779          qv  = prsccf/alinv/x*
 780      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 781          qsea= prsccf/alinv/x*
 782      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 783 c
 784          qu  = qv/3.0d0  + qsea/6.0d0
 785          qu  = qu*x
 786          ZUV=qu
 787          ZUB=qu
 788          qd  = qv/12.0d0 + qsea/6.0d0
 789          qd  = qd*x
 790          ZDV=qd
 791          ZDB=qd
 792          ZSB=qd
 793 c
 794          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 795             call WHIT4Q(x,mc2q2,cv,cs)
 796             qc = cv/alinv/2.0d0/PI + cs*alstpi
 797             qc  = qc*x
 798             ZCB=qc
 799          else
 800             qc = 0.0d0
 801             ZCB=qc
 802          endif
 803 c
 804          g   = WHIT4G(x,Q2)
 805          g   = g*x
 806          ZGL=g
 807 c
 808       else
 809 c over 100 GeV^2
 810 c
 811 c set scale s
 812          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
 813          prsccf = 1.0d0
 814          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
 815             s2=s**2
 816             s3=s2*s
 817             s4=s2**2
 818 c
 819 cccccc   WHIT4 quark (O100)
 820 c
 821       A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.619000d+00)
 822       A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.430000d+00)
 823       A2val= 2.837000d+00+s*( 6.470000d+00)+s2*( 4.090000d+00)
 824       Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.756000d-02)
 825       Cval = 1.728000d-01+s*(-2.479000d-02)+s2*( 1.446000d-01)
 826       A0sea= 1.188000d+00+s*(-1.396000d+00)+s2*( 8.710000d+00)
 827      $           +s3*(-2.542000d+01)+s4*( 2.492000d+01)
 828       B0sea=-2.448000d-01+s*(-4.190000d-01)+s2*( 1.007000d+00)
 829      $           +s3*(-2.689000d+00)+s4*( 2.517000d+00)
 830       BB0sea=1.942000d+00+s*(-6.040000d+00)+s2*( 5.030000d+01)
 831      $           +s3*(-1.478000d+02)+s4*( 1.481000d+02)
 832       C0sea= 5.420000d+00+s*( 6.110000d+00)+s2*(-5.380000d+01)
 833      $           +s3*( 1.632000d+02)+s4*(-1.716000d+02)
 834 c
 835          qv  = 1.0d0/alinv/x*
 836      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 837          qsea= 1.0d0/alinv/x*
 838      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 839 c
 840          qu  = qv/3.0d0  + qsea/6.0d0
 841          qu  = qu*x
 842          ZUV=qu
 843          ZUB=qu
 844          qd  = qv/12.0d0 + qsea/6.0d0
 845          qd  = qd*x
 846          ZDV=qd
 847          ZDB=qd
 848          ZSB=qd
 849          g   = WHIT4G(x,Q2)
 850          g   = g*x
 851          ZGL=g
 852 c
 853          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 854       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
 855      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
 856       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
 857      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
 858       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
 859      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
 860       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
 861      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
 862       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
 863      $           +s3*( 1.967000d+01)
 864       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
 865      $           +s3*( 4.630000d+01)
 866       Adcs =              s*(-2.821000d-02)+s2*(-2.649000d-04)
 867      $           +s3*( 7.040000d-03)
 868       B0dcs=-3.270000d-01+s*(-2.298000d-01)+s2*( 3.500000d-02)
 869       B1dcs= 1.254000d+00+s*( 8.780000d-01)+s2*( 2.086000d-01)
 870       Cdcs = 4.170000d+00+s*( 6.400000d-01)+s2*(-7.630000d+00)
 871      $           +s3*( 7.170000d+00)
 872 c
 873          dcv = 1.0d0/alinv/x*
 874      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
 875          dcs = 1.0d0/alinv/x*
 876      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
 877 c
 878            call WHIT4Q(x,mc*mc/Q2,cv,cs)
 879            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
 880            qc  = qc*x
 881            ZCB=qc
 882          else
 883            qc = 0.0d0
 884            ZCB=qc
 885          endif
 886       endif
 887 c
 888       return
 889       end
 890 c-------------------------------------------------------
 891       subroutine SFWHI5(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
 892 c-------------------------------------------------------
 893 c     WHIT5 parton distribution in the photon
 894 c
 895 c     INPUT:  integer ic  : if ic=0 then qc=0
 896 c                           else qc is calculated
 897 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
 898 c             DOUBLE PRECISION  x   : energy fraction
 899 c
 900 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
 901 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
 902 c             DOUBLE PRECISION  qc  : charm-quark dist.
 903 c             DOUBLE PRECISION  g   : gluon dist.
 904 c-------------------------------------------------------
 905 c     Modified by M.Tanaka on July 22, 1994.
 906 c     The bug pointed out by M.Drees is fixed.
 907 c-------------------------------------------------------
 908 c     Modified by I.Watanabe on July 22, 1994.
 909 c-------------------------------------------------------
 910       implicit none
 911       external WHIT5G
 912       double precision
 913      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
 914 c arg
 915       integer ic
 916       DOUBLE PRECISION Q2,x
 917       DOUBLE PRECISION qu,qd,qc,g
 918 c const
 919       DOUBLE PRECISION q42it,q52it,lam42,lam52
 920       DOUBLE PRECISION alinv,mc,PI
 921 c local
 922       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
 923       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
 924      $       A0sea,B0sea,BB0sea,C0sea
 925       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
 926       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
 927       DOUBLE PRECISION x1,x2,mc2q2
 928       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
 929       DOUBLE PRECISION WHIT5G
 930 c parameters
 931       parameter(lam42=0.16d0, lam52=0.091411319d0)
 932       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
 933       parameter(alinv=137.036d0, mc=1.5d0)
 934       parameter(pi=3.14159265358979323846d0)
 935       common /scale/ s,s2,s3,s4,prsccf
 936 c
 937 c begin
 938       x=ZX
 939       Q2=ZQ*ZQ
 940       ic=1
 941 c
 942       x1=1.0d0-x
 943       x2=x**2
 944       mc2q2=mc**2/Q2
 945 c
 946       if(Q2.lt.100.0d0) then
 947 c  under 100 GeV^2
 948 c
 949 c  set  scale s
 950          if(Q2.lt.4.0d0) then
 951 cccc  for under 4GeV^2 prescription
 952             s=  0.0d0
 953             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
 954             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
 955          else
 956             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
 957             prsccf = 1.0d0
 958             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
 959          endif
 960             s2=s**2
 961             s3=s2*s
 962             s4=s2**2
 963 c
 964 cccccc   WHIT5 quark (U100)
 965 c
 966       A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
 967       A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
 968       A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
 969      $           +s3*(-5.200000d+00)
 970       Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
 971       Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
 972      $           +s3*(-9.950000d-01)
 973       A0sea= 2.227000d+00+s*( 5.720000d+00)+s2*(-1.295000d+01)
 974      $           +s3*( 7.220000d+00)+s4*(-2.514000d-01)
 975       B0sea=-8.810000d-02+s*( 1.465000d-01)+s2*(-9.750000d-01)
 976      $           +s3*( 7.820000d-01)+s4*(-2.074000d-01)
 977       BB0sea=3.370000d+00+s*( 1.416000d+01)+s2*(-3.150000d+01)
 978      $           +s3*( 2.789000d+01)+s4*(-8.710000d+00)
 979       C0sea= 1.581000d+01+s*(-3.630000d+01)+s2*( 7.710000d+01)
 980      $           +s3*(-7.810000d+01)+s4*( 2.948000d+01)
 981 c
 982          qv  = prsccf/alinv/x*
 983      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
 984          qsea= prsccf/alinv/x*
 985      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
 986 c
 987          qu  = qv/3.0d0  + qsea/6.0d0
 988          qu  = qu*x
 989          ZUV=qu
 990          ZUB=qu
 991          qd  = qv/12.0d0 + qsea/6.0d0
 992          qd  = qd*x
 993          ZDV=qd
 994          ZDB=qd
 995          ZSB=qd
 996 c
 997          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
 998             call WHIT5Q(x,mc2q2,cv,cs)
 999             qc = cv/alinv/2.0d0/PI + cs*alstpi
1000             qc  = qc*x
1001             ZCB=qc
1002          else
1003             qc = 0.0d0
1004             ZCB=qc
1005          endif
1006 c
1007          g   = WHIT5G(x,Q2)
1008          g   = g*x
1009          ZGL=g
1010 c
1011       else
1012 c over 100 GeV^2
1013 c
1014 c set scale s
1015          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
1016          prsccf = 1.0d0
1017          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
1018             s2=s**2
1019             s3=s2*s
1020             s4=s2**2
1021 c
1022 cccccc   WHIT5 quark (O100)
1023 c
1024       A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.617000d+00)
1025       A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.417000d+00)
1026       A2val= 2.837000d+00+s*( 6.470000d+00)+s2*( 4.070000d+00)
1027       Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.750000d-02)
1028       Cval = 1.728000d-01+s*(-2.457000d-02)+s2*( 1.440000d-01)
1029       A0sea= 2.318000d+00+s*(-3.760000d+00)+s2*( 2.026000d+01)
1030      $           +s3*(-5.950000d+01)+s4*( 5.900000d+01)
1031       B0sea=-2.425000d-01+s*(-4.360000d-01)+s2*( 1.241000d+00)
1032      $           +s3*(-3.510000d+00)+s4*( 3.360000d+00)
1033       BB0sea=5.330000d+00+s*(-8.680000d+00)+s2*( 7.420000d+01)
1034      $           +s3*(-2.070000d+02)+s4*( 1.967000d+02)
1035       C0sea= 8.480000d+00+s*( 9.310000d+00)+s2*(-1.041000d+02)
1036      $           +s3*( 2.801000d+02)+s4*(-2.663000d+02)
1037 c
1038          qv  = 1.0d0/alinv/x*
1039      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1040          qsea= 1.0d0/alinv/x*
1041      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1042 c
1043          qu  = qv/3.0d0  + qsea/6.0d0
1044          qu  = qu*x
1045          ZUV=qu
1046          ZUB=qu
1047          qd  = qv/12.0d0 + qsea/6.0d0
1048          qd  = qd*x
1049          ZDV=qd
1050          ZDB=qd
1051          ZSB=qd
1052          g   = WHIT5G(x,Q2)
1053          g   = g*x
1054          ZGL=g
1055 c
1056          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1057       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
1058      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
1059       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
1060      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
1061       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
1062      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
1063       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
1064      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
1065       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
1066      $           +s3*( 1.967000d+01)
1067       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
1068      $           +s3*( 4.630000d+01)
1069       Adcs =              s*(-6.580000d-02)+s2*( 1.059000d-01)
1070      $           +s3*(-6.630000d-02)
1071       B0dcs=-2.750000d-01+s*(-4.760000d-01)+s2*( 1.191000d-01)
1072       B1dcs= 6.370000d+00+s*(-5.320000d+00)+s2*( 1.986000d+00)
1073       Cdcs = 3.400000d+00+s*( 3.750000d-01)+s2*(-8.790000d+00)
1074      $           +s3*( 1.001000d+01)
1075 c
1076          dcv = 1.0d0/alinv/x*
1077      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
1078          dcs = 1.0d0/alinv/x*
1079      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
1080 c
1081            call WHIT5Q(x,mc*mc/Q2,cv,cs)
1082            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
1083            qc  = qc*x
1084            ZCB=qc
1085          else
1086            qc = 0.0d0
1087            ZCB=qc
1088          endif
1089       endif
1090 c
1091       return
1092       end
1093 c-------------------------------------------------------
1094       subroutine SFWHI6(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
1095 c-------------------------------------------------------
1096 c     WHIT6 parton distribution in the photon
1097 c
1098 c     INPUT:  integer ic  : if ic=0 then qc=0
1099 c                           else qc is calculated
1100 c             DOUBLE PRECISION  Q2  : energy scale Q^2 (GeV^2)
1101 c             DOUBLE PRECISION  x   : energy fraction
1102 c
1103 c     OUTPUT: DOUBLE PRECISION  qu  : up-quark dist.
1104 c             DOUBLE PRECISION  qd  : down- or strange-quark dist.
1105 c             DOUBLE PRECISION  qc  : charm-quark dist.
1106 c             DOUBLE PRECISION  g   : gluon dist.
1107 c-------------------------------------------------------
1108 c     Modified by M.Tanaka on July 22, 1994.
1109 c     The bug pointed out by M.Drees is fixed.
1110 c-------------------------------------------------------
1111 c     Modified by I.Watanabe on July 22, 1994.
1112 c-------------------------------------------------------
1113       implicit none
1114       external WHIT6G
1115       double precision
1116      +       ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
1117 c arg
1118       integer ic
1119       DOUBLE PRECISION Q2,x
1120       DOUBLE PRECISION qu,qd,qc,g
1121 c const
1122       DOUBLE PRECISION q42it,q52it,lam42,lam52
1123       DOUBLE PRECISION alinv,mc,PI
1124 c local
1125       DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
1126       DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
1127      $       A0sea,B0sea,BB0sea,C0sea
1128       DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
1129       DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
1130       DOUBLE PRECISION x1,x2,mc2q2
1131       DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
1132       DOUBLE PRECISION WHIT6G
1133 c parameters
1134       parameter(lam42=0.16d0, lam52=0.091411319d0)
1135       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1136       parameter(alinv=137.036d0, mc=1.5d0)
1137       parameter(pi=3.14159265358979323846d0)
1138       common /scale/ s,s2,s3,s4,prsccf
1139 c
1140 c begin
1141       x=ZX
1142       Q2=ZQ*ZQ
1143       ic=1
1144 c
1145       x1=1.0d0-x
1146       x2=x**2
1147       mc2q2=mc**2/Q2
1148 c
1149       if(Q2.lt.100.0d0) then
1150 c  under 100 GeV^2
1151 c
1152 c  set  scale s
1153          if(Q2.lt.4.0d0) then
1154 cccc  for under 4GeV^2 prescription
1155             s=  0.0d0
1156             prsccf =  log(Q2/LAM42)/ log(Q42IT/LAM42)
1157             alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
1158          else
1159             s=   log(  log(Q2/LAM42)/ log(Q42IT/LAM42))
1160             prsccf = 1.0d0
1161             alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
1162          endif
1163             s2=s**2
1164             s3=s2*s
1165             s4=s2**2
1166 c
1167 cccccc   WHIT6 quark (U100)
1168 c
1169       A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
1170       A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
1171       A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
1172      $           +s3*(-5.200000d+00)
1173       Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
1174       Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
1175      $           +s3*(-9.950000d-01)
1176       A0sea= 3.180000d+00+s*( 8.690000d+00)+s2*(-2.287000d+01)
1177      $           +s3*( 1.896000d+01)+s4*(-5.140000d+00)
1178       B0sea=-1.003000d-01+s*( 1.603000d-01)+s2*(-1.037000d+00)
1179      $           +s3*( 9.440000d-01)+s4*(-2.915000d-01)
1180       BB0sea=5.690000d+00+s*( 1.867000d+01)+s2*(-4.670000d+01)
1181      $           +s3*( 5.050000d+01)+s4*(-1.835000d+01)
1182       C0sea= 2.149000d+01+s*(-5.650000d+01)+s2*( 1.293000d+02)
1183      $           +s3*(-1.459000d+02)+s4*( 5.750000d+01)
1184 c
1185          qv  = prsccf/alinv/x*
1186      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1187          qsea= prsccf/alinv/x*
1188      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1189 c
1190          qu  = qv/3.0d0  + qsea/6.0d0
1191          qu  = qu*x
1192          ZUV=qu
1193          ZUB=qu
1194          qd  = qv/12.0d0 + qsea/6.0d0
1195          qd  = qd*x
1196          ZDV=qd
1197          ZDB=qd
1198          ZSB=qd
1199 c
1200          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1201             call WHIT6Q(x,mc2q2,cv,cs)
1202             qc = cv/alinv/2.0d0/PI + cs*alstpi
1203             qc  = qc*x
1204             ZCB=qc
1205          else
1206             qc = 0.0d0
1207             ZCB=qc
1208          endif
1209 c
1210          g   = WHIT6G(x,Q2)
1211          g   = g*x
1212          ZGL=g
1213 c
1214       else
1215 c over 100 GeV^2
1216 c
1217 c set scale s
1218          s=   log(  log(Q2/LAM52)/ log(Q52IT/LAM52))
1219          prsccf = 1.0d0
1220          alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
1221             s2=s**2
1222             s3=s2*s
1223             s4=s2**2
1224 c
1225 cccccc   WHIT6 quark (O100)
1226 c
1227       A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.621000d+00)
1228       A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.439000d+00)
1229       A2val= 2.837000d+00+s*( 6.460000d+00)+s2*( 4.100000d+00)
1230       Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.758000d-02)
1231       Cval = 1.728000d-01+s*(-2.493000d-02)+s2*( 1.451000d-01)
1232       A0sea= 3.340000d+00+s*(-5.610000d+00)+s2*( 5.000000d+01)
1233      $           +s3*(-2.207000d+02)+s4*( 3.028000d+02)
1234       B0sea=-2.402000d-01+s*(-4.090000d-01)+s2*( 2.263000d+00)
1235      $           +s3*(-1.050000d+01)+s4*( 1.487000d+01)
1236       BB0sea=8.790000d+00+s*(-8.860000d+00)+s2*( 1.640000d+02)
1237      $           +s3*(-7.120000d+02)+s4*( 9.730000d+02)
1238       C0sea= 9.160000d+00+s*( 9.290000d+00)+s2*(-2.784000d+02)
1239      $           +s3*( 1.175000d+03)+s4*(-1.592000d+03)
1240 c
1241          qv  = 1.0d0/alinv/x*
1242      $         (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1243          qsea= 1.0d0/alinv/x*
1244      $         A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1245 c
1246          qu  = qv/3.0d0  + qsea/6.0d0
1247          qu  = qu*x
1248          ZUV=qu
1249          ZUB=qu
1250          qd  = qv/12.0d0 + qsea/6.0d0
1251          qd  = qd*x
1252          ZDV=qd
1253          ZDB=qd
1254          ZSB=qd
1255          g   = WHIT6G(x,Q2)
1256          g   = g*x
1257          ZGL=g
1258 c
1259          if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1260       A0dcv=              s*( 1.219000d-01)+s2*( 6.200000d+00)
1261      $           +s3*(-2.504000d+01)+s4*( 3.098000d+01)
1262       A1dcv=              s*( 1.913000d+00)+s2*(-7.690000d+01)
1263      $           +s3*( 3.180000d+02)+s4*(-3.920000d+02)
1264       A2dcv=              s*(-7.160000d+00)+s2*( 2.503000d+02)
1265      $           +s3*(-1.062000d+03)+s4*( 1.308000d+03)
1266       A3dcv=              s*( 3.190000d+00)+s2*(-2.301000d+02)
1267      $           +s3*( 1.012000d+03)+s4*(-1.250000d+03)
1268       Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
1269      $           +s3*( 1.967000d+01)
1270       Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
1271      $           +s3*( 4.630000d+01)
1272       Adcs =              s*(-4.990000d-02)+s2*( 1.026000d-01)
1273      $           +s3*(-7.870000d-02)
1274       B0dcs=-3.610000d-01+s*(-5.760000d-01)+s2*( 2.257000d-01)
1275       B1dcs= 7.680000d+00+s*(-8.830000d+00)+s2*( 3.880000d+00)
1276       Cdcs = 2.548000d+00+s*( 6.910000d-01)+s2*(-8.700000d+00)
1277      $           +s3*( 1.065000d+01)
1278 c
1279          dcv = 1.0d0/alinv/x*
1280      $         (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
1281          dcs = 1.0d0/alinv/x*
1282      $         Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
1283 c
1284            call WHIT6Q(x,mc*mc/Q2,cv,cs)
1285            qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
1286            qc  = qc*x
1287            ZCB=qc
1288          else
1289            qc = 0.0d0
1290            ZCB=qc
1291          endif
1292       endif
1293 c
1294       return
1295       end
1296 c
1297 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1298       DOUBLE PRECISION function WHIT1G(x,Q2)
1299 c               input: x,Q2
1300 c               output: clg
1301 c                        (gluon dist.)
1302 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1303       implicit none
1304 c arg
1305       DOUBLE PRECISION Q2,x
1306 c const
1307       DOUBLE PRECISION q42it,q52it,lam42,lam52
1308       DOUBLE PRECISION alinv
1309 c local
1310       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1311       DOUBLE PRECISION s,s2,s3,s4,prsccf
1312       DOUBLE PRECISION x1
1313 c parameters
1314       parameter(lam42=0.16d0, lam52=0.091411319d0)
1315       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1316       parameter(alinv=137.036d0)
1317       common /scale/ s,s2,s3,s4,prsccf
1318 c
1319 c begin
1320       x1=1.0d0-x
1321 c
1322       if(Q2.lt.100.0d0) then
1323 c  under 100 GeV^2
1324 c
1325 cccccc   WHIT1 gluon (U100)
1326 c
1327       A0g = 2.000000d+00+s*(-3.280000d+00)+s2*( 2.894000d+00)
1328      $          +s3*(-1.561000d+00)+s4*( 8.180000d-01)
1329       B0g =              s*(-7.610000d-01)+s2*(-4.900000d-02)
1330      $          +s3*( 4.460000d-01)
1331       C0g = 3.000000d+00+s*( 1.586000d+00)+s2*(-9.490000d-01)
1332      $          +s3*( 2.425000d+00)
1333       A1g =              s*( 4.610000d-01)+s2*( 1.041000d-01)
1334      $          +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1335       AA1g=              s*( 9.680000d-03)+s2*(-4.170000d-01)
1336      $          +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1337       B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1338      $          +s3*(-5.070000d-01)
1339       C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1340      $          +s3*( 8.090000d-01)
1341       else
1342 c over 100 GeV^2
1343 c
1344 cccccc   WHIT1 gluon (O100)
1345 c
1346       A0g = 7.840000d-01+s*(-2.238000d+00)+s2*( 1.617000d+01)
1347      $          +s3*(-6.250000d+01)+s4*( 8.390000d+01)
1348       B0g =-4.030000d-01+s*(-1.307000d+00)+s2*( 8.780000d+00)
1349      $          +s3*(-3.580000d+01)+s4*( 5.350000d+01)
1350       C0g = 4.450000d+00+s*( 1.027000d+00)+s2*( 4.460000d+01)
1351      $          +s3*(-1.600000d+02)+s4*( 1.816000d+02)
1352       A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1353      $          +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1354       AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1355      $          +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1356       B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1357      $          +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1358       C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1359      $          +s3*(-2.072000d+00)
1360       endif
1361 c
1362          WHIT1G = prsccf/alinv/x*
1363      $            ( A0g * x**B0g * x1**C0g
1364      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
1365 c
1366       return
1367       end
1368 c
1369 cccccccccccccccccccccccccccccc
1370 c   QPM calculation
1371       subroutine WHIT1Q(x,r,cv,cs)
1372 ccc INPUTS : x,r=mc^2/Q^2
1373 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
1374 ccc                  cv <-- cv / ( alpha / 2PI)
1375 ccc                  cs <-- cs / ( alpha_s/2PI)
1376 c
1377       implicit none
1378 c arg
1379       DOUBLE PRECISION x,r
1380       DOUBLE PRECISION cv,cs
1381 c CONST
1382       DOUBLE PRECISION ec,mc
1383       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1384 c  N=15 Gauss int. weights and points
1385       integer GN,i
1386       parameter(GN=15)
1387       DOUBLE PRECISION XG(GN), XW(GN)
1388       DATA (XG(i),i=1,GN)/6.003741d-03,
1389      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1390      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1391      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1392      $ 9.686367d-01, 9.939963d-01/
1393       DATA (XW(i),i=1,GN)/1.537662d-02,
1394      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1395      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1396      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1397      $ 3.518302d-02, 1.537662d-02/
1398 c local
1399       DOUBLE PRECISION sum, y,z,beta,w,WHIT1G,L
1400       parameter(L=4.0d0)
1401       DOUBLE PRECISION x1,rx,z1,rz
1402 c
1403 c begin
1404       x1=1.0d0-x
1405       rx=4.0d0*r*x
1406 c
1407 c direct
1408          beta=dsqrt(1.0d0-rx/x1)
1409          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1410      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1411      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1412          cv = 3.0d0*ec**2 * w / x
1413 c
1414 c resolved
1415       sum=0.0d0
1416       do 10 i=1,GN
1417          y= x+rx + (x1-rx)*XG(i)**L
1418          z=x/y
1419          z1=1.0d0-z
1420          rz=4.0d0*r*z
1421          beta=dsqrt(1.0d0-rz/z1)
1422          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1423      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1424      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1425 c
1426          sum= sum + w * WHIT1G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1427 c
1428    10 continue
1429 c
1430          cs = 0.5d0/x * (x1-rx) * sum
1431 c
1432       return
1433       end
1434 c
1435 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1436       DOUBLE PRECISION function WHIT2G(x,Q2)
1437 c               input: x,Q2
1438 c               output: clg
1439 c                        (gluon dist.)
1440 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1441       implicit none
1442 c arg
1443       DOUBLE PRECISION Q2,x
1444 c const
1445       DOUBLE PRECISION q42it,q52it,lam42,lam52
1446       DOUBLE PRECISION alinv
1447 c local
1448       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1449       DOUBLE PRECISION s,s2,s3,s4,prsccf
1450       DOUBLE PRECISION x1
1451 c parameters
1452       parameter(lam42=0.16d0, lam52=0.091411319d0)
1453       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1454       parameter(alinv=137.036d0)
1455       common /scale/ s,s2,s3,s4,prsccf
1456 c
1457 c begin
1458       x1=1.0d0-x
1459 c
1460       if(Q2.lt.100.0d0) then
1461 c  under 100 GeV^2
1462 c
1463 cccccc   WHIT2 gluon (U100)
1464 c
1465       A0g = 5.000000d+00+s*(-1.499000d+01)+s2*( 2.617000d+01)
1466      $          +s3*(-2.530000d+01)+s4*( 1.012000d+01)
1467       B0g =              s*(-9.370000d-01)+s2*( 4.100000d-01)
1468      $          +s3*( 3.390000d-02)
1469       C0g = 9.000000d+00+s*( 7.090000d-01)+s2*( 3.118000d+00)
1470      $          +s3*(-5.820000d-04)
1471       A1g =              s*( 4.610000d-01)+s2*( 1.041000d-01)
1472      $          +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1473       AA1g=              s*( 9.680000d-03)+s2*(-4.170000d-01)
1474      $          +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1475       B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1476      $          +s3*(-5.070000d-01)
1477       C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1478      $          +s3*( 8.090000d-01)
1479       else
1480 c over 100 GeV^2
1481 c
1482 cccccc   WHIT2 gluon (O100)
1483 c
1484       A0g = 1.095000d+00+s*(-2.388000d+00)+s2*( 9.190000d+00)
1485      $          +s3*(-3.032000d+01)+s4*( 3.480000d+01)
1486       B0g =-4.410000d-01+s*(-9.070000d-01)+s2*( 4.680000d+00)
1487      $          +s3*(-1.866000d+01)+s4*( 2.717000d+01)
1488       C0g = 1.099000d+01+s*( 4.710000d+00)+s2*( 2.801000d+01)
1489      $          +s3*(-1.279000d+02)+s4*( 1.640000d+02)
1490       A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1491      $          +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1492       AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1493      $          +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1494       B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1495      $          +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1496       C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1497      $          +s3*(-2.072000d+00)
1498       endif
1499 c
1500          WHIT2G = prsccf/alinv/x*
1501      $            ( A0g * x**B0g * x1**C0g
1502      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
1503 c
1504       return
1505       end
1506 c
1507 cccccccccccccccccccccccccccccc
1508 c   QPM calculation
1509       subroutine WHIT2Q(x,r,cv,cs)
1510 ccc INPUTS : x,r=mc^2/Q^2
1511 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
1512 ccc                  cv <-- cv / ( alpha / 2PI)
1513 ccc                  cs <-- cs / ( alpha_s/2PI)
1514 c
1515       implicit none
1516 c arg
1517       DOUBLE PRECISION x,r
1518       DOUBLE PRECISION cv,cs
1519 c CONST
1520       DOUBLE PRECISION ec,mc
1521       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1522 c  N=15 Gauss int. weights and points
1523       integer GN,i
1524       parameter(GN=15)
1525       DOUBLE PRECISION XG(GN), XW(GN)
1526       DATA (XG(i),i=1,GN)/6.003741d-03,
1527      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1528      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1529      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1530      $ 9.686367d-01, 9.939963d-01/
1531       DATA (XW(i),i=1,GN)/1.537662d-02,
1532      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1533      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1534      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1535      $ 3.518302d-02, 1.537662d-02/
1536 c local
1537       DOUBLE PRECISION sum, y,z,beta,w,WHIT2G,L
1538       parameter(L=4.0d0)
1539       DOUBLE PRECISION x1,rx,z1,rz
1540 c
1541 c begin
1542       x1=1.0d0-x
1543       rx=4.0d0*r*x
1544 c
1545 c direct
1546          beta=dsqrt(1.0d0-rx/x1)
1547          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1548      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1549      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1550          cv = 3.0d0*ec**2 * w / x
1551 c
1552 c resolved
1553       sum=0.0d0
1554       do 10 i=1,GN
1555          y= x+rx + (x1-rx)*XG(i)**L
1556          z=x/y
1557          z1=1.0d0-z
1558          rz=4.0d0*r*z
1559          beta=dsqrt(1.0d0-rz/z1)
1560          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1561      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1562      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1563 c
1564          sum= sum + w * WHIT2G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1565 c
1566    10 continue
1567 c
1568          cs = 0.5d0/x * (x1-rx) * sum
1569 c
1570       return
1571       end
1572 c
1573 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1574       DOUBLE PRECISION function WHIT3G(x,Q2)
1575 c               input: x,Q2
1576 c               output: clg
1577 c                        (gluon dist.)
1578 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1579       implicit none
1580 c arg
1581       DOUBLE PRECISION Q2,x
1582 c const
1583       DOUBLE PRECISION q42it,q52it,lam42,lam52
1584       DOUBLE PRECISION alinv
1585 c local
1586       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1587       DOUBLE PRECISION s,s2,s3,s4,prsccf
1588       DOUBLE PRECISION x1
1589 c parameters
1590       parameter(lam42=0.16d0, lam52=0.091411319d0)
1591       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1592       parameter(alinv=137.036d0)
1593       common /scale/ s,s2,s3,s4,prsccf
1594 c
1595 c begin
1596       x1=1.0d0-x
1597 c
1598       if(Q2.lt.100.0d0) then
1599 c  under 100 GeV^2
1600 c
1601 cccccc   WHIT3 gluon (U100)
1602 c
1603       A0g = 8.000000d+00+s*(-2.864000d+01)+s2*( 5.590000d+01)
1604      $          +s3*(-5.760000d+01)+s4*( 2.366000d+01)
1605       B0g =              s*(-9.870000d-01)+s2*( 5.100000d-01)
1606      $          +s3*(-6.670000d-02)
1607       C0g = 1.500000d+01+s*( 3.310000d-01)+s2*( 3.500000d+00)
1608      $          +s3*( 8.920000d-01)
1609       A1g =              s*( 4.610000d-01)+s2*( 1.041000d-01)
1610      $          +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1611       AA1g=              s*( 9.680000d-03)+s2*(-4.170000d-01)
1612      $          +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1613       B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1614      $          +s3*(-5.070000d-01)
1615       C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1616      $          +s3*( 8.090000d-01)
1617       else
1618 c over 100 GeV^2
1619 c
1620 cccccc   WHIT3 gluon (O100)
1621 c
1622       A0g = 1.270000d+00+s*(-2.817000d+00)+s2*( 5.740000d+00)
1623      $          +s3*(-1.327000d+01)+s4*( 1.268000d+01)
1624       B0g =-4.610000d-01+s*(-8.170000d-01)+s2*( 3.320000d+00)
1625      $          +s3*(-1.296000d+01)+s4*( 1.893000d+01)
1626       C0g = 1.721000d+01+s*( 1.257000d+00)+s2*( 5.050000d+01)
1627      $          +s3*(-2.761000d+02)+s4*( 4.900000d+02)
1628       A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1629      $          +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1630       AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1631      $          +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1632       B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1633      $          +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1634       C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1635      $          +s3*(-2.072000d+00)
1636       endif
1637 c
1638          WHIT3G = prsccf/alinv/x*
1639      $            ( A0g * x**B0g * x1**C0g
1640      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
1641 c
1642       return
1643       end
1644 c
1645 cccccccccccccccccccccccccccccc
1646 c   QPM calculation
1647       subroutine WHIT3Q(x,r,cv,cs)
1648 ccc INPUTS : x,r=mc^2/Q^2
1649 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
1650 ccc                  cv <-- cv / ( alpha / 2PI)
1651 ccc                  cs <-- cs / ( alpha_s/2PI)
1652 c
1653       implicit none
1654 c arg
1655       DOUBLE PRECISION x,r
1656       DOUBLE PRECISION cv,cs
1657 c CONST
1658       DOUBLE PRECISION ec,mc
1659       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1660 c  N=15 Gauss int. weights and points
1661       integer GN,i
1662       parameter(GN=15)
1663       DOUBLE PRECISION XG(GN), XW(GN)
1664       DATA (XG(i),i=1,GN)/6.003741d-03,
1665      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1666      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1667      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1668      $ 9.686367d-01, 9.939963d-01/
1669       DATA (XW(i),i=1,GN)/1.537662d-02,
1670      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1671      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1672      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1673      $ 3.518302d-02, 1.537662d-02/
1674 c local
1675       DOUBLE PRECISION sum, y,z,beta,w,WHIT3G,L
1676       parameter(L=4.0d0)
1677       DOUBLE PRECISION x1,rx,z1,rz
1678 c
1679 c begin
1680       x1=1.0d0-x
1681       rx=4.0d0*r*x
1682 c
1683 c direct
1684          beta=dsqrt(1.0d0-rx/x1)
1685          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1686      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1687      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1688          cv = 3.0d0*ec**2 * w / x
1689 c
1690 c resolved
1691       sum=0.0d0
1692       do 10 i=1,GN
1693          y= x+rx + (x1-rx)*XG(i)**L
1694          z=x/y
1695          z1=1.0d0-z
1696          rz=4.0d0*r*z
1697          beta=dsqrt(1.0d0-rz/z1)
1698          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1699      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1700      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1701 c
1702          sum= sum + w * WHIT3G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1703 c
1704    10 continue
1705 c
1706          cs = 0.5d0/x * (x1-rx) * sum
1707 c
1708       return
1709       end
1710 c
1711 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1712       DOUBLE PRECISION function WHIT4G(x,Q2)
1713 c               input: x,Q2
1714 c               output: clg
1715 c                        (gluon dist.)
1716 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1717       implicit none
1718 c arg
1719       DOUBLE PRECISION Q2,x
1720 c const
1721       DOUBLE PRECISION q42it,q52it,lam42,lam52
1722       DOUBLE PRECISION alinv
1723 c local
1724       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1725       DOUBLE PRECISION s,s2,s3,s4,prsccf
1726       DOUBLE PRECISION x1
1727 c parameters
1728       parameter(lam42=0.16d0, lam52=0.091411319d0)
1729       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1730       parameter(alinv=137.036d0)
1731       common /scale/ s,s2,s3,s4,prsccf
1732 c
1733 c begin
1734       x1=1.0d0-x
1735 c
1736       if(Q2.lt.100.0d0) then
1737 c  under 100 GeV^2
1738 c
1739 cccccc   WHIT4 gluon (U100)
1740 c
1741       A0g = 4.000000d+00+s*(-9.400000d+00)+s2*( 1.555000d+01)
1742      $          +s3*(-1.450000d+01)+s4*( 5.470000d+00)
1743       B0g =              s*(-1.142000d+00)+s2*( 1.034000d+00)
1744      $          +s3*(-4.410000d-01)
1745       C0g = 3.000000d+00+s*( 8.720000d-01)+s2*( 1.006000d+00)
1746      $          +s3*( 3.560000d-01)
1747       A1g =              s*( 6.020000d-01)+s2*( 5.090000d-01)
1748      $          +s3*(-2.054000d+00)+s4*( 1.392000d+00)
1749       AA1g=              s*(-9.220000d-02)+s2*(-1.899000d+00)
1750      $          +s3*( 4.180000d+00)+s4*(-2.494000d+00)
1751       B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
1752      $          +s3*( 1.116000d+00)
1753       C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
1754      $          +s3*( 7.330000d-01)+s4*(-7.620000d-01)
1755       else
1756 c over 100 GeV^2
1757 c
1758 cccccc   WHIT4 gluon (O100)
1759 c
1760       A0g = 1.384000d+00+s*(-2.455000d+00)+s2*( 8.940000d+00)
1761      $          +s3*(-2.906000d+01)+s4*( 3.710000d+01)
1762       B0g =-4.420000d-01+s*(-7.190000d-01)+s2*( 2.961000d+00)
1763      $          +s3*(-1.209000d+01)+s4*( 1.916000d+01)
1764       C0g = 4.210000d+00+s*( 2.524000d+00)+s2*( 1.003000d+01)
1765      $          +s3*(-1.827000d+01)+s4*( 2.162000d+00)
1766       A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
1767      $          +s3*( 7.260000d+00)+s4*(-1.839000d+01)
1768       AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
1769      $          +s3*(-6.660000d+00)+s4*( 1.923000d+01)
1770       B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
1771      $          +s3*( 1.075000d+01)+s4*(-1.993000d+01)
1772       C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
1773      $          +s3*(-2.057000d+00)
1774       endif
1775 c
1776          WHIT4G = prsccf/alinv/x*
1777      $            ( A0g * x**B0g * x1**C0g
1778      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
1779 c
1780       return
1781       end
1782 c
1783 cccccccccccccccccccccccccccccc
1784 c   QPM calculation
1785       subroutine WHIT4Q(x,r,cv,cs)
1786 ccc INPUTS : x,r=mc^2/Q^2
1787 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
1788 ccc                  cv <-- cv / ( alpha / 2PI)
1789 ccc                  cs <-- cs / ( alpha_s/2PI)
1790 c
1791       implicit none
1792 c arg
1793       DOUBLE PRECISION x,r
1794       DOUBLE PRECISION cv,cs
1795 c CONST
1796       DOUBLE PRECISION ec,mc
1797       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1798 c  N=15 Gauss int. weights and points
1799       integer GN,i
1800       parameter(GN=15)
1801       DOUBLE PRECISION XG(GN), XW(GN)
1802       DATA (XG(i),i=1,GN)/6.003741d-03,
1803      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1804      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1805      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1806      $ 9.686367d-01, 9.939963d-01/
1807       DATA (XW(i),i=1,GN)/1.537662d-02,
1808      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1809      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1810      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1811      $ 3.518302d-02, 1.537662d-02/
1812 c local
1813       DOUBLE PRECISION sum, y,z,beta,w,WHIT4G,L
1814       parameter(L=4.0d0)
1815       DOUBLE PRECISION x1,rx,z1,rz
1816 c
1817 c begin
1818       x1=1.0d0-x
1819       rx=4.0d0*r*x
1820 c
1821 c direct
1822          beta=dsqrt(1.0d0-rx/x1)
1823          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1824      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1825      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1826          cv = 3.0d0*ec**2 * w / x
1827 c
1828 c resolved
1829       sum=0.0d0
1830       do 10 i=1,GN
1831          y= x+rx + (x1-rx)*XG(i)**L
1832          z=x/y
1833          z1=1.0d0-z
1834          rz=4.0d0*r*z
1835          beta=dsqrt(1.0d0-rz/z1)
1836          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1837      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1838      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1839 c
1840          sum= sum + w * WHIT4G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1841 c
1842    10 continue
1843 c
1844          cs = 0.5d0/x * (x1-rx) * sum
1845 c
1846       return
1847       end
1848 c
1849 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1850       DOUBLE PRECISION function WHIT5G(x,Q2)
1851 c               input: x,Q2
1852 c               output: clg
1853 c                        (gluon dist.)
1854 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1855       implicit none
1856 c arg
1857       DOUBLE PRECISION Q2,x
1858 c const
1859       DOUBLE PRECISION q42it,q52it,lam42,lam52
1860       DOUBLE PRECISION alinv
1861 c local
1862       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1863       DOUBLE PRECISION s,s2,s3,s4,prsccf
1864       DOUBLE PRECISION x1
1865 c parameters
1866       parameter(lam42=0.16d0, lam52=0.091411319d0)
1867       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1868       parameter(alinv=137.036d0)
1869       common /scale/ s,s2,s3,s4,prsccf
1870 c
1871 c begin
1872       x1=1.0d0-x
1873 c
1874       if(Q2.lt.100.0d0) then
1875 c  under 100 GeV^2
1876 c
1877 cccccc   WHIT5 gluon (U100)
1878 c
1879       A0g = 1.000000d+01+s*(-3.400000d+01)+s2*( 6.900000d+01)
1880      $          +s3*(-7.530000d+01)+s4*( 3.230000d+01)
1881       B0g =              s*(-1.126000d+00)+s2*( 9.260000d-01)
1882      $          +s3*(-3.930000d-01)
1883       C0g = 9.000000d+00+s*( 4.810000d-01)+s2*( 3.200000d+00)
1884      $          +s3*(-3.470000d-01)
1885       A1g =              s*( 6.020000d-01)+s2*( 5.090000d-01)
1886      $          +s3*(-2.054000d+00)+s4*( 1.392000d+00)
1887       AA1g=              s*(-9.220000d-02)+s2*(-1.899000d+00)
1888      $          +s3*( 4.180000d+00)+s4*(-2.494000d+00)
1889       B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
1890      $          +s3*( 1.116000d+00)
1891       C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
1892      $          +s3*( 7.330000d-01)+s4*(-7.620000d-01)
1893       else
1894 c over 100 GeV^2
1895 c
1896 cccccc   WHIT5 gluon (O100)
1897 c
1898       A0g = 1.995000d+00+s*(-3.260000d+00)+s2*( 1.818000d+00)
1899      $          +s3*( 1.711000d+00)+s4*(-4.990000d+00)
1900       B0g =-4.660000d-01+s*(-6.100000d-01)+s2*( 1.691000d+00)
1901      $          +s3*(-6.680000d+00)+s4*( 1.019000d+01)
1902       C0g = 1.075000d+01+s*( 5.420000d+00)+s2*( 6.550000d+00)
1903      $          +s3*(-2.297000d+01)+s4*( 1.867000d+01)
1904       A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
1905      $          +s3*( 7.260000d+00)+s4*(-1.839000d+01)
1906       AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
1907      $          +s3*(-6.660000d+00)+s4*( 1.923000d+01)
1908       B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
1909      $          +s3*( 1.075000d+01)+s4*(-1.993000d+01)
1910       C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
1911      $          +s3*(-2.057000d+00)
1912       endif
1913 c
1914          WHIT5G = prsccf/alinv/x*
1915      $            ( A0g * x**B0g * x1**C0g
1916      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
1917 c
1918       return
1919       end
1920 c
1921 cccccccccccccccccccccccccccccc
1922 c   QPM calculation
1923       subroutine WHIT5Q(x,r,cv,cs)
1924 ccc INPUTS : x,r=mc^2/Q^2
1925 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
1926 ccc                  cv <-- cv / ( alpha / 2PI)
1927 ccc                  cs <-- cs / ( alpha_s/2PI)
1928 c
1929       implicit none
1930 c arg
1931       DOUBLE PRECISION x,r
1932       DOUBLE PRECISION cv,cs
1933 c CONST
1934       DOUBLE PRECISION ec,mc
1935       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1936 c  N=15 Gauss int. weights and points
1937       integer GN,i
1938       parameter(GN=15)
1939       DOUBLE PRECISION XG(GN), XW(GN)
1940       DATA (XG(i),i=1,GN)/6.003741d-03,
1941      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1942      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1943      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1944      $ 9.686367d-01, 9.939963d-01/
1945       DATA (XW(i),i=1,GN)/1.537662d-02,
1946      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1947      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1948      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1949      $ 3.518302d-02, 1.537662d-02/
1950 c local
1951       DOUBLE PRECISION sum, y,z,beta,w,WHIT5G,L
1952       parameter(L=4.0d0)
1953       DOUBLE PRECISION x1,rx,z1,rz
1954 c
1955 c begin
1956       x1=1.0d0-x
1957       rx=4.0d0*r*x
1958 c
1959 c direct
1960          beta=dsqrt(1.0d0-rx/x1)
1961          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1962      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1963      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1964          cv = 3.0d0*ec**2 * w / x
1965 c
1966 c resolved
1967       sum=0.0d0
1968       do 10 i=1,GN
1969          y= x+rx + (x1-rx)*XG(i)**L
1970          z=x/y
1971          z1=1.0d0-z
1972          rz=4.0d0*r*z
1973          beta=dsqrt(1.0d0-rz/z1)
1974          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1975      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1976      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
1977 c
1978          sum= sum + w * WHIT5G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1979 c
1980    10 continue
1981 c
1982          cs = 0.5d0/x * (x1-rx) * sum
1983 c
1984       return
1985       end
1986 c
1987 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1988       DOUBLE PRECISION function WHIT6G(x,Q2)
1989 c               input: x,Q2
1990 c               output: clg
1991 c                        (gluon dist.)
1992 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1993       implicit none
1994 c arg
1995       DOUBLE PRECISION Q2,x
1996 c const
1997       DOUBLE PRECISION q42it,q52it,lam42,lam52
1998       DOUBLE PRECISION alinv
1999 c local
2000       DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
2001       DOUBLE PRECISION s,s2,s3,s4,prsccf
2002       DOUBLE PRECISION x1
2003 c parameters
2004       parameter(lam42=0.16d0, lam52=0.091411319d0)
2005       parameter(Q42IT=4.0d0, Q52IT=100.0d0)
2006       parameter(alinv=137.036d0)
2007       common /scale/ s,s2,s3,s4,prsccf
2008 c
2009 c begin
2010       x1=1.0d0-x
2011 c
2012       if(Q2.lt.100.0d0) then
2013 c  under 100 GeV^2
2014 c
2015 cccccc   WHIT6 gluon (U100)
2016 c
2017       A0g = 1.600000d+01+s*(-6.100000d+01)+s2*( 1.278000d+02)
2018      $          +s3*(-1.399000d+02)+s4*( 5.990000d+01)
2019       B0g =              s*(-1.109000d+00)+s2*( 8.450000d-01)
2020      $          +s3*(-3.510000d-01)
2021       C0g = 1.500000d+01+s*( 1.596000d-01)+s2*( 4.180000d+00)
2022      $          +s3*(-1.765000d-01)
2023       A1g =              s*( 6.020000d-01)+s2*( 5.090000d-01)
2024      $          +s3*(-2.054000d+00)+s4*( 1.392000d+00)
2025       AA1g=              s*(-9.220000d-02)+s2*(-1.899000d+00)
2026      $          +s3*( 4.180000d+00)+s4*(-2.494000d+00)
2027       B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
2028      $          +s3*( 1.116000d+00)
2029       C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
2030      $          +s3*( 7.330000d-01)+s4*(-7.620000d-01)
2031       else
2032 c over 100 GeV^2
2033 c
2034 cccccc   WHIT6 gluon (O100)
2035 c
2036       A0g = 2.378000d+00+s*(-4.380000d+00)+s2*( 5.850000d-01)
2037      $          +s3*( 8.340000d+00)+s4*(-9.920000d+00)
2038       B0g =-4.790000d-01+s*(-6.070000d-01)+s2*( 1.458000d+00)
2039      $          +s3*(-6.030000d+00)+s4*( 9.330000d+00)
2040       C0g = 1.706000d+01+s*( 4.960000d+00)+s2*( 2.497000d+01)
2041      $          +s3*(-1.582000d+02)+s4*( 2.954000d+02)
2042       A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
2043      $          +s3*( 7.260000d+00)+s4*(-1.839000d+01)
2044       AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
2045      $          +s3*(-6.660000d+00)+s4*( 1.923000d+01)
2046       B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
2047      $          +s3*( 1.075000d+01)+s4*(-1.993000d+01)
2048       C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
2049      $          +s3*(-2.057000d+00)
2050       endif
2051 c
2052          WHIT6G = prsccf/alinv/x*
2053      $            ( A0g * x**B0g * x1**C0g
2054      $             +(A1g+AA1g*x) * x**B1g * x1**C1g )
2055 c
2056       return
2057       end
2058 c
2059 cccccccccccccccccccccccccccccc
2060 c   QPM calculation
2061       subroutine WHIT6Q(x,r,cv,cs)
2062 ccc INPUTS : x,r=mc^2/Q^2
2063 ccc OUTPUTS:   cv,cs (valence- and sea- charm quark dist)
2064 ccc                  cv <-- cv / ( alpha / 2PI)
2065 ccc                  cs <-- cs / ( alpha_s/2PI)
2066 c
2067       implicit none
2068 c arg
2069       DOUBLE PRECISION x,r
2070       DOUBLE PRECISION cv,cs
2071 c CONST
2072       DOUBLE PRECISION ec,mc
2073       parameter(ec=2.0d0/3.0d0,mc=1.5d0)
2074 c  N=15 Gauss int. weights and points
2075       integer GN,i
2076       parameter(GN=15)
2077       DOUBLE PRECISION XG(GN), XW(GN)
2078       DATA (XG(i),i=1,GN)/6.003741d-03,
2079      $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
2080      $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
2081      $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
2082      $ 9.686367d-01, 9.939963d-01/
2083       DATA (XW(i),i=1,GN)/1.537662d-02,
2084      $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
2085      $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
2086      $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
2087      $ 3.518302d-02, 1.537662d-02/
2088 c local
2089       DOUBLE PRECISION sum, y,z,beta,w,WHIT6G,L
2090       parameter(L=4.0d0)
2091       DOUBLE PRECISION x1,rx,z1,rz
2092 c
2093 c begin
2094       x1=1.0d0-x
2095       rx=4.0d0*r*x
2096 c
2097 c direct
2098          beta=dsqrt(1.0d0-rx/x1)
2099          w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
2100      $        +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
2101      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
2102          cv = 3.0d0*ec**2 * w / x
2103 c
2104 c resolved
2105       sum=0.0d0
2106       do 10 i=1,GN
2107          y= x+rx + (x1-rx)*XG(i)**L
2108          z=x/y
2109          z1=1.0d0-z
2110          rz=4.0d0*r*z
2111          beta=dsqrt(1.0d0-rz/z1)
2112          w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
2113      $        +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
2114      $         *  log( (1.0d0+beta)/(1.0d0-beta) ))
2115 c
2116          sum= sum + w * WHIT6G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
2117 c
2118    10 continue
2119 c
2120          cs = 0.5d0/x * (x1-rx) * sum
2121 c
2122       return
2123       end