cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     EqHaz01 Ver. A (2010/05/27)
c
c     The programs Haz01 can be used for generating synthetic catalog of earth
c     quakes. The synthetic catalog is generated based on random draw of earth
c     quake parameters (time, location, and magnitude) given the seismicity pa
c     rameters of seismic sources and their geometries. The intention is to us
c     e it as a part of set of programs to be used in probabilistic hazard ass
c     essment using Monte Carlo approach of Musson (1999 and 2000).
c
c     	* The program is developed using Compaq Visual Fortran 6.6 and Intel® 
c         Visual Fortran Compiler Professional Edition 11.1. Compatibility of 
c     	  the code with other compilers is not tested.
c     	* The program is initially developed for Canadian applications so th
c     	  e types of inputs are consistent with the inputs of NBCC 2005 haza
c     	  rd maps.
c     	* The program is capable of generating any number of catalogues for 
c     	  required time periods.
c     	* The version 1.0 of program is capable of generating synthetic cata
c     	  loges for areal sources. The case of linear sources is left for fu
c     	  ture developments.
c     
c     The program generates synthetic cataloges by random draw of earthquake p
c     arameters, including, time, location, depth, and magnitude as explained 
c     below.
c     1-Time: Generation of time is easily done by multiplication of uniformly
c       distributde random number (between zero and one) by the duration of ca
c       taloges.
c     2-Coordinates: Generation of coordinates is done in two steps:
c     	A-A grid of points with equal distance along NS and EW directions is
c     	  generated and stored in the list of zones' parameters.
c     	B-For each earthquake randomly one of the grid points is selected, t
c     	  he grid spacing is multiplied by a uniformly distributed random nu
c     	  mber between -0.5 and 0.5, and this number is added to the longitu
c     	  de of the grid point. repeatition of the same procedure for the gr
c     	  id point latitude provides the random coordinates for location.
c     	This two stage procedure is implemented to save the time to check if 
c     	a point lies in a polygon in each coordinate generation. Although th
c     	is time saving is compensated by allocating a part of the memory for 
c     	the storage of the coordinates of the grid points.
c     3-Depth: Generation of depth is simple, however, it could be done in two
c       ways depending on  the user's choice. The choice of the method can be 
c       different for different zones and the program is developed to handle i
c       t.
c     	* Choosing among (maximum of five) given depths with different proba
c     	  bilities of occurences. This option is applied if 'ndepth' value o 
c     	  f a source is an integer between one and five.
c     	* Choosing a depth from a uniformly distributed depth range. This op
c     	  tion is applied if 'ndepth' value of a source is set to zero. In t
c     	  his case two depths should be entered in input file.
c     4-Magnitude: The generation of magnitude is the most complicated one sin
c       ce five parameters with their probabilities are entering in the proces
c       s of random draw. The following parameters are entered for each source:
c     	* Beta, N0 and their probabilities (one probability for each pair an
c     	  d a maximum of five pairs).
c     	* Mmax and its probabilities (a maximum of five values are permitted
c     	  as given by nmmax. If nmmax is set to zero, an Mmax value from uni
c     	  form distribution will be chosen).
c     	* Mmin is the minimum accepted magnitude and is identical for all so
c     	  urces.
c     	* Recurrence relation choice among GR, truncated GR, and tapered GR.
c     	  The choice of relation is controlled by parameter 'irecurr' for ea
c     	  ch source.
c       Magnitude selection does ndraw=100000 iterations per source per catalo
c         g. For each iteration the program performs the following steps:
c     	A-Chooses a Beta and N0 pair among (maximum of five) possibilities c
c     	  onsidering the probabilities.
c     	B-Chooses an Mmax value among (maximum of five) possibilities consid
c     	  ering the probabilities or from continuous distribution.
c     	C-Based on the drawn values of Beta, N0, and Mmax, the program calcu
c     	  lates an m0 value, using truncated GR distribution, so that the nu
c     	  mber of earthquakes with M>=m0 for the period of the catalog will 
c     	  be ndraw=100000.
c     	D-The program passes m0, Mmax, and b (Beta) to selected subroutine (
c     	  GR, truncated GR, or tapered GR) for a random realization of 'M'.
c     	E-Compares the value of 'M' with minimum magnitude.
c
c Modifications:
c 2010/09/10:
c     -Changed the input format as suggested by G. Atkinson.
c     -Changed the format of sample input file to match with program.
c     -Developed the program "reformat" that reads old input files and reforma
c     ts them to new format.
c     -Added subroutine "checkinput" that checks sum of the weigts of seismici
c      city, mmax, and depth parameters. If weights of any parameters of a sou
c      rce don't sum up to one(1), that source will be disabled from generatin
c      g catalogues.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     EqHaz01 Ver. B (2011/11/20)
c
c     1-Program is cleaned up for reporting only useful outputs.
c     2-The program is modified to be capable of generating catalogs for multi
c      -segment fault zones. The designed algorithm allows the earthquakes to 
c      be generated on any ppoint on fault surface with equal probability.
c     3-Bug fixes are made to the sorting segments of the code and the modifie
c      d "sort" subroutine is changed with "sort2" from numerical recipes.
c     4-Added feature to mark the time stamp in output files. This is intended
c      for keeping track of catalogs processed by all EqHaz modules.
c     5-Adding a feature to report an echo file of input parameters.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     EqHaz01 Ver. C (2012/02/15)
c
c     1-The code is modified so that it accepts minimum magnitude for each sou
c      rce.
c     2-The program checks the minimum magnitude of each source and subcatalog
c      durations and signals the possibility of incomplete gatalog generation
c      and sugests proper parameter adjustments.
c     3-The generator of catalogs for fault sources is modified so that the ea
c      rthquake coordinates are driven by the choices of depth ranges.
c     4-The geometries of individual fault sources are fixed to have constant 
c      width, dip angle, and depth of upper corner for all segments.
c     5-The program is tested with some erronous inputs for checking the prope
c       r warning signals of program.
c     6-A new Capability is added to let the user to specify cumulative distri
c       bution function (CDF)of magnitudes in a file for each source. The only
c       restriction in this case is that the user can specify only one N0-beta
c       numbers and the beta is read but not used by program. (2012/03/08)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     EqHaz1_2012 Ver. c (2012/06/27)
c
c     1-A new capability is added to let the user to specify probability mass 
c       function (PMF) of magnitudes in a file for each source. The program wi
c       ll read N0-beta values, but will discard the beta and will generate N=
c       N0*(subcatalog duration) earthquakes with given PMF distribution.
c     2-The new version of program can be set to mix (or not mix) epistemic an
c       d aleatory uncertainties. When mixing is chosen every record of a zone
c       and subcatalog will adopt any possible combination of Mmax, depth, and
c       N0-beta combination. But when the program is in "not mix" mode one com
c       bination of these parameters will be drown and applied to all earthqua
c       kes of a zone in a subcatalog. When unmixing is chosen the program rep
c       orts another file named "*applied.par" with applied parameters for rec
c       ords of each source in each subcatalog.
c     3-This version reports the name of input file name and uncertainty flag
c       values in the output with purpose of allowing EqHaz2_2012c read (and u
c       se) these information.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     EqHaz1_2012 Ver. d (2012/08/03)
c
c     A new capability is added to this program to allow for generation of syn
c     thetic catalogs based on smooth seismicity approach. The smooth seismici
c     ty approach is applied through setting a flag in the input equal to one.
c     When smooth seismicity approach is adopted the program looks for two fil
c     es in the working directory: 1-A catalog file in CCSC format; and 2-A pa
c     rameter file that sets the values of the following parameters to each so
c     urce zone (note that name of these two files are given in EqHaz1 input):
c         *-Set of completeness magnitude and corresponding dates
c         *-Correlation distance
c         *-Width of areal rings in which events are counted
c     Temporarily a file is generated by the program with grid point coordinat
c     es inside each source zone and corresponding weight in that zone (this f
c     eature will be removed from the code).
c     As a result of huge amount of calculations, the program pauses for some
c     time (may be about a minute) and then quickly starts generating syntheti
c     c catalogs of earthquakes.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     This module does the definitions of seismic sources and input to program.
c     For each source (both areal and linear) a structure of components is gen
c     erated (the structure "source") and array of these structures are genera
c     ted ("zones" and "fault"). The definition of this module allows us to de
c     fine and declare these values once and implement them in any subroutine
c     when needed.
c     The module lists all the inputs that the program needs for running.
c
      module srcdef																
        structure /source/                                                      !Starting the structure of sources.
            character*8	 name                                                   !Source name. Maximum of eigth caracters are reserved.
            integer      typ                                                    !Source type flag (1:zone, 2:fault).
            integer      iuse                                                   !Source usage flag (1:use this source, 2:ignore this source).
            integer      iflg(10)                                               !Flags reserved for assigning properties for each source
            integer      nnode                                                  !Number of vertices of polygons of areal sources or node points linear sources.
            integer      ngrid                                                  !Number of grid points in polygons of areal sources or along linear sources.
            real         grd_space                                              !Grid points in polygons of areal sources or along linear sources.
            integer      nb_n0_val                                              !Number of Beta, N0, and their weight sets. A total of 5 sets are allowed.
            real         b_n0_val(5,3)                                          !Array of Beta, N0, and their weight sets.
            integer      nmmax                                                  !Number of Mmax values and their weights. A total of 5 sets are allowed. Zero means uniform distribution between two values.
            real         mmax(5,2)                                              !Mmax values and their weights. If nmmax=0 then two mmax values should be entered for uniform distribution range limits.
            real         mmin                                                   !Minimum magnitude considered in the catalog of a source, values smaller than this value will not appear in catalog.
            integer      irecurr                                                !Recurrence relation flag (1:GR, 2:Truncated GR, 3:Tapered GR).
            integer      ndepth                                                 !Number of Depth values and their weights. A total of 5 sets are allowed. Zero means uniform distribution between two depths.
            real         depth(5,2)                                             !Depth values and their weights. If ndepth=0 then two depth values should be entered for uniform distribution range limits.
            real         nodes(50,2),grids(90000,2)                             !Array of nodes (reserved 50 coordinates) and array of grid points (reserved 90000 coordinates). The geographical coordinates in degrees should be used (Order: Lat. long.).
            real         grids_w(90000)                                         !Array of weights for grid points of each source zone. (Added for incorporation of smooth seismicity approach 2012/07/29)
            real         geom(50,3)                                             !Array of geometrical information used for linear sources (width, dip, and depth of fault segments-Added on 2010/09/10)
            real         perturb(2)                                             !Array of horizontal and vertical perturbation values in kilometers (Added on 2010/9/28)
            real         extra(50,3)                                            !An array of extra information, for faults I use for storing lengths, projection width on on surface, and azimuths of fault segments (Added on 2010/09/30)
            character*30 file_recurr                                            !Name of a file with probability distribution (CDF) of events in a source. This works if irecurr=4 (Added on 2012/03/07)
            real         tab_recurr(10000,2)                                    !Table probability distribution (CDF) of events in a source. This works if irecurr=4 (Added on 2012/03/07)
            integer      n_tab_recurr                                           !Number of magnitude-probability pairs of probability distribution (CDF) of events in a source. This works if irecurr=4 (Added on 2012/03/07)
        end structure
        character*500    headline                                               !Header line in the input file.  Added on 2011/10/26
        character*12     real_clock(3)
        integer          numcat                                                 !Number of catalogues to be generated.
        real             numyears                                               !Number of years of each catalog.
        integer          nzones                                                 !Number of areal sources. A total of 50 is reserved. If more is requird change the number in following "record" declaration statement.
        integer          nfaults                                                !Number of linear sources. A total of 10 is reserved. If more is requird change the number in following "record" declaration statement.
        integer          nsource                                                !Total number of sources, including areas and lines (Added on 2010/09/28)
        integer          seed, epist, ismooth                                   !Seed value for randon number generator (0: random, none-zero values will enter to the program). Mixing epistemic/aleatory uncertainty flag (use 1 for mixing and 0 for not mixing [traditional]). Smooth seismicity flag (Yes: 1, No:0). (Last change 2012/07/06)
        
        
        record	/source/source(60)												
        record  /source/tmp_src
      end module srcdef


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     This module makes room for earthquake catalog storage and use in case sm
c     ooth seismicity approach is selected. (Last modified on 2012/07/06)
c
      module smooth_seis																
          structure /eq_cat/                                                    !Starting the structure of earthquake catalog.
              character (len=50)                  ::catfile, parfile            !Earthquake catalog input file name. Earthquake smoothing parameter file name
              character (len=512),dimension(200)  ::head                        !Earthquake catalog header lines
              integer                             ::nline_head, n_rec           !Number of earthquake catalog header lines, number of earthquake catalog record lines
              integer, dimension(1000000,5)       ::time                        !Date and time of each records of earthquake (year, month, day, hour, and minute)
              real, dimension(1000000,3)          ::coord                       !Earthquake coordinate
              real, dimension(1000000,9)          ::mag                         !Magnitudes in input catalog (MB, MN, ML, MS, MC, MD, MW, MM, and MZ)
              real, dimension(1000000)            ::Mpf                         !Preferred magnitude
              character (len=4),dimension(1000000)::Tm                          !Preferred magnitude type
              character (len=1),dimension(1000000)::dd                          !Depth designation
              integer, dimension(1000000,2)       ::czf                         !Catalog and zone flags
              real, dimension(1000000)            ::mf                          !Mw additive conversion factor [magnitude units]
              real, dimension(1000000)            ::Mw                          !calculated moment magnitude value
          end structure

          structure /smooth_src/
              real                                ::ring_w, corr_dis            !Ring width, correlation distance for smoothing
              integer                             ::n_comp                      !Number of pairs of year and magnitude completenesses
              real, dimension(10)                 ::comp_yr, comp_m             !Arrays of years and starting date of completeness
              real                                ::end_year                    !End year in real catalog considered for this source
              real                                ::beta                        !Beta value to be used for in Equation 10 in Weichert 1980 for rate calculation
              real                                ::bkground                    !This parameter is the fraction of seismicity of a source that goes to background
          end structure
          
          record /eq_cat/cat_smo
          record /smooth_src/source_smo(60)
          
      end module smooth_seis
      
      

      program               EqHaz1_2012d
      use                   srcdef                                              !Implements the module srcdef that includes definitions of sources and their parameters.
      use                   smooth_seis                                         !This module is used for reserving area and access to earthquake catalog values for smooth seismicity  approach (Last modified on 2012/07/06)
      real                  magcat(100000,10,60)                                !'magcat' stores the records of up to 100000 events for (up to) 60 sources per catalog. This array is refilled in the process of each catalogue generation and in the end its contents are stored in output file.
      real,allocatable ::   numyears_cat(:,:)                                   !The array 'numyears_cat' stores the records of all sources in 'magcat'. This array is implemented for sorting the records of a catalog chronologically.
      real,allocatable ::   srt(:),indx(:)                                      !These two arrays are added to help implementing sort2 subroutine (from Numerical recipes) in this program
      integer               num(60),cat_event_num,n_count                       !'num stores the number of records of a source in each catalog, cat_event_num and n_count store the total number of records in each catalog.
      character*50          input, output                                       !Variables for storing input and output file names.
