
      PROGRAM MAIN
 
C-----------------------------------------------------------------------
C  MAIN PROGRAM
C
C  SIMULATION OF EXTENSIVE AIR SHOWERS
C  PREPARES INITIALIZATIONS
C  GENERATES SHOWERS IN THE SHOWER LOOP
C  TREATES PARTICLES IN THE PARTICLE LOOP
C  PERFORMS PRINTING OF TABLES AT END OF SHOWER AND AT END OF RUN
C-----------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      parameter (xct=1)
      parameter (yct=2)
      parameter (zct=3)
      parameter (ctthet=4)
      parameter (ctphi=5)
      parameter (ctdiam=6)
      parameter (ctfoc=7)
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

*KEEP,BAL.
      COMMON /BAL/     EBAL
      DOUBLE PRECISION EBAL(10)
*KEEP,BUFFS.
      COMMON /BUFFS/   RUNH,RUNE,EVTH,EVTE,DATAB,LH
      INTEGER          MAXBUF,MAXLEN
      PARAMETER        (MAXBUF=39*7)
      PARAMETER        (MAXLEN=12)
      REAL             RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
     *                 RUNE(MAXBUF),DATAB(MAXBUF)
      INTEGER          LH
      CHARACTER*4      CRUNH,CRUNE,CEVTH,CEVTE
      EQUIVALENCE      (RUNH(1),CRUNH), (RUNE(1),CRUNE)
      EQUIVALENCE      (EVTH(1),CEVTH), (EVTE(1),CEVTE)
*KEEP,CHISTA.
      COMMON /CHISTA/  IHYCHI,IKACHI,IMUCHI,INNCHI,INUCHI,IPICHI
      INTEGER          IHYCHI(124),IKACHI(124),IMUCHI(124),
     *                 INNCHI(124),INUCHI(124),IPICHI(124)
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,CURVE.
      COMMON /CURVE/   CHAPAR,DEP,ERR,NSTP
      DOUBLE PRECISION CHAPAR(1100),DEP(1100),ERR(1100)
      INTEGER          NSTP
*KEEP,ELADPM.
      COMMON /ELADPM/  ELMEAN,ELMEAA,IELDPM,IELDPA
      DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
      INTEGER          IELDPM(37,13),IELDPA(37,13)
*KEEP,ELASTY.
      COMMON /ELASTY/  ELAST,IELIS,IELHM,IELNU,IELPI
      DOUBLE PRECISION ELAST
      INTEGER          IELIS(20),IELHM(20),IELNU(20),IELPI(20)
*KEEP,GENER.
      COMMON /GENER/   GEN,ALEVEL
      DOUBLE PRECISION GEN,ALEVEL
*KEEP,IRET.
      COMMON /IRET/    IRET1,IRET2
      INTEGER          IRET1,IRET2
*KEEP,ISTA.
      COMMON /ISTA/    IFINET,IFINNU,IFINKA,IFINPI,IFINHY
      INTEGER          IFINET,IFINNU,IFINKA,IFINPI,IFINHY
*KEEP,LONGI.
      COMMON /LONGI/   APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
     *                 NSTEP,LLONGI,FLGFIT
      DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
     *                 SPLONG(0:1040,9),THSTEP,THSTPI
      INTEGER          NSTEP
      LOGICAL          LLONGI,FLGFIT
*KEEP,MPARTI.
      COMMON /MPARTI/  MPARTO
      DOUBLE PRECISION MPARTO(10,25),MPHOTO(10),MPOSIT(10),MELECT(10),
     *                 MNU(10),MMUP(10),MMUM(10),MPI0(10),MPIP(10),
     *                 MPIM(10),MK0L(10),MKPL(10),MKMI(10),MNEUTR(10),
     *                 MPROTO(10),MPROTB(10),MK0S(10),MHYP(10),
     *                 MNEUTB(10),MDEUT(10),MTRIT(10),MALPHA(10),
     *                 MOTHER(10)
      EQUIVALENCE (MPARTO(1, 1),MPHOTO(1)), (MPARTO(1, 2),MPOSIT(1)),
     *            (MPARTO(1, 3),MELECT(1)), (MPARTO(1, 4),MNU(1))   ,
     *            (MPARTO(1, 5),MMUP(1))  , (MPARTO(1, 6),MMUM(1))  ,
     *            (MPARTO(1, 7),MPI0(1))  , (MPARTO(1, 8),MPIP(1))  ,
     *            (MPARTO(1, 9),MPIM(1))  , (MPARTO(1,10),MK0L(1))  ,
     *            (MPARTO(1,11),MKPL(1))  , (MPARTO(1,12),MKMI(1))  ,
     *            (MPARTO(1,13),MNEUTR(1)), (MPARTO(1,14),MPROTO(1)),
     *            (MPARTO(1,15),MPROTB(1)), (MPARTO(1,16),MK0S(1))  ,
     *            (MPARTO(1,18),MHYP(1))  , (MPARTO(1,19),MDEUT(1)) ,
     *            (MPARTO(1,20),MTRIT(1)) , (MPARTO(1,21),MALPHA(1)),
     *            (MPARTO(1,22),MOTHER(1)), (MPARTO(1,25),MNEUTB(1))
*KEEP,MULT.
      COMMON /MULT/    EKINL,MSMM,MULTMA,MULTOT
      DOUBLE PRECISION EKINL
      INTEGER          MSMM,MULTMA(37,13),MULTOT(37,13)
*KEEP,MUPART.
      COMMON /MUPART/  AMUPAR,BCUT,CMUON,FMUBRM,FMUORG
      DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11)
      LOGICAL          FMUBRM,FMUORG
*KEEP,NCOUNT.
      COMMON /NCOUNT/  NCOUN
      INTEGER          NCOUN(8)
*KEEP,NKGI.
      COMMON /NKGI/    SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST,
     *                 DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX,
     *                 OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT
      DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10),
     *                 ZSL(10),DIST(10),
     *                 DISX(-10:10),DISY(-10:10),
     *                 DISXY(-10:10,2),DISYX(-10:10,2),
     *                 DLAX (-10:10,2),DLAY (-10:10,2),
     *                 DLAXY(-10:10,2),DLAYX(-10:10,2),
     *                 OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10)
      INTEGER          IALT(2)
*KEEP,NKGS.
      COMMON /NKGS/    CZX,CZY,CZXY,CZYX,SAH,SL,ZNE
      DOUBLE PRECISION CZX(-10:10,2),CZY(-10:10,2),CZXY(-10:10,2),
     *                 CZYX(-10:10,2),SAH(10),SL(10),ZNE(10)
*KEEP,NPARTI.
      COMMON /NPARTI/  NPARTO,MUOND
      DOUBLE PRECISION NPARTO(10,25),NPHOTO(10),NPOSIT(10),NELECT(10),
     *                 NNU(10),NMUP(10),NMUM(10),NPI0(10),NPIP(10),
     *                 NPIM(10),NK0L(10),NKPL(10),NKMI(10),NNEUTR(10),
     *                 NPROTO(10),NPROTB(10),NK0S(10),NHYP(10),
     *                 NNEUTB(10),NDEUT(10),NTRIT(10),NALPHA(10),
     *                 NOTHER(10),MUOND
      EQUIVALENCE (NPARTO(1, 1),NPHOTO(1)), (NPARTO(1, 2),NPOSIT(1)),
     *            (NPARTO(1, 3),NELECT(1)), (NPARTO(1, 4),NNU(1))   ,
     *            (NPARTO(1, 5),NMUP(1))  , (NPARTO(1, 6),NMUM(1))  ,
     *            (NPARTO(1, 7),NPI0(1))  , (NPARTO(1, 8),NPIP(1))  ,
     *            (NPARTO(1, 9),NPIM(1))  , (NPARTO(1,10),NK0L(1))  ,
     *            (NPARTO(1,11),NKPL(1))  , (NPARTO(1,12),NKMI(1))  ,
     *            (NPARTO(1,13),NNEUTR(1)), (NPARTO(1,14),NPROTO(1)),
     *            (NPARTO(1,15),NPROTB(1)), (NPARTO(1,16),NK0S(1))  ,
     *            (NPARTO(1,18),NHYP(1))  , (NPARTO(1,19),NDEUT(1)) ,
     *            (NPARTO(1,20),NTRIT(1)) , (NPARTO(1,21),NALPHA(1)),
     *            (NPARTO(1,22),NOTHER(1)), (NPARTO(1,25),NNEUTB(1))
*KEEP,OBSPAR.
      COMMON /OBSPAR/  OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
     *                 THETPR,PHIPR,NOBSLV
      DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
     *                 THETAP,THETPR(2),PHIP,PHIPR(2)
      INTEGER          NOBSLV
*KEEP,PAM.
      COMMON /PAM/     PAMA,SIGNUM
      DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
*KEEP,PARPAR.
      COMMON /PARPAR/  CURPAR,SECPAR,PRMPAR,OUTPAR,C,
     *                 E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
      DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
     *                 C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
      INTEGER          ITYPE,LEVL
*KEEP,PARPAE.
      DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
      EQUIVALENCE      (CURPAR(2),GAMMA),  (CURPAR(3),COSTHE),
     *                 (CURPAR(4), PHI ),  (CURPAR(5), H    ),
     *                 (CURPAR(6), T   ),  (CURPAR(7), X    ),
     *                 (CURPAR(8), Y   ),  (CURPAR(9), CHI  ),
     *                 (CURPAR(10),BETA),  (CURPAR(11),GCM  ),
     *                 (CURPAR(12),ECM )
*KEEP,PBALA.
      COMMON /PBALA/   PBAL
      DOUBLE PRECISION PBAL(10)
*KEEP,PRIMSP.
      COMMON /PRIMSP/  PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX,ISPEC
      DOUBLE PRECISION PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX
      INTEGER          ISPEC
