LHAPDF/lhapdf-5.9.1/src/wrapcteq65-lite.f

   1 ! -*- F90 -*-
   2
   3
   4       subroutine CTEQ65evolve(x,Q,pdf)
   5       implicit real*8(a-h,o-z)
   6       include 'parmsetup.inc'
   7       character*16 name(nmxset)
   8       character*512 setpath
   9       integer nmem(nmxset),ndef(nmxset),mmem
  10       common/NAME/name,nmem,ndef,mmem
  11       double precision gridx(nmxgridx),gridq(nmxgridq)
  12       integer ngridx,ngridq,jx,jq
  13       real*8 pdf(-6:6)
  14       integer nset
  15       Character Line*80
  16       PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 0)
  17       PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
  18       Common                                                            &
  19      & / CtqPar1nhess65 /                                               &
  20      &  Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset),                 &
  21      &  UPD(0:nhess,MXPQX,nmxset)                                       &
  22      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
  23      & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset)             &
  24      & / QCDtable /  Alambda(nmxset), Nfl(nmxset), Iorder(nmxset)       &
  25      & / Masstbl / Amass(6,nmxset)
  26
  27       common/masses_LHA/cMass(nmxset),bMass(nmxset),tMass(nmxset)
  28
  29       data pi / 3.141592653589793d0 /
  30       save
  31 !
  32       call getnset(iset)
  33       call getnmem(iset,imem)
  34
  35       U =         X * CtLhCtq65Pdf(0,1,X,Q)
  36       D =         X * CtLhCtq65Pdf(0,2,X,Q)
  37       USEA =      X * CtLhCtq65Pdf(0,-1,X,Q)
  38       DSEA =      X * CtLhCtq65Pdf(0,-2,X,Q)
  39       STR =       X * CtLhCtq65Pdf(0,-3,X,Q)
  40       CHM =       X * CtLhCtq65Pdf(0,-4,X,Q)
  41       BOT =       X * CtLhCtq65Pdf(0,-5,X,Q)
  42       GLU  =      X * CtLhCtq65Pdf(0,0,X,Q)
  43 !
  44       pdf(0)  = glu
  45       pdf(1)  = d
  46       pdf(-1) = dsea
  47       pdf(2)  = u
  48       pdf(-2) = usea
  49       pdf(3)  = str
  50       pdf(-3) = str
  51       pdf(4)  = chm
  52       pdf(-4) = chm
  53       pdf(5)  = bot
  54       pdf(-5) = bot
  55       pdf(6)  = 0.0d0
  56       pdf(-6) = 0.0d0
  57       return
  58 !
  59                                           !like cteq65evolve, but allows
  60       entry CTEQ65cevolve(x,Q,pdf)
  61 !
  62       call getnset(iset)
  63       call getnmem(iset,imem)
  64
  65       U =         X * CtLhCtq65Pdf(0,1,X,Q)
  66       D =         X * CtLhCtq65Pdf(0,2,X,Q)
  67       USEA =      X * CtLhCtq65Pdf(0,-1,X,Q)
  68       DSEA =      X * CtLhCtq65Pdf(0,-2,X,Q)
  69       STR =       X * CtLhCtq65Pdf(0,3,X,Q)
  70       SBAR =      X * CtLhCtq65Pdf(0,-3,X,Q)
  71       CHM =       X * CtLhCtq65Pdf(0,4,X,Q)
  72       CBAR =      X * CtLhCtq65Pdf(0,-4,X,Q)
  73       BOT =       X * CtLhCtq65Pdf(0,5,X,Q)
  74       GLU  =      X * CtLhCtq65Pdf(0,0,X,Q)
  75 !
  76       pdf(0)  = glu
  77       pdf(1)  = d
  78       pdf(-1) = dsea
  79       pdf(2)  = u
  80       pdf(-2) = usea
  81       pdf(3)  = str
  82       pdf(-3) = sbar
  83       pdf(4)  = chm
  84       pdf(-4) = cbar
  85       pdf(5)  = bot
  86       pdf(-5) = bot
  87       pdf(6)  = 0.0d0
  88       pdf(-6) = 0.0d0
  89
  90       return
  91 !
  92                                           !like cteq65evolve, but allows
  93       entry CTEQ65sevolve(x,Q,pdf)
  94 !
  95       call getnset(iset)
  96       call getnmem(iset,imem)
  97
  98       U =         X * CtLhCtq65Pdf(0,1,X,Q)
  99       D =         X * CtLhCtq65Pdf(0,2,X,Q)
 100       USEA =      X * CtLhCtq65Pdf(0,-1,X,Q)
 101       DSEA =      X * CtLhCtq65Pdf(0,-2,X,Q)
 102       STR =       X * CtLhCtq65Pdf(0,3,X,Q)
 103       SBAR =      X * CtLhCtq65Pdf(0,-3,X,Q)
 104       CHM =       X * CtLhCtq65Pdf(0,-4,X,Q)
 105       BOT =       X * CtLhCtq65Pdf(0,5,X,Q)
 106       GLU  =      X * CtLhCtq65Pdf(0,0,X,Q)
 107 !
 108       pdf(0)  = glu
 109       pdf(1)  = d
 110       pdf(-1) = dsea
 111       pdf(2)  = u
 112       pdf(-2) = usea
 113       pdf(3)  = str
 114       pdf(-3) = sbar
 115       pdf(4)  = chm
 116       pdf(-4) = chm
 117       pdf(5)  = bot
 118       pdf(-5) = bot
 119       pdf(6)  = 0.0d0
 120       pdf(-6) = 0.0d0
 121
 122       return
 123 !
 124       entry CTEQ65getgrid(nset,ngridx,ngridq,gridx,gridq)
 125       do jx=0,nx(nset)
 126           gridx(jx+1)=xv(jx,nset)
 127       enddo
 128       do jq=0,nt(nset)
 129           qtemp = alambda(nset)*dexp(dexp(tv(jq,nset)))
 130           gridq(jq+1)=qtemp*qtemp
 131       enddo
 132       ngridx=nx(nset)
 133       ngridq=nt(nset)
 134       return
 135
 136       entry CTEQ65read(nset)
 137       call CtLhbldat1
 138       call CtLhbldat2
 139
 140       call LHct65set
 141
 142                                             !*** nmem+1=number of member
 143       read(1,*)nmem(nset),ndef(nset)
 144 !      if(nmem(nset) .gt. nhess) then
 145 !         print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
 146 !         stop
 147 !      endif
 148
 149       MxVal = 3
 150       Read  (1, '(A)') Line
 151       Read  (1, '(A)') Line
 152       Read  (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
 153       Iorder(nset) = Nint(Dr)
 154       Nfl(nset) = Nint(Fl)
 155       Alambda(nset) = Al(nset)
 156
 157       cMass(nset) = Amass(4,nset)
 158       bMass(nset) = Amass(5,nset)
 159       tMass(nset) = Amass(6,nset)
 160
 161       Read  (1, '(A)') Line
 162 !                                               This is the .pds (WKT) f
 163       Read  (1, *) N0, N0, N0, NfMx(nset), N0, N0
 164       Read  (1, '(A)') Line
 165       Read  (1, *) NX(nset),  NT(nset), N0, N0, N0
 166       Read  (1, '(A)') (Line,I=1,4)
 167       Read  (1, *)                                                      &
 168      &  QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
 169
 170       Read  (1, '(A)') Line
 171       Read  (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
 172       XV(0,nset)=0D0
 173 !
 174 !                  Since quark = anti-quark for nfl>2 at this stage,
 175 !                  we Read  out only the non-redundent data points
 176 !     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
 177
 178
 179       Nblk = (NX(nset)+1) * (NT(nset)+1)
 180       Npts =  Nblk  * (NfMx(nset)+1+MxVal)
 181
 182                                              !*** new version: allows nm
 183       do ihess = 0,nmem(nset)
 184
 185         Read  (1, '(A)') Line
 186         Read  (1, '(A)') Line
 187         Read  (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
 188
 189       enddo
 190       return
 191 !
 192       entry CTEQ66read(nset)
 193       call CtLhbldat1
 194       call CtLhbldat2
 195
 196       call LHct65set
 197
 198                                             !*** nmem+1=number of member
 199       read(1,*)nmem(nset),ndef(nset)
 200 !      if(nmem(nset) .gt. nhess) then
 201 !         print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
 202 !         stop
 203 !      endif
 204
 205       MxVal = 2
 206       Read  (1, '(A)') Line
 207       Read  (1, '(A)') Line
 208       Read  (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
 209       Iorder(nset) = Nint(Dr)
 210       Nfl(nset) = Nint(Fl)
 211       Alambda(nset) = Al(nset)
 212
 213       cMass(nset) = Amass(4,nset)
 214       bMass(nset) = Amass(5,nset)
 215       tMass(nset) = Amass(6,nset)
 216
 217       Read  (1, '(A)') Line
 218 !                                               This is the .pds (WKT) f
 219       Read  (1, *) N0, N0, N0, NfMx(nset), N0, N0
 220       Read  (1, '(A)') Line
 221       Read  (1, *) NX(nset),  NT(nset), N0, N0, N0
 222       Read  (1, '(A)') (Line,I=1,4)
 223       Read  (1, *)                                                      &
 224      &  QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
 225
 226       Read  (1, '(A)') Line
 227       Read  (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
 228       XV(0,nset)=0D0
 229 !
 230 !                  Since quark = anti-quark for nfl>2 at this stage,
 231 !                  we Read  out only the non-redundent data points
 232 !     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
 233
 234
 235       Nblk = (NX(nset)+1) * (NT(nset)+1)
 236       Npts =  Nblk  * (NfMx(nset)+1+MxVal)
 237
 238                                              !*** new version: allows nm
 239       do ihess = 0,nmem(nset)
 240
 241         Read  (1, '(A)') Line
 242         Read  (1, '(A)') Line
 243         Read  (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
 244
 245       enddo
 246       return
 247 !
 248                                              !like CTEQ65read, but c.ne.
 249       entry CTEQ65cread(nset)
 250       call CtLhbldat1
 251       call CtLhbldat2
 252       call LHct65set
 253
 254                                             !*** nmem+1=number of member
 255       read(1,*)nmem(nset),ndef(nset)
 256 !      if(nmem(nset) .gt. nhess) then
 257 !         print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
 258 !         stop
 259 !      endif
 260
 261                                                !one more species free
 262       MxVal = 4
 263       Read  (1, '(A)') Line
 264       Read  (1, '(A)') Line
 265       Read  (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
 266       Iorder(nset) = Nint(Dr)
 267       Nfl(nset) = Nint(Fl)
 268       Alambda(nset) = Al(nset)
 269
 270       cMass(nset) = Amass(4,nset)
 271       bMass(nset) = Amass(5,nset)
 272       tMass(nset) = Amass(6,nset)
 273
 274       Read  (1, '(A)') Line
 275 !                                               This is the .pds (WKT) f
 276       Read  (1, *) N0, N0, N0, NfMx(nset), N0, N0
 277       Read  (1, '(A)') Line
 278       Read  (1, *) NX(nset),  NT(nset), N0, N0, N0
 279       Read  (1, '(A)') (Line,I=1,4)
 280       Read  (1, *)                                                      &
 281      &  QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
 282
 283       Read  (1, '(A)') Line
 284       Read  (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
 285       XV(0,nset)=0D0
 286 !
 287 ! we Read only the non-redundent data points
 288
 289       Nblk = (NX(nset)+1) * (NT(nset)+1)
 290       Npts =  Nblk  * (NfMx(nset)+1+MxVal)
 291
 292                                              !*** new version: allows nm
 293       do ihess = 0,nmem(nset)
 294
 295         Read  (1, '(A)') Line
 296         Read  (1, '(A)') Line
 297         Read  (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
 298
 299       enddo
 300       return
 301 !
 302       entry CTEQ65alfa(alfas,Qalfa)
 303       alfas = pi*CtLhALPI(Qalfa)
 304       return
 305 !
 306       entry CTEQ65init(Eorder,Q2fit)
 307       return
 308 !
 309       entry CTEQ65pdf(mem)
 310 !        imem = mem
 311         call getnset(iset)
 312         call setnmem(iset,mem)
 313
 314 ! have to reopen stream 1
 315         call getsetpath(setpath)
 316         open(1,file=setpath(1:len_trim(setpath)),action='READ')
 317
 318         line = ''
 319         do while (line(1:3).ne.'==>')
 320           read(1,'(a)'),line
 321         enddo
 322 ! - backspace by one record
 323         backspace(1)
 324 ! - dummy read up to the member requested
 325         do n=0,mem-1
 326           read(1,'(a)'),line
 327           read(1,'(a)'),line
 328           Read  (1, *, IOSTAT=IRET) (UPD(0,I,iset), I=1,Npts)
 329         enddo
 330 !- read in the data of the requested member
 331           Read  (1, '(A)') Line
 332           Read  (1, '(A)') Line
 333           Read  (1, *, IOSTAT=IRET) (UPD(0,I,iset), I=1,Npts)
 334           close(1)
 335       return
 336 !
 337       END
 338
 339       subroutine LHct65set
 340       Implicit Double Precision (A-H,O-Z)
 341       common /ctq65co/ xlast, qlast, nxsave
 342       nxsave = -1000
 343       xlast = -2.
 344       qlast = -2.
 345       return
 346       END
 347
 348 !=======================================================================
 349       Function CtLhPartonX65 (imem,IPRTN, XX, QQ)
 350 !  Given the parton distribution function in the array U in
 351 !  COMMON / PEVLDT / , this routine interpolates to find
 352 !  the parton distribution at an arbitray point in x and q.
 353 !
 354       Implicit Double Precision (A-H,O-Z)
 355
 356       include 'parmsetup.inc'
 357       PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 0)
 358       PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
 359
 360       Common                                                            &
 361      & / CtqPar1nhess65 /                                               &
 362      &  Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset),                 &
 363      &  UPD(0:nhess,MXPQX,nmxset)                                       &
 364      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
 365      & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset)
 366       common /ctq65co/ xlast,qlast, nxsave
 367
 368         parameter(nqvec = 4)
 369
 370       Dimension fvec(4), fij(4)
 371       Dimension xvpow(0:mxx)
 372       Data OneP / 1.00001 /
 373                                 !**** choice of interpolation variable
 374       Data xpow / 0.3d0 /
 375       Save xvpow
 376       Data ixprint,iqprint/0,0/
 377       save ixprint,iqprint
 378
 379          save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
 380          save const1 , const2, const3, const4, const5, const6
 381          save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
 382
 383
 384       call getnset(iset)
 385 !      call getnmem(iset,imem)
 386
 387 ! store the powers used for interpolation on first call...
 388       if(nx(iset) .ne. nxsave) then
 389          nxsave = nx(iset)
 390          xvpow(0) = 0.D0
 391          do i = 1, nx(iset)
 392             xvpow(i) = xv(i,iset)**xpow
 393          enddo
 394       endif
 395
 396       X = XX
 397       Q = QQ
 398
 399       if((x.lt.xmin(iset)).or.(x.gt.1.d0)) then
 400         ixprint=ixprint+1
 401         if(ixprint.lt.11) print 98,x
 402         if(ixprint.eq.10) print *,                                      &
 403      &  'more warning messages like the last suppressed.'
 404       endif
 405    98 format(' WARNING:  X=',e12.5,' OUT OF RANGE')
 406       if((q.lt.qini(iset)).or.(q.gt.qmax(iset))) then
 407         iqprint=iqprint+1
 408         if(iqprint.lt.11) print 99,q
 409         if(iqprint.eq.10) print *,                                      &
 410      &  'more warning messages like the last suppressed.'
 411       endif
 412    99 format(' WARNING:  Q=',e12.5,' OUT OF RANGE')
 413
 414 ! skip the initialization in x if same as in the previous call.
 415         if(x .eq. xlast) goto 100
 416         xlast = x
 417
 418 !      -------------    find lower end of interval containing x, i.e.,
 419 !                       get jx such that xv(jx) .le. x .le. xv(jx+1)...
 420       JLx = -1
 421       JU = Nx(iset)+1
 422    11 If (JU-JLx .GT. 1) Then
 423          JM = (JU+JLx) / 2
 424          If (X .Ge. XV(JM,iset)) Then
 425             JLx = JM
 426          Else
 427             JU = JM
 428          Endif
 429          Goto 11
 430       Endif
 431 !                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
 432 !                           |---|---|---|...|---|-x-|---|...|---|---|
 433 !                     x     0  Xmin               x                 1
 434 !
 435       If     (JLx .LE. -1) Then
 436         Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX65 x=',x
 437         Stop
 438       ElseIf (JLx .Eq. 0) Then
 439          Jx = 0
 440       Elseif (JLx .LE. Nx(iset)-2) Then
 441
 442 !                For interior points, keep x in the middle, as shown abo
 443          Jx = JLx - 1
 444       Elseif (JLx.Eq.Nx(iset)-1 .or. x.LT.OneP) Then
 445
 446 !                  We tolerate a slight over-shoot of one (OneP=1.00001)
 447 !              perhaps due to roundoff or whatever, but not more than th
 448 !                                      Keep at least 4 points >= Jx
 449          Jx = JLx - 2
 450       Else
 451         Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX65 x=',x
 452         Stop
 453       Endif
 454 !          ---------- Note: JLx uniquely identifies the x-bin; Jx does n
 455
 456 !                       This is the variable to be interpolated in
 457       ss = x**xpow
 458
 459       If (JLx.Ge.2 .and. JLx.Le.Nx(iset)-2) Then
 460
 461 !     initiation work for "interior bins": store the lattice points in s
 462       svec1 = xvpow(jx)
 463       svec2 = xvpow(jx+1)
 464       svec3 = xvpow(jx+2)
 465       svec4 = xvpow(jx+3)
 466
 467       s12 = svec1 - svec2
 468       s13 = svec1 - svec3
 469       s23 = svec2 - svec3
 470       s24 = svec2 - svec4
 471       s34 = svec3 - svec4
 472
 473       sy2 = ss - svec2
 474       sy3 = ss - svec3
 475
 476 ! constants needed for interpolating in s at fixed t lattice points...
 477       const1 = s13/s23
 478       const2 = s12/s23
 479       const3 = s34/s23
 480       const4 = s24/s23
 481       s1213 = s12 + s13
 482       s2434 = s24 + s34
 483       sdet = s12*s34 - s1213*s2434
 484       tmp = sy2*sy3/sdet
 485       const5 = (s34*sy2-s2434*sy3)*tmp/s12
 486       const6 = (s1213*sy2-s12*sy3)*tmp/s34
 487
 488       EndIf
 489
 490   100      continue
 491
 492 ! skip the initialization in q if same as in the previous call.
 493         if(q .eq. qlast) goto 110
 494         qlast = q
 495
 496       tt = log(log(Q/Al(iset)))
 497
 498 !         --------------Now find lower end of interval containing Q, i.e
 499 !                          get jq such that qv(jq) .le. q .le. qv(jq+1).
 500       JLq = -1
 501       JU = NT(iset)+1
 502    12 If (JU-JLq .GT. 1) Then
 503          JM = (JU+JLq) / 2
 504          If (tt .GE. TV(JM,iset)) Then
 505             JLq = JM
 506          Else
 507             JU = JM
 508          Endif
 509          Goto 12
 510        Endif
 511
 512       If     (JLq .LE. 0) Then
 513          Jq = 0
 514       Elseif (JLq .LE. Nt(iset)-2) Then
 515 !                                  keep q in the middle, as shown above
 516          Jq = JLq - 1
 517       Else
 518 !                         JLq .GE. Nt-1 case:  Keep at least 4 points >=
 519         Jq = Nt(iset) - 3
 520
 521       Endif
 522 !                                   This is the interpolation variable i
 523
 524       If (JLq.GE.1 .and. JLq.LE.Nt(iset)-2) Then
 525 !                                        store the lattice points in t..
 526       tvec1 = Tv(jq,iset)
 527       tvec2 = Tv(jq+1,iset)
 528       tvec3 = Tv(jq+2,iset)
 529       tvec4 = Tv(jq+3,iset)
 530
 531       t12 = tvec1 - tvec2
 532       t13 = tvec1 - tvec3
 533       t23 = tvec2 - tvec3
 534       t24 = tvec2 - tvec4
 535       t34 = tvec3 - tvec4
 536
 537       ty2 = tt - tvec2
 538       ty3 = tt - tvec3
 539
 540       tmp1 = t12 + t13
 541       tmp2 = t24 + t34
 542
 543       tdet = t12*t34 - tmp1*tmp2
 544
 545       EndIf
 546
 547   110      continue
 548
 549 ! get the pdf function values at the lattice points...
 550 ! In this code, we store 10 flavors: u,ubar,d,dbar,s,sbar,c,cbar,b=bbar,
 551 ! hence Iprtn=5 (b) is obtained from -5 (bbar)
 552
 553       If (Iprtn .GE. 5) Then
 554          Ip = - Iprtn
 555       Else
 556          Ip = Iprtn
 557       EndIf
 558       jtmp = ((Ip + NfMx(iset))*(NT(iset)+1)+(jq-1))*(NX(iset)+1)+jx+1
 559
 560       Do it = 1, nqvec
 561
 562          J1  = jtmp + it*(NX(iset)+1)
 563
 564        If (Jx .Eq. 0) Then
 565 !                          For the first 4 x points, interpolate x^2*f(x
 566 !                           This applies to the two lowest bins JLx = 0,
 567 !            We cannot put the JLx.eq.1 bin into the "interior" section
 568 !                           (as we do for q), since Upd(J1) is undefined
 569          fij(1) = 0
 570          fij(2) = Upd(imem,J1+1,iset) * XV(1,iset)**2
 571          fij(3) = Upd(imem,J1+2,iset) * XV(2,iset)**2
 572          fij(4) = Upd(imem,J1+3,iset) * XV(3,iset)**2
 573 !
 574 !                 Use CtLhPolint which allows x to be anywhere w.r.t. th
 575
 576          Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
 577
 578          If (x .GT. 0D0)  fvec(it) =  Fx / x**2
 579 !                                              Pdf is undefined for x.eq
 580        ElseIf  (JLx .Eq. Nx(iset)-1) Then
 581 !                                                This is the highest x b
 582
 583 !** fix allow 4 consecutive elements with iset... mrw 19.9.2005
 584         fij(1) = Upd(imem,j1,iset)
 585         fij(2) = Upd(imem,j1+1,iset)
 586         fij(3) = Upd(imem,j1+2,iset)
 587         fij(4) = Upd(imem,j1+3,iset)
 588         Call CtLhPolint4 (XVpow(Nx(iset)-3), Fij(1), 4, ss, Fx, Dfx)
 589
 590         fvec(it) = Fx
 591
 592        Else
 593
 594 ! for all interior points, use Jon's in-line function
 595 ! This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
 596 ! (This is cubic spline interpolation, as used by cteq; it was
 597 ! changed to polint in previous Durham releases (jcp).)
 598          sf2 = Upd(imem,J1+1,iset)
 599          sf3 = Upd(imem,J1+2,iset)
 600
 601          Fvec(it) = (const5*(Upd(imem,J1,iset)                          &
 602      &                             - sf2*const1 + sf3*const2)           &
 603      &             + const6*(Upd(imem,J1+3,iset)                        &
 604      &                             + sf2*const3 - sf3*const4)           &
 605      &             + sf2*sy3 - sf3*sy2) / s23
 606
 607        Endif
 608
 609       enddo
 610 !                                   We now have the four values Fvec(1:4
 611 !     interpolate in t...
 612
 613       If (JLq .LE. 0) Then
 614 !                         1st Q-bin, as well as extrapolation to lower Q
 615         Call CtLhPolint4(TV(0,iset), Fvec(1), 4, tt, ff, Dfq)
 616
 617       ElseIf (JLq .GE. Nt(iset)-1) Then
 618 !                         Last Q-bin, as well as extrapolation to higher
 619         Call CtLhPolint4(TV(Nt(iset)-3,iset), Fvec(1), 4, tt, ff, Dfq)
 620       Else
 621 !                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
 622 !       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
 623 !                         the full range QV(0:Nt)  (in contrast to XV)
 624         tf2 = fvec(2)
 625         tf3 = fvec(3)
 626
 627         g1 = ( tf2*t13 - tf3*t12) / t23
 628         g4 = (-tf2*t34 + tf3*t24) / t23
 629
 630         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12                      &
 631      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
 632
 633         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
 634       EndIf
 635
 636       CtLhPartonX65 = ff
 637
 638       Return
 639       END
 640 !=======================================================================
 641       Function CtLhCtq65Pdf (imem,Iparton, X, Q)
 642       Implicit Double Precision (A-H,O-Z)
 643       include 'parmsetup.inc'
 644       Logical Warn
 645       Common                                                            &
 646      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
 647      & / QCDtable /  Alambda(nmxset), Nfl(nmxset), Iorder(nmxset)
 648
 649       Data Warn /.true./
 650       save Warn
 651       call getnset(iset)
 652
 653       If (X .lt. 0D0 .or. X .gt. 1D0) Then
 654         Print *, 'X out of range in CtLhCtq65Pdf: ', X
 655         Stop
 656       Endif
 657       If (Q .lt. Alambda(iset)) Then
 658         Print *, 'Q out of range in CtLhCtq65Pdf: ', Q
 659         Stop
 660       Endif
 661
 662 ! added to force pdf = 0.0 at x=1.0 exactly - mrw
 663       if(x .eq. 1.0d0) then
 664           CtLhCtq65Pdf = 0.0d0
 665           return
 666       endif
 667 !
 668       If ((Iparton .lt. -NfMx(iset) .or. Iparton .gt. NfMx(iset))) Then
 669          If (Warn) Then
 670 !        put a warning for calling extra flavor.
 671              Warn = .false.
 672              Print *, 'Warning: Iparton out of range in CtLhCtq65Pdf: ' &
 673      &              , Iparton
 674          Endif
 675          CtLhCtq65Pdf = 0D0
 676          Return
 677       Endif
 678
 679       CtLhCtq65Pdf = CtLhPartonX65 (imem,Iparton, X, Q)
 680       if(CtLhCtq65Pdf.lt.0.D0)  CtLhCtq65Pdf = 0.D0
 681
 682       Return
 683
 684 !                             ********************
 685       END