c
c     This section asks for input and output file names, opens the output 
c     file, and writes a header line with names of columns.
c
      write(*,*) 'Enter input file name: '
      read(*,*)	 input
      write(*,*) 'Enter output file stem e.g. hazrun, where the program'
      write(*,*) 'generate hazrun.out and hazrun.echo for catalog and '
      write(*,*) 'and implemented parameters files: '
      read(*,*)	 output
      open(unit=2,file=trim(output)//'.out',status='unknown')
c
c     This section performs some initialization jobs.
c     
      call date_and_time(date=real_clock(1), time=real_clock(2))
c     call date_and_time(real_clock(1), real_clock(2), 
c     1                   real_clock(3), date_time)
      call readinput(input)                                                     !Read the input file and store the values in the parameter structure 'srcdef'.
      call writepar(input,output)                                               !Makes a copy of input file called "stem.echo with time stamp on top for future access by user. Added on 2011/10/26
      call checkinput()                                                         !Added on 2010/09/10. Checks the sum of seismicity, mmax, and depth parameters' weights.
      call fillgrid()                                                           !Fill the grid points in source zones and store coordinates in corresponding arrays.
      if(ismooth/=0)then                                                        !Added on 2012/07/06 for reading earthquake catalog if smooth seismicity approach is selected
          write(*,*)'Calculating grid weights for smooth seismicity'
          call read_smooth_catfile()                                            !Added on 2012/07/06 for reading earthquake catalog if smooth seismicity approach is selected
          call read_smooth_parfile()                                            !Added on 2012/07/06 for reading earthquake catalog if smooth seismicity approach is selected
          call fillgrid_smo()                                                   !This subroutine calculates the weights of each grid point output from fillgrid based on smooth seismicity approach.
      endif                                                                     !Added on 2012/07/06 for reading earthquake catalog if smooth seismicity approach is selected
      do lm=1,nsource
          if(source(lm).typ.eq.2)cycle
          write(51,*)'Source # ',lm
          do ll=1,source(lm).ngrid
              write(51,*)source(lm).grids(ll,:),source(lm).grids_w(ll)
          enddo
      enddo
!      goto 100
      
      call randnumgen(rnd,seed)                                                 !Create an array of random numbers with the provided 'seed' values. 'randnumgen' is developed and applied this way for practical purposes.

      write(2,'(a34, a8)')'EqHaz1 run date (yyyymmdd)  :     ',                 !Adds date stamp of program run on top of catalog file. Added on 2011/10/26
     1                     trim(real_clock(1))
      write(2,'(a34,a10)')'EqHaz1 run time (hhmmss.sss):     ',                 !Adds time stamp of program run on top of catalog file. Added on 2011/10/26
     1                     trim(real_clock(2))
      write(2,'(a34,a50)')'EqHaz1 input file name:           ',                 !Adds name of input file name of current EqHaz1 run 2012/06/18
     1                     input
      write(2,'(a500)')headline                                                 !Adds header line of input file on top of catalog file. Added on 2011/10/26
      write(2,'(a28,i8)')'Number of subcatalogs:        ',numcat                !Adds number of catalogs on top of catalog file. Added on 2011/10/26
      write(2,'(a28,i8,2x,a5)')'Period of each subcatalog:    '                 !Adds time period of each catalog on top of catalog file. Added on 2011/10/26
     1                         ,int(numyears),'years'
      write(2,'(a55,i2)')'Mixing epistemic/aleatory uncertainties '//           !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                   '(Yes=1,No=0):  ',epist        
      write(2,'(a150)')' Subcatalog No| Time (years) |   Latitude   |   
     1Longitude  |   Depth(km)  |   Magnitude  |  Zone Number |  N0 per 
     2year |     Beta     |     Mmax     |'


      
      if(epist.eq.0)then                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          open(unit=3,file=trim(output)//'_applied.par',status='unknown'        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1        )
          write(3,'(a34,a12)')'EqHaz1 run date (yyyymmdd)  :     ',             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                         real_clock(1)                                    !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          write(3,'(a34,a12)')'EqHaz1 run time (hhmmss.sss):     ',             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                         real_clock(2)                                    !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          write(3,'(a34,a50)')'EqHaz1 input file name:           ',             !Adds name of input file name of current EqHaz1 run 2012/06/18
     1                     input
          write(3,'(a500)')headline                                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          write(3,'(a28,i8)')'Number of subcatalogs:        ',numcat            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          write(3,'(a28,i8,2x,a5)')'Period of each subcatalog:    '             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                             ,int(numyears),'years'                       !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          write(3,'(a90)')' Subcatalog No|  Zone Number |   Depth(km)  '        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                //'|  N0 per year |     Beta     |     Mmax     |'        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      endif                                                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
c
c     The main catalog generating program iterates on the cataloges and mak
c     es one catalog per cycle and writes its results on output file. There
c     are 'numcat' number of iterations.
c
      do i=1,numcat
          write(*,'(a10,i5)')'Subcat.#',i
          magcat=0.
          cat_event_num=0
          do j=1,nsource
c
c      This loop iterates on the sources and does the random generation
c      of earthquakes' magnitudes, locations, depths, and times and coun
c      ts the number of all events in the current catalog. The records o
c      f each zone are stored in different layers of array 'magcat'.
c
c       write(*,'(a16,i5,a11,i5)')'Processing Cat.#',i,' and Src. #',j
              if(source(j).iuse.eq.1)then
                  call magmaker(j,num(j),magcat(:,:,j))
                  call time_depth_coord_gen(j,num(j),magcat(:,:,j))
                  cat_event_num=cat_event_num+num(j)
                  if(epist.eq.0)then                                            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                      write(3,'(4x,i7,8x,i7,8x,f7.2,6x,f11.5,                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1                          6x,f7.2,8x,f7.2)')i,j,                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     2                          magcat(num(j),4,j),                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     3                          magcat(num(j),7,j),                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     4                          magcat(num(j),8,j),                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     5                          magcat(num(j),9,j)                              !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  endif                                                         !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              endif
          enddo
          allocate (numyears_cat(cat_event_num,10))
          allocate (srt(cat_event_num))
          allocate (indx(cat_event_num))
          n_count=0
          do j=1,nsource
c
c     This loop stores the records of all source zones in one table and
c     counts the records. The values of the new table are stored in the 
c     array 'numyears_cat'. This puts the records of all zones in one p
c     lace and allows for sorting.
c
              if(source(j).iuse.eq.1)then
                  numyears_cat(n_count+1:n_count+num(j),:)=
     1                                     magcat(1:num(j),:,j)
                  n_count=n_count+num(j)
              endif
          enddo
c
c     This stage sorts the records of all zones in a catalog chronologi
c     cally and then stores the values of current catalog in the output
c     file.
c
c     call sort(cat_event_num,10,1,numyears_cat,0)
          do ii=1,cat_event_num
              srt(ii)=numyears_cat(ii,1)
              indx(ii)=real(ii)
          enddo
          call sort2(cat_event_num,srt,indx)
          do jj=1,cat_event_num
              j=int(indx(jj))
c              write(2,'(11(g13.6,2x))')i,(numyears_cat(j,k),k=1,9)
              write(2,'(5x,i5,5x,f13.3,6x,f7.3,7x,f8.3,8x,f7.3,9x,f6.3
     1               ,10x,i3,9x,f9.1,7x,f7.3,9x,f5.2)')i
     2               ,(numyears_cat(j,k),k=1,5),int(numyears_cat(j,6))
     3               ,(numyears_cat(j,k),k=7,9)
              n=int(numyears_cat(j,6))
ccc              write(3)i,(numyears_cat(j,k),k=1,10),source(n)
          enddo
          deallocate (numyears_cat)
          deallocate (srt)
          deallocate (indx)
      enddo

      do i=1,nsource
          if(source(i).iflg(1).eq.1)then
              write(*,*)'*****WARNING*****'
              write(*,*)'Synthetic catalog for Source #',i,
     1                  'is incomplete. Increase mmin or',
     2                  'decrease subcatalog duration.'
          endif
      enddo

      close(unit=2)
      if(epist.eq.0)then
          close(unit=3)
      endif
100   end




c
c     This subroutine generates 'nmag' number of time (in years), depths, a
c     nd coordinates for current zone ('isource'). This subroutine stores t
c     he generated values in proper columns of array 'maglis'. This subrout
c     ine is called once per zone per catalogue.
c
      subroutine time_depth_coord_gen(isource,nmag,maglis)	
      use            srcdef
      real           maglis(100000,10)
c
c     The iteration is over 'nmag' number of events in source # 'isource'. 
c
c
c     	Depth generation segment for zones
c
      if(epist.eq.0)then                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          call randnumgen(rnd,0)                                                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          if(source(isource).ndepth.eq.0)then                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmin=source(isource).depth(1,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmax=source(isource).depth(2,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=rnd*(dmax-dmin)+dmin                                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          else                                                                  !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              do j=1,source(isource).ndepth                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  if(sum(source(isource).depth(1:j,2)).ge.rnd)exit              !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              enddo                                                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call randnumgen(rnd,0)                                            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=source(isource).depth(j,1)+                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1          2*(0.5-rnd)*source(isource).perturb(2)                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          endif                                                                 !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          maglis(1:nmag,4)=dep                                                  !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      endif                                                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      
      do i=1,nmag
c
c     Time generation segment
c
          call randnumgen(rnd,0)
          maglis(i,1)=rnd*numyears
c
c     Coordinate and depth generation segment
c
          if(source(isource).typ.eq.1)then                                      !Condition applied for incorporating linear sources (Added on 2010/09/30)
c         Coordinate generation segment for areal sources
              call find_zone_hypocenter(isource,
     1                     maglis(i,2),maglis(i,3),maglis(i,4))
          elseif(source(isource).typ.eq.2)then                                  !Condition applied for incorporating linear sources (Added on 2010/09/30)
c         Coordinate generation segment for linear sources
              call find_flt_on_surf(isource)
              call find_flt_hypocenter(isource,
     1                     maglis(i,2),maglis(i,3),maglis(i,4))
          endif                                                                 !Condition applied for incorporating linear sources (Added on 2010/09/30)
      enddo
      return
      end



c
c     This subroutine gets the serial number of the areal source as introdu
c     ced in input file and returns the latitude, longitude and depth of as
c     signed hypocenter. The assignment of epicenter is based on random sel
c     ection from the list of generated points for the zone. For longitude 
c     and latitude two values of perturbations. One is equal to [-1,+1] use
c     r defined perturbation length introduced by the user for each compone
c     nt in the input. The other is [-1,+1] times of initial grid interval 
c     in x and y directions. 
c
      subroutine find_zone_hypocenter(isource,xlat,xlon,dep)
      use					srcdef
      real  tmp_sum
      double precision	lat1,lon1,lat2,lon2,dist,azi
c
c     	Coordinate generation segment for zones
c
      call randnumgen(rnd,0)
!      igrid=int(source(isource).ngrid*rnd)+1                                    !This portion is changed to implement smooth seismicity approach for epicenter generation
      if(ismooth.eq.0)then                                                      !This portion is added to implement smooth seismicity approach for epicenter generation
          igrid=int(source(isource).ngrid*rnd)+1                                !This portion is added to implement smooth seismicity approach for epicenter generation
      else                                                                      !This portion is added to implement smooth seismicity approach for epicenter generation
          tmp_sum=0.                                                            !This portion is added to implement smooth seismicity approach for epicenter generation
          do i=1,source(isource).ngrid                                          !This portion is added to implement smooth seismicity approach for epicenter generation
              tmp_sum=tmp_sum+source(isource).grids_w(i)                        !This portion is added to implement smooth seismicity approach for epicenter generation
              if(tmp_sum.ge.rnd)exit                                            !This portion is added to implement smooth seismicity approach for epicenter generation
          enddo                                                                 !This portion is added to implement smooth seismicity approach for epicenter generation
          igrid=i                                                               !This portion is added to implement smooth seismicity approach for epicenter generation
      endif                                                                     !This portion is added to implement smooth seismicity approach for epicenter generation
      ygrid=source(isource).grids(igrid,1)
      xgrid=source(isource).grids(igrid,2)
      call radius_der(ygrid,r,r_der)
      call randnumgen(rnd,0)
      xperturb=rnd*source(isource).perturb(1)
      dx=rad2deg((source(isource).grd_space+xperturb)/
     1	(r*cos(deg2rad(ygrid))))
      call randnumgen(rnd,0)
      yperturb=rnd*source(isource).perturb(1)
      dy=rad2deg((source(isource).grd_space+yperturb)/
     1	(sqrt(r**2+r_der**2)))
      call randnumgen(rnd,0)
      xlat=ygrid+(0.5-rnd)*dy
      call randnumgen(rnd,0)
      xlon=xgrid+(0.5-rnd)*dx
c
c     	Depth generation segment for zones
c
      if(epist.eq.1)then                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          call randnumgen(rnd,0)                                                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          if(source(isource).ndepth.eq.0)then                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmin=source(isource).depth(1,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmax=source(isource).depth(2,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=rnd*(dmax-dmin)+dmin                                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          else                                                                  !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              do j=1,source(isource).ndepth                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  if(sum(source(isource).depth(1:j,2)).ge.rnd)exit              !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              enddo                                                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call randnumgen(rnd,0)                                            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=source(isource).depth(j,1)+                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1            2*(0.5-rnd)*source(isource).perturb(2)                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          endif                                                                 !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      endif                                                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012

      return
      end


c
c     This subroutine gets the serial number of the linear source as introd
c     uced in input file and returns the latitude, longitude and depth of a
c     ssigned hypocenter. The assignment of coordinates are based on random 
c     assignment of a point along the fault length and dip. Then each of th
c     e coordinates are perturbed based on the values given in the input fi
c     le. 
c
      subroutine find_flt_hypocenter(isource,xlat,xlon,dep)
      use					srcdef
      double precision	lat1,lon1,lat2,lon2,dist,azi
      call randnumgen(rnd,0)
      randlen=rnd*source(isource).extra(source(isource).nnode,1)
      complen=0.
      do n=1,source(isource).nnode-1
          complen=complen+source(isource).extra(n,1)
          if(complen.ge.randlen)exit
      enddo
      lat1=dble(source(isource).nodes(n,1))
      lon1=dble(source(isource).nodes(n,2))
      dist=dble(randlen-complen+source(isource).extra(n,1))
      azi =dble(source(isource).extra(n,3))
      call disazi2xy(lat1,lon1,dist,azi,lat2,lon2)

      if(epist.eq.1)then                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          call randnumgen(rnd,0)                                                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          if(source(isource).ndepth.eq.0)then                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmin=source(isource).depth(1,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dmax=source(isource).depth(2,1)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=rnd*(dmax-dmin)+dmin                                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          else                                                                  !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              do j=1,source(isource).ndepth                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  if(sum(source(isource).depth(1:j,2)).ge.rnd)exit              !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              enddo                                                             !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call randnumgen(rnd,0)                                            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              dep=source(isource).depth(j,1)+                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
     1          2*(0.5-rnd)*source(isource).perturb(2)                          !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          endif                                                                 !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      endif                                                                     !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      
      call correct_depth(isource,n,dep)
      if(source(isource).geom(n,2).eq.90.)then
          dist=0.
      elseif(source(isource).geom(n,2).eq.0.)then
          call randnumgen(rnd,0)
          dist=rnd*source(isource).geom(n,1)
      else
          dist=(dep-source(isource).geom(n,3))/
     1         tan(deg2rad(source(isource).geom(n,2)))
      endif

      azi=azi+90.
      call disazi2xy(lat2,lon2,dist,azi,lat1,lon1)
      xlat=real(lat1)
      xlon=real(lon1)
      return
      end	


      subroutine correct_depth(isource,n,dep)
      use					srcdef
      real                                      :: dep, depmin, depmax
      integer                                   :: n
      depmin=source(isource).geom(n,3)
      depmax=source(isource).geom(n,3)+source(isource).geom(n,1)*
     1                     sin(deg2rad(source(isource).geom(n,2)))
      if(dep.lt.depmin)dep=depmin
      if(dep.gt.depmax)dep=depmax
      return
      end



!c    This subroutine is pre Jan.25, 2012.
!c     This subroutine gets the serial number of the linear source as introd
!c     uced in input file and returns the latitude, longitude and depth of a
!c     ssigned hypocenter. The assignment of coordinates are based on random 
!c     assignment of a point along the fault length and dip. Then each of th
!c     e coordinates are perturbed based on the values given in the input fi
!c     le. 
!c
!	subroutine find_flt_hypocenter(isource,xlat,xlon,dep)
!	use					srcdef
!	double precision	lat1,lon1,lat2,lon2,dist,azi
!	call randnumgen(rnd,0)
!	randlen=rnd*source(isource).extra(source(isource).nnode,1)
!	complen=0.
!	do n=1,source(isource).nnode-1
!		complen=complen+source(isource).extra(n,1)
!		if(complen.ge.randlen)exit
!	enddo
!	lat1=dble(source(isource).nodes(n,1))
!	lon1=dble(source(isource).nodes(n,2))
!	dist=dble(randlen-complen+source(isource).extra(n,1))
!	azi =dble(source(isource).extra(n,3))
!	call disazi2xy(lat1,lon1,dist,azi,lat2,lon2)
!	call randnumgen(rnd,0)
!	dist=dble(source(isource).extra(n,2)*rnd)
!	azi=azi+90.
!	call disazi2xy(lat2,lon2,dist,azi,lat1,lon1)
!	if(source(isource).geom(n,2).ne.90.)then
!		dep=source(isource).geom(n,3)+real(dist)*
!	1		tan(deg2rad(source(isource).geom(n,2)))
!	else
!		call randnumgen(rnd,0)
!		dep=source(isource).geom(n,3)+rnd*(source(isource).geom(n,1))
!	endif
!	call randnumgen(rnd,0)
!	dep=dep+2.*(rnd-0.5)*source(isource).perturb(2)
!	call randnumgen(rnd,0)
!	dist=dble(source(isource).perturb(1)*rnd)
!	call randnumgen(rnd,0)
!	azi =dble(360.*rnd)
!	call disazi2xy(lat1,lon1,dist,azi,lat2,lon2)
!	xlat=real(lat2)
!	xlon=real(lon2)
!	return
!	end	
!
!
!
c
c     This subroutine finds the length of whole fault and fault segments, w
c     idth of segments' projection on surface, and azimuths of segments. Th
c     en it saves these values in an array called extra.(Added on 2010/10/01)
c
      subroutine find_flt_on_surf(isource)
      use					srcdef
      double precision	lat1,lon1,lat2,lon2,dist,azi
      source(isource).extra(source(isource).nnode,1)=0.0
      do i=1,source(isource).nnode-1
        lat1=dble(source(isource).nodes(i,1))			
        lon1=dble(source(isource).nodes(i,2))			
        lat2=dble(source(isource).nodes(i+1,1))			
        lon2=dble(source(isource).nodes(i+1,2))			
        call xy2disazi(lat1,lon1,lat2,lon2,dist,azi)
        source(isource).extra(i,1)=real(dist)                                   !Save length of fault segment
        source(isource).extra(i,2)=source(isource).geom(i,1)*
     1      cos(deg2rad(source(isource).geom(i,2)))                             !Save surface projection of fault segment width
        source(isource).extra(i,3)=real(azi)                                    !Save azimuth of fault segment width
        source(isource).extra(source(isource).nnode,1)=
     1          source(isource).extra(source(isource).nnode,1)+
     2          source(isource).extra(i,1)                                      !Save the total length of fault in the last node storage (note the number of fault segments is one unit less than node points, and the last node point storage is saved for total fault length
      enddo
      return
      end



c
c     This subroutine generates magnitudes for source# 'isource'. It also c
c     ounts the number of generated magnitudes and stores them in 'nmag' an
c     d writes the magnitudes in proper column in array 'maglis'. In additi
c     on it saves the currently selected Beta, N0, and Mx values in proper
c     columns of maglis.
c     This subroutine calls the subroutines 'b_N0_mmax_gen' and 'maggen' du
c     ring its application.
c     The idea behind making this subroutine is to generate 'ndraw' magnitu
c     des and throw away the small ones (We hardwired the value of 'ndraw' 
c     equal to 100000 and this explains why we chose this number of rows fo
c     r maglis and iterations for the main loop in this subroutine. This id
c     ea seems necessary because there could be conditions where no catalog
c     could be generated for a set of seismicity parameters and threshold m
c     agnitudes. So, this idea forces the algorithm to generate 100000 even
c     ts and remove the ones that are les than threshold level. To do this,
c     on each earthquake to be generated, a random draw of seismicity param
c     eters is taken (using the subroutine 'b_N0_mmax_gen') and a minimum m
c     agnitude based on truncated Gutenberg-Richter distribution for 100000 
c     events is calculated. Then this minimum magnitude along with the draw
c     n seismicity parameters is sent random magnitude generation subroutin
c     e ('maggen') and the generated magnitude is kept if it were above the
c     threshold magnitude. This whole process in repeated ndraw=100000 time
c     s and the resulting list of magnitudes is kept.
c
      subroutine magmaker(isource,nmag,maglis)
      use           srcdef
      real          maglis(100000,10)
      real          b,N0,m0,mx,xmx(1)                                           !Selected b, N0, m0, Mmax for current draw
      real          ndraw
      if(source(isource).irecurr.le.3)then
          ndraw=100000.
          nmag=0
          call b_N0_mmax_gen(isource,b,N0,mx)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          do i=1,int(ndraw)
              if(epist.eq.1) call b_N0_mmax_gen(isource,b,N0,mx)                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              m0=-log(ndraw/N0+exp(-log(10.)*b*mx))/(b*log(10.))
              if(source(isource).mmin.lt.m0)source(isource).iflg(1)=1           !This condition checks if the choice of mmin is enough for making complete list of events per zone
              call maggen(b,m0,mx,1,0,source(isource).irecurr,xmx)
              xm=xmx(1)
              if(xm.ge.source(isource).mmin)then
                  if(source(isource).irecurr.eq.1.and.xm.gt.mx)xm=mx
                  nmag=nmag+1
                  maglis(nmag,5)=xm
                  maglis(nmag,6)=isource
                  maglis(nmag,8)=b*log(10.)
                  maglis(nmag,7)=N0/numyears
                  maglis(nmag,9)=mx
                  maglis(nmag,10)=m0
              endif
          enddo
      elseif(source(isource).irecurr.eq.4)then                                  !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          ndraw=int(source(isource).b_n0_val(1,2)*numyears)                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          nmag=0                                                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          call b_N0_mmax_gen(isource,b,N0,mx)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          do i=1,int(ndraw)                                                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              if(epist.eq.1) call b_N0_mmax_gen(isource,b,N0,mx)                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call maggen_table_CDF(isource,xmx)                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              xm=xmx(1)                                                         !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              if(xm.ge.source(isource).mmin.and.                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
     1           xm.le.mx)then                                                  !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  nmag=nmag+1                                                   !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,5)=xm                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,6)=isource                                        !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,8)=b*log(10.)                                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,7)=N0/numyears                                    !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,9)=mx                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,10)=0                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              endif                                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          enddo                                                                 !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
      elseif(source(isource).irecurr.eq.5)then                                  !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          ndraw=int(source(isource).b_n0_val(1,2)*numyears)                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          nmag=0                                                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          call b_N0_mmax_gen(isource,b,N0,mx)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          do i=1,int(ndraw)                                                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              if(epist.eq.1) call b_N0_mmax_gen(isource,b,N0,mx)                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call maggen_table_PMF(isource,xmx)                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              xm=xmx(1)                                                         !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              if(xm.ge.source(isource).mmin.and.                                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
     1           xm.le.mx)then                                                  !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  nmag=nmag+1                                                   !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,5)=xm                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,6)=isource                                        !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,8)=b*log(10.)                                     !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,7)=N0/numyears                                    !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,9)=mx                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                  maglis(nmag,10)=0                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
              endif                                                             !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
          enddo                                                                 !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
      endif
      
      return
      end


