c====================================================================== Subroutine Structm(xin,qin,Upv,Dnv,Usea,Dsea,Str,Chm,Bot,Top,Glu) c====================================================================== c c Routine to Provide Identical PDFLIB outputs using CTEQ PDFs not c in the present CERNLIB release. c c Required: The appropriate *.ini file in your directory while running. c The .ini file is defined in the data statement below. c c Note: This completely ignores the normal PDFSET settings, you c can just leave your code as is and link in this routine. c In fact, this provides a dummy PDFSET so you can remove the c link to PDFLIB completely if you wish. In the dummy PDFSET c the value of alphas is set to the appropriate value for the c pdf being run, and is inserted into the PDFLIB common block c so that the PDFLIB alphas routine will work. c c implicit double precision (a-h, o-z) dimension qmss(3) logical ifirst data ifirst/.true./ character*60 inifile c ********************************************************************** c Change the following Data Statement To Use A Different PDF c c Now Available In This Format: c CTEQ3M Lambda(5 Flavors) = 0.158 *Old c CTEQ4A1 = 0.140 *Alphas-Scan c CTEQ4A2 = 0.169 *Alphas-Scan c CTEQ4M=CTEQ4A3 = 0.202 *Default c CTEQ4A4 = 0.239 *Alphas-Scan c CTEQ4A5 = 0.282 *Alphas-Scan c CTEQ4HJ = 0.206 *High-Pt Jet Fit c CTEQ4L = 0.181 *Leading-Order c ********************************************************************** Data inifile/'cteq.ini'/ if(ifirst) then ifirst = .false. qmss(1) = 1.6d0 qmss(2) = 5.0d0 qmss(3) = 175.0d0 open(unit=3,status='old',file=inifile) write(*,*) write(*,*) 'CTEQ Parameter File Opened',inifile xmin = 0.00001 qmax = 4000. call evlini(3,xmin,qmax,qmss,100,20,2) close(3) write(*,*) write(*,*) ' PDF evolution complete' write(*,*) endif x = xin qmaxn = qin Glu = x*pardis(0,x,qmaxn) Dsea = x*pardis(-2,x,qmaxn) Usea = x*pardis(-1,x,qmaxn) Dnv = x*pardis(2,x,qmaxn) - Dsea Upv = x*pardis(1,x,qmaxn) - Usea Str = x*pardis(3,x,qmaxn) Chm = x*pardis(4,x,qmaxn) Bot = x*pardis(5,x,qmaxn) Top = x*pardis(6,x,qmaxn) End c=================================================================== SUBROUTINE PDFSET(PARM,VAL) c=================================================================== c c dummy version to call cteq version of structm once in order to c get the appropriate alphas c DOUBLE PRECISION VAL(20) CHARACTER*20 PARM(20) DOUBLE PRECISION XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU c XX = 0.1 QQ = 100. CALL STRUCTM(XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU) RETURN END C======================================================================== Subroutine Evlini (Nini, Xmin,Qmax,Qmss, Nx,Nq,Lpr) C======================================================================== 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 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) csek 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 csek DATA D0, Ahdn / 0.0, 1.0 / DATA Ahdn / 1.0 / C ---------------------------- C Initialize pdf & Qcd parameters Dumm = SetPdf() Dumm = SetQcd() csek Print *, 'Reading input XXX.ini.....' C Nini is the input Xnnn.ini file Call iniread(Nini, qini, anfl, ordn, ordi, blam, Qm4, Qm5) csek If (Lpr .GE. 2) Print '(A/6F9.3, I5)', csek >' Nini, qini, anfl, ordn, ordi, blam, ifun=', csek > 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 csek Print '(A/A)', 'Initial distribution parameters:' csek > ,' N A1 A2 A3 A4' csek do 11 ifl = 2,-nf0,-1 csek 11 print '(1X, 5(1pE11.3))', (b(iex,ifl),iex=0,4) C Print *, ' ' Print *, 'Start evolution ....' Print *, ' ' C ---------------------------- Anout = 6 Call ParPdf (4, 'DUM', ANout, Jr) Call Evolve (FitIni, Irt) Return END C **************************** C **************************** 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 c============================================================================ FUNCTION FitIni (LP, XX) c============================================================================ C Initial Parton distr. at Q = Q0 C Ifun = valence sea gluon ! description C 1 Fa * F1 Fa * F1 Fa * F1 ! CTEQorg:(dbar,ubar) F1 C 2 Fa * F2 Fa * F2 Fa * F2 ! M-T C 3 Fa * F1 Fa * F1,3 Fa * F1 ! CTEQstd:(dbar-ubar) F3 C 4 Fa * F1 Fa * F1,3 Fa * F3 ! " w. gluon F3 C 5 Fa * F3 Fa * F** Fa * F3 ! MRS std** C 6 Fa * F6 Fa * F3,6 Fa * F6 ! ?? C 7 Fa * F1 Fa * F1,3 Fa * F7 ! CTEQstd w. gluon F7 C 8 Fa * F1 Fa * F1,3 Fa * F8 ! CTEQstd w. gluon F8 C 9 Fa * F1 Fa * F1 Fa * F1 ! CTEQstd:(dbar/ubar) F1 C 10 Fa * F1 Fa * F1 Fa * F1 ! CTEQstd:(dbar/ubar) F3 C 11 Fa * F1 Fa * F1,3 Fa * F1 +F11 ! distort glu C 12 Fa * F1 Fa * F1 Fa * F1 ! CTEQstd:(dbar/ubar) B-Wz C 13 Fa * F1 Fa * F1 Fa * F1 ! CTEQstd:(dbar/ubar)B-Wz+F3 C 14 Fa * F1 Fa * F1 Fa * F1 ! CTEQstd:(dbar/ubar)B-Wx+F3 C 101 Gf Fa * F Fa * F ! Gaussian Frg C 102 Pf Fa * F Fa * F ! Peterson Frg 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,A1,A2) = x**(A1-1.) * (1.-x)**A2 F1 (x,A3,A4) = 1. + (Exp(A3)-1.) * x**A4 ! CTEQstd F2 (x,A3) = (LOG(1.+1./x))**A3 ! M-T F3 (x,A3,A4) = 1.D0 + A3* x**0.5D0 + A4* x ! MRS F6 (x,x0,Del) = Exp( - (x-x0)**2/ 2D0/ Del**2 ) ! Gaussian around F7 (x,A3,A4) = 1.D0 + Exp(A3 + A4 * x) ! Exp1 F8 (x,A3,A4) = 1.D0 + A3* x + Exp(A4 * x) ! Exp2 F11(x,A5,A6) = A5*(1.-x)**A6 *(x - 1./(A6+2.)) / x ! MomConserve (Soper) F12(z,z0,Del) = > 1D0 / (1D0 + ((z - Exp(z0))/Exp(Del))**2) ! Breit-Wigner C c Gaussian Gf (x,d,x0) = exp(- (x-x0)**2 * d**2) C Peterson function Pf (x,e) = 1./ (1.- 1./x- e/(1.-x))**2 / 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 C Valence if (lp.eq.1.or.lp.eq.2) then FV = Fa(x,B(1,LP),B(2,LP)) if (ifun .eq. 1.or. ifun.eq.3 .or. ifun.eq.4 > .or. ifun.ge.7 .or. ifun.le.15) then FV = FV * F1(x, B(3,LP),B(4,LP)) ElseIf (Ifun .Eq. 2) Then FV = FV * F2(x, B(3,LP)) ElseIf (Ifun .Eq. 5) Then FV = FV * F3(x, B(3,LP),B(4,LP)) ElseIf (Ifun .Eq. 6) Then FV = FV * F6(x, B(3,LP),B(4,LP)) ElseIf (Ifun .Eq. 101) Then C Gaussian FV = Gf(x, B(3,Lp), B(4,Lp)) ElseIf (Ifun .Eq. 102) Then C Peterson for heavy quark frag. FV = Pf(x, B(3,Lp)) Else Print *, 'Ifun out of range in val. FitIni! Ifun = ', 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 C gluon,Sea KP = -ABS (LP) If (Ifun .Eq. 1) Then FS = Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) ElseIf (Ifun .Eq. 2) Then FS = Fa(x,B(1,KP),B(2,KP))* F2(x,B(3,KP))* B(0,KP) ElseIf (Ifun.eq.3) then C B(n,-1) : dbar - ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then dmu = Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) Fs = dpu / 2.D0 + (-Kp -1.5D0) * dmu else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.4) then if (KP.eq.-1 .or. KP.eq.-2) then dmu = Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) Fs = dpu / 2.D0 + (-Kp -1.5D0) * dmu elseif( KP.eq.0) then FS=Fa(x,B(1,KP),B(2,KP))* F3(x,B(3,KP),B(4,KP))* B(0,KP) else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.5) then C MRS form C Delta = B(0,-4); As=B(0,-3) If ( kp.eq.-1) then FS=Fa(x,B(1,-2),B(2,-2))* F3(x,B(3,-2),B(4,-2))*B(0,-2) > *(1.D0-B(0,-4))*0.2D0 > -.5D0*Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))*B(0,-1) ElseIf ( kp.eq.-2) then FS=Fa(x,B(1,-2),B(2,-2))* F3(x,B(3,-2),B(4,-2))*B(0,-2) > *(1.D0-B(0,-4))*0.2D0 > +.5D0*Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))*B(0,-1) ElseIf ( kp.eq.-3) then FS = Fa(x,B(1,KP),B(2,KP))* F3(x,B(3,KP),B(4,KP))*B(0,KP) > *(1.D0-B(0,-4))*0.1D0 ElseIf ( kp.eq.-4) then FS = Fa(x,B(1,KP),B(2,KP))* F3(x,B(3,KP),B(4,KP))*B(0,KP) > * B(0,-3) * 0.5D0 Elseif( kp.eq.0) then FS = Fa(x,B(1,KP),B(2,KP))* F3(x,B(3,KP),B(4,KP))*B(0,KP) endif ElseIf (Ifun.eq.6) then if (KP.eq.-1 .or. KP.eq.-2) then FS=Fa(x,B(1,-2),B(2,-2))* F6(x,B(3,-2),B(4,-2))* B(0,-2)/2. > + (-KP -1.5D0) > *Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) else FS=Fa(x,B(1,KP),B(2,KP))* F6(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.7) then if (KP.eq.-1 .or. KP.eq.-2) then FS=Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2)/2. > + (-KP -1.5D0) > *Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) elseif (KP.eq.0) then FS=Fa(x,B(1,KP),B(2,KP))* F7(x,B(3,KP),B(4,KP))* B(0,KP) else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.8) then if (KP.eq.-1 .or. KP.eq.-2) then FS=Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2)/2. > + (-KP -1.5D0) > *Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) elseif (KP.eq.0) then FS=Fa(x,B(1,KP),B(2,KP))* F8(x,B(3,KP),B(4,KP))* B(0,KP) else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.9) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0 and C (possibly) constant as x->1 dou = B(0,-1) * Fa(x,B(1,-1),B(2,-1))* F3(x,0D0,B(3,-1)) > + 1D0 + B(4,-1) * x dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) if (Kp.eq.-1) Then Fs = dpu / (dou + 1D0) else Fs = dpu * dou / (dou + 1D0) EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.10) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0 and C (possibly) constant as x->1 dou = B(0,-1) * Fa(x,B(1,-1),B(2,-1)) > + F3(x,B(3,-1), B(4,-1)) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) if (Kp.eq.-1) Then Fs = dpu / (dou + 1D0) else Fs = dpu * dou / (dou + 1D0) EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.11) then if (KP.eq.-1 .or. KP.eq.-2) then FS=Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2)/2. > + (-KP -1.5D0) > *Fa(x,B(1,-1),B(2,-1))* F3(x,B(3,-1),B(4,-1))* B(0,-1) elseif (KP.eq.0) then FS=( Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP)) > + F11(x,B(3,-4),B(4,-4)) ) * B(0,KP) If (FS .lt. 0) then print*,'x,FS:',x,FS FS=0 endif else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.12) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0, C peaks at z=B(1,-1) with width B(2,-1) in Z z = zfrmx(x) dou = 1D0 + B(0,-1) * F12 (z, B(1,-1), B(2,-1)) > * z**3 * (1D0 + B(3,-1)*z**2 + B(4,-1)*z**4) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) ubar = dpu / (dou + 1D0) if (Kp.eq.-1) Then Fs = ubar else Fs = ubar * dou EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.13) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0, C peaks at z=B(1,-1) with width B(2,-1) in Z z = x**(B(3,-1)) dou = 1D0 + B(4,-1)*x + > B(0,-1) * x**0.5 * F12(z,B(1,-1),B(2,-1)) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) ubar = dpu / (dou + 1D0) if (Kp.eq.-1) Then Fs = ubar else Fs = ubar * dou EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.14) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0, C peaks at x=B(1,-1) with width B(2,-1) dou = F3(x,B(3,-1),B(4,-1)) + > B(0,-1) *(F12(x,B(1,-1),B(2,-1))-F12(D0,B(1,-1),B(2,-1))) dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) ubar = dpu / (dou + 1D0) if (Kp.eq.-1) Then Fs = ubar else Fs = ubar * dou EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif ElseIf (Ifun.eq.15) then C B(n,-1) : dbar / ubar C B(n,-2) : dbar + ubar if (KP.eq.-1 .or. KP.eq.-2) then C do(ver)u is chosen to approach 1 at x=0, C and has a broad oscillating shape as per E866 ArgSin = x**2 *B(3,-1)*(1D0 + x *B(4,-1)) If (ArgSin .Le. 0 .or. ArgSin .Ge. 6.2D0) Then dou = 1.0 Else dou = 1D0 > + B(0,-1) * (1D0 + B(1,-1) * x**B(2,-1)) *sin (ArgSin) Endif dpu = Fa(x,B(1,-2),B(2,-2))* F1(x,B(3,-2),B(4,-2))* B(0,-2) ubar = dpu / (dou + 1D0) if (Kp.eq.-1) Then Fs = ubar else Fs = ubar * dou EndIf else FS=Fa(x,B(1,KP),B(2,KP))* F1(x,B(3,KP),B(4,KP))* B(0,KP) endif Else Print *, 'Ifun out of range in sea FitIni! Ifun = ', Ifun Stop Endif SEA = FS FitIni = VL + SEA Return C **************************** Entry FFN(Y) C Integrand for quark number sum rule C Revelant for valence quarks only IF ((Y.LT.SML).OR.(Y.GT.D1M)) THEN FFN = 0. Return ENDIF C Choose functional form if (ifun .eq. 1.or. ifun.eq.3 .or. ifun.eq.4 > .or. ifun.ge.7 .or. ifun.le.15) then Tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 2) then Tem = Fa(Y,A(1),A(2)) * F2(Y,A(3)) elseif (ifun .eq. 5) then Tem = Fa(Y,A(1),A(2)) * F3(Y,A(3),A(4)) elseif (ifun .eq. 6) then Tem = Fa(Y,A(1),A(2)) * F6(Y,A(3),A(4)) ElseIf (Ifun .Eq. 101) Then C Gaussian Tem = Gf(Y, A(3), A(4)) ElseIf (Ifun .Eq. 102) Then C Peterson for heavy quark frag. Tem = Pf(Y, A(3)) else Print *, 'Ifun out of range in FFN! Ifun = ', Ifun Stop endif ffn = tem Return C ------------- momentum integrand ENTRY FFM (Y) C used for valence, sea (incl. dbar+ubar) IF ((Y.LT.SML).OR.(Y.GT.D1M)) THEN FFM = 0. Return ENDIF C Choose functional form if (ifun .eq. 1.or. ifun.eq.3 .or. ifun.eq.4 > .or. ifun.ge.7 .or. ifun.le.15) then Tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 2) then Tem = Fa(Y,A(1),A(2)) * F2(Y,A(3)) elseif (ifun .eq. 5) then Tem = Fa(Y,A(1),A(2)) * F3(Y,A(3),A(4)) elseif (ifun .eq. 6) then Tem = Fa(Y,A(1),A(2)) * F6(Y,A(3),A(4)) else Print *, 'Ifun out of range in FFM! Ifun = ', Ifun Stop endif ffm = tem * Y Return C **************************** ENTRY FFg (Y) C Integrand of gluon momentum integral IF ((Y.LT.SML).OR.(Y.GT.D1M)) THEN FFg = 0. Return ENDIF C Choose functional form if (ifun .eq. 1) then tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 2) then tem = Fa(Y,A(1),A(2)) * F2(Y,A(3)) elseif (ifun .eq. 3) then tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 4) then tem = Fa(Y,A(1),A(2)) * F3(Y,A(3),A(4)) elseif (ifun .eq. 5) then tem = Fa(Y,A(1),A(2)) * F3(Y,A(3),A(4)) elseif (ifun .eq. 6) then Tem = Fa(Y,A(1),A(2)) * F6(Y,A(3),A(4)) elseif (ifun .eq. 7) then Tem = Fa(Y,A(1),A(2)) * F7(Y,A(3),A(4)) elseif (ifun .eq. 8) then Tem = Fa(Y,A(1),A(2)) * F8(Y,A(3),A(4)) elseif (ifun .eq. 9) then tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 10) then tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) elseif (ifun .eq. 11) then Tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) > + F11(Y,B(3,-4),B(4,-4)) If (Tem .lt. 0) then print*,'y,Tem:',y,Tem Tem=0 endif elseif (ifun.ge.12 .or. ifun.le.15) then tem = Fa(Y,A(1),A(2)) * F1(Y,A(3),A(4)) else Print *, 'Ifun out of range in FFg! Ifun = ', Ifun Stop endif ffg = Y * tem Return C ************************** End C 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) csek 2 Lines added INTEGER IKNL Data Iknl/0/ 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) csek 2 Lines added Integer IRR Data Irr/0/ 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) Integer IPRTN Print *, >'You just called a dummy test function PdfTst in PrtPdf!' IPRTN = 0 XX=0. QQ=0. 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) csek Dimension Fq(5), Df(5), QL(5) Dimension Fq(5), Df(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 csek CALL WARNR (IWRN1, NWRT, 'Q less than QINI in PARDIS call; Q SET csek >= QINI.', 'Q', Q, QINI, QMAX, 1) CALL WARNR(IWRN1,NWRT,'Q < QINI in Structm','Q',Q,QINI,QMAX,1) Q = QINI ElseIF (Q .GT. QMAX) THEN csek CALL WARNR(IWRN2,NWRT,'Q greater than QMAX in PARDIS call; csek >Extrapolation will be used.', 'Q', Q, QINI, QMAX, 1) CALL WARNR(IWRN2,NWRT,'Q > QMAX in Structm','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) csek Dimension Fq(5), Df(5), QL(5) Dimension Fq(5), Df(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 csek Line changed CALL WARNR(IWRN1,NWRT,'Q < QINI ','Q',Q,QINI,QMAX,1) Q = QINI ElseIF (Q .GT. QMAX) THEN csek Line changed CALL WARNR(IWRN2,NWRT,'Q > QMAX ','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 / csek DATA NA, IA, PUNY / 3, 5, 1D-30 / DATA PUNY / 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 in Xarray', > '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 / csek 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) csek DATA IRUN / 0 / csek 1 funny line to avoid the Vax compiler neff = neff 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 csek 1 line to avoid the Vax compiler iknl = iknl csek 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) csek DATA IRUN / 0 / csek 1 strange line to avoid the Vax compiler Neff = Neff 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) csek 3 funny lines to avoid the Vax compiler Y(1) = Y(1) TT = TT Nqrk = Nqrk 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 csek 4 more crazy statements to avoid the Vax compiler Qmx(1) = Qmx(1) Xmn(1) = Xmn(1) Nqrk(1) = Nqrk(1) Isch(1) = Isch(1) 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 csek 2 more silly statements to avoid the Vax compiler Mex(1) = Mex(1) Mpn(1) = Mpn(1) 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 / csek Data ist, lp, qsto, Aln2 / 0, -10, 1.2345, 0.6931 / Data ist, lp, qsto / 0, -10, 1.2345 / csek 2 more stupid statements to avoid Vax compiler Mpn(1) = Mpn(1) Mex(1) = Mex(1) 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) csek WRITE (NOUUT, 131) QINI, IPD0, IHDN, QMAX, IKNL, csek > XMIN, XCR, NX, NT, JT, KF csek 131 FORMAT ( / csek >' Current parameters and values are: '// csek >' Initiation parameters: Qini, Ipd0, Ihdn = ', F8.2, 2I8 // csek >' Maximum Q, Order of Alpha: Qmax, IKNL = ', 1PE10.2, I6// csek >' X- mesh parameters : Xmin, Xcr, Nx = ', 2(1PE10.2), I8// csek >' LnQ-mesh parameters : Nt, Jt = ', 2I8 // csek >' # 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 continue 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: ' ) csek READ(NIN, *, ERR=22) V1, V2, V3, V4, V5 READ(NIN, *) 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 csek DATA SML, IRUN, TEM, RER / 1E-5, 0, D1, 1E-3 / DATA TEM, RER / D1, 1E-3 / csek DATA E1, ZLOW, ZHIGH, IWRN1, IWRN2 / 1.00002, -10.0,1.00002, 2*0 / DATA ZLOW, ZHIGH, IWRN2 / -10.0,1.00002, 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 csek DATA IWRN1,IWRN2,IWRN3,IWRN4, HUGE, TINY / 4*0, 1.E35, 1.E-35 / DATA IWRN1, HUGE, TINY / 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 csek DATA LA, LB, TINY / 2 * .FALSE., 1.E-15 / DATA LA, LB / 2 * .FALSE. / 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) Common/W50512/QCDL4,QCDL5 LOGICAL SET IF (.NOT. SET) CALL LAMCWZ csek WRITE (NOUT,110) AL, NF, NORDER csek 110 FORMAT( csek 1 ' Lambda (MSBAR) =',G13.5,', NFL (total # of Flavors) =',I3, csek 2 ', Order (loops) =', I2) csek WRITE (NOUT,120) (I,VALQMS(I),I=1,NF) csek 120 FORMAT (3(' M', I1, '=', G13.5, :, ',')) QCDL4 = Alamf(4) QCDL5 = Alamf(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 csek 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 / csek DATA IW1, IW2, IW3, IW4, IR1, SML / 4*0, 0, 1.E-5 / DATA IR1, SML / 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 csek modified Data Nmax / 3 / IW = IWRN IV = IVAB IF (IW .EQ. 0) THEN csek PRINT '(1X,A/1X, 2A,I10 /A,I4)', MSG, NMVAR, ' = ', IV, csek > ' For all warning messages, check file unit #', NWRT IF (IACT .EQ. 1) THEN csek PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX csek 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 csek Print '(/A/)', '!!! Severe Warning, Too many errors !!!' csek Print '(/A/)', ' !!! Check The Error File !!!' Write (Nwrt, '(/A)') > 'Message suppressed from now on, extrapolation used.' Write (Nwrt,*) > 'But should your code be using this value of the parameter?' 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 / 3 / IW = IWRN VR = VARIAB IF (IW .EQ. 0) THEN csek PRINT '(1X, A/1X,2A,1PD16.7/A,I4)', MSG, NMVAR, ' = ', VR, csek > ' For all warning messages, check file unit #', NWRT IF (IACT .EQ. 1) THEN csek PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX csek 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 csek Print '(/A/)', '!!! Severe Warning, Too many errors !!!' csek Print '(/A/)', ' !!! Check The Error File !!!' Write (Nwrt, '(/A)') > 'Message suppressed from now on, extrapolation used ' Write (Nwrt,*) > 'But should your code be using this value of the parameter?' 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,IRT 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/ csek 2 dumb lines inserted next Err = err Irt = Irt 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 Zbrnt = B 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 csek DATA SML, HUGE, XX / 1.D-5, 1.D30, 0., 1.0, 2.0 / DATA SML, XX / 1.D-5, 0., 1.0, 2.0 / csek DATA IW1, IW2, IW3, IW4, IW5, IW6, IW7 / 7 * 0 / DATA IW1, IW3, IW5 / 3 * 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','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 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