cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc     EqHaz03 Ver. a (2011/11/22)
ccc
ccc	The programs Haz03 is developed to read GM catalog files generated by th 
ccc	e program 'EqHaz2' (at most 10 EqHaz2 outputs) and makes the following p
ccc   calculate the percent roducts:
ccc     1-Lists all records that exceed required rturn period motion level.
ccc     2-Calculates magnitude-distance contributions in the motions of above 
ccc       list (deaggregation)(2012/02/26).
ccc     3-Lists the records of peak motions in consecutive windows of required 
ccc       durations (eg. 475yr)(2012/02/26).
ccc     4-Sorts the pick motions of given windows and prepares their fractile 
ccc       table(2012/02/26).
ccc     5-Calculates magnitude-distance contributions in peak motions of windo
ccc       ws (second type of deaggregation)(2012/02/26).
ccc
ccc     * This program is developed using Intel Visual Fortran and the compati
ccc		  bility of this code with other compilers is not tested.
ccc     * The program asks for a control file. The control file should carry t
ccc       he information of ground motion catalog file(s) generated by EqHaz2,
ccc       flags for output types, return periods (in years) of motion for whic
ccc       h the percentiles are calculated, and output file stem name.
ccc     * The flags in input file make the program flexible for generating the 
ccc       type of required output. The flags correspond to sorted output and f
ccc       ractile output for user defined periods respectively. 
ccc       If one does not need any type of output, one has to set correspondin
ccc       g flag value equal to zero. For example setting all flags equal to 1
ccc       will result in generation of both types of files.
ccc     * For the calculation of pick motion percentiles for a given time peri
ccc       od, the program divides the whole duration of catalog to subwindows
ccc       of desired period and finds the maximum motion in every subwindow, t
ccc       hen it sorts the maximum motions in descending order and assigns a s
ccc       equence number to them. If there are N sub-windows the percentile as
ccc       sociated with i-th pick motion will be equal to 100*(N-i)/N.
ccc     * We chose the following convention for naming the output files:
ccc        - The sorted catalog file is called: stem_S?C?_sort-?.out   where 
ccc          "S?" specifies the site number in EqHaz2 input, "C?" specifies th
ccc          e catalog number in EqHaz2 input, and "-?" specifies GM catalog n
ccc          umber in EqHaz3 input.
ccc        - The file containing the fractiles of pick motions in a chosen per
ccc          iod of time: stem_S?C?_fract????-?.out where "S?", "C?", "-?" are
ccc          same as the sorted output, and "????" specifies the subwindow len
ccc          gth in years. 
ccc     * Both "out" files and "tab" files from EqHaz2 (type "a" and type "b")
ccc       can be processed with this program.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc     EqHaz03 Ver. a (2012/02/26)
ccc
ccc   1-A correcttion in code is applied for making efficient use of a loop an
ccc     d making the program much faster.
ccc   2-Deaggregation calculations for given levels of motions is added. The u
ccc     ser adds the minimum, maximum, and incremental magnitudes and distance
ccc     s in input, for which bin sizes and corresponding contributions are ca
ccc     lculated. These files carry the term "deagg" in their file names.
ccc   3-Another type of deaggregation calculation feature is added in which th
ccc     e deaggregation tables are prepared for extreme motions of subwindows 
ccc     of user defined durations (eg. 475 year). These extra information are
ccc     added in "fract" files.
ccc   4-The complete list of records of extreme motions of subwindows are repo
ccc     rted in "fract" files.
ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc     EqHaz03 Ver. 2012c (2012/06/26)
ccc
ccc   1-The code is changed and it is consistent with outpus from EqHaz2_2012c
ccc   2-This version has option for calculating mean hazard and hazard curves
ccc     for subcatalogs. The flag controling this feature accepts: 0 for not c
ccc     alculating hazard values for subcatalogs, 1 for calculating subcatalog
ccc     hazards, and 2 checking EqHaz2 output for mixing of epistemic and alea
ccc     tory uncertainties in which subcatalog hazards will be calculated if t
ccc     he uncertainties are not mixed.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



      module datastruct
          structure /catalog/
              character (len=50)                :: cat_name                     !Catalog file name
              character (len=300),dimension(100):: cat_head                     !Catalog file header lines
              real                              :: wgt                          !Weight of this catalog in hazard. The sum of weights of all catalogss should be 1.
              integer                           :: site_ser, cat_ser            !Site serial number and catalog serial number stored in "out"-files or "tab"-files generated by EqHaz2
              integer                           :: nline_hdr                    !Number of header lines in catalog file
              integer                           :: n_of_subcats,n_of_yrs        !Number of sub-catalogs and duration of each sub-catalog
              integer                           :: n_rec                        !Number of records in a catalog file
              integer, dimension(1000000)       :: subcat_no                    !Subcatalog number of a record
              integer                           :: isort,isubcat,iyears         !Flags for calculation of sorted GM, percentiles of maxs in subcatalog, and percentiles of maxs in in given years 
              integer                           :: ihaz_of_subcat               !Flag for calculation of hazard for subcats
              integer                           :: epist                        !Epistemic flag in EqHaz2 output suggesting mixing of epistemic and aleatory uncertainties (Yes=1 and No=0)
              real, dimension(1000000,2)        :: time                         !Time of a record in catalog
              real, dimension(1000000,3)        :: loc                          !Latitude, longitude, and depth of a record in catalog
              real, dimension(1000000)          :: mag                          !Magnitude of a record in catalog
              integer, dimension(1000000)       :: zone_no                      !Zone number of a record in catalog
              real, dimension(1000000,2)        :: N0beta                       !N0 and beta of a record in catalog
              real, dimension(1000000)          :: Mmax                         !Mmax of a record in catalog
              real, dimension(1000000,4)        :: distance                     !Distances of an arthquake in a record from site (Repi, Rhyp, Rcd, Rjb respectively)
              real, dimension(1000000,2)        :: GM                           !Ground motion parameter and applied multiple of standard deviations
              integer, dimension(1000000)       :: GM_no                        !GM model applied
              integer                           :: all_read                     !Flag specifying if catalog is read completely in arrays
              
              real, dimension(100000,4)         :: haz                          !This matrix carries the rate of exceedance for each percentil
              real, dimension(1000000)          :: indices                      !This array carries the indices used for later re-ordering
              real, dimension(1000000)          :: temp                         !This array carries the values of parameter for which sorting will be used
          end structure

          record /catalog/ cat(10)                                              !This line creates an array of 3 catalog structures
          
          real, dimension(100)                  :: yrs_for_percentile           !Return periods (in years) of motion for which the percentiles are calculated
          integer                               :: yrs_for_percentile_no        !Number of return periods for which the percentiles are calculated
          integer                               :: n_of_cats
          character (len=80)                    :: input,out_stem
          real                                  :: deagg_mmin,deagg_mmax
          real                                  :: deagg_mint, deagg_m
          real                                  :: deagg_dmin,deagg_dmax
          real                                  :: deagg_dint, deagg_d
          real                                  :: deagg_val
          real, dimension(1000000)              :: sigw
          integer                               :: sigw_n
          real                                  :: sigw_avg, sigw_std

      end module


      program pp
      use datastruct
      real, allocatable                         :: fractiles(:)
      integer                                   :: i
      character (len=5)                         :: c5
      
      write(*,*)'Enter input file name: '
      read(*,*)input
      write(*,*)'Reading control and input catalog file(s)'
      call read_input()
      do i=1, n_of_cats
          call find_sorted_indices(i)
          
          do j=1,yrs_for_percentile_no
              call make_deagg_output(i, yrs_for_percentile(j), 'out')
          enddo
          if(cat(i).isort /= 0)then
              write(*,*)'Preparing sorted file #',i
              call make_sorted_output_w(i, 'out')
          endif
          if(cat(i).iyears /= 0)then
              do j=1,yrs_for_percentile_no
                  write(*,*)'Preparing fractile file #',i,'for',
     1                       yrs_for_percentile(j)
                  call calc_haz(i,yrs_for_percentile(j),'out')
              enddo
          endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   Added block for calculating hazard of blocks (mainly intended for fractile
!!!   calculations).
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   This subsection generates mean hazard values for each subcatalog and lists
!!!   them in a fractile table (two columns). This section will be run if the va
!!!   lue of ihaz_of_subcat flag is set to 1 and skipped if the value of flag is
!!!   set to zero. If the value of ihaz_of_subcat is set equal to 2, EqHaz3 will
!!!   look for flag "epist" in EqHaz2 output and if the value of epist is equal
!!!   to 0 it means the epistemic and aleatory uncertainties are not mixed and E
!!!   qHaz3 will calculate mean hazaed for each subcatalog (will skip otherwise)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          if(cat(i).ihaz_of_subcat == 2)then
              cat(i).ihaz_of_subcat=1-cat(i).epist
          endif
          if(cat(i).ihaz_of_subcat /= 0)then
              open(unit=21,file='Fractiles.txt',status='unknown')
              write(21,'(a)')(cat(i).cat_head(1:cat(i).nline_hdr-1))
              write(21,'(a)')'Fractile tables:'
              do j=1,yrs_for_percentile_no
                  write(*,*)'Calculating subcatalog hazard for',
     1                       yrs_for_percentile(j)
                  write(21,*)
                  write(21,'(f9.1,x,a)')
     1                      yrs_for_percentile(j),'year period'
                  write(21,'(a)')'     %ile     |   GM level   '
                  n_last=0
                  allocate (fractiles(cat(i).n_of_subcats))
                  do k=1,cat(i).n_of_subcats
                      n_first=n_last+1
                      call find_subcat_indices(i,k,n_first,n_last)
                      call find_hazvalue(i,j,k,n_first,n_last,
     1                                   fractiles(k))
                  enddo
                  call sort(cat(i).n_of_subcats,fractiles)
                  write(21,'(1x,g13.5,2x,g13.5)')
     1                 ((((real(l-1)/real(cat(i).n_of_subcats))),
     2                 fractiles(l)),l=1,cat(i).n_of_subcats)
                  deallocate (fractiles)
              enddo
          endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   This subsection generates mean curves values for each subcatalog and list 
