]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.3.1/EVLCTEQ.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
pdfs for LHAPDF version 5.3.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / EVLCTEQ.f
CommitLineData
4e9e3152 1 SUBROUTINE CtLhALFSET (QS, ALFS)
2 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3 EXTERNAL CtLhRTALF
4 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
5 DATA ALAM, BLAM, ERR / 0.01, 10.0, 0.02 /
6 ALFST = ALFS
7 CALL CtLhParQcd (2, 'ORDR', ORDR, IR1)
8 JORD = Int(ORDR)
9 NEFF = LhCtNFL(QS)
10 EFLLN = CtLhQZBRNT (CtLhRTALF, ALAM, BLAM, ERR, IR2)
11 EFFLAM = QS / EXP (EFLLN)
12 CALL CtLhSETL1 (NEFF, EFFLAM)
13 END
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
18 LOGICAL SET
19 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
20 DATA IW1, IW2 / 2*0 /
21 IF(.NOT.SET) CALL CtLhLAMCWZ
22 NEFF = LhCtNFL(AMU)
23 ALM = ALAM(NEFF)
24 CtLhALPI = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT)
25 IF (IRT .EQ. 1) THEN
26 CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI', 'AMU', AMU,
27 > ALM, BIG, 1)
28 ELSEIF (IRT .EQ. 2) THEN
29 CALL CtLhWARNR (IW2, 'CtLhALPI > 3; Be aware!', 'CtLhALPI',
30 > CtLhALPI, D0, D1, 0)
31 ENDIF
32 RETURN
33 END
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)
38 IRT = 0
39 IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN
40 print *,
41 > 'Order out of range in CtLhALPQCD: IRDR = ', IRDR
42 STOP
43 ENDIF
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
46 RM2 = RML**2
47 IF (RM2 .LE. 1.) THEN
48 IRT = 1
49 CtLhALPQCD = 99.
50 RETURN
51 ENDIF
52 ALN = LOG (RM2)
53 AL = 4.d0/ B0 / ALN
54 IF (IRDR .GE. 2) AL = AL * (1.d0-B1*LOG(ALN) / ALN / B0**2)
55 IF (AL .GE. 3.) THEN
56 IRT = 2
57 ENDIF
58 CtLhALPQCD = AL
59 RETURN
60 END
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
65 LOGICAL 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'
69 CtLhAMHATF = 0
70 ELSE
71 CtLhAMHATF = AMHAT(I)
72 ENDIF
73 RETURN
74 END
75 FUNCTION CtLhDXDZ (Z)
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 /
79 ZZ = Z
80 X = CtLhXFRMZ (ZZ)
81 TEM = CtLhDZDX (X)
82 IF (TEM .NE. D0) THEN
83 TMP = D1 / TEM
84 Else
85 CALL CtLhWARNR(IWRN, 'CtLhDXDZ singular in CtLhDXDZ; set=HUGE',
86 > 'Z', Z, D0, D1, 0)
87 TMP = HUGE
88 EndIf
89 CtLhDXDZ = TMP
90 RETURN
91 END
92 SUBROUTINE CtLhEVLPAR (IACT, NAME, VALUE, IRET)
93 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
94 CHARACTER*(*) NAME
95 IRET = 1
96 IF (IACT .EQ. 0) THEN
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 ' /)
103 IRET = 4
104 ElseIF (IACT .EQ. 1) THEN
105 CALL CtLhEVLSET (NAME, VALUE, IRET)
106 Else
107 print *,'fatal evlpar'
108 stop
109 EndIf
110 RETURN
111 END
112 SUBROUTINE CtLhEVLSET (NAME, VALUE, IRET)
113 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
114 LOGICAL LSTX
115 CHARACTER*(*) NAME
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
123 IRET = 1
124 IF (NAME .EQ. 'QINI') THEN
125 IF (VALUE .LE. 0) GOTO 12
126 QINI = VALUE
127 ElseIF (NAME .EQ. 'IPD0') THEN
128 ITEM = NINT(VALUE)
129 IF (Item .Eq. 10 .or. Item .Eq. 11) GOTO 12
130 IPD0 = ITEM
131 ElseIF (NAME .EQ. 'IHDN') THEN
132 ITEM = NINT(VALUE)
133 IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12
134 IHDN = ITEM
135 ElseIF (NAME .EQ. 'QMAX') THEN
136 IF (VALUE .LE. QINI) GOTO 12
137 QMAX = VALUE
138 ElseIF (NAME .EQ. 'IKNL') THEN
139 ITMP = NINT(VALUE)
140 ITEM = ABS(ITMP)
141 IF (ITEM.NE.1.AND.ITEM.NE.2) GOTO 12
142 IKNL = ITMP
143 ElseIF (NAME .EQ. 'XCR') THEN
144 IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12
145 XCR = VALUE
146 LSTX = .FALSE.
147 ElseIF (NAME .EQ. 'XMIN') THEN
148 IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12
149 XMIN = VALUE
150 LSTX = .FALSE.
151 ElseIF (NAME .EQ. 'NX') THEN
152 ITEM = NINT(VALUE)
153 IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12
154 NX = ITEM
155 LSTX = .FALSE.
156 ElseIF (NAME .EQ. 'NT') THEN
157 ITEM = NINT(VALUE)
158 IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12
159 NT = ITEM
160 ElseIF (NAME .EQ. 'JT') THEN
161 ITEM = NINT(VALUE)
162 IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12
163 JT = ITEM
164 ElseIF (NAME .EQ. 'NFMX') THEN
165 ITEM = NINT(VALUE)
166 IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12
167 NfMx = ITEM
168 ElseIF (NAME .EQ. 'IPDMOD') THEN
169 ITEM = NINT(VALUE)
170 IF (Abs(Item) .Gt. 1) GOTO 12
171 IpdMod = ITEM
172 ElseIF (NAME .EQ. 'IPTN0') THEN
173 ITEM = NINT(VALUE)
174 IF (ABS(ITEM) .GT. MXF) GOTO 12
175 IPTN0 = ITEM
176 ElseIF (NAME .EQ. 'NUINI') THEN
177 ITEM = NINT(VALUE)
178 IF (ITEM .LE. 0) GOTO 12
179 NuIni = ITEM
180 Else
181 IRET = 0
182 EndIf
183 RETURN
184 12 IRET = 2
185 RETURN
186 END
187 SUBROUTINE CtLhEVOLVE (FINI, IRET)
188 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
189 include 'parmsetup.inc'
190 LOGICAL LSTX
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)
202 DIMENSION QRKP(MXF)
203 DIMENSION JI(-MXF : MXF+1)
204 EXTERNAL LhCtNSRHSP, LhCtNSRHSM, FINI
205 DATA DZER / 0.0 /
206 save nxsave, ntsave, jtsave, ngsave,
207 & xcrsave, xminsave, qinisave, qmaxsave, ientry
208 data ientry / 0 /
209c
210 call getnset(iset)
211c
212 ientry = ientry + 1
213 IRET = 0
214 IF (IHDN .LE. 4) THEN
215 MXVAL = 2
216 ElseIF (IHDN .LE. 6) THEN
217 MXVAL = 3
218 EndIf
219 IF (.NOT. LSTX) CALL CtLhXARRAY
220 CALL CtLhPARPDF (2, 'ALAM', AL, IR)
221 CALL CtLhQARRAY (NINI)
222 NFSN = NFMX + 1
223 KF = 2 * NFMX + 2
224 Nelmt = KF * (Nt+1) * (Nx+1)
225 DO 101 IFLV = -NFMX, NFMX+1
226 JFL = NFMX + IFLV
227 JI(IFLV) = JFL * (NT+1) * (NX+1)
228 101 CONTINUE
229 DO 31 IZ = 1, NX
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))
236 QRKP (IFLV) = A + B
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
240 331 CONTINUE
241 DO 332 IFLV = 1, NINI
242 UPD(JI( IFLV)+IZ+1,iset) =
243 > QRKP(IFLV) - UPD(JI(NFSN)+IZ+1,iset)/NINI
244 332 CONTINUE
245 31 CONTINUE
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
250 NITR = NTN (ICNT)
251 DT = DTN (ICNT)
252 TIN = TLN (ICNT)
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))
260 IF (IFLV .LE. MXVAL)
261 > CALL CtLhNSEVL (LhCtNSRHSM, IKNL, NX, NITR, JT, DT, TIN,
262 > NEFF, UPD(JI(-IFLV)+2,iset), UPD(JI(-IFLV)+1,iset))
263 DO 55 IS = 0, NITR
264 DO 56 IX = 0, NX
265 TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV),iset)
266 TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN),iset) / NEFF
267 TP = TP + TS
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
272 TM = MAX (TM, DZER)
273 TP = MAX (TP, TM)
274 EndIf
275 Else
276 TM = 0.
277 EndIf
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.
280 56 CONTINUE
281 55 CONTINUE
282333 CONTINUE
283 DO 334 IFLV = NEFF + 1, NFMX
284 DO 57 IS = 0, NITR
285 DO 58 IX = 0, NX
286 UPD(JI( IFLV) + IS*(NX+1) + IX + 1,iset) = 0
287 UPD(JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = 0
288 58 CONTINUE
289 57 CONTINUE
290 334 CONTINUE
291 88 CONTINUE
292 IF (NFMX .EQ. NEFF) GOTO 21
293 DO 335 IFLV = -NFMX, NFMX+1
294 JI(IFLV) = JI(IFLV) + NITR * (NX+1)
295 335 CONTINUE
296 CALL CtLhHQRK (NX, TT, NEFF+1, UPD(JI(0)+2,iset),
297 > UPD(JI(NEFF+1)+2,iset))
298 DO 32 IZ = 1, NX
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)
301 > + QRKP (NEFF+1)
302 VS00 = UPD (JI(NFSN)+IZ+1,iset) / (NEFF+1)
303 UPD(JI( NEFF+1) + IZ + 1,iset) = QRKP(NEFF+1) - VS00
304 DO 321 IFL = 1, NEFF
305 A = UPD(JI( IFL)+IZ+1,iset)
306 B = UPD(JI(-IFL)+IZ+1,iset)
307 QRKP(IFL) = A + B
308 UPD(JI( IFL)+IZ+1,iset) = QRKP(IFL) - VS00
309 IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1,iset) = A - B
310 321 CONTINUE
311 32 CONTINUE
312 21 CONTINUE
313 if(ientry .eq. 1) then
314 nxsave = nx
315 ntsave = nt
316 jtsave = jt
317 ngsave = ng
318 xcrsave = xcr
319 xminsave = xmin
320 qinisave = qini
321 qmaxsave = qmax
322 endif
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,
332 & qini, qmax, ientry
333669 format(1x,'evolve.f: nx,nt,jt,ng=',4i4,
334 & ' xcr,xmin=',2f9.6,
335 & ' qini, qmax',f7.4,1x,e12.5,' ientry=',i6)
336 nxsave = nx
337 ntsave = nt
338 jtsave = jt
339 ngsave = ng
340 qinisave = qini
341 qmaxsave = qmax
342 xcrsave = xcr
343 xminsave = xmin
344 endif
345 Return
346 End
347*
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)
351 PARAMETER (MX = 3)
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 /
355 IR = 0
356 X = XV
357 ERR = 0.
358 ANX = NX
359 CtLhFINTRP = 0.
360 IF (NX .LT. 1) THEN
361 CALL CtLhWARNI(IW1, 'Nx < 1, error in CtLhFINTRP.',
362 > 'NX', NX, 1, 256, 1)
363 IR = 1
364 RETURN
365 ELSE
366 MNX = MIN(NX+1, MX)
367 ENDIF
368 IF (DX .LE. 0) THEN
369 CALL CtLhWARNR(IW3, 'DX < 0, error in CtLhFINTRP.',
370 > 'DX', DX, D0, D1, 1)
371 IR = 2
372 RETURN
373 ENDIF
374 XM = X0 + DX * NX
375 IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN
376 CALL CtLhWARNR(IW5,
377 > 'X out of range in CtLhFINTRP, Extrapolation used.',
378 > 'X',X,X0,XM,1)
379 IR = 3
380 ENDIF
381 TX = (X - X0) / DX
382 IF (TX .LE. 1.) THEN
383 IX = 0
384 ELSEIF (TX .GE. ANX-1.) THEN
385 IX = NX - 2
386 ELSE
387 IX = Int(TX)
388 ENDIF
389 DDX = TX - IX
390 CALL CtLhRATINT (XX, FF(IX), MNX, DDX, TEM, ERR)
391 CtLhFINTRP = TEM
392 RETURN
393 END
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)
398 external f
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/
439 TOLABS=AERR
440 TOLREL=RERR
441
442 CtLhGausInt=0.
443 NLIMS=2
444 XLIMS(1)=XL
445 XLIMS(2)=XR
44610 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0
447 BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0
448 TVAL=0.
449 DO 15 I=1,3
45015 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
451 TVAL=TVAL*AA
452 DO 25 J=1,4
453 VAL=0.
454 DO 20 I=PTR(J),PTR(J)-1+NORD(J)
45520 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
456 VAL=VAL*AA
457 TOL=MAX(TOLABS,TOLREL*ABS(VAL))
458 IF (ABS(TVAL-VAL).LT.TOL) THEN
459 CtLhGausInt=CtLhGausInt+VAL
460 NLIMS=NLIMS-2
461 IF (NLIMS.NE.0) GO TO 10
462 RETURN
463 END IF
46425 TVAL=VAL
465 IF (NMAX.EQ.2) THEN
466 CtLhGausInt=VAL
467 RETURN
468 END IF
469 IF (NLIMS.GT.(NMAX-2)) THEN
470 write(*,50) CtLhGausInt,NMAX,BB-AA,BB+AA
471 RETURN
472 END IF
473 XLIMS(NLIMS+1)=BB
474 XLIMS(NLIMS+2)=BB+AA
475 XLIMS(NLIMS)=BB
476 NLIMS=NLIMS+2
477 GO TO 10
47850 FORMAT (' CtLhGausInt FAILS, CtLhGausInt,NMAX,XL,XR=',
479 > G15.7,I5,2G15.7)
480 END
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)
490 DZ = 1D0 / (NX-1)
491 DO 20 I = 1, NX-2
492 NP = NX - I + 1
493 TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2)
494 DO 30 KZ = 3, NP
495 IY = I + KZ - 1
496 W = XA(I,1) / XA(IY,1)
497 G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W)
498 30 CONTINUE
499 HTEM = CtLhSMPSNA (NP-2, DZ, G(3), ERR)
500 TEM1 = F(I) * ELY(I)
501 H(I) = TEM + HTEM + TEM1
502 20 CONTINUE
503 H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0))
504 H(NX) = 0
505 RETURN
506 END
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
512 11 CONTINUE
513 DO 230 IZ = 1, NX
514 IF (NX .GT. 1) THEN
515 F(IZ) = 0
516 GOTO 230
517 EndIf
518 230 CONTINUE
519 RETURN
520 END
521 SUBROUTINE CtLhINTEGR (NX, M, F, G, IR)
522 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
523 CHARACTER MSG*80
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 /
532 IRR = 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)
536 IRR = 1
537 EndIf
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)
541 IRR = 2
542 EndIf
543 G(NX) = 0D0
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)
546 IF (M .EQ. 0) THEN
547 G(NX-1) = TEM
548 Else
549 G(NX-1) = TEM * XA(NX-1, M)
550 EndIf
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)
554 IF (M .EQ. 0) THEN
555 G(I) = TEM
556 Else
557 G(I) = TEM * XA(I, M)
558 EndIf
559 10 CONTINUE
560 TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3)
561 IF (M .EQ. 0) THEN
562 G(1) = TEM
563 Else
564 G(1) = TEM * XA(1, M)
565 EndIf
566 IR = IRR
567 RETURN
568 END
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 /
575 IRT = 0
576 TRNF = TR * NFL
577 X = XX
578 IF (X .LE. 0. .OR. X .GE. 1.) THEN
579 CALL CtLhWARNR(IWRN, 'X out of range in CtLhKERNEL', 'X', X,
580 > D0, D1, 1)
581 IRT = 1
582 RETURN
583 EndIf
584 XI = 1./ X
585 X2 = X ** 2
586 XM1I = 1./ (1.- X)
587 XP1I = 1./ (1.+ X)
588 XLN = LOG (X)
589 XLN2 = XLN ** 2
590 XLN1M = LOG (1.- 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)
601 FG1 = CF * FGP * X
602 GF1 = 2.* TRNF * GFP
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)
616 > + FFM * 2.* SPEN2
617 FFCFG = 14./3.* (1.-X)
618 > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.)
619 > - FFM * SPEN2
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
624 > - 2.* X * XLN1M
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
629 > - PI2/6. + 0.5)
630 > + FGM * SPEN2
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.)
639 > + GFM * 2. * SPEN2
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
643 > - GGP * 20./9.
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.)
647 > + GGM * 2.* SPEN2
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.)
653 XG2 = XLG ** 2
654 FF2 = FF2 * X * (1.- X)
655 FG2 = FG2 * X / XG2
656 GF2 = GF2 * X / XG2
657 GG2 = GG2 * X * (1.- X)
658 RETURN
659 END
660 SUBROUTINE CtLhLAMCWZ
661 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
662 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
663 LOGICAL SET
664 CALL CtLhSETL1 (NF, AL)
665 END
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
670 LOGICAL SET
671 CHARACTER ONECH*(1)
672 ONECH = '0'
673 IASC0 = ICHAR(ONECH)
674 NAME = NNAME
675 LhCtNAMQCD=0
676 IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR.
677 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') )
678 2 LhCtNAMQCD=1
679 IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR.
680 1 (NAME .EQ. '# FL') )
681 2 LhCtNAMQCD=2
682 DO 10 I=1, 9
683 IF (NAME .EQ. 'M'//CHAR(I+IASC0))
684 1 LhCtNAMQCD=I+2
68510 CONTINUE
686 DO 20 I= 0, NF
687 IF (NAME .EQ. 'LAM'//CHAR(I+IASC0))
688 1 LhCtNAMQCD=I+13
68920 CONTINUE
690 IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') LhCtNAMQCD = 24
691 RETURN
692 END
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
697 LOGICAL SET
698 IF (.NOT. SET) CALL CtLhLAMCWZ
699 LhCtNFL = NF - NHQ
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
703 LhCtNFL = I
704 ELSE
705 GOTO 20
706 ENDIF
70710 CONTINUE
70820 RETURN
709 END
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)
719 external rhs
720 DDT = DT / JT
721 IF (NX .GT. MXX) THEN
722 WRITE (*,*) 'Nx =', NX, ' greater than Max pts in CtLhNSEVL.'
723 STOP 'Program stopped in CtLhNSEVL'
724 EndIf
725c ++ remove unused quantities (jcp)
726c ++ TMD = TIN + DT * NT / 2.
727c ++ AMU = EXP(TMD)
728c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
729c ++ TLAM = TMD - TEM
730 DO 9 IX = 1, NX
731 UN(IX, 0) = U0(IX)
732 9 CONTINUE
733 UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1)
734 TT = TIN
735 DO 10 IZ = 1, NX
736 Y0(IZ) = U0(IZ)
737 10 CONTINUE
738 DO 20 IS = 1, NT
739 DO 202 JS = 1, JT
740 IRND = (IS-1) * JT + JS
741 IF (IRND .EQ. 1) THEN
742 CALL RHS (TT, Neff, Y0, F0)
743 DO 250 IZ = 1, NX
744 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
745 250 CONTINUE
746 TT = TT + DDT
747 CALL RHS (TT, NEFF, Y0, F1)
748 DO 251 IZ = 1, NX
749 Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
750 251 CONTINUE
751 Else
752 CALL RHS (TT, NEFF, Y1, F1)
753 DO 252 IZ = 1, NX
754 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
755 252 CONTINUE
756 TT = TT + DDT
757 CALL RHS (TT, NEFF, YP, FP)
758 DO 253 IZ = 1, NX
759 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
760 F0(IZ) = F1(IZ)
761 253 CONTINUE
762 EndIf
763 202 CONTINUE
764 DO 260 IZ = 1, NX
765 UN (IZ, IS) = Y1(IZ)
766 260 CONTINUE
767 UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
768 20 CONTINUE
769 RETURN
770 END
771 SUBROUTINE LhCtNSRHSM (TT, NEFF, FI, FO)
772 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
773 LOGICAL LSTX
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)
785 S = EXP(TT)
786 Q = AL * EXP (S)
787 CPL = CtLhALPI(Q)
788 CPL2= CPL ** 2 / 2. * S
789 CPL = CPL * S
790 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
791 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
792 CALL CtLhHINTEG (NX, FI, WH)
793 DO 230 IZ = 1, NX
794 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
795 FO(IZ) = CPL * FO(IZ)
796 230 CONTINUE
797 IF (IKNL .EQ. 2) THEN
798 DZ = 1./ (NX - 1)
799 DO 21 IX = 1, NX-1
800 NP = NX - IX + 1
801 IS = NP
802 DO 31 KZ = 2, NP
803 IY = IX + KZ - 1
804 IT = NX - IY + 1
805 XY = ZZ (IS, IT)
806 G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX))
807 31 CONTINUE
808 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
809 TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2
810 FO(IX) = FO(IX) + TMP2
811 21 CONTINUE
812 EndIf
813 RETURN
814 END
815 SUBROUTINE LhCtNSRHSP (TT, NEFF, FI, FO)
816 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
817 LOGICAL LSTX
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)
829 S = EXP(TT)
830 Q = AL * EXP (S)
831 CPL = CtLhALPI(Q)
832 CPL2= CPL ** 2 / 2. * S
833 CPL = CPL * S
834 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
835 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
836 CALL CtLhHINTEG (NX, FI, WH)
837 DO 230 IZ = 1, NX
838 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
839 FO(IZ) = CPL * FO(IZ)
840 230 CONTINUE
841 IF (IKNL .EQ. 2) THEN
842 DZ = 1./ (NX - 1)
843 DO 21 IX = 1, NX-1
844 NP = NX - IX + 1
845 DO 31 KZ = 2, NP
846 IY = IX + KZ - 1
847 XY = ZZ (NX-IX+1, NX-IY+1)
848 G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX))
849 31 CONTINUE
850 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
851 TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + ZQQB)) * CPL2
852 FO(IX) = FO(IX) + TMP2
853 21 CONTINUE
854 EndIf
855 RETURN
856 END
857 FUNCTION CtLhPARDIS (IPRTN, XX, QQ)
858 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
859 include 'parmsetup.inc'
860 Character Msg*80
861 LOGICAL LSTX
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)
867 parameter(nqvec = 4)
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
880 data nxsave / 0 /
881 save xvpow, nxsave
882 save xlast, qlast
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
886c
887 call getnset(iset)
888c
889 if(nx .ne. nxsave) then
890 xvpow(0) = 0D0
891 do i = 1, nx
892 xvpow(i) = xv(i)**xpow
893 enddo
894 nxsave = nx
895 endif
896
897 X = XX
898 Q = QQ
899
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')
904
905c enforce threshold early to improve speed...
906 ii = iabs(IPRTN)
907 if(ii .ne. 0) then
908 if(QQ .lt. VALQMS(ii) ) then
909 ctlhpardis = 0.d0
910 return
911 endif
912 endif
913
914c force pardis = 0.0d0 at exactly =1.0d0 - added mrw 10/May/06
915 if(xx .eq. 1.0d0) then
916 ctlhpardis = 0.0d0
917 return
918 endif
919
920c skip the initialization in x if same as in the previous call.
921 if(x .eq. xlast) goto 100
922 xlast = x
923
924 JLx = -1
925 JU = Nx+1
926 11 If (JU-JLx .GT. 1) Then
927 JM = (JU+JLx) / 2
928 If (X .Ge. XV(JM)) Then
929 JLx = JM
930 Else
931 JU = JM
932 Endif
933 Goto 11
934 Endif
935 If (JLx .LE. -1) Then
936 Print '(A,1pE12.4)','Severe error: x <= 0 in ParDis x=', x
937 Stop
938 ElseIf (JLx .Eq. 0) Then
939 Jx = 0
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
943 Jx = JLx - 1
944 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
945 Jx = JLx - 2
946 Else
947 Print '(A,1pE12.4)','Severe error: x > 1 in ParDis x=', x
948 Stop
949 Endif
950 ss = x**xpow
951 If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
952 svec1 = xvpow(jx)
953 svec2 = xvpow(jx+1)
954 svec3 = xvpow(jx+2)
955 svec4 = xvpow(jx+3)
956 s12 = svec1 - svec2
957 s13 = svec1 - svec3
958 s23 = svec2 - svec3
959 s24 = svec2 - svec4
960 s34 = svec3 - svec4
961 sy2 = ss - svec2
962 sy3 = ss - svec3
963 const1 = s13/s23
964 const2 = s12/s23
965 const3 = s34/s23
966 const4 = s24/s23
967 s1213 = s12 + s13
968 s2434 = s24 + s34
969 sdet = s12*s34 - s1213*s2434
970 tmp = sy2*sy3/sdet
971 const5 = (s34*sy2-s2434*sy3)*tmp/s12
972 const6 = (s1213*sy2-s12*sy3)*tmp/s34
973 EndIf
974
975100 continue
976
977c skip the initialization in q if same as in the previous call.
978 if(q .eq. qlast) goto 110
979 qlast = q
980
981 tt = log(log(Q/Al))
982
983 JLq = -1
984 JU = NT+1
985 12 If (JU-JLq .GT. 1) Then
986 JM = (JU+JLq) / 2
987 If (Q .GE. QV(JM)) Then
988 JLq = JM
989 Else
990 JU = JM
991 Endif
992 Goto 12
993 Endif
994 If (JLq .LE. 0) Then
995 Jq = 0
996 If (JLq .LT. 0) Then
997 Msg = 'Q < Q0 in ParDis; extrapolation used!'
998 CALL CtLhWARNR (IWRN2, Msg, 'Q', Q, Qini, 1D0, 1)
999 EndIf
1000 Elseif (JLq .LE. Nt-2) Then
1001 Jq = JLq - 1
1002 Else
1003 Jq = Nt - 3
1004 If (JLq .GE. Nt) Then
1005 Msg = 'Q > Qmax in ParDis; extrapolation used!'
1006 CALL CtLhWARNR (IWRN3, Msg, 'Q', Q, Qmax, 1D0, 1)
1007 Endif
1008 Endif
1009
1010 If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
1011 tvec1 = Tv(jq)
1012 tvec2 = Tv(jq+1)
1013 tvec3 = Tv(jq+2)
1014 tvec4 = Tv(jq+3)
1015 t12 = tvec1 - tvec2
1016 t13 = tvec1 - tvec3
1017 t23 = tvec2 - tvec3
1018 t24 = tvec2 - tvec4
1019 t34 = tvec3 - tvec4
1020 ty2 = tt - tvec2
1021 ty3 = tt - tvec3
1022 tmp1 = t12 + t13
1023 tmp2 = t24 + t34
1024 tdet = t12*t34 - tmp1*tmp2
1025 EndIf
1026
1027110 continue
1028
1029 jtmp = ((IPRTN + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
1030 Do it = 1, nqvec
1031 J1 = jtmp + it*(NX+1)
1032 If (Jx .Eq. 0) Then
1033 fij(1) = 0
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)
1038
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)
1042 Fvec(it) = Fx
1043 Else
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
1051 Endif
1052 enddo
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)
1057 Else
1058 tf2 = fvec(2)
1059 tf3 = fvec(3)
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
1065 EndIf
1066 CtLhPARDIS = ff
1067 Return
1068 End
1069
1070 SUBROUTINE CtLhPARPDF (IACT, NAME, VALUE, IRET)
1071 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1072 CHARACTER NAME*(*), Uname*10
1073 LOGICAL START1
1074 DATA ILEVEL, LRET / 1, 1 /
1075 JRET = IRET
1076 CALL CtLhUPC (NAME, Ln, Uname)
1077 IF (IACT .EQ. 0 .OR. IACT .EQ. 4) then
1078c > IVALUE = NINT (VALUE) !tentatively remove this since it seems not to be used
1079 print *,'Fatal error: iact=',iact
1080 stop
1081 ENDIF
1082 START1 = (IACT .NE. 1) .AND. (IACT .NE. 2)
1083c prepare to remove this stuff, since I think IACT=1 or 2 always
1084 if(start1) then
1085 print *,'Fatal error: start1=',start1
1086 stop
1087 endif
1088 IF (START1) ILEVEL = 1
1089 GOTO (1, 2), ILEVEL
1090 1 START1 = .TRUE.
1091 ILEVEL = 0
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
1097 ILEVEL = ILEVEL + 1
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
1103 ILEVEL = ILEVEL + 1
1104 IF (.NOT. START1) GOTO 1
1105 IF (JRET .EQ. 0) GOTO 10
1106 GOTO 14
1107 10 CONTINUE
1108 11 CONTINUE
1109 12 CONTINUE
1110 13 CONTINUE
1111 14 CONTINUE
1112 15 CONTINUE
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'
1118 stop
1119 EndIf
1120 IRET= JRET
1121 RETURN
1122 END
1123 SUBROUTINE CtLhParQcd(IACT,NAME,VALUE,IRET)
1124 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1125 INTEGER IACT,IRET
1126 CHARACTER*(*) NAME
1127 IRET=1
1128 IF (IACT.EQ.0) THEN
1129 WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ',
1130 > '(i in 1 to 9), LAMi (i in 1 to NFL)'
1131 IRET=4
1132 ELSEIF (IACT.EQ.1) THEN
1133 CALL CtLhQCDSET (NAME,VALUE,IRET)
1134 ELSEIF (IACT.EQ.2) THEN
1135 CALL CtLhQCDGET (NAME,VALUE,IRET)
1136 ELSE
1137 IRET=3
1138 ENDIF
1139 RETURN
1140 END
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)
1147 PARAMETER (MX = 3)
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)
1154 SAVE
1155 DATA LA, LB / 2 * .FALSE. /
1156 LB = .TRUE.
1157 ENTRY CtLhTFF1(ZZ)
1158 LA = .TRUE.
1159 IF (LA .AND. .NOT.LB) THEN
1160 Z = ZZ
1161 X = CtLhXFRMZ (Z)
1162 Else
1163 X = XX
1164 EndIf
1165 IF (X .GE. D1) THEN
1166 CtLhPFF1 = 0
1167 RETURN
1168 ElseIF (X .GE. XMIN) THEN
1169 Z = CtLhZFRMX (X)
1170 TEM = CtLhFINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT)
1171 Else
1172 CALL CtLhPOLIN1 (XL, FF1(1), MX, X, TEM, ERR)
1173 EndIf
1174 IF (LA) THEN
1175 IF (LB) THEN
1176 CtLhPFF1 = TEM / (1.-X)
1177 LB =.FALSE.
1178 Else
1179 CtLhTFF1 = TEM / (1.-X) * CtLhDXDZ(Z)
1180 EndIf
1181 LA =.FALSE.
1182 Else
1183 IF (LB) THEN
1184c +++ something is wrong, since QFF1 and RFF1 are not used.
1185c +++ but this code appears to only be used for extrapolation
1186c +++ to small x, which is unreliable anyway, so ignore for now (jcp)
1187 QFF1 = TEM
1188 LB =.FALSE.
1189 Else
1190 RFF1 = TEM * X / (1.-X)
1191 EndIf
1192 EndIf
1193 RETURN
1194 ENTRY CtLhFNSP (XX)
1195 X = XX
1196 IF (X .GE. D1) THEN
1197 CtLhFNSP = 0.
1198 RETURN
1199 ElseIF (X .GE. XMIN) THEN
1200 Z = CtLhZFRMX (X)
1201 TEM = CtLhFINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT)
1202 Else
1203 CALL CtLhPOLIN1 (XL, PNSP(1), MX, X, TEM, ERR)
1204 EndIf
1205 CtLhFNSP = TEM / (1.- X)
1206 RETURN
1207 ENTRY CtLhFNSM (XX)
1208 X = XX
1209 IF (X .GE. D1) THEN
1210 CtLhFNSM = 0.
1211 RETURN
1212 ElseIF (X .GE. XMIN) THEN
1213 Z = CtLhZFRMX (X)
1214 TEM = CtLhFINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT)
1215 Else
1216 CALL CtLhPOLIN1 (XL, PNSM(1), MX, X, TEM, ERR)
1217 EndIf
1218 CtLhFNSM = TEM / (1.- X)
1219 RETURN
1220 ENTRY CtLhRGG1 (XX)
1221 X = XX
1222 IF (X .GE. D1) THEN
1223 CtLhRGG1= 0 !error corrected? (jcp)
1224 RETURN
1225 ElseIF (X .GE. XMIN) THEN
1226 Z = CtLhZFRMX (X)
1227 TEM = CtLhFINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT)
1228 Else
1229 CALL CtLhPOLIN1 (XL, GG1(1), MX, X, TEM, ERR)
1230 EndIf
1231 IF (LA) THEN
1232 CtLhRGG1 = TEM / X / (1.-X) !error corrected? (jcp)
1233 LA =.FALSE.
1234 Else
1235 IF (LB) THEN
1236 CtLhRGG1 = TEM / X !error corrected? (jcp)
1237 LB =.FALSE.
1238 Else
1239 CtLhRGG1 = TEM / (1.-X)
1240 EndIf
1241 EndIf
1242 RETURN
1243 ENTRY CtLhRFF2 (XX)
1244 X = XX
1245 IF (X .GE. D1) THEN
1246 CtLhRFF2 = 0 !error corrected? (jcp)
1247 RETURN
1248 ElseIF (X .GE. XMIN) THEN
1249 Z = CtLhZFRMX (X)
1250 TEM = CtLhFINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT)
1251 Else
1252 CALL CtLhPOLIN1 (XL, FF2(1), MX, X, TEM, ERR)
1253 EndIf
1254 IF (LA) THEN
1255 CtLhRFF2 = TEM / X / (1.-X) !error corrected? (jcp)
1256 LA =.FALSE.
1257 Else
1258 IF (LB) THEN
1259 CtLhRFF2 = TEM / X !error corrected? (jcp)
1260 LB =.FALSE.
1261 Else
1262 CtLhRFF2 = TEM / (1.-X)
1263 EndIf
1264 EndIf
1265 RETURN
1266 ENTRY CtLhRGG2 (XX)
1267 X = XX
1268 IF (X .GE. D1) THEN
1269 CtLhRGG2 = 0 !error corrected? (jcp)
1270 RETURN
1271 ElseIF (X .GE. XMIN) THEN
1272 Z = CtLhZFRMX (X)
1273 TEM = CtLhFINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT)
1274 Else
1275 CALL CtLhPOLIN1 (XL, GG2(1), MX, X, TEM, ERR)
1276 EndIf
1277 IF (LA) THEN
1278 CtLhRGG2 = TEM / X / (1.-X) !error corrected? (jcp)
1279 LA =.FALSE.
1280 Else
1281 IF (LB) THEN
1282 CtLhRGG2 = TEM / X !error corrected? (jcp)
1283 LB =.FALSE.
1284 Else
1285 CtLhRGG2 = TEM / (1.-X)
1286 EndIf
1287 EndIf
1288 RETURN
1289 END
1290 SUBROUTINE CtLhPOLIN1 (XA,YA,N,X,Y,DY)
1291 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1292 PARAMETER (NMAX=10)
1293 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1294 NS=1
1295 DIF=ABS(X-XA(1))
1296 DO 11 I=1,N
1297 DIFT=ABS(X-XA(I))
1298 IF (DIFT.LT.DIF) THEN
1299 NS=I
1300 DIF=DIFT
1301 ENDIF
1302 C(I)=YA(I)
1303 D(I)=YA(I)
130411 CONTINUE
1305 Y=YA(NS)
1306 NS=NS-1
1307 DO 13 M=1,N-1
1308 DO 12 I=1,N-M
1309 HO=XA(I)-X
1310 HP=XA(I+M)-X
1311 W=C(I+1)-D(I)
1312 DEN=HO-HP
1313 DEN=W/DEN
1314 D(I)=HP*DEN
1315 C(I)=HO*DEN
131612 CONTINUE
1317 IF (2*NS.LT.N-M)THEN
1318 DY=C(NS+1)
1319 ELSE
1320 DY=D(NS)
1321 NS=NS-1
1322 ENDIF
1323 Y=Y+DY
132413 CONTINUE
1325 RETURN
1326 END
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
1335 NCNT = 0
1336 IF (NT .GE. mxq) NT = mxq - 1
1337 S = LOG(QINI/AL)
1338 TINI = LOG(S)
1339 S = LOG(QMAX/AL)
1340 TMAX = LOG(S)
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)
1347 AFL = NfMx
1348 Call CtLhParQcd (1, 'NFL', AFL, Ir)
1349 Iordr = Nint (Ord)
1350 Ifl0 = Nint (Afl0)
1351 Call CtLhSetLam (Ifl0, Al0, Iordr)
1352 NG = NFMX - NINI + 1
1353 QIN = QINI
1354 QOUT = QINI
1355 S = LOG(QIN/AL)
1356 TIN = LOG(S)
1357 TLN(1) = TIN
1358 NTL(1) = 0
1359 QV(0) = QINI
1360 TV(0) = Tin
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)
1366 Else
1367 QOUN = QMAX
1368 EndIf
1369 IF (QOUN-QOUT .LE. 0.0001) THEN
1370 DT = 0
1371 NITR = 0
1372 Else
1373 QOUT = QOUN
1374 S = LOG(QOUT/AL)
1375 TOUT = LOG(S)
1376 TEM = TOUT - TIN
1377 NITR = INT (TEM / DT0) + 1
1378 DT = TEM / NITR
1379 EndIf
1380 DTN (ICNT) = DT
1381 NTN (ICNT) = NITR
1382 TLN (ICNT) = TIN
1383 NTL (ICNT+1) = NTL(ICNT) + NITR
1384 IF (NITR .NE. 0) THEN
1385 DO 205 I = 1, NITR
1386 TV (NTL(ICNT)+I) = TIN + DT * I
1387 S = EXP (TV(NTL(ICNT)+I))
1388 QV (NTL(ICNT)+I) = AL * EXP (S)
1389 205 CONTINUE
1390 EndIf
1391 QIN = QOUT
1392 TIN = TOUT
1393 20 CONTINUE
1394 NCNT = NCNT + 1
1395 NTP = NTL (NG + 1)
1396 ND = NTP - NT
1397 IF (NTP .GE. MXQ) THEN
1398 NT = MXQ - ND - NCNT
1399 GOTO 1
1400 EndIf
1401 NT = NTP
1402 RETURN
1403 END
1404 SUBROUTINE CtLhQCDGET(NAME,VALUE,IRET)
1405 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1406 CHARACTER*(*) NAME
1407 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1408 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1409 COMMON / LhCtCOMQMS / VALQMS(9)
1410 LOGICAL SET
1411 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1412 ICODE = LhCtNAMQCD(NAME)
1413 IRET = 1
1414 IF (ICODE .EQ. 1) THEN
1415 VALUE = AL
1416 ELSEIF (ICODE .EQ. 2) THEN
1417 VALUE = NF
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
1423 VALUE = NORDER
1424 ELSE
1425 IRET=0
1426 ENDIF
1427 END
1428 SUBROUTINE CtLhQCDSET (NAME,VALUE,IRET)
1429 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1430 CHARACTER*(*) NAME
1431 COMMON / LhCtCOMQMS / VALQMS(9)
1432 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1433 LOGICAL SET
1434 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1435 IVALUE = NINT(VALUE)
1436 ICODE = LhCtNAMQCD(NAME)
1437 IF (ICODE .EQ. 0) THEN
1438 IRET=0
1439c print *,'warning empty CtLhQCDSET call: NAME=',
1440c & NAME,' VALUE=',VALUE
1441 ELSE
1442 IRET = 1
1443 SET = .FALSE.
1444 IF (ICODE .EQ. 1) THEN
1445 IF (VALUE.LE.0) GOTO 12
1446 AL=VALUE
1447 ELSEIF (ICODE .EQ. 2) THEN
1448 IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12
1449 NF = IVALUE
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
1461 NORDER = IVALUE
1462 ENDIF
1463 IF (.NOT. SET) CALL CtLhLAMCWZ
1464 ENDIF
1465 RETURN
1466 12 IRET=2
1467 RETURN
1468 END
1469 FUNCTION CtLhQZBRNT(FUNC, X1, X2, TOLIN, IRT)
1470 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1471 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
1472 external func
1473 TOL = ABS(TOLIN)
1474 A=X1
1475 B=X2
1476 FA=FUNC(A)
1477 FB=FUNC(B)
1478 IF(FB*FA.GT.0.) THEN
1479 WRITE (*, *) 'Root must be bracketed for CtLhQZBRNT.'
1480 IRT = 1
1481 ENDIF
1482 FC=FB
1483 DO 11 ITER=1,ITMAX
1484 IF(FB*FC.GT.0.) THEN
1485 C=A
1486 FC=FA
1487 D=B-A
1488 E=D
1489 ENDIF
1490 IF(ABS(FC).LT.ABS(FB)) THEN
1491 A=B
1492 B=C
1493 C=A
1494 FA=FB
1495 FB=FC
1496 FC=FA
1497 ENDIF
1498 TOL1=2.*EPS*ABS(B)+0.5*TOL
1499 XM=.5*(C-B)
1500 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
1501 CtLhQZBRNT=B
1502 RETURN
1503 ENDIF
1504 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
1505 S=FB/FA
1506 IF(A.EQ.C) THEN
1507 P=2.*XM*S
1508 Q=1.-S
1509 ELSE
1510 Q=FA/FC
1511 R=FB/FC
1512 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
1513 Q=(Q-1.)*(R-1.)*(S-1.)
1514 ENDIF
1515 IF(P.GT.0.) Q=-Q
1516 P=ABS(P)
1517 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
1518 E=D
1519 D=P/Q
1520 ELSE
1521 D=XM
1522 E=D
1523 ENDIF
1524 ELSE
1525 D=XM
1526 E=D
1527 ENDIF
1528 A=B
1529 FA=FB
1530 IF(ABS(D) .GT. TOL1) THEN
1531 B=B+D
1532 ELSE
1533 B=B+SIGN(TOL1,XM)
1534 ENDIF
1535 FB=FUNC(B)
153611 CONTINUE
1537 WRITE (*, *) 'CtLhQZBRNT exceeding maximum iterations.'
1538 IRT = 2
1539 CtLhQZBRNT=B
1540 RETURN
1541 END
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)
1546 NS=1
1547 HH=ABS(X-XA(1))
1548 DO 11 I=1,N
1549 H=ABS(X-XA(I))
1550 IF (H.EQ.0.)THEN
1551 Y=YA(I)
1552 DY=0.0
1553 RETURN
1554 ELSE IF (H.LT.HH) THEN
1555 NS=I
1556 HH=H
1557 ENDIF
1558 C(I)=YA(I)
1559 D(I)=YA(I)+TINY
156011 CONTINUE
1561 Y=YA(NS)
1562 NS=NS-1
1563 DO 13 M=1,N-1
1564 DO 12 I=1,N-M
1565 W=C(I+1)-D(I)
1566 H=XA(I+M)-X
1567 T=(XA(I)-X)*D(I)/H
1568 DD=T-C(I+1)
1569 DD=W/DD
1570 D(I)=C(I+1)*DD
1571 C(I)=T*DD
157212 CONTINUE
1573 IF (2*NS.LT.N-M)THEN
1574 DY=C(NS+1)
1575 ELSE
1576 DY=D(NS)
1577 NS=NS-1
1578 ENDIF
1579 Y=Y+DY
158013 CONTINUE
1581 RETURN
1582 END
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)
1589 TEM1 = PI / ALFST
1590 TEM2 = 1. / CtLhALPQCD (JORD, NEFF, EFMULM, I)
1591 CtLhRTALF = TEM1 - TEM2
1592 END
1593 Subroutine CtLhbldat1
1594 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1595 include 'parmsetup.inc'
1596 LOGICAL LSTX
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)
1605 XMIN = .999999D-4
1606 XCR = 1.5
1607 JT = 1
1608 Return
1609 END
1610 SUBROUTINE CtLhSETL1 (NEF, VLAM)
1611 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1612 LOGICAL SET
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
1618 STOP
1619 ENDIF
1620 AMHAT(0) = 0.
1621 DO 5 N = 1, NF
1622 AMHAT(N) = VALQMS(N)
1623 5 CONTINUE
1624 ALAM(NEF) = VLAM
1625 DO 10 N = NEF, 1, -1
1626 CALL CtLhTRNLAM(NORDER, N, -1, IR1)
1627 10 CONTINUE
1628 DO 20 N = NEF, NF-1
1629 CALL CtLhTRNLAM(NORDER, N, 1, IR1)
1630 20 CONTINUE
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
1634 NHQ = NF - N
1635 GOTO 40
1636 ENDIF
1637 30 CONTINUE
1638 NHQ = NF
1639 40 CONTINUE
1640 DO 50, N = NF-NHQ, 1, -1
1641 AMHAT(N) = 0
1642 ALAM(N-1) = ALAM(N)
1643 50 CONTINUE
1644 AMN = ALAM(NF)
1645 DO 60, N = 0, NF-1
1646 IF (ALAM(N) .GT. AMN) AMN = ALAM(N)
1647 60 CONTINUE
1648 AMN = AMN * 1.0001
1649 AL = ALAM(NF)
1650 SET = .TRUE.
1651 RETURN
1652 END
1653 SUBROUTINE CtLhSETLAM (NEF, WLAM, IRDR)
1654 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1655 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1656 LOGICAL SET
1657 IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN
1658 WRITE(*,*)'NEF out of range in CtLhSETLAM: NEF NF=',NEF,NF
1659 STOP
1660 ENDIF
1661 VLAM = WLAM
1662 IF (IRDR .NE. NORDER) then
1663 PRINT *,'fatal error: wanted cnvl1'
1664 stop
1665 ENDIF
1666 CALL CtLhSETL1 (NEF, VLAM)
1667 END
1668 Subroutine CtLhbldat2
1669 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1670 COMMON / LhCtCOMQMS / VALQMS(9)
1671 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1672 LOGICAL SET
1673 AL = .226d0
1674 NF = 5
1675 NORDER = 2
1676 SET = .FALSE.
1677 VALQMS(1) = 0.
1678 VALQMS(2) = 0.
1679 VALQMS(3) = 0.2d0
1680 VALQMS(4) = 1.3d0
1681 VALQMS(5) = 4.5d0
1682 VALQMS(6) = 174.d0
1683 VALQMS(7) = 0.
1684 VALQMS(8) = 0.
1685 VALQMS(9) = 0.
1686 Return
1687 END
1688 FUNCTION CtLhSMPNOL (NX, DX, FN, ERR)
1689 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1690 DIMENSION FN(NX)
1691 MS = MOD(NX, 2)
1692 IF (NX .LE. 1 .OR. NX .GT. 1000) THEN
1693 PRINT *, 'NX =', NX, ' OUT OF RANGE IN CtLhSMPNOL!'
1694 STOP
1695 ELSEIF (NX .EQ. 2) THEN
1696 TEM = DX * FN(2)
1697 ELSEIF (NX .EQ. 3) THEN
1698 TEM = DX * FN(2) * 2.
1699 ELSE
1700 IF (MS .EQ. 0) THEN
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
1706 ELSE
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
1712 ENDIF
1713 ENDIF
1714 CtLhSMPNOL = TEM
1715 RETURN
1716 END
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)
1721 DIMENSION F(NX)
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',
1725 > DX, D0, D1, 0)
1726 CtLhSMPSNA = 0.
1727 RETURN
1728 ENDIF
1729 IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN
1730 CALL CtLhWARNI(IW1, 'NX out of range in CtLhSMPSNA', 'NX', NX,
1731 > 1, MAXX, 1)
1732 SIMP = 0.
1733 ELSEIF (NX .EQ. 1) THEN
1734 SIMP = 0.
1735 ELSEIF (NX .EQ. 2) THEN
1736 SIMP = (F(1) + F(2)) / 2.
1737 ERRD = (F(1) - F(2)) / 2.
1738 ELSE
1739 MS = MOD(NX, 2)
1740 IF (MS .EQ. 0) THEN
1741 ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24.
1742 NZ = NX - 1
1743 ELSE
1744 ADD = 0.
1745 NZ = NX
1746 ENDIF
1747 IF (NZ .EQ. 3) THEN
1748 SIMP = (F(1) + 4.* F(2) + F(3)) / 3.
1749 TRPZ = (F(1) + 2.* F(2) + F(3)) / 2.
1750 ELSE
1751 SE = F(2)
1752 SO = 0
1753 NM1 = NZ - 1
1754 DO 60 I = 4, NM1, 2
1755 IM1 = I - 1
1756 SE = SE + F(I)
1757 SO = SO + F(IM1)
1758 60 CONTINUE
1759 SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3.
1760 TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2.
1761 ENDIF
1762 ERRD = TRPZ - SIMP
1763 SIMP = SIMP + ADD
1764 ENDIF
1765 CtLhSMPSNA = SIMP * DX
1766 IF (ABS(SIMP) .GT. TINY) THEN
1767 ERR = ERRD / SIMP
1768 ELSE
1769 ERR = 0.
1770 ENDIF
1771 RETURN
1772 END
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)
1783 DATA D0 / 0.0 /
1784 JTT = 2 * JT
1785 DDT = DT / JTT
1786 IF (NX .GT. MXX) THEN
1787 WRITE (*,*) 'Nx =', NX, ' too many pts in CtLhSNEVL'
1788 STOP 'Program stopped in CtLhSNEVL'
1789 EndIf
1790c ++ remove unused quantities (jcp)
1791c ++ TMD = TIN + DT * NT / 2.
1792c ++ AMU = EXP(TMD)
1793c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
1794c ++ TLAM = TMD - TEM
1795 DO 9 IX = 1, NX
1796 US (IX, 0) = UI(IX)
1797 GS (IX, 0) = GI(IX)
1798 9 CONTINUE
1799 US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3)
1800 GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3)
1801 TT = TIN
1802 DO 10 IZ = 1, NX
1803 Y0(IZ) = UI(IZ)
1804 Z0(IZ) = GI(IZ)
1805 10 CONTINUE
1806 DO 20 IS = 1, NT
1807 DO 202 JS = 1, JTT
1808 IRND = (IS-1) * JTT + JS
1809 IF (IRND .EQ. 1) THEN
1810 CALL CtLhSNRHS (TT, NEFF, Y0,Z0, F0,G0)
1811 DO 250 IZ = 1, NX
1812 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
1813 Z0(IZ) = Z0(IZ) + DDT * G0(IZ)
1814 250 CONTINUE
1815 TT = TT + DDT
1816 CALL CtLhSNRHS (TT, NEFF, Y0, Z0, F1, G1)
1817 DO 251 IZ = 1, NX
1818 Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
1819 Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0
1820 251 CONTINUE
1821 Else
1822 CALL CtLhSNRHS (TT, NEFF, Y1, Z1, F1, G1)
1823 DO 252 IZ = 1, NX
1824 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
1825 ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0
1826 252 CONTINUE
1827 TT = TT + DDT
1828 CALL CtLhSNRHS (TT, NEFF, YP, ZP, FP, GP)
1829 DO 253 IZ = 1, NX
1830 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
1831 Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0
1832 F0(IZ) = F1(IZ)
1833 G0(IZ) = G1(IZ)
1834 253 CONTINUE
1835 EndIf
1836 202 CONTINUE
1837 DO 260 IX = 1, NX
1838 IF (IKNL .GT. 0) THEN
1839 US (IX, IS) = MAX(Y1(IX), D0)
1840 GS (IX, IS) = MAX(Z1(IX), D0)
1841 Else
1842 US (IX, IS) = Y1(IX)
1843 GS (IX, IS) = Z1(IX)
1844 EndIf
1845 260 CONTINUE
1846 US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
1847 GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3)
1848 20 CONTINUE
1849 RETURN
1850 END
1851 SUBROUTINE CtLhSNRHS (TT, NEFF, FI, GI, FO, GO)
1852 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1853 LOGICAL LSTX
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)
1866 S = EXP(TT)
1867 Q = AL * EXP (S)
1868 CPL = CtLhALPI(Q)
1869 CPL2= CPL ** 2 / 2. * S
1870 CPL = CPL * 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
1881 DO 230 IZ = 1, NX
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
1890 230 CONTINUE
1891 Else
1892 DO 240 IZ = 1, NX
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
1899 240 CONTINUE
1900 EndIf
1901 IF (IKNL .EQ. 2) THEN
1902 DZ = 1./(NX - 1)
1903 DO 21 I = 1, NX-1
1904 NP = NX - I + 1
1905 IS = NP
1906 g2(1)=0d0
1907 g3(1)=0d0
1908 DO 31 KZ = 2, NP
1909 IY = I + KZ - 1
1910 IT = NX - IY + 1
1911 XY = ZZ (IS, IT)
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)
1916 31 CONTINUE
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)
1925 TMF = TEM1 + TEM2
1926 TMG = TEM3 + TEM4
1927 FO(I) = FO(I) + TMF * CPL2
1928 GO(I) = GO(I) + TMG * CPL2
1929 21 CONTINUE
1930 EndIf
1931 RETURN
1932 END
1933 FUNCTION CtLhSPENC2 (X)
1934 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1935 EXTERNAL CtLhSPN2IN
1936 COMMON / LhCtSPENCC / XX
1937 DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 /
1938 XX = X
1939 TEM = CtLhGausInt(CtLhSPN2IN, XX, U1, AERR, RERR, ERR, IRT)
1940 CtLhSPENC2 = TEM + LOG (XX) ** 2 / 2.
1941 RETURN
1942 END
1943 FUNCTION CtLhSPN2IN (ZZ)
1944 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1945 COMMON / LhCtSPENCC / X
1946 Z = ZZ
1947 TEM = LOG (1.+ X - Z) / Z
1948 CtLhSPN2IN = TEM
1949 RETURN
1950 END
1951 SUBROUTINE CtLhSTUPKL (NFL)
1952 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1953 LOGICAL LSTX
1954 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1955 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1956 PARAMETER (MX = 3)
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)
1974 SAVE
1975 DATA AERR, RERR / 0.0, 0.02 /
1976 NNX = NX
1977 DZ = 1./ (NX - 1)
1978 DO 5 I0 = 1, MX
1979 XL(I0) = XV(I0)
1980 5 CONTINUE
1981 DO 10 I = 1, NX-1
1982 XZ = XV(I)
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)
1985 10 CONTINUE
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)
2006 RER = RERR * 4.
2007 AFF1(1) = CtLhGausInt(CtLhPFF1,D0,XV(1),AERR,RERR,ER1,IRT)
2008 DGG1 = NFL / 3.
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)
2013 AER = AFF1(1) * RER
2014 AFF2(1) = CtLhGausInt(CtLhRFF2, D0, XV(1), AER, RER, ER2, IRT)
2015 AER = AGG1(1) * RER
2016 AGG2(1) = CtLhGausInt(CtLhRGG2, D0, XV(1), AER, RER, ER4, IRT)
2017 DO 20 I2 = 2, NX-1
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)
2022 > +AFF2(I2-1)
2023 TEM = CtLhGausInt(CtLhRGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT)
2024 TMPG = TMPG + TEM
2025 AGG1(I2) = TMPG + DGG1
2026 AER = ABS(TEM * RER)
2027 AGG2(I2)=CtLhGausInt(CtLhRGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)
2028 > +AGG2(I2-1)
2029 ANSP(I2)=CtLhGausInt(CtLhFNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2030 > +ANSP(I2-1)
2031 ANSM(I2)=CtLhGausInt(CtLhFNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2032 > +ANSM(I2-1)
2033 20 CONTINUE
2034 ANSP(NX)=CtLhGausInt(CtLhFNSP,XV(NX-1),D1,AERR,RER,ERR,
2035 > IRT) + ANSP(NX-1)
2036 ANSM(NX)=CtLhGausInt(CtLhFNSM,XV(NX-1),D1,AERR,RER,ERR,
2037 > IRT) + ANSM(NX-1)
2038 TRNF = TR * NFL
2039 do i2=1,nx-1 !loop over x
2040 x=xv(i2)
2041c XI = 1./ X !unused - jcp
2042 X2 = X ** 2
2043 X3= x**3
2044 XLN = DLOG (X)
2045 XLN2 = XLN ** 2
2046 XLN1M = DLOG (1.- X)
2047 xLi2m=CtLhxLi(2,-x)
2048 xLi2=CtLhxLi(2,x)
2049 xLi3=CtLhxLi(3,x)
2050 xLi31m=CtLhxLi(3,1d0-x)
2051 xLi32=CtLhxLi(3,x2)
2052 xln1m2=xln1m*xln1m
2053 xln1p=dlog(1d0+x)
2054 x1m=1d0-x
2055 x1p=1d0+x
2056 x3m=3d0-x
2057 x3p=3d0+x
2058 wgfcft=
2059 > (9 + 4*Pi2 - 22*x + 13*x2 + 6*(3 - 4*x + x2)*xln1m +
2060 > 40*xln - 24*xLi2)/9.
2061 wgfcf2=
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
2065 wgfcfg=
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
2072 wfgcft=
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.
2078 wfgcgt=
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
2085 enddo !i2
2086 AGF2(nx)=0d0
2087 AFG2(nx)=0d0
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)
2091 DO 21 IX = 1, NX-1
2092 X = XV(IX)
2093 NP = NX - IX + 1
2094 IS = NP
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)
2099 DO 31 KZ = 2, NP
2100 IY = IX + KZ - 1
2101 IT = NX - IY + 1
2102 XY = X / XV(IY)
2103 XM1 = 1.- XY
2104 XG2 = (LOG(1./XM1) + 1.) ** 2
2105 Z = ZZ (IX, IY)
2106 TZ = (Z + DZ) / DZ
2107 IZ = Int(TZ)
2108 IZ = MAX (IZ, 0)
2109 IZ = MIN (IZ, NX-1)
2110 DT = TZ - IZ
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)
2123 31 CONTINUE
2124 21 CONTINUE
2125 RETURN
2126 END
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
2131 EXTERNAL CtLhZBRLAM
2132 DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 /
2133 DATA IR1, SML / 0, 1.E-5 /
2134 IRT = 0
2135 N = NF
2136 JRDR = IRDR
2137 JACT = IACT
2138 VLAM = ALAM(N)
2139 IF (JACT .GT. 0) THEN
2140 N1 = N + 1
2141 THMS = AMHAT(N1)
2142 ALM = LOG (THMS/VLAM)
2143 BLM = BLM0
2144 ELSE
2145 N1 = N -1
2146 THMS = AMHAT(N)
2147 ALM = ALM0
2148 THMS = MAX (THMS, SML)
2149 BLM = LOG (THMS/VLAM)
2150 ENDIF
2151 IF (VLAM .GE. 0.7 * THMS) THEN
2152 IF (JACT . EQ. 1) THEN
2153 AMHAT(N1) = 0
2154 ELSE
2155 AMHAT(N) = 0
2156 ENDIF
2157 IRT = 4
2158 ALAM(N1) = VLAM
2159 RETURN
2160 ENDIF
2161 IF (ALM .GE. BLM) THEN
2162 WRITE (*, *) 'CtLhTRNLAM has ALM >= BLM: ', ALM, BLM
2163 WRITE (*, *) 'I do not know how to continue'
2164 STOP
2165 ENDIF
2166 VMULM = THMS/VLAM
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'
2174 STOP
2175 ENDIF
2176 RETURN
2177 END
2178 SUBROUTINE CtLhUPC (A, La, UpA)
2179 CHARACTER A*(*), UpA*(*), C*(1)
2180 INTEGER I, La, Ld
2181 La = Len(A)
2182 Lb = Len(UpA)
2183 If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!'
2184 Ld = ICHAR('A')-ICHAR('a')
2185 DO 1 I = 1, Lb
2186 If (I .Le. La) Then
2187 c = A(I:I)
2188 IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN
2189 UpA (I:I) = CHAR(Ichar(c) + ld)
2190 Else
2191 UpA (I:I) = C
2192 ENDIF
2193 Else
2194 UpA (I:I) = ' '
2195 Endif
2196 1 CONTINUE
2197
2198 RETURN
2199 END
2200 SUBROUTINE CtLhWARNI (IWRN, MSG, NMVAR, IVAB,
2201 > IMIN, IMAX, IACT)
2202 CHARACTER*(*) MSG, NMVAR
2203 Save Iw
2204 Data Nmax / 100 /
2205 IW = IWRN
2206 IV = IVAB
2207
2208 IF (IW .EQ. 0) THEN
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
2212 ENDIF
2213 ENDIF
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'
2218 Endif
2219 IWRN = IW + 1
2220 RETURN
2221 END
2222 SUBROUTINE CtLhWARNR (IWRN, MSG, NMVAR, VARIAB,
2223 > VMIN, VMAX, IACT)
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
2227 Save Iw
2228 Data Nmax / 100 /
2229 IW = IWRN
2230 VR = VARIAB
2231 IF (IW .EQ. 0) THEN
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
2235 ENDIF
2236 ENDIF
2237 If (Iw .LT. Nmax) Then
2238 PRINT '(I5, 2A/1X,2A,1PD16.7)', IW, ' ', MSG,
2239 > NMVAR, ' = ', VR
2240 Elseif (Iw .Eq. Nmax) Then
2241 Print '(/A/)', 'CtLhWARNR Severe Warning: Too many errors'
2242 Endif
2243 IWRN = IW + 1
2244 RETURN
2245 END
2246 SUBROUTINE CtLhXARRAY
2247 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2248 LOGICAL LSTX
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)
2254 Character Msg*80
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 /
2263 DATA PUNY / 1D-30 /
2264 XV(0) = 0D0
2265 DZ = 1D0 / (NX-1)
2266 DO 10 I = 1, NX - 1
2267 Z = DZ * (I-1)
2268 ZV(I) = Z
2269 X = CtLhXFRMZ (Z)
2270 DXTZ(I) = CtLhDXDZ(Z) / X
2271 XV (I) = X
2272 XA(I, 1) = X
2273 XA(I, 0) = LOG (X)
2274 DO 20 L = L1, L2
2275 IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L
2276 20 CONTINUE
2277 10 CONTINUE
2278 XV(1) = Xmin
2279 XV(NX) = 1D0
2280 ZV(Nx) = 1D0
2281 DXTZ(NX) = CtLhDXDZ(1.D0)
2282 DO 21 L = L1, L2
2283 XA (NX, L) = 1D0
2284 21 CONTINUE
2285 XA (NX, 0) = 0D0
2286 DO 11 I = 1, NX-1
2287 ELY(I) = LOG(1D0 - XV(I))
2288 11 CONTINUE
2289 ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3)
2290 DO 17 IX = 1, NX
2291 ZZ (IX, IX) = 1.
2292 DO 17 IY = IX+1, NX
2293 XY = XV(IX) / XV(IY)
2294 ZZ (IX, IY) = CtLhZFRMX (XY)
2295 ZZ (NX-IX+1, NX-IY+1) = XY
2296 17 CONTINUE
2297 DO 30 I = 1, NX-1
2298 IF (I .NE. NX-1) THEN
2299 F11 = 1D0/XV(I)
2300 F21 = 1D0/XV(I+1)
2301 F31 = 1D0/XV(I+2)
2302 F13 = XV(I)
2303 F23 = XV(I+1)
2304 F33 = XV(I+2)
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)
2310 DET = PUNY
2311 EndIf
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)
2326 EndIf
2327 DO 51 J = 1, NDH
2328 DO 52 L = 1, NDG
2329 IF (I .EQ. 1) THEN
2330 GB(L,J,I) = G2(L,J)
2331 ElseIF (I .EQ. NX-1) THEN
2332 GB(L,J,I) = G1(L,J)
2333 Else
2334 GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0
2335 EndIf
2336 52 CONTINUE
2337 51 CONTINUE
2338 DO 35 MM = M1, M2
2339 DO 40 K = 1, NDG
2340 KK = K + MM - 2
2341 IF (KK .EQ. 0) THEN
2342 A(K) = XA(I+1, 0) - XA(I, 0)
2343 Else
2344 A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK)
2345 EndIf
2346 40 CONTINUE
2347 DO 41 J = 1, NDH
2348 TEM = 0
2349 DO 43 L = 1, NDG
2350 TEM = TEM + A(L) * GB(L,J,I)
2351 43 CONTINUE
2352 H(J,I,MM) = TEM
2353 41 CONTINUE
2354 35 CONTINUE
2355 DO 42 J = 1, NDG
2356 DO 44 L = 1, NDG
2357 G1(L,J) = G2(L,J+1)
2358 44 CONTINUE
2359 42 CONTINUE
2360 30 CONTINUE
2361 LSTX = .TRUE.
2362 RETURN
2363 END
2364 FUNCTION CtLhXFRMZ (Z)
2365 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2366 LOGICAL LSTX
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
2371 EXTERNAL CtLhZFXL
2372 DATA TEM, RER / D1, 1E-3 /
2373 DATA ZLOW, ZHIGH, IWRN2 / -10.0, 1.00002, 0 /
2374 EPS = TEM * RER
2375 ZA = Z
2376 IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN
2377 XLA = LOG (XMIN) * 1.5
2378 XLB = 0.00001
2379 TEM = CtLhZBRNT (CtLhZFXL, XLA, XLB, EPS, IRT)
2380 Else
2381 CALL CtLhWARNR (IWRN2, 'Z out of range in CtLhXFRMZ, X set=0.',
2382 > 'Z', Z, ZLOW, ZHIGH, 1)
2383 TEM = 0
2384 EndIf
2385 CtLhXFRMZ = EXP(TEM)
2386 RETURN
2387 END
2388 FUNCTION CtLhxLi(n,x)
2389 implicit NONE
2390 integer NCUT, i,n,m3
2391 real*8 CtLhxLi,Out,x,pi2by6,zeta3,c1,c2
2392 real*8 r,xt,L,xln1m
2393 parameter (m3=8)
2394 dimension c1(2:m3),c2(2:m3)
2395 data NCUT/27/
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/
2402 L=0.0
2403 i=0
2404 r=1.0
2405 if (abs(x).gt.r) then
2406 PRINT *,'Li: x out of range (-1,1) , x=',x
2407 STOP
2408 endif
2409 if (n.lt.0) then
2410 PRINT *,'Polylogarithm Li undefined for n=',n
2411 STOP
2412 elseif (n.eq.0) then
2413 Out=x/(1d0-x)
2414 elseif (n.eq.1) then
2415 Out=-dlog(1-x)
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
2420 do while(i.le.NCUT)
2421 i=i+1
2422 r=r*x
2423 L=L+r/i/i
2424 enddo
2425 Out=L
2426 elseif (x.eq.0) then
2427 Out=0d0
2428 elseif(x.gt.0.5) then !n.eq.2,x>0.5
2429 xt = 1.0-x
2430 L = pi2by6 - dlog(x)*dlog(xt)
2431 do while(i.le.NCUT)
2432 i=i+1
2433 r=r*xt
2434 L=L-r/i/i
2435 enddo
2436 Out=L
2437 elseif (x.lt.(-0.5)) then
2438 xt=-x/(1d0-x)
2439 L=-0.5*dlog(1-x)**2
2440 do while (i.le.NCUT)
2441 i=i+1
2442 r=r*xt
2443 L=L-r/i/i
2444 enddo
2445 Out=L
2446 endif
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)
2449 xt=(1d0-x)
2450 xln1m=dlog(xt)
2451 do i=2,m3
2452 L=L+(c1(i)+c2(i)*xln1m)*xt**i
2453 enddo
2454 Out=L
2455 else !n>3 or x=3,x<0.8
2456 do while(i.le.NCUT)
2457 i=i+1
2458 r=r*x
2459 L=L+r/dble(i)**dble(n)
2460 enddo
2461 Out=L
2462 endif
2463 CtLhxLi=Out
2464 End ! CtLhxLi
2465 FUNCTION CtLhZBRLAM (WLLN)
2466 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2467 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2468 WMULM = EXP (WLLN)
2469 TEM1 = 1./ CtLhALPQCD(JRDR, N1, WMULM, I)
2470 TEM2 = 1./ CtLhALPQCD(JRDR, N, VMULM, I)
2471 CtLhZBRLAM = TEM1 - TEM2
2472 END
2473 FUNCTION CtLhZBRNT(FUNC, X1, X2, TOL, IRT)
2474 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2475 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
2476 external func
2477 IRT = 0
2478 TOL = ABS(TOL)
2479 A=X1
2480 B=X2
2481 FA=FUNC(A)
2482 FB=FUNC(B)
2483 IF(FB*FA.GT.0.) THEN
2484 PRINT *, 'Root must be bracketed for CtLhZBRNT. Set = 0'
2485 IRT = 1
2486 CtLhZBRNT=0.
2487 RETURN
2488 ENDIF
2489 FC=FB
2490 DO 11 ITER=1,ITMAX
2491 IF(FB*FC.GT.0.) THEN
2492 C=A
2493 FC=FA
2494 D=B-A
2495 E=D
2496 ENDIF
2497 IF(ABS(FC).LT.ABS(FB)) THEN
2498 A=B
2499 B=C
2500 C=A
2501 FA=FB
2502 FB=FC
2503 FC=FA
2504 ENDIF
2505 TOL1=2.*EPS*ABS(B)+0.5*TOL
2506 XM=.5*(C-B)
2507 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
2508 CtLhZBRNT=B
2509 RETURN
2510 ENDIF
2511 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
2512 S=FB/FA
2513 IF(A.EQ.C) THEN
2514 P=2.*XM*S
2515 Q=1.-S
2516 ELSE
2517 Q=FA/FC
2518 R=FB/FC
2519 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
2520 Q=(Q-1.)*(R-1.)*(S-1.)
2521 ENDIF
2522 IF(P.GT.0.) Q=-Q
2523 P=ABS(P)
2524 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
2525 E=D
2526 D=P/Q
2527 ELSE
2528 D=XM
2529 E=D
2530 ENDIF
2531 ELSE
2532 D=XM
2533 E=D
2534 ENDIF
2535 A=B
2536 FA=FB
2537 IF(ABS(D) .GT. TOL1) THEN
2538 B=B+D
2539 ELSE
2540 B=B+SIGN(TOL1,XM)
2541 ENDIF
2542 FB=FUNC(B)
254311 CONTINUE
2544 PRINT *, 'CtLhZBRNT exceeding maximum iterations.'
2545 IRT = 2
2546 CtLhZBRNT=B
2547 RETURN
2548 END
2549 FUNCTION CtLhZFRMX (XX)
2550 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2551 LOGICAL LSTX
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)
2558 X = XX
2559 IF (X .GE. XMIN) THEN
2560 TEM = F(X) / F(D1)
2561 ElseIF (X .GE. D0) THEN
2562 X = MAX (X, TINY)
2563 TEM = F(X) / F(D1)
2564 Else
2565 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhZFRMX'
2566 > , 'X', X, TINY, HUGE, 1)
2567 TEM = 99.
2568 STOP
2569 EndIf
2570 CtLhZFRMX = TEM
2571 RETURN
2572 ENTRY CtLhDZDX (XX)
2573 X = XX
2574 IF (X .GE. XMIN) THEN
2575 TEM = D(X) / F(D1)
2576 ElseIF (X .GE. D0) THEN
2577 X = MAX (X, TINY)
2578 TEM = D(X) / F(D1)
2579 Else
2580 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhDZDX '
2581 > , 'X', X, TINY, HUGE, 1)
2582 TEM = 99.
2583 STOP
2584 EndIf
2585 CtLhDZDX = TEM
2586 RETURN
2587 END
2588 FUNCTION CtLhZFXL (XL)
2589 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2590 COMMON / LhCtINVERT / ZA
2591 X = EXP(XL)
2592 TT = CtLhZFRMX (X) - ZA
2593 CtLhZFXL = TT
2594 RETURN
2595 END
2596