1 subroutine MRSTevolve(x,Q,pdf)
2 implicit real*8(a-h,o-z)
3 include 'parmsetup.inc'
4 character*16 name(nmxset)
5 integer nmem(nmxset),ndef(nmxset),mmem
6 c integer member(nmxset)
8 common/NAME/name,nmem,ndef,mmem
9 parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26,
20 .,fc(nx,nqc),fb(nx,nqb)
22 .cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4),
23 .cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),
24 .cc6(0:nhess,nx,nq,4,4),cc8(0:nhess,nx,nq,4,4),
25 .ccc(0:nhess,nx,nqc,4,4),ccb(0:nhess,nx,nqb,4,4)
26 real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
27 data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
28 . 1d-4,2d-4,4d-4,6d-4,8d-4,
29 . 1d-3,2d-3,4d-3,6d-3,8d-3,
30 . 1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
31 . .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
32 . .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
33 . .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
35 data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
36 . 1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
37 . 1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
38 . 1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
39 . 1.8d6,3.2d6,5.6d6,1d7/
40 data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
52 call getnmem(iset,imem)
53 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
54 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
55 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
56 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
57 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
58 call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
62 call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
67 call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
89 read(1,*)nmem(nset),ndef(nset)
90 c print *,nmem(nset),ndef(nset)
95 read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
96 . f5(n,m),f7(n,m),f6(n,m),f8(n,m)
97 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
99 c write(*,*)'PDF set ',nm,' first element ',f1(1,1)
118 call jeppe1(nm,nx,nq,xxl,qql,f1,cc1)
119 call jeppe1(nm,nx,nq,xxl,qql,f2,cc2)
120 call jeppe1(nm,nx,nq,xxl,qql,f3,cc3)
121 call jeppe1(nm,nx,nq,xxl,qql,f4,cc4)
122 call jeppe1(nm,nx,nq,xxl,qql,f6,cc6)
123 call jeppe1(nm,nx,nq,xxl,qql,f8,cc8)
134 call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc)
142 call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb)
150 entry MRSTalfa(nflav,alfas,Qalfa)
153 c call setnmem(member(iset))
154 call alphamrs(nflav,alfas,Qalfa)
157 entry MRSTinit(Eorder,Q2fit)
161 c if(mem.eq.0) mem=ndef
165 call setnmem(iset,mem)
171 subroutine jeppe1(imem,nx,my,xx,yy,ff,cc)
172 implicit real*8(a-h,o-z)
173 parameter(nnx=49,mmy=37,nhess=30)
174 dimension xx(nx),yy(my),ff(nnx,mmy),
175 xff1(nnx,mmy),ff2(nnx,mmy),
176 xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
177 xcl(16),cc(0:nhess,nx,my,4,4),iwt(16,16)
179 data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
180 x 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
181 x -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
182 x 2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
183 x 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
184 x 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
185 x 0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
186 x 0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
187 x -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
188 x 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
189 x 9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
190 x -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
191 x 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
192 x 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
193 x -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
194 x 4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
199 ff1(1,m)=(ff(2,m)-ff(1,m))/dx
201 ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
203 ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
210 ff2(n,1)=(ff(n,2)-ff(n,1))/dy
212 ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
214 ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
221 ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
223 ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
225 ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
253 yy12(3)=ff12(n+1,m+1)
266 xxd=xxd+iwt(k,l)*z(k)
274 cc(imem,n,m,k,j)=cl(l)
282 subroutine jeppe2(i,x,y,nx,my,xx,yy,cc,z)
283 C-- G.W. 02/07/2007 Allow extrapolation to small x and large q.
284 implicit real*8(a-h,o-z)
286 dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)
291 if (n.gt.0.and.n.lt.nx.and.m.gt.0.and.m.lt.my) then
292 C-- Do usual interpolation.
293 t=(x-xx(n))/(xx(n+1)-xx(n))
294 u=(y-yy(m))/(yy(m+1)-yy(m))
297 z=t*z+((cc(i,n,m,l,4)*u+cc(i,n,m,l,3))*u
298 & +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
301 else if (n.eq.0.and.m.gt.0.and.m.lt.my) then
302 C-- Extrapolate to small x.
303 call jeppe3(i,xx(1),y,nx,my,xx,yy,cc,f0)
304 call jeppe3(i,xx(2),y,nx,my,xx,yy,cc,f1)
305 if (f0.gt.0.d0.and.f1.gt.0.d0) then
306 z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
308 z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
311 else if (n.gt.0.and.m.eq.my) then
312 C-- Extrapolate to large q.
313 call jeppe3(i,x,yy(my),nx,my,xx,yy,cc,f0)
314 call jeppe3(i,x,yy(my-1),nx,my,xx,yy,cc,f1)
315 if (f0.gt.0.d0.and.f1.gt.0.d0) then
316 z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
319 z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
322 else if (n.eq.0.and.m.eq.my) then
323 C-- Extrapolate to small x AND large q.
324 call jeppe3(i,xx(1),yy(my),nx,my,xx,yy,cc,f0)
325 call jeppe3(i,xx(1),yy(my-1),nx,my,xx,yy,cc,f1)
326 if (f0.gt.0.d0.and.f1.gt.0.d0) then
327 z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
330 z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
332 call jeppe3(i,xx(2),yy(my),nx,my,xx,yy,cc,f0)
333 call jeppe3(i,xx(2),yy(my-1),nx,my,xx,yy,cc,f1)
334 if (f0.gt.0.d0.and.f1.gt.0.d0) then
335 z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
338 z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
340 if (z0.gt.0.d0.and.z1.gt.0.d0) then
341 z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
343 z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
347 C-- Set parton distribution to zero otherwise.
355 C-- G.W. 02/07/2007 Copy of the original jeppe2,
356 C-- only used for extrapolation.
357 subroutine jeppe3(i,x,y,nx,my,xx,yy,cc,z)
358 implicit real*8(a-h,o-z)
360 dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)
363 t=(x-xx(n))/(xx(n+1)-xx(n))
364 u=(y-yy(m))/(yy(m+1)-yy(m))
367 z=t*z+((cc(i,n,m,l,4)*u+cc(i,n,m,l,3))*u
368 & +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
373 integer function locx(xx,nx,x)
374 implicit real*8(a-h,o-z)
376 c$$$ if(x.le.xx(1)) then
377 if(x.eq.xx(1)) then ! G.W. 02/07/2007
381 c$$$ if(x.ge.xx(nx)) then
382 if(x.eq.xx(nx)) then ! G.W. 02/07/2007
388 1 if((ju-jl).le.1) go to 2
400 real*8 function polderiv(x1,x2,x3,y1,y2,y3)
401 implicit real*8(a-h,o-z)
402 polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
403 .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))