LHAPDF/lhapdf5.2.2/wrapQCDNUM3.f

   1       subroutine QCDNUM3evolve(x,Q,pdf)
   2       implicit none
   3       character*64 gridname
   4       character*16 s1,dummy
   5       integer Eorder,index
   6       real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
   7       real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
   8       real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
   9       real*8 QALFAS,QPDFXQ,XFROMIX
  10       real*8 f(-6:6),pdf(-6:6)
  11       integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
  12       integer nx,nq
  13       integer nset,iset,isetlast
  14       data isetlast/-1/
  15       real*8 xmin,xmax,qmin,qmax,S
  16       save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
  17       save nx,xmin,xmax,nq,qmin,qmax,gridname
  18 *
  19       call getnset(iset)
  20       if (iset.ne.isetlast) then
  21         call get_pdfqcd(iset)
  22         isetlast = iset
  23       endif
  24 c
  25       Q2=Q*Q
  26       nf=3
  27       pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
  28       singlet=QPDFXQ('SINGLET',x,Q2,iflag)
  29       dm= QPDFXQ('DM',x,Q2,IFLAG)
  30       um= QPDFXQ('UM',x,Q2,IFLAG)
  31       dp= QPDFXQ('DP',x,Q2,IFLAG)
  32       up= QPDFXQ('UP',x,Q2,IFLAG)
  33       sp= QPDFXQ('SP',x,Q2,IFLAG)
  34       ub=0.5d0*(up-um+singlet/dble(nf))
  35       db=0.5d0*(dp-dm+singlet/dble(nf))
  36       sb=0.5d0*(sp+singlet/dble(nf))
  37
  38       pdf(1)=dm+db
  39       pdf(2)=um+ub
  40       pdf(3)=sb
  41       pdf(4)=0d0
  42       pdf(5)=0d0
  43       pdf(6)=0d0
  44       pdf(-1)=db
  45       pdf(-2)=ub
  46       pdf(-3)=sb
  47       pdf(-4)=0d0
  48       pdf(-5)=0d0
  49       pdf(-6)=0d0
  50       return
  51 *
  52       entry QCDNUM3alfa(alfas,Q)
  53       Q2=Q*Q
  54       nf=3
  55       alfas=QALFAS(Q2,Qlam,nf,iflag)
  56       return
  57 *
  58       entry QCDNUM3read(nset)
  59 c      read(1,*) s1
  60 c      print *,s1
  61 c         gridname='large.grid'
  62 c         nx=400
  63 c         xmin=1d-6
  64 c         xmax=1d0
  65 c         nq=112
  66 c         qmin=1d0
  67 c         qmax=1d10
  68 c         read(1,*) dummy
  69 *      else
  70          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
  71 c*      endif
  72 c      iset = nset
  73       return
  74 *
  75       entry QCDNUM3init(nset,Eorder,Q2fit)
  76       call qninit
  77       call QNISET('ORDER',Eorder+1)
  78       call grxdef(nx,xmin)
  79       call grqdef(nq,qmin,qmax)
  80       call grqinp(Q2fit,1)
  81       iq0=IQFROMQ(Q2fit)
  82 c      nf0= NFLGET(iq0)
  83       nf0=3
  84       CALL QNBOOK(2,'dm')
  85       CALL QNBOOK(3,'um')
  86       CALL QNBOOK(4,'dp')
  87       CALL QNBOOK(5,'up')
  88       CALL QNBOOK(6,'sp')
  89       CALL QNLSET('W1ANA',.TRUE.)
  90       CALL QNLSET('W2NUM',.TRUE.)
  91       CALL QNLSET('W2STF',.FALSE.)
  92       call qthres(1d10,2d10)
  93       if (index(gridname,'none').eq.1) then
  94          call qnfilw(0,0)
  95       else
  96          qnerr=-1
  97          open(unit=2,status='old',file=gridname,
  98      .        form='unformatted',err=1)
  99          call QNREAD(2,1,qnerr)
 100  1       close(2)
 101          if (qnerr.ne.0) then
 102             write(*,*) 'Grid file problem: ',gridname
 103             if (qnerr.lt.0) then
 104                write(*,*) 'Grid file does not exist'
 105                write(*,*) 'Calculating and creating grid file'
 106                call qnfilw(0,0)
 107                open(unit=2,status='unknown',file=gridname,
 108      .              form='unformatted')
 109                call QNDUMP(2)
 110                close(2)
 111             else
 112                write(*,*) 'Existing grid file is inconsistent'
 113                if (qnerr.eq.1)
 114      .              write(*,*) 'Defined grid different'
 115                if (qnerr.eq.2)
 116      .              write(*,*) 'Heavy quark weight table different'
 117                if (qnerr.eq.3)
 118      .              write(*,*) 'Charm mass different'
 119                if (qnerr.eq.4)
 120      .              write(*,*) 'Bottom mass different'
 121                stop
 122             endif
 123          endif
 124       endif
 125       return
 126 *
 127       entry QCDNUM3pdf(nset)
 128 c      print *,'entering QCDNUMpdf',nset
 129       call GetAlfas(nset,alfas0,scale0)
 130 c      print *,alfas0,scale0
 131       Q2=scale0*scale0
 132       CALL QNRSET('ALFAS',alfas0)
 133       CALL QNRSET('ALFQ0',Q2)
 134       DO ix = 1,nx
 135          xx = XFROMIX(ix)
 136 c        print *,'calling parmPDF',ix
 137          call parmPDF(nset,xx,f)
 138 c        if(ix.lt.6) print *,nset,xx,f
 139          singlet=0.d0
 140          do i=1,nf0
 141             singlet=singlet+f(i)+f(-i)
 142 c           print  *,i,singlet
 143          end do
 144          gl=f(0)
 145          dm=f(1)-f(-1)
 146          um=f(2)-f(-2)
 147          dp=f(1)+f(-1)-singlet/dble(nf0)
 148          up=f(2)+f(-2)-singlet/dble(nf0)
 149          sp=f(3)+f(-3)-singlet/dble(nf0)
 150          CALL QNPSET('SINGLET',ix,iq0,singlet)
 151          CALL QNPSET('GLUON',ix,iq0,gl)
 152          CALL QNPSET('DM',ix,iq0,DM)
 153          CALL QNPSET('UM',ix,iq0,UM)
 154          CALL QNPSET('DP',ix,iq0,DP)
 155          CALL QNPSET('UP',ix,iq0,UP)
 156          CALL QNPSET('SP',ix,iq0,SP)
 157       ENDDO
 158       CALL EVOLSG(iq0,1,nq)
 159       CALL EVOLNM('DM',iq0,1,nq)
 160       CALL EVOLNM('UM',iq0,1,nq)
 161 c
 162 c      CALL EVPLUS('DP',iq0,1,nq)
 163 c      CALL EVPLUS('UP',iq0,1,nq)
 164 c      CALL EVPLUS('SP',iq0,1,nq)
 165       CALL EVOLNP('DP',iq0,1,nq)
 166       CALL EVOLNP('UP',iq0,1,nq)
 167       CALL EVOLNP('SP',iq0,1,nq)
 168 c
 169       call getnset(iset)
 170       call save_pdfqcd(iset)
 171       return
 172 *
 173       end
 174