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