!!!   them in separate files (single column). This section will be run if the va
!!!   lue of ihaz_of_subcat flag is set to 1 and skipped if the value of flag is
!!!   set to zero. If the value of ihaz_of_subcat is set equal to 2, EqHaz3 will
!!!   look for flag "epist" in EqHaz2 output and if the value of epist is equal
!!!   to 0 it means the epistemic and aleatory uncertainties are not mixed and E
!!!   qHaz3 will calculate mean hazaed curves for each subcatalog (and will skip 
!!!   otherwise).
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          if(cat(i).ihaz_of_subcat /= 0)then
              write(*,*)'Preparing hazard curves for subcatalogs'
              n_last=0
              open(unit=23,file='hazcurves_all.txt',status='unknown')
              write(23,'(a)')'Hazard curves of all subcatalogs. '//
     1                       'First row of numbers are ground motions.'
              write(23,'(a)')'Second row of numbers shows rates of '//
     1                       'exceedances of first subcatalog and so on'
              do k=1,cat(i).n_of_subcats
                  write(c5,'(i5.5)')k
                  open(unit=22,file='hazcurve_subcat'//c5//'.txt',
     1                       status='unknown')
                  write(22,'(a,i4)')'Hazard curve for subcatalog #',k
                  write(22,'(a)')
     1                     '      GM level      | Rate of exceedance '
                  n_first=n_last+1
                  call find_subcat_indices(i,k,n_first,n_last)
                  call find_hazcurve(i,k,n_first,n_last)
                  close(unit=22)
              enddo
              close(unit=23)
          endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      enddo
      end



c
c     This subroutine finds the maximum of a subcatalog and makes an output
c     file with sorted values used for making GM rate of exceedance tables.
c
      subroutine calc_haz(n_c,n_year,c3)
      use datastruct
      real, dimension (100000)                  :: tmp1, tmp2
      real                                      :: n_year
      integer                                   :: n_c, n, n_
      character (len=2)                         :: c2
      character (len=3)                         :: c3, c3_
      character (len=5)                         :: c5
      character (len=7)                         :: c7
      cat(n_c).haz(:,3:4)=0.
      n=int(cat(n_c).n_of_subcats*cat(n_c).n_of_yrs/n_year)
      i=1
      do j=1,cat(n_c).n_rec
          if(cat(n_c).time(j,2) >  real(i*n_year))i=i+1
          if(cat(n_c).time(j,2) >= real((i-1)*n_year).and. 
     1       cat(n_c).time(j,2) <  real( i *  n_year).and. 
     2       cat(n_c).GM(j,1)   >= cat(n_c).haz(i,3))then
                  cat(n_c).haz(i,3)=cat(n_c).GM(j,1)
                  cat(n_c).haz(i,4)=real(j)
          endif
      enddo
      
      do i=1,n
          tmp1(i)=cat(n_c).haz(i,3)
          tmp2(i)=cat(n_c).haz(i,4)
      enddo
      call sort2(n,tmp1(1:n),tmp2(1:n))
      do i=1,n
          cat(n_c).haz(i,3)=tmp1(i)
          cat(n_c).haz(i,4)=tmp2(i)
      enddo
      
      write(c5,'(a1,i4)')'S',cat(n_c).site_ser
      if(cat(n_c).site_ser<1000)then
          write(c5,'(a1,i3,a1)')'S',cat(n_c).site_ser,' '
      endif
      if(cat(n_c).site_ser<100)then
          write(c5,'(a1,i2,a2)')'S',cat(n_c).site_ser,'  '
      endif
      if(cat(n_c).site_ser<10)then
          write(c5,'(a1,i1,a3)')'S',cat(n_c).site_ser,'   '
      endif
      write(c2,'(a1,i1)')'C',cat(n_c).cat_ser
      write(c7,'(i7)')int(n_year)
      if(n_year<1000000)write(c7,'(i6,a1)')int(n_year),' '
      if(n_year<100000)write(c7,'(i5,a2)')int(n_year),'  '
      if(n_year<10000)write(c7,'(i4,a3)')int(n_year),'   '
      if(n_year<1000)write(c7,'(i3,a4)')int(n_year),'    '
      if(n_year<100)write(c7,'(i2,a5)')int(n_year),'     '
      if(n_year<10)write(c7,'(i1,a6)')int(n_year),'      '
      write(c3_,'(a1,i2)')'-',n_c
      if(n_c<10)write(c3_,'(a1,i1,a1)')'-',n_c,' '
c
c     block3
c      
      open(unit=11,file=trim(out_stem)//'_'//trim(c5)//c2//'_fract'//
     1                  trim(c7)//'_b3'//trim(c3_)//'.'//trim(c3),
     2                  status='unknown')
      write(11,'(a39,f9.2,a14)')
     1 'List of fractiles of maximum motions in',n_year,' year windows.'
      write(11,'(a42,f9.2,a30)')
     1 'The maximum motions are the peak values of', n_year, 
     2 ' year consecutive time blocks.'
      write(11,53)trim(cat(n_c).cat_name)
53    format('This table is extracted from file:',a50)
      write(11,54)
