cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   EqHaz02 Ver. b (2011/11/22)
ccc
ccc	The programs Haz02 is developed to read catalog files generated by the p
ccc	rogram 'EqHaz1' (at most three EqHaz1 outputs) and calculate ground shak
ccc	ing level (PGA, PGV, and PSA) at a number of sites (at most 1000) for at
ccc   given ground motion prediction tables (at most 5). In this version of pr
ccc   ogram earthquake catalog files (EqHaz1 outputs) with records numbers lar
ccc   ger than 1000000 are not supported.
ccc
ccc		* This program is developed using Intel Visual Fortran and the compati
ccc		  bility of this code with other compilers is not tested.
ccc		* This program does the calculation of ground motion measures by inter 
ccc		  polating the expected values from provided tables.
ccc     * The program asks for a control file. The control file should carry t
ccc       he information of catalog file(s), ground motion table file(s) and t
ccc       he site(s). Also, some constraints on distance and magnitude ranges
ccc       should be included and the program will use ONLY the records of EqHa
ccc       z1 catalog file(s) that fall in the range.
ccc     * The outputs of the program are: a) ground motion catalogs for every 
ccc       combination of sites, earthquake catalogs (from EqHaz1), and ground
ccc       motion tables; b) ground motion cataloge for every combination site
ccc       and earthquake catalog, derived from combining catalogs of same sit
ccc       e/catalog-file set with different ground motions; c) hazard curce fi
ccc       les for every combined catalog file ("b" types in this paragraph).
ccc     * We used the following convention for naming these three file types (
ccc       ? stands for a number in the following):
ccc         - "a" type files: stem_S?C?G?.tab    where S? stands for site #, 
ccc           C? for catalog #, and G? for ground motion table#.
ccc         - "b" type files: stem_S?C?.out      where S? and C? are as above.
ccc         - "c" type files: stem_S?C?.haz      where S? and C? are as above.
ccc     * For each ground motion table file the following information should b
ccc       included in control file: a) file name; b)standard devition; c) mult
ccc       iplicative factor; and d) weighting factor.
ccc       Standard deviation is the value corresponding to ground motion model
ccc       at implemented frequency. Multiplicative factor is actually a number
ccc       which is added to the logarithm of calculated ground motion, and the
ccc       weight is a factor that assignd the importance of each specific grou
ccc       nd moion model and the sum of weights of all GM models should be one.
ccc       these weights are used for generating type "b" catalogs from type "a".
ccc     * For each earthquake catalog file, the name and weight of catalog file
ccc       should be specified in control file. Although the program checks for
ccc       the sum of weights of catalog files, but in this version of EqHaz2 we
ccc       didn't assign the program to join them and this feature is left for f
ccc       uture developments. The current advantage of multiple catalogs in thi
ccc       s version is that the user can generate all necessary GM catalogs in
ccc       one run (for example H and R models).
ccc     * For each site the user should include the coordinates of the site and
ccc       the frequency of required PSA in control file (99 and 89 are used for
ccc       PGA and PGV repectively) .
ccc     * The user needs to add the name of output file stem name and ranges of
ccc       magnitudes and distance for acceptable earthquake catalog records.
ccc     * All the files including EqHaz2b application, control file, earthquake
ccc       catalog file(s), ground motion table file(s) in one directory. Also i
ccc       t is necessary to have sufficient space for storing output files.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   EqHaz02 Ver. b (2011/12/15)
ccc
ccc   1-The program accepts magnitude type declaration in input and in case th
ccc     e type is MN, program converts it to Mw.
ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   EqHaz02 Ver. b (2012/02/20)
ccc
ccc   1-The program accepts different flags Rcd and Rjb calculation. Flag valu
ccc     es are: 1: apply Goda et al., 2010 formulae; 2: apply Atkinson 2011 fo
ccc     r west; 3: apply Atkinson 2011 for east; 4: apply Atkinson 2011 for ea
ccc     st if GMPE is AB06.
ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   EqHaz02 Ver. 2012c (2012/06/27)
ccc
ccc   1-The program works with EqHaz1_2012c output. The difference is only add
ccc     ition of two lines in header of EqHaz1_2012c output: 1-the name of inp
ccc     ut file to EqHaz1_2012c is added; 2-a flag is added showing if epistem
ccc     ic and aleatory uncertainties are mixed in EqHaz1_2012c output.
ccc   2-The program allows for mixing (or not mixing) epistemix and aleatory u
ccc     ncertainties by allowing different GMPEs apply to subcatalogs or fix i
ccc     t (note that the concept of subcatalogs is necessary when the user dec
ccc     ides not to mix uncertainties, although use of subcatalogs is always n
ccc     ecessary from practical point of view). This capability is controlled
ccc     by user in EqHaz1_2012c and it is set there. EqHaz2_2012c learns about 
ccc     this parameter from output of EqHaz1_2012c.
ccc   3-The program allows for mixing (or not mixing) uncertainties according
ccc     to what is done in EqHaz1_2012c.
ccc   4-Fault geometries are added to this version so the exact calculations o
ccc     f Rcd and Rjb for faults are possible (instead of appoximating them us
ccc     ing some empirical relationships). For calculating fracture area on fa
ccc     ult surface, Wells and Coppersmith equations are used. Although use of
ccc     this option makes the caculations more exact, but they make the progra
ccc     m a little slow and in the case of calculations for large numver of si
ccc     tes (map preparation) it may mean considerable time difference. When t
ccc     he program is assigned to do exact fault distance calculations, it fin
ccc     ds fault geometry information from EqHaz1_2012c input file, and in cas
ccc     e the file is not existing in working directory, the program will assu
ccc     me all sources are areal. The only thing the user has to do is to set 
ccc     the fault distance calculation flag to 1 in EqHaz2_2012c input and tha
ccc     t is all, the name of EqHaz1 input is read from its output.
ccc   5-An option is added for map preparation. when the program is in this mo
ccc     de, it does not generate lots of outputs. In the input file the user h
ccc     as to list the site coordinates for which mean hazard motions for requ
ccc     ired return period will be tabulated. Most of outputs and hazard curve
ccc     calculations are disables to make the program generate quick results.
ccc   6-Flag is added that allows the user to suppress generation of ground mo
ccc     tion tables for each individual GMPE (*.tab files). This helps the pro
ccc     gram to run much faster and not waste memory if these intermediate fil
ccc     es are unnecessary. To make the map preparation more practical the max
ccc     imum number of sites accepted by program is increased to 100000.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      module datastruct      
          structure /catalog/
              character (len=50)                :: catname, EqHaz1_inp          !Catalog file name. Corresponding EqHaz1 input file name
              character (len=80), dimension(50) :: cathead                      !Catalog file header lines
              real                              :: wgt                          !Weight of this catalog in hazard. The sum of weights of all catalogss should be 1.
              integer                           :: nline_hdr                    !Number of header lines in catalog file
              integer                           :: n_of_cats, n_of_yrs          !Number of sub-catalogs and duration of each sub-catalog
              integer                           :: n_rec                        !Number of records in a catalog file
              integer, dimension(1000000)       :: cat_no                       !Subcatalog number of a record
              real, dimension(1000000)          :: time                         !Time of a record in catalog
              real, dimension(1000000,3)        :: loc                          !Latitude, longitude, and depth of a record in catalog
              real, dimension(1000000)          :: mag                          !Magnitude of a record in catalog
              integer, dimension(1000000)       :: zone_no                      !Zone number of a record in catalog
              real, dimension(1000000,2)        :: N0beta                       !N0 and beta of a record in catalog
              real, dimension(1000000)          :: Mmax                         !Mmax of a record in catalog
              integer                           :: epist                        !Epistemic uncertainty mixing flag (use 1 for mixing and 0 for not mixing [traditional]). This change is for creating constant parameters for a zone in a subcatalog-April-2012
              integer                           :: allread                      !Flag specifying if catalog is read completely in arrays
          end structure
          
          structure /gmpe/
              character (len=50)                :: gmpname                      !GMPE file name
              character (len=80), dimension(50) :: gmphead                      !GMPE file header lines
              integer                           :: nline_hdr                    !Number of header lines in GMPE file
              real                              :: wgt                          !Weight of this GMP model in hazard. The sum of weights of all GMP models should be 1.
              real                              :: sigma, const                 !Standard deviation and constant value added to GM (note: in current version sigma is not used and standard deviations are introduced in GM table file)
              real                              :: sigma_model                  !Standard deviation considered for GM. These values are read from GM table
              integer                           :: n_mag, n_dis, n_freq         !Number of magnitudes, distances, and frequencies given in GMPE table
              integer                           :: disflag                      !Distance type flag (0:Repi, 1:Rhyp, 2:Rcd, and 3:Rjb)
              real, dimension(50)               :: mag, dis, freq, GMsig        !Arrays for keeping the values of magnitudes, distances, and frequencies given in tables
              real, dimension(50,50,50)         :: table                        !Array for keeping the GM values of tables
          end structure

          structure /site/
              real, dimension(2)                :: coord                        !Latitude and longitude of site
              real                              :: freq                         !Frequency for which GM is calculated
              real, dimension(2,2)              :: rng                          !Range of magnitudes and distances accepted for cataloging
          end structure

          structure /range_info/                                                !Realizations of components of this structure are mainly used for interpolation calculations
              real, dimension(1000000,4)        :: dist                         !Distances between event and site. col1: Repi, col2: Rhyp, col3: Rcd, and col4: Rjb
              integer, dimension(1000000,3)     :: flags                        !Three flags for each new calculation. For frequency, magnitude, and distance respectively
              integer                           :: nline, nline_inrange         !Number of lines and number of in-magnitude and in-distance range lines
          end structure
          
          structure /hazard/
              real, dimension(1000000,2)        :: gmpar                        !Calculated GM parameter and assigned multiple of standard deviation
              real, dimension(1000000,2)        :: res                          !reserved columns for future use
              integer                           :: nline                        !Number of lines 
          end structure
          
          record /catalog/ cat(6)                                               !This line creates an array of 3 catalog structures
          record /gmpe/ gmp(5)                                                  !This line creates an array of 10 GMPE model structures
          record /site/ sit(100000)                                             !This line creates an array of 100000 site structures
          record /range_info/ M_R_limit(6,5)                                    !This line creates an array of 3 range srtuctures for tracking in-range catalog records
          record /hazard/ haz(6,5)                                              !This line creates an array of 30 hazard scenario structures
          record /hazard/ hazall                                                !This line creates one hazard structure for final solution
          
          integer                                 ::cat_mag_typ
          integer                                 ::i_dist_conv
          integer                                 ::i_flt_dist
          integer                                 ::i_map
          integer                                 ::ngmpe, ncat, nsite
          integer                                 ::i_maketab
          integer                                 ::f_serial=0
          integer                                 ::epist                       !Epistemic uncertainty mixing flag (use 1 for mixing and 0 for not mixing [traditional]). This change is for creating constant parameters for a zone in a subcatalog-April-2012
