]>
Commit | Line | Data |
---|---|---|
4e9e3152 | 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 |