C 26 Oct 96 John Collins: Added Msg argument to call to WARNR around line 924 Function SetPdf () C Lead Group of ++++++ PdfPac.f ++++++ C (These general comments are included in the C lead subprogram in order to survive forsplit.) C This package is organized in 8 Groups C FileList: C 5setpdf: :setpdf datpdf pdfset C 5pdffns: :pdf pdfprt pdftst pardis prdis1 xf0 xf1 C 5runevl: :rnevlu chgevl prpdf C 5evolve: :evolve nsevl nsrhsm nsrhsp snevl snrhs kernel pff1 stupkl C 5evlaux: :evlwt wtupd evlrd1 hqrk fpin fpinms spnint spence C 5evlutl: :qarray xarray integr hinteg zfrmx C 5parpdf: :parpdf evlpar evlget evlgt1 evlin evlset C 5pdfmis: :pmsdis c2qx sqrk x2zconv amomen C EndFileList C ==================================================================== C GroupList: This is the lead Group with setup functions C setpdf datpdf pdfset C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5setpdf.f,v 5.9 96/10/20 13:42:08 wkt Exp $ C $Log: 5setpdf.f,v $ c Revision 5.9 96/10/20 13:42:08 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.6 96/06/02 23:36:31 wkt c minor c c Revision 5.5 96/06/02 23:21:16 wkt c Write and Read to file converted to 'Formatted' c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/20 15:40:23 wkt c Evlini and InRead removed; c ((They belong to FitPac; cause linking problem c x2zconv retained. c CC Function SetPdf () IMPLICIT DOUBLE PRECISION (A-H, O-Z) C Force link to DatPdf + other setup work External DatPDF SetPdf = 0. Return C **************************** END BLOCK DATA DATPDF C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C CHARACTER*10 NAMPDF, NMHDRN, NAMQRK CHARACTER*12 MRSFLN LOGICAL LSTX C PARAMETER (Z = 1D -10, ZZ = 1D -20) C MxF is the maximum number of quark flavors C MxAdF (= 1 gluon + 1 singlet quark + Additional fake flavors C to save space for storing info on evolution) PARAMETER (MxX = 525, MxQ = 25, MxF = 6, MxAdF = 6) C PARAMETER (MxPN = MxF * 2 + MxAdF) C JCC?? Changed to agree with other definitions. PARAMETER (MxPN = MxF * 2 + 2) PARAMETER (MxQX= MxQ * MxX, MxPQX = MxQX * MxPN) PARAMETER (MxPDF = 20, MXHDRN = 6) C COMMON / IOUNIT / NIN, NOUT, NWRT Common / PdCntrl/ LPrt, LDbg COMMON / XXARAY / XCR, XMIN, XV(0:MxX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MxQ),TV(0:MxQ), NT,JT,NG COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt COMMON / PARPRZ / XMNPRZ(9), QMNPRZ(9), ALF(9), > NfLPRZ(9), NfAL(9), MORD(9) COMMON / INITIA / AN(-MxF : MxF), AI(-MxF : MxF), > BI(-MxF : MxF), CI(-MxF : MxF) C COMMON / NAMPDF / NAMPDF(MxPDF), NMHDRN(-1 : MxHDRN), > NAMQRK(-MxF : MxF) COMMON / MRSDAT / MRSFLN (3) C Data LPrt, LDbg / 1, 0 / DATA QINI, QMAX, XMIN, XCR / 1.9, 1.001D2, 0.999D-3, 0.3 / DATA KF, IKNL, IPD0, IHDN / 10, 1, 1, 1 / DATA NX, NT, JT, LSTX / 40, 6, 1, .FALSE. / DATA (NfLPRZ(I), I=1,9) / 6,4,5,5,5,6,6,6,6/ DATA (XMNPRZ(I), I=1,9) / 3*1.d-4,6*1.D-5/ DATA (QMNPRZ(I), I=1,9) / 2.25, 2.0,3.2, 2*2.25, 4*2.0 / C Validity needs to be checked. DATA (MORD(I), I=1,9) / 2*1, 7*2/ DATA (ALF(I), I=1,9) / 0.2, 0.2, 0.3, 0.19, 0.215, $ 0.22, 0.225,0.2,0.2/ DATA (NfAL(I), I=1,9) / 3*4, 6*5 / C DATA MRSFLN / 'prmz:hmrse', 'prmz:hmrsb', 'GRID3' / DATA (NAMPDF(I),I=1,11) / 'Ehlq_1 ','Dk-Ow 1 ', 'DFLM_NLLA', >'KMRS B0 ','MRS92_D0 ','MT90_S1 ','CTEQ1M ', ' ', >' ', 'QCD_EVL_1','QCD_EVL_2' / C DATA (NMHDRN(I), I=-1,5)/ 'AntiProton', 'Neutron', 'Proton', > 'IsoScalar', 'Pi_Plus', 'Pi_Minus', 'K_Plus' / C DATA (NAMQRK(I), I=-6,6) / '-t_Quark', '-b_Quark', '-c_Quark', > '-s_Quark', '-d_Quark', '-u_Quark', 'Gluon', 'u_Quark', > 'd_Quark', 's_Quark', 'c_Quark', 'b_Quark', 't_Quark' / C C The following data correspond to EHLQ set 1 DATA AI / 7 * Z, 0.5, 0.4, 4 * Z / DATA BI / 6 * Z, 3.5, 2 * 1.51, 4 * 1./ DATA CI / 3 * 5, 3 * 8.54, 5.9, 3.5, 4.5, 4 * 5./ DATA AN / 3 *ZZ, .081, 2 *.182, 2.62, 1.78, 0.67, 4 *ZZ / C C **************************** END SUBROUTINE PDFSET (I1, IHDRN, ALAM, TPMS, QINI, QMAX, XMIN, > NU, HEADER, I2, I3, IRET, IRR, UserIn) C C ------------------------------ C Arguments: C C > signifies an input parameter, C < signifies an output parameter, C * indicates that the parameter is Input if I1=0 ; C Output if I1=1 ; C Inactive if I1=2 . C C > I1 = 0 if evolution calculation is to be carried out, C 1 if data table is to be read from an existing file; C (NU must be set to the Unit # of the file) C -1 same as 1, except data is read to alternate C common blocks for "second set" of PDF's C (NU must be set to the Unit # of the file.) C 2 if data table in memory is to be written to an C external file with no calculation (NU must be C specified) C C * Ihdrn = hadron ID code: proton is 1, anti-proton is -1; C C * Alam = QCD 'Lamda parameter' in the 'CWZ' scheme C (See Collins & Tung 'Calculating heavy ... ', C Fermilab-Pub-86/39-T for definition and C conversion graph to '1st order' and 'MS-bar' C Lamda values.) C C * Tpms = Mass of the Top-quark in GeV; C C * Qini = Q-initial (GeV), i.e. starting value of Q-evolution; C C * Qmax = Maxmum value of Q (GeV) to be used; C C * Xmin = Smallest value of X desired; C C > NU = Unit # to (from) which data is to be written (read); C C * Header= Informational header for the file C C > I2 (if I1 = 0 is chosen): code to choose initial distribu- C tion functions at Q = Qini (Ignored if I1 = 1, 2, -1) C C = 0 : Use initial distribution function C UserIn (Iparton, X) to be supplied by user; C UserIn must be declared EXTERNAL in the C calling program. C C 1 - 99: Use "canned" parton distribution sets as input. C The list of parametrizations is updated C periodically; it is given in the current version C of the PDF function. C C > I3 (if I1 = 0, 2 is chosen): switch to write to file C C = 0 : results of calculation not written to a file; C C 1 : results written to the file (with unit # NU) C C 2 : results written to the alternate common blocks C C < Iret : return error code C C = 0 : no error C C = 1 : Error opening file to read (check Irr) C C = 2 : Error condition in reading file (check Irr) C C = 10 : Xmin > 1. is not allowed C C = 11 : Xmin (< 1E-7) is too small for use of the C current formalism; Default Nx, Ke kept C C = 20 : I2 value is illegal C C = 22 : All parameters are unchanged from previous C run; calculation by-passed C C = 30 : I3 < 0 .OR. > 2 is illegal C C = 31 : Error condition in opening file for write C C = 32 : Write error from EVLWT (check Irr value) C C = 33 : Read error from EVLRD1 (check Irr value) C C = 40 : I1 < 0 .or. > 1 is illegal C C = 98 : Evolution calculation by-passed because C all relevant parameters remained the same C as before; data table in memory should C remain valid. C C < Irr = error code returned by Fortran when error condition C is encountered in reading from or writing C to files (check Fortran manual). C C ---------------------------- C C Additional parameters which can be 'dialed' by the user but which do C appear in the above argument list because they C are less likely to be changed: C C Physics parameters: C C NFL is the number of quark flavors C legal values: 0, 1, ... 6 ; default = 6; C C M1, M2, M3, M4, M5 are masses of quarks u, d, s, c, b respec. C C IKNL selects the evolution kernel on the RHS of the AP equation C = 1, 2 : 1st & 2nd order unpolarized case C -1, -2 : ...... polarized .. C GRID points in the variables X and T (=Log Q): C C NX (< MXX), default adjusted according to choise of Xmin; C NT (< MXQ), C (cf. program codes Parameter ... for values of Mxx and Mxq) C C The (integer) parameter JT (with default set at 1) subdivides C the step-size of the variable Q in the solution of the evolu- C tion equation. Using JT > 1 increases the numerical precision C of the results without increasing the size of the data block, C but, of course, at the expense of multiplying the execution C time. C C A similar parameter (JX) will be added in a future version to C control the X-mesh points used in actual calulation, C independent of the declared array size. This will allow us C to reduce the size of the data block without sacrificing C much numerical accuracy -- a desirable change for memory- C constrained applications. C C All parameters listed in this section can be changed ('set') by the C following fortran call to a general purpose input/output routine: C C Call ParPdf (1, 'NAME', VALUE, Iret) C C prior to invoking the PDFSET subroutine. C Here NAME is the (character) name of the variable as listed above, C VALUE is the (real) value of the variable to be set, and Iret is an C error return code (e.g. Iret=0 means 'variable NAME cannot be found'). C C Likewise, the current value of any variable can be recalled by the C following call: C C ParPdf (2, 'NAME', VALUE, Iret) C C For other use of PARPDF and further details consult the program listing C in the section EVLPAR inside PDFIT.FOR. C C For interactive use of these programs: C The entire list of adjustable parameters can be brought to the screen C and systematically changed in an interactive way by the Fortran call: C C Call ChgEvl (InputSet) C C where InputSet is an optional "setup subroutine" for any input initial C parton distribution function (at Q0) (denoted by UserIn in the above C argument list for PdfSet) which needs to be called before the C actual function evaluations (setting up common block of coefficients...etc.) C If needed, InputSet must be declared External in the calling program. C If not needed, it can be ignored. C The user will be prompted for possible inputs and other options. C C ---------------------------- C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C CHARACTER HEADER*78, FORMT*11 C LOGICAL LWT, LOP Common / IOUnit / Nin, Nout, Nwrt C EXTERNAL FPIN, UserIn C IRET = 0 IRR = 0 LWT = .FALSE. C ---------------------------- C Read only section: IF (I1 .EQ. 1) THEN C OPEN (NU, FORM='FORMATTED', STATUS='OLD',IOSTAT=IRR,ERR=10) INQUIRE (NU, OPENED=LOP, FORM=FORMT) IF (.NOT. LOP) > STOP 'Data table file not openned for read in PDFSET' C Read file CALL EVLRD (NU, HEADER, IRR) IF (IRR .NE. 0) THEN IRET = 2 RETURN EndIf C Recall physical parameters used in generating data in the file C CALL PARPDF (2, 'ALAM', ALAM, J1) CALL PARPDF (2, 'IHDN', HDRN, J2) CALL PARPDF (2, 'M6 ', TPMS, J3) CALL PARPDF (2, 'QINI', QINI, J4) CALL PARPDF (2, 'QMAX', QMAX, J5) CALL PARPDF (2, 'XMIN', XMIN, J6) IHDRN = NINT(HDRN) C Add Error testing codes (Js) if desired 10 IF (IRR .NE. 0) IRET = 1 C ElseIF (I1 .EQ. -1) THEN C INQUIRE (NU, OPENED=LOP, FORM=FORMT) IF (.NOT. LOP) > STOP 'Data table file not openned for read in PDFSET' C Read file CALL EVLRD1 (NU, HEADER, IRR) IF (IRR .NE. 0) THEN IRET = 33 RETURN EndIf C Recall physical parameters used in generating data in the file C CALL PARPDF (2, 'ALAM', ALAM, J1) CALL PARPDF (2, 'IHDN', HDRN, J2) CALL PARPDF (2, 'M6 ', TPMS, J3) CALL PARPDF (2, 'QINI', QINI, J4) CALL PARPDF (2, 'QMAX', QMAX, J5) CALL PARPDF (2, 'XMIN', XMIN, J6) IHDRN = NINT(HDRN) C Add Error testing codes (Js) if desired 12 IF (IRR .NE. 0) IRET = 1 C ---------------------------- C Main section to do evolution: C Set parameters for evolution ElseIF (I1 .EQ. 0) THEN C HDRN = DBLE(IHDRN) CALL PARPDF (1, 'ALAM', ALAM, J1) CALL PARPDF (1, 'IHDN', HDRN, J2) CALL PARPDF (1, 'M6 ', TPMS, J3) CALL PARPDF (1, 'QINI', QINI, J4) CALL PARPDF (1, 'QMAX', QMAX, J5) CALL PARPDF (1, 'XMIN', XMIN, J5) C Add Error testing code if desired C IF (XMIN .GT. 1D0) THEN IRET = 10 RETURN EndIf C C Call evolution routine according to the choice of initial distribution C IF (I2 .LT. 0) THEN IRET = 20 RETURN ElseIF (I2 .EQ. 0) THEN AI2 = 0. CALL PARPDF (1, 'IPD0', AI2, L1) IRR = 99 Aout = Nout Print *, 'Start evolution with parameters:' CALL PARPDF (4, 'ALL', Aout, L1) CALL EVOLVE (UserIn, IRR) C Add error checking here ElseIF (I2 .LE. 99) THEN AI2 = I2 CALL PARPDF (1, 'IPD0', AI2, L1) CALL EVOLVE (FPIN, IRR) IF (IRR .EQ. 98) IRET = 98 Else IRET = 20 PRINT *, 'I2 out of range in PDFSET; I2 =', I2 RETURN EndIf IF (IRR .EQ. 98) IRET = 22 C Write results to file if I3=1 IF (I3 .LT. 0 .OR. I3 .GT. 2) THEN IRET = 30 ElseIF (I3 .GT. 0) THEN LWT = .TRUE. EndIf C ---------------------------- ElseIF (I1 .EQ. 2) THEN LWT = .TRUE. Else IRET = 40 EndIf C IF (LWT) THEN INQUIRE (NU, OPENED=LOP, FORM=FORMT) IF (.NOT. LOP) > STOP 'Data table file not openned for write in PDFSET' C CALL EVLWT (NU, HEADER, IRR) IF (IRR .NE. 0) THEN PRINT *, 'Error Writing Data to File in PDFSET' IRET = 32 EndIf REWIND (NU) C IF (I3 .EQ. 2) THEN INQUIRE (NU, OPENED=LOP, FORM=FORMT) IF (.NOT. LOP) > STOP 'Data table file not openned for read in PDFSET' CALL EVLRD1 (NU, HEADER, IRR) IF (IRR .NE. 0) THEN PRINT *, > 'Error Reading Data to the Alternate Block in PDFSET' IRET = 33 EndIf EndIf 30 IF (IRR .NE. 0) IRET = 31 EndIf C RETURN C **************************** END FUNCTION PDF (Iset, Ihadron, Iparton, X, Q, IR) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: User Callable Parton distribution functions C 5pdffns: : pdf pdfprt pardis prdis1 PdfTst xf0 xf1 C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5pdffns.f,v 5.9 96/10/20 13:42:03 wkt Exp $ C $Log: 5pdffns.f,v $ c Revision 5.9 96/10/20 13:42:03 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.5 96/06/02 23:19:13 wkt c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/19 14:15:44 wkt c A lot of changes in PDF, PdfPrt, ... c C ==================================================================== C FUNCTION PDF (Iset, Ihadron, Iparton, X, Q, IR) C ------------------------------------------------- C Front-end function for parton distributions of pbar, n, p, D, C, Fe C all tranformed in terms of distribution of proton (PrtPdf) C Also checks to see if pdf < 0; if so issue warning and set pdf=0. C Revised 4/4/94 by HLL & WKT: C PDF retains its name and argument list for compatibility with all C existing programs which Call this function; PDF switches between C target hadrons; C PdfPrt is for proton target; it switches between different Iset's. C C Ihadron = -1, 0, 1, 2, 4, 6 : C pbar, n, p, D, C, Fe C 5 is "isoscalar-corrected iron" hence = D C In all cases, adjust the Iparton label C to convert to the corresponding proton distribution which is C given in Function PdfPrt; C --------------------------------------------------- C If any of the arguments are unphysical, it issues a warning. C Return error code IR: C 0 : O.K. C 1,2: Error from call to PdfPrt (proton parton distribution) C 3 : X < 0 or X > 1 : unphysical: warning + set PDF = 0; C 4 : Ihadron out of range; C 5 : PDF < 0., unphysical: warning + set PDF = 0; C 6 : Iparton out of range; warning + set PDF = 0; C ------------------------------------- IMPLICIT DOUBLE PRECISION (A-H, O-Z) Character Msg*80 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx DATA 1 HUGE, DUM, D0m, D1p 1 / 1D10, 0.D0, -1D-6, 1.000001D0 / 1 IW1, IW2 / 2*0 / C Ier = 0 C ------- Check x-range IF (X.LE.D0m .or. X.GT.D1p) Then Ier = 3 Call WARNR(IW1,NWRT,'X out of range in PDF.', 'X', X, > D0, D1, 1) TEM = DUM EndIf C Check bounds on Ihadron If (Ihadron.gt.6 .or. Ihadron.lt.-1 .or. Ihadron.eq.3) then Call WARNI(IW, NWRT, > 'Only Ihardon=-1,0,1,2,4,5,6 (pbar,n,p,D,*,) are active', > 'Ihadron', Ihadron, -1,6,1) Ir = 4 PDF=0.D0 Return Endif C Check bounds on Iparton Jp=Abs(Iparton) Neff = NFL(Q) c nfl(q) returns the number of `light' flavors at scale Q - effective If ( Jp .gt. NEFF) then C if Jp > Neff, then set PDF=0 and return Call WARNI(IW, NWRT, > 'Iparton out of range', > 'Iparton', Iparton, -Neff, Neff,1) Ier = 6 PDF = 0D0 Return Endif C --- Conversion of Ihadron to proton distributions, C if necessary ---- If (Jp.eq.1 .or. Jp.eq.2) Then C For u and d C Use Isospin symmetry n<->p == u<->d Ipartner=3-Jp If (Iparton.lt.0) Ipartner=-Ipartner If (Ihadron.eq.1) then C p: Tem= PdfPrt(Iset, Iparton, X, Q, Ir) Elseif (Ihadron.eq.-1) then C pbar: Tem= PdfPrt(Iset, -Iparton, X, Q, Ir) Elseif (Ihadron.eq.0) then C n: Tem= PdfPrt(Iset, Ipartner, X, Q, Ir) Elseif (Ihadron.eq.2 .or.Ihadron.eq.4 .or.Ihadron.eq.5) then C isoscalar: D, C, ... Tem=( PdfPrt(Iset, Iparton, X, Q, Ir) > +PdfPrt(Iset, Ipartner, X, Q, Ir) )/2.D0 Elseif (Ihadron.eq.6) then C Fe: Tem=( 26.D0* PdfPrt(Iset, Iparton, X, Q, Ir) > +30.D0* PdfPrt(Iset, Ipartner, X, Q, Ir) )/56.D0 Endif Else C For s,c,b,t Tem= PdfPrt(Iset, Iparton, X, Q, Ir) Endif C --- Make sure PDF >= 0 -------- C (unless Iknl<0 - polarized pdf) IF (TEM .LT. D0 .and. Iknl .ge. 1) Then IF (TEM .LT. D0m .AND. X .LE. 0.9D0) Then Call WARNR(IW2,NWRT,'PDF < 0; Set --> 0', 'PDF',TEM,D0,HUGE,1) WRITE (NWRT, '(A, 2I5, 2(1PE15.3))') > ' ISET, Iparton, X, Q = ', ISET, Iparton, X, Q Ier = 5 EndIf TEM = D0 EndIf C -------- Return function value and error code PDF = TEM IR = Ier RETURN C **************************** END FUNCTION PdfPrt (IPDF, LPRTN, XD, QD, IR) C C C This routine gives the parton ( IPARTN ) distribution function inside C the proton in a chosen evolved or parametrized form ( Ipdf ) C 2. For Ipdf = 10 or 11: steer to results from QCD evolution program; C 1 - 9: CTEQ parametrizations C 12 - 20: CTEQ tables or other forms C > 30 : other parametrizations: see codes below for details. C It gives the probability distribution, not the momemtum-weighted one. C C Return Code: IER = 0: No error c < 10: PDF fundtion error C = 10: Ipdf out of range c > 10: multi-error IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) C PARAMETER (MXPDF = 20, MXHDRN = 6, MXF = 6, MXPN = MXF*2+2) C CHARACTER MSG*75 C COMMON / IOUNIT / NIN, NOUT, NWRT C DATA IW1 / 0 / IER = 0 Irr = 0 C Iprtn = LPRTN x = XD Q = QD C ==== Begin of overall Ipdf IfBlock ===== If (Ipdf .eq. 0) Then C ----- Experimental set ------ Tmp = PdfTst (Iprtn, X, Q) ElseIf (Ipdf.ge.1 .and. Ipdf.le.3) then c 1-3 CTEQ3 Nset = Ipdf Tmp = Ctq3Pd(Nset, Iprtn, X, Q, irr) Elseif (Ipdf .le. 9) then c 4-9 CTEQ2 Nset = Ipdf - 3 Tmp = Ctq2Pd(Nset, Iprtn, X, Q, irr) Elseif (Ipdf .eq. 10) then C 10 evolve Tmp = ParDis (Iprtn, X, Q) Elseif (Ipdf .eq. 11) then C 11 evolve Tmp = PrDis1 (Iprtn, X, Q) Elseif (Ipdf.ge.21 .and. Ipdf.le.26) then c 21-26 CTEQ1 Nset = Ipdf-20 tmp = Ctq1Pd(Nset, Iprtn, X, Q, irr) Elseif (Ipdf .EQ. 31) then c 31 EHLQ1 Tmp = PSHK (Iprtn, X, Q, 1) Elseif (Ipdf .EQ. 32) then c 32 Duke-Owens 1 Tmp = PSHK (Iprtn, X, Q, 3) ElseIF (Ipdf .EQ. 33) then c 33 DFLM Tmp = DFLM (Iprtn, X, Q, 2) ElseIF (Ipdf .EQ. 34) then c 34 MT90-S1 Tmp = Pdxmt (1, Iprtn, X, Q, irr) ElseIf (Ipdf.GT.40 .and. Ipdf.LE.50) then c 41- MRS Nset = Ipdf - 40 Tmp = PDMRS (Iprtn, X, Q, Nset) Elseif (Ipdf.gt.60 .and. Ipdf.le.70) then c 61- GRV c 61:LO; 62:NLO(MSbar); 63:NLO(DIS) Nset = Ipdf - 60 Tmp = PDGRV (Iprtn, X, Q, Nset) Elseif (Ipdf.gt.80 .and. Ipdf.le.90) then c 81- CTEQ4 Nset = Ipdf - 80 Tmp = Ctq4Fn (Nset, Iprtn, X, Q) Else Ier = 10 MSG= >'Ipdf chosen is currently inactive. PdfPrt set equal to zero.' CALL WARNI (IW1, NWRT, MSG, 'Ipdf', Ipdf, 1, 10, 0) Tmp = 0. EndIf PdfPrt = Tmp Ir = Ier + Irr 10 RETURN C **************************** END FUNCTION PdfTst (IPRTN, XX, QQ) C Place-holder for a test pdf IMPLICIT DOUBLE PRECISION (A-H, O-Z) Print *, >'You just called a dummy test function PdfTst in PrtPdf!' PdfTst = 0 C RETURN C **************************** END FUNCTION PARDIS (IPRTN, XX, QQ) C C Given the parton distribution function in the array U in C COMMON / PEVLDT / , this routine fetches u(fl, x, q) at any value of C x and q using Mth-order polinomial (II=0) or rational fraction (II=1) C interpolation for x. It always uses quadratic polinomial interpolation C in ln ln (Q/lambda). C The calling program must ensure that 0 =< x =< 1 ; C If 0 =< x < Xmin, extrapolation is used and a warning is given C The calling program must ensure that Alambda < Q ; C If (Alambda < Q < Qini .or. Qmax < Q), C extrapolation is used and a warning is given C IMPLICIT DOUBLE PRECISION (A-H, O-Z) Character Msg*80 LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (Smll = 1D-9) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt C Dimension Fq(5), Df(5) C Dimension QLg(5) C M determines the order of the polynomial C II switches between the polint/ratint interpolation routines C They are fixed for this version Save Data M , II / 2, 0 / Data Iwrn1, Iwarn2, Iwarn3 / 3*0 / C C Check integrity of Upd C If (Irun .Eq. 0) Then C Print '(1pE13.5, 5E13.5)', (UPD(I), I=1,Nelmt-1) C Irun = 1 C Endif X = XX Q = QQ C M = Max (M, 1) C M = Min (M, 3) Md = M / 2 Amd= Md C IF (X .LT. Xmin-Smll) THEN Msg = '0 < X < Xmin in ParDis; extrapolation used!' CALL WARNR (IWRN3, NWRT, Msg, 'X', X, Xmin, 1D0, 1) EndIf C IF (Q .LT. QINI-Smll) THEN Msg = 'Q less than QINI in PARDIS call; Q SET = QINI.' CALL WARNR (IWRN1, NWRT, Msg, 'Q', Q, QINI, QMAX, 1) Q = QINI ElseIF (Q .GT. QMAX) THEN Msg = 'Q greater than QMAX in PARDIS call; ' Msg = Msg // 'Extrapolation will be used.' CALL WARNR(IWRN2, NWRT, Msg, 'Q', Q, QINI, QMAX, 1) EndIf C Find lower end of interval containing X JL = -1 JU = Nx+1 11 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (X .GT. XV(JM)) Then JL = JM Else JU = JM Endif Goto 11 Endif Jx = JL - (M-1)/2 If (Jx .LT. 0) Then Jx = 0 Elseif (Jx .GT. Nx-M) Then Jx = Nx - M Endif C Find the interval where Q lies JL = -1 JU = NT+1 12 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (Q .GT. QV(JM)) Then JL = JM Else JU = JM Endif Goto 12 Endif Jq = JL - (M-1)/2 If (Jq .LT. 0) Then Jq = 0 Elseif (Jq .GT. Nt-M) Then Jq = Nt - M Endif SS = LOG (Q/AL) TT = LOG (SS) C Find the off-set in the linear array Upd JFL = IPRTN + NfMx J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx C C Nt -| .......................... C | .......................... C Jq+M -| .....o ......o............ Iq=M+1 C | .........X................ C Jq -| .(J0)o ......o ........... Iq=1 C | ...... ................... C 0 --------|-------|-----------| C 0 Jx Jx+M Nx C Write (Nwrt, '(/ 10(1pE11.3)/)') X, (XV(I),I=Jx,Jx+M), Q Do 21 Iq = 1, M+1 J1 = J0 + (Nx+1)*(Iq-1) + 1 C Interpolate in x If (II .EQ. 0) Then Call Polint (XV(Jx), Upd(J1), M+1, X, Fq(Iq), Df(Iq)) Elseif (II .EQ. 1) Then Call RatInt (XV(Jx), Upd(J1), M+1, X, Fq(Iq), Df(Iq)) Else Print *, 'II out of range in Pardis; II = ', II Endif C QLg(Iq) = Log ( QV(Jq+Iq-1)/AL ) C Write (Nwrt, '(10(1pE11.3))') C >QV(Jq+Iq-1), (Upd(I), I=J1,J1+M), Fq(Iq), Df(Iq) 21 Continue C Interpolate in LnLnQ Call Polint (TV(Jq), Fq(1), M+1, TT, Ftmp, Ddf) C Or interpolate in LnQ C Call Polint (QLg(1), Fq(1), M+1, SS, Ftmp, Ddf) C Write (Nwrt, '(/ 10(1pE11.3))') (Fq(I), I=1,M+1), Ftmp, Ddf C Write (Nwrt, *) '----------' C PARDIS = Ftmp C RETURN C **************************** END FUNCTION PRDIS1 (IPRTN, XX, QQ) C C Copy of PARDIS for the alternate PDF set C Only the names of common blocks are changed C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C Character Msg*80 COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / X1RRAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / Q1RAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / Q1RAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / E1LPAR / AL, IKNL, IPD0, IHDN, NfMx COMMON / P1VLDT / UPD(MXPQX), KF, Nelmt C Dimension Fq(5), Df(5) C DATA Tiny, SMLL, M, II / 1D-10, 1.D-7, 2, 0 / C X = XX Q = QQ C M = Max (M, 1) C M = Min (M, 3) Md = M / 2 Amd= Md C IF (X .LT. Xmin) THEN If ( (Xmin-X) .gt. Tiny) Then Msg = '0 < X < Xmin in ParDis; extrapolation used!' CALL WARNR (IWRN3, NWRT, Msg, 'X', X, Xmin, 1D0, 1) Endif EndIf C IF (Q .LT. QINI) THEN If ( (QINI-Q) .gt. Tiny) Then Msg = 'Q less than QINI in PRDIS1 call; Q SET = QINI.' CALL WARNR (IWRN1, NWRT, Msg, 'Q', Q, QINI, QMAX, 1) Endif Q = QINI ElseIF (Q .GT. QMAX) THEN If ( (Qmax-Q) .gt. Tiny) Then Msg = 'Q greater than QMAX in PRDIS1 call; ' Msg = Msg // 'Extrapolation will be used.' CALL WARNR(IWRN2, NWRT, MSG, 'Q', Q, QINI, QMAX, 1) Endif EndIf C C Find lower end of interval containing X JL = -1 JU = Nx+1 11 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (X .GT. XV(JM)) Then JL = JM Else JU = JM Endif Goto 11 Endif Jx = JL - (M-1)/2 If (Jx .LT. 0) Then Jx = 0 Elseif (Jx .GT. Nx-M) Then Jx = Nx - M Endif C Find the interval where Q lies JL = -1 JU = NT+1 12 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (Q .GT. QV(JM)) Then JL = JM Else JU = JM Endif Goto 12 Endif Jq = JL - (M-1)/2 If (Jq .LT. 0) Then Jq = 0 Elseif (Jq .GT. Nt-M) Then Jq = Nt - M Endif SS = LOG (Q/AL) TT = LOG (SS) C Find the off-set in the linear array Upd JFL = IPRTN + (KF - 2) / 2 J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx C C Write (Nwrt, '(/ 10(1pE11.3)/)') X, (XV(I),I=Jx,Jx+M), Q Do 21 Iq = 1, M+1 J1 = J0 + (Nx+1)*(Iq-1) + 1 If (II .EQ. 0) Then Call Polint (XV(Jx), Upd(J1), M+1, X, Fq(Iq), Df(Iq)) Elseif (II .EQ. 1) Then Call RatInt (XV(Jx), Upd(J1), M+1, X, Fq(Iq), Df(Iq)) Else Print *, 'II out of range in PRDIS1; II = ', II Endif 21 Continue Call Polint (TV(Jq), Fq(1), M+1, TT, Ftmp, Ddf) PRDIS1 = Ftmp C RETURN C **************************** END FUNCTION XF0 (IX) C The following functions return (X, Q) lattice values C and the values of PDF at these lattice points; C The "0" set is for the standard set of evolved distr.; C The "1" set is for the 'alternate' set. C Value of X at site IX IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt Save DATA IW1, IW2 / 2 * 0 / IF (IX .LE. 0) THEN CALL WARNI(IW1,NWRT,'IX out of range in XF0', 'IX',IX,1,NX,1) XF0 = 0. ElseIF (IX .LE. NX) THEN XF0 = XV(IX) Else CALL WARNI(IW1,NWRT,'IX out of range in XF0', 'IX',IX,1,NX,1) XF0 = 0. EndIf RETURN C **************************** ENTRY QF0 (IQ) C Value of QQ at site IQ IF (IQ .LT. 0) THEN CALL WARNI(IW2,NWRT,'IQ out of range in QF0', 'IQ',IQ,0,NT,1) QF0 = 0. ElseIF (IQ .LE. NT) THEN QF0 = QV(IQ) Else CALL WARNI(IW2,NWRT,'IQ out of range in QF0', 'IQ',IQ,0,NT,1) QF0 = 0. EndIf RETURN C **************************** ENTRY PDF0 (IPRTN, IX, IQ) C Value of PDF at lattice point (IX, IQ) JFL = IPRTN + (KF - 2) / 2 J0 = JFL * (NT+1) * (NX+1) IUPD = J0 + IQ * (NX + 1) + IX + 1 PDF0 = UPD(IUPD) RETURN C **************************** END FUNCTION XF1 (IX) C Value of X at site IX - other set IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX = MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / X1ARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / Q1RAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / E1LPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / P1VLDT / UPD(MXPQX), KF, Nelmt Save DATA IW1, IW2 / 2 * 0 / IF (IX .LE. 0) THEN CALL WARNI(IW1,NWRT,'IX out of range in XF1', 'IX',IX,1,NX,1) XF1 = 0. ElseIF (IX .LE. NX) THEN XF1 = XV(IX) Else CALL WARNI(IW1,NWRT,'IX out of range in XF1', 'IX',IX,1,NX,1) XF1 = 0. EndIf RETURN C **************************** ENTRY QF1 (IQ) C Value of QQ for site IQ - alternate set IF (IQ .LE. 0) THEN CALL WARNI(IW2,NWRT,'IQ out of range in QF1', 'IQ',IQ,1,NT,1) QF1 = 0. ElseIF (IQ .LE. NT) THEN QF1 = XV(IQ) Else CALL WARNI(IW2,NWRT,'IQ out of range in QF1', 'IQ',IQ,1,NT,1) QF1 = 0. EndIf RETURN C **************************** ENTRY PDF1 (IPRTN, IX, IQ) C Value of PDF on the lattice - alternate set JFL = IPRTN + (KF - 2) / 2 J0 = JFL * (NT+1) * (NX+1) IUPD = J0 + IQ * (NX + 1) + IX + 1 PDF1 = UPD(IUPD) RETURN C **************************** END C SUBROUTINE RNEVLU (UinSet, UsIn) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: units used for interactively changing evolution parameters C 5runevl: : rnevlu chgevl prpdf C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5runevl.f,v 5.9 96/10/20 13:42:06 wkt Exp $ C $Log: 5runevl.f,v $ c Revision 5.9 96/10/20 13:42:06 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.7 96/06/03 01:46:09 wkt c setprzpdf removed (it belongs to the prz package) c c Revision 5.6 96/06/02 23:36:17 wkt c minor c c Revision 5.5 96/06/02 23:20:50 wkt c unFormatted --> Formatted c c Revision 5.2 95/07/19 15:17:09 wkt c HLL's setparpdf added c c Revision 5.1 95/07/19 14:17:48 wkt c c Additional argument (initial pdf function) added to PdfSet c C ==================================================================== C SUBROUTINE RNEVLU (UinSet, UsIn) C C An Interactive subprogram to control the running of the Evolution C program package. C C Usin is a dummy function name for a User defined initial parton C distribution function F(Iparton, x) to be used if he/she chooses C to perform QCD evolution based on his/her tailor made initial fn. C UinSet is an (optional) "setup" routine to be called to prepare the C way before UsIn is invoked (e.g. to set up needed coefficients.etc.) C If these functions are used, they must be declared as External in C the calling programs. IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D9=9D0, D6=6D0) CHARACTER FLNM*20, HD*78, msg*120 LOGICAL ASK C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / PARPRZ / XMNPRZ(9), QMNPRZ(9), ALF(9), > NFLPRZ(9), NFAL(9), MORD(9) C External UsIn, UinSet DATA IRUN, NDAT, NU, HD / 0, 55, 47, 'PDF DATA TABLE' / C ------------------------ I1 = 0 I2 = 0 I3 = 0 C Msg= 'Do the data array(S) for the QCD-evolved PDFs already' Msg= Msg//' exist in the regular and alternate common blocks' IF ( ASK(Msg)) RETURN C Msg= 'Do you wish to copy existing data from an external' Msg= Msg//' file to the memory' IF ( ASK(Msg)) THEN C I1 = 1 Msg = 'Input 1 or -1 for the regular or the alternate set.' 17 CALL INPUTI (I1, Msg, -1, 1, IR1) IF (I1 .EQ. 0) GOTO 17 C 1000 CALL INPUTC (FLNM, 'Specify Filename', L) OPEN (NU, FILE=FLNM, FORM='FORMATTED', STATUS='OLD', ERR=333) CALL PDFSET (I1, IH,AL,TM,QI,QM,XM, NU, HD,I2,I3,IR,IRET, UsIn) C IF (IR .EQ. 2) THEN WRITE (NOUT, *) 'Read Error! Fortran Iostat = ', IRET 7 WRITE (NOUT, 9) 9 FORMAT (' What to do? 1. Stop; 2. Proceed, ignore error; > 3. Start over again.' / '$Choose 1, 2, or 3: ') READ (NIN, *) K IF (K .EQ. 1) THEN STOP ElseIF (K .EQ. 3) THEN GOTO 1000 ElseIF (K .NE. 2) THEN WRITE (NOUT, *) 'Illegal choice, (1, 2, 3) only!' GOTO 7 EndIf EndIf C RETURN EndIf C ---------------------------- C Do the evolution I1 = 0 CALL PARPDF (2, 'IPD0', P0, J7) AFL = NFLPRZ (NINT(P0)) CALL INPT2R (P0, AFL, 'Input IPD0 and total NFL', D1, D9, > D1, D6, J4) IPD0 = P0 CALL PARPDF (1, 'IPD0', P0, J6) CALL PARPDF (1, 'XMIN', XMNPRZ(IPD0), J6) CALL PARPDF (1, 'QINI', QMNPRZ(IPD0), J6) CALL PARPDF (1, 'NFL', AFL , J6) If (Ipd0 .Eq. 0) Then Msg = 'For Ipd0=0 to work, you must have already prepared' Msg = Msg//' an input function ' Print *, Msg Msg = 'of (Iparton,X), with FuncName passed as argument to' Msg = Msg//' the RnEvlu routine.' Print *, Msg Msg = 'This function could have a "setup routine" whose' Msg = Msg//' name is also passed the same way.' Print *, Msg If (.not.Ask('Do you have Pdini ready')) > Stop 'Cannot proceed without Pdini; stopping program!' Endif 16 CALL CHGEVL (UinSet) CALL PARPDF (2, 'ORDER', ARDR, J5) CALL PARPDF (2, 'NFL', AFL, J5) CALL PARPDF (2, 'ALAM', ALM, J5) CALL INPT3R(AFL, ALM, ARDR, > 'Input NEFF & LAMDA(NEFF) for ORDER alpha:', > D0, D6, D0, D2, D1, D2, IR) CALL SETLAM (NINT(AFL), ALM, NINT(ARDR)) C CALL PARPDF (2, 'IPD0', P0, J1) CALL PARPDF (2, 'IHDN', HN, J1) CALL PARPDF (2, 'ALAM', AL, J1) CALL PARPDF (2, 'M6 ', TM, J1) CALL PARPDF (2, 'QINI', QI, J1) CALL PARPDF (2, 'QMAX', QM, J1) CALL PARPDF (2, 'XMIN', XM, J1) I2 = P0 IH = HN C CALL PDFSET (I1,IH,AL,TM,QI,QM,XM,NU,HD,I2,0,IR,IRET, UsIn) C IF (IR .EQ. 22) THEN Msg = 'No parameters have changed since last calculation;' Msg = Msg//' Change parameters and try again' IF ( ASK(Msg) ) GOTO 16 IRET = 0 EndIf C Msg = 'Do you wish to write the results to a file and/or' Msg = Msg//' the alternate common blocks' IF (ASK(Msg)) THEN C I1 = 2 1111 CALL INPUTC (FLNM, 'Specify New Filename.', L) OPEN (NU, FILE=FLNM, FORM='FORMATTED', STATUS='NEW', ERR=334) I3 = 1 CALL INPUTC (HD, 'Specify informational Header :', L) IF (ASK('Write data to the alternate common blocks')) I3 = 2 C CALL PDFSET (I1, IH,AL,TM,QI,QM,XM, NU, HD, I2,I3,IR,IRET, Usin) C IF (IR .EQ. 32 .OR. IR .EQ. 33) THEN WRITE (NOUT, *) 'Write, Read Error, Fortran Iostat = ', IRET 8 WRITE (NOUT, 9) READ (NIN, *) K IF (K .EQ. 1) THEN STOP ElseIF (K .EQ. 3) THEN GOTO 1111 ElseIF (K .NE. 2) THEN WRITE (NOUT, *) 'Illegal choice, (1, 2, 3) only!' GOTO 8 EndIf EndIf EndIf C RETURN 333 STOP 'Data table file cannot be openned for read in RNEVLU' 334 STOP 'Data table file cannot be openned for write in RNEVLU' C **************************** END SUBROUTINE CHGEVL (Inptn) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (LENNAM=4, LENREP=10, LPART1=4) CHARACTER NAME*(LENNAM), REPLY*(LENREP) C COMMON / IOUNIT / NIN, NOUT, NWRT C LOGICAL ASK C CALL PARPDF (4, NAME, DBLE(NOUT), IR) 1 WRITE (NOUT, 100) 100 FORMAT('$Type H(help), NC, VAL, OLD Name, > or NAME of parameter to change: ') READ (NIN, 150, END=8) REPLY CALL UPCASE (REPLY) NAME = REPLY (1:LPART1) C Analyze special cases: IF (NAME .EQ. '*') GOTO 1 IF (NAME .EQ. 'NC') GOTO 9 C IF (NAME(1:1) .EQ. 'H') THEN WRITE (NOUT, 101) 101 FORMAT (/' ''HELP'' gives this list and parameter names' > /' ''NC'' means no (more) changes, get on with work' > /' ''VAL'' gives a list of values of all parameters' > /' ''OLD name'' gives value of the named parameter') C WRITE (NOUT, 102) 102 FORMAT (/' List of parameter names:') CALL PARPDF (0, NAME, DBLE(NOUT), IRET) WRITE (NOUT, 103) GOTO 1 EndIf C IF (NAME .EQ. 'VAL') THEN CALL PARPDF (4, NAME, DBLE(NOUT), IRET) GOTO 1 EndIf C IF (NAME .EQ. 'OLD ') THEN NAME = REPLY(LPART1+1:) CALL PARPDF (2, NAME, VALUE, IRET) IF (IRET .EQ. 1) THEN WRITE (NOUT, 104) NAME, VALUE 104 FORMAT (1X, A, ' = ', G15.6) GOTO 1 EndIf C WRITE (NOUT, 105) NAME 105 FORMAT (/1X, A, ' is not a legal name. Try again'/) GOTO 1 EndIf C C If no special action, then change value of variable. C But first check if variable exists: CALL PARPDF (2, NAME, VALUE, IRET) 4 IF (IRET .EQ. 0) THEN WRITE (NOUT, 105) NAME GOTO 1 EndIf C IF (NAME .EQ. 'ISET') THEN IF (ASK ('Do you need to see the Codes for ISET') ) > CALL PRPDF (1, IR) ElseIF (NAME .EQ. 'IHDN') THEN IF (ASK ('Do you need to see the Codes for IHDN') ) > CALL PRPDF (2, IR) EndIf C 2 WRITE (NOUT, 106) NAME, VALUE 106 FORMAT (/ ' Current value of ', A, ' is:', G12.4 / > '$Type new value : ' ) READ (NIN, *, ERR=3, END=8) VALUE C IF (NAME .EQ. 'IPD0' .AND. NINT(VALUE) .EQ. 0) CALL INPTN C CALL PARPDF (1, NAME, VALUE, IRET) IF (IRET .EQ. 0) GOTO 4 IF (IRET .EQ. 1) GOTO 1 C Bad value arrives here: 3 WRITE (NOUT, 107) 107 FORMAT (/' Illegal value. Try again!'/) GOTO 2 C End of file: 8 CONTINUE C Good exit: C CALL XARRAY 9 CONTINUE C RETURN C Formats for input: 103 FORMAT (/) C Formats for output: 150 FORMAT (A) C C **************************** END SUBROUTINE PRPDF (IACT, IRET) C C Prints out menus for: C (1) Parton sets; C (2) Hadron; C (3) Parton label. C Input action code: C Iact = 0 : Print all codes ; C 1 : Print codes of PDF set; C 2 : hadron ID ; C 3 : parton ID ; C C Return Code (Output): C Iret = 0 : Input Iact out of range C 1 : successful try C PARAMETER (MXPDF = 20, MXHDRN = 6, MXF = 6, MXPN = MXF*2+2) C CHARACTER*10 NAMPDF, NMHDRN, NAMQRK C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / NAMPDF / NAMPDF(MXPDF), NMHDRN(-1 : MXHDRN), > NAMQRK(-MXF : MXF) C IRET = 0 C IF (IACT .EQ. 0 .OR. IACT .EQ. 1) THEN WRITE (NOUT, 1101) (I, NAMPDF(I), I = 1, 10) 1101 FORMAT ( ' Code for PDF set labels Ipdf, Ipd0' / > 5(1X, 2(I3, ' :', 5X, A, 10X) /) ) IRET = 1 EndIf C IF (IACT .EQ. 0 .OR. IACT .EQ. 3) THEN WRITE (NOUT, 1103) > (I, NAMQRK(I), -I, NAMQRK(-I), I = 1, 6), NAMQRK(0) 1103 FORMAT ( ' Code for Parton Flavor label Iquark' / > 6(1X, 2(I3, ' :', 5X, A, 10X) /), > 17X, ' 0 :', 5X, A /) IRET = 1 EndIf C IF (IACT .EQ. 0 .OR. IACT .EQ. 2) THEN I1 = -1 WRITE (NOUT, 1102) I1, NMHDRN(-1), (I, NMHDRN(I), I = 1, 5) 1102 FORMAT( ' Code for Hadron label (Ihdrn, etc..): ' / > 3(1X, 2(I3, ' :', 5X, A, 10X) /) ) IRET = 1 EndIf C RETURN C **************************** END C SUBROUTINE EVOLVE (FINI, IRET) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== CC GroupList: Main evolution modules C 5evolve: : evolve nsevl nsrhsm nsrhsp snevl snrhs kernel pff1 stupkl C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5evolve.f,v 5.9 96/10/20 13:41:54 wkt Exp $ C $Log: 5evolve.f,v $ c Revision 5.9 96/10/20 13:41:54 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.6 96/06/03 01:45:09 wkt c bug fix in evolve: missing End statment; Qarray changed. c c Revision 5.5 96/06/02 23:14:17 wkt c Confusing Kf = statement taken fro Qarray to Evolve. c Write and Read to file converted to 'Formatted'. c EvlWt and WtUpd and EvlRd1 modules moved to EvlAux package. c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/19 14:09:25 wkt c Printout of Upd matrix option added; use Ldbg to control print; c Prepare for Version6. c C ==================================================================== C SUBROUTINE EVOLVE (FINI, IRET) C C Input argument: FINI is a function C C FINI (LPARTN, X) C C where LPARTN = -6, ... 6 labels the parton flavor: C t-bar(-6), b-bar(-5), ... gluon(0), u(1), ... t(6) res. C C Output parameter: C C Iret = 0 : normal execution C C ---------------------------- C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / PEVLDT / UPD(MXPQX), Kf, Nelmt C COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) C DIMENSION QRKP(MXF) DIMENSION JI(-MXF : MXF+1) C EXTERNAL NSRHSP, NSRHSM, FINI DATA D0 / 0.0 / C 11 IRET = 0 C Set up number of "valence quarks" IF (IHDN .LE. 4) THEN MXVAL = 2 ElseIF (IHDN .LE. 6) THEN MXVAL = 3 EndIf C Set up X-mesh points and common parameters IF (.NOT. LSTX) CALL XARRAY DLX = 1D0 / NX J10 = NX / 10 C Set up Q-related arrays CALL PARPDF (2, 'ALAM', AL, IR) C Nini and NfMx determined by Qini and Qmax in Qarray CALL QARRAY (NINI) C Flavor # NFSN is for the singlet quark combination NFSN = NFMX + 1 C Kf is the range of the flavor index for the Upd data array C = quark flavors + 1 for singlet combination and 1 for gluon KF = 2 * NFMX + 2 C Total number data points to be stored in Upd Nelmt = Kf * (Nt+1) * (Nx+1) C ---------------------- C Various initiations: C C Offset, or Starting address -1, of Pdf(Iflv) in Upd for the first stage. C Will be updated each turn DO 101 IFLV = -NFMX, NFMX+1 JFL = NFMX + IFLV JI(IFLV) = JFL * (NT+1) * (NX+1) 101 CONTINUE C Define initial distributions for the evolution. C Input functions are defined on (1:Nx). C C Input gluon and quark distributions are inserted C in the Upd array directly; Output of each stage of C evolution serve as the input for the next stage. C 3 DO 31 IZ = 1, NX C Gluon: (Position Ji(0)+1 is reserved for x = 0) C UPD(JI(0)+IZ+1) = FINI (0, XV(IZ)) C Singlet quark C Input singlet quark distribution must be initiated for C each stage separately. C UPD(JI(NFSN)+IZ+1) = 0 C If no quark, bypass filling of quark arrays IF (NFMX .EQ. 0) GOTO 31 C DO 331 IFLV = 1, NINI A = FINI ( IFLV, XV(IZ)) B = FINI (-IFLV, XV(IZ)) QRKP (IFLV) = A + B C Acculumate singlet distribution UPD(JI(NFSN)+IZ+1) = UPD(JI(NFSN)+IZ+1) + QRKP (IFLV) C Initialize the "minus" non-singlet com- C bination (Q - Q-bar) in area for Q-bar UPD(JI(-IFLV)+IZ+1) = A - B 331 CONTINUE C The "plus" non-singlet combination is initialized C in the array area for the Quark of the same flavor DO 332 IFLV = 1, NINI UPD(JI( IFLV)+IZ+1) = QRKP(IFLV) - UPD(JI(NFSN)+IZ+1)/NINI 332 CONTINUE C 31 CONTINUE C ----------------------- C Start of the Q2- Evolution: C C Outer loop is by the effective number of flvr DO 21 NEFF = NINI, NFMX C Set up 1st and 2nd order kernel functions IF (IKNL .EQ. 2) CALL STUPKL (NEFF) C Singlet Calculation ICNT = NEFF - NINI + 1 C Skip if new quark mass = old IF (NTN(ICNT) .EQ. 0) GOTO 21 C Otherwise, recall iteration parameters NITR = NTN (ICNT) DT = DTN (ICNT) TIN = TLN (ICNT) C Perform evolution CALL SNEVL (IKNL, NX, NITR, JT, DT, TIN, NEFF, > UPD(JI(NFSN)+2), UPD(JI(0)+2), UPD(JI(NFSN)+1), UPD(JI(0)+1)) C C Non-singlet sector, one flavor a time. C C Skip this section if quark flavor = 0. IF (NEFF .EQ. 0) GOTO 88 C 5 DO 333 IFLV = 1, NEFF C First evolve the (Q+Q-bar) "PLUS NS" part CALL NSEVL (NSRHSP, IKNL, NX, NITR, JT, DT, TIN, NEFF, > UPD(JI( IFLV)+2), UPD(JI( IFLV)+1)) C C The (Q-Q-bar) "MINUS NS" evolution is needed C only for those flavors with valence IF (IFLV .LE. MXVAL) > CALL NSEVL (NSRHSM, IKNL, NX, NITR, JT, DT, TIN, NEFF, > UPD(JI(-IFLV)+2), UPD(JI(-IFLV)+1)) C C To obtain the real quark distribution functions, C combine the singlet piece to the two non-singlet pieces. C Enforce positivity conditions also at this stage. C DO 55 IS = 0, NITR DO 56 IX = 0, NX TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV)) TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN)) / NEFF TP = TP + TS IF (IKNL .GT. 0) TP = MAX (TP, D0) IF (IFLV .LE. MXVAL) THEN TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV)) IF (IKNL .GT. 0) THEN TM = MAX (TM, D0) TP = MAX (TP, TM) EndIf Else TM = 0. EndIf UPD (JI( IFLV) + IS*(NX+1) + IX + 1) = (TP + TM)/2. UPD (JI(-IFLV) + IS*(NX+1) + IX + 1) = (TP - TM)/2. 56 CONTINUE 55 CONTINUE 333 CONTINUE C Check results C DO 533 IFLV = 0, NEFF, 2 C WRITE (NOUT,'(/A, 2I5)') ' Neff, Iparton =', NEFF, IFLV C WRITE (44, '(/A, 2I5)') ' Neff, Iparton =', NEFF, IFLV C DO 533 IS = 1, NITR, 3 C QQ = QV(NTL(ICNT) + IS) - 0.00001 C WRITE (NOUT,'(/A, I5, 1PE15.3)') ' IQ, Q =', IS, QQ C WRITE (44, '(/A, I5, 1PE15.3)') ' IQ, Q =', IS, QQ C DO 533 IX = 1, NX, J10 C XX = XV(IX) + 0.00001 C TM = UPD (JI( IFLV) + IS*(NX+1) + IX + 1) C TL = DPDF(10, IHDN, IFLV, XX, QQ, JR) C TN = DPDF(IPD0, IHDN, IFLV, XX, QQ, JR) C WRITE (NOUT,'(I5, 4(1PE15.3))') IX, XX, TM, TL, TN C WRITE (44, '(I5, 4(1PE15.3))') IX, XX, TM, TL, TN C 533 CONTINUE C Heavy quarks above current threshold have zero distribution C DO 334 IFLV = NEFF + 1, NFMX DO 57 IS = 0, NITR DO 58 IX = 0, NX UPD(JI( IFLV) + IS*(NX+1) + IX + 1) = 0 UPD(JI(-IFLV) + IS*(NX+1) + IX + 1) = 0 58 CONTINUE 57 CONTINUE 334 CONTINUE 88 CONTINUE C ------------------------------ C Define initial parameters for next stage IF (NFMX .EQ. NEFF) GOTO 21 C New Offsets DO 335 IFLV = -NFMX, NFMX+1 JI(IFLV) = JI(IFLV) + NITR * (NX+1) 335 CONTINUE C New distributions: C gluon input functions are in place; C only non-singlet and singlet quark needs re-initiation C C Calculate initial heavy quark distribution due to C change of renormalization scheme across threshold: CALL HQRK (NX, TT, NEFF+1, UPD(JI(0)+2), UPD(JI(NEFF+1)+2)) C DO 32 IZ = 1, NX QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1) C New Singlet piece UPD (JI(NFSN)+IZ+1) = UPD (JI(NFSN)+IZ+1) + QRKP (NEFF+1) VS00 = UPD (JI(NFSN)+IZ+1) / (NEFF+1) C C "plus" non-singlet for the new flavor UPD(JI( NEFF+1) + IZ + 1) = QRKP(NEFF+1) - VS00 C C Calculate the non-singulet parts of the other quark distr. C C Change from the output of last stage of calcu- C lation due to two sources of change in Vs/Neff: C change of Neff and addition of new quark distr. DO 321 IFL = 1, NEFF A = UPD(JI( IFL)+IZ+1) B = UPD(JI(-IFL)+IZ+1) QRKP(IFL) = A + B C "plus" non-singlet for flavor IFL UPD(JI( IFL)+IZ+1) = QRKP(IFL) - VS00 C "minus" non-singlet for flavors with valence IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1) = A - B 321 CONTINUE C 32 CONTINUE C Return of Q-2 evolution loop to the next stage 21 CONTINUE C Conclusion of the full calculation Return C ********************** End C============================================================================== C Fini is "user defined". Here it is just a dummy place-holder. C============================================================================== C FUNCTION FINI (IPARTN, XX) C C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C FINI = 0. C C RETURN C **************** C END C============================================================================== C SUBROUTINE NSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN) C C IKNL determines to which order in Alpha is the C calculation to be carried out. C C Given the non-singlet parton distribution function U0 at some initial C QIN (Tt= 0) in the x interval (0, 1) covered by the array Iz = 1, NX, C this routine calculates the evoluted function U at Nt values of Tt at C intervals of Dt by numerically integrating the A-P equation using the C non-singlet kernel. C C Un(Ix, Tt) = Y(x,es) at the sites Ix=0,..,Nx (x = Ix * Dz); C Un(0, Tt) is obtained by quadratic extrapolation from Ix = 1, 2, 3 C for each Tt rather then by evolution because of possible C singular behavior at x = 0. C C Data is stored at Tt = Is * Dt, Is = 0, ... , Nt. C The function at Is = 0 is the input distribution. C C Numerical iteration is performed with finer grain Ddt = Dt/Jt. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT Common / PdCntrl/ LPrt, LDbg COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) C DIMENSION U0(NX), UN(0:NX, 0:NT) DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX) external rhs C If (Ldbg .Eq. 1) Then Write (Nwrt, '(A)') ' Non-Singlet:' N5 = Nx / 5 + 1 Endif DDT = DT / JT C IF (NX .GT. MXX) THEN WRITE (NOUT,*) 'Nx =', NX, ' greater than Max # of pts in NSEVL.' STOP 'Program stopped in NSEVL' EndIf C C Compute an effective first order lamda C to be used in checking of moment evl. TMD = TIN + DT * NT / 2. AMU = EXP(TMD) TEM = 6./ (33.- 2.* NEFF) / ALPI(AMU) TLAM = TMD - TEM C C Fill first rows of output array, for initial value of Q DO 9 IX = 1, NX UN(IX, 0) = U0(IX) 9 CONTINUE UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1) C Initiation TT = TIN DO 10 IZ = 1, NX Y0(IZ) = U0(IZ) 10 CONTINUE C loop in the Tt variable DO 20 IS = 1, NT C fine- grained iteration DO 202 JS = 1, JT C Irnd is the counter of the Q-iteration IRND = (IS-1) * JT + JS C Use Runge-Katta the first round IF (IRND .EQ. 1) THEN C CALL RHS (TT, Neff, Y0, F0) C DO 250 IZ = 1, NX Y0(IZ) = Y0(IZ) + DDT * F0(IZ) 250 CONTINUE C TT = TT + DDT C CALL RHS (TT, NEFF, Y0, F1) C DO 251 IZ = 1, NX Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0 251 CONTINUE C What follows is a combination of the 2-step method C plus the Adams Predictor-Corrector Algorithm Else C CALL RHS (TT, NEFF, Y1, F1) C Predictor DO 252 IZ = 1, NX YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0 252 CONTINUE C C Increment of Tt at this place is part of the formalism TT = TT + DDT C CALL RHS (TT, NEFF, YP, FP) C Corrector DO 253 IZ = 1, NX Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0 F0(IZ) = F1(IZ) 253 CONTINUE EndIf C 202 CONTINUE C Fill output array DO 260 IZ = 1, NX UN (IZ, IS) = Y1(IZ) 260 CONTINUE C C The value of the function at x=0 is obtained by extrapolation UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) C Print out for Debugging If (LDbg .Eq. 1) Then Write (Nwrt, '(A, 5(1pE12.3))') ' :', (Un(Iz,Is), Iz=1,Nx,N5) Endif C 20 CONTINUE C RETURN C **************************** END SUBROUTINE NSRHSM (TT, NEFF, FI, FO) C C Subroutine to compute the Right-Side of the Altarelli-Parisi Equation C This copy applies to the "NS-minus" piece -- (Qrk - Qrk-bar) C C See comments in NSRHSP for details IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) PARAMETER (MXQX = MXX*MXQ) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) COMMON / KRNL01 / AFF1 (MXX), AFF2 (MXX), AGG1 (MXX), AGG2 (MXX), > AFG2, AGF2, ANSP (MXX), ANSM (MXX), AQQB COMMON / KRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx DIMENSION G1(MXX), FI(NX), FO(NX), FFI(3*MXX) DIMENSION W0(MXX), W1(MXX), WH(MXX), W2(MXX), W3(MXX) C S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S C CALL NEWARRAY (NX, FI, FFI) CALL INTEGR (NX, 0, FI, W0, IR1) CALL NINTEGR (NX, 1, FFI, FI, W1, IR11) CALL NINTEGR (NX, 8, FFI, FI, W2, IR01) CALL NINTEGR (NX, 7, FFI, FI, W3, IR22) CALL HINTEG (NX, FI, WH) C DO 230 IZ = 1, NX FO(IZ) = 2.* FI(IZ) + 4./3.* (2D0*WH(IZ) - W0(IZ) > - W1(IZ) - W3(IZ) + W2(IZ)) FO(IZ) = CPL * FO(IZ) 230 CONTINUE C 2nd order calculation IF (IKNL .EQ. 2) THEN DZ = 1./ (NX - 1) DO 21 IX = 1, NX-1 X = XV(IX) C Number of points in the I-th integral NP = NX - IX + 1 IS = NP C Evaluate the integrand for the I-th integral DO 31 KZ = 2, NP IY = IX + KZ - 1 IT = NX - IY + 1 C XY = X / Y, already calculated in XARRAY XY = ZZ (IS, IT) G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX)) 31 CONTINUE C TEM1 = SMPNOL (NP, DZ, G1, ERR) C 2nd order contribution TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2 C C TWONE = TMP2 / FO(IX) C IF (TWONE .GE. 0.2) THEN C PRINT '(A, 4(1PE15.4))', C > 'Second order F big in NSRHS: X, Q, 2nd, 1st = ',X,Q,TMP2,FO(I C EndIf C 1st + 2nd order terms FO(IX) = FO(IX) + TMP2 C 21 CONTINUE C EndIf RETURN C **************************** END SUBROUTINE NSRHSP (TT, NEFF, FI, FO) C C Subroutine to compute the Right-Side of the Altarelli-Parisi Equation C This copy applies to the "NS-plus" piece involving (Qrk + Qrk-bar) C C IKNL = 1, 2 : 1st & 2nd order evolution of the unpolarized case; C -1, -2 : ..... polarized .. C C Nx is the number of mesh-points, Tt is the Log Q variable. C Y is the input distribution function defined on the mesh points C F is the RHS value (which is also = dY/dt) defined on the same mesh C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) PARAMETER (MXQX = MXX*MXQ) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) COMMON / KRNL01 / AFF1 (MXX), AFF2 (MXX), AGG1 (MXX), AGG2 (MXX), > AFG2, AGF2, ANSP (MXX), ANSM (MXX), AQQB COMMON / KRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx DIMENSION G1(MXX), FI(NX), FO(NX), FFI(3*MXX) DIMENSION W0(MXX), W1(MXX), WH(MXX), W2(MXX), W3(MXX) C S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S C C Print *, 'r', FI(2) CALL NEWARRAY (NX, FI, FFI) CALL INTEGR (NX, 0, FI, W0, IR1) CALL NINTEGR (NX, 1, FFI, FI, W1, IR11) CALL NINTEGR (NX, 8, FFI, FI, W2, IR01) CALL NINTEGR (NX, 7, FFI, FI, W3, IR22) CALL HINTEG (NX, FI, WH) C DO 230 IZ = 1, NX FO(IZ) = 2.* FI(IZ) + 4./3.* (2D0*WH(IZ) - W0(IZ) > - W1(IZ) - W3(IZ) + W2(IZ)) FO(IZ) = CPL * FO(IZ) 230 CONTINUE C 2nd order calculation IF (IKNL .EQ. 2) THEN DZ = 1./ (NX - 1) DO 21 IX = 1, NX-1 X = XV(IX) C Number of points in the I-th integral NP = NX - IX + 1 C Evaluate the integrand for the I-th integral DO 31 KZ = 2, NP IY = IX + KZ - 1 C XY = X / Y, already calculated in XARRAY XY = ZZ (NX-IX+1, NX-IY+1) G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX)) 31 CONTINUE C TEM1 = SMPNOL (NP, DZ, G1, ERR) C Assemble 2nd order contribution TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + AQQB)) * CPL2 C Sum 1st + 2nd order terms FO(IX) = FO(IX) + TMP2 C 21 CONTINUE C EndIf RETURN C **************************** END SUBROUTINE SNEVL(IKNL, NX, NT, JT, DT, TIN, NEFF, UI, GI, US, GS) C C This is the singlet counter-part of the NSEVL subroutine. Refer to C comments at the beginning of that program section. C C Input parton distributions are Gi (for gluon) and Ui (for singlet quark) C at Tt = 0; output distributions are Gs and Us C at Tt = IS*dt with IS = 1, 2, ... , Nt. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXQX= MXQ * MXX) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT Common / PdCntrl/ LPrt, LDbg COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) C DIMENSION UI(NX), US(0:NX, 0:NT) DIMENSION GI(NX), GS(0:NX, 0:NT) DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX) DIMENSION Z0(MXX), Z1(MXX), ZP(MXX), G0(MXX), G1(MXX), GP(MXX) DATA D0 / 0.0 / If (Ldbg .Eq. 1) Then Write (Nwrt, '(A)') 'Singlet:' N5 = Nx / 5 + 1 Endif C Faster evolution of the singlet sector C requires finer iteration to achieve the C same accuracy as in the non-singlet. JTT = 2 * JT DDT = DT / JTT C IF (NX .GT. MXX) THEN WRITE (NOUT,*) 'Nx =', NX, ' greater than Max # of pts in SNEVL.' STOP 'Program stopped in SNEVL' EndIf C Compute an effective first order lamda C to be used in checking of moment evl. TMD = TIN + DT * NT / 2. AMU = EXP(TMD) TEM = 6./ (33.- 2.* NEFF) / ALPI(AMU) TLAM = TMD - TEM C Initialization (see previous subroutine) DO 9 IX = 1, NX US (IX, 0) = UI(IX) GS (IX, 0) = GI(IX) 9 CONTINUE US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3) GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3) C TT = TIN DO 10 IZ = 1, NX Y0(IZ) = UI(IZ) Z0(IZ) = GI(IZ) 10 CONTINUE C loop in the Tt variable DO 20 IS = 1, NT C fine- grained iteration DO 202 JS = 1, JTT C Irnd is the counter for Q-iterations IRND = (IS-1) * JTT + JS C Use Runge-Katta the first round IF (IRND .EQ. 1) THEN C CALL SNRHS (TT, NEFF, Y0,Z0, F0,G0) C DO 250 IZ = 1, NX Y0(IZ) = Y0(IZ) + DDT * F0(IZ) Z0(IZ) = Z0(IZ) + DDT * G0(IZ) 250 CONTINUE C TT = TT + DDT C CALL SNRHS (TT, NEFF, Y0, Z0, F1, G1) C DO 251 IZ = 1, NX Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0 Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0 251 CONTINUE C What follows is a combination of the 2-step method C and the Adams Predictor-Corrector Algorithm Else C CALL SNRHS (TT, NEFF, Y1, Z1, F1, G1) C Predictor DO 252 IZ = 1, NX YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0 ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0 252 CONTINUE C Increment of Tt at this place is part of the formalism TT = TT + DDT C CALL SNRHS (TT, NEFF, YP, ZP, FP, GP) C Corrector DO 253 IZ = 1, NX Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0 Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0 F0(IZ) = F1(IZ) G0(IZ) = G1(IZ) 253 CONTINUE EndIf C 202 CONTINUE C Fill output array and restore factor of X, if necessary C For spin-averaged case, enforce positivity DO 260 IX = 1, NX IF (IKNL .GT. 0) THEN US (IX, IS) = MAX(Y1(IX), D0) GS (IX, IS) = MAX(Z1(IX), D0) Else US (IX, IS) = Y1(IX) GS (IX, IS) = Z1(IX) EndIf 260 CONTINUE C C The value of the function at x=0 is obtained by extrapolation C US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3) C C Print out for Debugging If (LDbg .Eq. 1) Then Write (Nwrt, '(A, 5(1pE12.3))') ' SQ:',(Us(Iz,Is), Iz=1,Nx,N5) Write (Nwrt, '(A, 5(1pE12.3))') ' G:',(Gs(Iz,Is), Iz=1,Nx,N5) Endif 20 CONTINUE C RETURN C **************************** END SUBROUTINE SNRHS (TT, NEFF, FI, GI, FO, GO) C C Subroutine to compute the Right-Side of the Altarelli-Parisi Equation C for the Singlet sector: C See comments in NSRHSP for notes on IKNL C C FI, Z are the input distributions for quark and gluon respectively; C FO, G are the output dY/dt, dZ/dt. C Nx is the number of mesh-points, Tt is the Log Q variable. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) PARAMETER (MXQX = MXX*MXQ) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) COMMON / KRNL01 / AFF1 (MXX), AFF2 (MXX), AGG1 (MXX), AGG2 (MXX), > AFG2, AGF2, ANSP (MXX), ANSM (MXX), AQQB COMMON / KRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX) DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), R4(MXX) DIMENSION R0(MXX), R1(MXX), RH(MXX),W2(MXX), R3(MXX) DIMENSION FII(3*MXX), GII(3*MXX), R5(MXX), R6(MXX), RH1(MXX) DIMENSION WH1(MXX), WH2(MXX) C S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S C CALL NEWARRAY (NX, FI, FII) CALL NEWARRAY (NX, GI, GII) C CALL NINTEGR (NX, 1, FII, FI, W1, IR11) CALL NINTEGR (NX, 6, FII, FI, W2, IR01) CALL NINTEGR (NX, 8, FII, FI, WH1, IR22) CALL NINTEGR (NX, 7, FII, FI, WH2, IR21) C CALL INTEGR (NX, 0, FI, W0, IR2) CALL INTEGR (NX, 0, GI, R0, IR5) CALL INTEGR (NX, 1, GI, R1, IR6) C CALL NINTEGR (NX, 2, GII, GI, R3, IR13) CALL NINTEGR (NX, 3, GII, GI, R4, IR23) CALL NINTEGR (NX, 4, GII, GI, R5, IR33) CALL NINTEGR (NX, 5, GII, GI, R6, IR42) CALL NINTEGR (NX, 8, GII, GI, RH1, IR43) C dudu CALL HINTEG (NX, FI, WH) CALL HINTEG (NX, GI, RH) C IF (IKNL .GT. 0) THEN DO 230 IZ = 1, NX FO(IZ) = 2D0*FI(IZ) + 4D0/3D0*(2D0*WH(IZ)-W0(IZ) > + WH1(IZ) - WH2(IZ) - W1(IZ)) > + NEFF * (R3(IZ) + R4(IZ)) FO(IZ) = FO(IZ) * CPL C Print *, ' r', FO(IZ) C GO(IZ) = 4D0 / 3D0 * (W0(IZ) + W2(IZ)) > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ) + 3D0 > * (2D0*RH(IZ) > - R0(IZ) - R1(IZ) - 2D0*R4(IZ) > + RH1(IZ) + R6(IZ) + 2D0*R5(IZ)) GO(IZ) = GO(IZ) * CPL C Print *, ' R', GO(IZ) C 230 CONTINUE Else DO 240 IZ = 1, NX FO(IZ) = NEFF * (-R0(IZ) + 2.* R1(IZ) ) > + 2.* FI(IZ) + 4./ 3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ) ) FO(IZ) = FO(IZ) * CPL C GO(IZ) = 4./ 3.* ( 2.* W0(IZ) - W1(IZ) ) >+ (33.- 2.* NEFF) / 6.* GI(IZ) + 6.*(RH(IZ) + R0(IZ) - 2.* R1(IZ)) GO(IZ) = GO(IZ) * CPL 240 CONTINUE EndIf C IF (IKNL .EQ. 2) THEN DZ = 1./(NX - 1) DO 21 I = 1, NX-1 X = XV(I) C Number of points in the I-th integral NP = NX - I + 1 IS = NP C Evaluate the integrand for the I-th integral G2(1) = FFG(IS, IS) * GI(I) G3(1) = GGF(IS, IS) * FI(I) DO 31 KZ = 2, NP IY = I + KZ - 1 IT = NX - IY + 1 C Quantities associated with kernel functions XY = ZZ (IS, IT) C 2nd order terms G1(KZ) = FFG(I, IY) * (FI(IY) - XY**2 *FI(I)) G4(KZ) = GGF(I, IY) * (GI(IY) - XY**2 *GI(I)) G2(KZ) = FFG(IS, IT) * GI(IY) G3(KZ) = GGF(IS, IT) * FI(IY) 31 CONTINUE C TEM1 = SMPNOL (NP, DZ, G1, ERR) TEM2 = SMPSNA (NP, DZ, G2, ERR) TEM3 = SMPSNA (NP, DZ, G3, ERR) TEM4 = SMPNOL (NP, DZ, G4, ERR) C PRINT '(A, F12.3/ A, 4(1PE15.3))', ' Q =', Q, ' Int', C > TEM1, TEM4, TEM2, TEM3 TEM1 = TEM1 - FI(I) * (AFF2(I) + AGF2) TEM4 = TEM4 - GI(I) * (AGG2(I) + AFG2) C PRINT '(A, 2(1PE15.3))', ' Dlx', TEM1, TEM4 TMF = TEM1 + TEM2 TMG = TEM3 + TEM4 FO(I) = FO(I) + TMF * CPL2 GO(I) = GO(I) + TMG * CPL2 C 21 CONTINUE EndIf RETURN C ************************* END SUBROUTINE KERNEL >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT) C New version with 'regularized' kernel functions which are smooth, hence C are more suitable for interpolation. C C Subroutine to calculate the values of the 1st and 2nd order evolution C kernel function at a given value of X. C Formulas used are from Furmanski & Petronzio, Phys.Lett. 97B, p437. C Notations are different than conventional. FG <--> GF, NFL * wrong place C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (PI = 3.141592653589793, PI2 = PI ** 2) PARAMETER (D0 = 0.0, D1 = 1.0) COMMON / IOUNIT / NIN, NOUT, NWRT C DATA CF, CG, TR, IWRN / 1.333333333333, 3.0, 0.5, 0 / C IRT = 0 TRNF = TR * NFL C X = XX IF (X .LE. 0. .OR. X .GE. 1.) THEN CALL WARNR(IWRN, NWRT, 'X out of range in KERNEL', 'X', X, > D0, D1, 1) IRT = 1 RETURN EndIf C XI = 1./ X X2 = X ** 2 XM1I = 1./ (1.- X) XP1I = 1./ (1.+ X) XLN = LOG (X) XLN2 = XLN ** 2 XLN1M = LOG (1.- X) SPEN2 = SPENC2 (X) C 1st order fns + x FFP = (1.+ X2) * XM1I FGP = (2.- 2.* X + X2) / X GFP = 1. - 2.* X + 2.* X2 GGP = XM1I + XI - 2. + X - X2 C 1st order fns - x FFM = (1.+ X2) * XP1I FGM = - (2.+ 2.* X + X2) / X GFM = 1. + 2.* X + 2.* X2 GGM = XP1I - XI - 2. - X - X2 C 1st order kernels with regularization factors, C and excluding subtractions FF1 = CF * FFP * (1.- X) FG1 = CF * FGP * X GF1 = 2.* TRNF * GFP GG1 = 2.* CG * GGP * X * (1.-X) C Begin 2nd order calculations C with the non-singlet pieces PCF2 = -2.* FFP *XLN*XLN1M - (3.*XM1I + 2.*X)*XLN > - (1.+X)/2.*XLN2 - 5.*(1.-X) PCFG = FFP * (XLN2 + 11.*XLN/3.+ 67./9.- PI**2 / 3.) > + 2.*(1.+X) * XLN + 40.* (1.-X) / 3. PCFT = (FFP * (- XLN - 5./3.) - 2.*(1.-X)) * 2./ 3. PQQB = 2.* FFM * SPEN2 + 2.*(1.+X)*XLN + 4.*(1.-X) PQQB = (CF**2-CF*CG/2.) * PQQB PQQ2 = CF**2 * PCF2 + CF*CG * PCFG / 2. + CF*TRNF * PCFT C "Regularize" PNSP = (PQQ2 + PQQB) * (1.-X) PNSM = (PQQ2 - PQQB) * (1.-X) C Now the singlet kernel matrix FFCF2 = - 1. + X + (1.- 3.*X) * XLN / 2. - (1.+ X) * XLN2 / 2. > - FFP * (3.* XLN / 2. + 2.* XLN * XLN1M) > + FFM * 2.* SPEN2 FFCFG = 14./3.* (1.-X) > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.) > - FFM * SPEN2 FFCFT = - 16./3. + 40./3.* X + (10.* X + 16./3.* X2 + 2.) * XLN > - 112./9.* X2 + 40./9./X - 2.* (1.+ X) * XLN2 > - FFP * (10./9. + 2./3. * XLN) C FGCF2 = - 5./2.- 7./2.* X + (2.+ 7./2.* X) * XLN + (X/2.-1.)*XLN2 > - 2.* X * XLN1M > - FGP * (3.* XLN1M + XLN1M ** 2) FGCFG = 28./9. + 65./18.* X + 44./9. * X2 - (12.+ 5.*X + 8./3.*X2) > * XLN + (4.+ X) * XLN2 + 2.* X * XLN1M > + FGP * (-2.*XLN*XLN1M + XLN2/2. + 11./3.*XLN1M + XLN1M**2 > - PI2/6. + 0.5) > + FGM * SPEN2 FGCFT = -4./3.* X - FGP * (20./9.+ 4./3.*XLN1M) C GFCFT = 4.- 9.*X + (-1.+ 4.*X)*XLN + (-1.+ 2.*X)*XLN2 + 4.*XLN1M > + GFP * (-4.*XLN*XLN1M + 4.*XLN + 2.*XLN2 - 4.*XLN1M > + 2.*XLN1M**2 - 2./3.* PI2 + 10.) GFCGT = 182./9.+ 14./9.*X + 40./9./X + (136./3.*X - 38./3.)*XLN > - 4.*XLN1M - (2.+ 8.*X)*XLN2 > + GFP * (-XLN2 + 44./3.*XLN - 2.*XLN1M**2 + 4.*XLN1M > + PI2/3. - 218./9.) > + GFM * 2. * SPEN2 C GGCFT = -16.+ 8.*X + 20./3.*X2 + 4./3./X + (-6.-10.*X)*XLN > - 2.* (1.+ X) * XLN2 GGCGT = 2.- 2.*X + 26./9.*X2 - 26./9./X - 4./3.*(1.+X)*XLN > - GGP * 20./9. GGCG2 = 27./2.*(1.-X) + 67./9.*(X2-XI) + 4.*(1.+X)*XLN2 > + (-25.+ 11.*X - 44.*X2)/3.*XLN > + GGP * (67./9.- 4.*XLN*XLN1M + XLN2 - PI2/3.) > + GGM * 2.* SPEN2 C 2nd order singlet kernel functions FF2 = CF * TRNF * FFCFT + CF ** 2 * FFCF2 + CF * CG * FFCFG FG2 = CF * TRNF * FGCFT + CF ** 2 * FGCF2 + CF * CG * FGCFG GF2 = CF * TRNF * GFCFT + CG * TRNF * GFCGT GG2 = CF * TRNF * GGCFT + CG ** 2 * GGCG2 + CG * TRNF * GGCGT C C Take out singular factors at both ends XLG = (LOG(1./(1.-X)) + 1.) XG2 = XLG ** 2 FF2 = FF2 * X * (1.- X) FG2 = FG2 * X / XG2 GF2 = GF2 * X / XG2 GG2 = GG2 * X * (1.- X) RETURN C **************************** END FUNCTION PFF1 (XX) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LA, LB, LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) PARAMETER (MX = 3) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / KRNL00 / DZ, XL(MX), NNX COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / KRN1ST / FF1(0:MXX), FG1(0:MXX), GF1(0:MXX), GG1(0:MXX), > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX), > PNSP(0:MXX), PNSM(0:MXX) SAVE DATA LA, LB / 2 * .FALSE. / C LB = .TRUE. ENTRY TFF1(ZZ) LA = .TRUE. GOTO 2 ENTRY QFF1 (XX) LB = .TRUE. ENTRY RFF1 (XX) C 2 IF (LA .AND. .NOT.LB) THEN Z = ZZ X = XFRMZ (Z) Else X = XX EndIf C If Xmin < X < 1. interpolate in Z (equally spaced) C If 0 < X < Xmin extrapolate in X (finite range) IF (X .GE. D1) THEN PFF1 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, FF1(1), MX, X, TEM, ERR) EndIf C IF (LA) THEN IF (LB) THEN PFF1 = TEM / (1.-X) LB =.FALSE. Else TFF1 = TEM / (1.-X) * DXDZ(Z) EndIf LA =.FALSE. Else IF (LB) THEN QFF1 = TEM LB =.FALSE. Else RFF1 = TEM * X / (1.-X) EndIf EndIf RETURN C ---------------------------- ENTRY FNSP (XX) X = XX IF (X .GE. D1) THEN FNSP = 0. RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, PNSP(1), MX, X, TEM, ERR) EndIf FNSP = TEM / (1.- X) RETURN C ---------------------------- ENTRY FNSM (XX) X = XX IF (X .GE. D1) THEN FNSM = 0. RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, PNSM(1), MX, X, TEM, ERR) EndIf FNSM = TEM / (1.- X) RETURN C ---------------------------- ENTRY PFG1 (XX) LA = .TRUE. ENTRY QFG1 (XX) LB = .TRUE. ENTRY RFG1 (XX) X = XX IF (X .GE. D1) THEN PFG1 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (FG1, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, FG1(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PFG1 = TEM LA =.FALSE. Else IF (LB) THEN QFG1 = TEM * (1.- X) LB =.FALSE. Else RFG1 = TEM * X EndIf EndIf RETURN C ---------------------------- ENTRY PGF1 (XX) LA = .TRUE. ENTRY QGF1 (XX) LB = .TRUE. ENTRY RGF1 (XX) X = XX IF (X .GE. D1) THEN PGF1= 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (GF1, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, GF1(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PGF1 = TEM / X LA =.FALSE. Else IF (LB) THEN QGF1 = TEM * (1.- X) / X LB =.FALSE. Else RGF1 = TEM EndIf EndIf RETURN C ---------------------------- ENTRY PGG1 (XX) LA = .TRUE. ENTRY QGG1 (XX) LB = .TRUE. ENTRY RGG1 (XX) X = XX IF (X .GE. D1) THEN PGG1= 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, GG1(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PGG1 = TEM / X / (1.-X) LA =.FALSE. Else IF (LB) THEN QGG1 = TEM / X LB =.FALSE. Else RGG1 = TEM / (1.-X) EndIf EndIf RETURN C ---------------------------- ENTRY PFF2 (XX) LA = .TRUE. ENTRY QFF2 (XX) LB = .TRUE. ENTRY RFF2 (XX) X = XX C XLG = LOG(X/(1.-X))**2 IF (X .GE. D1) THEN PFF2 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, FF2(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PFF2 = TEM / X / (1.-X) LA =.FALSE. Else IF (LB) THEN QFF2 = TEM / X LB =.FALSE. Else RFF2 = TEM / (1.-X) EndIf EndIf RETURN C ---------------------------- ENTRY PFG2 (XX) LA = .TRUE. ENTRY QFG2 (XX) LB = .TRUE. ENTRY RFG2 (XX) X = XX XM1 = 1.- X XLG = (LOG(1./XM1) + 1.)**2 IF (X .GE. D1) THEN PFG2 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (FG2, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, FG2(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PFG2 = TEM / X * XLG LA =.FALSE. Else IF (LB) THEN QFG2 = TEM / X * XLG * XM1 LB =.FALSE. Else RFG2 = TEM * XLG EndIf EndIf RETURN C ---------------------------- ENTRY PGF2 (XX) LA = .TRUE. ENTRY QGF2 (XX) LB = .TRUE. ENTRY RGF2 (XX) X = XX XM1 = 1.- X XLG = (LOG(1./XM1) + 1.) ** 2 IF (X .GE. D1) THEN PGF2 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (GF2, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, GF2(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PGF2 = TEM / X * XLG LA =.FALSE. Else IF (LB) THEN QGF2 = TEM / X * XLG * XM1 LB =.FALSE. Else RGF2 = TEM * XLG EndIf EndIf RETURN C ---------------------------- ENTRY PGG2 (XX) LA = .TRUE. ENTRY QGG2 (XX) LB = .TRUE. ENTRY RGG2 (XX) X = XX C XLG = LOG(X/(1.-X)) ** 2 IF (X .GE. D1) THEN PGG2 = 0 RETURN ElseIF (X .GE. XMIN) THEN Z = ZFRMX (X) TEM = FINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT) Else CALL POLINT (XL, GG2(1), MX, X, TEM, ERR) EndIf IF (LA) THEN PGG2 = TEM / X / (1.-X) LA =.FALSE. Else IF (LB) THEN QGG2 = TEM / X LB =.FALSE. Else RGG2 = TEM / (1.-X) EndIf EndIf RETURN C ---------------------------- ENTRY VFF1 (XX) X = XX TEM = (1.+ X**2) / (1.- X) VFF1 = TEM *4./3. RETURN C **************************** END SUBROUTINE STUPKL (NFL) C Set up the common block containing the arrays C for the first and second order kernels C C Also calculates the integrals and constants C needed to complete the [p(x)]sub+ integrals C Real calculation in done in the routine KERNEL which C follows strictly the Furmanski-Petronzio notation. C This routine converts their convention to my convention. C FG (mine) = GF (their) C GF (mine) = FG (their) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MX = 3) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / KRN1ST / FF1(0:MXX), FG1(0:MXX), GF1(0:MXX), GG1(0:MXX), > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX), > PNSP(0:MXX), PNSM(0:MXX) COMMON / KRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) COMMON / KRNL00 / DZ, XL(MX), NNX COMMON / KRNL01 / AFF1 (MXX), AFF2 (MXX), AGG1 (MXX), AGG2 (MXX), > AFG2, AGF2, ANSP (MXX), ANSM (MXX), AQQB EXTERNAL PFF1, RGG1, RFF2, RGG2, RGF2, RFG2 EXTERNAL FNSP, FNSM SAVE DATA AERR, RERR / 0.0, 0.02 / C I = 1, NX corresponds to Z = 0, 1 or X = Xmin, 1 C The point I = 0 corresponding to Z = -DZ. C The point X = 0. is singular, it cannot be reached. C C First past these quantities to FUNCTION PFF1 ... etc C via / KRNL00 / for interpolation and extrapolation purposes NNX = NX DZ = 1./ (NX - 1) DO 5 I0 = 1, MX XL(I0) = XV(I0) 5 CONTINUE C Calculation & switching FG <--> GF DO 10 I = 1, NX-1 XZ = XV(I) CALL KERNEL (XZ, FF1(I), GF1(I), FG1(I), GG1(I), PNSP(I), > PNSM(I), FF2(I), GF2(I), FG2(I), GG2(I), NFL, IRT) 10 CONTINUE C FF1(0) = FF1(1) * 3. - FF1(2) * 3. + FF1(3) FG1(0) = FG1(1) * 3. - FG1(2) * 3. + FG1(3) GF1(0) = GF1(1) * 3. - GF1(2) * 3. + GF1(3) GG1(0) = GG1(1) * 3. - GG1(2) * 3. + GG1(3) PNSP(0) = PNSP(1) * 3. - PNSP(2) * 3. + PNSP(3) PNSM(0) = PNSM(1) * 3. - PNSM(2) * 3. + PNSM(3) FF2(0) = FF2(1) * 3. - FF2(2) * 3. + FF2(3) FG2(0) = FG2(1) * 3. - FG2(2) * 3. + FG2(3) GF2(0) = GF2(1) * 3. - GF2(2) * 3. + GF2(3) GG2(0) = GG2(1) * 3. - GG2(2) * 3. + GG2(3) C FF1(NX) = FF1(NX-1) * 3. - FF1(NX-2) * 3. + FF1(NX-3) FG1(NX) = FG1(NX-1) * 3. - FG1(NX-2) * 3. + FG1(NX-3) GF1(NX) = GF1(NX-1) * 3. - GF1(NX-2) * 3. + GF1(NX-3) GG1(NX) = GG1(NX-1) * 3. - GG1(NX-2) * 3. + GG1(NX-3) PNSM(NX) = PNSM(NX-1) * 3. - PNSM(NX-2) * 3. + PNSM(NX-3) PNSP(NX) = PNSP(NX-1) * 3. - PNSP(NX-2) * 3. + PNSP(NX-3) FF2(NX) = FF2(NX-1) * 3. - FF2(NX-2) * 3. + FF2(NX-3) FG2(NX) = FG2(NX-1) * 3. - FG2(NX-2) * 3. + FG2(NX-3) GF2(NX) = GF2(NX-1) * 3. - GF2(NX-2) * 3. + GF2(NX-3) GG2(NX) = GG2(NX-1) * 3. - GG2(NX-2) * 3. + GG2(NX-3) C Int of PFF1, PNSP, and PNSM C Int of x * Kernel for FF2, GG1, and GG2 RER = RERR * 4. AFF1 (1) = GAUSS(PFF1, D0, XV(1), AERR, RERR, ER1, IRT) DGG1 = NFL / 3. TMPG = GAUSS(RGG1, D0, XV(1), AERR, RERR, ER3, IRT) AGG1 (1) = TMPG + DGG1 ANSM (1) = GAUSS(FNSM, D0, XV(1), AERR, RER, ER2, IRT) ANSP (1) = GAUSS(FNSP, D0, XV(1), AERR, RER, ER2, IRT) AER = AFF1(1) * RER AFF2 (1) = GAUSS(RFF2, D0, XV(1), AER, RER, ER2, IRT) AGF2 = GAUSS(RGF2, D0, XV(1), AER, RER, ER4, IRT) AER = AGG1(1) * RER AGG2 (1) = GAUSS(RGG2, D0, XV(1), AER, RER, ER4, IRT) AFG2 = GAUSS(RFG2, D0, XV(1), AER, RER, ER4, IRT) DO 20 I2 = 2, NX-1 TEM =GAUSS(PFF1,XV(I2-1),XV(I2),AERR,RERR,ER1,IRT) AFF1(I2) = TEM + AFF1(I2-1) AER = ABS(TEM * RER) AGF2 =GAUSS(RGF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT)+AGF2 AFF2(I2)=GAUSS(RFF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT)+AFF2(I2-1) TEM = GAUSS(RGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT) TMPG = TMPG + TEM AGG1(I2) = TMPG + DGG1 AER = ABS(TEM * RER) AFG2 =GAUSS(RFG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)+AFG2 AGG2(I2)=GAUSS(RGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)+AGG2(I2-1) ANSP(I2)=GAUSS(FNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)+ANSP(I2-1) ANSM(I2)=GAUSS(FNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)+ANSM(I2-1) C Compare with exact ans C X = XV(I2) C XLM1 = LOG(1./(1.-X)) C TEMF = (2.*XLM1 - X - X**2/2.) * 4./3. C TEMG = (XLM1 - X**2 + X**3/3. - X**4/4.) * 6. C PRINT '(/I6, 5(1PE14.3))', I2, X, AFF1(I2), TEMF, TMPG, TEMG 20 CONTINUE C Factor relevant to the full PNSP calculation ANSP(NX)=GAUSS(FNSP,XV(NX-1),D1,AERR,RER,ERR, IRT) + ANSP(NX-1) ANSM(NX)=GAUSS(FNSM,XV(NX-1),D1,AERR,RER,ERR, IRT) + ANSM(NX-1) AQQB = ANSP(NX) - ANSM(NX) C Factor relevant to the full FF; restore flavor factor AGF2 = GAUSS(RGF2, XV(NX-1), D1, AER, RER, ERR, IRT) + AGF2 C Factor relevant to the full GG; restore flavor factor AFG2 = GAUSS(RFG2, XV(NX-1), D1, AER, RER, ERR, IRT) + AFG2 C ---------------------------- C Set up kernel functions for direct use in 2nd-order calculations C C 1 . . . . . . . . . 1 . . . . . . . . . C . 1 . . . . . . . . . 1 . . . . . . . . C . . 1 . . . 1 . C . 1 FF2 . . 1 GG2 . C . 1 . . 1 . C . . . . . . C . FG2 1 . . GF2 1 . C . 1 . . 1 . C . . . . . . . . 1 . . . . . . . . . 1 . C . . . . . . . . . 1 . . . . . . . . . 1 C C The same applies to the array PNS containing PNSM \ PNSP DO 21 IX = 1, NX-1 X = XV(IX) C Number of points in the I-th integral NP = NX - IX + 1 IS = NP C Evaluate the FG & GF kernel along diagonal XG2 = (LOG(1./(1.-X)) + 1.) ** 2 FFG (IS, IS) = FG2(NX) * DXTZ(I) * XG2 GGF (IS, IS) = GF2(NX) * DXTZ(I) * XG2 PNS (IS, IS) =PNSM(NX) * DXTZ(I) DO 31 KZ = 2, NP IY = IX + KZ - 1 IT = NX - IY + 1 XY = X / XV(IY) XM1 = 1.- XY XG2 = (LOG(1./XM1) + 1.) ** 2 Z = ZZ (IX, IY) C Quantities needed for interpolation TZ = (Z + DZ) / DZ IZ = TZ IZ = MAX (IZ, 0) IZ = MIN (IZ, NX-1) DT = TZ - IZ C DXTZ is the Jacobian due to the change of variable dx/x -> dZ C 2nd order terms -- by interpolation TEM = (FF2(IZ) * (1.- DT) + FF2(IZ+1) * DT) / XM1 / XY FFG (IX, IY) = TEM * DXTZ(IY) TEM = (FG2(IZ) * (1.- DT) + FG2(IZ+1) * DT) * XG2 / XY FFG (IS, IT) = TEM * DXTZ(IY) TEM = (GF2(IZ) * (1.- DT) + GF2(IZ+1) * DT) * XG2 / XY GGF (IS, IT) = TEM * DXTZ(IY) TEM = (GG2(IZ) * (1.- DT) + GG2(IZ+1) * DT) / XM1 / XY GGF (IX, IY) = TEM * DXTZ(IY) TEM = (PNSP(IZ) * (1.- DT) + PNSP(IZ+1) * DT) / XM1 PNS (IX, IY) = TEM * DXTZ(IY) TEM = (PNSM(IZ) * (1.- DT) + PNSM(IZ+1) * DT) / XM1 PNS (IS, IT) = TEM * DXTZ(IY) 31 CONTINUE 21 CONTINUE RETURN C **************************** END C Subroutine EVLWT (NU, HEADER, IRR) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: Auxillary modules for evolution C 5evlaux: : evlwt evlrd wtupd evlrd1 hqrk fpin fpinms spnint spence C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5evlaux.f,v 5.9 96/10/20 13:26:33 wkt Exp $ C $Log: 5evlaux.f,v $ c Revision 5.9 96/10/20 13:26:33 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.8 96/10/15 23:43:18 wkt c minor bug fixes c c Revision 5.5 96/06/02 23:03:16 wkt c Write and Read to file converted to 'Formatted'. c EvlWt and WtUpd and EvlRd1 modules moved to EvlAux package. c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/19 14:08:19 wkt c comments only c C ==================================================================== C Subroutine EVLWT (NU, HEADER, IRR) C C Before calling this routine, the calling program must open an external C file for data transfer, provide a filename for the file, and C establish the equivalence of that file with unit=Nu, such as: C C Open (Nu, File='FILENAME', Form='FORMATTED', status='NEW', ... ) C C C The file is not automatically close because it is also used for memory C transfer (PEVLDT ---> PEVLD1) which required reading from this unit. C Input parameter: Nu = unit number for the external file to be C written to. (established in calling program) C C Header = informational header for the file C C Output parameter: Irr C C ---------------------------- IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER HEADER*78, Line*80 C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx COMMON / PEVLDT / UPD(MXPQX), KF, Nelmt CALL PARPDF (2, 'ALAM', AL, IR) CALL PARPDF (2, 'NFL', FL, IR) CALL PARPDF (2, 'ORDER', DR, IR) CALL PARPDF (2, 'M1', AM1, IR) CALL PARPDF (2, 'M2', AM2, IR) CALL PARPDF (2, 'M3', AM3, IR) CALL PARPDF (2, 'M4', AM4, IR) CALL PARPDF (2, 'M5', AM5, IR) CALL PARPDF (2, 'M6 ', AM6, IR) C Write (Nu, '(A)') Header Write (Nu, '(2x, A, 8x, A, 9x, A)') $ 'Ordr, Nfl, lambda', 'Qmass 1, 2, 3,', '4, 5, 6' Write (Nu, '(2F5.0, F8.4, 6F9.3)') $ Dr, Fl, Al, Am1, Am2, Am3, AM4, AM5, AM6 Write (Nu, '(3x, A)') 'IPD0, IHDN, IKNL, NfMx, KF, Nelmt' Write (Nu, '(5I6, I8)') IPD0, IHDN, IKNL, NfMx, KF, Nelmt Write (Nu, '(3x,A)') 'NX, NT, JT, NG, NTL(NG+1)' Write (Nu, '(5I5)') NX, NT, JT, NG, NTL(NG+1) Write (Nu, '(A)') '(NTL(I), NTN(I), TLN(I), DTN(I), I =1, NG)' Write (Nu, '(2I5, 1pE13.5, E13.5)') $ (NTL(I), NTN(I), TLN(I), DTN(I), I =1, NG) Write (Nu, '(A)') 'QINI, QMAX, (QV(I), TV(I), I =0, NT)' Write (Nu, '(1pE13.5, E13.5 / (2E13.5))') $ QINI, QMAX, (QV(I), TV(I), I =0, NT) Write (Nu, '(A)') 'XMIN, XCR, (XV(I), I =1, NX)' Write (Nu, '(1pE13.5, E13.5 / (6E13.5))') $ XMIN, XCR, (XV(I), I =1, NX) C C Since quark = anti-quark for nfl>2 at this stage, C we write out only the non-redundent data points C No of flavors = NfMx sea + 1 gluon + 2 valence Write (Nu, '(A)') 'Parton Distribution Table:' Npts = (NX+1) * (NT+1) * (NfMx+3) CALL WTUPD (UPD, Npts, NU, IR2) C IRR = IR2 C IF (IRR .NE. 0) PRINT *, 'Write error in EVLWT' C RETURN C ---------------------------- C Read external file ENTRY EVLRD (NU, HEADER, IRR) C C See comment lines of the Entry EVLWT section for detailed instructions. C C The call program must open the existing file to be read from with: C C Open (Nu, File='FILENAME', Form='UNFORMATTED', status='OLD', Err=... ) C C where the file named 'FILENAME' must exist and it has to be written C originally by EVLWT. (cf. above) c C HEADER is an output character variable. C C ---------------------------- C Read (Nu, '(A)') Header Read (Nu, '(A)') Line Read (Nu, *) Dr, Fl, Al, Am1, Am2, Am3, AM4, AM5, AM6 Read (Nu, '(A)') Line Read (Nu, *) IPD0, IHDN, IKNL, NfMx, KF, Nelmt Read (Nu, '(A)') Line Read (Nu, *) NX, NT, JT, NG, NTL(NG+1) Read (Nu, '(A)') Line Read (Nu, *) (NTL(I), NTN(I), TLN(I), DTN(I), I =1, NG) Read (Nu, '(A)') Line Read (Nu, *) QINI, QMAX, (QV(I), TV(I), I =0, NT) Read (Nu, '(A)') Line Read (Nu, *) XMIN, XCR, (XV(I), I =1, NX) C C Since quark = anti-quark for nfl>2 at this stage, C we Read out only the non-redundent data points C No of flavors = NfMx sea + 1 gluon + 2 valence Nblk = (NX+1) * (NT+1) Npts = Nblk * (NfMx+3) Read (Nu, '(A)') Line CALL RdUpd (UPD, Npts, NU, IRR) CLOSE (NU) IF (IRR .NE. 0) PRINT *, 'Read error in EVLRD' If (NfMx .GE. 3) Then C Refill Upd for s -> t Do 11 Nflv = 3, NfMx J0 = (-Nflv + NfMx) * Nblk C offset = NfMx sea + 1 gluon + (Nflv-1) less-massive quarks = NfMx+Nflv J1 = ( Nflv + NfMx) * Nblk Do 12 I = 1, Nblk Upd (J1 + I) = Upd (J0 + I) 12 Continue 11 Continue EndIf C To check read-in result against saved file. C Print '(A/(1pE13.5, 5E13.5))', 'Upd', (UPD(I), I=1,Nelmt-1) C C Set QCD parameters to current values C CALL PARPDF (1, 'ALAM', AL, IR) C CALL PARPDF (1, 'NFL', FL, IR) C CALL PARPDF (1, 'ORDER', DR, IR) C CALL PARPDF (1, 'M1', AM1, IR) C CALL PARPDF (1, 'M2', AM2, IR) C CALL PARPDF (1, 'M3', AM3, IR) C CALL PARPDF (1, 'M4', AM4, IR) C CALL PARPDF (1, 'M5', AM5, IR) C CALL PARPDF (1, 'M6 ', AM6, IR) C C PRINT '(/A/1X,A/)', ' Parameters from this data file are:',HEADER C CALL PARPDF(4, 'ALL', DBLE(NOUT), IR) C ------------------------------ Return C ***************** End SUBROUTINE WTUPD (UPD, NTL, NDAT, IRET) C C The I/O operation is made a stand-alone subprogram C so that fast execution is achieved with block C data transfer for the actual size of the array C (rather than the declared size in the main program). DOUBLE PRECISION UPD DIMENSION UPD (NTL) C WRITE (NDAT, '(1pE13.5, 5E13.5)', IOSTAT=IRET) UPD C RETURN C ---------------------------- C ENTRY RDUPD (UPD, NTL, NDAT, IRET) C READ (NDAT, *, IOSTAT=IRET) UPD C To check read-in result against saved file. C Print '(1pE13.5, 5E13.5)', UPD C RETURN C **************************** END SUBROUTINE EVLRD1 (NU, HEADER, IRR) C C READ DATA FROM FILE TO ALTERNATE COMMON BLOCKS C C See comment lines of the Entry EVLRD section for detailed instructions. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER HEADER*78 C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / X1RRAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / Q1RAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / Q1RAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / E1LPAR / AL, IKNL, IPD0, IHDN, NfMx COMMON / P1VLDT / UPD(MXPQX), KF, Nelmt C Read (Nu, '(A)') Header Read (Nu, '(A)') Line Read (Nu, *) Dr, Fl, Al, Am1, Am2, Am3, AM4, AM5, AM6 Read (Nu, '(A)') Line Read (Nu, *) IPD0, IHDN, IKNL, NfMx, KF, Nelmt Read (Nu, '(A)') Line Read (Nu, *) NX, NT, JT, NG, NTL(NG+1) Read (Nu, '(A)') Line Read (Nu, *) (NTL(I), NTN(I), TLN(I), DTN(I), I =1, NG) Read (Nu, '(A)') Line Read (Nu, *) QINI, QMAX, (QV(I), TV(I), I =0, NT) Read (Nu, '(A)') Line Read (Nu, *) XMIN, XCR, (XV(I), I =1, NX) C C Since quark = anti-quark for nfl>2 at this stage, C we Read out only the non-redundent data points C No of flavors = NfMx sea + 1 gluon + 2 valence Nblk = (NX+1) * (NT+1) Npts = Nblk * (NfMx+3) CALL RdUpd (UPD, Npts, NU, IRR) CLOSE (NU) IF (IRR .NE. 0) PRINT *, 'Read error in EVLRD1' If (NfMx .GE. 3) Then C Refill Upd for s -> t Do 11 Nflv = 3, NfMx J0 = (-Nflv + NfMx) * Nblk C offset = NfMx sea + 1 gluon + (Nflv-1) less-massive quarks = NfMx+Nflv J1 = ( Nflv + NfMx) * Nblk Do 12 I = 1, Nblk Upd (J1 + I) = Upd (J0 + I) 12 Continue 11 Continue EndIf C PRINT '(/A/1X,A/)',' EVL Parameters for this set are:', HEADER CALL EVLPAR (6, 'ALL', DBLE(NOUT), IR) PRINT '(/ A /)', ' QCD Parameters for this set are:' PRINT '(A, F7.3, 2F7.0)', ' ALAM, NFL, NRDR = ', AL, FL, DR PRINT '(A, 6F7.3)',' M1, M2, ..., M6 = ', AM1,AM2,AM3,AM4,AM5,AM6 C Compare two sets of QCD parameters CALL PARPDF (2, 'ALAM', BL, IR) CALL PARPDF (2, 'ORDER',RDR, IR) CALL PARPDF (2, 'M1', BM1, IR) CALL PARPDF (2, 'M2', BM2, IR) CALL PARPDF (2, 'M3', BM3, IR) CALL PARPDF (2, 'M4', BM4, IR) CALL PARPDF (2, 'M5', BM5, IR) CALL PARPDF (2, 'M6', BM6, IR) CALL PARPDF (2, 'NFL', BFL, IR) LF = NINT(2.*BFL+2.) C SML = 1.E-4 DIF = MAX(ABS(AM1-BM1), ABS(AM2-BM2), ABS(AM3-BM3), ABS(AM4-BM4), > ABS(AM5-BM5), ABS(AM6-BM6), ABS(RDR- DR), ABS(AL-BL)) IF (DIF .GT. SML .OR. KF .NE. LF) THEN PRINT *, 'Warning! Two PDF sets have different QCD parameters.' PRINT '(/A/)', ' Parameters for the regular set are:' CALL PARPDF(4, 'ALL', DBLE(NOUT), IR) EndIf C RETURN C **************************** END SUBROUTINE HQRK (NX, TT, NQRK, Y, F) C C Subroutine to compute the (heavy) quark distribution from (given) C gluon distribution, as the result of a change in the renormalization C scheme (from MS-bar to BPH) as the threshold for quark flavor Nqrk C is crossed. C C Nx is the number of mesh-points, Tt is the Log Q variable. C Y is the input g-distribution function defined on the mesh points C F is the outpur Qrk-distribution function defined on the same pts. C C If the threshold is chosen at the 'natural boundary' Mu = Mass(qrk), C then there is no renormalization and this routine returns zero. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) C DIMENSION Y(NX), F(NX) C DIMENSION W0(MXX), W1(MXX), W2(MXX), WH(MXX), WM(MXX) C C Returns zero, assuming 'natural boundary'. IF (NX .GT. 1) GOTO 11 C Q = EXP(TT) C AL = ALPI(Q) C AMS = AMASS(NQRK) C AMH = AMHATF(NQRK) C FAC = 2.* LOG (AMS / AMH) C FAB = AL / 4. C C CALL INTEGR (NX, 0, Y, W0, IR1) C CALL INTEGR (NX, 1, Y, W1, IR2) C CALL INTEGR (NX, 2, Y, W2, IR3) C 11 CONTINUE DO 230 IZ = 1, NX C Returns zero answer. IF (NX .GT. 1) THEN F(IZ) = 0 GOTO 230 EndIf C C F(IZ) = FAB * ( FAC * W0(IZ) C > + 2.* (FAC -1.) * (-W1(IZ) + W2(IZ) ) ) C F(Iz) = fab * ( Fac * (W0(Iz) - Y(Iz) / 2.) C > + 2.* (fac -1.) * (-W1(Iz) + W2(Iz) + Y(Iz) / 12.)) C 230 CONTINUE C RETURN C **************************** END FUNCTION FPIN (LPARTN, X) C Initial parton distribution at Q = Qini C from published PDF parametrizations IMPLICIT DOUBLE PRECISION (A-H, O-Z) C IR = 0 IR1 = 0 IR2 = 0 IR3 = 0 C CALL PARPDF (2, 'IPD0', SET, IR1) CALL PARPDF (2, 'IHDN', HDRN, IR2) CALL PARPDF (2, 'QINI', QINI, IR3) C IF (IR1 .NE. 1) PRINT *, 'Ir1 = ', IR1, ' in FPIN' IF (IR2 .NE. 1) PRINT *, 'Ir2 = ', IR2, ' in FPIN' IF (IR3 .NE. 1) PRINT *, 'Ir3 = ', IR3, ' in FPIN' C IH = HDRN IS = SET C FPIN = PDF (IS, IH, LPARTN, X, QINI, IR) C IF (IR .NE. 0) PRINT *, 'Ir = ', IR, ' in FPIN' RETURN C **************************** END FUNCTION FPINMS (LPARTN, X) C Initial parton distribution at Q = Qini C from published PDF parametrizations C C Conversion from DIS to MS-bar scheme is C performed before the function is passed on. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C IR = 0 IR1 = 0 IR2 = 0 IR3 = 0 C CALL PARPDF (2, 'IPD0', SET, IR1) CALL PARPDF (2, 'IHDN', HDRN, IR2) CALL PARPDF (2, 'QINI', QINI, IR3) C IF (IR1 .NE. 1) PRINT *, 'Ir1 = ', IR1, ' in FPINMS' IF (IR2 .NE. 1) PRINT *, 'Ir2 = ', IR2, ' in FPINMS' IF (IR3 .NE. 1) PRINT *, 'Ir3 = ', IR3, ' in FPINMS' C IH = HDRN IS = SET C MRS sets are in MS-Bar, the others in DIS scheme to begin with. IF (IS .GE. 5 .AND. IS .LE. 7) THEN IACT = -1 Else IACT = 1 EndIf TEM = PMSDIS (IS, IH, LPARTN, X, QINI, IACT) C FPINMS = TEM C IF (IR .NE. 0) PRINT *, 'Ir = ', IR, ' in FPINMS' RETURN C **************************** END FUNCTION SPNINT (ZZ) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C COMMON / SPENCC / X C Z = ZZ TEM = LOG (1.- Z) / Z SPNINT = TEM C RETURN C ---------------------------- ENTRY SPN2IN (ZZ) C Z = ZZ TEM = LOG (1.+ X - Z) / Z SPN2IN = TEM C RETURN C **************************** END FUNCTION SPENCE (X) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C EXTERNAL SPNINT, SPN2IN C COMMON / SPENCC / XX C DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 / C XX = X TEM = GAUSS(SPNINT, XX, U1, AERR, RERR, ERR, IRT) SPENCE = TEM C RETURN C ---------------------------- ENTRY SPENC2 (X) C XX = X TEM = GAUSS (SPN2IN, XX, U1, AERR, RERR, ERR, IRT) SPENC2 = TEM + LOG (XX) ** 2 / 2. C RETURN C **************************** END C SUBROUTINE QARRAY (NINI) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: routines to manipulate grids and integration for evolution C 5evlutl: :qarray xarray integr hinteg zfrmx C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5evlutl.f,v 5.9 96/10/20 13:41:47 wkt Exp $ C $Log: 5evlutl.f,v $ c Revision 5.9 96/10/20 13:41:47 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.6 96/06/03 01:44:04 wkt c bug fix in Qarray; argument NfMx removed (it is in Common) c c Revision 5.5 96/06/02 23:06:05 wkt c Confusing Kf = statement taken fro Qarray to Evolve. c Commons EVLPAC and PEVLDT modified. c Set Nfl and AL of QCDLIB to Max flavor # and Reset Lambda values in Qarray. c c Revision 5.1 95/07/19 14:06:54 wkt c kf mismatch fixed; XYARAY changed c C ==================================================================== C SUBROUTINE QARRAY (NINI) C C Given Qini, Qmax, and Nt, this routine go through the various C flavor thresholds; determines the step-sizes in-between each pair C of thresholds, computes a new NT, if necessary; and returns: C the variable NINI, NFMX as arguments; C the (NT+1)-dim Q- and T=lnlnQ/Lamda arrays in the common block QARAY1; C the Neff-dependent variable (see below) in QARAY2; and C KF in the common block EvlPac C C Note in particular that the value of NT may increase by 1 or 2 in a C call to this routine because of the logistics of setting up the steps C for the flavor thresholds. C C For given Neff, the variables TLN, DTN, NTL and NTN has the following C meaning: C Flavor thresholds C Nini Nf Nf+1 ....NFmx C |.. ........|<- DTN ->|<- DTN ->| .......................|..........| C <------------ NTL(Nf) steps -------------> C 0 NTN(Nf) NTN(Nf+1) NT C TLN(Nf) (TLN=lnln(Q/Lamda)) TNL(Nf+1) C TV(0).... TV(NTN(Nf)).....TV(NTN(Nf)+2).............TV(NTN(Nf+1))..TV(NT) C QV(0).... QV(NTN(Nf)).....QV(NTN(Nf)+2).............QV(NTN(Nf+1))..QV(NT) C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / QARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx NCNT = 0 IF (NT .GE. mxq) NT = mxq - 1 C Set up overall t=log(Q) parameters S = LOG(QINI/AL) TINI = LOG(S) S = LOG(QMAX/AL) TMAX = LOG(S) C Nt is the provisional # of mesh-points in t; C Dt0 is the approximate increment in T, C Actual values of (Nt, Dt) are determined later. 1 DT0 = (TMAX - TINI) / float(NT) C Determine effective Nflv at Qini and at Qmax NINI = NFL(QINI) NFMX = NFL(QMAX) Call ParQcd (2, 'ORDER', Ord, Ir) Call ParQcd (2, 'ALAM', Al0, Ir) Call ParQcd (2, 'NFL', Afl0, Ir) C Set the total number of Quark flavors C in the QCD package to NfMx AFL = NfMx Call ParQcd (1, 'NFL', AFL, Ir) C and restore the QCD coupling Iordr = Nint (Ord) Ifl0 = Nint (Afl0) Call SetLam (Ifl0, Al0, Iordr) C C Q2 evolution is carried out in stages separated by flavor thresholds C NG = NFMX - NINI + 1 C ----------------------- C Determine the threshold points and C Set up detailed t- mesh structure QIN = QINI QOUT = QINI S = LOG(QIN/AL) TIN = LOG(S) TLN(1) = TIN NTL(1) = 0 QV(0) = QINI TV(0) = Tin C DO 20 NEFF = NINI, NFMX C ICNT = NEFF - NINI + 1 IF (NEFF .LT. NFMX) THEN THRN = AMHATF (NEFF + 1) QOUN = MIN (QMAX, THRN) Else QOUN = QMAX EndIf C IF (QOUN-QOUT .LE. 0.0001) THEN DT = 0 NITR = 0 Else QOUT = QOUN S = LOG(QOUT/AL) TOUT = LOG(S) TEM = TOUT - TIN C Nitr = Number of iterations in this stage NITR = INT (TEM / DT0) + 1 DT = TEM / NITR EndIf C Save book-keeping data on array DTN (ICNT) = DT NTN (ICNT) = NITR TLN (ICNT) = TIN NTL (ICNT+1) = NTL(ICNT) + NITR C C QV is the physical Q-value for this lattice point. IF (NITR .NE. 0) THEN DO 205 I = 1, NITR TV (NTL(ICNT)+I) = TIN + DT * I S = EXP (TV(NTL(ICNT)+I)) QV (NTL(ICNT)+I) = AL * EXP (S) 205 CONTINUE EndIf C Initialize for next iteration QIN = QOUT TIN = TOUT C 20 CONTINUE C Redefine Nt, as the actual number of t-points; may have been C increased by 1 or 2 by the above algorithm C If NT > MXQ, reduce NT until it is within bounds NCNT = NCNT + 1 NTP = NTL (NG + 1) ND = NTP - NT IF (NTP .GE. MXQ) THEN NT = MXQ - ND - NCNT GOTO 1 EndIf NT = NTP C RETURN C **************************** END SUBROUTINE XARRAY C C Given NX, this routine fills the x-arrays of the common blocks C related to the x-variable for the use of various other routines. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C PARAMETER (D0 = 0.0, D10=10.0) PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C Character Msg*80 COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) COMMON / HINTEC / GH(NDG, MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ(MXX, MXX), ZV(0:MXX) C DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG) C DATA F12, F22, F32 / 1D0, 1D0, 1D0 / DATA (G1(I,NDH), G2(I,1), I=1,NDG) / 0.0,0.0,0.0,0.0,0.0,0.0 / DATA PUNY / 1D-30 / C C ---------------------------- C Change of variable X <--> Z C Range of Z is [0, 1] for X in [XM, 1] and I in [1, NX] C Use evenly spaced points in Z XV(0) = 0D0 DZ = 1D0 / (NX-1) DO 10 I = 1, NX - 1 Z = DZ * (I-1) ZV(I) = Z X = XFRMZ (Z) C DXDZ is the Jacobian for the change of variable DXTZ(I) = DXDZ(Z) / X C Fill x - arrays XV (I) = X XA(I, 1) = X XA(I, 0) = LOG (X) DO 20 L = L1, L2 IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L 20 CONTINUE 10 CONTINUE C Fill last point, I=Nx XV(NX) = 1D0 DXTZ(NX) = DXDZ(1.D0) DO 21 L = L1, L2 XA (NX, L) = 1D0 21 CONTINUE XA (NX, 0) = 0D0 C Fill ELY = Log (1. - x) DO 11 I = 1, NX-1 ELY(I) = LOG(1D0 - XV(I)) 11 CONTINUE C Log(1-x) is infinite at x=1, for the purpose of numerical C Calculations, the following extrapolated value is used C to avoid an artificial discontinuity. C ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3) C ---------------------------- C Matrix elements for 2nd order calculations C C 1 . . . . . . . . . C . 1 . . . . . . . . C . . 1 . C . 1 Z (X/Y) . C . 1 . C . . . C . X/Y 1 . C . 1 . C . . . . . . . . 1 . C . . . . . . . . . 1 DO 17 IX = 1, NX ZZ (IX, IX) = 1. DO 17 IY = IX+1, NX XY = XV(IX) / XV(IY) ZZ (IX, IY) = ZFRMX (XY) ZZ (NX-IX+1, NX-IY+1) = XY 17 CONTINUE C ------------------------------ C Start of x - loop to compute ceefficients C for integrals used in INTEG, AMOM, ... etc. DO 30 I = 1, NX-1 C "F" matrix {a(i)} --> {f(i)} IF (I .NE. NX-1) THEN F11 = 1D0/XV(I) F21 = 1D0/XV(I+1) F31 = 1D0/XV(I+2) F13 = XV(I) F23 = XV(I+1) F33 = XV(I+2) C Determinant for matrix C DET = F11*F22*F33 + F21*F32*F13 + F31*F12*F23 > - F31*F22*F13 - F21*F12*F33 - F11*F32*F23 IF (ABS(DET) .LT. PUNY) THEN Msg='Determinant close to zero; will be arbitrarily set to:' CALL WARNR(IWRN, NWRT, Msg, 'DET', PUNY, D0, D0, 0) DET = PUNY EndIf C Compute "G" matrix -- inverse of F C G1 is only needed from I=2 and on G2(1,2) = (F22*F33 - F23*F32) / DET G2(1,3) = (F32*F13 - F33*F12) / DET G2(1,4) = (F12*F23 - F13*F22) / DET C G2(2,2) = (F23*F31 - F21*F33) / DET G2(2,3) = (F33*F11 - F31*F13) / DET G2(2,4) = (F13*F21 - F11*F23) / DET C G2(3,2) = (F21*F32 - F22*F31) / DET G2(3,3) = (F31*F12 - F32*F11) / DET G2(3,4) = (F11*F22 - F12*F21) / DET C Compute coefficients for HINTEG B2 = LOG (XV(I+2)/XV(I)) B3 = XV(I) * (B2 - 1.) + XV(I+2) GH (1,I) = B2 * G2 (2,2) + B3 * G2 (3,2) GH (2,I) = B2 * G2 (2,3) + B3 * G2 (3,3) GH (3,I) = B2 * G2 (2,4) + B3 * G2 (3,4) EndIf C "G-bar" is the "average" of G1 and G2 DO 51 J = 1, NDH DO 52 L = 1, NDG C First interval IF (I .EQ. 1) THEN GB(L,J,I) = G2(L,J) C last interval ElseIF (I .EQ. NX-1) THEN GB(L,J,I) = G1(L,J) C intermidiate intervals Else GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0 EndIf 52 CONTINUE 51 CONTINUE C Compute "A" matrix DO 35 MM = M1, M2 DO 40 K = 1, NDG KK = K + MM - 2 IF (KK .EQ. 0) THEN A(K) = XA(I+1, 0) - XA(I, 0) Else A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK) EndIf 40 CONTINUE C Compute "H" matrix DO 41 J = 1, NDH TEM = 0 DO 43 L = 1, NDG TEM = TEM + A(L) * GB(L,J,I) 43 CONTINUE H(J,I,MM) = TEM 41 CONTINUE 35 CONTINUE C ------------------------------ C Initialize G1 matrix DO 42 J = 1, NDG DO 44 L = 1, NDG G1(L,J) = G2(L,J+1) 44 CONTINUE 42 CONTINUE C 30 CONTINUE C End of x - loop to calculate coefficients C ------------------------------ LSTX = .TRUE. RETURN C **************************** END SUBROUTINE NEWARRAY (NX, F, FF) C C In this routine the input parton distribution f(y) will be approximated C in the bin x_i, x_i+1 by making apolynomial approximation: C f(y)= a*y^2 + b*y +c C where the coefficients a,b,c are determined by knowing f(y) in C the points x_i-1, x_i, x_i+1. This relationship yields a 3x3 matrix C which has then to be inverted to find a,b,c! C The coeffiecints for the bins are C stored in the array FF as FF(i) = a, FF(i+1) = b, FF(i+2) = c C corresponding to the bin x_i, x_i+1. IMPLICIT DOUBLE PRECISION (A-H, O-Z) C CHARACTER MSG*80 C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) DIMENSION F(NX), FF(3*MXX) C C Initialize counter C K = 1 C C Loop to calculate a,b,c except first bin! C DO 10 I = NX-1, 2, -1 C DET =1D0/((XA(I+1,1) - XA(I-1,1))*(XA(I+1,1) - XA(I,1))) DET1 =1D0/((XA(I,1) - XA(I-1,1))*(XA(I+1,1) - XA(I,1))) DET2 =1D0/((XA(I,1) - XA(I-1,1))*(XA(I+1,1) - XA(I-1,1))) C FF(K) =F(I+1)*DET -F(I)*DET1 +F(I-1)*DET2 C FF(K+1) = - (XA(I-1,1) + XA(I,1))*F(I+1)*DET > + (XA(I+1,1) + XA(I-1,1))*F(I)*DET1 > - (XA(I,1) + XA(I+1,1))*F(I-1)*DET2 C FF(K+2) = + XA(I,1)*XA(I-1,1)*F(I+1)*DET > - XA(I-1,1)*XA(I+1,1)*F(I)*DET1 > + XA(I+1,1)*XA(I,1)*F(I-1)*DET2 C Print *, I, FF(K), F(I-1), F(I), F(I+1) K = K + 3 10 CONTINUE DET = 1D0/((XA(3,1) - XA(1,1))*(XA(3,1) - XA(2,1))) DET1 = 1D0/((XA(2,1) - XA(1,1))*(XA(3,1) - XA(2,1))) DET2 = 1D0/((XA(3,1) - XA(1,1))*(XA(2,1) - XA(1,1))) C C Fix a,b,c for the first bin C FF(K) = F(3)*DET - F(2)*DET1 + F(1)*DET2 C FF(K+1) = - (XA(1,1) + XA(2,1))*F(3)*DET > + (XA(3,1) + XA(1,1))*F(2)*DET1 > - (XA(3,1) + XA(2,1))*F(1)*DET2 C FF(K+2) = + XA(2,1)*XA(1,1)*F(3)*DET > - XA(1,1)*XA(3,1)*F(2)*DET1 > + XA(2,1)*XA(3,1)*F(1)*DET2 RETURN C ***************************************************** END SUBROUTINE INTEGR (NX, M, F, G, IR) C C Computes (x/y) ** M * F(y)dy/y integrated over the range [x, 1]; C Result is G. C Integand function F must be defined on an array of size NX which C covers the range [0, 1] of the variables x and y. C The output function G is returned on the same array. C The use of integration variable z defined by C C y = f(x) C C where f(x) can be any monotonic function. C C IR is an error return code: IR = 1 -- NX out of range C = 2 -- M out of range C IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER MSG*80 C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) C DIMENSION F(NX), G(NX) C DATA IWRN1, IWRN2 / 0, 0 / C IRR = 0 C IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN MSG = 'NX out of range in INTEGR call' CALL WARNI (IWRN1, NWRT, MSG, 'NX', NX, 0, MXX, 0) IRR = 1 EndIf C IF (M .LT. M1 .OR. M .GT. M2) THEN MSG ='Exponent M out of range in INTEGR' CALL WARNI (IWRN2, NWRT, MSG, 'M', M, M1, M2, 1) IRR = 2 EndIf C NX G(NX) = 0D0 C NX - 1 TEM = H(1, NX-1, -M) * F(NX-2) + H(2, NX-1, -M) * F(NX-1) > + H(3, NX-1, -M) * F(NX) IF (M .EQ. 0) THEN G(NX-1) = TEM Else G(NX-1) = TEM * XA(NX-1, M) EndIf C NX-2 : 2 DO 10 I = NX-2, 2, -1 TEM = TEM + H(1,I,-M)*F(I-1) + H(2,I,-M)*F(I) > + H(3,I,-M)*F(I+1) + H(4,I,-M)*F(I+2) IF (M .EQ. 0) THEN G(I) = TEM Else G(I) = TEM * XA(I, M) EndIf 10 CONTINUE C 1 TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3) IF (M .EQ. 0) THEN G(1) = TEM Else G(1) = TEM * XA(1, M) EndIf IR = IRR C RETURN C **************************** END SUBROUTINE NINTEGR (NX, M, FF, F, G, IR) C C Computes the integral of the more complicated pieces appearing C in the kernels of the ND case as found in .... over the range C [x, 1];The integrals are computed analytically within a bin and then C summed to give the final answer. C The input function F(y) has been replaced by FF an array containing C the constants a,b,c from the quadratic approximation of F(y). C The Result is G. C Integand function F must be defined on an array of size NX which C covers the range [0, 1] of the variables x and y. C The output function G is returned on the same array. C C IR is an error return code: IR = 1 -- NX out of range C IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER MSG*80 C PARAMETER (MXX = 525, MXQ = 25, MXF = 6, DEL = 0.0 ) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) C DIMENSION FF(3*MXX), G(NX), F(NX), HH(MXX), GH(MXX), FH(MXX) C DATA IWRN1, IWRN2 / 0, 0 / C IRR = 0 C IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN MSG = 'NX out of range in INTEGR call' CALL WARNI (IWRN1, NWRT, MSG, 'NX', NX, 0, MXX, 0) IRR = 1 EndIf C C Difference between x_min and delta is used to determine which analytic C expression for the integrals to use! T < 10^-10 is treated as C x_min = delta i.e T = 0 to avoid accuracy errors in the computation. C T > 0.98 x_min will be treated by using the analytic expression for C the integrals expanded to O(delta). This case also allows delta = 0 C i.e the diagonal case. C T = XA(1,1) - DEL T1 = 0.98*XA(1,1) IF (DEL .EQ. 0) then DIM = 0.0 ELSE DIM = 1D0/DEL ENDIF C NX G(NX) = 0 C Nx-1 -> 1 TEM = 0 K = 1 KK = 1 C C a = FF(K) b = FF(K+1) c = FF(K+2) from quadratic approximation C of F(X) knowing the function at xi , x(i-1), x(i+1) C The following are the analytical results of the integrals from C the kernels. C IF (M .EQ. 1) THEN IF (T .LT. 1D-10) THEN DO 5 IJ = NX-1, 1, -1 G(IJ) = 0 5 CONTINUE ENDIF IF (T .GT. T1) THEN DO 90 J = NX-1, 1, -1 TEM = TEM + (FF(K+1)+FF(K)*DEL)*(XA(J,-1) - XA(J+1,-1)) > + FF(K)*(XA(J+1,0) - XA(J,0)) > + (FF(K+2) + FF(K+1)*DEL)*(XA(J,-2) - XA(J+1,-2))/2. > + FF(K+2)*DEL*(XA(J,-3) - XA(J+1,-3))/3. G(J) = (XA(J,1)-DEL)*XA(J,1)*TEM K = K + 3 90 CONTINUE ELSE DO 10 I = NX-1, 1, -1 A2 = XA(I+1,1)*DIM - 1D0 A1 = XA(I,1)*DIM - 1D0 B1 = DEL*XA(I,-1) B2 = DEL*XA(I+1,-1) TEM = TEM + FF(K+2)*(B2 - B1) > + (FF(K+2) + FF(K+1)*DEL)*(XA(I,0) - XA(I+1,0)) > + (FF(K+2) + FF(K+1)*DEL + FF(K)*DEL**2) > *(LOG(A2) - LOG(A1)) G(I) = (A1+1D0)*(A1)*TEM K = K + 3 10 CONTINUE ENDIF ENDIF C IF (M .EQ. 2) THEN C Treatment of the case when x_min equals delta. IF (T .LT. 1D-10) THEN CALL INTEGR(NX, 1, F, FH, IR9) DO 20 II = NX-1, 1, -1 G(II) = FH(II) 20 CONTINUE ENDIF IF (T .GT. T1) THEN DO 91 JI = NX-1, 1, -1 DO 92 IJ = NX-1, JI, -1 TEM = TEM + (FF(K+1)*XA(JI,1) - 2.*FF(K)*XA(JI,2) + 2.*FF(K) > *XA(JI,1)*DEL)*(XA(IJ+1,0) - XA(IJ,0)) + FF(K)*XA(JI,1) > *(XA(IJ+1,1) - XA(IJ,1)) > + (FF(K)*XA(JI,3) - 2.*FF(K+1)*XA(JI,2) + FF(K+2)*XA(JI,1) > + 2.*FF(K+1)*XA(JI,1)*DEL - 4.*FF(K)*XA(JI,2)*DEL) > *(XA(IJ,-1) - XA(IJ+1,-1)) + (FF(K+1)*XA(JI,3)/2. - FF(K+2) > *XA(JI,2) + FF(K+2)*XA(JI,1)*DEL - 2.*FF(K+1)*XA(JI,2)*DEL > + FF(K)*XA(JI,3)*DEL)*(XA(IJ,-2) - XA(IJ+1,-2)) > + (FF(K+2)*XA(JI,3)/3. - 4.*FF(K+2)*XA(JI,2)*DEL/3. > + 2.*FF(K+1)*XA(JI,3)*DEL/3.)*(XA(IJ,-3) - XA(IJ+1,-3)) > + FF(K+2)*XA(JI,3)*DEL*(XA(IJ,-4) - XA(IJ+1,-4))/2. K = K + 3 92 CONTINUE G(JI) = TEM TEM = 0 K = 1 KK = 1 91 CONTINUE ELSE DO 12 J = NX-1, 1, -1 AA1 = XA(J,1)*DIM DO 23 I = NX-1, J, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM -1D0 TEM = TEM + AA1*(FF(K+2) + FF(K+1) > *DEL + FF(K)*DEL**2)*(1D0 - AA1)**2 > *(A1**-1 - A2**-1) > + FF(K)*XA(J,1)*(XA(I+1,1) - XA(I,1)) > + FF(K+2)*XA(J,1)*(XA(I,-1) - XA(I+1,-1))*AA1**2 > + AA1**2*(FF(K+1)*XA(J,1) + 2D0*FF(K+2)*AA1 > - 2D0*FF(K+2))*(XA(I+1,0) - XA(I,0)) > + AA1*(1D0 - AA1)*(FF(K+1)*DEL + 2D0*FF(K)*DEL**2 > + 2D0*FF(K+2)*AA1 + FF(K+1)*XA(J,1)) > *(LOG(A2) - LOG(A1)) K = K + 3 23 CONTINUE G(J) = TEM TEM = 0 K = 1 12 CONTINUE ENDIF ENDIF C IF (M .EQ. 3) THEN C IF (T .LT. 1D-10) THEN DO 6 IJ = NX-1, 1, -1 G(IJ) = FF(K+2)+FF(K+1)*XA(IJ,1)+FF(K)*XA(IJ,2) K = K + 3 6 CONTINUE ENDIF IF (T .GT. T1) THEN DO 93 J = NX-1, 1, -1 TEM = TEM + FF(K)*(XA(J,-1) - XA(J+1,-1)) + (FF(K+1)/2. > + FF(K)*DEL)*(XA(J,-2) - XA(J+1,-2)) + (FF(K+2)/3. > + 2.*FF(K+1)*DEL/3.)*(XA(J,-3) - XA(J+1,-3)) > + FF(K+2)*DEL*(XA(J,-4) - XA(J+1,-4))/2. G(J) = (XA(J,1)-DEL)*XA(J,2)*TEM K = K + 3 93 CONTINUE ELSE DO 13 I = NX-1, 1, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM -1D0 TEM = TEM + FF(K+2)*(DEL*XA(I,-1) - DEL*XA(I+1,-1)) > + (FF(K+1)*DEL + FF(K)*DEL**2 + FF(K+2)) > *((A1)**-1 - (A2)**-1) > + (2D0*FF(K+2) + FF(K+1)*DEL)*(XA(I+1,0) - XA(I,0) > + LOG(A1) - LOG(A2)) G(I) = (A1+1D0)**2*A1*TEM K = K + 3 13 CONTINUE ENDIF ENDIF C IF (M .EQ. 4) THEN IF (T .LT. 1D-10) THEN CALL INTEGR(NX, 0, F, GH, IR66) DO 66 II = NX-1, 1, -1 G(II) = GH(II) 66 CONTINUE ENDIF IF (T .GT. T1) THEN DO 95 JI = NX-1, 1, -1 DO 96 IJ = NX-1, JI, -1 TEM = TEM + FF(K)*(XA(IJ+1,2) - XA(IJ,2))/2. + (FF(K+1) > - 2.*FF(K)*XA(JI,1) + 2.*FF(K)*DEL) > *(XA(IJ+1,1) - XA(IJ,1)) + (FF(K+1)*XA(JI,2) - 2.*FF(K+2) > *(XA(JI,1) + DEL) - 4.*DEL*FF(K+1)*XA(JI,1) + 2.*FF(K) > *XA(JI,2)*DEL)*(XA(IJ,-1) - XA(IJ+1,-1)) + (FF(K+2) > *XA(JI,2)/2. - 2.*FF(K+2)*XA(JI,1)*DEL > + FF(K+1)*XA(JI,2)*DEL)*(XA(IJ,-2) - XA(IJ+1,-2)) > + 2.*FF(K+2)*XA(JI,2)*DEL* (XA(IJ,-3) - XA(IJ+1,-3))/3. > + (FF(K+2) - 2.*FF(K+1)*XA(JI,1) + FF(K)*XA(JI,2) > + 2.*FF(K+1)*DEL > - 4.*FF(K)*XA(JI,1)*DEL)*(XA(IJ+1,0) - XA(IJ,0)) K = K + 3 96 CONTINUE G(JI) = TEM TEM = 0 K = 1 95 CONTINUE ELSE DO 15 J = NX-1, 1, -1 AA1 = XA(J,1)*DIM DO 26 I = NX-1, J, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM -1D0 TEM = TEM + (FF(K+1) + 2D0*(DEL - XA(J,1))*FF(K)) > *(XA(I+1,1) - XA(I,1)) > + (FF(K+2) + FF(K+1)*DEL + FF(K)*DEL**2) > *(1D0 - AA1)**2*(A1**-1 > - A2**-1) + FF(K+2)*AA1**2 > *(XA(I+1,0) - XA(I,0)) > + (AA1 - 1D0)*(FF(K)*XA(J,1)*DEL - 2D0*FF(K+1)*DEL > - FF(K+2)*AA1 - 3D0*FF(K)*DEL**2 > - FF(K+2))*(LOG(A2) - LOG(A1)) + FF(K)*(XA(I+1,2) > - XA(I,2))/2D0 K = K + 3 26 CONTINUE G(J) = TEM TEM = 0 K = 1 15 CONTINUE ENDIF ENDIF C IF (M .EQ. 5) THEN IF (T .LT. 1D-10) THEN DO 7 IJ = NX-1, 1, -1 G(IJ) = FF(K+2)+FF(K+1)*XA(IJ,1)+FF(K)*XA(IJ,2) K = K + 3 7 CONTINUE ENDIF IF (T .GT. T1) THEN DO 97 II = NX-1, 1, -1 TEM = TEM + FF(K)*(XA(II+1,1) - XA(II,1)) + (FF(K+1) > + 2.*FF(K)*DEL)*(XA(II+1,0) - XA(II,0)) + (FF(K+2) > + 2.*FF(K+1)*DEL)*(XA(II,-1) - XA(II+1,-1)) + FF(K+2)*DEL > *(XA(II,-2) - XA(II+1,-2)) G(II) = (XA(II,1)-DEL)*TEM K = K + 3 97 CONTINUE ELSE DO 16 I = NX-1, 1, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM -1D0 TEM = TEM + FF(K)*(XA(I+1,1) - XA(I,1))*DEL > + (FF(K+2) + FF(K+1)*DEL + FF(K)*DEL**2) > *(A1**-1 - A2**-1) + (FF(K+1)*DEL > + FF(K)*2D0*DEL**2)*(LOG(A2) - LOG(A1)) G(I) = A1*TEM K = K + 3 16 CONTINUE ENDIF ENDIF C IF (M .EQ. 6) THEN IF (T .LT. 1D-10) THEN CALL INTEGR(NX, 0, F, GH, IR8) CALL INTEGR(NX, 1, F, FH, IR7) DO 88 II = NX-1, 1, -1 G(II) = GH(II) - FH(II) 88 CONTINUE ENDIF IF (T .GT. T1) THEN DO 98 JI = NX-1, 1, -1 DO 99 IJ = NX-1, JI, -1 TEM = TEM + FF(K)*(XA(IJ+1,2) - XA(IJ,2))/2. + (FF(K+1) > - 2.*FF(K)*XA(JI,1) + FF(K)*DEL)*(XA(IJ+1,1) - XA(IJ,1)) > + (FF(K+1)*XA(JI,2) - 2.*FF(K+2)*(XA(JI,1) - DEL/2.) > - 2.*FF(K+1)*XA(JI,1)*DEL + FF(K)*XA(JI,2)*DEL) > *(XA(IJ,-1) - XA(IJ+1,-1)) + (FF(K+2)*XA(JI,2)/2. > - FF(K+2)*XA(JI,1)*DEL + FF(K+1)*XA(JI,2)*DEL/2.) > *(XA(IJ,-2) - XA(IJ+1,-2)) + FF(K+2)*XA(JI,2)*DEL > * (XA(IJ,-3) - XA(IJ+1,-3))/3. + (FF(K+2) - 2.*FF(K+1) > *XA(JI,1) + FF(K)*XA(JI,2) + FF(K+1)*DEL - 2.*FF(K)*DEL > *XA(JI,1))*(XA(IJ+1,0) - XA(IJ,0)) K = K + 3 99 CONTINUE G(JI) = TEM TEM = 0 K = 1 98 CONTINUE ELSE DO 18 J = NX-1, 1, -1 AA1 = XA(J,1)*DIM DO 29 I = NX-1, J, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM -1D0 B1 = DEL*XA(I,-1) B2 = DEL*XA(I+1,-1) TEM = TEM + FF(K)*(XA(I+1,2) - XA(I,2))/2D0 > + (FF(K+1) + FF(K)*DEL*(1D0 - 2D0*AA1)) > *(XA(I+1,1) - XA(I,1)) + FF(K+2)*AA1**2 > *(B2 - B1) + (2D0*FF(K+2) > *AA1 - FF(K+1)*AA1*XA(J,1) > - FF(K+2)*AA1**2) > *(XA(I+1,0) - XA(I,0)) > + (AA1-1D0)**2*(FF(K+2) + FF(K+1)*DEL > + FF(K)*DEL**2)*(LOG(A2) - LOG(A1)) K = K + 3 29 CONTINUE G(J) = TEM TEM = 0 K = 1 18 CONTINUE ENDIF ENDIF C IF (M .EQ. 7) THEN IF (T .LT. 1D-10) THEN DO 8 II = NX-1, 1, -1 G(II) = 0 8 CONTINUE ENDIF IF (T .GT. T1) THEN DO 48 IJ = NX-1, 1, -1 TEM = TEM + FF(K)*(XA(IJ+1,1) - XA(IJ,1)) + (FF(K+1) + FF(K) > *DEL)*(XA(IJ+1,0) - XA(IJ,0)) + (FF(K+2) + FF(K+1)*DEL) > *(XA(IJ,-1) - XA(IJ+1,-1)) G(IJ) = (XA(IJ,1)-DEL)*TEM K = K + 3 48 CONTINUE ELSE DO 222 I = NX-1, 1, -1 A1 = XA(I,1)*DIM - 1D0 A2 = XA(I+1,1)*DIM - 1D0 TEM = TEM + DEL*FF(K)*(XA(I+1,1) - XA(I,1)) > + FF(K+2)*(XA(I,0) - XA(I+1,0)) > + (FF(K+2) + FF(K+1)*DEL + FF(K)*DEL**2) > *(LOG(A2) - LOG(A1)) G(I) = A1*TEM K = K + 3 222 CONTINUE ENDIF ENDIF C IF (M .EQ. 8) THEN IF (T .LT. 1D-10) THEN CALL HINTEG(NX, F, HH) DO 113 J = NX-1, 1, -1 G(J) = -HH(J) 113 CONTINUE ENDIF IF (T .GT. T1) THEN DO 83 LL = NX-1, 1, -1 DO 84 II = NX-1, LL, -1 TEM = TEM + FF(K)*(XA(II+1,2) - XA(II,2))/2D0 + (FF(K+1) > + FF(K)*DEL)*(XA(II+1,1) - XA(II,1)) + (FF(K+2) + FF(K+1) > *DEL)*(XA(II+1,0) - XA(II,0)) + DEL*XA(LL,1) > *(XA(II+1,-1) - XA(II,-1))*(FF(K+1) + FF(K)*XA(LL,1)) K = K + 3 84 CONTINUE G(LL) = - TEM - F(LL)*LOG(1D0 - DEL) TEM = 0 K = 1 83 CONTINUE ELSE DO 202 I = NX-1,1, -1 DO 112 J = NX-1,I, -1 A1 = XA(J,1)*DIM - 1D0 A2 = XA(J+1,1)*DIM - 1D0 TEM = TEM + (FF(K+1)+FF(K)*DEL)*(XA(J+1,1)-XA(J,1)) > + FF(K)*(XA(J+1,2) - XA(J,2))/2D0 > + (FF(K+1)*DEL - FF(KK+1)*XA(I,1) + FF(K)*DEL**2 > - FF(KK)*XA(I,2) + FF(K+2) - FF(KK+2)) > *(LOG(A2)-LOG(A1)) + (FF(KK)*XA(I,2) + FF(KK+1)*XA(I,1) > + FF(KK+2))*(XA(J+1,0) - XA(J,0)) K = K + 3 112 CONTINUE G(I) = - TEM - F(I)*LOG(1D0-DEL) TEM = 0 K = 1 KK = KK + 3 202 CONTINUE ENDIF ENDIF RETURN C ********************************************* END SUBROUTINE HINTEG (NX, F, H) C C Computes the integral [yF(y)-xF(x)]/(y-x) * dy/y over the range [x, 1]; C then add F(x) * Ln (1-x) to get Int 1/(1-x/y)(sub+)F(y)dy/y. C C The input function F must be specified on an array of size NX over C the range (0, 1] of the x variable. In order to allow a possible c singularity at x = 0, the first mesh-point is at x = 1/nx, not 0. C The output function H is given on the same array as above. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / HINTEC / GH(NDG, MXX) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) C DIMENSION F(NX), H(NX), G(MXX) C DZ = 1D0 / (NX-1) C Loop to calculate the required integral C |--|--|-----------|--|--| C Iy: I I+1 ... nx C Kz: 1 2 ... np C C ------------------------------ C Do each integral in turn. DO 20 I = 1, NX-2 C Number of points in the I-th integral NP = NX - I + 1 C Evaluate the first two bins of the integrals TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2) C Evaluate the integrand for the I-th integral DO 30 KZ = 3, NP IY = I + KZ - 1 C DXDZ is the Jacobian due to the change of variable X->Z C W = XA(I,1) / XA(IY,1) G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W) C 30 CONTINUE C HTEM = SMPSNA (NP-2, DZ, G(3), ERR) TEM1 = F(I) * ELY(I) H(I) = TEM + HTEM + TEM1 C 20 CONTINUE C H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0)) H(NX) = 0 C RETURN C **************************** END Function zfrmx(x) C Transformation x --> z for convenient and efficient solution of AP Eq. C (and for displaying the behavior of f(x,Q) vs. x!). C Motivation: (i) evenly spaced grids in z sample physics more evenly in C the range Xmin < x < 1; (ii) z spans the std range [0,1] C x : Xmin --------------- 1 C z : 0 --------------- 1 C Original version was essentially: z = A Ln x + B x + C: C great, but clumsy to invert; C This version use the simpler fn: z = A x**Dd + B C with Dd ~ 0.3 C (fractional power simulates the logarithm; and it inverts easily.) C With the boundary conditions shown, there is one adjustable parameter C in each case: Xcr in the original version; Dd in this version. IMPLICIT DOUBLE PRECISION (A-H, O-Z) Logical Lstx PARAMETER (MXX = 525) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX C Make sure parameters and coefficients are saved Save Xm, Dd, Dd0, A, B C For clarity, use inline functions C Optimizers will expand these ex (x) = x ** Dd dedx (x) = Dd * x **(Dd-1) xe (e) = e ** (1/Dd) dxde (e) = e **(1/Dd - 1D0) / Dd Data Xm, Dd, Dd0 / -99D0, 99D0, 0.30D0 / C Hard-wire the exponent Dd to 0.3 for the moment C May want to allow it to be fine-tuned by the parameter Xcr later If (Xm.Ne.Xmin .or. Dd.Ne.Xcr) Then Xcr = Dd0 C Remove above line when Xcr is activated Xm = Xmin Dd = Xcr exm = ex (Xm) ex1 = ex (1D0) A = 1D0 / (ex1 - exm) B = - A * exm EndIf Zx = A * ex(x) + B Zfrmx = Zx Return C --------------- Entry DzDx (x) tem = dedx (x) DzDx = tem * A Return C --------------- Entry xfrmz (z) C Transformation z --> x : Inverse of zfrmx (x) If (Xm.Ne.Xmin .or. Dd.Ne.Xcr) Then Xcr = Dd0 C Remove above line when Xcr is activated Xm = Xmin Dd = Xcr exm = ex (Xm) ex1 = ex (1D0) A = 1D0 / (ex1 - exm) B = - A * exm EndIf exz = (z - B) / A Xz = xe (exz) xfrmz = Xz Return C ----------------- Entry DxDz (z) exz = (z - B) / A tem = dxde (exz) dxdz = tem / A Return C ****************** End C SUBROUTINE PARPDF (IACT, NAME, VALUE, IRET) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: routines to input and output evolution parameters C 5parpdf: :parpdf evlpar evlget evlgt1 evlin evlset C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5parpdf.f,v 5.9 96/10/20 13:42:01 wkt Exp $ C $Log: 5parpdf.f,v $ c Revision 5.9 96/10/20 13:42:01 wkt c 1. Consolidated with Liang's Cteq4 version; c 2. Cleaned up numerous nuisances (variables not used, ..) which cause c compilation warnings c 3. ZfrmX and XfrmZ and ancillaries replaced with new version based on c fractional power rather than log + linear transformation; c 4. New Cteq, Mrs, Grv switches put in; but need their programs and c tables to run. c c Revision 5.6 96/06/02 23:35:47 wkt c minor c c Revision 5.5 96/06/02 23:16:05 wkt c Kf variable replaced by NfMx. (Kf = 2*NfMx + 1 glu + Madd for Sqrk ...etc.) c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/19 14:12:40 wkt c Parpdf updated to correct the UpCase problem c C ==================================================================== C SUBROUTINE PARPDF (IACT, NAME, VALUE, IRET) C C Actions: 0 type list of variables on unit VALUE. C 1 set variable with name NAME to VALUE, if C it exists, Else set IRET to 0. C 2 find value of variable. If it does not exist, C set IRET to 0. C 3 request values of all parameters from terminal. C 4 type list of all values on unit VALUE. c C IRET = 0 variable not found. C 1 successful search. C 2 variable found, but bad value. C 3 bad value for IACT. C 4 no variable search (i.e., IACT is 0, 3, or 4). C C If necessary, VALUE is converted to integer by NINT C C Use ILEVEL and START1 to start search for variable names close C to previous name to ensure effiency when reading in many values. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C CHARACTER NAME*(*), Uname*10 C LOGICAL START1 C COMMON / IOUNIT / NIN, NOUT, NWRT C DATA ILEVEL, LRET / 1, 1 / C JRET = IRET CALL UPC (NAME, Ln, Uname) IF (IACT .EQ. 0 .OR. IACT .EQ. 4) > IVALUE = NINT (VALUE) START1 = (IACT .NE. 1) .AND. (IACT .NE. 2) IF (START1) ILEVEL = 1 C GOTO (1, 2), ILEVEL C 1 START1 = .TRUE. ILEVEL = 0 CALL PARQCD (IACT, Uname(1:Ln), VALUE, JRET) IF (JRET .EQ. 1) GOTO 11 IF (JRET .EQ. 2) GOTO 12 IF (JRET .EQ. 3) GOTO 13 IF (JRET .GT. 4) GOTO 15 ILEVEL = ILEVEL + 1 2 CALL EVLPAR (IACT, Uname(1:Ln), VALUE, JRET) IF (JRET .EQ. 1) GOTO 11 IF (JRET .EQ. 2) GOTO 12 IF (JRET .EQ. 3) GOTO 13 IF (JRET .GT. 4) GOTO 15 ILEVEL = ILEVEL + 1 C IF (.NOT. START1) GOTO 1 C IF (JRET .EQ. 0) GOTO 10 C C Arrive here if IACT = 0, 3 and all is OK (IRET=4): 9 CONTINUE C WRITE (IVALUE, 100) GOTO 14 C Exits: 10 CONTINUE 11 CONTINUE 12 CONTINUE 13 CONTINUE 14 CONTINUE 15 CONTINUE C LRET is used for debugging purpose only IF (JRET .NE. 4) LRET = JRET IF (LRET.EQ.0 .OR. LRET.EQ.2 .OR. LRET.EQ.3) THEN PRINT *, 'Error in PARPDF: IRET, IACT, NAME, VALUE =', > LRET, IACT, NAME, VALUE PAUSE EndIf IRET= JRET RETURN C 100 FORMAT (/) C **************************** END SUBROUTINE EVLPAR (IACT, NAME, VALUE, IRET) C C For STANDARD codes on Iact and Iret, see SUBROUTINE PARPDF C C Additional options added for Version 7.2 and up of PDF package: C C IACT = 5 find value of variable from the alternate set. C 6 type list of all values from the alternate set. C c These options must be called from EVLPAR directly, not from the C front-end unit PARPDF because PARQCD cannot handle Iact > 4 C C NAME is assumed upper-case. C If necessary, VALUE is converted to integer by NINT C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C CHARACTER*(*) NAME C COMMON / IOUNIT / NIN, NOUT, NWRT C IRET = 1 IF (IACT .EQ. 0) THEN WRITE ( NINT (VALUE) , 101) 101 FORMAT (/ ' Initiation parameters: Qini, Ipd0, Ihdn ' / > ' Maximum Q, Order of Alpha: Qmax, IKNL ' / > ' X- mesh parameters : Xmin, Xcr, Nx ' / > ' LnQ-mesh parameters : Nt, Jt ' / > ' # of parton flavors : NfMx ' /) IRET = 4 ElseIF (IACT .EQ. 1) THEN CALL EVLSET (NAME, VALUE, IRET) ElseIF (IACT .EQ. 2) THEN CALL EVLGET (NAME, VALUE, IRET) ElseIF (IACT .EQ. 5) THEN CALL EVLGT1 (NAME, VALUE, IRET) ElseIF (IACT .EQ. 3) THEN CALL EVLIN IRET = 4 ElseIF (IACT .EQ. 4) THEN CALL EVLOUT ( NINT (VALUE) ) IRET = 4 ElseIF (IACT .EQ. 6) THEN CALL EVLOT1 ( NINT (VALUE) ) IRET = 4 Else IRET = 3 EndIf C RETURN C **************************** END SUBROUTINE EVLGET (NAME, VALUE, IRET) C C Gets VALUE of variable named NAME. C C IRET = 0 variable not found. C 1 success. C C NAME is assumed upper-case, and VALUE real. C If necessary, VALUE is converted to integer by NINT C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C CHARACTER*(*) NAME C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx C IRET = 1 C IF (NAME .EQ. 'QINI') THEN VALUE = QINI ElseIF (NAME .EQ. 'IPD0') THEN VALUE = IPD0 ElseIF (NAME .EQ. 'IHDN') THEN VALUE = IHDN ElseIF (NAME .EQ. 'QMAX') THEN VALUE = QMAX ElseIF (NAME .EQ. 'IKNL') THEN VALUE = IKNL ElseIF (NAME .EQ. 'XCR') THEN VALUE = XCR ElseIF (NAME .EQ. 'XMIN') THEN VALUE = XMIN ElseIF (NAME .EQ. 'NX') THEN VALUE = NX ElseIF (NAME .EQ. 'NT') THEN VALUE = NT ElseIF (NAME .EQ. 'JT') THEN VALUE = JT ElseIF (NAME .EQ. 'NFMX') THEN VALUE = NfMx Else IRET = 0 EndIf C RETURN C ____________________________ C ENTRY EVLOUT (NOUUT) C C Write current values of parameters to unit NOUUT C WRITE (NOUUT, 131) QINI, IPD0, IHDN, QMAX, IKNL, > XMIN, XCR, NX, NT, JT, NfMx C 131 FORMAT ( / >' Current parameters and values are: '// >' Initiation parameters: Qini, Ipd0, Ihdn = ', F8.2, 2I8 // >' Maximum Q, Order of Alpha: Qmax, IKNL = ', 1PE10.2, I6// >' X- mesh parameters : Xmin, Xcr, Nx = ', 2(1PE10.2), I8// >' LnQ-mesh parameters : Nt, Jt = ', 2I8 // >' # of parton flavors : NfMx = ', I8 /) C RETURN C **************************** END SUBROUTINE EVLGT1 (NAME, VALUE, IRET) C C COPY OF EVLGET for the alternate set C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C LOGICAL LSTX CHARACTER*(*) NAME C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / X1RRAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / Q1RAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / E1LPAR / AL, IKNL, IPD0, IHDN, NfMx C IRET = 1 C IF (NAME .EQ. 'QINI') THEN VALUE = QINI ElseIF (NAME .EQ. 'IPD0') THEN VALUE = IPD0 ElseIF (NAME .EQ. 'IHDN') THEN VALUE = IHDN ElseIF (NAME .EQ. 'QMAX') THEN VALUE = QMAX ElseIF (NAME .EQ. 'IKNL') THEN VALUE = IKNL ElseIF (NAME .EQ. 'XCR') THEN VALUE = XCR ElseIF (NAME .EQ. 'XMIN') THEN VALUE = XMIN ElseIF (NAME .EQ. 'NX') THEN VALUE = NX ElseIF (NAME .EQ. 'NT') THEN VALUE = NT ElseIF (NAME .EQ. 'JT') THEN VALUE = JT ElseIF (NAME .EQ. 'NFMX') THEN VALUE = NfMx Else IRET = 0 EndIf C RETURN C ____________________________ C ENTRY EVLOT1 (NOUUT) C C Write current values of parameters to unit NOUUT C WRITE (NOUUT, 131) QINI, IPD0, IHDN, QMAX, IKNL, > XMIN, XCR, NX, NT, JT, NfMx C 131 FORMAT ( / >' Current parameters and values are: '// >' Initiation parameters: Qini, Ipd0, Ihdn = ', F8.2, 2I8 // >' Maximum Q, Kernel ID : Qmax, IKNL = ', 1PE10.2, I6// >' X- mesh parameters : Xmin, Xcr, Nx = ', 2(1PE10.2), I8// >' LnQ-mesh parameters : Nt, Jt = ', 2I8 // >' # of parton flavors : NfMx = ', I8 /) C RETURN C **************************** END C SUBROUTINE EVLIN C C Solicits parameters in EVL calculations IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / IOUNIT / NIN, NOUT, NWRT CALL EVLGET ('QINI', V1, IR1) CALL EVLGET ('IPD0', V2, IR2) CALL EVLGET ('IHDN', V3, IR3) 1 WRITE (NOUT, 101) V1, V2, V3 101 FORMAT (/' Current values of parameters QINI, IPD0, IHDN: ', > 1PE10.2, 2I6 / '$Type new values: ' ) READ(NIN, *, ERR=302) V1, V2, V3 C CALL EVLSET ('QINI', V1, IRET1) CALL EVLSET ('IPD0', V2, IRET2) CALL EVLSET ('IHDN', V3, IRET3) goto 301 302 write(nout,120) goto 1 C 301 IF ((IRET1.NE.1) .OR. (IRET2.NE.1) .OR. (IRET3.NE.3)) THEN 20 WRITE (NOUT, 120) GOTO 1 EndIf CALL EVLGET ('QMAX', V1, IR1) CALL EVLGET ('NT', V2, IR2) CALL EVLGET ('JT', V3, IR3) C 2 WRITE (NOUT, 102) V1, V2, V3 102 FORMAT (/' Current values of parameters QMAX, NT, JT: ', > 1PE10.2, I6, I6 / '$Type new values: ' ) READ(NIN, *, ERR=304) V1, V2, V3 C CALL EVLSET ('QMAX', V1, IRET1) CALL EVLSET ('NT', V2, IRET2) CALL EVLSET ('JT', V3, IRET3) goto 303 304 write(nout,120) goto 2 C 303 IF ((IRET1.NE.1) .OR. (IRET2.NE.1) .OR. (IRET3.NE.1)) THEN 22 WRITE (NOUT, 120) GOTO 2 EndIf C CALL EVLGET ('XMIN', V1, IR1) CALL EVLGET ('XCR', V2, IR2) CALL EVLGET ('NX', V3, IR3) CALL EVLGET ('NFMX', V4, IR4) CALL EVLGET ('IKNL', V5, IR5) 3 WRITE (NOUT, 103) V1, V2, V3, V4, V5 103 FORMAT(/' Current values of parameters XMIN, XCR, NX,NFMX,IKNL: ', > 2(1PE12.3), 3I6 / '$Type new values: ' ) READ(NIN, *, ERR=22) V1, V2, V3, V4, V5 C CALL EVLSET ('XMIN', V1, IRET1) CALL EVLSET ('XCR', V2, IRET2) CALL EVLSET ('NX', V3, IRET3) CALL EVLSET ('NFMX', V4, IRET4) CALL EVLSET ('IKNL', V5, IRET5) C IF ( (IRET1 .NE. 1) .OR. (IRET2 .NE. 1) .OR. (IRET3 .NE. 1) > .OR. (IRET4 .NE. 1) .OR. (IRET5 .NE. 1) ) THEN 23 WRITE (NOUT, 120) GOTO 3 EndIf C 120 FORMAT(' Bad values, Try again!' /) C RETURN C **************************** END SUBROUTINE EVLSET (NAME, VALUE, IRET) C C Sets variable named NAME to VALUE. C C IRET = 0 variable not found. C 1 success. C 2 variable found, but bad value. C C NAME is assumed upper-case, and VALUE real. C If necessary, VALUE is converted to integer by NINT C IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX C CHARACTER*(*) NAME C PARAMETER (MXX = 525, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / QARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG COMMON / EVLPAC / AL, IKNL, IPD0, IHDN, NfMx C IRET = 1 IF (NAME .EQ. 'QINI') THEN IF (VALUE .LE. 0) GOTO 12 QINI = VALUE ElseIF (NAME .EQ. 'IPD0') THEN ITEM = NINT(VALUE) ITM = MOD (ITEM, 20) IF (ITM .LT. 0 .OR. ITM .GT. 9) GOTO 12 IPD0 = ITEM C ElseIF (NAME .EQ. 'IHDN') THEN ITEM = NINT(VALUE) IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12 IHDN = ITEM ElseIF (NAME .EQ. 'QMAX') THEN IF (VALUE .LE. QINI) GOTO 12 QMAX = VALUE ElseIF (NAME .EQ. 'IKNL') THEN ITMP = NINT(VALUE) ITEM = ABS(ITMP) IF (ITEM .NE. 1 .AND. ITEM .NE. 2) GOTO 12 IKNL = ITMP ElseIF (NAME .EQ. 'XCR') THEN IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12 XCR = VALUE LSTX = .FALSE. ElseIF (NAME .EQ. 'XMIN') THEN IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12 XMIN = VALUE LSTX = .FALSE. ElseIF (NAME .EQ. 'NX') THEN ITEM = NINT(VALUE) IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12 NX = ITEM LSTX = .FALSE. ElseIF (NAME .EQ. 'NT') THEN ITEM = NINT(VALUE) IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12 NT = ITEM ElseIF (NAME .EQ. 'JT') THEN ITEM = NINT(VALUE) IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12 JT = ITEM ElseIF (NAME .EQ. 'NFMX') THEN ITEM = NINT(VALUE) IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12 NfMx = ITEM Else IRET = 0 EndIf C RETURN C Error exit: 12 IRET = 2 C RETURN C **************************** END C FUNCTION PMSDIS (ISET, IHDRN, IPRTN, XX, QQ, IAT) C (These general comments are enclosed in the C lead subprogram in order to survive forsplit.) C ==================================================================== C GroupList: Misc. subprograms C DIS scheme <--> MS-bar scheme conversion :pmsdis c2qx sqrk C non-essential x <--> z conversion: x2zconv C Calculate moments: amomen C EndGroupList C ==================================================================== C $Header: /users/wkt/1hep/2pdf/evl/v5/RCS/5pdfmis.f,v 5.9 96/10/20 13:48:35 wkt Exp $ C $Log: 5pdfmis.f,v $ c Revision 5.9 96/10/20 13:48:35 wkt c 5dis2ms absorbed into 5pdfmis c c Revision 5.5 96/06/02 23:15:33 wkt c Commons EVLPAC and PEVLDT modified. c c Revision 5.1 95/07/19 14:11:41 wkt c comments only c C ==================================================================== C FUNCTION PMSDIS (ISET, IHDRN, IPRTN, XX, QQ, IAT) C Front-end function to DIS2MS which transforms between DIS and MS-bar C scheme distribution functions: C Also set up the common block MSTDIS for use by the integrand functions C On input, IACT = 1 : performs the DIS --> MS-BAR conversion C -1 : " DIS <-- MS-BAR " C 0 : no conversion C other: error condition IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / MSTDS1 / X, Q, NF, JSET, JHDRN, JPRTN, NCNT C EXTERNAL C2QX, QQY, QGY, GQY SAVE C DATA TINY, D1P, DUM, IST / 1D-20, 1.001, -999., -5 / DATA AERR, RERR, SML, D1M / 0.0, 0.01, 1.D-6, 0.999 / IACT = IAT IF (IACT .EQ. 0) THEN PMSDIS = PDF(ISET, IHDRN, IPRTN, XX, QQ, IR) PRINT *, 'IACT=0 in PMSDIS call, no convertion is done.' RETURN ElseIF (ABS(IACT) .GT. 1) THEN PRINT *, 'Illegal value of Iact in PMSDIS: IACT =', IACT PMSDIS = DUM RETURN EndIf IF (ISET .NE. IST .OR. IHDRN .NE. JHDRN) THEN IST = ISET JSET = ISET JHDRN = IHDRN IF(IST.LE.0 .OR. IST.GT.31 .OR. (IST.GT.11.AND.IST.LT.21)) THEN PRINT *, 'ISET = ', ISET, ' out of range in PMSDIS.' PMSDIS = DUM IR = 11 RETURN EndIf EndIf JPRTN = IPRTN Q = QQ NF = NFL(Q) X = XX IF (X .LE. TINY .OR. X .GT. D1P) THEN PRINT '(A,1PE13.4,A)', ' X =', X, ' out of range in PMSDIS.' PMSDIS = DUM IR = 13 RETURN ElseIF (X .GT. D1M) THEN PMSDIS = 0. RETURN EndIf AL = ALPI (Q) C Begin computation TEM0 = PDF (JSET, JHDRN, JPRTN, X, Q, IR0) IF (TEM0 .LT. -SML) > PRINT '(A, 3(1PE13.4))', 'PDF < 0 in PMSDIS', TEM0, X, Q IF (TEM0 .LT. SML) THEN PMSDIS = TEM0 RETURN EndIf AERR = ABS(TEM0 / 10.) IF (JPRTN .EQ. 0) THEN TMGQ1=-ADZINT (C2QX,D0,X,D0,RERR,ER2,IR, 2,1) * SQRK(X) TMGQ2=-ADZINT (GQY, X, D1, AERR, RERR, ER3, IR, 1,2) TMGG =-ADZINT (QGY, X,D1, AERR,RERR,ER1,IR, 1,2) * 2.*NF TEM = TEM0 + AL * (TMGQ2 - TMGQ1 + TMGG) * IACT Else TMQQ1= ADZINT (C2QX, D0,X,D0,RERR,ER1,IR, 2,1) * TEM0 TMQQ2= ADZINT (QQY, X, D1, AERR, RERR, ER2, IR, 1,2) TMQG = ADZINT (QGY, X, D1, AERR, RERR, ER3, IR, 1,2) TEM = TEM0 + AL * (TMQQ2 - TMQQ1 + TMQG) * IACT EndIf PMSDIS = TEM C RETURN C **************************** END FUNCTION C2QX (XX) C C Integrands for the convolution integral in DIS to MS-bar scheme C transformation IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D1=1D0,D2=2D0,D3=3D0,D4=4D0,D5=5D0,D6=6D0,D10=1D1) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / MSTDS1 / X, Q, NF, JSET, JHDRN, JPRTN, NCNT DATA 1 D0, DUM / 0.0, 0.0 / 1 IW1, IW2, IW3 / 3 * 0 / C Statement functions QRK(Y) = PDF (JSET, JHDRN, JPRTN, Y, Q, IRT) GLU(Y) = PDF (JSET, JHDRN, 0, Y, Q, IRT) C C2Q(X) = (2./3.) * (1.+ X**2) / (1.- X) * LOG(X/(1.-X)) > + 1./ (1.- X) - 2.- (4./3.) * X C C2G(X) = (1./4) * (X**2 + (1.- X)**2) * (LOG( X/(1.-X) ) + 1.) > - (3./2.) * X * (1.-X) C XV = XX IF (XV .LE. D0) THEN PRINT *, 'X < 0 out of range in C2QX, X = ', XV TEM = DUM ElseIF (XV .LT. D1) THEN TEM = C2Q(XV) Else PRINT *, 'X > 1 out of range in C2QX, X =', XV TEM = DUM EndIf C2QX = TEM RETURN C -------------------- ENTRY QQY (YY) C Quark to Quark piece Y = YY IF (Y .GT. X .AND. Y .LT. D1) THEN TEM = C2Q (Y) * (QRK(X/Y) - QRK(X)*Y) / Y ElseIF (Y .EQ. X) THEN TEM = C2Q (Y) * (-QRK(X)) Else TEM = 0. IF (Y .NE. D1) > CALL WARNR(IW1, NWRT, 'Y out of range in QQY', 'Y',Y,X,D1,I1) EndIf C QQY = TEM C NCNT = NCNT + 1 RETURN C ---------------------------- ENTRY GQY (YY) C Singlet quark to gluon Y = YY IF (Y .GT. X .AND. Y .LT. D1) THEN TEM = C2Q (Y) * (SQRK(X/Y) - Y *SQRK(X)) / Y ElseIF (Y .EQ. X) THEN TEM = C2Q (Y) * (-SQRK(X)) ElseIF (Y .EQ. D1) THEN TEM = 0. Else CALL WARNR (IW2, NWRT, 'Y out of range in GQY', 'Y',Y,X,D1,I1) TEM = 0. EndIf C GQY = TEM C NCNT = NCNT + 1 RETURN C ---------------------------- ENTRY QGY (YY) C Y = YY IF (Y .GT. X .AND. Y .LT. D1) THEN TEM = C2G (Y) * GLU(X/Y) / Y Else TEM = 0. IF (Y .LT. X .OR. Y .GT. D1) > CALL WARNR(IW3, NWRT, 'Y out of range in QGY', 'Y',Y,X,D1,I1) EndIf C QGY = TEM C NCNT = NCNT + 1 RETURN C **************************** END FUNCTION SQRK (Y) C Returns singlet quark PDF IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / MSTDS1 / X, Q, NF, JSET, JHDRN, JPRTN, NCNT TEM = 0 DO 201 IFL = 1, NF TEM = TEM + PDF(JSET, JHDRN, IFL, Y, Q, IR) > + PDF(JSET, JHDRN,-IFL, Y, Q, IR) 201 CONTINUE SQRK = TEM RETURN C **************************** END C FUNCTION AMOMEN (AM, NX, F) C Given the function F(x) defined in the interval (0,1] on the array C (1, NX), evenly distributed in the variable z where x = xfrmz (Z) C this routine returns its Amth moment defined as Int F(x) (x**Am) dx/x, C C The value of F(x) as x --> 0 is obtained by extrapolation in the C variable z from the first three points F1, F2, and F3. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 525) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) C COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) C DIMENSION F(NX), HH(NDH, MXX), A(NDG), DI(0:3) C IF (NX .LE. 10 .OR. NX .GT. MXX) THEN CALL WARNI(IWRN,NWRT,'NX out of range in AMOMEN subroutine ', > 'NX', NX, 10, MXX, 1) END IF C Compute "A" matrix DO 30 I = NX-1, 1, -1 DO 40 K = 1, NDG AK = AM - 2D0 + K IF (AK .EQ. 0D0) THEN A(K) = XA(I+1, 0) - XA(I, 0) Else A(K) = (XA(I+1,1)**AK - XA(I,1)**AK) / AK EndIf 40 CONTINUE C Compute "H" matrix DO 41 J = 1, NDH TEM = 0 DO 43 L = 1, NDG TEM = TEM + A(L) * GB(L,J,I) 43 CONTINUE HH(J,I) = TEM 41 CONTINUE 30 CONTINUE C Calculate integral from X= X(I=4) to 1.(I=Nx) C Patterned after INTEGR TEM1 = 0D0 DO 60 I = NX-1, 4, -1 TEM1 = TEM1 + HH(1,I)*F(I-1) + HH(2,I)*F(I) > + HH(3,I)*F(I+1) + HH(4,I)*F(I+2) 60 CONTINUE C Calculate contribution from the first 3 bins and C Use the first 3 bins to extrapolate to X = 0. TEM2 = 0 DO 61 I = 1, 3 DI(I) = HH(1,I)*F(I-1) + HH(2,I)*F(I) > + HH(3,I)*F(I+1) + HH(4,I)*F(I+2) TEM2 = TEM2 + DI(I) 61 CONTINUE C Quadratic extrapolation gives DI(0) = 3D0 * (DI(1) - DI(2)) + DI(3) C TEM = TEM1 + TEM2 + DI(0) C AMOMEN = TEM RETURN C **************************** END Subroutine X2Zconv (Nx, Xx, Zz, Xmn, Ixz) C Convert x to z: Ixz = 1 : x --> z ; C 2 : z --> x C Xmn is the Xmin parameter in the conversion C formula (along with Xcr which has 1.5 as default). C Xmn = 1e-3 or 1e-4 are good choices; C Z = (0,1) for x = (Xmn, 1); if x