c          integer                                 ::good_num
          integer, dimension(1000000,3)           ::m_ind, r_ind, f_ind         !Array of index limits for magnitude, distance and frequency
          real                                    ::period_map
          real, dimension(1000000)                ::good_gm                     !This is an array of ground motions (only) reported after all averagings
          real, dimension(1000000)                ::m_frac,r_frac,f_frac        !Fractions of magnitude, distance, and frequency to be added to lower limit values
          real, dimension(100,4)                  ::haz_curve                   !The calculated hazard curve is stored in this array
          character (len=50)                      ::input,out_stem
          character (len=80), dimension(1000,6,5) ::output
          character (len=80), dimension(1000,6)   ::output_w
          
      end module


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   This module is copied from EqHaz1 as its components will be used for rea
ccc   ding and storing EqHaz1 input. Also it has a couple of added lines from
ccc   the program fault_dist for carrying out asperity specification calculati
ccc   ons.
ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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(20000,2)                             !Array of nodes (reserved 50 coordinates) and array of grid points (reserved 20000 coordinates). The geographical coordinates in degrees should be used (Order: Lat. long.).
            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)
            integer      isasperity(50)                                         !An array specifying what segments have portions of asperities
            real         asperity(50,4,3)                                       !Coordinates of corners of asperity on a segment
        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                                            !Seed value for randon number generator (0: random, none-zero values will enter to the program). Epistemic uncertainty mixing flag (use 1 for mixing and 0 for not mixing [traditional])
        integer          icnt
!        integer          date_time(8)

        record	/source/source
        record  /source/tmp_src,src(60)
        
      end module srcdef

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc   This module is copied from the program fault_dist for incorporating para
ccc   meters for fault distance calculations.
ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      module fltdis
          real                                  :: Lf,Wf,l,w,x,y
          real                                  :: Ax,Ay,Bx,By
          real                                  :: eqlat,eqlon,eqdep
          real                                  :: eqmag
          real                                  :: Repi,Rhyp,Rjb,Rcd
          real                                  :: Rjb_all,Rcd_all
          real                                  :: plat,plon
          integer                               :: ctrl_nlines,ix
          integer                               :: input_nlines
          character (len=50)                    :: ctrl, srcfile
          character (len=50)                    :: output_file
          character (len=80), dimension (100)   :: ctrl_lines
          character (len=80), dimension (10)    :: input_lines
      end module fltdis

                          
      program EqHaz2
      use datastruct
      use srcdef
      use fltdis
      use dflib
      integer(2)				n_of_args

      icnt=0
      n_of_args=nargs()
      if(n_of_args.gt.1)then
          call getarg(1,input)
      else
          write(*,*) 'Enter input file name : '
          read(*,*) input
      endif
      
      write(*,'(a)')'Reading input and checking weights:'
      call read_input()
      call check_weights()
      if(i_map==1)then                                                          !Added/modified for allowing map generation
          open(unit=30,file=trim(out_stem)//'.map',status='unknown')            !Added/modified for allowing map generation
          write(30,'(f8.1,a)')period_map,' Year return period mean '//          !Added/modified for allowing map generation
     1                        'hazard at a grid of points:'                     !Added/modified for allowing map generation
          write(30,'(a)')'   Latitude      Longitude    '//                     !Added/modified for allowing map generation
     1                   '   Ground motion'                                     !Added/modified for allowing map generation
      endif                                                                     !Added/modified for allowing map generation
      do k=1,nsite
          do j=1,ncat
              if(i_flt_dist.eq.1)then                                           !Added for reading the contents of original EqHaz1 input file name
                  call readEQHAz1input(cat(j).EqHaz1_inp)                       !Added for reading the contents of original EqHaz1 input file name
              else                                                              !Added for reading the contents of original EqHaz1 input file name
                  src.typ=1                                                     !Added for reading the contents of original EqHaz1 input file name
              endif                                                             !Added for reading the contents of original EqHaz1 input file name
              do i=1,ngmpe
                  write(*,1)'Procesing site #',k,'Catalog #',j,'GM #',i
1                 format(a,i5,2x,a,i2,2x,a,i3)
                  f_serial=f_serial+1
                  call make_output_name(k,j,i)
                  call fill_dist(k,j,i)
                  call find_limits(k,j,i)
                  call calc_hazard(k,j,i)
                  if(i_maketab.eq.1.and.i_map.eq.0)then                         !Added/modified for allowing map generation
                      call make_output(k,j,i,'tab')                             !Added/modified for allowing map generation
                  endif                                                         !Added/modified for allowing map generation
              enddo
              write(*,1)'Calculating hazard for site #',k,'Catalog #',j
              f_serial=f_serial+1
              call make_output_w_name(k,j)
              call calc_hazall(k,j)
              if(i_map.eq.0)then                                                !Added/modified for allowing map generation
                  call make_output_w(k,j,'out')                                 !Added/modified for allowing map generation
                  call make_hazard_output_w(k,j,'haz')                          !Added/modified for allowing map generation
              else                                                              !Added/modified for allowing map generation
                  call make_output_w_map(k,j,'out')                             !Added/modified for allowing map generation
                  call make_hazard_map(k,j)                                     !Added/modified for allowing map generation
              endif                                                             !Added/modified for allowing map generation
          enddo
      enddo
      if(i_map.eq.1)close(unit=30)
      end

 
      
c
c	This subroutine makes an output file name based on output stem name,
c     site number, GM model number, and catalog number. Also, it assigns a
c     serial number to the name of file.
c
      subroutine make_output_name(n_s,n_c,n_g)
      use datastruct
      integer                                   :: n_s, n_c, n_g
      character (len=4)                         :: c4
      character (len=5)                         :: c5, c5_
      character (len=6)                         :: c6
      write(c4,'(a1,i2.2,a1)')'G',n_g,' '
      if(n_g <10)write(c4,'(a1,i1,a2)')'G',n_g,'  '
      write(c5,'(a1,i1,a3)')'C',n_c,'   '
      write(c6,'(a1,i4,a1)')'S',n_s,' '
      if(n_s <1000)write(c6,'(a1,i3,a2)')'S',n_s,'  '
      if(n_s <100)write(c6,'(a1,i2,a3)')'S',n_s,'   '
      if(n_s <10)write(c6,'(a1,i1,a4)')'S',n_s,'    '
      write(c5_,'(a1,i4.4)')'@',mod(f_serial,10000)
      
      output(n_s, n_c, n_g)=trim(out_stem)//'_'//trim(c6)//trim(c5)//
     1                      trim(c4)
     
      return
      end



      subroutine make_output_w_name(n_s,n_c)
      use datastruct
      integer                                   :: n_s, n_c
      character (len=5)                         :: c5, c5_
      character (len=6)                         :: c6

      write(c5,'(a1,i1,a1)')'C',n_c,' '
      write(c6,'(a1,i4,a1)')'S',n_s,' '
      if(n_s<1000)write(c6,'(a1,i3,a2)')'S',n_s,'  '
      if(n_s<100)write(c6,'(a1,i2,a3)')'S',n_s,'   '
      if(n_s<10)write(c6,'(a1,i1,a4)')'S',n_s,'    '
      write(c5_,'(a1,i4.4)')'@',mod(f_serial,10000)
      
      output_w(n_s, n_c)=trim(out_stem)//'_'//trim(c6)//trim(c5)
      
      return
      end



c
c     This subroutine checks the sum of the weights of GMPEs and catalogs t
c     o make sure they are properly weighted. In the case of an error by th
c     e user a report will be generated on screen and the program will stop
c     running.
c
      subroutine check_weights()
      use datastruct
      integer                                   :: flag=0
      real                                      :: sum
      sum=0.
      do i=1,ngmpe
          sum=sum+gmp(i).wgt
      enddo
      if (abs(sum-1).gt.0.001)then
          flag=1
          write(*,*)'Sum of weights of GMPEs is not equal to 1!'
      endif
      sum=0.
      do i=1,ncat
          sum=sum+cat(i).wgt
      enddo
      if (abs(sum-1).gt.0.001)then
          flag=1
          write(*,*)'Sum of weights of catalogs is not equal to 1!'
      endif
      if(flag.eq.1)then
          write(*,*)'Program aborted. Please check the input file.'
          call exit()
      endif
      return
      end



c
c     This subroutine reads the input file and it calls other reading subro
c     utines for reading the GMPE, catalog file. It only reads the informat
c     ion of sites as they are introduced in input file.
c
      subroutine read_input()
      use datastruct

      integer(2)                                :: delfiles 
      real, dimension(2,2)                      :: tmp
      character (len=50)                        :: tempinput
      tempinput='tempinput'
      call reformat(input,tempinput)

      open(unit=1,file=tempinput,status='old')
      read(1,*) cat_mag_typ, i_dist_conv, i_flt_dist
      read(1,*)i_map, period_map
      read(1,*) ngmpe, i_maketab
      do i=1,ngmpe
          read(1,*) gmp(i).gmpname
c          read(1,*) gmp(i).sigma, gmp(i).const
          read(1,*) gmp(i).const
          read(1,*) gmp(i).wgt
          call read_GMPE(i)
      enddo
      read(1,*) ncat
      do i=1,ncat
          read(1,*) cat(i).catname
          read(1,*) cat(i).wgt
          call read_catalog(i)
      enddo
      read(1,*)nsite
      read(1,*)tmp(:,:)
      do i=1,nsite
          read(1,*)sit(i).coord(:),sit(i).freq
          sit(i).rng=tmp
      enddo
      read(1,*)out_stem
      close(unit=1)
      
      delfiles=delfilesqq('tempinput')

      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=input,status='old')
      open(unit=12,file=tempinput,status='unknown')
      do while (not(eof(11)))
          flag=1
          read(11,'(a100)')line
          do i=1,10
              if(line(1:1).eq.c1(i))flag=0
          enddo
          if(flag.eq.1)then
              write(12,'(a)')line(1:80)
          endif
      enddo
      close(unit=11)
      close(unit=12)
      return
      end



