LHAPDF/lhapdf5.5.1/src/wrapmrst98.f

   1 ! -*- F90 -*-
   2
   3
   4       subroutine MRST98evolve(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       integer nset
  10       common/NAME/name,nmem,ndef,mmem
  11       parameter(nx=49,nq=37,ntenth=23,np=8,members=5)
  12       real*8 pdf(-6:6)
  13       real*8 f(0:members,np,nx,nq+1)
  14       real*8 qq(nq),xx(nx),xxin(nx),g(np),n0(np)
  15       data xxin/1d-5,2d-5,4d-5,6d-5,8d-5,                               &
  16      &              1d-4,2d-4,4d-4,6d-4,8d-4,                           &
  17      &              1d-3,2d-3,4d-3,6d-3,8d-3,                           &
  18      &              1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,               &
  19      &           .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,     &
  20      &           .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,     &
  21      &           .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,        &
  22      &           .8d0,.9d0,1d0/
  23       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,       &
  24      &        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,                          &
  25      &        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,          &
  26      &        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,              &
  27      &        1.8d6,3.2d6,5.6d6,1d7/
  28       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
  29       data n0/3,4,5,9,9,9,9,9/
  30       save
  31 !
  32       xsave=x
  33       qsq = q*q
  34       q2save=qsq
  35 !
  36       if(x.lt.xmin) x=xmin
  37       if(x.gt.xmax) x=xmax
  38       if(qsq.lt.qsqmin)        qsq=qsqmin
  39       if(qsq.gt.qsqmax)        qsq=qsqmax
  40 !
  41       xxx=x
  42       if(x.lt.xx(ntenth)) xxx=dlog10(x/xx(ntenth))+xx(ntenth)
  43       n=0
  44    70 n=n+1
  45       if(xxx.gt.xx(n+1)) goto 70
  46       a=(xxx-xx(n))/(xx(n+1)-xx(n))
  47       m=0
  48    80 m=m+1
  49       if(qsq.gt.qq(m+1)) goto 80
  50       b=(qsq-qq(m))/(qq(m+1)-qq(m))
  51       do 60 i=1,np
  52       g(i)= (1d0-a)*(1d0-b)*f(imem,i,n,m)+(1d0-a)*b*f(imem,i,n,m+1)     &
  53      &          +       a*(1d0-b)*f(imem,i,n+1,m)+a*b*f(imem,i,n+1,m+1)
  54       if(n.ge.ntenth) goto 65
  55       if(i.eq.5.or.i.eq.7) goto 65
  56           fac=(1d0-b)*f(imem,i,ntenth,m)+b*f(imem,i,ntenth,m+1)
  57            g(i)=fac*10d0**(g(i)-fac)
  58    65 continue
  59       g(i)=g(i)*(1d0-x)**n0(i)
  60    60 continue
  61       upv=g(1)
  62       dnv=g(2)
  63       usea=g(4)
  64       dsea=g(8)
  65       str=g(6)
  66       chm=g(5)
  67       glu=g(3)
  68       bot=g(7)
  69 !
  70       pdf(0)  = glu
  71       pdf(1)  = dnv+dsea
  72       pdf(-1) = dsea
  73       pdf(2)  = upv+usea
  74       pdf(-2) = usea
  75       pdf(3)  = str
  76       pdf(-3) = str
  77       pdf(4)  = chm
  78       pdf(-4) = chm
  79       pdf(5)  = bot
  80       pdf(-5) = bot
  81       pdf(6)  = 0.0d0
  82       pdf(-6) = 0.0d0
  83
  84       x=xsave
  85       qsq=q2save
  86       return
  87 !
  88       entry MRST98read(nset)
  89       read(1,*)nmem(nset),ndef(nset)
  90 ! - first resotre the xx array
  91       do j=1,nx
  92         xx(j)=xxin(j)
  93       enddo
  94 ! - next read in the data points
  95       do nm = 0,nmem(nset)
  96         do 20 n=1,nx-1
  97         do 20 m=1,nq
  98         read(1,50)f(nm,1,n,m),f(nm,2,n,m),f(nm,3,n,m),f(nm,4,n,m),      &
  99      &                  f(nm,5,n,m),f(nm,7,n,m),f(nm,6,n,m),f(nm,8,n,m)
 100 ! notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
 101         do 25 i=1,np
 102    25        f(nm,i,n,m)=f(nm,i,n,m)/(1d0-xx(n))**n0(i)
 103    20   continue
 104 !        write(*,*)'PDF set ',nm,' first element ',f(nm,1,1,1)
 105         do 31 j=1,ntenth-1
 106 !        xx(j)=dlog10(xx(j)/xx(ntenth))+xx(ntenth)
 107         do 31 i=1,8
 108         if(i.eq.5.or.i.eq.7) goto 31
 109         do 30 k=1,nq
 110    30   f(nm,i,j,k)=dlog10(f(nm,i,j,k)/f(nm,i,ntenth,k))                &
 111      &              +f(nm,i,ntenth,k)
 112    31   continue
 113    50   format(8f10.5)
 114         do 40 i=1,np
 115         do 40 m=1,nq
 116    40 f(nm,i,nx,m)=0d0
 117        enddo
 118       do 32 j=1,ntenth-1
 119         xx(j)=dlog10(xx(j)/xx(ntenth))+xx(ntenth)
 120    32 continue
 121       return
 122 !
 123       entry MRST98alfa(alfas,Qalfa)
 124         call alphamrs(5,alfas,Qalfa)
 125       return
 126 !
 127       entry MRST98init(Eorder,Q2fit)
 128       return
 129 !
 130       entry MRST98pdf(mem)
 131 !      if(mem.eq.0) mem=ndef
 132       imem = mem
 133 !      print *,imem
 134
 135       return
 136 !
 137       END