c=======================================================================
c
c  r e a d p a r t i c l e s . f
c  ----------------------------- 
c      write corsika shower data as readable ascii file, needs 
c      factor 3.5 of disk space as of binary particle data file.
c      This is impossible for 64bit simulations on all machines.
c      End of ascii print out may be incomplete because of
c      differencies in reading particle data of 32bit simulations
c      on 64bit machines. For 64bit corsika simulations only three
c      subblocks will be printed out (switch to 32bit simulations).
c-----------------------------------------------------------------------
c      gfortran -O0 -fbounds-check readparticles.f -o readparticles 
c      ifort -C -O0 -check bounds readparticles.f -o readparticles
c-----------------------------------------------------------------------
c           RUNH = 211285.2812500000000000;
c           EVTH = 217433.0781250000000000;
c           LONG =  52815.2968750000000000;
c           EVTE =   3397.3918457031250000;
c           RUNE =   3301.3325195312500000;
c-----------------------------------------------------------------------
c      input-files:
c           unit=*: number of showers and file name(s):
c           ------------------------------------------------
c                       1          'total_number_of_showers'
c                       1          'total_number_of_files'
c            '/lxdata/d2lx14/joe/DAT045216'
c                       1          '_showers_of_this_file_'
c           ------------------------------------------------
c           unit=3: current corsika particle data file.
c     output-files: 
c           unit=*: protocol output.
c           unit=9: ascii file named DATnnnnnn.ascithin or ....ascistnd.
c-----------------------------------------------------------------------
c                                     juergen.oehlschlaeger@kit.edu
c-----------------------------------------------------------------------
  
      program readparticles
 
      implicit double precision (a-h,o-z), integer (i-n) 

      parameter (lenthin=6552,lenstnd=5733) 

      character cout*120,crunh*4,cevte*4,qpatext(1:200)*19
      character cdata(100)*120,cdat*120,cblk*120
      double precision aatm(5),batm(5),catm(5)
      double precision parmas(0:101),phead(30)
      dimension qpartic(0:400),jpartic(0:201,32)
      real pdata(lenthin),qdata(936),prunh,pevte
      equivalence (crunh,prunh),(cevte,pevte)
      common /atmos/aatm,batm,catm
      common /integ/lobs,lsho,ishu,irec,isho,isub,jpartic
      common /utabl/pdata,parmas,phead,cpi180,c180pi,qpartic
      data crunh/'RUNH'/,cevte/'EVTE'/

