This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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