pythiaeRHIC/pythia/pyofsh.f

0001
0002 C***********************************************************************
0003
0004 C...PYOFSH
0005 C...Calculates partial width and differential cross-section maxima
0006 C...of channels/processes not allowed on mass-shell, and selects
0007 C...masses in such channels/processes.
0008
0009       SUBROUTINE PYOFSH(MOFSH,KFMO,KFD1,KFD2,PMMO,RET1,RET2)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Commonblocks.
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0019       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0020       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0021       COMMON/PYINT1/MINT(400),VINT(400)
0022       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0023       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0024       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,/PYINT1/,
0025      &/PYINT2/,/PYINT5/
0026 C...Local arrays.
0027       DIMENSION KFD(2),MBW(2),PMD(2),PGD(2),PMG(2),PML(2),PMU(2),
0028      &PMH(2),ATL(2),ATU(2),ATH(2),RMG(2),INX1(100),XPT1(100),
0029      &FPT1(100),INX2(100),XPT2(100),FPT2(100),WDTP(0:400),
0030      &WDTE(0:400,0:5)
0031
0032 C...Find if particles equal, maximum mass, matrix elements, etc.
0033       MINT(51)=0
0034       ISUB=MINT(1)
0035       KFD(1)=IABS(KFD1)
0036       KFD(2)=IABS(KFD2)
0037       MEQL=0
0038       IF(KFD(1).EQ.KFD(2)) MEQL=1
0039       MLM=0
0040       IF(MOFSH.GE.2.AND.MEQL.EQ.1) MLM=INT(1.5D0+PYR(0))
0041       IF(MOFSH.LE.2.OR.MOFSH.EQ.5) THEN
0042         NOFF=44
0043         PMMX=PMMO
0044       ELSE
0045         NOFF=40
0046         PMMX=VINT(1)
0047         IF(CKIN(2).GT.CKIN(1)) PMMX=MIN(CKIN(2),VINT(1))
0048       ENDIF
0049       MMED=0
0050       IF((KFMO.EQ.25.OR.KFMO.EQ.35.OR.KFMO.EQ.36).AND.MEQL.EQ.1.AND.
0051      &(KFD(1).EQ.23.OR.KFD(1).EQ.24)) MMED=1
0052       IF((KFMO.EQ.32.OR.IABS(KFMO).EQ.34).AND.(KFD(1).EQ.23.OR.
0053      &KFD(1).EQ.24).AND.(KFD(2).EQ.23.OR.KFD(2).EQ.24)) MMED=2
0054       IF((KFMO.EQ.32.OR.IABS(KFMO).EQ.34).AND.(KFD(2).EQ.25.OR.
0055      &KFD(2).EQ.35.OR.KFD(2).EQ.36)) MMED=3
0056       LOOP=1
0057
0058 C...Find where Breit-Wigners are required, else select discrete masses.
0059   100 DO 110 I=1,2
0060         KFCA=PYCOMP(KFD(I))
0061         IF(KFCA.GT.0) THEN
0062           PMD(I)=PMAS(KFCA,1)
0063           PGD(I)=PMAS(KFCA,2)
0064         ELSE
0065           PMD(I)=0D0
0066           PGD(I)=0D0
0067         ENDIF
0068         IF(MSTP(42).LE.0.OR.PGD(I).LT.PARP(41)) THEN
0069           MBW(I)=0
0070           PMG(I)=PMD(I)
0071           RMG(I)=(PMG(I)/PMMX)**2
0072         ELSE
0073           MBW(I)=1
0074         ENDIF
0075   110 CONTINUE
0076
0077 C...Find allowed mass range and Breit-Wigner parameters.
0078       DO 120 I=1,2
0079         IF(MOFSH.EQ.1.AND.LOOP.EQ.1.AND.MBW(I).EQ.1) THEN
0080           PML(I)=PARP(42)
0081           PMU(I)=PMMX-PARP(42)
0082           IF(MBW(3-I).EQ.0) PMU(I)=MIN(PMU(I),PMMX-PMD(3-I))
0083           IF(PMU(I).LT.PML(I)+PARJ(64)) MBW(I)=-1
0084         ELSEIF(MBW(I).EQ.1.AND.MOFSH.NE.5) THEN
0085           ILM=I
0086           IF(MLM.EQ.2) ILM=3-I
0087           PML(I)=MAX(CKIN(NOFF+2*ILM-1),PARP(42))
0088           IF(MBW(3-I).EQ.0) THEN
0089             PMU(I)=PMMX-PMD(3-I)
0090           ELSE
0091             PMU(I)=PMMX-MAX(CKIN(NOFF+5-2*ILM),PARP(42))
0092           ENDIF
0093           IF(CKIN(NOFF+2*ILM).GT.CKIN(NOFF+2*ILM-1)) PMU(I)=
0094      &    MIN(PMU(I),CKIN(NOFF+2*ILM))
0095           IF(I.EQ.MLM) PMU(I)=MIN(PMU(I),0.5D0*PMMX)
0096           IF(MEQL.EQ.0) PMH(I)=MIN(PMU(I),0.5D0*PMMX)
0097           IF(PMU(I).LT.PML(I)+PARJ(64)) MBW(I)=-1
0098           IF(MBW(I).EQ.1) THEN
0099             ATL(I)=ATAN((PML(I)**2-PMD(I)**2)/(PMD(I)*PGD(I)))
0100             ATU(I)=ATAN((PMU(I)**2-PMD(I)**2)/(PMD(I)*PGD(I)))
0101             IF(MEQL.EQ.0) ATH(I)=ATAN((PMH(I)**2-PMD(I)**2)/(PMD(I)*
0102      &      PGD(I)))
0103           ENDIF
0104         ELSEIF(MBW(I).EQ.1.AND.MOFSH.EQ.5) THEN
0105           ILM=I
0106           IF(MLM.EQ.2) ILM=3-I
0107           PML(I)=MAX(CKIN(48+I),PARP(42))
0108           PMU(I)=PMMX-MAX(CKIN(51-I),PARP(42))
0109           IF(MBW(3-I).EQ.0) PMU(I)=MIN(PMU(I),PMMX-PMD(3-I))
0110           IF(I.EQ.MLM) PMU(I)=MIN(PMU(I),0.5D0*PMMX)
0111           IF(MEQL.EQ.0) PMH(I)=MIN(PMU(I),0.5D0*PMMX)
0112           IF(PMU(I).LT.PML(I)+PARJ(64)) MBW(I)=-1
0113           IF(MBW(I).EQ.1) THEN
0114             ATL(I)=ATAN((PML(I)**2-PMD(I)**2)/(PMD(I)*PGD(I)))
0115             ATU(I)=ATAN((PMU(I)**2-PMD(I)**2)/(PMD(I)*PGD(I)))
0116             IF(MEQL.EQ.0) ATH(I)=ATAN((PMH(I)**2-PMD(I)**2)/(PMD(I)*
0117      &      PGD(I)))
0118           ENDIF
0119         ENDIF
0120   120 CONTINUE
0121       IF(MBW(1).LT.0.OR.MBW(2).LT.0.OR.(MBW(1).EQ.0.AND.MBW(2).EQ.0))
0122      &THEN
0123         CALL PYERRM(3,'(PYOFSH:) no allowed decay product masses')
0124         MINT(51)=1
0125         RETURN
0126       ENDIF
0127
0128 C...Calculation of partial width of resonance.
0129       IF(MOFSH.EQ.1) THEN
0130
0131 C..If only one integration, pick that to be the inner.
0132         IF(MBW(1).EQ.0) THEN
0133           PM2=PMD(1)
0134           PMD(1)=PMD(2)
0135           PGD(1)=PGD(2)
0136           PML(1)=PML(2)
0137           PMU(1)=PMU(2)
0138         ELSEIF(MBW(2).EQ.0) THEN
0139           PM2=PMD(2)
0140         ENDIF
0141
0142 C...Start outer loop of integration.
0143         IF(MBW(1).EQ.1.AND.MBW(2).EQ.1) THEN
0144           ATL2=ATAN((PML(2)**2-PMD(2)**2)/(PMD(2)*PGD(2)))
0145           ATU2=ATAN((PMU(2)**2-PMD(2)**2)/(PMD(2)*PGD(2)))
0146           NPT2=1
0147           XPT2(1)=1D0
0148           INX2(1)=0
0149           FMAX2=0D0
0150         ENDIF
0151   130   IF(MBW(1).EQ.1.AND.MBW(2).EQ.1) THEN
0152           PM2S=PMD(2)**2+PMD(2)*PGD(2)*TAN(ATL2+XPT2(NPT2)*(ATU2-ATL2))
0153           PM2=MIN(PMU(2),MAX(PML(2),SQRT(MAX(0D0,PM2S))))
0154         ENDIF
0155         RM2=(PM2/PMMX)**2
0156
0157 C...Start inner loop of integration.
0158         PML1=PML(1)
0159         PMU1=MIN(PMU(1),PMMX-PM2)
0160         IF(MEQL.EQ.1) PMU1=MIN(PMU1,PM2)
0161         ATL1=ATAN((PML1**2-PMD(1)**2)/(PMD(1)*PGD(1)))
0162         ATU1=ATAN((PMU1**2-PMD(1)**2)/(PMD(1)*PGD(1)))
0163         IF(PML1+PARJ(64).GE.PMU1.OR.ATL1+1D-7.GE.ATU1) THEN
0164           FUNC2=0D0
0165           GOTO 180
0166         ENDIF
0167         NPT1=1
0168         XPT1(1)=1D0
0169         INX1(1)=0
0170         FMAX1=0D0
0171   140   PM1S=PMD(1)**2+PMD(1)*PGD(1)*TAN(ATL1+XPT1(NPT1)*(ATU1-ATL1))
0172         PM1=MIN(PMU1,MAX(PML1,SQRT(MAX(0D0,PM1S))))
0173         RM1=(PM1/PMMX)**2
0174
0175 C...Evaluate function value - inner loop.
0176         FUNC1=SQRT(MAX(0D0,(1D0-RM1-RM2)**2-4D0*RM1*RM2))
0177         IF(MMED.EQ.1) FUNC1=FUNC1*((1D0-RM1-RM2)**2+8D0*RM1*RM2)
0178         IF(MMED.EQ.2) FUNC1=FUNC1**3*(1D0+10D0*RM1+10D0*RM2+RM1**2+
0179      &  RM2**2+10D0*RM1*RM2)
0180         IF(FUNC1.GT.FMAX1) FMAX1=FUNC1
0181         FPT1(NPT1)=FUNC1
0182
0183 C...Go to next position in inner loop.
0184         IF(NPT1.EQ.1) THEN
0185           NPT1=NPT1+1
0186           XPT1(NPT1)=0D0
0187           INX1(NPT1)=1
0188           GOTO 140
0189         ELSEIF(NPT1.LE.8) THEN
0190           NPT1=NPT1+1
0191           IF(NPT1.LE.4.OR.NPT1.EQ.6) ISH1=1
0192           ISH1=ISH1+1
0193           XPT1(NPT1)=0.5D0*(XPT1(ISH1)+XPT1(INX1(ISH1)))
0194           INX1(NPT1)=INX1(ISH1)
0195           INX1(ISH1)=NPT1
0196           GOTO 140
0197         ELSEIF(NPT1.LT.100) THEN
0198           ISN1=ISH1
0199   150     ISH1=ISH1+1
0200           IF(ISH1.GT.NPT1) ISH1=2
0201           IF(ISH1.EQ.ISN1) GOTO 160
0202           DFPT1=ABS(FPT1(ISH1)-FPT1(INX1(ISH1)))
0203           IF(DFPT1.LT.PARP(43)*FMAX1) GOTO 150
0204           NPT1=NPT1+1
0205           XPT1(NPT1)=0.5D0*(XPT1(ISH1)+XPT1(INX1(ISH1)))
0206           INX1(NPT1)=INX1(ISH1)
0207           INX1(ISH1)=NPT1
0208           GOTO 140
0209         ENDIF
0210
0211 C...Calculate integral over inner loop.
0212   160   FSUM1=0D0
0213         DO 170 IPT1=2,NPT1
0214           FSUM1=FSUM1+0.5D0*(FPT1(IPT1)+FPT1(INX1(IPT1)))*
0215      &    (XPT1(INX1(IPT1))-XPT1(IPT1))
0216   170   CONTINUE
0217         FUNC2=FSUM1*(ATU1-ATL1)/PARU(1)
0218   180   IF(MBW(1).EQ.1.AND.MBW(2).EQ.1) THEN
0219           IF(FUNC2.GT.FMAX2) FMAX2=FUNC2
0220           FPT2(NPT2)=FUNC2
0221
0222 C...Go to next position in outer loop.
0223           IF(NPT2.EQ.1) THEN
0224             NPT2=NPT2+1
0225             XPT2(NPT2)=0D0
0226             INX2(NPT2)=1
0227             GOTO 130
0228           ELSEIF(NPT2.LE.8) THEN
0229             NPT2=NPT2+1
0230             IF(NPT2.LE.4.OR.NPT2.EQ.6) ISH2=1
0231             ISH2=ISH2+1
0232             XPT2(NPT2)=0.5D0*(XPT2(ISH2)+XPT2(INX2(ISH2)))
0233             INX2(NPT2)=INX2(ISH2)
0234             INX2(ISH2)=NPT2
0235             GOTO 130
0236           ELSEIF(NPT2.LT.100) THEN
0237             ISN2=ISH2
0238   190       ISH2=ISH2+1
0239             IF(ISH2.GT.NPT2) ISH2=2
0240             IF(ISH2.EQ.ISN2) GOTO 200
0241             DFPT2=ABS(FPT2(ISH2)-FPT2(INX2(ISH2)))
0242             IF(DFPT2.LT.PARP(43)*FMAX2) GOTO 190
0243             NPT2=NPT2+1
0244             XPT2(NPT2)=0.5D0*(XPT2(ISH2)+XPT2(INX2(ISH2)))
0245             INX2(NPT2)=INX2(ISH2)
0246             INX2(ISH2)=NPT2
0247             GOTO 130
0248           ENDIF
0249
0250 C...Calculate integral over outer loop.
0251   200     FSUM2=0D0
0252           DO 210 IPT2=2,NPT2
0253             FSUM2=FSUM2+0.5D0*(FPT2(IPT2)+FPT2(INX2(IPT2)))*
0254      &      (XPT2(INX2(IPT2))-XPT2(IPT2))
0255   210     CONTINUE
0256           FSUM2=FSUM2*(ATU2-ATL2)/PARU(1)
0257           IF(MEQL.EQ.1) FSUM2=2D0*FSUM2
0258         ELSE
0259           FSUM2=FUNC2
0260         ENDIF
0261
0262 C...Save result; second integration for user-selected mass range.
0263         IF(LOOP.EQ.1) WIDW=FSUM2
0264         WID2=FSUM2
0265         IF(LOOP.EQ.1.AND.(CKIN(46).GE.CKIN(45).OR.CKIN(48).GE.CKIN(47)
0266      &  .OR.MAX(CKIN(45),CKIN(47)).GE.1.01D0*PARP(42))) THEN
0267           LOOP=2
0268           GOTO 100
0269         ENDIF
0270         RET1=WIDW
0271         RET2=WID2/WIDW
0272
0273 C...Select two decay product masses of a resonance.
0274       ELSEIF(MOFSH.EQ.2.OR.MOFSH.EQ.5) THEN
0275   220   DO 230 I=1,2
0276           IF(MBW(I).EQ.0) GOTO 230
0277           PMBW=PMD(I)**2+PMD(I)*PGD(I)*TAN(ATL(I)+PYR(0)*
0278      &    (ATU(I)-ATL(I)))
0279           PMG(I)=MIN(PMU(I),MAX(PML(I),SQRT(MAX(0D0,PMBW))))
0280           RMG(I)=(PMG(I)/PMMX)**2
0281   230   CONTINUE
0282         IF((MEQL.EQ.1.AND.PMG(MAX(1,MLM)).GT.PMG(MIN(2,3-MLM))).OR.
0283      &  PMG(1)+PMG(2)+PARJ(64).GT.PMMX) GOTO 220
0284
0285 C...Weight with matrix element (if none known, use beta factor).
0286         FLAM=SQRT(MAX(0D0,(1D0-RMG(1)-RMG(2))**2-4D0*RMG(1)*RMG(2)))
0287         IF(MMED.EQ.1) THEN
0288           WTBE=FLAM*((1D0-RMG(1)-RMG(2))**2+8D0*RMG(1)*RMG(2))
0289         ELSEIF(MMED.EQ.2) THEN
0290           WTBE=FLAM**3*(1D0+10D0*RMG(1)+10D0*RMG(2)+RMG(1)**2+
0291      &    RMG(2)**2+10D0*RMG(1)*RMG(2))
0292         ELSEIF(MMED.EQ.3) THEN
0293           WTBE=FLAM*(RMG(1)+FLAM**2/12D0)
0294         ELSE
0295           WTBE=FLAM
0296         ENDIF
0297         IF(WTBE.LT.PYR(0)) GOTO 220
0298         RET1=PMG(1)
0299         RET2=PMG(2)
0300
0301 C...Find suitable set of masses for initialization of 2 -> 2 processes.
0302       ELSEIF(MOFSH.EQ.3) THEN
0303         IF(MBW(1).NE.0.AND.MBW(2).EQ.0) THEN
0304           PMG(1)=MIN(PMD(1),0.5D0*(PML(1)+PMU(1)))
0305           PMG(2)=PMD(2)
0306         ELSEIF(MBW(2).NE.0.AND.MBW(1).EQ.0) THEN
0307           PMG(1)=PMD(1)
0308           PMG(2)=MIN(PMD(2),0.5D0*(PML(2)+PMU(2)))
0309         ELSE
0310           IDIV=-1
0311   240     IDIV=IDIV+1
0312           PMG(1)=MIN(PMD(1),0.1D0*(IDIV*PML(1)+(10-IDIV)*PMU(1)))
0313           PMG(2)=MIN(PMD(2),0.1D0*(IDIV*PML(2)+(10-IDIV)*PMU(2)))
0314           IF(IDIV.LE.9.AND.PMG(1)+PMG(2).GT.0.9D0*PMMX) GOTO 240
0315         ENDIF
0316         RET1=PMG(1)
0317         RET2=PMG(2)
0318
0319 C...Evaluate importance of excluded tails of Breit-Wigners.
0320         IF(MEQL.EQ.0.AND.MBW(1).EQ.1.AND.MBW(2).EQ.1.AND.PMD(1)+PMD(2)
0321      &  .GT.PMMX.AND.PMH(1).GT.PML(1).AND.PMH(2).GT.PML(2)) MEQL=2
0322         IF(MEQL.LE.1) THEN
0323           VINT(80)=1D0
0324           DO 250 I=1,2
0325             IF(MBW(I).NE.0) VINT(80)=VINT(80)*1.25D0*(ATU(I)-ATL(I))/
0326      &      PARU(1)
0327   250     CONTINUE
0328         ELSE
0329           VINT(80)=(1.25D0/PARU(1))**2*MAX((ATU(1)-ATL(1))*
0330      &    (ATH(2)-ATL(2)),(ATH(1)-ATL(1))*(ATU(2)-ATL(2)))
0331         ENDIF
0332         IF((ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.30.OR.ISUB.EQ.35).AND.
0333      &  MSTP(43).NE.2) VINT(80)=2D0*VINT(80)
0334         IF(ISUB.EQ.22.AND.MSTP(43).NE.2) VINT(80)=4D0*VINT(80)
0335         IF(MEQL.GE.1) VINT(80)=2D0*VINT(80)
0336
0337 C...Pick one particle to be the lighter (if improves efficiency).
0338       ELSEIF(MOFSH.EQ.4) THEN
0339         IF(MEQL.EQ.0.AND.MBW(1).EQ.1.AND.MBW(2).EQ.1.AND.PMD(1)+PMD(2)
0340      &  .GT.PMMX.AND.PMH(1).GT.PML(1).AND.PMH(2).GT.PML(2)) MEQL=2
0341   260   IF(MEQL.EQ.2) MLM=INT(1.5D0+PYR(0))
0342
0343 C...Select two masses according to Breit-Wigner + flat in s + 1/s.
0344         DO 270 I=1,2
0345           IF(MBW(I).EQ.0) GOTO 270
0346           PMV=PMU(I)
0347           IF(MEQL.EQ.2.AND.I.EQ.MLM) PMV=PMH(I)
0348           ATV=ATU(I)
0349           IF(MEQL.EQ.2.AND.I.EQ.MLM) ATV=ATH(I)
0350           RBR=PYR(0)
0351           IF((ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.ISUB.EQ.30.OR.
0352      &    ISUB.EQ.35).AND.MSTP(43).NE.2) RBR=2D0*RBR
0353           IF(RBR.LT.0.8D0) THEN
0354             PMSR=PMD(I)**2+PMD(I)*PGD(I)*TAN(ATL(I)+PYR(0)*(ATV-ATL(I)))
0355             PMG(I)=MIN(PMV,MAX(PML(I),SQRT(MAX(0D0,PMSR))))
0356           ELSEIF(RBR.LT.0.9D0) THEN
0357             PMG(I)=SQRT(MAX(0D0,PML(I)**2+PYR(0)*(PMV**2-PML(I)**2)))
0358           ELSEIF(RBR.LT.1.5D0) THEN
0359             PMG(I)=PML(I)*(PMV/PML(I))**PYR(0)
0360           ELSE
0361             PMG(I)=SQRT(MAX(0D0,PML(I)**2*PMV**2/(PML(I)**2+PYR(0)*
0362      &      (PMV**2-PML(I)**2))))
0363           ENDIF
0364   270   CONTINUE
0365         IF((MEQL.GE.1.AND.PMG(MAX(1,MLM)).GT.PMG(MIN(2,3-MLM))).OR.
0366      &  PMG(1)+PMG(2)+PARJ(64).GT.PMMX) THEN
0367           IF(MINT(48).EQ.1.AND.MSTP(171).EQ.0) THEN
0368             NGEN(0,1)=NGEN(0,1)+1
0369             NGEN(MINT(1),1)=NGEN(MINT(1),1)+1
0370             GOTO 260
0371           ELSE
0372             MINT(51)=1
0373             RETURN
0374           ENDIF
0375         ENDIF
0376         RET1=PMG(1)
0377         RET2=PMG(2)
0378
0379 C...Give weight for selected mass distribution.
0380         VINT(80)=1D0
0381         DO 280 I=1,2
0382           IF(MBW(I).EQ.0) GOTO 280
0383           PMV=PMU(I)
0384           IF(MEQL.EQ.2.AND.I.EQ.MLM) PMV=PMH(I)
0385           ATV=ATU(I)
0386           IF(MEQL.EQ.2.AND.I.EQ.MLM) ATV=ATH(I)
0387           F0=PMD(I)*PGD(I)/((PMG(I)**2-PMD(I)**2)**2+
0388      &    (PMD(I)*PGD(I))**2)/PARU(1)
0389           F1=1D0
0390           F2=1D0/PMG(I)**2
0391           F3=1D0/PMG(I)**4
0392           FI0=(ATV-ATL(I))/PARU(1)
0393           FI1=PMV**2-PML(I)**2
0394           FI2=2D0*LOG(PMV/PML(I))
0395           FI3=1D0/PML(I)**2-1D0/PMV**2
0396           IF((ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.ISUB.EQ.30.OR.
0397      &    ISUB.EQ.35).AND.MSTP(43).NE.2) THEN
0398             VINT(80)=VINT(80)*20D0/(8D0+(FI0/F0)*(F1/FI1+6D0*F2/FI2+
0399      &      5D0*F3/FI3))
0400           ELSE
0401             VINT(80)=VINT(80)*10D0/(8D0+(FI0/F0)*(F1/FI1+F2/FI2))
0402           ENDIF
0403           VINT(80)=VINT(80)*FI0
0404   280   CONTINUE
0405         IF(MEQL.GE.1) VINT(80)=2D0*VINT(80)
0406       ENDIF
0407
0408       RETURN
0409       END