c--initialize names of particles----------------------------------------
      data (qpatext(i),i=1,60)/
     c' gamma             ',' positron          ',' electron          ',
     c' _stck_in_         ',' muon+             ',' muon-             ',
     c' pi0               ',' pi+               ',' pi-               ',
     c' K0long            ',' K+                ',' K-                ',
     c' neutron           ',' proton            ',' anti proton       ',
     c' K0short           ','                   ',' Lambda            ',
     c' Sigma+            ',' Sigma0            ',' Sigma-            ',
     c' Xi0               ',' Xi-               ',' Omega-            ',
     c' anti neutron      ',' anti Lambda       ',' anti Sigma-       ',
     c' anti Sigma0       ',' anti Sigma+       ',' anti Xi0          ',
     c' anti Xi+          ',' anti Omega+       ',8*'                 ',
     c' Deuteron          ',' Tritium           ',' Helium3           ',
     c' Helium4           ',' addi muon+        ',' addi muon-        ',
     c 3*'                ',' omega             ',' rho0              ',
     c' rho+              ',' rho-              ',' Delta++           ',
     c' Delta+            ',' Delta0            ',' Delta-            ',
     c' anti Delta--      ',' anti Delta-       ',' anti Delta0       '/
      data (qpatext(i),i=61,129)/
     c' anti Delta+       ',' K*0               ',' K*+               ',
     c' K*-               ',' anti K*0          ',' electron neutrino ',
     c' anti elec neutrino',' muon neutrino     ',' anti muon neutrino',
     c'                   ',' eta=>2*gamma      ',' eta=>3*pi0        ',
     c' eta=>pi+pi-pi0    ',' eta=>pi+pi-gamma  ',' addi muon+        ',
     c' addi muon-        ','                   ',37*'                ',
     c'                   ',' D0                ',' D+                ',
     c' anti D-           ',' anti D0           ',' Ds+               ',
     c' anti Ds-          ',' eta c             ',' D*0               ',
     c' D*+               ',' anti D*-          ',' anti D*0          ',
     c' D*s+              ',' anti D*s-         ','                   '/
      data (qpatext(i),i=130,200)/
     c' J/psi             ',' tau +             ',' tau -             ',
     c' tau neutrino      ',' anti tau neutrino ','                   ',
     c'                   ',' Lambda c +        ',' Xi c +            ',
     c' Xi c 0            ',' Sigma c ++        ',' Sigma c +         ',
     c' Sigma c 0         ',' Xi c prime +      ',' Xi c prime 0      ',
     c' Omega c 0         ','                   ','                   ',
     c'                   ',' anti Lambda c -   ',' anti Xi c -       ',
     c' anti Xi c 0       ',' anti Sigma c --   ',' anti Sigma c -    ',
     c' anti Sigma c 0    ',' anti Xi c prime - ',' anti Xi c prime 0 ',
     c' anti Omega c 0    ','                   ','                   ',
     c'                   ',' Sigma c * ++      ',' Sigma c * +       ',
     c' Sigma c * 0       ',7*'                 ',' anti Sigma c * -- ',
     c' anti Sigma c * -  ',' anti Sigma c * 0  ',25*'                ',
     c' Cherenkov photon  ','                   '/
      qpatext(176) = ' Stnd              '
      qpatext(196) = ' Thin              '

c--initialize some quantities-------------------------------------------
      cblk='                                                  '
      cdat=cblk
      cpi = 4.d0 * atan(1.d0)
      cpi180 = cpi/180.d0
      c180pi = 180.d0/cpi
      isho = 0
      iend = 0
      irec = 0
      irun = 0
 
c--read run parameters including file names-----------------------------
      read(*,*,end=498,err=498) lsho
      read(*,*,end=498,err=498) nfil
      ! write(*,'(22x,''total number of showers:'',i7)') lsho
      ! write(*,'(22x,''total number of files  :'',i7)') nfil
      if ( nfil .gt. 100 ) goto 498
      do  ifil=1,nfil
         read(*,*,end=498,err=498) cdat
         read(*,*,end=498,err=498) nfsh
         cdata(ifil) = cdat
      enddo
      idat = index(cdata(1),'DAT')
      read(cdata(1)(idat+3:idat+8),'(i6)') irun
 
c----------one or more showers in big disk file-------------------------
      write(*,'(40('' _''),/,12x,''readparticles.f:'')')
      do  444  ifil=1,nfil
      if (ifil.gt.1) close(unit=3)
**    if (ifil.gt.1) close(unit=9)
      ilost = 0
      iswit = 0
      irec = 0
      idat = index(cdata(ifil),'DAT')
      ilen = index(cdata(ifil),' ') - 1
      write(*,'(/,8x,a)') cdata(ifil)(1:ilen)
* - - - - - - read data record with lenstnd words - - - -
      open(unit=3,file=cdata(ifil),status='old',
     +     form='unformatted',access='sequential')
      ishift = 2
      itype = 2
      read(unit=3,err=496,end=432) (pdata(i),i=1,lenstnd)
      close(unit=3)
      write(*,'(/,8x,''pdata(1:6)'',f11.2,i8.6,''.00'',f11.2,
     +   2f9.4,f12.2)') pdata(1),int(pdata(2)),(pdata(i),i=3,6)
      write(*,'(22x,''RUNH'',6x,''runnr'',7x,''date'',5x,''version'',
     +   4x,''nobsl'',4x,''obslev'')')
