1 SUBROUTINE CtLhALFSET (QS, ALFS)
2 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
5 DATA ALAM, BLAM, ERR / 0.01, 10.0, 0.02 /
7 CALL CtLhParQcd (2, 'ORDR', ORDR, IR1)
10 EFLLN = CtLhQZBRNT (CtLhRTALF, ALAM, BLAM, ERR, IR2)
11 EFFLAM = QS / EXP (EFLLN)
12 CALL CtLhSETL1 (NEFF, EFFLAM)
14 FUNCTION CtLhALPI (AMU)
15 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
16 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
17 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
19 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
21 IF(.NOT.SET) CALL CtLhLAMCWZ
24 CtLhALPI = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT)
26 CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI', 'AMU', AMU,
28 ELSEIF (IRT .EQ. 2) THEN
29 CALL CtLhWARNR (IW2, 'CtLhALPI > 3; Be aware!', 'CtLhALPI',
30 > CtLhALPI, D0, D1, 0)
34 FUNCTION CtLhALPQCD (IRDR, NF, RML, IRT)
35 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
36 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
37 PARAMETER (CG = 3.0d0, TR = 0.5d0, CF = 4.0d0/3.0d0)
39 IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN
41 > 'Order out of range in CtLhALPQCD: IRDR = ', IRDR
44 B0 = (11.d0*CG - 2.* NF) / 3.d0
45 B1 = (34.d0*CG**2 - 10.d0*CG*NF - 6.d0*CF*NF) / 3.d0
54 IF (IRDR .GE. 2) AL = AL * (1.d0-B1*LOG(ALN) / ALN / B0**2)
61 FUNCTION CtLhAMHATF(I)
62 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
63 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
64 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
66 IF (.NOT.SET) CALL CtLhLAMCWZ
67 IF ((I.LE.0).OR.(I.GT.9)) THEN
68 print *,'warning I OUT OF RANGE IN CtLhAMHATF'
76 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
77 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
78 DATA HUGE, IWRN / 1.E20, 0 /
85 CALL CtLhWARNR(IWRN, 'CtLhDXDZ singular in CtLhDXDZ; set=HUGE',
92 SUBROUTINE CtLhEVLPAR (IACT, NAME, VALUE, IRET)
93 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
97 WRITE ( NINT(VALUE) , 101)
98 101 FORMAT (/ ' Initiation parameters: Qini, Ipd0, Ihdn ' /
99 > ' Maximum Q, Order of Alpha: Qmax, IKNL ' /
100 > ' X- mesh parameters : Xmin, Xcr, Nx ' /
101 > ' LnQ-mesh parameters : Nt, Jt ' /
102 > ' # of parton flavors : NfMx ' /)
104 ElseIF (IACT .EQ. 1) THEN
105 CALL CtLhEVLSET (NAME, VALUE, IRET)
107 print *,'fatal evlpar'
112 SUBROUTINE CtLhEVLSET (NAME, VALUE, IRET)
113 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
116 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
117 PARAMETER (MXPN = MXF * 2 + 2)
118 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
119 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
120 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
121 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
122 > / PdfSwh / Iset, IpdMod, Iptn0, NuIni
124 IF (NAME .EQ. 'QINI') THEN
125 IF (VALUE .LE. 0) GOTO 12
127 ElseIF (NAME .EQ. 'IPD0') THEN
129 IF (Item .Eq. 10 .or. Item .Eq. 11) GOTO 12
131 ElseIF (NAME .EQ. 'IHDN') THEN
133 IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12
135 ElseIF (NAME .EQ. 'QMAX') THEN
136 IF (VALUE .LE. QINI) GOTO 12
138 ElseIF (NAME .EQ. 'IKNL') THEN
141 IF (ITEM.NE.1.AND.ITEM.NE.2) GOTO 12
143 ElseIF (NAME .EQ. 'XCR') THEN
144 IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12
147 ElseIF (NAME .EQ. 'XMIN') THEN
148 IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12
151 ElseIF (NAME .EQ. 'NX') THEN
153 IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12
156 ElseIF (NAME .EQ. 'NT') THEN
158 IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12
160 ElseIF (NAME .EQ. 'JT') THEN
162 IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12
164 ElseIF (NAME .EQ. 'NFMX') THEN
166 IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12
168 ElseIF (NAME .EQ. 'IPDMOD') THEN
170 IF (Abs(Item) .Gt. 1) GOTO 12
172 ElseIF (NAME .EQ. 'IPTN0') THEN
174 IF (ABS(ITEM) .GT. MXF) GOTO 12
176 ElseIF (NAME .EQ. 'NUINI') THEN
178 IF (ITEM .LE. 0) GOTO 12
187 SUBROUTINE CtLhEVOLVE (FINI, IRET)
188 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
189 include 'parmsetup.inc'
191 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
192 PARAMETER (MXPN = MXF * 2 + 2)
193 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
194 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
195 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
196 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
197 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
198 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
199 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
200 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
201 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
203 DIMENSION JI(-MXF : MXF+1)
204 EXTERNAL LhCtNSRHSP, LhCtNSRHSM, FINI
206 save nxsave, ntsave, jtsave, ngsave,
207 & xcrsave, xminsave, qinisave, qmaxsave, ientry
214 IF (IHDN .LE. 4) THEN
216 ElseIF (IHDN .LE. 6) THEN
219 IF (.NOT. LSTX) CALL CtLhXARRAY
220 CALL CtLhPARPDF (2, 'ALAM', AL, IR)
221 CALL CtLhQARRAY (NINI)
224 Nelmt = KF * (Nt+1) * (Nx+1)
225 DO 101 IFLV = -NFMX, NFMX+1
227 JI(IFLV) = JFL * (NT+1) * (NX+1)
230 UPD(JI(0)+IZ+1,iset) = FINI (0, XV(IZ))
231 UPD(JI(NFSN)+IZ+1,iset) = 0
232 IF (NFMX .EQ. 0) GOTO 31
233 DO 331 IFLV = 1, NINI
234 A = FINI ( IFLV, XV(IZ))
235 B = FINI (-IFLV, XV(IZ))
237 UPD(JI(NFSN)+IZ+1,iset) =
238 > UPD(JI(NFSN)+IZ+1,iset) + QRKP (IFLV)
239 UPD(JI(-IFLV)+IZ+1,iset) = A - B
241 DO 332 IFLV = 1, NINI
242 UPD(JI( IFLV)+IZ+1,iset) =
243 > QRKP(IFLV) - UPD(JI(NFSN)+IZ+1,iset)/NINI
246 DO 21 NEFF = NINI, NFMX
247 IF (IKNL .EQ. 2) CALL CtLhSTUPKL (NEFF)
248 ICNT = NEFF - NINI + 1
249 IF (NTN(ICNT) .EQ. 0) GOTO 21
253 CALL CtLhSNEVL (IKNL, NX, NITR, JT, DT, TIN, NEFF,
254 > UPD(JI(NFSN)+2,iset), UPD(JI(0)+2,iset),
255 > UPD(JI(NFSN)+1,iset), UPD(JI(0)+1,iset))
256 IF (NEFF .EQ. 0) GOTO 88
257 DO 333 IFLV = 1, NEFF
258 CALL CtLhNSEVL (LhCtNSRHSP, IKNL, NX, NITR, JT, DT, TIN,
259 > NEFF, UPD(JI( IFLV)+2,iset), UPD(JI( IFLV)+1,iset))
261 > CALL CtLhNSEVL (LhCtNSRHSM, IKNL, NX, NITR, JT, DT, TIN,
262 > NEFF, UPD(JI(-IFLV)+2,iset), UPD(JI(-IFLV)+1,iset))
265 TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV),iset)
266 TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN),iset) / NEFF
268 IF (IKNL .GT. 0) TP = MAX (TP, DZER)
269 IF (IFLV .LE. MXVAL) THEN
270 TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV),iset)
271 IF (IKNL .GT. 0) THEN
278 UPD (JI( IFLV) + IS*(NX+1) + IX + 1,iset) = (TP + TM)/2.
279 UPD (JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = (TP - TM)/2.
283 DO 334 IFLV = NEFF + 1, NFMX
286 UPD(JI( IFLV) + IS*(NX+1) + IX + 1,iset) = 0
287 UPD(JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = 0
292 IF (NFMX .EQ. NEFF) GOTO 21
293 DO 335 IFLV = -NFMX, NFMX+1
294 JI(IFLV) = JI(IFLV) + NITR * (NX+1)
296 CALL CtLhHQRK (NX, TT, NEFF+1, UPD(JI(0)+2,iset),
297 > UPD(JI(NEFF+1)+2,iset))
299 QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1,iset)
300 UPD (JI(NFSN)+IZ+1,iset) = UPD (JI(NFSN)+IZ+1,iset)
302 VS00 = UPD (JI(NFSN)+IZ+1,iset) / (NEFF+1)
303 UPD(JI( NEFF+1) + IZ + 1,iset) = QRKP(NEFF+1) - VS00
305 A = UPD(JI( IFL)+IZ+1,iset)
306 B = UPD(JI(-IFL)+IZ+1,iset)
308 UPD(JI( IFL)+IZ+1,iset) = QRKP(IFL) - VS00
309 IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1,iset) = A - B
313 if(ientry .eq. 1) then
323 if((nx .ne. nxsave) .or.
324 & (nt .ne. ntsave) .or.
325 & (jt .ne. jtsave) .or.
326 & (ng .ne. ngsave) .or.
327 & (xcr .ne. xcrsave) .or.
328 & (xmin .ne. xminsave) .or.
329 & (qini .ne. qinisave) .or.
330 & (qmax .ne. qmaxsave)) then
331 write(6,669) nx, nt, jt, ng, xcr, xmin,
333 669 format(1x,'evolve.f: nx,nt,jt,ng=',4i4,
334 & ' xcr,xmin=',2f9.6,
335 & ' qini, qmax',f7.4,1x,e12.5,' ientry=',i6)
348 FUNCTION CtLhFINTRP (FF, X0, DX, NX, XV, ERR, IR)
349 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
350 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
352 DIMENSION FF (0:NX), XX(MX)
353 DATA SML, XX / 1.D-5, 0., 1.0, 2.0 /
354 DATA IW1, IW3, IW5 / 3 * 0 /
361 CALL CtLhWARNI(IW1, 'Nx < 1, error in CtLhFINTRP.',
362 > 'NX', NX, 1, 256, 1)
369 CALL CtLhWARNR(IW3, 'DX < 0, error in CtLhFINTRP.',
370 > 'DX', DX, D0, D1, 1)
375 IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN
377 > 'X out of range in CtLhFINTRP, Extrapolation used.',
384 ELSEIF (TX .GE. ANX-1.) THEN
390 CALL CtLhRATINT (XX, FF(IX), MNX, DDX, TEM, ERR)
394 FUNCTION CtLhGausInt(F,XL,XR,AERR,RERR,ERR,IRT)
395 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
396 DIMENSION XLIMS(100), R(93), W(93)
397 INTEGER PTR(4),NORD(4)
399 DATA PTR,NORD/4,10,22,46, 6,12,24,48/
400 DATA R/.2386191860,.6612093865,.9324695142,
401 1 .1252334085,.3678314990,.5873179543,.7699026742,.9041172563,
402 1 .9815606342,.0640568929,.1911188675,.3150426797,.4337935076,
403 1 .5454214714,.6480936519,.7401241916,.8200019860,.8864155270,
404 1 .9382745520,.9747285560,.9951872200,.0323801710,.0970046992,
405 1 .1612223561,.2247637903,.2873624873,.3487558863,.4086864820,
406 1 .4669029048,.5231609747,.5772247261,.6288673968,.6778723796,
407 1 .7240341309,.7671590325,.8070662040,.8435882616,.8765720203,
408 1 .9058791367,.9313866907,.9529877032,.9705915925,.9841245837,
409 1 .9935301723,.9987710073,.0162767488,.0488129851,.0812974955,
410 1 .1136958501,.1459737146,.1780968824,.2100313105,.2417431561,
411 1 .2731988126,.3043649444,.3352085229,.3656968614,.3957976498,
412 1 .4254789884,.4547094222,.4834579739,.5116941772,.5393881083,
413 1 .5665104186,.5930323648,.6189258401,.6441634037,.6687183100,
414 1 .6925645366,.7156768123,.7380306437,.7596023411,.7803690438,
415 1 .8003087441,.8194003107,.8376235112,.8549590334,.8713885059,
416 1 .8868945174,.9014606353,.9150714231,.9277124567,.9393703398,
417 1 .9500327178,.9596882914,.9683268285,.9759391746,.9825172636,
418 1 .9880541263,.9925439003,.9959818430,.9983643759,.9996895039/
419 DATA W/.4679139346,.3607615730,.1713244924,
420 1 .2491470458,.2334925365,.2031674267,.1600783285,.1069393260,
421 1 .0471753364,.1279381953,.1258374563,.1216704729,.1155056681,
422 1 .1074442701,.0976186521,.0861901615,.0733464814,.0592985849,
423 1 .0442774388,.0285313886,.0123412298,.0647376968,.0644661644,
424 1 .0639242386,.0631141923,.0620394232,.0607044392,.0591148397,
425 1 .0572772921,.0551995037,.0528901894,.0503590356,.0476166585,
426 1 .0446745609,.0415450829,.0382413511,.0347772226,.0311672278,
427 1 .0274265097,.0235707608,.0196161605,.0155793157,.0114772346,
428 1 .0073275539,.0031533461,.0325506145,.0325161187,.0324471637,
429 1 .0323438226,.0322062048,.0320344562,.0318287589,.0315893308,
430 1 .0313164256,.0310103326,.0306713761,.0302999154,.0298963441,
431 1 .0294610900,.0289946142,.0284974111,.0279700076,.0274129627,
432 1 .0268268667,.0262123407,.0255700360,.0249006332,.0242048418,
433 1 .0234833991,.0227370697,.0219666444,.0211729399,.0203567972,
434 1 .0195190811,.0186606796,.0177825023,.0168854799,.0159705629,
435 1 .0150387210,.0140909418,.0131282296,.0121516047,.0111621020,
436 1 .0101607705,.0091486712,.0081268769,.0070964708,.0060585455,
437 1 .0050142027,.0039645543,.0029107318,.0018539608,.0007967921/
438 DATA TOLABS,TOLREL,NMAX/1.E-35,5.E-4,100/
446 10 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0
447 BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0
450 15 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
454 DO 20 I=PTR(J),PTR(J)-1+NORD(J)
455 20 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
457 TOL=MAX(TOLABS,TOLREL*ABS(VAL))
458 IF (ABS(TVAL-VAL).LT.TOL) THEN
459 CtLhGausInt=CtLhGausInt+VAL
461 IF (NLIMS.NE.0) GO TO 10
469 IF (NLIMS.GT.(NMAX-2)) THEN
470 write(*,50) CtLhGausInt,NMAX,BB-AA,BB+AA
478 50 FORMAT (' CtLhGausInt FAILS, CtLhGausInt,NMAX,XL,XR=',
481 SUBROUTINE CtLhHINTEG (NX, F, H)
482 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
483 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
484 PARAMETER (MXPN = MXF * 2 + 2)
485 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
486 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
487 COMMON / LhCtHINTEC / GH(NDG, MXX)
488 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
489 DIMENSION F(NX), H(NX), G(MXX)
493 TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2)
496 W = XA(I,1) / XA(IY,1)
497 G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W)
499 HTEM = CtLhSMPSNA (NP-2, DZ, G(3), ERR)
501 H(I) = TEM + HTEM + TEM1
503 H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0))
507 SUBROUTINE CtLhHQRK (NX, TT, NQRK, Y, F)
508 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
509 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
510 DIMENSION Y(NX), F(NX)
511 IF (NX .GT. 1) GOTO 11
521 SUBROUTINE CtLhINTEGR (NX, M, F, G, IR)
522 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
524 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
525 PARAMETER (MXPN = MXF * 2 + 2)
526 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
527 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
528 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
529 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
530 DIMENSION F(NX), G(NX)
531 DATA IWRN1, IWRN2 / 0, 0 /
533 IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN
534 MSG = 'NX out of range in CtLhINTEGR call'
535 CALL CtLhWARNI (IWRN1, MSG, 'NX', NX, 0, MXX, 0)
538 IF (M .LT. M1 .OR. M .GT. M2) THEN
539 MSG ='Exponent M out of range in CtLhINTEGR'
540 CALL CtLhWARNI (IWRN2, MSG, 'M', M, M1, M2, 1)
544 TEM = H(1, NX-1, -M) * F(NX-2) + H(2, NX-1, -M) * F(NX-1)
545 > + H(3, NX-1, -M) * F(NX)
549 G(NX-1) = TEM * XA(NX-1, M)
551 DO 10 I = NX-2, 2, -1
552 TEM = TEM + H(1,I,-M)*F(I-1) + H(2,I,-M)*F(I)
553 > + H(3,I,-M)*F(I+1) + H(4,I,-M)*F(I+2)
557 G(I) = TEM * XA(I, M)
560 TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3)
564 G(1) = TEM * XA(1, M)
569 SUBROUTINE CtLhKERNEL
570 >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT)
571 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
572 PARAMETER (PI = 3.141592653589793d0, PI2 = PI ** 2)
573 PARAMETER (D0 = 0.0, D1 = 1.0)
574 DATA CF, CG, TR, IWRN / 1.33333333333333d0, 3.0d0, 0.5d0, 0 /
578 IF (X .LE. 0. .OR. X .GE. 1.) THEN
579 CALL CtLhWARNR(IWRN, 'X out of range in CtLhKERNEL', 'X', X,
591 SPEN2 = CtLhSPENC2 (X)
592 FFP = (1.+ X2) * XM1I
593 FGP = (2.- 2.* X + X2) / X
594 GFP = 1. - 2.* X + 2.* X2
595 GGP = XM1I + XI - 2. + X - X2
596 FFM = (1.+ X2) * XP1I
597 FGM = - (2.+ 2.* X + X2) / X
598 GFM = 1. + 2.* X + 2.* X2
599 GGM = XP1I - XI - 2. - X - X2
600 FF1 = CF * FFP * (1.- X)
603 GG1 = 2.* CG * GGP * X * (1.-X)
604 PCF2 = -2.* FFP *XLN*XLN1M - (3.*XM1I + 2.*X)*XLN
605 > - (1.+X)/2.*XLN2 - 5.*(1.-X)
606 PCFG = FFP * (XLN2 + 11.*XLN/3.+ 67./9.- PI**2 / 3.)
607 > + 2.*(1.+X) * XLN + 40.* (1.-X) / 3.
608 PCFT = (FFP * (- XLN - 5./3.) - 2.*(1.-X)) * 2./ 3.
609 PQQB = 2.* FFM * SPEN2 + 2.*(1.+X)*XLN + 4.*(1.-X)
610 PQQB = (CF**2-CF*CG/2.) * PQQB
611 PQQ2 = CF**2 * PCF2 + CF*CG * PCFG / 2. + CF*TRNF * PCFT
612 PNSP = (PQQ2 + PQQB) * (1.-X)
613 PNSM = (PQQ2 - PQQB) * (1.-X)
614 FFCF2 = - 1. + X + (1.- 3.*X) * XLN / 2. - (1.+ X) * XLN2 / 2.
615 > - FFP * (3.* XLN / 2. + 2.* XLN * XLN1M)
617 FFCFG = 14./3.* (1.-X)
618 > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.)
620 FFCFT = - 16./3. + 40./3.* X + (10.* X + 16./3.* X2 + 2.) * XLN
621 > - 112./9.* X2 + 40./9./X - 2.* (1.+ X) * XLN2
622 > - FFP * (10./9. + 2./3. * XLN)
623 FGCF2 = - 5./2.- 7./2.* X + (2.+ 7./2.* X) * XLN + (X/2.-1.)*XLN2
625 > - FGP * (3.* XLN1M + XLN1M ** 2)
626 FGCFG = 28./9. + 65./18.* X + 44./9. * X2 - (12.+ 5.*X + 8./3.*X2)
627 > * XLN + (4.+ X) * XLN2 + 2.* X * XLN1M
628 > + FGP * (-2.*XLN*XLN1M + XLN2/2. + 11./3.*XLN1M + XLN1M**2
631 FGCFT = -4./3.* X - FGP * (20./9.+ 4./3.*XLN1M)
632 GFCFT = 4.- 9.*X + (-1.+ 4.*X)*XLN + (-1.+ 2.*X)*XLN2 + 4.*XLN1M
633 > + GFP * (-4.*XLN*XLN1M + 4.*XLN + 2.*XLN2 - 4.*XLN1M
634 > + 2.*XLN1M**2 - 2./3.* PI2 + 10.)
635 GFCGT = 182./9.+ 14./9.*X + 40./9./X + (136./3.*X - 38./3.)*XLN
636 > - 4.*XLN1M - (2.+ 8.*X)*XLN2
637 > + GFP * (-XLN2 + 44./3.*XLN - 2.*XLN1M**2 + 4.*XLN1M
638 > + PI2/3. - 218./9.)
640 GGCFT = -16.+ 8.*X + 20./3.*X2 + 4./3./X + (-6.-10.*X)*XLN
641 > - 2.* (1.+ X) * XLN2
642 GGCGT = 2.- 2.*X + 26./9.*X2 - 26./9./X - 4./3.*(1.+X)*XLN
644 GGCG2 = 27./2.*(1.-X) + 67./9.*(X2-XI) + 4.*(1.+X)*XLN2
645 > + (-25.+ 11.*X - 44.*X2)/3.*XLN
646 > + GGP * (67./9.- 4.*XLN*XLN1M + XLN2 - PI2/3.)
648 FF2 = CF * TRNF * FFCFT + CF ** 2 * FFCF2 + CF * CG * FFCFG
649 FG2 = CF * TRNF * FGCFT + CF ** 2 * FGCF2 + CF * CG * FGCFG
650 GF2 = CF * TRNF * GFCFT + CG * TRNF * GFCGT
651 GG2 = CF * TRNF * GGCFT + CG ** 2 * GGCG2 + CG * TRNF * GGCGT
652 XLG = (LOG(1./(1.-X)) + 1.)
654 FF2 = FF2 * X * (1.- X)
657 GG2 = GG2 * X * (1.- X)
660 SUBROUTINE CtLhLAMCWZ
661 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
662 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
664 CALL CtLhSETL1 (NF, AL)
666 FUNCTION LhCtNAMQCD(NNAME)
667 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
668 CHARACTER NNAME*(*), NAME*8
669 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
676 IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR.
677 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') )
679 IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR.
680 1 (NAME .EQ. '# FL') )
683 IF (NAME .EQ. 'M'//CHAR(I+IASC0))
687 IF (NAME .EQ. 'LAM'//CHAR(I+IASC0))
690 IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') LhCtNAMQCD = 24
693 FUNCTION LhCtNFL(AMU)
694 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
695 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
696 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
698 IF (.NOT. SET) CALL CtLhLAMCWZ
700 IF ((LhCtNFL .EQ. NF) .OR. (AMU .LE. AMN)) GOTO 20
701 DO 10 I = NF - NHQ + 1, NF
702 IF (AMU .GE. AMHAT(I)) THEN
710 SUBROUTINE CtLhNSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN)
711 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
712 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
713 PARAMETER (MXPN = MXF * 2 + 2)
714 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
715 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
716 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
717 DIMENSION U0(NX), UN(0:NX, 0:NT)
718 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
721 IF (NX .GT. MXX) THEN
722 WRITE (*,*) 'Nx =', NX, ' greater than Max pts in CtLhNSEVL.'
723 STOP 'Program stopped in CtLhNSEVL'
725 c ++ remove unused quantities (jcp)
726 c ++ TMD = TIN + DT * NT / 2.
728 c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
729 c ++ TLAM = TMD - TEM
733 UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1)
740 IRND = (IS-1) * JT + JS
741 IF (IRND .EQ. 1) THEN
742 CALL RHS (TT, Neff, Y0, F0)
744 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
747 CALL RHS (TT, NEFF, Y0, F1)
749 Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
752 CALL RHS (TT, NEFF, Y1, F1)
754 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
757 CALL RHS (TT, NEFF, YP, FP)
759 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
767 UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
771 SUBROUTINE LhCtNSRHSM (TT, NEFF, FI, FO)
772 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
774 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
775 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
776 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
777 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
778 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
779 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
780 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
781 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
782 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
783 DIMENSION G1(MXX), FI(NX), FO(NX)
784 DIMENSION W0(MXX), W1(MXX), WH(MXX)
788 CPL2= CPL ** 2 / 2. * S
790 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
791 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
792 CALL CtLhHINTEG (NX, FI, WH)
794 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
795 FO(IZ) = CPL * FO(IZ)
797 IF (IKNL .EQ. 2) THEN
806 G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX))
808 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
809 TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2
810 FO(IX) = FO(IX) + TMP2
815 SUBROUTINE LhCtNSRHSP (TT, NEFF, FI, FO)
816 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
818 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
819 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
820 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
821 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
822 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
823 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
824 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
825 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
826 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
827 DIMENSION G1(MXX), FI(NX), FO(NX)
828 DIMENSION W0(MXX), W1(MXX), WH(MXX)
832 CPL2= CPL ** 2 / 2. * S
834 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
835 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
836 CALL CtLhHINTEG (NX, FI, WH)
838 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
839 FO(IZ) = CPL * FO(IZ)
841 IF (IKNL .EQ. 2) THEN
847 XY = ZZ (NX-IX+1, NX-IY+1)
848 G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX))
850 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
851 TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + ZQQB)) * CPL2
852 FO(IX) = FO(IX) + TMP2
857 FUNCTION CtLhPARDIS (IPRTN, XX, QQ)
858 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
859 include 'parmsetup.inc'
862 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
863 PARAMETER (MXPN = MXF * 2 + 2)
864 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
865 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
866 PARAMETER (Smll = 1D-9)
868 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
869 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
870 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
871 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
872 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
873 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
874 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
875 COMMON / LhCtCOMQMS / VALQMS(9)
876 dimension fvec(4), fij(4)
877 dimension xvpow(0:mxx)
878 Data Iwrn1, Iwrn2, Iwrn3, OneP / 3*0, 1.00001 /
879 data xpow / 0.3d0 / !**** choice of interpolation variable
883 save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
884 save const1 , const2, const3, const4, const5, const6
885 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
889 if(nx .ne. nxsave) then
892 xvpow(i) = xv(i)**xpow
900 if((x.lt.xmin).or.(x.gt.1.d0)) print 98,x
901 98 format(' WARNING: X=',e12.5,' OUT OF RANGE')
902 if((q.lt.qini).or.(q.gt.qmax)) print 99,q
903 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE')
905 c enforce threshold early to improve speed...
908 if(QQ .lt. VALQMS(ii) ) then
914 c force pardis = 0.0d0 at exactly =1.0d0 - added mrw 10/May/06
915 if(xx .eq. 1.0d0) then
920 c skip the initialization in x if same as in the previous call.
921 if(x .eq. xlast) goto 100
926 11 If (JU-JLx .GT. 1) Then
928 If (X .Ge. XV(JM)) Then
935 If (JLx .LE. -1) Then
936 Print '(A,1pE12.4)','Severe error: x <= 0 in ParDis x=', x
938 ElseIf (JLx .Eq. 0) Then
940 Msg = '0 < X < Xmin in ParDis; extrapolation used!'
941 CALL CtLhWARNR (IWRN1, Msg, 'X', X, Xmin, 1D0, 1)
942 Elseif (JLx .LE. Nx-2) Then
944 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
947 Print '(A,1pE12.4)','Severe error: x > 1 in ParDis x=', x
951 If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
969 sdet = s12*s34 - s1213*s2434
971 const5 = (s34*sy2-s2434*sy3)*tmp/s12
972 const6 = (s1213*sy2-s12*sy3)*tmp/s34
977 c skip the initialization in q if same as in the previous call.
978 if(q .eq. qlast) goto 110
985 12 If (JU-JLq .GT. 1) Then
987 If (Q .GE. QV(JM)) Then
997 Msg = 'Q < Q0 in ParDis; extrapolation used!'
998 CALL CtLhWARNR (IWRN2, Msg, 'Q', Q, Qini, 1D0, 1)
1000 Elseif (JLq .LE. Nt-2) Then
1004 If (JLq .GE. Nt) Then
1005 Msg = 'Q > Qmax in ParDis; extrapolation used!'
1006 CALL CtLhWARNR (IWRN3, Msg, 'Q', Q, Qmax, 1D0, 1)
1010 If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
1024 tdet = t12*t34 - tmp1*tmp2
1029 jtmp = ((IPRTN + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
1031 J1 = jtmp + it*(NX+1)
1034 fij(2) = Upd(J1+1,iset) * Xa(1,2)
1035 fij(3) = Upd(J1+2,iset) * Xa(2,2)
1036 fij(4) = Upd(J1+3,iset) * Xa(3,2)
1037 Call CtLhPolint4 (XVpow(0), Fij(1), 4, ss, Fx, Dfx)
1039 If (x .GT. 0D0) Fvec(it) = Fx / x**2
1040 ElseIf (JLx .Eq. Nx-1) Then
1041 Call CtLhPolint4 (XVpow(Nx-3), Upd(J1,iset), 4, ss, Fx, Dfx)
1044 sf2 = Upd(J1+1,iset)
1045 sf3 = Upd(J1+2,iset)
1046 Fvec(it) = (const5*(Upd(J1,iset)
1047 & - sf2*const1 + sf3*const2)
1048 & + const6*(Upd(J1+3,iset)
1049 & + sf2*const3 - sf3*const4)
1050 & + sf2*sy3 - sf3*sy2) / s23
1053 If (JLq .LE. 0) Then
1054 Call CtLhPolint4 (TV(0), Fvec(1), 4, tt, ff, Dfq)
1055 ElseIf (JLq .GE. Nt-1) Then
1056 Call CtLhPolint4 (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
1060 g1 = ( tf2*t13 - tf3*t12) / t23
1061 g4 = (-tf2*t34 + tf3*t24) / t23
1062 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
1063 & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
1064 ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
1070 SUBROUTINE CtLhPARPDF (IACT, NAME, VALUE, IRET)
1071 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1072 CHARACTER NAME*(*), Uname*10
1074 DATA ILEVEL, LRET / 1, 1 /
1076 CALL CtLhUPC (NAME, Ln, Uname)
1077 IF (IACT .EQ. 0 .OR. IACT .EQ. 4) then
1078 c > IVALUE = NINT (VALUE) !tentatively remove this since it seems not to be used
1079 print *,'Fatal error: iact=',iact
1082 START1 = (IACT .NE. 1) .AND. (IACT .NE. 2)
1083 c prepare to remove this stuff, since I think IACT=1 or 2 always
1085 print *,'Fatal error: start1=',start1
1088 IF (START1) ILEVEL = 1
1092 CALL CtLhParQcd (IACT, Uname(1:Ln), VALUE, JRET)
1093 IF (JRET .EQ. 1) GOTO 11
1094 IF (JRET .EQ. 2) GOTO 12
1095 IF (JRET .EQ. 3) GOTO 13
1096 IF (JRET .GT. 4) GOTO 15
1098 2 CALL CtLhEVLPAR (IACT, Uname(1:Ln), VALUE, JRET)
1099 IF (JRET .EQ. 1) GOTO 11
1100 IF (JRET .EQ. 2) GOTO 12
1101 IF (JRET .EQ. 3) GOTO 13
1102 IF (JRET .GT. 4) GOTO 15
1104 IF (.NOT. START1) GOTO 1
1105 IF (JRET .EQ. 0) GOTO 10
1113 IF (JRET .NE. 4) LRET = JRET
1114 IF (LRET.EQ.0 .OR. LRET.EQ.2 .OR. LRET.EQ.3) THEN
1115 PRINT *, 'Error in CtLhPARPDF: IRET, IACT, NAME, VALUE =',
1116 > LRET, IACT, NAME, VALUE
1117 PRINT *, 'fatal error in CtLhparpdf'
1123 SUBROUTINE CtLhParQcd(IACT,NAME,VALUE,IRET)
1124 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1129 WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ',
1130 > '(i in 1 to 9), LAMi (i in 1 to NFL)'
1132 ELSEIF (IACT.EQ.1) THEN
1133 CALL CtLhQCDSET (NAME,VALUE,IRET)
1134 ELSEIF (IACT.EQ.2) THEN
1135 CALL CtLhQCDGET (NAME,VALUE,IRET)
1141 FUNCTION CtLhPFF1 (XX)
1142 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1143 LOGICAL LA, LB, LSTX
1144 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1145 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1146 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1148 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1149 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1150 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1151 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1152 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1153 > PNSP(0:MXX), PNSM(0:MXX)
1155 DATA LA, LB / 2 * .FALSE. /
1159 IF (LA .AND. .NOT.LB) THEN
1168 ElseIF (X .GE. XMIN) THEN
1170 TEM = CtLhFINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT)
1172 CALL CtLhPOLIN1 (XL, FF1(1), MX, X, TEM, ERR)
1176 CtLhPFF1 = TEM / (1.-X)
1179 CtLhTFF1 = TEM / (1.-X) * CtLhDXDZ(Z)
1184 c +++ something is wrong, since QFF1 and RFF1 are not used.
1185 c +++ but this code appears to only be used for extrapolation
1186 c +++ to small x, which is unreliable anyway, so ignore for now (jcp)
1190 RFF1 = TEM * X / (1.-X)
1199 ElseIF (X .GE. XMIN) THEN
1201 TEM = CtLhFINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT)
1203 CALL CtLhPOLIN1 (XL, PNSP(1), MX, X, TEM, ERR)
1205 CtLhFNSP = TEM / (1.- X)
1212 ElseIF (X .GE. XMIN) THEN
1214 TEM = CtLhFINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT)
1216 CALL CtLhPOLIN1 (XL, PNSM(1), MX, X, TEM, ERR)
1218 CtLhFNSM = TEM / (1.- X)
1223 CtLhRGG1= 0 !error corrected? (jcp)
1225 ElseIF (X .GE. XMIN) THEN
1227 TEM = CtLhFINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT)
1229 CALL CtLhPOLIN1 (XL, GG1(1), MX, X, TEM, ERR)
1232 CtLhRGG1 = TEM / X / (1.-X) !error corrected? (jcp)
1236 CtLhRGG1 = TEM / X !error corrected? (jcp)
1239 CtLhRGG1 = TEM / (1.-X)
1246 CtLhRFF2 = 0 !error corrected? (jcp)
1248 ElseIF (X .GE. XMIN) THEN
1250 TEM = CtLhFINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT)
1252 CALL CtLhPOLIN1 (XL, FF2(1), MX, X, TEM, ERR)
1255 CtLhRFF2 = TEM / X / (1.-X) !error corrected? (jcp)
1259 CtLhRFF2 = TEM / X !error corrected? (jcp)
1262 CtLhRFF2 = TEM / (1.-X)
1269 CtLhRGG2 = 0 !error corrected? (jcp)
1271 ElseIF (X .GE. XMIN) THEN
1273 TEM = CtLhFINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT)
1275 CALL CtLhPOLIN1 (XL, GG2(1), MX, X, TEM, ERR)
1278 CtLhRGG2 = TEM / X / (1.-X) !error corrected? (jcp)
1282 CtLhRGG2 = TEM / X !error corrected? (jcp)
1285 CtLhRGG2 = TEM / (1.-X)
1290 SUBROUTINE CtLhPOLIN1 (XA,YA,N,X,Y,DY)
1291 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1293 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1298 IF (DIFT.LT.DIF) THEN
1317 IF (2*NS.LT.N-M)THEN
1327 SUBROUTINE CtLhQARRAY (NINI)
1328 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1329 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1330 PARAMETER (MXPN = MXF * 2 + 2)
1331 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
1332 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1333 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
1334 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1336 IF (NT .GE. mxq) NT = mxq - 1
1341 1 DT0 = (TMAX - TINI) / float(NT)
1342 NINI = LhCtNFL(QINI)
1343 NFMX = LhCtNFL(QMAX)
1344 Call CtLhParQcd (2, 'ORDER', Ord, Ir)
1345 Call CtLhParQcd (2, 'ALAM', Al0, Ir)
1346 Call CtLhParQcd (2, 'NFL', Afl0, Ir)
1348 Call CtLhParQcd (1, 'NFL', AFL, Ir)
1351 Call CtLhSetLam (Ifl0, Al0, Iordr)
1352 NG = NFMX - NINI + 1
1361 DO 20 NEFF = NINI, NFMX
1362 ICNT = NEFF - NINI + 1
1363 IF (NEFF .LT. NFMX) THEN
1364 THRN = CtLhAMHATF (NEFF + 1)
1365 QOUN = MIN (QMAX, THRN)
1369 IF (QOUN-QOUT .LE. 0.0001) THEN
1377 NITR = INT (TEM / DT0) + 1
1383 NTL (ICNT+1) = NTL(ICNT) + NITR
1384 IF (NITR .NE. 0) THEN
1386 TV (NTL(ICNT)+I) = TIN + DT * I
1387 S = EXP (TV(NTL(ICNT)+I))
1388 QV (NTL(ICNT)+I) = AL * EXP (S)
1397 IF (NTP .GE. MXQ) THEN
1398 NT = MXQ - ND - NCNT
1404 SUBROUTINE CtLhQCDGET(NAME,VALUE,IRET)
1405 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1407 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1408 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1409 COMMON / LhCtCOMQMS / VALQMS(9)
1411 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1412 ICODE = LhCtNAMQCD(NAME)
1414 IF (ICODE .EQ. 1) THEN
1416 ELSEIF (ICODE .EQ. 2) THEN
1418 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 12)) THEN
1419 VALUE = VALQMS(ICODE - 2)
1420 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1421 VALUE = ALAM(ICODE - 13)
1422 ELSEIF (ICODE .EQ. 24) THEN
1428 SUBROUTINE CtLhQCDSET (NAME,VALUE,IRET)
1429 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1431 COMMON / LhCtCOMQMS / VALQMS(9)
1432 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1434 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1435 IVALUE = NINT(VALUE)
1436 ICODE = LhCtNAMQCD(NAME)
1437 IF (ICODE .EQ. 0) THEN
1439 c print *,'warning empty CtLhQCDSET call: NAME=',
1440 c & NAME,' VALUE=',VALUE
1444 IF (ICODE .EQ. 1) THEN
1445 IF (VALUE.LE.0) GOTO 12
1447 ELSEIF (ICODE .EQ. 2) THEN
1448 IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12
1450 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 11)) THEN
1451 IF (VALUE .LT. 0) GOTO 12
1452 Scle = Min (Value , VALQMS(ICODE - 2))
1453 AlfScle = CtLhALPI(Scle) * Pi
1454 VALQMS(ICODE - 2) = VALUE
1455 Call CtLhAlfSet (Scle, AlfScle)
1456 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1457 IF (VALUE .LE. 0) GOTO 12
1458 CALL CtLhSETL1 (ICODE-13, VALUE)
1459 ELSEIF (ICODE .EQ. 24) THEN
1460 IF ((IVALUE .LT. 1) .OR. (IVALUE .GT. 2)) GOTO 12
1463 IF (.NOT. SET) CALL CtLhLAMCWZ
1469 FUNCTION CtLhQZBRNT(FUNC, X1, X2, TOLIN, IRT)
1470 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1471 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
1478 IF(FB*FA.GT.0.) THEN
1479 WRITE (*, *) 'Root must be bracketed for CtLhQZBRNT.'
1484 IF(FB*FC.GT.0.) THEN
1490 IF(ABS(FC).LT.ABS(FB)) THEN
1498 TOL1=2.*EPS*ABS(B)+0.5*TOL
1500 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
1504 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
1512 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
1513 Q=(Q-1.)*(R-1.)*(S-1.)
1517 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
1530 IF(ABS(D) .GT. TOL1) THEN
1537 WRITE (*, *) 'CtLhQZBRNT exceeding maximum iterations.'
1542 SUBROUTINE CtLhRATINT(XA,YA,N,X,Y,DY)
1543 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1544 PARAMETER (NMAX=10,TINY=1.E-25)
1545 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1554 ELSE IF (H.LT.HH) THEN
1573 IF (2*NS.LT.N-M)THEN
1583 FUNCTION CtLhRTALF (EFLLN)
1584 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1585 include 'parmsetup.inc'
1586 PARAMETER (PI = 3.141592653589793d0)
1587 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
1588 EFMULM = EXP (EFLLN)
1590 TEM2 = 1. / CtLhALPQCD (JORD, NEFF, EFMULM, I)
1591 CtLhRTALF = TEM1 - TEM2
1593 Subroutine CtLhbldat1
1594 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1595 include 'parmsetup.inc'
1597 PARAMETER (MXX = 105, MXQ = 25, MxF = 6)
1598 PARAMETER (MxPN = MxF * 2 + 2)
1599 PARAMETER (MxQX= MXQ * MXX, MxPQX = MxQX * MxPN)
1600 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1601 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1602 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1603 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
1604 PARAMETER (NF0 = 4, Nshp = 8,NEX = Nshp+2)
1610 SUBROUTINE CtLhSETL1 (NEF, VLAM)
1611 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1613 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1614 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1615 COMMON / LhCtCOMQMS / VALQMS(9)
1616 IF (NEF .LT. 0 .OR. NEF .GT. NF) THEN
1617 WRITE(*,*)'NEF out of range in CtLhSETL1: NEF NF =',NEF,NF
1622 AMHAT(N) = VALQMS(N)
1625 DO 10 N = NEF, 1, -1
1626 CALL CtLhTRNLAM(NORDER, N, -1, IR1)
1629 CALL CtLhTRNLAM(NORDER, N, 1, IR1)
1631 DO 30, N = NF, 1, -1
1632 IF ((ALAM(N) .GE. 0.7 * AMHAT(N))
1633 > .OR. (ALAM(N-1) .GE. 0.7 * AMHAT(N)))THEN
1640 DO 50, N = NF-NHQ, 1, -1
1646 IF (ALAM(N) .GT. AMN) AMN = ALAM(N)
1653 SUBROUTINE CtLhSETLAM (NEF, WLAM, IRDR)
1654 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1655 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1657 IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN
1658 WRITE(*,*)'NEF out of range in CtLhSETLAM: NEF NF=',NEF,NF
1662 IF (IRDR .NE. NORDER) then
1663 PRINT *,'fatal error: wanted cnvl1'
1666 CALL CtLhSETL1 (NEF, VLAM)
1668 Subroutine CtLhbldat2
1669 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1670 COMMON / LhCtCOMQMS / VALQMS(9)
1671 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1688 FUNCTION CtLhSMPNOL (NX, DX, FN, ERR)
1689 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1692 IF (NX .LE. 1 .OR. NX .GT. 1000) THEN
1693 PRINT *, 'NX =', NX, ' OUT OF RANGE IN CtLhSMPNOL!'
1695 ELSEIF (NX .EQ. 2) THEN
1697 ELSEIF (NX .EQ. 3) THEN
1698 TEM = DX * FN(2) * 2.
1701 TEM = DX * (23.* FN(2) - 16.* FN(3) + 5.* FN(4)) / 12.
1702 TMP = DX * (3.* FN(2) - FN(3)) / 2.
1703 ERR = ABS(TEM - TMP)
1704 TEM = TEM + CtLhSMPSNA (NX-1, DX, FN(2), ER1)
1705 ERR = ABS(ER1) + ERR
1707 TEM = DX * (8.* FN(2) - 4.* FN(3) + 8.* FN(4)) / 3.
1708 TMP = DX * (3.* FN(2) + 2.* FN(3) + 3.* FN(4)) / 2.
1709 ERR = ABS(TEM - TMP)
1710 TEM = TEM + CtLhSMPSNA (NX-4, DX, FN(5), ER1)
1711 ERR = ABS(ER1) + ERR
1717 FUNCTION CtLhSMPSNA (NX, DX, F, ERR)
1718 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1719 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1720 PARAMETER (MAXX = 1000)
1722 DATA IW1, IW2, TINY / 2*0, 1.E-35 /
1723 IF (DX .LE. 0.) THEN
1724 CALL CtLhWARNR(IW2,'DX cannot be < 0. in CtLhSMPSNA', 'DX',
1729 IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN
1730 CALL CtLhWARNI(IW1, 'NX out of range in CtLhSMPSNA', 'NX', NX,
1733 ELSEIF (NX .EQ. 1) THEN
1735 ELSEIF (NX .EQ. 2) THEN
1736 SIMP = (F(1) + F(2)) / 2.
1737 ERRD = (F(1) - F(2)) / 2.
1741 ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24.
1748 SIMP = (F(1) + 4.* F(2) + F(3)) / 3.
1749 TRPZ = (F(1) + 2.* F(2) + F(3)) / 2.
1759 SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3.
1760 TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2.
1765 CtLhSMPSNA = SIMP * DX
1766 IF (ABS(SIMP) .GT. TINY) THEN
1773 SUBROUTINE CtLhSNEVL(IKNL,NX,NT,JT,DT,TIN,NEFF,UI,GI,US,GS)
1774 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1775 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1776 PARAMETER (MXQX= MXQ * MXX)
1777 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1778 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1779 DIMENSION UI(NX), US(0:NX, 0:NT)
1780 DIMENSION GI(NX), GS(0:NX, 0:NT)
1781 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
1782 DIMENSION Z0(MXX), Z1(MXX), ZP(MXX), G0(MXX), G1(MXX), GP(MXX)
1786 IF (NX .GT. MXX) THEN
1787 WRITE (*,*) 'Nx =', NX, ' too many pts in CtLhSNEVL'
1788 STOP 'Program stopped in CtLhSNEVL'
1790 c ++ remove unused quantities (jcp)
1791 c ++ TMD = TIN + DT * NT / 2.
1793 c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
1794 c ++ TLAM = TMD - TEM
1799 US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3)
1800 GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3)
1808 IRND = (IS-1) * JTT + JS
1809 IF (IRND .EQ. 1) THEN
1810 CALL CtLhSNRHS (TT, NEFF, Y0,Z0, F0,G0)
1812 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
1813 Z0(IZ) = Z0(IZ) + DDT * G0(IZ)
1816 CALL CtLhSNRHS (TT, NEFF, Y0, Z0, F1, G1)
1818 Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
1819 Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0
1822 CALL CtLhSNRHS (TT, NEFF, Y1, Z1, F1, G1)
1824 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
1825 ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0
1828 CALL CtLhSNRHS (TT, NEFF, YP, ZP, FP, GP)
1830 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
1831 Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0
1838 IF (IKNL .GT. 0) THEN
1839 US (IX, IS) = MAX(Y1(IX), D0)
1840 GS (IX, IS) = MAX(Z1(IX), D0)
1842 US (IX, IS) = Y1(IX)
1843 GS (IX, IS) = Z1(IX)
1846 US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
1847 GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3)
1851 SUBROUTINE CtLhSNRHS (TT, NEFF, FI, GI, FO, GO)
1852 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1854 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1855 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1856 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1857 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1858 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1859 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1860 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1861 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1862 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1863 DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX)
1864 DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), WM(MXX)
1865 DIMENSION R0(MXX), R1(MXX), R2(MXX), RH(MXX), RM(MXX)
1869 CPL2= CPL ** 2 / 2. * S
1871 CALL CtLhINTEGR (NX,-1, FI, WM, IR1)
1872 CALL CtLhINTEGR (NX, 0, FI, W0, IR2)
1873 CALL CtLhINTEGR (NX, 1, FI, W1, IR3)
1874 CALL CtLhINTEGR (NX,-1, GI, RM, IR4)
1875 CALL CtLhINTEGR (NX, 0, GI, R0, IR5)
1876 CALL CtLhINTEGR (NX, 1, GI, R1, IR6)
1877 CALL CtLhINTEGR (NX, 2, GI, R2, IR7)
1878 CALL CtLhHINTEG (NX, FI, WH)
1879 CALL CtLhHINTEG (NX, GI, RH)
1880 IF (IKNL .GT. 0) THEN
1882 FO(IZ) = ( 2D0 * FI(IZ)
1883 > + 4D0 / 3D0 * ( 2D0 * WH(IZ) - W0(IZ) - W1(IZ) ))
1884 > + NEFF * ( R0(IZ) - 2D0 * R1(IZ) + 2D0 * R2(IZ) )
1885 FO(IZ) = FO(IZ) * CPL
1886 GO(IZ) = 4D0 / 3D0 * ( 2D0 * WM(IZ) - 2D0 * W0(IZ) + W1(IZ) )
1887 > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ)
1888 > + 6D0 * (RH(IZ) + RM(IZ) - 2D0 * R0(IZ) + R1(IZ) - R2(IZ))
1889 GO(IZ) = GO(IZ) * CPL
1893 FO(IZ) = NEFF * (-R0(IZ) + 2.* R1(IZ) )
1894 > + 2.* FI(IZ) + 4./ 3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ) )
1895 FO(IZ) = FO(IZ) * CPL
1896 GO(IZ) = 4./ 3.* ( 2.* W0(IZ) - W1(IZ) )
1897 >+ (33.- 2.* NEFF) / 6.* GI(IZ) + 6.*(RH(IZ) + R0(IZ) - 2.* R1(IZ))
1898 GO(IZ) = GO(IZ) * CPL
1901 IF (IKNL .EQ. 2) THEN
1912 G1(KZ) = FFG(I, IY) * (FI(IY) - XY**2 *FI(I))
1913 G4(KZ) = GGF(I, IY) * (GI(IY) - XY**2 *GI(I))
1914 G2(KZ) = FFG(IS,IT) * (GI(IY) - xy*GI(I)) !FG
1915 G3(KZ) = GGF(IS,IT) * (FI(IY) - XY*FI(I)) !GF (usual notations)
1917 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
1918 TEM2 = CtLhSMPSNA (NP, DZ, G2, ERR)
1919 TEM3 = CtLhSMPSNA (NP, DZ, G3, ERR)
1920 TEM4 = CtLhSMPNOL (NP, DZ, G4, ERR)
1921 TEM1 = TEM1 - FI(I) * (AFF2(I) + ZGF2)
1922 TEM4 = TEM4 - GI(I) * (AGG2(I) + ZFG2)
1923 tem2 = tem2 + GI(I)*AFG2(I)
1924 tem3= tem3 + FI(I)*AGF2(I)
1927 FO(I) = FO(I) + TMF * CPL2
1928 GO(I) = GO(I) + TMG * CPL2
1933 FUNCTION CtLhSPENC2 (X)
1934 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1936 COMMON / LhCtSPENCC / XX
1937 DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 /
1939 TEM = CtLhGausInt(CtLhSPN2IN, XX, U1, AERR, RERR, ERR, IRT)
1940 CtLhSPENC2 = TEM + LOG (XX) ** 2 / 2.
1943 FUNCTION CtLhSPN2IN (ZZ)
1944 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1945 COMMON / LhCtSPENCC / X
1947 TEM = LOG (1.+ X - Z) / Z
1951 SUBROUTINE CtLhSTUPKL (NFL)
1952 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1954 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1955 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1957 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1958 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1959 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1960 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1961 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1962 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1963 > PNSP(0:MXX), PNSM(0:MXX)
1964 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1965 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1966 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1967 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1968 EXTERNAL CtLhPFF1, CtLhRGG1, CtLhRFF2, CtLhRGG2
1969 EXTERNAL CtLhFNSP, CtLhFNSM
1970 dimension aff1(mxx),agg1(mxx)
1971 PARAMETER (PI = 3.141592653589793d0, PI2 = PI**2)
1972 DATA CF, CG, TR / 1.33333333333333d0, 3.0, 0.5 /
1973 data zeta3/1.20205690315959d0/ ! zeta(3.0)
1975 DATA AERR, RERR / 0.0, 0.02 /
1983 CALL CtLhKERNEL (XZ, FF1(I), GF1(I), FG1(I), GG1(I), PNSP(I),
1984 > PNSM(I), FF2(I), GF2(I), FG2(I), GG2(I), NFL, IRT)
1986 FF1(0) = FF1(1) * 3. - FF1(2) * 3. + FF1(3)
1987 FG1(0) = FG1(1) * 3. - FG1(2) * 3. + FG1(3)
1988 GF1(0) = GF1(1) * 3. - GF1(2) * 3. + GF1(3)
1989 GG1(0) = GG1(1) * 3. - GG1(2) * 3. + GG1(3)
1990 PNSP(0) = PNSP(1) * 3. - PNSP(2) * 3. + PNSP(3)
1991 PNSM(0) = PNSM(1) * 3. - PNSM(2) * 3. + PNSM(3)
1992 FF2(0) = FF2(1) * 3. - FF2(2) * 3. + FF2(3)
1993 FG2(0) = FG2(1) * 3. - FG2(2) * 3. + FG2(3)
1994 GF2(0) = GF2(1) * 3. - GF2(2) * 3. + GF2(3)
1995 GG2(0) = GG2(1) * 3. - GG2(2) * 3. + GG2(3)
1996 FF1(NX) = FF1(NX-1) * 3. - FF1(NX-2) * 3. + FF1(NX-3)
1997 FG1(NX) = FG1(NX-1) * 3. - FG1(NX-2) * 3. + FG1(NX-3)
1998 GF1(NX) = GF1(NX-1) * 3. - GF1(NX-2) * 3. + GF1(NX-3)
1999 GG1(NX) = GG1(NX-1) * 3. - GG1(NX-2) * 3. + GG1(NX-3)
2000 PNSM(NX) = PNSM(NX-1) * 3. - PNSM(NX-2) * 3. + PNSM(NX-3)
2001 PNSP(NX) = PNSP(NX-1) * 3. - PNSP(NX-2) * 3. + PNSP(NX-3)
2002 FF2(NX) = FF2(NX-1) * 3. - FF2(NX-2) * 3. + FF2(NX-3)
2003 FG2(NX) = FG2(NX-1) * 3. - FG2(NX-2) * 3. + FG2(NX-3)
2004 GF2(NX) = GF2(NX-1) * 3. - GF2(NX-2) * 3. + GF2(NX-3)
2005 GG2(NX) = GG2(NX-1) * 3. - GG2(NX-2) * 3. + GG2(NX-3)
2007 AFF1(1) = CtLhGausInt(CtLhPFF1,D0,XV(1),AERR,RERR,ER1,IRT)
2009 TMPG = CtLhGausInt(CtLhRGG1,D0,XV(1),AERR,RERR,ER3,IRT)
2010 AGG1(1) = TMPG + DGG1
2011 ANSM(1) = CtLhGausInt(CtLhFNSM,D0,XV(1),AERR,RER,ER2,IRT)
2012 ANSP(1) = CtLhGausInt(CtLhFNSP,D0,XV(1),AERR,RER,ER2,IRT)
2014 AFF2(1) = CtLhGausInt(CtLhRFF2, D0, XV(1), AER, RER, ER2, IRT)
2016 AGG2(1) = CtLhGausInt(CtLhRGG2, D0, XV(1), AER, RER, ER4, IRT)
2018 TEM =CtLhGausInt(CtLhPFF1,XV(I2-1),XV(I2),AERR,RERR,ER1,IRT)
2019 AFF1(I2) = TEM + AFF1(I2-1)
2020 AER = ABS(TEM * RER)
2021 AFF2(I2)=CtLhGausInt(CtLhRFF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT)
2023 TEM = CtLhGausInt(CtLhRGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT)
2025 AGG1(I2) = TMPG + DGG1
2026 AER = ABS(TEM * RER)
2027 AGG2(I2)=CtLhGausInt(CtLhRGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)
2029 ANSP(I2)=CtLhGausInt(CtLhFNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2031 ANSM(I2)=CtLhGausInt(CtLhFNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2034 ANSP(NX)=CtLhGausInt(CtLhFNSP,XV(NX-1),D1,AERR,RER,ERR,
2036 ANSM(NX)=CtLhGausInt(CtLhFNSM,XV(NX-1),D1,AERR,RER,ERR,
2039 do i2=1,nx-1 !loop over x
2041 c XI = 1./ X !unused - jcp
2046 XLN1M = DLOG (1.- X)
2050 xLi31m=CtLhxLi(3,1d0-x)
2059 > (9 + 4*Pi2 - 22*x + 13*x2 + 6*(3 - 4*x + x2)*xln1m +
2060 > 40*xln - 24*xLi2)/9.
2062 > (6*(2*(-9 + Pi2) + 3*x*(5 + x)) +4*(3 +2*Pi2+3*x*(-3 + 2*x))*
2063 > xln1m + 6*x3m*x1m*xln1m2 - 6*(x*(8 + 3*x) + 4*xln1m2)*
2064 > xln - 3*(-4 + x)*x*xln2)/12 - 2*(3 + 2*xln1m)*xLi2 - 4*xLi31m
2066 > (3637-186*Pi2-x*(3198+72*Pi2+x*(231 + 208*x)))/108.- xln +
2067 > (3*xln1m*(-33 - 4*Pi2 + (50 - 17*x)*x - 3*x3m*x1m*xln1m) +
2068 > 2*(x*(198 + x*(27+8*x))+9*xln1m*(3 - 4*x + x2 + 2*xln1m))*
2069 > xln - 9*x*(4 + x)*xln2)/18- x1p*x3p*xln*xln1p-
2070 > (x1p*x3p - 4*xln)*xLi2m + (31d0/3d0 +4*xln1m- 4*xln)*xLi2 +
2071 > 4*xLi31m + 12*xLi3 - 2*xLi32 - 10*zeta3
2073 > (18 - 81*x + 6*Pi2*x + 123*x2 - 6*Pi2*x2 - 60*x3 +
2074 > 4*Pi2*x3 - 6*(-2 + 3*x - 3*x2 + 2*x3)*xln1m2 -33*x*xln +
2075 > 15*x2*xln - 24*x3*xln - 9*x*xln2 + 9*x2*xln2 -
2076 > 12*x3*xln2 - 12*x1m*xln1m*(-1 + 2*x2 + 2*xln - x*xln +
2077 > 2*x2*xln) - 24*xLi2)/9.
2079 > (2*(-67 + 2*Pi2 + x*(64 + x*(-91 + 3*Pi2 + 94*x)) +
2080 > x1m*(7+x*(-5+16*x))*xln1m -3*x1m*(2+ x*(-1+2*x))*xln1m2 -
2081 > 20*xln - 3*x*xln*(13 + 16*x*x1p - 3*x1p*xln) +
2082 > 6*x1p*(2+x+2*x2)*xln*xln1p+6*x1p*(2+x+2*x2)*xLi2m))/9.
2083 AGF2(I2) = CF*TRNF*WGFCFT+CF**2* WGFCF2+CF*CG*WGFCFG
2084 AFG2(I2) = CF*TRNF*WFGCFT +CG*TRNF*WFGCGT
2088 ZGF2=-28./27.*Cf**2+94./27.*Cf*Cg -52./27.*Cf*TrNf
2089 ZFG2= 37./27.*Cf*TrNf + 35./54.*Cg*TrNf
2090 ZQQB=1.43862321154902*(Cf**2-0.5*Cf*Cg)
2095 XG2 = (LOG(1./(1.-X)) + 1.) ** 2
2096 FFG (IS, IS) = FG2(NX) * DXTZ(I) * XG2
2097 GGF (IS, IS) = GF2(NX) * DXTZ(I) * XG2
2098 PNS (IS, IS) =PNSM(NX) * DXTZ(I)
2104 XG2 = (LOG(1./XM1) + 1.) ** 2
2111 TEM = (FF2(IZ) * (1.- DT) + FF2(IZ+1) * DT) / XM1 / XY
2112 FFG (IX, IY) = TEM * DXTZ(IY)
2113 TEM = (FG2(IZ) * (1.- DT) + FG2(IZ+1) * DT) * XG2 / XY
2114 FFG (IS, IT) = TEM * DXTZ(IY)
2115 TEM = (GF2(IZ) * (1.- DT) + GF2(IZ+1) * DT) * XG2 / XY
2116 GGF (IS, IT) = TEM * DXTZ(IY)
2117 TEM = (GG2(IZ) * (1.- DT) + GG2(IZ+1) * DT) / XM1 / XY
2118 GGF (IX, IY) = TEM * DXTZ(IY)
2119 TEM = (PNSP(IZ) * (1.- DT) + PNSP(IZ+1) * DT) / XM1
2120 PNS (IX, IY) = TEM * DXTZ(IY)
2121 TEM = (PNSM(IZ) * (1.- DT) + PNSM(IZ+1) * DT) / XM1
2122 PNS (IS, IT) = TEM * DXTZ(IY)
2127 SUBROUTINE CtLhTRNLAM (IRDR, NF, IACT, IRT)
2128 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2129 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
2130 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2132 DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 /
2133 DATA IR1, SML / 0, 1.E-5 /
2139 IF (JACT .GT. 0) THEN
2142 ALM = LOG (THMS/VLAM)
2148 THMS = MAX (THMS, SML)
2149 BLM = LOG (THMS/VLAM)
2151 IF (VLAM .GE. 0.7 * THMS) THEN
2152 IF (JACT . EQ. 1) THEN
2161 IF (ALM .GE. BLM) THEN
2162 WRITE (*, *) 'CtLhTRNLAM has ALM >= BLM: ', ALM, BLM
2163 WRITE (*, *) 'I do not know how to continue'
2167 ERR = RERR * LOG (VMULM)
2168 WLLN = CtLhQZBRNT (CtLhZBRLAM, ALM, BLM, ERR, IR1)
2169 ALAM(N1) = THMS / EXP (WLLN)
2170 IF (IR1 .NE. 0) THEN
2171 WRITE (*, *) 'CtLhQZBRNT failed in CtLhTRNLAM; ',
2172 > 'NF, VLAM =', NF, VLAM
2173 WRITE (*, *) 'I do not know how to continue'
2178 SUBROUTINE CtLhUPC (A, La, UpA)
2179 CHARACTER A*(*), UpA*(*), C*(1)
2183 If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!'
2184 Ld = ICHAR('A')-ICHAR('a')
2188 IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN
2189 UpA (I:I) = CHAR(Ichar(c) + ld)
2200 SUBROUTINE CtLhWARNI (IWRN, MSG, NMVAR, IVAB,
2202 CHARACTER*(*) MSG, NMVAR
2209 PRINT '(1X,A/1X, 2A,I10)', MSG, NMVAR,' = ', IV
2210 IF (IACT .EQ. 1) THEN
2211 PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX
2214 If (Iw .LT. Nmax) Then
2215 PRINT '(1X,A/1X, 2A,I10)', MSG, NMVAR,' = ', IV
2216 Elseif (Iw .Eq. Nmax) Then
2217 Print '(/A/)', 'CtLhWARNI Severe Warning: Too many errors'
2222 SUBROUTINE CtLhWARNR (IWRN, MSG, NMVAR, VARIAB,
2224 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2225 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2226 CHARACTER*(*) MSG, NMVAR
2232 PRINT '(1X, A/1X,2A,1PD16.7)', MSG, NMVAR, ' = ', VR
2233 IF (IACT .EQ. 1) THEN
2234 PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX
2237 If (Iw .LT. Nmax) Then
2238 PRINT '(I5, 2A/1X,2A,1PD16.7)', IW, ' ', MSG,
2240 Elseif (Iw .Eq. Nmax) Then
2241 Print '(/A/)', 'CtLhWARNR Severe Warning: Too many errors'
2246 SUBROUTINE CtLhXARRAY
2247 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2249 PARAMETER (D0 = 0.0, D10=10.0)
2250 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
2251 PARAMETER (MXPN = MXF * 2 + 2)
2252 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
2253 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
2255 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
2256 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
2257 COMMON / LhCtHINTEC / GH(NDG, MXX)
2258 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2259 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
2260 DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG)
2261 DATA F12, F22, F32 / 1D0, 1D0, 1D0 /
2262 DATA (G1(I,NDH), G2(I,1), I=1,NDG) / 0.0,0.0,0.0,0.0,0.0,0.0 /
2270 DXTZ(I) = CtLhDXDZ(Z) / X
2275 IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L
2281 DXTZ(NX) = CtLhDXDZ(1.D0)
2287 ELY(I) = LOG(1D0 - XV(I))
2289 ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3)
2293 XY = XV(IX) / XV(IY)
2294 ZZ (IX, IY) = CtLhZFRMX (XY)
2295 ZZ (NX-IX+1, NX-IY+1) = XY
2298 IF (I .NE. NX-1) THEN
2305 DET = F11*F22*F33 + F21*F32*F13 + F31*F12*F23
2306 > - F31*F22*F13 - F21*F12*F33 - F11*F32*F23
2307 IF (ABS(DET) .LT. PUNY) THEN
2308 Msg='Determinant close to zero; will be arbitrarily set to:'
2309 CALL CtLhWARNR(IWRN, Msg, 'DET', PUNY, D0, D0, 0)
2312 G2(1,2) = (F22*F33 - F23*F32) / DET
2313 G2(1,3) = (F32*F13 - F33*F12) / DET
2314 G2(1,4) = (F12*F23 - F13*F22) / DET
2315 G2(2,2) = (F23*F31 - F21*F33) / DET
2316 G2(2,3) = (F33*F11 - F31*F13) / DET
2317 G2(2,4) = (F13*F21 - F11*F23) / DET
2318 G2(3,2) = (F21*F32 - F22*F31) / DET
2319 G2(3,3) = (F31*F12 - F32*F11) / DET
2320 G2(3,4) = (F11*F22 - F12*F21) / DET
2321 B2 = LOG (XV(I+2)/XV(I))
2322 B3 = XV(I) * (B2 - 1.) + XV(I+2)
2323 GH (1,I) = B2 * G2 (2,2) + B3 * G2 (3,2)
2324 GH (2,I) = B2 * G2 (2,3) + B3 * G2 (3,3)
2325 GH (3,I) = B2 * G2 (2,4) + B3 * G2 (3,4)
2331 ElseIF (I .EQ. NX-1) THEN
2334 GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0
2342 A(K) = XA(I+1, 0) - XA(I, 0)
2344 A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK)
2350 TEM = TEM + A(L) * GB(L,J,I)
2364 FUNCTION CtLhXFRMZ (Z)
2365 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2367 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2368 PARAMETER (MXX = 105)
2369 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2370 COMMON / LhCtINVERT / ZA
2372 DATA TEM, RER / D1, 1E-3 /
2373 DATA ZLOW, ZHIGH, IWRN2 / -10.0, 1.00002, 0 /
2376 IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN
2377 XLA = LOG (XMIN) * 1.5
2379 TEM = CtLhZBRNT (CtLhZFXL, XLA, XLB, EPS, IRT)
2381 CALL CtLhWARNR (IWRN2, 'Z out of range in CtLhXFRMZ, X set=0.',
2382 > 'Z', Z, ZLOW, ZHIGH, 1)
2385 CtLhXFRMZ = EXP(TEM)
2388 FUNCTION CtLhxLi(n,x)
2390 integer NCUT, i,n,m3
2391 real*8 CtLhxLi,Out,x,pi2by6,zeta3,c1,c2
2394 dimension c1(2:m3),c2(2:m3)
2396 data c1/0.75,-0.5833333333333333d0,0.454861111111111d0,
2397 > -0.3680555555555555d0,0.3073611111111111d0,
2398 > -0.2630555555555555d0,0.2294880243764172d0/
2399 data c2/-0.5d0,0.5d0,-0.4583333333333333d0,0.416666666666666d0,
2400 > -0.3805555555555555d0,0.35d0,-0.3241071428571428d0/
2401 data zeta3,pi2by6 /1.20205690315959d0,1.64493406684823d0/
2405 if (abs(x).gt.r) then
2406 PRINT *,'Li: x out of range (-1,1) , x=',x
2410 PRINT *,'Polylogarithm Li undefined for n=',n
2412 elseif (n.eq.0) then
2414 elseif (n.eq.1) then
2416 elseif (n.eq.2) then
2417 !Calculate dilogarithm
2418 !separately for x<0.5 and x>0.5
2419 if (x.ge.(-0.5).and.x.le.0.5) then
2426 elseif (x.eq.0) then
2428 elseif(x.gt.0.5) then !n.eq.2,x>0.5
2430 L = pi2by6 - dlog(x)*dlog(xt)
2437 elseif (x.lt.(-0.5)) then
2440 do while (i.le.NCUT)
2447 elseif (n.eq.3.and.x.ge.0.8) then !use the expansion of Li3 near x=1
2448 L=zeta3+pi2by6*dlog(x)
2452 L=L+(c1(i)+c2(i)*xln1m)*xt**i
2455 else !n>3 or x=3,x<0.8
2459 L=L+r/dble(i)**dble(n)
2465 FUNCTION CtLhZBRLAM (WLLN)
2466 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2467 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2469 TEM1 = 1./ CtLhALPQCD(JRDR, N1, WMULM, I)
2470 TEM2 = 1./ CtLhALPQCD(JRDR, N, VMULM, I)
2471 CtLhZBRLAM = TEM1 - TEM2
2473 FUNCTION CtLhZBRNT(FUNC, X1, X2, TOL, IRT)
2474 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2475 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
2483 IF(FB*FA.GT.0.) THEN
2484 PRINT *, 'Root must be bracketed for CtLhZBRNT. Set = 0'
2491 IF(FB*FC.GT.0.) THEN
2497 IF(ABS(FC).LT.ABS(FB)) THEN
2505 TOL1=2.*EPS*ABS(B)+0.5*TOL
2507 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
2511 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
2519 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
2520 Q=(Q-1.)*(R-1.)*(S-1.)
2524 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
2537 IF(ABS(D) .GT. TOL1) THEN
2544 PRINT *, 'CtLhZBRNT exceeding maximum iterations.'
2549 FUNCTION CtLhZFRMX (XX)
2550 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2552 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2553 PARAMETER (MXX = 105)
2554 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2555 DATA IWRN1, HUGE, TINY / 0, 1.E35, 1.E-35 /
2556 F(X) = (XCR-XMIN) * LOG (X/XMIN) + LOG (XCR/XMIN) * (X-XMIN)
2557 D(X) = (XCR-XMIN) / X + LOG (XCR/XMIN)
2559 IF (X .GE. XMIN) THEN
2561 ElseIF (X .GE. D0) THEN
2565 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhZFRMX'
2566 > , 'X', X, TINY, HUGE, 1)
2574 IF (X .GE. XMIN) THEN
2576 ElseIF (X .GE. D0) THEN
2580 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhDZDX '
2581 > , 'X', X, TINY, HUGE, 1)
2588 FUNCTION CtLhZFXL (XL)
2589 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2590 COMMON / LhCtINVERT / ZA
2592 TT = CtLhZFRMX (X) - ZA