LHAPDF/lhapdf5.3.1/wrapacfgpg.f

   1       subroutine ACFGPevolvep(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       real*4 par,par2,calc,celc
   6       common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
   7       character*16 name(nmxset)
   8       integer nmem(nmxset),ndef(nmxset),mmem
   9       common/NAME/name,nmem,ndef,mmem
  10       integer nset
  11
  12       save
  13
  14       if(imem.eq.1) then
  15         call ACFGP1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  16
  17       elseif(imem.eq.2) then
  18         call ACFGP2(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  19
  20       elseif(imem.eq.3.or.imem.eq.0) then
  21         call SFAFG1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
  22
  23       else
  24         CONTINUE
  25       endif
  26
  27       pdf(-6)= 0.0d0
  28       pdf(6)= 0.0d0
  29       pdf(-5)= 0.0d0
  30       pdf(5 )= 0.0d0
  31       pdf(-4)= chm
  32       pdf(4 )= chm
  33       pdf(-3)= str
  34       pdf(3 )= str
  35       pdf(-2)= usea
  36       pdf(2 )= upv
  37       pdf(-1)= dsea
  38       pdf(1 )= dnv
  39       pdf(0 )= glu
  40
  41       return
  42 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  43       entry ACFGPread(nset)
  44       read(1,*)nmem(nset),ndef(nset)
  45       do iset = 1,3
  46         read(1,*)(par(j,iset),j=1,4)
  47         read(1,*)(par(j,iset),j=5,8)
  48         read(1,*)(par(j,iset),j=9,12)
  49         read(1,*)(par(j,iset),j=13,16)
  50         read(1,*)(par(j,iset),j=17,20)
  51         read(1,*)(par(j,iset),j=21,24)
  52         read(1,*)(par(j,iset),j=25,28)
  53         read(1,*)(par(j,iset),j=29,30)
  54         do j=1,32
  55           do k=1,20
  56             read(1,*)(calc(i,k,j,iset),i=1,4)
  57             read(1,*)(calc(i,k,j,iset),i=5,8)
  58           enddo
  59         enddo
  60       enddo
  61 c last one
  62       read(1,*)(par2(j),j=1,4)
  63       read(1,*)(par2(j),j=5,8)
  64       read(1,*)(par2(j),j=9,12)
  65       read(1,*)(par2(j),j=13,16)
  66       read(1,*)(par2(j),j=17,20)
  67       read(1,*)(par2(j),j=21,24)
  68       read(1,*)(par2(j),j=25,28)
  69       read(1,*)(par2(j),j=29,30)
  70       do j=1,32
  71         do k=1,20
  72           read(1,*)(celc(i,k,j),i=1,4)
  73           read(1,*)(celc(i,k,j),i=5,8)
  74         enddo
  75       enddo
  76
  77       return
  78
  79
  80 c
  81 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  82       entry ACFGPalfa(alfas,qalfa)
  83         call getnset(iset)
  84         call GetOrderAsM(iset,iord)
  85         call Getlam4M(iset,imem,qcdl4)
  86         call Getlam5M(iset,imem,qcdl5)
  87         call aspdflib(alfas,Qalfa,iord,qcdl5)
  88
  89       return
  90 c
  91 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  92       entry ACFGPinit(Eorder,Q2fit)
  93       return
  94 c
  95 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  96       entry ACFGPpdf(mem)
  97       imem = mem
  98       return
  99 c
 100  1000 format(5e13.5)
 101       end
 102 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 103       SUBROUTINE ACFGP1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
 104 C
 105 C     INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
 106 C     QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
 107 C
 108 C     THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
 109 C     CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
 110 C
 111 C     THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
 112 C     BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
 113 C
 114 C     THE PROGRAM WORKS FOR  2. GEV**2 < Q2 <5.5E+5   AND .00137 < X < .9986
 115 C
 116 C     THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
 117 C
 118 C     THE VALUE OF LAMBDA-MSB IS 200 MEV
 119 C
 120 C     THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
 121 C                                  X*U=X*U(X,Q2)
 122 C                                  X*D= ...
 123 C                                  X*S= ...
 124 C                                  X*C= ...  (MASSLESS CHARM WITH     C(X,2.)=0)
 125 C                                 X*CM= ...  (MASSIVE CHARM WITH MC=1.5 GEV )
 126 C                              X*GLU=GLUON(X,Q2)*X
 127 C
 128 C
 129 C                    F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
 130 C                    F2C=  "        "         "     WITH MASSIVE CHARM
 131 C
 132       double precision
 133      +       DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
 134       REAL    X, Q, UV, DV, UB, DB, SB, CB, BB, GL
 135       REAL       Q2
 136       common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
 137       DIMENSION XPDF(7),CALCO(8,20,32)
 138       COMMON/W5051I7/CALCO
 139       EXTERNAL AFCPLU
 140       DATA ZERO/0.0/
 141 C----------------------------------------------------------------------
 142        DATA ISTART/0/
 143        SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
 144 C
 145       IF (ISTART.EQ.0) THEN
 146         ISTART=1
 147         DO 10 K=1,32
 148         DO 10 I=1,20
 149         DO 10 M=1,8
 150    10   CALCO(M,I,K) = CALC(M,I,K,1)
 151            OWLAM=PAR(1,1)
 152            OWLAM2=OWLAM**2
 153            Q02=PAR(30,1)
 154            FLAV=PAR(25,1)
 155            DELTA=PAR(29,1)
 156            CALL WATE32
 157          ENDIF
 158 C
 159       X = DX
 160       Q = DQ
 161       Q2 = Q*Q
 162       IDQ2=2
 163       SB=0.
 164       IF(Q2-Q02) 1,1,2
 165     2 IF(IDQ2-1) 1,1,3
 166     3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
 167     1 CONTINUE
 168       CALL AURGAM(8,0,X,SB,XPDF(7))
 169       CALL AURGAM(7,0,X,SB,SING)
 170       CALL AURGAM(4,0,X,SB,DPLUSNS)
 171       CALL AURGAM(3,0,X,SB,CPLUSNS)
 172       CALL AURGAM(5,0,X,SB,UPLUSNS)
 173       CALL AURGAM(6,0,X,SB,SPLUSNS)
 174       XPDF(3) = CPLUSNS
 175       XPDF(4) = DPLUSNS
 176       XPDF(5) = UPLUSNS
 177       XPDF(6) = SPLUSNS
 178       XPDF(1) = SING
 179 C
 180       ADD = XPDF(1)/FLAV
 181       UPLUS=XPDF(5)+ADD
 182       DPLUS=-XPDF(4)+ADD
 183       SPLUS=-XPDF(6)+ADD
 184       CPLUS=-XPDF(3)+ADD
 185       UB=UPLUS*0.5
 186       UV=UB
 187       DB=DPLUS*0.5
 188       DV=DB
 189       SB=SPLUS*0.5
 190       CB=CPLUS*0.5
 191       SING=XPDF(1)
 192       GLU=XPDF(7)
 193       GL=GLU
 194 *
 195       DUV=MAX(ZERO,UV)
 196       DDV=MAX(ZERO,DV)
 197       DUB=MAX(ZERO,UB)
 198       DDB=MAX(ZERO,DB)
 199       DSB=MAX(ZERO,SB)
 200       DCB=MAX(ZERO,CB)
 201       DGL=MAX(ZERO,GL)
 202 C
 203       RETURN
 204       END
 205 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 206       SUBROUTINE ACFGP2(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
 207 C
 208 C     INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
 209 C     QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
 210 C
 211 C     THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
 212 C     CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
 213 C
 214 C     THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
 215 C     BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
 216 C
 217 C     THE PROGRAM WORKS FOR  2. GEV**2 < Q2 <5.5E+5   AND .00137 < X < .9986
 218 C
 219 C     THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
 220 C
 221 C     THE VALUE OF LAMBDA-MSB IS 200 MEV
 222 C
 223 C     THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
 224 C                                  X*U=X*U(X,Q2)
 225 C                                  X*D= ...
 226 C                                  X*S= ...
 227 C                                  X*C= ...  (MASSLESS CHARM WITH     C(X,2.)=0)
 228 C                                 X*CM= ...  (MASSIVE CHARM WITH MC=1.5 GEV )
 229 C                              X*GLU=GLUON(X,Q2)*X
 230 C
 231 C
 232 C                    F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
 233 C                    F2C=  "        "         "     WITH MASSIVE CHARM
 234 C
 235       double precision
 236      +       DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
 237       REAL    X, Q, UV, DV, UB, DB, SB, CB, BB, GL
 238       REAL       Q2
 239       common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
 240       DIMENSION XPDF(7),CALCO(8,20,32)
 241       COMMON/W5051I7/CALCO
 242       EXTERNAL AFCPLU
 243       DATA ZERO/0.0/
 244 C----------------------------------------------------------------------
 245        DATA ISTART/0/
 246        SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
 247 C
 248       IF (ISTART.EQ.0) THEN
 249         ISTART=1
 250         DO 10 K=1,32
 251         DO 10 I=1,20
 252         DO 10 M=1,8
 253    10   CALCO(M,I,K) = CALC(M,I,K,2)
 254            OWLAM=PAR(1,2)
 255            OWLAM2=OWLAM**2
 256            Q02=PAR(30,2)
 257            FLAV=PAR(25,2)
 258            DELTA=PAR(29,2)
 259            CALL WATE32
 260          ENDIF
 261 C
 262       X = DX
 263       Q = DQ
 264       Q2 = Q*Q
 265       IDQ2=2
 266       SB=0.
 267       IF(Q2-Q02) 1,1,2
 268     2 IF(IDQ2-1) 1,1,3
 269     3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
 270     1 CONTINUE
 271       CALL AURGAM(8,0,X,SB,XPDF(7))
 272       CALL AURGAM(7,0,X,SB,SING)
 273       CALL AURGAM(4,0,X,SB,DPLUSNS)
 274       CALL AURGAM(3,0,X,SB,CPLUSNS)
 275       CALL AURGAM(5,0,X,SB,UPLUSNS)
 276       CALL AURGAM(6,0,X,SB,SPLUSNS)
 277       XPDF(3) = CPLUSNS
 278       XPDF(4) = DPLUSNS
 279       XPDF(5) = UPLUSNS
 280       XPDF(6) = SPLUSNS
 281       XPDF(1) = SING
 282 C
 283       ADD = XPDF(1)/FLAV
 284       UPLUS=XPDF(5)+ADD
 285       DPLUS=-XPDF(4)+ADD
 286       SPLUS=-XPDF(6)+ADD
 287       CPLUS=-XPDF(3)+ADD
 288       UB=UPLUS*0.5
 289       UV=UB
 290       DB=DPLUS*0.5
 291       DV=DB
 292       SB=SPLUS*0.5
 293       CB=CPLUS*0.5
 294       SING=XPDF(1)
 295       GLU=XPDF(7)
 296       GL=GLU
 297 C... get parton density with massive charm
 298       CPLUM=AFCPLU(X,Q2)
 299       CB=CPLUM*0.5
 300 *
 301       DUV=MAX(ZERO,UV)
 302       DDV=MAX(ZERO,DV)
 303       DUB=MAX(ZERO,UB)
 304       DDB=MAX(ZERO,DB)
 305       DSB=MAX(ZERO,SB)
 306       DCB=MAX(ZERO,CB)
 307       DGL=MAX(ZERO,GL)
 308 C
 309       RETURN
 310       END
 311 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 312       SUBROUTINE SFAFG1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
 313 C
 314 C**************************************************************************
 315 C                       ( 1st of February 1994)
 316 C     This is an interpolation program which reads the files GRPOL and
 317 C     GRVDM and gives the quark and gluon distributions in real photon
 318 C     as functions of x and Q**2.
 319 C
 320 C     The Q**2 evolution is a BLL evolution (MSbar scheme) with Nf=4
 321 C     and LAMBDA(MSbar)=.200 Gev.
 322 C
 323 C     A massless charm distribution is generated for Q**2 > 2 Gev**2.
 324 C
 325 C     The distributions are the sum of a pointlike part (PL) and of a
 326 C     Vdm part (VDM):
 327 C                     dist=PL + KA*VDM
 328 C     KA is a factor which can be adjusted ( the default value is KA=1.0).
 329 C     The file GRPOL contains the pointlike part of the distributions.
 330 C     The file GRVDM contains the vdm part (A precise definition of this
 331 C     latter is given in the paper "PARTON DISTRIBUTIONS IN THE PHOTON",
 332 C     Preprint LPTHE Orsay 93-37, by P.Aurenche,M.Fontannaz and J.Ph.Guillet).
 333 C
 334 C     The output of the program is written in the file GETOUT with the
 335 C     following conventions
 336 C                              UPLUS=x(u+ubar)
 337 C                              DPLUS=x(d+dbar)
 338 C                              SPLUS=x(s+sbar)
 339 C                              CPLUS=x(c+cbar)
 340 C                              SING =UPLUS+DPLUS+SPLUS+CPLUS
 341 C                              GLU  =x*g
 342 C
 343 C      The interpolation is valid for     2. < Q**2 < 5.5E+5 Gev**2,
 344 C                             and for   .0015<  x   < .99
 345 C
 346 C      The program also gives the structure function F2:
 347 C                        F2 = q*Cq + g*Cg + Cgam
 348 C      Cq and Cg are the Wilson coeficients and Cgam is the direct term.
 349 C
 350 C      Although the charm quark evolution is massless, the direct term
 351 C      Cgam includes the effects due to the charm quark mass. The charm
 352 C      quark threshold is therefore correctly described at the lowest
 353 C      ordre in alphastrong (Details are given in the preprint).
 354 C
 355 C
 356 C**************************************************************************
 357 C
 358       double precision
 359      +       DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
 360       REAL    X, Q, UV, DV, UB, DB, SB, CB, BB, GL
 361       REAL       Q2
 362       DIMENSION XPDF(7)
 363       common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
 364       DIMENSION CALCO(8,20,32)
 365       DIMENSION CELCO(8,20,32)
 366       COMMON/W5051IA/CALCO
 367       COMMON/W5051IB/CELCO
 368       EXTERNAL AFCPLU
 369       DATA ZERO/0.0/
 370 C----------------------------------------------------------------------
 371        DATA ISTART/0/
 372        SAVE ISTART,OWLAM2,Q02,FLAV,KA, /W5051IA/, /W5051IB/
 373 C
 374       IF (ISTART.EQ.0) THEN
 375         ISTART=1
 376         DO 10 K=1,32
 377         DO 10 I=1,20
 378         DO 10 M=1,8
 379         CALCO(M,I,K) = CALC(M,I,K,3)
 380    10   CELCO(M,I,K) = CELC(M,I,K)
 381            OWLAM=PAR(1,3)
 382            OWLAM2=OWLAM**2
 383            Q02=PAR(30,3)
 384            FLAV=PAR(25,3)
 385            DELTA=PAR(29,3)
 386            CALL WATE32
 387            KA=1.0
 388          ENDIF
 389 C
 390       X = DX
 391       Q = DQ
 392       Q2 = Q*Q
 393       IDQ2=2
 394       SB=0.
 395       IF(Q2-Q02) 1,1,2
 396     2 IF(IDQ2-1) 1,1,3
 397     3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
 398     1 CONTINUE
 399       CALL AFGINT(8,0,X,SB,XPDF(7))
 400       CALL AFGINT(7,0,X,SB,SING)
 401       CALL AFGINT(4,0,X,SB,DPLUSNS)
 402       CALL AFGINT(3,0,X,SB,CPLUSNS)
 403       CALL AFGINT(5,0,X,SB,UPLUSNS)
 404       CALL AFGINT(6,0,X,SB,SPLUSNS)
 405       XPDF(3) = CPLUSNS
 406       XPDF(4) = DPLUSNS
 407       XPDF(5) = UPLUSNS
 408       XPDF(6) = SPLUSNS
 409       XPDF(1) = SING
 410 C
 411       ADD = XPDF(1)/FLAV
 412       UPLUS= XPDF(5)+ADD
 413       DPLUS=-XPDF(4)+ADD
 414       SPLUS=-XPDF(6)+ADD
 415       CPLUS=-XPDF(3)+ADD
 416       SING=XPDF(1)
 417       GLU=XPDF(7)
 418       GL=GLU
 419 *
 420       CALL AFGIN2(8,0,X,SB,XPDF(7))
 421       CALL AFGIN2(7,0,X,SB,SING)
 422       CALL AFGIN2(4,0,X,SB,DPLUSNS)
 423       CALL AFGIN2(3,0,X,SB,CPLUSNS)
 424       CALL AFGIN2(5,0,X,SB,UPLUSNS)
 425       CALL AFGIN2(6,0,X,SB,SPLUSNS)
 426       XPDF(3) = CPLUSNS
 427       XPDF(4) = DPLUSNS
 428       XPDF(5) = UPLUSNS
 429       XPDF(6) = SPLUSNS
 430       XPDF(1) = SING
 431 C
 432       ADD2 = XPDF(1)/FLAV
 433       UPLU2= XPDF(5)+ADD2
 434       DPLU2=-XPDF(4)+ADD2
 435       SPLU2=-XPDF(6)+ADD2
 436       CPLU2=-XPDF(3)+ADD2
 437       SING2=XPDF(1)
 438       GLU2=XPDF(7)
 439       UB=UPLUS+UPLU2*KA
 440       UV=UB
 441       DB=DPLUS+DPLU2*KA
 442       DV=DB
 443       SB=SPLUS+SPLU2*KA
 444       CB=CPLUS+CPLU2*KA
 445       SING=SING+SING2*KA
 446       GL=GLU+GLU2*KA
 447 *
 448       DUV=MAX(ZERO,UV)
 449       DDV=MAX(ZERO,DV)
 450       DUB=MAX(ZERO,UB)
 451       DDB=MAX(ZERO,DB)
 452       DSB=MAX(ZERO,SB)
 453       DCB=MAX(ZERO,CB)
 454       DGL=MAX(ZERO,GL)
 455 C
 456       RETURN
 457       END
 458 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 459       SUBROUTINE WATE32
 460 C  32 POINT GAUSSIAN QUADRATURE ROUTINE
 461       double precision
 462      +       X(16),W(16)
 463       double precision
 464      +               XI(32),WI(32),XX(33)
 465       COMMON/W5051I9/XI,WI,XX,NTERMS
 466       NTERMS=32
 467       X(1)=0.048307665687738316235D0
 468       X(2)=0.144471961582796493485D0
 469       X(3)=0.239287362252137074545D0
 470       X(4)=0.331868602282127649780D0
 471       X(5)=0.421351276130635345364D0
 472       X(6)=0.506899908932229390024D0
 473       X(7)=0.587715757240762329041D0
 474       X(8)=0.663044266930215200975D0
 475       X(9)=0.732182118740289680387D0
 476       X(10)=0.794483795967942406963D0
 477       X(11)=0.849367613732569970134D0
 478       X(12)=0.896321155766052123965D0
 479       X(13)=0.934906075937739689171D0
 480       X(14)=0.964762255587506430774D0
 481       X(15)=0.985611511545268335400D0
 482       X(16)=0.997263861849481563545D0
 483       W(1)=0.096540088514727800567D0
 484       W(2)=0.095638720079274859419D0
 485       W(3)=0.093844399080804565639D0
 486       W(4)=0.091173878695763884713D0
 487       W(5)=0.087652093004403811143D0
 488       W(6)=0.083311924226946755222D0
 489       W(7)=0.078193895787070306472D0
 490       W(8)=0.072345794108848506225D0
 491       W(9)=0.065822222776361846838D0
 492       W(10)=0.058684093478535547145D0
 493       W(11)=0.050998059262376176196D0
 494       W(12)=0.042835898022226680657D0
 495       W(13)=0.034273862913021433103D0
 496       W(14)=0.025392065309262059456D0
 497       W(15)=0.016274394730905670605D0
 498       W(16)=0.007018610009470096600D0
 499       NTERMH = NTERMS/2
 500       DO 1 I=1,NTERMH
 501       XI(I)=-X(17-I)
 502       WI(I)=W(17-I)
 503       XI(I+16)=X(I)
 504       WI(I+16)=W(I)
 505     1 CONTINUE
 506       DO 2 I=1,NTERMS
 507     2 XX(I)=0.5D0*(XI(I)+1.0D0)
 508       XX(33)=1.0D0
 509       RETURN
 510       END
 511 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 512       SUBROUTINE AURGAM(I,NDRV,X,S,ANS)
 513       DIMENSION F1(32),F2(32),F3(32)
 514       DIMENSION AF(3),AS(3)
 515       DIMENSION CALCO(8,20,32)
 516       COMMON/W5051I7/CALCO
 517       DATA DELTA/0.8000E-01/
 518       ANS=0.
 519       IF(X.GT.0.9985) RETURN
 520       N=3
 521       IS=S/DELTA+1
 522       IF(IS.GE.17) IS=17
 523       IS1=IS+1
 524       IS2=IS1+1
 525       DO 1 L=1,32
 526       KL=L+32*NDRV
 527       F1(L)=CALCO(I,IS,KL)
 528       F2(L)=CALCO(I,IS1,KL)
 529       F3(L)=CALCO(I,IS2,KL)
 530     1 CONTINUE
 531       AF(1)=AFGETFV(X,F1)
 532       AF(2)=AFGETFV(X,F2)
 533       AF(3)=AFGETFV(X,F3)
 534       AS(1)=(IS-1)*DELTA
 535       AS(2)=AS(1)+DELTA
 536       AS(3)=AS(2)+DELTA
 537       CALL AFPOLIN(AS,AF,N,S,AANS,DY)
 538       ANS=AANS
 539       RETURN
 540       END
 541 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 542       SUBROUTINE AFGINT(I,NDRV,X,S,ANS)
 543       DIMENSION F1(32),F2(32),F3(32)
 544       DIMENSION AF(3),AS(3)
 545       DIMENSION CALCO(8,20,32)
 546       COMMON/W5051IA/CALCO
 547       DATA DELTA/0.8000E-01/
 548       ANS=0.
 549       IF(X.GT.0.9985) RETURN
 550       N=3
 551       IS=S/DELTA+1
 552 C     IF(IS.GE.17) IS=17
 553       IS1=IS+1
 554       IS2=IS1+1
 555       DO 1 L=1,32
 556       KL=L+32*NDRV
 557       F1(L)=CALCO(I,IS,KL)
 558       F2(L)=CALCO(I,IS1,KL)
 559       F3(L)=CALCO(I,IS2,KL)
 560     1 CONTINUE
 561       AF(1)=AFGETFV(X,F1)
 562       AF(2)=AFGETFV(X,F2)
 563       AF(3)=AFGETFV(X,F3)
 564       AS(1)=(IS-1)*DELTA
 565       AS(2)=AS(1)+DELTA
 566       AS(3)=AS(2)+DELTA
 567       CALL AFPOLIN(AS,AF,N,S,AANS,DY)
 568       ANS=AANS
 569       RETURN
 570       END
 571 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 572       SUBROUTINE AFGIN2(I,NDRV,X,S,ANS)
 573       DIMENSION F1(32),F2(32),F3(32)
 574       DIMENSION AF(3),AS(3)
 575       DIMENSION CELCO(8,20,32)
 576       COMMON/W5051IB/CELCO
 577       DATA DELTA/0.8000E-01/
 578       ANS=0.
 579       IF(X.GT.0.9985) RETURN
 580       N=3
 581       IS=S/DELTA+1
 582 C     IF(IS.GE.17) IS=17
 583       IS1=IS+1
 584       IS2=IS1+1
 585       DO 1 L=1,32
 586       KL=L+32*NDRV
 587       F1(L)=CELCO(I,IS,KL)
 588       F2(L)=CELCO(I,IS1,KL)
 589       F3(L)=CELCO(I,IS2,KL)
 590     1 CONTINUE
 591       AF(1)=AFGETFV(X,F1)
 592       AF(2)=AFGETFV(X,F2)
 593       AF(3)=AFGETFV(X,F3)
 594       AS(1)=(IS-1)*DELTA
 595       AS(2)=AS(1)+DELTA
 596       AS(3)=AS(2)+DELTA
 597       CALL AFPOLIN(AS,AF,N,S,AANS,DY)
 598       ANS=AANS
 599       RETURN
 600       END
 601 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 602       FUNCTION AFCPLU(X,Q2)
 603       CMS=1.5**2
 604       BETS=1-4.*CMS*X/(1.-X)/Q2
 605       IF(BETS.LE..0) THEN
 606          AFCPLU=.0
 607          RETURN
 608       ENDIF
 609       BETA=SQRT(BETS)
 610       CPLU=(8.*X*(1.-X)-1.-4.*CMS*X*(1.-X)/Q2)*BETA
 611       CAU=X**2+(1.-X)**2+4.*CMS*X*(1.-3.*X)/Q2-8.*CMS**2*X**2/Q2**2
 612       CPLU=CPLU+CAU* LOG((1.+BETA)/(1.-BETA))
 613       AFCPLU=3.*(4./9.)*CPLU*X/(3.1415*137.)
 614   1   RETURN
 615       END
 616 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 617        FUNCTION AFGETFV(X,FVL)
 618 C  NOUVEAU PROGRAMME D'INTERPOLATION UTILISANT UNE ROUTINE DE MATH. RECIPES
 619        DIMENSION FVL(32)
 620        double precision
 621      +                XI(32),WI(32),XX(33)
 622        COMMON/W5051I9/XI,WI,XX,NTERMS
 623        DIMENSION A(4),B(4)
 624        N=4
 625        EPS=1.E-7
 626        XAM=XX(1)-EPS
 627        XAP=XX(1)+EPS
 628 C      IF(X.LT.XAM) PRINT*,' X = ',X
 629        IF(X.GT.XAM.AND.X.LT.XAP) GO TO 50
 630        GO TO 80
 631    50  Y=FVL(1)
 632        GO TO 77
 633    80  IF(X.LT.XX(2)) GO TO 51
 634        IF(X.GT.XX(30)) GO TO 61
 635        DO 1 I=3,30
 636        IF(X.GT.XX(I)) GO TO 1
 637        A(1)=XX(I-2)
 638        A(2)=XX(I-1)
 639        A(3)=XX(I)
 640        A(4)=XX(I+1)
 641        B(1)=FVL(I-2)
 642        B(2)=FVL(I-1)
 643        B(3)=FVL(I)
 644        B(4)=FVL(I+1)
 645        GO TO 70
 646    1   CONTINUE
 647   61   A(1)=XX(29)
 648        A(2)=XX(30)
 649        A(3)=XX(31)
 650        A(4)=XX(32)
 651        B(1)=FVL(29)
 652        B(2)=FVL(30)
 653        B(3)=FVL(31)
 654        B(4)=FVL(32)
 655        GO  TO 70
 656   51   A(1)=XX(1)
 657        A(2)=XX(2)
 658        A(3)=XX(3)
 659        A(4)=XX(4)
 660        B(1)=FVL(1)
 661        B(2)=FVL(2)
 662        B(3)=FVL(3)
 663        B(4)=FVL(4)
 664 C 70   IF(X.GT..2.AND.X.LT..8) THEN
 665 C            CALL AFPOLIN(A,B,N,X,Y,DY)
 666 C      ELSE
 667 C            CALL AFRATIN(A,B,N,X,Y,DY)
 668 C      ENDIF
 669   70   CONTINUE
 670              CALL AFPOLIN(A,B,N,X,Y,DY)
 671   77   AFGETFV=Y
 672        RETURN
 673        END
 674 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 675       SUBROUTINE AFPOLIN(XA,YA,N,X,Y,DY)
 676       PARAMETER (NMAX=10)
 677       DIMENSION XA(NMAX),YA(NMAX),C(NMAX),D(NMAX)
 678       Y=0.
 679       IF(N.GT.NMAX) RETURN
 680       NS=1
 681       DIF=ABS(X-XA(1))
 682       DO 11 I=1,N
 683         DIFT=ABS(X-XA(I))
 684         IF (DIFT.LT.DIF) THEN
 685           NS=I
 686           DIF=DIFT
 687         ENDIF
 688         C(I)=YA(I)
 689         D(I)=YA(I)
 690 11    CONTINUE
 691       Y=YA(NS)
 692       NS=NS-1
 693       DO 13 M=1,N-1
 694         DO 12 I=1,N-M
 695           HO=XA(I)-X
 696           HP=XA(I+M)-X
 697           W=C(I+1)-D(I)
 698           DEN=HO-HP
 699 C         IF(DEN.EQ.0.)PAUSE
 700           DEN=W/DEN
 701           D(I)=HP*DEN
 702           C(I)=HO*DEN
 703 12      CONTINUE
 704         IF (2*NS.LT.N-M)THEN
 705           DY=C(NS+1)
 706         ELSE
 707           DY=D(NS)
 708           NS=NS-1
 709         ENDIF
 710         Y=Y+DY
 711 13    CONTINUE
 712       RETURN
 713       END
 714 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc