GEANT321/peanut/prepre.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:22:03  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 28/03/94  01.30.59  by  S.Giani
  11 *-- Author :
  12 #if defined(CERNLIB_HPUX)
  13 $OPTIMIZE OFF
  14 #endif
  15 *$ CREATE PREPRE.FOR
  16 *COPY PREPRE
  17 *
  18 *=== prepre ===========================================================*
  19 *
  20       SUBROUTINE PREPRE ( WEE, NNEUT, NPROT, NHOLE, ICYCL )
  21
  22 #include "geant321/dblprc.inc"
  23 #include "geant321/dimpar.inc"
  24 #include "geant321/iounit.inc"
  25 *
  26 *----------------------------------------------------------------------*
  27 *----------------------------------------------------------------------*
  28 *
  29 #include "geant321/balanc.inc"
  30 #include "geant321/eva0.inc"
  31 #include "geant321/fheavy.inc"
  32 #include "geant321/finuc.inc"
  33 #include "geant321/nucdat.inc"
  34 #include "geant321/nucgeo.inc"
  35 #include "geant321/parevt.inc"
  36 #include "geant321/parnuc.inc"
  37 #include "geant321/part.inc"
  38 #include "geant321/resnuc.inc"
  39 *
  40       REAL RNDM(2)
  41       COMMON / FKCOSP / C1ST (3), C2ND (3), LEMISS
  42       LOGICAL LEMISS
  43       COMMON / FKCMCY / NPCYCL (20,2), IEVT, LOUT
  44       COMMON / FKPLOC / IABCOU
  45 *
  46       DIMENSION ERNCM  (2), EPSMX  (2), AMMAFT (2), AMNAFT (2),
  47      &          ZAFTR  (2), EEXMNM (2), EEXDLT (2), EEXANW (2),
  48      &          EPSEEX (2), EPSDLT (2), EPSANY (2), EPSFIX (29),
  49      &          C (3), UMULEG (0:2), ACOLEG (0:2), NPART (2),
  50      &          NPRFIX (2)
  51       LOGICAL LSPREJ, LNUCTS
  52 *
  53       NPART (1) = NNEUT
  54       NPART (2) = NPROT
  55       IF ( LHLFIX ) THEN
  56          NHLFIX = NHLEXP
  57          NHOLE  = NHOLE - NHLFIX
  58          EHLFIX = 0.D+00
  59          DO 50 IHOLE = 1, NHLFIX
  60             EHLFIX = EHLFIX + HOLEXP (IHOLE)
  61    50    CONTINUE
  62          IF ( NHOLE .GT. 0 ) THEN
  63             RHOIMP = ( RHOIMP * ( NHOLE + NHLFIX ) - RHOEXP ) / NHOLE
  64             EKFIMP = ( EKFIMP * ( NHOLE + NHLFIX ) - EKFEXP ) / NHOLE
  65          ELSE
  66             RHOIMP = RHOAVE
  67             EKFIMP = EKFAVE
  68          END IF
  69       ELSE
  70          NHLFIX = 0
  71          EHLFIX = 0.D+00
  72       END IF
  73       NHINI  = NHOLE
  74       IF ( PTRES .LT. ANGLGB ) THEN
  75          UMO2  = ERES * ERES
  76          UMO   = SQRT (UMO2)
  77          EKRES = 0.D+00
  78          GAMCM = 1.D+00
  79          ETAX  = 0.D+00
  80          ETAY  = 0.D+00
  81          ETAZ  = 0.D+00
  82          PXORI = 0.D+00
  83          PYORI = 0.D+00
  84          PZORI = 0.D+00
  85          PCMORI = 0.D+00
  86          CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
  87       ELSE
  88          UMO2  = ( ERES - PTRES ) * ( ERES + PTRES )
  89          UMO   = SQRT (UMO2)
  90          EKRES = ERES - UMO
  91          GAMCM = ERES  / UMO
  92          ETACM = PTRES / UMO
  93          ETAX  = PXRES / UMO
  94          ETAY  = PYRES / UMO
  95          ETAZ  = PZRES / UMO
  96          ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
  97          PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
  98          PCMSX = PXORI - ETAX * PHELP
  99          PCMSY = PYORI - ETAY * PHELP
 100          PCMSZ = PZORI - ETAZ * PHELP
 101          PCMORI = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
 102          CXAXCM = PCMSX / PCMORI
 103          CYAXCM = PCMSY / PCMORI
 104          CZAXCM = PCMSZ / PCMORI
 105       END IF
 106       BBTAR = IBTAR
 107       ZZTAR = ICHTAR
 108       TVCMS  = UMO - AMNRES
 109       EHLFIX = MIN ( TVCMS, EHLFIX )
 110       TVEFF  = TVCMS - EHLFIX
 111       IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
 112       TEMNU = SQRT ( TVCMS / ANOW * ALEVEL  )
 113       NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
 114       NPTOT  = NPART (1) + NPART (2)
 115       IF ( NHOLE + NPTOT + NHLFIX .GE. NEXMX .OR.
 116      &     NPTOT .GT. NINT (0.5D+00 * ANOW) ) GO TO 7000
 117       ANPTOT = NPTOT
 118       ANPROT = NPART (2)
 119       ANNEUT = NPART (1)
 120       AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
 121      &       + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
 122       IAAFT = IBRES-1
 123       IZAFT = ICRES
 124       AAFTR     = IAAFT
 125       ZAFTR (1) = IZAFT
 126       AMMAFT(1) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(1))
 127       AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC + ELBNDE (IZAFT)
 128       CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
 129       BNDGAV = BNENRG (2)
 130       EPSMX  (1) = UMO - AMNRES - BNDGAV - EHLFIX
 131       EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
 132       EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
 133       IZAFT = ICRES - 1
 134       ZAFTR (2) = IZAFT
 135       AMMAFT(2) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(2))
 136       AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC + ELBNDE (IZAFT)
 137       CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
 138       BNDGAV = BNENRG (1)
 139       EPSMX  (2) = UMO - AMNRES - BNDGAV - EHLFIX
 140       IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU ) GO TO 7000
 141       EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
 142       EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
 143       IF ( NP .LE. NP0 ) THEN
 144          IF ( KPORI .EQ. 8 ) THEN
 145             EEXANW (1) = EEXDLT (1)
 146             EPSANY (1) = MIN ( EPSDLT (1), EPSMX (1) )
 147             EEXANW (2) = 0.D+00
 148             EPSANY (2) = EPSMX  (2)
 149          ELSE IF ( KPORI .EQ. 1 ) THEN
 150             EEXANW (1) = 0.D+00
 151             EPSANY (1) = EPSMX  (1)
 152             EEXANW (2) = EEXDLT (2)
 153             EPSANY (2) = MIN ( EPSDLT (2), EPSMX (2) )
 154          ELSE
 155             EEXANW (1) = 0.D+00
 156             EPSANY (1) = EPSMX (1)
 157             EEXANW (2) = 0.D+00
 158             EPSANY (2) = EPSMX (2)
 159          END IF
 160       ELSE
 161          EEXANW (1) = 0.D+00
 162          EPSANY (1) = EPSMX (1)
 163          EEXANW (2) = 0.D+00
 164          EPSANY (2) = EPSMX (2)
 165       END IF
 166       SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
 167       NEMISS = 0
 168       MNUCTS = 0
 169  1000 CONTINUE
 170          NPRFIX (1) = 0
 171          NPRFIX (2) = 0
 172          JEMFIX = 0
 173          EPRFIX = 0.D+00
 174          IF ( .NOT. LEMISS ) ICYCL = ICYCL + 1
 175          IF ( NPTOT .LE. 0 ) GO TO 4600
 176          IF ( NNUCTS .GT. MNUCTS ) THEN
 177             MNUCTS = MNUCTS + 1
 178             JNUCTS = INUCTS (MNUCTS)
 179             JEMIS0 = 2 - ABS (KPNUCL(JNUCTS)) / 8
 180             JEMIS1 = JEMIS0
 181             JDEMIS = 1
 182             EKFSAV = EKFIMP
 183             RHOSAV = RHOIMP
 184             TMPRHO = 0.001D+00 * RHOAVE
 185             RHOIMP = MAX ( RHNUCL (JNUCTS), TMPRHO )
 186             TMPEKF = 0.001D+00 * EKFAVE
 187             EKFIMP = MAX ( EKFNUC (JNUCTS), TMPEKF )
 188             LNUCTS = .TRUE.
 189          ELSE
 190             LNUCTS = .FALSE.
 191             CALL GRNDM(RNDM,1)
 192             IF ( RNDM (1) .LT. ANNEUT / ANPTOT ) THEN
 193                JEMIS0 = 1
 194                JEMIS1 = 2
 195                JDEMIS = 1
 196             ELSE
 197                JEMIS0 = 2
 198                JEMIS1 = 1
 199                JDEMIS = -1
 200             END IF
 201          END IF
 202          NEXTOT = NPTOT + NHOLE
 203          RHORED = ACOLL / BBTAR
 204          IF ( LGDHPR ) THEN
 205             EKFCNN = ( ZNOW * EKFCEN (1) + ( ANOW - ZNOW ) * EKFCEN (2)
 206      &             ) / ANOW
 207             RHOWEI = 0.5D+00
 208             IF ( LNUCTS ) THEN
 209                AEFREF = EKFIMP
 210                AEFRAV = EKFIMP
 211                RHONUC = RHOIMP
 212             ELSE IF ( NPTOT .LE. 2 .AND. ICYCL .EQ. 1 ) THEN
 213                AHLTOT = MAX ( NHOLE, 1 )
 214                WEIGH1 = NHINI / AHLTOT
 215                WEIGH2 = 1.D+00 - WEIGH1
 216                AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
 217                RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
 218                AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
 219                RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
 220      &                * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
 221                AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
 222      &                / ( AHLTOT + 1.D+00 )
 223             ELSE
 224                AHLTOT = MAX ( NHOLE, 1 )
 225                WEIGH1 = NHINI / AHLTOT
 226                WEIGH2 = 1.D+00 - WEIGH1
 227                AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
 228                RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
 229                AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
 230                RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
 231      &                * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
 232                AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
 233      &                / ( AHLTOT + 1.D+00 )
 234             END IF
 235             CLAMDI = 1.0D+00
 236          ELSE
 237             EKFCNN = AEFRMX
 238             AEFREF = 0.5D+00 * AEFRMX
 239             AEFRAV = AEFREF
 240             RHONUC = 0.5D+00**1.5D+00 * RHONU0
 241             CLAMDI = 0.5D+00
 242          END IF
 243          DO 4500 JEMISS = JEMIS0, JEMIS1, JDEMIS
 244             IF ( NPART (JEMISS) .LE. 0 ) GO TO 4500
 245             IEMISS = 3 - JEMISS
 246             RATEC0 = 2.D+00 * AMNUCL (IEMISS) / PLABRC**3 / PISQ
 247             BNDGEN = AMNAFT (JEMISS) + AMNUCL (IEMISS) - AMNRES
 248             IF ( .NOT. LNUCTS .AND. EPRFIX .GT. ANGLGB ) THEN
 249                TVEFF  = TVEFF - EPRFIX
 250                NPTOT  = NPTOT - NPRFIX (1) - NPRFIX (2)
 251                NPART (1) = NPART (1) - NPRFIX (1)
 252                NPART (2) = NPART (2) - NPRFIX (2)
 253                JEMFIX = JEMISS
 254                ERNCM0 = ERNCM (JEMISS)
 255                EPSMX0 = EPSMX (JEMISS)
 256                BNDGAV = BNENRG (IEMISS)
 257                EPSMX (JEMISS) = UMO - AMNRES - BNDGAV - EHLFIX - EPRFIX
 258             ELSE
 259                JEMFIX = 0
 260             END IF
 261             BNDGAV = BNENRG (IEMISS)
 262             EPSMIN = BNDGEN - BNDGAV
 263             BNDHLP = BNENRG (IEMISS)
 264             IF ( EPSMX (JEMISS) - EPSMIN .LT. TEMNU ) GO TO 4500
 265             IF ( EHLFIX + EPRFIX .LE. 0.5D+00 * EEXANW (JEMISS) ) THEN
 266                UUMIN  = TVEFF - 0.5D+00 * ( EPSANY (JEMISS)
 267      &                + EPSMX (JEMISS) ) - BNDGAV
 268             ELSE
 269                UUMIN  = TVEFF - EPSMX (JEMISS) - BNDGAV
 270             END IF
 271             UUMAX  = TVEFF - BNDGEN
 272             BNDUSE = BNENRG (IEMISS)
 273             IF ( LNUCTS ) THEN
 274                EKNNUC = ENNUC  (JNUCTS) - AMNUCL (IEMISS) + BNDGEN
 275      &                - RSTNUC (JNUCTS)
 276                EPSHLP = BNDHLP - BNDGEN + EPSMIN
 277                EKNNUC = EKNNUC + EPSHLP
 278                IF ( EKNNUC .LE. EPSMIN ) THEN
 279                   GO TO 4500
 280                ELSE IF ( EKNNUC .LE. 0.D+00 ) THEN
 281                   DELTAE = 0.5D+00 * ( - EPSMIN - EKNNUC )
 282                   EKNNUC = EKNNUC + DELTAE
 283                ELSE
 284                   DELTAE = 0.D+00
 285                END IF
 286                ENNUC (JNUCTS) = EKNNUC + AMNUCL (IEMISS)
 287                PNUCCO = SQRT ( EKNNUC * ( EKNNUC + 2.D+00
 288      &                * AMNUCL (IEMISS) ) )
 289                PHELP = PNUCCO / PNUCL (JNUCTS)
 290                PCMSX = PXNUCL (JNUCTS) * PHELP
 291                PCMSY = PYNUCL (JNUCTS) * PHELP
 292                PCMSZ = PZNUCL (JNUCTS) * PHELP
 293                ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
 294                PHELP = ENNUC (JNUCTS) - ETAPCM / ( GAMCM + 1.D+00 )
 295                PCMSX = PCMSX - ETAX * PHELP
 296                PCMSY = PCMSY - ETAY * PHELP
 297                PCMSZ = PCMSZ - ETAZ * PHELP
 298                PCMST = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
 299                EPS = GAMCM * ENNUC (JNUCTS) - ETAPCM - AMNUCL (IEMISS)
 300                IF ( EPS .LE. EPSMIN ) THEN
 301                   GO TO 4500
 302                ELSE IF ( EPS .LT. ANGLGB * ENNUC (JNUCTS) ) THEN
 303                   GO TO 4500
 304                END IF
 305                EPS   = EPS - DELTAE
 306                C (1) = PCMSX / PCMST
 307                C (2) = PCMSY / PCMST
 308                C (3) = PCMSZ / PCMST
 309                EPS = MIN ( EPS, HLFHLF * ( EPSMX (JEMISS)
 310      &             + EPSANY (JEMISS) ) )
 311                IF ( JEMISS .EQ. 2 ) THEN
 312                   FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
 313                   CCOUL  = DOST ( 3, ZAFTR (JEMISS) )
 314                   ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
 315      &                   / RMASS (IBRES-1)
 316                END IF
 317                FLAMMX = 1.D+00
 318                NTENT  = 1
 319                ISAMPL = 0
 320                LSPREJ = .FALSE.
 321                GO TO 3000
 322             END IF
 323             IF ( UUMIN .GE. UUMAX ) GO TO 4500
 324             EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
 325      &             + AMNAFT (JEMISS) )**2 ) / AMNAFT (JEMISS)
 326             EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
 327             BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
 328             EKEWLL = EPSWLL - AMNUCL (IEMISS)
 329             IF ( JEMISS .EQ. 1 ) THEN
 330                AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
 331                BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
 332      &                / AALPHA * 1.D-03
 333                TMP07  = 0.07D+00
 334                EPSEFF = MIN ( EPSINV, TMP07 )
 335                SIGINV = 0.1D+00 * XINNEU ( EPSEFF, ZAFTR (JEMISS),
 336      &                  BBETAA )
 337                SFTYMX = 1.1D+00
 338                GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) / EKFCEN (IEMISS)
 339                IF ( EKEWLL .LE. 0.04D+00 ) THEN
 340                   SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
 341      &                   -1.86D+00 + 0.09415D+03 * EKEWLL
 342      &                   + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
 343      &                   * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
 344      &                   + 0.13D+03 * EKEWLL )**2 )
 345                   IF ( EKEWLL .LT. 0.02D+00 ) THEN
 346                      SIGMAN = 0.3333333333333333D+00 * SIGMAP
 347                   ELSE
 348                      SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 349      &                      + 42.9D+00
 350                   END IF
 351                ELSE
 352                   SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
 353      &                   + 82.2D+00
 354                   SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 355      &                   + 42.9D+00
 356                END IF
 357                SIGMAP = 0.1D+00 * SIGMAP
 358                SIGMAN = 0.1D+00 * SIGMAN
 359                SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
 360      &                + ( ZNOW - ANPROT ) * SIGMAP )
 361      &                / ( ANOW - ANPTOT )
 362             ELSE
 363                EPSCOU = UMO - AMNAFT (JEMISS) - AMNUCL (IEMISS)
 364                FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
 365                CCOUL  = DOST ( 3, ZAFTR (JEMISS) )
 366                ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
 367      &                / RMASS (IBRES-1)
 368                TMP07  = 0.07D+00
 369                EPSEFF = MIN ( EPSINV, TMP07 )
 370                SIGINV = 1.D-01 * XINPRO ( EPSEFF, ZAFTR (JEMISS),
 371      &                                    ETHRES )
 372                SFTYMX = 1.2D+00
 373                GSNGLV = 1.5D+00 * ZNOW / EKFCEN (IEMISS)
 374                IF ( EKEWLL .LE. 0.04D+00 ) THEN
 375                   SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
 376      &                   -1.86D+00 + 0.09415D+03 * EKEWLL
 377      &                   + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
 378      &                   * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
 379      &                   + 0.13D+03 * EKEWLL )**2 )
 380                   IF ( EKEWLL .LT. 0.02D+00 ) THEN
 381                      SIGMAP = 0.3333333333333333D+00 * SIGMAN
 382                   ELSE
 383                      SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 384      &                      + 42.9D+00
 385                   END IF
 386                ELSE
 387                   SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
 388      &                   + 82.2D+00
 389                   SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 390      &                   + 42.9D+00
 391                END IF
 392                SIGMAP = 0.1D+00 * SIGMAP
 393                SIGMAN = 0.1D+00 * SIGMAN
 394                SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
 395      &                + ( ZNOW - ANPROT ) * SIGMAP )
 396      &                / ( ANOW - ANPTOT )
 397             END IF
 398             ZITA   = AEFRAV / EKEWLL
 399             IF ( ZITA .LE. 0.5D+00 ) THEN
 400                PZITA = 1.D+00 - 1.4D+00 * ZITA
 401             ELSE
 402                PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
 403      &               * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
 404             END IF
 405             RATEC  = SIGINV * RATEC0 * EPSINV / GSNGLV
 406             IF ( RATEC .LE. 0.D+00 ) GO TO 4500
 407             ALAMDC = BETWLL / RATEC
 408             ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (TWOTWO
 409      &             * PI * PLABRC * BETWLL * EPSWLL / AMNUSQ (IEMISS),
 410      &             RZNUCL )
 411             AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
 412             RATEI  = BETWLL / ALAMDI
 413             FLAMMX = SFTYMX * RATEC / ( RATEC + CLAMDI * RATEI )
 414             IF ( FLAMMX .LT. 0.D+00 ) THEN
 415                GO TO 4500
 416             END IF
 417             LSPREJ = .NOT. ( NHOLE .GT. 2 .OR. NPTOT .GT. 4 )
 418             IF ( LSPREJ ) THEN
 419                IF ( NPTOT .EQ. 1 ) THEN
 420                   IF ( NHOLE .EQ. 0 ) THEN
 421                      EPS = EPSMX (JEMISS)
 422                      FLAMMX = 1.D+00
 423                      NTENT  = 1
 424                      ISAMPL = 0
 425                      LSPREJ = .FALSE.
 426                      GO TO 3000
 427                   ELSE IF ( NHOLE .EQ. 1 ) THEN
 428                      ISAMPL = 1
 429                   ELSE IF ( NHOLE .EQ. 2 ) THEN
 430                      ISAMPL = 2
 431                   ELSE
 432                      ISAMPL = 7
 433                   END IF
 434                ELSE IF ( NPTOT .EQ. 2 ) THEN
 435                   IF ( NHOLE .EQ. 0 ) THEN
 436                      ISAMPL = 7
 437                   ELSE IF ( NHOLE .EQ. 1 ) THEN
 438                      ISAMPL = 3
 439                   ELSE IF ( NHOLE .EQ. 2 ) THEN
 440                      ISAMPL = 4
 441                   ELSE
 442                      ISAMPL = 7
 443                   END IF
 444                ELSE IF ( NPTOT .EQ. 3 ) THEN
 445                   IF ( NHOLE .EQ. 0 ) THEN
 446                      ISAMPL = 7
 447                   ELSE IF ( NHOLE .EQ. 1 ) THEN
 448                      ISAMPL = 5
 449                   ELSE IF ( NHOLE .EQ. 2 ) THEN
 450                      ISAMPL = 6
 451                   ELSE
 452                      ISAMPL = 7
 453                   END IF
 454                ELSE
 455                   ISAMPL = 7
 456                END IF
 457             ELSE
 458                ISAMPL = 7
 459             END IF
 460             GO TO ( 2100, 2200, 2300, 2400, 2500, 2600, 2700 ), ISAMPL
 461  2100       CONTINUE
 462                LSPREJ = .FALSE.
 463                AEFRPA = AEFREF
 464                UUMAX  = MIN ( AEFRPA, UUMAX )
 465                AITOT  = ( AEFRPA**3 - ( AEFRPA + UUMIN - UUMAX )**3 )
 466      &                / AEFRPA**3
 467                IF ( UUMAX .GT. TVEFF ) THEN
 468                   AIPRO = MIN ( AITOT, ONEONE ) * FLAMMX
 469                   AITOT = AITOT * FLAMMX
 470                ELSE
 471                   AITOT = AITOT * FLAMMX
 472                   AIPRO = AITOT
 473                END IF
 474             GO TO 2800
 475  2200       CONTINUE
 476                LSPREJ = .FALSE.
 477                AEFRPA = AEFREF
 478                UUMAX  = MIN ( TWOTWO * AEFRPA, UUMAX )
 479                IF ( TVEFF .LE. AEFRPA ) THEN
 480                   EDENOM = TVEFF**2 * 0.5D+00
 481                ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
 482                   EDENOM = ( TVEFF**2 - 2.D+00 * ( TVEFF - AEFRPA )**2 )
 483      &                   * 0.5D+00
 484                ELSE
 485                   EDENOM = AEFRPA**2
 486                END IF
 487                UUMN2 = UUMIN**2
 488                IF ( UUMAX .GT. AEFRPA ) THEN
 489                   UUDIV = AEFRPA
 490                   UUDV2 = UUDIV**2
 491                   AIPR1 = 0.5D+00 * UUDV2
 492                   AIPR2 = 0.5D+00 * ( UUMAX - UUDIV ) * ( 3.D+00 * UUDIV
 493      &                  - UUMAX )
 494                   IF ( UUMIN .LE. UUDIV ) THEN
 495                      AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
 496                      AITO2 = AIPR2
 497                      AIHLP = 0.D+00
 498                   ELSE
 499                      AITO1 = 0.D+00
 500                      AIHLP = 0.5D+00 * ( UUMIN - UUDIV ) * ( 3.D+00
 501      &                     * UUDIV - UUMIN )
 502                      AITO2 = AIPR2 - AIHLP
 503                   END IF
 504                ELSE
 505                   UUDIV = UUMAX
 506                   UUDV2 = UUDIV**2
 507                   AIPR1 = 0.5D+00 * UUDV2
 508                   AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
 509                   AITO2 = 0.D+00
 510                   AIPR2 = 0.D+00
 511                END IF
 512                AITOT = AITO1 + AITO2
 513                AITO1 = AITO1 / AITOT
 514                AITO2 = AITO2 / AITOT
 515                AITOT = AITOT * FLAMMX / EDENOM
 516                IF ( UUMAX .GT. TVEFF ) THEN
 517                   AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM, ONEONE )
 518      &                  * FLAMMX
 519                ELSE
 520                   AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
 521                END IF
 522             GO TO 2800
 523  2300       CONTINUE
 524                LSPREJ = .FALSE.
 525                AEFRPA = AEFREF
 526                IF ( TVEFF .LE. AEFRPA ) THEN
 527                   EDENOM = 0.25D+00 * TVEFF * TVEFF
 528                ELSE
 529                   EDENOM = 0.25D+00 * AEFRPA * ( 2.D+00 * TVEFF
 530      &                   - AEFRPA )
 531                END IF
 532                FCHLP = 0.25D+00 * NPART (JEMISS)
 533                IF ( UUMAX .LE. AEFRPA ) THEN
 534                   UUDIV = UUMAX
 535                   AIPR1 = FCHLP * UUDIV**2
 536                   AITO1 = AIPR1 - FCHLP * UUMIN**2
 537                   AITO2 = 0.D+00
 538                   AIPR2 = 0.D+00
 539                ELSE
 540                   UUDIV = AEFRPA
 541                   AIPR1 = FCHLP * UUDIV**2
 542                   AIPR2 = FCHLP * ( UUMAX - UUDIV ) * UUDIV
 543                   IF ( UUMIN .GT. UUDIV ) THEN
 544                      AITO1 = 0.D+00
 545                      AIHLP = FCHLP * ( UUMIN - UUDIV ) * UUDIV
 546                      AITO2 = AIPR2 - AIHLP
 547                   ELSE
 548                      AIHLP = 0.D+00
 549                      AITO1 = AIPR1 - FCHLP * UUMIN**2
 550                      AITO2 = AIPR2
 551                   END IF
 552                END IF
 553                AITOT = AITO1 + AITO2
 554                AITO1 = AITO1 / AITOT
 555                AITO2 = AITO2 / AITOT
 556                IF ( UUMAX .GT. TVEFF ) THEN
 557                   DDNPAR = NPART(JEMISS)
 558                   AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
 559      &                    DDNPAR ) * FLAMMX
 560                ELSE
 561                   AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
 562                END IF
 563             GO TO 2800
 564  2400       CONTINUE
 565                LSPREJ = .TRUE.
 566                AEFRPA = AEFREF
 567                IF ( TVEFF .LE. AEFRPA ) THEN
 568                   LSPREJ = .FALSE.
 569                   FEMAX  = 1.D+00
 570                   ISAMPL = 7
 571                ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
 572                   FEMAX = TVEFF**3 / ( TVEFF**3 - 2.D+00 * ( TVEFF
 573      &                  - AEFRPA )**3 )
 574                ELSE
 575                   FEMAX = 0.1666666666666667D+00 * TVEFF**3
 576      &                  / ( AEFRPA**2 * ( TVEFF - AEFRPA ) )
 577                END IF
 578                UUMNR  = ( UUMIN / TVEFF )**(NEXTOT-1)
 579                UUMXR  = ( UUMAX / TVEFF )**(NEXTOT-1)
 580                AITOT  = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
 581      &                * FEMAX
 582                AIPRO  = NPART (JEMISS) * FLAMMX * UUMXR * FEMAX
 583             GO TO 2800
 584  2500       CONTINUE
 585                LSPREJ = .FALSE.
 586                AEFRPA = AEFREF
 587                IF ( TVEFF .LE. AEFRPA ) THEN
 588                   EDENOM = TVEFF**3 / 36.D+00
 589                ELSE
 590                   EDENOM = ( TVEFF**3 - ( TVEFF - AEFRPA )**3 )
 591      &                   / 26.D+00
 592                END IF
 593                IF ( UUMAX .GT. AEFRPA ) THEN
 594                   UUDIV = AEFRPA
 595                   UUDV3 = UUDIV**3
 596                   UUMN3 = UUMIN**3
 597                   FCHLP = NPART(JEMISS) / 36.D+00
 598                   AIPR1 = FCHLP * UUDV3
 599                   AIPR2 = 3.D+00 * FCHLP * UUMAX * UUDIV * ( UUMAX
 600      &                  - UUDIV )
 601                   IF ( UUMIN .LE. UUDIV ) THEN
 602                      AITO1 = AIPR1 - FCHLP * UUMN3
 603                      AITO2 = AIPR2
 604                      AIHLP = 0.D+00
 605                   ELSE
 606                      AITO1 = 0.D+00
 607                      AIHLP = 3.D+00 * FCHLP * UUMIN * UUDIV * ( UUMIN
 608      &                     - UUDIV )
 609                      AITO2 = AIPR2 - AIHLP
 610                   END IF
 611                ELSE
 612                   UUDIV = UUMAX
 613                   UUDV3 = UUDIV**3
 614                   UUMN3 = UUMIN**3
 615                   FCHLP = NPART(JEMISS) / 36.D+00
 616                   AIPR1 = FCHLP * UUDV3
 617                   AITO1 = AIPR1 - FCHLP * UUMN3
 618                   AIPR2 = 0.D+00
 619                   AITO2 = 0.D+00
 620                END IF
 621                AITOT = AITO1 + AITO2
 622                AITO1 = AITO1 / AITOT
 623                AITO2 = AITO2 / AITOT
 624                AITOT = AITOT * FLAMMX / EDENOM
 625                IF ( UUMAX .GT. TVEFF ) THEN
 626                   DDNPAR = NPART(JEMISS)
 627                   AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
 628      &                    DDNPAR ) * FLAMMX
 629                ELSE
 630                   AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
 631                END IF
 632             GO TO 2800
 633  2600       CONTINUE
 634                AEFRPA = AEFREF
 635                IF ( TVEFF .LE. AEFRPA ) THEN
 636                   EDENOM = TVEFF**4 / 288.D+00
 637                ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
 638                   EDENOM = ( TVEFF**4 - 2.D+00 * ( TVEFF - AEFRPA )**4 )
 639      &                   / 288.D+00
 640                ELSE
 641                   EDENOM = AEFRPA**2 * ( 0.5D+00 * TVEFF**2 - TVEFF
 642      &                   * AEFRPA + 7.D+00 / 12.D+00 * AEFRPA**2 )
 643      &                   / 12.D+00
 644                END IF
 645                FCHLP = NPART (JEMISS) / 288.D+00
 646                IF ( UUMAX .GT. 2.D+00 * AEFRPA ) THEN
 647                   LSPREJ = .TRUE.
 648                   UUDIV = AEFRPA
 649                   UUDI2 = 2.D+00 * AEFRPA
 650                   UUDV4 = UUDIV**4
 651                   UUMN4 = UUMIN**4
 652                   UUD24 = UUDI2**4
 653                   AIPR1 = FCHLP * UUDV4
 654                   AIPR2 = FCHLP * ( UUD24 - UUDV4 )
 655                   AIPR3 = 12.D+00 * FCHLP * AEFRPA**2 * UUMAX
 656      &                  * ( UUMAX - UUDI2 )
 657                   IF ( UUMIN .LE. UUDIV ) THEN
 658                      AITO1 = AIPR1 - FCHLP * UUMN4
 659                      AITO2 = AIPR2
 660                      AITO3 = AIPR3
 661                      AIHLP = 0.D+00
 662                      AIHL2 = 0.D+00
 663                   ELSE IF ( UUMIN .LE. UUDI2 ) THEN
 664                      AITO1 = 0.D+00
 665                      AIHLP = FCHLP * ( UUMN4 - UUDV4 )
 666                      AITO2 = AIPR2 - AIHLP
 667                      AITO3 = AIPR3
 668                      AIHL2 = 0.D+00
 669                   ELSE
 670                      AITO1 = 0.D+00
 671                      AITO2 = 0.D+00
 672                      AIHL2 = 12.D+00 * FCHLP * AEFRPA**2 * UUMIN
 673      &                     * ( UUMIN - UUDI2 )
 674                      AITO3 = AIPR3 - AIHL2
 675                   END IF
 676                ELSE IF ( UUMAX .GT. AEFRPA ) THEN
 677                   LSPREJ = .TRUE.
 678                   UUDIV = AEFRPA
 679                   UUDI2 = UUMAX
 680                   UUDV4 = UUDIV**4
 681                   UUMN4 = UUMIN**4
 682                   UUD24 = UUDI2**4
 683                   AIPR1 = FCHLP * UUDV4
 684                   AIPR2 = FCHLP * ( UUD24 - UUDV4 )
 685                   IF ( UUMIN .LE. UUDIV ) THEN
 686                      AITO1 = AIPR1 - FCHLP * UUMN4
 687                      AITO2 = AIPR2
 688                      AIHLP = 0.D+00
 689                   ELSE
 690                      AITO1 = 0.D+00
 691                      AIHLP = FCHLP * ( UUMN4 - UUDV4 )
 692                      AITO2 = AIPR2 - AIHLP
 693                   END IF
 694                   AIPR3 = 0.D+00
 695                   AITO3 = 0.D+00
 696                ELSE
 697                   LSPREJ = .FALSE.
 698                   UUDI2 = UUMAX
 699                   UUDIV = UUMAX
 700                   UUDV4 = UUDIV**4
 701                   UUMN4 = UUMIN**4
 702                   AIPR1 = FCHLP * UUDV4
 703                   AITO1 = AIPR1 - FCHLP * UUMN4
 704                   AITO2 = 0.D+00
 705                   AIPR2 = 0.D+00
 706                   AITO3 = 0.D+00
 707                   AIPR3 = 0.D+00
 708                END IF
 709                AITOT = AITO1 + AITO2 + AITO3
 710                AITO1 = AITO1 / AITOT
 711                AITO2 = AITO2 / AITOT
 712                AITO3 = AITO3 / AITOT
 713                AITOT = AITOT * FLAMMX / EDENOM
 714                AIPRO = ( AIPR1 + AIPR2 + AIPR3 ) * FLAMMX / EDENOM
 715             GO TO 2800
 716  2700       CONTINUE
 717                LSPREJ = .FALSE.
 718                TVLEFF = LOG(TVEFF)
 719                IF(UUMAX.LE.0.) THEN
 720                   UULMXR=-100
 721                ELSE
 722                   UULMXR = (LOG(UUMAX)-TVLEFF)*(NEXTOT-1)
 723                ENDIF
 724                IF(UULMXR.LT.-88) THEN
 725                   UUMNR = 0.
 726                   UUMXR = 0.
 727                ELSE
 728                   UUMXR = EXP(UULMXR)
 729                   IF(UUMIN.LE.0.) THEN
 730                      UULMNR=-100
 731                   ELSE
 732                      UULMNR = (LOG(UUMIN)-TVLEFF)*(NEXTOT-1)
 733                   ENDIF
 734                   IF(UULMNR.LT.-88) THEN
 735                      UUMNR = 0.
 736                   ELSE
 737                      UUMNR = EXP(UULMNR)
 738                   ENDIF
 739                ENDIF
 740                AITOT  = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
 741                AIPRO  = NPART (JEMISS) * FLAMMX * MIN ( UUMXR, ONEONE )
 742  2800       CONTINUE
 743             NTENT  = INT  (AIPRO)
 744             CALL GRNDM(RNDM,1)
 745             RNTENT = RNDM (1)
 746             IF ( RNTENT .LT. AIPRO - NTENT ) NTENT = NTENT + 1
 747  3000       CONTINUE
 748             ITACC = 0
 749             DO 4100 IT = 1, NTENT
 750                GO TO ( 3100, 3200, 3300, 3400, 3500, 3600, 3700 ),
 751      &            ISAMPL
 752                GO TO 4000
 753  3100          CONTINUE
 754                   CALL GRNDM(RNDM,1)
 755                   UURND = AEFRPA * ( 1.D+00 - AITOT / FLAMMX
 756      &                  * RNDM (1) )**0.3333333333333333D+00
 757                   UURND = AEFRPA - UURND + UUMIN
 758                GO TO 3800
 759  3200          CONTINUE
 760                   CALL GRNDM(RNDM,1)
 761                   RNDEPS = RNDM (1)
 762                   IF ( RNDEPS .LT. AITO1 ) THEN
 763                      RNDEPS = RNDEPS / AITO1
 764                      UURND  = SQRT ( RNDEPS * ( UUDV2 - UUMN2 ) + UUMN2)
 765                   ELSE
 766                      RNDEPS = RNDEPS - AITO1
 767                      RNDEPS = AIHLP + RNDEPS * AITOT * EDENOM / FLAMMX
 768                      UURND  = 2.D+00 * UUDIV - SQRT ( UUDIV**2
 769      &                      - 2.D+00 * RNDEPS )
 770                   END IF
 771                GO TO 3800
 772  3300          CONTINUE
 773                   CALL GRNDM(RNDM,1)
 774                   RNDEPS = RNDM (1)
 775                   IF ( RNDEPS .LT. AITO1 ) THEN
 776                      RNDEPS = RNDEPS * AITOT
 777                      UURND  = RNDEPS / FCHLP + UUMIN**2
 778                      UURND  = SQRT (UURND)
 779                   ELSE
 780                      RNDEPS = ( RNDEPS - AITO1 ) * AITOT + AIHLP
 781                      UURND  = RNDEPS / FCHLP / UUDIV + UUDIV
 782                   END IF
 783                GO TO 3800
 784  3400          CONTINUE
 785                   CALL GRNDM(RNDM,1)
 786                   UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
 787      &                  / FEMAX / NPART(JEMISS) + UUMNR )**
 788      &                  ( 1.D+00 / (NEXTOT-1) )
 789                   IF ( UURND .LE. AEFRPA ) THEN
 790                      LSPREJ = .FALSE.
 791                      FREJE  = 1.D+00
 792                   ELSE IF ( UURND .LE. 2.D+00 * AEFRPA ) THEN
 793                      FREJE  = ( UURND**2 - 2.D+00 * ( UURND - AEFRPA )
 794      &                      **2 ) / UURND**2
 795                   ELSE
 796                      FREJE  = 2.D+00 * AEFRPA**2 / UURND**2
 797                   END IF
 798                GO TO 3800
 799  3500          CONTINUE
 800                   CALL GRNDM(RNDM,1)
 801                   RNDEPS = RNDM (1)
 802                   IF ( RNDEPS .LT. AITO1 ) THEN
 803                      RNDEPS = RNDEPS / AITO1
 804                      UURND  = RNDEPS * ( UUDV3 - UUMN3 ) + UUMN3
 805                      UURND  = UURND**0.333333333333333D+00
 806                      LSPREJ = .FALSE.
 807                   ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
 808                      RNDEPS = RNDEPS - AITO1
 809                      RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
 810                      UURND  = 0.5D+00 * ( UUDIV + SQRT ( UUDIV**2
 811      &                      + 1.333333333333333D+00 * RNDEPS / FCHLP
 812      &                      / UUDIV ) )
 813                      LSPREJ = .FALSE.
 814                   END IF
 815                GO TO 3800
 816  3600          CONTINUE
 817                   CALL GRNDM(RNDM,1)
 818                   RNDEPS = RNDM (1)
 819                   IF ( RNDEPS .LT. AITO1 ) THEN
 820                      RNDEPS = RNDEPS / AITO1
 821                      UURND  = RNDEPS * ( UUDV4 - UUMN4 ) + UUMN4
 822                      UURND  = UURND**0.25D+00
 823                      LSPREJ = .FALSE.
 824                   ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
 825                      RNDEPS = RNDEPS - AITO1
 826                      RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
 827                      UURND  = RNDEPS / FCHLP + UUDV4
 828                      UURND  = UURND**0.25D+00
 829                      LSPREJ = .TRUE.
 830                      FREJE  = 1.D+00 - 2.D+00 * ( 1.D+00 - UUDIV / UURND
 831      &                      )**3
 832                   ELSE
 833                      RNDEPS = RNDEPS - AITO1 - AITO2
 834                      RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHL2
 835                      UURND  = AEFRPA + SQRT ( AEFRPA**2 + RNDEPS /
 836      &                      ( 12.D+00 * FCHLP * AEFRPA**2 ) )
 837                      LSPREJ = .FALSE.
 838                   END IF
 839                GO TO 3800
 840  3700          CONTINUE
 841                   CALL GRNDM(RNDM,1)
 842                   UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
 843      &                  / NPART(JEMISS) + UUMNR )**
 844      &                  ( 1.D+00 / (NEXTOT-1) )
 845  3800          CONTINUE
 846                IF ( LSPREJ ) THEN
 847                   CALL GRNDM(RNDM,1)
 848                   IF ( RNDM (1) .GE. FREJE ) GO TO 4100
 849                END IF
 850                EPS    = TVEFF - UURND - BNDGAV
 851  4000          CONTINUE
 852                IF ( LPRFIX ) THEN
 853                   ITACC = ITACC + 1
 854                   EPSFIX (ITACC) = EPS - EPSMIN + BNDGEN
 855                END IF
 856                DEPSEX = EPS - EPSEEX (JEMISS)
 857                IF ( DEPSEX .GT. 0.D+00 ) THEN
 858                   IF ( EEXDLT (JEMISS) .LT. EEXMNM (JEMISS) .AND.
 859      &               DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS) - EPSEEX
 860      &              (JEMISS) ) ) THEN
 861                      DEPSEX = EPS - EPSDLT (JEMISS)
 862                      IF ( DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS)
 863      &                  - EPSANY (JEMISS) ) ) THEN
 864                         EPS = EPSANY (JEMISS)
 865                      ELSE
 866                         EPS = EPSDLT (JEMISS)
 867                      END IF
 868                   ELSE
 869                      EPS = EPSEEX (JEMISS)
 870                   END IF
 871                   IF ( EPS .LE. EPSMIN ) GO TO 4100
 872                END IF
 873                AMNHLP = UMO - EPS + EPSMIN - AMNUCL (IEMISS)
 874                ERNCM (JEMISS) = 0.5D+00 * ( UMO2 + AMNHLP**2
 875      &                        - AMNUSQ (IEMISS) ) / UMO
 876                EPS    = UMO - ERNCM (JEMISS) - AMNUCL (IEMISS)
 877                EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
 878      &                + AMNHLP )**2 ) / AMNHLP
 879                EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
 880                BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
 881                EKEWLL = EPSWLL - AMNUCL (IEMISS)
 882                EPSCOU = UMO - AMNHLP - AMNUCL (IEMISS)
 883                IF ( JEMISS .EQ. 1 ) THEN
 884                   GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) * SQRT ( ( EKEWLL
 885      &                   + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
 886      &                   / EKFCEN (IEMISS)
 887                   AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
 888                   BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
 889      &                / AALPHA * 1.D-03
 890                   SIGINV = 0.1D+00 * XINNEU ( EPSINV, ZAFTR (JEMISS),
 891      &                                        BBETAA )
 892                   IF ( EKEWLL .LE. 0.04D+00 ) THEN
 893                      SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
 894      &                      -1.86D+00 + 0.09415D+03 * EKEWLL
 895      &                      + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
 896      &                      * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
 897      &                      + 0.13D+03 * EKEWLL )**2 )
 898                      IF ( EKEWLL .LT. 0.02D+00 ) THEN
 899                         SIGMAN = 0.3333333333333333D+00 * SIGMAP
 900                      ELSE
 901                         SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00
 902      &                         / BETWLL + 42.9D+00
 903                      END IF
 904                   ELSE
 905                      SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
 906      &                      + 82.2D+00
 907                      SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 908      &                      + 42.9D+00
 909                   END IF
 910                   SIGMAP = 0.1D+00 * SIGMAP
 911                   SIGMAN = 0.1D+00 * SIGMAN
 912                   SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
 913      &                   + ( ZNOW - ANPROT ) * SIGMAP )
 914      &                   / ( ANOW - ANPTOT )
 915                ELSE
 916                   GSNGLV = 1.5D+00 * ZNOW * SQRT ( ( EKEWLL
 917      &                   + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
 918      &                   / EKFCEN (IEMISS)
 919                   SIGINV = 1.D-01 * XINPRO ( EPSINV, ZAFTR (JEMISS),
 920      &                                       ETHRES )
 921                   IF ( EKEWLL .LE. 0.04D+00 ) THEN
 922                      SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
 923      &                      -1.86D+00 + 0.09415D+03 * EKEWLL
 924      &                      + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
 925      &                      * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
 926      &                      + 0.13D+03 * EKEWLL )**2 )
 927                      IF ( EKEWLL .LT. 0.02D+00 ) THEN
 928                         SIGMAP = 0.3333333333333333D+00 * SIGMAN
 929                      ELSE
 930                         SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00
 931      &                         / BETWLL + 42.9D+00
 932                      END IF
 933                   ELSE
 934                      SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
 935      &                      + 82.2D+00
 936                      SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
 937      &                      + 42.9D+00
 938                   END IF
 939                   SIGMAP = 0.1D+00 * SIGMAP
 940                   SIGMAN = 0.1D+00 * SIGMAN
 941                   SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
 942      &                   + ( ZNOW - ANPROT ) * SIGMAP )
 943      &                   / ( ANOW - ANPTOT )
 944                END IF
 945                ZITA   = AEFRAV / EKEWLL
 946                IF ( ZITA .LE. 0.5D+00 ) THEN
 947                   PZITA = 1.D+00 - 1.4D+00 * ZITA
 948                ELSE
 949                   PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
 950      &                  * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
 951                END IF
 952                RATEC  = SIGINV * RATEC0 * EPSINV / GSNGLV
 953                IF ( RATEC .LE. 0.D+00 ) GO TO 4100
 954                ALAMDC = BETWLL / RATEC
 955                ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (
 956      &                  TWOTWO * PI * PLABRC * BETWLL * EPSWLL
 957      &                / AMNUSQ (IEMISS), RZNUCL )
 958                AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
 959                RATEI  = BETWLL / ALAMDI
 960                FLAMDA = RATEC / ( RATEC + CLAMDI * RATEI ) / FLAMMX
 961                CALL GRNDM(RNDM,1)
 962                IF ( RNDM (1) .LT. FLAMDA ) GO TO 4200
 963  4100       CONTINUE
 964             IF ( JEMISS .EQ. JEMIS1 ) GO TO 4500
 965             IF ( ITACC .LE. NPART (JEMISS) ) THEN
 966                NPRFIX (JEMISS) = NPRFIX (JEMISS) + ITACC
 967                DO 4150 IT = 1, ITACC
 968                   EPRFIX = EPRFIX + EPSFIX (IT)
 969  4150          CONTINUE
 970             ELSE
 971                NPRFIX (JEMISS) = NPART (JEMISS)
 972                ACCEP = 0.D+00
 973                DO 4160 IT = 1, ITACC
 974                   PRACC  = ( NPART (JEMISS) - ACCEP )
 975      &                   / ( ITACC - IT + 1 )
 976                   CALL GRNDM(RNDM,1)
 977                   IF ( RNDM (1) .LT. PRACC ) THEN
 978                      EPRFIX = EPRFIX + EPSFIX (IT)
 979                      ACCEP  = ACCEP  + 1.D+00
 980                   END IF
 981  4160          CONTINUE
 982             END IF
 983             GO TO 4500
 984  4200       CONTINUE
 985             ZNOW = ZAFTR (JEMISS)
 986             ANOW = AAFTR
 987             IBRES = IBRES - 1
 988             ICRES = ICRES - JEMISS + 1
 989             AMMRES = AMMAFT (JEMISS)
 990             AMNRES = AMNAFT (JEMISS)
 991             EEPCM  = EPS + AMNUCL (IEMISS)
 992             ERESCM = UMO - EEPCM
 993             IF ( LNUCTS ) THEN
 994             ELSE IF ( ICYCL .EQ. 1 .AND. LEMISS ) THEN
 995                C (1) = C2ND (1)
 996                C (2) = C2ND (2)
 997                C (3) = C2ND (3)
 998             ELSE IF ( ICYCL .LE. 1 .AND. PCMORI .LE. ANGLGB ) THEN
 999                C (1) = C1ST (1)
