LHAPDF/lhapdf5.3.1/wrapa02m.f

   1       subroutine A02Mevolve(xb,Q,fset)
   2       implicit none
   3       include 'parmsetup.inc'
   4       integer npdf,npar,kschem,i,k,n,m,kx,nxbb
   5       integer nxb,nq,np,nvar,pdfmem,imem
   6       parameter(nxb=99,nq=20,np=9,nvar=17)
   7       integer nexp(0:np)
   8       data nexp / 0, 3, 4, 5, 5, 5, 5, 5, 5, 5 /
   9       character*16 name(nmxset)
  10       integer nmem(nmxset),ndef(nmxset),mem
  11       common/NAME/name,nmem,ndef,mem
  12 c      integer iset,iimem
  13 c      common/SET/iset,iimem
  14       integer nset
  15       real*8 f(0:nvar,nxb,nq+1,0:np),xx(nxb)
  16       real*8 fsp(nxb),bs(nxb),cs(nxb),ds(nxb)
  17       real*8 bsp(0:nvar,nxb,nq+1,0:np),csp(0:nvar,nxb,nq+1,0:np)
  18      +      ,dsp(0:nvar,nxb,nq+1,0:np)
  19       real*8 x1,xd,del,dels,delx,delx1,xlog1
  20       real*8 pdfs(0:np),fset(-6:6),alfas
  21       real*8 x,Q,xb,q2,xmin,xmax,qsq,qsqmin,qsqmax,b,ss
  22       real*8 aa,f0,fp,fm
  23       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
  24
  25       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb,xx
  26      +,fsp,bsp,csp,dsp
  27 c      save
  28
  29       q2=Q*Q
  30       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
  31       if(xb.lt.xmin.or.xb.gt.xmax)       print 98,x
  32   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ')
  33   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ')
  34
  35       x=max(xb,xmin)
  36       x=min(x,xmax)
  37       qsq=max(q2,qsqmin)
  38       qsq=min(qsq,qsqmax)
  39
  40       if (x.gt.x1) then
  41         xd=(1d0-x1)**2-(1d0-x)**2
  42         n=int(xd/delx1)+nxbb
  43       else
  44         xd=dlog(x)-xlog1
  45         n=nxbb+int(xd/DELX)-1
  46       end if
  47       aa=x-xx(n)
  48
  49       ss=dlog(dlog(qsq/0.04d0))-dlog(dlog(qsqmin/0.04d0))
  50       m=int(ss/dels)+1
  51       b=ss/dels-dble(m)+1.d0
  52
  53       k=pdfmem
  54
  55       do i=0,npdf
  56         f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i)
  57      +              + aa**3*dsp(k,n,m,i)
  58         fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i)
  59      +                + aa**3*dsp(k,n,m+1,i)
  60         if (m.ge.2) then
  61           fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i)
  62      +                   +aa**3*dsp(k,n,m-1,i)
  63           pdfs(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
  64         else
  65           pdfs(i)= f0*(1d0-b) + fp*b
  66         end if
  67         pdfs(i) = pdfs(i)*(1d0-x)**dble(nexp(i))
  68       end do
  69
  70       fset(-6)=pdfs(9)
  71       fset(-5)=pdfs(8)
  72       fset(-4)=pdfs(7)
  73       fset(-3)=pdfs(5)
  74 c--reversed mrs 7/7/04 due
  75       fset(-1)=pdfs(6)
  76       fset(-2)=pdfs(4)
  77       fset(0)=pdfs(3)
  78 c--reversed mrs 7/7/04 due
  79       fset(2)=pdfs(1)+pdfs(4)
  80       fset(1)=pdfs(2)+pdfs(6)
  81       fset(3)=pdfs(5)
  82       fset(4)=pdfs(7)
  83       fset(5)=pdfs(8)
  84       fset(6)=pdfs(9)
  85       return
  86 *
  87       entry A02Malfa(alfas,Q)
  88       q2=Q*Q
  89       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
  90       qsq=max(q2,qsqmin)
  91       qsq=min(qsq,qsqmax)
  92
  93       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
  94       m=int(ss/dels)+1
  95       b=ss/dels-m+1
  96
  97       k=pdfmem
  98       alfas=(1d0-b)*f(k,1,m,0)+b*f(k,1,m+1,0)
  99       return
 100 *
 101       entry A02Mread(nset)
 102 c following fix because members are 0-nvar
 103 c      nmem = nvar + 1
 104
 105       nmem(nset) = nvar
 106       ndef(nset) = 0
 107
 108
 109       read(1,*) npdf
 110 c      print *,'number of members',npdf
 111       npar=nvar
 112       do k=0,npar
 113          do n=1,nxb-1
 114             do m=1,nq
 115                read(1,*) (f(k,n,m,i),i=0,npdf)
 116 c              print 100,(f(k,n,m,i),i=0,npdf)
 117             enddo
 118          enddo
 119       enddo
 120   100 format (13f11.5)
 121
 122       nxbb=nxb/2
 123       x1=0.3d0
 124       xlog1=dlog(x1)
 125       delx=(dlog(x1)-dlog(xmin))/dble(nxbb-1)
 126       DELX1=(1.d0-x1)**2/dble(nxbb+1)
 127
 128 *...X GRID
 129       do kx=1,nxbb
 130         xx(kx)=dexp(dlog(xmin)+delx*dble(kx-1))
 131       end do
 132       do kx=nxbb+1,nxb-1
 133         xx(kx)=1.d0-dsqrt(dabs((1.d0-x1)**2-delx1*dble(kx-nxbb)))
 134       end do
 135       xx(nxb)=1d0
 136
 137       do k=0,npar
 138       do i=0,npdf
 139         do m=1,nq
 140           if (i.ne.0) then
 141             f(k,nxb,m,i)=0d0
 142           else
 143             f(k,nxb,m,i)=f(k,nxb-1,m,i)
 144           end if
 145           do n=1,nxb-1
 146             f(k,n,m,i)=f(k,n,m,i)/(1d0-xx(n))**nexp(i)
 147           end do
 148           do n=1,nxb
 149             fsp(n)=f(k,n,m,i)
 150           end do
 151           call a02mspline (nxb,xx,fsp,bs,cs,ds)
 152           do n=1,nxb
 153             bsp(k,n,m,i)=bs(n)
 154             csp(k,n,m,i)=cs(n)
 155             dsp(k,n,m,i)=ds(n)
 156           end do
 157         end do
 158       end do
 159       end do
 160
 161       return
 162 *
 163       entry A02Minit
 164
 165       dels=(dlog(dlog(qsqmax/0.04d0))-
 166      +      dlog(dlog(qsqmin/0.04d0)))/dble(nq-1)
 167
 168
 169       return
 170 *
 171       entry A02Mpdf(imem)
 172       pdfmem=imem
 173       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
 174          write(*,*) 'A02M PDF set:'
 175          write(*,*) 'PDF member out of range:'
 176          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
 177          stop
 178       endif
 179       return
 180       end
 181
 182 * ---------------------------------------------------------------------
 183       SUBROUTINE A02MSPLINE(N,X,Y,B,C,D)
 184 * ---------------------------------------------------------------------
 185 * CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
 186 * INTERPOLATION SUBROUTINES ARE TAKEN FROM
 187 * G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
 188 * COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
 189 *
 190       IMPLICIT REAL*8(A-H,O-Z)
 191
 192       DIMENSION X(N), Y(N), B(N), C(N), D(N)
 193 *
 194       NM1=N-1
 195       IF(N.LT.2) RETURN
 196       IF(N.LT.3) GO TO 250
 197       D(1)=X(2)-X(1)
 198       C(2)=(Y(2)-Y(1))/D(1)
 199       DO 210 K=2,NM1
 200          D(K)=X(K+1)-X(K)
 201          B(K)=2.0D0*(D(K-1)+D(K))
 202          C(K+1)=(Y(K+1)-Y(K))/D(K)
 203          C(K)=C(K+1)-C(K)
 204   210 CONTINUE
 205       B(1)=-D(1)
 206       B(N)=-D(N-1)
 207       C(1)=0.0D0
 208       C(N)=0.0D0
 209       IF(N.EQ.3) GO TO 215
 210       C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
 211       C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
 212       C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
 213       C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
 214  215  CONTINUE
 215       DO 220 K=2,N
 216          T=D(K-1)/B(K-1)
 217          B(K)=B(K)-T*D(K-1)
 218          C(K)=C(K)-T*C(K-1)
 219  220  CONTINUE
 220       C(N)=C(N)/B(N)
 221       DO 230 IB=1,NM1
 222          K=N-IB
 223          C(K)=(C(K)-D(K)*C(K+1))/B(K)
 224  230  CONTINUE
 225       B(N)=(Y(N)-Y(NM1))/D(NM1)
 226      1     +D(NM1)*(C(NM1)+2.0D0*C(N))
 227       DO 240 K=1,NM1
 228          B(K)=(Y(K+1)-Y(K))/D(K)
 229      1        -D(K)*(C(K+1)+2.0D0*C(K))
 230          D(K)=(C(K+1)-C(K))/D(K)
 231          C(K)=3.0D0*C(K)
 232  240  CONTINUE
 233       C(N)=3.0D0*C(N)
 234       D(N)=D(N-1)
 235       RETURN
 236  250  CONTINUE
 237       B(1)=(Y(2)-Y(1))/(X(2)-X(1))
 238       C(1)=0.0D0
 239       D(1)=0.0D0
 240       B(2)=B(1)
 241       C(2)=0.0D0
 242       D(2)=0.0D0
 243       RETURN
 244       END