c
c     This subroutine does a random draw of seismicity parameters (b,N0,mx)
c     from the parameter lists of source zones considering given probabilit
c     ies. 
c
      subroutine b_N0_mmax_gen(isource,b,N0,mx)
      use           srcdef
      real          b,N0,mx                                                     !Selected Mmax for current draw
      real          mxmin,mxmax
c
c     This section does a random draw of N0 and Beta values.
c
      call randnumgen(rnd,0)
      do i=1,source(isource).nb_n0_val
          if(sum(source(isource).b_n0_val(1:i,3)).ge.rnd)exit		
      enddo
      b=source(isource).b_n0_val(i,1)/log(10.)                                  !Conversion of beta to b
      N0=numyears*source(isource).b_n0_val(i,2)                                 !Conversion of N0 to Nyears
c
c     This section does a random draw of Mmax values. If the nmmax value of
c     a source is equal to zero, the drawn Mmax value will be taken from a 
c     uniform distribution between two values.
c
      call randnumgen(rnd,0)
      if(source(isource).nmmax.eq.0)then
          mxmin=source(isource).mmax(1,1)
          mxmax=source(isource).mmax(2,1)
          mx=(mxmax-mxmin)*rnd+mxmin
      else
          do i=1,source(isource).nmmax
              if(sum(source(isource).mmax(1:i,2)).ge.rnd)exit		
          enddo
          mx=source(isource).mmax(i,1)
      endif

      return
      end