*KEEP,RANDPA.
      COMMON /RANDPA/  FAC,U1,U2,RD,NSEQ,ISEED,KNOR
      DOUBLE PRECISION FAC,U1,U2
      REAL             RD(3000)
      INTEGER          ISEED(103,10),NSEQ
      LOGICAL          KNOR
*KEEP,RECORD.
      COMMON /RECORD/  IRECOR
      INTEGER          IRECOR
*KEEP,REJECT.
      COMMON /REJECT/  AVNREJ,
     *                 ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
     *                 FNPRIM
      DOUBLE PRECISION AVNREJ(10)
      REAL             ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
     *                 CUTLN,EONCUT
      LOGICAL          FNPRIM
*KEEP,RESON.
      COMMON /RESON/   RDRES,RESRAN,IRESPAR
      REAL             RDRES(2),RESRAN(1000)
      INTEGER          IRESPAR
 
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEEP,STACKF.
      COMMON /STACKF/  STACK,STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
      INTEGER          MAXSTK
      PARAMETER        (MAXSTK = 12*340*2)
      DOUBLE PRECISION STACK(MAXSTK)
      INTEGER          STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
*KEEP,STATI.
      COMMON /STATI/   SABIN,SBBIN,INBIN,IPBIN,IKBIN,IHBIN
      DOUBLE PRECISION SABIN(37),SBBIN(37)
      INTEGER          INBIN(37),IPBIN(37),IKBIN(37),IHBIN(37)
*KEEP,THNVAR.
      COMMON /THNVAR/  STACKINT,INT_ICOUNT,THINNING
      INTEGER          MAXICOUNT
      PARAMETER        (MAXICOUNT=20000)
      DOUBLE PRECISION STACKINT(MAXICOUNT,13)
      INTEGER          INT_ICOUNT
      LOGICAL          THINNING
*KEEP,VERS.
      COMMON /VERS/    VERNUM,MVDATE,VERDAT
      DOUBLE PRECISION VERNUM
      INTEGER          MVDATE
      CHARACTER*18     VERDAT
*KEEP,CEREN1.
      COMMON /CEREN1/  CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
     *                 CERSIZ,LCERFI
      DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
      REAL             CERSIZ
      LOGICAL          LCERFI
*KEEP,CEREN2.
      COMMON /CEREN2/  PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
     *                 DCERX,DCERY,ACERX,ACERY,
     *                 XCMAX,YCMAX,EPSX,EPSY,
     *                 DCERXI,DCERYI,FCERX,FCERY,
     *                 XSCATT,YSCATT,CERXOS,CERYOS,
     *                 NCERX,NCERY,ICERML
      REAL             PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
     *                 DCERX,DCERY,ACERX,ACERY,
     *                 XCMAX,YCMAX,EPSX,EPSY,
     *                 DCERXI,DCERYI,FCERX,FCERY,
     *                 XSCATT,YSCATT,CERXOS(20),CERYOS(20)
      INTEGER          NCERX,NCERY,ICERML
*KEEP,CEREN3.
      COMMON /CEREN3/  CERCNT,DATAB2,LHCER
      INTEGER          MAXBF2
      PARAMETER        (MAXBF2 = 39 * 7)
      DOUBLE PRECISION CERCNT
      REAL             DATAB2(MAXBF2)
      INTEGER          LHCER
*KEEP,CEREN4.
      COMMON /CEREN4/  NRECER
      INTEGER          NRECER
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*keep,certel.
      common /certel/  cormxd,cord,coralp,ctpars,omega,
     +                 photn,photnp,phpt,pht,vphot,
     +                 vchi,veta,vzeta,vchim,vetam,vzetam,
     +                 lambda,mu,nu,nctels,ncph
      double precision cormxd,cord,coralp,ctpars(20,7),omega(20,3,3),
     +                 photn(3),photnp(3),phpt(3),pht,vphot(3),
     +                 vchi(3),veta(3),vzeta(3),vchim,vetam,vzetam,
     +                 lambda,mu,nu
      integer          nctels,ncph(5)
      double precision xg,yg,zg,xgp,ygp,zgp,up,vp,wp,xpcut,ypcut,zpcut
      equivalence (photn(1) ,xg)   ,(photn(2) ,yg)   ,(photn(3) ,zg)  ,
     +            (photnp(1),xgp)  ,(photnp(2),ygp)  ,(photnp(3),zgp),
     +            (phpt(1)  ,xpcut),(phpt(2)  ,ypcut),(phpt(3)  ,zpcut),
     +            (vphot(1) ,up)   ,(vphot(2) ,vp)   ,(vphot(3) ,wp)    
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     Modificacion hecha por Aitor (5-feb-98)
      common /aitor/   aitoth
      double precision aitoth
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c     Angles for the "spinning" of a particle around the 
c     main axis of the CT
      common /spinang/ spinxi
      double precision spinxi
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


*KEND.
 
      INTEGER JNBIN(37),JPBIN(37),JKBIN(37),JHBIN(37)
      DOUBLE PRECISION CHI2,FPARAM(6)
      DOUBLE PRECISION MPART2(10,25),MPHOT2(10),MPOSI2(10),MELEC2(10),
     *                 MNU2(10),MMUP2(10),MMUM2(10),MPI02(10),MPIP2(10),
     *                 MPIM2(10),MK0L2(10),MKPL2(10),MKMI2(10),
     *                 MNETR2(10),MPROT2(10),MPRTB2(10),MK0S2(10),
     *                 MHYP2(10),MNETB2(10),MDEUT2(10),MTRIT2(10),
     *                 MALPH2(10),MOTH2(10)
      EQUIVALENCE (MPART2(1, 1),MPHOT2(1)), (MPART2(1, 2),MPOSI2(1)),
     *            (MPART2(1, 3),MELEC2(1)), (MPART2(1, 4),MNU2(1))  ,
     *            (MPART2(1, 5),MMUP2(1)) , (MPART2(1, 6),MMUM2(1)) ,
     *            (MPART2(1, 7),MPI02(1)) , (MPART2(1, 8),MPIP2(1)) ,
     *            (MPART2(1, 9),MPIM2(1)) , (MPART2(1,10),MK0L2(1)) ,
     *            (MPART2(1,11),MKPL2(1)) , (MPART2(1,12),MKMI2(1)) ,
     *            (MPART2(1,13),MNETR2(1)), (MPART2(1,14),MPROT2(1)),
     *            (MPART2(1,15),MPRTB2(1)), (MPART2(1,16),MK0S2(1)) ,
     *            (MPART2(1,18),MHYP2(1)) , (MPART2(1,19),MDEUT2(1)),
     *            (MPART2(1,20),MTRIT2(1)), (MPART2(1,21),MALPH2(1)),
     *            (MPART2(1,22),MOTH2 (1)), (MPART2(1,25),MNETB2(1))
C  VARIABLES BEING USED FOR RUNTIME
      REAL     TDIFF
      INTEGER  ILEFTA,ILEFTB,TIME
      EXTERNAL TIME

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      double precision ctdiams(20)
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      double precision theprim, phiprim
      double precision spinthe, spinphi
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


C-----------------------------------------------------------------------
 
 
      CERELE = 0.D0
      CERHAD = 0.D0
      NRECER = 0
C  INITIALIZE AND READ RUN STEERING CARDS
      CALL START
 
      IF ( CERSIZ .LE. 0. ) THEN
        ICRSIZ = 0
      ELSE
        ICRSIZ = 1
      ENDIF
 
C  RESET COUNTER FOR WORDS WRITTEN TO TAPE
      IRECOR = 0
 
C  RESET COUNTER FOR AVERAGE HIGHT OF 1ST INTERACTION
      CHISUM = 0.D0
      CHISM2 = 0.D0
 
C  SET ARRAYS FOR SCALES OF KINETIC ENERGY-INTERACTION TABLE
      SABIN(1) = 0.D0
      SBBIN(1) = 0.1D0
      DO 13  J = 2,37
        SABIN(J) = 10.D0**((J-4.D0)/3.D0)
        SBBIN(J) = 10.D0**((J-3.D0)/3.D0)
  13  CONTINUE
 
C  CHECK AND SET PRIMARY PARAMETERS
      CALL INPRM

      do 161 i=1,nctels
        ctdiams(i) = ctpars(i,ctdiam)
 161  continue

C  INITIALIZE NKG ROUTINES
      CALL ININKG
 
C  RESET COUNTERS FOR NUCLEON, PION AND KAON TABLE FOR ALL SHOWERS
C  RESET ENERGY-MULTIPLICITY & ENERGY-ELASTICITY MATRIX FOR ALL SHOWERS
      DO 17  J = 1,37
        JNBIN(J)  = 0
        JPBIN(J)  = 0
        JKBIN(J)  = 0
        JHBIN(J)  = 0
        ELMEAA(J) = 0.D0
        DO 17  L = 1,13
          MULTOT(J,L) = 0
          IELDPA(J,L) = 0
  17  CONTINUE
 
C  RESET OTHER ARRAYS FOR ALL SHOWERS
      DO 99  J = 1,20
        IELNU(J) = 0
        IELPI(J) = 0
        IELIS(J) = 0
        IELHM(J) = 0
  99  CONTINUE
 
C  RESET ARRAYS FOR INTERACTION LENGTH STATISTICS
      DO 90  J = 1,124
        IHYCHI(J) = 0
        IKACHI(J) = 0
        IMUCHI(J) = 0
        INUCHI(J) = 0
        IPICHI(J) = 0
        INNCHI(J) = 0
  90  CONTINUE
 
C  RESET ARRAY FOR MEAN VALUES AND STANDARD DEVIATION
      DO 477  K = 1,25
      DO 477  J = 1,10
        MPARTO(J,K) = 0.D0
        MPART2(J,K) = 0.D0
 477  CONTINUE
 
C  RESET ARRAYS FOR LONGITUDINAL DISTRIBUTION
      IF ( LLONGI ) THEN
        DO 478  K = 1,9
          DO 4781  J = 0,NSTEP
            APLONG(J,K) = 0.D0
            SPLONG(J,K) = 0.D0
 4781     CONTINUE
 478    CONTINUE
      ENDIF
 
C  STEERING OF PRINTOUT OF RANDOM GENERATOR SEEDS
      IPROUT = MIN(100,NSHOW/20)
      IPROUT = MAX(1,IPROUT)
 
C  TIME AT BEGINNING
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c      ILEFTA = TIME()
      ILEFTA = 0
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     Modified by C. Bigongiari 2001 Jan 16 
C 
C
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      print *,'JCIO::========================================'
      print *,'JCIO:: Initializing JCIO system for advanced'
      print *,'JCIO:: saving of data.'
      print *,'JCIO::========================================'
C
Cc- initialize jcio system
C      
      call jcinitio(dsn,nrrun)
Cc- create file run######
C      call jcstartrun(runh)
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C- Modified JCSTARTRUN creates cer###### and dat###### files !
C
C     ###### is the RUN number !
C

      call jcstartrun(RUNH)

C
C- write Run Header on cer and dat files  
C 
      CALL TOBUF(RUNH,0)
      IF ( LCERFI ) CALL TOBUFC(RUNH,0)

C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C-----------------------------------------------------------------------
C  LOOP OVER SHOWERS

      DO 2  ISHW = 1,NSHOW

 
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c   Next block of code has been modified, and comes from INPRM 
c----------------------------------------------------------------------
C  SCATTERING OF CENTER OF CHERENKOV ARRAY RELATIVE TO SHOWER AXIS
c>> Actually, XSCATT and YSCATT should be RminSCATT and RmaxSCATT
        ICERML = MIN(MAX(ICERML,1),20)
        XSCATT = ABS(XSCATT)
        YSCATT = ABS(YSCATT)
        WRITE(MONIOU,5225)ICERML,XSCATT,YSCATT
 5225   FORMAT(' ** USING EACH SHOWER SEVERAL TIMES:'/ 
     +     ' USE EACH EVENT ',I2,' TIMES'/
     +     ' THE EVENTS ARE SCATTERED RANDOMLY IN A SECTOR OF RADII:'/
     +     '   Rmin = ',F10.0,'   Rmax = ',F10.0)
        DO 4438 I=1,ICERML
 5226     CALL RMMAR( RD,2,3 )
          CERXOS(I) = 2.0*YSCATT*(RD(1)-0.5)
          CERYOS(I) = 2.0*YSCATT*(RD(2)-0.5)
          R=SQRT(CERXOS(I)**2+CERYOS(I)**2)
          IF ((R.LT.XSCATT).OR.(R.GT.YSCATT)) GOTO 5226
          WRITE(MONIOU,4437) I,CERXOS(I),CERYOS(I)
 4437     FORMAT('    CORE OF EVENT ',I2,'  AT  ',2F12.2)
 4438   CONTINUE
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c   Next block of code comes from INPRM 
c----------------------------------------------------------------------
        EVTH(98) = FLOAT(ICERML)
        DO  480 I=1,20
          EVTH( 98+I) = CERXOS(I)
          EVTH(118+I) = CERYOS(I)
 480    CONTINUE
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

        SHOWNO = SHOWNO + 1
        I      = ISHW
        IF ( ISHW .LE. MAXPRT ) THEN
          FPRINT = .TRUE.
        ELSE
          FPRINT = .FALSE.
        ENDIF
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Cc   Create cer######,dat######,sta###### files
Cc------------------------------------------------------------
C        call jcnewshower
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
C  RESET COUNTERS
        DO 447  K = 1,25
        DO 447  J = 1,10
          NPARTO(J,K) = 0.D0
 447    CONTINUE
        MUOND = 0.D0
 
C  RESET ARRAY FOR LONGITUDINAL DISTRIBUTION PER SHOWER
        IF ( LLONGI ) THEN
          DO 479  K = 1,9
            DO 4791  J = 0,NSTEP
              PLONG(J,K) = 0.D0
 4791       CONTINUE
 479      CONTINUE
        ENDIF
 
        NRECS = 0
        NBLKS = 0
        DO 922  KKK = 1,10
          AVNREJ(KKK) = 0.D0
 922    CONTINUE
        IRESPAR = 0
 
 
C  FIRST INTERACTION DATA
        FIRSTI = .TRUE.
        IFINET = 0
        IFINNU = 0
        IFINKA = 0
        IFINPI = 0
        IFINHY = 0
        ELAST  = 0.D0
 
C  RESET COUNTERS FOR NUCLEON, PION AND KAON TABLE FOR SHOWER
C  RESET ENERGY-MULTIPLICITY & ENERGY-ELASTICITY MATRIX FOR SHOWER
        DO 11  J = 1,37
          INBIN(J) = 0
          IPBIN(J) = 0
          IKBIN(J) = 0
          IHBIN(J) = 0
          ELMEAN(J) = 0.D0
          DO 11  L = 1,13
            MULTMA(J,L) = 0
            IELDPM(J,L) = 0
  11    CONTINUE
 
        DO 12  J = 1,10
          PBAL(J) = 0.D0
          EBAL(J) = 0.D0
  12    CONTINUE
 
C  INITIALIZE PARTICLE STACK
        CALL ISTACK
C  RESET STACKINT
        DO J=1,MAXICOUNT
          DO K=1,MAXLEN
            STACKINT(J,K) = 0.D0
          ENDDO
        ENDDO
 
C  INITIALIZE EVENT HEADER AND END FOR EACH EVENT
        DO 2123  L = 2,43
          EVTH(L) = 0.
 2123   CONTINUE
        DO 123  L = 2,MAXBUF
          EVTE(L) = 0.
 123    CONTINUE
 
C  SHOWER BEGIN PRINTOUT
        IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,105) SHOWNO
 105    FORMAT ('1',10('='),' SHOWER NO ',I10,' ',47('=')/)
 
C  RANDOM GENERATOR STATUS AT BEGINNING OF SHOWER CALCULATION
        EVTH(13) = NSEQ
        DO 45  L = 1,NSEQ
          CALL RMMAQ( ISEED(1,L), L, 'R' )
C  SEED
          EVTH(11+L*3) = ISEED(1,L)
C  NUMBER OF CALLS
          EVTH(12+L*3) = MOD ( ISEED(2,L), 1000000 )
C  NUMBER OF MILLIONS
          EVTH(13+L*3) = ISEED(3,L)*1000 + INT( ISEED(2,L)/1000000 )
  45    CONTINUE
        IF ( FPRINT  .OR.  DEBUG  .OR.  MOD(ISHW-1,IPROUT).EQ.0 ) THEN
          CALL PRTIME(TTIME)
          WRITE(MONIOU,158) SHOWNO,(L,(ISEED(J,L),J=1,3),L=1,NSEQ)
 158      FORMAT(' AND RANDOM NUMBER GENERATOR AT BEGIN OF EVENT :',I8,
     *            /,(' SEQUENCE = ',I2,'  SEED = ',I9 ,'  CALLS = ',I9,
     *               '  BILLIONS = ',I9))
        ENDIF
C  RESET KNOR
        KNOR = .TRUE.
 