* - - - - - - check on reading 32bit simulation on 64bit machine:
      if ( 211285.2 .lt. pdata(1) .and.
     +     pdata(1) .lt. 211285.4 ) then
         ishift = 0
      elseif ( 211285.2 .lt. pdata(2) .and.
     +     pdata(2) .lt. 211285.4 ) then
         ishift = -1
         do  i=1,936-1
            pdata(i) = pdata(i+1)
         enddo
      elseif ( 2.0202 .lt. pdata(3) .and. pdata(3) .lt. 9.9999 ) then 
         ! check version number instead of testing (273+1) and (312+1).
         ishift = 1
         do  i=936,2,-1
            pdata(i) = pdata(i-1)
         enddo
         pdata(1) = prunh ! 211285.281 
      endif
* - - - - - - detect `standard` simulation instead of `thinning`:
      if ( 217433.0 .lt. pdata(273+1) .and.
     +                   pdata(273+1) .lt. 217433.2 ) then
         itype = 0
         cout = cdata(ifil)(idat:ilen)//'.ascistnd'
         lenrec = lenstnd
      elseif ( 217433.0 .lt. pdata(312+1) .and.
     +                       pdata(312+1) .lt. 217433.2 ) then
         itype = 1
         cout = cdata(ifil)(idat:ilen)//'.ascithin'
         lenrec = lenthin
      else
         write(*,*) '    ______________________________________________'
         write(*,*) '    ERROR: this corsika syntax should not occur.'
         goto 498
      endif
      lenblk = lenrec / 21
      lenpar = lenblk / 39 
      write(*,'(/,6x,''runnr ='',i7.6,5x,''log10(E) ='',f8.4,
     +   5x,''theta ='',f6.2,5x,''phi ='',f7.2)') irun,
     +   log10(1.d9*pdata(lenblk+4)),57.29578*pdata(lenblk+11),
     +   57.29578*pdata(lenblk+12)
      do  i=1,936 ! keep three subblocks including some particles.
         qdata(i) = pdata(i)
      enddo
      mthin = 0
      do  i=575,936,7
         if ( int(qdata(i)/1000.) .eq. 8888. ) mthin = mthin + 1
      enddo
      write(*,'(/,6x,''primary ='',i5,5x,''obslev(m) ='',f9.2,5x,
     +   ''date/vers ='',i9,''/'',f8.5)') int(pdata(3+lenblk)),
     +   pdata(48+lenblk)*1.d-2,20000000+ int(pdata(3)),pdata(4)
      if ( qdata(274)+qdata(547) .lt. 1. ) then
         write(*,'(/,12x,''`thinning run`'',/)')
      else
         if ( mthin .eq. 0 ) then
            write(*,'(/,12x''`standard run`'',/)')
         else
            write(*,'(/,12x,''`multithin run`'',/)')
         endif
      endif
      write(*,'(27(''  -''),/)')
      if ( ifil .eq. nfil ) write(*,'(30x,a)') cdata(ifil)(idat:ilen)
      if ( ishift .eq. -1 ) goto 443 
* - - - - - - read data record with lenrec words - - - -
**    optionally check `write(9,` for ascii tabular output.  
**    write(*,'(20x,a)') cout(1:ilen+9)
**    open(unit=9,file=cout,status='unknown',form='formatted')
      open(unit=3,file=cdata(ifil),status='old',
     +     form='unformatted',access='sequential')
  431 continue
      read(unit=3,err=496,end=432) (pdata(i),i=1,lenrec)
      irec = irec + 1
* - - - - - - check on original 64bit simulation:
      pdatone = pdata(1)
      if ( irec .gt. 1 .and. ( 3.21d-41 .lt. pdatone
     +                 .and. pdatone .lt. 3.23d-41 ) ) goto 443
      if (mod(irec,4000) .eq. 0) write(*,'(30x,''irec ='',i8)') irec