c
c     This subroutine fills equal distance grid points in a source zone. Th
c     e spacing between successive grid points along EW and NS is 'grd_spac
c     e'. This solution considers a spherical earth with 6378.137 km radius.
c     This subroutine is not implemented in 2010/05/27 version of program b 
c     ut it is reserved for future improvements or modifications.
c
      subroutine fillgrid_()
      use               srcdef
      r=6378.137
      do i=1,nsource
          if(source(i).typ.eq.2) cycle                                          !This line added to disable grid generation for faults (Added on 2010/09/30)
          call zoneRange(i,xmin,xmax,ymin,ymax)
          source(i).ngrid=0
          y=ymin
          do while (y.le.ymax)	
              x=xmin
              dx=rad2deg(source(i).grd_space/(r*cos(deg2rad(y))))
              dy=rad2deg(source(i).grd_space/r)
              do while (x.le.xmax)
                  call pnpoly(x,y,
     1                        source(i).nodes(1:source(i).nnode,2),
     2                        source(i).nodes(1:source(i).nnode,1),
     3                        source(i).nnode,inout)
                  if(inout.gt.-1)then
                      source(i).ngrid=source(i).ngrid+1
                      source(i).grids(source(i).ngrid,1)=y
                      source(i).grids(source(i).ngrid,2)=x
                  endif
                  x=x+dx
              enddo
              y=y+dy
          enddo
      enddo
      return
      end



c
c     This subroutine fills equal distance grid points in a source zone. Th
c     e spacing between successive grid points along EW and NS is 'grd_spac
c     e'. This solution considers the earth flattening with different polar
c     and equatorial radii. This subroutine is more complicated than 'fillg
c     rid_' and it is implemented in the program.
c
      subroutine fillgrid()
      use				srcdef
      do i=1,nsource
          if(source(i).typ.eq.2) cycle                                          !This line added to disable grid generation for faults (Added on 2010/09/30)
          call zoneRange(i,xmin,xmax,ymin,ymax)
          source(i).ngrid=0
          y=ymin
          do while (y.le.ymax)	
              x=xmin
              call radius_der(y,r,r_der)
              dx=rad2deg(source(i).grd_space/(r*cos(deg2rad(y))))
              dy=rad2deg(source(i).grd_space/(sqrt(r**2+r_der**2)))
              do while (x.le.xmax)
              call pnpoly(x,y,
     1					source(i).nodes(1:source(i).nnode,2),
     2					source(i).nodes(1:source(i).nnode,1),
     3					source(i).nnode,inout)
              if(inout.gt.-1)then
                  source(i).ngrid=source(i).ngrid+1
                  source(i).grids(source(i).ngrid,1)=y
                  source(i).grids(source(i).ngrid,2)=x
              endif
              x=x+dx
              enddo
          y=y+dy
          enddo
      enddo
      return
      end



c
c     This subroutine extracts the minimum and maximum limits of longitudes
c     and latitudes of a source zone.
c
      subroutine zoneRange(n,xmin,xmax,ymin,ymax)
      use				srcdef
      ymax=maxval(source(n).nodes(1:source(n).nnode,1))
      xmax=maxval(source(n).nodes(1:source(n).nnode,2))
      ymin=minval(source(n).nodes(1:source(n).nnode,1))
      xmin=minval(source(n).nodes(1:source(n).nnode,2))
      return
      end



c
c     Conversion from degrees to radians.
c     
      function deg2rad(x)
      pi=4.*atan(1.)
      deg2rad=x*pi/180.
      return
      end



c
c     Conversion from radians to degrees.
c
      function rad2deg(x)
      pi=4.*atan(1.)
      rad2deg=x*180./pi
      return
      end



c
c     This subroutine calculates the radius and the derivative of radius wi
c     th respect to polar angle changes of a point in latitude 'x'. This su
c     broutine considers a flattened earth with semi-major axis 'a' and sem
c     i-minor axis 'b'.
c
      subroutine radius_der(x,r,r_der)                                          !x is latitude in degrees, r is the earth's radius, r_der is derivative of r
      a=6378.1370
      b=6356.7523
      fi=deg2rad(x)
      r=sqrt(((cos(fi)*a**2)**2+(sin(fi)*b**2)**2)
     1     /((cos(fi)*a)**2+(sin(fi)*b)**2))
      r_der=0.5*((a**2*b**2*(b**2-a**2)*sin(2*fi))
     1		/(((a*cos(fi))**2+(b*sin(fi))**2)**2))/r
      return
      end


c
c     Subroutine "checkinput" checks sum of the weigts of seismicity, mmax,
c     and depth parameters. If weights of any parameters of a source don't
c     sum up to one(1), that source will be disabled from generating catalo
c     gues. Added on 2010/09/10.
c     
      subroutine checkinput()
      use srcdef
      epsilon=0.0001

      do i=1,nsource
          if(source(i).iuse.ne.0)then
              if(source(i).nb_n0_val.ne.0)then
                  w_of_seism=0.
                  do j=1,source(i).nb_n0_val
                  	  w_of_seism=w_of_seism+source(i).b_n0_val(j,3)
                  enddo
                  if(abs(w_of_seism-1).gt.epsilon)then
                  	  write(*,'(a18,a9,i2,a19)')'Seismicity weights',
     1				  ' of zone:',i,' don''t sum up to 1.'
                  	  source(i).iuse=0
                  endif
              endif

              if(source(i).nmmax.ne.0)then
                  w_of_mmax=0.
                  do j=1,source(i).nmmax
                  	  w_of_mmax=w_of_mmax+source(i).mmax(j,2)
                  enddo
                  if(abs(w_of_mmax-1).gt.epsilon)then
                  	  write(*,'(a21,i2,a19)')'Mmax weights of zone:',i,
     1				  ' don''t sum up to 1.'
                  	  source(i).iuse=0
                  endif
              endif

              if(source(i).ndepth.ne.0)then
                  w_of_depth=0.
                  do j=1,source(i).ndepth
                  	  w_of_depth=w_of_depth+source(i).depth(j,2)
                  enddo
                  if(abs(w_of_depth-1).gt.epsilon)then
                  	  write(*,'(a22,i2,a19)')'Depth weights of zone:',i,
     1				  ' don''t sum up to 1.'
                  	  source(i).iuse=0
                  endif
              endif
          endif		
      enddo
      
      return
      end


