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