c
c	This subrouotine reads the GMPE tables and stores the values in arrays
c	table(50,50,50), mag(50), dis(50), freq(50), and GMsig(50). The array 
c	'table' carries the values of PSAs, PGA, and PGV values. The array mag
c	, dis, and freq carry list of magnitudes, distances, and frequencies.  
c	These arrays reserves a maximum of 30 values for these three parameter
c	s. The numbers of these parameters are stores in the variables: n_mag, 
c	n_dis, and n_freq. These values are returned to the main and other sub
c     -programs through the implementation datastruct module.
c	
      subroutine read_GMPE(n)
      use datastruct
      real                                      :: tmp
      character (len=80)                        :: c80
      open(unit=2,file=gmp(n).gmpname,status='old')
      mag= 0.
      dis= 0.
      freq= 0.
      GMsig=0.
      gmp(n).nline_hdr=1
      do i=1,gmp(n).nline_hdr
        read(2,'(a80)')gmp(n).gmphead(i)
      enddo
      read(2,*)gmp(n).n_mag,gmp(n).n_dis,gmp(n).n_freq,gmp(n).disflag           !# of magnitudes, # of distances, and # of frequencies
      read(2,*)(gmp(n).freq(k),k=1,gmp(n).n_freq)
      read(2,*)(gmp(n).GMsig(i),i=1,gmp(n).n_freq)                              !Standard deviations of GM in base 10 logarithmic units
      do i=1,gmp(n).n_mag
          do j=1,gmp(n).n_dis
            read(2,*)gmp(n).mag(i),gmp(n).dis(j),
     1              (gmp(n).table(k,j,i),k=1,gmp(n).n_freq)
            gmp(n).dis(j)=log10(gmp(n).dis(j))
          enddo
      enddo
      close(unit=2)
      
      return
      end



c
c	This subroutine reads the catalog file generated by program 'haz01'.
c     The information are stored in an element of a structure through the
c     implementation of datastruct module.
c
      subroutine read_catalog(n)
      use datastruct
      integer                                   :: nl=0

!      cat(n).nline_hdr=6                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      cat(n).nline_hdr=8                                                        !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      open(unit=3,file=cat(n).catname,status='old')
      do i=1,cat(n).nline_hdr
          read(3,'(a80)')cat(n).cathead(i)
      enddo
      read(cat(n).cathead(3)(35:),*)cat(n).EqHaz1_inp                           !Added for reading the name of original EqHaz1 input file name
      read(cat(n).cathead(7)(56:),*)cat(n).epist                                !Added for reading the applied epistemic flag value
      do while ( not( eof(3)))
          nl=nl+1
          read(3,*)cat(n).cat_no(nl),cat(n).time(nl),
     1             cat(n).loc(nl,:),cat(n).mag(nl),cat(n).zone_no(nl),
     2             cat(n).N0beta(nl,:),cat(n).mmax(nl)
      enddo
      close(unit=3)
      cat(n).n_rec=nl
      read(cat(n).cathead(5)(28:),*)cat(n).n_of_cats
      read(cat(n).cathead(6)(28:),*)cat(n).n_of_yrs
      if(cat_mag_typ /= 0)call correct_magnitude(n)
      return
      end



      subroutine correct_magnitude(n)
      use datastruct
      real                                      :: tmp
      
      if(cat_mag_typ ==1)then
          do i=1,cat(n).n_rec
              if(cat(n).mag(i) <= 5.5)then
                  tmp=0.98*cat(n).mag(i)-0.39
                  cat(n).mag(i)=tmp
              else
                  tmp=2.715-0.277*cat(n).mag(i)+0.127*cat(n).mag(i)**2
                  cat(n).mag(i)=tmp
              endif
          enddo
      endif
      return
      end



c
c	This subroutine allows to skip 'lines' records in the unit # 'iunit'.
c
      subroutine skip_lines(iunit,lines)
      integer                                   :: iunit, lines
      if(lines.ge.1)then
          do i=1,lines
              read(iunit,*)
          enddo
      endif
      return
      end
      
      
      