C  GET FULL RANDOM GENERATOR STATUS (103 WORDS PER SEQUENCE)
CC      DO 495  L = 1,NSEQ
CC        CALL RMMAQ( ISEED(1,L), L, 'RV' )
CC        WRITE(MONIOU,658) L,(ISEED(J,L),J=1,103)
CC658     FORMAT ( ' FULL RANDOM NUMBER GENERATOR STATUS ',
CC   *             'FOR SEQUENCE ',I2,/(' ',10I11))
CC495   CONTINUE
 
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> *** ATENTION *** ATENTION *** ATENTION *** ATENTION *** ATENTION >>
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>>                                                                  >>
c>> In the next block (between this ATENTION comments) CORSIKA makes >>
c>> three things, in this order:                                     >>
c>>                                                                  >>
c>>   i. Set ANGLES OF INCIDENCE (different distributions of theta   >>
c>>      for gammas -flat- and hadrons -standard.                    >>
c>>  ii. Set HEIGHT for start at THICK0 (normally = 0 => 112.8 Km)   >>
c>> iii. Set ENERGY of the primary.                                  >>
c>>                                                                  >>
c>> (The original order was ii., iii. and i.)                        >>
c>>                                                                  >>
c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  GET PRIMARY ANGLES OF INCIDENCE
        IF ( FIXINC ) THEN

          THETAP = THETPR(1)
          PHIP   = PHIPR(1)
          PRMPAR(3) = COS(THETAP)

        ELSE

          if ( prmpar(1).eq.1 ) then

C>> GAMMAS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C NOTE!! We will use a FLAT distribution for THETA: 
C Then, next block (original block) must be commented.
c The modificated code follows this block
c
            CALL RMMAR( RD,3,1 )
            CT1 = THETPR(1)
            CT2 = THETPR(2)
            THETAP = RD(2)*(CT2 - CT1) + CT1
            CTT  = COS(THETAP)
            PRMPAR(3) = CTT

          else

C>> HADRONS AND ELECTRONS (AND ANY OTHER BUT GAMMAS) >>>>>>>>>>>>>>>>
c  Choose angles at random with equal flux for all directions
c  with horizontal detector array (see: O.C. Allkofer & P.K.F. Grieder,
c  Cosmic Rays on Earth, in: Physics Data 25/1, H.Behrens & G.Ebel Ed.,
c  (Fachinformationszentrum Karlsruhe, Germany, 1983) chpt. 1.1.2)
c
            CALL RMMAR( RD,3,1 )
            CT1 = SIN(THETPR(1))**2
            CT2 = SIN(THETPR(2))**2
            CTT  = SQRT( 1.D0 - RD(2)*(CT2 - CT1) - CT1 )
            PRMPAR(3) = CTT
            THETAP = ACOS(CTT)

          endif

          PHIP = RD(1) * ( PHIPR(2) - PHIPR(1) ) + PHIPR(1)

        ENDIF
        PRMPAR(4) = PHIP


c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  DEFINE HEIGHT FOR START AT THICK0 (IN G/CM**2) (112.8 KM FOR THICK0=0)
        PRMPAR(5) = HEIGH(THICK0)
 
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  GET PRIMARY ENERGY INTO PRMPAR(2)
        IF ( ISPEC .EQ. 0 ) THEN
          PRMPAR(2) = LLIMIT
        ELSE
          CALL RMMAR( RD,1,1 )
          IF ( PSLOPE .NE. -1.D0 ) THEN
            PRMPAR(2) = ( RD(1)*UL + ( 1.D0-RD(1) )*LL )**SLEX
          ELSE
            PRMPAR(2) = LLIMIT * LL**RD(1)
          ENDIF
        ENDIF

c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> *** ATENTION *** ATENTION *** ATENTION *** ATENTION *** ATENTION >>
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> Modification: this is no longer needed >>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> (Superseeded by Sphere algorithm, see cerenkov.f) >>>>>>>>>>>>>>>>>
c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
cc Btw, we now modify the "shadow area" of the telescopes,
cc to cover the angle theta. 
c        do 160 i=1,nctels
c          ctpars(i,ctdiam) = ctdiams(i)/cos(thetap)
c          write (MONIOU,*) 
c     *        'New region for CT',i,' = ',ctpars(i,ctdiam)
c 160    continue
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c Here we calculate the angles spinphi and spinthe, which are the 
c phi and theta angles of the particle, with respect to the direction
c where the CT is pointing to. spinthe is the angular displacement
c of the new direction respect to the original (CT); spinphi=0 means
c that the new direction is towards the zenith, spinphi=+-180 means 
c towards the horizont.
c See the document "simulation.tex" 

c First, save the "CT" orientation 
c (moved from a couple of lines below, marked with [*])
        EVTH(11) = THETAP
        EVTH(12) = PHIP
        
        CALL RMMAR( RD,3,1 )

c Then, calculate the new direction relative to the CT direction
        spinphi = RD(1)*PI
        spinthe = RD(2)*spinxi*pi/180

c And then, RE-calculate the GLOBAL direction in CORSIKA
c We use formulae for spherical triangles

        if ( THETAP .EQ. 0.0 ) then                   
           theprim = spinthe 
           phiprim = spinphi
        elseif ( spinthe .eq. 0.0 ) then 
           theprim = THETAP 
           phiprim = 0.0
        else 
           theprim = acos( cos(THETAP)*cos(spinthe)+
     $                  sin(THETAP)*sin(spinthe)*cos(spinphi) )
           phiprim = asin( sin(spinthe)*sin(spinphi)/sin(theprim) )
        endif 
        THETAP = theprim
        EVTH(140) = spinthe
        if (RD(3).gt.0.5) then
          PHIP = PHIP - phiprim 
          EVTH(141) = -spinphi
        else
          PHIP = PHIP + phiprim 
          EVTH(141) = spinphi
        endif
        PRMPAR(3) = COS(THETAP)
        PRMPAR(4) = PHIP

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> Modification (HZA trick) cancelled >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C     Modificacion hecha por Aitor
c         aitoth = THETAP
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

C  CALCULATE COORDINATE CORRECTION FOR TOP OF ATMOSPHERE
        CALL COORIN( PRMPAR(5) )
 
C  COUNTER FOR PARTICLE OUTPUT
        LH = 0
C  COUNTER FOR CERENKOV OUTPUT
        IF ( LCERFI ) LHCER = 0
C  CALCULATE BUNCH SIZE FOR CERENKOV PHOTONS IF NOT SET IN DATAC
        IF ( ICRSIZ .EQ. 0 ) THEN
          CALL GETBUS( NINT(PRMPAR(1)),SNGL(PRMPAR(2)),SNGL(PRMPAR(3)),
     *                  CERSIZ )
          WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS CALCULATED TO=',CERSIZ
        ENDIF
C  GET GAMMA FACTOR FROM ENERGY
C  FOR GAMMAS PRMPAR(2) STAYS = ENERGY
        IF ( PRMPAR(1) .NE. 1.D0 )
     *              PRMPAR(2) = PRMPAR(2) / PAMA(NINT(PRMPAR(1)))
 
C  SET PRIMARY TO CURRENT PARTICLE
        DO 3  J = 1,8
          CURPAR(J) = PRMPAR(J)
          NCOUN(J)  = 0
  3     CONTINUE
C  SET WEIGHT
 
C  CALCULATE FIRST INTERACTION POINT IF HADRONIC
        GEN = 0.D0
 
        H = HEIGH(THICK0)
        CALL BOX2
        IF ( FIX1I ) THEN
          CHI = THICK(FIXHEI) / PRMPAR(3)
          H = FIXHEI
          FDECAY = .FALSE.
        ELSE
          H = HEIGH ( CHI*PRMPAR(3) + THICK0 )
        ENDIF
        CHISUM = CHISUM + CHI
        CHISM2 = CHISM2 + CHI**2
        ALEVEL = H
C  INITIALIZE COORDINATE CORRECTIONS FOR HADRONIC PRIMARIES
C  FOR EM PRIMARIES IT IS DONE IN EGS
        HH = MAX( H, 0.D0 )
        IF ( CURPAR(1) .GT. 3.D0 ) CALL COORIN( HH )
 
        IF ( FMUADD ) THEN
          IF ( CURPAR(1) .EQ. 5  .OR.  CURPAR(1) .EQ. 6) THEN
            DO J = 1,MAXLEN
              AMUPAR(J) = CURPAR(J)
            ENDDO
            AMUPAR(5) = PRMPAR(5)
            IF(DEBUG)WRITE(MDEBUG,*)'MAIN  : MUON STORED IN AMUPAR'
            FMUORG = .TRUE.
          ENDIF
        ENDIF
 
C  SET TARGET FLAG IF SELECTED FOR FIRST INTERACTION
        IF ( N1STTR .GT. 0 ) THEN
          FIXTAR  = .TRUE.
          FDECAY  = .FALSE.
          EVTH(6) = REAL(N1STTR)
        ELSE
          FIXTAR  = .FALSE.
          EVTH(6) = 0.0
        ENDIF
 
C  INITIALIZE ARRAYS FOR NKG FOR EACH SHOWER
        IF ( FNKG ) CALL STANKG
 
C  STORE FIRST PARTICLE IN HEADER AND PRINT IT OUT
        EVTH( 2) = REAL(SHOWNO)
        EVTH( 3) = CURPAR(1)
        IF ( CURPAR(1) .EQ. 1.D0 ) THEN
C  PRIMARY ENERGY FOR PHOTONS
          E00    = GAMMA
          E00PN  = GAMMA
          INUCL  = 1
        ELSE
          E00    = GAMMA * PAMA(NINT(CURPAR(1)))
          INUCL  = INT(MAX(1.D0,CURPAR(1)/100.D0))
          E00PN  = E00 / INUCL
        ENDIF
        EVTH(147) = 0.

        IF ( FEGS ) THEN
C  PARAMETER FOR ELECTRON AND PHOTON REJECT (CONVERT ENERGY TO MEV)
          EONCUT = .5E-9*SQRT(E00*1000.D0)
          CUTLN  = LOG(EONCUT)
        ENDIF
        EVTH( 4) = E00
        EVTH( 5) = THICK0
        EVTH( 7) = H
        PTOT0    = SQRT( E00**2 - PAMA(NINT(CURPAR(1)))**2 )
        PTOT0N   = PTOT0 / INUCL
        ST       = SQRT(1.D0-COSTHE**2)
        EVTH( 8) = PTOT0 * ST * COS(PHI)
        EVTH( 9) = PTOT0 * ST * SIN(PHI)
        EVTH(10) = PTOT0 * COSTHE
     

c
c [*] one block from here sent above
c

        EVTH(85) = CERSIZ
 
        IF ( CURPAR(1) .GT. 3.D0 ) THEN
          IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,102) (CURPAR(J),J = 1,8)
 102      FORMAT (/' PRIMARY PARAMETERS AT FIRST INTERACTION POINT'/
     *               16X,1P,8E10.3)
        ELSE
          IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,132)
 132      FORMAT (/' PRIMARY PARTICLE IS ELECTROMAGNETIC')
        ENDIF
 
C  WRITE EVENT HEADER INTO BUFFER
C  FOR EM PARTICLES EVTH IS WRITTEN TO BUFFER IN EGS (IF ACTIVE)
        IF ( EVTH(3) .GT. 3.0  .OR.  .NOT. FEGS ) THEN
          CALL TOBUF ( EVTH,0 )
          IF ( LCERFI ) CALL TOBUFC( EVTH,0 ) 
        ENDIF
 
C  PRINT HEADER FOR HIGH ENERGY PARTICLES
        IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,103)
 103    FORMAT(/'                   TYPE      GAMMA   COSTHETA ',
     *          '    PHI     HEIGHT     TIME      X-CM      Y-CM   ',
     *          '    GEN      LEVEL  E ON STACK'/)
        NOPART = 0

        IF ( CURPAR(1) .LE. 3.D0  .OR.
     *      (CURPAR(1) .GE. 5.D0  .AND.  CURPAR(1) .LE. 7.D0) ) THEN
C  GIVE PARTICLE TO EGS OR NKG IF ELECTROMAGNETIC
C  AND TAKE THEN NEXT PARTICLE FROM STACK
C  FLAG FOR NO PRIMARY INTERACTION IS SET FOR ALL BUT ELM. PRIMARIES
          IF ( CURPAR(1) .LE. 3.D0 ) THEN
            FNPRIM = .FALSE.
          ELSE
            FNPRIM = .TRUE.
            H = PRMPAR(5)
          ENDIF
          CALL BOX3
          IF ( FEGS ) THEN
            CHISUM = CHISUM + THICK(DBLE(EVTH(7)))
            CHISM2 = CHISM2 + THICK(DBLE(EVTH(7)))**2
          ENDIF
          FIRSTI = .FALSE.
          GOTO 4
 
        ELSE
