LHAPDF/lhapdf5.3.1/inputPDF.f

   1
   2       subroutine weightPDF(x)
   3       implicit none
   4       integer nset,imem,nfmax
   5       real*8 x,Q
   6       nset=1
   7       call weightPDFM(nset,x)
   8       return
   9
  10       entry GetNF(nfmax)
  11       nset=1
  12       call GetNFM(nset,nfmax)
  13       return
  14
  15       entry GetThreshold(imem,Q)
  16       nset=1
  17       call GetThresholdM(nset,imem,Q)
  18       return
  19
  20       end
  21
  22       subroutine parmPDF(nset,x,pdf)
  23       implicit none
  24       include 'parmsetup.inc'
  25       integer nset
  26       character*16 name(nmxset)
  27       integer nmem(nmxset),ndef(nmxset),mem
  28       common/NAME/name,nmem,ndef,mem
  29       real*4 xp(40),x4,fgdis4
  30       character*16 s1,s2
  31       integer i,j,imem,nop,parmN,Nfunc(nmxset),Fw(nmxset),nfmax,M
  32       real*8 x,N,b0,Poly,pdf(-6:6),Fparm(nopmax),F(nofmax)
  33       real*8 Ccoef(nmxset,-6:6,nofmax),Fpow(nmxset,nofmax)
  34       real*8 Q,Treshold(nmxset,-6:6)
  35 c      data Treshold/39*0d0/
  36       integer Fmap(nmxset,nofmax,npfmax)
  37       integer Ftype(nmxset,nofmax),Fn(nmxset,nofmax),Ctype(nmxset,-6:6)
  38       integer lhasilent
  39       common/lhasilent/lhasilent
  40       logical first
  41       data first/.true./
  42       save Nfunc,Fn,Fw,Fpow,Fmap,Ccoef,Fparm,Ftype,Ctype,first,Treshold
  43
  44       do i=1,Nfunc(nset)
  45          if (Ftype(nset,i).eq.1) then
  46             Poly=1.0
  47
  48             do j=4,Fn(nset,i)
  49          Poly=Poly+Fparm(Fmap(nset,i,j))*x**(float(j-3)/Fpow(nset,i))
  50             enddo
  51             Poly=Fparm(Fmap(nset,i,1))*Poly
  52        F(i)=x**Fparm(Fmap(nset,i,2))*(1.0-x)**Fparm(Fmap(nset,i,3))*Poly
  53          endif
  54          if (Ftype(nset,i).eq.2) then
  55             if (x.lt.0.9999999) then
  56       Poly=Fparm(Fmap(nset,i,2))*log(x)+Fparm(Fmap(nset,i,3))*log(1.0-x)
  57      .      +Fparm(Fmap(nset,i,4))*x
  58      .      +Fparm(Fmap(nset,i,6))*log(1.0+x*exp(Fparm(Fmap(nset,i,5))))
  59          F(i)=Fparm(Fmap(nset,i,1))*exp(Poly)
  60             else
  61                F(i)=0d0
  62             endif
  63          endif
  64          if (Ftype(nset,i).eq.101) then
  65        Poly=exp(Fparm(Fmap(nset,i,1)))*x**(Fparm(Fmap(nset,i,2))-1)
  66      .       *(1d0-x)**Fparm(Fmap(nset,i,3))
  67        Poly=
  68      . Poly+(1d0+Fparm(Fmap(nset,i,4))*x)*(1d0-x)**Fparm(Fmap(nset,i,5))
  69             b0=10d0
  70             if (Poly.gt.b0) then
  71                F(i)=Poly
  72             elseif (Poly.lt.-b0) then
  73                F(i)=0d0
  74             else
  75                F(i)=Poly+log(1d0+exp(-b0*Poly)-exp(-b0))/b0
  76             endif
  77          endif
  78 c - to add the mrst2004 gluon convolution
  79          if (Ftype(nset,i).eq.201) then
  80            xp(2) = Fparm(Fmap(nset,i,1))
  81            xp(3) = Fparm(Fmap(nset,i,2))
  82            xp(23)= Fparm(Fmap(nset,i,3))
  83            xp(16)= Fparm(Fmap(nset,i,4))
  84            xp(5) = Fparm(Fmap(nset,i,5))
  85            xp(40)= Fparm(Fmap(nset,i,6))
  86            xp(24)= Fparm(Fmap(nset,i,7))
  87            xp(20)= Fparm(Fmap(nset,i,8))
  88            xp(4) = Fparm(Fmap(nset,i,9))
  89            xp(1) = Fparm(Fmap(nset,i,10))
  90            x4 = x
  91            call gconv(x4,xp,fgdis4)
  92           F(i) = FGDIS4
  93          endif
  94 c--
  95       enddo
  96       do i=-6,6
  97          pdf(i)=0.0
  98          if (Ctype(nset,i).gt.0) then
  99             if (Ctype(nset,i).eq.1) then
 100                do j=1,Nfunc(nset)
 101                   pdf(i)=pdf(i)+Ccoef(nset,i,j)*F(j)
 102 c                 print *,i,j,Ccoef(i,j),F(j)
 103                enddo
 104             endif
 105             if (Ctype(nset,i).eq.101) then
 106                if (i.eq.-2) then
 107        pdf(i)=F(int(Ccoef(nset,i,1)))/(F(int(Ccoef(nset,i,2)))+1d0)
 108                endif
 109                if (i.eq.-1) then
 110       pdf(i)=
 111      .F(int(Ccoef(nset,i,1)))/(1d0/F(int(Ccoef(nset,i,2)))+1d0)
 112                endif
 113                if (i.eq.1) then
 114                   pdf(i)=F(1)+pdf(-1)
 115                endif
 116                if (i.eq.2) then
 117                   pdf(i)=F(2)+pdf(-2)
 118                endif
 119             endif
 120          endif
 121       enddo
 122 c      print *,pdf
 123       return
 124 *
 125       entry weightPDFM(nset,x)
 126       if (Fw(nset).ge.0) then
 127          x=Fparm(Fw(nset))
 128       else
 129          call numberPDF(nop)
 130          x=1.0/float(nop)
 131       endif
 132       return
 133 *
 134 c      entry GetParmPDFM(nset,imem,x)
 135       entry GetParmPDF(nset,imem,x)
 136       x=Fparm(imem)
 137       return
 138 *
 139       entry GetNfM(nset,nfmax)
 140       if(name(nset).eq.'MRST'.or.
 141      +   name(nset).eq.'EVLCTEQ'.or.
 142      +   name(nset)(1:6).eq.'QCDNUM') then
 143         nfmax=0
 144         do i=1,6
 145            if (Treshold(nset,-i).ge.0d0) nfmax=nfmax+1
 146            if (Treshold(nset,i).ge.0d0) nfmax=nfmax+1
 147         enddo
 148         nfmax=nfmax/2
 149       else
 150         nfmax = -1
 151       endif
 152       return
 153 *
 154       entry GetThresholdM(nset,imem,Q)
 155       Q=Treshold(nset,imem)
 156       return
 157 *
 158 c      entry InitEvolvePDFM(nset,imem)
 159       entry InitEvolvePDF(nset,imem)
 160 c      print*, 'calling listPDF',nset,imem
 161       call listPDF(nset,imem,Fparm)
 162 c      print *,Fparm
 163       return
 164 *
 165 c      entry initInputPDFM(nset)
 166       entry initInputPDF(nset)
 167
 168       if(first) then
 169         do i=1,nmxset
 170           do j=-6,6
 171              Treshold(i,j)=0.0d0
 172           enddo
 173         enddo
 174         first=.false.
 175       endif
 176       read(1,*) s1,Fw(nset),Nfunc(nset)
 177 c      print *,s1,Fw,Nfunc
 178       if(lhasilent.eq.0) then
 179         write(*,*) 'Parametrization: ',s1
 180         write(*,*)
 181       endif
 182       do i=1,Nfunc(nset)
 183          Ftype(nset,i)=-1
 184          read(1,*) s1,s2
 185          if (index(s2,'x-taylor').eq.1) then
 186             Ftype(nset,i)=1
 187             read(1,*) FPow(nset,i),Fn(nset,i)
 188             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
 189          endif
 190          if (index(s2,'log-pade').eq.1) then
 191             Ftype(nset,i)=2
 192             FPow(nset,i)=0d0
 193             read(1,*) Fn(nset,i)
 194             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
 195          endif
 196          if (index(s2,'cteq6-ratio').eq.1) then
 197             Ftype(nset,i)=101
 198             Fpow(nset,i)=0d0
 199             Fn(nset,i)=5
 200             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
 201          endif
 202          if (index(s2,'convol').eq.1) then
 203             Ftype(nset,i)=201
 204             read(1,*) Fn(nset,i)
 205             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
 206          endif
 207          if (Ftype(nset,i).lt.0) then
 208             write(*,*) 'File description error:'
 209             write(*,*) 'Unknown functional ',s2
 210             stop
 211          endif
 212       enddo
 213       read(1,*) s1
 214       do i=-6,6
 215          Ctype(nset,i)=-1
 216          read(1,*) s1,s2
 217 c        print *,s1,s2
 218          if (index(s2,'none').eq.1) then
 219             Ctype(nset,i)=0
 220             Treshold(nset,i)=-1d0
 221          endif
 222          if (index(s2,'treshold').eq.1) then
 223             Ctype(nset,i)=0
 224             read(1,*) Treshold(nset,i)
 225          endif
 226          if (index(s2,'composite').eq.1) then
 227             Ctype(nset,i)=1
 228             Treshold(nset,i)=0d0
 229             read(1,*) (Ccoef(nset,i,j),j=1,Nfunc(nset))
 230          endif
 231          if (index(s2,'cteq6-ratio').eq.1) then
 232             Ctype(nset,i)=101
 233             Treshold(nset,i)=0d0
 234             read(1,*) (Ccoef(nset,i,j),j=1,3)
 235          endif
 236          if (Ctype(nset,i).lt.0) then
 237             write(*,*) 'File description error:'
 238             write(*,*) 'Unknown composit type ',s2
 239             stop
 240          endif
 241       enddo
 242       if (Fw(nset).ge.0) then
 243          write(*,*) '***********************************************'
 244          write(*,*) '* Note that this is a weigthed PDF set.       *'
 245          write(*,*) '* See manual for proper use.                  *'
 246          write(*,*) '***********************************************'
 247       endif
 248       return
 249 *
 250       end
 251
 252       FUNCTION BETA_LHA(X1,X2)
 253       CALL GAMMA_LHA(X1,G1,IER)
 254       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
 255       CALL GAMMA_LHA(X2,G2,IER)
 256       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
 257       X3=X1+X2
 258       CALL GAMMA_LHA(X3,G3,IER)
 259       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
 260       BETA_LHA=G1*G2/G3
 261       RETURN
 262       END
 263
 264       SUBROUTINE GAMMA_LHA(XX,GX,IER)
 265       IF(XX-34.5)6,6,4
 266     4 IER=2
 267       GX=1.E38
 268       RETURN
 269     6 X=XX
 270       ERR=1.0E-6
 271       IER=0
 272       GX=1.0
 273       IF(X-2.0)50,50,15
 274    10 IF(X-2.0)110,110,15
 275    15 X=X-1.0
 276       GX=GX*X
 277       GO TO 10
 278    50 IF(X-1.0)60,120,110
 279 C        SEE IF X IS NEAR NEGATIVE INTEGER OR ZERO
 280    60 IF(X-ERR)62,62,80
 281    62 K=X
 282       Y=FLOAT(K)-X
 283       IF(ABS(Y)-ERR)130,130,64
 284    64 IF(1.0-Y-ERR)130,130,70
 285 C        X NOT NEAR A NEGATIVE INTEGER OR ZERO
 286    70 IF(X-1.0)80,80,110
 287    80 GX=GX/X
 288       X=X+1.0
 289       GO TO 70
 290   110 Y=X-1.0
 291       GY=1.0+Y*(-0.5771017+Y*(+0.9858540+Y*(-0.8764218+Y*(+0.8328212+
 292      1Y*(-0.5684729+Y*(+0.2548205+Y*(-0.05149930)))))))
 293       GX=GX*GY
 294   120 RETURN
 295   130 IER=1
 296       RETURN
 297       END
 298
 299       FUNCTION  ALPHA(T,AL)
 300       COMMON/AINPUT/IORD,QSCT,QSDT
 301       COMMON/PARAM/PARA(40)
 302       DATA PI/3.14159/
 303       DATA TOL/.0005/
 304       ITH=0
 305       TT=T
 306       qsctt=qsct/4.
 307       qsdtt=qsdt/4.
 308 c      AL=para(1)
 309       AL2=AL*AL
 310       FLAV=4.
 311       QS=AL2*EXP(T)
 312
 313       if(qs.lt.0.5d0) then   !!  running stops below 0.5
 314           qs=0.5d0
 315           t=alog(qs/al2)
 316           tt=t
 317       endif
 318
 319       IF(QS.gt.QSCTT) GO        TO 12
 320       IF(QS.lt.QSDTT) GO        TO 312
 321    11 CONTINUE
 322       B0=11-2.*FLAV/3.
 323       IF(IORD)1,1,2
 324 c     IF(IORD)2,2,2     !TAKE CARE !!
 325     1 CONTINUE
 326       ALPHA=4.*PI/B0/T
 327       RETURN
 328     2 CONTINUE
 329       X1=4.*PI/B0
 330       B1=102.-38.*FLAV/3.
 331       X2=B1/B0**2
 332       AS=X1/T*(1.-X2*aLOG(T)/T)
 333     5 CONTINUE
 334       F=-T+X1/AS-X2*aLOG(X1/AS+X2)
 335       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
 336       AS2=AS-F/FP
 337       DEL=ABS(F/FP/AS)
 338       IF(DEL-TOL)3,3,4
 339     3 CONTINUE
 340       ALPHA=AS2
 341       IF(ITH.EQ.0) RETURN
 342       GO TO (13,14,15) ITH
 343     4 CONTINUE
 344       AS=AS2
 345       GO TO 5
 346    12 ITH=1
 347       T=aLOG(QSCTT/AL2)
 348       GO TO 11
 349    13 ALFQC4=ALPHA
 350       FLAV=5.
 351       ITH=2
 352       GO TO 11
 353    14 ALFQC5=ALPHA
 354       ITH=3
 355       T=TT
 356       GO TO 11
 357    15 ALFQS5=ALPHA
 358       ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
 359       ALPHA=1./ALFINV
 360       RETURN
 361
 362   311 CONTINUE
 363       B0=11-2.*FLAV/3.
 364       IF(IORD)31,31,32
 365 c     IF(IORD)32,32,32  !TAKE CARE !!
 366    31 CONTINUE
 367       ALPHA=4.*PI/B0/T
 368       RETURN
 369    32 CONTINUE
 370       X1=4.*PI/B0
 371       B1=102.-38.*FLAV/3.
 372       X2=B1/B0**2
 373       AS=X1/T*(1.-X2*aLOG(T)/T)
 374    35 CONTINUE
 375       F=-T+X1/AS-X2*aLOG(X1/AS+X2)
 376       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
 377       AS2=AS-F/FP
 378       DEL=ABS(F/FP/AS)
 379       IF(DEL-TOL)33,33,34
 380    33 CONTINUE
 381       ALPHA=AS2
 382       IF(ITH.EQ.0) RETURN
 383       GO TO (313,314,315) ITH
 384    34 CONTINUE
 385       AS=AS2
 386       GO TO 35
 387   312 ITH=1
 388       T=aLOG(QSDTT/AL2)
 389       GO TO 311
 390   313 ALFQC4=ALPHA
 391       FLAV=3.
 392       ITH=2
 393       GO TO 311
 394   314 ALFQC3=ALPHA
 395       ITH=3
 396       T=TT
 397       GO TO 311
 398   315 ALFQS3=ALPHA
 399       ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
 400       ALPHA=1./ALFINV
 401       RETURN
 402       END
 403
 404       SUBROUTINE WATE96
 405 C*******************************************************************
 406 C*****                                                         *****
 407 C***** THE X(I) AND W(I) ARE THE DIRECT OUTPUT FROM A PROGRAM  *****
 408 C***** USING NAG ROUTINE D01BCF TO CALCULATE THE               *****
 409 C***** GAUSS-LEGENDRE WEIGHTS FOR 96 POINT INTEGRATION.        *****
 410 C***** THEY AGREE TO TYPICALLY 14 DECIMAL PLACES WITH THE      *****
 411 C***** TABLE IN ABRAMOWITZ & STEGUN, PAGE 919.                 *****
 412 C*****                                                         *****
 413 C***** ---->   PETER HARRIMAN, APRIL 3RD 1990.                 *****
 414 C*****                                                         *****
 415 C*******************************************************************
 416       DIMENSION X(48),W(48)
 417       COMMON/GAUS96/XI(96),WI(96),nterms,XX(97)
 418       NTERMS=96
 419
 420       X( 1)=   0.01627674484960183561
 421       X( 2)=   0.04881298513604856015
 422       X( 3)=   0.08129749546442434360
 423       X( 4)=   0.11369585011066471632
 424       X( 5)=   0.14597371465489567682
 425       X( 6)=   0.17809688236761733390
 426       X( 7)=   0.21003131046056591064
 427       X( 8)=   0.24174315616383866556
 428       X( 9)=   0.27319881259104774468
 429       X(10)=   0.30436494435449495954
 430       X(11)=   0.33520852289262397655
 431       X(12)=   0.36569686147231213885
 432       X(13)=   0.39579764982890709712
 433       X(14)=   0.42547898840729897474
 434       X(15)=   0.45470942216774136446
 435       X(16)=   0.48345797392059470382
 436       X(17)=   0.51169417715466604391
 437       X(18)=   0.53938810832435567233
 438       X(19)=   0.56651041856139533470
 439       X(20)=   0.59303236477757022282
 440       X(21)=   0.61892584012546672523
 441       X(22)=   0.64416340378496526886
 442       X(23)=   0.66871831004391424358
 443       X(24)=   0.69256453664216964528
 444       X(25)=   0.71567681234896561582
 445       X(26)=   0.73803064374439816819
 446       X(27)=   0.75960234117664555964
 447       X(28)=   0.78036904386743123629
 448       X(29)=   0.80030874413913884180
 449       X(30)=   0.81940031073792957139
 450       X(31)=   0.83762351122818502758
 451       X(32)=   0.85495903343459936363
 452       X(33)=   0.87138850590929436968
 453       X(34)=   0.88689451740241818933
 454       X(35)=   0.90146063531585023110
 455       X(36)=   0.91507142312089592706
 456       X(37)=   0.92771245672230655266
 457       X(38)=   0.93937033975275308073
 458       X(39)=   0.95003271778443564022
 459       X(40)=   0.95968829144874048809
 460       X(41)=   0.96832682846326217918
 461       X(42)=   0.97593917458513455843
 462       X(43)=   0.98251726356301274934
 463       X(44)=   0.98805412632962202890
 464       X(45)=   0.99254390032376081654
 465       X(46)=   0.99598184298720747465
 466       X(47)=   0.99836437586317963722
 467       X(48)=   0.99968950388322870559
 468       W( 1)=   0.03255061449236316962
 469       W( 2)=   0.03251611871386883307
 470       W( 3)=   0.03244716371406427668
 471       W( 4)=   0.03234382256857594104
 472       W( 5)=   0.03220620479403026124
 473       W( 6)=   0.03203445623199267876
 474       W( 7)=   0.03182875889441101874
 475       W( 8)=   0.03158933077072719007
 476       W( 9)=   0.03131642559686137819
 477       W(10)=   0.03101033258631386231
 478       W(11)=   0.03067137612366917839
 479       W(12)=   0.03029991542082762553
 480       W(13)=   0.02989634413632842385
 481       W(14)=   0.02946108995816795100
 482       W(15)=   0.02899461415055528410
 483       W(16)=   0.02849741106508543861
 484       W(17)=   0.02797000761684838950
 485       W(18)=   0.02741296272602931385
 486       W(19)=   0.02682686672559184485
 487       W(20)=   0.02621234073567250055
 488       W(21)=   0.02557003600534944960
 489       W(22)=   0.02490063322248370695
 490       W(23)=   0.02420484179236479915
 491       W(24)=   0.02348339908592633665
 492       W(25)=   0.02273706965832950717
 493       W(26)=   0.02196664443874448477
 494       W(27)=   0.02117293989219144572
 495       W(28)=   0.02035679715433347898
 496       W(29)=   0.01951908114014518992
 497       W(30)=   0.01866067962741165898
 498       W(31)=   0.01778250231604547316
 499       W(32)=   0.01688547986424539715
 500       W(33)=   0.01597056290256253144
 501       W(34)=   0.01503872102699521608
 502       W(35)=   0.01409094177231515264
 503       W(36)=   0.01312822956696188190
 504       W(37)=   0.01215160467108866759
 505       W(38)=   0.01116210209983888144
 506       W(39)=   0.01016077053500880978
 507       W(40)=   0.00914867123078384552
 508       W(41)=   0.00812687692569928101
 509       W(42)=   0.00709647079115442616
 510       W(43)=   0.00605854550423662775
 511       W(44)=   0.00501420274292825661
 512       W(45)=   0.00396455433844564804
 513       W(46)=   0.00291073181793626202
 514       W(47)=   0.00185396078894924657
 515       W(48)=   0.00079679206555731759
 516       DO 1 I=1,48
 517       XI(I)=-X(49-I)
 518       WI(I)=W(49-I)
 519       XI(I+48)=X(I)
 520       WI(I+48)=W(I)
 521     1 CONTINUE
 522       DO 2 I=1,96
 523     2 XX(I)=0.5*(XI(I)+1.)
 524       XX(97)=1.0
 525       EXPON=1.0
 526       DO 3 I=1,96
 527       YI=2.*(0.5*(1.+XI(I)))**EXPON-1.
 528       WI(I)=WI(I)/(1.+XI(I))*(1.+YI)*EXPON
 529       XI(I)=YI
 530       XX(I)=0.5*(1.+YI)
 531     3 CONTINUE
 532       RETURN
 533       END
 534
 535       subroutine gconv(x,xp,fgdis)
 536       COMMON/AINPUT/IORD,QSCT,QSDT
 537       common/GAUS96/XI(96),WI(96),NTERMS,XX(97)
 538       dimension xp(40)
 539       logical first
 540       data first/.true./
 541       if(first) then
 542          call wate96
 543          first=.false.
 544       endif
 545       PI = 3.14159
 546       PI2 = PI*PI
 547       iord = 1
 548       qsdt=8.18    !!  This is the value of 4m_c^2
 549       qsct=74.0    !!  This is the value of 4m_b^2
 550       cf = 4./3.
 551       eta4 = xp(40)
 552       T=alog(1/xp(1)**2)
 553 c      AL = 0.550/(4.* 3.14159)
 554       AL=ALPHA(T,xp(1))/(4.* pi)
 555       rx=sqrt(x)
 556       FF1=BETA_LHA(XP(2),XP(3)+1.)+XP(16)*BETA_LHA(XP(2)+1.,XP(3)+1.)+
 557      XXP(23)*BETA_LHA(XP(2)+0.5,XP(3)+1.)
 558       FF2=BETA_LHA(XP(2)+1.,XP(3)+1.)+
 559      XXP(16)*BETA_LHA(XP(2)+2.,XP(3)+1.)+
 560      XXP(23)*BETA_LHA(XP(2)+1.5,XP(3)+1.)
 561       FF3=BETA_LHA(XP(5),ETA4+1.)+XP(20)*BETA_LHA(XP(5)+1.,ETA4+1.)+
 562      XXP(24)*BETA_LHA(XP(5)+0.5,ETA4+1.)
 563       FF4=BETA_LHA(XP(5)+1.,ETA4+1.)+XP(20)*BETA_LHA(XP(5)+2.,ETA4+1.)+
 564      XXP(24)*BETA_LHA(XP(5)+1.5,ETA4+1.)
 565       COEFU=2.*XP(4)/FF1
 566       COEFD=   XP(4)/FF3
 567 c          print *,'coefu ',coefu
 568 c          print *,'coefd ',coefd
 569       UV=coefu*X**XP(2)*(1.-X)**XP(3)*(1.+XP(16)*X+XP(23)*SQRT(X))
 570       DV=coefd*X**XP(5)*(1.-X)**ETA4*(1.+XP(20)*X+XP(24)*SQRT(X))
 571       FGDIS=al*CF*(-9.-2.*PI2/3.+alog(1.-x)*(-3.+
 572      .2.*alog(1.-x)))*(UV+DV)
 573
 574       DO 23 M=1,NTERMS
 575       Y=0.5*(1.-X)*XI(M)+0.5*(1.+X)
 576       XY=X/Y
 577       UVXY=coefu*XY**XP(2)*(1.-XY)**XP(3)*(1.+XP(16)*XY
 578      .+XP(23)*SQRT(XY))
 579       DVXY=coefd*XY**XP(5)*(1.-XY)**ETA4*(1.+XP(20)*XY
 580      .+XP(24)*SQRT(XY))
 581       AL1=ALOG(1.-Y)
 582       C22=CF*(6.+4.*Y-2.*(1.+Y*Y)/(1.-Y)*ALOG(Y)-2.*(1.+Y)*ALOG(1.-Y))
 583       C23=CF*(-3.+4.*ALOG(1.-Y))/(1.-Y)
 584
 585       FGDIS=FGDIS+.5*(1.-X)*WI(M)*al*
 586      .(C22*(uvxy+dvxy)+C23*(uvxy+dvxy-uv-dv))
 587
 588    23 CONTINUE
 589
 590       return
 591       end