* - - - - - - - - - - -
      if ( ishift .eq. 1 ) then
* - - - - - - shift data corresp. to a 32bit simul. on a 64bit machine:
         do  i=lenrec,2,-1
            pdata(i) = pdata(i-1)
         enddo
         if ( irec .gt. 1 ) then
            pdata(1) = 177177.1
            ipart = pdata(lenpar+1) * 1.000001e-3
            if ( ipart .ge. 1 .and. ipart .le. 3 )
     +         pdata(1) = 1000. + mod(pdata(lenpar+1),1000.)
            if ( ipart .eq. 0 ) pdata(1) = 1091. ! fixed generation.
            if ( ipart .eq. 5 ) pdata(1) = 70000. + pdata(lenpar+1)
            if ( ipart .eq. 6 ) pdata(1) = 70000. + pdata(lenpar+1)
            if ( ipart .eq. 75 ) pdata(1) = 177177.1
            if ( ipart .eq. 76 ) pdata(1) = 177177.1
            if ( 3301.2 .lt. pdata(lenblk+1) .and.
     +         pdata(lenblk+1) .lt. 3301.4 ) pdata(1) = pevte ! 3397.39185
            if ( pdata(1) .eq. 177177.1 ) then
               ilost = ilost + 1
            else
               iswit = iswit + 1
            endif
         else ! irec = 1.
            pdata(1) = prunh ! 211285.281
         endif
      elseif ( ishift .eq. -1 ) then
* - - - - - - shift data corresp. to a 64bit simul. on a 32bit machine:
         do  i=1,lenrec-1
            pdata(i) = pdata(i+1)
         enddo
         pdata(lenrec) = -pdata(lenrec-lenpar) 
         ! check negative elements carefully.
      endif
* - - - - - - - - - - -
      call blwork(iend,lenrec)
      if ( iend .le. 0 ) goto 431
  432 continue
* - - - - - - end of corsika particle data file reached: 
**    write(*,*) '       isub =',isub,' (rune)    irec =',irec,
**   +      '    irwc =',819*irec
      write(*,'(30x,''irec ='',i8)') irec
      if ( 0 .lt. isub .and. isub .lt. 21 ) then
         do  i=1+isub*39,819
**          write(9,'(1p,8e14.6)') (0.,l=0,lenpar-1) 
         enddo
      endif
  443 continue
* - - - - - - distinguish contents of pdata vector:
      if ( ishift .eq. 1 ) then
         write(*,*) '    ______________________________________________'
         write(*,*) '    WARNING: will be better using a 32bit machine.'
         if ( ilost+iswit .gt. 0 ) then
            if ( ilost .ge. 0 )
     +         write(*,*) '    minimal lost particles =',ilost
            if ( iswit .ge. 0 )
     +         write(*,*) '    max switched particles =',iswit
         endif
      elseif ( ishift .eq. -1 ) then
**      write(*,'(5x,''irec=1.01'',f11.1,5f15.1)') (pdata(i),i=1,6)
**      write(*,'(10x,''1.02'',f11.1,5f15.1)') 
**   +       (pdata(i),i=lenblk+1,lenblk+6)
**      write(*,'(10x,''1.03'',f11.1,5f15.1)')
**   +       (pdata(i),i=lenblk*2+1,lenblk*2+5)
**      write(*,*) '             ',(pdata(i),i=lenblk*2+lenpar+1,
**   +                                         lenblk*2+lenpar+5)
        write(*,*) '    irec=1.01',(pdata(i),i=1,6)
        write(*,*) '         1.02',(pdata(i),i=lenblk+1,lenblk+6)
        write(*,*) '         1.03',(pdata(i),i=lenblk*2+1,lenblk*2+5)
        write(*,*) '             ',(pdata(i),i=lenblk*2+lenpar+1,
     +                                         lenblk*2+lenpar+5)
        write(*,*) '    _______________________________________________'
        write(*,*) '    fortran: impossible to read a 64bit simulation;'
        write(*,*) '        use new 32bit particle data files.'
        open(unit=9,file=cout,status='unknown',form='formatted')
        do  l=1,lenblk*3,lenpar ! to check content of first record.
           write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
        enddo
        close(unit=9)
      elseif ( ishift .eq. 0 ) then