1000                C (2) = C1ST (2)
1001                C (3) = C1ST (3)
1002             ELSE IF ( ICYCL .LE. 1 ) THEN
1003                PCMSXT = - AMNUCL (IEMISS) * ETAX
1004                PCMSYT = - AMNUCL (IEMISS) * ETAY
1005                PCMSZT = - AMNUCL (IEMISS) * ETAZ
1006                APSHAR = MAX ( ANPTOT + NEMISS - 1, ONEONE )
1007                PCMSX  = ( CXAXCM * PCMORI + PCMSXT ) / APSHAR - PCMSXT
1008                PCMSY  = ( CYAXCM * PCMORI + PCMSYT ) / APSHAR - PCMSYT
1009                PCMSZ  = ( CZAXCM * PCMORI + PCMSZT ) / APSHAR - PCMSZT
1010 *
1011                PCMSTT = PCMSXT**2 + PCMSYT**2 + PCMSZT**2
1012                PCMSPT = PCMSX**2  + PCMSY**2  + PCMSZ**2
1013                ECMSTT = SQRT ( AMNUSQ (IEMISS) + PCMSTT )
1014                ECMSPT = SQRT ( AMNUSQ (IEMISS) + PCMSPT )
1015 *
1016                UMMO2  = 2.D+00 * AMNUSQ (IEMISS) + 2.D+00 * ECMSTT
1017      &                * ECMSPT - 2.D+00 * ( PCMSXT * PCMSX + PCMSYT
1018      &                * PCMSY + PCMSZT * PCMSZ )
1019                UMMO   = SQRT ( UMMO2 )
1020                GAMCCM = ( ECMSTT + ECMSPT ) / UMMO
1021                ETAXM  = ( PCMSX + PCMSXT ) / UMMO
1022                ETAYM  = ( PCMSY + PCMSYT ) / UMMO
1023                ETAZM  = ( PCMSZ + PCMSZT ) / UMMO
1024                ETAPCM = ETAXM * PCMSX + ETAYM * PCMSY + ETAZM * PCMSZ
1025                PHELP  = ECMSPT - ETAPCM / ( GAMCCM + 1.D+00 )
1026                PCCMSX = PCMSX - ETAXM * PHELP
1027                PCCMSY = PCMSY - ETAYM * PHELP
1028                PCCMSZ = PCMSZ - ETAZM * PHELP
1029                PCCMSP = SQRT ( PCCMSX**2 + PCCMSY**2 + PCCMSZ**2 )
1030                CALL RACO (C(1),C(2),C(3))
1031                SCAPRO = PCCMSX * C (1) + PCCMSY * C (2) + PCCMSZ * C (3)
1032                IF ( SCAPRO .LT. 0.D+00 ) THEN
1033                   C (1) = - C (1)
1034                   C (2) = - C (2)
1035                   C (3) = - C (3)
1036                END IF
1037                PCCMSX = C(1) * PCCMSP
1038                PCCMSY = C(2) * PCCMSP
1039                PCCMSZ = C(3) * PCCMSP
1040                ETAPCM = ETAXM * PCCMSX + ETAYM * PCCMSY + ETAZM * PCCMSZ
1041                PHELP  = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
1042                PCMSX  = PCCMSX + ETAXM * PHELP
1043                PCMSY  = PCCMSY + ETAYM * PHELP
1044                PCMSZ  = PCCMSZ + ETAZM * PHELP
1045                PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
1046                C (1)  = PCMSX / PCMSPT
1047                C (2)  = PCMSY / PCMSPT
1048                C (3)  = PCMSZ / PCMSPT
1049                ETAPCM = - ETAPCM
1050                PHELP  = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
1051                PCMSX  = - PCCMSX + ETAXM * PHELP
1052                PCMSY  = - PCCMSY + ETAYM * PHELP
1053                PCMSZ  = - PCCMSZ + ETAZM * PHELP
1054                PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
1055                C2ND (1) = PCMSX / PCMSPT
1056                C2ND (2) = PCMSY / PCMSPT
1057                C2ND (3) = PCMSZ / PCMSPT
1058             ELSE
1059                UMULEG (0) = 1.0D+00
1060                UMULEG (1) = 0.6666666666666667D+00
1061                UMULEG (2) = 0.25D+00
1062                ACOLEG (0) = 0.5D+00
1063                ACOLEG (1) = 1.5D+00 * UMULEG (1)**ICYCL
1064                ACOLEG (2) = 2.5D+00 * UMULEG (2)**ICYCL
1065                COSTHE = COSLEG ( ACOLEG )
1066                SINTHE = SQRT ( ( 1.D+00 - COSTHE ) * ( 1.D+00 + COSTHE )
1067      &                )
1068  4300          CONTINUE
1069                   CALL GRNDM(RNDM,2)
1070                   RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
1071                   RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
1072                   RPHI12 = RPHI1 * RPHI1
1073                   RPHI22 = RPHI2 * RPHI2
1074                   RSQ = RPHI12 + RPHI22
1075                IF ( RSQ .GT. 1.D+00 ) GO TO 4300
1076                SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
1077                COSPHI = ( RPHI12 - RPHI22 ) / RSQ
1078                SINT02 = ( 1.D+00 - CZAXCM ) * ( 1.D+00 + CZAXCM )
1079                IF ( SINT02 .LT. ANGLSQ ) THEN
1080                   C (1) = COSPHI * SINTHE
1081                   C (2) = SINPHI * SINTHE
1082                   C (3) = CZAXCM * COSTHE
1083                ELSE
1084                   SINTH0 = SQRT ( SINT02 )
1085                   UPRIME = SINTHE * COSPHI
1086                   VPRIME = SINTHE * SINPHI
1087                   COSPH0 = CXAXCM / SINTH0
1088                   SINPH0 = CYAXCM / SINTH0
1089                   C (1) = UPRIME * COSPH0 * CZAXCM - VPRIME * SINPH0
1090      &                  + COSTHE * CXAXCM
1091                   C (2) = UPRIME * SINPH0 * CZAXCM + VPRIME * COSPH0
1092      &                  + COSTHE * CYAXCM
1093                   C (3) = - UPRIME * SINTH0 + COSTHE * CZAXCM
1094                END IF
1095             END IF
1096             PCMS  = SQRT ( ( EEPCM - AMNUCL (IEMISS) ) * ( EEPCM
1097      &            + AMNUCL (IEMISS) ) )
1098             UMO2   = ( ERESCM - PCMS ) * ( ERESCM + PCMS )
1099             UMO    = SQRT (UMO2)
1100             PCMSX = PCMS * C (1)
1101             PCMSY = PCMS * C (2)
1102             PCMSZ = PCMS * C (3)
1103             ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
1104             NP = NP + 1
1105             KPART (NP) = 1 - 7 * ( JEMISS - 2 )
1106             TKI   (NP) = GAMCM * EEPCM + ETAPCM - AMNUCL (IEMISS)
1107             PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEPCM
1108             WEI   (NP) = WEE
1109             PLBPX = PCMSX + ETAX * PHELP
1110             PLBPY = PCMSY + ETAY * PHELP
1111             PLBPZ = PCMSZ + ETAZ * PHELP
1112             PHELP = SQRT ( PLBPX * PLBPX + PLBPY * PLBPY
1113      &            + PLBPZ * PLBPZ )
1114             CXR   (NP) = PLBPX / PHELP
1115             CYR   (NP) = PLBPY / PHELP
1116             CZR   (NP) = PLBPZ / PHELP
1117             PLR   (NP) = PHELP
1118             IGREYP = IGREYP + JEMISS - 1
1119             IGREYN = IGREYN + 2 - JEMISS
1120             IF ( JEMISS .EQ. 1 ) THEN
1121                ISTORE = IGREYN
1122             ELSE
1123                ISTORE = IGREYP
1124             END IF
1125             IF ( LNUCTS ) THEN
1126                NPCYCL (ISTORE,JEMISS) = -ICYCL
1127             ELSE
1128                NPCYCL (ISTORE,JEMISS) = ICYCL
1129             END IF
1130             EINTR  = EINTR  + TKI (NP) + AMNUCL (IEMISS)
1131             PXINTR = PXINTR + CXR (NP) * PLR (NP)
1132             PYINTR = PYINTR + CYR (NP) * PLR (NP)
1133             PZINTR = PZINTR + CZR (NP) * PLR (NP)
1134             IBINTR = IBINTR + IBAR ( KPART (NP) )
1135             ICINTR = ICINTR + ICH  ( KPART (NP) )
1136             ERES  = GAMCM * ERESCM - ETAPCM
1137             EKRES = ERES - UMO
1138             PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERESCM
1139             PXRES = - PCMSX + ETAX * PHELP
1140             PYRES = - PCMSY + ETAY * PHELP
1141             PZRES = - PCMSZ + ETAZ * PHELP
1142             PTRES2= PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
1143             PTRES = SQRT (PTRES2)
1144             TVCMS = UMO - AMNRES
1145             IF ( TVCMS .LT. 0.D+00 ) THEN
1146                IF ( TVCMS .LT. - GAMMIN ) THEN
1147                   WRITE (LUNOUT,*)' PREPRE: TVCMS',TVCMS
1148                   WRITE (LUNERR,*)' PREPRE: TVCMS',TVCMS
1149                   WRITE (77,*)' ^^^ PREPRE: TVCMS',TVCMS
1150                END IF
1151                TVCMS = 0.D+00
1152             END IF
1153             EHLFIX = MIN ( TVCMS, EHLFIX )
1154             TVEFF = TVCMS - EHLFIX
1155             IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
1156             GAMCM = ERES  / UMO
1157             ETACM = PTRES / UMO
1158             ETAX  = PXRES / UMO
1159             ETAY  = PYRES / UMO
1160             ETAZ  = PZRES / UMO
1161             ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
1162             PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
1163             PCMSX = PXORI - ETAX * PHELP
1164             PCMSY = PYORI - ETAY * PHELP
1165             PCMSZ = PZORI - ETAZ * PHELP
1166             PCMORI = PCMSX**2 + PCMSY**2 + PCMSZ**2
1167             IF ( PCMORI .LE. ANGLGB ) THEN
1168                PCMORI = 0.D+00
1169                CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
1170             ELSE
1171                PCMORI = SQRT ( PCMORI )
1172                CXAXCM = PCMSX / PCMORI
1173                CYAXCM = PCMSY / PCMORI
1174                CZAXCM = PCMSZ / PCMORI
1175             END IF
1176             IAAFT = IBRES-1
1177             IZAFT = ICRES
1178             AAFTR = IAAFT
1179             ZAFTR (1) = IZAFT
1180             AMMAFT(1) = AAFTR * AMUAMU + 1.D-03
1181      &                * FKENER (AAFTR,ZAFTR(1))
1182             AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC
1183      &                + ELBNDE (IZAFT)
1184             CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
1185             BNDGAV = BNENRG (2)
1186             EPSMX  (1) = UMO - AMNRES - BNDGAV - EHLFIX
1187             EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
1188             EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
1189             IZAFT = ICRES - 1
1190             ZAFTR (2) = IZAFT
1191             AMMAFT(2) = AAFTR * AMUAMU + 1.D-03
1192      &                * FKENER (AAFTR,ZAFTR(2))
1193             AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC
1194      &                + ELBNDE (IZAFT)
1195             CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
1196             BNDGAV = BNENRG (1)
1197             EPSMX  (2) = UMO - AMNRES - BNDGAV - EHLFIX
1198             EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
1199             EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
1200             EEXANW (1) = 0.D+00
1201             EPSANY (1) = EPSMX (1)
1202             EEXANW (2) = 0.D+00
1203             EPSANY (2) = EPSMX (2)
1204             NPART (JEMISS) = NPART (JEMISS) - 1
1205             IF ( JEMFIX .GT. 0 ) THEN
1206                NPART (1) = NPART (1) + NPRFIX (1)
1207                NPART (2) = NPART (2) + NPRFIX (2)
1208             END IF
1209             NEMISS = NEMISS + 1
1210             LEMISS = .TRUE.
1211             TEMNU = SQRT ( TVCMS / ANOW * ALEVEL  )
1212             IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU )
1213      &         GO TO 7000
1214             SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
1215             NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
1216             ANPROT =  NPART (2)
1217             ANNEUT =  NPART (1)
1218             NPTOT  = NPART (1) + NPART (2)
1219             ANPTOT = NPTOT
1220             GO TO 6000
1221  4500    CONTINUE
1222  4600    CONTINUE
1223          IF ( LNUCTS ) THEN
1224             LEMISS = .TRUE.
1225             GO TO 6000
1226          ELSE
1227             LEMISS = .FALSE.
1228             NPART (1) = NPART (1) + NPRFIX (1)
1229             NPART (2) = NPART (2) + NPRFIX (2)
1230             IF ( JEMFIX .GT. 0 ) THEN
1231                TVEFF = TVCMS - EHLFIX
1232                EPSMX (JEMFIX) = EPSMX0
1233                ERNCM (JEMFIX) = ERNCM0
1234             END IF
1235          END IF
1236  5000    CONTINUE
1237          ANPROT =  NPART (2)
1238          ANNEUT =  NPART (1)
1239          IF ( NPTOT .GT. 0 ) THEN
1240             PNPROT = ( ZNOW - ANPROT ) * ( 3.D+00 * ANNEUT + ANPROT )
1241      &             / ( ANPROT * ( ZNOW - ANPROT + 3.D+00 * ( ANOW
1242      &             - ANNEUT - ZNOW ) ) + ANNEUT * ( 3.D+00 * ( ZNOW
1243      &             - ANPROT ) + ANOW - ANNEUT - ZNOW ) )
1244          ELSE
1245             PNPROT = ZNOW / ANOW
1246          END IF
1247          CALL GRNDM(RNDM,1)
1248          IF ( RNDM (1) .LT. PNPROT ) THEN
1249             NPART (2) = NPART (2) + 1
1250             ZCOLL = ZCOLL - 1.D+00
1251          ELSE
1252             NPART (1) = NPART (1) + 1
1253          END IF
1254          NPTOT  = NPART (1) + NPART (2)
1255          ANPTOT = NPTOT
1256          NHOLE  = NHOLE + 1
1257          ACOLL  = ACOLL - 1.D+00
1258          AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
1259      &          + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
1260  6000    CONTINUE
1261          IF ( LNUCTS ) THEN
1262             RHOIMP = RHOSAV
1263             EKFIMP = EKFSAV
1264          END IF
1265       IF ( NHOLE + NHLFIX + NPTOT .LT. NEXMX .AND.
1266      &     NPTOT .LE. NINT (0.5D+00 * ANOW) ) GO TO 1000
1267  7000 CONTINUE
1268 *=== End of subroutine prepre =========================================*
1269       RETURN
1270       END
1271 #if defined(CERNLIB_HPUX)
1272 $OPTIMIZE ON
1273 #endif