54    format('With the following header:')
      do i=1,cat(n_c).nline_hdr-1
          write(11,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(11,*)
      write(11,55)
55    format('Block 1:',/,'     %ile     |   GM level   ')
      do i=n,1,-1
          write(11,56)real((i-1.)/n)*100.,
     1                cat(n_c).haz(i,3)
56        format(1x,g13.5,2x,g13.5)
      enddo
      close(unit=11)
c
c     block2
c      
      open(unit=11,file=trim(out_stem)//'_'//trim(c5)//c2//'_fract'//
     1                  trim(c7)//'_b2'//trim(c3_)//'.'//trim(c3),
     2                  status='unknown')
      write(11,'(a37,f9.2,a14)')
     1   'List of records of maximum motions in',n_year,' year windows.'
      write(11,'(a42,f9.2,a30)')
     1 'The maximum motions are the peak values of', n_year, 
     2 ' year consecutive time blocks.'
      write(11,53)trim(cat(n_c).cat_name)
      write(11,54)
      do i=1,cat(n_c).nline_hdr-1
          write(11,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(11,*)
      write(11,57)cat(n_c).cat_head(cat(n_c).nline_hdr)
57    format(a300)
      do i=n,1,-1
          write(11,58) cat(n_c).subcat_no(int(cat(n_c).haz(i,4))),
     1                 cat(n_c).time(int(cat(n_c).haz(i,4)),1),
     2                 cat(n_c).loc(int(cat(n_c).haz(i,4)),1:3),
     3                 cat(n_c).mag(int(cat(n_c).haz(i,4))),
     4                 cat(n_c).zone_no(int(cat(n_c).haz(i,4))),
     5                 cat(n_c).N0beta(int(cat(n_c).haz(i,4)),1:2),
     6                 cat(n_c).Mmax(int(cat(n_c).haz(i,4))),
     7                 cat(n_c).distance(int(cat(n_c).haz(i,4)),1:4),
     8                 cat(n_c).GM(int(cat(n_c).haz(i,4)),1:2),
     9                 cat(n_c).GM_no(int(cat(n_c).haz(i,4)))
58        format(16(g13.6,2x),3(2x,i1,2x))
      enddo
      close(unit=11)
c
c     block1
c      
      open(unit=11,file=trim(out_stem)//'_'//trim(c5)//c2//'_fract'//
     1                  trim(c7)//'_b1'//trim(c3_)//'.'//trim(c3),
     2                  status='unknown')
      write(11,'(a60)')
     1   'Table of contributions for preparation of deaggregation plot'
      write(11,'(a21,f9.2,a14)')
     1   'of maximum motions in',n_year,' year windows.'
      write(11,'(a42,f9.2,a30)')
     1 'The maximum motions are the peak values of', n_year, 
     2 ' year consecutive time blocks.'
      write(11,53)trim(cat(n_c).cat_name)
      write(11,54)
      do i=1,cat(n_c).nline_hdr-1
          write(11,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(11,*)
      write(11,'(a100)')'    Magnitude       Distance[km]     Contributi
     1on   Avg.times sigma  Stdv.times sigma    # of records '
      deagg_m=deagg_mmin
      do while ((deagg_m - deagg_mmax)<-0.0001)
          deagg_d=deagg_dmin
          do while ((deagg_d - deagg_dmax)<-0.0001)
              deagg_val=0
              sigw_n=0
              do i=n,1,-1
                  n_=int(cat(n_c).haz(i,4))
                  if(cat(n_c).mag(n_)>=deagg_m.and.
     1               cat(n_c).mag(n_)< deagg_m+deagg_mint.and.
     2               cat(n_c).distance(n_,1)>=deagg_d.and.
     3               cat(n_c).distance(n_,1)< deagg_d+deagg_dint)then
                         deagg_val=deagg_val+1
                         sigw_n=sigw_n+1
                         sigw(sigw_n)=cat(n_c).GM(n_,2)
                  endif
              enddo
              call get_avg_std()
              write(11,12)deagg_m+deagg_mint/2., deagg_d+deagg_dint/2.,
     1              deagg_val/real(n), sigw_avg, sigw_std, sigw_n
12            format(5(4x,f9.4,4x),4x,i6)
              deagg_d=deagg_d + deagg_dint
          enddo
          deagg_m=deagg_m + deagg_mint
      enddo
      close(unit=11)

      return
      end
      


c
c     This subroutine finds the maximum of a subcatalog and makes an output
c     file with sorted values used for making GM rate of exceedance tables.
c
      subroutine calc_haz_subcat(n_c,c3)
      use datastruct
      real, dimension (100000)                  :: tmp1, tmp2
      integer                                   :: n_c, n
      character (len=2)                         :: c2
      character (len=3)                         :: c3, c3_
      character (len=5)                         :: c5
      
      cat(n_c).haz(:,1:2)=0.
      
      do i=1,cat(n_c).n_of_subcats
          do j=1,cat(n_c).n_rec
              if(cat(n_c).subcat_no(j) == i .and. 
     1           cat(n_c).GM(j,1) >= cat(n_c).haz(i,1))then
                      cat(n_c).haz(i,1)=cat(n_c).GM(j,1)
                      cat(n_c).haz(i,2)=real(j)
              endif
          enddo
      enddo
      
      n=cat(n_c).n_of_subcats
      do i=1,n
          tmp1(i)=cat(n_c).haz(i,1)
          tmp2(i)=cat(n_c).haz(i,2)
      enddo
      call sort2(n,tmp1(1:n),tmp2(1:n))
      do i=1,n
          cat(n_c).haz(i,1)=tmp1(i)
          cat(n_c).haz(i,2)=tmp2(i)
      enddo

      write(c5,'(a1,i4)')'S',cat(n_c).site_ser
      if(cat(n_c).site_ser<1000)then
          write(c5,'(a1,i3,a1)')'S',cat(n_c).site_ser,' '
      endif
      if(cat(n_c).site_ser<100)then
          write(c5,'(a1,i2,a2)')'S',cat(n_c).site_ser,'  '
      endif
      if(cat(n_c).site_ser<10)then
          write(c5,'(a1,i1,a3)')'S',cat(n_c).site_ser,'   '
      endif
      write(c2,'(a1,i1)')'C',cat(n_c).cat_ser
      write(c3_,'(a1,i2)')'-',n_c
      if(n_c<10)write(c3_,'(a1,i1,a1)')'-',n_c,' '

      open(unit=11,file=trim(out_stem)//'_'//trim(c5)//c2//'_fract'//
     1            trim(c3_)//'.'//trim(c3),status='unknown')

      write(11,51)cat(n_c).n_of_yrs
51    format('This file lists the percentiles of motions with 1/',i6)
      write(11,52)
52    format('occurrences per year.')
      write(11,53)trim(cat(n_c).cat_name)
53    format('This table is extracted from file:',a50)
      write(11,54)
54    format('With the following header:')
      do i=1,cat(n_c).nline_hdr-1
          write(11,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(11,55)
55    format('     %ile     |   GM level   ')
      do i=cat(n_c).n_of_subcats,1,-1
          write(11,56)real((i-1.)/cat(n_c).n_of_subcats)*100,
     1                cat(n_c).haz(i,1)
56        format(1x,g13.5,2x,g13.5)
      enddo
      
      
      return
      end



      subroutine read_input()
      use datastruct
      integer(2)                                :: delfiles 
      call reformat(input,'tempinput')
      open(unit=1,file='tempinput',status='unknown')
      read(1,*)n_of_cats
      do i=1,n_of_cats
          read(1,*)cat(i).cat_name,cat(i).wgt
          read(1,*)cat(i).isort,cat(i).iyears,cat(i).ihaz_of_subcat
!          read(1,*)cat(i).isort,cat(i).isubcat,cat(i).iyears
!          read(1,*)cat(i).yrs_for_percentile
      enddo
      read(1,*)yrs_for_percentile_no
      read(1,*)(yrs_for_percentile(i),i=1,yrs_for_percentile_no)
      read(1,*)out_stem
      read(1,*)deagg_mmin,deagg_mmax,deagg_mint
      read(1,*)deagg_dmin,deagg_dmax,deagg_dint
      close(unit=1)
      do i=1,n_of_cats
          call read_cat_file(i)
      enddo
      delfiles=delfilesqq('tempinput')
      return
      end



      subroutine read_cat_file(n_c)
      use datastruct
      integer                                   :: n_c,n_tmp,m_tmp,o
      pausems=10
      cat(n_c).n_rec=0
      n_tmp=0
      m_tmp=0
      open(unit=2,file=cat(n_c).cat_name,status='unknown')
      do while (not(eof(2)))
          n_tmp=n_tmp+1
          read(2,'(a300)')cat(n_c).cat_head(n_tmp)
          if(n_tmp>=100)exit
      enddo

      cat(n_c).nline_hdr=find_line_num(cat(n_c).cat_head,
     1                                 'Subcatalog No|',n_tmp)
      m_tmp=find_line_num(cat(n_c).cat_head,
     1                    'Site number:',
     2                    cat(n_c).nline_hdr)
      read(cat(n_c).cat_head(m_tmp)(20:),*)cat(n_c).site_ser
      m_tmp=find_line_num(cat(n_c).cat_head,
     1                    'Catalog number:',
     2                    cat(n_c).nline_hdr)
      read(cat(n_c).cat_head(m_tmp)(20:),*)cat(n_c).cat_ser
      m_tmp=find_line_num(cat(n_c).cat_head,
     1                    'Number of subcatalogs:',
     2                    cat(n_c).nline_hdr)
      read(cat(n_c).cat_head(m_tmp)(40:),*)cat(n_c).n_of_subcats
      m_tmp=find_line_num(cat(n_c).cat_head,
     1                    'Period of each subcatalog:',
     2                    cat(n_c).nline_hdr)
      read(cat(n_c).cat_head(m_tmp)(40:),*)cat(n_c).n_of_yrs
      m_tmp=find_line_num(cat(n_c).cat_head,
     1                    'Mixing epistemic/aleatory uncertainties',
     2                    cat(n_c).nline_hdr)
      read(cat(n_c).cat_head(m_tmp)(65:),*)cat(n_c).epist

      rewind(unit=2)
      call skip_lines(2,cat(n_c).nline_hdr)
!      read(2,'(a180)')
      do while (not(eof(2)))
          cat(n_c).n_rec=cat(n_c).n_rec+1
          read(2,*)cat(n_c).subcat_no(cat(n_c).n_rec),
     1             cat(n_c).time(cat(n_c).n_rec,1),
     2             cat(n_c).loc(cat(n_c).n_rec,1:3),
     3             cat(n_c).mag(cat(n_c).n_rec),
     4             cat(n_c).zone_no(cat(n_c).n_rec),
     5             cat(n_c).N0beta(cat(n_c).n_rec,1:2),
     6             cat(n_c).Mmax(cat(n_c).n_rec),
     7             cat(n_c).distance(cat(n_c).n_rec,1:4),
     8             cat(n_c).GM(cat(n_c).n_rec,1:2),
     9             cat(n_c).GM_no(cat(n_c).n_rec)
          cat(n_c).time(cat(n_c).n_rec,2)=
     1    cat(n_c).time(cat(n_c).n_rec,1)+
     2    real(
     3    cat(n_c).subcat_no(cat(n_c).n_rec)-1.)*cat(n_c).n_of_yrs
      enddo
      close(unit=2)
      return
      end
      
      
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   This subroutine finds line number of header (array of nline_chk char   !!!
!!!   acter strings str1) that contains the string str2.                     !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function find_line_num(str1,str2,nline_chk)
      character (len=*)                         :: str1(nline_chk)
      character (len=*)                         :: str2
      integer                                   :: nline_chk, nline_now
      integer                                   :: find_line_num
      find_line_num=-1
      do i=1,nline_chk
          nline_now=i_strcmp(str1(i),str2)
          if(nline_now>0)then
              find_line_num=i
              return
          endif
      enddo
      write(*,*)'didnt find "',str2,'"',find_line_num
      write(*,*)'Program terminated!, check inputs.'
      call exit()
      return
      end



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   This subroutine finds the column number in string str1 that str2 str   !!!
!!!   ing sequence matches. If the length of str2 is larger than str1 or t   !!!
!!!   he squence is not existing a value of -1 is returned.                  !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function i_strcmp(str1,str2)
      implicit none
      integer                                    :: i,n,n1,n2,i_strcmp
      character (len=*)                          :: str1
      character (len=*)                          :: str2
      i_strcmp=-1
      n=len(str1)-len(str2)
      if(n < 0)then
          write(*,*)'String:', trim(str2), 'is longer than line length'
          return
      endif
      do i=1,n+1
          if(str1(i:i+len(str2)-1) == str2)then
              i_strcmp=i
              return
          endif
      enddo
      return
      end



c
c	This subroutine rewrites a new input file for the program to read wit
c	h the comments removed. If the first characters in the beginning of a 
c	line are from the following list the line will be considered as a com
c	ment line. the list is !,@,#,$,%,^,&,*,(,).
c
	subroutine reformat(input,tempinput)
	character (len=1)                         :: c1(10)
	character (len=80)                        :: input, tempinput
	character (len=100)                       :: line
	integer                                   :: flag
	data c1/'!','@','#','$','%','^','&','*','(',')'/
	open(unit=11,file=input,status='old')
	open(unit=12,file=tempinput,status='unknown')
	do while (not(eof(11)))
		flag=1
		read(11,'(a100)')line
		do i=1,10
			if(line(1:1).eq.c1(i))flag=0
		enddo
		if(flag.eq.1)then
			write(12,'(a)')line(1:80)
		endif
	enddo
	close(unit=11)
	close(unit=12)
	return
	end
	
	
	
c
c	This subroutine allows to skip 'lines' records in the unit # 'iunit'.
c
	subroutine skip_lines(iunit,lines)
	integer                                   :: iunit, lines
	if(lines.ge.1)then
		  do i=1,lines
	        read(iunit,*)
		  enddo
	endif
	return
	end
      
      
      
      subroutine find_sorted_indices(n_c)
      use datastruct
      integer                                   :: n_c, nl
      nl=cat(n_c).n_rec
      do i=1,nl
          cat(n_c).indices(i)=real(i)
          cat(n_c).temp(i)=cat(n_c).GM(i,1)
      enddo
      call sort2(nl,cat(n_c).temp(1:nl),cat(n_c).indices(1:nl))
      return
      end



      SUBROUTINE sort2(n,arr,brr)
      INTEGER n,M,NSTACK
      REAL arr(n),brr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,b,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          b=brr(j)
          do 11 i=j-1,l,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
            brr(i+1)=brr(i)
11        continue
          i=l-1
2         arr(i+1)=a
          brr(i+1)=b
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        temp=brr(k)
        brr(k)=brr(l+1)
        brr(l+1)=temp
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
          temp=brr(l)
          brr(l)=brr(ir)
          brr(ir)=temp
        endif
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
          temp=brr(l+1)
          brr(l+1)=brr(ir)
          brr(ir)=temp
        endif
        if(arr(l).gt.arr(l+1))then
          temp=arr(l)
          arr(l)=arr(l+1)
          arr(l+1)=temp
          temp=brr(l)
          brr(l)=brr(l+1)
          brr(l+1)=temp
        endif
        i=l+1
        j=ir
        a=arr(l+1)
        b=brr(l+1)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        temp=brr(i)
        brr(i)=brr(j)
        brr(j)=temp
        goto 3
5       arr(l+1)=arr(j)
        arr(j)=a
        brr(l+1)=brr(j)
        brr(j)=b
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END
	
	
	
      subroutine make_sorted_output_w(n_c, c3)
      use datastruct
      integer                                   :: n_c
      integer                                   :: n
      character (len=2)                         :: c2
      character (len=3)                         :: c3,c3_
      character (len=4)                         :: c4
      character (len=5)                         :: c5
      
      write(c5,'(a1,i4)')'S',cat(n_c).site_ser
      if(cat(n_c).site_ser<1000)then
          write(c5,'(a1,i3,a1)')'S',cat(n_c).site_ser,' '
      endif
      if(cat(n_c).site_ser<100)then
          write(c5,'(a1,i2,a2)')'S',cat(n_c).site_ser,'  '
      endif
      if(cat(n_c).site_ser<10)then
          write(c5,'(a1,i1,a3)')'S',cat(n_c).site_ser,'   '
      endif
      write(c2,'(a1,i1)')'C',cat(n_c).cat_ser
      write(c3_,'(a1,i2)')'-',n_c
      if(n_c<10)write(c3_,'(a1,i1,a1)')'-',n_c,' '

      open(unit=10,file=trim(out_stem)//'_'//trim(c5)//c2//'_sort'//
     1            trim(c3_)//'.'//trim(c3),status='unknown')
      do i=1,cat(n_c).nline_hdr
          write(10,'(a300)')cat(n_c).cat_head(i)
      enddo
      do i=cat(n_c).n_rec,1,-1
          n=int(cat(n_c).indices(i))
          write(10,'(16(g13.6,2x),(7x,i2))')
     1          cat(n_c).subcat_no(n),   cat(n_c).time(n,1),
     2          cat(n_c).loc(n,:),       cat(n_c).mag(n),
     3          cat(n_c).zone_no(n),     cat(n_c).N0beta(n,:),
     4          cat(n_c).Mmax(n),        cat(n_c).distance(n,:),
     5          cat(n_c).GM(n,:),        cat(n_c).GM_no(n)
      enddo
      close(unit=10)
      return
      end
      
      
      
      subroutine make_deagg_output(n_c, n_year, c3)
      use datastruct
      real                                      :: n_year
      integer                                   :: n_c, n_bigs, n
      character (len=2)                         :: c2
      character (len=3)                         :: c3, c3_
      character (len=5)                         :: c5
      character (len=7)                         :: c7
      
      n_bigs=int(real(cat(n_c).n_of_subcats*cat(n_c).n_of_yrs)/
     1           real(n_year))
      write(c5,'(a1,i4)')'S',cat(n_c).site_ser
      if(cat(n_c).site_ser<1000)then
          write(c5,'(a1,i3,a1)')'S',cat(n_c).site_ser,' '
      endif
      if(cat(n_c).site_ser<100)then
          write(c5,'(a1,i2,a2)')'S',cat(n_c).site_ser,'  '
      endif
      if(cat(n_c).site_ser<10)then
          write(c5,'(a1,i1,a3)')'S',cat(n_c).site_ser,'   '
      endif
      write(c2,'(a1,i1)')'C',cat(n_c).cat_ser
      write(c7,'(i7)')int(n_year)
      if(n_year<1000000)write(c7,'(i6,a1)')int(n_year),' '
      if(n_year<100000)write(c7,'(i5,a2)')int(n_year),'  '
      if(n_year<10000)write(c7,'(i4,a3)')int(n_year),'   '
      if(n_year<1000)write(c7,'(i3,a4)')int(n_year),'    '
      if(n_year<100)write(c7,'(i2,a5)')int(n_year),'     '
      if(n_year<10)write(c7,'(i1,a6)')int(n_year),'      '
      write(c3_,'(a1,i2)')'-',n_c
      if(n_c<10)write(c3_,'(a1,i1,a1)')'-',n_c,' '

      open(unit=10,file=trim(out_stem)//'_'//trim(c5)//c2//'_deagg'//
     1                  trim(c7)//'_b2'//trim(c3_)//'.'//trim(c3),
     2                  status='unknown')
      write(10,'(a29,1x,f9.2,1x,a25)')
     1           'List of records of motions in', n_year,
     2           'year return period level.'
      write(10,*)
      do i=1,cat(n_c).nline_hdr-1
          write(10,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(10,*)
      write(10,'(a300)')cat(n_c).cat_head(cat(n_c).nline_hdr)
      do i=cat(n_c).n_rec,cat(n_c).n_rec-n_bigs+1,-1
          n=int(cat(n_c).indices(i))
          write(10,'(16(g13.6,2x),(7x,i2))')
     1          cat(n_c).subcat_no(n),   cat(n_c).time(n,1),
     2          cat(n_c).loc(n,:),       cat(n_c).mag(n),
     3          cat(n_c).zone_no(n),     cat(n_c).N0beta(n,:),
     4          cat(n_c).Mmax(n),        cat(n_c).distance(n,:),
     5          cat(n_c).GM(n,:),        cat(n_c).GM_no(n)
      enddo
      close(unit=10)
      
      open(unit=10,file=trim(out_stem)//'_'//trim(c5)//c2//'_deagg'//
     1                  trim(c7)//'_b1'//trim(c3_)//'.'//trim(c3),
     2                  status='unknown')
      write(10,'(a56,1x,f9.2,1x,a33)')
     1        'Table of contributions for making deaggregation plots of'
     2        , n_year, 'year return period level motions.'
      write(10,*)
      do i=1,cat(n_c).nline_hdr-1
          write(10,'(a300)')cat(n_c).cat_head(i)
      enddo
      write(10,*)
      write(10,'(a100)')'    Magnitude       Distance[km]     Contributi
     1on   Avg.times sigma  Stdv.times sigma    # of records '
      deagg_m=deagg_mmin
      do while ((deagg_m - deagg_mmax)<-0.0001)
          deagg_d=deagg_dmin
          do while ((deagg_d - deagg_dmax)<-0.0001)
              deagg_val=0
              sigw_n=0
              do i=cat(n_c).n_rec,cat(n_c).n_rec-n_bigs+1,-1
                  n=int(cat(n_c).indices(i))
                  if(cat(n_c).mag(n)>=deagg_m.and.
     1               cat(n_c).mag(n)< deagg_m+deagg_mint.and.
     2               cat(n_c).distance(n,1)>=deagg_d.and.
     3               cat(n_c).distance(n,1)< deagg_d+deagg_dint)then
                         deagg_val=deagg_val+1
                         sigw_n=sigw_n+1
                         sigw(sigw_n)=cat(n_c).GM(n,2)
                  endif
              enddo
              call get_avg_std()
              write(10,11)deagg_m+deagg_mint/2., deagg_d+deagg_dint/2.,
     1            deagg_val/real(n_bigs),sigw_avg,sigw_std,sigw_n
11            format(5(4x,f9.4,4x),4x,i6)
              deagg_d=deagg_d + deagg_dint
          enddo
          deagg_m=deagg_m + deagg_mint
      enddo
      close(unit=10)

      return
      end
      
      
      
      subroutine get_avg_std()
      use datastruct
      real                            :: x2b, xb
      x2b=0.
      xb =0.
      if(sigw_n == 0)then
          sigw_avg=0.
          sigw_std=0.
          return
      endif
      do i=1,sigw_n
          xb =xb + sigw(i)/real(sigw_n)
          x2b=x2b+(sigw(i)**2.)/real(sigw_n)
      enddo
      sigw_avg=xb
      sigw_std=sqrt(x2b-xb**2.)
      return
      end      
      
      
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   Added block for calculating hazard of blocks (mainly intended for fractile
!!!   calculations).
!!!   This subroutine finds the first and the last indices opf subcatalog "subca
!!!   t_num in the catalog n_c.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine find_subcat_indices(n_c,subcat_num,n_first,n_last)
      use datastruct
      integer                                   :: n_c, subcat_num
      integer                                   :: n_first, n_last
      do i=n_first,cat(n_c).n_rec
          if(subcat_num /= cat(n_c).subcat_no(i))exit
          n_last=i
      enddo
      return
      end      
      
      
      
      subroutine find_hazvalue(n_c,nyr,subcat_num,
     1                         n_first,n_last,hazvalue)
      use datastruct
      real, allocatable                         :: temp_arr(:)
      integer                                   :: n, n_first, n_last
      integer                                   :: n_c, nyr, subcat_num
      npoint=n_last-n_first+1
      allocate (temp_arr(npoint))
      do i=1,npoint
          temp_arr(i)=cat(n_c).GM(i+n_first-1,1)
      enddo
      n=int(real(cat(n_c).n_of_yrs)/real(yrs_for_percentile(nyr)))
      if(n == 0)n=1
      call sort(npoint,temp_arr) 
      hazvalue=temp_arr(1+npoint-n)
      deallocate(temp_arr)
      return
      end
      
      
      
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   Added block for calculating hazard of blocks (mainly intended for fractile
!!!   calculations).
!!!   This subroutine writes the hazard curves of each subcatalog in separate fi
!!!   le (unit 98) as well as all hazard curves in one file (unit 97- note that
!!!   the block carries recurrence rates of each subcatalog in a row, in other w
!!!   ords the recurrence block is transposed. The first row shows the motion le
!!!   vels).
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine find_hazcurve(n_c,subcat_num,n_first,n_last)
      use datastruct
      real, dimension(61,2)                     :: hazcurve
      real, allocatable                         :: temp_arr(:)
      integer                                   :: n, n_first, n_last
      integer                                   :: n_c, subcat_num
      hazcurve=0.
      npoint=n_last-n_first+1
      allocate (temp_arr(npoint))
      do i=1,npoint
          temp_arr(i)=cat(n_c).GM(i+n_first-1,1)
      enddo
      call sort(npoint,temp_arr) 
      do i=1,61
          hazcurve(i,1)=0.1*10.**(0.1*real(i-1))
      enddo
      do i=1,61
          do j=npoint,1,-1
              hazcurve(i,2)=(npoint-j)/real(cat(n_c).n_of_yrs)
              if(temp_arr(j)<hazcurve(i,1))exit
          enddo
          if(j == npoint)exit
      enddo
      do k=1,i-1
          write(22,'(2(3x,g14.6,4x))')(hazcurve(k,l),l=1,2)
      enddo
      write(22,'(2(3x,g14.6,4x))')temp_arr(npoint),
     1                            1./real(cat(n_c).n_of_yrs)
      
      if(subcat_num==1)then
          write(23,'(61(2x,f9.4))')(hazcurve(k,1),k=1,50)
      endif
      write(23,'(61(2x,f9.6))')(hazcurve(k,2),k=1,i)
      deallocate(temp_arr)
      return
      end
      
      
      
      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
          do j=l+1,ir
              a=arr(j)
              do i=j-1,l,-1
                  if(arr(i).le.a)goto 2
                  arr(i+1)=arr(i)
              enddo
              i=l-1
2             arr(i+1)=a
          enddo
          if(jstack.eq.0)return
          ir=istack(jstack)
          l=istack(jstack-1)
          jstack=jstack-2
      else
          k=(l+ir)/2
          temp=arr(k)
          arr(k)=arr(l+1)
          arr(l+1)=temp
          if(arr(l).gt.arr(ir))then
              temp=arr(l)
              arr(l)=arr(ir)
              arr(ir)=temp
          endif
          if(arr(l+1).gt.arr(ir))then
              temp=arr(l+1)
              arr(l+1)=arr(ir)
              arr(ir)=temp
          endif
          if(arr(l).gt.arr(l+1))then
              temp=arr(l)
              arr(l)=arr(l+1)
              arr(l+1)=temp
          endif
          i=l+1
          j=ir
          a=arr(l+1)
3         continue
              i=i+1
          if(arr(i).lt.a)goto 3
4         continue
              j=j-1
          if(arr(j).gt.a)goto 4
          if(j.lt.i)goto 5
          temp=arr(i)
          arr(i)=arr(j)
          arr(j)=temp
          goto 3
5         arr(l+1)=arr(j)
          arr(j)=a
          jstack=jstack+2
          if(jstack.gt.NSTACK)pause
          if(ir-i+1.ge.j-l)then
              istack(jstack)=ir
              istack(jstack-1)=i
              ir=j-1
          else
              istack(jstack)=j-1
              istack(jstack-1)=l
              l=i
          endif
      endif
      goto 1
      END