C  HADRONIC PARTICLES
          FNPRIM = .TRUE.

C  FILL LONGITUDINAL DISTRIBUTION FOR THE PRIMARY PARTICLE
C  THE PARTICLE IS TRACKED FROM THICK0 DOWN TO THICK0+CHI*PRMPAR(3)
C  COUNT THE PARTICLES FOR THE LONGITUDINAL DEVELOPMENT
          IF ( LLONGI ) THEN
            LPCT1 = INT( THICK0 * THSTPI )
            LPCT2 = INT( (THICK0 + PRMPAR(3)*CHI) * THSTPI )
            LPCT2 = MIN(NSTEP,LPCT2)
C  GAMMAS, ELECTRONS AND POSITRONS ARE NOT TRANSPORTED HERE, SEE EGS
C  MUONS ARE TRANSPORTED IN MUTRAC
C  HADRONS
            IF     ( ITYPE .GE. 7 .AND. ITYPE .LE. 41 ) THEN
              DO 5004 L = LPCT1,LPCT2
                PLONG(L,6) = PLONG(L,6) + 1.
 5004         CONTINUE
C  CHARGED HADRONS
              IF ( SIGNUM(ITYPE) .NE. 0.D0 ) THEN
                DO 5005 L = LPCT1,LPCT2
                  PLONG(L,7) = PLONG(L,7) + 1.
 5005           CONTINUE
              ENDIF
C  NUCLEI
            ELSEIF ( ITYPE .GT. 100 ) THEN
              DO 5006 L = LPCT1,LPCT2
                PLONG(L,8) = PLONG(L,8) + 1.
 5006         CONTINUE
            ENDIF
          ENDIF

C  CHECK OBSERVATION LEVEL PASSAGE AND UPDATE PARTICLE COORDINATES
          HNEW = H
C  FOR UPDATE WE NEED THE START ALTITUDE H
          H = HEIGH(THICK0)
          DO  251  J = 1,NOBSLV
C  JUMP INTO NORMAL PARTICLE TREATMENT FOR HADRONS
            IF ( HNEW .GT. OBSLEV(J) ) THEN
              H = HNEW
              GOTO 6
            ENDIF
            IF ( H .LT. OBSLEV(J) ) GOTO 251
C  REMEMBER NUMBER OF LEVEL FOR OUTPUT
            LEVL   = J
            CALL UPDATE( OBSLEV(J),THCKOB(J),J )
            IF (DEBUG) WRITE(MDEBUG,256) J,IRET1,IRET2
 256        FORMAT(' MAIN  : LEVEL ',I5,' IRET1,2=',2I5)
C  IF PARTICLE IS NOT CUTTED, BRING IT TO OUTPUT
            IF ( IRET2 .EQ. 0 ) THEN
              CALL OUTPUT
            ENDIF
 251      CONTINUE
          IF (DEBUG) WRITE(MDEBUG,*)
     *       'MAIN  : PRIMARY REACHED LOWEST OBSERVATION LEVEL'
          GOTO 4
        ENDIF
 

C-----------------------------------------------------------------------
C  NORMAL CYCLE
  7     CONTINUE
 
C  IF ENERGY TOO SMALL TAKE NEXT PARTICLE
        IF ( GAMMA .LE. 1.D0 ) THEN
          IF ( CURPAR(1) .NE. 1.D0 ) THEN
            IF ( CURPAR(1).EQ.5.D0 .OR. CURPAR(1).EQ.6.D0 )
     *                                        FMUORG = .FALSE.
            GOTO 4
          ENDIF
C  SPECIAL TREATMENT FOR PHOTONS
          ITYPE = 1
          CHI   = 0.D0
          GOTO 5
        ENDIF
 
C  DETERMINE PLACE OF NEXT INTERACTION
        CALL BOX2
 
C  CHECK PASSAGE THROUGH OBSERVATION LEVELS AND TRACK PARTICLES TO THE
C  PLACE OF INTERACTION
  5     CONTINUE
        IRET1 = 0
        CALL BOX3
        IF ( IRET1 .NE. 0 ) GOTO 4
 
  6     CONTINUE
        IRET1 = 0
        MSMM  = 0
 
C  INCREMENT PARTICLE GENERATION AND PROCESS NUCLEAR INTERACTION
        GEN   = GEN + 1.D0
C  INITIALIZE INTERMEDIATE STACK FOR ONE REACTION
        CALL TSTINI
        CALL NUCINT
C  TRANSFER INTERMEDIATE STACK FOR ONE REACTION
        CALL TSTEND
 
C  ENERGY - MULTIPLICITY STATISTICS
        IF ( EKINL .LE. 0.1D0 ) THEN
          MEN = 1
        ELSE
          MEN = 4.D0 + 3.D0 * LOG10(EKINL)
          MEN = MIN( MEN, 37 )
        ENDIF
        IF ( MSMM .LE. 1 ) THEN
          MMU = 1
        ELSE
          MMU = 1.D0 + 3.D0 * LOG10(DBLE(MSMM))
          MMU = MIN( MMU, 13 )
        ENDIF
        MULTMA(MEN,MMU) = MULTMA(MEN,MMU) + 1
        MULTOT(MEN,MMU) = MULTOT(MEN,MMU) + 1
        IF ( DEBUG ) WRITE(MDEBUG,*) 'MAIN  : EKINL,MSMM=',
     *                                  SNGL(EKINL),MSMM
 
        IF ( IRET1 .EQ. 0 ) THEN
          IF ( DEBUG ) WRITE(MDEBUG,666) (CURPAR(II),II=1,11)
 666      FORMAT(' MAIN  : CURPAR=',1P,11E10.3)
          GOTO 7
        ENDIF
 
C  GET NEXT PARTICLE FROM STACK, IF IRET=1 ALL PARTICLES ARE DONE
  4     CONTINUE
        IRET1 = 0
        CALL FSTACK
        IF ( FMUADD ) THEN
          IF ( (CURPAR(1) .EQ. 5  .OR.  CURPAR(1) .EQ. 6)
     *         .AND.  IRET1 .EQ. 0  .AND.  .NOT. FMUORG ) THEN
            DO J = 1,MAXLEN
              AMUPAR(J) = CURPAR(J)
            ENDDO
            IF(DEBUG)WRITE(MDEBUG,*)'MAIN  : MUON STORED IN AMUPAR'
            FMUORG = .TRUE.
          ENDIF
        ENDIF
 
        IF ( IRET1 .EQ. 0 ) GOTO 7
 
C-----------------------------------------------------------------------
C  FINISH SHOWER AND PRINT INFORMATION
        CALL OUTEND
 