c
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(arr,n)
      dimension arr(n,1)
      pi=4*atan(1.)
          do i=1,n
              u1=0.
              u2=0.
              do while (u1.eq.0.or.u2.eq.0.)
                  call rand_num_gen(u1,0)
                  call rand_num_gen(u2,0)
              enddo
              arr(i,1)=sqrt(-2*log(u1))*cos(2*pi*u2)
          enddo
      return
      end



      function boxmul()
      real                                      :: boxmul, pi, u1, u2
      pi=4*atan(1.)
      u1=0.
      u2=0.
      do while (u1.eq.0.or.u2.eq.0.)
          call rand_num_gen(u1,0)
          call rand_num_gen(u2,0)
      enddo
      boxmul=sqrt(-2*log(u1))*cos(2*pi*u2)
      return
      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 rand_num_gen(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	This subroutine makes a list of epicentral, hypocentral, closest dist
c	ances for a list of earthquakes in catalog 'cat' from the site 'sit'.
c	This list of distances is for all events catalog 'cat' from one site 
c     'sit'. The index of catalog is 'n_c' and the index of site is 'n_s'. 
c     There are a total of 'cat(n_c).n_rec' records.
c     The calculated distances are all in km.
c
      subroutine fill_dist(n_s,n_c,n_g)
      use datastruct
      use srcdef
      use fltdis
      integer                                   :: n_s, n_c, n_g
      integer                                   :: i
      M_R_limit(n_c,n_g).dist=0.      
      M_R_limit(n_c,n_g).flags=0
      M_R_limit(n_c,n_g).nline=cat(n_c).n_rec
      do i=1,M_R_limit(n_c,n_g).nline
          if(src(cat(n_c).zone_no(i)).typ == 2.and.i_flt_dist == 1)then
              source=src(cat(n_c).zone_no(i))
              eqlat=cat(n_c).loc(i,1)
              eqlon=cat(n_c).loc(i,2)
              eqdep=cat(n_c).loc(i,3)
              eqmag=cat(n_c).mag(i)
              plat=sit(n_s).coord(1)
              plon=sit(n_s).coord(2)
              call find_distances()
              M_R_limit(n_c,n_g).dist(i,1)=Repi
              M_R_limit(n_c,n_g).dist(i,2)=Rhyp
              M_R_limit(n_c,n_g).dist(i,3)=Rcd
              M_R_limit(n_c,n_g).dist(i,4)=Rjb
          else
              M_R_limit(n_c,n_g).dist(i,1)=dis(sit(n_s).coord(1),
     1                                     sit(n_s).coord(2),
     2                                     cat(n_c).loc(i,1),
     3                                     cat(n_c).loc(i,2))
              M_R_limit(n_c,n_g).dist(i,2)=sqrt(M_R_limit(n_c,n_g).
     1                           dist(i,1)**2+cat(n_c).loc(i,3)**2)
              M_R_limit(n_c,n_g).dist(i,3)=Rhyp2Rrup(cat(n_c).mag(i),
     1                                     M_R_limit(n_c,n_g).dist(i,2),
     2                                     cat(n_c).loc(i,3),
     3                                     i_dist_conv)
              M_R_limit(n_c,n_g).dist(i,4)=Repi2Rjb(cat(n_c).mag(i),
     1                                     M_R_limit(n_c,n_g).dist(i,1),
     2                                     i_dist_conv)
          endif
      enddo
      return
      end



c
c	This subroutine estimated the rupture distance of a fault from a site
c	in kilometers given the hypocentral distance and magnitude of the ear
c	thquake. This approximation is based on study by Goda et al., 2010, a
c     nd Atkinson 2011. 
c	
      function Rhyp2Rrup(M,R,h,i_dist_conv)
      real                                      :: M, R, h
      real, dimension(7)                        :: a
      real                                      :: flt_l, flt_w
      real                                      :: Dtop, Ztor, Dsurf
      integer                                   :: i_dist_conv
      data a/0.363,-0.0437,0.978,-0.0256,-1.213,-0.149,0.104/
      if(i_dist_conv.eq.1)then
          dR=(1-exp(-(a(1)+a(2)*M)*R**(a(3)+a(4)*M)))
     1         *exp(a(5)+a(6)*M+a(7)*M**2)
          Rhyp2Rrup=R-dR
      elseif(i_dist_conv.eq.2)then
          flt_l=10.**(-2.44+0.59*M)
          flt_w=10.**(-1.01+0.32*M)
          Dtop=h-0.5*flt_w
          Ztor=amax1(Dtop,0.)
          Dsurf=amax1((sqrt(R**2.-h**2.)-0.3*flt_l),0.1)
          Rhyp2Rrup=sqrt(Dsurf*Dsurf+Ztor*Ztor)
      elseif(i_dist_conv.eq.3)then
          flt_l=0.6*10.**(-2.44+0.59*M)
          flt_w=0.6*10.**(-1.01+0.32*M)
          Dtop=amax1((21.-2.5*M),(h-0.5*flt_w))
          Ztor=amax1(Dtop,0.)
          Dsurf=amax1((sqrt(R**2.-h**2.)-0.3*flt_l),0.1)
          Rhyp2Rrup=sqrt(Dsurf*Dsurf+Ztor*Ztor)
      else
          write(*,*)'Incorrect distance conversion flag-program stopped'
          call exit()
      endif
      return
      end



!c
!c	This subroutine estimated the rupture distance of a fault from a site
!c	in kilometers given the hypocentral distance and magnitude of the ear
!c	thquake. This approximation is based on study by Goda et al., 2010, a
!c     nd Atkinson 2011. 
!c	
!      function Rhyp2Rrup(M,R,h,i_dist_conv)
!      real                                      :: M, R, h
!      real, dimension(7)                        :: a
!      real                                      :: flt_l, flt_w
!      real                                      :: Dtop, Ztor, Dsurf
!      integer                                   :: i_dist_conv
!      data a/0.363,-0.0437,0.978,-0.0256,-1.213,-0.149,0.104/
!      if(i_dist_conv.eq.1)then
!          dR=(1-exp(-(a(1)+a(2)*M)*R**(a(3)+a(4)*M)))
!     1         *exp(a(5)+a(6)*M+a(7)*M**2)
!          Rhyp2Rrup=R-dR
!      elseif(i_dist_conv.eq.2)then
!          flt_l=10.**(-2.44+0.59*M)
!          flt_w=10.**(-1.01+0.32*M)
!          Dtop=h-0.5*flt_w
!          Ztor=amax1(Dtop,0.)
!          Dsurf=amax1((sqrt(R**2.-h**2.)-0.3*flt_l),0.1)
!          Rhyp2Rrup=sqrt(Dsurf*Dsurf+Ztor*Ztor)
!      elseif(i_dist_conv.eq.3)then
!          flt_l=0.6*10.**(-2.44+0.59*M)
!          flt_w=0.6*10.**(-1.01+0.32*M)
!          Dtop=h-0.5*flt_w
!          Ztor=amax1(Dtop,0.)
!          Dsurf=amax1((sqrt(R**2.-h**2.)-0.3*flt_l),0.1)
!          Rhyp2Rrup=sqrt(Dsurf*Dsurf+Ztor*Ztor)
!      elseif(i_dist_conv.eq.4)then
!          flt_l=0.6*10.**(-2.44+0.59*M)
!          flt_w=0.6*10.**(-1.01+0.32*M)
!          Dtop=amax1((21.-2.5*M),(h-0.5*flt_w))
!          Ztor=amax1(Dtop,0.)
!          Dsurf=amax1((sqrt(R**2.-h**2.)-0.3*flt_l),0.1)
!          Rhyp2Rrup=sqrt(Dsurf*Dsurf+Ztor*Ztor)
!      else
!          write(*,*)'Incorrect distance conversion flag-program stopped'
!          call exit()
!      endif
!      return
!      end



c
c	This subroutine estimated the Joyner-Boore distance of a fault from a 
c	site in kilometers given the epicentral distance and magnitude of the 
c	earthquake. This approximation is based on study by Goda et al., 2010
c     , and Atkinson 2011. 
c	
      function Repi2Rjb(M,R,i_dist_conv)
      real                                      :: M, R
      real, dimension(7)                        :: a
      real                                      :: flt_l
      integer                                   :: i_dist_conv
      data a/0.458,-0.0549,1.046,-0.0361,-1.297,-0.138,0.105/
      
      if(i_dist_conv.eq.1)then
          dR=(1-exp(-(a(1)+a(2)*M)*R**(a(3)+a(4)*M)))
     1         *exp(a(5)+a(6)*M+a(7)*M**2)
          Repi2Rjb=R-dR
      elseif(i_dist_conv.eq.2)then
          flt_l=10.**(-2.44+0.59*M)
          Repi2Rjb=amax1((R-0.3*flt_l),0.1)
      elseif(i_dist_conv.eq.3)then
!      elseif(i_dist_conv.eq.3.or.i_dist_conv.eq.4)then
          flt_l=0.6*10.**(-2.44+0.59*M)
          Repi2Rjb=amax1((R-0.3*flt_l),0.1)
      else
          write(*,*)'Incorrect distance conversion flag-program stopped'
          call exit()
      endif
      return
      end



c
c	This function calculates the distance between two points on earth in
c	kilometers given the coordinates of those points in degrees. It is as
c	sumed that earth is spherical with 6378.137 km radius.
c
      function dis(lat1,lon1,lat2,lon2)
      real                                      :: lat1,lon1,lat2,lon2
      real                                      :: dis
      r=6378.137
      y1=deg2rad(lat1)
      x1=deg2rad(lon1)
      y2=deg2rad(lat2)
      x2=deg2rad(lon2)
      dis=r*acos(cos(y1)*cos(y2)*cos(x1-x2)+sin(y1)*sin(y2))
      return
      end
      


c
c	This subroutine converts the scaler 'x' from degrees to radians.
c
      function deg2rad(x)
      pi=4.*atan(1.)
      deg2rad=x*pi/180.
      return
      end



c
c	This subroutine converts the scaler 'x' from radians to degrees.
c
      function rad2deg(x)
      pi=4.*atan(1.)
      rad2deg=x*180./pi
      return
      end
      
      

c
c     This subroutine links between previously made code "get_limits" and t
c     he EqHaz2b version with its structured data. This subroutine makes th
c     e flags, finds indices of each magnitude, distance, and frequency val
c     ues from GM table for being used for interpolation.
c
      subroutine find_limits(n_s, n_c, n_g)
      use datastruct
      integer                                   :: n_s, n_c,n_g
      integer                                   :: i
      M_R_limit(n_c,n_g).nline=cat(n_c).n_rec
      do i=1,M_R_limit(n_c,n_g).nline
          call get_limits(cat(n_c).mag(i),
     1                    M_R_limit(n_c,n_g).dist(i,gmp(n_g).disflag+1),
     2                    sit(n_s).freq,
     3                    gmp(n_g).mag,gmp(n_g).dis,gmp(n_g).freq,
     4                    gmp(n_g).n_mag,gmp(n_g).n_dis,gmp(n_g).n_freq,
     5                    m_frac(i),r_frac(i),f_frac(i),
     6                    m_ind(i,1),m_ind(i,2),
     7                    r_ind(i,1),r_ind(i,2),
     8                    f_ind(i,1),f_ind(i,2),
     9                    M_R_limit(n_c,n_g).flags(i,1:3))
      enddo
      return
      end
      
      
      
c
c     This subroutine gets the value of questioned magnitude, distance, and
c     frequency along with the available values of same quantities in GM ta
c     ble, and returns indices of the bounding values in m-r-f space. Also
c     it returns the values of fractional changes on each vector of m-r-f s
c     pace. These fractional values can be used for interpolation.
c     This subroutine is kept this form as it was working well in EqHaz2a v
c     ersion and its results were verified.
c
      subroutine get_limits(amag,r,freq,xmag,rlog,f,nm,nd,nf,fracm,fracd
     1                     ,fracf,jml,jmu,jdl,jdu,jfl,jfu,iflag)
      real, dimension(30)             :: xmag, rlog, f                          !xmag, rlog, and f arrays taken from table
      integer                         :: nm, nd, nf
      real                            :: fracm, fracd, fracf
      real                            :: sigma_model
      integer                         :: jml, jmu, jdl, jdu, jfl, jfu 
      integer, dimension(1,3)         :: iflag
      
      iflag=4
      
c     Find bounding jf values.
      jfl = 1
      jfu = 1
      fracf=0.
      do j=1,nf-2
          if(abs(freq-f(j)) .lt. 0.01) then
              jfl = j
              jfu = j
              iflag(1,1)=0
              exit
          endif
      enddo
      if(abs(freq-f(nf-1)) .lt. 0.01) then
          jfl = nf-1
          jfu = nf-1
          iflag(1,1)=2
      endif
      if(abs(freq-f(nf)) .lt. 0.01) then
          jfl = nf
          jfu = nf
          iflag(1,1)=3
      endif
      if(iflag(1,1).eq.4.and.freq.lt.f(nf-2).and.freq.ge.f(1))then
          do j = 1,nf-3
              if (freq .ge. f(j) .and. freq .lt. f(j+1))then
                  jfl = j
                  jfu = j + 1
                  fracf = (alog10 (freq) - alog10 (f(jfl)) ) /
     1                    (alog10 (f(jfl+1)) - alog10 (f(jfl)))
              endif
          enddo
          iflag(1,1)=1
      elseif(iflag(1,1).eq.4.and.freq.gt.f(nf-2))then
          iflag(1,1)=5
      endif

c     Find bounding jm values.
      jml = 1
      jmu = 1
      fracm=0.
      do j=1,nm
          if(abs(amag-xmag(j)) .lt. 0.01) then
              jml = j
              jmu = j
              iflag(1,2)=0
              exit
          endif
      enddo
      if(iflag(1,2).eq.4.and.amag.lt.xmag(nm).and.amag.ge.xmag(1))then
          do j = 1,nm-1
              if (amag .ge. xmag(j) .and. amag .lt. xmag(j+1))then
                  jml = j
                  jmu = j + 1
                  fracm = (amag - xmag(jml)) / (xmag(jml+1)-xmag(jml))
              endif
          enddo
          iflag(1,2)=1
      elseif(iflag(1,2).eq.4.and.amag.gt.xmag(nm))then
          jml = nm                                                              !These two lines are added to apply the effect of magnitude saturation
          jmu = nm                                                              !These two lines are added to apply the effect of magnitude saturation
          iflag(1,2)=5
      endif

c     find bounding distance values.
      rl = alog10(r)
      jdl = 1
      jdu = 1
      fracd=0.
      do j = 1, nd
          if(abs(rl-rlog(j)).lt.0.01)then
              jdl = j
              jdu = j
              iflag(1,3)=0
              exit
          endif
      enddo
      if(iflag(1,3).eq.4.and.rl.lt.rlog(nd).and.rl.ge.rlog(1))then
          do j=1,nd-1
              if(rl .ge. rlog(j) .and. rl .lt. rlog(j+1))then
                 jdl=j
                 jdu=j + 1
                 fracd = (rl - rlog(jdl)) / (rlog(jdl+1)-rlog(jdl))
              endif
          enddo
          iflag(1,3)=1
      elseif(iflag(1,3).eq.4.and.rl.gt.rlog(nd))then
          jdl = nd
          jdu = nd
          iflag(1,3)=5
      endif

      return
      end



      subroutine calc_hazard(n_s, n_c, n_g)
      use datastruct
      integer                                   :: n_s, n_c, n_g
      integer                                   :: i, j, num
      real                                      :: arl, aru
      real                                      :: dis_, mag_
      real, dimension(2)                        :: avg
      num=0
      M_R_limit(n_c,n_g).nline=cat(n_c).n_rec
      call boxmul_arr(haz(n_c,n_g).gmpar(1:cat(n_c).n_rec,2),
     1                cat(n_c).n_rec)
     
      do i=1,M_R_limit(n_c,n_g).nline
          haz(n_c,n_g).gmpar(i,1)=-20.
c          dis_=M_R_limit(n_c,n_g).dist(i,2)                                     
          dis_=M_R_limit(n_c,n_g).dist(i,gmp(n_g).disflag+1)                    !This makes the dis_ to be hypocentral distance and will be used ac a measure for accepting this record in catalog
          mag_=cat(n_c).mag(i)
          if(mag_.lt.sit(n_s).rng(1,1))cycle
          if(mag_.gt.sit(n_s).rng(2,1))cycle
          if(dis_.lt.sit(n_s).rng(1,2))cycle
          if(dis_.gt.sit(n_s).rng(2,2))cycle
          num=num+1
          do j=f_ind(i,1),f_ind(i,2)
              arl = gmp(n_g).table(j,r_ind(i,1),m_ind(i,1))+m_frac(i)*
     1             (gmp(n_g).table(j,r_ind(i,1),m_ind(i,2))-
     2              gmp(n_g).table(j,r_ind(i,1),m_ind(i,1)))
              aru = gmp(n_g).table(j,r_ind(i,2),m_ind(i,1))+m_frac(i)*
     1             (gmp(n_g).table(j,r_ind(i,2),m_ind(i,2))-
     2              gmp(n_g).table(j,r_ind(i,2),m_ind(i,1)))
              avg(j) = arl + r_frac(i)*(aru-arl)
          enddo
          haz(n_c,n_g).gmpar(i,1) = avg(f_ind(i,1)) + 
     1                 f_frac(i) * (avg(f_ind(i,2)) - avg(f_ind(i,1)))
          if(M_R_limit(n_c,n_g).flags(i,3).eq.5)then
              haz(n_c,n_g).gmpar(i,1)=haz(n_c,n_g).gmpar(i,1)+
     1                                gmp(n_g).dis(gmp(n_g).n_dis)-
     2        log10(M_R_limit(n_c,n_g).dist(i,gmp(n_g).disflag+1))
          endif
          gmp(n_g).sigma_model=gmp(n_g).GMsig(f_ind(i,1))+f_frac(i)*
     1                        (gmp(n_g).GMsig(f_ind(i,2))-
     2                         gmp(n_g).GMsig(f_ind(i,1)))
          
!          haz(n_c,n_g).gmpar(i,1)=haz(n_c,n_g).gmpar(i,1)
!     1                           +haz(n_c,n_g).gmpar(i,2)*
!     2                            (gmp(n_g).sigma+gmp(n_g).sigma_model)
!     3                           +gmp(n_g).const
          haz(n_c,n_g).gmpar(i,1)=haz(n_c,n_g).gmpar(i,1)
     1                           +haz(n_c,n_g).gmpar(i,2)*
     2                            gmp(n_g).sigma_model
     3                           +gmp(n_g).const
          haz(n_c,n_g).gmpar(i,1)=10**(haz(n_c,n_g).gmpar(i,1))
      enddo
      M_R_limit(n_c,n_g).nline_inrange=num
      return
      end



      subroutine calc_hazall(n_s, n_c)
      use datastruct
      integer                                   :: n_s, n_c, n_subcat           !This change is for creating constant parameters for a zone in a subcatalog-April-2012
      integer                                   :: i, j, ngmpe_pick
      real                                      :: rands
      real                                      :: sum
      if(cat(n_c).epist.eq.0)then
      n_subcat=0                                                                !This change is for creating constant parameters for a zone in a subcatalog-April-2012
          do i=1,cat(n_c).n_rec
              if(n_subcat.ne.cat(n_c).cat_no(i))then                            !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  n_subcat=cat(n_c).cat_no(i)                                   !This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  sum=0                                                         !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  call rand_num_gen(rands,0)                                    !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  do j=1,ngmpe                                                  !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                      ngmpe_pick=j                                              !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                      sum=sum+gmp(j).wgt                                        !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                      if(sum.gt.rands)exit                                      !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  enddo                                                         !Indented: 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
              hazall.gmpar(i,:)=haz(n_c,ngmpe_pick).gmpar(i,:)
              hazall.res(i,1)  =real(ngmpe_pick)
              hazall.res(i,2)  =rands
          enddo
      else
          do i=1,cat(n_c).n_rec
              sum=0                                                             !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
              call rand_num_gen(rands,0)                                        !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
              do j=1,ngmpe                                                      !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  ngmpe_pick=j                                                  !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  sum=sum+gmp(j).wgt                                            !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
                  if(sum.gt.rands)exit                                          !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
              enddo                                                             !Indented: This change is for creating constant parameters for a zone in a subcatalog-April-2012
              hazall.gmpar(i,:)=haz(n_c,ngmpe_pick).gmpar(i,:)
              hazall.res(i,1)  =real(ngmpe_pick)
              hazall.res(i,2)  =rands
          enddo
      endif
      return
      end



      subroutine make_output(n_s, n_c, n_g, c3)
      use datastruct
      integer                                   :: n_s, n_c, n_g
      real                                      :: dis_, mag_
      character (len=3)                         :: c3
      character (len=15)                        :: GMtype
      call findGMtype(sit(n_s).freq,GMtype)

      open(unit=10,file=trim(output(n_s, n_c, n_g))//'.'//c3
     1            ,status='unknown')
      write(10,7)n_s
7     format('Site number:           ',i4)
      write(10,8)n_c
8     format('Catalog number:        ',i4)
      write(10,9)n_g
9     format('GM model number:       ',i4)
      write(10,10)sit(n_s).coord(1)
10    format('Site latitude:         ',f9.3)
      write(10,11)sit(n_s).coord(2)
11    format('Site longitude:        ',f9.3)
      write(10,12)GMtype
12    format('Ground motion:         ',a14)
      write(10,13)gmp(n_g).sigma_model
13    format('Applied GM sigma:      ',f5.2)
      write(10,14)10**gmp(n_g).const
14    format('Multiplication factor: ',f5.2)
      write(10,15)gmp(n_g).gmpname
15    format('GM table file name:    ',a50)
      write(10,16)cat(n_c).catname
16    format('Catalog file name:     ',a50)
      write(10,17)
17    format('Catalog file header:')
      write(10,18)(cat(n_c).cathead(1:7))
18    format(3x,'|--->',3x,a80)
      write(10,*)
      write(10,41)GMtype
41    format(' Subcatalog No| Time (years) |   Latitude   |   Longitude 
     1 |   Depth(km)  |   Magnitude  |  Zone Number |  N0 per year |    
     2 Beta     |     Mmax     |Distance(Epi.)|Distance(Hyp.)|Distance(R
     3cd.)|Distance(Rjb.)|',a15,' Times sigma  | Flags: F M D |')
      do i=1,M_R_limit(n_c,n_g).nline
          dis_=M_R_limit(n_c,n_g).dist(i,gmp(n_g).disflag+1)
          mag_=cat(n_c).mag(i)
          if(mag_.lt.sit(n_s).rng(1,1))cycle
          if(mag_.gt.sit(n_s).rng(2,1))cycle
          if(dis_.lt.sit(n_s).rng(1,2))cycle
          if(dis_.gt.sit(n_s).rng(2,2))cycle
          write(10,'(16(g13.6,2x),3(2x,i1,2x))')
     1          cat(n_c).cat_no(i),      cat(n_c).time(i),
     2          cat(n_c).loc(i,:),       cat(n_c).mag(i),
     3          cat(n_c).zone_no(i),     cat(n_c).N0beta(i,:),
     4          cat(n_c).Mmax(i),        M_R_limit(n_c,n_g).dist(i,1:4),
     5          haz(n_c,n_g).gmpar(i,:), M_R_limit(n_c,n_g).flags(i,:)
      enddo
      close(unit=10)
      return
      end



      subroutine findGMtype(freq,GMtype)
      character*15 GMtype
      if(freq.eq.99.)then
          GMtype=' PGA(cm/sec^2)|'
      elseif(freq.eq.89.)then
          GMtype='  PGV(cm/sec) |'
      else
          write(GMtype,'(a6,f5.2,a4)')' PSA (',freq,'Hz)|'
      endif
      return
      end



      subroutine make_output_w(n_s, n_c, c3)
      use datastruct
      integer                                   :: n_s, n_c
      integer                                   :: good_num
      real                                      :: dis_, mag_
      character (len=3)                         :: c3
      character (len=15)                        :: GMtype
      good_num=0
      call findGMtype(sit(n_s).freq,GMtype)
      open(unit=10,file=trim(output_w(n_s, n_c))//'.'//c3
     1            ,status='unknown')
      write(10,7)n_s
7     format('Site number:           ',i4)
      write(10,8)n_c
8     format('Catalog number:        ',i4)
      write(10,10)sit(n_s).coord(1)
10    format('Site latitude:         ',f9.3)
      write(10,11)sit(n_s).coord(2)
11    format('Site longitude:        ',f9.3)
      write(10,12)GMtype
12    format('Ground motion:         ',a14)
      write(10,21)ngmpe
21    format('There are',i3,' GMPE tables with parameters: ')
      do i=1,ngmpe
          write(10,22)i
22        format('           GM model number:       ',i5)
          write(10,13)gmp(i).gmpname
13        format('           GM table file name:    ',a50)
          write(10,14)gmp(i).sigma_model
14        format('           Applied GM sigma:      ',f5.2)
          write(10,15)10**gmp(i).const
15        format('           Multiplication factor: ',f5.2)
      enddo
      write(10,*)
      write(10,16)cat(n_c).catname
16    format('Catalog file name:     ',a50)
      write(10,17)
17    format('Catalog file header:')
      write(10,18)(cat(n_c).cathead(1:7))
18    format(3x,'|--->',3x,a80)
      write(10,*)
      write(10,41)GMtype

41    format(' Subcatalog No| Time (years) |   Latitude   |   Longitude 
     1 |   Depth(km)  |   Magnitude  |  Zone Number |  N0 per year |    
     2 Beta     |     Mmax     |Distance(Epi.)|Distance(Hyp.)|Distance(R
     3cd.)|Distance(Rjb.)|',a15,' Times sigma  | GMPE model # |')
      do i=1,M_R_limit(n_c,1).nline
          dis_=M_R_limit(n_c,1).dist(i,gmp(1).disflag+1)
!          dis_=M_R_limit(n_c,1).dist(i,2)
          mag_=cat(n_c).mag(i)
          if(mag_.lt.sit(n_s).rng(1,1))cycle
          if(mag_.gt.sit(n_s).rng(2,1))cycle
          if(dis_.lt.sit(n_s).rng(1,2))cycle
          if(dis_.gt.sit(n_s).rng(2,2))cycle
c          write(10,'(15(g13.6,2x),3(2x,i1,2x))')
          write(10,'(16(g13.6,2x),(7x,i2))')
     1          cat(n_c).cat_no(i),      cat(n_c).time(i),
     2          cat(n_c).loc(i,:),       cat(n_c).mag(i),
     3          cat(n_c).zone_no(i),     cat(n_c).N0beta(i,:),
     4          cat(n_c).Mmax(i),        M_R_limit(n_c,1).dist(i,1:4),
     5          hazall.gmpar(i,:),       int(hazall.res(i,1))
          good_num=good_num+1
          good_gm(good_num)=hazall.gmpar(i,1)
      enddo
      close(unit=10)
      return
      end



      subroutine make_output_w_map(n_s, n_c, c3)
      use datastruct
      integer                                   :: n_s, n_c
      integer                                   :: good_num
      real                                      :: dis_, mag_
      character (len=3)                         :: c3
      character (len=15)                        :: GMtype
      good_num=0
      do i=1,M_R_limit(n_c,1).nline
          dis_=M_R_limit(n_c,1).dist(i,gmp(1).disflag+1)
          mag_=cat(n_c).mag(i)
          if(mag_.lt.sit(n_s).rng(1,1))cycle
          if(mag_.gt.sit(n_s).rng(2,1))cycle
          if(dis_.lt.sit(n_s).rng(1,2))cycle
          if(dis_.gt.sit(n_s).rng(2,2))cycle
          good_num=good_num+1
          good_gm(good_num)=hazall.gmpar(i,1)
      enddo
      return
      end



      subroutine make_hazard_output_w(n_s, n_c, c3)
      use datastruct
      integer                                   :: n_s, n_c
      integer                                   :: n, i, j
      real                                      :: dis_, mag_
      real                                      :: min_gm, max_gm
      character (len=3)                         :: c3
      character (len=15)                        :: GMtype
      call findGMtype(sit(n_s).freq,GMtype)
      call find_min_max(good_gm(1:M_R_limit(n_c,1).nline_inrange),
     1                  M_R_limit(n_c,1).nline_inrange,min_gm,max_gm)
!      write(*,*)M_R_limit(n_c,1).nline_inrange,min_gm,max_gm
      call calc_haz_curve(n_s, n_c)

      open(unit=10,file=trim(output_w(n_s, n_c))//'.'//c3
     1            ,status='unknown')
      write(10,7)n_s
7     format('Site number:           ',i4)
      write(10,8)n_c
8     format('Catalog number:        ',i4)
      write(10,10)sit(n_s).coord(1)
10    format('Site latitude:         ',f9.3)
      write(10,11)sit(n_s).coord(2)
11    format('Site longitude:        ',f9.3)
      write(10,12)GMtype
12    format('Ground motion:         ',a14)
      write(10,21)ngmpe
21    format('There are',i3,' GMPE tables with parameters: ')
      do i=1,ngmpe
          write(10,22)i
22        format('           GM model number:       ',i5)
          write(10,13)gmp(i).gmpname
13        format('           GM table file name:    ',a50)
          write(10,14)gmp(i).sigma_model
14        format('           Applied GM sigma:      ',f5.2)
          write(10,15)10**gmp(i).const
15        format('           Multiplication factor: ',f5.2)
      enddo
      write(10,*)
      write(10,16)cat(n_c).catname
16    format('Catalog file name:     ',a50)
      write(10,17)
17    format('Catalog file header:')
      write(10,18)(cat(n_c).cathead(1:7))
18    format(3x,'|--->',3x,a80)
      write(10,*)
      write(10,30)min_gm
30    format('Minimum ground motion: ',f9.3)
      write(10,31)max_gm
31    format('Maximum ground motion: ',f9.3)
      write(10,*)
      write(10,32)
32    format('      GM level      |  No of exceedance  | Rate of exceeda
     1nce | Prob. of exceedance /yr')
      write(10,33)
33    format(90('-'))
      do i=1,61
          if(haz_curve(i,1)<max_gm)then
              write(10,34)(haz_curve(i,j),j=1,4)
          endif
      enddo
c      write(10,34)(haz_curve(90,j),j=1,4)
34    format(4(3x,g14.6,4x))
      
      close(unit=10)
      return
      end



      subroutine make_hazard_map(n_s, n_c)
      use datastruct
      integer                                   :: n_s, n_c
      integer                                   :: n, i, j
      real                                      :: gm

      call sort(M_R_limit(n_c,1).nline_inrange,
     1          good_gm(1:M_R_limit(n_c,1).nline_inrange))
      n=int(cat(n_c).n_of_cats*cat(n_c).n_of_yrs/period_map)
      gm=good_gm(M_R_limit(n_c,1).nline_inrange-n+1)
      write(30,*)sit(n_s).coord,gm

      return
      end
      
      

      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,l,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=l-1
2         arr(i+1)=a
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
        if(arr(l).gt.arr(ir))then
          temp=arr(l)

          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(l+1))then
          temp=arr(l)
          arr(l)=arr(l+1)
          arr(l+1)=temp
        endif
        i=l+1
        j=ir
        a=arr(l+1)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l+1)=arr(j)
        arr(j)=a
        jstack=jstack+2

        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END



      subroutine calc_haz_curve(n_s, n_c)
      use datastruct
      integer                                   :: i, j, n_occur
      integer                                   :: n_s, n_c
      real                                      :: x_min, x_max
      do i=1,61
          haz_curve(i,1)=10.**(-2+0.1*real(i-1))
          n_occur=0
          do j=1,M_R_limit(n_c, 1).nline_inrange
              if(good_gm(j)>= haz_curve(i,1))n_occur=n_occur+1
          enddo
          haz_curve(i,2)=real(n_occur)
          haz_curve(i,3)=haz_curve(i,2)/
     1                   (cat(n_c).n_of_cats*cat(n_c).n_of_yrs)
          haz_curve(i,4)=1-exp(-haz_curve(i,3))
      enddo
      call find_min_max(good_gm,M_R_limit(n_c, 1).nline_inrange,
     1                  x_min,x_max)
      haz_curve(90,1)=x_max
      haz_curve(90,2)=1.
      haz_curve(90,3)=1./(cat(n_c).n_of_cats*cat(n_c).n_of_yrs)
      haz_curve(90,4)=1-exp(-haz_curve(90,3))
      return
      end



      subroutine find_min_max(arr,n,x_min,x_max)
      real, dimension(*)                        :: arr
      real                                      :: x_min, x_max
      x_min=arr(1)
      x_max=arr(1)
      do i=2,n
          if(arr(i)<x_min)x_min=arr(i)
          if(arr(i)>x_max)x_max=arr(i)
      enddo
      return
      end
      
      
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

cc             The following subroutines are kept for future use             cc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	v is an eight component matrix corresponding to values of the functio
c	n at corners of a box surrounding the point in question. p0 is the co
c	ordinates of initial point of box. dp is the vector of sizes of box s
c	ides, p is the coordinate of the study point, and value is the value 
c	at p.
c	This subroutine is not implemented in this version of program but res
c	erved for possible use in the future.
c
      subroutine interpolate3d(v,p0,dp,p,vlu)
      real					v(2,2,2),p0(3),dp(3),p(3)
      x=(p(1)-p0(1))/dp(1)
      y=(p(2)-p0(2))/dp(2)
      z=(p(3)-p0(3))/dp(3)
      vlu=v(1,1,1)*(1-x)*(1-y)*(1-z)+v(2,1,1)*x*(1-y)*(1-z)+
     1    v(1,2,1)*(1-x)*y*(1-z)    +v(1,1,2)*(1-x)*(1-y)*z+
     2    v(2,1,2)*x*(1-y)*z        +v(1,2,2)*(1-x)*y*z+
     3    v(2,2,1)*x*y*(1-z)        +v(2,2,2)*x*y*z
      return
      end



c
c     Function maxfinder finds the maximum value stored in 1-D array arr. The
c     parameter nelement is the number of elelemts (or size) of matrix arr.
c
      function maxfinder(arr,nelement)
      integer            :: maxfinder, arr(nelement)
      maxfinder=0
      do i=1,nelement
          maxfinder=max(arr(i),maxfinder)
      enddo
      return
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

cc                   Section of subprograms from fault_dist                   cc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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 readEQHAz1input(input)
      use srcdef
      integer(2)				delfiles 
      logical                   ex
      character*50 input,tempinput
      
      inquire(file=trim(input),exist=ex)
      if(not(ex))then
          src(1:60).typ=1
          write(*,*)'EqHaz1 input file: ',trim(input),' not found.'
          write(*,*)'All sources will be assumed to be areal.'
          return
      endif
      
      call getheadline(input,headline)                                          !Reading the header line of input. Added on 2011/10/26

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

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

      izone=0
      ifault=0
      nsource=nzones+nfaults
      if(nsource.ge.1)then
          do i=1,nsource
            read(21,*)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(21,*)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(21,*)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(21,*)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(21,*)(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(21,*)(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(21,*)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(21,*)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(21,*)(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(21,*)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(21,*)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(21,*)(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(21,*)(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(21,*)(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)
            src(i)=tmp_src

          enddo
      endif
      close(unit=21)
      delfiles=delfilesqq(tempinput)
      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=23,file=input,status='old')
	read(23,'(a500)')headline
	close(unit=23)
	return
	end



      subroutine find_distances()
      use srcdef
      use fltdis
      call find_flt_on_surf()                                                   !!!!!!This subroutine finds surface length, width, and azimuth of segments
      call find_eqk_on_fault(eqlat,eqlon,eqdep,eqmag,                           !!!!!
     1                       Lf,Wf,l,w,x,y)
      call find_asperity(Lf,Wf,l,w,x,y,Ax,Ay,Bx,By)
      call find_asperity_segments(Ax,Ay,Bx,By)                                  !!!!!
      Rjb_all=20000.
      Rcd_all=20000.
      do ix=1,source.nnode-1                                                     !!!!!
          if(source.isasperity(ix).eq.1)then                                     !!!!!
              call find_Rcd_Rjb(ix,plat,plon,phei,Rcd,Rjb)                       !!!!!
              if(Rjb.lt.Rjb_all)Rjb_all=Rjb
              if(Rcd.lt.Rcd_all)Rcd_all=Rcd
          endif
      enddo
      call find_Repi_Rhyp(plat,plon,eqlat,eqlon,eqdep,Repi,Rhyp)
      return
      end



      subroutine find_Repi_Rhyp(plat,plon,eqlat,eqlon,eqdep,Repi,Rhyp)
      implicit double precision (a-z)
      real                                      :: plat, plon
      real                                      :: eqlat,eqlon,eqdep
      real                                      :: Repi,Rhyp
      lat1=dble(plat)
      lon1=dble(plon)
      lat2=dble(eqlat)
      lon2=dble(eqlon)
      call xy2disazi(lat1,lon1,lat2,lon2,dist,azi)
      Repi=real(dist)
      Rhyp=sqrt(Repi**2+eqdep**2)
      return
      end



      subroutine find_Rcd_Rjb(n,plat,plon,phei,Rcd,Rjb)                         !!!!!This subroutine has its isource values cleaned
      use srcdef
      implicit double precision (a-z)
      dimension vec(3),mat(3,3)
      real                                      :: plat, plon, phei
      real                                      :: Rjb, Rcd
      integer                                   :: n
      asp_lat=rad(dble(source.asperity(n,1,1)))                                 !!!!!
      asp_lon=rad(dble(source.asperity(n,1,2)))                                 !!!!!
      asp_dep=dble(source.asperity(n,1,3))                                      !!!!!
      call find_asp_length(n,asp_len,asp_azi)                                   !!!!!
      call find_asp_width(n,asp_wid,asp_wid_image,asp_dip)                      !!!!!

      lat=rad(dble(plat))
      lon=rad(dble(plon))

      vec(1)=dcos(lat)*dcos(lon)
      vec(2)=dcos(lat)*dsin(lon)
      vec(3)=dsin(lat)
      call rot_matrix(-asp_lon,mat)
      call mult_rot_matrix(mat,vec,'z')
      call rot_matrix(+asp_lat,mat)
      call mult_rot_matrix(mat,vec,'y')
      call rot_matrix(+asp_azi-rad(dble(90.)),mat)
      call mult_rot_matrix(mat,vec,'x')
      call vec2xy(vec,lat,lon)
      call find_Rjb(lat,lon,asp_len,asp_wid_image,Rjb)
      
      vec(1)=6378.137*lon
      vec(2)=6378.137*lat
      vec(3)=dble(phei)-asp_dep
      call rot_matrix(-asp_dip,mat)
      call mult_rot_matrix(mat,vec,'x')
      call find_Rcd(vec(1),vec(2),vec(3),asp_len,asp_wid,Rcd)

      return
      end



c
c     This subroutine gets the coordinates of Ax,..,By as well as fault geometry
c     and calculates the coordinates of corners of rectangular areas on each seg
c     ment corresponding to the activated asperity.
c
      subroutine find_asperity_segments(Ax,Ay,Bx,By)                            !!!!!This subroutine has its isource values cleaned
      use srcdef
      implicit double precision (a-z)
      dimension vec(3),mat(3,3)
      real Ax,Ay,Bx,By,sum_len
      integer i, j, isegsum, iflag                                              !!!!!

      source.isasperity=0                                                       !!!!!
      source.asperity=0.                                                        !!!!!
      
      sum_len=0.
      do i=1,source.nnode-1                                                     !!!!!
          source.asperity(i,1,1)=0.                                             !!!!!
          source.asperity(i,1,2)=Ay                                             !!!!!
          source.asperity(i,2,1)=source.extra(i,1)                              !!!!!
          source.asperity(i,2,2)=Ay                                             !!!!!
          source.asperity(i,3,1)=source.extra(i,1)                              !!!!!
          source.asperity(i,3,2)=By                                             !!!!!
          source.asperity(i,4,1)=0                                              !!!!!
          source.asperity(i,4,2)=By                                             !!!!!
          if(Ax.ge.sum_len.and.
     1    Ax.le.(sum_len+source.extra(i,1)))then                                !!!!!
              source.isasperity(i)=1                                            !!!!!
              source.asperity(i,1,1)=Ax-sum_len                                 !!!!!
              source.asperity(i,1,2)=Ay                                         !!!!!
              source.asperity(i,4,1)=Ax-sum_len                                 !!!!!
              source.asperity(i,4,2)=By                                         !!!!!
          endif
          if(Bx.ge.sum_len.and.
     1    Bx.le.(sum_len+source.extra(i,1)))then                                !!!!!
              source.isasperity(i)=1                                            !!!!!
              source.asperity(i,2,1)=Bx-sum_len                                 !!!!!
              source.asperity(i,2,2)=Ay                                         !!!!!
              source.asperity(i,3,1)=Bx-sum_len                                 !!!!!
              source.asperity(i,3,2)=By                                         !!!!!
          endif
          sum_len=sum_len+source.extra(i,1)                                     !!!!!
      enddo
      
      isegsum=0
      do i=1,source.nnode-1                                                     !!!!!
          isegsum=isegsum+source.isasperity(i)                                  !!!!!
      enddo
      if(isegsum.eq.2)then
          iflag=0
          do i=1,source.nnode-1                                                 !!!!!
              if(source.isasperity(i).eq.0.and.iflag.eq.1)then                  !!!!!
                  source.isasperity(i)=1                                        !!!!!
              elseif(source.isasperity(i).eq.1)then                             !!!!!
                  iflag=1-iflag
              endif
          enddo
      endif

      do i=1,source.nnode-1                                                     !!!!!
          if(source.isasperity(i).eq.1)then                                     !!!!!
              do j=1,4
                  source.asperity(i,j,3)=                                       !!!!!
     1            source.asperity(i,j,2)*                                       !!!!!
     2               sin(deg2rad(source.geom(i,2)))-                            !!!!!!!!Probable error correction
     3                           source.geom(i,3)                               !!!!!!!!Probable erroe correction
                  vec(1)=
     1            dcos(dis2angle(dble(source.asperity(i,j,1))))                 !!!!!
     2           *dcos(dis2angle(dble(source.asperity(i,j,2)                    !!!!!
     3           *cos(deg2rad(source.geom(i,2))))))                             !!!!!
                  vec(2)=
     1            dsin(dis2angle(dble(source.asperity(i,j,1))))                 !!!!!
     2           *dcos(dis2angle(dble(source.asperity(i,j,2)                    !!!!!
     3           *cos(deg2rad(source.geom(i,2))))))                             !!!!!
                  vec(3)=
     1            dsin(dis2angle(dble(source.asperity(i,j,2)                    !!!!!
     2           *cos(deg2rad(source.geom(i,2))))))                             !!!!!

                  call rot_matrix
     1                 (+rad(dble(90-source.extra(i,3))),mat)                   !!!!!
                  call mult_rot_matrix(mat,vec,'x')
                  call rot_matrix
     1                 (-rad(dble(source.nodes(i,1))),mat)                      !!!!!
                  call mult_rot_matrix(mat,vec,'y')
                  call rot_matrix
     1                 (+rad(dble(source.nodes(i,2))),mat)                      !!!!!
                  call mult_rot_matrix(mat,vec,'z')
                  call vec2xy(vec,lat,lon)
                  source.asperity(i,j,1)=real(deg(lat))                         !!!!!
                  source.asperity(i,j,2)=real(deg(lon))                         !!!!!
              enddo
          endif
      enddo
      
      return
      end
      
      
            
c
c     This subroutine gets the geometry of multisegment fault, coordinates and m
c     agnitude of an earthquake that happened on fault surface, and finds the va
c     lues of Lf, Wf, l, w, x, and y of the source and hypocenter. These values 
c     are used later by other subroutines to derive the coordinates of corners o
c     f rectangular areas on each segment corresponding to activated asperity. F
c     or the calculations of asperity sizes Wells and Coppersmith equations are 
c     used. 
      subroutine find_eqk_on_fault(eqlat,eqlon,eqdep,eqmag,                     !!!!!This subroutine has its isource values cleaned
     1                             Lf,Wf,l,w,x,y)
      use srcdef
      real                                      :: Lf,Wf,l,w,x,y
      Lf=source.extra(source.nnode,1)                                           !!!!!
      Wf=source.geom(1,1)                                                       !!!!!
      l=10**(-2.44+0.59*eqmag)
      w=10**(-1.01+0.32*eqmag)
      call find_hypo(eqlat,eqlon,eqdep,x,y)                                     !!!!!
      return
      end



      subroutine find_asperity(Lf,Wf,l,w,x,y,Ax,Ay,Bx,By)
c     This subroutine takes the inputs about the following parameter:
c        Lf: Rectangular fault length
c        Wf: Rectangular fault width
c        l:  Length of asperity
c        w:  Width of asperity
c        x and y: Distance of the hypocenter (asperity center) along length and
c                 width of fault
c     respect to fault and returns the coordinates of corners of asperity on 
c     fault surface with corner of the fault. Ax, Ay, Bx, and By are the coordin
c     ates of asperity corners:
c
c     Fault:
c       (0,0) -------------------------------------------------------(Lf,0)
c             |                                                     |
c             |                        (Ax,Ay)  l                   |
c             |                            |---------|              |
c             |                            |         |              | 
c             |                          w |    *    |              |
c             |                            |  (x,y)  |              |
c             |                            |---------|              |
c             |                                     (Bx,By)         |
c             |                                                     |
c             -------------------------------------------------------
c      (0,-Wf)                                                      (Lf,-Wf)
c     This subroutine follows Gail Atkinson's suggested approach:
c     1-For asperity size larger than fault size the whole fault will be active.
c     2-If width of asperity gets larger than fault width, the maximum width wil
c       l be the fault width and the length of asperity will grow to conserve th
c       e asperity area.
c     3-If length of asperity gets larger than fault length, the maximum length 
c       will be the fault length and the width of asperity will grow to conserve 
c       the asperity area.
c     4-If the asperity sides fall outside the fault borders, the whole asperity
c       will be translated so that its edges will match the edges of the fault a
c       rea.
      real                                      :: Lf,Wf,l,w,x,y
      real                                      :: Ax,Ay,Bx,By
      real                                      :: leff,weff                    !Effective length and width of asperity
      leff=l
      weff=w
      
      if((leff*weff)>=(Lf*Wf))then
          Ax=0
          Ay=0
          Bx=Lf
          By=-Wf
          return
      endif
      
      if(weff>=Wf)then
          weff=Wf
          leff=l*w/Wf
          Ax=x-leff/2.
          Ay=0.
          Bx=x+leff/2.
          By=-weff
          if(Ax<0.)then
              Ax=0.
              Bx=leff
          endif
          if(Bx>Lf)then
              Ax=Lf-leff
              Bx=Lf
          endif
          return
      endif

      if(leff>=Lf)then
          leff=Lf
          weff=l*w/Lf
          Ax=0.
          Ay=y+weff/2.
          Bx=leff
          By=y-weff/2.
          if(Ay>0)then
              Ay=0.
              By=-weff
          endif
          if(By<-Wf)then
              Ay=-Wf+weff
              By=-Wf
          endif
          return
      endif
      
      Ax=x-leff/2.
      Ay=y+weff/2.
      Bx=x+leff/2.
      By=y-weff/2.
      if(Ax<0.)then
          Bx=Bx-Ax
          Ax=0.
      endif
      if(Ay>0.)then
          By=By-Ay
          Ay=0.
      endif
      if(Bx>Lf)then
          Ax=Ax-(Bx-Lf)
          Bx=Lf
      endif
      if(By<-Wf)then
          Ay=Ay-(By+Wf)
          By=-Wf
      endif
      
      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()                                             !!!!!This subroutine has its isource values cleaned
	use					srcdef
	double precision	lat1,lon1,lat2,lon2,dist,azi
      source.extra(source.nnode,1)=0.0                                          !!!!!
      do i=1,source.nnode-1                                                     !!!!!
          lat1=dble(source.nodes(i,1))                                          !!!!!
          lon1=dble(source.nodes(i,2))                                          !!!!!
          lat2=dble(source.nodes(i+1,1))                                        !!!!!
          lon2=dble(source.nodes(i+1,2))                                        !!!!!
          call xy2disazi(lat1,lon1,lat2,lon2,dist,azi)
          source.extra(i,1)=real(dist)                                          !!!!!!Save length of fault segment
          source.extra(i,2)=source.geom(i,1)*                                   !!!!!
     1                      cos(deg2rad(source.geom(i,2)))                      !!!!!!Save surface projection of fault segment width
          source.extra(i,3)=real(azi)                                           !!!!!!Save azimuth of fault segment width
              source.extra(source.nnode,1)=                                     !!!!!
     1        source.extra(source.nnode,1)+                                     !!!!!
     2        source.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



      subroutine find_Rcd(xlon,xlat,xdep,len,wid,Rcd)
      implicit double precision (a-z)
      real                                      :: Rcd
      x=0
      if(xlon.le.0)x=xlon
      if(xlon.ge.len)x=xlon-len
      y=0
      if(xlat.ge.0)y=xlat
      if(xlat.le.(-wid))y=xlat+wid
      Rcd=real(dsqrt(x**2+y**2+xdep**2))
      return
      end


      subroutine find_Rjb(lat,lon,len,wid,Rjb)
      implicit double precision (a-z)
      real                                      :: Rjb
      xlon=6378.137*lon
      xlat=6378.137*lat
      x=0.
      if(xlon.le.0.)x=xlon
      if(xlon.ge.len)x=xlon-len
      y=0.
      if(xlat.ge.0.)y=xlat
      if(xlat.le.(-wid))y=xlat+wid
      Rjb=real(dsqrt(x**2.+y**2.))
      return
      end


      subroutine find_asp_length(n,asp_len,asp_azi)                             !!!!!This subroutine has its isource values cleaned
      use srcdef
      implicit double precision (a-z)
      integer                                   :: n                            !!!!!
      lat1=dble(source.asperity(n,1,1))                                         !!!!!
      lon1=dble(source.asperity(n,1,2))                                         !!!!!
      lat2=dble(source.asperity(n,2,1))                                         !!!!!
      lon2=dble(source.asperity(n,2,2))                                         !!!!!
	call xy2disazi(lat1,lon1,lat2,lon2,asp_len,azimuth)
	asp_azi=rad(azimuth)
	return
	end



      subroutine find_asp_width(n,asp_wid,asp_wid_image,asp_dip)                !!!!!This subroutine has its isource values cleaned
      use srcdef
      implicit double precision (a-z)
      integer                                   :: n                            !!!!!
      lat1=dble(source.asperity(n,1,1))                                         !!!!!
      lon1=dble(source.asperity(n,1,2))                                         !!!!!
      lat2=dble(source.asperity(n,4,1))                                         !!!!!
      lon2=dble(source.asperity(n,4,2))                                         !!!!!
	call xy2disazi(lat1,lon1,lat2,lon2,asp_wid_image,azimuth)
      delta_depth=dabs(dble(source.asperity(n,4,3)-                             !!!!!
     1                      source.asperity(n,1,3)))                            !!!!!
      if(asp_wid_image.le.0.01)then
          asp_dip=2.*datan(dble(1.))
          asp_wid=delta_depth
      else
          asp_dip=datan(delta_depth/asp_wid_image)
          asp_wid=dsqrt(delta_depth**2+asp_wid_image**2)
      endif
      return
      end
      
      
            
      subroutine find_hypo(eqlat,eqlon,eqdep,x,y)                               !!!!!This subroutine has its isource values cleaned
      use srcdef
      implicit double precision (a-z)
      dimension vec(3),mat(3,3)
      real eqlat,eqlon,eqdep,x,y
      integer i
      x=0.
      y=0.
      do i=1,source.nnode-1                                                     !!!!!
          vec(1)=dcos((rad(dble(eqlat))))*dcos(rad(dble(eqlon)))
          vec(2)=dcos(rad(dble(eqlat)))*dsin(rad(dble(eqlon)))
          vec(3)=dsin(rad(dble(eqlat)))
!          write(*,*)'Eqlat=',eqlat,'EqLon=',eqlon
          call rot_matrix(-rad(dble(source.nodes(i,2))),mat)                    !!!!!
          call mult_rot_matrix(mat,vec,'z')
          call rot_matrix(+rad(dble(source.nodes(i,1))),mat)                    !!!!!
          call mult_rot_matrix(mat,vec,'y')
          call rot_matrix(+rad(dble(source.extra(i,3)-90)),mat)                 !!!!!
          call mult_rot_matrix(mat,vec,'x')
          call vec2xy(vec,lat,lon)
          if(6378.137*lon.ge.0.and.
     1       6378.137*lon.le.source.extra(i,1))then                             !!!!!
             x=x+6378.137*lon
             if(source.geom(i,2).eq.90.)then                                    !!!!!
                 y=source.geom(i,3)-eqdep                                       !!!!!
             else
                 y=-abs(real(lat)*6378.137)/
     1              cos(deg2rad(source.geom(i,2)))                              !!!!!
             endif
             exit
          else
          x=x+source.extra(i,1)                                                 !!!!!
          endif
      enddo
      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 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 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 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	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	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 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	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



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