c        write(*,*) '    ______________________________________________'
c        write(*,*) '    WARNING: last record of the ascii output file '
c        write(*,*) '        may be incompletely filled with zeroes as '
c        write(*,*) '        compared to the exactly last record of the'
c        write(*,*) '        corsika particle data file.'
cc       write(*,*) '    INFO: check record length 22940',
cc   +             ' of `standard` simulation'
cc       write(*,*) '        and other record length 26216',
cc   +             ' of `thinning` simulation.'
cc       write(*,*) '        qdata(1)___qdata(2)___qdata(3)',
cc   +              '___qdata(4)___qdata(5)'
cc       write(*,*) '      ',qdata(1),qdata(2),
cc   +                       qdata(3),qdata(4),qdata(5)
cc       write(*,*) '        qdata(1)__qdata(274)__qdata(313)',
cc   +              '__qdata(547)__qdata(625)'
cc       write(*,*) '      ',qdata(1),qdata(274),
cc   +                       qdata(313),qdata(547),qdata(625)
      endif
  444 continue
  445 continue
      close(unit=3)
      goto 499
 
c--end of data----------------------------------------------------------
  496 continue
      write(*,*) '       irec =',irec
      inull = 0
      do  l=1,5733,7
         if ( pdata(l) .eq. 0.d0 ) inull = inull + 1
      enddo
      do  i=1,5733,273
         write(*,'(f12.0,1p,3e14.6,''   i='',i4)') (pdata(l),l=i,i+3),i
      enddo
      write(*,*)
     +'    ____________________________________________________________'
      write(*,*) '    ERROR: simulation type of corsika is `standard`.'
      write(*,'(9x,''estimated number of particles'',f15.0,
     +   '' (as 1 sh.)'')') 741.d0+(819.d0*(irec-1))-inull  
      write(*,*) '        subblocks EVTE and RUNE missing!'
      goto 499 
  497 continue
      write(*,*) '       irec =',irec
      inull = 0
      do  l=1,6552,8
         if ( pdata(l) .eq. 0.d0 ) inull = inull + 1
      enddo
      do  i=1,6552,312
         write(*,'(f12.0,1p,3e14.6,''   i='',i4)') (pdata(l),l=i,i+3),i
      enddo
      write(*,*) '    _________________________________________________'
      write(*,*) '    ERROR: simulation type of corsika is `thinning`.'
      write(*,'(9x,''estimated number of particles'',f15.0,
     +   '' (as 1 sh.)'')') 741.d0+(819.d0*(irec-1))-inull  
      write(*,*) '        subblocks EVTE and RUNE missing!'
      goto 499 
  498 continue
      write(*,*) '    _________________________________________________'
      write(*,*) '    ERROR: missing some input parameters.'
  499 continue
c--print out tabulars---------------------------------------------------
      npot = 16 ! logarithmic bins per decade.
      fpot = 10.d0**(1.d0/npot) ! factor per bin.
      if ( npot .lt. 0 ) then
      write(*,'(10x,''photons'')')
      write(*,'(10x,''-------'')')
      do  jp=-4,4
         ja = npot * (jp+4)
         write(*,'(/,1p,e7.0,1x,0p,15f8.5)') 10.d0**jp,(fpot**i,i=1,15)
         write(*,'(16i8)') (jpartic(i,1),i=ja,ja+npot-1)
      enddo
      write(*,'(10x,''positrons'')')
      write(*,'(10x,''---------'')')
      do  jp=-4,4
         ja = npot * (jp+4)
         write(*,'(/,1p,e7.0,1x,0p,15f8.5)') 10.d0**jp,(fpot**i,i=1,15)
         write(*,'(16i8)') (jpartic(i,2),i=ja,ja+npot-1)
      enddo
      write(*,'(10x,''electrons'')')
      write(*,'(10x,''---------'')')
      do  jp=-4,4
         ja = npot * (jp+4)
         write(*,'(/,1p,e7.0,1x,0p,15f8.5)') 10.d0**jp,(fpot**i,i=1,15)
         write(*,'(16i8)') (jpartic(i,3),i=ja,ja+npot-1)
      enddo
      endif
