]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.2.2/EVLCTEQ.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / EVLCTEQ.f
CommitLineData
3c5d1739 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 QST = QS
7 ALFST = ALFS
8 CALL CtLhParQcd (2, 'ORDR', ORDR, IR1)
9 JORD = ORDR
10 NEFF = LhCtNFL(QS)
11 EFLLN = CtLhQZBRNT (CtLhRTALF, ALAM, BLAM, ERR, IR2)
12 EFFLAM = QS / EXP (EFLLN)
13 CALL CtLhSETL1 (NEFF, EFFLAM)
14 END
15 FUNCTION CtLhALPI (AMU)
16 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
17 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
18 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
19 LOGICAL SET
20 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
21 DATA IW1, IW2 / 2*0 /
22 IF(.NOT.SET) CALL CtLhLAMCWZ
23 NEFF = LhCtNFL(AMU)
24 ALM = ALAM(NEFF)
25 CtLhALPI = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT)
26 IF (IRT .EQ. 1) THEN
27 CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI', 'AMU', AMU,
28 > ALM, BIG, 1)
29 ELSEIF (IRT .EQ. 2) THEN
30 CALL CtLhWARNR (IW2, 'CtLhALPI > 3; Be aware!', 'CtLhALPI',
31 > CtLhALPI, D0, D1, 0)
32 ENDIF
33 RETURN
34 END
35 FUNCTION CtLhALPQCD (IRDR, NF, RML, IRT)
36 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
37 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
38 PARAMETER (CG = 3.0d0, TR = 0.5d0, CF = 4.0d0/3.0d0)
39 IRT = 0
40 IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN
41 print *,
42 > 'Order out of range in CtLhALPQCD: IRDR = ', IRDR
43 STOP
44 ENDIF
45 B0 = (11.d0*CG - 2.* NF) / 3.d0
46 B1 = (34.d0*CG**2 - 10.d0*CG*NF - 6.d0*CF*NF) / 3.d0
47 RM2 = RML**2
48 IF (RM2 .LE. 1.) THEN
49 IRT = 1
50 CtLhALPQCD = 99.
51 RETURN
52 ENDIF
53 ALN = LOG (RM2)
54 AL = 4.d0/ B0 / ALN
55 IF (IRDR .GE. 2) AL = AL * (1.d0-B1*LOG(ALN) / ALN / B0**2)
56 IF (AL .GE. 3.) THEN
57 IRT = 2
58 ENDIF
59 CtLhALPQCD = AL
60 RETURN
61 END
62 FUNCTION CtLhAMHATF(I)
63 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
64 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
65 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
66 LOGICAL SET
67 IF (.NOT.SET) CALL CtLhLAMCWZ
68 IF ((I.LE.0).OR.(I.GT.9)) THEN
69 print *,'warning I OUT OF RANGE IN CtLhAMHATF'
70 CtLhAMHATF = 0
71 ELSE
72 CtLhAMHATF = AMHAT(I)
73 ENDIF
74 RETURN
75 END
76 FUNCTION CtLhDXDZ (Z)
77 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
78 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
79 DATA HUGE, IWRN / 1.E20, 0 /
80 ZZ = Z
81 X = CtLhXFRMZ (ZZ)
82 TEM = CtLhDZDX (X)
83 IF (TEM .NE. D0) THEN
84 TMP = D1 / TEM
85 Else
86 CALL CtLhWARNR(IWRN, 'CtLhDXDZ singular in CtLhDXDZ; set=HUGE',
87 > 'Z', Z, D0, D1, 0)
88 TMP = HUGE
89 EndIf
90 CtLhDXDZ = TMP
91 RETURN
92 END
93 SUBROUTINE CtLhEVLPAR (IACT, NAME, VALUE, IRET)
94 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
95 CHARACTER*(*) NAME
96 IRET = 1
97 IF (IACT .EQ. 0) THEN
98 WRITE ( NINT(VALUE) , 101)
99 101 FORMAT (/ ' Initiation parameters: Qini, Ipd0, Ihdn ' /
100 > ' Maximum Q, Order of Alpha: Qmax, IKNL ' /
101 > ' X- mesh parameters : Xmin, Xcr, Nx ' /
102 > ' LnQ-mesh parameters : Nt, Jt ' /
103 > ' # of parton flavors : NfMx ' /)
104 IRET = 4
105 ElseIF (IACT .EQ. 1) THEN
106 CALL CtLhEVLSET (NAME, VALUE, IRET)
107 Else
108 print *,'fatal evlpar'
109 stop
110 EndIf
111 RETURN
112 END
113 SUBROUTINE CtLhEVLSET (NAME, VALUE, IRET)
114 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
115 LOGICAL LSTX
116 CHARACTER*(*) NAME
117 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
118 PARAMETER (MXPN = MXF * 2 + 2)
119 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
120 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
121 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
122 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
123 > / PdfSwh / Iset, IpdMod, Iptn0, NuIni
124 IRET = 1
125 IF (NAME .EQ. 'QINI') THEN
126 IF (VALUE .LE. 0) GOTO 12
127 QINI = VALUE
128 ElseIF (NAME .EQ. 'IPD0') THEN
129 ITEM = NINT(VALUE)
130 IF (Item .Eq. 10 .or. Item .Eq. 11) GOTO 12
131 IPD0 = ITEM
132 ElseIF (NAME .EQ. 'IHDN') THEN
133 ITEM = NINT(VALUE)
134 IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12
135 IHDN = ITEM
136 ElseIF (NAME .EQ. 'QMAX') THEN
137 IF (VALUE .LE. QINI) GOTO 12
138 QMAX = VALUE
139 ElseIF (NAME .EQ. 'IKNL') THEN
140 ITMP = NINT(VALUE)
141 ITEM = ABS(ITMP)
142 IF (ITEM.NE.1.AND.ITEM.NE.2) GOTO 12
143 IKNL = ITMP
144 ElseIF (NAME .EQ. 'XCR') THEN
145 IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12
146 XCR = VALUE
147 LSTX = .FALSE.
148 ElseIF (NAME .EQ. 'XMIN') THEN
149 IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12
150 XMIN = VALUE
151 LSTX = .FALSE.
152 ElseIF (NAME .EQ. 'NX') THEN
153 ITEM = NINT(VALUE)
154 IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12
155 NX = ITEM
156 LSTX = .FALSE.
157 ElseIF (NAME .EQ. 'NT') THEN
158 ITEM = NINT(VALUE)
159 IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12
160 NT = ITEM
161 ElseIF (NAME .EQ. 'JT') THEN
162 ITEM = NINT(VALUE)
163 IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12
164 JT = ITEM
165 ElseIF (NAME .EQ. 'NFMX') THEN
166 ITEM = NINT(VALUE)
167 IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12
168 NfMx = ITEM
169 ElseIF (NAME .EQ. 'IPDMOD') THEN
170 ITEM = NINT(VALUE)
171 IF (Abs(Item) .Gt. 1) GOTO 12
172 IpdMod = ITEM
173 ElseIF (NAME .EQ. 'IPTN0') THEN
174 ITEM = NINT(VALUE)
175 IF (ABS(ITEM) .GT. MXF) GOTO 12
176 IPTN0 = ITEM
177 ElseIF (NAME .EQ. 'NUINI') THEN
178 ITEM = NINT(VALUE)
179 IF (ITEM .LE. 0) GOTO 12
180 NuIni = ITEM
181 Else
182 IRET = 0
183 EndIf
184 RETURN
185 12 IRET = 2
186 RETURN
187 END
188 SUBROUTINE CtLhEVOLVE (FINI, IRET)
189 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
190 include 'parmsetup.inc'
191 LOGICAL LSTX
192 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
193 PARAMETER (MXPN = MXF * 2 + 2)
194 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
195 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
196 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
197 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
198 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
199 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
200 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
201 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
202 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
203 DIMENSION QRKP(MXF)
204 DIMENSION JI(-MXF : MXF+1)
205 EXTERNAL LhCtNSRHSP, LhCtNSRHSM, FINI
206 DATA DZER / 0.0 /
207 save nxsave, ntsave, jtsave, ngsave,
208 & xcrsave, xminsave, qinisave, qmaxsave, ientry, ishow
209 data ientry / 0 /
210 dimension xvsave(0:MXX), qvsave(0:MXQ), tvsave(0:MXQ)
211c
212 call getnset(iset)
213c
214 ientry = ientry + 1
215 if(ientry .lt. 100) then
216666 format(1x,'enter evolve',i3)
217 elseif(ientry .eq. 100) then
218667 format(1x,'enter evolve',i3,' further suppressed')
219 endif
220 ishow = 0 !set to no display
221 if(ientry .eq. 1) then
222 ishow = 1 !turn display on
223 endif
224 11 IRET = 0
225 IF (IHDN .LE. 4) THEN
226 MXVAL = 2
227 ElseIF (IHDN .LE. 6) THEN
228 MXVAL = 3
229 EndIf
230 IF (.NOT. LSTX) CALL CtLhXARRAY
231 DLX = 1.D0 / NX
232 J10 = NX / 10
233 CALL CtLhPARPDF (2, 'ALAM', AL, IR)
234 CALL CtLhQARRAY (NINI)
235 NFSN = NFMX + 1
236 KF = 2 * NFMX + 2
237 Nelmt = KF * (Nt+1) * (Nx+1)
238 DO 101 IFLV = -NFMX, NFMX+1
239 JFL = NFMX + IFLV
240 JI(IFLV) = JFL * (NT+1) * (NX+1)
241 101 CONTINUE
242 3 DO 31 IZ = 1, NX
243 UPD(JI(0)+IZ+1,iset) = FINI (0, XV(IZ))
244 UPD(JI(NFSN)+IZ+1,iset) = 0
245 IF (NFMX .EQ. 0) GOTO 31
246 DO 331 IFLV = 1, NINI
247 A = FINI ( IFLV, XV(IZ))
248 B = FINI (-IFLV, XV(IZ))
249 QRKP (IFLV) = A + B
250 UPD(JI(NFSN)+IZ+1,iset) =
251 > UPD(JI(NFSN)+IZ+1,iset) + QRKP (IFLV)
252 UPD(JI(-IFLV)+IZ+1,iset) = A - B
253 331 CONTINUE
254 DO 332 IFLV = 1, NINI
255 UPD(JI( IFLV)+IZ+1,iset) =
256 > QRKP(IFLV) - UPD(JI(NFSN)+IZ+1,iset)/NINI
257 332 CONTINUE
258 31 CONTINUE
259 DO 21 NEFF = NINI, NFMX
260 IF (IKNL .EQ. 2) CALL CtLhSTUPKL (NEFF)
261 ICNT = NEFF - NINI + 1
262 IF (NTN(ICNT) .EQ. 0) GOTO 21
263 NITR = NTN (ICNT)
264 DT = DTN (ICNT)
265 TIN = TLN (ICNT)
266 CALL CtLhSNEVL (IKNL, NX, NITR, JT, DT, TIN, NEFF,
267 > UPD(JI(NFSN)+2,iset), UPD(JI(0)+2,iset),
268 > UPD(JI(NFSN)+1,iset), UPD(JI(0)+1,iset))
269 IF (NEFF .EQ. 0) GOTO 88
270 5 DO 333 IFLV = 1, NEFF
271 CALL CtLhNSEVL (LhCtNSRHSP, IKNL, NX, NITR, JT, DT, TIN,
272 > NEFF, UPD(JI( IFLV)+2,iset), UPD(JI( IFLV)+1,iset))
273 IF (IFLV .LE. MXVAL)
274 > CALL CtLhNSEVL (LhCtNSRHSM, IKNL, NX, NITR, JT, DT, TIN,
275 > NEFF, UPD(JI(-IFLV)+2,iset), UPD(JI(-IFLV)+1,iset))
276 DO 55 IS = 0, NITR
277 DO 56 IX = 0, NX
278 TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV),iset)
279 TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN),iset) / NEFF
280 TP = TP + TS
281 IF (IKNL .GT. 0) TP = MAX (TP, DZER)
282 IF (IFLV .LE. MXVAL) THEN
283 TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV),iset)
284 IF (IKNL .GT. 0) THEN
285 TM = MAX (TM, DZER)
286 TP = MAX (TP, TM)
287 EndIf
288 Else
289 TM = 0.
290 EndIf
291 UPD (JI( IFLV) + IS*(NX+1) + IX + 1,iset) = (TP + TM)/2.
292 UPD (JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = (TP - TM)/2.
293 56 CONTINUE
294 55 CONTINUE
295333 CONTINUE
296 DO 334 IFLV = NEFF + 1, NFMX
297 DO 57 IS = 0, NITR
298 DO 58 IX = 0, NX
299 UPD(JI( IFLV) + IS*(NX+1) + IX + 1,iset) = 0
300 UPD(JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = 0
301 58 CONTINUE
302 57 CONTINUE
303 334 CONTINUE
304 88 CONTINUE
305 IF (NFMX .EQ. NEFF) GOTO 21
306 DO 335 IFLV = -NFMX, NFMX+1
307 JI(IFLV) = JI(IFLV) + NITR * (NX+1)
308 335 CONTINUE
309 CALL CtLhHQRK (NX, TT, NEFF+1, UPD(JI(0)+2,iset),
310 > UPD(JI(NEFF+1)+2,iset))
311 DO 32 IZ = 1, NX
312 QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1,iset)
313 UPD (JI(NFSN)+IZ+1,iset) = UPD (JI(NFSN)+IZ+1,iset)
314 > + QRKP (NEFF+1)
315 VS00 = UPD (JI(NFSN)+IZ+1,iset) / (NEFF+1)
316 UPD(JI( NEFF+1) + IZ + 1,iset) = QRKP(NEFF+1) - VS00
317 DO 321 IFL = 1, NEFF
318 A = UPD(JI( IFL)+IZ+1,iset)
319 B = UPD(JI(-IFL)+IZ+1,iset)
320 QRKP(IFL) = A + B
321 UPD(JI( IFL)+IZ+1,iset) = QRKP(IFL) - VS00
322 IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1,iset) = A - B
323 321 CONTINUE
324 32 CONTINUE
325 21 CONTINUE
326 if(ientry .eq. 1) then
327 nxsave = nx
328 ntsave = nt
329 jtsave = jt
330 ngsave = ng
331 xcrsave = xcr
332 xminsave = xmin
333 qinisave = qini
334 qmaxsave = qmax
335 endif
336 if((nx .ne. nxsave) .or.
337 & (nt .ne. ntsave) .or.
338 & (jt .ne. jtsave) .or.
339 & (ng .ne. ngsave) .or.
340 & (xcr .ne. xcrsave) .or.
341 & (xmin .ne. xminsave) .or.
342 & (qini .ne. qinisave) .or.
343 & (qmax .ne. qmaxsave)) then
344 write(6,669) nx, nt, jt, ng, xcr, xmin,
345 & qini, qmax, ientry
346669 format(1x,'evolve.f: nx,nt,jt,ng=',4i4,
347 & ' xcr,xmin=',2f9.6,
348 & ' qini, qmax',f7.4,1x,e12.5,' ientry=',i6)
349 nxsave = nx
350 ntsave = nt
351 jtsave = jt
352 ngsave = ng
353 qinisave = qini
354 qmaxsave = qmax
355 xcrsave = xcr
356 xminsave = xmin
357 endif
358 ixshow = 0
359 do i = 0, mxx
360 if(xvsave(i) .ne. xv(i)) then
361 ixshow = 1
362 xvsave(i) = xv(i)
363 endif
364 enddo
365 iqshow = 0
366 itshow = 0
367 do i = 0, mxq
368 if(qvsave(i) .ne. qv(i)) then
369 iqshow = 1
370 qvsave(i) = qv(i)
371 endif
372 if(tvsave(i) .ne. tv(i)) then
373 itshow = 1
374 tvsave(i) = tv(i)
375 endif
376 enddo
377 Return
378 End
379 FUNCTION CtLhFINTRP (FF, X0, DX, NX, XV, ERR, IR)
380 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
381 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
382 PARAMETER (MX = 3)
383 DIMENSION FF (0:NX), XX(MX)
384 DATA SML, XX / 1.D-5, 0., 1.0, 2.0 /
385 DATA IW1, IW3, IW5 / 3 * 0 /
386 IR = 0
387 X = XV
388 ERR = 0.
389 ANX = NX
390 CtLhFINTRP = 0.
391 IF (NX .LT. 1) THEN
392 CALL CtLhWARNI(IW1, 'Nx < 1, error in CtLhFINTRP.',
393 > 'NX', NX, 1, 256, 1)
394 IR = 1
395 RETURN
396 ELSE
397 MNX = MIN(NX+1, MX)
398 ENDIF
399 IF (DX .LE. 0) THEN
400 CALL CtLhWARNR(IW3, 'DX < 0, error in CtLhFINTRP.',
401 > 'DX', DX, D0, D1, 1)
402 IR = 2
403 RETURN
404 ENDIF
405 XM = X0 + DX * NX
406 IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN
407 CALL CtLhWARNR(IW5,
408 > 'X out of range in CtLhFINTRP, Extrapolation used.',
409 > 'X',X,X0,XM,1)
410 IR = 3
411 ENDIF
412 TX = (X - X0) / DX
413 IF (TX .LE. 1.) THEN
414 IX = 0
415 ELSEIF (TX .GE. ANX-1.) THEN
416 IX = NX - 2
417 ELSE
418 IX = TX
419 ENDIF
420 DDX = TX - IX
421 CALL CtLhRATINT (XX, FF(IX), MNX, DDX, TEM, ERR)
422 CtLhFINTRP = TEM
423 RETURN
424 END
425 FUNCTION CtLhGausInt(F,XL,XR,AERR,RERR,ERR,IRT)
426 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
427 DIMENSION XLIMS(100), R(93), W(93)
428 INTEGER PTR(4),NORD(4)
429 external f
430 DATA PTR,NORD/4,10,22,46, 6,12,24,48/
431 DATA R/.2386191860,.6612093865,.9324695142,
432 1 .1252334085,.3678314990,.5873179543,.7699026742,.9041172563,
433 1 .9815606342,.0640568929,.1911188675,.3150426797,.4337935076,
434 1 .5454214714,.6480936519,.7401241916,.8200019860,.8864155270,
435 1 .9382745520,.9747285560,.9951872200,.0323801710,.0970046992,
436 1 .1612223561,.2247637903,.2873624873,.3487558863,.4086864820,
437 1 .4669029048,.5231609747,.5772247261,.6288673968,.6778723796,
438 1 .7240341309,.7671590325,.8070662040,.8435882616,.8765720203,
439 1 .9058791367,.9313866907,.9529877032,.9705915925,.9841245837,
440 1 .9935301723,.9987710073,.0162767488,.0488129851,.0812974955,
441 1 .1136958501,.1459737146,.1780968824,.2100313105,.2417431561,
442 1 .2731988126,.3043649444,.3352085229,.3656968614,.3957976498,
443 1 .4254789884,.4547094222,.4834579739,.5116941772,.5393881083,
444 1 .5665104186,.5930323648,.6189258401,.6441634037,.6687183100,
445 1 .6925645366,.7156768123,.7380306437,.7596023411,.7803690438,
446 1 .8003087441,.8194003107,.8376235112,.8549590334,.8713885059,
447 1 .8868945174,.9014606353,.9150714231,.9277124567,.9393703398,
448 1 .9500327178,.9596882914,.9683268285,.9759391746,.9825172636,
449 1 .9880541263,.9925439003,.9959818430,.9983643759,.9996895039/
450 DATA W/.4679139346,.3607615730,.1713244924,
451 1 .2491470458,.2334925365,.2031674267,.1600783285,.1069393260,
452 1 .0471753364,.1279381953,.1258374563,.1216704729,.1155056681,
453 1 .1074442701,.0976186521,.0861901615,.0733464814,.0592985849,
454 1 .0442774388,.0285313886,.0123412298,.0647376968,.0644661644,
455 1 .0639242386,.0631141923,.0620394232,.0607044392,.0591148397,
456 1 .0572772921,.0551995037,.0528901894,.0503590356,.0476166585,
457 1 .0446745609,.0415450829,.0382413511,.0347772226,.0311672278,
458 1 .0274265097,.0235707608,.0196161605,.0155793157,.0114772346,
459 1 .0073275539,.0031533461,.0325506145,.0325161187,.0324471637,
460 1 .0323438226,.0322062048,.0320344562,.0318287589,.0315893308,
461 1 .0313164256,.0310103326,.0306713761,.0302999154,.0298963441,
462 1 .0294610900,.0289946142,.0284974111,.0279700076,.0274129627,
463 1 .0268268667,.0262123407,.0255700360,.0249006332,.0242048418,
464 1 .0234833991,.0227370697,.0219666444,.0211729399,.0203567972,
465 1 .0195190811,.0186606796,.0177825023,.0168854799,.0159705629,
466 1 .0150387210,.0140909418,.0131282296,.0121516047,.0111621020,
467 1 .0101607705,.0091486712,.0081268769,.0070964708,.0060585455,
468 1 .0050142027,.0039645543,.0029107318,.0018539608,.0007967921/
469 DATA TOLABS,TOLREL,NMAX/1.E-35,5.E-4,100/
470 TOLABS=AERR
471 TOLREL=RERR
472
473 CtLhGausInt=0.
474 NLIMS=2
475 XLIMS(1)=XL
476 XLIMS(2)=XR
47710 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0
478 BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0
479 TVAL=0.
480 DO 15 I=1,3
48115 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
482 TVAL=TVAL*AA
483 DO 25 J=1,4
484 VAL=0.
485 DO 20 I=PTR(J),PTR(J)-1+NORD(J)
48620 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
487 VAL=VAL*AA
488 TOL=MAX(TOLABS,TOLREL*ABS(VAL))
489 IF (ABS(TVAL-VAL).LT.TOL) THEN
490 CtLhGausInt=CtLhGausInt+VAL
491 NLIMS=NLIMS-2
492 IF (NLIMS.NE.0) GO TO 10
493 RETURN
494 END IF
49525 TVAL=VAL
496 IF (NMAX.EQ.2) THEN
497 CtLhGausInt=VAL
498 RETURN
499 END IF
500 IF (NLIMS.GT.(NMAX-2)) THEN
501 write(*,50) CtLhGausInt,NMAX,BB-AA,BB+AA
502 RETURN
503 END IF
504 XLIMS(NLIMS+1)=BB
505 XLIMS(NLIMS+2)=BB+AA
506 XLIMS(NLIMS)=BB
507 NLIMS=NLIMS+2
508 GO TO 10
50950 FORMAT (' CtLhGausInt FAILS, CtLhGausInt,NMAX,XL,XR=',
510 > G15.7,I5,2G15.7)
511 END
512 SUBROUTINE CtLhHINTEG (NX, F, H)
513 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
514 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
515 PARAMETER (MXPN = MXF * 2 + 2)
516 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
517 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
518 COMMON / LhCtHINTEC / GH(NDG, MXX)
519 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
520 DIMENSION F(NX), H(NX), G(MXX)
521 DZ = 1D0 / (NX-1)
522 DO 20 I = 1, NX-2
523 NP = NX - I + 1
524 TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2)
525 DO 30 KZ = 3, NP
526 IY = I + KZ - 1
527 W = XA(I,1) / XA(IY,1)
528 G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W)
529 30 CONTINUE
530 HTEM = CtLhSMPSNA (NP-2, DZ, G(3), ERR)
531 TEM1 = F(I) * ELY(I)
532 H(I) = TEM + HTEM + TEM1
533 20 CONTINUE
534 H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0))
535 H(NX) = 0
536 RETURN
537 END
538 SUBROUTINE CtLhHQRK (NX, TT, NQRK, Y, F)
539 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
540 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
541 DIMENSION Y(NX), F(NX)
542 IF (NX .GT. 1) GOTO 11
543 11 CONTINUE
544 DO 230 IZ = 1, NX
545 IF (NX .GT. 1) THEN
546 F(IZ) = 0
547 GOTO 230
548 EndIf
549 230 CONTINUE
550 RETURN
551 END
552 SUBROUTINE CtLhINTEGR (NX, M, F, G, IR)
553 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
554 CHARACTER MSG*80
555 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
556 PARAMETER (MXPN = MXF * 2 + 2)
557 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
558 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
559 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
560 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
561 DIMENSION F(NX), G(NX)
562 DATA IWRN1, IWRN2 / 0, 0 /
563 IRR = 0
564 IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN
565 MSG = 'NX out of range in CtLhINTEGR call'
566 CALL CtLhWARNI (IWRN1, MSG, 'NX', NX, 0, MXX, 0)
567 IRR = 1
568 EndIf
569 IF (M .LT. M1 .OR. M .GT. M2) THEN
570 MSG ='Exponent M out of range in CtLhINTEGR'
571 CALL CtLhWARNI (IWRN2, MSG, 'M', M, M1, M2, 1)
572 IRR = 2
573 EndIf
574 G(NX) = 0D0
575 TEM = H(1, NX-1, -M) * F(NX-2) + H(2, NX-1, -M) * F(NX-1)
576 > + H(3, NX-1, -M) * F(NX)
577 IF (M .EQ. 0) THEN
578 G(NX-1) = TEM
579 Else
580 G(NX-1) = TEM * XA(NX-1, M)
581 EndIf
582 DO 10 I = NX-2, 2, -1
583 TEM = TEM + H(1,I,-M)*F(I-1) + H(2,I,-M)*F(I)
584 > + H(3,I,-M)*F(I+1) + H(4,I,-M)*F(I+2)
585 IF (M .EQ. 0) THEN
586 G(I) = TEM
587 Else
588 G(I) = TEM * XA(I, M)
589 EndIf
590 10 CONTINUE
591 TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3)
592 IF (M .EQ. 0) THEN
593 G(1) = TEM
594 Else
595 G(1) = TEM * XA(1, M)
596 EndIf
597 IR = IRR
598 RETURN
599 END
600 SUBROUTINE CtLhKERNEL
601 >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT)
602 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
603 PARAMETER (PI = 3.141592653589793d0, PI2 = PI ** 2)
604 PARAMETER (D0 = 0.0, D1 = 1.0)
605 DATA CF, CG, TR, IWRN / 1.33333333333333d0, 3.0d0, 0.5d0, 0 /
606 IRT = 0
607 TRNF = TR * NFL
608 X = XX
609 IF (X .LE. 0. .OR. X .GE. 1.) THEN
610 CALL CtLhWARNR(IWRN, 'X out of range in CtLhKERNEL', 'X', X,
611 > D0, D1, 1)
612 IRT = 1
613 RETURN
614 EndIf
615 XI = 1./ X
616 X2 = X ** 2
617 XM1I = 1./ (1.- X)
618 XP1I = 1./ (1.+ X)
619 XLN = LOG (X)
620 XLN2 = XLN ** 2
621 XLN1M = LOG (1.- X)
622 SPEN2 = CtLhSPENC2 (X)
623 FFP = (1.+ X2) * XM1I
624 FGP = (2.- 2.* X + X2) / X
625 GFP = 1. - 2.* X + 2.* X2
626 GGP = XM1I + XI - 2. + X - X2
627 FFM = (1.+ X2) * XP1I
628 FGM = - (2.+ 2.* X + X2) / X
629 GFM = 1. + 2.* X + 2.* X2
630 GGM = XP1I - XI - 2. - X - X2
631 FF1 = CF * FFP * (1.- X)
632 FG1 = CF * FGP * X
633 GF1 = 2.* TRNF * GFP
634 GG1 = 2.* CG * GGP * X * (1.-X)
635 PCF2 = -2.* FFP *XLN*XLN1M - (3.*XM1I + 2.*X)*XLN
636 > - (1.+X)/2.*XLN2 - 5.*(1.-X)
637 PCFG = FFP * (XLN2 + 11.*XLN/3.+ 67./9.- PI**2 / 3.)
638 > + 2.*(1.+X) * XLN + 40.* (1.-X) / 3.
639 PCFT = (FFP * (- XLN - 5./3.) - 2.*(1.-X)) * 2./ 3.
640 PQQB = 2.* FFM * SPEN2 + 2.*(1.+X)*XLN + 4.*(1.-X)
641 PQQB = (CF**2-CF*CG/2.) * PQQB
642 PQQ2 = CF**2 * PCF2 + CF*CG * PCFG / 2. + CF*TRNF * PCFT
643 PNSP = (PQQ2 + PQQB) * (1.-X)
644 PNSM = (PQQ2 - PQQB) * (1.-X)
645 FFCF2 = - 1. + X + (1.- 3.*X) * XLN / 2. - (1.+ X) * XLN2 / 2.
646 > - FFP * (3.* XLN / 2. + 2.* XLN * XLN1M)
647 > + FFM * 2.* SPEN2
648 FFCFG = 14./3.* (1.-X)
649 > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.)
650 > - FFM * SPEN2
651 FFCFT = - 16./3. + 40./3.* X + (10.* X + 16./3.* X2 + 2.) * XLN
652 > - 112./9.* X2 + 40./9./X - 2.* (1.+ X) * XLN2
653 > - FFP * (10./9. + 2./3. * XLN)
654 FGCF2 = - 5./2.- 7./2.* X + (2.+ 7./2.* X) * XLN + (X/2.-1.)*XLN2
655 > - 2.* X * XLN1M
656 > - FGP * (3.* XLN1M + XLN1M ** 2)
657 FGCFG = 28./9. + 65./18.* X + 44./9. * X2 - (12.+ 5.*X + 8./3.*X2)
658 > * XLN + (4.+ X) * XLN2 + 2.* X * XLN1M
659 > + FGP * (-2.*XLN*XLN1M + XLN2/2. + 11./3.*XLN1M + XLN1M**2
660 > - PI2/6. + 0.5)
661 > + FGM * SPEN2
662 FGCFT = -4./3.* X - FGP * (20./9.+ 4./3.*XLN1M)
663 GFCFT = 4.- 9.*X + (-1.+ 4.*X)*XLN + (-1.+ 2.*X)*XLN2 + 4.*XLN1M
664 > + GFP * (-4.*XLN*XLN1M + 4.*XLN + 2.*XLN2 - 4.*XLN1M
665 > + 2.*XLN1M**2 - 2./3.* PI2 + 10.)
666 GFCGT = 182./9.+ 14./9.*X + 40./9./X + (136./3.*X - 38./3.)*XLN
667 > - 4.*XLN1M - (2.+ 8.*X)*XLN2
668 > + GFP * (-XLN2 + 44./3.*XLN - 2.*XLN1M**2 + 4.*XLN1M
669 > + PI2/3. - 218./9.)
670 > + GFM * 2. * SPEN2
671 GGCFT = -16.+ 8.*X + 20./3.*X2 + 4./3./X + (-6.-10.*X)*XLN
672 > - 2.* (1.+ X) * XLN2
673 GGCGT = 2.- 2.*X + 26./9.*X2 - 26./9./X - 4./3.*(1.+X)*XLN
674 > - GGP * 20./9.
675 GGCG2 = 27./2.*(1.-X) + 67./9.*(X2-XI) + 4.*(1.+X)*XLN2
676 > + (-25.+ 11.*X - 44.*X2)/3.*XLN
677 > + GGP * (67./9.- 4.*XLN*XLN1M + XLN2 - PI2/3.)
678 > + GGM * 2.* SPEN2
679 FF2 = CF * TRNF * FFCFT + CF ** 2 * FFCF2 + CF * CG * FFCFG
680 FG2 = CF * TRNF * FGCFT + CF ** 2 * FGCF2 + CF * CG * FGCFG
681 GF2 = CF * TRNF * GFCFT + CG * TRNF * GFCGT
682 GG2 = CF * TRNF * GGCFT + CG ** 2 * GGCG2 + CG * TRNF * GGCGT
683 XLG = (LOG(1./(1.-X)) + 1.)
684 XG2 = XLG ** 2
685 FF2 = FF2 * X * (1.- X)
686 FG2 = FG2 * X / XG2
687 GF2 = GF2 * X / XG2
688 GG2 = GG2 * X * (1.- X)
689 RETURN
690 END
691 SUBROUTINE CtLhLAMCWZ
692 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
693 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
694 LOGICAL SET
695 CALL CtLhSETL1 (NF, AL)
696 END
697 FUNCTION LhCtNAMQCD(NNAME)
698 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
699 CHARACTER NNAME*(*), NAME*8
700 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
701 LOGICAL SET
702 CHARACTER ONECH*(1)
703 ONECH = '0'
704 IASC0 = ICHAR(ONECH)
705 NAME = NNAME
706 LhCtNAMQCD=0
707 IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR.
708 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') )
709 2 LhCtNAMQCD=1
710 IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR.
711 1 (NAME .EQ. '# FL') )
712 2 LhCtNAMQCD=2
713 DO 10 I=1, 9
714 IF (NAME .EQ. 'M'//CHAR(I+IASC0))
715 1 LhCtNAMQCD=I+2
71610 CONTINUE
717 DO 20 I= 0, NF
718 IF (NAME .EQ. 'LAM'//CHAR(I+IASC0))
719 1 LhCtNAMQCD=I+13
72020 CONTINUE
721 IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') LhCtNAMQCD = 24
722 RETURN
723 END
724 FUNCTION LhCtNFL(AMU)
725 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
726 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
727 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
728 LOGICAL SET
729 IF (.NOT. SET) CALL CtLhLAMCWZ
730 LhCtNFL = NF - NHQ
731 IF ((LhCtNFL .EQ. NF) .OR. (AMU .LE. AMN)) GOTO 20
732 DO 10 I = NF - NHQ + 1, NF
733 IF (AMU .GE. AMHAT(I)) THEN
734 LhCtNFL = I
735 ELSE
736 GOTO 20
737 ENDIF
73810 CONTINUE
73920 RETURN
740 END
741 SUBROUTINE CtLhNSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN)
742 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
743 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
744 PARAMETER (MXPN = MXF * 2 + 2)
745 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
746 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
747 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
748 DIMENSION U0(NX), UN(0:NX, 0:NT)
749 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
750 external rhs
751 DDT = DT / JT
752 IF (NX .GT. MXX) THEN
753 WRITE (*,*) 'Nx =', NX, ' greater than Max pts in CtLhNSEVL.'
754 STOP 'Program stopped in CtLhNSEVL'
755 EndIf
756 TMD = TIN + DT * NT / 2.
757 AMU = EXP(TMD)
758 TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
759 TLAM = TMD - TEM
760 DO 9 IX = 1, NX
761 UN(IX, 0) = U0(IX)
762 9 CONTINUE
763 UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1)
764 TT = TIN
765 DO 10 IZ = 1, NX
766 Y0(IZ) = U0(IZ)
767 10 CONTINUE
768 DO 20 IS = 1, NT
769 DO 202 JS = 1, JT
770 IRND = (IS-1) * JT + JS
771 IF (IRND .EQ. 1) THEN
772 CALL RHS (TT, Neff, Y0, F0)
773 DO 250 IZ = 1, NX
774 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
775 250 CONTINUE
776 TT = TT + DDT
777 CALL RHS (TT, NEFF, Y0, F1)
778 DO 251 IZ = 1, NX
779 Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
780 251 CONTINUE
781 Else
782 CALL RHS (TT, NEFF, Y1, F1)
783 DO 252 IZ = 1, NX
784 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
785 252 CONTINUE
786 TT = TT + DDT
787 CALL RHS (TT, NEFF, YP, FP)
788 DO 253 IZ = 1, NX
789 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
790 F0(IZ) = F1(IZ)
791 253 CONTINUE
792 EndIf
793 202 CONTINUE
794 DO 260 IZ = 1, NX
795 UN (IZ, IS) = Y1(IZ)
796 260 CONTINUE
797 UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
798 20 CONTINUE
799 RETURN
800 END
801 SUBROUTINE LhCtNSRHSM (TT, NEFF, FI, FO)
802 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
803 LOGICAL LSTX
804 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
805 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
806 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
807 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
808 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
809 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
810 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
811 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
812 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
813 DIMENSION G1(MXX), FI(NX), FO(NX)
814 DIMENSION W0(MXX), W1(MXX), WH(MXX)
815 S = EXP(TT)
816 Q = AL * EXP (S)
817 CPL = CtLhALPI(Q)
818 CPL2= CPL ** 2 / 2. * S
819 CPL = CPL * S
820 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
821 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
822 CALL CtLhHINTEG (NX, FI, WH)
823 DO 230 IZ = 1, NX
824 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
825 FO(IZ) = CPL * FO(IZ)
826 230 CONTINUE
827 IF (IKNL .EQ. 2) THEN
828 DZ = 1./ (NX - 1)
829 DO 21 IX = 1, NX-1
830 X = XV(IX)
831 NP = NX - IX + 1
832 IS = NP
833 DO 31 KZ = 2, NP
834 IY = IX + KZ - 1
835 IT = NX - IY + 1
836 XY = ZZ (IS, IT)
837 G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX))
838 31 CONTINUE
839 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
840 TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2
841 FO(IX) = FO(IX) + TMP2
842 21 CONTINUE
843 EndIf
844 RETURN
845 END
846 SUBROUTINE LhCtNSRHSP (TT, NEFF, FI, FO)
847 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
848 LOGICAL LSTX
849 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
850 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
851 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
852 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
853 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
854 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
855 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
856 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
857 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
858 DIMENSION G1(MXX), FI(NX), FO(NX)
859 DIMENSION W0(MXX), W1(MXX), WH(MXX)
860 S = EXP(TT)
861 Q = AL * EXP (S)
862 CPL = CtLhALPI(Q)
863 CPL2= CPL ** 2 / 2. * S
864 CPL = CPL * S
865 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
866 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
867 CALL CtLhHINTEG (NX, FI, WH)
868 DO 230 IZ = 1, NX
869 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
870 FO(IZ) = CPL * FO(IZ)
871 230 CONTINUE
872 IF (IKNL .EQ. 2) THEN
873 DZ = 1./ (NX - 1)
874 DO 21 IX = 1, NX-1
875 X = XV(IX)
876 NP = NX - IX + 1
877 DO 31 KZ = 2, NP
878 IY = IX + KZ - 1
879 XY = ZZ (NX-IX+1, NX-IY+1)
880 G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX))
881 31 CONTINUE
882 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
883 TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + ZQQB)) * CPL2
884 FO(IX) = FO(IX) + TMP2
885 21 CONTINUE
886 EndIf
887 RETURN
888 END
889 FUNCTION CtLhPARDIS (IPRTN, XX, QQ)
890 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
891 include 'parmsetup.inc'
892 Character Msg*80
893 LOGICAL LSTX
894 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
895 PARAMETER (MXPN = MXF * 2 + 2)
896 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
897 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
898 PARAMETER (Smll = 1D-9)
899 parameter(nqvec = 4)
900 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
901 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
902 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
903 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
904 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
905 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
906 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
907 COMMON / LhCtCOMQMS / VALQMS(9)
908 dimension fvec(4), fij(4)
909 dimension xvpow(0:mxx)
910 Data Iwrn1, Iwrn2, Iwrn3, OneP / 3*0, 1.00001 /
911 data xpow / 0.3d0 / !**** choice of interpolation variable
912 data nxsave / 0 /
913 save xvpow, nxsave
914 save xlast, qlast
915 save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
916 save const1 , const2, const3, const4, const5, const6
917 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
918c
919 call getnset(iset)
920c
921 if(nx .ne. nxsave) then
922 xvpow(0) = 0D0
923 do i = 1, nx
924 xvpow(i) = xv(i)**xpow
925 enddo
926 nxsave = nx
927 endif
928
929 X = XX
930 Q = QQ
931
932c enforce threshold early to improve speed...
933 ii = iabs(IPRTN)
934 if(ii .ne. 0) then
935 if(QQ .lt. VALQMS(ii) ) then
936 ctlhpardis = 0.d0
937 return
938 endif
939 endif
940
941c force pardis = 0.0d0 at exactly =1.0d0 - added mrw 10/May/06
942 if(xx .eq. 1.0d0) then
943 ctlhpardis = 0.0d0
944 return
945 endif
946
947c skip the initialization in x if same as in the previous call.
948 if(x .eq. xlast) goto 100
949 xlast = x
950
951 JLx = -1
952 JU = Nx+1
953 11 If (JU-JLx .GT. 1) Then
954 JM = (JU+JLx) / 2
955 If (X .Ge. XV(JM)) Then
956 JLx = JM
957 Else
958 JU = JM
959 Endif
960 Goto 11
961 Endif
962 If (JLx .LE. -1) Then
963 Print '(A,1pE12.4)','Severe error: x <= 0 in ParDis x=', x
964 Stop
965 ElseIf (JLx .Eq. 0) Then
966 Jx = 0
967 Msg = '0 < X < Xmin in ParDis; extrapolation used!'
968 CALL CtLhWARNR (IWRN1, Msg, 'X', X, Xmin, 1D0, 1)
969 Elseif (JLx .LE. Nx-2) Then
970 Jx = JLx - 1
971 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
972 Jx = JLx - 2
973 Else
974 Print '(A,1pE12.4)','Severe error: x > 1 in ParDis x=', x
975 Stop
976 Endif
977 ss = x**xpow
978 If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
979 svec1 = xvpow(jx)
980 svec2 = xvpow(jx+1)
981 svec3 = xvpow(jx+2)
982 svec4 = xvpow(jx+3)
983 s12 = svec1 - svec2
984 s13 = svec1 - svec3
985 s23 = svec2 - svec3
986 s24 = svec2 - svec4
987 s34 = svec3 - svec4
988 sy2 = ss - svec2
989 sy3 = ss - svec3
990 const1 = s13/s23
991 const2 = s12/s23
992 const3 = s34/s23
993 const4 = s24/s23
994 s1213 = s12 + s13
995 s2434 = s24 + s34
996 sdet = s12*s34 - s1213*s2434
997 tmp = sy2*sy3/sdet
998 const5 = (s34*sy2-s2434*sy3)*tmp/s12
999 const6 = (s1213*sy2-s12*sy3)*tmp/s34
1000 EndIf
1001
1002100 continue
1003
1004c skip the initialization in q if same as in the previous call.
1005 if(q .eq. qlast) goto 110
1006 qlast = q
1007
1008 tt = log(log(Q/Al))
1009
1010 JLq = -1
1011 JU = NT+1
1012 12 If (JU-JLq .GT. 1) Then
1013 JM = (JU+JLq) / 2
1014 If (Q .GE. QV(JM)) Then
1015 JLq = JM
1016 Else
1017 JU = JM
1018 Endif
1019 Goto 12
1020 Endif
1021 If (JLq .LE. 0) Then
1022 Jq = 0
1023 If (JLq .LT. 0) Then
1024 Msg = 'Q < Q0 in ParDis; extrapolation used!'
1025 CALL CtLhWARNR (IWRN2, Msg, 'Q', Q, Qini, 1D0, 1)
1026 EndIf
1027 Elseif (JLq .LE. Nt-2) Then
1028 Jq = JLq - 1
1029 Else
1030 Jq = Nt - 3
1031 If (JLq .GE. Nt) Then
1032 Msg = 'Q > Qmax in ParDis; extrapolation used!'
1033 CALL CtLhWARNR (IWRN3, Msg, 'Q', Q, Qmax, 1D0, 1)
1034 Endif
1035 Endif
1036
1037 If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
1038 tvec1 = Tv(jq)
1039 tvec2 = Tv(jq+1)
1040 tvec3 = Tv(jq+2)
1041 tvec4 = Tv(jq+3)
1042 t12 = tvec1 - tvec2
1043 t13 = tvec1 - tvec3
1044 t23 = tvec2 - tvec3
1045 t24 = tvec2 - tvec4
1046 t34 = tvec3 - tvec4
1047 ty2 = tt - tvec2
1048 ty3 = tt - tvec3
1049 tmp1 = t12 + t13
1050 tmp2 = t24 + t34
1051 tdet = t12*t34 - tmp1*tmp2
1052 EndIf
1053
1054110 continue
1055
1056 jtmp = ((IPRTN + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
1057 Do it = 1, nqvec
1058 J1 = jtmp + it*(NX+1)
1059 If (Jx .Eq. 0) Then
1060 fij(1) = 0
1061 fij(2) = Upd(J1+1,iset) * Xa(1,2)
1062 fij(3) = Upd(J1+2,iset) * Xa(2,2)
1063 fij(4) = Upd(J1+3,iset) * Xa(3,2)
1064 Call CtLhPolint4 (XVpow(0), Fij(1), 4, ss, Fx, Dfx)
1065
1066 If (x .GT. 0D0) Fvec(it) = Fx / x**2
1067 ElseIf (JLx .Eq. Nx-1) Then
1068 Call CtLhPolint4 (XVpow(Nx-3), Upd(J1,iset), 4, ss, Fx, Dfx)
1069 Fvec(it) = Fx
1070 Else
1071 sf2 = Upd(J1+1,iset)
1072 sf3 = Upd(J1+2,iset)
1073 Fvec(it) = (const5*(Upd(J1,iset)
1074 & - sf2*const1 + sf3*const2)
1075 & + const6*(Upd(J1+3,iset)
1076 & + sf2*const3 - sf3*const4)
1077 & + sf2*sy3 - sf3*sy2) / s23
1078 Endif
1079 enddo
1080 If (JLq .LE. 0) Then
1081 Call CtLhPolint4 (TV(0), Fvec(1), 4, tt, ff, Dfq)
1082 ElseIf (JLq .GE. Nt-1) Then
1083 Call CtLhPolint4 (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
1084 Else
1085 tf2 = fvec(2)
1086 tf3 = fvec(3)
1087 g1 = ( tf2*t13 - tf3*t12) / t23
1088 g4 = (-tf2*t34 + tf3*t24) / t23
1089 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
1090 & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
1091 ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
1092 EndIf
1093 CtLhPARDIS = ff
1094 Return
1095 End
1096
1097 SUBROUTINE CtLhPARPDF (IACT, NAME, VALUE, IRET)
1098 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1099 CHARACTER NAME*(*), Uname*10
1100 LOGICAL START1
1101 DATA ILEVEL, LRET / 1, 1 /
1102 JRET = IRET
1103 CALL CtLhUPC (NAME, Ln, Uname)
1104 IF (IACT .EQ. 0 .OR. IACT .EQ. 4) then
1105c > IVALUE = NINT (VALUE) !tentatively remove this since it seems not to be used
1106 print *,'Fatal error: iact=',iact
1107 stop
1108 ENDIF
1109 START1 = (IACT .NE. 1) .AND. (IACT .NE. 2)
1110c prepare to remove this stuff, since I think IACT=1 or 2 always
1111 if(start1) then
1112 print *,'Fatal error: start1=',start1
1113 stop
1114 endif
1115 IF (START1) ILEVEL = 1
1116 GOTO (1, 2), ILEVEL
1117 1 START1 = .TRUE.
1118 ILEVEL = 0
1119 CALL CtLhParQcd (IACT, Uname(1:Ln), VALUE, JRET)
1120 IF (JRET .EQ. 1) GOTO 11
1121 IF (JRET .EQ. 2) GOTO 12
1122 IF (JRET .EQ. 3) GOTO 13
1123 IF (JRET .GT. 4) GOTO 15
1124 ILEVEL = ILEVEL + 1
1125 2 CALL CtLhEVLPAR (IACT, Uname(1:Ln), VALUE, JRET)
1126 IF (JRET .EQ. 1) GOTO 11
1127 IF (JRET .EQ. 2) GOTO 12
1128 IF (JRET .EQ. 3) GOTO 13
1129 IF (JRET .GT. 4) GOTO 15
1130 ILEVEL = ILEVEL + 1
1131 IF (.NOT. START1) GOTO 1
1132 IF (JRET .EQ. 0) GOTO 10
1133 9 CONTINUE
1134 GOTO 14
1135 10 CONTINUE
1136 11 CONTINUE
1137 12 CONTINUE
1138 13 CONTINUE
1139 14 CONTINUE
1140 15 CONTINUE
1141 IF (JRET .NE. 4) LRET = JRET
1142 IF (LRET.EQ.0 .OR. LRET.EQ.2 .OR. LRET.EQ.3) THEN
1143 PRINT *, 'Error in CtLhPARPDF: IRET, IACT, NAME, VALUE =',
1144 > LRET, IACT, NAME, VALUE
1145 PRINT *, 'fatal error in CtLhparpdf'
1146 stop
1147 EndIf
1148 IRET= JRET
1149 RETURN
1150 100 FORMAT (/)
1151 END
1152 SUBROUTINE CtLhParQcd(IACT,NAME,VALUE,IRET)
1153 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1154 INTEGER IACT,IRET
1155 CHARACTER*(*) NAME
1156 IRET=1
1157 IF (IACT.EQ.0) THEN
1158 WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ',
1159 > '(i in 1 to 9), LAMi (i in 1 to NFL)'
1160 IRET=4
1161 ELSEIF (IACT.EQ.1) THEN
1162 CALL CtLhQCDSET (NAME,VALUE,IRET)
1163 ELSEIF (IACT.EQ.2) THEN
1164 CALL CtLhQCDGET (NAME,VALUE,IRET)
1165 ELSE
1166 IRET=3
1167 ENDIF
1168 RETURN
1169 END
1170 FUNCTION CtLhPFF1 (XX)
1171 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1172 LOGICAL LA, LB, LSTX
1173 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1174 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1175 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1176 PARAMETER (MX = 3)
1177 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1178 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1179 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1180 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1181 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1182 > PNSP(0:MXX), PNSM(0:MXX)
1183 SAVE
1184 DATA LA, LB / 2 * .FALSE. /
1185 LB = .TRUE.
1186 ENTRY CtLhTFF1(ZZ)
1187 LA = .TRUE.
1188 2 IF (LA .AND. .NOT.LB) THEN
1189 Z = ZZ
1190 X = CtLhXFRMZ (Z)
1191 Else
1192 X = XX
1193 EndIf
1194 IF (X .GE. D1) THEN
1195 CtLhPFF1 = 0
1196 RETURN
1197 ElseIF (X .GE. XMIN) THEN
1198 Z = CtLhZFRMX (X)
1199 TEM = CtLhFINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT)
1200 Else
1201 CALL CtLhPOLIN1 (XL, FF1(1), MX, X, TEM, ERR)
1202 EndIf
1203 IF (LA) THEN
1204 IF (LB) THEN
1205 CtLhPFF1 = TEM / (1.-X)
1206 LB =.FALSE.
1207 Else
1208 CtLhTFF1 = TEM / (1.-X) * CtLhDXDZ(Z)
1209 EndIf
1210 LA =.FALSE.
1211 Else
1212 IF (LB) THEN
1213 QFF1 = TEM
1214 LB =.FALSE.
1215 Else
1216 RFF1 = TEM * X / (1.-X)
1217 EndIf
1218 EndIf
1219 RETURN
1220 ENTRY CtLhFNSP (XX)
1221 X = XX
1222 IF (X .GE. D1) THEN
1223 CtLhFNSP = 0.
1224 RETURN
1225 ElseIF (X .GE. XMIN) THEN
1226 Z = CtLhZFRMX (X)
1227 TEM = CtLhFINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT)
1228 Else
1229 CALL CtLhPOLIN1 (XL, PNSP(1), MX, X, TEM, ERR)
1230 EndIf
1231 CtLhFNSP = TEM / (1.- X)
1232 RETURN
1233 ENTRY CtLhFNSM (XX)
1234 X = XX
1235 IF (X .GE. D1) THEN
1236 CtLhFNSM = 0.
1237 RETURN
1238 ElseIF (X .GE. XMIN) THEN
1239 Z = CtLhZFRMX (X)
1240 TEM = CtLhFINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT)
1241 Else
1242 CALL CtLhPOLIN1 (XL, PNSM(1), MX, X, TEM, ERR)
1243 EndIf
1244 CtLhFNSM = TEM / (1.- X)
1245 RETURN
1246 ENTRY CtLhRGG1 (XX)
1247 X = XX
1248 IF (X .GE. D1) THEN
1249 PGG1= 0
1250 RETURN
1251 ElseIF (X .GE. XMIN) THEN
1252 Z = CtLhZFRMX (X)
1253 TEM = CtLhFINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT)
1254 Else
1255 CALL CtLhPOLIN1 (XL, GG1(1), MX, X, TEM, ERR)
1256 EndIf
1257 IF (LA) THEN
1258 PGG1 = TEM / X / (1.-X)
1259 LA =.FALSE.
1260 Else
1261 IF (LB) THEN
1262 QGG1 = TEM / X
1263 LB =.FALSE.
1264 Else
1265 CtLhRGG1 = TEM / (1.-X)
1266 EndIf
1267 EndIf
1268 RETURN
1269 ENTRY CtLhRFF2 (XX)
1270 X = XX
1271 IF (X .GE. D1) THEN
1272 PFF2 = 0
1273 RETURN
1274 ElseIF (X .GE. XMIN) THEN
1275 Z = CtLhZFRMX (X)
1276 TEM = CtLhFINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT)
1277 Else
1278 CALL CtLhPOLIN1 (XL, FF2(1), MX, X, TEM, ERR)
1279 EndIf
1280 IF (LA) THEN
1281 PFF2 = TEM / X / (1.-X)
1282 LA =.FALSE.
1283 Else
1284 IF (LB) THEN
1285 QFF2 = TEM / X
1286 LB =.FALSE.
1287 Else
1288 CtLhRFF2 = TEM / (1.-X)
1289 EndIf
1290 EndIf
1291 RETURN
1292 ENTRY CtLhRGG2 (XX)
1293 X = XX
1294 IF (X .GE. D1) THEN
1295 PGG2 = 0
1296 RETURN
1297 ElseIF (X .GE. XMIN) THEN
1298 Z = CtLhZFRMX (X)
1299 TEM = CtLhFINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT)
1300 Else
1301 CALL CtLhPOLIN1 (XL, GG2(1), MX, X, TEM, ERR)
1302 EndIf
1303 IF (LA) THEN
1304 PGG2 = TEM / X / (1.-X)
1305 LA =.FALSE.
1306 Else
1307 IF (LB) THEN
1308 QGG2 = TEM / X
1309 LB =.FALSE.
1310 Else
1311 CtLhRGG2 = TEM / (1.-X)
1312 EndIf
1313 EndIf
1314 RETURN
1315 END
1316 SUBROUTINE CtLhPOLIN1 (XA,YA,N,X,Y,DY)
1317 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1318 PARAMETER (NMAX=10)
1319 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1320 NS=1
1321 DIF=ABS(X-XA(1))
1322 DO 11 I=1,N
1323 DIFT=ABS(X-XA(I))
1324 IF (DIFT.LT.DIF) THEN
1325 NS=I
1326 DIF=DIFT
1327 ENDIF
1328 C(I)=YA(I)
1329 D(I)=YA(I)
133011 CONTINUE
1331 Y=YA(NS)
1332 NS=NS-1
1333 DO 13 M=1,N-1
1334 DO 12 I=1,N-M
1335 HO=XA(I)-X
1336 HP=XA(I+M)-X
1337 W=C(I+1)-D(I)
1338 DEN=HO-HP
1339 DEN=W/DEN
1340 D(I)=HP*DEN
1341 C(I)=HO*DEN
134212 CONTINUE
1343 IF (2*NS.LT.N-M)THEN
1344 DY=C(NS+1)
1345 ELSE
1346 DY=D(NS)
1347 NS=NS-1
1348 ENDIF
1349 Y=Y+DY
135013 CONTINUE
1351 RETURN
1352 END
1353 SUBROUTINE CtLhQARRAY (NINI)
1354 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1355 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1356 PARAMETER (MXPN = MXF * 2 + 2)
1357 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
1358 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1359 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
1360 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1361 NCNT = 0
1362 IF (NT .GE. mxq) NT = mxq - 1
1363 S = LOG(QINI/AL)
1364 TINI = LOG(S)
1365 S = LOG(QMAX/AL)
1366 TMAX = LOG(S)
1367 1 DT0 = (TMAX - TINI) / float(NT)
1368 NINI = LhCtNFL(QINI)
1369 NFMX = LhCtNFL(QMAX)
1370 Call CtLhParQcd (2, 'ORDER', Ord, Ir)
1371 Call CtLhParQcd (2, 'ALAM', Al0, Ir)
1372 Call CtLhParQcd (2, 'NFL', Afl0, Ir)
1373 AFL = NfMx
1374 Call CtLhParQcd (1, 'NFL', AFL, Ir)
1375 Iordr = Nint (Ord)
1376 Ifl0 = Nint (Afl0)
1377 Call CtLhSetLam (Ifl0, Al0, Iordr)
1378 NG = NFMX - NINI + 1
1379 QIN = QINI
1380 QOUT = QINI
1381 S = LOG(QIN/AL)
1382 TIN = LOG(S)
1383 TLN(1) = TIN
1384 NTL(1) = 0
1385 QV(0) = QINI
1386 TV(0) = Tin
1387 DO 20 NEFF = NINI, NFMX
1388 ICNT = NEFF - NINI + 1
1389 IF (NEFF .LT. NFMX) THEN
1390 THRN = CtLhAMHATF (NEFF + 1)
1391 QOUN = MIN (QMAX, THRN)
1392 Else
1393 QOUN = QMAX
1394 EndIf
1395 IF (QOUN-QOUT .LE. 0.0001) THEN
1396 DT = 0
1397 NITR = 0
1398 Else
1399 QOUT = QOUN
1400 S = LOG(QOUT/AL)
1401 TOUT = LOG(S)
1402 TEM = TOUT - TIN
1403 NITR = INT (TEM / DT0) + 1
1404 DT = TEM / NITR
1405 EndIf
1406 DTN (ICNT) = DT
1407 NTN (ICNT) = NITR
1408 TLN (ICNT) = TIN
1409 NTL (ICNT+1) = NTL(ICNT) + NITR
1410 IF (NITR .NE. 0) THEN
1411 DO 205 I = 1, NITR
1412 TV (NTL(ICNT)+I) = TIN + DT * I
1413 S = EXP (TV(NTL(ICNT)+I))
1414 QV (NTL(ICNT)+I) = AL * EXP (S)
1415 205 CONTINUE
1416 EndIf
1417 QIN = QOUT
1418 TIN = TOUT
1419 20 CONTINUE
1420 NCNT = NCNT + 1
1421 NTP = NTL (NG + 1)
1422 ND = NTP - NT
1423 IF (NTP .GE. MXQ) THEN
1424 NT = MXQ - ND - NCNT
1425 GOTO 1
1426 EndIf
1427 NT = NTP
1428 RETURN
1429 END
1430 SUBROUTINE CtLhQCDGET(NAME,VALUE,IRET)
1431 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1432 CHARACTER*(*) NAME
1433 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1434 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1435 COMMON / LhCtCOMQMS / VALQMS(9)
1436 LOGICAL SET
1437 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1438 ICODE = LhCtNAMQCD(NAME)
1439 IRET = 1
1440 IF (ICODE .EQ. 1) THEN
1441 VALUE = AL
1442 ELSEIF (ICODE .EQ. 2) THEN
1443 VALUE = NF
1444 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 12)) THEN
1445 VALUE = VALQMS(ICODE - 2)
1446 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1447 VALUE = ALAM(ICODE - 13)
1448 ELSEIF (ICODE .EQ. 24) THEN
1449 VALUE = NORDER
1450 ELSE
1451 IRET=0
1452 ENDIF
1453 END
1454 SUBROUTINE CtLhQCDSET (NAME,VALUE,IRET)
1455 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1456 CHARACTER*(*) NAME
1457 COMMON / LhCtCOMQMS / VALQMS(9)
1458 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1459 LOGICAL SET
1460 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1461 IVALUE = NINT(VALUE)
1462 ICODE = LhCtNAMQCD(NAME)
1463 IF (ICODE .EQ. 0) THEN
1464 IRET=0
1465c print *,'warning empty CtLhQCDSET call: NAME=',
1466c & NAME,' VALUE=',VALUE
1467 ELSE
1468 IRET = 1
1469 SET = .FALSE.
1470 IF (ICODE .EQ. 1) THEN
1471 IF (VALUE.LE.0) GOTO 12
1472 AL=VALUE
1473 ELSEIF (ICODE .EQ. 2) THEN
1474 IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12
1475 NF = IVALUE
1476 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 11)) THEN
1477 IF (VALUE .LT. 0) GOTO 12
1478 Scle = Min (Value , VALQMS(ICODE - 2))
1479 AlfScle = CtLhALPI(Scle) * Pi
1480 VALQMS(ICODE - 2) = VALUE
1481 Call CtLhAlfSet (Scle, AlfScle)
1482 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1483 IF (VALUE .LE. 0) GOTO 12
1484 CALL CtLhSETL1 (ICODE-13, VALUE)
1485 ELSEIF (ICODE .EQ. 24) THEN
1486 IF ((IVALUE .LT. 1) .OR. (IVALUE .GT. 2)) GOTO 12
1487 NORDER = IVALUE
1488 ENDIF
1489 IF (.NOT. SET) CALL CtLhLAMCWZ
1490 ENDIF
1491 RETURN
1492 12 IRET=2
1493 RETURN
1494 END
1495 FUNCTION CtLhQZBRNT(FUNC, X1, X2, TOLIN, IRT)
1496 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1497 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
1498 external func
1499 TOL = ABS(TOLIN)
1500 A=X1
1501 B=X2
1502 FA=FUNC(A)
1503 FB=FUNC(B)
1504 IF(FB*FA.GT.0.) THEN
1505 WRITE (*, *) 'Root must be bracketed for CtLhQZBRNT.'
1506 IRT = 1
1507 ENDIF
1508 FC=FB
1509 DO 11 ITER=1,ITMAX
1510 IF(FB*FC.GT.0.) THEN
1511 C=A
1512 FC=FA
1513 D=B-A
1514 E=D
1515 ENDIF
1516 IF(ABS(FC).LT.ABS(FB)) THEN
1517 A=B
1518 B=C
1519 C=A
1520 FA=FB
1521 FB=FC
1522 FC=FA
1523 ENDIF
1524 TOL1=2.*EPS*ABS(B)+0.5*TOL
1525 XM=.5*(C-B)
1526 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
1527 CtLhQZBRNT=B
1528 RETURN
1529 ENDIF
1530 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
1531 S=FB/FA
1532 IF(A.EQ.C) THEN
1533 P=2.*XM*S
1534 Q=1.-S
1535 ELSE
1536 Q=FA/FC
1537 R=FB/FC
1538 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
1539 Q=(Q-1.)*(R-1.)*(S-1.)
1540 ENDIF
1541 IF(P.GT.0.) Q=-Q
1542 P=ABS(P)
1543 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
1544 E=D
1545 D=P/Q
1546 ELSE
1547 D=XM
1548 E=D
1549 ENDIF
1550 ELSE
1551 D=XM
1552 E=D
1553 ENDIF
1554 A=B
1555 FA=FB
1556 IF(ABS(D) .GT. TOL1) THEN
1557 B=B+D
1558 ELSE
1559 B=B+SIGN(TOL1,XM)
1560 ENDIF
1561 FB=FUNC(B)
156211 CONTINUE
1563 WRITE (*, *) 'CtLhQZBRNT exceeding maximum iterations.'
1564 IRT = 2
1565 CtLhQZBRNT=B
1566 RETURN
1567 END
1568 SUBROUTINE CtLhRATINT(XA,YA,N,X,Y,DY)
1569 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1570 PARAMETER (NMAX=10,TINY=1.E-25)
1571 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1572 NS=1
1573 HH=ABS(X-XA(1))
1574 DO 11 I=1,N
1575 H=ABS(X-XA(I))
1576 IF (H.EQ.0.)THEN
1577 Y=YA(I)
1578 DY=0.0
1579 RETURN
1580 ELSE IF (H.LT.HH) THEN
1581 NS=I
1582 HH=H
1583 ENDIF
1584 C(I)=YA(I)
1585 D(I)=YA(I)+TINY
158611 CONTINUE
1587 Y=YA(NS)
1588 NS=NS-1
1589 DO 13 M=1,N-1
1590 DO 12 I=1,N-M
1591 W=C(I+1)-D(I)
1592 H=XA(I+M)-X
1593 T=(XA(I)-X)*D(I)/H
1594 DD=T-C(I+1)
1595 DD=W/DD
1596 D(I)=C(I+1)*DD
1597 C(I)=T*DD
159812 CONTINUE
1599 IF (2*NS.LT.N-M)THEN
1600 DY=C(NS+1)
1601 ELSE
1602 DY=D(NS)
1603 NS=NS-1
1604 ENDIF
1605 Y=Y+DY
160613 CONTINUE
1607 RETURN
1608 END
1609 FUNCTION CtLhRTALF (EFLLN)
1610 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1611 include 'parmsetup.inc'
1612 PARAMETER (PI = 3.141592653589793d0)
1613 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
1614 EFMULM = EXP (EFLLN)
1615 TEM1 = PI / ALFST
1616 TEM2 = 1. / CtLhALPQCD (JORD, NEFF, EFMULM, I)
1617 CtLhRTALF = TEM1 - TEM2
1618 END
1619 Subroutine CtLhbldat1
1620 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1621 include 'parmsetup.inc'
1622 LOGICAL LSTX
1623 PARAMETER (MXX = 105, MXQ = 25, MxF = 6)
1624 PARAMETER (MxPN = MxF * 2 + 2)
1625 PARAMETER (MxQX= MXQ * MXX, MxPQX = MxQX * MxPN)
1626 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1627 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1628 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1629 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
1630 PARAMETER (NF0 = 4, Nshp = 8,NEX = Nshp+2)
1631 XMIN = .999999D-4
1632 XCR = 1.5
1633 JT = 1
1634 Return
1635 END
1636 SUBROUTINE CtLhSETL1 (NEF, VLAM)
1637 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1638 LOGICAL SET
1639 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1640 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1641 COMMON / LhCtCOMQMS / VALQMS(9)
1642 IF (NEF .LT. 0 .OR. NEF .GT. NF) THEN
1643 WRITE(*,*)'NEF out of range in CtLhSETL1: NEF NF =',NEF,NF
1644 STOP
1645 ENDIF
1646 AMHAT(0) = 0.
1647 DO 5 N = 1, NF
1648 AMHAT(N) = VALQMS(N)
1649 5 CONTINUE
1650 ALAM(NEF) = VLAM
1651 DO 10 N = NEF, 1, -1
1652 CALL CtLhTRNLAM(NORDER, N, -1, IR1)
1653 10 CONTINUE
1654 DO 20 N = NEF, NF-1
1655 CALL CtLhTRNLAM(NORDER, N, 1, IR1)
1656 20 CONTINUE
1657 DO 30, N = NF, 1, -1
1658 IF ((ALAM(N) .GE. 0.7 * AMHAT(N))
1659 > .OR. (ALAM(N-1) .GE. 0.7 * AMHAT(N)))THEN
1660 NHQ = NF - N
1661 GOTO 40
1662 ENDIF
1663 30 CONTINUE
1664 NHQ = NF
1665 40 CONTINUE
1666 DO 50, N = NF-NHQ, 1, -1
1667 AMHAT(N) = 0
1668 ALAM(N-1) = ALAM(N)
1669 50 CONTINUE
1670 AMN = ALAM(NF)
1671 DO 60, N = 0, NF-1
1672 IF (ALAM(N) .GT. AMN) AMN = ALAM(N)
1673 60 CONTINUE
1674 AMN = AMN * 1.0001
1675 AL = ALAM(NF)
1676 SET = .TRUE.
1677 RETURN
1678 END
1679 SUBROUTINE CtLhSETLAM (NEF, WLAM, IRDR)
1680 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1681 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1682 LOGICAL SET
1683 IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN
1684 WRITE(*,*)'NEF out of range in CtLhSETLAM: NEF NF=',NEF,NF
1685 STOP
1686 ENDIF
1687 VLAM = WLAM
1688 IF (IRDR .NE. NORDER) then
1689 PRINT *,'fatal error: wanted cnvl1'
1690 stop
1691 ENDIF
1692 CALL CtLhSETL1 (NEF, VLAM)
1693 END
1694 Subroutine CtLhbldat2
1695 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1696 COMMON / LhCtCOMQMS / VALQMS(9)
1697 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1698 LOGICAL SET
1699 AL = .226d0
1700 NF = 5
1701 NORDER = 2
1702 SET = .FALSE.
1703 VALQMS(1) = 0.
1704 VALQMS(2) = 0.
1705 VALQMS(3) = 0.2d0
1706 VALQMS(4) = 1.3d0
1707 VALQMS(5) = 4.5d0
1708 VALQMS(6) = 174.d0
1709 VALQMS(7) = 0.
1710 VALQMS(8) = 0.
1711 VALQMS(9) = 0.
1712 Return
1713 END
1714 FUNCTION CtLhSMPNOL (NX, DX, FN, ERR)
1715 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1716 DIMENSION FN(NX)
1717 MS = MOD(NX, 2)
1718 IF (NX .LE. 1 .OR. NX .GT. 1000) THEN
1719 PRINT *, 'NX =', NX, ' OUT OF RANGE IN CtLhSMPNOL!'
1720 STOP
1721 ELSEIF (NX .EQ. 2) THEN
1722 TEM = DX * FN(2)
1723 ELSEIF (NX .EQ. 3) THEN
1724 TEM = DX * FN(2) * 2.
1725 ELSE
1726 IF (MS .EQ. 0) THEN
1727 TEM = DX * (23.* FN(2) - 16.* FN(3) + 5.* FN(4)) / 12.
1728 TMP = DX * (3.* FN(2) - FN(3)) / 2.
1729 ERR = ABS(TEM - TMP)
1730 TEM = TEM + CtLhSMPSNA (NX-1, DX, FN(2), ER1)
1731 ERR = ABS(ER1) + ERR
1732 ELSE
1733 TEM = DX * (8.* FN(2) - 4.* FN(3) + 8.* FN(4)) / 3.
1734 TMP = DX * (3.* FN(2) + 2.* FN(3) + 3.* FN(4)) / 2.
1735 ERR = ABS(TEM - TMP)
1736 TEM = TEM + CtLhSMPSNA (NX-4, DX, FN(5), ER1)
1737 ERR = ABS(ER1) + ERR
1738 ENDIF
1739 ENDIF
1740 CtLhSMPNOL = TEM
1741 RETURN
1742 END
1743 FUNCTION CtLhSMPSNA (NX, DX, F, ERR)
1744 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1745 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1746 PARAMETER (MAXX = 1000)
1747 DIMENSION F(NX)
1748 DATA IW1, IW2, TINY / 2*0, 1.E-35 /
1749 IF (DX .LE. 0.) THEN
1750 CALL CtLhWARNR(IW2,'DX cannot be < 0. in CtLhSMPSNA', 'DX',
1751 > DX, D0, D1, 0)
1752 CtLhSMPSNA = 0.
1753 RETURN
1754 ENDIF
1755 IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN
1756 CALL CtLhWARNI(IW1, 'NX out of range in CtLhSMPSNA', 'NX', NX,
1757 > 1, MAXX, 1)
1758 SIMP = 0.
1759 ELSEIF (NX .EQ. 1) THEN
1760 SIMP = 0.
1761 ELSEIF (NX .EQ. 2) THEN
1762 SIMP = (F(1) + F(2)) / 2.
1763 ERRD = (F(1) - F(2)) / 2.
1764 ELSE
1765 MS = MOD(NX, 2)
1766 IF (MS .EQ. 0) THEN
1767 ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24.
1768 NZ = NX - 1
1769 ELSE
1770 ADD = 0.
1771 NZ = NX
1772 ENDIF
1773 IF (NZ .EQ. 3) THEN
1774 SIMP = (F(1) + 4.* F(2) + F(3)) / 3.
1775 TRPZ = (F(1) + 2.* F(2) + F(3)) / 2.
1776 ELSE
1777 SE = F(2)
1778 SO = 0
1779 NM1 = NZ - 1
1780 DO 60 I = 4, NM1, 2
1781 IM1 = I - 1
1782 SE = SE + F(I)
1783 SO = SO + F(IM1)
1784 60 CONTINUE
1785 SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3.
1786 TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2.
1787 ENDIF
1788 ERRD = TRPZ - SIMP
1789 SIMP = SIMP + ADD
1790 ENDIF
1791 CtLhSMPSNA = SIMP * DX
1792 IF (ABS(SIMP) .GT. TINY) THEN
1793 ERR = ERRD / SIMP
1794 ELSE
1795 ERR = 0.
1796 ENDIF
1797 RETURN
1798 END
1799 SUBROUTINE CtLhSNEVL(IKNL,NX,NT,JT,DT,TIN,NEFF,UI,GI,US,GS)
1800 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1801 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1802 PARAMETER (MXQX= MXQ * MXX)
1803 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1804 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1805 DIMENSION UI(NX), US(0:NX, 0:NT)
1806 DIMENSION GI(NX), GS(0:NX, 0:NT)
1807 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
1808 DIMENSION Z0(MXX), Z1(MXX), ZP(MXX), G0(MXX), G1(MXX), GP(MXX)
1809 DATA D0 / 0.0 /
1810 JTT = 2 * JT
1811 DDT = DT / JTT
1812 IF (NX .GT. MXX) THEN
1813 WRITE (*,*) 'Nx =', NX, ' too many pts in CtLhSNEVL'
1814 STOP 'Program stopped in CtLhSNEVL'
1815 EndIf
1816 TMD = TIN + DT * NT / 2.
1817 AMU = EXP(TMD)
1818 TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
1819 TLAM = TMD - TEM
1820 DO 9 IX = 1, NX
1821 US (IX, 0) = UI(IX)
1822 GS (IX, 0) = GI(IX)
1823 9 CONTINUE
1824 US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3)
1825 GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3)
1826 TT = TIN
1827 DO 10 IZ = 1, NX
1828 Y0(IZ) = UI(IZ)
1829 Z0(IZ) = GI(IZ)
1830 10 CONTINUE
1831 DO 20 IS = 1, NT
1832 DO 202 JS = 1, JTT
1833 IRND = (IS-1) * JTT + JS
1834 IF (IRND .EQ. 1) THEN
1835 CALL CtLhSNRHS (TT, NEFF, Y0,Z0, F0,G0)
1836 DO 250 IZ = 1, NX
1837 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
1838 Z0(IZ) = Z0(IZ) + DDT * G0(IZ)
1839 250 CONTINUE
1840 TT = TT + DDT
1841 CALL CtLhSNRHS (TT, NEFF, Y0, Z0, F1, G1)
1842 DO 251 IZ = 1, NX
1843 Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
1844 Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0
1845 251 CONTINUE
1846 Else
1847 CALL CtLhSNRHS (TT, NEFF, Y1, Z1, F1, G1)
1848 DO 252 IZ = 1, NX
1849 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
1850 ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0
1851 252 CONTINUE
1852 TT = TT + DDT
1853 CALL CtLhSNRHS (TT, NEFF, YP, ZP, FP, GP)
1854 DO 253 IZ = 1, NX
1855 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
1856 Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0
1857 F0(IZ) = F1(IZ)
1858 G0(IZ) = G1(IZ)
1859 253 CONTINUE
1860 EndIf
1861 202 CONTINUE
1862 DO 260 IX = 1, NX
1863 IF (IKNL .GT. 0) THEN
1864 US (IX, IS) = MAX(Y1(IX), D0)
1865 GS (IX, IS) = MAX(Z1(IX), D0)
1866 Else
1867 US (IX, IS) = Y1(IX)
1868 GS (IX, IS) = Z1(IX)
1869 EndIf
1870 260 CONTINUE
1871 US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
1872 GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3)
1873 20 CONTINUE
1874 RETURN
1875 END
1876 SUBROUTINE CtLhSNRHS (TT, NEFF, FI, GI, FO, GO)
1877 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1878 LOGICAL LSTX
1879 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1880 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1881 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1882 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1883 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1884 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1885 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1886 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1887 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1888 DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX)
1889 DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), WM(MXX)
1890 DIMENSION R0(MXX), R1(MXX), R2(MXX), RH(MXX), RM(MXX)
1891 S = EXP(TT)
1892 Q = AL * EXP (S)
1893 CPL = CtLhALPI(Q)
1894 CPL2= CPL ** 2 / 2. * S
1895 CPL = CPL * S
1896 CALL CtLhINTEGR (NX,-1, FI, WM, IR1)
1897 CALL CtLhINTEGR (NX, 0, FI, W0, IR2)
1898 CALL CtLhINTEGR (NX, 1, FI, W1, IR3)
1899 CALL CtLhINTEGR (NX,-1, GI, RM, IR4)
1900 CALL CtLhINTEGR (NX, 0, GI, R0, IR5)
1901 CALL CtLhINTEGR (NX, 1, GI, R1, IR6)
1902 CALL CtLhINTEGR (NX, 2, GI, R2, IR7)
1903 CALL CtLhHINTEG (NX, FI, WH)
1904 CALL CtLhHINTEG (NX, GI, RH)
1905 IF (IKNL .GT. 0) THEN
1906 DO 230 IZ = 1, NX
1907 FO(IZ) = ( 2D0 * FI(IZ)
1908 > + 4D0 / 3D0 * ( 2D0 * WH(IZ) - W0(IZ) - W1(IZ) ))
1909 > + NEFF * ( R0(IZ) - 2D0 * R1(IZ) + 2D0 * R2(IZ) )
1910 FO(IZ) = FO(IZ) * CPL
1911 GO(IZ) = 4D0 / 3D0 * ( 2D0 * WM(IZ) - 2D0 * W0(IZ) + W1(IZ) )
1912 > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ)
1913 > + 6D0 * (RH(IZ) + RM(IZ) - 2D0 * R0(IZ) + R1(IZ) - R2(IZ))
1914 GO(IZ) = GO(IZ) * CPL
1915 230 CONTINUE
1916 Else
1917 DO 240 IZ = 1, NX
1918 FO(IZ) = NEFF * (-R0(IZ) + 2.* R1(IZ) )
1919 > + 2.* FI(IZ) + 4./ 3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ) )
1920 FO(IZ) = FO(IZ) * CPL
1921 GO(IZ) = 4./ 3.* ( 2.* W0(IZ) - W1(IZ) )
1922 >+ (33.- 2.* NEFF) / 6.* GI(IZ) + 6.*(RH(IZ) + R0(IZ) - 2.* R1(IZ))
1923 GO(IZ) = GO(IZ) * CPL
1924 240 CONTINUE
1925 EndIf
1926 IF (IKNL .EQ. 2) THEN
1927 DZ = 1./(NX - 1)
1928 DO 21 I = 1, NX-1
1929 X = XV(I)
1930 NP = NX - I + 1
1931 IS = NP
1932 g2(1)=0d0
1933 g3(1)=0d0
1934 DO 31 KZ = 2, NP
1935 IY = I + KZ - 1
1936 IT = NX - IY + 1
1937 XY = ZZ (IS, IT)
1938 G1(KZ) = FFG(I, IY) * (FI(IY) - XY**2 *FI(I))
1939 G4(KZ) = GGF(I, IY) * (GI(IY) - XY**2 *GI(I))
1940 G2(KZ) = FFG(IS,IT) * (GI(IY) - xy*GI(I)) !FG
1941 G3(KZ) = GGF(IS,IT) * (FI(IY) - XY*FI(I)) !GF (usual notations)
1942 31 CONTINUE
1943 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
1944 TEM2 = CtLhSMPSNA (NP, DZ, G2, ERR)
1945 TEM3 = CtLhSMPSNA (NP, DZ, G3, ERR)
1946 TEM4 = CtLhSMPNOL (NP, DZ, G4, ERR)
1947 TEM1 = TEM1 - FI(I) * (AFF2(I) + ZGF2)
1948 TEM4 = TEM4 - GI(I) * (AGG2(I) + ZFG2)
1949 tem2 = tem2 + GI(I)*AFG2(I)
1950 tem3= tem3 + FI(I)*AGF2(I)
1951 TMF = TEM1 + TEM2
1952 TMG = TEM3 + TEM4
1953 FO(I) = FO(I) + TMF * CPL2
1954 GO(I) = GO(I) + TMG * CPL2
1955 21 CONTINUE
1956 EndIf
1957 RETURN
1958 END
1959 FUNCTION CtLhSPENC2 (X)
1960 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1961 EXTERNAL CtLhSPN2IN
1962 COMMON / LhCtSPENCC / XX
1963 DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 /
1964 XX = X
1965 TEM = CtLhGausInt(CtLhSPN2IN, XX, U1, AERR, RERR, ERR, IRT)
1966 CtLhSPENC2 = TEM + LOG (XX) ** 2 / 2.
1967 RETURN
1968 END
1969 FUNCTION CtLhSPN2IN (ZZ)
1970 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1971 COMMON / LhCtSPENCC / X
1972 Z = ZZ
1973 TEM = LOG (1.+ X - Z) / Z
1974 CtLhSPN2IN = TEM
1975 RETURN
1976 END
1977 SUBROUTINE CtLhSTUPKL (NFL)
1978 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1979 LOGICAL LSTX
1980 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1981 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1982 PARAMETER (MX = 3)
1983 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1984 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1985 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1986 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1987 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1988 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1989 > PNSP(0:MXX), PNSM(0:MXX)
1990 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1991 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1992 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1993 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1994 EXTERNAL CtLhPFF1, CtLhRGG1, CtLhRFF2, CtLhRGG2
1995 EXTERNAL CtLhFNSP, CtLhFNSM
1996 dimension aff1(mxx),agg1(mxx)
1997 PARAMETER (PI = 3.141592653589793d0, PI2 = PI**2)
1998 DATA CF, CG, TR / 1.33333333333333d0, 3.0, 0.5 /
1999 data zeta3/1.20205690315959d0/ ! zeta(3.0)
2000 SAVE
2001 DATA AERR, RERR / 0.0, 0.02 /
2002 NNX = NX
2003 DZ = 1./ (NX - 1)
2004 DO 5 I0 = 1, MX
2005 XL(I0) = XV(I0)
2006 5 CONTINUE
2007 DO 10 I = 1, NX-1
2008 XZ = XV(I)
2009 CALL CtLhKERNEL (XZ, FF1(I), GF1(I), FG1(I), GG1(I), PNSP(I),
2010 > PNSM(I), FF2(I), GF2(I), FG2(I), GG2(I), NFL, IRT)
2011 10 CONTINUE
2012 FF1(0) = FF1(1) * 3. - FF1(2) * 3. + FF1(3)
2013 FG1(0) = FG1(1) * 3. - FG1(2) * 3. + FG1(3)
2014 GF1(0) = GF1(1) * 3. - GF1(2) * 3. + GF1(3)
2015 GG1(0) = GG1(1) * 3. - GG1(2) * 3. + GG1(3)
2016 PNSP(0) = PNSP(1) * 3. - PNSP(2) * 3. + PNSP(3)
2017 PNSM(0) = PNSM(1) * 3. - PNSM(2) * 3. + PNSM(3)
2018 FF2(0) = FF2(1) * 3. - FF2(2) * 3. + FF2(3)
2019 FG2(0) = FG2(1) * 3. - FG2(2) * 3. + FG2(3)
2020 GF2(0) = GF2(1) * 3. - GF2(2) * 3. + GF2(3)
2021 GG2(0) = GG2(1) * 3. - GG2(2) * 3. + GG2(3)
2022 FF1(NX) = FF1(NX-1) * 3. - FF1(NX-2) * 3. + FF1(NX-3)
2023 FG1(NX) = FG1(NX-1) * 3. - FG1(NX-2) * 3. + FG1(NX-3)
2024 GF1(NX) = GF1(NX-1) * 3. - GF1(NX-2) * 3. + GF1(NX-3)
2025 GG1(NX) = GG1(NX-1) * 3. - GG1(NX-2) * 3. + GG1(NX-3)
2026 PNSM(NX) = PNSM(NX-1) * 3. - PNSM(NX-2) * 3. + PNSM(NX-3)
2027 PNSP(NX) = PNSP(NX-1) * 3. - PNSP(NX-2) * 3. + PNSP(NX-3)
2028 FF2(NX) = FF2(NX-1) * 3. - FF2(NX-2) * 3. + FF2(NX-3)
2029 FG2(NX) = FG2(NX-1) * 3. - FG2(NX-2) * 3. + FG2(NX-3)
2030 GF2(NX) = GF2(NX-1) * 3. - GF2(NX-2) * 3. + GF2(NX-3)
2031 GG2(NX) = GG2(NX-1) * 3. - GG2(NX-2) * 3. + GG2(NX-3)
2032 RER = RERR * 4.
2033 AFF1(1) = CtLhGausInt(CtLhPFF1,D0,XV(1),AERR,RERR,ER1,IRT)
2034 DGG1 = NFL / 3.
2035 TMPG = CtLhGausInt(CtLhRGG1,D0,XV(1),AERR,RERR,ER3,IRT)
2036 AGG1(1) = TMPG + DGG1
2037 ANSM(1) = CtLhGausInt(CtLhFNSM,D0,XV(1),AERR,RER,ER2,IRT)
2038 ANSP(1) = CtLhGausInt(CtLhFNSP,D0,XV(1),AERR,RER,ER2,IRT)
2039 AER = AFF1(1) * RER
2040 AFF2(1) = CtLhGausInt(CtLhRFF2, D0, XV(1), AER, RER, ER2, IRT)
2041 AER = AGG1(1) * RER
2042 AGG2(1) = CtLhGausInt(CtLhRGG2, D0, XV(1), AER, RER, ER4, IRT)
2043 DO 20 I2 = 2, NX-1
2044 TEM =CtLhGausInt(CtLhPFF1,XV(I2-1),XV(I2),AERR,RERR,ER1,IRT)
2045 AFF1(I2) = TEM + AFF1(I2-1)
2046 AER = ABS(TEM * RER)
2047 AFF2(I2)=CtLhGausInt(CtLhRFF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT)
2048 > +AFF2(I2-1)
2049 TEM = CtLhGausInt(CtLhRGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT)
2050 TMPG = TMPG + TEM
2051 AGG1(I2) = TMPG + DGG1
2052 AER = ABS(TEM * RER)
2053 AGG2(I2)=CtLhGausInt(CtLhRGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)
2054 > +AGG2(I2-1)
2055 ANSP(I2)=CtLhGausInt(CtLhFNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2056 > +ANSP(I2-1)
2057 ANSM(I2)=CtLhGausInt(CtLhFNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2058 > +ANSM(I2-1)
2059 20 CONTINUE
2060 ANSP(NX)=CtLhGausInt(CtLhFNSP,XV(NX-1),D1,AERR,RER,ERR,
2061 > IRT) + ANSP(NX-1)
2062 ANSM(NX)=CtLhGausInt(CtLhFNSM,XV(NX-1),D1,AERR,RER,ERR,
2063 > IRT) + ANSM(NX-1)
2064 TRNF = TR * NFL
2065 do i2=1,nx-1 !loop over x
2066 x=xv(i2)
2067 XI = 1./ X !auxiliary definitions
2068 X2 = X ** 2
2069 X3= x**3
2070 XLN = DLOG (X)
2071 XLN2 = XLN ** 2
2072 XLN1M = DLOG (1.- X)
2073 xLi2m=CtLhxLi(2,-x)
2074 xLi2=CtLhxLi(2,x)
2075 xLi3=CtLhxLi(3,x)
2076 xLi31m=CtLhxLi(3,1d0-x)
2077 xLi32=CtLhxLi(3,x2)
2078 xln1m2=xln1m*xln1m
2079 xln1p=dlog(1d0+x)
2080 x1m=1d0-x
2081 x1p=1d0+x
2082 x3m=3d0-x
2083 x3p=3d0+x
2084 wgfcft=
2085 > (9 + 4*Pi2 - 22*x + 13*x2 + 6*(3 - 4*x + x2)*xln1m +
2086 > 40*xln - 24*xLi2)/9.
2087 wgfcf2=
2088 > (6*(2*(-9 + Pi2) + 3*x*(5 + x)) +4*(3 +2*Pi2+3*x*(-3 + 2*x))*
2089 > xln1m + 6*x3m*x1m*xln1m2 - 6*(x*(8 + 3*x) + 4*xln1m2)*
2090 > xln - 3*(-4 + x)*x*xln2)/12 - 2*(3 + 2*xln1m)*xLi2 - 4*xLi31m
2091 wgfcfg=
2092 > (3637-186*Pi2-x*(3198+72*Pi2+x*(231 + 208*x)))/108.- xln +
2093 > (3*xln1m*(-33 - 4*Pi2 + (50 - 17*x)*x - 3*x3m*x1m*xln1m) +
2094 > 2*(x*(198 + x*(27+8*x))+9*xln1m*(3 - 4*x + x2 + 2*xln1m))*
2095 > xln - 9*x*(4 + x)*xln2)/18- x1p*x3p*xln*xln1p-
2096 > (x1p*x3p - 4*xln)*xLi2m + (31d0/3d0 +4*xln1m- 4*xln)*xLi2 +
2097 > 4*xLi31m + 12*xLi3 - 2*xLi32 - 10*zeta3
2098 wfgcft=
2099 > (18 - 81*x + 6*Pi2*x + 123*x2 - 6*Pi2*x2 - 60*x3 +
2100 > 4*Pi2*x3 - 6*(-2 + 3*x - 3*x2 + 2*x3)*xln1m2 -33*x*xln +
2101 > 15*x2*xln - 24*x3*xln - 9*x*xln2 + 9*x2*xln2 -
2102 > 12*x3*xln2 - 12*x1m*xln1m*(-1 + 2*x2 + 2*xln - x*xln +
2103 > 2*x2*xln) - 24*xLi2)/9.
2104 wfgcgt=
2105 > (2*(-67 + 2*Pi2 + x*(64 + x*(-91 + 3*Pi2 + 94*x)) +
2106 > x1m*(7+x*(-5+16*x))*xln1m -3*x1m*(2+ x*(-1+2*x))*xln1m2 -
2107 > 20*xln - 3*x*xln*(13 + 16*x*x1p - 3*x1p*xln) +
2108 > 6*x1p*(2+x+2*x2)*xln*xln1p+6*x1p*(2+x+2*x2)*xLi2m))/9.
2109 AGF2(I2) = CF*TRNF*WGFCFT+CF**2* WGFCF2+CF*CG*WGFCFG
2110 AFG2(I2) = CF*TRNF*WFGCFT +CG*TRNF*WFGCGT
2111 enddo !i2
2112 AGF2(nx)=0d0
2113 AFG2(nx)=0d0
2114 ZGF2=-28./27.*Cf**2+94./27.*Cf*Cg -52./27.*Cf*TrNf
2115 ZFG2= 37./27.*Cf*TrNf + 35./54.*Cg*TrNf
2116 ZQQB=1.43862321154902*(Cf**2-0.5*Cf*Cg)
2117 DO 21 IX = 1, NX-1
2118 X = XV(IX)
2119 NP = NX - IX + 1
2120 IS = NP
2121 XG2 = (LOG(1./(1.-X)) + 1.) ** 2
2122 FFG (IS, IS) = FG2(NX) * DXTZ(I) * XG2
2123 GGF (IS, IS) = GF2(NX) * DXTZ(I) * XG2
2124 PNS (IS, IS) =PNSM(NX) * DXTZ(I)
2125 DO 31 KZ = 2, NP
2126 IY = IX + KZ - 1
2127 IT = NX - IY + 1
2128 XY = X / XV(IY)
2129 XM1 = 1.- XY
2130 XG2 = (LOG(1./XM1) + 1.) ** 2
2131 Z = ZZ (IX, IY)
2132 TZ = (Z + DZ) / DZ
2133 IZ = TZ
2134 IZ = MAX (IZ, 0)
2135 IZ = MIN (IZ, NX-1)
2136 DT = TZ - IZ
2137 TEM = (FF2(IZ) * (1.- DT) + FF2(IZ+1) * DT) / XM1 / XY
2138 FFG (IX, IY) = TEM * DXTZ(IY)
2139 TEM = (FG2(IZ) * (1.- DT) + FG2(IZ+1) * DT) * XG2 / XY
2140 FFG (IS, IT) = TEM * DXTZ(IY)
2141 TEM = (GF2(IZ) * (1.- DT) + GF2(IZ+1) * DT) * XG2 / XY
2142 GGF (IS, IT) = TEM * DXTZ(IY)
2143 TEM = (GG2(IZ) * (1.- DT) + GG2(IZ+1) * DT) / XM1 / XY
2144 GGF (IX, IY) = TEM * DXTZ(IY)
2145 TEM = (PNSP(IZ) * (1.- DT) + PNSP(IZ+1) * DT) / XM1
2146 PNS (IX, IY) = TEM * DXTZ(IY)
2147 TEM = (PNSM(IZ) * (1.- DT) + PNSM(IZ+1) * DT) / XM1
2148 PNS (IS, IT) = TEM * DXTZ(IY)
2149 31 CONTINUE
2150 21 CONTINUE
2151 RETURN
2152 END
2153 SUBROUTINE CtLhTRNLAM (IRDR, NF, IACT, IRT)
2154 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2155 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
2156 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2157 EXTERNAL CtLhZBRLAM
2158 DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 /
2159 DATA IR1, SML / 0, 1.E-5 /
2160 IRT = 0
2161 N = NF
2162 JRDR = IRDR
2163 JACT = IACT
2164 VLAM = ALAM(N)
2165 IF (JACT .GT. 0) THEN
2166 N1 = N + 1
2167 THMS = AMHAT(N1)
2168 ALM = LOG (THMS/VLAM)
2169 BLM = BLM0
2170 ELSE
2171 N1 = N -1
2172 THMS = AMHAT(N)
2173 ALM = ALM0
2174 THMS = MAX (THMS, SML)
2175 BLM = LOG (THMS/VLAM)
2176 ENDIF
2177 IF (VLAM .GE. 0.7 * THMS) THEN
2178 IF (JACT . EQ. 1) THEN
2179 AMHAT(N1) = 0
2180 ELSE
2181 AMHAT(N) = 0
2182 ENDIF
2183 IRT = 4
2184 ALAM(N1) = VLAM
2185 RETURN
2186 ENDIF
2187 IF (ALM .GE. BLM) THEN
2188 WRITE (*, *) 'CtLhTRNLAM has ALM >= BLM: ', ALM, BLM
2189 WRITE (*, *) 'I do not know how to continue'
2190 STOP
2191 ENDIF
2192 VMULM = THMS/VLAM
2193 ERR = RERR * LOG (VMULM)
2194 WLLN = CtLhQZBRNT (CtLhZBRLAM, ALM, BLM, ERR, IR1)
2195 ALAM(N1) = THMS / EXP (WLLN)
2196 IF (IR1 .NE. 0) THEN
2197 WRITE (*, *) 'CtLhQZBRNT failed in CtLhTRNLAM; ',
2198 > 'NF, VLAM =', NF, VLAM
2199 WRITE (*, *) 'I do not know how to continue'
2200 STOP
2201 ENDIF
2202 RETURN
2203 END
2204 SUBROUTINE CtLhUPC (A, La, UpA)
2205 CHARACTER A*(*), UpA*(*), C*(1)
2206 INTEGER I, La, Ld
2207 La = Len(A)
2208 Lb = Len(UpA)
2209 If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!'
2210 Ld = ICHAR('A')-ICHAR('a')
2211 DO 1 I = 1, Lb
2212 If (I .Le. La) Then
2213 c = A(I:I)
2214 IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN
2215 UpA (I:I) = CHAR(Ichar(c) + ld)
2216 Else
2217 UpA (I:I) = C
2218 ENDIF
2219 Else
2220 UpA (I:I) = ' '
2221 Endif
2222 1 CONTINUE
2223
2224 RETURN
2225 END
2226 SUBROUTINE CtLhWARNI (IWRN, MSG, NMVAR, IVAB,
2227 > IMIN, IMAX, IACT)
2228 CHARACTER*(*) MSG, NMVAR
2229 Save Iw
2230 Data Nmax / 100 /
2231 IW = IWRN
2232 IV = IVAB
2233
2234 IF (IW .EQ. 0) THEN
2235 PRINT '(1X,A/1X, 2A,I10 /A,I4)', MSG, NMVAR, ' = ', IV
2236 IF (IACT .EQ. 1) THEN
2237 PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX
2238 ENDIF
2239 ENDIF
2240 If (Iw .LT. Nmax) Then
2241 PRINT '(1X,A/1X,I10,A, I10)', MSG, NMVAR, ' = ', IV
2242 Elseif (Iw .Eq. Nmax) Then
2243 Print '(/A/)', 'CtLhWARNI Severe Warning: Too many errors'
2244 Endif
2245 IWRN = IW + 1
2246 RETURN
2247 END
2248 SUBROUTINE CtLhWARNR (IWRN, MSG, NMVAR, VARIAB,
2249 > VMIN, VMAX, IACT)
2250 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2251 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2252 CHARACTER*(*) MSG, NMVAR
2253 Save Iw
2254 Data Nmax / 100 /
2255 IW = IWRN
2256 VR = VARIAB
2257 IF (IW .EQ. 0) THEN
2258 PRINT '(1X, A/1X,2A,1PD16.7/A,I4)', MSG, NMVAR, ' = ', VR
2259 IF (IACT .EQ. 1) THEN
2260 PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX
2261 ENDIF
2262 ENDIF
2263 If (Iw .LT. Nmax) Then
2264 PRINT '(I5, 2A/1X,2A,I10,1PD16.7)', IW, ' ', MSG,
2265 > NMVAR, ' = ', VR
2266 Elseif (Iw .Eq. Nmax) Then
2267 Print '(/A/)', 'CtLhWARNR Severe Warning: Too many errors'
2268 Endif
2269 IWRN = IW + 1
2270 RETURN
2271 END
2272 SUBROUTINE CtLhXARRAY
2273 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2274 LOGICAL LSTX
2275 PARAMETER (D0 = 0.0, D10=10.0)
2276 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
2277 PARAMETER (MXPN = MXF * 2 + 2)
2278 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
2279 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
2280 Character Msg*80
2281 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
2282 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
2283 COMMON / LhCtHINTEC / GH(NDG, MXX)
2284 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2285 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
2286 DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG)
2287 DATA F12, F22, F32 / 1D0, 1D0, 1D0 /
2288 DATA (G1(I,NDH), G2(I,1), I=1,NDG) / 0.0,0.0,0.0,0.0,0.0,0.0 /
2289 DATA PUNY / 1D-30 /
2290 XV(0) = 0D0
2291 DZ = 1D0 / (NX-1)
2292 DO 10 I = 1, NX - 1
2293 Z = DZ * (I-1)
2294 ZV(I) = Z
2295 X = CtLhXFRMZ (Z)
2296 DXTZ(I) = CtLhDXDZ(Z) / X
2297 XV (I) = X
2298 XA(I, 1) = X
2299 XA(I, 0) = LOG (X)
2300 DO 20 L = L1, L2
2301 IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L
2302 20 CONTINUE
2303 10 CONTINUE
2304 XV(1) = Xmin
2305 XV(NX) = 1D0
2306 ZV(Nx) = 1D0
2307 DXTZ(NX) = CtLhDXDZ(1.D0)
2308 DO 21 L = L1, L2
2309 XA (NX, L) = 1D0
2310 21 CONTINUE
2311 XA (NX, 0) = 0D0
2312 DO 11 I = 1, NX-1
2313 ELY(I) = LOG(1D0 - XV(I))
2314 11 CONTINUE
2315 ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3)
2316 DO 17 IX = 1, NX
2317 ZZ (IX, IX) = 1.
2318 DO 17 IY = IX+1, NX
2319 XY = XV(IX) / XV(IY)
2320 ZZ (IX, IY) = CtLhZFRMX (XY)
2321 ZZ (NX-IX+1, NX-IY+1) = XY
2322 17 CONTINUE
2323 DO 30 I = 1, NX-1
2324 IF (I .NE. NX-1) THEN
2325 F11 = 1D0/XV(I)
2326 F21 = 1D0/XV(I+1)
2327 F31 = 1D0/XV(I+2)
2328 F13 = XV(I)
2329 F23 = XV(I+1)
2330 F33 = XV(I+2)
2331 DET = F11*F22*F33 + F21*F32*F13 + F31*F12*F23
2332 > - F31*F22*F13 - F21*F12*F33 - F11*F32*F23
2333 IF (ABS(DET) .LT. PUNY) THEN
2334 Msg='Determinant close to zero; will be arbitrarily set to:'
2335 CALL CtLhWARNR(IWRN, Msg, 'DET', PUNY, D0, D0, 0)
2336 DET = PUNY
2337 EndIf
2338 G2(1,2) = (F22*F33 - F23*F32) / DET
2339 G2(1,3) = (F32*F13 - F33*F12) / DET
2340 G2(1,4) = (F12*F23 - F13*F22) / DET
2341 G2(2,2) = (F23*F31 - F21*F33) / DET
2342 G2(2,3) = (F33*F11 - F31*F13) / DET
2343 G2(2,4) = (F13*F21 - F11*F23) / DET
2344 G2(3,2) = (F21*F32 - F22*F31) / DET
2345 G2(3,3) = (F31*F12 - F32*F11) / DET
2346 G2(3,4) = (F11*F22 - F12*F21) / DET
2347 B2 = LOG (XV(I+2)/XV(I))
2348 B3 = XV(I) * (B2 - 1.) + XV(I+2)
2349 GH (1,I) = B2 * G2 (2,2) + B3 * G2 (3,2)
2350 GH (2,I) = B2 * G2 (2,3) + B3 * G2 (3,3)
2351 GH (3,I) = B2 * G2 (2,4) + B3 * G2 (3,4)
2352 EndIf
2353 DO 51 J = 1, NDH
2354 DO 52 L = 1, NDG
2355 IF (I .EQ. 1) THEN
2356 GB(L,J,I) = G2(L,J)
2357 ElseIF (I .EQ. NX-1) THEN
2358 GB(L,J,I) = G1(L,J)
2359 Else
2360 GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0
2361 EndIf
2362 52 CONTINUE
2363 51 CONTINUE
2364 DO 35 MM = M1, M2
2365 DO 40 K = 1, NDG
2366 KK = K + MM - 2
2367 IF (KK .EQ. 0) THEN
2368 A(K) = XA(I+1, 0) - XA(I, 0)
2369 Else
2370 A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK)
2371 EndIf
2372 40 CONTINUE
2373 DO 41 J = 1, NDH
2374 TEM = 0
2375 DO 43 L = 1, NDG
2376 TEM = TEM + A(L) * GB(L,J,I)
2377 43 CONTINUE
2378 H(J,I,MM) = TEM
2379 41 CONTINUE
2380 35 CONTINUE
2381 DO 42 J = 1, NDG
2382 DO 44 L = 1, NDG
2383 G1(L,J) = G2(L,J+1)
2384 44 CONTINUE
2385 42 CONTINUE
2386 30 CONTINUE
2387 LSTX = .TRUE.
2388 RETURN
2389 END
2390 FUNCTION CtLhXFRMZ (Z)
2391 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2392 LOGICAL LSTX
2393 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2394 PARAMETER (MXX = 105)
2395 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2396 COMMON / LhCtINVERT / ZA
2397 EXTERNAL CtLhZFXL
2398 DATA TEM, RER / D1, 1E-3 /
2399 DATA ZLOW, ZHIGH, IWRN2 / -10.0, 1.00002, 0 /
2400 7 EPS = TEM * RER
2401 ZA = Z
2402 IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN
2403 XLA = LOG (XMIN) * 1.5
2404 XLB = 0.00001
2405 TEM = CtLhZBRNT (CtLhZFXL, XLA, XLB, EPS, IRT)
2406 Else
2407 CALL CtLhWARNR (IWRN2, 'Z out of range in CtLhXFRMZ, X set=0.',
2408 > 'Z', Z, ZLOW, ZHIGH, 1)
2409 TEM = 0
2410 EndIf
2411 CtLhXFRMZ = EXP(TEM)
2412 RETURN
2413 END
2414 FUNCTION CtLhxLi(n,x)
2415 implicit NONE
2416 integer NCUT, i,n,m3
2417 real*8 CtLhxLi,Out,x,pi2by6,zeta3,c1,c2
2418 real*8 r,xt,L,xln1m
2419 parameter (m3=8)
2420 dimension c1(2:m3),c2(2:m3)
2421 data NCUT/27/
2422 data c1/0.75,-0.5833333333333333d0,0.454861111111111d0,
2423 > -0.3680555555555555d0,0.3073611111111111d0,
2424 > -0.2630555555555555d0,0.2294880243764172d0/
2425 data c2/-0.5d0,0.5d0,-0.4583333333333333d0,0.416666666666666d0,
2426 > -0.3805555555555555d0,0.35d0,-0.3241071428571428d0/
2427 data zeta3,pi2by6 /1.20205690315959d0,1.64493406684823d0/
2428 L=0.0
2429 i=0
2430 r=1.0
2431 if (abs(x).gt.r) then
2432 PRINT *,'Li: x out of range (-1,1) , x=',x
2433 STOP
2434 endif
2435 if (n.lt.0) then
2436 PRINT *,'Polylogarithm Li undefined for n=',n
2437 STOP
2438 elseif (n.eq.0) then
2439 Out=x/(1d0-x)
2440 elseif (n.eq.1) then
2441 Out=-dlog(1-x)
2442 elseif (n.eq.2) then
2443 !Calculate dilogarithm
2444 !separately for x<0.5 and x>0.5
2445 if (x.ge.(-0.5).and.x.le.0.5) then
2446 do while(i.le.NCUT)
2447 i=i+1
2448 r=r*x
2449 L=L+r/i/i
2450 enddo
2451 Out=L
2452 elseif (x.eq.0) then
2453 Out=0d0
2454 elseif(x.gt.0.5) then !n.eq.2,x>0.5
2455 xt = 1.0-x
2456 L = pi2by6 - dlog(x)*dlog(xt)
2457 do while(i.le.NCUT)
2458 i=i+1
2459 r=r*xt
2460 L=L-r/i/i
2461 enddo
2462 Out=L
2463 elseif (x.lt.(-0.5)) then
2464 xt=-x/(1d0-x)
2465 L=-0.5*dlog(1-x)**2
2466 do while (i.le.NCUT)
2467 i=i+1
2468 r=r*xt
2469 L=L-r/i/i
2470 enddo
2471 Out=L
2472 endif
2473 elseif (n.eq.3.and.x.ge.0.8) then !use the expansion of Li3 near x=1
2474 L=zeta3+pi2by6*dlog(x)
2475 xt=(1d0-x)
2476 xln1m=dlog(xt)
2477 do i=2,m3
2478 L=L+(c1(i)+c2(i)*xln1m)*xt**i
2479 enddo
2480 Out=L
2481 else !n>3 or x=3,x<0.8
2482 do while(i.le.NCUT)
2483 i=i+1
2484 r=r*x
2485 L=L+r/dble(i)**dble(n)
2486 enddo
2487 Out=L
2488 endif
2489 CtLhxLi=Out
2490 End ! CtLhxLi
2491 FUNCTION CtLhZBRLAM (WLLN)
2492 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2493 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2494 WMULM = EXP (WLLN)
2495 TEM1 = 1./ CtLhALPQCD(JRDR, N1, WMULM, I)
2496 TEM2 = 1./ CtLhALPQCD(JRDR, N, VMULM, I)
2497 CtLhZBRLAM = TEM1 - TEM2
2498 END
2499 FUNCTION CtLhZBRNT(FUNC, X1, X2, TOL, IRT)
2500 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2501 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
2502 external func
2503 IRT = 0
2504 TOL = ABS(TOL)
2505 A=X1
2506 B=X2
2507 FA=FUNC(A)
2508 FB=FUNC(B)
2509 IF(FB*FA.GT.0.) THEN
2510 PRINT *, 'Root must be bracketed for CtLhZBRNT. Set = 0'
2511 IRT = 1
2512 CtLhZBRNT=0.
2513 RETURN
2514 ENDIF
2515 FC=FB
2516 DO 11 ITER=1,ITMAX
2517 IF(FB*FC.GT.0.) THEN
2518 C=A
2519 FC=FA
2520 D=B-A
2521 E=D
2522 ENDIF
2523 IF(ABS(FC).LT.ABS(FB)) THEN
2524 A=B
2525 B=C
2526 C=A
2527 FA=FB
2528 FB=FC
2529 FC=FA
2530 ENDIF
2531 TOL1=2.*EPS*ABS(B)+0.5*TOL
2532 XM=.5*(C-B)
2533 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
2534 CtLhZBRNT=B
2535 RETURN
2536 ENDIF
2537 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
2538 S=FB/FA
2539 IF(A.EQ.C) THEN
2540 P=2.*XM*S
2541 Q=1.-S
2542 ELSE
2543 Q=FA/FC
2544 R=FB/FC
2545 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
2546 Q=(Q-1.)*(R-1.)*(S-1.)
2547 ENDIF
2548 IF(P.GT.0.) Q=-Q
2549 P=ABS(P)
2550 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
2551 E=D
2552 D=P/Q
2553 ELSE
2554 D=XM
2555 E=D
2556 ENDIF
2557 ELSE
2558 D=XM
2559 E=D
2560 ENDIF
2561 A=B
2562 FA=FB
2563 IF(ABS(D) .GT. TOL1) THEN
2564 B=B+D
2565 ELSE
2566 B=B+SIGN(TOL1,XM)
2567 ENDIF
2568 FB=FUNC(B)
256911 CONTINUE
2570 PRINT *, 'CtLhZBRNT exceeding maximum iterations.'
2571 IRT = 2
2572 CtLhZBRNT=B
2573 RETURN
2574 END
2575 FUNCTION CtLhZFRMX (XX)
2576 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2577 LOGICAL LSTX
2578 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2579 PARAMETER (MXX = 105)
2580 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2581 DATA IWRN1, HUGE, TINY / 0, 1.E35, 1.E-35 /
2582 F(X) = (XCR-XMIN) * LOG (X/XMIN) + LOG (XCR/XMIN) * (X-XMIN)
2583 D(X) = (XCR-XMIN) / X + LOG (XCR/XMIN)
2584 X = XX
2585 IF (X .GE. XMIN) THEN
2586 TEM = F(X) / F(D1)
2587 ElseIF (X .GE. D0) THEN
2588 X = MAX (X, TINY)
2589 TEM = F(X) / F(D1)
2590 Else
2591 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhZFRMX'
2592 > , 'X', X, TINY, HUGE, 1)
2593 TEM = 99.
2594 STOP
2595 EndIf
2596 CtLhZFRMX = TEM
2597 RETURN
2598 ENTRY CtLhDZDX (XX)
2599 X = XX
2600 IF (X .GE. XMIN) THEN
2601 TEM = D(X) / F(D1)
2602 ElseIF (X .GE. D0) THEN
2603 X = MAX (X, TINY)
2604 TEM = D(X) / F(D1)
2605 Else
2606 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhDZDX '
2607 > , 'X', X, TINY, HUGE, 1)
2608 TEM = 99.
2609 STOP
2610 EndIf
2611 CtLhDZDX = TEM
2612 RETURN
2613 END
2614 FUNCTION CtLhZFXL (XL)
2615 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2616 COMMON / LhCtINVERT / ZA
2617 X = EXP(XL)
2618 TT = CtLhZFRMX (X) - ZA
2619 CtLhZFXL = TT
2620 RETURN
2621 END
2622