PROGRAM SeismicityGRD_2_EQC !A utility program for generating one (or more) synthetic seismic catalogs of shallow earthquakes, ! expressed in Peter Bird's .EQC catalog format, based on a long-term forecast of shallow ! seismicity in Peter Bird's .GRD format. !Note that this is a simple Poissonian forecast in which the events have no clustering ! in time (i.e., no aftershock sequences). !By Peter Bird, UCLA, 2017.10.23 and revised 2018.01.05, ! re-using some code from my application PseudoCSEP. !========================= IMSL option ================================================================== !USE Numerical_Libraries ! needed for RNSET, RNUN, RNUNF, & ANORIN: ! (1) CALL RNSET(234579) ! initialize random-number generator. ! (2a) CALL RNUN(length, vector) ! fills "vector" with random numbers "p" out to "length" ! (2b) p = RNUNF() ! returns a single random number (less efficient). !N.B. The results of these calls are uniformly distributed (0.0, 1.0), and can be converted to !values from a standard normal distribution (in approximately (-4, +4)) by the "probit" function. ! (aka "inverse cumulative distribution function"). !The probit function is called "ANORIN" in the IMSL library: Probit(p) = ANORIN(p) = 2**(1/2) * ERFI(2*p - 1) !======Lines above commented-out when compiling on MEGAMIND, which does not have access to IMSL.========= !=========================== Intel MKL-VSL option ====================================================== !Using Intel's Math Kernel Library (MKL); specifically, the Vector Statistical Library (VSL) part. !The following normal-distribution generators are what I will need: ! TYPE SPECIFICATIONS: ! INTEGER :: status ! TYPE(VSL_STREAM_STATE) :: my_stream ! REAL, DIMENSION(1) :: MKL_VSL_random ! .... ! INITIALIZATION: ! status = vslnewstream( stream = my_stream, brng = VSL_BRNG_MCG31, seed = 1) ! IF (status /= VSL_STATUS_OK) THEN ... {stop} ! status == 0 means OK ! BOXCAR/UNIFORM DISTRIBUTION (between 0. and 1.): ! status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = MKL_VSL_random, a = 0.0, b = 1.0 ) ! IF (status /= VSL_STATUS_OK) THEN ... {stop} ! status == 0 means OK ! GAUSSIAN DISTRIBUTION (with mean of zero and standard deviation of 1.): ! status = vsRngGaussian( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = MKL_VSL_random, a = 0.0, sigma = 1.0 ) ! IF (status /= VSL_STATUS_OK) THEN ... {stop} ! status == 0 means OK !INCLUDE 'mkl_vsl.f90' <== while logically necessary to include this, it was not syntactically correct ! to declare a MODULE inside a PROGRAM. So I saved this file in my own folder, ! and added it explicitly as a component file for the project. !USE mkl_vsl_type ! for uniform-distribution and normal-distribution generators USE mkl_vsl ! for uniform-distribution and normal-distribution generators !======================================================================================================== IMPLICIT NONE CHARACTER*2 :: c2 CHARACTER*3 :: c3 CHARACTER*4 :: c4 CHARACTER*5 :: c5 CHARACTER*7 :: c7 CHARACTER*8 :: c8 CHARACTER*58 :: build_EQC_line, spare_EQC_line CHARACTER*58, DIMENSION(:), ALLOCATABLE :: EQC_line_store CHARACTER*127 :: GRD_filename ! SeismicityGRD file name (input) CHARACTER*127 :: EQC_filename ! to be constructed, for output INTEGER, PARAMETER :: m = 10 ! m = number of randomly-generated simulated earthquake catalogs for this forecast. INTEGER, PARAMETER :: Factorial_limit = 170 ! Length of R8 Factorial(0:Factorial_limit) and LOG_Factorial(0:Factorial_limit) vectors; ! should be larger than the greatest number of EQs in any one cell (actual, or simulated). ! HOWEVER, if set too high, you will get immediate overflow abend on start-up. ! You can determine the highest integer that a particular Fortran compiler will ! support by trial-and-error, as any error occurs right at start of the run. ! My own experiments (using Intel Visual Fortran Composer XE 2013, with Math Kernel Library) ! showed that a value of 170 is acceptable, but that a value of 171 causes an abend. INTEGER, PARAMETER :: comb_limit = 140 ! Maximum value for the integer number of earthquakes in one spatial cell !(and thus, independent variable in Poisson PDF and CDF) during simulations of ! virtual catalogs. Must be larger than the maximum number of earthquakes expected ! in the most-active cell of any forecast being scored. ! One hard limit is that comb_limit is not allowed to exceed Factorial_limit; ! this will be checked during execution, and SeismicityGRD_2_EQC will stop. ! Another problem is that, depending on the input data, ! SOME RUNS MAY EXPERIENCE floating-overflow ABENDS ! while computing lamba**omega, as omega goes to its upper limit of comb_limit, ! if comb_limit is too high. Therefore, comb_limit should be less than Factorial_limit, ! but it is not possible to say exactly how much less; it depends on the problem, ! because the input data for the problem determine the value(s) of lambda. ! Here we suggest a value based on some limited experimentation with long catalogs... ! Please let me (Peter Bird) know if this separation is not enough for your problem! ! NOTE: If your usage of SeismicityGRD_2_EQC is encountering problems with the limits above, try: ! *raising your magnitude threshold, OR ! *decreasing the size (d_lon, d_lat) of the spatial cells in your seismicity .GRD file, OR ! *shortening your simulated catalog duration (number of years), ! as any/all of these will reduce the numbers of earthquakes considered. INTEGER :: col, day_integer, depth_km_integer, hour_integer, i, iBin, ios, j, k, lines, & & minute_integer, month_integer, nGRDrows, nGRDcols, N, N_hat, & & omega, row, seconds_integer, trimmed_length, unmasked_cells, year1, year999, year_integer INTEGER*2, DIMENSION(:,:,:), ALLOCATABLE :: sim_N_hat_3D_matrix !---Parameters needed for statistical package Intel MKL/VSL:------------------------------- INTEGER(4), PARAMETER :: I4_one = 1 ! needed as input to some MKL subprograms, if default INTEGER ever changes to INTEGER(8). INTEGER(4) :: status ! specifying (4) in case the default INTEGER type ever changes TYPE(VSL_STREAM_STATE) :: my_stream !--------------------------------------------------------------------------------------- LOGICAL :: using_mask = .FALSE. ! I might program in this capability later(?) LOGICAL*1, DIMENSION(:,:), ALLOCATABLE :: cell_in_use ! same size and shape as forecast grid; will be re-sampled from mask grid REAL :: depth_km, GBs_requested, grid_lat, grid_lon, & & lat_max, lat_min, lon_max, lon_min, & & mean, max_EQ_depth_km, min_EQ_depth_km, new_EQC_duration_years, sigma REAL, DIMENSION(:), ALLOCATABLE :: N_hat_vector ! (k = 1, ..., m is the simulation #) REAL, DIMENSION(:,:), ALLOCATABLE :: dt_da_factor ! (row, col); converts from EQ rates in #/m^2/s to EQ #'s, in one cell, in test time window. REAL, DIMENSION(m) :: sim_temporary ! holds probabilities from which EQ #s will be deduced !----------------------------------------------------------------------------------------- !Dummy vectors (of 1 entry each) are needed to sync with MKL subprograms: REAL, DIMENSION(1) :: VSCdfNorm_in, VSCdfNorm_out, VSCdfNormInv_in, VSCdfNormInv_out !--------------------------------------------------------------------------------------- DOUBLE PRECISION, PARAMETER :: degrees_per_radian = 57.2957795130823D0 DOUBLE PRECISION, PARAMETER :: radians_per_degree = 0.0174532925199433D0 DOUBLE PRECISION, PARAMETER :: s_per_year = 3.1557600D7 ! assuming 365.25 days/year, on average. DOUBLE PRECISION :: R_m = 6371.0D3 ! radius of Earth DOUBLE PRECISION :: area_m2, beta, corner_magnitude, corner_moment, d_lat, d_lon, d_magnitude, dayX, & & factor1, factor2, factor3, fraction, hourX, lambda, latitude, lat_factor, longitude, lost_probability, & & magnitude, minuteX, mislocation_km, mislocation_degrees_EW, mislocation_degrees_NS, monthX, & & omega_R8, partial_day, partial_hour, partial_minute, partial_month, partial_second, partial_year, & & secondsX, sum, spare_R8, threshold_magnitude, threshold_moment, total_probability, variance, yearX DOUBLE PRECISION, DIMENSION(0:comb_limit ) :: probability_of_integer, cumulative_probability DOUBLE PRECISION, DIMENSION(0:Factorial_limit) :: Factorial, LOG_Factorial ! We will precompute these, to save time. DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: forecast ! DOUBLE PRECISION because of small exponents, like E-22! DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: s_within_window DOUBLE PRECISION, DIMENSION(m) :: random_R8 DOUBLE PRECISION, DIMENSION(0:50) :: magnitude_list, moment_list, G_list !--------------------------------------------------------------------------------------- !Dummy vectors (of 1 entry each) are needed to sync with MKL subprograms: DOUBLE PRECISION, DIMENSION(1) :: VDCdfNorm_in, VDCdfNorm_out !--------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------- !To get results of random-number generation -------------------------------------------- !REAL :: RNUNF ! IMSL function; or (alternatively)... REAL, DIMENSION(1) :: MKL_VSL_random ! under MKL's VSL package option !----------------------------------------------------------------------------------------- !============================================================================================================ !-------------------------- begin executable code ----------------------------------------------------------- WRITE (*, "(' ')") WRITE (*, "(' PROGRAM SeismicityGRD_2_EQC')") WRITE (*, "(' ')") WRITE (*, "(' A utility program for generating one (or many)')") WRITE (*, "(' synthetic seismic catalogs of shallow earthquakes,')") WRITE (*, "(' expressed in Peter Bird''s .EQC catalog format,')") WRITE (*, "(' based on a long-term forecast of shallow')") WRITE (*, "(' seismicity in Peter Bird''s .GRD format.')") WRITE (*, "(' Note that this is a simple Poissonian forecast')") WRITE (*, "(' in which the events have no clustering')") WRITE (*, "(' in time (i.e., no aftershock sequences).')") WRITE (*, "(' By Peter Bird, UCLA, 2017.10.23, and revised 2018.01.05,')") WRITE (*, "(' re-using some code from my application PseudoCSEP.')") WRITE (*, "(' ')") CALL Pause() ! initialize random-number generator of CALL RNUN(length, vector) or X = RNUNF(). !----Initialize IMSL version of random numbers (below): ----------------------------------- !CALL RNSET(234579) ! seeding the random number generator !------------------------------------------------------------------------------------------ IF (comb_limit > Factorial_limit) THEN WRITE (*, "(' ERROR: Parameter comb_limit (',I4,') exceeds Factorial_limit (',I4').')") comb_limit, Factorial_limit CALL Pause() STOP END IF !----Initialize Intel MKL/VSL version of random numbers (below): status = VslNewStream( stream = my_stream, brng = VSL_BRNG_MCG31, seed = 1) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When initializing my_stream with Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !------------------------------------------------------------------------------------------- !Precompute DOUBLE PRECISION, DIMENSION(0:Factorial_limit) :: Factorial and LOG_Factorial, to save time: Factorial(0) = 1.0D0 LOG_Factorial(0) = 0.0D0 DO i = 1, Factorial_limit Factorial(i) = Factorial(i-1) * i ! WARNING: This line will cause overflow abends if Factorial_limit is too high! LOG_Factorial(i) = LOG_Factorial(i-1) + LOG(1.0D0 * i) !WRITE (*, *) i, Factorial(i), LOG_Factorial(i) ! <=== un-comment this line to see the progression toward an overflow abend END DO !------------------------------------------------------------------------------------------- !======================================================================================================== !Locate and open the input SeismicityGRD file: WRITE (*, *) WRITE (*, "(' Enter the filename of your Seismicity.GRD file: ')") READ (*, "(A)") GRD_filename OPEN (UNIT = 1, FILE = TRIM(GRD_filename), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This file is not present (in the current folder.)')") CALL Pause() STOP END IF trimmed_length = LEN_TRIM(GRD_filename) - 4 ! (not counting the ".GRD" part) !read headers of this GRD file: READ (1, *) lon_min, d_lon, lon_max REWIND (1) READ (1, "(106X, F6.3)", IOSTAT = ios) threshold_magnitude IF (ios == 0) THEN WRITE (*, "(' Threshold magnitude of this forecast = ',F6.3)") threshold_magnitude ELSE WRITE (*, "(' NOT able to read the threshold magnitude from the first header line;')") WRITE (*, "(' Please enter the threshold magnitude for this forecast: ')") READ (*, *) threshold_magnitude END IF threshold_moment = Moment(threshold_magnitude) !Finish reading the Seismicity.GRD, and memorize it: READ (1, *) lat_min, d_lat, lat_max nGRDrows = NINT(1.0 + ((lat_max - lat_min) / d_lat)) nGRDcols = NINT(1.0 + ((lon_max - lon_min) / d_lat)) WRITE (*, "(' This forecast has ',I6,' rows and ',I6,' columns.')") nGRDrows, nGRDcols ALLOCATE ( forecast(nGRDrows, nGRDcols) ) READ (1, *, IOSTAT = ios)((forecast(i, j), j = 1, nGRDcols), i = 1, nGRDrows) IF (ios == 0) THEN WRITE (*, "(' Forecast grid was successfully read.')") ELSE WRITE (*, "(' ERROR while attempting to read the forecast grid.')") CALL Pause() STOP END IF CLOSE(1) !Describe the frequence/magnitude/moment relation for simulated seismicity: WRITE (*, *) WRITE (*, "(' Simulated seismicity will have a frequency/magnitude/moment relation')") WRITE (*, "(' which follows the tapered Gutenberg-Richter relation, from')") WRITE (*, "(' equation (9) of Bird & Kagan [2004, BSSA].')") WRITE (*, "(' Enter the beta value to be used (suggestion: 0.63): ')") READ (*, *) beta IF ((beta <= 0.0D0).OR.(beta >= 1.0D0)) THEN WRITE (*, "(' ERROR: Beta must be BETWEEN 0.0 and 1.0.')") CALL Pause() STOP END IF WRITE (*, "(' Enter the corner magnitude to be used: ')") READ (*, *) corner_magnitude corner_moment = Moment(corner_magnitude) !Pre-build the frequency/magnitude/moment relation: magnitude_list(0) = threshold_magnitude moment_list(0) = threshold_moment G_list(0) = 1.0D0 d_magnitude = (9.5D0 - threshold_magnitude) / 50.0D0 ! bin-width, in magnitude units DO i = 1, 50 magnitude = threshold_magnitude + i * d_magnitude magnitude_list(i) = magnitude moment_list(i) = Moment(magnitude) G_list(i) = ((moment_list(i) / threshold_moment)**(-beta)) * & & EXP((threshold_moment - moment_list(i)) / corner_moment) END DO !Determine temporal extent (first & last years) for desired virtual catalog EQC: !convert EQ rates in SI units of #/m^2/s into absolute (real) numbers per test time-window: WRITE (*, *) WRITE (*, "(' Please define the length of the desired synthetic catalog:')") WRITE (*, "(' Enter the first year# (as an integer): ')") READ (*, *) year1 WRITE (*, "(' Enter the last year# (as an integer): ')") READ (*, *) year999 IF (year999 >= year1) THEN new_EQC_duration_years = 1 + (year999 - year1) ELSE WRITE (*, "(' ERROR: Last year # must be >= first year #.')") CALL Pause() STOP END IF !Describe depth range for earthquakes: WRITE (*, *) WRITE (*, "(' Enter minimum depth for earthquakes, in km: ')") READ (*, *) min_EQ_depth_km WRITE (*, "(' Enter maximum depth for earthquakes, in km: ')") READ (*, *) max_EQ_depth_km IF (max_EQ_depth_km < min_EQ_depth_km) THEN WRITE (*, "(' ERROR: Maximum must be larger than minimum!')") CALL Pause() STOP END IF !Describe typical (median) amount of earthquake mislocation: WRITE (*, *) WRITE (*, "(' Enter the median amount of earthquake mislocation, in kilometers: ')") READ (*, *) mislocation_km IF (mislocation_km < 0.0D0) THEN WRITE (*, "(' ERROR: Negative mislocation error is unphysical.')") CALL Pause() STOP END IF !Precompute conversion factor, from seismicity rate density (SI) to # of events in cell: ALLOCATE ( dt_da_factor(nGRDrows, nGRDcols) ) !first, dt = length of test time-window: dt_da_factor = new_EQC_duration_years * s_per_year !next, include # of square meters (neglecting edge effects) DO row = 1, NGRDrows latitude = lat_max - (row - 1) * d_lat area_m2 = Area_of_cell(d_lon, d_lat, latitude, R_m) dt_da_factor(row, 1:nGRDcols) = dt_da_factor(row, 1:nGRDcols) * area_m2 END DO GBs_requested = 2. * m * nGRDrows * nGRDcols / (1024.**3) IF (GBs_requested >= 1.0D0) THEN WRITE (*, *) WRITE (*, "(' ALLOCATing ',F8.3,' GB array sim_N_hat_3D_matrix;')") GBs_requested CALL Pause() END IF ALLOCATE ( sim_N_hat_3D_matrix(m, nGRDrows, nGRDcols) ) ! May be BIG! (e.g., For m = 1000, and 0.1-degree global grid: 12.07 GB) sim_N_hat_3D_matrix = 0 ALLOCATE ( N_hat_vector(m) ) !Conduct (k = 1, ..., m) simulations of hypothetical expected catalogs, based on each input forecast model: WRITE (*, *) WRITE (*, "(' Computing ',I5,' simulated catalogs for this forecast model:')") m WRITE (*, "(' -----------------------------------------------------------')") WRITE (*, "(' Forecast model = ', A)") TRIM(GRD_filename(1:50)) IF (using_mask) THEN WRITE (*, "(' ERROR: Cannot use masking because it has not been programmed yet.')") CALL Pause() STOP WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I10,' unmasked cells...')") m, unmasked_cells ELSE ALLOCATE ( cell_in_use(nGRDrows, nGRDcols) ) cell_in_use = .TRUE. ! whole matrix WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I4,' rows and ',I4,' columns...')") m, nGRDrows, nGRDcols END IF DO row = 1, nGRDrows DO col = 1, nGRDcols lambda = forecast(row, col) * dt_da_factor(row, col) ! Number of earthquakes expected by this forecast, in this particular bin. IF (cell_in_use(row, col)) THEN !Note that lambda is not an integer, but a real number, e.g., 0.00000123 or 37.443 IF (NINT(REAL(lambda + 3.0 * SQRT(lambda))) > comb_limit) THEN WRITE (*, *) WRITE (*, "(' ERROR: comb_limit = ',I5,' is too low when lambda = ',F10.2)") comb_limit, lambda WRITE (*, "(' Increase this parameter and recompile, or else ...')") WRITE (*, "(' If you don''t want to recompile, or if comb_limit is already ~140,')") WRITE (*, "(' (and higher values run more risk of floating-overflow abends!),')") WRITE (*, "(' there are some other kinds of work-arounds that you could try:')") WRITE (*, "(' *raise your magnitude threshold, OR')") WRITE (*, "(' *decrease the size of the spatial cells in your seismicity .GRD file, OR')") WRITE (*, "(' *shorten your simulated catalog duration (number of years),')") WRITE (*, "(' as any/all of these will reduce the #s of earthquakes simulated in one cell.')") WRITE (*, *) CALL Pause() STOP END IF !Compute probabilities of different integer outcomes, given this expectation: probability_of_integer = 0.0D0 ! whole REAL*8 vector initialized total_probability = 0.0D0 ! REAL*8 scalar initialized create_comb_2: DO omega = 0, comb_limit !--- IMSL version: ------------------------------------------------------------------------------ !probability_of_integer(omega) = (lambda ** omega) * DEXP(-lambda) / DFAC(omega) !Using DFAC(n) = n! from IMSL Library. !--- MKL does not have factorial; using my own: ------------------------------------------------- omega_R8 = omega ! convert INTEGER to REAL*8, changing method of exponentiation in next statement: IF (omega == 0) THEN factor1 = 1.0D0 ! which should be implied by the general formula; just making sure! ELSE IF (omega == 1) THEN factor1 = lambda ! which should be implied by the general formula; just making sure! ELSE factor1 = lambda ** omega_R8 END IF factor2 = DEXP(-lambda) factor3 = 1.0D0 / Factorial(omega) probability_of_integer(omega) = factor1 * factor2 * factor3 ! where my Factorial() is a pre-computed vector with a much higher limit on the argument! !------------------------------------------------------------------------------------------------ total_probability = total_probability + probability_of_integer(omega) END DO create_comb_2 IF ((total_probability < 0.99D0).OR.(total_probability > 1.01D0)) THEN WRITE (*, "(' ERROR: total_probability = ',ES12.4,' when lambda = ',ES12.4,' and comb_limit = ',I6)") total_probability, lambda, comb_limit CALL Pause() STOP END IF !Compute cumulative probability of different integer outcomes: cumulative_probability(0) = probability_of_integer(0) DO omega = 1, comb_limit cumulative_probability(omega) = cumulative_probability(omega - 1) + probability_of_integer(omega) END DO lost_probability = 1.0D0 - cumulative_probability(comb_limit) IF (lost_probability > 0.0D0) THEN IF (lost_probability >= 0.01D0) THEN WRITE (*, "(' ERROR: lost_probability = ', ES12.4, ' is unacceptable.')") lost_probability CALL Pause() STOP END IF cumulative_probability(0:comb_limit) = cumulative_probability(0:comb_limit) + lost_probability cumulative_probability(comb_limit) = 1.0D0 ! (This should be implied by the line above; just making sure.) END IF !sample this Poisson distribution by reverse look-up (starting at low end, which is most probable): !Get one random number (0, 1) per simulation (for this bin and forecast): !--- IMSL method: -------------------------------------------------------------------------------------- !CALL RNUN(m, sim_temporary) ! Some day, I should upgrade this to a REAL*8 function version... !random_R8 = sim_temporary ! REAL*8(1:m) vector = REAL(1:m) vector !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = m, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- DO k = 1, m ! (simulation #) search: DO omega = 0, comb_limit ! which upper limit we do not expect to reach! (Most bins will produce 0, 1, 2, ... ?) IF (random_R8(k) <= cumulative_probability(omega)) THEN N_hat = omega sim_N_hat_3D_matrix(k, row, col) = omega ! sampled INTEGER # of EQs for this bin, simulation, [& forecast] EXIT search END IF END DO search ! omega = 0, comb_limit IF (N_hat == comb_limit) WRITE (*, "(' WARNING: N_hat hit compiled-in upper limit of ',I4,'/cell.')") comb_limit END DO ! k = 1, m (simulation number) !score all the (k = 1, ..., m) simulations for this forecast #j (using all bins) ELSE sim_N_hat_3D_matrix(1:m, row, col) = NINT(lambda) ! don't bother to perturb rates that will not be used END IF ! cell_in_use(row, col), or not END DO ! col = 1, nGRDcols END DO ! row = 1, nGRDrows !Count total earthquakes in each simulation of this forecast: WRITE (*, "(' Totalling virtual earthquakes for each of these simulations...')") sum = 0.0D0 ! just initializing the total # EQs over ALL m simulations of this forecast: DO k = 1, m N = 0 ! just initializing the # in one particular virtual catalog... DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN N = N + sim_N_hat_3D_matrix(k, row, col) END IF END DO END DO N_hat_vector(k) = N sum = sum + N END DO !Output of EQC file(s): one, or many, depending on value of m: DO k = 1, m !Output EQC filename: WRITE (c3, "(I3)") MIN(k, 999) IF (c3(1:1) == ' ') c3(1:1) = '0' IF (c3(2:2) == ' ') c3(2:2) = '0' EQC_filename = GRD_filename(1:trimmed_length) // '_' // c3 // ".EQC" OPEN (UNIT = 2, FILE = TRIM(EQC_filename)) !Create text-lines of output file, initially without regard for their order: N_hat = N_hat_vector(k) ALLOCATE ( EQC_line_store(N_hat) ) EQC_line_store = ' ' ! just in case... ALLOCATE ( s_within_window(N_hat) ) s_within_window = 0.0D0 ! just in case... !Build and store lines of the EQC file: lines = 0 ! initially DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col).AND.(sim_N_hat_3D_matrix(k, row, col) > 0)) THEN grid_lat = lat_max - (row - 1) * d_lat grid_lon = lon_min + (col - 1) * d_lon n = sim_N_hat_3D_matrix(k, row, col) DO i = 1, n lines = lines + 1 ! within the new EQC file !Get one REAL*8 random number (0.0D0, 1.0D0) per earthquake, to be s_within_window: !--- IMSL method: -------------------------------------------------------------------------------------- !CALL DRNUN(1, random_R8) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF s_within_window(lines) = random_R8(1) !---------------------------------------------------------------------------------------------------------------------------------- !Build the line of output yearX = year1 + s_within_window(lines) * new_EQC_duration_years year_integer = yearX ! truncating any fraction partial_year = yearX - year_integer ! should be in range of 0.0 ~ 0.99999 monthX = 1.0D0 + partial_year * 12.0D0 month_integer = monthX ! truncating any fraction partial_month = monthX - month_integer dayX = 1.0D0 + partial_month * 30.0D0 ! FUDGE; treating all months as 30 days long! IF (month_integer == 2) THEN ! February IF (MOD(year_integer, 4) == 0) THEN ! leap year dayX = 1.0D0 + partial_month * 29.0D0 ELSE ! other year dayX = 1.0D0 + partial_month * 28.0D0 END IF ELSE IF (month_integer == 9) THEN ! September dayX = 1.0D0 + partial_month * 30.0D0 ELSE IF (month_integer == 4) THEN ! April dayX = 1.0D0 + partial_month * 30.0D0 ELSE IF (month_integer == 6) THEN ! June dayX = 1.0D0 + partial_month * 30.0D0 ELSE IF (month_integer == 11) THEN ! November dayX = 1.0D0 + partial_month * 30.0D0 ELSE ! months with 31 days dayX = 1.0D0 + partial_month * 31.0D0 END IF day_integer = dayX ! truncating any fraction partial_day = dayX - day_integer hourX = partial_day * 24.0D0 hour_integer = hourX ! truncating any fraction partial_hour = hourX - hour_integer minuteX = partial_hour * 60.0D0 minute_integer = minuteX ! truncating any fraction partial_minute = minuteX - minute_integer secondsX = partial_minute * 60.0D0 seconds_integer = secondsX ! truncating any fraction partial_second = secondsX - seconds_integer build_EQC_line(1:9) = "GRD_2_EQC" WRITE (c5, "(I5)") year_integer ! which might have a leading minus sign? build_EQC_line(10:14) = c5 build_EQC_line(15:15) = '.' WRITE (c2, "(I2)") month_integer IF (c2(1:1) == ' ') c2(1:1) = '0' build_EQC_line(16:17) = c2 build_EQC_line(18:18) = '.' WRITE (c2, "(I2)") day_integer IF (c2(1:1) == ' ') c2(1:1) = '0' build_EQC_line(19:20) = c2 build_EQC_line(21:21) = ' ' WRITE (c2, "(I2)") hour_integer IF (c2(1:1) == ' ') c2(1:1) = '0' build_EQC_line(22:23) = c2 build_EQC_line(24:24) = ':' WRITE (c2, "(I2)") minute_integer IF (c2(1:1) == ' ') c2(1:1) = '0' build_EQC_line(25:26) = c2 build_EQC_line(27:27) = ':' WRITE (c4, "(F4.1)") secondsX IF (c4(1:1) == ' ') c4(1:1) = '0' build_EQC_line(28:31) = c4 build_EQC_line(32:32) = ' ' longitude = grid_lon ! before any perturbation latitude = grid_lat !First, perturb EQ location within the Seismicity.GRD cell, using boxcar distributions: !Get one REAL*8 random number (0.0D0, 1.0D0) per earthquake, to perturb the longitude: !--- IMSL method: -------------------------------------------------------------------------------------- !CALL DRNUN(1, random_R8) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- longitude = longitude + (random_R8(1) - 0.5D0) * d_lon !Get one REAL*8 more random number (0.0D0, 1.0D0) per earthquake, to perturb the latitude !--- IMSL method: -------------------------------------------------------------------------------------- !CALL DRNUN(1, random_R8) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- latitude = latitude + (random_R8(1) - 0.5D0) * d_lat !Second, add the mislocation perturbation, using Gaussian (normal) distributions: !---------------------------------------------------------------------------------------------------------------------------------- ! GAUSSIAN DISTRIBUTION (with mean of zero and standard deviation of 1.): status = vdRngGaussian( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, sigma = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- mislocation_degrees_EW = (((mislocation_km * 1000.0D0) / (R_m * COS(latitude * radians_per_degree))) * degrees_per_radian) / SQRT(2.0D0) longitude = longitude + mislocation_degrees_EW * random_R8(1) !---------------------------------------------------------------------------------------------------------------------------------- ! GAUSSIAN DISTRIBUTION (with mean of zero and standard deviation of 1.): status = vdRngGaussian( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, sigma = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- mislocation_degrees_NS = (((mislocation_km * 1000.0D0) / R_m) * degrees_per_radian) / SQRT(2.0D0) latitude = latitude + mislocation_degrees_NS * random_R8(1) WRITE (c8, "(F8.3)") longitude build_EQC_line(33:40) = c8 build_EQC_line(41:41) = ' ' WRITE (c7, "(F7.3)") latitude build_EQC_line(42:48) = c7 build_EQC_line(49:49) = ' ' !Get one more REAL*8 random number (0.0D0, 1.0D0) per earthquake, to determine its depth !--- IMSL method: -------------------------------------------------------------------------------------- !CALL DRNUN(1, random_R8) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- depth_km = min_EQ_depth_km + random_R8(1) * (max_EQ_depth_km - min_EQ_depth_km) depth_km_integer = NINT(depth_km) WRITE (c3, "(I3)") depth_km_integer build_EQC_line(50:52) = c3 build_EQC_line(53:53) = ' ' !Get one more REAL*8 random number (0.0D0, 1.0D0) per earthquake, to determine its magnitude ! (by reverse look-up in the frequency-magnitude relation, relative to G values): !--- IMSL method: -------------------------------------------------------------------------------------- !CALL DRNUN(1, random_R8) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vdRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = random_R8, a = 0.0D0, b = 1.0D0 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When requesting random_R8 from Intel MKL/VSL, status = ',I12)") status CALL Pause() STOP END IF !---------------------------------------------------------------------------------------------------------------------------------- scanning: DO iBin = 1, 50 ! index of high-side of bin; index of low-side will be 0:49 IF ((random_R8(1) >= G_list(iBin)).AND.(random_R8(1) <= G_list(iBin-1))) THEN ! Note that G_list is a descending function of magnitude !This is the right bin, which contains the G-value fraction = (random_R8(1) - G_list(iBin)) / (G_list(iBin-1) - G_list(iBin)) ! from right to left on magnitude axis magnitude = magnitude_list(iBin) - fraction * d_magnitude EXIT scanning END IF END DO scanning WRITE (c5, "(F5.2)") magnitude build_EQC_line(54:58) = c5 EQC_line_store(lines) = build_EQC_line END DO ! i = 1, n (EQs in this cell) END IF ! 1 or more EQs in this cell, in this simulation END DO ! col = 1, nGRDcols END DO ! row = 1, nGRDrows !Sort lines of the EQC file: DO i = 1, (N_hat - 1) DO j = (i + 1), N_hat IF (s_within_window(i) > s_within_window(j)) THEN ! Swap them! spare_EQC_line = EQC_line_store(i) EQC_line_store(i) = EQC_line_store(j) EQC_line_store(j) = spare_EQC_line spare_R8 = s_within_window(i) s_within_window(i) = s_within_window(j) s_within_window(j) = spare_R8 END IF END DO END DO !Write out lines of the EQC file DO i = 1, N_hat WRITE (2, "(A)") EQC_line_store(i)(1:58) END DO DEALLOCATE ( s_within_window ) DEALLOCATE ( EQC_line_store ) CLOSE (2) END DO ! k = 1, m (outputting multiple files?) !Print summary statistics: mean = sum / m ! REAL*8 = REAL*8 / INTEGER (not integer-division, which truncates) variance = 0.0D0 ! just initializing... DO k = 1, m variance = variance + (N_hat_vector(k) - mean)**2 END DO IF (m > 1) THEN sigma = DSQRT(variance / (m - 1)) WRITE (*, "(' Mean = ',F10.2,' and sigma = ',F10.2)") mean, sigma ELSE WRITE (*, "(' Earthquake count = ',F10.2)") mean END IF WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause CONTAINS REAL*8 FUNCTION Area_of_cell(d_lon, d_lat, center_lat, R_m) !Returns area (in m**2) of a grid cell with width d_lon (degrees) !and height of d_lat (degrees), with central-point latitude of center_lat (degrees), !on a planet of radius R_m (meters). !Uses the exact (calculus) solution, not the approximate ~d_lat*d_lan*COS(center_lat)... USE DSphere ! provided by Peter Bird in file DSphere.f90 IMPLICIT NONE REAL*8, INTENT(IN) :: d_lon, d_lat, center_lat, R_m REAL*8 :: theta1, theta2 theta1 = (90.0D0 - (center_lat + d_lat/2.0D0)) * radians_per_degree ! integration limits, measured theta2 = (90.0D0 - (center_lat - d_lat/2.0D0)) * radians_per_degree ! from N pole, in radians theta1 = MAX(0.0D0, MIN(Pi, theta1)) ! just for good luck theta2 = MAX(0.0D0, MIN(Pi, theta2)) Area_of_cell = Two_Pi * (R_m**2) * (d_lon / 360.0D0) * (COS(theta1) - COS(theta2)) ! sic; !last term is equal to [ -COS(theta2) - (-COS(theta1) ]. END FUNCTION Area_of_cell DOUBLE PRECISION FUNCTION Moment(magnitude) !Per equation (1) of Bird & Kagan [2004, BSSA] DOUBLE PRECISION, INTENT(IN) :: magnitude DOUBLE PRECISION :: argument argument = 1.5D0 * magnitude + 9.05D0 Moment = 10.0D0 ** argument END FUNCTION Moment SUBROUTINE Pause () IMPLICIT NONE CHARACTER*1 c1 WRITE (*,"(' Press [Enter] to continue...'\)") READ (*,"(A)") c1 END SUBROUTINE Pause END PROGRAM SeismicityGRD_2_EQC