LHAPDF/lhapdf5.3.1/wrapa02.f

   1       subroutine A02evolve(xb,Q,fset)
   2       implicit none
   3       include 'parmsetup.inc'
   4       integer npdf,npar,kschem,i,k,n,m,nxbb
   5       integer nxb,nq,np,nvar,pdfmem,imem
   6       parameter(nxb=99,nq=20,np=9,nvar=15)
   7       character*16 name(nmxset)
   8       integer nmem(nmxset),ndef(nmxset),mem
   9       common/NAME/name,nmem,ndef,mem
  10       integer nset
  11       real*4 f(0:nvar,nxb,nq+1,0:np)
  12       real*8 x1,xd,del,dels,delx,delx1,xlog1
  13       real*8 pdfs(np),fset(-6:6),alfas
  14       real*8 x,Q,xb,q2,xmin,xmax,qsq,qsqmin,qsqmax,a,b,ss
  15       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
  16
  17       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb
  18
  19       q2=Q*Q
  20       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
  21       if(xb.lt.xmin.or.xb.gt.xmax)       print 98,x
  22   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ')
  23   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ')
  24
  25       x=max(xb,xmin)
  26       x=min(x,xmax)
  27       qsq=max(q2,qsqmin)
  28       qsq=min(qsq,qsqmax)
  29
  30       if (x.gt.x1) then
  31         xd=(1-x1)**2-(1-x)**2
  32         n=int(xd/delx1)+nxbb
  33         a=xd/delx1-n+nxbb
  34       else
  35         xd=log(x)-xlog1
  36         n=nxbb+int(xd/DELX)-1
  37         a=xd/delx-n+nxbb
  38       end if
  39
  40       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
  41       m=int(ss/dels)+1
  42       b=ss/dels-m+1
  43
  44       k=pdfmem
  45
  46       do i=1,npdf
  47         if (n.gt.1.and.m.gt.1.and.n.ne.49) then
  48         pdfs(i)= f(k,n,m,i)*(1+a*b-a**2-b**2)
  49      +         + f(k,n+1,m+1,i)*a*b
  50      +         + f(k,n+1,m,i)*a*(a-2*b+1)/2.
  51      +         + f(k,n,m+1,i)*b*(b-2*a+1)/2.
  52      +         +  f(k,n-1,m,i)*a*(a-1)/2.
  53      +         + f(k,n,m-1,i)*b*(b-1)/2.
  54         else
  55         pdfs(i)= (1-a)*(1-b)*f(k,n,m,i)
  56      .         + (1-a)*b*f(k,n,m+1,i)
  57      .         + a*(1-b)*f(k,n+1,m,i)
  58      .         + a*b*f(k,n+1,m+1,i)
  59         end if
  60       end do
  61
  62
  63       fset(-6)=pdfs(9)
  64       fset(-5)=pdfs(8)
  65       fset(-4)=pdfs(7)
  66       fset(-3)=pdfs(5)
  67 c--reversed mrs 7/7/04 due
  68       fset(-1)=pdfs(6)
  69       fset(-2)=pdfs(4)
  70       fset(0)=pdfs(3)
  71 c--reversed mrs 7/7/04 due
  72       fset(2)=pdfs(1)+pdfs(4)
  73       fset(1)=pdfs(2)+pdfs(6)
  74       fset(3)=pdfs(5)
  75       fset(4)=pdfs(7)
  76       fset(5)=pdfs(8)
  77       fset(6)=pdfs(9)
  78       return
  79 *
  80       entry A02alfa(alfas,Q)
  81       q2=Q*Q
  82       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
  83       qsq=max(q2,qsqmin)
  84       qsq=min(qsq,qsqmax)
  85
  86       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
  87       m=int(ss/dels)+1
  88       b=ss/dels-m+1
  89
  90       k=pdfmem
  91       alfas=(1d0-b)*f(k,1,m,0)+b*f(k,1,m+1,0)
  92       return
  93 *
  94       entry A02read(nset)
  95 c following fix because members are 0-nvar
  96 c      nmem = nvar + 1
  97       nmem(nset) = nvar
  98       ndef(nset) = 0
  99       read(1,*) npdf
 100       npar=nvar
 101       do k=0,npar
 102          do n=1,nxb-1
 103             do m=1,nq
 104                read(1,100) (f(k,n,m,i),i=0,npdf)
 105             enddo
 106          enddo
 107       enddo
 108   100 format (13f11.5)
 109       return
 110 *
 111       entry A02init
 112
 113       dels=(log(log(qsqmax/0.04))-log(log(qsqmin/0.04)))/(nq-1)
 114
 115       nxbb=nxb/2
 116       x1=0.3
 117       xlog1=log(x1)
 118       delx=(log(x1)-log(xmin))/(nxbb-1)
 119       DELX1=(1-x1)**2/(nxbb+1)
 120
 121       do i=1,npdf
 122          do m=1,nq
 123             do k=1,npar
 124                f(k,nxb,m,i)=0d0
 125             end do
 126          end do
 127          do m=1,nq
 128             do k=1,npar
 129                f(k,nxb,m,0)=f(k,nxb-1,m,0)
 130             end do
 131          end do
 132       end do
 133       return
 134 *
 135       entry A02pdf(imem)
 136       pdfmem=imem
 137       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
 138          write(*,*) 'A02 PDF set:'
 139          write(*,*) 'PDF member out of range:'
 140          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
 141          stop
 142       endif
 143       return
 144       end