* - - - - - - print total counts of particles (also with weights):
      write(*,'(2x)')
      write(*,'(/,6x,''code name'',18x,''particles '',a)')
     +   cdata(nfil)(idat:ilen)
      write(*,'(6x,''---- ----'',18x,''--------- '',16(''-''))')
      do  l=1,50 ! 200
         if ( qpartic(l) .gt. 0. ) write(*,'(i10,a19,f13.0)')
     +      l,qpatext(l),qpartic(l)
      enddo
      write(*,'(6x,''---- ---- ----'',13x,''--------- '')')
      if ( lenblk .eq. 312 ) then ! prints particles incl. weights.
        do  l=1,50 ! 200
          if ( qpartic(l+200) .gt. 0. ) write(*,'(i10,a19,f14.1)')
     +       l,qpatext(l),qpartic(l+200)
        enddo
        write(*,'(6x,''---- ---- ----'',13x,''--------- '')')
      endif
      write(*,'(2x)')
      stop
      end
c=======================================================================
c
c     block data initialization
c
c=======================================================================
      block data blinit

      implicit double precision (a-h,o-z), integer (i-n)

      parameter (lenthin=6552)

      double precision aatm(5),batm(5),catm(5)
      double precision parmas(0:101),phead(30)
      dimension qpartic(0:400),jpartic(0:201,32)
      real pdata(lenthin)
      common /atmos/aatm,batm,catm
      common /integ/lobs,lsho,ishu,irec,isho,isub,jpartic
      common /utabl/pdata,parmas,phead,cpi180,c180pi,qpartic
      data aatm /-186.5562d0,  -94.919d0, 0.61289d0, 0.d0, .01128292d0/
      data batm /1222.6562d0,1144.9069d0, 1305.5948d0, 540.1778d0,0.d0/
      data catm /994186.38d0,878153.55d0,636143.04d0,772170.16d0,1.d-9/
      data parmas/
     + 2*0.0000d0 , 0.000511d0 , 0.000511d0 , 0.000000d0 , 0.105658d0,
     + 0.105658d0 , 0.134973d0 , 0.139568d0 , 0.139568d0 , 0.497671d0,
     + 0.493646d0 , 0.493646d0 , 0.939566d0 , 0.938272d0 , 0.938272d0,
     + 0.497671d0 , 0.5488d0   , 1.11563d0  , 1.18937d0  , 1.192550d0,
     + 1.19743d0  , 1.3149d0   , 1.32132d0  , 1.67243d0  , 0.939566d0,
     + 1.11563d0  , 1.18937d0  , 1.19255d0  , 1.19743d0  , 1.3149d0,
     + 1.32132d0  , 1.67243d0  , 1.7841d0   , 1.7841d0   , 1.8693d0,
     + 1.8693d0   , 1.8645d0   , 1.8645d0   , 1.9693d0   , 1.9693d0,
     + 2.2852d0   , 80.6d0     , 80.6d0     , 91.161d0   , 1.8770d0,
     + 2.817d0    , 3.755d0    , 0.0d0      , 0.0d0      , 0.0000d0,
     + 0.7669d0   , 0.7681d0   , 0.7681d0   , 1.2309d0   , 1.2323d0,
     + 1.2336d0   , 1.2349d0   , 1.2309d0   , 1.2323d0   , 1.2336d0,
     + 1.2349d0   , 0.0d0      , 0.0d0      , 0.0d0      , 0.0000d0,
     + 0.0d0      , 0.0d0      , 0.0d0      , 0.0d0      , 0.0000d0,
     + 0.5488d0   , 0.5488d0   , 0.5488d0   , 0.5488d0   , 0.105658d0,
     + 0.105658d0 , 0.0d0      , 0.0d0      , 0.0d0      , 22*0.0d0/
      end
