LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / EVLCTEQ.f
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)
211 c
212        call getnset(iset)
213 c       
214         ientry = ientry + 1
215         if(ientry .lt. 100) then
216 666        format(1x,'enter evolve',i3)
217         elseif(ientry .eq. 100) then
218 667        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
295 333      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
346 669        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
477 10      AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0
478         BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0
479         TVAL=0.
480         DO 15 I=1,3
481 15      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)
486 20      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
495 25      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
509 50      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
716 10       CONTINUE
717       DO 20 I= 0, NF
718          IF (NAME .EQ. 'LAM'//CHAR(I+IASC0))
719      1             LhCtNAMQCD=I+13
720 20       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
738 10       CONTINUE
739 20    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
918 c
919       call getnset(iset)
920 c
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
932 c 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
941 c 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         
947 c 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
1002 100     continue
1003
1004 c 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
1054 110     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
1105 c     >    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)
1110 c 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)
1330 11    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
1342 12      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
1350 13    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
1465 c     print *,'warning empty CtLhQCDSET call: NAME=',
1466 c     &          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)
1562 11    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
1586 11    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
1598 12      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
1606 13    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)
2569 11    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