PROGRAM Kagan_2009_GJI_I_scores ! A program for scoring earthquake-forecast grids by the methods of: ! ! Kagan, Yan Y. [2009] Testing long-term earthquake forecasts: ! likelihood methods and error diagrams, ! Geophys. J. Int., v. 177, pages 532-542. ! ! Some advantages of these methods are that they: ! -are insensitive to the grid used to cover the Earth; ! -are insensitive to changes in the overall seismicity rate; ! -do not modify the locations or magnitudes of test earthquakes; ! -do not require simulation of virtual catalogs; ! -return relative quality measures, not just "pass" or "fail;" and ! -indicate relative specificity of forecasts as well as relative success. ! ! This program by Peter Bird, UCLA, May 2013. Masking added October 2014. ! 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). ! 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 at or exceeding threshold) presented at gridded ! longitude/latitude points, 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.) ! ! *Optionally, a "mask.grd" file (in Peter Bird's .GRD format) ! can be read to prevent the use of certain forecast cells in ! computations of both I_0 and I_1. An entry of "1" in this ! .grd file indicates masking; an entry of "0" indicates normal ! processing of that cell. IMPLICIT NONE CHARACTER*1 :: new_eq_tenths CHARACTER*2 :: new_eq_month, new_eq_day, new_eq_hour, new_eq_minute, & & new_eq_second CHARACTER*132 :: EQC_file, forecast_GRD_file, log_file, mask_file, text_line INTEGER :: all_cells, i, i_mask, ios, j, j_mask, mask_nRows, mask_nColumns, masked_cells, max_depth_km_int, N, nRows, nColumns, & & new_eq_year, new_eq_depth_int, unmasked_cells INTEGER*1, DIMENSION(:,:), ALLOCATABLE :: mask_grd_I1 ! (when using_mask = .TRUE.) LOGICAL :: global_mask, report_scores, use_this_EQ, using_mask LOGICAL*1, DIMENSION(:,:), ALLOCATABLE :: cell_in_use ! (when using_mask = .TRUE.) REAL*8, PARAMETER :: Pi = 3.141592653589793D0 REAL*8, PARAMETER :: radians_per_degree = 0.01745329252D0 REAL*8, PARAMETER :: R = 6371000.D0 ! Earth radius in m; not really necessary. REAL*8 :: bottom_lat, box_lat_max, box_lat_min, box_lon_max, box_lon_min, & & d_lat, d_lon, & & Easting, & & factor, & & grid_center_lon, grid_lat_max, grid_lat_min, grid_lon_max, grid_lon_min, & & lat, lon, & & mask_grid_lat_min, mask_d_lat, mask_grid_lat_max, & & mask_grid_lon_min, mask_d_lon, mask_grid_lon_max, & & new_eq_Elon, new_eq_Nlat, new_eq_mag, & & strip_area_steradian, square_m, & & test_lon, threshold, this_EQ_score, top_lat, & & uniform_forecast REAL*8, DIMENSION(:, :), ALLOCATABLE :: area, forecast, nu_events DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 ! N.B. "DOUBLE PRECISION" variables were already REAL*8 when the rest of the program was REAL*4. DOUBLE PRECISION :: I_0, I_1, total_area, total_events WRITE (*, *) WRITE (*, "(' PROGRAM Kagan_2009_GJI_I_scores')") WRITE (*, *) WRITE (*, "(' A program for scoring earthquake-forecast grids by the methods of:')") WRITE (*, *) WRITE (*, "(' Kagan, Yan Y. [2009] Testing long-term earthquake forecasts:')") WRITE (*, "(' likelihood methods and error diagrams,')") WRITE (*, "(' Geophys. J. Int., v. 177, pages 532-542.')") WRITE (*, *) WRITE (*, "(' Some advantages of these methods are that they:')") WRITE (*, "(' -are insensitive to the grid used to cover the Earth;')") WRITE (*, "(' -are insensitive to changes in the overall seismicity rate;')") WRITE (*, "(' -do not modify the locations or magnitudes of test earthquakes;')") WRITE (*, "(' -do not require simulation of virtual catalogs;')") WRITE (*, "(' -return relative quality measures, not just ""pass"" or ""fail;"" and')") WRITE (*, "(' -indicate relative specificity of forecasts as well as relative success.')") WRITE (*, *) WRITE (*, "(' This program by Peter Bird, UCLA, May 2013 & March 2014.' & & /' Masking added October 2014.' & & /' Converted to REAL*8 and 64-bit in July 2015.')") WRITE (*, *) WRITE (*, "(' -------------------------------------------------------------------------')") WRITE (*, *) log_file = ' ' CALL DPrompt_for_String('Enter name for (new) log file of this run:', log_file, log_file) OPEN (UNIT = 21, FILE = TRIM(log_file)) ! Absolute OPEN; overwrites any existing file(?) WRITE (21, "('PROGRAM Kagan_2009_GJI_I_scores')") WRITE (21, *) WRITE (21, "('A program for scoring earthquake-forecast grids by the methods of:')") WRITE (21, *) WRITE (21, "('Kagan, Yan Y. [2009] Testing long-term earthquake forecasts:')") WRITE (21, "(' likelihood methods and error diagrams,')") WRITE (21, "(' Geophys. J. Int., v. 177, pages 532-542.')") WRITE (21, *) WRITE (21, "('Some advantages of these methods are that they:')") WRITE (21, "('-are insensitive to the grid used to cover the Earth;')") WRITE (21, "('-are insensitive to changes in the overall seismicity rate;')") WRITE (21, "('-do not modify the locations or magnitudes of test earthquakes;')") WRITE (21, "('-do not require simulation of virtual catalogs;')") WRITE (21, "('-return relative quality measures, not just ""pass"" or ""fail;"" and')") WRITE (21, "('-indicate relative specificity of forecasts as well as relative success.')") WRITE (21, *) WRITE (21, "('This program by Peter Bird, UCLA, May 2013 & March 2014.' & & /'Masking added October 2014.' & & /'Converted to REAL*8 and 64-bit in July 2015.')") WRITE (21, *) WRITE (21, "('-------------------------------------------------------------------------')") 10 WRITE (*, *) WRITE (*, "(' REQUIRED INPUT: AT LEAST ONE FORECAST:')") WRITE (*, "(' *A stationary forecast of seismicity rate with a single magnitude')") WRITE (*, "(' bin (everything at or 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 (21, *) WRITE (21, "('REQUIRED INPUT: AT LEAST ONE FORECAST:')") WRITE (21, "('*A stationary forecast of seismicity rate with a single magnitude')") WRITE (21, "(' bin (everything at or exceeding threshold) presented at gridded')") WRITE (21, "(' longitude/latitude points, according to the Long_Term_Seismicity')") WRITE (21, "(' variant of the .grd (Gridded Data) format documented at:')") WRITE (21, "(' http://peterbird.name/guide/grd_format.htm')") WRITE (21, "(' (NOTE that the Long_Term_Seismicity variant of the .grd format')") WRITE (21, "(' inserts the threshold magnitude into the first line, and also')") WRITE (21, "(' adds 2 lines with the forecast integral after the end of the')") WRITE (21, "(' long number list.)')") WRITE (*, *) forecast_GRD_file = ' ' CALL DPrompt_for_String('Enter name of forecast .GRD file: ', forecast_GRD_file, forecast_GRD_file) WRITE (21, *) WRITE (21, "('Enter name of forecast .GRD file: ')") WRITE (21, "(A)") TRIM(forecast_GRD_file) OPEN (UNIT = 1, FILE = TRIM(forecast_GRD_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found in current folder/directory. Try again...')") CALL Pause() GO TO 10 END IF WRITE (*, *) WRITE (21, *) WRITE (*, "(' Here are the first 2 lines of this file:')") WRITE (21, "('Here are the first 2 lines of this file:')") READ (1, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (21, "(A)") TRIM(text_line) READ (1, "(A)") text_line WRITE (*, "(' ',A)") TRIM(text_line) WRITE (21, "(A)") TRIM(text_line) WRITE (*, *) WRITE (21, *) threshold = 5.767D0 CALL DPrompt_for_Real('What is the threshold magnitude of the forecast?', threshold, threshold) WRITE (21, "('What is the threshold magnitude of the forecast?')") WRITE (21, *) threshold CLOSE (1) WRITE (*, *) max_depth_km_int = 70 CALL DPrompt_for_Integer('What is the maximum depth (in km, using INTEGER) for this forecast?', max_depth_km_int, max_depth_km_int) WRITE (21, *) WRITE (21, "('What is the maximum depth (in km, using INTEGER) for this forecast?')") WRITE (21, *) max_depth_km_int OPEN (UNIT = 1, FILE = TRIM(forecast_GRD_file), STATUS = "OLD", IOSTAT = ios) READ (1, *) grid_lon_min, d_lon, grid_lon_max grid_center_lon = (grid_lon_min + grid_lon_max) / 2.0D0 READ (1, *) grid_lat_min, d_lat, grid_lat_max nRows = NINT(1 + (grid_lat_max - grid_lat_min) / d_lat) nColumns = NINT(1 + (grid_lon_max - grid_lon_min) / d_lon) ALLOCATE ( forecast(nRows, nColumns ) ) forecast = 0.0D0 ! whole matrix ALLOCATE ( area(nRows, nColumns) ) area = 0.0D0 ! whole matrix ALLOCATE ( nu_events(nRows, nColumns) ) nu_events = 0.0D0 ! whole matrix READ (1, *)((forecast(i, j), j = 1, nColumns), i = 1, nRows) CLOSE (1) ! forecast_GRD file !Consider possible use of a MASK.GRD file(?) 30 WRITE (*, *) WRITE (*, "(' *Optionally, a ""mask.grd"" file (in Peter Bird''s .GRD format)')") WRITE (*, "(' can be read to block the use of certain forecast cells in')") WRITE (*, "(' computations of both I_0 and I_1. An entry of ""1"" in this')") WRITE (*, "(' .grd file indicates masking; an entry of ""0"" indicates normal')") WRITE (*, "(' processing of that cell.')") WRITE (21, *) WRITE (21, "('*Optionally, a ""mask.grd"" file (in Peter Bird''s .GRD format)')") WRITE (21, "(' can be read to block the use of certain forecast cells in')") WRITE (21, "(' computations of both I_0 and I_1. An entry of ""1"" in this')") WRITE (21, "(' .grd file indicates masking; an entry of ""0"" indicates normal')") WRITE (21, "(' processing of that cell.')") CALL DPrompt_for_Logical('Do you wish to read in a MASK.GRD file now?:', .FALSE., using_mask) WRITE (21, "('Do you wish to read in a MASK.GRD file now?:', L7)") using_mask IF (using_mask) THEN mask_file = ' ' CALL DPrompt_for_String('Name of MASK.GRD file: ', mask_file, mask_file) WRITE (21, "('Name of MASK.GRD file:', A)") TRIM(mask_file) OPEN (UNIT = 3, FILE = TRIM(mask_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Mask file ',A,' not found (in active folder).')") TRIM(mask_file) WRITE (21, "('ERROR: Mask file ',A,' not found (in active folder).')") TRIM(mask_file) CALL Pause() GO TO 30 END IF READ (3, *) mask_grid_lon_min, mask_d_lon, mask_grid_lon_max READ (3, *) 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.9D0) IF (global_mask) THEN WRITE (*, "(' This mask has global longitude range, so longitude cycling is permitted.')") WRITE (21, "('This mask has global longitude range, so longitude cycling is permitted.')") ELSE WRITE (*, "(' This mask has only a local longitude range, so queries must stay in-bounds.')") WRITE (21, "('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 (3, *, 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 (21, "('ERROR in reading MASK.GRD file; check format and size.')") CALL Pause() STOP ELSE WRITE (*, "(' Mask file ',A,' was read.')") TRIM(mask_file) WRITE (21, "('Mask file ',A,' was read.')") TRIM(mask_file) CALL Pause() END IF CLOSE (3) !Translate between grids of different sizes, from (mask_nRows, mask_nColumns) to (nRows, nColumns): ALLOCATE ( cell_in_use(nRows, nColumns) ) DO i = 1, nRows lat = grid_lat_max - (i-1) * d_lat i_mask = 1 + NINT((mask_grid_lat_max - lat) / mask_d_lat) lat = MIN(MAX(lat, -90.0D0), 90.0D0) ! preventing silly out-of-range abends at 90.0000000001 DO j = 1, nColumns lon = grid_lon_min + (j-1) * d_lon j_mask = 1 + NINT((lon - mask_grid_lon_min) / mask_d_lon) IF ((i_mask < 1).OR.(i_mask > mask_nRows)) THEN WRITE (*, "(' ERROR: Attempt to read mask_grid_I1(',I6,', ',I6,') is out-of-range!')") i_mask, j_mask WRITE (21, "('ERROR: Attempt to read mask_grid_I1(',I6,', ',I6,') is out-of-range!')") i_mask, j_mask CALL Pause() STOP END IF IF (global_mask) THEN IF (j_mask < 1) THEN test_lon = lon + 360.0D0 j_mask = 1 + NINT((test_lon - mask_grid_lon_min) / mask_d_lon) END IF IF (j_mask > mask_nColumns) THEN test_lon = lon - 360.0D0 j_mask = 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 (21, "('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 = nRows * nColumns unmasked_cells = 0 DO i = 1, nRows DO j = 1, nColumns IF (cell_in_use(i, j)) unmasked_cells = unmasked_cells + 1 END DO END DO masked_cells = all_cells - unmasked_cells WRITE (*, "(' Referring to the seismicity-forecast grid,')") WRITE (21, "('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 (21, "('Out of ',I10,' total cells, ',I10,' were masked and ',I10,' were not.')") all_cells, masked_cells, unmasked_cells CALL Pause() END IF ! using_mask !Decide limits to forecast box in (lon, lat) space: IF (ABS((grid_lon_max - grid_lon_min) - 360.0D0) < 0.01D0) THEN ! grid includes redundant extra column (last = first) box_lon_min = grid_lon_min box_lon_max = grid_lon_max ELSE ! allow for half-width of cells in first and last column: box_lon_min = grid_lon_min - 0.5D0 * d_lon box_lon_max = grid_lon_max + 0.5D0 * d_lon END IF IF (grid_lat_max == 90.0D0) THEN ! one cannot go beyond the pole box_lat_max = 90.0D0 ELSE box_lat_max = grid_lat_max + 0.5D0 * d_lat END IF IF (grid_lat_min == -90.0D0) THEN ! one cannot go beyond the pole box_lat_min = -90.0D0 ELSE box_lat_min = grid_lat_min - 0.5D0 * d_lat END IF !Compute cell areas, in m^2: DO i = 1, nRows ! N to S: top_lat = MIN((grid_lat_max + 0.5D0 * d_lat - (i - 1) * d_lat), 90.0D0) bottom_lat = MAX((top_lat - d_lat), -90.0D0) strip_area_steradian = 2.0D0 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) & & - 2.0D0 * 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) / (nColumns - 1) ELSE ! no redundant column; each column gets full weight square_m = strip_area_steradian * (R**2) * ((box_lon_max - box_lon_min)/360.0D0) / nColumns END IF DO j = 1, nColumns IF (box_lon_min == grid_lon_min) THEN ! redundant extra column (last = first) was included: IF ((j == 1).OR.(j == nColumns)) THEN factor = 0.5D0 ELSE factor = 1.0D0 END IF ELSE ! no redundant column; every column gets full weight: factor = 1.0D0 END IF area(i, j) = factor * square_m END DO END DO !Next, sum areas and normalize them all: IF (using_mask) THEN ! complex case: total_area = 0.0D0 DO i = 1, nRows DO j = 1, nColumns IF (cell_in_use(i, j)) total_area = total_area + area(i, j) END DO END DO IF (total_area > 0.0D0) THEN area = area / total_area ! whole matrix; note that sum of values == 1.0D0 inside mask, but more than 1.0D0 for whole matrix. ELSE WRITE (*, "(' ERROR: Mask covers entire grid; nothing left to score!')") WRITE (21, "('ERROR: Mask covers entire grid; nothing left to score!')") CALL Pause() STOP END IF ELSE ! simple case total_area = 0.0D0 DO i = 1, nRows DO j = 1, nColumns total_area = total_area + area(i, j) END DO END DO area = area / total_area ! whole matrix END IF !Compute normalized numbers of events per cell (nu), !which is the same as probability that the next earthquake falls in that cell: nu_events = forecast * area ! whole array; for one unit of elapsed time IF (using_mask) THEN ! complex case: total_events = 0.0D0 DO i = 1, nRows DO j = 1, nColumns IF (cell_in_use(i, j)) total_events = total_events + nu_events(i, j) END DO END DO IF (total_events > 0.0D0) THEN nu_events = nu_events / total_events ! whole array; now represents fraction of total forecast events that are not masked ELSE WRITE (*, "(' ERROR: No events are forecast within the un-masked portion of the forecast.')") WRITE (21, "('ERROR: No events are forecast within the un-masked portion of the forecast.')") CALL Pause() STOP END IF ELSE ! simple case: total_events = 0.0D0 DO i = 1, nRows DO j = 1, nColumns total_events = total_events + nu_events(i, j) END DO END DO nu_events = nu_events / total_events ! whole array; now represents fraction of total forecast events END IF WRITE (*, *) WRITE (*, "(' Computing forecast specificity (I_0)...')") !COMPUTATION OF SPECIFICITY (I_0): IF (using_mask) THEN ! complex case: I_0 = 0.0D0 ! initializing before sum DO i = 1, nRows DO j = 1, nColumns IF (cell_in_use(i, j)) I_0 = I_0 + nu_events(i, j) * DLog2(nu_events(i, j) / area(i, j)) END DO END DO ELSE ! simple case I_0 = 0.0D0 ! initializing before sum DO i = 1, nRows DO j = 1, nColumns I_0 = I_0 + nu_events(i, j) * DLog2(nu_events(i, j) / area(i, j)) END DO END DO END IF WRITE (*, *) WRITE (*, "(' -------------------------------------------------')") WRITE (*, "(' I_0 = ',F10.4,' (specificity of forecast)')") I_0 WRITE (*, "(' -------------------------------------------------')") WRITE (21, *) WRITE (21, "('-------------------------------------------------')") WRITE (21, "('I_0 = ',F10.4,' (specificity of forecast)')") I_0 WRITE (21, "('-------------------------------------------------')") CALL Pause() WRITE (21, *) WRITE (21, "('REQUIRED INPUT: EARTHQUAKE CATALOG for the test time-window.')") WRITE (21, "('*A seismic catalog for (exactly) the test time-window, which is')") WRITE (21, "(' either global or spatially larger than the forecast region;')") WRITE (21, "(' which extends at least as deep as the depth limit of forecast(s);')") WRITE (21, "(' which is complete down to the threshold magnitude of the')") WRITE (21, "(' forecast(s).')") WRITE (21, "(' The catalog must be in Peter Bird''s .eqc = EarthQuake Catalog')") WRITE (21, "(' format (as output by program Seismicity) which is documented in:')") WRITE (21, "(' Bird, P., and Y. Y. Kagan [2004] Plate-tectonic analysis')") WRITE (21, "(' of shallow seismicity: Apparent boundary width, beta,')") WRITE (21, "(' corner magnitude, coupled lithosphere thickness, and')") WRITE (21, "(' coupling in seven tectonic settings,')") WRITE (21, "(' Bull. Seismol. Soc. Am., 94(6), 2380-2399,')") WRITE (21, "(' plus electronic supplement (eqc_format.pdf).')") WRITE (21, "(' NOTE that the .EQC catalog file should be trimmed to the')") WRITE (21, "(' exact time window of the intended tests.')") WRITE (21, "(' NOTE that the optional focal-mechanism information at the')") WRITE (21, "(' end of each line in the .eqc file is not currently used.')") WRITE (21, *) WRITE (21, "('Enter name of EQC file: ')") 110 WRITE (*, *) 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).')") 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 (*, *) EQC_file = ' ' CALL DPrompt_for_String('Enter name of (existing) EQC file: ', EQC_file, EQC_file) OPEN (UNIT = 2, FILE = TRIM(EQC_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found in current folder/directory. Try again...')") CALL Pause() GO TO 110 END IF WRITE (21, "(A)") TRIM(EQC_file) WRITE (*, *) WRITE (21, *) CALL DPrompt_for_Logical("Do you want an output file listing score per test EQ?", .FALSE., report_scores) WRITE (21, "(' Do you want an output file listing score per test EQ?: ',L5)") report_scores IF (report_scores) THEN OPEN (UNIT = 3, FILE = 'KaganIs_per-EQ.dat') WRITE ( *, "(' Writing KaganIs_per-EQ.dat...')") WRITE (21, "(' Writing KaganIs_per-EQ.dat...')") END IF WRITE (*, *) WRITE (21, *) WRITE (*, "(' Computing forecast success (I_1)...')") WRITE (21, "(' Computing forecast success (I_1)...')") !Compute total forecast rate and equivalent uniform epicentroidal rate density: IF (using_mask) THEN ! complex case: total_events = 0.0D0 ! just initializing DO i = 1, nRows DO j = 1, nColumns IF (cell_in_use(i, j)) total_events = total_events + forecast(i, j) * area(i, j) * & & total_area * 3.1557600D7 * 100.0D0 !giving result in epicentroids/century, !where we note that matrix "area" has already been normalized to sum of 1.0 (within the mask), !while total_area represents the part of of the forecast area within the mask, in m^2. END DO END DO ELSE ! simple case total_events = 0.0D0 ! just initializing DO i = 1, nRows DO j = 1, nColumns total_events = total_events + forecast(i, j) * area(i, j) * & & total_area * 3.1557600D7 * 100.0D0 !giving result in epicentroids/century, !where we note that matrix "area" has already been normalized to sum of 1.0, !while total_area represents the entire forecast grid area. END DO END DO END IF uniform_forecast = total_events / (total_area * 3.1557600D7 * 100.0D0) !in units of eqs/m^2/2, as in the forecast GRD file. N = 0 I_1 = 0.0D0 ! initializing, before sum listening: DO READ (2, 120, IOSTAT = ios) & & new_eq_year, new_eq_month, new_eq_day, & & new_eq_hour, new_eq_minute, new_eq_second, new_eq_tenths, & & new_eq_Elon, new_eq_Nlat, & & new_eq_depth_int, new_eq_mag !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 listening IF (new_eq_Nlat >= box_lat_min) THEN IF (new_eq_Nlat <= box_lat_max) THEN Easting = new_eq_Elon - box_lon_min IF (Easting < 0.0D0) Easting = Easting + 360.0D0 IF (Easting < 0.0D0) Easting = Easting + 360.0D0 IF (Easting >= 360.0D0) Easting = Easting - 360.0D0 IF (Easting >= 360.0D0) Easting = Easting - 360.0D0 IF (Easting <= (box_lon_max - box_lon_min)) THEN IF (new_eq_depth_int <= max_depth_km_int) THEN IF (new_eq_mag >= threshold) THEN IF (using_mask) THEN !determine grid coordinates to see if this cell is in use: i = NINT(1 + (grid_lat_max - new_eq_Nlat) / d_lat) IF ((new_eq_Elon - grid_center_lon) > +180.0D0) new_eq_Elon = new_eq_Elon - 360.0D0 IF ((new_eq_Elon - grid_center_lon) > +180.0D0) new_eq_Elon = new_eq_Elon - 360.0D0 IF ((new_eq_Elon - grid_center_lon) < -180.0D0) new_eq_Elon = new_eq_Elon + 360.0D0 IF ((new_eq_Elon - grid_center_lon) < -180.0D0) new_eq_Elon = new_eq_Elon + 360.0D0 j = NINT(1 + (new_eq_Elon - grid_lon_min) / d_lon) j = MAX(1, MIN(j, nColumns)) use_this_EQ = cell_in_use(i, j) ELSE ! accept everything! use_this_EQ = .TRUE. END IF IF (use_this_EQ) THEN !Located an eligable test earthquake! N = N + 1 i = NINT(1 + (grid_lat_max - new_eq_Nlat) / d_lat) IF ((new_eq_Elon - grid_center_lon) > +180.0D0) new_eq_Elon = new_eq_Elon - 360.0D0 IF ((new_eq_Elon - grid_center_lon) > +180.0D0) new_eq_Elon = new_eq_Elon - 360.0D0 IF ((new_eq_Elon - grid_center_lon) < -180.0D0) new_eq_Elon = new_eq_Elon + 360.0D0 IF ((new_eq_Elon - grid_center_lon) < -180.0D0) new_eq_Elon = new_eq_Elon + 360.0D0 j = NINT(1 + (new_eq_Elon - grid_lon_min) / d_lon) j = MAX(1, MIN(j, nColumns)) this_EQ_score = DLog2(forecast(i, j) / uniform_forecast) I_1 = I_1 + DLog2(forecast(i, j) / uniform_forecast) IF (report_scores) THEN WRITE (3, 220) this_EQ_score, & & new_eq_year, new_eq_month, new_eq_day, & & new_eq_hour, new_eq_minute, new_eq_second, new_eq_tenths, & & new_eq_Elon, new_eq_Nlat, & & new_eq_depth_int, new_eq_mag 220 FORMAT (F8.4, 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) END IF ! report_scores END IF ! use_this_EQ END IF ! EQ is at or above threshold END IF ! EQ is within depth limit END IF ! Easting is within limits END IF ! EQ latitude is within second limit END IF ! EQ latitude is within first limit END DO listening CLOSE (2) ! EQC_file IF (report_scores) CLOSE(3) ! per-EQ log file WRITE (*, *) WRITE (*, "(' ',I6,' relevant earthquakes found in test catalog.')") N WRITE (21, *) WRITE (21, "(I6,' relevant earthquakes found in test catalog.')") N I_1 = I_1 / N WRITE (*, *) WRITE (*, "(' -------------------------------------------------')") WRITE (*, "(' I_1 = ',F10.4,' (success of forecast)')") I_1 WRITE (*, "(' -------------------------------------------------')") WRITE (21, *) WRITE (21, "('-------------------------------------------------')") WRITE (21, "('I_1 = ',F10.4,' (success of forecast)')") I_1 WRITE (21, "('-------------------------------------------------')") CALL Pause() WRITE (21, *) WRITE (21, "('Work completed.')") CLOSE (21) ! log_file WRITE (*, *) WRITE (*, "(' Work completed.')") !CALL Pause() CONTAINS CHARACTER*10 FUNCTION DASCII10(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*8, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL*8 :: x_log DOUBLE PRECISION :: y IF (x == 0.0D0) THEN DASCII10=' 0' RETURN ELSE IF (x > 0.0D0) THEN punt = (x >= 999999999.5D0).OR.(x < 0.0000100D0) ELSE ! x < 0.0 punt = (x <= -99999999.5D0).OR.(x > -0.000100D0) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,D10.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 DASCII10 = temp10 ELSE ! can represent without exponential notation x_log = DLOG10(ABS(x)) zeros = DInt_Below(x_log) - 3 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0D0) 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.0D0) 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 DASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION DASCII10 INTEGER FUNCTION DInt_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER :: i REAL*8 :: y i = INT(x) IF (x >= 0.0D0) THEN DInt_Below = i ELSE ! x < 0. y = 1.0D0*i IF (y <= x) THEN DInt_Below = i ELSE ! most commonly DInt_Below = i - 1 END IF END IF END FUNCTION DInt_Below SUBROUTINE Pause () IMPLICIT NONE CHARACTER*1 c1 WRITE (*,"(' Press [Enter] to continue...'\)") READ (*,"(A)") c1 END SUBROUTINE Pause SUBROUTINE DPrompt_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". ! 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 occasional 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 DPrompt_for_Integer SUBROUTINE DPrompt_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". ! 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 DPrompt_for_Logical SUBROUTINE DPrompt_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*8, INTENT(IN) :: default REAL*8, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL*8 :: 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 = DASCII10(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 DPrompt_for_Real SUBROUTINE DPrompt_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 DPrompt_for_String REAL FUNCTION DLog2(x) IMPLICIT NONE REAL*8, INTENT(IN) :: x IF (x > 0.0D0) THEN DLog2 = DLOG10(x) / 0.301029995D0 ELSE WRITE (*, "(' ERROR in DLog2: Log of ', ES12.5,' is undefined.')") x WRITE (*, "(' All cells in the forecast must have positive rate density!')") WRITE (*, "(' (Even cells that may be masked must have a positive value.)')") CALL Pause() STOP END IF END FUNCTION DLog2 END PROGRAM Kagan_2009_GJI_I_scores