c=======================================================================
c
c     analyze the contents of all 21 subblocks in a record
c
c=======================================================================
c
c        primary    energy    runnumb   simnrsh     #hadr    #muons
c        #gammas     #elec   #elecnkg    obslv    theta    phi
c         h1wkm    h1wgcmq  ( emymax    ehasum    ehamax    nha )
c
c=======================================================================
 
      subroutine blwork(iend,lenrec)

      implicit double precision (a-h,o-z), integer (i-n)

      parameter (lenthin=6552)
 
      double precision aatm(5),batm(5),catm(5)
      double precision parmas(0:101),phead(30)
      dimension qpartic(0:400),jpartic(0:201,32)
      real pdata(lenthin)
      common /integ/lobs,lsho,ishu,irec,isho,isub,jpartic
      common /utabl/pdata,parmas,phead,cpi180,c180pi,qpartic
      isub = 0
      lenblk = lenrec / 21
      lenpar = lenblk / 39

c-----------loop over subblocks-----------------------------------------
      do  948  lia=1,lenrec,lenblk
      isub = isub + 1
      if ( 211285.2 .le. pdata(lia).and.pdata(lia) .le. 211285.4 ) then
c----------------subblock run header------------------------------------
         lobs = int(pdata(lia+4))
         if (isho.le.0) then
            phead(6) = lobs
            do  908  i=1,lobs
               phead(6+i) = 1.e-2 * pdata(lia+4+i)
  908       continue
            if (lobs.ge.1.and.lobs.le.10) then
               if (lobs.gt.1) write(*,
     +            '(6x,''observation level:'',f10.2,'' meter'')')
     +            phead(6+lobs)
            else
               write(*,'(6x,''nr of observation level undefined.'')')
               iend = iend + 9876543
               goto 949
            endif
         endif
c - - - - - - - - - - - observation level (meter) - - - - - - - - - - -
         phead(17) = 1.e-2 * pdata(lia+lobs)
c - - - - - - - - - - - run number  - - - - - - - - - - - - - - - - - -
         phead(18) = pdata(lia+1)
         do  909  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
  909    continue
      elseif ( 217433.0 .le. pdata(lia).and.pdata(lia) .le. 217433.2 )
     +   then
c----------------subblock event header----------------------------------
         if (iend.lt.0) goto 948
         isho = isho + 1
c- - - - - - - - - - simulated shower number - - - - - - - - - - - - -
         phead(1) = pdata(lia+1)
c - - - - - - - - - - height of first interaction (km)  - - - - - - - -
         phead(2) = 1.e-5 * pdata(lia+6)
c - - - - - - - - - - height of first interaction (grams/cm2)- - - - -
         pdata(lia+4) = thickgr(dble(pdata(lia+6)))
         phead(3) = pdata(lia+4)
c - - - - - - - - - - primary particle code - - - - - - - - - - - - - -
         phead(4) = pdata(lia+2)
c - - - - - - - - - - primary particle energy in gev  - - - - - - - - -
         phead(5) = pdata(lia+3)
c - - - - - - - - - - - phi angle this shower in radian - - - - - - - -
         phead(19) = pdata(lia+11)
c - - - - - - - - - - theta angle this shower in radian - - - - - - - -
         phead(20) = pdata(lia+10)
         timev0 = ( phead(2) * 1.e+5 - phead(6+lobs)*100. ) / 
     /      cos(phead(20)) / 29.9792458e0
         pdata(lia+37) = timev0
         do  912  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
  912    continue
c----------------subblock longi information-----------------------------
      elseif (52815.2 .le. pdata(lia).and.pdata(lia) .le. 52815.4) then
         do  915  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
  915    continue