c
c     This subroutine opens the input file, reads its contents, and closes 
c     it. This subroutine will treat any line with its first character fro
c     m the following list as comment and any number of comment lines are
c     accepted in any place in input file. Character list: !@#$%^&*()
c
      subroutine readinput(input)
      use srcdef
      use smooth_seis
      integer(2)				delfiles 
      character*50 input,tempinput
      
      call getheadline(input,headline)                                          !Reading the header line of input. Added on 2011/10/26

      tempinput='tempinput.txt'
      call reformat(input,tempinput)
      open(unit=1,file=tempinput,status='old')

      read(1,*)numcat,numyears                                                  !Read number of cataloges and the number of years of each catalog
      read(1,*)nzones,nfaults                                                   !Read the number of zones (areal sources) and the number of faults (linear sources)
      read(1,*)seed, epist, ismooth                                             !Read the seed of random number generator. It sould be an integer. Zero means random seed from computer.
      if(ismooth.eq.1)read(1,*)cat_smo.catfile, cat_smo.parfile

      izone=0
      ifault=0
      nsource=nzones+nfaults
      if(nsource.ge.1)then
          do i=1,nsource
            read(1,*)tmp_src.name,tmp_src.iuse,tmp_src.irecurr,                 !Read source name, usage flag (iuse=0 => don't use this source), and 
     1               tmp_src.mmin                                               !recuurence relation type for catalog (1:GR, 2:Truncated GR, 3:Tapered GR), minimum magnitude
            if(tmp_src.irecurr.eq.4.or.tmp_src.irecurr.eq.5)then                !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                read(1,*)tmp_src.file_recurr                                    !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
                call read_file_recurr()                                         !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
            endif                                                               !Addition for reading probability distribution (CDF) of events in a source from file. This works if irecurr=4 (Added on 2012/03/07)
            read(1,*)tmp_src.typ,tmp_src.nnode,tmp_src.grd_space,               !Read source type (areal or linear),# of nodes, spacing between grid points,
     1               tmp_src.perturb(1),tmp_src.perturb(2)                      !horizontal and vertical perturbations in km (Added on 2010/09/28).

            tmp_src.b_n0_val=0.0                                                !Reset b_n0_val of tmp_src to zero and  next read source "i" beta and N0 values
            read(1,*)tmp_src.nb_n0_val                                          !Read number of source "i" beta-N0-weight sets
            if(tmp_src.nb_n0_val.gt.5.or.tmp_src.nb_n0_val.lt.0)then            !If the number of values is out of [0,5] range report error and ignore the source
                write(*,*)'Error with "b-value" of: ',tmp_src.name,i
                tmp_src.iuse=0
            elseif(tmp_src.nb_n0_val.eq.0)then                                  !If nb_n0_val=0 means linear distribution of beta and N0 values
                do j=1,2                                                        !Range of beta and N0 values (with uniform distribution in the range)
                    read(1,*)(tmp_src.b_n0_val(j,k),k=1,2)                      !Read i-th source, j-th set, beta(k=1) and N0(k=2) values
                enddo	
            else                                                                !If nb_n0_val is in [1,5] then:
                do j=1,tmp_src.nb_n0_val
                    read(1,*)(tmp_src.b_n0_val(j,k),k=1,3)                      !Read nb_n0_val beta values
                enddo
            endif

            tmp_src.mmax=0.0                                                    !Reset mmax of tmp_src to zero and next read source "i" Mmax values
            read(1,*)tmp_src.nmmax                                              !Read number of source "i" Mmax-weight sets
            if(tmp_src.nmmax.gt.5.or.tmp_src.nmmax.lt.0)then                    !If the number of values is out of [0,5] range report error and ignore the source
                write(*,*)'Error with "Mmax" of : ',tmp_src.name,i
                tmp_src.iuse=0
            elseif(tmp_src.nmmax.eq.0)then                                      !If nmmax=0 means linear distribution of Mmax values
                do j=1,2
                    read(1,*)tmp_src.mmax(j,1)                                  !Read two Mmax values (with uniform distribution in the range)
                enddo
            else                                                                !If nmmax is in [1,5] then:
                do j=1,tmp_src.nmmax
                    read(1,*)(tmp_src.mmax(j,k),k=1,2)                          !Read Mmax(k=1) and Mmax probability(k=2) values
                enddo
            endif

            tmp_src.depth=0.0                                                   !Reset depth of tmp_src to zero and next read source "i" depth values
            read(1,*)tmp_src.ndepth                                             !Read number of source "i" depth-weight sets
            if(tmp_src.ndepth.gt.5.or.tmp_src.ndepth.lt.0)then                  !If the number of values is out of [0,5] range report error and ignore the source
                write(*,*)'Error with "Depth" of : ',tmp_src.name,i
                tmp_src.iuse=0
            elseif(tmp_src.ndepth.eq.0)then                                     !If ndepth=0 means linear distribution of depth values
                do j=1,2
                    read(1,*)tmp_src.depth(j,1)                                 !Read range of depth values (with uniform distribution in the range)
                enddo
            else                                                                !If ndepth is in [1,5] then:
                do j=1,tmp_src.ndepth
                    read(1,*)(tmp_src.depth(j,k),k=1,2)                         !Read depth(k=1) and depth probability(k=2) values
                enddo
            endif

            tmp_src.geom=0.0
            if(tmp_src.typ.eq.2)then                                            !Read fault geometry if the source is linear source
                read(1,*)(tmp_src.geom(1,j),j=1,3)
!               do j=2,tmp_src.nnode
!                   do k=1,3
!                       tmp_src.geom(j,k)=tmp_src.geom(1,k)
!                   enddo
!               enddo
                if(tmp_src.nnode.gt.2)then
                    do k=2,tmp_src.nnode-1
                        tmp_src.geom(k,:)=tmp_src.geom(1,:)
                    enddo
                endif
            endif

            tmp_src.nodes=0.0
            do j=1,tmp_src.nnode                                                !Read source "i" node coordinates
                read(1,*)(tmp_src.nodes(j,k),k=1,2)
            enddo
            tmp_src_iflg=0

            source(i)=tmp_src                                                  !Assign read source parameters to source(i) (Added on 2010/09/29)

        enddo
      endif
      close(unit=1)
      delfiles=delfilesqq(tempinput)
      return
      end



c     
c     This subroutine reads the table of cumulative distribution from file
c     with sorted values of magnitude probability pairs. This subroutine is
c     added for the cases where the magnitude frequency distribution functi
c     on is not following predefined ones: GR, truncated GR, and tapered GR
c     . This case is called when irecurr is set to 4 for a source and it is
c     mainly intended for characteristic sources.
c
      subroutine read_file_recurr()
      use srcdef
      open(unit=5, file=tmp_src.file_recurr,status='old')
      read(5,*)tmp_src.n_tab_recurr
      do i=1,tmp_src.n_tab_recurr
          read(5,*)(tmp_src.tab_recurr(i,j),j=1,2)
      enddo
      close(unit=5)
      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*1			c1(10)
      character*50		    input,tempinput
      character*100		    line
      integer				flag
      data c1/'!','@','#','$','%','^','&','*','(',')'/
      open(unit=11,file=trim(input),status='old')
      open(unit=12,file=trim(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,'(a100)')line
          endif
      enddo
      close(unit=11)
      close(unit=12)
      return
      end



c
c     This this subroutine reads and returns the header line of the input f
c     ile. Added on 2011/10/26.
c
      subroutine getheadline(input,headline)
      character*50		input
      character*500     headline
      open(unit=13,file=input,status='old')
      read(13,'(a500)')headline
      close(unit=13)
      return
      end



c
c     This this subroutine writes input parameters on echo file
c
      subroutine writepar(input,output)
      use srcdef
      character*50      input,output
      character*500     lines(10000)
      integer           nline
      lines=''
      nline=0
      open(unit=14,file=input,status='old')
      open(unit=15,file=trim(output)//'.echo',status='unknown')
      do while (not(eof(14)))
          nline=nline+1
          read(14,'(a500)')lines(nline)
      enddo
      close(unit=14)
      write(15,'(a34,a12)')'EqHaz1 run date (yyyymmdd)  :     ',
     1                      real_clock(1)
      write(15,'(a34,a12)')'EqHaz1 run time (hhmmss.sss):     ',
     1                      real_clock(2)
      do i=1,nline
          write(15,'(a500)')lines(i)
      enddo
      close(unit=15)
      return
      end


c
c     The subroutine maggen_table generates a magnitude from given CDF tabl
c     e values. It refers to corresponding source CDF file readings.
c
      subroutine maggen_table_CDF(isource,xm)
      use srcdef
      dimension         xm(1)
      real              rnd
      call randnumgen(rnd,0)
      do i=1,source(isource).n_tab_recurr-1
          x1=source(isource).tab_recurr(i,1)
          y1=source(isource).tab_recurr(i,2)
          x2=source(isource).tab_recurr(i+1,1)
          y2=source(isource).tab_recurr(i+1,2)
          if(rnd.ge.y1.and.rnd.lt.y2)exit
      enddo
      
      xm(1)=(x2-x1)*(rnd-y1)/(y2-y1)+x1
      return
      end



c
c     The subroutine maggen_table generates a magnitude from given PMF tabl
c     e values. It refers to corresponding source PMF file readings.
c
      subroutine maggen_table_PMF(isource,xm)
      use srcdef
      dimension         xm(1)
      real              rnd, prob_sum
      prob_sum=0.
      call randnumgen(rnd,0)
      do i=1,source(isource).n_tab_recurr
          prob_sum=prob_sum + source(isource).tab_recurr(i,2)
          if(rnd.le.prob_sum)exit
      enddo
      xm(1)=source(isource).tab_recurr(i,1)
      return
      end



c
c     The subroutine maggen runs subroutines GR, TR, and TP. The parameters
c     are: (i) b the Gutenberg-Richter slop of the source; (ii) m0 the comp
c     leteness magnitude of catalogue; (iii) mmax the maximum magnitude for
c     truncated and tapered G-R distribution; (iv) npoint the number of mag
c     nitudes desired to be generated; (v) iseed is the seed value for rand
c     om number generator, it is supposed to be an integer number set by us
c     er. If it is set to an identical  non-zero value in different calls, 
c     the random number engine will generate the same sequence of pseudo-ra
c     ndom numbers. If it is set to zero, the sequences will differ from on
c     e set to another; (vi) imethod is a flag accepting the values 1, 2, a
c     nd 3 for calling GR, TR, and TP subroutines respectively; (vii) is a
c     size npoint vector that the generated magnitudes are stored in it and
c     sent back to calling program.
c
      subroutine maggen(b,m0,mmax,npoint,iseed,imethod,xm)
      dimension xm(npoint)
      real m0,mmax
      select case (imethod)
          case (1)
              call GR(b,m0,mmax,npoint,iseed,xm)
          case (2)
              call TR(b,m0,mmax,npoint,iseed,xm)
          case (3)
              call Tp(b,m0,mmax,npoint,iseed,xm)
          case default
              write(*,*)'Improper imethod value selected: ',imethod
      end select
      end



c
c     This subroutine generates magnitudes of n events following tapered GR
c     distribution with parameters b,m0,mmax and writes these magnitudes in
c     array XM. Thecalculation is based on random number rejection method u
c     sing probabiluty distribution function of tapered GR in Frohlich 2007
c     The paper by Frohlich presented cumulative sum of tapered GR for mome
c     nts. Karen used equation "m=c*log(M0)+d" to get a functional form for
c     density distribution function by differentiation. So the value of "c"
c     in this subrouutine is taken equal to 2/3 based on comparison of Froh
c     lich's equation given for GR distribution of moments and the one deri
c     ved by writing GR equation for magnitudes. Be careful on use of 'b' v
c     alue as it may mean slightly different thing in tapered GR than other
c     distribution.
c
      subroutine Tp(b,m0,mmax,n,iseed,XM)
      real m,m0,mmax
      dimension XM(n)
      c=2./3.
      u=10.**(m0*b)
      v=10.**(m0/c)
      w=10.**(mmax/c)
      pmax=log(10.)*(b+v/(w*c))
      i=1
          do while (i.le.n)
              call randnumgen(rnd,iseed)
              m=m0+(2*mmax-m0)*rnd
              p=u*alog(10.)*10.**(-b*m)*(b+10.**(m/c)/(w*c))
     1		  *exp((v-10.**(m/c))/w)
              call randnumgen(rnd2,iseed)
              rnd2=rnd2*pmax
              if(rnd2.le.p)then
                  XM(i)=m
                  i=i+1
              endif
          end do
      return
      end	
      
      
      
c
c     This subroutine generates magnitudes of n events following truncated 
c     GR distribution with parameters b,m0,mmax and writes these magnitudes 
c     in array XM. The calculation is based on the application of uniformly 
c     distributed random number to the inverse cumulative distribution func
c     tion.
c
      subroutine TR(b,m0,mmax,n,iseed,XM)
      real m0,mmax
      dimension XM(n)
      do i=1,n
          call randnumgen(rnd,iseed)
          XM(i)=m0-(log10(1-rnd*(1-10**(-b*(mmax-m0)))))/b
      enddo
      return
      end	
      
      
      
c
c     This subroutine generates magnitudes of n events following ordinary G
c     R distribution with lower cutoff for parameters b,m0,mmax and writes  
c     these magnitudesin array XM. The calculation is based on the applicat 
c     ion of uniformly distributed random number to the inverse cumulative 
c     distribution function.
c
      subroutine GR(b,m0,mmax,n,iseed,XM)
      real m0,mmax
      dimension XM(n)
      x=mmax
          do 10 i=1,n
              call randnumgen(rnd,iseed)
              XM(i)=m0-log10(1-rnd)/b
10		continue
      return
      end



c
c     This subroutine is not implemented in this version of program. It is 
c     reserved for future improvements or modifications.
c     The subroutine boxmul generates size "n" array "arr" and fills it wit
c     h standard normal values (mean=0 and standard deviation =1).
c
      subroutine boxmul(arr,n)
      dimension arr(n)
      pi=4*atan(1.)
          do i=1,n
              u1=0.
              u2=0.
              do while (u1.eq.0.or.u2.eq.0.)
                  call randnumgen(u1,0)
                  call randnumgen(u2,0)
              enddo
              arr(i)=sqrt(-2*log(u1))*cos(2*pi*u2)
          enddo
      return
      end



c
c     Modification of W. Randolph Franklin's code based on Shimrat 1962 for
c     checking if a point lies in a ploygon. px and py are x and y coordina
c     tes of the point in question and xx and yy are one dimensional x and 
c     y coordinate arrays of polygon, n is the number of polygon nodes, and 
c     inout is a flag with values 1, 0, and -1 for the point inside, on, an
c     d outside polygon respectively. The original code is downloaded from:
c     www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
c
      subroutine pnpoly(px,py,xx,yy,n,inout)
      real x(n),y(n),xx(n),yy(n)
      logical mx,my,nx,ny
6     do 1 i=1,n
      x(i)=xx(i)-px
1     y(i)=yy(i)-py
      inout=-1 
      do 2 i=1,n 
      j=1+mod(i,n) 
      mx=x(i).ge.0.0 
      nx=x(j).ge.0.0 
      my=y(i).ge.0.0 
      ny=y(j).ge.0.0 
      if(.not.((my.or.ny).and.(mx.or.nx)).or.(mx.and.nx)) go to 2 
      if(.not.(my.and.ny.and.(mx.or.nx).and..not.(mx.and.nx))) go to 3 
      inout=-inout 
      go to 2 
3     if((y(i)*x(j)-x(i)*y(j))/(x(j)-x(i))) 2,4,5 
4     inout=0 
      return 
5     inout=-inout 
2     continue 
      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



c
c     This subroutine generates a pseudo-random number upon its call and re
c     turns it in argument 'rnd'. The argument 'itype' introduces the seed 
c     of random number generator. This argument can be any integer number a
c     nd if it is set to zero, the seed value will be initialized from comp
c     uter's internal clock. However, the seed can be set only once, in the
c     first call of the subroutine and setting 'itype' in successive calls 
c     will not affect the seed value while the program is running. This fea
c     ture is implemented for practical purposes.
c     This subroutine generates 1000000 uniformly distributed pseudo-random 
c     numbers and stores them in array 'ranlist'. It stores the value of la
c     st number's seed too. Then upon each call of the subroutine the succe
c     ssive numbers from 'ranlist' are returned to the calling function. Up
c     on the 1000000th call, the subroutineuses the stored seed to generate 
c     another 1000000 random numbers and repeats the same procedure. The va
c     riables of this subroutine are chosen to be static, so their values w
c     ill remain untouched while the program is running and only the succes
c     sive calls of this subroutine changes their values.
c
      subroutine randnumgen(rnd,itype)
      logical, static     ::		visit
      real, static     ::		ranlist(1000000)
      integer, static     ::		count
      integer,static     ::		tmp(2)

      if(not(visit))then
          visit=not(visit)
          count=0
          tmp(1)=itype
          tmp(2)=itype
      endif

      if(mod(count,1000000).eq.0)then
          if(tmp(1).eq.0)then
              call random_seed()
          else
              call random_seed(put=tmp(1:2))
          endif
          call random_number(ranlist)
          call random_seed(get=tmp(1:2))
          count=0
      endif

          count=count+1
          rnd=ranlist(count)
      return
      end
      
      
      
c
c     The subsequent subroutines and functions are used for generating area
c     around linear sources so that the program becomes capable of generati
c     ng catalogs for linear sources. Note that EqHaz1 considers the geomet
c     ry of multisegment faults by accepting their node points, width, dept
c     h of the top of the fault, and dip angle.
c
      subroutine line2zone(nodes,geom,nnode)
      real nodes(50,2),geom(50,3),tmp(50,2)
      double precision p1(2),p2(2),q1(2),q2(2)
      double precision x1(2),x2(2),x3(2),x(2),width
      double precision rad
      nloop=nnode-2
      if(nloop.eq.0)then
          x1(1)=dble(nodes(1,1))		
          x1(2)=dble(nodes(1,2))		
          x2(1)=dble(nodes(2,1))		
          x2(2)=dble(nodes(2,2))
          width=dble(geom(1,1))*dcos(dble(rad(dble(geom(1,2)))))
          call find_conj_points(x1,x2,p1,p2,width)
          nodes(3,1)=p2(1)
          nodes(3,2)=p2(2)
          nodes(4,1)=p1(1)
          nodes(4,2)=p1(2)
          nnode=4
          return
      endif
      n=0                                                                       !n is the node number in temporary array
      do i=1,nloop
          x1(1)=dble(nodes(i,1))
          x1(2)=dble(nodes(i,2))
          x2(1)=dble(nodes(i+1,1))
          x2(2)=dble(nodes(i+1,2))
          x3(1)=dble(nodes(i+2,1))
          x3(2)=dble(nodes(i+2,2))
          if(i.eq.1)then
              width=dble(geom(1,1))*dcos(dble(rad(dble(geom(1,2)))))
              call find_conj_points(x1,x2,p1,p2,width)
              n=n+1
              tmp(n,1)=p1(1)
              tmp(n,2)=p1(2)
          endif
          n=n+1
          width=dble(geom(i,1))*dcos(dble(rad(dble(geom(i,2)))))
          call find_conj_points(x1,x2,p1,p2,width)
          width=dble(geom(i+1,1))*dcos(dble(rad(dble(geom(i+1,2)))))
          call find_conj_points(x2,x3,q1,q2,width)
          call intersect(p1,p2,q1,q2,x(1),x(2))
          tmp(n,1)=x(1)
          tmp(n,2)=x(2)
          if(i.eq.nloop)then
              n=n+1
              tmp(n,1)=q2(1)
              tmp(n,2)=q2(2)
          endif
      enddo
      do i=1,n
          nodes(nnode+i,1)=tmp(n+1-i,1)
          nodes(nnode+i,2)=tmp(n+1-i,2)
      enddo
      nnode=nnode+n
      return
      end


      subroutine find_conj_points(p1,p2,cp1,cp2,width)
      double precision p1(2),p2(2),cp1(2),cp2(2)
      double precision azi,azi_perp,width,dist
      call xy2disazi(p1(1),p1(2),p2(1),p2(2),dist,azi)
      azi_perp=azi+dble(90.)
      call disazi2xy(p1(1),p1(2),width,azi_perp,cp1(1),cp1(2))
      call disazi2xy(p2(1),p2(2),width,azi_perp,cp2(1),cp2(2))
      return
      end


c     The subroutine intersect gets the coordinates of two pairs of points (in d
c     egrees) and generates the coordinates of a point resulting from intersecti
c     ng two lines containing the pairs of points. 
      subroutine intersect(p1,p2,q1,q2,xlat,xlon)
      double precision p1(2),p2(2),q1(2),q2(2)
      double precision vp1(3),vp2(3),vq1(3),vq2(3)
      double precision vp12(3),vq12(3),vec(3)
      double precision xlat,xlon,xlat1,xlon1,xlat2,xlon2
      double precision dist,dist1,dist2,azi,azi1,azi2
      double precision deg, rad

      do i=1,2
          p1(i)=rad(p1(i))
          p2(i)=rad(p2(i))
          q1(i)=rad(q1(i))
          q2(i)=rad(q2(i))
      enddo

      call xy2vec(p1(1),p1(2),vp1)
      call xy2vec(p2(1),p2(2),vp2)
      call xy2vec(q1(1),q1(2),vq1)
      call xy2vec(q2(1),q2(2),vq2)

      call mulvec(vp1,vp2,vp12)
      call mulvec(vq1,vq2,vq12)
      call mulvec(vp12,vq12,vec)

      call vec2xy(vec,xlat1,xlon1)
      call vec2xy(-vec,xlat2,xlon2) 

      xlat1=deg(xlat1)
      xlon1=deg(xlon1)
      xlat2=deg(xlat2)
      xlon2=deg(xlon2)
      do i=1,2
          p1(i)=deg(p1(i))
          p2(i)=deg(p2(i))
          q1(i)=deg(q1(i))
          q2(i)=deg(q2(i))
      enddo

      call xy2disazi(p1(1),p1(2),xlat1,xlon1,dist1,azi)
      call xy2disazi(p1(1),p1(2),xlat2,xlon2,dist2,azi)
      xlat=xlat1
      xlon=xlon1

      if(dist1.gt.dist2)then
          xlat=xlat2
          xlon=xlon2
      endif

      call xy2disazi(p1(1),p1(2),p2(1),p2(2),dist,azi1)
      call xy2disazi(q1(1),q1(2),q2(1),q2(2),dist,azi2)
      if(abs(azi1-azi2).lt.dble(1.))then
          xlat=p2(1)
          xlon=p2(2)
      endif
      end


      subroutine mulvec(v1,v2,v)
      double precision v(3),v1(3),v2(3)
      v(1)=v1(2)*v2(3)-v1(3)*v2(2)
      v(2)=v1(3)*v2(1)-v1(1)*v2(3)
      v(3)=v1(1)*v2(2)-v1(2)*v2(1)
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     The main program runs the coordinate finding subroutines and transfers the
c     points to the proper location that gives the coordinates of the fault.
      subroutine array(xlat1,xlon1,xlat2,xlon2,dis,ndis,num,iflag)
      implicit double precision (a-h,o-z)
      dimension xy(10000,2),xy_all(100000,2),rotmat(3,3),vec(3)
      re=6378.137					!Radius of earth
      pi=dble(4)*datan(dble(1))
c      write(*,*)'Enter latitude and longitude of fault first corner:'
c      read(*,*)xlat1,xlon1
c      write(*,*)'Enter latitude and longitude of fault second corner:'
c      read(*,*)xlat2,xlon2
c      write(*,*)'Enter distance and distance sub-dvisions: '
c      read (*,*)dis,ndis
c      write(*,*)'Enter the number of points per track: '
c      read(*,*)num
c      write(*,*)'Enter (1) for equal azimuth or'
c      write(*,*)'      (2) for equal distance points:'
c      read(*,*)iflag
      open(unit=1,file='array.txt',status='unknown')
      call xy2disazi(xlat1,xlon1,xlat2,xlon2,flt_len,strike)
      do i=1,ndis
          dist=dble(i)*dis/dble(ndis)
          if(iflag.eq.1)then
              call track_eq_angle(flt_len,dist,num,xy)
          elseif(iflag.eq.2)then
              call track_eq_dist(flt_len,dist,num,xy)
          endif
          xy_all((i-1)*num+1:i*num,1)=xy(1:num,1)
          xy_all((i-1)*num+1:i*num,2)=xy(1:num,2)
      enddo
      do i=1,ndis*num
          call xy2vec(rad(xy_all(i,1)),rad(xy_all(i,2)),vec)
          call rot_matrix(flt_len/(2d0*re),rotmat)
          call mult_rot_matrix(rotmat,vec,'z')
          call rot_matrix(pi/2d0,rotmat)
          call mult_rot_matrix(rotmat,vec,'x')
          call rot_matrix(-rad(strike),rotmat)
          call mult_rot_matrix(rotmat,vec,'x')
          call rot_matrix(-rad(xlat1),rotmat)
          call mult_rot_matrix(rotmat,vec,'y')
          call rot_matrix(rad(xlon1),rotmat)
          call mult_rot_matrix(rotmat,vec,'z')
          call vec2xy(vec,xy_all(i,1),xy_all(i,2))
          xy_all(i,1)=deg(xy_all(i,1))
          xy_all(i,2)=deg(xy_all(i,2))
          write(1,'(2(e18.10,2x))')(xy_all(i,j),j=1,2)
      enddo
      close(unit=1)
      end
      

c     This subroutine finds the coordinates of num points on a track from fault 
c     with length flt_len at distance dist in array coord. The fault is oriented
c     in EW direction and centered in the geographic origin.
      subroutine track_eq_dist(flt_len,dist,num,coord)
      implicit double precision (a-h,o-z)
      dimension coord(10000,2),p1(2),p2(2),c1(2),c2(2)
      re=dble(6378.137)					!Radius of earth
      pi=dble(4)*datan(dble(1))
      call track_len(flt_len,dist,track,h_circle,h_line)
      call find_points(flt_len,dist,p1,p2,c1,c2)
      do i=1,num
          dist_on_track=dble(i-1)*track/dble(num)
          call segment_finder(dist_on_track,h_circle,h_line,
     1						nseg,dist_on_seg)
          call find_point(flt_len,dist,p1,p2,c1,c2,nseg,dist_on_seg,y,x)
          coord(i,1)=y
          coord(i,2)=x
      enddo
      return
      end


c     This subroutine finds the coordinates of a point on track, given line leng
c     th, distance of track from line, points c1, c2, p1, p2, segment number , a
c     nd the distance on the segment.
      subroutine find_point(flt_len,dist,p1,p2,c1,c2,
     1						nseg,dist_on_seg,y,x)
      implicit double precision (a-h,o-z)
      dimension coord(10000,2),p1(2),p2(2),c1(2),c2(2)
      re=dble(6378.137)					!Radius of earth
      pi=dble(4)*datan(dble(1))
      if(nseg.eq.1)then
      y=deg(p1(1))
      x=deg(p1(2))+deg(dist_on_seg/(re*dcos(dist/re)))
      return
      endif
      if(nseg.eq.2)then
      azi=deg(dist_on_seg/(re*dsin(dist/re)))
      call disazi2xy(deg(c2(1)),deg(c2(2)),dist,azi,y,x)
      return
      endif
      if(nseg.eq.3)then
      y=deg(p2(1))
      x=deg(p2(2))-deg(dist_on_seg/(re*dcos(dist/re)))
      return
      endif
      if(nseg.eq.4)then
      azi=deg(dist_on_seg/(re*dsin(dist/re))+pi)
      call disazi2xy(deg(c1(1)),deg(c1(2)),dist,azi,y,x)
      return
      endif
      return
      end


c     This subroutine finds the segment number and distance on that segment, giv
c     en distance clockwise from p1 and length of line and semicircles. (See fig
c     ure 3).
      subroutine segment_finder(dist_on_track,h_circle,h_line,
     1					nseg,dist_on_seg)
      implicit double precision (a-h,o-z)
      nseg=1
      if(dist_on_track.gt.h_line)nseg=2
      if(dist_on_track.gt.(h_line+h_circle))nseg=3
      if(dist_on_track.gt.(2d0*h_line+h_circle))nseg=4
      if(nseg.eq.1)dist_on_seg=dist_on_track
      if(nseg.eq.2)dist_on_seg=dist_on_track-h_line
      if(nseg.eq.3)dist_on_seg=dist_on_track-h_line-h_circle
      if(nseg.eq.4)dist_on_seg=dist_on_track-2d0*h_line-h_circle
      return
      end


c     This subroutine finds the coordinates of the points c1, c2, p1, p2 in radi
c     ans. (See Figure03).
      subroutine find_points(flt_len,dist,p1,p2,c1,c2)
      implicit double precision (a-h,o-z)
      dimension p1(2),p2(2),c1(2),c2(2)
      re=dble(6378.137)					!Radius of earth
      pi=dble(4)*datan(dble(1))
      c1(1)=0d0							!Latitude (rad)
      c1(2)=-flt_len/(2*re)				!Longitude (rad)
      c2(1)=0d0							!Latitude (rad)
      c2(2)=+flt_len/(2*re)				!Longitude (rad)
      p1(1)=dist/re						!Latitude (rad)
      p1(2)=-flt_len/(2*re)				!Longitude (rad)
      p2(1)=-dist/re					!Latitude (rad)
      p2(2)=+flt_len/(2*re)				!Longitude (rad)
      return
      end


c     This subroutine finds the length of the track on a spherical earth.
      subroutine track_len(flt_len,dist,track,h_circle,h_line)
      implicit double precision (a-h,o-z)
      re=dble(6378.137)					!Radius of earth
      pi=dble(4)*datan(dble(1))
      h_line=flt_len*dcos(dist/re)
      h_circle=pi*re*dsin(dist/re)
      track=2d0*(h_line+h_circle)
      return
      end


c     This subroutine fills the coordinates of points (in degrees) along a track
c     , given a fault length (in km), distance (in km), and number of points on
c     track. The fault is assumed to be in EW direction centered on geographical
c     origin. This subroutine finds the coordinates in equal angles relative to 
c     the fault center.
      subroutine track_eq_angle(flt_len,dist,num,coord)
      implicit double precision (a-h,o-z)
      dimension coord(10000,2)
      re=dble(6378.137)					!Radius of earth
      pi=dble(4)*datan(dble(1))
      teta=find_critic_teta(flt_len,dist)
      do i=1,num
          azi=dble(i-1)*360d0/dble(num)
          b=dis2angle(dist)
          bb=calc_bb(azi)
          c=dis2angle(flt_len/2d0)
          if(i_zone(azi,teta).eq.1)then
              a=a_in_zone1(bb,b)
          else
              a=a_in_zone2(bb,b,c)
          endif
          call disazi2xy(dble(0),dble(0),a*re,azi,coord(i,1),coord(i,2))
      enddo
      end


c     This function calculates the angle between line and the connector of the 
c     center of the line to the point (BB in radians) given the azimuth of the
c     point (in radians).
      function calc_bb(azi)
      implicit double precision (a-h,o-z)
      if(azi.ge.0d0.and.azi.lt.90d0)calc_bb=90d0-azi
      if(azi.ge.90d0.and.azi.lt.180d0)calc_bb=azi-90d0
      if(azi.ge.180d0.and.azi.lt.270d0)calc_bb=270d0-azi
      if(azi.ge.270d0.and.azi.lt.360d0)calc_bb=azi-270d0
      calc_bb=rad(calc_bb)
      return
      end
      

c     This function gives the distance of a point from the center of a line (in 
c     radians) based on its distabce from tip of the line (b in radians), the a
c     ngle between line and the connector of the center of the line to the poin
c     t (bb in radians), and the half-length of line (c). This subroutine does 
c     the calculations for points in zone 2. (See Figure 2).
      function a_in_zone2(bb,b,c)
      implicit double precision (a-h,o-z)
      pi=dble(4)*datan(dble(1))
      if(bb.eq.90d0.or.c.eq.0d0)then
          a_in_zone2=b
          return
      endif
      term1=dcos(b)/dsqrt(dcos(bb)**2*dsin(c)**2+dcos(c)**2)
      term2=dcos(c)/(dcos(bb)*dsin(c))
      if(term1.gt.1d0)then
          term1=1d0
      elseif(term1.lt.-1d0)then
          term1=-1d0
      endif
      a_in_zone2=pi-dasin(term1)-datan(term2)
      return
      end


c     This function gives the distance of a point from the center of a line (in 
c     radians) based on its distance from the line (b in radians) and the angle
c     between line and the connector of the center of the line to the point (BB 
c      in radians). This subroutine does the calculations for points in zone 1.
c     (See Figure01).
      function a_in_zone1(bb,b)
      double precision a_in_zone1,b,bb,pi
      pi=dble(4)*datan(dble(1))
      if(dsin(bb).eq.0d0)then
          a_in_zone1=b
          return
      endif
      if(dsin(b)/dabs(dsin(bb)).gt.1d0)then
          a_in_zone1=pi/2d0
          return
      endif
      if(dsin(b)/dabs(dsin(bb)).lt.-1d0)then
          a_in_zone1=-pi/2d0
          return
      endif
      a_in_zone1=dasin(dsin(b)/dabs(dsin(bb)))
      return
      end


c     This function decides if a point with given azimuth and critical angle is
c     in zone 1 or zone 2. Zone 1 is the region in which a point is along a para
c     llel and zone 2 is the region in which a point is along a circle. Check th
c     e explanations for function find_critic_teta.
      function i_zone(azi,teta)
      double precision azi,teta
      i_zone=1
      if(azi.ge.teta)i_zone=2
      if(azi.ge.(180d0-teta))i_zone=1
      if(azi.ge.(180d0+teta))i_zone=2
      if(azi.ge.(360d0-teta))i_zone=1
      return
      end


c     This function finds the critical azimuth associated with an EW line (with
c     its center on origine) and a distance. The critical angle is the one below
c     which equidistant points are along a parallel and beyond that they are on
c     circle. The fault length and distanc are in km and the function value is i
c     n degrees. 
      function find_critic_teta(flt_len,dist)
      implicit double precision (a-z)
      pi=dble(4)*datan(dble(1))
      b=dis2angle(dist)
      c=dis2angle(flt_len/dble(2))
      a=dacos(dcos(b)*dcos(c))
      find_teta=(pi/2.d0)-dasin(dsin(b)/dsqrt(1-dcos(b)**2*dcos(c)**2))
      find_critic_teta=deg(find_teta)
      return
      end


c     This subroutine finds coordinates of the point given its azimuth (azi in d
c     egrees) and distance (in km) from a point. The coordinates of points are 
c     in degrees.
      subroutine disazi2xy(lat1,lon1,dist,azi,lat2,lon2)
      implicit double precision (a-z)
      dimension vec(3),mat(3,3)
c     write(*,*)'**',lat1,lon1,dist,azi
      lat1=rad(lat1)
      lon1=rad(lon1)
      azi=rad(azi)
      vec=dble(0)
      vec(1)=dcos(dis2angle(dist))
      vec(3)=dsin(dis2angle(dist))
      call rot_matrix(-azi,mat)
      call mult_rot_matrix(mat,vec,'x')
      call rot_matrix(-lat1,mat)
      call mult_rot_matrix(mat,vec,'y')
      call rot_matrix(lon1,mat)
      call mult_rot_matrix(mat,vec,'z')
      call vec2xy(vec,lat2,lon2) 
      lat1=deg(lat1)
      lon1=deg(lon1)
      lat2=deg(lat2)
      lon2=deg(lon2)
      azi=deg(azi)
c     write(*,*)'**',lat1,lon1,dist,azi,lat2,lon2
      return
      end


c     This subroutine finds the distance (dist in km) and azimuth (azi in degree
c     s) of a second point relative to first point given the coordinates of two
c     points in degrees. 
      subroutine xy2disazi(lat1,lon1,lat2,lon2,dist,azi)
      implicit double precision (a-z)
      dimension vec2(3),v_tmp(3),mat(3,3)
      re=6378.137					!Radius of earth
      lat1=rad(lat1)
      lon1=rad(lon1)
      lat2=rad(lat2)
      lon2=rad(lon2)
      arg=dcos(lat1)*dcos(lat2)*dcos(lon1-lon2)+dsin(lat1)*dsin(lat2)
      if(arg.lt.dble(-1))arg=dble(-1)
      if(arg.gt.dble(+1))arg=dble(+1)
      dist=re*dacos(arg)
      call xy2vec(lat2,lon2,vec2)
      call rot_matrix(-lon1,mat)
      call mult_rot_matrix(mat,vec2,'z')
      call rot_matrix(+lat1,mat)
      call mult_rot_matrix(mat,vec2,'y')
      call findazi(vec2,azi)
      lat1=deg(lat1)
      lon1=deg(lon1)
      lat2=deg(lat2)
      lon2=deg(lon2)
      azi=deg(azi)
      return
      end


c     This subroutine calculates the azimuth of a point from geographical origin
c     (0,0) given unit vector components of the point. The azimuth is in radians
c     and it is the angle betweenprime meridian and half-great circle connecting
c     the origin and the point.
      subroutine findazi(vec,azi)
      double precision vec(3),azi,pi
      pi=dble(4)*datan(dble(1))
      if(vec(3).eq.dble(0))then
          if(vec(2).eq.dble(0.))azi=0
          if(vec(2).gt.dble(0.))azi=pi/2.
          if(vec(2).lt.dble(0.))azi=3.*pi/2.
          return
      endif
      azi=datan(vec(2)/vec(3))
      if(vec(3).lt.dble(0.))azi=azi+pi
      if(azi.lt.dble(0.))azi=azi+2*pi
      return
      end
      

c     This subroutine converts the x, y, and z of a unit vector to latitude and
c     longitude values. Latitude and longitude are in radians.
      subroutine vec2xy(vec,lat,lon) 
      implicit double precision (a-z)
      dimension vec(3)
      pi=dble(4)*datan(dble(1))
      arg=0.
      if(vec(1).ne.0.or.vec(2).ne.0.or.vec(3).ne.0.)then
          arg=vec(3)/dsqrt(vec(1)**2+vec(2)**2+vec(3)**2)
      endif
          if(arg.gt.dble(1.))arg=dble(1.)
          if(arg.lt.dble(-1.))arg=dble(-1.)
      lat=dasin(arg)

      arg=0.
      if(vec(1).ne.0.or.vec(2).ne.0.)then
          arg=vec(2)/dsqrt(vec(1)**2+vec(2)**2)
      endif
          if(arg.gt.dble(1.))arg=dble(1.)
          if(arg.lt.dble(-1.))arg=dble(-1.)
      lon=dasin(arg)
      if(lon.lt.dble(0.).and.vec(1).lt.dble(0.))lon=-pi-lon
      if(lon.ge.dble(0.).and.vec(1).lt.dble(0.))lon= pi-lon
      return	
      end


c     This subroutine converts latitude and longitude of a vector to x, y, and z
c     components of a unit vector. The latitude and longitude are in radians.
      subroutine xy2vec(lat,lon,vec)
      double precision lat,lon,vec(3)
      vec(1)=dcos(lat)*dcos(lon)
      vec(2)=dcos(lat)*dsin(lon)
      vec(3)=dsin(lat)
c     write(*,*)'XX',lat,lon,vec
      return
      end


c     This subroutine derives the vector resulted from rotating the vector vec a
c     round one of x, y, or z axes using rotation matrix mat. Note that three su
c     broutines mult_rot_matrix, shift_vec, and rot_matrix are working doing rot
c     ations, by playing a trick. Instead of deriving rotation matrix for each a
c     xis, they switch components of the vector, do the rotation around z-axis a
c     nd switch the components of the vector back to its initial state.
c     The rotation matrix is 'mat', the vector is 'vec', and the axis of rotatio
c     n is 'c' that takes one of the values 'x', 'y', or 'z'.
      subroutine mult_rot_matrix(mat,vec,c)
      double precision mat(3,3),vec(3),vec_tmp(3)
      character*1 c
      if(c.eq.'x')call shift_vec(vec,1)
      if(c.eq.'y')call shift_vec(vec,2)
      call mult_mat_vec(mat,vec,vec_tmp)
      vec=vec_tmp
      if(c.eq.'x')call shift_vec(vec,2)
      if(c.eq.'y')call shift_vec(vec,1)
      return	
      end


c     Subroutine for shifting components of a vector, one unit up or down.
      subroutine shift_vec(vec,n)
      double precision vec(3),tmp(3)
      if(n.eq.-1)n=2
      do i=1,3
          tmp(i)=vec(max(mod(i+n,3),mod(i+n,4)))
      enddo
      vec=tmp
      return
      end


c     Subroutine for derivation of rotation matrix in xy plane around z-axis in
c     counter-clockwise direction. Angle of rotation is alpha in radians from po
c     sitive x to positive y.
      subroutine rot_matrix(alpha,mat)
      double precision alpha,mat(3,3)
      mat=dble(0)
      mat(1,1)= dcos(alpha)
      mat(1,2)=-dsin(alpha)
      mat(2,1)= dsin(alpha)
      mat(2,2)= dcos(alpha)
      mat(3,3)= dble(1)
      return
      end


c     Subroutine for multiplying a vector (vec_in) by a matrix (mat) and store
c     the results in another vector (vec_out).
      subroutine mult_mat_vec(mat,vec_in,vec_out)
      double precision mat(3,3),vec_in(3),vec_out(3)
      do i=1,3
          vec_out(i)=dble(0)
          do j=1,3
              vec_out(i)=vec_out(i)+mat(i,j)*vec_in(j)
          enddo
      enddo
      return
      end


c     Function to convert degrees to radians.	
      function rad(x)
      double precision rad,pi,x
      pi=dble(4)*datan(dble(1.))
      rad=x*pi/dble(180)
      return
      end


c     Function to convert radians to degrees.
      function deg(x)
      double precision deg,pi,x
      pi=dble(4)*datan(dble(1.))
      deg=x*dble(180)/pi
      return
      end


c     Function to convert distance on earth to angle in radians.
      function dis2angle(dis)
      double precision dis2angle,dis,re
      re=dble(6378.137)                                                         !Earth's radius
      dis2angle=dis/re
      return
      end


      subroutine read_smooth_parfile()
      use                   srcdef
      use                   smooth_seis                                         !This module is used for reserving area and access to earthquake catalog values for smooth seismicity  approach (Last modified on 2012/07/06)
      integer(2)                                  ::delfiles
      integer                                     ::i,j
      character (len=50)                          :: cat_smo_par_tmp
      write(*,*)'Parametr file for smooth seismicity: ', cat_smo.parfile
      cat_smo_par_tmp='cat_smo_par.tmp'
c      write(*,*)'Number of sources is:', nsource
c      write(*,*)'I am 1..............'
      call reformat(cat_smo.parfile,cat_smo_par_tmp)
c      write(*,*)'I am 2..............'
      open(unit=17, file=cat_smo_par_tmp,status='unknown')
c      write(*,*)'I am 3..............'
      do i=1, nsource
c          write(*,*)'I am 4..............'
          read(17,*)
          read(17,*)source_smo(i).n_comp
c          write(*,*)'I am 5..............'
          do j=1,source_smo(i).n_comp
c              write(*,*)'I am 6..............'
              read(17,*)source_smo(i).comp_m(j),
     1                  source_smo(i).comp_yr(j)
          enddo
c          write(*,*)'I am 7..............'
!          read(17,*)source_smo(i).end_year
!          read(17,*)source_smo(i).beta
          read(17,*)source_smo(i).ring_w, source_smo(i).corr_dis,
     1              source_smo(i).bkground     
      enddo
c      write(*,*)'I am 8..............'
      close(unit=17)
      delfiles=delfilesqq(cat_smo_par_tmp)
      return
      end



      subroutine read_smooth_catfile()
      use                   smooth_seis                                         !This module is used for reserving area and access to earthquake catalog values for smooth seismicity  approach (Last modified on 2012/07/06)
      integer                                     ::nline
      character (Len=512),allocatable             ::content(:)                  !Temporary array of strings for reading the catalog file
      write(*,*)'Catalog file for smooth seismicity: ',cat_smo.catfile
      open(unit=16, file=trim(cat_smo.catfile),status='unknown')
      nline=0
      do while(not(eof(16)))
          read(16,*)
          nline=nline+1
      enddo
      rewind(unit=16)
      allocate(content(nline))
      do i=1,nline
          read(16,'(a512)')content(i)
      enddo
      do i=1,200
          if(i_strcmp(content(i),'yr   mo dy hr mi')>0)exit
      enddo
      cat_smo.nline_head=i+1
      cat_smo.n_rec=nline-cat_smo.nline_head
      do i=1,cat_smo.n_rec
          read(content(cat_smo.nline_head+i),*)(cat_smo.time(i,j),j=1,5)
     1                                    ,(cat_smo.coord(i,j),j=1,2)
     2                                    ,(cat_smo.mag(i,j),j=1,9)
     3                                    , cat_smo.Mpf(i),cat_smo.Tm(i)
     4                                    , cat_smo.coord(i,3)
     5                                    , cat_smo.dd(i)
     6                                    ,(cat_smo.czf(i,j),j=1,2)
     7                                    , cat_smo.mf(i), cat_smo.Mw(i)
      enddo
      close(unit=16)
      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



      subroutine fillgrid_smo()                                                 !This subroutine calculates the weights of each grid point output from fillgrid based on smooth seismicity approach.
      use                   srcdef
      use                   smooth_seis                                         !This module is used for reserving area and access to earthquake catalog values for smooth seismicity  approach (Last modified on 2012/07/06)
      integer                                    ::i, j, k, l, m, n_ring
      integer                                    ::inrange
      real                                       ::Wei_fac
      real                                       ::x, dR, pi
      real                                       ::tmp_dist, tmp_sum
      real                                       ::numer, denom
      real, dimension(100)                       ::eq_in_ring
      double precision                           ::dist, azi
      double precision                           ::x1, y1, x2, y2
      pi=4.*atan(1.)
      do i=1,nsource
          if(source(i).typ.eq.2) cycle                                          !This line added to disable grid generation for faults (Added on 2010/09/30)
!          call Weichert_fac(i,Wei_fac)
          call find_n_ring(i,n_ring)
          tmp_sum=0.
          do j=1,source(i).ngrid
              eq_in_ring=0.
              cat_loop: do k=1,cat_smo.n_rec
!                  x1=dble(source(i).grids(j,2))
!                  y1=dble(source(i).grids(j,1))
!                  x2=dble(cat_smo.coord(k,2))
!                  y2=dble(cat_smo.coord(k,1))
!                  write(53,'(4(D10.4,2x)i5,2x,a2)')y1,x1,y2,x2,k,'in'
!                  call xy2disazi(y1,x1,y2,x2,dist,azi)
!                  write(53,'(6(D10.4,2x))')y1,x1,y2,x2,dist, azi
                  call xy2disazi(dble(source(i).grids(j,1)),
     1                           dble(source(i).grids(j,2)),
     2                           dble(cat_smo.coord(k,1)),
     3                           dble(cat_smo.coord(k,2)),
     4                           dist,azi)
                  inrange=0
                  tmp_dist=real(dist)
                  if(tmp_dist.gt.(n_ring*source_smo(i).ring_w))cycle
                  do m=1,source_smo(i).n_comp
                      if(cat_smo.time(k,1).ge.source_smo(i).comp_yr(m).
     1                    and.cat_smo.Mw(k).ge.source_smo(i).comp_m(m))
     2                    then
                          inrange=1
                          exit
                      endif
                  enddo
                  if(inrange.eq.1)then
                      do l=1,n_ring
                          if(tmp_dist.le.(real(l)*source_smo(i).ring_w))
     1                        then
                              eq_in_ring(l)=eq_in_ring(l)+1.
                              cycle cat_loop
                          endif
                      enddo
                  endif
              enddo cat_loop
              numer=0.
              denom=0.
              dR=source_smo(i).ring_w
              do l=1,n_ring
                  x=dR*sqrt((real(l)-1.)**2.+real(l)-0.5)
                  numer=numer+eq_in_ring(l)*
     1                  exp(-(x**2.)/(source_smo(i).corr_dis**2.))/
     2                  (pi*(2*real(l)-1)*source_smo(i).ring_w)                 !Needs correction but it is unimportant as ratios matter, it should be "(pi*(2*real(l)-1)*source_smo(i).ring_w**2)
                  denom=denom+exp(-(x**2.)/(source_smo(i).corr_dis**2.))
              enddo
              source(i).grids_w(j)=numer/denom
              tmp_sum=tmp_sum+numer/denom
          enddo
          if(tmp_sum.eq.0.)then
              source(i).iuse=0
          endif
          source(i).grids_w=source(i).grids_w/tmp_sum
          if(source_smo(i).bkground.ne.0.)call set_background(i)
      enddo
      return
      end



      subroutine set_background(n)
      use                   srcdef
      use                   smooth_seis
      integer                                    :: i, n
      real                                       :: added_w
      if(source_smo(n).bkground.ge.1.or.
     1                 source_smo(n).bkground.le.0)then
          write(*,*)'Outside range background seismicity for zone#',n
          write(*,*)'No background seismicity applied.',
     1               source_smo(n).bkground
          return
      endif
      if(source(n).ngrid.eq.0)return
      added_w=source_smo(n).bkground/
     1        ((1.-source_smo(n).bkground)*real(source(n).ngrid))
      do i=1,source(n).ngrid
          source(n).grids_w(i)=(source(n).grids_w(i)+added_w)/
     1    (1.+source_smo(n).bkground/(1.-source_smo(n).bkground))
      enddo
      return
      end



      subroutine find_n_ring(n,n_ring)
      use                   smooth_seis
      integer                                    ::n, n_ring
      real                                       ::c_ring
      c_ring=3.0*source_smo(n).corr_dis/source_smo(n).ring_w
      n_ring=int(c_ring)
      if((c_ring-real(n_ring))>0.01)n_ring=n_ring+1
      return
      end



      subroutine Weichert_fac(n,Wei_fac)
      use                   smooth_seis
      integer                                    ::i, j, n
      real                                       ::Wei_fac, numer, denom
      numer=0.
      denom=0.
      Wei_fac=0.
      do i=1,source_smo(n).n_comp
          numer=numer+exp(-source_smo(n).beta*source_smo(n).comp_m(i))
          denom=denom+(source_smo(n).end_year-source_smo(n).comp_yr(i))
     1               *exp(-source_smo(n).beta*source_smo(n).comp_m(i))
      enddo
      Wei_fac=numer/denom
      return
      end      