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