LHAPDF/lhapdf5.3.1/wrapcteq5.f

   1       subroutine CTEQ5evolve(x,Q,pdf)
   2       implicit real*8(a-h,o-z)
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6       common/NAME/name,nmem,ndef,mmem
   7       integer nset
   8       real*8 pdf(-6:6)
   9       Character Line*80
  10       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
  11       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
  12       PARAMETER (M= 2, M1 = M + 1)
  13       Common
  14      > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
  15      > / CtqPar2 / Nx, Nt, NfMx
  16      > / XQrange / Qini, Qmax, Xmin
  17      > / QCDtable /  Alambda, Nfl, Iorder
  18      > / Masstbl / Amass(6)
  19       data pi / 3.141592653589793d0 /
  20       save
  21 c
  22       U =         X * CtLhCtq5Pdf(1,X,Q)
  23       D =         X * CtLhCtq5Pdf(2,X,Q)
  24       USEA =      X * CtLhCtq5Pdf(-1,X,Q)
  25       DSEA =      X * CtLhCtq5Pdf(-2,X,Q)
  26       STR =       X * CtLhCtq5Pdf(3,X,Q)
  27       CHM =       X * CtLhCtq5Pdf(4,X,Q)
  28       BOT =       X * CtLhCtq5Pdf(5,X,Q)
  29       GLU  =      X * CtLhCtq5Pdf(0,X,Q)
  30       UPV=U-USEA
  31       DNV=D-DSEA
  32 c
  33       pdf(0)  = glu
  34       pdf(1)  = dnv+dsea
  35       pdf(-1) = dsea
  36       pdf(2)  = upv+usea
  37       pdf(-2) = usea
  38       pdf(3)  = str
  39       pdf(-3) = str
  40       pdf(4)  = chm
  41       pdf(-4) = chm
  42       pdf(5)  = bot
  43       pdf(-5) = bot
  44       pdf(6)  = 0.0d0
  45       pdf(-6) = 0.0d0
  46
  47       return
  48 *
  49       entry CTEQ5read(nset)
  50
  51       call CtLhbldat1
  52       call CtLhbldat2
  53
  54       read(1,*)nmem(nset),ndef(nset)
  55       Read  (1, '(A)') Line
  56       Read  (1, '(A)') Line
  57       Read  (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
  58       Iorder = Nint(Dr)
  59       Nfl = Nint(Fl)
  60       Alambda = Al
  61
  62       Read  (1, '(A)') Line
  63       Read  (1, *) NX,  NT, NfMx
  64
  65       Read  (1, '(A)') Line
  66       Read  (1, *) QINI, QMAX, (QL(I), I =0, NT)
  67
  68       Read  (1, '(A)') Line
  69       Read  (1, *) XMIN, (XV(I), I =0, NX)
  70
  71       Do 11 Iq = 0, NT
  72          QL(Iq) = Log (QL(Iq) /Al)
  73    11 Continue
  74 C
  75 C                  Since quark = anti-quark for nfl>2 at this stage,
  76 C                  we Read  out only the non-redundent data points
  77 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
  78
  79       Nblk = (NX+1) * (NT+1)
  80       Npts =  Nblk  * (NfMx+3)
  81
  82       Read  (1, '(A)') Line
  83       Read  (1, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
  84       Return
  85 *
  86
  87       entry CTEQ5alfa(alfas,Qalfa)
  88
  89        alfas = pi*CtLhALPI(Qalfa)
  90
  91       return
  92 *
  93       entry CTEQ5init(Eorder,Q2fit)
  94
  95       return
  96 *
  97       entry CTEQ5pdf(mem)               !error corrected (jcp)
  98 c        imem = mem
  99         call getnset(iset)
 100         call setnmem(iset,mem)
 101       return
 102 *
 103       end
 104
 105       FUNCTION CtLhPartonX5 (IPRTN, X, Q)
 106 C
 107 C   Given the parton distribution function in the array UPD in
 108 C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
 109 C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
 110 C
 111       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 112 C
 113       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
 114       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
 115       PARAMETER (M= 2, M1 = M + 1)
 116 C
 117       Common
 118      > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
 119      > / CtqPar2 / Nx, Nt, NfMx
 120      > / XQrange / Qini, Qmax, Xmin
 121 C
 122       Dimension Fq(M1), Df(M1)
 123
 124 c note, this routine doesn't have the features of not recalculating
 125 c the x or q point values if they have not changed since the last call.
 126 c that makes it slower than cteq6, but since cteq5 is already obsolete,
 127 c speed is not so important for it.
 128 C                                                 Work with Log (Q)
 129       QG  = LOG (Q/AL)
 130
 131 C                           Find lower end of interval containing X
 132       JL = -1
 133       JU = Nx+1
 134  11   If (JU-JL .GT. 1) Then
 135          JM = (JU+JL) / 2
 136          If (X .GT. XV(JM)) Then
 137             JL = JM
 138          Else
 139             JU = JM
 140          Endif
 141          Goto 11
 142       Endif
 143
 144       Jx = JL - (M-1)/2
 145 c crude treatment if outside the defined interval.
 146       If (Jx .LT. 0) then
 147          Jx = 0
 148       Elseif (Jx .GT. Nx-M) Then
 149          Jx = Nx - M
 150       Endif
 151 C                                    Find the interval where Q lies
 152       JL = -1
 153       JU = NT+1
 154  12   If (JU-JL .GT. 1) Then
 155          JM = (JU+JL) / 2
 156          If (QG .GT. QL(JM)) Then
 157             JL = JM
 158          Else
 159             JU = JM
 160          Endif
 161          Goto 12
 162       Endif
 163
 164       Jq = JL - (M-1)/2
 165       If (Jq .LT. 0) Then
 166          Jq = 0
 167 c         If (Q .lt. Qini)  Print '(A, 2(1pE12.4))',
 168 c     >     ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
 169       Elseif (Jq .GT. Nt-M) Then
 170          Jq = Nt - M
 171 c         If (Q .gt. Qmax)  Print '(A, 2(1pE12.4))',
 172 c     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
 173       Endif
 174
 175       If (Iprtn .GE. 3) Then
 176          Ip = - Iprtn
 177       Else
 178          Ip = Iprtn
 179       EndIf
 180 C                             Find the off-set in the linear array Upd
 181       JFL = Ip + NfMx
 182       J0  = (JFL * (NT+1) + Jq) * (NX+1) + Jx
 183 C
 184 C                                           Now interpolate in x for M1 Q's
 185       Do 21 Iq = 1, M1
 186          J1 = J0 + (Nx+1)*(Iq-1) + 1
 187          Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
 188  21   Continue
 189 C                                          Finish off by interpolating in Q
 190       Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)
 191
 192       CtLhPartonX5 = Ftmp
 193 C
 194       RETURN
 195
 196 C                        ****************************
 197       END
 198
 199       Function CtLhCtq5Pdf (Iparton, X, Q)
 200       Implicit Double Precision (A-H,O-Z)
 201       Logical Warn
 202       Common
 203      > / CtqPar2 / Nx, Nt, NfMx
 204      > / QCDtable /  Alambda, Nfl, Iorder
 205
 206       Data Warn /.true./
 207       save Warn
 208
 209       If (X .lt. 0D0 .or. X .gt. 1D0) Then
 210         Print *, 'X out of range in CtLhCtq5Pdf: ', X
 211         Stop
 212       Endif
 213       If (Q .lt. Alambda) Then
 214         Print *, 'Q out of range in CtLhCtq5Pdf: ', Q
 215         Stop
 216       Endif
 217       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
 218          If (Warn) Then
 219 C        put a warning for calling extra flavor.
 220              Warn = .false.
 221              Print *, 'Warning: Iparton out of range in CtLhCtq5Pdf: '
 222      >              , Iparton
 223          Endif
 224          CtLhCtq5Pdf = 0D0
 225          Return
 226       Endif
 227
 228       CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
 229       if(CtLhCtq5Pdf.lt.0.D0)  CtLhCtq5Pdf = 0.D0
 230
 231       Return
 232
 233 C                             ********************
 234       End