LHAPDF/lhapdf5.3.1/wrapgrv.f

   1       subroutine GRVevolve(xin,qin,pdf)
   2       implicit real*8 (a-h,o-z)
   3       include 'parmsetup.inc'
   4       PARAMETER(ngrid=2)
   5       PARAMETER (NPART=6, NX=68, NQ=27, NARG=2)
   6       character*16 name(nmxset)
   7       integer nmem(nmxset),ndef(nmxset),mmem
   8       common/NAME/name,nmem,ndef,mmem
   9       DIMENSION XXUVF(0:ngrid,NX,NQ), XXDVF(0:ngrid,NX,NQ),
  10      +          XXDEF(0:ngrid,NX,NQ), XXUDF(0:ngrid,NX,NQ),
  11      1          XXSF(0:ngrid,NX,NQ),   XXGF(0:ngrid,NX,NQ),
  12      +          XUVF(NX,NQ), XDVF(NX,NQ),
  13      +          XDEF(NX,NQ), XUDF(NX,NQ),
  14      1          XSF(NX,NQ),   XGF(NX,NQ),
  15      +          PARTON (NPART,NQ,NX-1),
  16      2          QS(NQ), XB(NX), XT(NARG), NA(NARG), ARRF(NX+NQ)
  17       CHARACTER*80 LINE
  18       dimension pdf(-6:6)
  19       save
  20       x=xin
  21       q2=qin*qin
  22 *...CHECK OF X AND Q2 VALUES :
  23       IF ( (X.LT.0.99D-9) .OR. (X.GT.1.D0) ) THEN
  24          WRITE(6,91)
  25   91     FORMAT (2X,'PARTON INTERPOLATION: X OUT OF RANGE')
  26          STOP
  27       ENDIF
  28       IF ( (Q2.LT.0.799) .OR. (Q2.GT.1.01E6) ) THEN
  29          WRITE(6,92)
  30   92     FORMAT (2X,'PARTON INTERPOLATION: Q2 OUT OF RANGE')
  31          STOP
  32       ENDIF
  33 *...INTERPOLATION :
  34       DO IQ=1,NQ
  35       DO IX=1,NX
  36       xuvf(ix,iq)=xxuvf(imem,ix,iq)
  37       xdvf(ix,iq)=xxdvf(imem,ix,iq)
  38       xdef(ix,iq)=xxdef(imem,ix,iq)
  39       xudf(ix,iq)=xxudf(imem,ix,iq)
  40       xsf(ix,iq)=xxsf(imem,ix,iq)
  41       xgf(ix,iq)=xxgf(imem,ix,iq)
  42       enddo
  43       enddo
  44       XT(1) = DLOG(X)
  45       XT(2) = DLOG(Q2)
  46       X1 = 1.- X
  47       XV = X**0.5
  48       XS = X**(-0.2)
  49       UV = FINT(NARG,XT,NA,ARRF,XUVF) * X1**3 * XV
  50       DV = FINT(NARG,XT,NA,ARRF,XDVF) * X1**4 * XV
  51       DE = FINT(NARG,XT,NA,ARRF,XDEF) * X1**7 * XV
  52       UD = FINT(NARG,XT,NA,ARRF,XUDF) * X1**7 * XS
  53       US = 0.5 * (UD - DE)
  54       DS = 0.5 * (UD + DE)
  55       SS = FINT(NARG,XT,NA,ARRF,XSF)  * X1**7 * XS
  56       GL = FINT(NARG,XT,NA,ARRF,XGF)  * X1**5 * XS
  57 *
  58
  59       pdf(-6) = 0.0d0
  60        pdf(6) = 0.0d0
  61       pdf(-5) = 0.0d0
  62        pdf(5) = 0.0d0
  63       pdf(-4) = 0.0d0
  64        pdf(4) = 0.0d0
  65       pdf(-3) = ss
  66        pdf(3) = ss
  67       pdf(-2) = us
  68        pdf(2) = uv+us
  69       pdf(-1) = ds
  70        pdf(1) = dv+ds
  71        pdf(0) = gl
  72       return
  73 c
  74 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  75       entry GRVread(nset)
  76 c
  77 c      print *,'calling grvread'
  78       read(1,*)nmem(nset),ndef(nset)
  79
  80       read(1,93)xb
  81 c      print *,xb
  82       read(1,93)qs
  83 c      print *,qs
  84   93  format(8e8.2)
  85 c
  86       do ng=0,nmem(nset)
  87 c
  88       READ(1,89) LINE
  89   89  FORMAT(A80)
  90 c      print *,line
  91       DO 15 M = 1, NX-1
  92       DO 15 N = 1, NQ
  93       READ(1,90) PARTON(1,N,M), PARTON(2,N,M), PARTON(3,N,M),
  94      1           PARTON(4,N,M), PARTON(5,N,M), PARTON(6,N,M)
  95   90  FORMAT (6(1PE10.3))
  96   15  CONTINUE
  97 *
  98 *....ARRAYS FOR THE INTERPOLATION SUBROUTINE :
  99       DO 10 IQ = 1, NQ
 100       DO 20 IX = 1, NX-1
 101         XB0V = XB(IX)**0.5
 102         XB0S = XB(IX)**(-0.2)
 103         XB1 = 1.-XB(IX)
 104         xXUVF(ng,IX,IQ) = PARTON(1,IQ,IX) / (XB1**3 * XB0V)
 105         xXDVF(ng,IX,IQ) = PARTON(2,IQ,IX) / (XB1**4 * XB0V)
 106         xXDEF(ng,IX,IQ) = PARTON(3,IQ,IX) / (XB1**7 * XB0V)
 107         xXUDF(ng,IX,IQ) = PARTON(4,IQ,IX) / (XB1**7 * XB0S)
 108         xXSF(ng,IX,IQ)  = PARTON(5,IQ,IX) / (XB1**7 * XB0S)
 109         xXGF(ng,IX,IQ)  = PARTON(6,IQ,IX) / (XB1**5 * XB0S)
 110   20  CONTINUE
 111         xXUVF(ng,NX,IQ) = 0.E0
 112         xXDVF(ng,NX,IQ) = 0.E0
 113         xXDEF(ng,NX,IQ) = 0.E0
 114         xXUDF(ng,NX,IQ) = 0.E0
 115         xXSF(ng,NX,IQ)  = 0.E0
 116         xXGF(ng,NX,IQ)  = 0.E0
 117   10  CONTINUE
 118       NA(1) = NX
 119       NA(2) = NQ
 120       DO 30 IX = 1, NX
 121         ARRF(IX) = DLOG(XB(IX))
 122   30  CONTINUE
 123       DO 40 IQ = 1, NQ
 124         ARRF(NX+IQ) = DLOG(QS(IQ))
 125   40  CONTINUE
 126
 127       enddo
 128
 129       return
 130 c
 131 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 132       entry GRValfa(alfas,qalfa)
 133       call getnset(iset)
 134       q2alfa = qalfa*qalfa
 135       call GetOrderAsM(iset,iord)
 136       nord=iord+1
 137       alfas=grvals(q2alfa,nord)
 138       alfas = 4.0d0*3.14159d0*alfas
 139       return
 140 c
 141 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 142       entry GRVinit(Eorder,Q2fit)
 143       return
 144 c
 145 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 146       entry GRVpdf(mem)
 147       imem = mem
 148       return
 149 c
 150  1000 format(5e13.5)
 151       end
 152 c
 153       FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
 154 *********************************************************************
 155 *                                                                   *
 156 *   THE INTERPOLATION ROUTINE (CERN LIBRARY ROUTINE E104)           *
 157 *                                                                   *
 158 *********************************************************************
 159       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 160       DIMENSION ARG(5),NENT(5),ENT(10),TABLE(10)
 161       DIMENSION D(5),NCOMB(5),IENT(5)
 162       KD=1
 163       M=1
 164       JA=1
 165          DO 5 I=1,NARG
 166       NCOMB(I)=1
 167       JB=JA-1+NENT(I)
 168          DO 2 J=JA,JB
 169       IF (ARG(I).LE.ENT(J)) GO TO 3
 170     2 CONTINUE
 171       J=JB
 172     3 IF (J.NE.JA) GO TO 4
 173       J=J+1
 174     4 JR=J-1
 175       D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
 176       IENT(I)=J-JA
 177       KD=KD+IENT(I)*M
 178       M=M*NENT(I)
 179     5 JA=JB+1
 180       FINT=0.
 181    10 FAC=1.
 182       IADR=KD
 183       IFADR=1
 184          DO 15 I=1,NARG
 185       IF (NCOMB(I).EQ.0) GO TO 12
 186       FAC=FAC*(1.-D(I))
 187       GO TO 15
 188    12 FAC=FAC*D(I)
 189       IADR=IADR-IFADR
 190    15 IFADR=IFADR*NENT(I)
 191       FINT=FINT+FAC*TABLE(IADR)
 192       IL=NARG
 193    40 IF (NCOMB(IL).EQ.0) GO TO 80
 194       NCOMB(IL)=0
 195       IF (IL.EQ.NARG) GO TO 10
 196       IL=IL+1
 197          DO 50  K=IL,NARG
 198    50 NCOMB(K)=1
 199       GO TO 10
 200    80 IL=IL-1
 201       IF(IL.NE.0) GO TO 40
 202       RETURN
 203       END
 204
 205 c      FUNCTION ALPHAS (Q2, NAORD)
 206       FUNCTION grvals (Q2, NAORD)
 207 *********************************************************************
 208 *                                                                   *
 209 *   THE ALPHA_S ROUTINE.                                            *
 210 *                                                                   *
 211 *   INPUT :  Q2    =  scale in GeV**2  (not too low, of course);    *
 212 *            NAORD =  1 (LO),  2 (NLO).                             *
 213 *                                                                   *
 214 *   OUTPUT:  alphas_s/(4 pi) for use with the GRV(98) partons.      *
 215 *                                                                   *
 216 *******************************************************i*************
 217 *
 218       IMPLICIT DOUBLE PRECISION (A - Z)
 219       INTEGER NF, K, I, NAORD
 220       DIMENSION LAMBDAL (3:6),  LAMBDAN (3:6), Q2THR (3)
 221 *
 222 *...HEAVY QUARK THRESHOLDS AND LAMBDA VALUES :
 223       DATA Q2THR   /  1.960,  20.25,  30625. /
 224       DATA LAMBDAL / 0.2041, 0.1750, 0.1320, 0.0665 /
 225       DATA LAMBDAN / 0.2994, 0.2460, 0.1677, 0.0678 /
 226 *
 227 *...DETERMINATION OF THE APPROPRIATE NUMBER OF FLAVOURS :
 228       NF = 3
 229       DO 10 K = 1, 3
 230       IF (Q2 .GT. Q2THR (K)) THEN
 231          NF = NF + 1
 232       ELSE
 233           GO TO 20
 234        END IF
 235   10   CONTINUE
 236 *
 237 *...LO ALPHA_S AND BETA FUNCTION FOR NLO CALCULATION :
 238   20   B0 = 11.- 2./3.* NF
 239        B1 = 102.- 38./3.* NF
 240        B10 = B1 / (B0*B0)
 241        IF (NAORD .EQ. 1) THEN
 242          LAM2 = LAMBDAL (NF) * LAMBDAL (NF)
 243          ALP  = 1./(B0 * DLOG (Q2/LAM2))
 244          GO TO 1
 245        ELSE IF (NAORD .EQ. 2) then
 246          LAM2 = LAMBDAN (NF) * LAMBDAN (NF)
 247          B1 = 102.- 38./3.* NF
 248          B10 = B1 / (B0*B0)
 249        ELSE
 250          WRITE (6,91)
 251   91     FORMAT ('INVALID CHOICE FOR ORDER IN ALPHA_S')
 252          STOP
 253        END IF
 254 *
 255 *...START VALUE FOR NLO ITERATION :
 256        LQ2 = DLOG (Q2 / LAM2)
 257        ALP = 1./(B0*LQ2) * (1.- B10*DLOG(LQ2)/LQ2)
 258 *
 259 *...EXACT NLO VALUE, FOUND VIA NEWTON PROCEDURE :
 260        DO 2 I = 1, 6
 261        XL  = DLOG (1./(B0*ALP) + B10)
 262        XLP = DLOG (1./(B0*ALP*1.01) + B10)
 263        XLM = DLOG (1./(B0*ALP*0.99) + B10)
 264        Y  = LQ2 - 1./ (B0*ALP) + B10 * XL
 265        Y1 = (- 1./ (B0*ALP*1.01) + B10 * XLP
 266      1       + 1./ (B0*ALP*0.99) - B10 * XLP) / (0.02D0*ALP)
 267        ALP = ALP - Y/Y1
 268   2    CONTINUE
 269 *
 270 *...OUTPUT :
 271 c  1    ALPHAS = ALP
 272   1    grvals = ALP
 273        RETURN
 274        END
 275
 276
 277