*       IF ( DEBUG ) WRITE(MDEBUG,442) NPARTO
*442    FORMAT(' MAIN  : NPARTO='/(' ',10F10.0))
 
        IF ( FPRINT .OR. DEBUG ) THEN
          IOBSLV = MIN( 5, NOBSLV )
          WRITE(MONIOU,54) (K,K=1,IOBSLV)
  54      FORMAT (/' PARTICLES AT DETECTOR LEVEL :'/
     *             ' FOR LEVEL         ', 5I13)
          WRITE(MONIOU,55) (OBSLEV(K),K=1,IOBSLV)
  55      FORMAT ( ' HEIGHT IN CM        ',1P, 5E13.3/)
          WRITE(MONIOU,776) 'PROTONS      ',(NPROTO(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'ANTIPROTONS  ',(NPROTB(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'NEUTRONS     ',(NNEUTR(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'ANTINEUTRONS ',(NNEUTB(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'PHOTONS      ',(NPHOTO(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'ELECTRONS    ',(NELECT(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'POSITRONS    ',(NPOSIT(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'NEUTRINOS    ',(NNU   (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'MU -         ',(NMUM  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'MU +         ',(NMUP  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'PI 0         ',(NPI0  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'PI -         ',(NPIM  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'PI +         ',(NPIP  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'K0L          ',(NK0L  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'K0S          ',(NK0S  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'K -          ',(NKMI  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'K +          ',(NKPL  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'STR. BARYONS ',(NHYP  (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'DEUTERONS    ',(NDEUT (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'TRITONS      ',(NTRIT (K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'ALPHAS       ',(NALPHA(K),K=1,IOBSLV)
          WRITE(MONIOU,776) 'OTHER PARTIC.',(NOTHER(K),K=1,IOBSLV)
          WRITE(MONIOU,*)
          WRITE(MONIOU,776) 'DECAYED MUONS',MUOND
 776      FORMAT(' NO OF ',A13, '= ',5F13.0)
 
          IF ( NOBSLV .GT. 5 ) THEN
            IOBSLV =  NOBSLV
            WRITE(MONIOU,54) (K,K=6,IOBSLV)
            WRITE(MONIOU,55) (OBSLEV(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'PROTONS      ',(NPROTO(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'ANTIPROTONS  ',(NPROTB(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'NEUTRONS     ',(NNEUTR(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'ANTINEUTRONS ',(NNEUTB(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'PHOTONS      ',(NPHOTO(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'ELECTRONS    ',(NELECT(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'POSITRONS    ',(NPOSIT(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'NEUTRINOS    ',(NNU   (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'MU -         ',(NMUM  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'MU +         ',(NMUP  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'PI 0         ',(NPI0  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'PI -         ',(NPIM  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'PI +         ',(NPIP  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'K0L          ',(NK0L  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'K0S          ',(NK0S  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'K -          ',(NKMI  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'K +          ',(NKPL  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'STR. BARYONS ',(NHYP  (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'DEUTERONS    ',(NDEUT (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'TRITONS      ',(NTRIT (K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'ALPHAS       ',(NALPHA(K),K=6,IOBSLV)
            WRITE(MONIOU,776) 'OTHER PARTIC.',(NOTHER(K),K=6,IOBSLV)
            WRITE(MONIOU,*)
          ENDIF
        ENDIF
 
C  ADD UP FOR MEAN VALUES
        DO 779  K = 1,25
        DO 779  J = 1,10
          MPARTO(J,K) = MPARTO(J,K) + NPARTO(J,K)
          MPART2(J,K) = MPART2(J,K) + NPARTO(J,K)**2
 779    CONTINUE
        EVTE(2) = SHOWNO
        DO 335  K = 1,NOBSLV
          EVTE(3) = EVTE(3) + NPHOTO(K)
          EVTE(4) = EVTE(4) + NELECT(K) + NPOSIT(K)
          EVTE(5) = EVTE(5) + NPROTO(K) + NPROTB(K) + NNEUTR(K) +
     *              NNEUTB(K) + NPI0(K) + NPIM(K) + NPIP(K) + NK0L(K) +
     *              NK0S(K) + NKMI(K) + NKPL(K) + NHYP(K) +
     *              NDEUT(K) + NTRIT(K) + NALPHA(K) + NOTHER(K)
          EVTE(6) = EVTE(6) + NMUP(K) + NMUM(K)
 335    CONTINUE
        EVTE(7)   = NOPART
 
        IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,110)
     *                  IFINNU,IFINPI,IFINET,IFINKA,IFINHY,
     *                  IFINNU+IFINPI+IFINET+IFINKA+IFINHY,ELAST
 110    FORMAT(/' NO OF NUCLEONS  PRODUCED IN FIRST INTERACTION =',I10/
     *          ' NO OF PIONS     PRODUCED IN FIRST INTERACTION =',I10/
     *          ' NO OF ETAS      PRODUCED IN FIRST INTERACTION =',I10/
     *          ' NO OF KAONS     PRODUCED IN FIRST INTERACTION =',I10/
     *          ' NO OF S.BARYONS PRODUCED IN FIRST INTERACTION =',I10/
     *          ' TOTAL MULTIPLICITY       OF FIRST INTERACTION =',I10/
     *        ' ELASTICITY               OF FIRST INTERACTION =',F10.4)
 
C  PRINT OUT NKG RESULT FOR ONE SHOWER IF SELECTED
        IF ( FNKG ) CALL AVAGE
 
        IF ( LLONGI ) THEN
C  TREAT LONGITUDINAL DISTRIBUTIONS
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c   calculated here again, 'cos it's rewrite I dont know where
          LPCT1 = INT( THICK0 * THSTPI )
          LPCT2 = INT( (THICK0 + PRMPAR(3)*CHI) * THSTPI )
          LPCT2 = MIN(NSTEP,LPCT2)
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          DO 980  J = LPCT1,NSTEP
C  ADD ELECTRONS, POSITRONS, MUONS AND NUCLEI TO THE CHARGED PARTICLES
            PLONG(J,7) = PLONG(J,7) + PLONG(J,2) + PLONG(J,3)
     *                 + PLONG(J,4) + PLONG(J,5) + PLONG(J,8)
C  ADD UP FOR MEAN VALUES OF LONGITUDINAL DISTRIBUTION
            DO 979  K = 1,9
              APLONG(J,K) = APLONG(J,K) + PLONG(J,K)
              SPLONG(J,K) = SPLONG(J,K) + PLONG(J,K)**2
 979        CONTINUE
 980      CONTINUE
 
C  PRINT LONGITUDINAL DISTRIBUTIONS PER SHOWER
          IF ( FPRINT .OR. DEBUG )  WRITE(MONIOU,910) THSTEP,
     *      'GAMMAS','POSITRONS','ELECTRONS','MU-','MU+','HADRONS',
     *      'CHARGED','NUCLEI','CERENKOV',
     *      (J*THSTEP,(PLONG(J,K),K=1,9),J=LPCT1,NSTEP)
 910      FORMAT(/' ---------- LONGITUDINAL DISTRIBUTION IN STEPS OF ',
     *        F5.0,' G/CM**2 ----------------'/
     *        '  DEPTH ',3A12,3A11,A12,A11,A12/
     *        (' ',F6.0,F13.0,2F12.0,3F11.0,F12.0,F11.0,1P,E12.5,0P) )
CJOK  ADAPTED FOR HEAT CALCULATION
C910      FORMAT(/
C    *    ' LONGITUDINAL DISTRIBUTION IN STEPS OF ',F5.0,' G/CM**2'
C    *      /' ',92('=')/'  DEPTH',8A10,A12/1P
C    *      (' ',0P,F6.0,1P,9E11.4))
CJOK
 
          IF ( FLGFIT ) THEN
C  PERFORM FIT TO THE LONGITUDINAL DISTRIBUTION OF ALL CHARGED PARTICLES
C  IF EGS IS SELECTED THIS IS THE DISTRIBUTION WHICH IS TO BE TAKEN
            IF ( FEGS ) THEN
              DO 930 J=0,NSTEP-LPCT1
                DEP(J+1)    = (J+LPCT1)*THSTEP
                CHAPAR(J+1) = PLONG(J+LPCT1,7)
 930          CONTINUE
              NSTP = NSTEP + 1 - LPCT1
              WRITE(MONIOU,8229) 'ALL CHARGED PARTICLES'
 8229         FORMAT(/' FIT OF THE CURVE   ',
     *      ' N(T) = P1*((T-P2)/(P3-P2))**((P3-T)/(P4+P5*T+P6*T**2))'/
     *      ' TO LONGITUDINAL DISTRIBUTION OF ',A35)
C  IF NKG IS SELECTED ONLY THE ELECTRON DISTRIBUTION IS AVAILABLE
            ELSEIF ( FNKG ) THEN
              DEP(1)    = 0.D0
              CHAPAR(1) = 0.D0
              DO 931 J = 1,IALT(1)
                DEP(J+1)    = TLEV(J)
                CHAPAR(J+1) = SL(J)
 931          CONTINUE
              NSTP = IALT(1) + 1
              WRITE(MONIOU,8229) 'NKG ELECTRONS'
C  IF NONE IS SELECTED IT DOES NOT REALLY MAKE SENSE TO FIT
C  BUT LET'S TAKE THEN ALL CHARGED WHICH ARE MUONS AND HADRONS
            ELSE
              DO 932 J=0,NSTEP-LPCT1
                DEP(J+1)    = (J+LPCT1)*THSTEP
                CHAPAR(J+1) = PLONG(J+LPCT1,7)
 932          CONTINUE
              NSTP = NSTEP + 1 - LPCT1
              WRITE(MONIOU,8229) 'MUONS AND CHARGED HADRONS'
            ENDIF
            IF ( NSTP .GT. 6 ) THEN
C  THERE ARE MORE THAN 6 STEP VALUES, A FIT SHOULD BE POSSIBLE.
C  DO THE FIT: NPAR AND FPARAM GIVE THE NUMBER OF PARAMETERS USED
C  AND THE FINAL VALUES FOR THE PARAMETERS. CHISQ GIVES THE CHI**2/DOF
C  FOR THE FIT.
              CALL LONGFT(FPARAM,CHI2)
              WRITE(MONIOU,8230) FPARAM,CHI2,CHI2/SQRT(FPARAM(1))*100.D0
 8230         FORMAT(' PARAMETERS         = ',1P,6E12.4,0P/
     *               ' CHI**2/DOF         = ',F10.1/
     *               ' AV. DEVIATION IN % = ',F10.4)
C  STORE RESULT IN END EVENT BLOCK
              DO 933 K = 1,6
                EVTE(255+K) = FPARAM(K)
 933          CONTINUE
              EVTE(262) = CHI2
            ELSE
              WRITE(MONIOU,*) 'NO LONGI. FIT POSSIBLE, ',
     *          ' NSTP = ',NSTP,'  TOO SMALL.'
              DO 934 K = 1,6
                EVTE(255+K) = 0.
 934          CONTINUE
              EVTE(262) = 0.
            ENDIF
          ENDIF
        ENDIF
 
C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     Modified by C. Bigongiari 2001 Jan 16 
C 
C
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c   Saves statistics to sta###### file
C        call jcstadata(EVTH,EVTE,
C     +      NPROTO,NPROTB,NNEUTR,NNEUTB,NPHOTO,NELECT,NPOSIT,
C     +      NNU   ,NMUM  ,NMUP  ,NPI0  ,NPIM  ,NPIP  ,NK0L  ,
C     +      NK0S  ,NKMI  ,NKPL  ,NHYP  ,NDEUT ,NTRIT ,NALPHA,
C     +      NOTHER,IFINNU,IFINPI,IFINET,IFINKA,IFINHY,
C     +      CERELE,CERHAD,PLONG,LPCT1,NSTEP,THSTEP)
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  WRITE SHOWER END TO OUTPUT BUFFER
Cc        CALL TOBUF( EVTE,0 )
C        CALL TOBUF( EVTE,1 )
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C        IF ( LCERFI ) THEN
C          CALL OUTND2
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Cc          CALL TOBUFC( EVTE,0 )
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C        ENDIF

C
C  WRITE SHOWER END TO OUTPUT BUFFER
        CALL TOBUF( EVTE,0 )
        IF ( LCERFI ) THEN
          CALL OUTND2
          CALL TOBUFC( EVTE,0 )
        ENDIF

C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
 
        IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,*)
     *                 'CERENKOV PH. FROM ELECTRONS = ',SNGL(CERELE),
     *                 '  CERENKOV PH. FROM HADRONS = ',SNGL(CERHAD)
        CERELE = 0.D0
        CERHAD = 0.D0
        NRECER = 0
 
        IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,210) SHOWNO
 210    FORMAT(/'   END OF SHOWER NO ',I10)
 
        DO 19  J = 1,37
          JNBIN(J) = JNBIN(J) + INBIN(J)
          JPBIN(J) = JPBIN(J) + IPBIN(J)
          JKBIN(J) = JKBIN(J) + IKBIN(J)
          JHBIN(J) = JHBIN(J) + IHBIN(J)
  19    CONTINUE


  2   CONTINUE

C  END OF SHOWER LOOP
 
C-----------------------------------------------------------------------
 992  CONTINUE
 
C  RESET NUMBER OF SHOWERS TO CORRECT VALUE
      ISHW = I
 
      RUNE(3) = REAL(ISHW)
C  WRITE RUN END TO OUTPUT BUFFER AND FINISH OUTPUT
C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     Modified by C. Bigongiari 2001 Jan 16 
C
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Cc      CALL TOBUF ( RUNE,1 )
C      call jcendrun(rune)
Cc      IF ( LCERFI ) CALL TOBUFC( RUNE,1 )
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C- write Run End      
C

      CALL TOBUF ( RUNE,1 )
      IF ( LCERFI ) CALL TOBUFC( RUNE,1)

C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C   
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  TIME SINCE BEGINNING
c      ILEFTB = TIME()
      ILEFTB = 1
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      TDIFF  = ILEFTB - ILEFTA
 
C  MEAN VALUE FOR FIRST INTERACTION ALTITUDE (G/CM**2)
      IF ( ISHW .GT. 1 ) THEN
        CHISM2 = SQRT( ABS(CHISM2-CHISUM**2/ISHW) / (ISHW-1) )
        CHISUM = CHISUM / ISHW
      ELSE
        CHISM2 = 0.D0
      ENDIF
 
C  OUTPUTS FOR ALL SHOWERS
      WRITE(MONIOU,201) ISHW,TDIFF,TDIFF/ISHW,IRECOR,IRECOR/ISHW,
     *                   CHISUM,CHISM2
 201  FORMAT('1',10('='),' RUN SUMMARY ',56('=')//
     *  ' NUMBER OF GENERATED EVENTS = ',I10,/
     *  ' TOTAL TIME USED            = ',E10.3,' SEC'/
     *  ' TIME PER EVENT             = ',E10.3,' SEC'/
     *  ' TOTAL SPACE ON PATAPE USED = ',I10,' WORDS'/
     *  ' SPACE PER EVENT ON PATAPE  = ',I10,' WORDS'/
     *  ' AVERAGE HEIGHT OF 1ST INT. = ',F10.3,' +-',F10.3,' G/CM**2'/)
 
C  ENERGY - MULTIPLICITY MATRIX FOR ALL SHOWERS
      WRITE(MONIOU,209) (K,K=1,13),
     *  (J,(MULTOT(J,K),K=1,13),10**((J-4.)/3.),10**((J-3.)/3.),J=1,37),
     *   1,(INT(10**((K-1.)/3.)+1),K = 2,13),
     *   2,(INT(10**((K   )/3.)  ),K = 2,13)
 209  FORMAT(//' ENERGY - MULTIPLICITY MATRIX FOR ALL SHOWERS'/
     *       ' ENERGY RUNS VERTICALLY, MULTIPLICITY HORIZONTALLY'//,
     *       ' ',6X,5I9,3I8,5I7'   ENERGY RANGE (GEV)'/
     *       37(/' ',I4,1X,I10,4I9,3I8,5I7,1X,1P,2E10.1,0P)//
     *       ' MULT. ',5I9,3I8,5I7,4X,'LOWER BIN LIMIT'/
     *       ' RANGE ',5I9,3I8,5I7,4X,'UPPER BIN LIMIT')
 
C  GET MEAN OF ELASTICITY FOR ENERGY BINS
      DO 3377  J = 1,37
        NELMEA = 0
        DO 3378  K = 1,10
          NELMEA = NELMEA + IELDPA(J,K)
 3378   CONTINUE
        IF ( NELMEA .NE. 0 ) ELMEAA(J) = ELMEAA(J) / NELMEA
 3377 CONTINUE
 
C  PRINT ENERGY - ELASTICITY MATRIX FOR ALL SHOWERS
      WRITE(MONIOU,408) (K,K=1,10),  (J,(IELDPA(J,K),K=1,10),
     *  ELMEAA(J),10**((J-4.D0)/3.D0),10**((J-3.)/3.D0),J=1,37),
     *      ((K-1)*0.1D0,K=1,10),(K*0.1D0,K=1,10)
 408  FORMAT (//' ENERGY - ELASTICITY MATRIX FOR ALL SHOWERS'/
     *          ' ENERGY RUNS VERTICALLY, ELASTICITY HORIZONTALLY'//
     *          ' ',5X,10I9,'   MEAN EL.   ENERGY RANGE (GEV)'/
     *          37(/' ',I4,1X,10I9,2X,1P,E10.3,2E10.1,0P)//
     *          ' ELA. ',10F9.2,5X,'LOWER BIN LIMIT'/
     *          ' RANGE',10F9.2,5X,'UPPER BIN LIMIT')
 
      WRITE(MONIOU,204)
 204  FORMAT (//' INTERACTIONS PER KINETIC ENERGY INTERVAL FOR ALL ',
     *    'SHOWERS'//'   BIN    LOWER LIMIT    UPPER LIMIT     ',
     *    'NUCLEON     PIONS     KAONS S.BARYONS      TOTAL'/
     *    12X,'IN GEV',9X,'IN GEV',7X,
     *     '  EVENTS    EVENTS    EVENTS    EVENTS    '//)
      WRITE(MONIOU,207) (I,SABIN(I),SBBIN(I),JNBIN(I),JPBIN(I),JKBIN(I)
     *            ,JHBIN(I),JNBIN(I)+JPBIN(I)+JKBIN(I)+JHBIN(I),I=1,37)
 207  FORMAT(' ',I5,1P,2E15.4,0P,I12,3I10,I11)
 
      IF ( .NOT.GHEISH ) THEN
C  PRINT ELASTICITY STATISTICS
        WRITE(MONIOU,89) (I,(I-1)*.05,I*.05,
     *                   IELIS(I),IELHM(I),IELNU(I),IELPI(I),I = 1,20)
  89    FORMAT (//' ELASTICITY STATISTICS '//
     *          ' BIN   LOW  HIGH EDGE   FOR ISOBARS     HEAVY MESONS',
     *          '  SINGLE NUCLEONS        AND PIONS'/
     *         (' ',I3,'  ',F4.2,'  ',F4.2,'  ',4I17))
      ENDIF
 
C  CALCULATE MEAN VALUES AND STANDARD DEVIATIONS OF PARTICLE NUMBERS
      IF ( ISHW .GT. 1 ) THEN
        DO 879  K = 1,25
        DO 879  J = 1,NOBSLV
          MPART2(J,K) = SQRT( abs(MPART2(J,K)-MPARTO(J,K)**2/ISHW)
     *                                                  /(ISHW-1) )
          MPARTO(J,K) = MPARTO(J,K)/ISHW
 879    CONTINUE
      ELSE
        DO 880  K = 1,25
        DO 880  J = 1,NOBSLV
          MPART2(J,K) = 0.D0
 880    CONTINUE
      ENDIF
 
C  PRINT MEAN VALUES AND STANDARD DEVIATIONS OF PARTICLE NUMBERS
      IOBSLV = MIN( 3, NOBSLV )
      WRITE(MONIOU,854) (K,K=1,IOBSLV)
 854  FORMAT (/ ' AVERAGE NUMBER OF PARTICLES PER EVENT :'/
     *          ' FROM LEVEL NUMBER ', 3(10X,I10,10X) )
      WRITE(MONIOU,855) (OBSLEV(K),K=1,IOBSLV)
 855  FORMAT (  ' HEIGHT IN CM',1P,3(20X,E10.3)/)
 
      WRITE(MONIOU,778)'PROTONS     ',(MPROTO(K),MPROT2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'NEUTRONS    ',(MNEUTR(K),MNETR2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'PHOTONS     ',(MPHOTO(K),MPHOT2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'ELECTRONS   ',(MELECT(K),MELEC2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'POSITRONS   ',(MPOSIT(K),MPOSI2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'NEUTRINOS   ',(MNU   (K),MNU2  (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'MU -        ',(MMUM  (K),MMUM2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'MU +        ',(MMUP  (K),MMUP2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'PI 0        ',(MPI0  (K),MPI02 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'PI -        ',(MPIM  (K),MPIM2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'PI +        ',(MPIP  (K),MPIP2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'K0L         ',(MK0L  (K),MK0L2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'K0S         ',(MK0S  (K),MK0S2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'K -         ',(MKMI  (K),MKMI2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'K +         ',(MKPL  (K),MKPL2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'STR. BARYONS',(MHYP  (K),MHYP2 (K),K=1,IOBSLV)
      WRITE(MONIOU,778)'DEUTERONS   ',(MDEUT (K),MDEUT2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'TRITONS     ',(MTRIT (K),MTRIT2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'ALPHAS      ',(MALPHA(K),MALPH2(K),K=1,IOBSLV)
      WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=1,IOBSLV)
      WRITE(MONIOU,*)
 778  FORMAT(' NO OF ',A12,' = ',3(F13.1,' +-',F13.1,' '))
 
      IF ( NOBSLV .GT. 3 ) THEN
        IOBSLV = MIN( 6, NOBSLV )
        WRITE(MONIOU,854) (K,K=4,IOBSLV)
        WRITE(MONIOU,855) (OBSLEV(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'PROTONS     ',(MPROTO(K),MPROT2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'NEUTRONS    ',(MNEUTR(K),MNETR2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'PHOTONS     ',(MPHOTO(K),MPHOT2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'ELECTRONS   ',(MELECT(K),MELEC2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'POSITRONS   ',(MPOSIT(K),MPOSI2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'NEUTRINOS   ',(MNU   (K),MNU2  (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'MU -        ',(MMUM  (K),MMUM2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'MU +        ',(MMUP  (K),MMUP2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'PI 0        ',(MPI0  (K),MPI02 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'PI -        ',(MPIM  (K),MPIM2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'PI +        ',(MPIP  (K),MPIP2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'K0L         ',(MK0L  (K),MK0L2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'K0S         ',(MK0S  (K),MK0S2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'K -         ',(MKMI  (K),MKMI2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'K +         ',(MKPL  (K),MKPL2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'STR. BARYONS',(MHYP  (K),MHYP2 (K),K=4,IOBSLV)
        WRITE(MONIOU,778)'DEUTERONS   ',(MDEUT (K),MDEUT2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'TRITONS     ',(MTRIT (K),MTRIT2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'ALPHAS      ',(MALPHA(K),MALPH2(K),K=4,IOBSLV)
        WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=4,IOBSLV)
        WRITE(MONIOU,*)
 
        IF ( NOBSLV .GT. 6 ) THEN
          IOBSLV = MIN( 9, NOBSLV )
          WRITE(MONIOU,854) (K,K=7,IOBSLV)
          WRITE(MONIOU,855) (OBSLEV(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'PROTONS     ',(MPROTO(K),MPROT2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'NEUTRONS    ',(MNEUTR(K),MNETR2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'PHOTONS     ',(MPHOTO(K),MPHOT2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'ELECTRONS   ',(MELECT(K),MELEC2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'POSITRONS   ',(MPOSIT(K),MPOSI2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'NEUTRINOS   ',(MNU   (K),MNU2  (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'MU -        ',(MMUM  (K),MMUM2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'MU +        ',(MMUP  (K),MMUP2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'PI 0        ',(MPI0  (K),MPI02 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'PI -        ',(MPIM  (K),MPIM2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'PI +        ',(MPIP  (K),MPIP2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'K0L         ',(MK0L  (K),MK0L2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'K0S         ',(MK0S  (K),MK0S2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'K -         ',(MKMI  (K),MKMI2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'K +         ',(MKPL  (K),MKPL2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'STR. BARYONS',(MHYP  (K),MHYP2 (K),K=7,IOBSLV)
        WRITE(MONIOU,778)'DEUTERONS   ',(MDEUT (K),MDEUT2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'TRITONS     ',(MTRIT (K),MTRIT2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'ALPHAS      ',(MALPHA(K),MALPH2(K),K=7,IOBSLV)
        WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=7,IOBSLV)
          WRITE(MONIOU,*)
 
          IF ( NOBSLV .GT. 9 ) THEN
            IOBSLV = MIN( 10, NOBSLV )
            WRITE(MONIOU,854) (K,K=9,IOBSLV)
            WRITE(MONIOU,855) (OBSLEV(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'PROTONS     ',(MPROTO(K),MPROT2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'NEUTRONS    ',(MNEUTR(K),MNETR2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'PHOTONS     ',(MPHOTO(K),MPHOT2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'ELECTRONS   ',(MELECT(K),MELEC2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'POSITRONS   ',(MPOSIT(K),MPOSI2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'NEUTRINOS   ',(MNU   (K),MNU2  (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'MU -        ',(MMUM  (K),MMUM2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'MU +        ',(MMUP  (K),MMUP2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'PI 0        ',(MPI0  (K),MPI02 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'PI -        ',(MPIM  (K),MPIM2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'PI +        ',(MPIP  (K),MPIP2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'K0L         ',(MK0L  (K),MK0L2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'K0S         ',(MK0S  (K),MK0S2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'K -         ',(MKMI  (K),MKMI2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'K +         ',(MKPL  (K),MKPL2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'STR. BARYONS',(MHYP  (K),MHYP2 (K),K=9,IOBSLV)
        WRITE(MONIOU,778)'DEUTERONS   ',(MDEUT (K),MDEUT2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'TRITONS     ',(MTRIT (K),MTRIT2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'ALPHAS      ',(MALPHA(K),MALPH2(K),K=9,IOBSLV)
        WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=9,IOBSLV)
            WRITE(MONIOU,*)
          ENDIF
 
        ENDIF
      ENDIF
 
C  PRINT OUT NKG RESULT FOR ALL SHOWERS IF SELECTED
      IF ( FNKG ) CALL MITAGE
 
C  CALCULATE MEAN VALUES AND SIGMAS OF LONGITUDINAL DISTRIBUTION
      IF ( LLONGI ) THEN
        IF ( ISHW .GT. 1 ) THEN
          DO 790  K = 1,9
            DO 789  J = LPCT1,NSTEP
              SPLONG(J,K) = SQRT( abs(SPLONG(J,K)-APLONG(J,K)**2/ISHW)
     *                                                   /(ISHW-1) )
              APLONG(J,K) = APLONG(J,K)/ISHW
 789        CONTINUE
 790      CONTINUE
        ELSE
          DO 990  K = 1,9
            DO 989  J = LPCT1,NSTEP
              SPLONG(J,K) = 0.D0
 989        CONTINUE
 990      CONTINUE
        ENDIF
 
C  PRINT AVERAGE LONGITUDINAL DISTRIBUTIONS
        WRITE(MONIOU,911) THSTEP,
     *     'GAMMAS ','POSITRONS','ELECTRONS','MU-  ','MU+  ',
     *     (J*THSTEP,(APLONG(J,K),SPLONG(J,K),K=1,5),J=LPCT1,NSTEP)
 911    FORMAT(/' AVERAGE LONGITUDINAL DISTRIBUTION IN STEPS OF ',
     *      F5.0,' G/CM**2 '/' ',131('=')/
     *      ' DEPTH',8X,3(A10,16X),A9,15X,A9 //
     *     (' ',F5.0,2X,1P,E11.4,'+-',E11.4,0P,1X,F12.0,'+-',F11.0,
     *                1X,F12.0,'+-',E11.4,1X,F11.1,'+-',F10.1,
     *                1X,F11.1,'+-',F10.1 ))
        WRITE(MONIOU,912) THSTEP,
     *     'HADRONS','CHARGED','NUCLEI','CERENKOV',
     *     (J*THSTEP,(APLONG(J,K),SPLONG(J,K),K=6,9),J=LPCT1,NSTEP)
 912    FORMAT(/' AVERAGE LONGITUDINAL DISTRIBUTION IN STEPS OF ',
     *      F5.0,' G/CM**2 '/' ',115('=')/
     *      ' DEPTH',8X,A9,16X,A10,16X,A9,21X,A9 //
     *     (' ',F5.0,1X,F11.1,'+-',F11.1,1X,F12.0,'+-',F12.0,
     *                2X,F10.1,'+-',F10.1,1X,1P,E16.6,'+-',E16.6,0P))
      ENDIF
 
      IF ( FLGFIT ) THEN
C  PERFORM FIT TO THE LONGITUDINAL DISTRIBUTION OF ALL CHARGED PARTICLES
C  IF EGS IS SELECTED THIS IS THE DISTRIBUTION WHICH IS TO BE TAKEN
        IF ( FEGS ) THEN
          DO 730 J=0,NSTEP-LPCT1
            DEP(J+1)    = (J+LPCT1)*THSTEP
            CHAPAR(J+1) = APLONG(J+LPCT1,7)
 730      CONTINUE
          NSTP = NSTEP + 1 - LPCT1
          WRITE(MONIOU,8229) 'AVERAGE ALL CHARGED PARTICLES'
C  IF NKG IS SELECTED ONLY THE ELECTRON DISTRIBUTION IS AVAILABLE
        ELSEIF ( FNKG ) THEN
          DEP(1)    = 0.D0
          CHAPAR(1) = 0.D0
          DO 731 J = 1,IALT(1)
            DEP(J+1)    = TLEV(J)
            CHAPAR(J+1) = SEL(J)/ISHW
 731      CONTINUE
          NSTP = IALT(1) + 1
          WRITE(MONIOU,8229) 'AVERAGE NKG ELECTRONS'
C  IF NONE IS SELECTED IT DOES NOT REALLY MAKE SENSE TO FIT
C  BUT LET'S TAKE THEN ALL CHARGED WHICH ARE MUONS AND HADRONS
        ELSE
          DO 732 J=0,NSTEP-LPCT1
            DEP(J+1)    = (J+LPCT1)*THSTEP
            CHAPAR(J+1) = APLONG(J+LPCT1,7)
 732      CONTINUE
          NSTP = NSTEP + 1 - LPCT1
          WRITE(MONIOU,8229) 'AVERAGE MUONS AND CHARGED HADRONS'
        ENDIF
        IF ( NSTP .GT. 6 ) THEN
C  THERE ARE MORE THAN 6 STEP VALUES, A FIT SHOULD BE POSSIBLE.
C  DO THE FIT: NPAR AND FPARAM GIVE THE NUMBER OF PARAMETERS USED
C  AND THE FINAL VALUES FOR THE PARAMETERS. CHISQ GIVES THE CHI**2/DOF
C  FOR THE FIT.
          CALL LONGFT(FPARAM,CHI2)
          WRITE(MONIOU,8230) FPARAM,CHI2,CHI2/SQRT(FPARAM(1))*100.D0
        ELSE
          WRITE(MONIOU,*) 'NO LONGI. FIT POSSIBLE, ',
     *                      ' NSTP = ',NSTP,'  TOO SMALL.'
        ENDIF
      ENDIF
 
 
C  CONTROL PRINT OUTPUT OF CONSTANTS
      IF ( DEBUG ) THEN
        CALL STAEND
        WRITE(MDEBUG,*) 'MAIN  : STAEND CALLED'
      ENDIF
 
C
CBC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     Modified by C. Bigongiari 2001 Jan 16 
C
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      call jcenddata(runh,rune)
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
      WRITE(MONIOU,*)' '
      CALL PRTIME(TTIME)
      WRITE(MONIOU,101)
 101  FORMAT (/' ',10('='),' END OF RUN ',67('='))
 
C  CLOSE ALL OPEN UNITS
      IF ( MONIOU .NE. 6 ) CLOSE( MONIOU )
      IF ( MDEBUG .NE. 6 ) CLOSE( MDEBUG )
      CLOSE( EXST )

C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     Modified by C. Bigongiari 2001 Jan 16 
C
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Cc      CLOSE( PATAPE )
Cc      IF ( LCERFI ) CLOSE( CETAPE )
Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C- close files       
C
      
      CLOSE( PATAPE )
      IF ( LCERFI ) CLOSE( CETAPE )
C      CALL JCENDRUN( RUNE ) 

C
CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C    
 
      STOP
      END
