LHAPDF/lhapdf5.3.1/wrapmrstqed.f

   1       subroutine MRSTqedevolve(x,Q,pdf,photon)
   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)
   7       integer nset,iset
   8       common/NAME/name,nmem,ndef,mmem
   9       parameter(nx=49,nq=37,np=9,nqc0=2,nqb0=11,nqc=35,nqb=26,
  10      .nhess=30)
  11       real*8 pdf(-6:6),photon,phot
  12       real*8 f1(nx,nq)
  13      .,f2(nx,nq)
  14      .,f3(nx,nq)
  15      .,f4(nx,nq)
  16      .,f5(nx,nq)
  17      .,f6(nx,nq)
  18      .,f7(nx,nq)
  19      .,f8(nx,nq)
  20      .,f9(nx,nq)
  21      .,fc(nx,nqc),fb(nx,nqb)
  22       real*8 qq(nq),xx(nx),
  23      .cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4),
  24      .cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),
  25      .cc6(0:nhess,nx,nq,4,4),cc8(0:nhess,nx,nq,4,4),
  26      .cc9(0:nhess,nx,nq,4,4),
  27      .ccc(0:nhess,nx,nqc,4,4),ccb(0:nhess,nx,nqb,4,4)
  28       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
  29       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
  30      .        1d-4,2d-4,4d-4,6d-4,8d-4,
  31      .        1d-3,2d-3,4d-3,6d-3,8d-3,
  32      .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
  33      .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
  34      .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
  35      .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
  36      .     .8d0,.9d0,1d0/
  37       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
  38      .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
  39      .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
  40      .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
  41      .        1.8d6,3.2d6,5.6d6,1d7/
  42       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
  43       save
  44
  45       xsave=x
  46       qsq = q*q
  47       q2save=qsq
  48
  49       xlog=dlog(x)
  50       qsqlog=dlog(qsq)
  51
  52       call getnset(iset)
  53 c      imem=member(iset)
  54       call getnmem(iset,imem)
  55       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
  56       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
  57       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
  58       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
  59       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
  60       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
  61       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc9,photon)
  62
  63       chm=0.d0
  64       if(qsq.gt.emc2) then
  65       call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
  66       endif
  67
  68       bot=0.d0
  69       if(qsq.gt.emb2) then
  70       call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
  71       endif
  72
  73       pdf(0)  = glu
  74       pdf(1)  = dnv+dsea
  75       pdf(-1) = dsea
  76       pdf(2)  = upv+usea
  77       pdf(-2) = usea
  78       pdf(3)  = str
  79       pdf(-3) = str
  80       pdf(4)  = chm
  81       pdf(-4) = chm
  82       pdf(5)  = bot
  83       pdf(-5) = bot
  84       pdf(6) = 0.0d0
  85       pdf(-6) = 0.0d0
  86
  87       x=xsave
  88       qsq=q2save
  89       return
  90 *
  91       entry MRSTqedread(nset)
  92       read(1,*)nmem(nset),ndef(nset)
  93 c      print *,nmem(nset),ndef(nset)
  94 c      do nm = 0,nmem-1
  95       do nm = 0,nmem(nset)
  96         do 20 n=1,nx-1
  97         do 20 m=1,nq
  98         read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
  99      .            f5(n,m),f7(n,m),f6(n,m),f8(n,m),f9(n,m)
 100 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea +9 - photon
 101   20  continue
 102 c      write(*,*)'PDF set ',nm,' first element ',f1(1,1)
 103       do 40 m=1,nq
 104       f1(nx,m)=0.d0
 105       f2(nx,m)=0.d0
 106       f3(nx,m)=0.d0
 107       f4(nx,m)=0.d0
 108       f5(nx,m)=0.d0
 109       f6(nx,m)=0.d0
 110       f7(nx,m)=0.d0
 111       f8(nx,m)=0.d0
 112       f9(nx,m)=0.d0
 113   40  continue
 114
 115       do n=1,nx
 116       xxl(n)=dlog(xx(n))
 117       enddo
 118       do m=1,nq
 119       qql(m)=dlog(qq(m))
 120       enddo
 121
 122       call jeppe1(nm,nx,nq,xxl,qql,f1,cc1)
 123       call jeppe1(nm,nx,nq,xxl,qql,f2,cc2)
 124       call jeppe1(nm,nx,nq,xxl,qql,f3,cc3)
 125       call jeppe1(nm,nx,nq,xxl,qql,f4,cc4)
 126       call jeppe1(nm,nx,nq,xxl,qql,f6,cc6)
 127       call jeppe1(nm,nx,nq,xxl,qql,f8,cc8)
 128       call jeppe1(nm,nx,nq,xxl,qql,f9,cc9)
 129
 130       emc2=2.045
 131       emb2=18.5
 132
 133       do 44 m=1,nqc
 134       qqlc(m)=qql(m+nqc0)
 135       do 44 n=1,nx
 136       fc(n,m)=f5(n,m+nqc0)
 137    44 continue
 138       qqlc(1)=dlog(emc2)
 139       call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc)
 140
 141       do 45 m=1,nqb
 142       qqlb(m)=qql(m+nqb0)
 143       do 45 n=1,nx
 144       fb(n,m)=f7(n,m+nqb0)
 145    45 continue
 146       qqlb(1)=dlog(emb2)
 147       call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb)
 148
 149
 150       enddo
 151   50  format(9f10.5)
 152       return
 153 *
 154       end