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