Subroutine Evlini (Nini, Xmin,Qmax,Qmss, Nx,Nq,Lpr) C !! Warning: C !! The input nnn.ini file must be open and associated with Unit# = Nini C Explanations: C * Read in nnn.ini file (Unit # = Nini) for input evolu parameters from fit nnn C * Do evolution in the range (Xmin, 1) and (Qini, Qmax) C (Qini is given in .ini) C * Qmss(4:6) are heavy quark masses to be used for charm, bottom, and top. C * Nx and Nq are # of grids in (x,Q) - tradeoff between accuracy/efficiency C * Set up Pdf(Iset, Ihadn, Iparton, x, Q, Irt) C Pdf's to be used as Iset = 10 C * Lpr : print level -- for diagnostics and debugging C To simplify the logistics, the total # of flavors will be set always = 6 C To generate PDF's with less than 6 flavors, simply decouple the unwanted C heavy flavors by setting their masses beyond Qmax. IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (NF0 = 4, nshp = 4,NEX = nshp+2) CHARACTER head*78 PARAMETER (MXQ = 25, MxQrk = 6) PARAMETER (MXX = 105, MXF = 6) External FitIni C Dimension Qmss(3) Common /FitIpt/ BP(0:NEX, -NF0:2) Common /InpFin/ A(Nshp), B(0:NEX, -NF0:2), Ifun DATA D0, Ahdn / 0.0, 1.0 / C ---------------------------- C Initialize pdf & Qcd parameters Dumm = SetPdf() Dumm = SetQcd() Print *, 'Reading input XXX.ini.....' C Nini is the input Xnnn.ini file Call iniread(Nini, qini, anfl, ordn, ordi, blam, Qm4, Qm5) If (Lpr .GE. 2) Print '(A/6F9.3, I5)', >' Nini, qini, anfl, ordn, ordi, blam, ifun=', > float(Nini), qini, anfl, ordn, ordi, blam, ifun C ---------------------------- C Set up evolution parameters C Will always generate 6 total flavors Tnfl = 6.0 CALL PARQCD(1, 'NFL', Tnfl, JR) CALL PARPDF(1, 'IHDN', Ahdn, JR) C From subroutine arguments Anx = Nx Ant = Nq if(Qmss(1).ne.Qm4 .or. Qmss(2).ne.Qm5) then print *,'WARNING:M4,M5 not quite match' write (*,'(A4,2f10.4)')'M4:',Qmss(1),Qm4 write (*,'(A4,2f10.4)')'M5:',Qmss(2),Qm5 endif Call parqcd(1, 'M4', Qmss(1), jr) Call parqcd(1, 'M5', Qmss(2), jr) Call parqcd(1, 'M6', Qmss(3), jr) Call parpdf(1, 'QMAX', qmax, jr) Call parpdf(1, 'XMIN', xmin, jr) Call parpdf(1, 'NX', anx, jr) Call parpdf(1, 'NT', ant, jr) C From .ini CALL PARQCD(1, 'ORDR', ORDN, JR) CALL PARPDF(1, 'IKNL', ORDI, JR) CALL PARPDF(1, 'QINI', QINI, IR) C Set Lambda for Anfl effective flavors C without affecting total number of flavors Neff = Nint (Anfl) IordEvl = Nint (Ordi) Call SetLam (Neff, Blam, IordEvl) c Fill B() matrix (doing the equivalent job of Inptn in fitpac.) DO 10 IFL = -NF0, 2 DO 10 IEX = 0, NEX B(IEX, IFL) = BP(IEX, IFL) 10 CONTINUE Print '(A/A)', 'Initial distribution parameters:' > ,' N A1 A2 A3 A4' do 11 ifl = 2,-nf0,-1 11 print '(1X, 5(1pE11.3))', (b(iex,ifl),iex=0,4) C Print *, 'Start evolution ....' C ---------------------------- Anout = 6 Call ParPdf (4, 'DUM', ANout, Jr) Call Evolve (FitIni, Irt) Return C **************************** END subroutine iniread(Nini,qini,anfl,ordn,ordi,blam,Qm4,Qm5) implicit double precision (a-h,o-z) character*80 lin80 parameter(nf0=4,nshp=4,nex=nshp+2) Common /InpFin/ A(Nshp), B(0:NEX, -NF0:2), Ifun Common /FitIpt/ BP(0:NEX, -NF0:2) read(Nini,'(a80)')lin80 read(Nini,'(a80)')lin80 read(Nini,*)ordn,ordi,anfl,blam,qini,afun,Qm4,Qm5 ifun = nint(afun) icn=0 5 read(Nini,'(a80)')lin80 if(lin80(1:5).eq.' Coef')icn=1 if(icn.eq.0) goto 5 read(Nini,'(a80)')lin80 do 10 ifl=2,-4,-1 10 read(Nini,'(6f12.3)')(bp(iex,ifl),iex=0,5) return end SUBROUTINE PARPDF (IACT, NAME, VALUE, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER NAME*(*), Uname*10 LOGICAL START1 COMMON / IOUNIT / NIN, NOUT, NWRT DATA ILEVEL, LRET / 1, 1 / 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 GOTO (1, 2), ILEVEL 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 IF (.NOT. START1) GOTO 1 IF (JRET .EQ. 0) GOTO 10 9 CONTINUE GOTO 14 10 CONTINUE 11 CONTINUE 12 CONTINUE 13 CONTINUE 14 CONTINUE 15 CONTINUE 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 100 FORMAT (/) END FUNCTION PDF (Iset, Ihadron, Iparton, X, Q, IR) 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 DATA 1 HUGE, DUM, D0m, D1p 1 / 1E10, 0.0, -1E-6, 1.000001 / 1 IW1, IW2 / 2*0 / Ier = 0 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 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. Return Endif Jp=Abs(Iparton) If (Jp.eq.1 .or. Jp.eq.2) Then Ipartner=3-Jp If (Iparton.lt.0) Ipartner=-Ipartner If (Ihadron.eq.1) then Tem= PdfPrt(Iset, Iparton, X, Q, Ir) Elseif (Ihadron.eq.-1) then Tem= PdfPrt(Iset, -Iparton, X, Q, Ir) Elseif (Ihadron.eq.0) then Tem= PdfPrt(Iset, Ipartner, X, Q, Ir) Elseif (Ihadron.eq.2 .or.Ihadron.eq.4 .or.Ihadron.eq.5) then Tem=( PdfPrt(Iset, Iparton, X, Q, Ir) > +PdfPrt(Iset, Ipartner, X, Q, Ir) )/2. Elseif (Ihadron.eq.6) then Tem=( 26.* PdfPrt(Iset, Iparton, X, Q, Ir) > +30.* PdfPrt(Iset, Ipartner, X, Q, Ir) )/56. Endif Else Tem= PdfPrt(Iset, Iparton, X, Q, Ir) Endif IF (TEM .LT. D0 .and. Iknl .ge. 1) Then IF (TEM .LT. D0m .AND. X .LE. 0.9) 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 PDF = TEM IR = Ier RETURN END FUNCTION PdfPrt (IPDF, LPRTN, XD, QD, IR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXPDF = 20, MXHDRN = 6, MXF = 6, MXPN = MXF*2+2) CHARACTER MSG*75 COMMON / IOUNIT / NIN, NOUT, NWRT DATA IW1 / 0 / IER = 0 Iprtn = LPRTN x = XD Q = QD If (Ipdf .eq. 0) Then Tmp = PdfTst (Iprtn, X, Q) ElseIf (Ipdf.ge.1 .and. Ipdf.le.3) then Nset = Ipdf Tmp = Ctq3Pd(Nset, Iprtn, X, Q, ir) Elseif (Ipdf .le. 9) then Nset = Ipdf - 3 Tmp = Ctq2Pd(Nset, Iprtn, X, Q, ir) Elseif (Ipdf .eq. 10) then Tmp = ParDis (Iprtn, X, Q) Elseif (Ipdf .eq. 11) then Tmp = PrDis1 (Iprtn, X, Q) Elseif (Ipdf.ge.21 .and. Ipdf.le.26) then Nset = Ipdf-20 tmp = Ctq1Pd(Nset, Iprtn, X, Q, ir) Else Ier = 1 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 END Function SetPdf () IMPLICIT DOUBLE PRECISION (A-H, O-Z) External DatPDF SetPDF = 0. Return END SUBROUTINE EVOLVE (FINI, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER*(*) HEADER PARAMETER (MXX = 105, 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) 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, KF COMMON / PEVLDT / UPD(MXPQX) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / VARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) DIMENSION QRKP(MXF) DIMENSION JI(-MXF : MXF+1) EXTERNAL NSRHSP, NSRHSM, FINI DATA D0 / 0.0 / 11 IRET = 0 IF (IHDN .LE. 4) THEN MXVAL = 2 ElseIF (IHDN .LE. 6) THEN MXVAL = 3 EndIf IF (.NOT. LSTX) CALL XARRAY DLX = 1D0 / NX J10 = NX / 10 CALL PARPDF (2, 'ALAM', AL, IR) CALL QARRAY (NINI, NFMX) NFSN = NFMX + 1 DO 101 IFLV = -NFMX, NFMX+1 JFL = NFMX + IFLV JI(IFLV) = JFL * (NT+1) * (NX+1) 101 CONTINUE 3 DO 31 IZ = 1, NX UPD(JI(0)+IZ+1) = FINI (0, XV(IZ)) UPD(JI(NFSN)+IZ+1) = 0 IF (NFMX .EQ. 0) GOTO 31 DO 331 IFLV = 1, NINI A = FINI ( IFLV, XV(IZ)) B = FINI (-IFLV, XV(IZ)) QRKP (IFLV) = A + B UPD(JI(NFSN)+IZ+1) = UPD(JI(NFSN)+IZ+1) + QRKP (IFLV) UPD(JI(-IFLV)+IZ+1) = A - B 331 CONTINUE DO 332 IFLV = 1, NINI UPD(JI( IFLV)+IZ+1) = QRKP(IFLV) - UPD(JI(NFSN)+IZ+1)/NINI 332 CONTINUE 31 CONTINUE DO 21 NEFF = NINI, NFMX IF (IKNL .EQ. 2) CALL STUPKL (NEFF) ICNT = NEFF - NINI + 1 IF (NTN(ICNT) .EQ. 0) GOTO 21 NITR = NTN (ICNT) DT = DTN (ICNT) TIN = TLN (ICNT) 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)) IF (NEFF .EQ. 0) GOTO 88 5 DO 333 IFLV = 1, NEFF CALL NSEVL (NSRHSP, IKNL, NX, NITR, JT, DT, TIN, NEFF, > UPD(JI( IFLV)+2), UPD(JI( IFLV)+1)) IF (IFLV .LE. MXVAL) > CALL NSEVL (NSRHSM, IKNL, NX, NITR, JT, DT, TIN, NEFF, > UPD(JI(-IFLV)+2), UPD(JI(-IFLV)+1)) 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 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 IF (NFMX .EQ. NEFF) GOTO 21 DO 335 IFLV = -NFMX, NFMX+1 JI(IFLV) = JI(IFLV) + NITR * (NX+1) 335 CONTINUE CALL HQRK (NX, TT, NEFF+1, UPD(JI(0)+2), UPD(JI(NEFF+1)+2)) DO 32 IZ = 1, NX QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1) UPD (JI(NFSN)+IZ+1) = UPD (JI(NFSN)+IZ+1) + QRKP (NEFF+1) VS00 = UPD (JI(NFSN)+IZ+1) / (NEFF+1) UPD(JI( NEFF+1) + IZ + 1) = QRKP(NEFF+1) - VS00 DO 321 IFL = 1, NEFF A = UPD(JI( IFL)+IZ+1) B = UPD(JI(-IFL)+IZ+1) QRKP(IFL) = A + B UPD(JI( IFL)+IZ+1) = QRKP(IFL) - VS00 IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1) = A - B 321 CONTINUE 32 CONTINUE 21 CONTINUE GOTO 900 ENTRY EVLWT (NU, HEADER, IRR) 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) WRITE (NU, IOSTAT=IR1) > HEADER, AL, FL, DR, AM1, AM2, AM3, AM4, AM5, AM6, > IPD0, IHDN, IKNL, KF, QINI, QMAX, NT, JT, NG, XMIN, XCR, NX, > (TLN(I), NTL(I), DTN(I), NTN(I), I =1, NG), NTL(NG+1), > (XV(I), I =1, NX), (QV(I), I =0, NT), (TV(I), I =0, NT) CALL WTUPD (UPD, (NX+1)*(NT+1)*KF+1, NU, IR2) IRR = IR1 + IR2 IF (IRR .NE. 0) PRINT *, 'Write error in EVLWT' RETURN ENTRY EVLRD (NU, HEADER, IRR) READ (NU, IOSTAT=IR1) > HEADER, AL, FL, DR, AM1, AM2, AM3, AM4, AM5, AM6, > IPD0, IHDN, IKNL, KF, QINI, QMAX, NT, JT, NG, XMIN, XCR, NX, > (TLN(I), NTL(I), DTN(I), NTN(I), I =1, NG), NTL(NG+1), > (XV(I), I =1, NX), (QV(I), I =0, NT), (TV(I), I =0, NT) CALL RDUPD (UPD, (NX+1)*(NT+1)*KF+1, NU, IR2) IRR = IR1 + IR2 CLOSE (NU) IF (IRR .NE. 0) PRINT *, 'Read error in EVLRD' CALL PARPDF (1, 'ALAM', AL, IR) CALL PARPDF (1, 'NFL', FL, IR) CALL PARPDF (1, 'ORDER', DR, IR) CALL PARPDF (1, 'M1', AM1, IR) CALL PARPDF (1, 'M2', AM2, IR) CALL PARPDF (1, 'M3', AM3, IR) CALL PARPDF (1, 'M4', AM4, IR) CALL PARPDF (1, 'M5', AM5, IR) CALL PARPDF (1, 'M6 ', AM6, IR) PRINT '(/A/1X,A/)', ' Parameters from this data file are:',HEADER CALL PARPDF(4, 'ALL', DBLE(NOUT), IR) 900 CONTINUE RETURN END FUNCTION PdfTst (IPRTN, XX, QQ) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Print *, >'You just called a dummy test function PdfTst in PrtPdf!' PdfTst = 0 RETURN END FUNCTION Ctq3Pd (Iset, Iparton, X, Q, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Ifl = Iparton JFL = ABS(Ifl) IF (Ifl.Eq.1 .or. Ifl.Eq.2) THEN VL = Ctq3df(Iset, Ifl, X, Q, Irt) ELSE VL = 0. ENDIF SEA = Ctq3df (Iset, -JFL, X, Q, Irt) Ctq3pd = (VL + SEA) / X Return END FUNCTION Ctq2Pd (Iset, Iparton, X, Q, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Ifl = Iparton JFL = ABS(Ifl) IF (Ifl.Eq.1 .or. Ifl.Eq.2) THEN VL = Ctq2df(Iset, Ifl, X, Q, Irt) ELSE VL = 0. ENDIF SEA = Ctq2df (Iset, -JFL, X, Q, Irt) Ctq2pd = (VL + SEA) / X Return END FUNCTION PARDIS (IPRTN, XX, QQ) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF COMMON / PEVLDT / UPD(MXPQX) Dimension Fq(5), Df(5), QL(5) DATA SMLL, M, II / 1.D-7, 2, 0 / X = XX Q = QQ Md = M / 2 Amd= Md IF (X .LT. -SMLL .OR. X .GT. 1D0+SMLL) THEN PRINT *, 'X = ',X, 'Out of range in the PARDIS function call!' STOP EndIf IF (Q .LT. QINI) THEN CALL WARNR (IWRN1, NWRT, 'Q less than QINI in PARDIS call; Q SET >= QINI.', 'Q', Q, QINI, QMAX, 1) Q = QINI ElseIF (Q .GT. QMAX) THEN CALL WARNR(IWRN2,NWRT,'Q greater than QMAX in PARDIS call; >Extrapolation will be used.', 'Q', Q, QINI, QMAX, 1) EndIf 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 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) JFL = IPRTN + (KF - 2) / 2 J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx 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 Pardis; II = ', II Endif 21 Continue Call Polint (TV(Jq), Fq(1), M+1, TT, Ftmp, Ddf) PARDIS = Ftmp RETURN END FUNCTION PRDIS1 (IPRTN, XX, QQ) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF COMMON / P1VLDT / UPD(MXPQX) Dimension Fq(5), Df(5), QL(5) DATA SMLL, M, II / 1.D-7, 2, 0 / X = XX Q = QQ Md = M / 2 Amd= Md IF (X .LT. -SMLL .OR. X .GT. 1D0+SMLL) THEN PRINT *, 'X = ',X, 'Out of range in the PRDIS1 function call!' STOP EndIf IF (Q .LT. QINI) THEN CALL WARNR (IWRN1, NWRT, 'Q less than QINI in PRDIS1 call; Q SET >= QINI.', 'Q', Q, QINI, QMAX, 1) Q = QINI ElseIF (Q .GT. QMAX) THEN CALL WARNR(IWRN2,NWRT,'Q greater than QMAX in PRDIS1 call; >Extrapolation will be used.', 'Q', Q, QINI, QMAX, 1) EndIf 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 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) JFL = IPRTN + (KF - 2) / 2 J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx 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 RETURN END FUNCTION Ctq1Pd (Iset, Iparton, X, Q, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Ifl = Iparton JFL = ABS(Ifl) IF (Ifl .LE. 0) THEN VL = 0 ELSEIF (Ifl .LE. 2) THEN VL = CtqPdf(Iset, Ifl, X, Q, Irt) ELSE VL = 0 ENDIF SEA = CtqPdf (Iset, -JFL, X, Q, Irt) Ctq1Pd = (VL + SEA) / X Return END SUBROUTINE EVLPAR (IACT, NAME, VALUE, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER*(*) NAME COMMON / IOUNIT / NIN, NOUT, NWRT 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 : Kf ' /) 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 RETURN END SUBROUTINE XARRAY IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0 = 0.0, D10=10.0) PARAMETER (MXX = 105, 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) 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) DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG) 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 NA, IA, PUNY / 3, 5, 1D-30 / XV(0) = 0D0 DZ = 1D0 / (NX-1) DO 10 I = 1, NX - 1 Z = DZ * (I-1) X = XFRMZ (Z) DXTZ(I) = DXDZ(Z) / X 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 XV(NX) = 1D0 DXTZ(NX) = DXDZ(1.D0) DO 21 L = L1, L2 XA (NX, L) = 1D0 21 CONTINUE XA (NX, 0) = 0D0 DO 11 I = 1, NX-1 ELY(I) = LOG(1D0 - XV(I)) 11 CONTINUE ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3) 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 DO 30 I = 1, NX-1 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) 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 CALL WARNR(IWRN, NWRT, 'Determinant close to zero; > will be arbitrarily set to:', 'DET', PUNY, D0, D0, 0) DET = PUNY EndIf G2(1,2) = (F22*F33 - F23*F32) / DET G2(1,3) = (F32*F13 - F33*F12) / DET G2(1,4) = (F12*F23 - F13*F22) / DET G2(2,2) = (F23*F31 - F21*F33) / DET G2(2,3) = (F33*F11 - F31*F13) / DET G2(2,4) = (F13*F21 - F11*F23) / DET G2(3,2) = (F21*F32 - F22*F31) / DET G2(3,3) = (F31*F12 - F32*F11) / DET G2(3,4) = (F11*F22 - F12*F21) / DET 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 DO 51 J = 1, NDH DO 52 L = 1, NDG IF (I .EQ. 1) THEN GB(L,J,I) = G2(L,J) ElseIF (I .EQ. NX-1) THEN GB(L,J,I) = G1(L,J) Else GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0 EndIf 52 CONTINUE 51 CONTINUE 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 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 DO 42 J = 1, NDG DO 44 L = 1, NDG G1(L,J) = G2(L,J+1) 44 CONTINUE 42 CONTINUE 30 CONTINUE LSTX = .TRUE. RETURN END SUBROUTINE QARRAY (NINI, NFMX) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF NCNT = 0 IF (NT .GE. mxq) NT = mxq - 1 S = LOG(QINI/AL) TINI = LOG(S) S = LOG(QMAX/AL) TMAX = LOG(S) 1 DT0 = (TMAX - TINI) / float(NT) NINI = NFL(QINI) NFMX = NFL(QMAX) KF = 2 * NFMX + 2 NG = NFMX - NINI + 1 QIN = QINI QOUT = QINI S = LOG(QIN/AL) TIN = LOG(S) TLN(1) = TIN NTL(1) = 0 QV(0) = QINI TV(0) = Tin DO 20 NEFF = NINI, NFMX ICNT = NEFF - NINI + 1 IF (NEFF .LT. NFMX) THEN THRN = AMHATF (NEFF + 1) QOUN = MIN (QMAX, THRN) Else QOUN = QMAX EndIf IF (QOUN-QOUT .LE. 0.0001) THEN DT = 0 NITR = 0 Else QOUT = QOUN S = LOG(QOUT/AL) TOUT = LOG(S) TEM = TOUT - TIN NITR = INT (TEM / DT0) + 1 DT = TEM / NITR EndIf DTN (ICNT) = DT NTN (ICNT) = NITR TLN (ICNT) = TIN NTL (ICNT+1) = NTL(ICNT) + NITR 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 QIN = QOUT TIN = TOUT 20 CONTINUE NCNT = NCNT + 1 NTP = NTL (NG + 1) ND = NTP - NT IF (NTP .GE. MXQ) THEN NT = MXQ - ND - NCNT GOTO 1 EndIf NT = NTP RETURN END SUBROUTINE STUPKL (NFL) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MX = 3) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ (MXX, 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, TFF1, RFF1, VFF1 EXTERNAL FNSP, FNSM SAVE DATA AERR, RERR / 0.0, 0.02 / NNX = NX DZ = 1./ (NX - 1) DO 5 I0 = 1, MX XL(I0) = XV(I0) 5 CONTINUE 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 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) 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) 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) 20 CONTINUE 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) AGF2 = GAUSS(RGF2, XV(NX-1), D1, AER, RER, ERR, IRT) + AGF2 AFG2 = GAUSS(RFG2, XV(NX-1), D1, AER, RER, ERR, IRT) + AFG2 DO 21 IX = 1, NX-1 X = XV(IX) NP = NX - IX + 1 IS = NP 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) TZ = (Z + DZ) / DZ IZ = TZ IZ = MAX (IZ, 0) IZ = MIN (IZ, NX-1) DT = TZ - IZ 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 END SUBROUTINE SNEVL(IKNL, NX, NT, JT, DT, TIN, NEFF, UI, GI, US, GS) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 105, 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) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / ANOMS1 / GG, GF, FG, FF, GPL, GMI COMMON / ANOMS2 / SGG, SGF, SFG, SFF, TGG, TGF, TFG, TFF COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) 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 / ESF(T) = 4./ (33.- 2.*NEFF) * LOG((T - TLAM)/(TIN - TLAM)) N5 = (NX - 1) / 5 JTT = 2 * JT DDT = DT / JTT IF (NX .GT. MXX) THEN WRITE (NOUT,*) 'Nx =', NX, ' greater than Max # of pts in SNEVL.' STOP 'Program stopped in SNEVL' EndIf TMD = TIN + DT * NT / 2. AMU = EXP(TMD) TEM = 6./ (33.- 2.* NEFF) / ALPI(AMU) TLAM = TMD - TEM 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) TT = TIN DO 10 IZ = 1, NX Y0(IZ) = UI(IZ) Z0(IZ) = GI(IZ) 10 CONTINUE DO 20 IS = 1, NT DO 202 JS = 1, JTT IRND = (IS-1) * JTT + JS IF (IRND .EQ. 1) THEN CALL SNRHS (TT, NEFF, Y0, Z0, F0, G0) DO 250 IZ = 1, NX Y0(IZ) = Y0(IZ) + DDT * F0(IZ) Z0(IZ) = Z0(IZ) + DDT * G0(IZ) 250 CONTINUE TT = TT + DDT CALL SNRHS (TT, NEFF, Y0, Z0, F1, G1) 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 Else CALL SNRHS (TT, NEFF, Y1, Z1, F1, G1) 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 TT = TT + DDT CALL SNRHS (TT, NEFF, YP, ZP, FP, GP) 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 202 CONTINUE 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 US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3) 20 CONTINUE RETURN END SUBROUTINE NSRHSP (TT, NEFF, FI, FO) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ (MXX, 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, KF DIMENSION G1(MXX), FI(NX), FO(NX) DIMENSION W0(MXX), W1(MXX), WH(MXX) DATA IRUN / 0 / S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S CALL INTEGR (NX, 0, FI, W0, IR1) CALL INTEGR (NX, 1, FI, W1, IR2) CALL HINTEG (NX, FI, WH) DO 230 IZ = 1, NX FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ)) FO(IZ) = CPL * FO(IZ) 230 CONTINUE IF (IKNL .EQ. 2) THEN DZ = 1./ (NX - 1) DO 21 IX = 1, NX-1 X = XV(IX) NP = NX - IX + 1 DO 31 KZ = 2, NP IY = IX + KZ - 1 XY = ZZ (NX-IX+1, NX-IY+1) G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX)) 31 CONTINUE TEM1 = SMPNOL (NP, DZ, G1, ERR) TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + AQQB)) * CPL2 FO(IX) = FO(IX) + TMP2 21 CONTINUE EndIf RETURN END SUBROUTINE NSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 105, 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) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / ANOMS1 / GG, GF, FG, FF, GPL, GMI COMMON / ANOMS2 / SGG, SGF, SFG, SFF, TGG, TGF, TFG, TFF COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) DIMENSION U0(NX), UN(0:NX, 0:NT) DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX) external rhs ESF(T) = 4./ (33.- 2.*NEFF) * LOG((T - TLAM)/(TIN - TLAM)) N5 = (NX - 1) / 5 DDT = DT / JT IF (NX .GT. MXX) THEN WRITE (NOUT,*) 'Nx =', NX, ' greater than Max # of pts in NSEVL.' STOP 'Program stopped in NSEVL' EndIf TMD = TIN + DT * NT / 2. AMU = EXP(TMD) TEM = 6./ (33.- 2.* NEFF) / ALPI(AMU) TLAM = TMD - TEM DO 9 IX = 1, NX UN(IX, 0) = U0(IX) 9 CONTINUE UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1) TT = TIN DO 10 IZ = 1, NX Y0(IZ) = U0(IZ) 10 CONTINUE DO 20 IS = 1, NT DO 202 JS = 1, JT IRND = (IS-1) * JT + JS IF (IRND .EQ. 1) THEN CALL RHS (TT, Neff, Y0, F0) DO 250 IZ = 1, NX Y0(IZ) = Y0(IZ) + DDT * F0(IZ) 250 CONTINUE TT = TT + DDT CALL RHS (TT, NEFF, Y0, F1) DO 251 IZ = 1, NX Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0 251 CONTINUE Else CALL RHS (TT, NEFF, Y1, F1) DO 252 IZ = 1, NX YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0 252 CONTINUE TT = TT + DDT CALL RHS (TT, NEFF, YP, FP) DO 253 IZ = 1, NX Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0 F0(IZ) = F1(IZ) 253 CONTINUE EndIf 202 CONTINUE DO 260 IZ = 1, NX UN (IZ, IS) = Y1(IZ) 260 CONTINUE UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) 20 CONTINUE RETURN END SUBROUTINE NSRHSM (TT, NEFF, FI, FO) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ (MXX, 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, KF DIMENSION G1(MXX), FI(NX), FO(NX) DIMENSION W0(MXX), W1(MXX), WH(MXX) DATA IRUN / 0 / S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S CALL INTEGR (NX, 0, FI, W0, IR1) CALL INTEGR (NX, 1, FI, W1, IR2) CALL HINTEG (NX, FI, WH) DO 230 IZ = 1, NX FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ)) FO(IZ) = CPL * FO(IZ) 230 CONTINUE IF (IKNL .EQ. 2) THEN DZ = 1./ (NX - 1) DO 21 IX = 1, NX-1 X = XV(IX) NP = NX - IX + 1 IS = NP DO 31 KZ = 2, NP IY = IX + KZ - 1 IT = NX - IY + 1 XY = ZZ (IS, IT) G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX)) 31 CONTINUE TEM1 = SMPNOL (NP, DZ, G1, ERR) TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2 FO(IX) = FO(IX) + TMP2 21 CONTINUE EndIf RETURN END SUBROUTINE HQRK (NX, TT, NQRK, Y, F) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 105, MXQ = 25, MXF = 6) DIMENSION Y(NX), F(NX) IF (NX .GT. 1) GOTO 11 11 CONTINUE DO 230 IZ = 1, NX IF (NX .GT. 1) THEN F(IZ) = 0 GOTO 230 EndIf 230 CONTINUE RETURN END SUBROUTINE WTUPD (UPD, NTL, NDAT, IRET) DOUBLE PRECISION UPD DIMENSION UPD (NTL) WRITE (NDAT, IOSTAT=IRET) UPD RETURN ENTRY RDUPD (UPD, NTL, NDAT, IRET) READ (NDAT, IOSTAT=IRET) UPD RETURN END FUNCTION Ctq3df (Iset, Iprtn, XX, QQ, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) Parameter (Nst = 3) DIMENSION > Iord(Nst), Isch(Nst), Nqrk(Nst),Alm(Nst) > , Vlm(4:6,Nst), Qms(4:6, Nst) > , Xmn(Nst), Qmn(Nst), Qmx(Nst) DATA > Isch(1), Iord(1), Nqrk(1), Alm(1) / 1, 2, 6, .239 / > (Vlm(I,1), I=4,6) / .239, .158, .063 / > (Qms(I,1), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(1), Qmn(1), Qmx(1) / 1.E-6, 1.60, 1.E4 / DATA > Isch(2), Iord(2), Nqrk(2), Alm(2) / 1, 1, 6, .177 / > (Vlm(I,2), I=4,6) / .177, .132, .066 / > (Qms(I,2), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(2), Qmn(2), Qmx(2) / 1.E-6, 1.60, 1.E4 / DATA > Isch(3), Iord(3), Nqrk(3), Alm(3) / 1, 2, 6, .247 / > (Vlm(I,3), I=4,6) / .247, .164, .066 / > (Qms(I,3), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(3), Qmn(3), Qmx(3) / 1.E-6, 1.60, 1.E4 / Data Ist, Lp, Qsto / 0, -10, 1.2345 / save Ist, Lp, Qsto save SB, SB2, SB3 X = XX Irt = 0 if(Iset.eq.Ist .and. Qsto.eq.QQ) then if (Iprtn.eq.Lp) goto 100 if (Iprtn.ge.-3 .and. Lp.ge.-3) goto 501 endif Ip = abs(Iprtn) If (Ip .GE. 4) then If (QQ .LE. Qms(Ip, Iset)) Then Ctq3df = 0.0 Return Endif Qi = Qms(ip, Iset) Else Qi = Qmn(Iset) Endif Alam = Alm (Iset) SBL = LOG(QQ/Alam) / LOG(Qi/Alam) SB = LOG (SBL) SB2 = SB*SB SB3 = SB2*SB 501 Iflv = 3 - Iprtn Goto (1,2,3, 311) Iset 1 Goto(11,12,13,14,15,16,17,18,19)Iflv 11 A0=Exp(-0.7266E+00-0.1584E+01*SB +0.1259E+01*SB2-0.4305E-01*SB3) A1= 0.5285E+00-0.3721E+00*SB +0.5150E+00*SB2-0.1697E+00*SB3 A2= 0.4075E+01+0.8282E+00*SB -0.4496E+00*SB2+0.2107E+00*SB3 A3= 0.3279E+01+0.5066E+01*SB -0.9134E+01*SB2+0.2897E+01*SB3 A4= 0.4399E+00-0.5888E+00*SB +0.4802E+00*SB2-0.1664E+00*SB3 A5= 0.3678E+00-0.8929E+00*SB +0.1592E+01*SB2-0.5713E+00*SB3 goto 100 12 A0=Exp( 0.2259E+00+0.1237E+00*SB +0.3035E+00*SB2-0.2935E+00*SB3) A1= 0.5085E+00+0.1651E-01*SB -0.3592E-01*SB2+0.2782E-01*SB3 A2= 0.3732E+01+0.4901E+00*SB +0.2218E+00*SB2-0.1116E+00*SB3 A3= 0.7011E+01-0.6620E+01*SB +0.2557E+01*SB2-0.1360E+00*SB3 A4= 0.8969E+00-0.2429E+00*SB +0.1811E+00*SB2-0.6888E-01*SB3 A5= 0.8636E-01+0.2558E+00*SB -0.3082E+00*SB2+0.2535E+00*SB3 goto 100 13 A0=Exp(-0.2318E+00-0.9779E+00*SB -0.3783E+00*SB2+0.1037E-01*SB3) A1=-0.2916E+00+0.1754E+00*SB -0.1884E+00*SB2+0.6116E-01*SB3 A2= 0.5349E+01+0.7460E+00*SB +0.2319E+00*SB2-0.2622E+00*SB3 A3= 0.6920E+01-0.3454E+01*SB +0.2027E+01*SB2-0.7626E+00*SB3 A4= 0.1013E+01+0.1423E+00*SB -0.1798E+00*SB2+0.1872E-01*SB3 A5=-0.5465E-01+0.2303E+01*SB -0.9584E+00*SB2+0.3098E+00*SB3 goto 100 14 A0=Exp(-0.2906E+01-0.1069E+00*SB -0.1055E+01*SB2+0.2496E+00*SB3) A1=-0.2875E+00+0.6571E-01*SB -0.1987E-01*SB2-0.1800E-02*SB3 A2= 0.9854E+01-0.2715E+00*SB -0.7407E+00*SB2+0.2888E+00*SB3 A3= 0.1583E+02-0.7687E+01*SB +0.3428E+01*SB2-0.3327E+00*SB3 A4= 0.9763E+00+0.7599E-01*SB -0.2128E+00*SB2+0.6852E-01*SB3 A5=-0.8444E-02+0.9434E+00*SB +0.4152E+00*SB2-0.1481E+00*SB3 goto 100 15 A0=Exp(-0.2328E+01-0.3061E+01*SB +0.3620E+01*SB2-0.1602E+01*SB3) A1=-0.3358E+00+0.3198E+00*SB -0.4210E+00*SB2+0.1571E+00*SB3 A2= 0.8478E+01-0.3112E+01*SB +0.5243E+01*SB2-0.2255E+01*SB3 A3= 0.1971E+02+0.3389E+00*SB -0.5268E+01*SB2+0.2099E+01*SB3 A4= 0.1128E+01-0.4701E+00*SB +0.7779E+00*SB2-0.3506E+00*SB3 A5=-0.4708E+00+0.3341E+01*SB -0.3375E+01*SB2+0.1353E+01*SB3 goto 100 16 A0=Exp(-0.3780E+01+0.2499E+01*SB -0.4962E+01*SB2+0.1936E+01*SB3) A1=-0.2639E+00-0.1575E+00*SB +0.3584E+00*SB2-0.1646E+00*SB3 A2= 0.8082E+01+0.2794E+01*SB -0.5438E+01*SB2+0.2321E+01*SB3 A3= 0.1811E+02-0.2000E+02*SB +0.1951E+02*SB2-0.6904E+01*SB3 A4= 0.9822E+00+0.4972E+00*SB -0.8690E+00*SB2+0.3415E+00*SB3 A5= 0.1772E+00-0.6078E+00*SB +0.3341E+01*SB2-0.1473E+01*SB3 goto 100 17 A0=SB** 0.1122E+01*Exp(-0.4232E+01-0.1808E+01*SB +0.5348E+00*SB2) A1=-0.2824E+00+0.5846E+00*SB -0.7230E+00*SB2+0.2419E+00*SB3 A2= 0.5683E+01-0.2948E+01*SB +0.5916E+01*SB2-0.2560E+01*SB3 A3= 0.2051E+01+0.4795E+01*SB -0.4271E+01*SB2+0.4174E+00*SB3 A4= 0.1737E+00+0.1717E+01*SB -0.1978E+01*SB2+0.6643E+00*SB3 A5= 0.8689E+00+0.3500E+01*SB -0.3283E+01*SB2+0.1026E+01*SB3 goto 100 18 A0=SB** 0.9906E+00*Exp(-0.1496E+01-0.6576E+01*SB +0.1569E+01*SB2) A1=-0.2140E+00-0.6419E-01*SB -0.2741E-02*SB2+0.3185E-02*SB3 A2= 0.5781E+01+0.1049E+00*SB -0.3930E+00*SB2+0.5174E+00*SB3 A3=-0.9420E+00+0.5511E+00*SB +0.8817E+00*SB2+0.1903E+01*SB3 A4= 0.2418E-01+0.4232E-01*SB -0.1244E-01*SB2-0.2365E-01*SB3 A5= 0.7664E+00+0.1794E+01*SB -0.4917E+00*SB2-0.1284E+00*SB3 goto 100 19 A0=SB** 0.1000E+01*Exp(-0.8460E+01+0.1154E+01*SB +0.8838E+01*SB2) A1=-0.4316E-01-0.2976E+00*SB +0.3174E+00*SB2-0.1429E+01*SB3 A2= 0.4910E+01+0.2273E+01*SB +0.5631E+01*SB2-0.1994E+02*SB3 A3= 0.1190E+02-0.2000E+02*SB -0.2000E+02*SB2+0.1292E+02*SB3 A4= 0.5771E+00-0.2552E+00*SB +0.7510E+00*SB2+0.6923E+00*SB3 A5= 0.4402E+01-0.1627E+01*SB -0.2085E+01*SB2-0.6737E+01*SB3 goto 100 2 Goto(21,22,23,24,25,26,27,28,29)Iflv 21 A0=Exp( 0.1141E+00+0.4764E+00*SB -0.1745E+01*SB2+0.7728E+00*SB3) A1= 0.4275E+00-0.1290E+00*SB +0.3609E+00*SB2-0.1689E+00*SB3 A2= 0.3000E+01+0.2946E+01*SB -0.4117E+01*SB2+0.1989E+01*SB3 A3=-0.1302E+01+0.2322E+01*SB -0.4258E+01*SB2+0.2109E+01*SB3 A4= 0.2586E+01-0.1920E+00*SB -0.3754E+00*SB2+0.2731E+00*SB3 A5=-0.2251E+00-0.5374E+00*SB +0.2245E+01*SB2-0.1034E+01*SB3 goto 100 22 A0=Exp( 0.1907E+00+0.4205E-01*SB +0.2752E+00*SB2-0.3171E+00*SB3) A1= 0.4611E+00+0.2331E-01*SB -0.3403E-01*SB2+0.3174E-01*SB3 A2= 0.3504E+01+0.5739E+00*SB +0.2676E+00*SB2-0.1553E+00*SB3 A3= 0.7452E+01-0.6742E+01*SB +0.2849E+01*SB2-0.1964E+00*SB3 A4= 0.1116E+01-0.3435E+00*SB +0.2865E+00*SB2-0.1288E+00*SB3 A5= 0.6659E-01+0.2714E+00*SB -0.2688E+00*SB2+0.2763E+00*SB3 goto 100 23 A0=Exp(-0.7631E+00-0.7241E+00*SB -0.1170E+01*SB2+0.5343E+00*SB3) A1=-0.3573E+00+0.3469E+00*SB -0.3396E+00*SB2+0.9188E-01*SB3 A2= 0.5604E+01+0.7458E+00*SB -0.5082E+00*SB2+0.1844E+00*SB3 A3= 0.1549E+02-0.1809E+02*SB +0.1162E+02*SB2-0.3483E+01*SB3 A4= 0.9881E+00+0.1364E+00*SB -0.4421E+00*SB2+0.2051E+00*SB3 A5=-0.9505E-01+0.3259E+01*SB -0.1547E+01*SB2+0.2918E+00*SB3 goto 100 24 A0=Exp(-0.2740E+01-0.7987E-01*SB -0.9015E+00*SB2-0.9872E-01*SB3) A1=-0.3909E+00+0.1244E+00*SB -0.4487E-01*SB2+0.1277E-01*SB3 A2= 0.9163E+01+0.2823E+00*SB -0.7720E+00*SB2-0.9360E-02*SB3 A3= 0.1080E+02-0.3915E+01*SB -0.1153E+01*SB2+0.2649E+01*SB3 A4= 0.9894E+00-0.1647E+00*SB -0.9426E-02*SB2+0.2945E-02*SB3 A5=-0.3395E+00+0.6998E+00*SB +0.7000E+00*SB2-0.6730E-01*SB3 goto 100 25 A0=Exp(-0.2449E+01-0.3513E+01*SB +0.4529E+01*SB2-0.2031E+01*SB3) A1=-0.4050E+00+0.3411E+00*SB -0.3669E+00*SB2+0.1109E+00*SB3 A2= 0.7470E+01-0.2982E+01*SB +0.5503E+01*SB2-0.2419E+01*SB3 A3= 0.1503E+02+0.1638E+01*SB -0.8772E+01*SB2+0.3852E+01*SB3 A4= 0.1137E+01-0.1006E+01*SB +0.1485E+01*SB2-0.6389E+00*SB3 A5=-0.5299E+00+0.3160E+01*SB -0.3104E+01*SB2+0.1219E+01*SB3 goto 100 26 A0=Exp(-0.3640E+01+0.1250E+01*SB -0.2914E+01*SB2+0.8390E+00*SB3) A1=-0.3595E+00-0.5259E-01*SB +0.3122E+00*SB2-0.1642E+00*SB3 A2= 0.7305E+01+0.9727E+00*SB -0.9788E+00*SB2-0.5193E-01*SB3 A3= 0.1198E+02-0.1799E+02*SB +0.2614E+02*SB2-0.1091E+02*SB3 A4= 0.9882E+00-0.6101E+00*SB +0.9737E+00*SB2-0.4935E+00*SB3 A5=-0.1186E+00-0.3231E+00*SB +0.3074E+01*SB2-0.1274E+01*SB3 goto 100 27 A0=SB** 0.1122E+01*Exp(-0.3718E+01-0.1335E+01*SB +0.1651E-01*SB2) A1=-0.4719E+00+0.7509E+00*SB -0.8420E+00*SB2+0.2901E+00*SB3 A2= 0.6194E+01-0.1641E+01*SB +0.4907E+01*SB2-0.2523E+01*SB3 A3= 0.4426E+01-0.4270E+01*SB +0.6581E+01*SB2-0.3474E+01*SB3 A4= 0.2683E+00+0.9876E+00*SB -0.7612E+00*SB2+0.1780E+00*SB3 A5=-0.4547E+00+0.4410E+01*SB -0.3712E+01*SB2+0.1245E+01*SB3 goto 100 28 A0=SB** 0.9838E+00*Exp(-0.2548E+01-0.7660E+01*SB +0.3702E+01*SB2) A1=-0.3122E+00-0.2120E+00*SB +0.5716E+00*SB2-0.3773E+00*SB3 A2= 0.6257E+01-0.8214E-01*SB -0.2537E+01*SB2+0.2981E+01*SB3 A3=-0.6723E+00+0.2131E+01*SB +0.9599E+01*SB2-0.7910E+01*SB3 A4= 0.9169E-01+0.4295E-01*SB -0.5017E+00*SB2+0.3811E+00*SB3 A5= 0.2402E+00+0.2656E+01*SB -0.1586E+01*SB2+0.2880E+00*SB3 goto 100 29 A0=SB** 0.1001E+01*Exp(-0.6934E+01+0.3050E+01*SB -0.6943E+00*SB2) A1=-0.1713E+00-0.5167E+00*SB +0.1241E+01*SB2-0.1703E+01*SB3 A2= 0.6169E+01+0.3023E+01*SB -0.1972E+02*SB2+0.1069E+02*SB3 A3= 0.4439E+01-0.1746E+02*SB +0.1225E+02*SB2+0.8350E+00*SB3 A4= 0.5458E+00-0.4586E+00*SB +0.9089E+00*SB2-0.4049E+00*SB3 A5= 0.3207E+01-0.3362E+01*SB +0.5877E+01*SB2-0.7659E+01*SB3 goto 100 3 Goto(31,32,33,34,35,36,37,38,39)Iflv 31 A0=Exp( 0.3961E+00+0.4914E+00*SB -0.1728E+01*SB2+0.7257E+00*SB3) A1= 0.4162E+00-0.1419E+00*SB +0.3680E+00*SB2-0.1618E+00*SB3 A2= 0.3248E+01+0.3028E+01*SB -0.4307E+01*SB2+0.1920E+01*SB3 A3=-0.1100E+01+0.2184E+01*SB -0.3820E+01*SB2+0.1717E+01*SB3 A4= 0.2082E+01-0.2756E+00*SB +0.3043E+00*SB2-0.1260E+00*SB3 A5=-0.4822E+00-0.5706E+00*SB +0.2243E+01*SB2-0.9760E+00*SB3 goto 100 32 A0=Exp( 0.2148E+00+0.5814E-01*SB +0.2734E+00*SB2-0.2902E+00*SB3) A1= 0.4810E+00+0.1657E-01*SB -0.3800E-01*SB2+0.3125E-01*SB3 A2= 0.3509E+01+0.3923E+00*SB +0.4010E+00*SB2-0.1932E+00*SB3 A3= 0.7055E+01-0.6552E+01*SB +0.3466E+01*SB2-0.5657E+00*SB3 A4= 0.1061E+01-0.3453E+00*SB +0.4089E+00*SB2-0.1817E+00*SB3 A5= 0.8687E-01+0.2548E+00*SB -0.2967E+00*SB2+0.2647E+00*SB3 goto 100 33 A0=Exp(-0.4665E+00-0.7554E+00*SB -0.3323E+00*SB2-0.2734E-04*SB3) A1=-0.3359E+00+0.2395E+00*SB -0.2377E+00*SB2+0.7059E-01*SB3 A2= 0.5451E+01+0.6086E+00*SB +0.8606E-01*SB2-0.1425E+00*SB3 A3= 0.1026E+02-0.9352E+01*SB +0.4879E+01*SB2-0.1150E+01*SB3 A4= 0.9935E+00-0.5017E-01*SB -0.1707E-01*SB2-0.1464E-02*SB3 A5=-0.4160E-01+0.2305E+01*SB -0.1063E+01*SB2+0.3211E+00*SB3 goto 100 34 A0=Exp(-0.3323E+01+0.2296E+00*SB -0.1109E+01*SB2+0.2223E+00*SB3) A1=-0.3410E+00+0.8847E-01*SB -0.1111E-01*SB2-0.5927E-02*SB3 A2= 0.9753E+01-0.5182E+00*SB -0.4670E+00*SB2+0.1921E+00*SB3 A3= 0.1977E+02-0.1600E+02*SB +0.9481E+01*SB2-0.1864E+01*SB3 A4= 0.9818E+00+0.2839E-02*SB -0.1188E+00*SB2+0.3584E-01*SB3 A5=-0.7934E-01+0.1004E+01*SB +0.3704E+00*SB2-0.1220E+00*SB3 goto 100 35 A0=Exp(-0.2714E+01-0.2868E+01*SB +0.3700E+01*SB2-0.1671E+01*SB3) A1=-0.3893E+00+0.3341E+00*SB -0.3897E+00*SB2+0.1420E+00*SB3 A2= 0.8359E+01-0.3267E+01*SB +0.5327E+01*SB2-0.2245E+01*SB3 A3= 0.2359E+02-0.5669E+01*SB -0.4602E+01*SB2+0.3153E+01*SB3 A4= 0.1106E+01-0.4745E+00*SB +0.7739E+00*SB2-0.3417E+00*SB3 A5=-0.5557E+00+0.3433E+01*SB -0.3390E+01*SB2+0.1354E+01*SB3 goto 100 36 A0=Exp(-0.3985E+01+0.2855E+01*SB -0.5208E+01*SB2+0.1937E+01*SB3) A1=-0.3337E+00-0.1150E+00*SB +0.3691E+00*SB2-0.1709E+00*SB3 A2= 0.7968E+01+0.3641E+01*SB -0.6599E+01*SB2+0.2642E+01*SB3 A3= 0.1873E+02-0.1999E+02*SB +0.1734E+02*SB2-0.5813E+01*SB3 A4= 0.9731E+00+0.5082E+00*SB -0.8780E+00*SB2+0.3231E+00*SB3 A5=-0.5542E-01-0.4189E+00*SB +0.3309E+01*SB2-0.1439E+01*SB3 goto 100 37 A0=SB** 0.1105E+01*Exp(-0.3952E+01-0.1901E+01*SB +0.5137E+00*SB2) A1=-0.3543E+00+0.6055E+00*SB -0.6941E+00*SB2+0.2278E+00*SB3 A2= 0.5955E+01-0.2629E+01*SB +0.5337E+01*SB2-0.2300E+01*SB3 A3= 0.1933E+01+0.4882E+01*SB -0.3810E+01*SB2+0.2290E+00*SB3 A4= 0.1806E+00+0.1655E+01*SB -0.1893E+01*SB2+0.6395E+00*SB3 A5= 0.4790E+00+0.3612E+01*SB -0.3152E+01*SB2+0.9684E+00*SB3 goto 100 38 A0=SB** 0.9818E+00*Exp(-0.1825E+01-0.7464E+01*SB +0.2143E+01*SB2) A1=-0.2604E+00-0.1400E+00*SB +0.1702E+00*SB2-0.8476E-01*SB3 A2= 0.6005E+01+0.6275E+00*SB -0.2535E+01*SB2+0.2219E+01*SB3 A3=-0.9067E+00+0.1149E+01*SB +0.1974E+01*SB2+0.4716E+01*SB3 A4= 0.3915E-01+0.5945E-01*SB -0.9844E-01*SB2+0.2783E-01*SB3 A5= 0.5500E+00+0.1994E+01*SB -0.6727E+00*SB2-0.1510E+00*SB3 goto 100 39 A0=SB** 0.1002E+01*Exp(-0.8553E+01+0.3793E+00*SB +0.9998E+01*SB2) A1=-0.5870E-01-0.2792E+00*SB +0.6526E+00*SB2-0.1984E+01*SB3 A2= 0.4716E+01+0.4473E+00*SB +0.1128E+02*SB2-0.1937E+02*SB3 A3= 0.1289E+02-0.1742E+02*SB -0.1983E+02*SB2-0.9274E+00*SB3 A4= 0.5647E+00-0.2732E+00*SB +0.1074E+01*SB2+0.5981E+00*SB3 A5= 0.4390E+01-0.1262E+01*SB -0.9026E+00*SB2-0.9394E+01*SB3 goto 100 311 stop 'This option is not currently supported.' 100 Ctq3df = A0 *(x**A1) *((D1-x)**A2) *(D1+A3*(x**A4)) $ *(log(D1+D1/x))**A5 if(Ctq3df.lt.D0) then Ctq3df = D0 Irt=1 endif Ist = Iset Lp = Iprtn Qsto = QQ Return ENTRY Wlamd3 (Iset, Iorder, Neff) Iorder = Iord (Iset) Wlamd3 = VLM (Neff, Iset) RETURN END FUNCTION Ctq2df (Iset, Iprtn, XX, QQ, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (Nex = 5, MxFl = 6, Npn = 3, Nst = 30, Nexpt=20) Parameter (Nst4 = Nst*4) DIMENSION > Iord(Nst), Isch(Nst), Nqrk(Nst),Alm(Nst) > , Vlm(4:6,Nst), Qms(4:6, Nst) > , Xmn(Nst), Qmn(Nst), Qmx(Nst), Nexp(Nexpt) > , Mex(Nst), Mpn(Nst), ExpN(Nexpt, Nst), ExpNor(Nexpt) DATA > Isch(1), Iord(1), Nqrk(1), Alm(1) / 1, 2, 6, .213 / > (Vlm(I,1), I=4,6) / .213, .139, .053 / > (Qms(I,1), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(1), Qmn(1), Qmx(1) / 1.E-5, 1.60, 1.E3 / > Mex(1), Mpn(1), Nexp(1) / 5, 3, 10 / > (ExpN(I,1), I=1,10 ) > / 0.990,1.012,1.022,0.980,1.062,0.870,0.843 > ,0.815,0.974,1.029 / DATA > Isch(2), Iord(2), Nqrk(2), Alm(2) / 1, 2, 6, .208 / > (Vlm(I,2), I=4,6) / .208, .135, .051 / > (Qms(I,2), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(2), Qmn(2), Qmx(2) / 1.E-5, 1.60, 1.E3 / > Mex(2), Mpn(2), Nexp(2) / 5, 3, 10 / > (ExpN(I,2), I=1,10 ) > / 0.992,1.017,1.023,0.982,1.079,0.879,0.845 > ,0.814,0.984,1.036 / DATA > Isch(3), Iord(3), Nqrk(3), Alm(3) / 1, 2, 6, .208 / > (Vlm(I,3), I=4,6) / .208, .135, .051 / > (Qms(I,3), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(3), Qmn(3), Qmx(3) / 1.E-5, 1.60, 1.E3 / > Mex(3), Mpn(3), Nexp(3) / 5, 3, 10 / > (ExpN(I,3), I=1,10 ) > / 0.989,1.014,1.021,0.977,1.099,0.812,0.789 > ,0.822,0.909,0.950 / DATA > Isch(4), Iord(4), Nqrk(4), Alm(4) / 1, 2, 6, .322 / > (Vlm(I,4), I=4,6) / .322, .220, .088 / > (Qms(I,4), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(4), Qmn(4), Qmx(4) / 1.E-5, 1.60, 1.E3 / > Mex(4), Mpn(4), Nexp(4) / 5, 3, 10 / > (ExpN(I,4), I=1,10 ) > / 0.990,1.009,1.022,0.978,1.050,1.027,0.955 > ,0.848,0.985,1.053 / DATA > Isch(5), Iord(5), Nqrk(5), Alm(5) / 0, 1, 6, .190 / > (Vlm(I,5), I=4,6) / .190, .143, .072 / > (Qms(I,5), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(5), Qmn(5), Qmx(5) / 1.E-5, 1.60, 1.E3 / > Mex(5), Mpn(5), Nexp(5) / 5, 3, 10 / > (ExpN(I,5), I=1,10 ) > / 0.990,1.018,1.022,0.982,0.725,0.931,0.855 > ,0.795,0.965,1.007 / DATA > Isch(6), Iord(6), Nqrk(6), Alm(6) / 2, 2, 6, .235 / > (Vlm(I,6), I=4,6) / .235, .155, .060 / > (Qms(I,6), I=4,6) / 1.60, 5.00, 180.0 / > Xmn(6), Qmn(6), Qmx(6) / 1.E-5, 1.60, 1.E3 / > Mex(6), Mpn(6), Nexp(6) / 5, 3, 10 / > (ExpN(I,6), I=1,10 ) > / 0.989,1.019,1.021,0.978,0.961,0.958,0.904 > ,0.857,0.965,1.022 / Data Ist, Lp, Qsto / 0, -10, 1.2345 / save Ist, Lp, Qsto save SB, SB2, SB3 X = XX Irt = 0 if(Iset.eq.Ist .and. Qsto.eq.QQ) then if (Iprtn.eq.Lp) goto 100 if (Iprtn.ge.-3 .and. Lp.ge.-3) goto 501 endif Ip = abs(Iprtn) If (Ip .GE. 4) then If (QQ .LE. Qms(Ip, Iset)) Then Ctq2df = 0.0 Return Endif Qi = Qms(ip, Iset) Else Qi = Qmn(Iset) Endif Alam = Alm (Iset) SBL = LOG(QQ/Alam) / LOG(Qi/Alam) SB = LOG (SBL) SB2 = SB*SB SB3 = SB2*SB 501 Iflv = 3 - Iprtn Goto (1,2,3,4,5,6, 311) Iset 1 Goto(11,12,13,14,15,16,17,18,19)Iflv 11 A0=Exp( 0.2143E+00+0.8417E+00*SB -0.2451E+01*SB2+0.9875E+00*SB3) A1= 0.5209E+00-0.2384E+00*SB +0.5086E+00*SB2-0.2123E+00*SB3 A2= 0.3178E+01+0.5258E+01*SB -0.8102E+01*SB2+0.3334E+01*SB3 A3=-0.8537E+00+0.5921E+01*SB -0.1007E+02*SB2+0.4146E+01*SB3 A4= 0.1821E+01+0.2822E-01*SB +0.1662E+00*SB2-0.1058E+00*SB3 A5= 0.0000E+00-0.1090E+01*SB +0.3136E+01*SB2-0.1301E+01*SB3 goto 100 12 A0=Exp(-0.1314E+01-0.1342E-01*SB +0.1136E+00*SB2-0.1557E+00*SB3) A1= 0.2780E+00+0.2558E-01*SB +0.4467E-02*SB2-0.2472E-02*SB3 A2= 0.3672E+01+0.5324E+00*SB +0.3531E-01*SB2+0.7928E-03*SB3 A3= 0.2957E+02-0.2000E+02*SB +0.5929E+01*SB2+0.3390E+00*SB3 A4= 0.8069E+00-0.2877E+00*SB +0.3574E-01*SB2+0.5622E-02*SB3 A5= 0.0000E+00+0.2287E+00*SB -0.4052E-01*SB2+0.5589E-01*SB3 goto 100 13 A0=Exp(-0.1059E+00-0.1461E+01*SB -0.2544E+00*SB2+0.4526E-01*SB3) A1=-0.2578E+00+0.1385E+00*SB -0.1383E+00*SB2+0.3811E-01*SB3 A2= 0.5195E+01+0.9648E+00*SB -0.2103E+00*SB2-0.6701E-01*SB3 A3= 0.5131E+01+0.2151E+01*SB -0.2880E+01*SB2+0.6608E+00*SB3 A4= 0.1118E+01+0.2636E+00*SB -0.5140E+00*SB2+0.1613E+00*SB3 A5= 0.0000E+00+0.2456E+01*SB -0.8741E+00*SB2+0.2136E+00*SB3 goto 100 14 A0=Exp(-0.2732E+00-0.3523E+01*SB +0.3657E+01*SB2-0.1415E+01*SB3) A1=-0.3807E+00+0.1211E+00*SB -0.1231E+00*SB2+0.3753E-01*SB3 A2= 0.9698E+01-0.2596E+01*SB +0.2412E+01*SB2-0.9257E+00*SB3 A3=-0.6165E+00+0.1120E+01*SB -0.1708E+01*SB2+0.6383E+00*SB3 A4= 0.7292E-01-0.1339E+00*SB +0.2104E+00*SB2-0.7987E-01*SB3 A5=-0.1370E+01+0.2452E+01*SB -0.1804E+01*SB2+0.6459E+00*SB3 goto 100 15 A0=Exp(-0.2319E+01-0.3182E+01*SB +0.3572E+01*SB2-0.1431E+01*SB3) A1=-0.2622E+00+0.3085E+00*SB -0.4394E+00*SB2+0.1496E+00*SB3 A2= 0.9481E+01-0.3627E+01*SB +0.5640E+01*SB2-0.2265E+01*SB3 A3= 0.5000E+02-0.1851E+02*SB +0.2640E+01*SB2-0.6001E+00*SB3 A4= 0.1566E+01-0.7375E+00*SB +0.8736E+00*SB2-0.3449E+00*SB3 A5=-0.7983E-01+0.3236E+01*SB -0.3373E+01*SB2+0.1236E+01*SB3 goto 100 16 A0=Exp(-0.1855E+01-0.5302E+01*SB +0.8433E+00*SB2-0.1236E+00*SB3) A1=-0.4000E-02-0.1345E+01*SB +0.1192E+01*SB2-0.3039E+00*SB3 A2= 0.6870E+01+0.1246E+01*SB -0.8968E+00*SB2-0.9791E-01*SB3 A3= 0.0000E+00+0.4616E+01*SB +0.1026E+02*SB2+0.2844E+02*SB3 A4= 0.1000E-02+0.4098E+00*SB -0.4250E+00*SB2+0.1100E+00*SB3 A5= 0.0000E+00-0.2151E+01*SB +0.2991E+01*SB2-0.7717E+00*SB3 goto 100 17 A0=SB** 0.7722E+00*Exp(-0.7241E+01-0.7885E-01*SB -0.1124E+01*SB2) A1=-0.3971E+00+0.9132E+00*SB -0.1175E+01*SB2+0.3573E+00*SB3 A2= 0.6367E+01-0.6565E+01*SB +0.8114E+01*SB2-0.2666E+01*SB3 A3= 0.2878E+02-0.2000E+02*SB +0.7000E+00*SB2+0.3000E+02*SB3 A4= 0.1010E+00-0.4592E+00*SB +0.5877E+00*SB2-0.1472E+00*SB3 A5= 0.1749E+00+0.3875E+01*SB -0.3768E+01*SB2+0.1316E+01*SB3 goto 100 18 A0=SB** 0.1299E+00*Exp(-0.4868E+01-0.4339E+01*SB +0.7080E+00*SB2) A1=-0.1705E+00-0.3381E+00*SB +0.5287E+00*SB2-0.2644E+00*SB3 A2= 0.5610E+01-0.1365E+01*SB +0.1835E+01*SB2-0.5655E+00*SB3 A3=-0.1001E+01+0.3044E+01*SB +0.2680E+01*SB2+0.1426E+02*SB3 A4= 0.3814E-02+0.3430E+00*SB -0.6926E+00*SB2+0.3486E+00*SB3 A5= 0.1156E+01+0.2016E+01*SB -0.1674E+01*SB2+0.5981E+00*SB3 goto 100 19 A0=SB** 0.9819E+00*Exp(-0.7859E+01+0.6819E+00*SB -0.3386E+01*SB2) A1=-0.1055E+00-0.1413E+01*SB +0.3451E+01*SB2-0.2466E+01*SB3 A2= 0.4055E+01+0.8107E+01*SB -0.1576E+02*SB2+0.8094E+01*SB3 A3= 0.3799E+01+0.9616E+01*SB -0.1984E+02*SB2+0.2641E+02*SB3 A4= 0.3619E+00-0.8627E+00*SB -0.9390E-01*SB2+0.9196E+00*SB3 A5= 0.3779E+01-0.6073E+01*SB +0.9999E+01*SB2-0.4304E+01*SB3 goto 100 2 Goto(21,22,23,24,25,26,27,28,29)Iflv 21 A0=Exp( 0.2790E+00+0.7294E+00*SB -0.2202E+01*SB2+0.8599E+00*SB3) A1= 0.5380E+00-0.2261E+00*SB +0.4636E+00*SB2-0.1871E+00*SB3 A2= 0.3259E+01+0.2141E+01*SB -0.2947E+01*SB2+0.1245E+01*SB3 A3=-0.8390E+00+0.1448E+01*SB -0.2331E+01*SB2+0.8658E+00*SB3 A4= 0.1847E+01-0.3943E+01*SB +0.5998E+01*SB2-0.2191E+01*SB3 A5= 0.0000E+00-0.9719E+00*SB +0.2830E+01*SB2-0.1137E+01*SB3 goto 100 22 A0=Exp(-0.1318E+01+0.2328E-01*SB +0.5179E-01*SB2-0.1305E+00*SB3) A1= 0.2760E+00+0.4429E-01*SB -0.2626E-01*SB2+0.7143E-02*SB3 A2= 0.3660E+01+0.5232E+00*SB +0.5491E-01*SB2-0.4115E-02*SB3 A3= 0.2910E+02-0.2000E+02*SB +0.6631E+01*SB2-0.3050E-01*SB3 A4= 0.8010E+00-0.2688E+00*SB +0.1051E-01*SB2+0.1195E-01*SB3 A5= 0.0000E+00+0.2887E+00*SB -0.1398E+00*SB2+0.8194E-01*SB3 goto 100 23 A0=Exp(-0.1623E+01-0.7232E+00*SB +0.1889E+00*SB2+0.1140E+00*SB3) A1=-0.5000E+00+0.8611E-01*SB +0.2203E-01*SB2-0.1401E-01*SB3 A2= 0.3821E+01+0.8976E+00*SB +0.1400E+00*SB2-0.9163E-01*SB3 A3= 0.5809E+01-0.5060E+01*SB +0.3808E+00*SB2+0.2519E+00*SB3 A4= 0.4500E+00-0.5121E+00*SB +0.1979E+00*SB2-0.2705E-01*SB3 A5= 0.0000E+00+0.1210E+01*SB -0.2921E+00*SB2+0.1240E+00*SB3 goto 100 24 A0=Exp(-0.6986E-01-0.5954E+00*SB -0.1582E+01*SB2+0.5104E+00*SB3) A1=-0.8461E+00+0.2127E+00*SB +0.9425E-01*SB2-0.5264E-01*SB3 A2= 0.1200E+02+0.1659E+01*SB -0.5354E+01*SB2+0.1795E+01*SB3 A3= 0.2958E+02+0.3000E+02*SB +0.3000E+02*SB2-0.1965E+02*SB3 A4= 0.4000E+01-0.4865E+00*SB +0.9460E+00*SB2+0.3432E+00*SB3 A5=-0.3378E+01+0.1656E+01*SB +0.1123E+01*SB2-0.4667E+00*SB3 goto 100 25 A0=Exp(-0.1929E+01-0.2626E+01*SB +0.2926E+01*SB2-0.1297E+01*SB3) A1=-0.6627E+00+0.4561E+00*SB -0.3818E+00*SB2+0.1239E+00*SB3 A2= 0.9506E+01-0.2724E+01*SB +0.4283E+01*SB2-0.1804E+01*SB3 A3= 0.1897E+02+0.1642E+01*SB -0.8390E+01*SB2+0.3894E+01*SB3 A4= 0.1024E+01-0.1786E+00*SB +0.4535E+00*SB2-0.2075E+00*SB3 A5=-0.1746E+01+0.3572E+01*SB -0.2908E+01*SB2+0.1093E+01*SB3 goto 100 26 A0=Exp(-0.4913E+00-0.6866E+01*SB +0.1432E+01*SB2-0.1749E+00*SB3) A1=-0.1157E+00-0.1567E+01*SB +0.1439E+01*SB2-0.3724E+00*SB3 A2= 0.7730E+01+0.9748E+00*SB -0.1157E+01*SB2-0.8358E-02*SB3 A3=-0.6050E+00+0.1835E+01*SB +0.3788E+01*SB2+0.3000E+02*SB3 A4= 0.1620E-08+0.4590E+00*SB -0.4070E+00*SB2+0.8900E-01*SB3 A5=-0.7048E+00-0.2505E+01*SB +0.4000E+01*SB2-0.1161E+01*SB3 goto 100 27 A0=SB** 0.7393E+00*Exp(-0.6518E+01-0.3998E+00*SB -0.1111E+01*SB2) A1=-0.6482E+00+0.1125E+01*SB -0.1290E+01*SB2+0.3940E+00*SB3 A2= 0.8487E+01-0.9235E+01*SB +0.9353E+01*SB2-0.2913E+01*SB3 A3= 0.2265E+02-0.1999E+02*SB +0.4105E+01*SB2+0.2144E+02*SB3 A4= 0.8990E-01-0.4372E+00*SB +0.5941E+00*SB2-0.1469E+00*SB3 A5=-0.9690E+00+0.5068E+01*SB -0.4368E+01*SB2+0.1503E+01*SB3 goto 100 28 A0=SB** 0.9880E+00*Exp(-0.7180E+01-0.2494E+01*SB +0.3561E-01*SB2) A1=-0.4301E+00-0.2611E+00*SB +0.3914E+00*SB2-0.1638E+00*SB3 A2= 0.5137E+01+0.1506E+01*SB -0.9588E+00*SB2-0.1596E+00*SB3 A3= 0.1483E+02+0.2998E+02*SB +0.2357E+02*SB2-0.9353E+01*SB3 A4= 0.2426E+00+0.1371E+00*SB -0.3791E+00*SB2+0.1948E+00*SB3 A5= 0.1463E+01+0.1907E+00*SB +0.3557E+00*SB2+0.2097E-01*SB3 goto 100 29 A0=SB** 0.1005E+01*Exp(-0.5255E+01-0.9866E-01*SB -0.2737E+01*SB2) A1=-0.3140E+00-0.2055E+00*SB +0.5594E+00*SB2-0.2960E+00*SB3 A2= 0.9227E+01-0.4569E+01*SB -0.9724E+01*SB2+0.1026E+02*SB3 A3= 0.1131E+02-0.1972E+02*SB -0.1107E+02*SB2+0.2311E+02*SB3 A4= 0.1488E+01+0.1737E+01*SB +0.4323E+01*SB2-0.9925E+01*SB3 A5= 0.1895E+01-0.7350E+00*SB +0.3780E+01*SB2-0.1408E+01*SB3 goto 100 3 Goto(31,32,33,34,35,36,37,38,39)Iflv 31 A0=Exp(-0.7913E+00-0.2789E+01*SB -0.7289E-01*SB2+0.1770E+00*SB3) A1= 0.4942E+00-0.7886E-01*SB +0.9057E-01*SB2-0.5259E-01*SB3 A2= 0.3727E+01+0.1089E+01*SB -0.1004E+01*SB2+0.4345E+00*SB3 A3= 0.1944E+01+0.7846E+01*SB +0.7984E+01*SB2+0.5548E+01*SB3 A4= 0.2940E-02+0.8428E-04*SB +0.1266E+00*SB2-0.3517E-01*SB3 A5=-0.1060E+00-0.1192E-01*SB +0.1130E+01*SB2-0.4527E+00*SB3 goto 100 32 A0=Exp(-0.1344E+01+0.7859E-02*SB +0.4623E-01*SB2-0.1273E+00*SB3) A1= 0.2760E+00+0.4201E-01*SB -0.1795E-01*SB2+0.3212E-02*SB3 A2= 0.3660E+01+0.5247E+00*SB +0.4405E-01*SB2+0.1391E-02*SB3 A3= 0.2981E+02-0.2000E+02*SB +0.6566E+01*SB2+0.2479E-01*SB3 A4= 0.7950E+00-0.2732E+00*SB +0.2470E-01*SB2+0.6157E-02*SB3 A5= 0.0000E+00+0.2793E+00*SB -0.9197E-01*SB2+0.5953E-01*SB3 goto 100 33 A0=Exp( 0.9746E+00-0.3252E+01*SB +0.1664E+01*SB2-0.6410E+00*SB3) A1=-0.5271E-02-0.3198E+00*SB +0.1279E+00*SB2-0.1256E-02*SB3 A2= 0.5740E+01-0.3139E+01*SB +0.3841E+01*SB2-0.1415E+01*SB3 A3= 0.7161E-01-0.4363E+01*SB +0.4925E+01*SB2-0.1614E+01*SB3 A4= 0.1860E+01+0.1342E+01*SB -0.2234E+01*SB2+0.1047E+01*SB3 A5= 0.7409E-01+0.2390E+01*SB -0.1457E+01*SB2+0.5853E+00*SB3 goto 100 34 A0=Exp(-0.8454E+00-0.3334E+01*SB +0.3591E+01*SB2-0.1485E+01*SB3) A1=-0.2826E-02-0.2810E+00*SB -0.3809E-01*SB2+0.6585E-01*SB3 A2= 0.9139E+01-0.2811E+01*SB +0.4730E+01*SB2-0.2157E+01*SB3 A3=-0.3120E+00+0.1217E+01*SB -0.1726E+01*SB2+0.6220E+00*SB3 A4= 0.1793E-01-0.4608E-01*SB +0.5294E-01*SB2-0.1709E-01*SB3 A5=-0.1471E+00+0.1104E+01*SB -0.1358E+01*SB2+0.7200E+00*SB3 goto 100 35 A0=Exp(-0.1398E+01-0.3536E+01*SB +0.3849E+01*SB2-0.1549E+01*SB3) A1=-0.1332E-01-0.2155E-01*SB -0.3404E+00*SB2+0.1569E+00*SB3 A2= 0.9981E+01-0.3499E+01*SB +0.5448E+01*SB2-0.2198E+01*SB3 A3= 0.3736E+02-0.2000E+02*SB +0.6675E+01*SB2-0.7276E+00*SB3 A4= 0.1705E+01-0.1013E+01*SB +0.1122E+01*SB2-0.4057E+00*SB3 A5=-0.1189E-01+0.2698E+01*SB -0.3429E+01*SB2+0.1389E+01*SB3 goto 100 36 A0=Exp(-0.2979E+01-0.6085E+01*SB +0.2428E+01*SB2-0.6482E+00*SB3) A1=-0.1372E+00-0.1281E+00*SB +0.1587E+00*SB2-0.9637E-01*SB3 A2= 0.7009E+01-0.1609E+01*SB +0.2765E+01*SB2-0.1177E+01*SB3 A3= 0.1308E+01+0.9583E+01*SB +0.2360E+02*SB2+0.2999E+02*SB3 A4= 0.2509E-01+0.2106E+00*SB -0.4405E+00*SB2+0.2075E+00*SB3 A5=-0.2069E-01+0.1971E+01*SB -0.1615E+01*SB2+0.6039E+00*SB3 goto 100 37 A0=SB** 0.8072E+00*Exp(-0.6920E+01-0.5031E+00*SB -0.9965E+00*SB2) A1=-0.2118E+00+0.7930E+00*SB -0.1101E+01*SB2+0.3302E+00*SB3 A2= 0.8039E+01-0.7170E+01*SB +0.8657E+01*SB2-0.2893E+01*SB3 A3= 0.2926E+02-0.1993E+02*SB +0.1841E+01*SB2+0.2996E+02*SB3 A4= 0.1339E+00-0.5531E+00*SB +0.6505E+00*SB2-0.1595E+00*SB3 A5= 0.7439E+00+0.3307E+01*SB -0.3284E+01*SB2+0.1152E+01*SB3 goto 100 38 A0=SB** 0.9925E+00*Exp(-0.2190E+01-0.3393E+01*SB -0.8631E+00*SB2) A1=-0.1261E+00-0.2368E+00*SB +0.4143E+00*SB2-0.1577E+00*SB3 A2= 0.4585E+01+0.5227E+01*SB -0.3248E+01*SB2-0.2599E+00*SB3 A3=-0.1094E+01+0.4927E+00*SB -0.9921E+00*SB2+0.3138E+01*SB3 A4= 0.1396E+00+0.2562E+00*SB +0.1844E+00*SB2-0.1599E+00*SB3 A5= 0.8621E+00+0.4715E+00*SB +0.2547E+01*SB2-0.8429E+00*SB3 goto 100 39 A0=SB** 0.1016E+01*Exp(-0.5397E+01-0.1979E+01*SB -0.2441E+00*SB2) A1=-0.1426E+00-0.2861E+00*SB +0.7434E+00*SB2-0.5214E+00*SB3 A2= 0.6363E+01+0.4028E+00*SB -0.8356E+01*SB2+0.6814E+01*SB3 A3=-0.2526E+00+0.2425E+01*SB -0.1407E+02*SB2+0.3000E+02*SB3 A4= 0.1125E+00-0.1089E+01*SB +0.9977E+01*SB2+0.1000E+02*SB3 A5= 0.2669E+01-0.6366E+00*SB +0.4355E+01*SB2-0.2919E+01*SB3 goto 100 4 Goto(41,42,43,44,45,46,47,48,49)Iflv 41 A0=Exp( 0.3760E+00+0.5491E+00*SB -0.1845E+01*SB2+0.6803E+00*SB3) A1= 0.5650E+00-0.1953E+00*SB +0.3761E+00*SB2-0.1419E+00*SB3 A2= 0.3464E+01+0.3817E+01*SB -0.5384E+01*SB2+0.2057E+01*SB3 A3=-0.5850E+00+0.5566E+01*SB -0.9000E+01*SB2+0.3433E+01*SB3 A4= 0.2322E+01-0.1431E+00*SB +0.3901E+00*SB2-0.1678E+00*SB3 A5= 0.0000E+00-0.7370E+00*SB +0.2310E+01*SB2-0.8743E+00*SB3 goto 100 42 A0=Exp(-0.1324E+01+0.1169E-01*SB +0.1969E-01*SB2-0.7583E-01*SB3) A1= 0.2890E+00+0.5832E-01*SB -0.2921E-01*SB2+0.4701E-02*SB3 A2= 0.3580E+01+0.5291E+00*SB -0.5662E-02*SB2+0.2746E-01*SB3 A3= 0.3021E+02-0.1999E+02*SB +0.6250E+01*SB2-0.3035E+00*SB3 A4= 0.7990E+00-0.2531E+00*SB +0.5556E-02*SB2+0.8272E-02*SB3 A5= 0.0000E+00+0.3674E+00*SB -0.1383E+00*SB2+0.4665E-01*SB3 goto 100 43 A0=Exp(-0.1920E+00-0.7015E+00*SB -0.9113E+00*SB2+0.2352E+00*SB3) A1=-0.2120E+00+0.1133E-01*SB -0.1553E-01*SB2+0.2822E-02*SB3 A2= 0.4549E+01+0.1250E+01*SB -0.4647E+00*SB2+0.9617E-01*SB3 A3= 0.1197E+02-0.4156E+01*SB +0.1413E+00*SB2+0.1607E+00*SB3 A4= 0.1616E+01+0.1082E+00*SB -0.6651E+00*SB2+0.2356E+00*SB3 A5= 0.0000E+00+0.1824E+01*SB -0.2063E+00*SB2+0.1148E-01*SB3 goto 100 44 A0=Exp(-0.1388E+01-0.7408E+00*SB -0.6454E+00*SB2+0.2373E+00*SB3) A1=-0.2928E+00-0.1726E-01*SB +0.4033E-01*SB2-0.2514E-01*SB3 A2= 0.9975E+01-0.2048E+01*SB -0.6060E+00*SB2+0.5225E+00*SB3 A3= 0.2687E+02-0.4683E+01*SB -0.1999E+02*SB2+0.1188E+02*SB3 A4= 0.4000E+01-0.6773E+00*SB +0.4301E+00*SB2+0.4524E+00*SB3 A5=-0.7164E+00+0.7488E+00*SB +0.5766E+00*SB2-0.2609E+00*SB3 goto 100 45 A0=Exp(-0.2272E+01-0.2998E+01*SB +0.3282E+01*SB2-0.1203E+01*SB3) A1=-0.2062E+00+0.3320E+00*SB -0.5074E+00*SB2+0.1655E+00*SB3 A2= 0.9667E+01-0.3497E+01*SB +0.5271E+01*SB2-0.1984E+01*SB3 A3= 0.4996E+02-0.3241E+01*SB -0.1425E+02*SB2+0.3849E+01*SB3 A4= 0.1619E+01-0.5354E+00*SB +0.5753E+00*SB2-0.2238E+00*SB3 A5= 0.8755E-01+0.3195E+01*SB -0.3496E+01*SB2+0.1197E+01*SB3 goto 100 46 A0=Exp(-0.1864E+01-0.5258E+01*SB +0.1034E+01*SB2-0.1550E+00*SB3) A1= 0.1000E-02-0.1090E+01*SB +0.8345E+00*SB2-0.1887E+00*SB3 A2= 0.6898E+01-0.4951E+00*SB +0.4279E+00*SB2-0.2727E+00*SB3 A3= 0.0000E+00+0.4322E+01*SB +0.8181E+01*SB2+0.2309E+02*SB3 A4= 0.1000E-02+0.3550E+00*SB -0.3220E+00*SB2+0.7294E-01*SB3 A5= 0.0000E+00-0.1347E+01*SB +0.1896E+01*SB2-0.4491E+00*SB3 goto 100 47 A0=SB** 0.7528E+00*Exp(-0.7684E+01+0.6791E-01*SB -0.9094E+00*SB2) A1=-0.3732E+00+0.8408E+00*SB -0.1020E+01*SB2+0.3046E+00*SB3 A2= 0.4984E+01-0.5534E+01*SB +0.6418E+01*SB2-0.1856E+01*SB3 A3= 0.3761E+02-0.1999E+02*SB -0.3358E+01*SB2+0.2999E+02*SB3 A4= 0.1161E+00-0.4680E+00*SB +0.5567E+00*SB2-0.1633E+00*SB3 A5= 0.3028E+00+0.3339E+01*SB -0.3004E+01*SB2+0.9160E+00*SB3 goto 100 48 A0=SB** 0.1011E+01*Exp(-0.7217E+01-0.2288E+01*SB +0.3450E+00*SB2) A1=-0.1955E+00-0.3371E+00*SB +0.5111E+00*SB2-0.2210E+00*SB3 A2= 0.4302E+01-0.1214E+01*SB +0.3104E+01*SB2-0.1408E+01*SB3 A3= 0.1487E+02+0.1549E+02*SB +0.2875E+02*SB2-0.1922E+02*SB3 A4= 0.8935E-02+0.3571E+00*SB -0.6668E+00*SB2+0.3037E+00*SB3 A5= 0.1570E+01+0.7105E+00*SB -0.6070E+00*SB2+0.3796E+00*SB3 goto 100 49 A0=SB** 0.9986E+00*Exp(-0.5847E+01-0.2798E+00*SB -0.9882E+00*SB2) A1=-0.2154E+00-0.8282E-01*SB +0.3611E-01*SB2+0.2623E-01*SB3 A2= 0.3250E+01+0.9635E+01*SB -0.1274E+02*SB2+0.4453E+01*SB3 A3=-0.2594E+01+0.9097E+01*SB +0.1581E+02*SB2-0.9123E+01*SB3 A4= 0.1768E+01-0.2749E+01*SB +0.9999E+01*SB2+0.9995E+01*SB3 A5= 0.2521E+01-0.1802E-01*SB +0.4820E+00*SB2+0.2004E+00*SB3 goto 100 5 Goto(51,52,53,54,55,56,57,58,59)Iflv 51 A0=Exp( 0.7248E-01+0.3941E+00*SB -0.1772E+01*SB2+0.7629E+00*SB3) A1= 0.4964E+00-0.1224E+00*SB +0.3646E+00*SB2-0.1685E+00*SB3 A2= 0.3000E+01+0.2780E+01*SB -0.4028E+01*SB2+0.1816E+01*SB3 A3=-0.1064E+01+0.3062E+01*SB -0.5927E+01*SB2+0.2785E+01*SB3 A4= 0.3193E+01+0.1499E+01*SB -0.2765E+01*SB2+0.1019E+01*SB3 A5= 0.1524E-01-0.4541E+00*SB +0.2281E+01*SB2-0.1033E+01*SB3 goto 100 52 A0=Exp(-0.1794E+01-0.2055E+00*SB -0.3350E-01*SB2-0.5084E-01*SB3) A1= 0.1748E+00+0.4637E-01*SB -0.2048E-01*SB2+0.2596E-02*SB3 A2= 0.3321E+01+0.6253E+00*SB +0.2148E-01*SB2+0.1288E-01*SB3 A3= 0.4355E+02-0.2000E+02*SB +0.5486E+01*SB2+0.1536E+00*SB3 A4= 0.9586E+00-0.3217E+00*SB +0.4458E-01*SB2-0.1404E-03*SB3 A5=-0.6595E-02+0.3499E+00*SB -0.7048E-01*SB2+0.2619E-01*SB3 goto 100 53 A0=Exp(-0.6194E+00-0.2643E+00*SB -0.1875E+01*SB2+0.6011E+00*SB3) A1=-0.2600E+00+0.8704E-01*SB -0.7375E-01*SB2+0.1876E-01*SB3 A2= 0.4620E+01+0.1578E+01*SB -0.8411E+00*SB2+0.1527E+00*SB3 A3= 0.1604E+02-0.1230E+02*SB +0.6939E+01*SB2-0.2012E+01*SB3 A4= 0.1255E+01+0.4769E+00*SB -0.9915E+00*SB2+0.3439E+00*SB3 A5= 0.1116E-02+0.2409E+01*SB -0.4442E+00*SB2+0.3431E-01*SB3 goto 100 54 A0=Exp(-0.1571E+01-0.1905E+00*SB -0.8672E+00*SB2+0.2070E+00*SB3) A1=-0.3266E+00+0.6428E-01*SB -0.8694E-01*SB2+0.1778E-01*SB3 A2= 0.8921E+01-0.5010E+00*SB -0.9658E+00*SB2+0.3893E+00*SB3 A3= 0.1329E+02+0.4652E+01*SB -0.2000E+02*SB2+0.1001E+02*SB3 A4= 0.3283E+01-0.3400E+00*SB -0.1957E+00*SB2+0.8063E+00*SB3 A5=-0.5701E+00+0.4042E+00*SB +0.5239E+00*SB2-0.1665E+00*SB3 goto 100 55 A0=Exp(-0.2281E+01-0.2768E+01*SB +0.3137E+01*SB2-0.1278E+01*SB3) A1=-0.2624E+00+0.4142E+00*SB -0.5936E+00*SB2+0.1937E+00*SB3 A2= 0.9438E+01-0.3179E+01*SB +0.5107E+01*SB2-0.2179E+01*SB3 A3= 0.5000E+02-0.1802E+02*SB -0.7515E+01*SB2+0.2991E+01*SB3 A4= 0.1809E+01-0.9121E+00*SB +0.8854E+00*SB2-0.3582E+00*SB3 A5= 0.4056E-01+0.3033E+01*SB -0.3431E+01*SB2+0.1253E+01*SB3 goto 100 56 A0=Exp(-0.2318E+01-0.4104E+01*SB -0.1502E+00*SB2+0.1693E+00*SB3) A1=-0.2251E-01-0.1101E+01*SB +0.1037E+01*SB2-0.3290E+00*SB3 A2= 0.6989E+01+0.1794E+01*SB -0.1811E+01*SB2+0.3061E+00*SB3 A3= 0.7972E+00+0.7806E+01*SB +0.1869E+02*SB2+0.2999E+02*SB3 A4= 0.4795E-01+0.1622E+00*SB -0.3977E+00*SB2+0.1920E+00*SB3 A5=-0.5275E-01-0.2616E+01*SB +0.3076E+01*SB2-0.7425E+00*SB3 goto 100 57 A0=SB** 0.8431E+00*Exp(-0.6539E+01-0.1875E+00*SB -0.1346E+01*SB2) A1=-0.4970E+00+0.9062E+00*SB -0.1169E+01*SB2+0.3703E+00*SB3 A2= 0.4939E+01-0.2995E+01*SB +0.4483E+01*SB2-0.1704E+01*SB3 A3= 0.3113E+02-0.1997E+02*SB +0.1540E+01*SB2+0.3000E+02*SB3 A4= 0.1349E+00-0.5418E+00*SB +0.6142E+00*SB2-0.1360E+00*SB3 A5=-0.8590E+00+0.3956E+01*SB -0.3612E+01*SB2+0.1401E+01*SB3 goto 100 58 A0=SB** 0.2639E-01*Exp(-0.2099E+01-0.2681E+01*SB +0.2925E+00*SB2) A1=-0.2243E+00-0.5343E-01*SB -0.1953E-01*SB2+0.1586E-01*SB3 A2= 0.4294E+01+0.1102E+01*SB -0.1822E+00*SB2-0.2481E+00*SB3 A3=-0.9998E+00+0.8275E-01*SB +0.5494E+00*SB2-0.1982E+00*SB3 A4= 0.5904E-04+0.9222E-01*SB -0.9293E-01*SB2+0.9159E-01*SB3 A5= 0.2657E+00+0.1770E+01*SB -0.7111E+00*SB2+0.2525E+00*SB3 goto 100 59 A0=SB** 0.1009E+01*Exp(-0.7032E+01+0.4562E+01*SB -0.9081E+01*SB2) A1=-0.1412E+00-0.5076E+00*SB +0.9513E+00*SB2-0.4326E+00*SB3 A2= 0.5385E+01+0.3023E+01*SB -0.1162E+02*SB2+0.7006E+01*SB3 A3= 0.4997E+01-0.1600E+02*SB +0.1342E+02*SB2+0.1197E+02*SB3 A4= 0.5825E+00+0.3994E+00*SB -0.1255E+01*SB2+0.6486E+00*SB3 A5= 0.3365E+01-0.4026E+01*SB +0.8385E+01*SB2-0.2260E+01*SB3 goto 100 6 Goto(61,62,63,64,65,66,67,68,69)Iflv 61 A0=Exp( 0.1590E+00+0.5580E+00*SB -0.1838E+01*SB2+0.7018E+00*SB3) A1= 0.5110E+00-0.1625E+00*SB +0.3547E+00*SB2-0.1412E+00*SB3 A2= 0.3158E+01+0.3962E+01*SB -0.5866E+01*SB2+0.2375E+01*SB3 A3=-0.6000E+00+0.6144E+01*SB -0.1056E+02*SB2+0.4345E+01*SB3 A4= 0.2306E+01-0.4669E-01*SB +0.2711E+00*SB2-0.1640E+00*SB3 A5= 0.0000E+00-0.6638E+00*SB +0.2239E+01*SB2-0.8843E+00*SB3 goto 100 62 A0=Exp(-0.1182E+01+0.1449E+00*SB +0.2753E-01*SB2-0.1009E+00*SB3) A1= 0.2540E+00+0.2686E-01*SB -0.1546E-01*SB2+0.5396E-02*SB3 A2= 0.3442E+01+0.5576E+00*SB +0.1937E-01*SB2+0.6696E-02*SB3 A3= 0.2545E+02-0.2000E+02*SB +0.7355E+01*SB2-0.7058E+00*SB3 A4= 0.9170E+00-0.3090E+00*SB +0.1705E-01*SB2+0.8534E-02*SB3 A5= 0.0000E+00+0.1449E+00*SB -0.7821E-01*SB2+0.6405E-01*SB3 goto 100 63 A0=Exp(-0.3410E+00-0.9613E+00*SB -0.4969E+00*SB2+0.9360E-01*SB3) A1=-0.2400E+00+0.1473E+00*SB -0.1593E+00*SB2+0.4538E-01*SB3 A2= 0.4841E+01+0.9311E+00*SB +0.1601E-03*SB2-0.1331E+00*SB3 A3= 0.7427E+01-0.1397E+01*SB +0.1489E+00*SB2-0.2848E+00*SB3 A4= 0.9600E+00+0.3697E+00*SB -0.4246E+00*SB2+0.1032E+00*SB3 A5= 0.0000E+00+0.2484E+01*SB -0.9908E+00*SB2+0.2568E+00*SB3 goto 100 64 A0=Exp( 0.1176E+00-0.3418E+01*SB +0.3529E+01*SB2-0.1367E+01*SB3) A1=-0.3654E+00+0.1914E+00*SB -0.2192E+00*SB2+0.6933E-01*SB3 A2= 0.1099E+02-0.4281E+01*SB +0.3729E+01*SB2-0.1254E+01*SB3 A3=-0.7514E+00+0.7696E+00*SB -0.1134E+01*SB2+0.4245E+00*SB3 A4= 0.7690E-01-0.6558E-01*SB +0.8726E-01*SB2-0.3345E-01*SB3 A5=-0.1447E+01+0.2617E+01*SB -0.2094E+01*SB2+0.7536E+00*SB3 goto 100 65 A0=Exp(-0.2412E+01-0.2522E+01*SB +0.3126E+01*SB2-0.1305E+01*SB3) A1=-0.2353E+00+0.3118E+00*SB -0.4864E+00*SB2+0.1689E+00*SB3 A2= 0.9017E+01-0.2437E+01*SB +0.4659E+01*SB2-0.2044E+01*SB3 A3= 0.5000E+02-0.1158E+02*SB -0.9260E+01*SB2+0.2847E+01*SB3 A4= 0.1726E+01-0.6849E+00*SB +0.7864E+00*SB2-0.3300E+00*SB3 A5= 0.5080E-01+0.2858E+01*SB -0.3297E+01*SB2+0.1246E+01*SB3 goto 100 66 A0=Exp(-0.1966E+01-0.4405E+01*SB +0.2436E+00*SB2+0.4576E-01*SB3) A1=-0.4000E-02-0.1229E+01*SB +0.1118E+01*SB2-0.2988E+00*SB3 A2= 0.6902E+01+0.1266E+01*SB -0.1068E+01*SB2+0.3062E-01*SB3 A3= 0.0000E+00+0.3987E+01*SB +0.9389E+01*SB2+0.1881E+02*SB3 A4= 0.1000E-02+0.3528E+00*SB -0.4201E+00*SB2+0.1248E+00*SB3 A5= 0.0000E+00-0.2149E+01*SB +0.2925E+01*SB2-0.7609E+00*SB3 goto 100 67 A0=SB** 0.7561E+00*Exp(-0.6960E+01+0.5634E-01*SB -0.1170E+01*SB2) A1=-0.4232E+00+0.9269E+00*SB -0.1161E+01*SB2+0.3470E+00*SB3 A2= 0.6057E+01-0.5790E+01*SB +0.7352E+01*SB2-0.2435E+01*SB3 A3= 0.2941E+02-0.1999E+02*SB -0.8345E+00*SB2+0.3000E+02*SB3 A4= 0.1069E+00-0.4620E+00*SB +0.5614E+00*SB2-0.1336E+00*SB3 A5=-0.1865E+00+0.3953E+01*SB -0.3791E+01*SB2+0.1315E+01*SB3 goto 100 68 A0=SB** 0.5661E-02*Exp(-0.2123E+01-0.3026E+01*SB +0.1912E+00*SB2) A1=-0.2011E+00-0.1338E-01*SB -0.3974E-01*SB2+0.1948E-01*SB3 A2= 0.4906E+01+0.1740E+01*SB -0.1387E+01*SB2+0.1263E+00*SB3 A3=-0.1000E+01+0.5767E-01*SB +0.6377E+00*SB2+0.4736E-01*SB3 A4= 0.5927E-04+0.1039E+00*SB -0.9797E-01*SB2+0.6881E-01*SB3 A5= 0.4017E+00+0.1981E+01*SB -0.7758E+00*SB2+0.2916E+00*SB3 goto 100 69 A0=SB** 0.1008E+01*Exp(-0.7211E+01+0.3273E+01*SB -0.6979E+01*SB2) A1=-0.1026E+00-0.4948E+00*SB +0.1188E+01*SB2-0.8016E+00*SB3 A2= 0.5397E+01+0.2135E+01*SB -0.9531E+01*SB2+0.6115E+01*SB3 A3= 0.4966E+01-0.1111E+02*SB +0.4732E+01*SB2+0.1568E+02*SB3 A4= 0.5345E+00-0.1935E+00*SB +0.5816E+00*SB2-0.6794E+00*SB3 A5= 0.3569E+01-0.3477E+01*SB +0.8756E+01*SB2-0.4139E+01*SB3 goto 100 311 stop 'This option is not currently supported.' 100 Ctq2df = A0 *(x**A1) *((1.-x)**A2) *(1.+A3*(x**A4)) $ *(log(1.+1./x))**A5 if(Ctq2df.lt.0.0) then Ctq2df = 0.0 Irt=1 endif Ist = Iset Lp = Iprtn Qsto = QQ Return ENTRY Wlamd2 (Iset, Iorder, Neff) Iorder = IORD (Iset) Wlamd2 = VLM (Neff, Iset) RETURN Entry PrCtq2 > (Iset, Iordr, Ischeme, MxFlv, > Alam4, Alam5, Alam6, Amas4, Amas5, Amas6, > Xmin, Qini, Qmax, ExpNor) Iordr = Iord (Iset) Ischeme= Isch (Iset) MxFlv = Nqrk (Iset) Alam4 = Vlm(4,Iset) Alam5 = Vlm(5,Iset) Alam6 = Vlm(6,Iset) Amas4 = Qms(4,Iset) Amas5 = Qms(5,Iset) Amas6 = Qms(6,Iset) Xmin = Xmn (Iset) Qini = Qmn (Iset) Qmax = Qmx (Iset) Do 201 Iexp = 1, Nexp(Iset) 201 ExpNor(Iexp) = ExpN(Iexp, Iset) Return END FUNCTION CtqPdf (Iset, Iprtn, XX, QQ, Irt) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (Nex = 5, MxFl = 6, Npn = 3, Nst = 30, Nexpt=20) Parameter (Nst4 = Nst*4) DIMENSION > Iord(Nst), Isch(Nst), Nqrk(Nst),Alm(Nst) > , Vlm(4:6,Nst), Qms(4:6, Nst) > , Xmn(Nst), Qmn(Nst), Qmx(Nst), Nexp(Nexpt) > , Mex(Nst), Mpn(Nst), ExpN(Nexpt, Nst), ExpNor(Nexpt) DATA > Isch(1), Iord(1), Nqrk(1), Alm(1) / 1, 2, 6, .152 / > (Vlm(I,1), I=4,6) / .231, .152, .059 / > (Qms(I,1), I=4,6) / 1.50, 5.00, 180.0 / > Xmn(1), Qmn(1), Qmx(1) / 1.E-5, 2.00, 1.E3 / > Mex(1), Mpn(1), Nexp(1) / 5, 3, 8 / > (ExpN(I, 1), I=1,8) > / 0.989, 1.00, 1.02, 0.978, 1.10, 0.972, 0.987, 0.846 / DATA > Isch(2), Iord(2), Nqrk(2), Alm(2) / 1, 2, 6, .152 / > (Vlm(I,2), I=4,6) / .231, .152, .059 / > (Qms(I,2), I=4,6) / 1.50, 5.00, 180.0 / > Xmn(2), Qmn(2), Qmx(2) / 1.E-5, 2.00, 1.E3 / > Mex(2), Mpn(2), Nexp(2) / 5, 3, 8 / > (ExpN(I, 2), I=1,8 ) > / 0.989, 1.00, 1.02, 0.984, 1.05, 0.891, 0.923, 0.824 / DATA > Isch(3), Iord(3), Nqrk(3), Alm(3) / 1, 2, 6, .220 / > (Vlm(I,3), I=4,6) / .322, .220, .088 / > (Qms(I,3), I=4,6) / 1.50, 5.00, 180.0 / > Xmn(3), Qmn(3), Qmx(3) / 1.E-5, 2.00, 1.E3 / > Mex(3), Mpn(3), Nexp(3) / 5, 3, 8 / > (ExpN(I, 3), I=1,8 ) > / 0.985, 1.00, 1.01, 0.977, 1.07, 1.31, 1.19, 1.09 / DATA > Isch(4), Iord(4), Nqrk(4), Alm(4) / 2, 2, 6, .164 / > (Vlm(I,4), I=4,6) / .247, .164, .064 / > (Qms(I,4), I=4,6) / 1.50, 5.00, 180.0 / > Xmn(4), Qmn(4), Qmx(4) / 1.E-5, 2.00, 1.E3 / > Mex(4), Mpn(4), Nexp(4) / 5, 3, 8 / > (ExpN(I, 4), I=1,8 ) > / 0.983, 1.00, 1.01, 0.975, 0.964, 1.23, 1.00, 1.12 / DATA > Isch(5), Iord(5), Nqrk(5), Alm(5) / 1, 1, 6, .125 / > (Vlm(I,5), I=4,6) / .168, .125, .063 / > (Qms(I,5), I=4,6) / 1.50, 5.00, 180.0 / > Xmn(5), Qmn(5), Qmx(5) / 1.E-5, 2.00, 1.E3 / > Mex(5), Mpn(5), Nexp(5) / 5, 3, 8 / > (ExpN(I, 5), I=1,8 ) > / 0.982, 1.01, 1.00, 0.972, 0.840, 0.959, 0.930, 0.861 / Data ist, lp, qsto, Aln2 / 0, -10, 1.2345, 0.6931 / X = XX if(iset.eq.ist.and.iprtn.eq.lp.and.qsto.eq.qq) goto 100 Irt = 0 ip = abs(iprtn) If (Ip.GE.5.and.QQ.LE.Qms(Ip, Iset)) Then CtqPdf = 0.0 Return Endif if(ip.ge.5) then Qi = qms(ip,iset) else Qi = Qmn (Iset) endif Alam = Alm (Iset) sta = log(qq/alam) stb = log(qi/alam) sbl = sta/stb SB = LOG (SBL) SB2 = SB*SB SB3 = SB2*SB iflv = 3 - iprtn Goto (1, 2, 3, 4, 5), Iset 1 Goto(11,12,13,14,15,16,17,18,19)Iflv 11 A0=0.3636E+01*(1.0 + 0.3122E+00*SB+0.1396E+00*SB2+0.4251E+00*SB3) A1=0.6930E+00-.2574E-01*SB+0.1047E+00*SB2-.2794E-01*SB3 A2=0.3195E+01+0.4045E+00*SB-.3737E+00*SB2-.1677E+00*SB3 A3=0.1009E+00*(1.0 -.1784E+01*SB+0.6263E+00*SB2+0.7337E-01*SB3) $ -1.0 A4=0.2910E+00-.2793E+00*SB+0.6155E-01*SB2+0.5150E-02*SB3 A5=0.0000E+00+0.3185E+00*SB+0.1953E+00*SB2+0.4184E-01*SB3 goto 100 12 A0=0.2851E+00*(1.0 + 0.3617E+00*SB-.4526E+00*SB2+0.5787E-01*SB3) A1=0.2690E+00+0.1104E-01*SB+0.1888E-01*SB2-.1031E-01*SB3 A2=0.3766E+01+0.7850E+00*SB-.3053E+00*SB2+0.1822E+00*SB3 A3=0.2865E+02*(1.0 -.9774E+00*SB+0.5958E+00*SB2-.1234E+00*SB3) $ -1.0 A4=0.8230E+00-.3612E+00*SB+0.5520E-01*SB2+0.1571E-01*SB3 A5=0.0000E+00+0.2145E-01*SB+0.2289E+00*SB2-.4947E-01*SB3 goto 100 13 A0=0.2716E+01*(1.0 -.2092E+01*SB+0.1500E+01*SB2-.3703E+00*SB3) A1=-.3100E-01-.7963E+00*SB+0.1129E+01*SB2-.4191E+00*SB3 A2=0.8015E+01+0.1168E+01*SB-.1625E+01*SB2-.1130E+01*SB3 A3=0.4813E+02*(1.0 -.4951E+00*SB-.8715E+00*SB2+0.5893E+00*SB3) $ -1.0 A4=0.2773E+01-.6329E+00*SB-.1048E+01*SB2+0.1418E+00*SB3 A5=0.0000E+00+0.5048E+00*SB+0.2390E+01*SB2-.4159E+00*SB3 goto 100 14 A0=0.3085E+00*(1.0 + 0.9422E+00*SB-.2606E+01*SB2+0.1364E+01*SB3) A1=0.5000E-02-.6433E+00*SB+0.4980E+00*SB2-.1780E+00*SB3 A2=0.7490E+01+0.9112E+00*SB-.2047E+01*SB2+0.1456E+01*SB3 A3=0.1145E-01*(1.0 + 0.4610E+01*SB+0.1699E+01*SB2+0.1296E+00*SB3) $ -1.0 A4=0.6030E+00-.8081E+00*SB+0.9410E+00*SB2-.4458E+00*SB3 A5=0.0000E+00-.1736E+01*SB+0.2863E+01*SB2-.1268E+01*SB3 goto 100 15 A0=0.1324E+00*(1.0 -.1050E+01*SB+0.4844E+00*SB2-.1043E+00*SB3) A1=-.1580E+00+0.1672E+00*SB-.4100E+00*SB2+0.1793E+00*SB3 A2=0.8559E+01-.7351E-01*SB+0.5898E+00*SB2-.2655E+00*SB3 A3=0.2378E+02*(1.0 -.1108E+00*SB-.1646E-01*SB2+0.1129E-01*SB3) $ -1.0 A4=0.1477E+01+0.3312E-01*SB-.2191E+00*SB2+0.9588E-01*SB3 A5=0.0000E+00+0.1850E+01*SB-.1481E+01*SB2+0.6222E+00*SB3 goto 100 16 A0=0.3208E+00*(1.0 -.4755E+00*SB-.4003E+00*SB2+0.2300E+00*SB3) A1=-.3200E-01-.3357E+00*SB+0.3222E-01*SB2+0.5011E-01*SB3 A2=0.1164E+02+0.1048E+01*SB-.1097E+01*SB2-.4431E+00*SB3 A3=0.5065E+02*(1.0 + 0.2484E+00*SB-.9235E+00*SB2+0.1935E+00*SB3) $ -1.0 A4=0.3300E+01-.6785E+00*SB+0.5337E+00*SB2-.4035E+00*SB3 A5=0.0000E+00-.2496E+00*SB+0.3903E+00*SB2+0.1392E+00*SB3 goto 100 17 A0=0.7967E-06*(1.0 + 0.1587E+01*SB+0.1812E+02*SB2-.1333E+02*SB3) $ *sqrt(sta - stb) A1=0.1096E+01-.1236E+01*SB+0.1014E+02*SB2+0.1940E+01*SB3 A2=0.4366E+00+0.1197E+02*SB-.5471E+00*SB2-.5427E+01*SB3 A3=0.4650E+03*(1.0 + 0.1310E+02*SB-.1918E+02*SB2+0.6791E+01*SB3) $ -1.0 A4=-.8486E+00+0.7457E+00*SB-.1083E+02*SB2-.1210E+01*SB3 A5=0.3494E+01-.3511E+01*SB-.1766E+01*SB2+0.3442E+01*SB3 goto 100 18 A0=0.1713E-03*(1.0 + 0.2562E+02*SB-.2988E+02*SB2+0.4798E+01*SB3) $ *sqrt(sta - stb) A1=-.5276E-01+0.4105E+00*SB-.1079E+01*SB2+0.6278E+00*SB3 A2=0.4515E+01+0.8369E+01*SB-.1192E+02*SB2+0.3403E+01*SB3 A3=0.1756E+01*(1.0 + 0.1325E+02*SB-.2997E+02*SB2+0.1758E+02*SB3) $ -1.0 A4=0.3557E-01+0.4159E+01*SB-.6947E+01*SB2+0.2982E+01*SB3 A5=0.2551E+01+0.2168E+01*SB-.5119E+01*SB2+0.3739E+01*SB3 goto 100 19 A0=0.7510E-04*(1.0 + 0.2836E+02*SB-.3000E+02*SB2-.2979E+02*SB3) $ *sqrt(sta - stb) A1=-.1855E+00+0.4543E+00*SB-.1448E+01*SB2+0.2009E-01*SB3 A2=0.6775E+01-.4210E+01*SB-.1221E+01*SB2+0.1199E+02*SB3 A3=0.1070E+01*(1.0 + 0.8356E+01*SB-.2992E+02*SB2+0.2433E+02*SB3) $ -1.0 A4=-.4601E-01+0.4248E+01*SB-.1736E+01*SB2+0.1187E+02*SB3 A5=0.2771E+01+0.1382E+01*SB-.4797E+01*SB2+0.1273E+01*SB3 goto 100 2 Goto(21,22,23,24,25,26,27,28,29)Iflv 21 A0=0.1828E+01*(1.0 -.8698E+00*SB+0.2906E+00*SB2-.2003E-01*SB3) A1=0.6060E+00+0.8595E-01*SB-.4934E-01*SB2+0.2221E-01*SB3 A2=0.3454E+01-.3115E+00*SB+0.1321E+01*SB2-.3490E+00*SB3 A3=0.2616E+00*(1.0 -.1670E+01*SB+0.2333E+01*SB2+0.7730E-01*SB3) $ -1.0 A4=0.8920E+00-.8500E-02*SB+0.4960E+00*SB2-.4045E-01*SB3 A5=0.0000E+00+0.1091E+01*SB-.1613E+00*SB2+0.3773E-01*SB3 goto 100 22 A0=0.2885E+00*(1.0 + 0.3388E+00*SB-.4550E+00*SB2+0.6005E-01*SB3) A1=0.2730E+00+0.1198E-01*SB+0.1880E-01*SB2-.1077E-01*SB3 A2=0.3736E+01+0.7687E+00*SB-.2731E+00*SB2+0.1638E+00*SB3 A3=0.2741E+02*(1.0 -.9585E+00*SB+0.5925E+00*SB2-.1239E+00*SB3) $ -1.0 A4=0.8040E+00-.3546E+00*SB+0.6123E-01*SB2+0.1086E-01*SB3 A5=0.0000E+00+0.4277E-01*SB+0.2187E+00*SB2-.4646E-01*SB3 goto 100 23 A0=0.8416E-01*(1.0 -.1996E+01*SB+0.1903E+01*SB2-.6722E+00*SB3) A1=-.4790E+00-.5459E+00*SB+0.1638E+01*SB2-.4342E+00*SB3 A2=0.5071E+01+0.1470E+01*SB-.2401E+01*SB2+0.1273E+01*SB3 A3=0.2847E+02*(1.0 + 0.1124E+00*SB-.1338E+01*SB2+0.7115E+00*SB3) $ -1.0 A4=0.4990E+00-.7208E+00*SB+0.3333E-03*SB2-.2354E+00*SB3 A5=0.0000E+00-.4480E+00*SB+0.3720E+01*SB2-.1838E+01*SB3 goto 100 24 A0=0.4378E+00*(1.0 -.1244E+01*SB+0.3278E+01*SB2-.2098E+01*SB3) A1=0.3500E-01-.1298E+01*SB+0.1229E+01*SB2-.3665E+00*SB3 A2=0.6781E+01+0.4078E+01*SB-.9711E+00*SB2-.1536E+01*SB3 A3=0.1527E-03*(1.0 + 0.1430E+02*SB+0.3000E+02*SB2+0.2771E+02*SB3) $ -1.0 A4=0.3060E+00+0.1011E+01*SB-.2045E+01*SB2+0.9422E+00*SB3 A5=0.0000E+00-.3205E+01*SB+0.2683E+01*SB2-.1746E+00*SB3 goto 100 25 A0=0.7413E-01*(1.0 + 0.1291E+01*SB-.2667E+01*SB2+0.1076E+01*SB3) A1=-.2730E+00-.1206E+00*SB+0.1828E+00*SB2-.1001E+00*SB3 A2=0.7719E+01+0.1537E+01*SB-.6410E+00*SB2-.3920E-01*SB3 A3=0.1799E+02*(1.0 -.1334E+01*SB+0.1916E+01*SB2-.8878E+00*SB3) $ -1.0 A4=0.1167E+01-.9176E-01*SB+0.5132E+00*SB2-.3460E+00*SB3 A5=0.0000E+00-.5023E+00*SB+0.1951E+01*SB2-.8427E+00*SB3 goto 100 26 A0=0.6551E+00*(1.0 -.5968E-01*SB+0.5621E-02*SB2-.2074E+00*SB3) A1=0.2800E-01-.1138E+01*SB+0.1178E+01*SB2-.4425E+00*SB3 A2=0.7553E+01+0.3996E+01*SB-.4448E+01*SB2+0.1673E+01*SB3 A3=0.9264E-01*(1.0 -.1760E+01*SB+0.1634E+01*SB2-.4067E+00*SB3) $ -1.0 A4=0.1970E+00+0.5256E+00*SB-.9775E+00*SB2+0.4488E+00*SB3 A5=0.0000E+00-.3668E+01*SB+0.4757E+01*SB2-.1717E+01*SB3 goto 100 27 A0=0.1486E-03*(1.0 + 0.2107E+01*SB-.1056E+02*SB2+0.1403E+02*SB3) $ * sqrt(sta - stb) A1=0.2115E+00-.1702E+01*SB+0.2571E+01*SB2-.1177E+01*SB3 A2=0.3533E+01+0.1367E+01*SB-.3397E+01*SB2+0.6260E+01*SB3 A3=0.1096E+02*(1.0 + 0.9213E+01*SB-.2020E+02*SB2+0.1084E+02*SB3) $ -1.0 A4=0.7041E+00-.7236E+00*SB+0.2766E-01*SB2+0.7352E+00*SB3 A5=0.3904E+01-.4398E+01*SB+0.7056E+01*SB2-.3722E+01*SB3 goto 100 28 A0=0.1201E-03*(1.0 + 0.5408E+01*SB-.1489E+02*SB2+0.1667E+02*SB3) $ * sqrt(sta - stb) A1=0.1420E-01-.1525E+01*SB+0.2408E+01*SB2-.1154E+01*SB3 A2=0.4254E+01+0.2836E+01*SB-.6018E+00*SB2+0.4133E+00*SB3 A3=0.5696E+01*(1.0 + 0.9451E+01*SB-.2029E+02*SB2+0.1033E+02*SB3) $ -1.0 A4=0.4775E+00-.6695E+00*SB+0.2747E+00*SB2-.1051E+00*SB3 A5=0.3330E+01-.5133E+01*SB+0.6921E+01*SB2-.3283E+01*SB3 goto 100 29 A0=0.7697E-04*(1.0 + 0.2801E+02*SB-.1901E+02*SB2-.2880E+02*SB3) $ *sqrt(sta - stb) A1=-.2249E+00+0.4432E+00*SB-.1454E+01*SB2+0.3509E-01*SB3 A2=0.6642E+01-.2702E+01*SB+0.8229E+01*SB2+0.8243E+01*SB3 A3=0.1146E+01*(1.0 + 0.8104E+01*SB-.2998E+02*SB2+0.2812E+02*SB3) $ -1.0 A4=-.6421E-01+0.4246E+01*SB-.2908E+01*SB2+0.9686E-02*SB3 A5=0.2606E+01+0.1261E+01*SB-.4933E+01*SB2+0.3476E+00*SB3 goto 100 3 Goto(31,32,33,34,35,36,37,38,39)Iflv 31 A0=0.3777E+01*(1.0 + 0.6986E+00*SB-.20655E+01*SB2+.10334E+01*SB3) A1=0.7100E+00+.2880E-01*SB-.7930E-01*SB2+0.5600E-01*SB3 A2=0.3259E+01+0.1508E+01*SB-.3932E+01*SB2+0.20613E+01*SB3 A3=0.1304E+00*(1.0 -.2016E+00*SB-.30015E+01*SB2+0.19118E+01*SB3) $ -1.0 A4=0.2890E+00-0.4311E+00*SB+0.7387E+00*SB2-.3697E+00*SB3 A5=0.0000E+00+0.4320E+00*SB+0.2449E+00*SB2-0.6670E-01*SB3 goto 100 32 A0=0.2780E+00*(1.0 + 0.4355E+00*SB-0.4584E+00*SB2+0.4390E-01*SB3) A1=0.2760E+00+0.1420E-01*SB+0.1480E-01*SB2-.9800E-02*SB3 A2=0.3710E+01+0.8250E+00*SB-.3581E+00*SB2+0.1978E+00*SB3 A3=0.2928E+02*(1.0 -.10154E+01*SB+0.6037E+00*SB2-.1175E+00*SB3) $ -1.0 A4=0.8070E+00-.3575E+00*SB+0.4920E-01*SB2+0.1584E-01*SB3 A5=0.0000E+00+0.1860E-01*SB+0.2080E+00*SB2-.450E-01*SB3 goto 100 33 A0=0.2924E+01*(1.0 -.18916E+01*SB+0.1191E+01*SB2-.2492E+00*SB3) A1=0.0000E+00-.9167E+00*SB+0.11147E+01*SB2-.3329E+00*SB3 A2=0.8529E+01+0.7080E+00*SB-.11345E+01*SB2-.10563E+01*SB3 A3=0.1420E+03*(1.0 -.15346E+01*SB+0.7261E+00*SB2-.5730E-01*SB3) $ -1.0 A4=0.3396E+01-.11541E+01*SB-.8834E+00*SB2+0.2430E+00*SB3 A5=0.0000E+00+0.1645E+00*SB+0.19041E+01*SB2+0.1474E+00*SB3 goto 100 34 A0=0.3471E+00*(1.0- 0.1753E+00*SB-.9189E+00*SB2+0.6211E+00*SB3) A1=0.1900E-01-.4579E+00*SB+0.2112E+00*SB2-.6180E-01*SB3 A2=0.7301E+01-.17308E+01*SB+.13666E+01*SB2-.6400E-02*SB3 A3=0.1853E-04*(1.0 -.18260E+02*SB-.2872E+02*SB2-.23456E+02*SB3) $ -1.0 A4=0.4400E+00-.4672E+00*SB+0.6532E+00*SB2-.3222E+00*SB3 A5=0.0000E+00-.4679E+00*SB+0.10741E+01*SB2-.5663E+00*SB3 goto 100 35 A0=0.1702E+00*(1.0 -.1041E+01*SB+0.4064E+00*SB2-.5888E-01*SB3) A1=-.9300E-01-.4742E-01*SB-.1959E+00*SB2+0.1039E+00*SB3 A2=0.9119E+01-.7331E-01*SB+0.3506E+00*SB2-.2081E+00*SB3 A3=0.2981E+02*(1.0 -.1912E+00*SB-.8947E-02*SB2+0.8805E-02*SB3) $ -1.0 A4=0.1668E+01-.6678E-02*SB-.2894E+00*SB2+0.1221E+00*SB3 A5=0.0000E+00+0.1245E+01*SB-.7843E+00*SB2+0.3724E+00*SB3 goto 100 36 A0=0.3910E+00*(1.0 -.1103E+01*SB+0.5383E+00*SB2-.1083E+00*SB3) A1=-.1400E-01-.2471E+00*SB-.8042E-01*SB2+0.7193E-01*SB3 A2=0.9812E+01-.4860E+01*SB+0.5958E+01*SB2-.2342E+01*SB3 A3=0.3749E+00*(1.0 -.3569E+01*SB+0.5456E+01*SB2-.2344E+01*SB3) $ -1.0 A4=0.4940E+00+0.2772E+00*SB-.2732E+00*SB2+0.6466E-01*SB3 A5=0.0000E+00+0.3927E+00*SB-.3216E+00*SB2+0.2164E+00*SB3 goto 100 37 A0=0.3815E-02*(1.0 + 0.2039E+02*SB-.2834E+02*SB2+0.1070E+02*SB3) $ * sqrt(sta - stb) A1=-.2789E-01-.7345E-03*SB-.3251E+00*SB2+0.1946E+00*SB3 A2=0.3223E+01-.4268E+00*SB+0.4387E+01*SB2-.2401E+01*SB3 A3=0.3338E-01*(1.0 -.1163E+02*SB+0.2995E+02*SB2-.1471E+02*SB3) $ -1.0 A4=0.3646E+00-.5767E+00*SB+0.6088E+00*SB2-.2514E+00*SB3 A5=0.1200E+01+0.2178E+00*SB-.4230E+00*SB2+0.4739E+00*SB3 goto 100 38 A0=0.1666E-02*(1.0 + 0.9518E+01*SB-.4715E+01*SB2-.1060E+01*SB3) $ * sqrt(sta - stb) A1=-.1231E+00+0.1656E+00*SB-.5219E+00*SB2+0.2750E+00*SB3 A2=0.3693E+01+0.4922E+01*SB-.1200E+02*SB2+0.7929E+01*SB3 A3=0.1778E+00*(1.0 + 0.3036E+01*SB-.1184E+02*SB2+0.7940E+01*SB3) $ -1.0 A4=0.5353E+00-.1401E+01*SB+0.1970E+01*SB2-.9405E+00*SB3 A5=0.1590E+01+0.1025E+01*SB-.2318E+01*SB2+0.1380E+01*SB3 goto 100 39 A0=0.4319E-03*(1.0 + 0.1100E+02*SB-.9520E+00*SB2+0.1434E+02*SB3) $ * sqrt(sta - stb) A1=-.2512E+00+0.3554E+00*SB-.4120E+00*SB2+0.1328E+00*SB3 A2=0.4764E+01-.3513E+00*SB+0.1199E+02*SB2-.8290E+01*SB3 A3=0.8458E-01*(1.0 + 0.2618E+01*SB+0.4407E+01*SB2+0.2991E+02*SB3) $ -1.0 A4=0.3991E+00-.1363E+01*SB+0.1526E+01*SB2-.3179E+01*SB3 A5=0.1981E+01+0.1496E+01*SB-.1501E+01*SB2+0.3880E+01*SB3 goto 100 4 Goto(41,42,43,44,45,46,47,48,49)Iflv 41 A0=0.1634E+01*(1.0 -.8336E+00*SB+0.1640E+00*SB2+0.1530E+00*SB3) A1=0.5790E+00+0.8587E-01*SB-.6087E-01*SB2+0.1361E-01*SB3 A2=0.2839E+01+0.3720E+00*SB+0.5264E+00*SB2+0.3538E-01*SB3 A3=0.1095E+00*(1.0 -.4830E+00*SB+0.3708E+01*SB2-.6165E+00*SB3) $ -1.0 A4=0.8010E+00-.1432E+00*SB+0.1442E+01*SB2-.1286E+01*SB3 A5=0.0000E+00+0.1035E+01*SB-.5910E-01*SB2-.1982E+00*SB3 goto 100 42 A0=0.3535E+00*(1.0 + 0.4352E+00*SB-.2095E+00*SB2-.8455E-02*SB3) A1=0.2660E+00-.4096E-03*SB+0.1502E-01*SB2-.1163E-01*SB3 A2=0.3514E+01+0.8219E+00*SB-.2330E+00*SB2+0.1055E+00*SB3 A3=0.2200E+02*(1.0 -.9716E+00*SB+0.4552E+00*SB2-.8202E-01*SB3) $ -1.0 A4=0.9000E+00-.3207E+00*SB-.4808E-01*SB2+0.3492E-01*SB3 A5=0.0000E+00-.6273E-01*SB+0.1497E+00*SB2-.5683E-01*SB3 goto 100 43 A0=0.2743E+01*(1.0 -.2027E+01*SB+0.1517E+01*SB2-.4145E+00*SB3) A1=0.7000E-02-.9431E+00*SB+0.1231E+01*SB2-.4834E+00*SB3 A2=0.8200E+01+0.1827E+01*SB-.3453E+01*SB2+0.6763E+00*SB3 A3=0.4975E+02*(1.0 -.2329E+00*SB-.1245E+01*SB2+0.7194E+00*SB3) $ -1.0 A4=0.2387E+01-.4077E+00*SB-.5542E+00*SB2-.9677E-02*SB3 A5=0.0000E+00+0.2702E+00*SB+0.2389E+01*SB2-.8274E+00*SB3 goto 100 44 A0=0.2015E+00*(1.0 -.2133E+00*SB-.6770E+00*SB2+0.5011E+00*SB3) A1=-.7700E-01-.7104E-01*SB-.3720E+00*SB2+0.2159E+00*SB3 A2=0.8008E+01-.2049E+01*SB+0.1800E+01*SB2-.4660E+00*SB3 A3=0.2923E-05*(1.0 + 0.2327E+02*SB+0.1500E+02*SB2+0.2633E+02*SB3) $ -1.0 A4=0.9020E+00-.9191E+00*SB+0.1104E+01*SB2-.5863E+00*SB3 A5=0.0000E+00+0.5840E+00*SB-.8720E+00*SB2+0.4234E+00*SB3 goto 100 45 A0=0.9117E-01*(1.0 -.4089E+00*SB-.4361E+00*SB2+0.2512E+00*SB3) A1=-.2370E+00+0.2492E+00*SB-.3267E+00*SB2+0.1055E+00*SB3 A2=0.8447E+01+0.6009E+00*SB+0.1003E+01*SB2-.1287E+01*SB3 A3=0.3106E+02*(1.0 -.3901E-01*SB+0.1443E+00*SB2-.3433E+00*SB3) $ -1.0 A4=0.1629E+01+0.7855E-01*SB-.1573E+00*SB2-.8595E-01*SB3 A5=0.0000E+00+0.1558E+01*SB-.6295E+00*SB2+0.1847E+00*SB3 goto 100 46 A0=0.3997E+00*(1.0 -.1046E+01*SB+0.6194E+00*SB2-.1342E+00*SB3) A1=0.2000E-02-.2544E+00*SB-.1958E+00*SB2+0.1458E+00*SB3 A2=0.9613E+01-.3919E+01*SB+0.9573E+01*SB2-.5623E+01*SB3 A3=0.3620E+00*(1.0 -.1858E+01*SB+0.8312E+01*SB2-.5900E+01*SB3) $ -1.0 A4=0.3840E+00+0.3572E+00*SB-.1191E+01*SB2+0.7310E+00*SB3 A5=0.0000E+00+0.3351E+00*SB-.7709E+00*SB2+0.4296E+00*SB3 goto 100 47 A0=0.2156E-03*(1.0 + 0.2879E+02*SB-.2310E+02*SB2+0.9812E+01*SB3) $ * sqrt(sta - stb) A1=0.9086E-01-.1250E+00*SB-.7373E-01*SB2-.2201E-01*SB3 A2=0.3588E+01+0.4518E+01*SB-.8930E-01*SB2+0.9163E-02*SB3 A3=0.5216E+01*(1.0 + 0.5912E+00*SB-.4111E+00*SB2+0.7330E+00*SB3) $ -1.0 A4=0.3145E+00+0.1233E+01*SB-.7478E+00*SB2+0.4657E+00*SB3 A5=0.2723E+01-.4110E+00*SB+0.4868E-01*SB2-.3075E+00*SB3 goto 100 48 A0=0.7476E-03*(1.0 + 0.1454E+02*SB-.2509E+02*SB2+0.1184E+02*SB3) $ * sqrt(sta - stb) A1=-.1955E-01-.1712E+00*SB-.1686E+00*SB2+0.2339E+00*SB3 A2=0.4616E+01-.6859E+00*SB-.3959E+01*SB2+0.5530E+01*SB3 A3=0.9881E+01*(1.0 -.1239E+02*SB+0.2721E+02*SB2-.1850E+02*SB3) $ -1.0 A4=0.1200E+02-.1133E+02*SB+0.8138E+01*SB2+0.1199E+02*SB3 A5=0.2226E+01-.5738E+00*SB+0.5239E+00*SB2+0.3825E+00*SB3 goto 100 49 A0=0.8392E-06*(1.0 + 0.1844E+02*SB-.1110E+02*SB2-.2504E+02*SB3) $ * sqrt(sta - stb) A1=0.2127E+00-.5602E+00*SB+0.4777E+01*SB2-.1014E+02*SB3 A2=0.1229E+01+0.7495E+01*SB-.5024E+01*SB2-.1200E+02*SB3 A3=0.2868E+02*(1.0 + 0.7634E+01*SB-.2916E+02*SB2+0.2953E+02*SB3) $ -1.0 A4=0.5970E+00+0.1138E+01*SB-.1439E+01*SB2-.1966E+01*SB3 A5=0.6429E+01-.6673E+00*SB+0.7008E+01*SB2-.1157E+02*SB3 goto 100 5 Goto(51,52,53,54,55,56,57,58,59)Iflv 51 A0= 1.791*(1.0 -0.449*SB-0.445*SB2+ 0.401*SB3) A1= 0.608+ 0.069*SB+ 0.005*SB2-0.037*SB3 A2= 3.470-0.375*SB+ 2.267*SB2-1.261*SB3 A3= 0.315*(1.0 -2.628*SB+ 6.481*SB2-3.834*SB3)-1.0 A4= 1.007-0.732*SB+ 1.490*SB2-0.966*SB3 A5= 0.000+ 0.741*SB+ 0.563*SB2-0.525*SB3 goto 100 52 A0= 0.513*(1.0 + 0.032*SB-0.120*SB2+ 0.013*SB3) A1= 0.276+ 0.052*SB+ 0.000*SB2-0.006*SB3 A2= 3.579+ 0.763*SB-0.135*SB2+ 0.083*SB3 A3= 17.993*(1.0 -0.725*SB+ 0.241*SB2-0.020*SB3)-1.0 A4= 1.120-0.357*SB+ 0.008*SB2+ 0.028*SB3 A5= 0.000+ 0.311*SB+ 0.029*SB2-0.010*SB3 goto 100 53 A0= 2.710*(1.0 -1.773*SB+ 0.970*SB2-0.149*SB3) A1= -0.010-1.636*SB+ 2.087*SB2-0.637*SB3 A2= 7.174+ 2.102*SB-2.209*SB2-0.420*SB3 A3= 29.904*(1.0 -0.756*SB-0.506*SB2+ 0.605*SB3)-1.0 A4= 2.572-0.437*SB-0.968*SB2+ 0.243*SB3 A5= 0.000-1.776*SB+ 4.266*SB2-0.335*SB3 goto 100 54 A0= 0.278*(1.0 - 1.022*SB+ 0.6457*SB2-0.1824*SB3) A1= 0.0862*SB-0.8657*SB2+ 0.4185*SB3 A2= 11.000-1.2809*SB+ 1.2516*SB2+0.061*SB3 A3= 37.338*(1.0 - 0.9404*SB+ 0.2517*SB2+0.1364*SB3)-1.0 A4= 1.960- 0.3385*SB-0.3422*SB2+0.3653*SB3 A5= 0.000+1.424*SB-2.7503*SB2+ 1.2226*SB3 goto 100 55 A0= 0.154*(1.0 -0.659*SB+ 0.005*SB2+ 0.061*SB3) A1= -0.128+ 0.279*SB-0.786*SB2+ 0.363*SB3 A2= 8.649+ 0.071*SB+ 0.351*SB2-0.051*SB3 A3= 43.685*(1.0 -0.603*SB+ 0.037*SB2+ 0.134*SB3)-1.0 A4= 2.238-0.338*SB-0.199*SB2+ 0.157*SB3 A5= 0.000+ 1.681*SB-2.068*SB2+ 0.975*SB3 goto 100 56 A0= 0.372*(1.0 -1.939*SB+ 1.504*SB2-0.440*SB3) A1= 0.009+ 0.610*SB-1.387*SB2+ 0.579*SB3 A2= 10.273-4.833*SB+ 6.583*SB2-2.633*SB3 A3= 0.160*(1.0 + 10.325*SB-2.027*SB2+ 1.571*SB3)-1.0 A4= 0.819-1.660*SB+ 1.845*SB2-0.829*SB3 A5= 0.000+ 3.558*SB-3.940*SB2+ 1.302*SB3 goto 100 57 A0= (7.5242E-5)*(1.0+22.0905*SB+7.1209*SB2-8.303*SB3)* $ sqrt(sta - stb) A1= 0.125-0.3027*SB+0.1564*SB2-0.091*SB3 A2= 2.0388+1.2161*SB+11.5296*SB2-8.0659*SB3 A3= 14.849*(1.0 -2.556*SB+3.5268*SB2-1.6353*SB3)-1.0 A4= 0.3061-0.0901*SB+0.953*SB2-0.4871*SB3 A5= 2.7352+0.1811*SB-0.5167*SB2+0.0543*SB3 goto 100 58 A0= (3.751E-4)*(1.0 + 21.5993*SB+3.1379*SB2-18.8328*SB3)* $ sqrt(sta - stb) A1= -0.0256-0.7717*SB+ 1.1499*SB2-0.5037*SB3 A2= 4.9241+4.0107*SB-4.7012*SB2+0.1097*SB3 A3= 2.842*(1.0 -2.2184*SB+ 2.0293*SB2-0.6907*SB3)-1.0 A4= -0.1352+ 0.8753*SB-1.2626*SB2+ 0.667*SB3 A5= 1.5627-0.4917*SB+ 1.5927*SB2-0.351*SB3 goto 100 59 A0=(2.725E-4)*(1.0 + 18.8497*SB-26.5797*SB2-29.0774*SB3)* $ sqrt(sta - stb) A1= -0.2204-1.0048*SB+0.9415*SB2-0.4274*SB3 A2= 11.034-9.8362*SB-11.1034*SB2-9.1977*SB3 A3= 2.084*(1.0 -2.881*SB+1.2778*SB2-2.9328*SB3)-1.0 A4= -0.0872+ 0.200*SB-1.6187*SB2-1.6058*SB3 A5= 0.8684+4.7047*SB-1.4614*SB2-5.2309*SB3 goto 100 100 CtqPdf = A0*(x**A1)*((1.-x)**A2)*(1.+A3*(x**A4)) $ *(log(1.+1./x))**A5 if(ctqpdf.lt.0.0) ctqpdf = 0.0 Ist = Iset Lp = Iprtn Qsto = QQ Return ENTRY WLAMBD (ISET, IORDER) IORDER = IORD (ISET) WLAMBD = ALM (ISET) RETURN Entry PrCtq1 > (Iset, Iordr, Ischeme, MxFlv, > Alam4, Alam5, Alam6, Amas4, Amas5, Amas6, > Xmin, Qini, Qmax, ExpNor) Iordr = Iord (Iset) Ischeme= Isch (Iset) MxFlv = Nqrk (Iset) Alam4 = Vlm(4,Iset) Alam5 = Vlm(5,Iset) Alam6 = Vlm(6,Iset) Amas4 = Qms(4,Iset) Amas5 = Qms(5,Iset) Amas6 = Qms(6,Iset) Xmin = Xmn (Iset) Qini = Qmn (Iset) Qmax = Qmx (Iset) Do 101 Iexp = 1, Nexp(Iset) ExpNor(Iexp) = ExpN(Iexp, Iset) 101 Continue Return END SUBROUTINE EVLSET (NAME, VALUE, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER*(*) NAME PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF 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 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. 'KF') THEN ITEM = NINT(VALUE) IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12 KF = ITEM Else IRET = 0 EndIf RETURN 12 IRET = 2 RETURN END SUBROUTINE EVLGET (NAME, VALUE, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER*(*) NAME PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF IRET = 1 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. 'KF') THEN VALUE = KF Else IRET = 0 EndIf RETURN ENTRY EVLOUT (NOUUT) WRITE (NOUUT, 131) QINI, IPD0, IHDN, QMAX, IKNL, > XMIN, XCR, NX, NT, JT, KF 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 : Kf (2 * Nf + 1) = ', I8 /) RETURN END SUBROUTINE EVLGT1 (NAME, VALUE, IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX CHARACTER*(*) NAME PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPN = MXF * 2 + 2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) 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, KF IRET = 1 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. 'KF') THEN VALUE = KF Else IRET = 0 EndIf RETURN ENTRY EVLOT1 (NOUUT) WRITE (NOUUT, 131) QINI, IPD0, IHDN, QMAX, IKNL, > XMIN, XCR, NX, NT, JT, KF 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 : Kf (2 * Nf + 1) = ', I8 /) RETURN END SUBROUTINE EVLIN 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 CALL EVLSET ('QINI', V1, IRET1) CALL EVLSET ('IPD0', V2, IRET2) CALL EVLSET ('IHDN', V3, IRET3) goto 301 302 write(nout,120) goto 1 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) 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 CALL EVLSET ('QMAX', V1, IRET1) CALL EVLSET ('NT', V2, IRET2) CALL EVLSET ('JT', V3, IRET3) goto 303 304 write(nout,120) goto 2 303 IF ((IRET1.NE.1) .OR. (IRET2.NE.1) .OR. (IRET3.NE.1)) THEN 22 WRITE (NOUT, 120) GOTO 2 EndIf CALL EVLGET ('XMIN', V1, IR1) CALL EVLGET ('XCR', V2, IR2) CALL EVLGET ('NX', V3, IR3) CALL EVLGET ('KF', 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, KF, IKNL: ', > 2(1PE12.3), 3I6 / '$Type new values: ' ) READ(NIN, *, ERR=22) V1, V2, V3, V4, V5 CALL EVLSET ('XMIN', V1, IRET1) CALL EVLSET ('XCR', V2, IRET2) CALL EVLSET ('NX', V3, IRET3) CALL EVLSET ('KF', V4, IRET4) CALL EVLSET ('IKNL', V5, IRET5) 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 120 FORMAT(' Bad values, Try again!' /) RETURN END FUNCTION XFRMZ (Z) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 105) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / INVERT / ZZ EXTERNAL ZFXL DATA SML, IRUN, TEM, RER / 1E-5, 0, D1, 1E-3 / DATA E1, ZLOW, ZHIGH, IWRN1, IWRN2 / 1.00002, -10.0,1.00002, 2*0 / 7 EPS = TEM * RER ZZ = Z IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN XLA = LOG (XMIN) * 1.5 XLB = 0.00001 TEM = ZBRNT (ZFXL, XLA, XLB, EPS, IRT) Else CALL WARNR (IWRN2, NWRT, 'Z out of range in XFRMZ, X set=0.', > 'Z', Z, ZLOW, ZHIGH, 1) TEM = 0 EndIf XFRMZ = EXP(TEM) RETURN END FUNCTION DXDZ (Z) 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 DATA HUGE, IWRN / 1.E20, 0 / ZZ = Z X = XFRMZ (ZZ) TEM = DZDX (X) IF (TEM .NE. D0) THEN TMP = D1 / TEM Else CALL WARNR(IWRN, NWRT, 'DXDZ is singular in DXDZ; DXDZ set=HUGE', > 'Z', Z, D0, D1, 0) TMP = HUGE EndIf DXDZ = TMP RETURN END FUNCTION ZFRMX (XX) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXX = 105) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX DATA IWRN1,IWRN2,IWRN3,IWRN4, HUGE, TINY / 4*0, 1.E35, 1.E-35 / F(X) = (XCR-XMIN) * LOG (X/XMIN) + LOG (XCR/XMIN) * (X-XMIN) D(X) = (XCR-XMIN) / X + LOG (XCR/XMIN) X = XX IF (X .GE. XMIN) THEN TEM = F(X) / F(D1) ElseIF (X .GE. D0) THEN X = MAX (X, TINY) TEM = F(X) / F(D1) Else CALL WARNR(IWRN1, NWRT, 'X out of range in ZFRMX, Z set to 99.' > , 'X', X, TINY, HUGE, 1) TEM = 99. STOP EndIf ZFRMX = TEM RETURN ENTRY DZDX (XX) X = XX IF (X .GE. XMIN) THEN TEM = D(X) / F(D1) ElseIF (X .GE. D0) THEN X = MAX (X, TINY) TEM = D(X) / F(D1) Else CALL WARNR(IWRN1, NWRT, 'X out of range in DZDX, Z set to 99.' > , 'X', X, TINY, HUGE, 1) TEM = 99. STOP EndIf DZDX = TEM RETURN END SUBROUTINE KERNEL >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT) 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 DATA CF, CG, TR, IWRN / 1.333333333333, 3.0, 0.5, 0 / IRT = 0 TRNF = TR * NFL 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 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) FFP = (1.+ X2) * XM1I FGP = (2.- 2.* X + X2) / X GFP = 1. - 2.* X + 2.* X2 GGP = XM1I + XI - 2. + X - X2 FFM = (1.+ X2) * XP1I FGM = - (2.+ 2.* X + X2) / X GFM = 1. + 2.* X + 2.* X2 GGM = XP1I - XI - 2. - X - X2 FF1 = CF * FFP * (1.- X) FG1 = CF * FGP * X GF1 = 2.* TRNF * GFP GG1 = 2.* CG * GGP * X * (1.-X) 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 PNSP = (PQQ2 + PQQB) * (1.-X) PNSM = (PQQ2 - PQQB) * (1.-X) 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) 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) 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 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 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 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 END FUNCTION PFF1 (XX) 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 = 105, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) PARAMETER (MX = 3) 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, TINY / 2 * .FALSE., 1.E-15 / LB = .TRUE. ENTRY TFF1(ZZ) LA = .TRUE. GOTO 2 ENTRY QFF1 (XX) LB = .TRUE. ENTRY RFF1 (XX) 2 IF (LA .AND. .NOT.LB) THEN Z = ZZ X = XFRMZ (Z) Else X = XX EndIf 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 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 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 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 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 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 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 ENTRY PFF2 (XX) LA = .TRUE. ENTRY QFF2 (XX) LB = .TRUE. ENTRY RFF2 (XX) X = XX 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 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 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 ENTRY PGG2 (XX) LA = .TRUE. ENTRY QGG2 (XX) LB = .TRUE. ENTRY RGG2 (XX) X = XX 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 ENTRY VFF1 (XX) X = XX TEM = (1.+ X**2) / (1.- X) VFF1 = TEM *4./3. RETURN END SUBROUTINE SNRHS (TT, NEFF, FI, GI, FO, GO) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL LSTX PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) COMMON / XXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX COMMON / XYARAY / ZZ (MXX, 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, KF DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX) DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), WM(MXX) DIMENSION R0(MXX), R1(MXX), R2(MXX), RH(MXX), RM(MXX) S = EXP(TT) Q = AL * EXP (S) CPL = ALPI(Q) CPL2= CPL ** 2 / 2. * S CPL = CPL * S CALL INTEGR (NX,-1, FI, WM, IR1) CALL INTEGR (NX, 0, FI, W0, IR2) CALL INTEGR (NX, 1, FI, W1, IR3) CALL INTEGR (NX,-1, GI, RM, IR4) CALL INTEGR (NX, 0, GI, R0, IR5) CALL INTEGR (NX, 1, GI, R1, IR6) CALL INTEGR (NX, 2, GI, R2, IR7) CALL HINTEG (NX, FI, WH) CALL HINTEG (NX, GI, RH) IF (IKNL .GT. 0) THEN DO 230 IZ = 1, NX FO(IZ) = ( 2D0 * FI(IZ) > + 4D0 / 3D0 * ( 2D0 * WH(IZ) - W0(IZ) - W1(IZ) )) > + NEFF * ( R0(IZ) - 2D0 * R1(IZ) + 2D0 * R2(IZ) ) FO(IZ) = FO(IZ) * CPL GO(IZ) = 4D0 / 3D0 * ( 2D0 * WM(IZ) - 2D0 * W0(IZ) + W1(IZ) ) > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ) > + 6D0 * (RH(IZ) + RM(IZ) - 2D0 * R0(IZ) + R1(IZ) - R2(IZ)) GO(IZ) = GO(IZ) * CPL 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 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 IF (IKNL .EQ. 2) THEN DZ = 1./(NX - 1) DO 21 I = 1, NX-1 X = XV(I) NP = NX - I + 1 IS = NP 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 XY = ZZ (IS, IT) 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 TEM1 = SMPNOL (NP, DZ, G1, ERR) TEM2 = SMPSNA (NP, DZ, G2, ERR) TEM3 = SMPSNA (NP, DZ, G3, ERR) TEM4 = SMPNOL (NP, DZ, G4, ERR) TEM1 = TEM1 - FI(I) * (AFF2(I) + AGF2) TEM4 = TEM4 - GI(I) * (AGG2(I) + AFG2) TMF = TEM1 + TEM2 TMG = TEM3 + TEM4 FO(I) = FO(I) + TMF * CPL2 GO(I) = GO(I) + TMG * CPL2 21 CONTINUE EndIf RETURN END SUBROUTINE INTEGR (NX, M, F, G, IR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER MSG*80 PARAMETER (MXX = 105, 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) 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), G(NX) DATA IWRN1, IWRN2 / 0, 0 / IRR = 0 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 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 G(NX) = 0D0 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 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 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 RETURN END SUBROUTINE HINTEG (NX, F, H) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXX = 105, 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) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / HINTEC / GH(NDG, MXX) COMMON / VARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) DIMENSION F(NX), H(NX), G(MXX) DZ = 1D0 / (NX-1) DO 20 I = 1, NX-2 NP = NX - I + 1 TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2) DO 30 KZ = 3, NP IY = I + KZ - 1 W = XA(I,1) / XA(IY,1) G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W) 30 CONTINUE HTEM = SMPSNA (NP-2, DZ, G(3), ERR) TEM1 = F(I) * ELY(I) H(I) = TEM + HTEM + TEM1 20 CONTINUE H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0)) H(NX) = 0 RETURN END FUNCTION ZFXL (XL) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / INVERT / ZZ X = EXP(XL) TT = ZFRMX (X) - ZZ ZFXL = TT RETURN END FUNCTION SPENCE (X) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL SPNINT, SPN2IN COMMON / SPENCC / XX DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 / XX = X TEM = GAUSS(SPNINT, XX, U1, AERR, RERR, ERR, IRT) SPENCE = TEM RETURN ENTRY SPENC2 (X) XX = X TEM = GAUSS (SPN2IN, XX, U1, AERR, RERR, ERR, IRT) SPENC2 = TEM + LOG (XX) ** 2 / 2. RETURN END FUNCTION SPNINT (ZZ) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / SPENCC / X Z = ZZ TEM = LOG (1.- Z) / Z SPNINT = TEM RETURN ENTRY SPN2IN (ZZ) Z = ZZ TEM = LOG (1.+ X - Z) / Z SPN2IN = TEM RETURN END Function SetQCD () IMPLICIT DOUBLE PRECISION (A-H, O-Z) External DatQCD SetQCD = 0. Return END SUBROUTINE PARQCD(IACT,NAME,VALUE,IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) INTEGER IACT,IRET CHARACTER*(*) NAME IRET=1 IF (IACT.EQ.0) THEN WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ', > '(i in 1 to 9), LAMi (i in 1 to NFL)' IRET=4 ELSEIF (IACT.EQ.1) THEN CALL QCDSET (NAME,VALUE,IRET) ELSEIF (IACT.EQ.2) THEN CALL QCDGET (NAME,VALUE,IRET) ELSEIF (IACT.EQ.3) THEN CALL QCDIN IRET=4 ELSEIF (IACT.EQ.4) THEN CALL QCDOUT(NINT(VALUE)) IRET=4 ELSE IRET=3 ENDIF RETURN END SUBROUTINE SETLAM (NEF, WLAM, IRDR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / QCDPAR / AL, NF, NORDER, SET LOGICAL SET IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN WRITE(NOUT,*)'NEF out of range in SETLAM, NEF, NF=', NEF,NF STOP ENDIF VLAM = WLAM IF (IRDR .NE. NORDER) CALL CNVL1 (IRDR, NORDER, NEF, VLAM) CALL SETL1 (NEF, VLAM) END FUNCTION NFL(AMU) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON /QCDPAR/ AL, NF, NORDER, SET LOGICAL SET IF (.NOT. SET) CALL LAMCWZ NFL = NF - NHQ IF ((NFL .EQ. NF) .OR. (AMU .LE. AMN)) GOTO 20 DO 10 I = NF - NHQ + 1, NF IF (AMU .GE. AMHAT(I)) THEN NFL = I ELSE GOTO 20 ENDIF 10 CONTINUE 20 RETURN END FUNCTION AMHATF(I) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / QCDPAR / AL, NF, NORDER, SET COMMON / IOUNIT / NIN, NOUT, NWRT LOGICAL SET IF (.NOT.SET) CALL LAMCWZ IF ((I.LE.0).OR.(I.GT.9)) THEN WRITE (NOUT,*) 'I IS OUT OF RANGE IN AMHATF' AMHATF = 0 ELSE AMHATF = AMHAT(I) ENDIF RETURN END FUNCTION ALPI (AMU) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / QCDPAR / AL, NF, NORDER, SET LOGICAL SET PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15) DATA IW1, IW2 / 2*0 / IF(.NOT.SET) CALL LAMCWZ NEFF = NFL(AMU) ALM = ALAM(NEFF) ALPI = ALPQCD (NORDER, NEFF, AMU/ALM, IRT) IF (IRT .EQ. 1) THEN CALL QWARN (IW1, NWRT, 'AMU < ALAM in ALPI', 'MU', AMU, > ALM, BIG, 1) WRITE (NWRT, '(A,I4,F15.3)') 'NEFF, LAMDA = ', NEFF, ALM ELSEIF (IRT .EQ. 2) THEN CALL QWARN (IW2, NWRT, 'ALPI > 1; Be aware!', 'ALPI', ALPI, > D0, D1, 0) WRITE (NWRT, '(A,I4,2F15.3)') 'NF, LAM, MU= ', NEFF, ALM, AMU ENDIF RETURN END SUBROUTINE QCDSET (NAME,VALUE,IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER*(*) NAME COMMON / COMQMS / VALMAS(9) COMMON / QCDPAR / AL, NF, NORDER, SET LOGICAL SET PARAMETER (PI=3.1415927, EULER=0.57721566) IVALUE = NINT(VALUE) ICODE = NAMQCD(NAME) IF (ICODE .EQ. 0) THEN IRET=0 ELSE IRET = 1 SET = .FALSE. IF (ICODE .EQ. 1) THEN IF (VALUE.LE.0) GOTO 12 AL=VALUE ELSEIF (ICODE .EQ. 2) THEN IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12 NF = IVALUE ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 11)) THEN IF (VALUE .LT. 0) GOTO 12 Scle = Min (Value , VALMAS(ICODE - 2)) AlfScle = Alpi(Scle) * Pi VALMAS(ICODE - 2) = VALUE Call AlfSet (Scle, AlfScle) ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN IF (VALUE .LE. 0) GOTO 12 CALL SETL1 (ICODE-13, VALUE) ELSEIF (ICODE .EQ. 24) THEN IF ((IVALUE .LT. 1) .OR. (IVALUE .GT. 2)) GOTO 12 NORDER = IVALUE ENDIF IF (.NOT. SET) CALL LAMCWZ ENDIF RETURN 12 IRET=2 RETURN END SUBROUTINE QCDGET(NAME,VALUE,IRET) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER*(*) NAME COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / QCDPAR / AL, NF, NORDER, SET COMMON / COMQMS / VALQMS(9) LOGICAL SET PARAMETER (PI=3.1415927, EULER=0.57721566) ICODE = NAMQCD(NAME) IRET = 1 IF (ICODE .EQ. 1) THEN VALUE = AL ELSEIF (ICODE .EQ. 2) THEN VALUE = NF ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 12)) THEN VALUE = VALQMS(ICODE - 2) ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN VALUE = ALAM(ICODE - 13) ELSEIF (ICODE .EQ. 24) THEN VALUE = NORDER ELSE IRET=0 ENDIF END SUBROUTINE QCDIN IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /IOUNIT/ NIN, NOUT, NWRT DIMENSION VALMAS(9) CHARACTER ONECH*(1) ONECH = '0' IASC0 = ICHAR(ONECH) CALL QCDGET ('LAM',ALAM,IRET1) CALL QCDGET ('NFL',ANF,IRET2) CALL QCDGET ('ORDER',ORDER,IRET3) NF = NINT(ANF) NORDER = NINT(ORDER) 1 WRITE (NOUT, *) 'LambdaMSBAR, # Flavors, loop order ?' READ (NIN,*, IOSTAT = IRET) ALAM, NF, NORDER ORDER = NORDER ANF = NF IF (IRET .LT. 0) GOTO 22 IF (IRET .EQ. 0) THEN CALL QCDSET ('LAM',ALAM,IRET1) CALL QCDSET ('NFL',ANF,IRET2) CALL QCDSET ('ORDER',ORDER,IRET3) ENDIF IF ((IRET.NE.0) .OR. (IRET1.NE.1) .OR. (IRET2.NE.1) > .OR. (IRET3.NE.1)) THEN WRITE (NOUT, *) 'Bad value(s), try again.' GOTO 1 ENDIF DO 20, I = 1, NF CALL QCDGET('M'//CHAR(I+IASC0),VALMAS(I),IRET1) 10 WRITE (NOUT, '(1X,A,I2,A)') 'Mass of Quark', I, '?' READ (NIN,*, IOSTAT=IRET) VALMAS(I) IF (IRET .LT. 0) GOTO 22 IF (IRET .EQ. 0) > CALL QCDSET('M'//CHAR(I+IASC0),VALMAS(I),IRET1) IF ((IRET .NE. 0) .OR. (IRET1 .NE. 1)) THEN WRITE (NOUT, *) 'Bad value, try again.' GOTO 10 ENDIF 20 CONTINUE RETURN 22 WRITE (NOUT, *) 'END OF FILE ON INPUT' WRITE (NOUT, *) RETURN END SUBROUTINE QCDOUT(NOUT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /QCDPAR/ AL, NF, NORDER, SET COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON /COMQMS/ VALQMS(9) LOGICAL SET IF (.NOT. SET) CALL LAMCWZ WRITE (NOUT,110) AL, NF, NORDER 110 FORMAT( 1 ' Lambda (MSBAR) =',G13.5,', NFL (total # of Flavors) =',I3, 2 ', Order (loops) =', I2) WRITE (NOUT,120) (I,VALQMS(I),I=1,NF) 120 FORMAT (3(' M', I1, '=', G13.5, :, ',')) IF (NHQ .GT. 0) 1 WRITE (NOUT,130) (I, ALAMF(I), I = NF-NHQ, NF) 130 FORMAT (: ' ! Effective lambda given number of light quarks:'/ > (2(' ! ', I1, ' quarks => lambda = ', G13.5 : '; ')) ) RETURN END SUBROUTINE CNVL1 (IRDR, JRDR, NF, VLAM) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL ZCNVLM COMMON / LAMCNV / AMU, ULAM, NFL, IRD, JRD COMMON / IOUNIT / NIN, NOUT, NWRT DATA ALM, BLM, ERR, AMUMIN / 0.001, 2.0, 0.02, 1.5 / IRD = IRDR JRD = JRDR ULAM = VLAM CALL PARQCD(2, 'NFL', ANF, IRT) NTL = NFLTOT() IF (NF .GT. NTL) THEN WRITE (NOUT, *) ' NF .GT. NTOTAL in CNVL1; set NF = NTOTAL' WRITE (NOUT, *) ' NF, NTOTAL = ', NF, NTL NF = NTL ENDIF NFL = NF AMU = AMHATF(NF) AMU = MAX (AMU, AMUMIN) VLM1 = QZBRNT (ZCNVLM, ALM, BLM, ERR, IR1) IF (NF .LT. NTL) THEN AMU = AMHATF(NF+1) AMU = MAX (AMU, AMUMIN) VLM2 = QZBRNT(ZCNVLM, ALM, BLM, ERR, IR2) ELSE VLM2 = VLM1 ENDIF VLAM = (VLM1 + VLM2) / 2 RETURN END SUBROUTINE SETL1 (NEF, VLAM) IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL SET COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / QCDPAR / AL, NF, NORDER, SET COMMON / COMQMS / QMS(9) COMMON / IOUNIT / NIN, NOUT, NWRT IF (NEF .LT. 0 .OR. NEF .GT. NF) THEN WRITE(NOUT,*)'NEF out of range in SETL1, NEF, NF =',NEF,NF STOP ENDIF AMHAT(0) = 0. DO 5 N = 1, NF AMHAT(N) = QMS(N) 5 CONTINUE ALAM(NEF) = VLAM DO 10 N = NEF, 1, -1 CALL TRNLAM(NORDER, N, -1, IR1) 10 CONTINUE DO 20 N = NEF, NF-1 CALL TRNLAM(NORDER, N, 1, IR1) 20 CONTINUE DO 30, N = NF, 1, -1 IF ((ALAM(N) .GE. 0.7 * AMHAT(N)) > .OR. (ALAM(N-1) .GE. 0.7 * AMHAT(N)))THEN NHQ = NF - N GOTO 40 ENDIF 30 CONTINUE NHQ = NF 40 CONTINUE DO 50, N = NF-NHQ, 1, -1 AMHAT(N) = 0 ALAM(N-1) = ALAM(N) 50 CONTINUE AMN = ALAM(NF) DO 60, N = 0, NF-1 IF (ALAM(N) .GT. AMN) AMN = ALAM(N) 60 CONTINUE AMN = AMN * 1.0001 AL = ALAM(NF) SET = .TRUE. RETURN END SUBROUTINE LAMCWZ IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / QCDPAR / AL, NF, NORDER, SET LOGICAL SET CALL SETL1 (NF, AL) END FUNCTION ALPQCD (IRDR, NF, RML, IRT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15) PARAMETER (CG = 3.0, TR = 0.5, CF = 4.0/3.0) COMMON / IOUNIT / NIN, NOUT, NWRT DATA IW1 / 0/ IRT = 0 IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN WRITE(NOUT, *) 'Order parameter out of range in ALPQCD; > IRDR = ', IRDR IRT = 3 STOP ENDIF B0 = (11.* CG - 2.* NF) / 3. B1 = (34.* CG**2 - 10.* CG * NF - 6.* CF * NF) / 3. RM2 = RML ** 2 IF (RM2 .LE. 1.) THEN IRT = 1 ALPQCD = 99 RETURN ENDIF ALN = LOG (RM2) AL = 4./ B0 / ALN IF (IRDR .GE. 2) AL = AL * (1.- B1 * LOG(ALN) / ALN / B0**2) IF (AL .GE. 1.) THEN IRT = 2 ENDIF ALPQCD = AL RETURN END SUBROUTINE QWARN (IWRN, NWRT1, MSG, NMVAR, VARIAB, > VMIN, VMAX, IACT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /IOUNIT/ NIN, NOUT, NWRT PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) CHARACTER*(*) MSG, NMVAR IW = IWRN VR = VARIAB WRITE (NWRT1,'(I5, 3X,A/ 1X,A,'' = '',1PD16.7)') IW, MSG, > NMVAR, VR IF (IW .EQ. 0) THEN WRITE (NOUT, '(1X, A/1X, A,'' = '', 1PD16.7/A,I4)') > MSG, NMVAR, VR, > ' Complete set of warning messages on file unit #', NWRT1 IF (IACT .EQ. 1) THEN WRITE (NOUT,'(1X,A/2(1PD15.3))')'The limits are: ', VMIN,VMAX WRITE (NWRT1,'(1X,A/2(1PD15.3))')'The limits are: ', VMIN,VMAX ENDIF ENDIF IWRN = IW + 1 RETURN END FUNCTION NAMQCD(NNAME) IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER NNAME*(*), NAME*8 COMMON /IOUNIT/ NIN, NOUT, NWRT COMMON /QCDPAR/ AL, NF, NORDER, SET LOGICAL SET CHARACTER ONECH*(1) ONECH = '0' IASC0 = ICHAR(ONECH) NAME = NNAME NAMQCD=0 IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR. 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') ) 2 NAMQCD=1 IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR. 1 (NAME .EQ. '# FL') ) 2 NAMQCD=2 DO 10 I=1, 9 IF (NAME .EQ. 'M'//CHAR(I+IASC0)) 1 NAMQCD=I+2 10 CONTINUE DO 20 I= 0, NF IF (NAME .EQ. 'LAM'//CHAR(I+IASC0)) 1 NAMQCD=I+13 20 CONTINUE IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') NAMQCD = 24 RETURN END SUBROUTINE ALFSET (QS, ALFS) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL RTALF COMMON / RTALFC / ALFST, JORD, NEFF DATA ALAM, BLAM, ERR / 0.01, 10.0, 0.02 / QST = QS ALFST = ALFS CALL PARQCD (2, 'ORDR', ORDR, IR1) JORD = ORDR NEFF = NFL(QS) EFLLN = QZBRNT (RTALF, ALAM, BLAM, ERR, IR2) EFFLAM = QS / EXP (EFLLN) CALL SETL1 (NEFF, EFFLAM) END FUNCTION ALAMF(N) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / QCDPAR / AL, NF, NORDER, SET COMMON / IOUNIT / NIN, NOUT, NWRT LOGICAL SET IF (.NOT.SET) CALL LAMCWZ IF ((N.LT.0) .OR. (N.GT.9)) THEN WRITE (NOUT, *) ' N IS OUT OF RANGE IN ALAMF' ALAMF=0. ELSE ALAMF = ALAM(MAX(N, NF-NHQ)) ENDIF RETURN END FUNCTION NFLTOT() IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /QCDPAR/ AL, NF, NORDER, SET LOGICAL SET NFLTOT=NF RETURN END FUNCTION ZCNVLM (VLAM) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / LAMCNV / AMU, ULAM, NFL, IRD, JRD ZCNVLM= ALPQCD (IRD,NFL,AMU/ULAM,I) - ALPQCD (JRD,NFL,AMU/VLAM,I) END FUNCTION QZBRNT(FUNC, X1, X2, TOLIN, IRT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /IOUNIT/ NIN, NOUT, NWRT PARAMETER (ITMAX = 1000, EPS = 3.E-12) external func TOL = ABS(TOLIN) A=X1 B=X2 FA=FUNC(A) FB=FUNC(B) IF(FB*FA.GT.0.) THEN WRITE (NOUT, *) 'Root must be bracketed for QZBRNT.' IRT = 1 ENDIF FC=FB DO 11 ITER=1,ITMAX IF(FB*FC.GT.0.) THEN C=A FC=FA D=B-A E=D ENDIF IF(ABS(FC).LT.ABS(FB)) THEN A=B B=C C=A FA=FB FB=FC FC=FA ENDIF TOL1=2.*EPS*ABS(B)+0.5*TOL XM=.5*(C-B) IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN QZBRNT=B RETURN ENDIF IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN S=FB/FA IF(A.EQ.C) THEN P=2.*XM*S Q=1.-S ELSE Q=FA/FC R=FB/FC P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.)) Q=(Q-1.)*(R-1.)*(S-1.) ENDIF IF(P.GT.0.) Q=-Q P=ABS(P) IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN E=D D=P/Q ELSE D=XM E=D ENDIF ELSE D=XM E=D ENDIF A=B FA=FB IF(ABS(D) .GT. TOL1) THEN B=B+D ELSE B=B+SIGN(TOL1,XM) ENDIF FB=FUNC(B) 11 CONTINUE WRITE (NOUT, *) 'QZBRNT exceeding maximum iterations.' IRT = 2 QZBRNT=B RETURN END SUBROUTINE TRNLAM (IRDR, NF, IACT, IRT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / IOUNIT / NIN, NOUT, NWRT COMMON / CWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ COMMON / TRNCOM / VMULM, JRDR, N, N1 EXTERNAL ZBRLAM DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 / DATA IW1, IW2, IW3, IW4, IR1, SML / 4*0, 0, 1.E-5 / IRT = 0 N = NF JRDR = IRDR JACT = IACT VLAM = ALAM(N) IF (JACT .GT. 0) THEN N1 = N + 1 THMS = AMHAT(N1) ALM = LOG (THMS/VLAM) BLM = BLM0 ELSE N1 = N -1 THMS = AMHAT(N) ALM = ALM0 THMS = MAX (THMS, SML) BLM = LOG (THMS/VLAM) ENDIF IF (VLAM .GE. 0.7 * THMS) THEN IF (JACT . EQ. 1) THEN AMHAT(N1) = 0 ELSE AMHAT(N) = 0 ENDIF IRT = 4 ALAM(N1) = VLAM RETURN ENDIF IF (ALM .GE. BLM) THEN WRITE (NOUT, *) 'TRNLAM has ALM >= BLM: ', ALM, BLM WRITE (NOUT, *) 'I do not know how to continue' STOP ENDIF VMULM = THMS/VLAM ERR = RERR * LOG (VMULM) WLLN = QZBRNT (ZBRLAM, ALM, BLM, ERR, IR1) ALAM(N1) = THMS / EXP (WLLN) IF (IR1 .NE. 0) THEN WRITE (NOUT, *) 'QZBRNT failed to find VLAM in TRNLAM; ', > 'NF, VLAM =', NF, VLAM WRITE (NOUT, *) 'I do not know how to continue' STOP ENDIF RETURN END FUNCTION RTALF (EFLLN) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (PI = 3.141592653589) COMMON / RTALFC / ALFST, JORD, NEFF EFMULM = EXP (EFLLN) TEM1 = PI / ALFST TEM2 = 1. / ALPQCD (JORD, NEFF, EFMULM, I) RTALF = TEM1 - TEM2 END FUNCTION ZBRLAM (WLLN) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON / TRNCOM / VMULM, JRDR, N, N1 WMULM = EXP (WLLN) TEM1 = 1./ ALPQCD(JRDR, N1, WMULM, I) TEM2 = 1./ ALPQCD(JRDR, N, VMULM, I) ZBRLAM = TEM1 - TEM2 END SUBROUTINE TRMSTR(STRING,ILEN) CHARACTER STRING*(*),SPACE*1 DATA SPACE/' '/ ILEN=0 IF (STRING.EQ.SPACE) RETURN 1 IF (STRING(1:1).NE.SPACE) GOTO 2 STRING=STRING(2:) GOTO 1 2 CONTINUE DO 3 I=LEN(STRING),1,-1 IF (STRING(I:I).NE.SPACE) THEN ILEN=I RETURN END IF 3 CONTINUE END SUBROUTINE WARNI (IWRN, NWRT, MSG, NMVAR, IVAB, > IMIN, IMAX, IACT) CHARACTER*(*) MSG, NMVAR Data Nmax / 1000 / IW = IWRN IV = IVAB IF (IW .EQ. 0) THEN PRINT '(1X,A/1X, 2A,I10 /A,I4)', MSG, NMVAR, ' = ', IV, > ' For all warning messages, check file unit #', NWRT IF (IACT .EQ. 1) THEN PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX WRITE (NWRT,'(A/2I10)') ' The limits are: ', IMIN, IMAX ENDIF ENDIF If (Iw .LT. Nmax) Then WRITE (NWRT,'(1X,A/1X,2A, I10)') MSG, NMVAR, ' = ', IV Elseif (Iw .Eq. Nmax) Then Print '(/A/)', '!!! Severe Warning, Too many errors !!!' Print '(/A/)', ' !!! Check The Error File !!!' Write (Nwrt, '(//A//)') > 'Too many warnings, Message suppressed !!' Endif IWRN = IW + 1 RETURN END SUBROUTINE UpC (A, La, UpA) CHARACTER A*(*), UpA*(*), C*(1) INTEGER I, La, Ld La = Len(A) Lb = Len(UpA) If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!' Ld = ICHAR('A')-ICHAR('a') DO 1 I = 1, Lb If (I .Le. La) Then c = A(I:I) IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN UpA (I:I) = CHAR(Ichar(c) + ld) Else UpA (I:I) = C ENDIF Else UpA (I:I) = ' ' Endif 1 CONTINUE RETURN END SUBROUTINE WARNR (IWRN, NWRT, MSG, NMVAR, VARIAB, > VMIN, VMAX, IACT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) CHARACTER*(*) MSG, NMVAR Data Nmax / 1000 / IW = IWRN VR = VARIAB IF (IW .EQ. 0) THEN PRINT '(1X, A/1X,2A,1PD16.7/A,I4)', MSG, NMVAR, ' = ', VR, > ' For all warning messages, check file unit #', NWRT IF (IACT .EQ. 1) THEN PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX WRITE (NWRT,'(A/2(1PE15.4))') ' The limits are: ', VMIN, VMAX ENDIF ENDIF If (Iw .LT. Nmax) Then WRITE (NWRT,'(I5, 2A/1X,2A,1PD16.7)') IW, ' ', MSG, > NMVAR, ' = ', VR Elseif (Iw .Eq. Nmax) Then Print '(/A/)', '!!! Severe Warning, Too many errors !!!' Print '(/A/)', ' !!! Check The Error File !!!' Write (Nwrt, '(//A//)') > '!! Too many warnings, Message suppressed from now on !!' Endif IWRN = IW + 1 RETURN END SUBROUTINE POLINT (XA,YA,N,X,Y,DY) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (NMAX=10) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END SUBROUTINE RATINT(XA,YA,N,X,Y,DY) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (NMAX=10,TINY=1.E-25) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 HH=ABS(X-XA(1)) DO 11 I=1,N H=ABS(X-XA(I)) IF (H.EQ.0.)THEN Y=YA(I) DY=0.0 RETURN ELSE IF (H.LT.HH) THEN NS=I HH=H ENDIF C(I)=YA(I) D(I)=YA(I)+TINY 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M W=C(I+1)-D(I) H=XA(I+M)-X T=(XA(I)-X)*D(I)/H DD=T-C(I+1) IF(DD.EQ.0.)PAUSE DD=W/DD D(I)=C(I+1)*DD C(I)=T*DD 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END FUNCTION GAUSS(F,XL,XR,AERR,RERR,ERR,IRT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION XLIMS(100), R(93), W(93) INTEGER PTR(4),NORD(4),NIN,NOUT,NWRT external f COMMON/IOUNIT/NIN,NOUT,NWRT DATA PTR,NORD/4,10,22,46, 6,12,24,48/ DATA R/.2386191860,.6612093865,.9324695142, 1 .1252334085,.3678314990,.5873179543,.7699026742,.9041172563, 1 .9815606342,.0640568929,.1911188675,.3150426797,.4337935076, 1 .5454214714,.6480936519,.7401241916,.8200019860,.8864155270, 1 .9382745520,.9747285560,.9951872200,.0323801710,.0970046992, 1 .1612223561,.2247637903,.2873624873,.3487558863,.4086864820, 1 .4669029048,.5231609747,.5772247261,.6288673968,.6778723796, 1 .7240341309,.7671590325,.8070662040,.8435882616,.8765720203, 1 .9058791367,.9313866907,.9529877032,.9705915925,.9841245837, 1 .9935301723,.9987710073,.0162767488,.0488129851,.0812974955, 1 .1136958501,.1459737146,.1780968824,.2100313105,.2417431561, 1 .2731988126,.3043649444,.3352085229,.3656968614,.3957976498, 1 .4254789884,.4547094222,.4834579739,.5116941772,.5393881083, 1 .5665104186,.5930323648,.6189258401,.6441634037,.6687183100, 1 .6925645366,.7156768123,.7380306437,.7596023411,.7803690438, 1 .8003087441,.8194003107,.8376235112,.8549590334,.8713885059, 1 .8868945174,.9014606353,.9150714231,.9277124567,.9393703398, 1 .9500327178,.9596882914,.9683268285,.9759391746,.9825172636, 1 .9880541263,.9925439003,.9959818430,.9983643759,.9996895039/ DATA W/.4679139346,.3607615730,.1713244924, 1 .2491470458,.2334925365,.2031674267,.1600783285,.1069393260, 1 .0471753364,.1279381953,.1258374563,.1216704729,.1155056681, 1 .1074442701,.0976186521,.0861901615,.0733464814,.0592985849, 1 .0442774388,.0285313886,.0123412298,.0647376968,.0644661644, 1 .0639242386,.0631141923,.0620394232,.0607044392,.0591148397, 1 .0572772921,.0551995037,.0528901894,.0503590356,.0476166585, 1 .0446745609,.0415450829,.0382413511,.0347772226,.0311672278, 1 .0274265097,.0235707608,.0196161605,.0155793157,.0114772346, 1 .0073275539,.0031533461,.0325506145,.0325161187,.0324471637, 1 .0323438226,.0322062048,.0320344562,.0318287589,.0315893308, 1 .0313164256,.0310103326,.0306713761,.0302999154,.0298963441, 1 .0294610900,.0289946142,.0284974111,.0279700076,.0274129627, 1 .0268268667,.0262123407,.0255700360,.0249006332,.0242048418, 1 .0234833991,.0227370697,.0219666444,.0211729399,.0203567972, 1 .0195190811,.0186606796,.0177825023,.0168854799,.0159705629, 1 .0150387210,.0140909418,.0131282296,.0121516047,.0111621020, 1 .0101607705,.0091486712,.0081268769,.0070964708,.0060585455, 1 .0050142027,.0039645543,.0029107318,.0018539608,.0007967921/ DATA TOLABS,TOLREL,NMAX/1.E-35,5.E-4,100/ TOLABS=AERR TOLREL=RERR GAUSS=0. NLIMS=2 XLIMS(1)=XL XLIMS(2)=XR 10 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0 BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0 TVAL=0. DO 15 I=1,3 15 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I))) TVAL=TVAL*AA DO 25 J=1,4 VAL=0. DO 20 I=PTR(J),PTR(J)-1+NORD(J) 20 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I))) VAL=VAL*AA TOL=MAX(TOLABS,TOLREL*ABS(VAL)) IF (ABS(TVAL-VAL).LT.TOL) THEN GAUSS=GAUSS+VAL NLIMS=NLIMS-2 IF (NLIMS.NE.0) GO TO 10 RETURN END IF 25 TVAL=VAL IF (NMAX.EQ.2) THEN GAUSS=VAL RETURN END IF IF (NLIMS.GT.(NMAX-2)) THEN WRITE(NOUT,50) GAUSS,NMAX,BB-AA,BB+AA RETURN END IF XLIMS(NLIMS+1)=BB XLIMS(NLIMS+2)=BB+AA XLIMS(NLIMS)=BB NLIMS=NLIMS+2 GO TO 10 50 FORMAT (' GAUSS FAILS, GAUSS,NMAX,XL,XR=',G15.7,I5,2G15.7) END FUNCTION SMPNOL (NX, DX, FN, ERR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION FN(NX) MS = MOD(NX, 2) IF (NX .LE. 1 .OR. NX .GT. 1000) THEN PRINT *, 'NX =', NX, ' OUT OF RANGE IN SMPNOL!' STOP ELSEIF (NX .EQ. 2) THEN TEM = DX * FN(2) ELSEIF (NX .EQ. 3) THEN TEM = DX * FN(2) * 2. ELSE IF (MS .EQ. 0) THEN TEM = DX * (23.* FN(2) - 16.* FN(3) + 5.* FN(4)) / 12. TMP = DX * (3.* FN(2) - FN(3)) / 2. ERR = ABS(TEM - TMP) TEM = TEM + SMPSNA (NX-1, DX, FN(2), ER1) ERR = ABS(ER1) + ERR ELSE TEM = DX * (8.* FN(2) - 4.* FN(3) + 8.* FN(4)) / 3. TMP = DX * (3.* FN(2) + 2.* FN(3) + 3.* FN(4)) / 2. ERR = ABS(TEM - TMP) TEM = TEM + SMPSNA (NX-4, DX, FN(5), ER1) ERR = ABS(ER1) + ERR ENDIF ENDIF SMPNOL = TEM RETURN END FUNCTION ZBRNT(FUNC, X1, X2, TOL, IRT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (ITMAX = 1000, EPS = 3.E-12) external func IRT = 0 TOL = ABS(TOL) A=X1 B=X2 FA=FUNC(A) FB=FUNC(B) IF(FB*FA.GT.0.) THEN PRINT *, 'Root must be bracketed for ZBRNT.' IRT = 1 RETURN ENDIF FC=FB DO 11 ITER=1,ITMAX IF(FB*FC.GT.0.) THEN C=A FC=FA D=B-A E=D ENDIF IF(ABS(FC).LT.ABS(FB)) THEN A=B B=C C=A FA=FB FB=FC FC=FA ENDIF TOL1=2.*EPS*ABS(B)+0.5*TOL XM=.5*(C-B) IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN ZBRNT=B RETURN ENDIF IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN S=FB/FA IF(A.EQ.C) THEN P=2.*XM*S Q=1.-S ELSE Q=FA/FC R=FB/FC P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.)) Q=(Q-1.)*(R-1.)*(S-1.) ENDIF IF(P.GT.0.) Q=-Q P=ABS(P) IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN E=D D=P/Q ELSE D=XM E=D ENDIF ELSE D=XM E=D ENDIF A=B FA=FB IF(ABS(D) .GT. TOL1) THEN B=B+D ELSE B=B+SIGN(TOL1,XM) ENDIF FB=FUNC(B) 11 CONTINUE PRINT *, 'ZBRNT exceeding maximum iterations.' IRT = 2 ZBRNT=B RETURN END FUNCTION FINTRP (FF, X0, DX, NX, XV, ERR, IR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MX = 3) DIMENSION FF (0:NX), XX(MX) COMMON / IOUNIT / NIN, NOUT, NWRT DATA SML, HUGE, XX / 1.D-5, 1.D30, 0., 1.0, 2.0 / DATA IW1, IW2, IW3, IW4, IW5, IW6, IW7 / 7 * 0 / IR = 0 X = XV ERR = 0. ANX = NX FINTRP = 0. IF (NX .LT. 1) THEN CALL WARNI(IW1, NWRT, 'Nx < 1, error in FINTRP.', > 'NX', NX, 1, 256, 1) IR = 1 RETURN ELSE MNX = MIN(NX+1, MX) ENDIF IF (DX .LE. 0) THEN CALL WARNR(IW3, NWRT, 'DX < 0, error in FINTRP.', > 'DX', DX, D0, D1, 1) IR = 2 RETURN ENDIF XM = X0 + DX * NX IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN CALL WARNR(IW5,NWRT,'X out of range in FINTRP, > Extrapolation used.','X',X,X0,XM,1) IR = 3 ENDIF TX = (X - X0) / DX IF (TX .LE. 1.) THEN IX = 0 ELSEIF (TX .GE. ANX-1.) THEN IX = NX - 2 ELSE IX = TX ENDIF DDX = TX - IX CALL RATINT (XX, FF(IX), MNX, DDX, TEM, ERR) FINTRP = TEM RETURN END FUNCTION SMPSNA (NX, DX, F, ERR) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MAXX = 1000) COMMON / IOUNIT / NIN, NOUT, NWRT DIMENSION F(NX) DATA IW1, IW2, TINY / 2*0, 1.E-35 / IF (DX .LE. 0.) THEN CALL WARNR(IW2,NWRT,'DX cannot be < 0. in SMPSNA', 'DX', DX, > D0, D1, 0) SMPSNA = 0. RETURN ENDIF IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN CALL WARNI(IW1, NWRT, 'NX out of range in SMPSNA', 'NX', NX, > 1, MAXX, 1) SIMP = 0. ELSEIF (NX .EQ. 1) THEN SIMP = 0. ELSEIF (NX .EQ. 2) THEN SIMP = (F(1) + F(2)) / 2. ERRD = (F(1) - F(2)) / 2. ELSE MS = MOD(NX, 2) IF (MS .EQ. 0) THEN ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24. NZ = NX - 1 ELSE ADD = 0. NZ = NX ENDIF IF (NZ .EQ. 3) THEN SIMP = (F(1) + 4.* F(2) + F(3)) / 3. TRPZ = (F(1) + 2.* F(2) + F(3)) / 2. ELSE SE = F(2) SO = 0 NM1 = NZ - 1 DO 60 I = 4, NM1, 2 IM1 = I - 1 SE = SE + F(I) SO = SO + F(IM1) 60 CONTINUE SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3. TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2. ENDIF ERRD = TRPZ - SIMP SIMP = SIMP + ADD ENDIF SMPSNA = SIMP * DX IF (ABS(SIMP) .GT. TINY) THEN ERR = ERRD / SIMP ELSE ERR = 0. ENDIF RETURN END FUNCTION FitIni (LP, XX) Implicit Double Precision (A-H, O-Z) Parameter (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) Parameter (NF0 = 4, Nshp = 4, NCP = 4, NEX = Nshp+2) Common > /InpFin/ A(Nshp), B(0:NEX, -NF0:2), Ifun DATA D1M, SML / 0.99999999, 1.0D-10 / Fa(X,L) = X**(B(1,L)-1.) *(1.-X)**B(2,L) F1(X,L) = 1.+ (Exp(B(3,L)) - 1.) * x**B(4,L) F2(x,L) = (LOG(1.+ 1./x)) **B(3,L) F3(x,L) = x**(B(1,L) - 1.) *((1.-x)**B(2,L)) $ *(1. + B(3,L)*(x**.5) + B(4,L)*x)*B(0,L) F5(x,L) = 1. + B(3,L)*(x**.5) + B(4,L)*x IF (LP .LT. -NF0 .OR. LP .GT. NF0) Then FitIni = 0. Return EndIf X = XX IF (x .LE. SML .OR. x .GE. D1M) Then FitIni = 0. Return EndIf if (lp.eq.1.or.lp.eq.2) then if (ifun .eq. 1.or. ifun.eq.3 .or. ifun.eq.4) then FV = Fa(x, LP) * F1(x, LP) ElseIf (Ifun .Eq. 2) Then FV = Fa(x, LP) * F2(x, LP) ElseIf (Ifun .Eq. 5) Then FV = Fa(x, LP) * F5(x, LP) Else Print *, 'Ifun = 1...5 in FitIni! Not ', Ifun Stop Endif else fv = 0. endif IF (LP .LE. 0) Then VL = 0. ElseIf (LP .LE. 2) Then VL = FV * B(0, LP) Else VL = 0. EndIf KP = -ABS (LP) If (Ifun .Eq. 1) Then FS = Fa(x, KP) * F1(x, KP)* B(0, KP) ElseIf (Ifun .Eq. 2) Then FS = Fa(x, KP) * F2(x, KP)* B(0, KP) ElseIf (Ifun.eq.3) then if(kp.eq.-1 .or. kp.eq.-2) then FS = .5D0 * fa(x,-2) * f1(x, -2) * b(0,-2) > + ( -kp -1.5D0) * F3(x,-1) Else FS = Fa(x, KP) * F1(x, KP)* B(0, KP) endif ElseIf (Ifun.eq.4) then if(kp.eq.-1 .or. kp.eq.-2) then FS = .5D0 * fa(x,-2) * f1(x, -2) * b(0,-2) > + ( -kp -1.5D0) * F3(x,-1) Elseif( kp.eq.0) then FS = F3(x,kp) Else FS = Fa(x, KP) * F1(x, KP)* B(0, KP) endif ElseIf (Ifun.eq.5) then If ( kp.eq.-1) then FS = F3(x,-2)*(1.D0-B(0,-4))*0.2D0-.5D0*F3(x,-1) ElseIf ( kp.eq.-2) then FS = F3(x,-2)*(1.D0-B(0,-4))*0.2D0+.5D0*F3(x,-1) ElseIf ( kp.eq.-3) then FS = F3(x,kp)*(1.D0-B(0,-4))*0.1D0 ElseIf ( kp.eq.-4) then FS = F3(x,kp)* B(0,-3) * 0.5D0 Elseif( kp.eq.0) then FS = F3(x,kp) endif Endif SEA = FS FitIni = VL + SEA Return End BLOCK DATA DATPDF IMPLICIT DOUBLE PRECISION (A-H, O-Z) CHARACTER*10 NAMPDF, NMHDRN, NAMQRK CHARACTER*12 MRSFLN LOGICAL LSTX PARAMETER (Z = 1D -10, ZZ = 1D -20) PARAMETER (MXX = 105, MXQ = 25, MXF = 6, MXHDRN = 6) PARAMETER (MXPN = MXF*2+2) PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) PARAMETER (MXPDF = 20) 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, KF COMMON / PEVLDT / UPD(MXPQX) 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) COMMON / NAMPDF / NAMPDF(MXPDF), NMHDRN(-1 : MXHDRN), > NAMQRK(-MXF : MXF) COMMON / MRSDAT / MRSFLN (3) DATA QINI, QMAX, XMIN, XCR / 1.9, 1.001D2, 0.999D-3, 1.5 / 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 / 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 / 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' / DATA (NMHDRN(I), I=-1,5)/ 'AntiProton', 'Dummy', 'Proton', > 'Neutron', 'Pi_Plus', 'Pi_Minus', 'K_Plus' / 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' / 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 / END BLOCK DATA DATQCD IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /COMQCH/ VALQCH(9) COMMON /COMQMS/ VALQMS(9) COMMON /QCDPAR/ AL, NF, NORDER, SET COMMON /COMALP/ ALPHA LOGICAL SET DATA AL, NF, NORDER, SET / .20, 5, 2, .FALSE. / DATA VALQCH/ 0.66666667, -0.33333333, > -0.33333333, 0.66666667, > -0.33333333, 0.66666667, > 3*0./ DATA VALQMS/ 2*0., .2, 1.6, 5., 180., 3*0./ DATA ALPHA/ 7.29927E-3 / END