PROGRAM pseudoCSEP ! A collection of tests of gridded, stationary forecast(s) of ! earthquake rates, which are based on tests performed (first) ! in the RELM = Regional Earthquake Likelihood Models experiment: ! *Field, E. H. [2007] Overview of the Working Group for the ! Development of Regional Earthquake Likelihood Models ! (RELM), in S. E. Hough & K. B. Olsen (editors), ! Special Issue on: Regional Earthquake Likelihood Models, ! Seismol. Res. Lett., 78(1), 7-16. ! *Schorlemmer, D., M. C. Gerstenberger, S. Wiemer, D. D. Jackson, ! and D. A. Rhoades [2007] Earthquake likelihood model ! testing, Seismol. Res. Lett., 78(1), 17-29. ! *Schorlemmer, D., and M. C. Gerstenberger [2007] ! RELM Testing Center, Seismol. Res. Lett., 78(1), 30-36. ! ! and (second) in the CSEP = Collaboratory for the Study of ! Earthquake Predictability experiment: ! *Schorlemmer, D., J. D. Zechar, M. J. Werner, E. H. Field, ! D. D. Jackson, T. H. Jordan, and the RELM Working Group ! [2010] First results of the Regional Earthquake Likelihood ! Models experiment, Pure Appl. Geoph., 167, 859-876. ! *Zechar, J. D. [2008] Methods of evaluating earthquake predictions, ! Ph.D. dissertation, U. Southern Cal., Los Angeles, 148 pages. ! *Zechar, J. D., and T. H. Jordan [2008] Testing alarm-based earthquake ! predictions, Geophys. J. Int., 172, 715-724. ! *Zechar, J. D., M. C. Gerstenberger, and D. A. Rhoades [2010] ! Likelihood-based tests for evaluating space-rate-magnitude ! earthquake s, Bull. Seismol. Soc. Am., 100(3), 1184-1195. ! ! These tests are only a subset of the RELM/CSEP tests because ! they only apply to forecasts which are stationary during the ! time window of the forecast. Also, I use only one magnitude ! bin, which is "everything above the threshold magnitude." ! (Therefore, no M-test appears in this program.) ! The spatial area is assumed to be rectangular in (lon, lat) space ! (but more complex shapes can be handled by masking). ! Furthermore, these tests do not decluster the catalog or attempt ! to distinguish between mainshocks and aftershocks; ! therefore they should only be used on forecasts of all earthquakes. ! As in the RELM and CSEP tests, the depth range covers "shallow" ! seismicity, but the user is free to define what this means. ! By Peter Bird, UCLA, 2010.10 & 2013.05, based on above references. ! Option to mask part of the forecast, and options for various ! different kinds of ASS tests, added 2015.02.27. ! Minor revisions for better numerical precision 2018.01.08. ! INPUT consists of (at least) 2 files: ! ------------------------------------- ! *A seismic catalog for (exactly) the test time-window, which is ! either global or spatially larger than the forecast region; ! which extends at least as deep as the depth limit of forecast(s); ! which is complete at least down to the threshold magnitude of the ! forecast(s); and which also includes smaller events (whose ! magnitudes may be boosted in some modifications). ! The catalog must in Peter Bird's .eqc = EarthQuake Catalog ! format (as output by program Seismicity) which is documented in: ! Bird, P., and Y. Y. Kagan [2004] Plate-tectonic analysis ! of shallow seismicity: Apparent boundary width, beta, ! corner magnitude, coupled lithosphere thickness, and ! coupling in seven tectonic settings, ! Bull. Seismol. Soc. Am., 94(6), 2380-2399, ! plus electronic supplement (eqc_format.pdf). ! NOTE that the .EQC catalog file should be trimmed to the ! exact time window of the intended tests. ! NOTE that the optional focal-mechanism information at the ! end of each line in the .eqc file is not currently used. ! ! *A stationary forecast of seismicity rate with a single magnitude ! bin (everything exceeding threshold) presented using gridded ! longitude/latitude cells, according to the Long_Term_Seismicity ! variant of the .grd (Gridded Data) format documented at: ! http://peterbird.name/guide/grd_format.htm ! (NOTE that the Long_Term_Seismicity variant of the .grd format ! inserts the threshold magnitude into the first line, and also ! adds 2 lines with the forecast integral after the end of the ! long number list.) ! ! *Additional forecasts may be included, allowing for the R-test. ! ! *Optionally, a mask file may be read to block scoring in some ! cells of the forecast (e.g., if the test catalog is incomplete ! there); the same mask is used for all competing forecasts. ! OUTPUT is to the screen, and to a log file on Fortran unit 3. ! ----------------------- !========================= 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 !The following parameters directly control execution time; higher values mean more time; and higher m gives better precision. INTEGER, PARAMETER :: m = 1000 ! m = number of randomly-generated "simulated" earthquake catalogs for each forecast. INTEGER, PARAMETER :: s = 1000 ! s = number of sets of minor "modifications" of the actual earthquake catalog (if allowed). !In my opinion, this is terrible notation, and the symbols s & m should have been swapped; it is ! far too easy to jump to the conclusion that "m" goes with "modifications" and "s" goes with "simulations" ! whereas actually the reverse pairing is correct! !Note that Schorlemmer et al. [2007, SRL] used 10000 for both m and s (if I understand their paper correctly). !I have often used less because of both memory and execution-speed limits on my desktop computer. !Also, I doubt whether modifications (which require large s) add ANY value to the scoring; more likely, they degrade it! ! So, I am not willing to burn a lot of computer time on this. !CAUTION: Before increasing m > 1000, consider effects of m on storage as well as on execution time. ! For example, INTEGER*2, DIMENSION(m, nGRDrows, nGRDcols) :: sim_N_hat_3D_matrix ! already occupies 497 MB when m = 1000 and the forecasts use a global grid with 0.5-degree resolution: ! and it occupies 12.13 GB when m = 1000 and the forecasts use a global grid with 0.1-degree resolution. ! (Because sim_N_hat is re-used for each forecast, the number of forecasts being scored is irrelevant.) ! The storage of these integers is quite important for execution efficiency, because ! --to avoid recreating the Poisson cumulative distribution of each bin m times over-- ! we have to "sample" all the results (k = 1, ..., m) for that bin at once; ! on the other hand, we have to score each simulation (k = 1, ..., m) across all bins for each simulation. INTEGER, PARAMETER :: Factorial_limit = 170 ! Length of REAL*8 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 pseudoCSEP 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! INTEGER, PARAMETER :: Poisson_limit = 10000 ! Upper limit on integer argument to Poisson distribution, when using the normal distribution ! as an approximation to the Poisson distribution for large arguments. ! Must be larger than comb_limit above. ! There is no real upper limit on this PARAMETER, but really high values waste execution time & memory. ! NOTE: If your usage of pseudoCSEP is encountering problems with the limits above, try raising your magnitude threshold ! or shortening your test-time-window duration, as either of these will reduce the numbers of earthquakes considered. !---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 !----------------------------------------------------------------------------------------- CHARACTER*1 :: c1 CHARACTER*2 :: c2a, c2b, c2c, c2d, c2e CHARACTER*9 :: comment CHARACTER*13 :: c13 CHARACTER*127 :: mask_GRD_pathfile, Molchan_pathfile, & & new_EQC_pathfile, new_EQC_line_1, new_EQC_last_line, & & old_EQC_pathfile, old_EQC_line_1, old_EQC_last_line, & & pathfile, pseudoCSEP_logfile, & & reference_forecast_pathfile, & & text_line CHARACTER*1, DIMENSION(:), ALLOCATABLE :: new_eq_tenths CHARACTER*2, DIMENSION(:), ALLOCATABLE :: new_eq_month, new_eq_day, new_eq_hour, new_eq_minute, new_eq_second CHARACTER*127, DIMENSION(:), ALLOCATABLE :: GRD_pathfile INTEGER :: all_cells, ASS_type_index, col, denominator, i, i_mask, ii, iEq, ios, iYear, j, j_mask, jj, k, & & mask_nRows, mask_nColumns, masked_cells, maxDepthKm, & & N, n_done, N_hat, N_RI, N_true, ncStem, new_EQC_lines, nForecasts, nGRDcols, nGRDrows, numerator, & & old_EQC_lines, omega, q, row, s_in_use, unmasked_cells INTEGER, DIMENSION(:), ALLOCATABLE :: N_tilde ! a set of alternative EQ counts based on s modifications of the catalog INTEGER, DIMENSION(:), ALLOCATABLE :: mod_depth_int, new_eq_year, new_eq_depth_int, old_eq_depth_int INTEGER, DIMENSION(:,:), ALLOCATABLE :: actual_n_in_bin ! actual EQ count, per spatial bin INTEGER, DIMENSION(:,:), ALLOCATABLE :: mod_n_in_bin ! EQ count, per spatial bin, under one particular random modification of catalog INTEGER, DIMENSION(:,:), ALLOCATABLE :: old_n_in_bin ! EQ count during pre-test relative-intensity calibration period INTEGER, DIMENSION(:,:), ALLOCATABLE :: N_hat_matrix ! Number of simulated earthquakes, with scripts (k, j): ! 1 <= k <= m is the subscript identifying which simulation of this forecast; ! 1 <= j <= nForecasts is the superscript identifying the forecast. INTEGER, DIMENSION(:,:), ALLOCATABLE :: scaled_N_hat_matrix ! Number of simulated earthquakes, with scripts (k, j): ! 1 <= k <= m is the subscript identifying which simulation of this SCALED forecast; ! 1 <= j <= nForecasts is the superscript identifying the SCALED forecast. INTEGER*2, DIMENSION(:,:,:), ALLOCATABLE :: sim_N_hat_3D_matrix ! BIG working array, used temporarily during simulation process INTEGER*1, DIMENSION(:,:), ALLOCATABLE :: mask_grd_I1 ! (when using_mask = .TRUE.) LOGICAL :: any_rejection, exact_Poisson, global_mask, & & isEQCglobal, isEQCframed, isEQCdeeper, isEQCcomplete, isEQCmagBuffered, isTimeWindowAligned, isTimeWindowPrevious, & & not_warned_yet, randomize_test_EQs, tau_0p2_expectation_defined, using_mask LOGICAL, DIMENSION(:), ALLOCATABLE :: selected, old_selected ! will mark EQs which should be counted. LOGICAL*1, DIMENSION(:,:), ALLOCATABLE :: cell_in_use ! same size and shape as forecast grid(s); will be re-sampled from mask grid REAL :: amplification, & & box_lat_max, box_lat_min, box_lon_max, box_lon_min, & & contour_d_P, d_lat, d_lon, depth_sigma_km, & & EW_or_NS_sigma_km, expectation, & & GBs_requested, grid_lat_max, grid_lat_min, grid_lon_max, grid_lon_min, & & L, lat, lat_factor, latitude, lon, & & mag_offset, mag_sigma, & & mask_grid_lon_min, mask_d_lon, mask_grid_lon_max, & & mask_grid_lat_min, mask_d_lat, mask_grid_lat_max, mean, & & N_tilde_mean, N_tilde_sigma, new_EQC_duration_years, nu, & & sigma, sigmas, standard_deviation, & & t_d_lat, t_d_lon, t_lat_max, t_lat_min, t_lon_max, t_lon_min, temp, test_lon, threshold !--------------------------------------------------------------------------------------- !Dummy vectors (of 1 entry each) are needed to sync with MKL subprogramsL REAL, DIMENSION(1) :: VSCdfNorm_in, VSCdfNorm_out, VSCdfNormInv_in, VSCdfNormInv_out !--------------------------------------------------------------------------------------- REAL, DIMENSION(m) :: sim_temporary ! holds probabilities from which EQ #s will be deduced REAL, DIMENSION(:), ALLOCATABLE :: real_vector_s ! holds results of a set of catalog modifications REAL, DIMENSION(:), ALLOCATABLE :: ASS_score, ASS_quantile, & & L_sup_j, L_tilde_mean, L_tilde_sigma, & & mod_temporary, mod_Elon, mod_Nlat, mod_mag, & & new_eq_Elon, new_eq_Nlat, new_eq_mag, N_sup_j, & & old_eq_Elon, old_eq_Nlat, old_eq_mag, & & scaled_L_sup_j, zeta REAL, DIMENSION(:,:), ALLOCATABLE :: del_sub_q_sup_j, gamma_sub_q_sup_j ! both are quantile scores for particular ! simulated catalogs (#q) based on particular ! forecasts (#j). REAL, DIMENSION(:,:), ALLOCATABLE :: d_P ! Based on the current forecast, what is the conditional probability that ! the next qualifying earthquake (in the space/depth/magnitude window of the test) ! will occur in this bin (#i = 1, nGRDrows, #j = 1, nGRDcols)? N.B. Sum of entries should be 1. REAL, DIMENSION(:,:), ALLOCATABLE :: d_tau ! Based on pre-test seismicity, what is the conditional probability that ! the next qualifying earthquake (in the space/depth/magnitude window of the test) ! will occur in this bin (#i = 1, nGRDrows, #j = 1, nGRDcols)? N.B. Sum of entries should be 1. ! NOTE: For the ASS-test. d_tau is the prior null hypothesis; d_P is the test hypothesis. REAL, DIMENSION(:,:), ALLOCATABLE :: alpha_matrix ! (nForecasts by nForecasts) REAL, DIMENSION(:,:), ALLOCATABLE :: R_square REAL, DIMENSION(:,:), ALLOCATABLE :: PoissonCDF ! (0:Poisson_limit, j = 1:nForecasts) REAL, DIMENSION(:,:), ALLOCATABLE :: dt_da_factor ! converts from EQ rates in #/m^2/s to EQ #'s in test time window. REAL, DIMENSION(:,:), ALLOCATABLE :: L_tilde_matrix ! log-likelihoods, with scripts (q, j): ! 1 <= q <= s is the subscript identifying which modification of the catalog; ! 1 <= j <= nForecasts is the superscript identifying the forecast. REAL, DIMENSION(:,:), ALLOCATABLE :: Molchan ! Molchan curve for one particular forecast, with 2 dimensionless coordinates ! (1st subscript 1 = nu; 1st subscript 2 = tau) ! (2nd subscript 0 for upper-left corner, or 1-N for number of actual earthquakes) REAL, DIMENSION(:,:), ALLOCATABLE :: scaled_L_hat_matrix ! (k = 1, ..., m is the simulation #; j = 1, nForecasts is the forecast #) REAL, DIMENSION(:,:,:), ALLOCATABLE :: L_hat_matrix ! log-likelihoods, with scripts (k, j, jj): ! 1 <= k <= m is the subscript identifying which simulation of this forecast; ! 1 <= j <= nForecasts is the superscript identifying the forecast ! used to create the simulation. ! 1 <= jj <= nForecasts is the forecast scored against this simulation. REAL, DIMENSION(:,:,:), ALLOCATABLE :: R_tilde_3D_matrix ! subscript q = 1, ..., s for modification of actual catalog; ! two superscripts for R^ij == L^i - L^j (i, j = 1, ..., nForecasts). DOUBLE PRECISION :: argument, factor1, factor2, factor3, L_sum, lambda, lost_probability, N_per_year, & & omega_R8, root_lambda, sum, tau, temp8, total_probability, variance !--------------------------------------------------------------------------------------- !Dummy vectors (of 1 entry each) are needed to sync with MKL subprogramsL DOUBLE PRECISION, DIMENSION(1) :: VDCdfNorm_in, VDCdfNorm_out !--------------------------------------------------------------------------------------- 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 :: reference_forecast ! IF (ASS_type_index == 3) DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: forecast, scaled_forecast ! DOUBLE PRECISION because of small exponents, like E-22! !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 ----------------------------------------------------------- ! 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 (',I3,') exceeds Factorial_limit (',I3').')") comb_limit, Factorial_limit WRITE (3, "(' ERROR: PARAMETER comb_limit (',I3,') exceeds Factorial_limit (',I3').')") comb_limit, Factorial_limit CALL Pause() STOP END IF IF (comb_limit > 140) THEN WRITE (*, *) WRITE (*, "(' WARNING: PARAMETER comb_limit (',I3,') > 140 may lead to floating-overflow abends!')") comb_limit WRITE (*, *) WRITE (3, *) WRITE (3, "(' WARNING: PARAMETER comb_limit (',I3,') > 140 may lead to floating-overflow abends!')") comb_limit WRITE (3, *) CALL Pause() END IF IF (Poisson_limit <= comb_limit) THEN WRITE (*, "(' ERROR: PARAMETER Poisson_limit (',I6,') must exceed comb_limit (',I3').')") Poisson_limit, comb_limit WRITE (3, "(' ERROR: PARAMETER Poisson_limit (',I6,') must exceed comb_limit (',I3').')") Poisson_limit, comb_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)") 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 LOG_Factorial(i) = LOG_Factorial(i-1) + LOG(1.0D0 * i) !WRITE (*, *) i, Factorial(i), LOG_Factorial(i) ! <=== un-comment to see progression toward overflow END DO !------------------------------------------------------------------------------------------- WRITE (*, "(' --== PROGRAM pseudoCSEP ==--')") WRITE (*, *) WRITE (*, "(' A collection of tests of gridded, stationary forecast(s)')") WRITE (*, "(' of earthquake rates, which are based on the tests performed')") WRITE (*, "(' (first) in the RELM = Regional Earthquake Likelihood Models')") WRITE (*, "(' experiment:')") WRITE (*, "(' *Field, E. H. [2007] Overview of the Working Group for the')") WRITE (*, "(' Development of Regional Earthquake Likelihood Models')") WRITE (*, "(' (RELM), in S. E. Hough & K. B. Olsen (editors),')") WRITE (*, "(' Special Issue on: Regional Earthquake Likelihood Models,')") WRITE (*, "(' Seismol. Res. Lett., 78(1), 7-16.')") WRITE (*, "(' *Schorlemmer, D., M. C. Gerstenberger, S. Wiemer, D. D. Jackson,')") WRITE (*, "(' and D. A. Rhoades [2007] Earthquake likelihood model')") WRITE (*, "(' testing, Seismol. Res. Lett., 78(1), 17-29.')") WRITE (*, "(' *Schorlemmer, D., and M. C. Gerstenberger [2007] ')") WRITE (*, "(' RELM Testing Center, Seismol. Res. Lett., 78(1), 30-36.')") WRITE (*, *) CALL Pause() WRITE (*, *) WRITE (*, "(' and (second) in the CSEP = Collaboratory for the Study of Earthquake')") WRITE (*, "(' Predictability experiment:')") WRITE (*, "(' *Schorlemmer, D., J. D. Zechar, M. J. Werner, E. H. Field,')") WRITE (*, "(' D. D. Jackson, T. H. Jordan, and the RELM Working Group')") WRITE (*, "(' [2010] First results of the Regional Earthquake Likelihood')") WRITE (*, "(' Models experiment, Pure Appl. Geoph., 167, 859-876.')") WRITE (*, "(' *Zechar, J. D. [2008] Methods of evaluating earthquake predictions,')") WRITE (*, "(' Ph.D. dissertation, U. Southern Cal., Los Angeles, 148 pages.')") WRITE (*, "(' *Zechar, J. D., and T. H. Jordan [2008] Testing alarm-based earthquake')") WRITE (*, "(' predictions, Geophys. J. Int., 172, 715-724.')") WRITE (*, "(' *Zechar, J. D., M. C. Gerstenberger, and D. A. Rhoades [2010]')") WRITE (*, "(' Likelihood-based tests for evaluating space-rate-magnitude')") WRITE (*, "(' earthquake forecasts, Bull. Seismol. Soc. Am., 100(3), 1184-1195.')") WRITE (*, *) CALL Pause() WRITE (*, *) WRITE (*, "(' These tests are only a subset of the RELM/CSEP tests because')") WRITE (*, "(' they only apply to forecasts which are stationary during the')") WRITE (*, "(' time window of the forecast. Also, I use only one magnitude')") WRITE (*, "(' bin, which is ""everything above the threshold magnitude.""')") WRITE (*, "(' (Therefore, no M-test appears in this program.)')") WRITE (*, "(' The spatial area is assumed to be rectangular in (lon, lat) space')") WRITE (*, "(' (but more complex shapes can be handled by masking).')") WRITE (*, "(' Furthermore, these tests do not decluster the catalog or attempt')") WRITE (*, "(' to distinguish between mainshocks and aftershocks;')") WRITE (*, "(' therefore they should only be used on forecasts of all earthquakes.')") WRITE (*, "(' As in the RELM and CSEP tests, the depth range covers ""shallow""')") WRITE (*, "(' seismicity, but the user is free to define what this means.')") WRITE (*, *) WRITE (*, "(' By Peter Bird, UCLA, 2010.10 & 2013.05, based on above references.')") WRITE (*, "(' Option to mask part of the forecast, and options for various')") WRITE (*, "(' different kinds of ASS tests, added 2015.02.27.')") WRITE (*, "(' Minor revisions for better numerical precision 2018.01.08.')") WRITE (*, *) CALL Pause() !======================================== begin INPUTS ===================================================== WRITE (*, *) pseudoCSEP_logfile = "pseudoCSEP_log.txt" CALL Prompt_for_String("Filename for new pseudoCSEP log file?", pseudoCSEP_logfile, pseudoCSEP_logfile) OPEN (UNIT = 3, FILE = TRIM(pseudoCSEP_logfile)) ! unconditional OPEN; will overwrite any existing file of same name! !N.B. By choosing different, customized logfile names, you can run multiple copies of pseudoCSEP at the same time. WRITE (3, "(' --== PROGRAM pseudoCSEP ==--')") WRITE (3, "(' ')") WRITE (3, "(' A collection of tests of gridded, stationary forecast(s)')") WRITE (3, "(' of earthquake rates, which are based on the tests performed')") WRITE (3, "(' (first) in the RELM = Regional Earthquake Likelihood Models')") WRITE (3, "(' experiment:')") WRITE (3, "(' *Field, E. H. [2007] Overview of the Working Group for the')") WRITE (3, "(' Development of Regional Earthquake Likelihood Models')") WRITE (3, "(' (RELM), in S. E. Hough & K. B. Olsen (editors),')") WRITE (3, "(' Special Issue on: Regional Earthquake Likelihood Models,')") WRITE (3, "(' Seismol. Res. Lett., 78(1), 7-16.')") WRITE (3, "(' *Schorlemmer, D., M. C. Gerstenberger, S. Wiemer, D. D. Jackson,')") WRITE (3, "(' and D. A. Rhoades [2007] Earthquake likelihood model')") WRITE (3, "(' testing, Seismol. Res. Lett., 78(1), 17-29.')") WRITE (3, "(' *Schorlemmer, D., and M. C. Gerstenberger [2007] ')") WRITE (3, "(' RELM Testing Center, Seismol. Res. Lett., 78(1), 30-36.')") WRITE (3, "(' and (second) in the CSEP = Collaboratory for the Study of Earthquake')") WRITE (3, "(' Predictability experiment:')") WRITE (3, "(' *Schorlemmer, D., J. D. Zechar, M. J. Werner, E. H. Field,')") WRITE (3, "(' D. D. Jackson, T. H. Jordan, and the RELM Working Group')") WRITE (3, "(' [2010] First results of the Regional Earthquake Likelihood')") WRITE (3, "(' Models experiment, Pure Appl. Geoph., 167, 859-876.')") WRITE (3, "(' *Zechar, J. D. [2008] Methods of evaluating earthquake predictions,')") WRITE (3, "(' Ph.D. dissertation, U. Southern Cal., Los Angeles, 148 pages.')") WRITE (3, "(' *Zechar, J. D., and T. H. Jordan [2008] Testing alarm-based earthquake')") WRITE (3, "(' predictions, Geophys. J. Int., 172, 715-724.')") WRITE (3, "(' *Zechar, J. D., M. C. Gerstenberger, and D. A. Rhoades [2010]')") WRITE (3, "(' Likelihood-based tests for evaluating space-rate-magnitude')") WRITE (3, "(' earthquake forecasts, Bull. Seismol. Soc. Am., 100(3), 1184-1195.')") WRITE (3, "(' These tests are only a subset of the RELM/CSEP tests because')") WRITE (3, "(' they only apply to forecasts which are stationary during the')") WRITE (3, "(' time window of the forecast. Also, I use only one magnitude')") WRITE (3, "(' bin, which is ""everything above the threshold magnitude.""')") WRITE (3, "(' (Therefore, no M-test appears in this program.)')") WRITE (3, "(' The spatial area is assumed to be rectangular in (lon, lat) space')") WRITE (3, "(' (but more complex shapes can be handled by masking).')") WRITE (3, "(' Furthermore, these tests do not decluster the catalog or attempt')") WRITE (3, "(' to distinguish between mainshocks and aftershocks;')") WRITE (3, "(' therefore they should only be used on forecasts of all earthquakes.')") WRITE (3, "(' As in the RELM and CSEP tests, the depth range covers ""shallow""')") WRITE (3, "(' seismicity, but the user is free to define what this means.')") WRITE (3, *) WRITE (3, "(' By Peter Bird, UCLA, 2010.10, based on above references.')") WRITE (3, "(' Option to mask part of the forecast, and options for various')") WRITE (3, "(' different kinds of ASS tests, added 2015.02.25.')") WRITE (3, "(' Minor revisions for better numerical precision 2018.01.08.')") 110 WRITE (*, *) WRITE (3, *) WRITE (*, "(' REQUIRED INPUT: EARTHQUAKE CATALOG for the test time-window.')") WRITE (*, "(' *A seismic catalog for (exactly) the test time-window, which is')") WRITE (*, "(' either global or spatially larger than the forecast region;')") WRITE (*, "(' which extends at least as deep as the depth limit of forecast(s);')") WRITE (*, "(' which is complete down to the threshold magnitude of the')") WRITE (*, "(' forecast(s); and which also includes smaller events (whose')") WRITE (*, "(' magnitudes may be boosted in some modifications).')") WRITE (*, "(' The catalog must be in Peter Bird''s .eqc = EarthQuake Catalog')") WRITE (*, "(' format (as output by program Seismicity) which is documented in:')") WRITE (*, "(' Bird, P., and Y. Y. Kagan [2004] Plate-tectonic analysis')") WRITE (*, "(' of shallow seismicity: Apparent boundary width, beta,')") WRITE (*, "(' corner magnitude, coupled lithosphere thickness, and')") WRITE (*, "(' coupling in seven tectonic settings,')") WRITE (*, "(' Bull. Seismol. Soc. Am., 94(6), 2380-2399,')") WRITE (*, "(' plus electronic supplement (eqc_format.pdf).')") WRITE (*, "(' NOTE that the .EQC catalog file should be trimmed to the')") WRITE (*, "(' exact time window of the intended tests.')") WRITE (*, "(' NOTE that the optional focal-mechanism information at the')") WRITE (*, "(' end of each line in the .eqc file is not currently used.')") WRITE (*, *) WRITE (3, "(' REQUIRED INPUT: EARTHQUAKE CATALOG for the test time-window.')") WRITE (3, "(' *A seismic catalog for (exactly) the test time-window, which is')") WRITE (3, "(' either global or spatially larger than the forecast region;')") WRITE (3, "(' which extends at least as deep as the depth limit of forecast(s);')") WRITE (3, "(' which is complete down to the threshold magnitude of the')") WRITE (3, "(' forecast(s); and which also includes smaller events (whose')") WRITE (3, "(' magnitudes may be boosted in some modifications).')") WRITE (3, "(' The catalog must be in Peter Bird''s .eqc = EarthQuake Catalog')") WRITE (3, "(' format (as output by program Seismicity) which is documented in:')") WRITE (3, "(' Bird, P., and Y. Y. Kagan [2004] Plate-tectonic analysis')") WRITE (3, "(' of shallow seismicity: Apparent boundary width, beta,')") WRITE (3, "(' corner magnitude, coupled lithosphere thickness, and')") WRITE (3, "(' coupling in seven tectonic settings,')") WRITE (3, "(' Bull. Seismol. Soc. Am., 94(6), 2380-2399,')") WRITE (3, "(' plus electronic supplement (eqc_format.pdf).')") WRITE (3, "(' NOTE that the .EQC catalog file should be trimmed to the')") WRITE (3, "(' exact time window of the intended tests.')") WRITE (3, "(' NOTE that the optional focal-mechanism information at the')") WRITE (3, "(' end of each line in the .eqc file is not currently used.')") WRITE (3, *) CALL Prompt_for_String ('[Drive:][\path\]filename of the new .EQC file:', '', new_EQC_pathfile) WRITE (3, "(' [Drive:][\path\]filename of the new .EQC file:')") WRITE (3, "(' ',A)") TRIM(new_EQC_pathfile) OPEN (UNIT = 1, FILE = TRIM(new_EQC_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: ',A,' could not be found (in current folder).')") TRIM(new_EQC_pathfile) WRITE (*,"(' Please check [Drive:][\path\]filename and/or move file, and try again...')") CALL Pause() GO TO 110 END IF !Count lines which include a valid year number: new_EQC_lines = 0 not_warned_yet = .TRUE. count_lines: DO READ (1, "(9X,I5)", IOSTAT = ios) iYear IF (ios == -1) EXIT count_lines ! EOF new_EQC_lines = new_EQC_lines + 1 IF (((iYear < 1800).OR.(iYear > 2099)).AND.not_warned_yet) THEN WRITE (*, *) WRITE (*, "(' WARNING: After reading ',I7,' lines, met bad(?) year number: ',I5)") new_EQC_lines, iYear WRITE (*, "(' This warning is non-fatal, and execution will continue...')") WRITE (3, *) WRITE (3, "(' WARNING: After reading ',I7,' lines, met bad(?) year number: ',I5)") new_EQC_lines, iYear WRITE (3, "(' This warning is non-fatal, and execution will continue...')") CALL Pause() WRITE (*, *) WRITE (3, *) not_warned_yet = .FALSE. END IF BACKSPACE(1) READ (1, "(A)") text_line IF (new_EQC_lines == 1) new_EQC_line_1 = TRIM(text_line) new_EQC_last_line = TRIM(text_line) END DO count_lines CLOSE(1) WRITE (*, "(' This catalog file includes ',I10,' events.')") new_EQC_lines WRITE (3, "(' This catalog file includes ',I10,' events.')") new_EQC_lines WRITE (*, *) WRITE (*, *) WRITE (*, "(' pseudoCSEP includes an option to add Gaussian random noise to ')") WRITE (*, "(' magnitudes and/or epicenters/epicentroids and/or depths of quakes')") WRITE (*, "(' in this test catalog.')") WRITE (*, "(' Such randomization was used in the RELM test, but is now less common')") WRITE (*, "(' in routine CSEP testing; it is reserved for sensitivity studies.')") WRITE (3, *) WRITE (3, "(' pseudoCSEP includes an option to add Gaussian random noise to ')") WRITE (3, "(' magnitudes and/or epicenters/epicentroids and/or depths of quakes')") WRITE (3, "(' in this test catalog.')") WRITE (3, "(' Such randomization was used in the RELM test, but is now less common')") WRITE (3, "(' in routine CSEP testing; it is reserved for sensitivity studies.')") CALL Prompt_for_Logical('Do you want to randomize test earthquake data?',.FALSE.,randomize_test_EQs) WRITE (3, "(' Do you want to randomize test earthquake data?')") WRITE (3, "(' ',L5)") randomize_test_EQs IF (randomize_test_EQs) THEN s_in_use = s ! INTEGER, PARAMETER, set at the top of this code ELSE s_in_use = 1 END IF IF (.NOT.ALLOCATED(N_tilde)) ALLOCATE ( N_tilde(s_in_use) ) ! (testing, because we may loop back here from below) IF (.NOT.ALLOCATED(real_vector_s)) ALLOCATE ( real_vector_s(s_in_use) ) WRITE (*, *) WRITE (3, *) CALL Prompt_for_Logical ('Is this catalog global in its areal coverage?', .FALSE., isEQCglobal) WRITE (3, "(' Is this catalog global in its areal coverage?')") WRITE (3, "(L5)") isEQCglobal IF (isEQCglobal) THEN isEQCframed = .TRUE. ELSE CALL Prompt_for_Logical ('Does this catalog include a lateral buffer zone (frame)?', .FALSE., isEQCframed) WRITE (3, "(' Does this catalog include a lateral buffer zone (frame)?')") WRITE (3, "(L5)") isEQCframed IF ((.NOT.isEQCframed).AND.randomize_test_EQs) THEN WRITE (*, "(' Lateral framing permits some mislocated earthquakes to be recovered.')") WRITE (*, "(' Please expand the catalog area and try again.')") WRITE (*, "(' If this is not possible, then enter zero for epicenter location uncertainty (below).')") WRITE (3, "(' Lateral framing permits some mislocated earthquakes to be recovered.')") WRITE (3, "(' Please expand the catalog area and try again.')") WRITE (3, "(' If this is not possible, then enter zero for epicenter location uncertainty (below).')") CALL Pause() GO TO 110 END IF END IF CALL Prompt_for_Logical ('Does this catalog extend deeper than the forecast(s)?', .FALSE., isEQCdeeper) WRITE (3, "(' Does this catalog extend deeper than the forecast(s)?')") WRITE (3, "(L5)") isEQCdeeper IF ((.NOT.isEQCdeeper).AND.randomize_test_EQs) THEN WRITE (*, "(' Vertical framing permits some mislocated earthquakes to be recovered.')") WRITE (*, "(' Please expand the catalog depth limit and try again.')") WRITE (*, "(' If this is not possible, then enter zero for hypocenter depth uncertainty (below).')") WRITE (3, "(' Vertical framing permits some mislocated earthquakes to be recovered.')") WRITE (3, "(' Please expand the catalog depth limit and try again.')") WRITE (3, "(' If this is not possible, then enter zero for hypocenter depth uncertainty (below).')") CALL Pause() GO TO 110 END IF CALL Prompt_for_Logical ('Is this catalog complete above the threshold of the forecast(s)?', .FALSE., isEQCcomplete) WRITE (3, "(' Is this catalog complete above the threshold of the forecast(s)?')") WRITE (3, "(L5)") isEQCcomplete IF (.NOT.isEQCcomplete) THEN WRITE (*, "(' Please lower the catalog magnitude threshold, or change catalogs, and try again.')") WRITE (*, "(' It could be seriously misleading to test forecast(s) against an incomplete catalog.')") WRITE (3, "(' Please lower the catalog magnitude threshold, or change catalogs, and try again.')") WRITE (3, "(' It could be seriously misleading to test forecast(s) against an incomplete catalog.')") CALL Pause() GO TO 110 END IF CALL Prompt_for_Logical ('Does this catalog include (some) events below threshold of the forecast(s)?', .FALSE., isEQCmagBuffered) WRITE (3, "(' Does this catalog include (some) events below threshold of the forecast(s)?')") WRITE (3, "(L5)") isEQCmagBuffered IF ((.NOT.isEQCmagBuffered).AND.randomize_test_EQs) THEN WRITE (*, "(' Magnitude framing permits some mis-measured earthquakes to be recovered.')") WRITE (*, "(' Please lower the catalog magnitude limit and try again.')") WRITE (*, "(' If this is not possible, then enter zero for magnitude uncertainty (below).')") WRITE (3, "(' Magnitude framing permits some mis-measured earthquakes to be recovered.')") WRITE (3, "(' Please lower the catalog magnitude limit and try again.')") WRITE (3, "(' If this is not possible, then enter zero for magnitude uncertainty (below).')") CALL Pause() GO TO 110 END IF WRITE (*, "(' ')") WRITE (*, "(' Here are the first and last lines of the .EQC file:')") WRITE (3, "(' Here are the first and last lines of the .EQC file:')") WRITE (*, "(' ',A)") TRIM(new_EQC_line_1) WRITE (3, "(' ',A)") TRIM(new_EQC_line_1) WRITE (*, "(' ',A)") TRIM(new_EQC_last_line) WRITE (3, "(' ',A)") TRIM(new_EQC_last_line) CALL Prompt_for_Real ('What is the duration of this catalog, in decimal years?', 5.0, new_EQC_duration_years) WRITE (3, "(' What is the duration of this catalog, in decimal years?')") WRITE (3, "(F10.2)") new_EQC_duration_years CALL Prompt_for_Logical ('Is its time-window EXACTLY the same as that of the forecast?',.FALSE.,isTimeWindowAligned) WRITE (3, "(' Is its time-window EXACTLY the same as that of the forecast?')") WRITE (3, "(L5)") isTimeWindowAligned IF (.NOT.isTimeWindowAligned) THEN WRITE (*, "(' Please adjust the catalog time-window, or change catalogs, and try again.')") WRITE (*, "(' It would be seriously misleading to test forecast(s) against a catalog of incorrect length!')") WRITE (3, "(' Please adjust the catalog time-window, or change catalogs, and try again.')") WRITE (3, "(' It would be seriously misleading to test forecast(s) against a catalog of incorrect length!')") CALL Pause() GO TO 110 END IF IF (randomize_test_EQs) THEN !Quiz the user about the uncertainties in the catalog (to be used for catalog modifications below): WRITE (*, *) WRITE (3, *) WRITE (*, "(' Next, please describe the UNCERTAINTIES of numbers in this catalog:')") WRITE (3, "(' Next, please describe the UNCERTAINTIES of numbers in this catalog:')") WRITE (*, *) WRITE (*, "(' Magnitude uncertainty is dimensionless and typically less than 1.0.')") WRITE (*, "(' For example, results of Kagan [2003, PEPI] suggest that magnitudes in the')") WRITE (*, "(' Centroid Moment Tensor catalog have uncertainties of about 0.0884.')") WRITE (3, "(' Magnitude uncertainty is dimensionless and typically less than 1.0.')") WRITE (3, "(' For example, results of Kagan [2003, PEPI] suggest that magnitudes in the')") WRITE (3, "(' Centroid Moment Tensor catalog have uncertainties of about 0.0884.')") CALL Prompt_for_Real ('What is the uncertainty/sigma/S.D. of magnitudes in this catalog?', 0.0884, mag_sigma) WRITE (3, "(' What is the uncertainty/sigma/S.D. of magnitudes in this catalog?')") WRITE (3, "(F10.4)") mag_sigma WRITE (*, *) WRITE (3, *) WRITE (*, "(' Lateral/positional uncertainty of hypocenters is often greater than the nominal')") WRITE (*, "(' uncertainty quoted in the catalog, due to systematic errors in velocity model.')") WRITE (*, "(' Results of Kagan [2003, PEPI] suggest about 17.7 km for the CMT catalog.')") WRITE (3, "(' Lateral/positional uncertainty of hypocenters is often greater than the nominal')") WRITE (3, "(' uncertainty quoted in the catalog, due to systematic errors in velocity model.')") WRITE (3, "(' Results of Kagan [2003, PEPI] suggest about 17.7 km for the CMT catalog.')") CALL Prompt_for_Real ('What is the lateral (E-W, or N-S) uncertainty/sigma/S.D. of epicenters, in km?', 17.7, EW_or_NS_sigma_km) WRITE (3, "(' What is the lateral (E-W, or N-S) uncertainty/sigma/S.D. of epicenters, in km?')") WRITE (3, "(F10.3)") EW_or_NS_sigma_km WRITE (*, *) WRITE (3, *) WRITE (*, "(' Uncertainty in hypocental/hypocentroid depths is MUCH GREATER than nominal')") WRITE (*, "(' uncertainty quoted in the catalog, due to systematic errors in velocity model.')") WRITE (*, "(' Results of Bird & Liu [2007, SRL] suggest ~25 km for offshore regions in CMT.')") WRITE (3, "(' Uncertainty in hypocental/hypocentroid depths is MUCH GREATER than nominal')") WRITE (3, "(' uncertainty quoted in the catalog, due to systematic errors in velocity model.')") WRITE (3, "(' Results of Bird & Liu [2007, SRL] suggest ~25 km for offshore regions in CMT.')") CALL Prompt_for_Real ('What is the uncertainty/sigma/S.D. of hypocentral depths, in km?', 25., depth_sigma_km) WRITE (3, "(' What is the uncertainty/sigma/S.D. of hypocentral depths, in km?')") WRITE (3, "(F10.3)") depth_sigma_km ELSE mag_sigma = 0.0 EW_or_NS_sigma_km = 0.0 depth_sigma_km = 0.0 END IF !Allocate storage, and then re-read the catalog to memorize it. ALLOCATE ( new_eq_year(new_EQC_lines) ) ALLOCATE ( new_eq_month(new_EQC_lines) ) ALLOCATE ( new_eq_day(new_EQC_lines) ) ALLOCATE ( new_eq_hour(new_EQC_lines) ) ALLOCATE ( new_eq_minute(new_EQC_lines) ) ALLOCATE ( new_eq_second(new_EQC_lines) ) ALLOCATE ( new_eq_tenths(new_EQC_lines) ) ALLOCATE ( new_eq_Elon(new_EQC_lines) ) ALLOCATE ( new_eq_Nlat(new_EQC_lines) ) ALLOCATE ( new_eq_depth_int(new_EQC_lines) ) ALLOCATE ( new_eq_mag(new_EQC_lines) ) ALLOCATE ( mod_temporary(new_EQC_lines) ) ! temporary vector, used to modify the catalog ALLOCATE ( mod_Elon(new_EQC_lines) ) ALLOCATE ( mod_Nlat(new_EQC_lines) ) ALLOCATE ( mod_depth_int(new_EQC_lines) ) ALLOCATE ( mod_mag(new_EQC_lines) ) OPEN (UNIT = 1, FILE = TRIM(new_EQC_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) next_line: DO iEq = 1, new_EQC_lines READ (1, 120, IOSTAT = ios) & & new_eq_year(iEq), new_eq_month(iEq), new_eq_day(iEq), & & new_eq_hour(iEq), new_eq_minute(iEq), new_eq_second(iEq), new_eq_tenths(iEq), & & new_eq_Elon(iEq), new_eq_Nlat(iEq), & & new_eq_depth_int(iEq), new_eq_mag(iEq) !Note: Any appended data, such as focal mechanism, is ignored. 120 FORMAT (9X, & & I5,'.',A2,'.',A2, 1X, & ! using I5 for negative years (B.C.) & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) IF (ios /= 0) EXIT next_line END DO next_line ! iEq = 1, new_EQC_lines CLOSE(1) ! end reading of new (test-time-window) seismic catalog !----------------------------------------------------------------------------------------------- ! begin reading of forecast(s): WRITE (*, *) WRITE (3, *) WRITE (*, "(' REQUIRED INPUT: AT LEAST ONE FORECAST:')") WRITE (*, "(' *A stationary forecast of seismicity rate with a single magnitude')") WRITE (*, "(' bin (everything exceeding threshold) presented at gridded')") WRITE (*, "(' longitude/latitude points, according to the Long_Term_Seismicity')") WRITE (*, "(' variant of the .grd (Gridded Data) format documented at:')") WRITE (*, "(' http://peterbird.name/guide/grd_format.htm')") WRITE (*, "(' (NOTE that the Long_Term_Seismicity variant of the .grd format')") WRITE (*, "(' inserts the threshold magnitude into the first line, and also')") WRITE (*, "(' adds 2 lines with the forecast integral after the end of the')") WRITE (*, "(' long number list.)')") WRITE (3, "(' REQUIRED INPUT: AT LEAST ONE FORECAST:')") WRITE (3, "(' *A stationary forecast of seismicity rate with a single magnitude')") WRITE (3, "(' bin (everything exceeding threshold) presented in gridded')") WRITE (3, "(' longitude/latitude cells, according to the Long_Term_Seismicity')") WRITE (3, "(' variant of the .grd (Gridded Data) format documented at:')") WRITE (3, "(' http://peterbird.name/guide/grd_format.htm')") WRITE (3, "(' (NOTE that the Long_Term_Seismicity variant of the .grd format')") WRITE (3, "(' inserts the threshold magnitude into the first line, and also')") WRITE (3, "(' adds 2 lines with the forecast integral after the end of the')") WRITE (3, "(' long number list.)')") WRITE (*, *) WRITE (3, *) CALL Prompt_for_Integer('How many forecasts will be input for this run?',1,nForecasts) WRITE (3, "(' How many forecasts will be input for this run?')") WRITE (3, *) nForecasts !Note: Schorlemmer et al. [2007, SRL] identify models by j, but have no symbol corresponding to nForecasts. ALLOCATE ( GRD_pathfile(nForecasts) ) ALLOCATE ( L_tilde_mean(nForecasts) ) ALLOCATE ( L_tilde_sigma(nForecasts) ) ALLOCATE ( zeta(nForecasts) ) ! results of the S-test ALLOCATE ( PoissonCDF(0:Poisson_limit, nForecasts) ) ALLOCATE ( L_tilde_matrix(s_in_use, nForecasts) ) ! (q, j), where 1 <= q <= s, and 1 <= j <= nForecasts. 200 IF (nForecasts == 1) THEN CALL Prompt_for_String('Enter [Drive:][\path\]filename (.GRD) for this forecast:','',pathfile) WRITE (3, "(' Enter [Drive:][\path\]filename (.GRD) for this forecast:')") WRITE (3, "(' ',A)") TRIM(pathfile) ELSE WRITE (*, "(' CAUTION: All forecasts submitted in this one run must share the:')") WRITE (*, "(' *SAME forecast time-window;')") WRITE (*, "(' *SAME spatial grid, and area covered;')") WRITE (*, "(' *SAME depth limit; and')") WRITE (*, "(' *SAME threshold (minimum) magnitude!')") WRITE (3, "(' CAUTION: All forecasts submitted in this one run must share the:')") WRITE (3, "(' *SAME forecast time-window;')") WRITE (3, "(' *SAME spatial grid, and area covered;')") WRITE (3, "(' *SAME depth limit; and')") WRITE (3, "(' *SAME threshold (minimum) magnitude!')") CALL Pause() CALL Prompt_for_String('Enter [Drive:][\path\]filename (.GRD) for the first forecast:','',pathfile) WRITE (3, "(' Enter [Drive:][\path\]filename (.GRD) for the first forecast:')") WRITE (3, "(' ',A)") TRIM(pathfile) END IF GRD_pathfile(1) = pathfile OPEN (UNIT = 2, FILE = TRIM(GRD_pathfile(1)), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Forecast .GRD file ',A,' not found (in current folder).')") TRIM(GRD_pathfile(1)) CALL Pause() GO TO 200 END IF WRITE (*, *) WRITE (3, *) WRITE (*, "(' Here are the first 2 lines of this file:')") WRITE (3, "(' Here are the first 2 lines of this file:')") READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) WRITE (*, *) WRITE (3, *) IF (nForecasts == 1) THEN CALL Prompt_for_Real('What is the threshold magnitude of the forecast?',5.767,threshold) WRITE (3, "(' What is the threshold magnitude of the forecast?')") WRITE (3, *) threshold ELSE CALL Prompt_for_Real('What is the threshold magnitude of all the forecasts?',5.767,threshold) WRITE (3, "(' What is the threshold magnitude of all the forecasts?')") WRITE (3, *) threshold END IF WRITE (*, "(' ')") CALL Prompt_for_Integer('What (integer) gives the forecast depth limit, in km?',70,maxDepthKm) WRITE (3, "(' What (integer) gives the forecast depth limit, in km?')") WRITE (3, *) maxDepthKm CLOSE(2) OPEN (UNIT = 2, FILE = TRIM(GRD_pathfile(1)), STATUS = "OLD", IOSTAT = ios) READ (2, *) grid_lon_min, d_lon, grid_lon_max READ (2, *) grid_lat_min, d_lat, grid_lat_max nGRDcols = NINT(1 + (grid_lon_max - grid_lon_min)/d_lon) nGRDrows = NINT(1 + (grid_lat_max - grid_lat_min)/d_lat) IF (ABS((grid_lon_max - grid_lon_min) - 360.0) < 0.01) THEN ! Earth-girdling grid includes redundant extra column (last = first) box_lon_min = grid_lon_min box_lon_max = grid_lon_max ELSE ! NOT girdling, or does not include extra column: allow for half-width of cells in first and last column: box_lon_min = grid_lon_min - (0.5 * d_lon) box_lon_max = grid_lon_max + (0.5 * d_lon) END IF IF (grid_lat_max == 90.0) THEN ! one cannot go beyond the pole box_lat_max = 90.0 ELSE box_lat_max = grid_lat_max + (0.5 * d_lat) END IF IF (grid_lat_min == -90.0) THEN ! one cannot go beyond the pole box_lat_min = -90.0 ELSE box_lat_min = grid_lat_min - (0.5 * d_lat) END IF ALLOCATE ( forecast(nForecasts, nGRDrows, nGRDcols) ) ! To be filled immediately with competing forecasts. ALLOCATE ( d_tau(nGRDrows, nGRDcols) ) ! Null hypothesis for the ASS-test; to be filled when "old" pre-test-window .EQC catalog is read, below. READ (2, *) ((forecast(1, row, col), col = 1, nGRDcols), row = 1, nGRDrows) WRITE (*, "(' Here are the last 2 lines of this file:')") WRITE (3, "(' Here are the last 2 lines of this file:')") READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) CLOSE(2) IF (nForecasts > 1) THEN DO j = 2, nForecasts WRITE (*, "(' ')") WRITE (3, "(' ')") WRITE (*, "(' Forecast #',I3)") j WRITE (3, "(' Forecast #',I3)") j 210 CALL Prompt_for_String('Enter [Drive:][\path\]filename (.GRD) for next forecast:','',pathfile) WRITE (3, "(' Enter [Drive:][\path\]filename (.GRD) for next forecast:')") WRITE (3, "(' ',A)") TRIM(pathfile) GRD_pathfile(j) = pathfile OPEN (UNIT = 2, FILE = TRIM(GRD_pathfile(j)), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Forecast .GRD file ',A,' not found (in current folder).')") TRIM(GRD_pathfile(j)) CALL Pause() GO TO 210 END IF READ (2, *) ! grid_lon_min, d_lon, grid_lon_max (MUST MATCH!) READ (2, *) ! grid_lat_min, d_lat, grid_lat_max (MUST MATCH!) READ (2, *) ((forecast(j, row, col), col = 1, nGRDcols), row = 1, nGRDrows) WRITE (*, "(' Here are the last 2 lines of this file:')") WRITE (3, "(' Here are the last 2 lines of this file:')") READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) READ (2, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) CLOSE(2) END DO END IF ! nForecasts > 1 ! end reading of forecast(s) !----------------------------------------------------------------------------------------------- ! Optionally, read a mask file to block scoring in certain cells: WRITE (*, *) WRITE (*, *) WRITE (*, "(' Optionally, you may now enter an additional .GRD file to mask out certain cells')") WRITE (*, "(' of the forecast, so that they will not be considered in any of the scoring.')") WRITE (*, "(' For example:')") WRITE (*, "(' *If scoring against an historical catalog which is very incomplete')") WRITE (*, "(' in offshore regions, you could mask out all cells that are underwater.')") WRITE (*, "(' *To score forecasts only within a certain country, you could mask out')") WRITE (*, "(' foreign countries that appear at the edges of the forecast grid.')") WRITE (*, "(' ')") WRITE (*, "(' A masking file must be in .GRD format, with value 0 for any cell that is to')") WRITE (*, "(' be scored in the usual way, and non-zero (e.g., 1) for any cell to be masked.')") WRITE (*, "(' The masking .GRD file can be larger than the forecast .GRD file; however,')") WRITE (*, "(' it is not allowed to be smaller.')") WRITE (*, "(' The masking .GRD file may use different grid-spacings than the forecast.')") WRITE (*, "(' ')") WRITE (3, *) WRITE (3, *) WRITE (3, "(' Optionally, you may now enter an additional .GRD file to mask out certain cells')") WRITE (3, "(' of the forecast, so that they will not be considered in any of the scoring.')") WRITE (3, "(' For example:')") WRITE (3, "(' *If scoring against an historical catalog which is very incomplete')") WRITE (3, "(' in offshore regions, you could mask out all cells that are underwater.')") WRITE (3, "(' *To score forecasts only within a certain country, you could mask out')") WRITE (3, "(' foreign countries that appear at the edges of the forecast grid.')") WRITE (3, "(' ')") WRITE (3, "(' A masking file must be in .GRD format, with value 0 for any cell that is to')") WRITE (3, "(' be scored in the usual way, and non-zero (e.g., 1) for any cell to be masked.')") WRITE (3, "(' The masking .GRD file can be larger than the forecast .GRD file; however,')") WRITE (3, "(' it is not allowed to be smaller.')") WRITE (3, "(' The masking .GRD file may use different grid-spacings than the forecast.')") WRITE (3, "(' ')") 260 CALL Prompt_for_Logical('Do you wish to use a mask file?:', .FALSE., using_mask) WRITE (3, "(' Do you wish to use a mask file?:', L6)") using_mask IF (using_mask) THEN CALL Prompt_for_String('Enter [Drive:][\path\]filename of masking .GRD file:', 'underwater5.grd', mask_GRD_pathfile) OPEN (UNIT = 2, FILE = TRIM(mask_GRD_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This masking GRD file was not found.')") WRITE (3, "(' ERROR: This masking GRD file was not found.')") CALL Pause() GO TO 260 END IF WRITE (3, "(' Enter [Drive:][\path\]filename of masking .GRD file:')") WRITE (3, "(' ', A)") TRIM(mask_GRD_pathfile) WRITE (3, *) READ (2, *) mask_grid_lon_min, mask_d_lon, mask_grid_lon_max READ (2, *) mask_grid_lat_min, mask_d_lat, mask_grid_lat_max mask_nRows = NINT(1 + (mask_grid_lat_max - mask_grid_lat_min) / mask_d_lat) mask_nColumns = NINT(1 + (mask_grid_lon_max - mask_grid_lon_min) / mask_d_lon) global_mask = ((mask_nColumns * mask_d_lon) >= 359.9) IF (global_mask) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' This mask has global longitude range, so longitude cycling is permitted.')") WRITE (3, "(' This mask has global longitude range, so longitude cycling is permitted.')") ELSE WRITE (*, *) WRITE (3, *) WRITE (*, "(' This mask has only a local longitude range, so queries must stay in-bounds.')") WRITE (3, "(' This mask has only a local longitude range, so queries must stay in-bounds.')") END IF ALLOCATE ( mask_grd_I1(mask_nRows, mask_nColumns) ) READ (2, *, IOSTAT = ios) ((mask_grd_I1(i, j), j = 1, mask_nColumns), i = 1, mask_nRows) IF (ios /= 0) THEN WRITE (*, "(' ERROR in reading MASK.GRD file; check format and size.')") WRITE (3, "(' ERROR in reading MASK.GRD file; check format and size.')") CALL Pause() STOP ELSE WRITE (*, "(' Mask file ',A,' was read.')") TRIM(mask_GRD_pathfile) WRITE (3, "(' Mask file ',A,' was read.')") TRIM(mask_GRD_pathfile) CALL Pause() END IF CLOSE (2) !Translate between grids of different sizes, from (mask_nRows, mask_nColumns) to (nGRDrows, nGRDcols): ALLOCATE ( cell_in_use(nGRDrows, nGRDcols) ) DO i = 1, nGRDrows lat = grid_lat_max - (i-1) * d_lat i_mask = 1 + NINT((mask_grid_lat_max - lat) / mask_d_lat) IF (i_mask == 0) THEN ! fix latitude that is just-barely-out-of-range IF ((lat - mask_grid_lat_max) <= mask_d_lat) THEN i_mask = 1 END IF END IF ! i_mask == 0 (possibly fixable) IF (i_mask == (mask_nRows + 1)) THEN ! fix latitude that is just-barely-out-of-range IF ((mask_grid_lat_min - lat) <= mask_d_lat) THEN i_mask = mask_nRows END IF END IF ! i_mask == (mask_nRows + 1) (possibly fixable) IF ((i_mask < 1).OR.(i_mask > mask_nRows)) THEN WRITE (*, "(' ERROR: Attempt to read mask_grid_I1(row ',I6,', any column) is out-of-range!')") i_mask WRITE (3, "(' ERROR: Attempt to read mask_grid_I1(row ',I6,', any column) is out-of-range!')") i_mask CALL Pause() STOP END IF DO j = 1, nGRDcols lon = grid_lon_min + (j-1) * d_lon j_mask = 1 + NINT((lon - mask_grid_lon_min) / mask_d_lon) IF (global_mask) THEN IF (j_mask < 1) THEN test_lon = lon + 360.0 j_mask = MAX(1, MIN(mask_nColumns , 1 + NINT((test_lon - mask_grid_lon_min) / mask_d_lon))) ELSE IF (j_mask > mask_nColumns) THEN test_lon = lon - 360.0 j_mask = MAX(1, MIN(mask_nColumns, 1 + NINT((test_lon - mask_grid_lon_min) / mask_d_lon))) END IF ELSE ! local mask IF ((j_mask < 1).OR.(j_mask > mask_nColumns)) THEN WRITE (*, "(' ERROR: Attempt to read mask_grid_I1(',I6,', ',I6,') is out-of-range!')") i_mask, j_mask WRITE (3, "(' ERROR: Attempt to read mask_grid_I1(',I6,', ',I6,') is out-of-range!')") i_mask, j_mask CALL Pause() STOP END IF END IF !****************************************************** cell_in_use(i, j) = (mask_grd_I1(i_mask, j_mask) == 0) !****************************************************** END DO END DO !Report summary statistics: all_cells = nGRDrows * nGRDcols unmasked_cells = 0 DO i = 1, nGRDrows DO j = 1, nGRDcols IF (cell_in_use(i, j)) unmasked_cells = unmasked_cells + 1 END DO END DO masked_cells = all_cells - unmasked_cells WRITE (*, *) WRITE (3, *) WRITE (*, "(' Referring to the seismicity-forecast grid,')") WRITE (3, "(' Referring to the seismicity-forecast grid,')") WRITE (*, "(' Out of ',I10,' total cells, ',I10,' were masked and ',I10,' were not.')") all_cells, masked_cells, unmasked_cells WRITE (3, "(' Out of ',I10,' total cells, ',I10,' were masked and ',I10,' were not.')") all_cells, masked_cells, unmasked_cells CALL Pause() ELSE ! .NOT.using_mask, so set up dummy cell_in_use array which is always .TRUE. ALLOCATE ( cell_in_use(nGRDrows, nGRDcols) ) cell_in_use = .TRUE. ! whole matrix (and, re-used if there are multiple forecasts) masked_cells = 0 unmasked_cells = nGRDrows * nGRDcols END IF ! using_mask, or not !----------------------------------------------------------------------------------------------- ! Decide which version of the Area Skill Score (ASS) test will be run: WRITE (*, *) WRITE (*, "(' The final test that might be run is the Area Skill Score (ASS) test')") WRITE (*, "(' described by Zechar & Jordan [2008]. Please choose the reference')") WRITE (*, "(' model (which you are hoping to beat) for this test:')") WRITE (*, "(' (0) [do not perform any kind of ASS test]')") WRITE (*, "(' (1) uniform seismicity across the forecast spatial window')") WRITE (*, "(' (2) prior Relative Intensity (RI), from a pre-test catalog')") WRITE (*, "(' (3) any model in an equal-size, aligned .GRD file')") WRITE (*, "(' -----------------------------------------------------------------')") WRITE (*, "(' (Hints: If you don''t have a prior [pre-test] catalog .EQC, )')") ! N.B. Right ')'s appear misaligned due to double apostrophes; WRITE (*, "(' ( or an equal-size, aligned .GRD file then select: 1 . )')") ! however, the alignment should be correct after compilation. WRITE (*, "(' ( Also, I don''t consider RI a viable model, because )')") WRITE (*, "(' ( it always fails a likelihood test (e.g., L, S, or I_1) )')") WRITE (*, "(' ( as soon as there is a single ''surprise'' earthquake )')") WRITE (*, "(' ( outside the set of cells which have had EQs before. )')") WRITE (*, "(' -----------------------------------------------------------------')") 300 CALL Prompt_for_Integer('Enter integer index for reference model:', 0, ASS_type_index) IF ((ASS_type_index < 0).OR.(ASS_type_index > 3)) THEN WRITE (*, "(' ERROR: Please choose an integer from the list above.')") CALL Pause() GO TO 300 END IF WRITE (3, *) WRITE (3, "(' The final test that might be run is the Area Skill Score (ASS) test')") WRITE (3, "(' described by Zechar & Jordan [2008]. Please choose the reference')") WRITE (3, "(' model (which you are hoping to beat) for this test:')") WRITE (3, "(' (0) [do not perform any kind of ASS test]')") WRITE (3, "(' (1) uniform seismicity across the forecast spatial window')") WRITE (3, "(' (2) prior relative intensity (RI), from a pre-test catalog')") WRITE (3, "(' (3) any model in an equal-size, aligned .GRD file')") WRITE (3, "(' -----------------------------------------------------------------')") WRITE (3, "(' (Hints: If you don''t have a prior [pre-test] catalog .EQC, )')") WRITE (3, "(' ( or an equal-size, aligned .GRD file then select: 1 . )')") WRITE (3, "(' ( Also, I don''t consider RI a viable model, because )')") WRITE (3, "(' ( it always fails a likelihood test (e.g., L, S, or I_1) )')") WRITE (3, "(' ( as soon as there is a single ''surprise'' earthquake )')") WRITE (3, "(' ( outside the set of cells which have had EQs before. )')") WRITE (3, "(' -----------------------------------------------------------------')") WRITE (3, "(' Enter integer index for reference model:', I3)") ASS_type_index IF (ASS_type_index == 2) THEN !For use in ASS-test, read an older (pre-test-time-window) seismic catalog (.EQC file) !which covers (at least) the same geographic area, depth range, and magnitude range !as the forecast(s) being tested in this run. The magnitude threshold may be lower. !Note that overlaps in area, depth, and magnitude are not necessary and will not be used. !Catalog duration should be as long as possible, but not including the test time window! WRITE (*, *) WRITE (3, *) WRITE (*, "(' For use in ASS-test, read an older (pre-test-time-window) seismic catalog')") WRITE (*, "(' which covers (at least) the same geographic area, depth, and magnitude')") WRITE (*, "(' ranges as the forecast(s) being tested in this run.')") WRITE (*, "(' The magnitude threshold MAY BE LOWER than that of the forecast,')") WRITE (*, "(' as this can be desirable to get the best estimate of prior RI;')") WRITE (*, "(' however, the EQC file supplied should be complete, so it should')") WRITE (*, "(' be truncated from below at (or above) its completeness magnitude.')") WRITE (*, "(' Note that overlaps in area, depth, and magnitude are not necessary.')") WRITE (3, "(' For use in ASS-test, read an older (pre-test-time-window) seismic catalog')") WRITE (3, "(' which covers (at least) the same geographic area, depth, and magnitude')") WRITE (3, "(' ranges as the forecast(s) being tested in this run.')") WRITE (3, "(' The magnitude threshold MAY BE LOWER than that of the forecast,')") WRITE (3, "(' as this can be desirable to get the best estimate of prior RI;')") WRITE (3, "(' however, the EQC file supplied should be complete, so it should')") WRITE (3, "(' be truncated from below at (or above) its completeness magnitude.')") WRITE (3, "(' Note that overlaps in area, depth, and magnitude are not necessary.')") WRITE (*, *) WRITE (3, *) 310 CALL Prompt_for_String ('[Drive:][\path\]filename of the old .EQC file:', '', old_EQC_pathfile) WRITE (3, "(' [Drive:][\path\]filename of the old .EQC file:')") WRITE (3, "(' ',A)") TRIM(old_EQC_pathfile) OPEN (UNIT = 1, FILE = TRIM(old_EQC_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: ',A,' could not be found (in current folder).')") TRIM(old_EQC_pathfile) WRITE (*,"(' Please check [Drive:][\path\]filename and/or move file, and try again...')") CALL Pause() GO TO 310 END IF !Count lines which include a valid year number: old_EQC_lines = 0 old_count_lines: DO READ (1, "(9X,I5)", IOSTAT = ios) iYear IF (ios == -1) EXIT old_count_lines ! EOF IF ((ios /= 0).OR.(iYear < 1800).OR.(iYear > 2099)) THEN WRITE (*,"(' After reading ',I7,' lines, met bad(?) year number: ',I5)") old_EQC_lines, iYear CALL Pause() ELSE old_EQC_lines = old_EQC_lines + 1 END IF BACKSPACE(1) READ (1, "(A)") text_line IF (old_EQC_lines == 1) old_EQC_line_1 = TRIM(text_line) old_EQC_last_line = TRIM(text_line) END DO old_count_lines CLOSE(1) WRITE (*, "(' This catalog file includes ',I10,' events.')") old_EQC_lines WRITE (3, "(' This catalog file includes ',I10,' events.')") old_EQC_lines WRITE (*, *) WRITE (3, *) WRITE (*, "(' Here are the first and last lines of the old .EQC file:')") WRITE (3, "(' Here are the first and last lines of the old .EQC file:')") WRITE (*, "(' ',A)") TRIM(old_EQC_line_1) WRITE (3, "(' ',A)") TRIM(old_EQC_line_1) WRITE (*, "(' ',A)") TRIM(old_EQC_last_line) WRITE (3, "(' ',A)") TRIM(old_EQC_last_line) CALL Prompt_for_Logical ('Is its time-window entirely BEFORE that of the forecast?',.FALSE.,isTimeWindowPrevious) WRITE (3, "(' Is its time-window entirely BEFORE that of the forecast?')") WRITE (3, "(L5)") isTimeWindowPrevious IF (.NOT.isTimeWindowPrevious) THEN WRITE (*, "(' Please change the catalog time-window, or change catalogs, and try again.')") WRITE (*, "(' It would be seriously misleading to use test-window seismicity as ""prior"" relative intensity.')") WRITE (3, "(' Please change the catalog time-window, or change catalogs, and try again.')") WRITE (3, "(' It would be seriously misleading to use test-window seismicity as ""prior"" relative intensity.')") CALL Pause() GO TO 310 END IF WRITE (3, "(' Is this catalog global in its areal coverage?')") WRITE (3, "(L5)") isEQCglobal IF (isEQCglobal) THEN isEQCframed = .TRUE. ELSE CALL Prompt_for_Logical ('Does this catalog cover the whole spatial window of the test?', .FALSE., isEQCframed) WRITE (3, "(' Does this catalog cover the whole spatial window of the test?')") WRITE (3, "(L5)") isEQCframed IF (.NOT.isEQCframed) THEN WRITE (*, "(' Please expand the older catalog''s area and try again.')") WRITE (3, "(' Please expand the older catalog''s area and try again.')") CALL Pause() GO TO 310 END IF END IF CALL Prompt_for_Logical ('Does this catalog extend (at least) as deep as the forecast(s)?', .FALSE., isEQCdeeper) WRITE (3, "(' Does this catalog extend (at least) as deep as the forecast(s)?')") WRITE (3, "(L5)") isEQCdeeper IF (.NOT.isEQCdeeper) THEN WRITE (*, "(' Please expand the depth limit of the old reference catalog, and try again.')") WRITE (3, "(' Please expand the depth limit of the old reference catalog, and try again.')") CALL Pause() GO TO 310 END IF CALL Prompt_for_Logical ('Is this catalog complete above the threshold of the forecast(s)?', .FALSE., isEQCcomplete) WRITE (3, "(' Is this catalog complete above the threshold of the forecast(s)?')") WRITE (3, "(L5)") isEQCcomplete IF (.NOT.isEQCcomplete) THEN WRITE (*, "(' Please lower the catalog magnitude threshold, or change catalogs, and try again.')") WRITE (*, "(' It could be seriously misleading to test use incomplete seismicity as ""prior"" relative intensity.')") WRITE (3, "(' Please lower the catalog magnitude threshold, or change catalogs, and try again.')") WRITE (3, "(' It could be seriously misleading to test use incomplete seismicity as ""prior"" relative intensity.')") CALL Pause() GO TO 310 END IF ALLOCATE ( old_eq_Elon(old_EQC_lines) ) ALLOCATE ( old_eq_Nlat(old_EQC_lines) ) ALLOCATE ( old_eq_depth_int(old_EQC_lines) ) ALLOCATE ( old_eq_mag(old_EQC_lines) ) OPEN (UNIT = 1, FILE = TRIM(old_EQC_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) next_line_2: DO iEq = 1, old_EQC_lines READ (1, 320, IOSTAT = ios) & & i, c2a, c2b, & & c2c, c2d, c2e, c1, & & old_eq_Elon(iEq), old_eq_Nlat(iEq), & & old_eq_depth_int(iEq), old_eq_mag(iEq) !Note: Any appended data, such as focal mechanism, is ignored. 320 FORMAT (9X, & ! N.B. This format is the same as 120. & I5,'.',A2,'.',A2, 1X, & ! using I5 for negative years (B.C.) & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) IF (ios /= 0) EXIT next_line_2 END DO next_line_2 ! iEq = 1, old_EQC_lines CLOSE(1) ELSE IF (ASS_type_index == 3) THEN ! (reading any equal-size, aligned .GRD file) 350 WRITE (*, *) WRITE (*, "(' Please supply the reference forecast that you wish to use in the ASS test.')") WRITE (*, "(' You will need to input an equal-sized, aligned .GRD file giving')") WRITE (*, "(' epicentroid rate densities (per unit area & time), rather than EQ counts.')") WRITE (3, *) WRITE (3, "(' Please supply the reference forecast that you wish to use in the ASS test.')") WRITE (3, "(' You will need to input an equal-sized, aligned .GRD file giving')") WRITE (3, "(' epicentroid rate densities (per unit area & time), rather than EQ counts.')") CALL Prompt_for_String('Enter [Drive:][\path\]filename (.GRD) for this forecast:','',reference_forecast_pathfile) WRITE (3, "(' Enter [Drive:][\path\]filename (.GRD) for this forecast:')") WRITE (3, "(' ',A)") TRIM(reference_forecast_pathfile) OPEN (UNIT = 12, FILE = TRIM(reference_forecast_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Reference forecast .GRD file ',A,' not found (in current folder).')") TRIM(reference_forecast_pathfile) CALL Pause() GO TO 350 END IF WRITE (*, *) WRITE (3, *) WRITE (*, "(' Here are the first 2 lines of this file:')") WRITE (3, "(' Here are the first 2 lines of this file:')") READ (12, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) READ (12, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (3, "(' ',A)") TRIM(text_line) CLOSE (12) OPEN (UNIT = 12, FILE = TRIM(reference_forecast_pathfile), STATUS = "OLD", IOSTAT = ios) READ (12, *) t_lon_min, t_d_lon, t_lon_max READ (12, *) t_lat_min, t_d_lat, t_lat_max IF ((t_lon_min /= grid_lon_min).OR.(t_d_lon /= d_lon).OR.(t_lon_max /= grid_lon_max).OR. & & (t_lat_min /= grid_lat_min).OR.(t_d_lat /= d_lat).OR.(t_lat_max /= grid_lat_max)) THEN WRITE (*, "(' ERROR: This reference forecast is NOT equal-sized and aligned!')") WRITE (*, "(' All 6 grid descriptors in the first two lines must match!')") WRITE (3, "(' ERROR: This reference forecast is NOT equal-sized and aligned!')") WRITE (3, "(' All 6 grid descriptors in the first two lines must match!')") CALL Pause() GO TO 350 END IF ! grids are not aligned ALLOCATE ( reference_forecast(nGRDrows, nGRDcols) ) READ (12, *)((reference_forecast(i, j), j = 1, nGRDcols), i = 1, nGRDrows) WRITE (*, "(' Reference forecast has been read.')") WRITE (3, "(' Reference forecast has been read.')") CLOSE (12) END IF ! ASS_type_index == 2 (so reading pre-test-time-window .EQC and storing it), or == 3 (reading any equal-size, aligned .GRD file) !======================================== end INPUTS ========================================================== !======================================== begin OUTPUTS ======================================================= WRITE (*, *) WRITE (3, *) WRITE (*, "(' INPUT section concluded; now beginning SCORING and OUTPUT:')") WRITE (*, "(' ---------------------------------------------------------')") WRITE (3, "(' INPUT section concluded; now beginning SCORING and OUTPUT:')") WRITE (3, "(' ---------------------------------------------------------')") !convert EQ rates in SI units of #/m^2/s into absolute (real) numbers per test time-window: ALLOCATE ( dt_da_factor(nGRDrows, nGRDcols) ) !first, dt = length of test time-window: dt_da_factor = new_EQC_duration_years * 3.15576E7 !next, include # of square meters (neglecting edge effects) DO row = 1, NGRDrows latitude = grid_lat_max - (row - 1) * d_lat lat_factor = COS(latitude / 57.29578) dt_da_factor(row, 1:nGRDcols) = 1.2364E10 * d_lat * d_lon * lat_factor * dt_da_factor(row, 1:nGRDcols) END DO !last, trim E & W edges for reduced area, if necessary: IF (box_lon_min == grid_lon_min) THEN dt_da_factor(1:nGRDrows, 1 ) = dt_da_factor(1:nGRDrows, 1 ) / 2.0 dt_da_factor(1:nGRDrows, nGRDcols) = dt_da_factor(1:nGRDrows, nGRDcols) / 2.0 END IF !Count number of qualifying events in actual new (test-time-window) catalog (and list them, if they are not too many): ALLOCATE ( selected(new_EQC_lines) ) ! LOGICAL: does this EQ count? ALLOCATE ( actual_n_in_bin(nGRDrows, nGRDcols) ) CALL Count_Events (new_EQC_lines, new_eq_Elon, new_eq_Nlat, new_eq_depth_int, new_eq_mag, & ! critical parts of catalog (excluding dates & times) & box_lon_min, box_lon_max, box_lat_min, box_lat_max, threshold, maxDepthKm, & ! box selection criteria (maxDepthKm = integer) & using_mask, cell_in_use, grid_lon_min, d_lon, grid_lon_max, grid_lat_min, d_lat, grid_lat_max, nGRDrows, nGRDcols, & ! mask description & N, selected) ! results N_true = N ! save under unique name CALL Assign_to_Bins (new_EQC_lines, selected, new_eq_Elon, new_eq_Nlat, & ! inputs describing catalog & d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max, nGRDrows, nGRDcols, & !inputs describing GRD & using_mask, cell_in_use, & ! inputs describing mask & actual_n_in_bin) ! output; will NOT be overwritten by modifications below, which go to array mod_n_in_bin. IF (using_mask) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' New (test-time-window) catalog contains ',I6,' EQs')") N WRITE (*, "(' in box, unmasked, above threshold, & not-too-deep.')") WRITE (3, "(' New (test-time-window) catalog contains ',I6,' EQs')") N WRITE (3, "(' in box, unmasked, above threshold, & not-too-deep.')") ELSE ! .NOT. using_mask WRITE (*, *) WRITE (3, *) WRITE (*, "(' New (test-time-window) catalog contains ',I6,' EQs')") N WRITE (*, "(' in (lon,lat) box, above threshold, & not-too-deep.')") WRITE (3, "(' New (test-time-window) catalog contains ',I6,' EQs')") N WRITE (3, "(' in (lon,lat) box, above threshold, & not-too-deep.')") END IF IF (N <= 40) THEN DO iEq = 1, new_EQC_lines IF (selected(iEq)) THEN WRITE (*, 120) & & new_eq_year(iEq), new_eq_month(iEq), new_eq_day(iEq), & & new_eq_hour(iEq), new_eq_minute(iEq), new_eq_second(iEq), new_eq_tenths(iEq), & & new_eq_Elon(iEq), new_eq_Nlat(iEq), & & new_eq_depth_int(iEq), new_eq_mag(iEq) WRITE (3, 120) & & new_eq_year(iEq), new_eq_month(iEq), new_eq_day(iEq), & & new_eq_hour(iEq), new_eq_minute(iEq), new_eq_second(iEq), new_eq_tenths(iEq), & & new_eq_Elon(iEq), new_eq_Nlat(iEq), & & new_eq_depth_int(iEq), new_eq_mag(iEq) !Note: Any appended data, such as focal mechanism, is not written because it was not stored. END IF END DO END IF !Integrate each forecast across all its spatial bins to get the total number of earthquakes expected. ALLOCATE ( N_sup_j(nForecasts) ) ! expectation N^j = SUM lambda_i^j; this array is REAL WRITE (*, *) WRITE (3, *) DO j = 1, nForecasts CALL Integrate_forecast(forecast, j, & & nGRDrows, d_lat, grid_lat_max, & & nGRDcols, box_lon_min, grid_lon_min, d_lon, box_lon_max, & & cell_in_use, & & N_per_year) ! output N_per_year_is DOUBLE PRECISION N_sup_j(j) = N_per_year * new_EQC_duration_years IF (using_mask) THEN WRITE (*, "(' After masking, forecast #', I2,' is for ',F10.2,' EQs in ',F6.2,' years.')") j, N_sup_j(j), new_EQC_duration_years WRITE (3, "(' After masking, forecast #', I2,' is for ',F10.2,' EQs in ',F6.2,' years.')") j, N_sup_j(j), new_EQC_duration_years ELSE WRITE (*, "(' Forecast #', I2,' is for ',F10.2,' EQs in ',F6.2,' years.')") j, N_sup_j(j), new_EQC_duration_years WRITE (3, "(' Forecast #', I2,' is for ',F10.2,' EQs in ',F6.2,' years.')") j, N_sup_j(j), new_EQC_duration_years END IF END DO !Compute scaled forecasts (for the S-test) which all have expected N = actual catalog N (not modified). ALLOCATE ( scaled_forecast(nForecasts, nGRDrows, nGRDcols) ) DO j = 1, nForecasts IF (N_sup_j(j) > 0.0) THEN amplification = (1.0 * N_true) / N_sup_j(j) DO row = 1, nGRDrows DO col = 1, nGRDcols scaled_forecast(j, row, col) = amplification * forecast(j, row, col) END DO END DO ELSE WRITE (*, "(' ERROR: Cannot scale forecast #',I2,' because it has an expectation of zero!')") j CALL Pause() STOP END IF END DO !Compute Poisson cumulative distribution function, PoissonCDF(0:comb_limit, 1:nForecasts), for each forecast, centered on N_sup_j value. DO j = 1, nForecasts lambda = N_sup_j(j) ! Number of earthquakes expected by this forecast, in whole box (or its unmasked part) and time-window. !Note that lambda is not an integer, but a REAL number, e.g., 37.443 or 1023.237 IF (NINT(REAL(lambda + 3.0 * SQRT(lambda))) > Poisson_limit) THEN WRITE (*, *) WRITE (*, "(' ERROR: Poisson_limit = ',I5,' is too low when lambda = ',F10.2)") Poisson_limit, lambda WRITE (*, "(' Increase this parameter and recompile, or else ...')") WRITE (*, *) CALL Pause() STOP END IF exact_Poisson = ((lambda + 3.0 * SQRT(lambda)) < comb_limit) ! This is ONE test that must be passed. exact_Poisson = exact_Poisson .AND. ((Factorial_limit * DLOG10(lambda)) < 300.D0) ! We must ALSO check that lambda**omega will not overflow the REAL*8 limit of about 1E+304 IF (exact_Poisson) THEN !Compute probabilities of different integer outcomes, given this expectation: probability_of_integer = 0.0D0 ! whole DP vector (0:comb_limit) initialized total_probability = 0.0D0 ! DP scalar initialized create_comb: 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 ! REAL*8 <--- INTEGER factor1 = lambda**omega_R8 factor2 = DEXP(-lambda) factor3 = 1.0D0 / Factorial(omega) probability_of_integer(omega) = factor1 * factor2 * factor3 !------------------------------------------------------------------------------------ total_probability = total_probability + probability_of_integer(omega) ! scalar DP sum END DO create_comb probability_of_integer((comb_limit + 1):Poisson_limit) = 0.0D0 IF ((total_probability < 0.99D0).OR.(total_probability > 1.01D0)) THEN WRITE (*, "(' ERROR: total_probability = ',ES12.4,' when lambda = ',ES12.4,' and Factorial_limit = ',I6)") total_probability, lambda, Factorial_limit CALL Pause() STOP END IF !Compute cumulative probability of different integer outcomes: PoissonCDF(0, j) = probability_of_integer(0) DO omega = 1, Poisson_limit PoissonCDF(omega, j) = PoissonCDF(omega - 1, j) + probability_of_integer(omega) END DO ELSE ! use normal distribution as a proxy for the Poisson distribution: root_lambda = DSQRT(lambda) DO omega = 0, Poisson_limit argument = (omega - lambda) / root_lambda IF (argument < -5.9D0) THEN PoissonCDF(omega, j) = 0.0D0 ELSE IF (argument > 5.9D0) THEN PoissonCDF(omega, j) = 1.0D0 ELSE ! limiting the function calls to -5.9 < argument < +5.9 to avoid under/overflows. !---- IMSL version: -------------------------------------------------------------------------------------- !PoissonCDF(omega, j) = DNORDF(argument) !Using the standard (normal) Gaussian distribution function, with input argument = # of sigmas from mean. !----- MKL version: -------------------------------------------------------------------------------------- VDCdfNorm_in(1) = argument CALL VDcdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) PoissonCDF(omega, j) = VDCdfNorm_out(1) !--------------------------------------------------------------------------------------------------------- END IF END DO END IF END DO ! j = 1, nForecasts !Compute log-likelihood for each forecast, using actual catalog: ALLOCATE ( L_sup_j(nForecasts) ) WRITE (*, *) WRITE (3, *) DO j = 1, nForecasts CALL LogLikelihood (j, nForecasts, forecast, dt_da_factor, & ! inputs describing forecast & nGRDrows, nGRDcols, & & cell_in_use, & ! describes masking (or, all .TRUE. if no mask) & actual_n_in_bin, & ! input with actual counts of qualifying earthquakes & L ) ! output: REAL*4 log-likelihood L (using natural logarithms) L_sup_j(j) = L WRITE (*, "(' Log-likelihood of forecast #',I2,' = ',F10.2,', based on original catalog.')") j, L WRITE (3, "(' Log-likelihood of forecast #',I2,' = ',F10.2,', based on original catalog.')") j, L END DO !Compute log-likelihood for each scaled forecast (for the S-test), using actual catalog: ALLOCATE ( scaled_L_sup_j(nForecasts) ) WRITE (*, *) WRITE (3, *) DO j = 1, nForecasts CALL LogLikelihood (j, nForecasts, scaled_forecast, dt_da_factor, & ! inputs describing forecast & nGRDrows, nGRDcols, & & cell_in_use, & ! describes masking (or, all .TRUE. if no mask) & actual_n_in_bin, & ! input with actual counts of qualifying earthquakes & L ) ! output: log-likelihood (using natural logarithms) scaled_L_sup_j(j) = L WRITE (*, "(' Log-likelihood of scaled forecast #',I2,' = ',F10.2,', based on original catalog.')") j, L WRITE (3, "(' Log-likelihood of scaled forecast #',I2,' = ',F10.2,', based on original catalog.')") j, L END DO !Generate "s_in_use" different modifications of the catalog, and count events and compute per-forecast likelihoods in each: ALLOCATE ( mod_n_in_bin(nGRDrows, nGRDcols) ) ! Note that this will be re-used for each random mod of catalog. ALLOCATE ( R_tilde_3D_matrix(s_in_use, nForecasts, nForecasts) ) mag_offset = 0.5 * mag_sigma**2 * 2.303 ! N.B. 2.303 is actually b * ln(10), where b is from Gutenberg-Richter plot. ! Here I have just assumed b = 1 for simplicity. ! This offset formula is due to (10) in Rhoades & Dowrick [2000]; ! this term is needed to avoid artifically inflating the number ! of countable EQs when we add Gaussian random noise to their magnitudes! !NOTE that Zechar et al. [2010, BSSA] used a Laplace distribution to modify the magnitudes; perhaps this ! is an alternative method that avoids the need for my term "mag_offset"? (However, they don't discuss this.) IF (randomize_test_EQs) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' Computing ', I5, ' modifications of the catalog...')") s_in_use WRITE (3, "(' Computing ', I5, ' modifications of the catalog...')") s_in_use ELSE WRITE (*, *) WRITE (3, *) WRITE (*, "(' Computing tests with no random perturbations of test catalog...')") WRITE (3, "(' Computing tests with no random perturbations of test catalog...')") END IF DO q = 1, s_in_use ! Loop to create "s_in_use" different modifications of the observed catalog, count N_tilde's, and score L_tilde_matrix's. !Randomly modify the magnitude: !Get one random number in range (0., 1.) per earthquake: !--- IMSL method: ---------------------------------------------------------------------- !CALL RNUN(new_EQC_lines, mod_temporary) ! get one random number (0, 1) per earthquake !--- MKL method: ------------------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = new_EQC_lines, r = mod_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------------------- DO iEq = 1, new_EQC_lines !--- IMSL version: ------------------------------------------------------------------------------ !sigmas = ANORIN(mod_temporary(iEq)) ! convert to # of sigmas, in a standard normal distribution. !--- MKL version: ------------------------------------------------------------------------------- VSCdfNormInv_in(1) = mod_temporary(iEq) CALL VSCdfNormInv(I4_one, VSCdfNormInv_in, VSCdfNormInv_out) sigmas = VSCdfNormInv_out(1) !------------------------------------------------------------------------------------------------ mod_mag(iEq) = new_eq_mag(iEq) + mag_sigma * sigmas - mag_offset END DO !Randomly modify the depth: !Get one random number in range (0., 1.) per earthquake: !--- IMSL method: ---------------------------------------------------------------------- !CALL RNUN(new_EQC_lines, mod_temporary) ! get one random number (0, 1) per earthquake !--- MKL method: ------------------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = new_EQC_lines, r = mod_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------------------- DO iEq = 1, new_EQC_lines !--- IMSL version: ------------------------------------------------------------------------------ !sigmas = ANORIN(mod_temporary(iEq)) ! convert to # of sigmas, in a standard normal distribution. !--- MKL version: ------------------------------------------------------------------------------- VSCdfNormInv_in(1) = mod_temporary(iEq) CALL VSCdfNormInv(I4_one, VSCdfNormInv_in, VSCdfNormInv_out) sigmas = VSCdfNormInv_out(1) !------------------------------------------------------------------------------------------------ mod_depth_int(iEq) = MAX(0, NINT(new_eq_depth_int(iEq) + depth_sigma_km * sigmas)) END DO !Randomly modify the latitude: !Get one random number in range (0., 1.) per earthquake: !--- IMSL method: ---------------------------------------------------------------------- !CALL RNUN(new_EQC_lines, mod_temporary) ! get one random number (0, 1) per earthquake !--- MKL method: ------------------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = new_EQC_lines, r = mod_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------------------- DO iEq = 1, new_EQC_lines !--- IMSL version: ------------------------------------------------------------------------------ !sigmas = ANORIN(mod_temporary(iEq)) ! convert to # of sigmas, in a standard normal distribution. !--- MKL version: ------------------------------------------------------------------------------- VSCdfNormInv_in(1) = mod_temporary(iEq) CALL VSCdfNormInv(I4_one, VSCdfNormInv_in, VSCdfNormInv_out) sigmas = VSCdfNormInv_out(1) !------------------------------------------------------------------------------------------------ mod_Nlat(iEq) = new_eq_Nlat(iEq) + EW_or_NS_sigma_km * sigmas/111.2 END DO !Randomly modify the longitude: !Get one random number in range (0., 1.) per earthquake: !--- IMSL method: ---------------------------------------------------------------------- !CALL RNUN(new_EQC_lines, mod_temporary) ! get one random number (0, 1) per earthquake !--- MKL method: ------------------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = new_EQC_lines, r = mod_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------------------- DO iEq = 1, new_EQC_lines !--- IMSL version: ------------------------------------------------------------------------------ !sigmas = ANORIN(mod_temporary(iEq)) ! convert to # of sigmas, in a standard normal distribution. !--- MKL version: ------------------------------------------------------------------------------- VSCdfNormInv_in(1) = mod_temporary(iEq) CALL VSCdfNormInv(I4_one, VSCdfNormInv_in, VSCdfNormInv_out) sigmas = VSCdfNormInv_out(1) !------------------------------------------------------------------------------------------------ mod_Elon(iEq) = new_eq_Elon(iEq) + EW_or_NS_sigma_km * sigmas/(111.2 * COS(new_eq_Nlat(iEq)/57.296)) END DO CALL Count_Events (new_EQC_lines, mod_Elon, mod_Nlat, mod_depth_int, mod_mag, & ! critical parts of modified catalog (excluding dates & times) & box_lon_min, box_lon_max, box_lat_min, box_lat_max, threshold, maxDepthKm, & ! selection criteria (maxDepthKm = integer) & using_mask, cell_in_use, grid_lon_min, d_lon, grid_lon_max, grid_lat_min, d_lat, grid_lat_max, nGRDrows, nGRDcols, & ! mask description & N, selected) ! results N_tilde(q) = N CALL Assign_to_Bins (new_EQC_lines, selected, mod_Elon, mod_Nlat, & ! inputs describing catalog & d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max, nGRDrows, nGRDcols, & !inputs describing GRD & using_mask, cell_in_use, & ! inputs describing mask & mod_n_in_bin) ! output DO j = 1, nForecasts CALL LogLikelihood (j, nForecasts, forecast, dt_da_factor, & ! inputs describing forecast & nGRDrows, nGRDcols, & & cell_in_use, & ! describes masking (or, all .TRUE. if no mask) & mod_n_in_bin, & ! input with actual counts of qualifying earthquakes & L ) ! output: log-likelihood (using natural logarithms) L_tilde_matrix(q, j) = L END DO !Expand into the R_tilde_3D_matrix: DO row = 1, nForecasts DO col = 1, nForecasts R_tilde_3D_matrix(q, row, col) = L_tilde_matrix(q, row) - L_tilde_matrix(q, col) END DO END DO END DO ! q = 1, s_in_use (different modifications of the catalog) !compute mean and sigma of N_tilde(q = 1, ... s): real_vector_s(1:s_in_use) = N_tilde(1:s_in_use) ! convert integers to reals CALL Moments (real_vector_s, s_in_use, & ! inputs & N_tilde_mean, N_tilde_sigma) ! outputs IF (randomize_test_EQs) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' After ',I6,' modifications of the catalog,'/' N_tilde_mean = ',F10.1'; N_tilde_sigma = ',F10.2)") s_in_use, N_tilde_mean, N_tilde_sigma WRITE (3, "(' After ',I6,' modifications of the catalog,'/' N_tilde_mean = ',F10.1'; N_tilde_sigma = ',F10.2)") s_in_use, N_tilde_mean, N_tilde_sigma ELSE WRITE (*, *) WRITE (3, *) WRITE (*, "(' With unchanged test catalog,'/' N_tilde_mean = ',F10.1'; N_tilde_sigma = ',F10.2)") N_tilde_mean, N_tilde_sigma WRITE (3, "(' With unchanged test catalog,'/' N_tilde_mean = ',F10.1'; N_tilde_sigma = ',F10.2)") N_tilde_mean, N_tilde_sigma END IF !compute means and standard deviations of L_tilde_matrix(q = 1, ... s, j = 1, ... nForecasts), for all j. DO j = 1, nForecasts real_vector_s(1:s_in_use) = L_tilde_matrix(1:s_in_use, j) ! convert integers to reals CALL Moments (real_vector_s, s_in_use, & ! inputs & L_tilde_mean(j), L_tilde_sigma(j)) ! outputs WRITE (*, "(' L_tilde_mean(',I2,') = ',F10.1'; L_tilde_sigma(',I2,') = ',F10.2)") j, L_tilde_mean(j), j, L_tilde_sigma(j) WRITE (3, "(' L_tilde_mean(',I2,') = ',F10.1'; L_tilde_sigma(',I2,') = ',F10.2)") j, L_tilde_mean(j), j, L_tilde_sigma(j) END DO !sort each column in L_tilde_matrix, from lowest (#1) to highest (#s) DO j = 1, nForecasts real_vector_s(1:s_in_use) = L_tilde_matrix(1:s_in_use, j) ! convert integers to reals CALL Sort (s_in_use, & ! input (length of vector) & real_vector_s) ! INTENT(INOUT) L_tilde_matrix(1:s_in_use, j) = real_vector_s(1:s_in_use) END DO IF (nForecasts > 1) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' Log-likelihood ratios R^ij = L^i - L^j (based on mean of all modifications):')") WRITE (*, "(' ----------------------------------------------------------------------------')") WRITE (*, "(' ',' j =',99I10)") (col, col = 1, nForecasts) WRITE (*, "(' ',14X,99(1X,A9))") (GRD_pathfile(col)(1:9), col = 1, nForecasts) WRITE (3, "(' Log-likelihood ratios R^ij = L^i - L^j (based on mean of all modifications):')") WRITE (3, "(' ----------------------------------------------------------------------------')") WRITE (3, "(' ',' j =',99I10)") (col, col = 1, nForecasts) WRITE (3, "(' ',14X,99(1X,A9))") (GRD_pathfile(col)(1:9), col = 1, nForecasts) ALLOCATE ( R_square(nForecasts, nForecasts) ) DO row = 1, nForecasts DO col = 1, nForecasts R_square(row, col) = L_tilde_mean(row) - L_tilde_mean(col) END DO WRITE (*, "(' ','i=',I2,1X,A9,99F10.1)") row, GRD_pathfile(row)(1:9), (R_square(row, col), col = 1, nForecasts) WRITE (3, "(' ','i=',I2,1X,A9,99F10.1)") row, GRD_pathfile(row)(1:9), (R_square(row, col), col = 1, nForecasts) END DO END IF !FIRST Simulation group (of 2): GBs_requested = 2. * m * nGRDrows * nGRDcols / (1024.**3) WRITE (*, *) WRITE (*, "(' ALLOCATing ',F8.3,' GB array sim_N_hat_3D_matrix;')") GBs_requested WRITE (*, "(' if computer has insufficient memory, this program will crash...')") WRITE (3, *) WRITE (3, "(' ALLOCATing ',F8.3,' GB array sim_N_hat_3D_matrix;')") GBs_requested WRITE (3, "(' if computer has insufficient memory, this program will crash...')") CALL Pause() ALLOCATE ( sim_N_hat_3D_matrix(m, nGRDrows, nGRDcols) ) ! BIG! For m = 1000, and 0.1-degree global grid: 12.07 GB sim_N_hat_3D_matrix = 0 ALLOCATE ( N_hat_matrix(m, nForecasts) ) ALLOCATE ( L_hat_matrix(m, nForecasts, nForecasts) ) !Conduct (k = 1, ..., m) simulations of hypothetical expected catalogs, based on each input forecast model: WRITE (*, *) WRITE (3, *) WRITE (*, "(' Computing ',I5,' simulated catalogs for each forecast model:')") m WRITE (*, "(' -----------------------------------------------------------')") WRITE (3, "(' Computing ',I5,' simulated catalogs for each forecast model:')") m WRITE (3, "(' -----------------------------------------------------------')") DO j = 1, nForecasts ! Note that this subscript identifies which forecast was the ASSUMED model for the simulation; ! further below we will also compute log-likelihoods (L_hat_matrix) for alternate forecast jj. !Note that all of sim_N_hat_3D_matrix is overwritten and reused for each new forecast. WRITE (*, "(' Forecast model j = ',I2,': ',A)") j, TRIM(GRD_pathfile(j)(1:50)) WRITE (3, "(' Forecast model j = ',I2,': ',A)") j, TRIM(GRD_pathfile(j)(1:50)) IF (using_mask) THEN WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I10,' unmasked cells...')") m, unmasked_cells WRITE (3, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I10,' unmasked cells...')") m, unmasked_cells ELSE WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I4,' rows and ',I4,' columns...')") m, nGRDrows, nGRDcols WRITE (3, "(' 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(j, 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 (*, "(' ERROR: comb_limit = ',I5,' is too low when lambda = ',F10.2)") comb_limit, lambda WRITE (*, "(' Increase parameter, or reduce test-time-window length so that this does not happen.')") WRITE (3, "(' ERROR: comb_limit = ',I5,' is too low when lambda = ',F10.2)") comb_limit, lambda WRITE (3, "(' Increase parameter, or reduce test-time-window length so that this does not happen.')") CALL Pause() STOP END IF !Compute probabilities of different integer outcomes, given this expectation: probability_of_integer = 0.0D0 ! whole DP vector initialized total_probability = 0.0D0 ! DP 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: ------------------------------------------------- probability_of_integer(omega) = (lambda ** omega) * DEXP(-lambda) / Factorial(omega) ! and note that 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 cumulative_probability(0:comb_limit) = cumulative_probability(0:comb_limit) + lost_probability cumulative_probability(comb_limit) = 1.0D0 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) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = m, r = sim_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------- 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 (sim_temporary(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...')") WRITE (3, "(' 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 forecast... 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_matrix(k, j) = N sum = sum + N END DO mean = sum / m variance = 0.0D0 ! just initializing... DO k = 1, m variance = variance + (N_hat_matrix(k, j) - mean)**2 END DO sigma = DSQRT(variance / (m - 1)) WRITE (*, "(' Mean = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (3, "(' Mean = ',F10.2,' and sigma = 'F10.2)") mean, sigma !Compute log-likelihood L_hat_matrix for each of the (k = 1, ..., m) simulations of this forecast #j, !and also for all other forecasts #jj: WRITE (*, "(' Computing ',I2,' log-likelihoods for each of these simulations...')") nForecasts WRITE (3, "(' Computing ',I2,' log-likelihoods for each of these simulations...')") nForecasts DO jj = 1, nForecasts ! Which forecast is tested? (Not necessarily the one that underlay the simulation!) DO k = 1, m ! each simulation of forecast #j L_sum = 0.0D0 ! just initializing sum DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN lambda = forecast(jj, row, col) * dt_da_factor(row, col) N_hat = sim_N_hat_3D_matrix(k, row, col) ! convert INT(2) --> INT(4); also, single access into huge array! IF (lambda == 0.0) THEN ! no EQs forecast in this bin IF (N_hat == 0) THEN ! no problem ! p = 1, log(p) = 0; no change to L ELSE L_sum = -999999.D0 ! temporary substitute for -infinity; forecast will be rejected END IF ELSE ! lamba > 0 ; OK to take log(lambda) IF (N_hat > Factorial_limit) THEN WRITE (*, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") N_hat, Factorial_limit WRITE (*, "(' Reduce length of test-time-window so that this does not happen.')") WRITE (3, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") N_hat, Factorial_limit WRITE (3, "(' Reduce length of test-time-window so that this does not happen.')") CALL Pause() STOP END IF !--- IMSL version: --------------------------------------------------------------------------------------------------- !L = L -lambda + N_hat * LOG(lambda) - LOG(DFAC(N_hat)) !N.B. DFAC() is the factorial function from IMSL; DFAC(n) = n! !--- MKL does not have a factorial; using my own: -------------------------------------------------------------------- L_sum = L_sum -lambda + N_hat * LOG(lambda) - LOG_Factorial(N_hat) !--------------------------------------------------------------------------------------------------------------------- !N.B. lambda > 0, so LOG(lambda) will not cause any crashes !N.B. 0! == 1, so last term will not cause any crashes. END IF END IF ! cell_in_use(row, col) END DO ! col = 1, nGRDcols END DO ! row = 1, nGRDrows L = L_sum ! R4 = R8; OK, as magnitude is not a problem; L_sum was used to avoid loss of tiny contributions to sum. L_hat_matrix(k, j, jj) = L ! where (k = 1, ..., m) is the simulation #, ! (j = 1, ..., nForecasts) is the basis for the simulation; and ! (jj = 1, ..., nForecasts) is another(?) forecast being scored. END DO ! k = 1, m END DO ! jj = 1, nForecasts (the one being tested, not (usually) the one which was the basis of the simulation sim_temporary(1:m) = L_hat_matrix(1:m, j, j) ! Here, reporting only the diagonal results ( jj == j ). CALL Moments (sim_temporary, m, & ! inputs & mean, sigma) ! outputs WRITE (*, "(' Diagonal (self-scored) mean log-likelihood = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (3, "(' Diagonal (self-scored) mean log-likelihood = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (*, *) ! provide blank line to separate different forecasts WRITE (3, *) ! provide blank line to separate different forecasts END DO ! j = 1, nForecasts (This is the forecast which was the basis for the simulation.) !SECOND Simulation group (of 2): !Conduct (k = 1, ..., m) simulations of hypothetical expected catalogs, based on each SCALED forecast model: WRITE (*, *) WRITE (3, *) WRITE (*, "(' Computing ',I5,' simulated catalogs for each scaled forecast model:')") m WRITE (*, "(' ------------------------------------------------------------------')") WRITE (3, "(' Computing ',I5,' simulated catalogs for each scaled forecast model:')") m WRITE (3, "(' ------------------------------------------------------------------')") ALLOCATE ( scaled_N_hat_matrix(m, nForecasts) ) ALLOCATE ( scaled_L_hat_matrix(m, nForecasts) ) DO j = 1, nForecasts ! (which, in this block of code, are always SCALED versions of forecasts) !Note that all of sim_N_hat_3D_matrix is overwritten and reused for each new scaled forecast. WRITE (*, "(' Scaled version of forecast model j = ',I2,': ',A)") j, TRIM(GRD_pathfile(j)(1:50)) WRITE (3, "(' Scaled version of forecast model j = ',I2,': ',A)") j, TRIM(GRD_pathfile(j)(1:50)) IF (using_mask) THEN WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I10,' unmasked cells...')") m, unmasked_cells WRITE (3, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I10,' unmasked cells...')") m, unmasked_cells ELSE WRITE (*, "(' Sampling ',I5,' EQ #s from Poisson for each of ',I4,' rows and ',I4,' columns...')") m, nGRDrows, nGRDcols WRITE (3, "(' 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 = scaled_forecast(j, 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., 37.443 or 1023.237 IF (NINT(REAL(lambda + 3.0 * SQRT(lambda))) > comb_limit) THEN WRITE (*, "(' ERROR: comb_limit = ',I5,' is too low when lambda = ',F10.2)") comb_limit, lambda WRITE (*, "(' Increase parameter, or reduce test-time-window length so that this does not happen.')") WRITE (3, "(' ERROR: comb_limit = ',I5,' is too low when lambda = ',F10.2)") comb_limit, lambda WRITE (3, "(' Increase parameter, or reduce test-time-window length so that this does not happen.')") CALL Pause() STOP END IF !Compute probabilities of different integer outcomes, given this expectation: probability_of_integer = 0.0D0 ! whole DP vector initialized total_probability = 0.0D0 ! DP scalar initialized create_comb_3: DO omega = 0, comb_limit ! This is the limit of what FAC can handle without over/underflows. !--- IMSL version: -------------------------------------------------------------------------------------- !probability_of_integer(omega) = (lambda ** omega) * DEXP(-lambda) / DFAC(omega) !Using DFAC(n) = n! from IMSL Library. !--- MKL does not have a factorial; using my own: ------------------------------------------------------- probability_of_integer(omega) = (lambda ** omega) * DEXP(-lambda) / Factorial(omega) ! where Factorial() is a pre-computed REAL*8 vector. !-------------------------------------------------------------------------------------------------------- total_probability = total_probability + probability_of_integer(omega) END DO create_comb_3 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 WRITE (3, "(' 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 cumulative_probability(0:comb_limit) = cumulative_probability(0:comb_limit) + lost_probability cumulative_probability(comb_limit) = 1.0D0 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) !--- MKL method: ------------------------------------------------------------------------------------------------------------------ status = vsRngUniform( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = m, r = sim_temporary, a = 0.0, b = 1.0 ) !---------------------------------------------------------------------------------------------------------------------------------- DO k = 1, m ! (simulation #) search_2: DO omega = 0, comb_limit ! which upper limit we do not expect to reach! (Most bins will produce 0, 1, 2, ... ?) IF (sim_temporary(k) <= cumulative_probability(omega)) THEN N_hat = omega sim_N_hat_3D_matrix(k, row, col) = omega ! sampled # of EQs for this bin, simulation, & forecast EXIT search_2 END IF END DO search_2 ! 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) ELSE ! .NOT.cell_in_use(row, col) sim_N_hat_3D_matrix(1:m, row, col) = NINT(lambda) 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 scaled simulations...')") WRITE (3, "(' Totalling virtual earthquakes for each of these scaled simulations...')") sum = 0.0D0 ! just initializing the total # EQs over ALL m simulations of this SCALED forecast: DO k = 1, m N = 0 ! just initializing the # in one particular forecast... 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 scaled_N_hat_matrix(k, j) = N sum = sum + N END DO mean = sum / m variance = 0.0D0 ! just initializing... DO k = 1, m variance = variance + (scaled_N_hat_matrix(k, j) - mean)**2 END DO sigma = DSQRT(variance / (m - 1)) WRITE (*, "(' Mean = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (3, "(' Mean = ',F10.2,' and sigma = 'F10.2)") mean, sigma !Compute log-likelihood scaled_L_hat_matrix for each of the (k = 1, ..., m) simulations of this scaled forecast #j: WRITE (*, "(' Computing log-likelihood for each of these simulations...')") WRITE (3, "(' Computing log-likelihood for each of these simulations...')") DO k = 1, m ! each simulation of forecast #j L_sum = 0.0D0 ! just initializing sum DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN lambda = scaled_forecast(j, row, col) * dt_da_factor(row, col) N_hat = sim_N_hat_3D_matrix(k, row, col) ! convert INT(2) --> INT(4); also, single access into huge array! IF (lambda == 0.0) THEN ! no EQs forecast in this bin IF (N_hat == 0) THEN ! no problem ! p = 1, log(p) = 0; no change to L ELSE L_sum = -999999. ! temporary substitute for -infinity; forecast will be rejected END IF ELSE ! lamba > 0 ; OK to take log(lambda) IF (N_hat > Factorial_limit) THEN WRITE (*, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") N_hat, Factorial_limit WRITE (*, "(' Reduce length of test-time-window so that this does not happen.')") WRITE (3, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") N_hat, Factorial_limit WRITE (3, "(' Reduce length of test-time-window so that this does not happen.')") CALL Pause() STOP END IF !--- IMSL version: ----------------------------------------------------------------------------------------------------- !L = L -lambda + N_hat * LOG(lambda) - LOG(DFAC(N_hat)) !N.B. DFAC() is the factorial function from IMSL; DFAC(n) = n!. !--- MKL does not have factorial; using my own: ------------------------------------------------------------------------ L_sum = L_sum -lambda + N_hat * LOG(lambda) - LOG_Factorial(N_hat) !----------------------------------------------------------------------------------------------------------------------- !N.B. lambda > 0, so LOG(lambda) will not cause any crashes !N.B. 0! == 1, so last term will not cause any crashes. END IF END IF ! cell_in_use(row, col) END DO ! columns END DO ! rows L = L_sum ! R4 = R8 scaled_L_hat_matrix(k, j) = L ! where (k = 1, ..., m) is the simulation #, j is the SCALED forecast # END DO ! k = 1, m sim_temporary(1:m) = scaled_L_hat_matrix(1:m, j) CALL Moments (sim_temporary, m, & ! inputs & mean, sigma) ! outputs WRITE (*, "(' Mean log-likelihood = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (3, "(' Mean log-likelihood = ',F10.2,' and sigma = 'F10.2)") mean, sigma WRITE (*, *) ! provide blank line to separate different forecasts WRITE (3, *) ! provide blank line to separate different forecasts END DO ! j = 1, nForecasts (all are SCALED forecasts in this block of code) IF (ASS_type_index > 0) THEN !Computations for the Area Skill Score (ASS) test: WRITE (*, *) WRITE (3, *) WRITE (*, "(' ------------------------------------------------------------------------')") WRITE (3, "(' ------------------------------------------------------------------------')") WRITE (*, *) WRITE (3, *) WRITE (*, "(' Computing Molchan diagram(s) and Area Skill Score(s)...')") WRITE (3, "(' Computing Molchan diagram(s) and Area Skill Score(s)...')") IF (ASS_type_index == 1) THEN ! spatially-uniform seismicity sum = 0.0D0 DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN d_tau(row, col) = dt_da_factor(row, col) ! note that (uniform) seismicity density factor is not needed; it would just cancel out sum = sum + dt_da_factor(row, col) ELSE d_tau(row, col) = 0.0 END IF END DO END DO d_tau = d_tau / sum ! whole grid; converting to conditional probability for next qualified test EQ ELSE IF (ASS_type_index == 2) THEN ! prior Relative Intensity (RI), from the pre-test-time-window .EQC file !Count number of qualifying (i.e., in-box, unmasked, large, shallow) events in old (pre-test-time-window) catalog (to establish prior relative intensity): ALLOCATE ( old_selected(old_EQC_lines) ) ! LOGICAL: does this EQ count? ALLOCATE ( old_n_in_bin(nGRDrows, nGRDcols) ) CALL Count_Events (old_EQC_lines, old_eq_Elon, old_eq_Nlat, old_eq_depth_int, old_eq_mag, & ! critical parts of catalog (excluding dates & times) & box_lon_min, box_lon_max, box_lat_min, box_lat_max, 0.0, maxDepthKm, & ! selection criteria (maxDepthKm = integer) & using_mask, cell_in_use, grid_lon_min, d_lon, grid_lon_max, grid_lat_min, d_lat, grid_lat_max, nGRDrows, nGRDcols, & ! mask description & N, old_selected) ! results N_RI = N ! save under unique name IF (using_mask) THEN WRITE (*, "(' Prior Relative Intensity is based on ',I6,' earthquakes')") N WRITE (*, "(' in the unmasked forecast area.')") WRITE (3, "(' Prior Relative Intensity is based on ',I6,' earthquakes')") N WRITE (3, "(' in the unmasked forecast area.')") ELSE WRITE (*, "(' Prior Relative Intensity is based on ',I6,' earthquakes')") N WRITE (*, "(' in the forecast area.')") WRITE (3, "(' Prior Relative Intensity is based on ',I6,' earthquakes')") N WRITE (3, "(' in the forecast area.')") END IF CALL Assign_to_Bins (old_EQC_lines, old_selected, old_eq_Elon, old_eq_Nlat, & ! inputs describing catalog & d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max, nGRDrows, nGRDcols, & !inputs describing GRD & using_mask, cell_in_use, & ! inputs describing mask & old_n_in_bin) ! output WRITE (*, "(' Computing Relative Intensity (RI) across (lon, lat) space for')") WRITE (*, "(' pre-test-time-window period.')") WRITE (3, "(' Computing Relative Intensity (RI) across (lon, lat) space for')") WRITE (3, "(' pre-test-time-window period.')") ! Compute per-cell quanta of "cost" (or, "fraction of space-time") d_tau, using RI as the reference model: d_tau = 0.0 ! whole grid sum = 0.0D0 DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN d_tau(row, col) = (1.0 * old_n_in_bin(row, col)) / (1.0 * N_RI) sum = sum + d_tau(row, col) ELSE d_tau(row, col) = 0.0 END IF END DO END DO d_tau = d_tau / sum ! whole grid; converting to conditional probability for next qualified test EQ ELSE IF (ASS_type_index == 3) THEN sum = 0.0D0 DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN d_tau(row, col) = reference_forecast(row, col) * dt_da_factor(row, col) ! converting to EQ counts/forecast-window sum = sum + d_tau(row, col) ELSE d_tau(row, col) = 0.0 END IF END DO END DO d_tau = d_tau / sum ! whole grid; converting to conditional probability for next qualified test EQ DEALLOCATE (reference_forecast) END IF ! ASS_type_index == 1, or == 2 (determines how d_tau is computed), or == 3 (using reference .GRD file) ALLOCATE ( d_P(nGRDrows, nGRDcols) ) ! Conditional probability, per bin, and based on current forecast model, ! that the next qualifying earthquake will occur in this spatial bin. ! Note that this array is re-used by each hypothesis in turn. ALLOCATE ( Molchan(1:2, 0:N_true) ) ! (tau, nu) for upper left corner (2nd subscript = 0) and then for each actual EQ, in order as: (tau, nu) ! Note that this array is re-used by each hypothesis in turn. ! If you want to plot the curve(s), then dump this to a different file for each forecast. ALLOCATE ( ASS_score(nForecasts), ASS_quantile(nForecasts) ) ! only results saved after this block of code is done ncStem = LEN_TRIM(pseudoCSEP_logfile) - 4 ! filename stem for MolchanNN files (omitting ".txt") DO j = 1, nForecasts WRITE (*, "(' Computing Molchan diagram for forecast #',I2,'...')") j !Compute per-cell quanta d_P of conditional probability (for the next qualifying test EQ), based on this forecast: d_P = 0.0 ! whole grid sum = 0.0D0 DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN d_P(row, col) = forecast(j, row, col) * dt_da_factor(row, col) / N_sup_j(j) ! Note that N_sup_j is REAL. sum = sum + d_P(row, col) ELSE d_P(row, col) = 0.0 END IF END DO END DO IF ((sum < 0.998D0).OR.(sum > 1.002D0)) THEN WRITE (*, "(' ERROR: Sum of per-bin conditional probabilities d_P should be 1, but is ',F10.5,' for forecast #',I2)") sum, j WRITE (3, "(' ERROR: Sum of per-bin conditional probabilities d_P should be 1, but is ',F10.5,' for forecast #',I2)") sum, j CALL Pause() STOP END IF d_P = d_P / sum ! normalizing to total conditional probability of 1.0D0 !Create list of tau values (Molchan(1, ...)); one per actual earthquake, without regard to order: Molchan = 0.0 ! just to simplify inspection of values during debugging, if needed. Molchan(2, 0) = 1.0 ! upper left-hand corner always has (tau == 0, nu == 1) by definition !Note: Following two loops on bins, and 3rd loop on EQs in a bin, essentially amount to looping on EQs, n = 1, N_true n_done = 0 ! just initializing count of actual EQs that have been processed DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN IF (actual_n_in_bin(row, col) > 0) THEN ! this array saved from very early processing of actual test EQC ! must compute tau (good for all actual EQs in this bin) tau = 0.0D0 ! just initializing before sum contour_d_P = d_P(row, col) DO ii = 1, nGRDrows DO jj = 1, nGRDcols IF (cell_in_use(ii, jj)) THEN IF (d_P(ii, jj) >= contour_d_P) THEN ! this bin is within the contour tau = tau + d_tau(ii, jj) ! "cost" of including this bin within the contour END IF END IF END DO END DO DO n = 1, actual_n_in_bin(row, col) ! processing each EQ separately !tau has been computed; use it for any and all actual EQs in this bin n_done = n_done + 1 ! should end up == N_true when all loops are completed Molchan(1, n_done) = tau ! R4 = R8 (tau is of order 1. I just used a R8 accumulator to avoid losing tiny contributions.) END DO ! n = 1, actual_n_in_bin(row, col) END IF ! actual_n_in_bin > 0 END IF ! cell_in_use(row, col) END DO ! columns END DO ! rows IF (n_done /= N_true) THEN WRITE (*, "(' ERROR in code logic; n_done = ',I6,', but N_true = ',I6)") n_done, N_true WRITE (3, "(' ERROR in code logic; n_done = ',I6,', but N_true = ',I6)") n_done, N_true CALL Pause() STOP END IF !sort the tau values in Molchan(1, ...) from smallest to largest: (Skipping Molchan(1, 0) == 0.0) DO ii = 1, (N_true - 1) DO jj = (ii + 1), N_true !Note that ii < jj. IF (Molchan(1, ii) > Molchan(1, jj)) THEN !swap them temp = Molchan(1, ii) Molchan(1, ii) = Molchan(1, jj) Molchan(1, jj) = temp END IF END DO END DO !assign nu values (Molchan(2, ...)) to the sorted Molchan points. (Molchan(2, 0) == 1 already.) DO n = 1, N_true nu = (1.0 * N_true - n) / (1.0 * N_true) Molchan(2, n) = nu END DO !************************************************************** ! At this point, Molchan curves are output for plotting. !************************************************************** WRITE (c2a, "(I2)") j IF (c2a(1:1) == ' ') c2a(1:1) = '0' Molchan_pathfile = pseudoCSEP_logfile(1:ncStem) // "_Molchan" // c2a // ".dat" OPEN (UNIT = 9, FILE = TRIM(Molchan_pathfile)) ! unconditional OPEN; overwrites any existing. DO n = 0, N_true WRITE (9, *) Molchan(1, n), Molchan(2, n) END DO CLOSE (9) !Integrate the area above and to the right of the Molchan curve: sum = 0.0D0 ! just initializing DO n = 1, N_true sum = sum + (Molchan(1, n) - Molchan(1, n-1)) * (1.0 - ((Molchan(2, n-1) + Molchan(2, n))/2.0)) END DO IF (Molchan(1, N_true) < 1.0) sum = sum + (1.0 - Molchan(1, N_true)) * 1.0! Add rectangle to right of end point. ASS_score(j) = sum ! result saved for report below !Describe signifigance of this score, compared to expected distribution of scores under the prior RI null hypothesis: variance = 1.0D0 / (12.0D0 * N_true) ! per Zechar & Jordan [2008, GJI] standard_deviation = DSQRT(variance) sigmas = (ASS_score(j) - 0.5000) / standard_deviation !--- IMSL version: ------------------------------------------------------------------------------- !ASS_quantile(j) = ANORDF(sigmas) ! using standard Gaussian distribution function from IMSL; result saved for report below !--- MKL version: -------------------------------------------------------------------------------- VSCdfNorm_in(1) = sigmas CALL VSCdfNorm(I4_one, VSCdfNorm_in, VSCdfNorm_out) ASS_quantile(j) = VSCdfNorm_out(1) !------------------------------------------------------------------------------------------------- END DO ! j = 1, nForecasts, computing ASS-test metrics END IF ! ASS_type_index > 0 !OUTPUT FORMAL RESULTS: WRITE (*, *) WRITE (*, *) WRITE (3, *) WRITE (3, *) WRITE (*, "(' ==============================================================================')") WRITE (*, "(' === Results of (this unofficial version of) CSEP test algorithms: ========')") WRITE (3, "(' ==============================================================================')") WRITE (3, "(' === Results of (this unofficial version of) CSEP test algorithms: ========')") WRITE (*, *) WRITE (3, *) IF (using_mask) THEN WRITE (*, "(' Masking GRD file ',A,' was applied.')") TRIM(mask_GRD_pathfile) WRITE (3, "(' Masking GRD file ',A,' was applied.')") TRIM(mask_GRD_pathfile) ELSE WRITE (*, "(' No masking file was applied.')") WRITE (3, "(' No masking file was applied.')") END IF WRITE (*, *) WRITE (3, *) WRITE (*, "(' L-test for data consistency:')") WRITE (*, "(' ----------------------------')") WRITE (3, "(' L-test for data consistency:')") WRITE (3, "(' ----------------------------')") ALLOCATE ( gamma_sub_q_sup_j(s_in_use, nForecasts) ) DO j = 1, nForecasts DO q = 1, s_in_use ! numerator = 0 ! just initializing, before sum DO k = 1, m ! scanning over all simulations based on this forecast: IF (L_hat_matrix(k, j, j) <= L_tilde_matrix(q, j)) numerator = numerator + 1 END DO denominator = m ! (number of simulations of each forecast) gamma_sub_q_sup_j(q, j) = (1.0 * numerator) / (1.0 * denominator) END DO ! q = 1, s_in_use (modification # of catalog) real_vector_s(1:s_in_use) = gamma_sub_q_sup_j(1:s_in_use, j) CALL Moments (real_vector_s, s_in_use, & ! inputs & mean, sigma) ! outputs IF (mean <= 0.05) THEN ! Note: This limiting value from section 3.3 of Schorlemmer et al. [2010, Pure Appl. Geophys.] WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED.')") j, GRD_pathfile(j)(1:10), mean WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED.')") j, GRD_pathfile(j)(1:10), mean ELSE WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), mean WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), mean END IF END DO ! j = 1, nForecasts (L-tests) WRITE (*, *) WRITE (3, *) WRITE (*, "(' N-test for total number of forecast earthquakes:')") WRITE (*, "(' ------------------------------------------------')") WRITE (3, "(' N-test for total number of forecast earthquakes:')") WRITE (3, "(' ------------------------------------------------')") ALLOCATE ( del_sub_q_sup_j(s_in_use, nForecasts) ) DO j = 1, nForecasts DO q = 1, s_in_use ! numerator = 0 ! just initializing, before sum DO k = 1, m ! scanning over all simulations based on this forecast: IF (N_hat_matrix(k, j) <= N_tilde(q)) numerator = numerator + 1 END DO denominator = m ! (number of simulations of each forecast) del_sub_q_sup_j(q, j) = (1.0 * numerator) / (1.0 * denominator) END DO ! q = 1, s_in_use (modification # of catalog) real_vector_s(1:s_in_use) = del_sub_q_sup_j(1:s_in_use, j) CALL Moments (real_vector_s, s_in_use, & ! inputs & mean, sigma) ! outputs IF (mean < 0.025) THEN ! Note: This limiting value from section 3.2 of Schorlemmer et al. [2010, Pure Appl. Geophys.]. WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED (forecast too many EQs).')") j, GRD_pathfile(j)(1:10), mean WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED (forecast too many EQs).')") j, GRD_pathfile(j)(1:10), mean ELSE IF (mean <= 0.975) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), mean WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), mean ELSE ! mean > 0.975: WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED (forecast too few EQs).')") j, GRD_pathfile(j)(1:10), mean WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED (forecast too few EQs).')") j, GRD_pathfile(j)(1:10), mean END IF END DO ! j = 1, nForecasts (N-tests) WRITE (*, *) WRITE (3, *) WRITE (*, "(' R-test for comparison of forecasts:')") WRITE (*, "(' -----------------------------------')") WRITE (3, "(' R-test for comparison of forecasts:')") WRITE (3, "(' -----------------------------------')") ALLOCATE ( alpha_matrix(nForecasts, nForecasts) ) DO row = 1, nForecasts ! indicates which forecast was basis for simulated catalogs; also called j elsewhere DO col = 1, nForecasts ! indicates which forecast is tested against simulations based on (usually) another. DO q = 1, s_in_use ! numerator = 0 ! just initializing, before sum DO k = 1, m ! scanning over all simulations based on this forecast: IF ((L_hat_matrix(k, row, row) - L_hat_matrix(k, row, col)) <= R_tilde_3D_matrix(q, row, col)) numerator = numerator + 1 END DO denominator = m ! (number of simulations based on each forecast) real_vector_s(q) = (1.0 * numerator) / (1.0 * denominator) END DO ! q = 1, s_in_use (modification # of catalog) CALL Moments (real_vector_s, s_in_use, & ! inputs & mean, sigma) ! outputs alpha_matrix(row, col) = mean END DO ! col = 1, nForecasts (column of alpha matrix) END DO ! row = 1, nForecasts (row of alpha matrix) !output the alpha matrix just computed: WRITE (*, "(' ',' j =',99I10)") (col, col = 1, nForecasts) WRITE (3, "(' ',' j =',99I10)") (col, col = 1, nForecasts) WRITE (*, "(' ',24X,99(1X,A9))") (GRD_pathfile(col)(1:9), col = 1, nForecasts) WRITE (3, "(' ',24X,99(1X,A9))") (GRD_pathfile(col)(1:9), col = 1, nForecasts) DO row = 1, nForecasts any_rejection = .FALSE. ! just initializing, for each row forecast DO col = 1, nForecasts IF (alpha_matrix(row, col) <= 0.05) any_rejection = .TRUE. ! Note: This limiting value from section 3.5 of Schorlemmer et al. [2010, Pure Appl. Geophys.]. END DO IF (any_rejection) THEN comment = "REJECTED." ELSE comment = "SURVIVES." END IF WRITE (*, "(' ','i=',I2,1X,A9,1X,A9,99F10.5)") row, GRD_pathfile(row)(1:9), comment, (alpha_matrix(row, col), col = 1, nForecasts) WRITE (3, "(' ','i=',I2,1X,A9,1X,A9,99F10.5)") row, GRD_pathfile(row)(1:9), comment, (alpha_matrix(row, col), col = 1, nForecasts) END DO WRITE (*, *) WRITE (3, *) WRITE (*, "(' S-test for spatial pattern:')") WRITE (*, "(' ---------------------------')") WRITE (3, "(' S-test for spatial pattern:')") WRITE (3, "(' ---------------------------')") DO j = 1, nForecasts ! (which are all scaled versions in this test) numerator = 0 ! just initializing, before sum DO k = 1, m ! scanning over all simulations based on this forecast: IF (scaled_L_hat_matrix(k, j) <= scaled_L_sup_j(j)) numerator = numerator + 1 END DO denominator = m ! (number of simulations of each scaled forecast) zeta(j) = (1.0 * numerator) / (1.0 * denominator) IF (zeta(j) <= 0.05) THEN ! Note: This limiting value is somewhat unclear in Zechar et al. [2010; BSSA]; ! contrast implied 0.025 in first paragraph of Results section (p. 1188) ! with later explicit statement of criterion 0.05 in caption to Figure 4. ! The idea that a result exactly matching the numerical criterion value fails ! is mentioned in first paragraph of Results (p. 1188), and reinforced in Table 1. ! I have chosen <= 0.05 for failure based on Figure 4, and also to ! enhance the analogies between the S-test and the older L-test. WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED.')") j, GRD_pathfile(j)(1:10), zeta(j) WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = FAILED.')") j, GRD_pathfile(j)(1:10), zeta(j) ELSE WRITE (*, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), zeta(j) WRITE (3, "(' j = ',I2,': ',A10,': ',F10.5,' = PASSED.')") j, GRD_pathfile(j)(1:10), zeta(j) END IF END DO ! j = 1, [ scaled ] nForecasts (S-test results) IF (ASS_type_index > 0) THEN WRITE (*, *) WRITE (3, *) WRITE (*, "(' ASS-test for spatial pattern:')") WRITE (*, "(' -----------------------------')") WRITE (3, "(' ASS-test for spatial pattern:')") WRITE (3, "(' -----------------------------')") WRITE (*, "(' j: Forecast : ASS quantile comment')") WRITE (3, "(' j: Forecast : ASS quantile comment')") DO j = 1, nForecasts ! IF (ASS_quantile(j) <= 0.05) THEN IF (ASS_type_index == 1) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 2) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 3) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = INFERIOR to reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) END IF ELSE IF (ASS_quantile(j) < 0.95) THEN IF (ASS_type_index == 1) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 2) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 3) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = not distinguishable from reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) END IF ELSE ! quantile >= 0.95 IF (ASS_type_index == 1) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to spatially-uniform seismicity.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 2) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to prior Relative Intensity (RI).')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) ELSE IF (ASS_type_index == 3) THEN WRITE (*, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) WRITE (3, "(' j = ',I2,': ',A10,': ',2F10.5,' = SUPERIOR to reference forecast.')") j, GRD_pathfile(j)(1:10), ASS_score(j), ASS_quantile(j) END IF END IF END DO ! j = 1, nForecasts (ASS-test results) END IF ! ASS_type_index > 0 WRITE (*, *) WRITE (3, *) WRITE (*, "(' ==============================================================================')") WRITE (3, "(' ==============================================================================')") !======================================== end OUTPUTS ===================================================== DEALLOCATE ( sim_N_hat_3D_matrix ) ! freeing up several GB which are probably not needed, even for debugging. WRITE (*, *) WRITE (3, *) WRITE (*, "(' This concludes the pseudoCSEP scoring program (for now)...')") WRITE (3, "(' This concludes the pseudoCSEP scoring program (for now)...')") CLOSE (3) ! end the log file CALL Pause() CONTAINS CHARACTER*10 FUNCTION ASCII10(x) ! Returns a right-justified 10-byte (or shorter) version of a ! floating-point number, in "human" format, with no more ! than 4 significant digits. IMPLICIT NONE REAL, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL :: x_log DOUBLE PRECISION :: y IF (x == 0.0) THEN ASCII10=' 0' RETURN ELSE IF (x > 0.0) THEN punt = (x >= 999999999.5).OR.(x < 0.0000100) ELSE ! x < 0.0 punt = (x <= -99999999.5).OR.(x > -0.000100) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,E10.3)') x !consider possible improvements, from left to right: IF (temp10(3:6) == '.000') THEN ! right-shift 4 spaces over it temp20(7:10) = temp10(7:10) temp20(5:6) = temp10(1:2) temp20(1:4) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(5:6) == '00') THEN ! right-shift 2 spaces over it temp20(7:10) = temp10(7:10) temp20(3:6) = temp10(1:4) temp20(1:2) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(6:6) == '0') THEN ! right-shift 1 space over it temp20(7:10) = temp10(7:10) temp20(2:6) = temp10(1:5) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(8:8) == '+') THEN ! right-shift over + sign in exponent temp20(9:10) = temp10(9:10) temp20(2:8) = temp10(1:7) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(9:9) == '0') THEN ! right-shift over leading 0 in exponent temp20(10:10) = temp10(10:10) temp20(2:9) = temp10(1:8) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF ASCII10 = temp10 ELSE ! can represent without exponential notation x_log = LOG10(ABS(x)) zeros = Int_Below(x_log) - 3 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0) y = -y WRITE (temp20,"(F20.9)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp20(19:20) == '01') temp20(19:20) = '00' !Find first important byte from right; change 0 -> ' ' k10 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 20, 12, -1 IF (temp20(j:j) == '0') THEN temp20(j:j) = ' ' ELSE k10 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0) THEN IF (temp20(10:11) == ' .') temp20(10:11) = '0.' ELSE ! x < 0.0 IF (k10 <= 18) THEN IF (temp20(9:11) == ' -.') temp20(9:11) = '-0.' END IF END IF k1 = k10 - 9 ASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION ASCII10 SUBROUTINE Assign_to_Bins (EQC_lines, selected, eq_Elon, eq_Nlat, & ! inputs describing catalog & d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max, nGRDrows, nGRDcols, & ! inputs describing GRD & using_mask, cell_in_use, & ! inputs describing mask & n_in_bin) ! output IMPLICIT NONE INTEGER, INTENT(IN) :: EQC_lines ! # of EQs in catalog (including outside-box, too-deep, too-small, ...) LOGICAL, DIMENSION(:), INTENT(IN) :: selected ! for being in-box, shallow, large-enough, ... REAL, DIMENSION(:), INTENT(IN) :: eq_Elon, eq_Nlat ! but on some CALLS these will be mod_Elon, mod_Nlat REAL, INTENT(IN) :: d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max ! describes outer perimeter of union of GRD cells INTEGER, INTENT(IN) :: nGRDrows, nGRDcols ! describes GRD LOGICAL, INTENT(IN) :: using_mask LOGICAL*1, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: cell_in_use INTEGER, DIMENSION(nGRDrows, nGRDcols), INTENT(OUT) :: n_in_bin REAL :: real_row, real_col INTEGER :: row, col REAL :: Easting, Easting_col1 n_in_bin = 0 ! whole array; just initializing before sum: !Allocate selected EQs to their proper bins: DO iEq = 1, EQC_lines IF (selected(iEq)) THEN CALL LonLat_2_RowCol(eq_Elon(iEq), eq_Nlat(iEq), & ! most-variable inputs & grid_lon_min, d_lon, nGRDcols, box_lon_min, & & grid_lat_max, d_lat, nGRDrows, & ! non-varying inputs that describe the grid & row, col) ! output IF ((row >= 1).AND.(row <= nGRDrows).AND.(col >= 1).AND.(col <= nGRDcols)) THEN IF (using_mask) THEN IF (cell_in_use(row, col)) THEN n_in_bin(row, col) = n_in_bin(row, col) + 1 END IF ELSE ! not using mask; no need to consult cell_in_use = .TRUE. throughout n_in_bin(row, col) = n_in_bin(row, col) + 1 END IF ! using_mask END IF ! (row, col) is in the legal range END IF ! this EQ was pre-selected END DO ! iEq = 1, EQC_lines END SUBROUTINE Assign_to_Bins SUBROUTINE Count_Events (EQC_lines, eq_Elon, eq_Nlat, eq_depth_int, eq_mag, & ! critical parts of catalog (excluding dates & times) & box_lon_min, box_lon_max, box_lat_min, box_lat_max, threshold, maxDepthKm, & ! box selection criteria (maxDepthKm = integer) & using_mask, cell_in_use, grid_lon_min, d_lon, grid_lon_max, grid_lat_min, d_lat, grid_lat_max, nGRDrows, nGRDcols, & ! mask description & N, selected) ! results IMPLICIT NONE INTEGER, INTENT(IN) :: EQC_lines ! number of events in catalog REAL, DIMENSION(:), INTENT(IN) :: eq_Elon, eq_Nlat ! E longitude and N latitude of each EQ, in degrees INTEGER, DIMENSION(:), INTENT(IN) :: eq_depth_int ! EQ depths, in km, rounded to nearest integer REAL, DIMENSION(:), INTENT(IN) :: eq_mag ! EQ magnitudes REAL, INTENT(IN) :: box_lon_min, box_lon_max, box_lat_min, box_lat_max ! outline of rectangle defining the forecast region REAL, INTENT(IN) :: threshold ! minimum magnitude for this event count INTEGER, INTENT(IN) :: maxDepthKm ! maximum depth, in km, rounded to integer, for this event-count LOGICAL, INTENT(IN) :: using_mask LOGICAL*1, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: cell_in_use ! All .TRUE. if .NOT.using_mask REAL, INTENT(IN) :: grid_lon_min, d_lon, grid_lon_max, grid_lat_min, d_lat, grid_lat_max ! referring to forecast grid; for computing (i, j) of cell_in_use(i, j) INTEGER, INTENT(IN) :: nGRDrows, nGRDcols ! see comment above INTEGER, INTENT(OUT) :: N ! number of events satisfying all criteria LOGICAL, DIMENSION(:), INTENT(OUT) :: selected ! flags will be set = T for EQs that are counted toward total N. INTEGER :: col, iEq, row selected = .FALSE. ! default; values will be changed for those EQs which are selected N = 0 ! just initializing this counter DO iEq = 1, EQC_lines IF (eq_mag(iEq) >= threshold) THEN IF (eq_depth_int(iEq) <= maxDepthKm) THEN IF (eq_Nlat(iEq) >= box_lat_min) THEN IF (eq_Nlat(iEq) <= box_lat_max) THEN !NOTE: Being in the box is a necessary condition, but not sufficient. Assigned (row, col) are critical: CALL LonLat_2_RowCol(eq_Elon(iEq), eq_Nlat(iEq), & ! most-variable inputs & grid_lon_min, d_lon, nGRDcols, box_lon_min, & & grid_lat_max, d_lat, nGRDrows, & ! non-varying inputs that describe the grid & row, col) ! output IF ((row >= 1).AND.(row <= nGRDrows).AND.(col >= 1).AND.(col <= nGRDcols)) THEN IF (using_mask) THEN IF (cell_in_use(row, col)) THEN selected(iEq) = .TRUE. N = N + 1 END IF ! this cell is in use ELSE ! not using_mask: no further tests needed selected(iEq) = .TRUE. N = N + 1 END IF ! using_mask, or not END IF ! assigned (row, col) are both in the legal range END IF ! eq_Nlat(iEq) <= box_lat_max END IF ! eq_Nlat(iEq) >= box_lat_min END IF ! eq_depth_int(iEq) <= maxDepthKm END IF ! eq_mag(iEq) >= threshold END DO ! iEq = 1, EQC_lines END SUBROUTINE Count_Events INTEGER FUNCTION Int_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL, INTENT(IN) :: x INTEGER :: i REAL :: y i = INT(x) IF (x >= 0.) THEN Int_Below = i ELSE ! x < 0. y = 1.*i IF (y <= x) THEN Int_Below = i ELSE ! most commonly Int_Below = i - 1 END IF END IF END FUNCTION Int_Below SUBROUTINE Integrate_forecast(seismicity, j, & & nGRDrows, d_lat, grid_lat_max, & & nGRDcols, box_lon_min, grid_lon_min, d_lon, box_lon_max, & & cell_in_use, & & N_per_year) IMPLICIT NONE DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN):: seismicity INTEGER, INTENT(IN) :: j ! superscript identifying the forecast within mega-array seismicity. INTEGER, INTENT(IN) :: nGRDrows, nGRDcols ! dimensions of each forecast grid REAL, INTENT(IN) :: d_lat, grid_lat_max, box_lon_min, grid_lon_min, d_lon, box_lon_max LOGICAL*1, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: cell_in_use DOUBLE PRECISION, INTENT(OUT) :: N_per_year DOUBLE PRECISION :: N_per_day INTEGER :: row, col REAL, PARAMETER :: Pi = 3.1415927 REAL, PARAMETER :: radians_per_degree = 0.0174533 REAL :: bottom_lat, factor, longitude, square_m, strip_area_steradian, top_lat REAL :: R = 6371000. ! radius of Earth, in m DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 INTEGER :: s_per_day = 24 * 60 * 60 N_per_day = 0.0D0 DO row = 1, nGRDrows top_lat = MIN((grid_lat_max + 0.5 * d_lat - (row - 1) * d_lat), 90.0) bottom_lat = MAX((top_lat - d_lat), -90.0) strip_area_steradian = 2 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) & & - 2 * Pi * (unity - DCOS((90.0D0 - top_lat ) * radians_per_degree)) IF (box_lon_min == grid_lon_min) THEN ! redundant extra column (last = first) was included: square_m = strip_area_steradian * (R**2) / (nGRDcols - 1) ELSE ! no redundant column; each column gets full weight square_m = strip_area_steradian * (R**2) * ((box_lon_max - box_lon_min)/360.0) / nGRDcols END IF DO col = 1, nGRDcols !central longitude = grid_lon_min + (col-1) * d_lon IF (box_lon_min == grid_lon_min) THEN ! redundant extra column (last = first) was included: IF ((col == 1).OR.(col == nGRDcols)) THEN factor = 0.5 ELSE factor = 1.0 END IF ELSE ! no redundant column; every column gets full weight: factor = 1.0 END IF IF (cell_in_use(row, col)) THEN ! If .NOT.using_mask, then all values are .TRUE. N_per_day = N_per_day + seismicity(j, row, col) * s_per_day * square_m * factor END IF END DO END DO N_per_year = N_per_day * 365.25 ! mean length of years; 1 in 4 is a leap year. END SUBROUTINE Integrate_forecast SUBROUTINE LogLikelihood (j, nForecasts, forecast, dt_da_factor, & ! inputs describing forecast & nGRDrows, nGRDcols, & & cell_in_use, & ! describes masking (or, all .TRUE. if no mask) & actual_n_in_bin, & ! input with actual counts of qualifying earthquakes & L ) ! output: log-likelihood (using natural logarithms) IMPLICIT NONE INTEGER, INTENT(IN) :: j, nForecasts ! identifies forecast; 1 <= j <= nForecasts DOUBLE PRECISION, DIMENSION(nForecasts, nGRDrows, nGRDcols), INTENT(IN) :: forecast REAL, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: dt_da_factor INTEGER, INTENT(IN) :: nGRDrows, nGRDcols LOGICAL*1, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: cell_in_use INTEGER, DIMENSION(nGRDrows, nGRDcols), INTENT(IN) :: actual_n_in_bin REAL, INTENT(OUT) :: L INTEGER :: row, col LOGICAL, SAVE :: dumped = .FALSE. DOUBLE PRECISION :: L_sum, lambda, term !NOTE that the global vector LOG_Factorial is used; this must be pre-computed! Also global Factorial_limit is used. L_sum = 0.0D0 ! just initializing sum (which will include lots of tiny entries from quiet cells, and some order-1 entries from active ones. DO row = 1, nGRDrows DO col = 1, nGRDcols IF (cell_in_use(row, col)) THEN ! note that all values are .TRUE. if no mask is in place lambda = forecast(j, row, col) * dt_da_factor(row, col) ! expected earthquake count in this cell (R8 = R8 * R4) IF (lambda == 0.0D0) THEN ! no EQs forecast in this bin IF (actual_n_in_bin(row, col) == 0) THEN ! no problem ! p = 1, log(p) = 0; no change to L ELSE L_sum = -999999.D0 ! temporary substitute for -infinity; forecast will be rejected END IF ELSE ! lamba > 0 ; OK to take log(lambda) IF (actual_n_in_bin(row, col) > Factorial_limit) THEN WRITE (*, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") actual_n_in_bin(row, col), Factorial_limit WRITE (*, "(' Reduce length of test-time-window so that this does not happen.')") WRITE (3, "(' ERROR: ',I6,' EQs in one spatial cell exceeds limit of ', I6)") actual_n_in_bin(row, col), Factorial_limit WRITE (3, "(' Reduce length of test-time-window so that this does not happen.')") CALL Pause() STOP END IF !--- IMSL version: --------------------------------------------------------------------------------- !L = L -lambda + n_in_bin(row, col) * LOG(lambda) - LOG(DFAC(n_in_bin(row, col))) !N.B. DFAC() is the factorial function from IMSL; DFAC(n) = n! !--- MKL does not have a factorial function; using my own: ---------------------------------------- term = -lambda + actual_n_in_bin(row, col) * DLOG(lambda) - LOG_Factorial(actual_n_in_bin(row, col)) L_sum = L_sum + term !-------------------------------------------------------------------------------------------------- !N.B. lambda > 0, so LOG(lambda) will not cause any crashes !N.B. 0! == 1, so last term will not cause any crashes. END IF END IF ! cell_in_use(row, col) END DO END DO L = L_sum ! R*4 = R*8 (OK, because typical values are something like -200 to -1200). ! We just used R8 during the lengthy sum to avoid losing the many tiny contributions from empty cells. END SUBROUTINE LogLikelihood SUBROUTINE LonLat_2_RowCol(lon, lat, & ! most-variable inputs: EQ epicenter or epicentroid coordinates (degrees_East, degrees_North). & grid_lon_min, d_lon, nGRDcols, box_lon_min, & & grid_lat_max, d_lat, nGRDrows, & ! non-varying inputs that describe the grid & row, col) ! outputs !Caution: Outputs (row, col) are not guaranteed to be in legal range; calling program must check this. IMPLICIT NONE REAL, INTENT(IN) :: lon, lat REAL, INTENT(IN) :: grid_lon_min, d_lon, box_lon_min, grid_lat_max, d_lat INTEGER, INTENT(IN) :: nGRDrows, nGRDcols INTEGER, INTENT(OUT) :: row, col REAL*8 :: row_R8, col_R8, Easting, Easting_col1 REAL*8, PARAMETER :: quantum = 0.0001D0 ! small perturbation to all EQ locations, to avoid falling "exactly" on cell boundaries. !(because in finite-precision numerical computing, such locations are not EXACTLY on the boundary, ! and fall to one side or the other, in ways that are hard to understand and debug!) row_R8 = 1.0D0 + quantum + (DBLE(grid_lat_max) - DBLE(lat)) / DBLE(d_lat) ! avoiding truncation of rational fractions row = IDNINT(row_R8) ! NOTE that this gives "arbitrary" result when row_R8 == (n + 0.5D0); however, rounding of output will at least be CONSISTENT. Easting = DBLE(lon) - DBLE(box_lon_min) IF (Easting >= 360.0D0) Easting = Easting - 360.0D0 IF (Easting >= 360.0D0) Easting = Easting - 360.0D0 IF (Easting < 0.0D0) Easting = Easting + 360.0D0 IF (Easting < 0.0D0) Easting = Easting + 360.0D0 !Note: Easting should now be in range from 0. to 359.99, measured Eward, relative to W side of (lon,lat) box. Easting_col1 = DBLE(grid_lon_min) - DBLE(box_lon_min) ! may be 0.0D0, or DBLE(d_lon)/2.0D0 col_R8 = 1.0D0 + quantum + (Easting - Easting_col1) / DBLE(d_lon) col = IDNINT(col_R8) ! NOTE that this gives "arbitrary" result when col_R8 == (n + 0.5D0); however, rounding of output will at least be CONSISTENT. END SUBROUTINE LonLat_2_RowCol SUBROUTINE Moments (vector, length, & ! inputs & mean, sigma) ! outputs IMPLICIT NONE INTEGER, INTENT(IN) :: length REAL, DIMENSION(length), INTENT(IN) :: vector REAL, INTENT(OUT) :: mean, sigma DOUBLE PRECISION :: sum, variance INTEGER :: i sum = 0.0D0 DO i = 1, length sum = sum + vector(i) END DO mean = sum / length variance = 0.0D0 DO i = 1, length variance = variance + (vector(i) - mean)**2 END DO IF (length > 1) THEN sigma = DSQRT(variance / (length - 1)) ELSE sigma = 0.0 END IF END SUBROUTINE Moments SUBROUTINE Pause () IMPLICIT NONE CHARACTER*1 c1 WRITE (*,"(' Press [Enter] to continue...'\)") READ (*,"(A)") c1 END SUBROUTINE Pause SUBROUTINE Prompt_for_Integer (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with an integer value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! This also occurs IF (mt_flashby), without waiting for user. ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text INTEGER, INTENT(IN) :: default INTEGER, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, trial, written LOGICAL :: finished WRITE (suggested,"(I11)") default suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an integer of 9 digits or less.')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE Prompt_for_Integer SUBROUTINE Prompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! The same happens IF (mt_flashby), without waiting for the user. ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE Prompt_for_Logical SUBROUTINE Prompt_for_Real (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a real value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text REAL, INTENT(IN) :: default REAL, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use ASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = ASCII10(default) !(2) Be sure that the number contains some sign that it is floating-point, not integer: IF (INDEX(suggested, '.') == 0) THEN IF ((INDEX(suggested, 'E') == 0).AND.(INDEX(suggested, 'e') == 0).AND. & & (INDEX(suggested, 'D') == 0).AND.(INDEX(suggested, 'd') == 0)) THEN suggested = ADJUSTL(suggested) point = LEN_TRIM(suggested) + 1 suggested(point:point) = '.' END IF END IF !------------------------------------------------------------------------------------ suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an real number using 11 characters (or less).')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE Prompt_for_Real SUBROUTINE Prompt_for_String (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a character-string value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text CHARACTER*(*), INTENT(IN) :: default CHARACTER*(*), INTENT(OUT) :: answer CHARACTER*80 trial INTEGER :: blank_at, default_bytes, leftover, & & prompt_bytes, written prompt_bytes = LEN_TRIM(prompt_text) default_bytes = LEN_TRIM(default) written = 0 leftover = 79 - prompt_bytes - 4 ! unless changed below DO WHILE ((prompt_bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) leftover = 79 - (blank_at - (written+1) + 1) - 4 written = blank_at END DO IF (leftover >= default_bytes) THEN WRITE (*,"(' ',A,' [',A,']')") prompt_text(written+1:prompt_bytes), TRIM(default) ELSE WRITE (*,"(' ',A)") prompt_text(written+1:prompt_bytes) WRITE (*,"(' [',A,']')") TRIM(default) END IF WRITE (*,"(' ?: '\)") READ (*,"(A)") trial IF (LEN_TRIM(trial) == 0) THEN answer = TRIM(default) ELSE answer = TRIM(trial) END IF END SUBROUTINE Prompt_for_String SUBROUTINE Sort (length, & ! input & vector) ! INTENT(INOUT) !Sorts from lowest (#1) to highest (#length) real value. IMPLICIT NONE INTEGER, INTENT(IN) :: length REAL, DIMENSION(length), INTENT(INOUT) :: vector INTEGER :: i, j REAL :: temp DO i = 1, (length - 1) DO j = (i + 1), length IF (vector(j) < vector(i)) THEN ! swap temp = vector(i) vector(i) = vector(j) vector(j) = temp END IF END DO END DO END SUBROUTINE Sort END PROGRAM pseudoCSEP