c----------------subblock event end-------------------------------------
      elseif ( 3397.3 .le. pdata(lia).and.pdata(lia) .le. 3397.5 ) then 
         if (iend.lt.0) then ! ignore first showers.
            iend = iend + 1
            isho = iend
            goto 948
         endif
         do  913  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
  913    continue
**       write(*,*) '       isub =',isub,' (evte)    irec =',irec,
**   +      '    isho =',isho
c----------------subblock run end---------------------------------------
      elseif ( 3301.3 .le. pdata(lia).and.pdata(lia) .le. 3301.5 ) then
         iend = iend + iend + isho
         do  914  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
  914    continue
         goto 949
      else
c-----------subblock with particle data---------------------------------
         do  917  l=lia,lia+lenblk-1,lenpar
**          write(9,'(1p,8e14.6)') (pdata(i+l),i=0,lenpar-1)
            icode = int(1.e-3*(pdata(l)*1.0000001))
            if ( 0 .lt. icode .and. icode .le. 32 ) then
               equad = parmas(icode)*parmas(icode)
     +               + pdata(l+1)*pdata(l+1)
     +               + pdata(l+2)*pdata(l+2) + pdata(l+3)*pdata(l+3)
               if ( equad .gt. 0 ) then 
                  epart = sqrt(equad) 
                  eplog = log10(epart) + 4. ! multiply by 1.e4
                  jlog = eplog * 16. ! 16 log bins per decade
                  if ( jlog .lt. 0 ) jlog = 0
                  if ( jlog .gt. 201 ) jlog = 201 
**                if ( jlog .eq. 15 .and. icode .eq. 2 ) 
**   +               write(*,*) sngl(epart), sngl(eplog), icode, jlog
                  jpartic(jlog,icode) = jpartic(jlog,icode) + 1 
               endif  
            endif
            if ( icode .gt. 50 ) then
               if ( icode .eq. 75 ) then
                  icode = 45
               else if ( icode .eq. 76 ) then
                  icode = 46
               else if ( icode .eq. 201 ) then
                  icode = 41
               else if ( icode .eq. 301 ) then
                  icode = 42
               else if ( icode .eq. 302 ) then
                  icode = 43
               else if ( icode .eq. 402 ) then
                  icode = 44
               else if ( icode .gt. 200 ) then
                  icode = 4
               else
                  if ( icode .gt. 96 ) then
                     icode = 49 ! ignore 85,86,95,96.
                  endif
               endif
            endif
            if ( icode .gt. 0 ) then
               qpartic(icode) = qpartic(icode) + 1.
               qpartic(icode+100) = qpartic(icode+100)+pdata(l+lenpar-1)
            endif
  917    continue
      endif
      if (iend.gt.0) goto 949
  948 continue
 
c-----------all records analyzed or end of showers reached--------------
  949 continue
      return
      end
c=======================================================================
 
      double precision function thickgr( arg )
 
c-----------------------------------------------------------------------
c  calculates thickness (g/cm**2) of atmosphere depending on height (cm)
c  argument:    arg    = height in cm
c-----------------------------------------------------------------------
 
      double precision aatm(5),batm(5),catm(5),arg
      common /atmos/aatm,batm,catm
 
      if     ( arg .lt. 4.d5 ) then
         thickgr = aatm(1) + batm(1) * exp( -arg/catm(1) )
      elseif ( arg .lt. 1.d6 ) then
         thickgr = aatm(2) + batm(2) * exp( -arg/catm(2) )
      elseif ( arg .lt. 4.d6 ) then
         thickgr = aatm(3) + batm(3) * exp( -arg/catm(3) )
      elseif ( arg .lt. 1.d7 ) then
         thickgr = aatm(4) + batm(4) * exp( -arg/catm(4) )
      else
         thickgr = aatm(5) - arg * catm(5)
      endif
 
      return
      end
