PROGRAM SHIFT_GSRM2x !Program to create several versions (x = a, b, c, d, e, f) !of a long-term, indefinite-window global shallow seismicity map !based on the Global Strain Rate Map, version 2.1 (GSRM2.1) !created in 2013-2014 by Corne Kreemer. !The methods of scaling strain-rate to seismicity progress !from simple to complex in the various versions, with the !more complex versions using the set of methods known as !Seismic Hazard Inferred From Tectonics (SHIFT), as presented !by [Bird & Liu, 2007, Seismol. Res. Lett.] and adapted to !global strain-rate maps by Bird et al. [2010, Seismol. Res. Lett.]. !By Peter Bird, UCLA, April 2013-April 2014, for the GEM project. !Output will be in Peter Bird's .GRD format, in units of epicenters/m^2/s. !========================= 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 !======================================================================================================== USE Sphere ! Peter Bird's module of geometry-on-sphere utility programs. IMPLICIT NONE INTEGER, PARAMETER :: last_version = 6 ! Number of model versions programmed to date; ! see user-prompt menu below at "SELECT DESIRED MODEL" ! for the meaning of each version. !---Parameters needed for statistical package Intel MKL/VSL: INTEGER :: status TYPE(VSL_STREAM_STATE) :: my_stream !----------------------------------------------------------------------------------------- !Definitions specific to GSRM2.1 input file "GSRM_average_strain_v2.1.txt" of Kreemer [p.c., 2014], !and to its internal (filled-out with zeros) representation GSRM2_grid: INTEGER :: GSRM2_rows, GSRM2_columns ! to be computed from data below: REAL, PARAMETER:: GSRM2_lat_min = -89.800, GSRM2_lat_max = 89.800, GSRM2_d_lat = 0.200, & & GSRM2_lon_min = -179.875, GSRM2_lon_max = 179.875, GSRM2_d_lon = 0.250, & & GSRM2_central_lon = 0.0 ! in decimal degrees of East longitude & North latitude. REAL, PARAMETER:: strain_rate_units = 1.E-9 / 3.15576E7 ! "units of 10^-9/year" REAL, DIMENSION(:, :, :), ALLOCATABLE :: GSRM2_grid !Note that longitude range in this dataset is from GSRM_lon_min to (GSRM_lon_min + 360.0). !Each strain-rate tensor in this grid (and each epicenter_rate_density point computed from it) !represents the average in a rectangular region centered on the grid point. REAL, DIMENSION(:, :), ALLOCATABLE :: epicenter_rate_density REAL, DIMENSION(:, :), ALLOCATABLE :: mean_G ! mean elastic shear modulus used in seismic-moment-rate calculation INTEGER, DIMENSION(:, :), ALLOCATABLE :: effective_cz_meters ! includes retrospective ! adjustment factor (to correct for wrong dips of ! many dip-slip faults), so not equal to ! from Table 5 of Bird & Kagan [2004]. INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: jT5column ! identifies column(s) of Table 5 of Bird & Kagan [2004] used for coupled-thickness REAL, DIMENSION(:, :, :), ALLOCATABLE :: fractionT5column ! identifies fractions of cell moment rate based on each column of Table 5 (above) !----------------------------------------------------------------------------------------- !Data file underwater.grd (containing either 1 or 0 at each grid point) is intended to !have the same resolution and phase as the GSRM2.1 data grid: INTEGER :: underwater_rows, underwater_columns ! to be computed from data below: REAL :: underwater_lat_min = -89.800, underwater_lat_max = 89.800, underwater_d_lat = 0.200, & & underwater_lon_min = -179.875, underwater_lon_max = 179.875, underwater_d_lon = 0.250, & & underwater_central_lon = 0.0 ! in decimal degrees of East longitude & North latitude. INTEGER(KIND = 2), DIMENSION(:, :), ALLOCATABLE :: underwater_INT_GRD LOGICAL*1, DIMENSION(:, :), ALLOCATABLE :: underwater_TF_GRD !----------------------------------------------------------------------------------------- !PB2002 plate-boundary model: REAL, PARAMETER :: PB2002_ministep_m = 3000.0 INTEGER, PARAMETER :: PB2002_list_max = 100000 ! depends on previous PARAMETER PB2002_ministep_m !and on total length of plate boundaries: 260,487 km. REAL, DIMENSION(:, :), ALLOCATABLE :: PB2002_uvecs REAL, DIMENSION(:), ALLOCATABLE :: PB2002_elevation INTEGER*1, DIMENSION(:), ALLOCATABLE :: PB2002_code !----------------------------------------------------------------------------------------- CHARACTER*1 :: dip_c1, eq_tenths CHARACTER*1, DIMENSION(20) :: alphabet CHARACTER*2 :: eq_month, eq_day, eq_hour, eq_minute, eq_second CHARACTER*3 :: c3 CHARACTER*21 :: appended_data CHARACTER*132 :: DAT_file, EQC_file, first_line, GRD_file, last_line, line, new_EQC_file INTEGER :: best_azimuth, best_code, bins, & & boundary_calibration_count, & & code_now, & & dCol_3_sigma, dRow_3_sigma, & & eq_year, eq_depth_int, & & i, ii, iN, internal_points, intraplate_calibration_count, & & ios, iS, isolation_count, it, iTarget, iteration, iTest, iVersion, iZone, & & j, jBeside, jE, jj, jTarget, jW, & & k, & & line_count, & & nAzim, nColumns, nElevation_m, neighbor_count, new_assigns, nGridPoints, nProblems, & & nRows, nTensorParts, nTZ_columns, nTZ_rows, & & old_percent_done, & & PB2002_list_highest, PB2002_list_i, percent_done, & & support, & & T5_column, top_zone, total_EQ_count_above_catalog_threshold INTEGER, DIMENSION(2) :: T5_column_list INTEGER, DIMENSION(4) :: best_support ! subscript is (positive) zone code INTEGER, DIMENSION(4) :: boundary_calibration_count_in_zone ! subscript is (positive) zone code INTEGER, DIMENSION(4) :: neighbors_in_zone ! subscript is (positive) zone code INTEGER*1, DIMENSION(:, :), ALLOCATABLE :: GSRM2_closest_PB2002_code INTEGER*1, DIMENSION(:, :), ALLOCATABLE :: smoothed_GSRM2_tectonic_zone_I1 INTEGER*1, DIMENSION(:, :), ALLOCATABLE :: snapshot_of_zones INTEGER*1, DIMENSION(:, :), ALLOCATABLE :: tectonic_zones_I1 INTEGER*2, DIMENSION(:), ALLOCATABLE :: PB2002_seaward_azimuth INTEGER*2, DIMENSION(:, :), ALLOCATABLE :: GSRM2_closest_PB2002_seaward_azimuth LOGICAL :: eligable_receiver, NS_this_time, smooth, test1, test2, test3, test4, test5 LOGICAL*1, DIMENSION(:, :), ALLOCATABLE :: GSRM2_raw_SUB REAL, PARAMETER :: s_per_year = 31557600.0 ! 365.25 * 24 * 60* 60 REAL, PARAMETER :: Earth_radius_m = 6371000.0 ! 6371 km, but expressed in m REAL, PARAMETER :: GSRM2_Intraplate_corner = 9.0 !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 REAL :: arc_radians, arc_radian2, azimuth_radians, & & best_elevation, best_velocity_mmpa, biggest_eRate, biggest_n_per_year, & & bottom_lat, boundary_scaling_factor, & & catalog_threshold_magnitude, catalog_years, coefficient, & & distance_m, dx, dy, & & E_Sup, effective_cz_denominator, effective_cz_numerator, eq_Elon, eq_mag, eq_Nlat, & & e1_rate, e1h, e2_rate, e2h, e3_rate, err, EW_fraction, exx, exy, eyy, & & forecast_threshold_magnitude, fraction, & & getLon, GSRM2_Intraplate_corner_Nm, & & integral_converted, intraplate_forecast_over_catalog_factor, & & km_long, & & lat, lat1, lat2, lon, lon1, lon2, & & M_rate_per_m2, mean_G_denominator, mean_G_numerator, ministep_radians, & & N_Sup, NS_fraction, & & percent_boundary_area, percent_boundary_eqs, & & percent_intraplate_area, percent_intraplate_eqs, polar_cap_area_m2, & & R2, & & S_Sup, scalar_strainrate, seismicity, sigma_m, sphere_area_m2, standard, strip_area_steradian, & & test_distance_m, test_radians2, tester, & & this_Sup, three_sigma_m, top_lat, total_cell_seismicity, total_moment_rate_SI, travel_m, & & TZ_d_lat, TZ_d_lon, TZ_lat_max, TZ_lat_min, TZ_lon_max, TZ_lon_min, & & velocity_mmpa, & & W_Sup REAL, DIMENSION(2) :: e1h_list, e2h_list, err_list, velocity_factor REAL, DIMENSION(3) :: mid_uvec, omega_uvec, uvec1, uvec2 REAL, DIMENSION(4) :: boundary_scaling_factor_in_zone REAL, DIMENSION(:), ALLOCATABLE :: cell_area_m2 REAL, DIMENSION(:), ALLOCATABLE :: PB2002_velocity_mmpa REAL, DIMENSION(:, :), ALLOCATABLE :: GSRM2_closest_PB2002_radian2 REAL, DIMENSION(:, :), ALLOCATABLE :: GSRM2_closest_PB2002_elevation REAL, DIMENSION(:, :), ALLOCATABLE :: GSRM2_closest_PB2002_velocity_mmpa REAL, DIMENSION(:, :, :), ALLOCATABLE :: smoothed_GSRM2_grid DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 DOUBLE PRECISION :: catalog_threshold_moment_NM, center_normal_rate, & & forecast_threshold_moment_Nm, dilatation_rate, & & integral, intraplate_area_m2, intraplate_seismicity_SI, & & plate_boundary_area_m2, & & radius_rate, & & second_invariant, & & tentative_boundary_eqs_per_catalog, tentative_boundary_eqs_per_s, & & total_area_m2 DOUBLE PRECISION, DIMENSION(4) :: plate_boundary_area_m2_in_zone DOUBLE PRECISION, DIMENSION(4) :: tentative_boundary_eqs_per_s_in_zone DOUBLE PRECISION, DIMENSION(4) :: tentative_boundary_eqs_per_catalog_in_zone !**************************************************************************************** !Definitions specific to parameters imported from Table 5 of Bird & Kagan [2004]: CHARACTER*6, DIMENSION(10) :: T5_column_header REAL :: CMT_duration_s REAL, DIMENSION(10) :: T5_CMT_pure_events, T5_CMT_pure_event_rate, & & T5_CMT_threshold_moment, & & T5_beta, & & T5_corner_magnitude, T5_corner_moment, & & T5_tGR_model_moment_rate, & & T5_length_km, T5_length_m, & & T5_mean_velocity_mmpa, T5_mean_velocity_mps, & & T5_assumed_dip_degrees, T5_assumed_dip_radians, & & T5_assumed_mu_GPa, T5_assumed_mu_Pa, & & T5_lineIntegral_Nps, & & T5_coupledThickness_km, T5_coupledThickness_m, & & T5_assumed_lithosphere_km, T5_assumed_lithosphere_m, & & T5_coupling !GPBT5 ! Empirical coefficients from selected rows of Table 5 of Bird & Kagan [2004, BSSA]: ! T5_column = 1 2 3 4 5 6 7 8 9 10 ! CRB CTF CCB OSR OSR OTF OTF OTF OCB SUB ! normal other slow medium fast T5_column_header = (/"CRB ","CTF ","CCB ","OSRnor","OSRoth","OTFslo","OTFmed","OTFfas","OCB ","SUB "/) T5_CMT_pure_events = (/ 285.9, 198.5, 259.4, 424.3, 77.0, 398.0, 406.9, 376.6, 117.7, 2052.8/) T5_CMT_threshold_moment = (/ 1.13E17, 3.5E17, 3.5E17, 1.13E17, 1.13E17, 2.0E17, 2.0E17, 2.0E17, 3.5E17, 3.5E17/) T5_beta = (/ 0.65, 0.65, 0.62, 0.92, 0.82, 0.64, 0.65, 0.73, 0.53, 0.64/) T5_corner_magnitude = (/ 7.64, 8.01, 8.46, 5.86, 7.39, 8.14, 6.55, 6.63, 8.04, 9.58/) T5_tGR_model_moment_rate = (/ 1.67E12, 3.8E12, 1.06E13, 6.7E11, 1.9E11, 6.7E12, 9.4E11, 9.0E11, 4.6E12, 2.85E14/) T5_length_km = (/ 18126., 19375., 12516., 61807., 61807., 27220., 10351., 6331., 13236., 38990./) T5_mean_velocity_mmpa = (/ 18.95, 21.54, 18.16, 46.40, 7.58, 20.68, 57.53, 97.11, 19.22, 61.48/) T5_assumed_dip_degrees = (/ 55., 73., 20., 55., 55., 73., 73., 73., 20., 14./) T5_assumed_mu_GPa = (/ 27.7, 27.7, 27.7, 25.7, 25.7, 25.7, 25.7, 25.7, 49., 49./) T5_lineIntegral_Nps = (/ 5.5E8, 4.4E8, 6.0E8, 5.0E9, 4.7E8, 5.2E8, 5.3E8, 5.5E8, 1.2E9, 1.58E10/) T5_coupledThickness_km = (/ 3.0, 8.6, 18., 0.13, 0.40, 13., 1.8, 1.6, 3.8, 18./) T5_assumed_lithosphere_km =(/ 6., 12., 13., 8., 8., 14., 14., 14., 14., 26./) T5_coupling = (/ 0.50, 0.72, 1.00, 0.016, 0.05, 0.93, 0.13, 0.11, 0.27, 0.69/) ! and, some readily-derived quantitites: CMT_duration_s = 25.7474 * s_per_year ! 1977.01.01 through 2002.09.30 <== Do NOT update! Must match Bird & Kagan [2004] calibration study. DO i = 1, 10 T5_CMT_pure_event_rate(i) = T5_CMT_pure_events(i) / CMT_duration_s T5_corner_moment(i) = Moment(T5_corner_magnitude(i)) T5_length_m(i) = T5_length_km(i) * 1000.0 T5_mean_velocity_mps(i) = T5_mean_velocity_mmpa(i) * 0.001 / s_per_year T5_assumed_dip_radians(i) = T5_assumed_dip_degrees(i) * radians_per_degree T5_assumed_mu_Pa(i) = T5_assumed_mu_GPa(i) * 1.0E9 T5_coupledThickness_m(i) = T5_coupledThickness_km(i) * 1000.0 T5_assumed_lithosphere_m(i) = T5_assumed_lithosphere_km(i) * 1000.0 END DO !**************************************************************************************** GSRM2_Intraplate_corner_Nm = Moment(GSRM2_Intraplate_corner) alphabet = (/'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t' /) !----------------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' PROGRAM SHIFT_GSRM2x')") WRITE (*, *) WRITE (*, "(' Program to create several versions (x = a, b, c, d, e, f)')") WRITE (*, "(' of a long-term, indefinite-window global shallow seismicity map')") WRITE (*, "(' based on the Global Strain Rate Map, version 2.1 (GSRM2.1)')") WRITE (*, "(' created in 2013-2014 by Corne Kreemer.')") WRITE (*, "(' ')") WRITE (*, "(' The methods of scaling strain-rate to seismicity progress')") WRITE (*, "(' from simple to complex in the various versions, with the')") WRITE (*, "(' more complex versions using the set of methods known as')") WRITE (*, "(' Seismic Hazard Inferred From Tectonics (SHIFT), as presented')") WRITE (*, "(' by [Bird & Liu, 2007, Seismol. Res. Lett.] and adapted to')") WRITE (*, "(' global strain-rate maps by Bird et al. [2010, Seismol. Res. Lett.].')") WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird, UCLA, April 2013-April 2014, for the GEM project.')") WRITE (*, "(' Output will be in Peter Bird''s .GRD format, in units of epicenters/m^2/s.')") WRITE (*, *) CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' Now, opening log file SHIFT_GSRM2x_log.txt.')") OPEN (UNIT = 21, FILE = "SHIFT_GSRM2x_log.txt") ! absolute OPEN; overwrites any existing WRITE (21, "('PROGRAM SHIFT_GSRM2x')") WRITE (21, *) WRITE (21, "('Program to create several versions (x = a, b, c, d, e, f,)')") WRITE (21, "('of a long-term, indefinite-window global shallow seismicity map')") WRITE (21, "('based on the Global Strain Rate Map, version 2.1 (GSRM2.1)')") WRITE (21, "('created in 2013-2014 by Corne Kreemer.')") WRITE (21, *) WRITE (21, "('The methods of scaling strain-rate to seismicity progress')") WRITE (21, "('from simple to complex in the various versions, with the')") WRITE (21, "('more complex versions using the set of methods known as')") WRITE (21, "('Seismic Hazard Inferred From Tectonics (SHIFT), as presented')") WRITE (21, "('by [Bird & Liu, 2007, Seismol. Res. Lett.] and adapted to')") WRITE (21, "('global strain-rate maps by Bird et al. [2010, Seismol. Res. Lett.].')") WRITE (21, *) WRITE (21, "('By Peter Bird, UCLA, April 2013-March 2014, for the GEM project.')") WRITE (21, "('Output will be in Peter Bird''s .GRD format, in units of epicenters/m^2/s.')") WRITE (21, *) WRITE (21, "('Now, opening log file SHIFT_GSRM2x_log.txt.')") nColumns = NINT(360.0 / GSRM2_d_lon) nRows = 1 + NINT((GSRM2_lat_max - GSRM2_lat_min) / GSRM2_d_lat) ALLOCATE ( GSRM2_grid(nRows, nColumns, 3 ) ) GSRM2_grid = 0.0 ALLOCATE ( epicenter_rate_density(nRows, nColumns) ) epicenter_rate_density = 0.0 ALLOCATE ( effective_cz_meters(nRows, nColumns) ) effective_cz_meters = 0 ! And, zero values that remain unchanged will mark intraplate cells. ALLOCATE ( mean_G(nRows, nColumns) ) mean_G = 0.0 ! which should be replaced for all active cells. ALLOCATE ( jT5column(2, nRows, nColumns) ) jT5column = 0 ! (for the many entries that will never be used) ALLOCATE ( fractionT5column(2, nRows, nColumns) ) fractionT5column = 0.0 ! (for the many entries that will never be used) nGridPoints = nRows * nColumns WRITE (*, *) WRITE (*, "(' Creating global grids of ',I6,' rows * ',I6,' columns = ',I12,' points')") nRows, nColumns, nGridPoints WRITE (*, "(' for input strain rate tensors, output epicentral rate densities, & effective cz.')") WRITE (21, *) WRITE (21, "('Creating global grids of ',I6,' rows * ',I6,' columns = ',I12,' points')") nRows, nColumns, nGridPoints WRITE (21, "('for input strain rate tensors, output epicentral rate densities, & mean cz.')") !Compute areas of grid cells and check the total: ALLOCATE ( cell_area_m2(nRows) ) R2 = Earth_radius_m**2 total_area_m2 = 0.0D0 DO i = 1, nRows ! N to S top_lat = GSRM2_lat_max + 0.5 * GSRM2_d_lat - (i-1) * GSRM2_d_lat bottom_lat = top_lat - GSRM2_d_lat 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)) cell_area_m2(i) = strip_area_steradian * R2 / nColumns total_area_m2 = total_area_m2 + nColumns * cell_area_m2(i) END DO WRITE (*, "(' ')") WRITE (*, "(' Total grid area = ',ES12.5,' square meters,')") total_area_m2 WRITE (21, *) WRITE (21, "('Total grid area = ',ES12.5,' square meters,')") total_area_m2 sphere_area_m2 = 4 * Pi * R2 WRITE (*, "(' compared to ideal sphere area of ',ES12.5,' square meters.')") sphere_area_m2 WRITE (21, "('compared to ideal sphere area of ',ES12.5,' square meters.')") sphere_area_m2 polar_cap_area_m2 = 2 * Pi * R2 * (unity - DCOS(0.5D0 * GSRM2_d_lat * radians_per_degree)) WRITE (*, "(' Difference is due to 2 polar caps of ',ES10.2,' square meters each.')") polar_cap_area_m2 WRITE (21, "('Difference is due to 2 polar caps of ',ES10.2,' square meters each.')") polar_cap_area_m2 !Read in Corne Kreemer's GSRM2 strain rates: OPEN (UNIT = 1, FILE = "GSRM_average_strain_v2.1.txt", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: GSRM_average_strain_v2.1.txt not found. Fix, and try again...')") CALL Pause() STOP END IF DO i = 1, 25 ! lines of headers READ (1, *) ! read (& discard) header END DO scanning: DO READ (1, *, IOSTAT = ios) lat, lon, exx, eyy, exy IF (ios < 0) EXIT scanning ! hit the EOF exx = exx * strain_rate_units eyy = eyy * strain_rate_units exy = exy * strain_rate_units i = 1 + NINT((GSRM2_lat_max - lat) / GSRM2_d_lat) j = 1 + NINT((lon - GSRM2_lon_min) / GSRM2_d_lon) GSRM2_grid(i, j, 1) = exx ! = e_E-W GSRM2_grid(i, j, 2) = eyy ! = e_N-S GSRM2_grid(i, j, 3) = exy ! = e_NE-SW END DO scanning CLOSE (1) ! GSRM_average_strain.txt, from Corne Kreemer !Based on original GSM data (which is unchanged if iVersion == 1), !Compute total areas of straining plate boundaries, and of intraplate regions: plate_boundary_area_m2 = 0.0D0 intraplate_area_m2 = 0.0D0 DO i = 1, nRows DO j = 1, nColumns IF ((GSRM2_grid(i,j,1)/=0.0).OR.(GSRM2_grid(i,j,2)/=0.0).OR.(GSRM2_grid(i,j,3)/=0.0)) THEN plate_boundary_area_m2 = plate_boundary_area_m2 + cell_area_m2(i) ELSE intraplate_area_m2 = intraplate_area_m2 + cell_area_m2(i) END IF END DO END DO percent_boundary_area = 100. * plate_boundary_area_m2 / sphere_area_m2 percent_intraplate_area = 100. * intraplate_area_m2 / sphere_area_m2 WRITE (*, "(' ')") WRITE (*, "(' According to the GSRM2 data:')") WRITE (*, "(' Plate boundary area = ',ES12.5,' m^2 (',F6.3,'%).')") plate_boundary_area_m2, percent_boundary_area WRITE (*, "(' Intraplate area = ',ES12.5,' m^2 (',F6.3,'%).')") intraplate_area_m2, percent_intraplate_area WRITE (*, "(' and we can assume the 2 polar caps are intraplate, as well.')") WRITE (21, *) WRITE (21, "('According to the GSRM2 data:')") WRITE (21, "('Plate boundary area = ',ES12.5,' m^2 (',F6.3,'%).')") plate_boundary_area_m2, percent_boundary_area WRITE (21, "('Intraplate area = ',ES12.5,' m^2 (',F6.3,'%).')") intraplate_area_m2, percent_intraplate_area WRITE (21, "('and we can assume the 2 polar caps are intraplate, as well.')") !Get seismic catalog file for calibration WRITE (*, *) 10 WRITE (*, "(' Enter name of (global, shallow) seismic catalog .EQC file for calibration:')") READ (*, "(A)") EQC_file OPEN (UNIT = 2, FILE = TRIM(EQC_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ',A,' not found. Fix, and try again...')") TRIM(EQC_file) CALL Pause() GO TO 10 END IF WRITE (21, "('Enter name of (global, shallow) seismic catalog .EQC file for calibration:')") WRITE (21, "(A)") EQC_file READ (2, "(A)") first_line quaking: DO READ (2, "(A)", IOSTAT = ios) line IF (ios == 0) THEN last_line = line ELSE EXIT quaking END IF END DO quaking CLOSE (2) !Get duration of seismic catalog: WRITE (*, *) WRITE (*, "(' Here are the first and last earthquakes in the catalog: ')") WRITE (*, "(' ',A)") TRIM(first_line) WRITE (*, "(' ',A)") TRIM(last_line) WRITE (*, "(' What is the duration of this catalog, in years? '\)") READ (*, *) catalog_years WRITE (21, *) WRITE (21, "('Here are the first and last earthquakes in the catalog: ')") WRITE (21, "(A)") TRIM(first_line) WRITE (21, "(A)") TRIM(last_line) WRITE (21, "('What is the duration of this catalog, in years? '\)") WRITE (21, *) catalog_years !Get threshold magnitude of seismic catalog: WRITE (*, *) WRITE (*, "(' What is the completeness-threshold magnitude of this catalog? '\)") READ (*, *) catalog_threshold_magnitude WRITE (21, *) WRITE (21, "('What is the completeness-threshold magnitude of this catalog? '\)") WRITE (21, *) catalog_threshold_magnitude catalog_threshold_moment_Nm = Moment(catalog_threshold_magnitude) WRITE (*, "(' Catalog threshold (minimum) scalar moment =',ES11.3,' N m')") catalog_threshold_moment_Nm WRITE (21, "('Catalog threshold (minimum) scalar moment =',ES11.3,' N m')") catalog_threshold_moment_Nm WRITE (*, "(' This threshold will be used for calibration of the forecast,')") WRITE (*, "(' regardless of the requested threshold for forecast reporting.')") WRITE (21, "('This threshold will be used for calibration of the forecast,')") WRITE (21, "('regardless of the requested threshold for forecast reporting.')") !-------------------------------------------------------------------------------- !Select model iVersion: WRITE (*, *) WRITE (*, "(' PLEASE SELECT DESIRED MODEL VERSION FROM THE FOLLOWING LIST:')") WRITE (*, "(' 1(a): Uniform (calibrated) coupled thickness in all boundaries [2 DOF].')") WRITE (*, "(' 2(b): Add isotropic smoothing of offshore non-SUB boundaries [2 DOF].')") WRITE (*, "(' 3(c): Add seaward smoothing of SUB boundaries to outer rises [2 DOF].')") WRITE (*, "(' 4(d): Add separate calibrations in tectonic zones #0~4 [5 DOF].')") WRITE (*, "(' 5(e): Add classification & partitioning of strain rate tensors [5 DOF].')") WRITE (*, "(' 6(f): Add velocity-dependence of SUB and CCB coupling [5 DOF].')") 20 WRITE (*, "(' CHOOSE INTEGER CODE: '\)") READ (*, *) iVersion IF ((iVersion < 1).OR.(iVersion > last_version)) THEN WRITE (*, "(' ERROR: Illegal integer entered. Try again...')") CALL Pause() GO TO 20 END IF WRITE (21, *) WRITE (21, "('PLEASE SELECT DESIRED MODEL VERSION FROM THE FOLLOWING LIST:')") WRITE (21, "(' 1(a): Uniform (calibrated) coupled thickness in all boundaries [2 DOF].')") WRITE (21, "(' 2(b): Add isotropic smoothing of offshore non-SUB boundaries [2 DOF].')") WRITE (21, "(' 3(c): Add seaward smoothing of SUB boundaries to outer rises [2 DOF].')") WRITE (21, "(' 4(d): Add separate calibrations in tectonic zones #0~4 [5 DOF].')") WRITE (21, "(' 5(e): Add classification & partitioning of strain rate tensors [5 DOF].')") WRITE (21, "(' 6(f): Add velocity-dependence of SUB and CCB coupling [5 DOF].')") WRITE (21, "('CHOOSE INTEGER CODE: '\)") WRITE (21, *) iVersion !Select forecast_threshold_magnitude: WRITE (*, *) WRITE (*, "(' Enter desired magnitude threshold for output forecast: '\)") READ (*, *) forecast_threshold_magnitude WRITE (21, *) WRITE (21, "('Enter desired magnitude threshold for output forecast: '\)") WRITE (21, *) forecast_threshold_magnitude forecast_threshold_moment_Nm = Moment(forecast_threshold_magnitude) WRITE (*, "(' Forecast threshold (minimum) scalar moment =',ES11.3,' N m')") forecast_threshold_moment_Nm WRITE (21, "('Forecast threshold (minimum) scalar moment =',ES11.3,' N m')") forecast_threshold_moment_Nm !Set GRD_file name for output of seismicity (EQ/m^2/s) results: WRITE (*, *) WRITE (*, "(' Create a .GRD file-name for the output: ')") READ (*, "(A)") GRD_file WRITE (21, *) WRITE (21, "('Create a .GRD file-name for the output: ')") WRITE (21, "(A)") GRD_file !-------------------------------------------------------------------------------- IF (iVersion == 1) THEN !Count shallow earthquakes in calibration period (both boundary & intraplate), !at or above the catalog (not forecast) threshold magnitude, !and prepare to save catalog earthquakes designated "intraplate" to a new .EQC file: boundary_calibration_count = 0 intraplate_calibration_count = 0 OPEN (UNIT = 2, FILE = TRIM(EQC_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios) OPEN (UNIT = 20, FILE = "SHIFT_GSRM2a_intraplate.eqc") ticking: DO READ (2, 100, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data ! Note: PAD = "YES" permits READ with no appended data 100 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, & & A) IF (ios /= 0) EXIT ticking IF (eq_mag >= catalog_threshold_magnitude) THEN i = 1 + NINT((GSRM2_lat_max - eq_Nlat) / GSRM2_d_lat) IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 j = 1 + NINT((eq_Elon - GSRM2_lon_min) / GSRM2_d_lon) IF ((GSRM2_grid(i,j,1)/=0.0).OR.(GSRM2_grid(i,j,2)/=0.0).OR.(GSRM2_grid(i,j,3)/=0.0)) THEN boundary_calibration_count = boundary_calibration_count + 1 ELSE intraplate_calibration_count = intraplate_calibration_count + 1 WRITE (20, 100) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data END IF END IF END DO ticking CLOSE (2) ! EQC_file, for calibration CLOSE (20) ! new EQC file with events (above catalog threshold) selected as "intraplate" percent_boundary_eqs = 100. * boundary_calibration_count / & & (boundary_calibration_count + intraplate_calibration_count) percent_intraplate_eqs = 100. * intraplate_calibration_count / & & (boundary_calibration_count + intraplate_calibration_count) WRITE (*, *) WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate boundaries: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, boundary_calibration_count, percent_boundary_eqs WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs intraplate_seismicity_SI = intraplate_calibration_count / & & (intraplate_area_m2 * catalog_years * s_per_year) WRITE (*, "(' Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI WRITE (21, *) WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate boundaries: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, boundary_calibration_count, percent_boundary_eqs WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs WRITE (21, "('Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI !First-pass analysis of boundaries: Treat all as subduction zones !(N.B. This will be checked and corrected by a common factor, later on.) tentative_boundary_eqs_per_s = 0.0D0 DO i = 1, nRows DO j = 1, nColumns exx = GSRM2_grid(i, j, 1) eyy = GSRM2_grid(i, j, 2) exy = GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value T5_column = 10 ! SUB CALL SHIFT_continuum_seismicity (catalog_threshold_moment_Nm, & ! INTENT(IN) & e1h, e2h, err, T5_column, & ! INTENT(IN) & seismicity) ! INTENT(OUT) tentative_boundary_eqs_per_s = tentative_boundary_eqs_per_s + seismicity * cell_area_m2(i) !Note that this summation uses the catalog (not forecast) threshold magnitude. !Also, it does not include the underlying "floor" of intraplate_seismicity_SI which is everywhere. END DO END DO tentative_boundary_eqs_per_catalog = tentative_boundary_eqs_per_s * catalog_years * s_per_year boundary_scaling_factor = ( boundary_calibration_count - & & (plate_boundary_area_m2 * intraplate_seismicity_SI * catalog_years * s_per_year) ) / & & tentative_boundary_eqs_per_catalog WRITE (*, *) WRITE (*, "(' Global plate-boundary seismicity scaling = ',F7.5,' (relative to SUB).')") boundary_scaling_factor WRITE (21, *) WRITE (21, "('Global plate-boundary seismicity scaling = ',F7.5,' (relative to SUB).')") boundary_scaling_factor !Scale the plate-boundary seismicity, and impose intraplate value as floor (both at forecast threshold): !First, determine how much intraplate seismicity shifts at forecast threshold (relative to catalog): !Use intraplate corner magnitude of 9.0 and beta = 0.63 inferred by Bird et al. [2010, SRL]. intraplate_forecast_over_catalog_factor = (( forecast_threshold_moment_Nm / catalog_threshold_moment_Nm )**(-0.63)) * & & EXP(( catalog_threshold_moment_Nm - forecast_threshold_moment_Nm ) / & & Moment(9.0) ) !using equation (7) in Bird et al. [2010, SRL]. intraplate_seismicity_SI = intraplate_seismicity_SI * intraplate_forecast_over_catalog_factor !Then, recompute all plate-boundary seismicities at forecast threshold, and scale: DO i = 1, nRows DO j = 1, nColumns exx = GSRM2_grid(i, j, 1) eyy = GSRM2_grid(i, j, 2) exy = GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value T5_column = 10 ! SUB CALL SHIFT_continuum_seismicity (forecast_threshold_moment_Nm, & ! INTENT(IN) & e1h, e2h, err, T5_column, & ! INTENT(IN) & seismicity) ! INTENT(OUT) epicenter_rate_density(i, j) = seismicity * boundary_scaling_factor + intraplate_seismicity_SI IF ((exx /= 0.0).OR.(eyy /= 0.0).OR.(exy /= 0.0)) THEN effective_cz_meters(i, j) = NINT(boundary_scaling_factor * T5_coupledThickness_m(T5_column)) mean_G(i, j) = T5_assumed_mu_Pa(T5_column) END IF END DO END DO !----------------------------------------------------------------------------- ELSE ! iVersion > 1; strain-rates must be smoothed in the offshore regions !locate, read, and memorize underwater.grd WRITE (21, *) WRITE (21, "('File underwater.grd must be available. Please check this.')") 200 WRITE (*, *) WRITE (*, "(' File underwater.grd must be available. Please check this.')") CALL Pause() OPEN (UNIT = 1, FILE = "underwater.grd", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Copy file underwater.grd was not found in active folder.')") CALL Pause() GO TO 200 END IF READ (1, *) underwater_lon_min, underwater_d_lon, underwater_lon_max READ (1, *) underwater_lat_min, underwater_d_lat, underwater_lat_max underwater_rows = 1 + NINT((underwater_lat_max - underwater_lat_min) / underwater_d_lat) underwater_columns = 1 + NINT((underwater_lon_max - underwater_lon_min) / underwater_d_lon) ALLOCATE ( underwater_INT_GRD(underwater_rows, underwater_columns) ) ALLOCATE ( underwater_TF_GRD(underwater_rows, underwater_columns) ) READ (1, *)((underwater_INT_GRD(i, j), j = 1, underwater_columns), i = 1, underwater_rows) DO i = 1, underwater_rows DO j = 1, underwater_columns underwater_TF_GRD(i, j) = (underwater_INT_GRD(i, j) == 1) END DO END DO CLOSE(1) !locate and open PB2002_steps.dat WRITE (21, *) WRITE (21, "('File PB2002_steps.dat must be available. Please check this.')") 210 WRITE (*, *) WRITE (*, "(' File PB2002_steps.dat must be available. Please check this.')") IF (iVersion >= 4) THEN ! also remind user to provide tectonic zones .GRD: WRITE (*, *) WRITE (*, "(' File PB2002_tectonic_zones.grd must also be available. Please check this.')") END IF CALL Pause() OPEN (UNIT = 1, FILE = "PB2002_steps.dat", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Copy file PB2002_steps.dat was not found in active folder.')") CALL Pause() GO TO 210 END IF !memorize contents as a list of categorized, step-interior uvecs: ALLOCATE ( PB2002_uvecs(3, PB2002_list_max) ) PB2002_uvecs = 0.0 ALLOCATE ( PB2002_code(PB2002_list_max) ) PB2002_code = 0 ALLOCATE ( PB2002_elevation(PB2002_list_max) ) PB2002_elevation = 0.0 ALLOCATE ( PB2002_seaward_azimuth(PB2002_list_max) ) ! only relevant to SUB boundaries PB2002_seaward_azimuth = 0 ! INTEGER*2 ALLOCATE ( PB2002_velocity_mmpa(PB2002_list_max) ) PB2002_velocity_mmpa = 0.0 line_count = 0 PB2002_list_i = 0 bounding: DO READ (1, "(8X,A1,2X,F9.3,F8.3,F9.3,F8.3,F6.1,I4,F6.1,18X,I7,6X,A3)", IOSTAT = ios) & dip_c1, lon1, lat1, lon2, lat2, km_long, nAzim, velocity_mmpa, nElevation_m, c3 IF (dip_c1 == '\') THEN nAzim = nAzim -90 ELSE IF (dip_c1 == '/') THEN nAzim = nAzim + 90 ELSE nAzim = 0 ! (and it will not be used) END IF IF (ios /= 0) EXIT bounding line_count = line_count + 1 !N.B. the following are in order of mean , to minimize effects of possible inappropriate ! rounding in the ratio smoothed_GSRM2_code(i, j) = NINT(smoothed_GSRM2_utility(i, j, 1) / & ! & smoothed_GSRM2_utility(i, j, 2)) ! All values quoted in comments are from Table 5 of Bird & Kagan [2004, BSSA]. IF (c3 == "OSR") THEN ! = 0.13 km (for normal EQs) code_now = 1 ELSE IF (c3 == "CRB") THEN ! = 3 km code_now = 2 ELSE IF (c3 == "OCB") THEN ! = 3.8 km code_now = 3 ELSE IF (c3 == "OTF") THEN ! = 1.6~13 km, depending on speed code_now = 4 ELSE IF (c3 == "CTF") THEN ! = 8.6 km code_now = 5 ELSE IF (c3 == "CCB") THEN ! = 18 km code_now = 6 ELSE IF (c3 == "SUB") THEN ! = 18 km, but less or more depending on speed (fast = 2 x slow). code_now = 7 ELSE WRITE (*, "(' ERROR: Illegal c3 code in PB2002_steps.dat, line ',I6)") line_count CALL Pause() STOP END IF internal_points = MAX(1, (NINT(km_long * 1000.0 / PB2002_ministep_m) - 1)) CALL LonLat_2_Uvec(lon1, lat1, uvec1) CALL LonLat_2_Uvec(lon2, lat2, uvec2) arc_radians = Arc(uvec1, uvec2) azimuth_radians = Relative_Compass (from_uvec = uvec1, to_uvec = uvec2) DO j = 1, internal_points fraction = j / (internal_points + 1.0) CALL Turn_To (azimuth_radians = azimuth_radians, & & base_uvec = uvec1, & & far_radians = fraction * arc_radians, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = mid_uvec) PB2002_list_i = PB2002_list_i + 1 IF (PB2002_list_i > PB2002_list_max) THEN WRITE (*, "(' ERROR: Increase PARAMETER PB2002_list_max and recompile.')") CALL Pause() STOP END IF PB2002_uvecs(1:3, PB2002_list_i) = mid_uvec(1:3) ! memorize internal point PB2002_code(PB2002_list_i) = code_now ! 1 = CCB, 2 = CRB, etc. PB2002_elevation(PB2002_list_i) = nElevation_m ! converting INT to REAL PB2002_seaward_azimuth(PB2002_list_i) = nAzim ! INTEGER*2 PB2002_velocity_mmpa(PB2002_list_i) = velocity_mmpa END DO END DO bounding PB2002_list_highest = PB2002_list_i WRITE (*, "(' PB2002_steps.dat has been converted to a list of ',I6,' internal points.')") PB2002_list_highest WRITE (21, "('PB2002_steps.dat has been converted to a list of ',I6,' internal points.')") PB2002_list_highest CLOSE (1) ! PB2002_steps.dat !Find PB2002 plate boundary nearest to each active GSRM2 grid cell... WRITE (*, "(' Finding PB2002 plate boundary nearest to each active GSRM2 grid cell...')") WRITE (21, "('Finding PB2002 plate boundary nearest to each active GSRM2 grid cell...')") ALLOCATE ( GSRM2_closest_PB2002_radian2(nRows, nColumns) ) GSRM2_closest_PB2002_radian2 = 1.0 ! "not close at all!" ALLOCATE ( GSRM2_closest_PB2002_code(nRows, nColumns) ) GSRM2_closest_PB2002_code = 0 ! and this will remain zero for inactive GSRM2 cells ALLOCATE ( GSRM2_closest_PB2002_elevation(nRows, nColumns) ) GSRM2_closest_PB2002_elevation = 0.0 ALLOCATE ( GSRM2_closest_PB2002_velocity_mmpa(nRows, nColumns) ) GSRM2_closest_PB2002_velocity_mmpa = 0.0 ALLOCATE ( GSRM2_closest_PB2002_seaward_azimuth(nRows, nColumns) ) GSRM2_closest_PB2002_seaward_azimuth = 0 ! INTEGER*2; set to 0 for all non-SUB boundaries. ALLOCATE ( GSRM2_raw_SUB(nRows, nColumns) ) GSRM2_raw_SUB = .FALSE. old_percent_done = 0 test_radians2 = (200000. / Earth_radius_m)*2 ! 200 km limit on distance from trench (for raw_SUB designation) DO i = 1, nRows lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat DO j = 1, nColumns IF ((GSRM2_grid(i,j,1)/=0.0).OR.(GSRM2_grid(i,j,2)/=0.0).OR.(GSRM2_grid(i,j,3)/=0.0)) THEN !conduct a search: lon = GSRM2_lon_min + (j-1) * GSRM2_d_lon CALL LonLat_2_Uvec(lon, lat, uvec1) standard = 1.0 ! not close at all! DO k = 1, PB2002_list_highest uvec2(1:3) = PB2002_uvecs(1:3, k) arc_radian2 = (uvec1(1) - uvec2(1))**2 + & & (uvec1(2) - uvec2(2))**2 + & & (uvec1(3) - uvec2(3))**2 IF (arc_radian2 < standard) THEN standard = arc_radian2 best_code = PB2002_code(k) best_elevation = PB2002_elevation(k) best_azimuth = PB2002_seaward_azimuth(k) best_velocity_mmpa = PB2002_velocity_mmpa(k) END IF END DO GSRM2_closest_PB2002_radian2(i, j) = standard GSRM2_closest_PB2002_code(i, j) = best_code GSRM2_closest_PB2002_elevation(i, j) = best_elevation GSRM2_closest_PB2002_velocity_mmpa(i, j) = best_velocity_mmpa GSRM2_closest_PB2002_seaward_azimuth(i, j) = best_azimuth ! INTEGER*2, and 0 for non-SUB GSRM2_raw_SUB(i, j) = (best_code == 7).AND.(standard < test_radians2) ! IF (TRUE), then ! a candidate for later sandblasting, IF ((iVersion >= 3).AND. ! ({not cancelled below, during isotropic smoothing}) END IF END DO percent_done = (100. * i) / nRows IF (percent_done > old_percent_done) THEN WRITE (*, "(' ',I3,'% done...')") percent_done old_percent_done = percent_done END IF END DO WRITE (*, "(' Search completed.')") WRITE (21, "('Search completed.')") !Create isotropically-smoothed version of GSRM2 strain rates about offshore non-SUB boundaries: WRITE (*, *) WRITE (*, "(' Isotropically smoothing rates adjacent to offshore non-SUB plate boundaries...')") WRITE (21, *) WRITE (21, "('Isotropically smoothing rates adjacent to offshore non-SUB plate boundaries...')") ALLOCATE ( smoothed_GSRM2_grid(nRows, nColumns, 3) ) smoothed_GSRM2_grid = 0.0 test_distance_m = 186000. ! = sup(OSR, OTF, OCB) half-widths, from last row of Table 1 of Bird & Kagan [2004, BSSA]. DO i = 1, nRows lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat DO j = 1, nColumns lon = GSRM2_lon_min + (j-1) * GSRM2_d_lon IF ((GSRM2_grid(i,j,1)/=0.0).OR.(GSRM2_grid(i,j,2)/=0.0).OR.(GSRM2_grid(i,j,3)/=0.0)) THEN !test1 = nearest plate boundary is underwater: test1 = (GSRM2_closest_PB2002_elevation(i, j) < 0.0) !Kludge: Apply isotropic smoothing to new plate boundary offshore Baja California, ! which was not present in PB2002, but added for GSRM2: IF ((lat > 22.5).AND.(lat < 32.7).AND.(lon > -120.).AND.(lon < -110.).AND. & & (lat < (34. - lon - (-120.)))) test1 = .TRUE. !Kludge: Apply isotropic smoothing to revised plate boundary offshore South Africa, ! which is in a new location in GSRM2 compared to PB2002. IF ((lat > -35.).AND.(lat < -25.9).AND.(lon > 35.).AND.(lon < 38.)) test1 = .TRUE. !test2 = nearest plate boundary is NOT type SUB test2 = (GSRM2_closest_PB2002_code(i, j) /= 7) !Kludge: Apply isotropic smoothing to PA\SC & PA\SL subduction zones, because ! both are represented by a single line of cells in GSRM2: IF ((lat > -63.).AND.(lat < -52.1).AND.(lon > -79.).AND.(lon < -56.)) test2 = .TRUE. !Kludge: Apply isotropic smoothing to SS\SB & SS\WL subduction zones, because ! both are represented by a single line of cells in GSRM2: IF ((lat > -9.5).AND.(lat < -5.5).AND.(lon > 148.5).AND.(lon < 154.5)) test2 = .TRUE. !test3 is whether distance is acceptable: test3 = ((SQRT(GSRM2_closest_PB2002_radian2(i, j)) * Earth_radius_m) <= test_distance_m) !Kludge: Apply isotropic smoothing to new plate boundary offshore Baja California, ! which was not present in PB2002, but added for GSRM2: IF ((lat > 22.5).AND.(lat < 32.7).AND.(lon > -120.).AND.(lon < -110.).AND. & & (lat < (34. - lon - (-120.)))) test3 = .TRUE. !Kludge: Apply isotropic smoothing to revised plate boundary offshore South Africa, ! which is in a new location in GSRM2 compared to PB2002. IF ((lat > -35.).AND.(lat < -25.9).AND.(lon > 35.).AND.(lon < 38.)) test3 = .TRUE. !test4 is whether cell is adjacent to rigid plate interior on any of its 4 sides: IF ((i == 1).OR.(i == nRows)) THEN test4 = .TRUE. ELSE ! a normal row, with neighbors to N and S iN = i - 1 iS = i + 1 jW = j - 1 IF (j == 1) jW = nColumns jE = j + 1 IF (j == nColumns) jE = 1 test4 = (GSRM2_closest_PB2002_code(iN, j) == 0) .OR. & & (GSRM2_closest_PB2002_code(i, jW) == 0) .OR. & & (GSRM2_closest_PB2002_code(i, jE) == 0) .OR. & & (GSRM2_closest_PB2002_code(iS, j) == 0) END IF !test5 is whether donor cell is underwater: test5 = underwater_TF_GRD(i, j) !Special exemption from the no-smoothing-on-land rule, for the New Britain area, because: ! (a) with only 2 GPS sites on island, rigidity cannot be demonstrated; and ! (b) if treated as "intraplate" it collects ~40 "intraplate" EQs in 1977-2004 and seriously biases the global rate! IF ((lon > 143.).AND.(lon < 156.).AND.(lat > -10.).AND.(lat < -1.)) THEN test5 = .TRUE. END IF !all tests must be passed, or there is no smoothing: smooth = test1.AND.test2.AND.test3.AND.test4.AND.test5 !but note that some quanta of the smoothed strain-rate could still be returned if receiver cells are not underwater. IF (smooth) THEN ! redistribute strain-rate tensor (and weighted code) into smoothed_GSRM2_grid GSRM2_raw_SUB(i, j) = .FALSE. ! Whether SUB or not, it is not longer "raw_SUB" after isotropic smoothing; ! so sandblasting should be blocked (IF (iVersion >= 3)). !locate source cell center: lat1 = GSRM2_lat_max - (i-1) * GSRM2_d_lat lon1 = GSRM2_lon_min + (j-1) * GSRM2_d_lon CALL LonLat_2_Uvec(lon1, lat1, uvec1) !determine extent of smoothing: IF (GSRM2_closest_PB2002_code(i, j) == 1) THEN ! OSR sigma_m = 66000. ! last row, Table 1, Bird & Kagan [2004, BSSA] ELSE IF (GSRM2_closest_PB2002_code(i, j) == 2) THEN ! CRB sigma_m = 77000. ELSE IF (GSRM2_closest_PB2002_code(i, j) == 3) THEN ! OCB sigma_m = 93000. ELSE IF (GSRM2_closest_PB2002_code(i, j) == 4) THEN ! OTF sigma_m = 64000. ELSE IF (GSRM2_closest_PB2002_code(i, j) == 5) THEN ! CTF sigma_m = 128500. ELSE IF (GSRM2_closest_PB2002_code(i, j) == 6) THEN ! CCB sigma_m = 94500. ELSE IF (GSRM2_closest_PB2002_code(i, j) == 7) THEN ! SUB !N.B. This occurs rarely; only where I add special local "kludge" code ! to permit local isotropic smoothing of a subduction zone which was ! represented by a single row of deforming cells in GSRM2. ! The sigma value is the mean of the forearc & backarc sides of the ! distributions shown in Table 1 of Bird & Kagan [2004, BSSA]. sigma_m = 88750. ELSE WRITE (*, "(' ERROR: At row ',I6,' & column ',I6,' GSRM2_closest_PB2002_code = ',I2,'(illegal).')") & & i, j, GSRM2_closest_PB2002_code(i, j) WRITE (21, "('ERROR: At row ',I6,' & column ',I6,' GSRM2_closest_PB2002_code = ',I2,'(illegal).')") & & i, j, GSRM2_closest_PB2002_code(i, j) CALL Pause() STOP END IF three_sigma_m = 3.0 * sigma_m dRow_3_sigma = NINT(three_sigma_m / (GSRM2_d_lat * radians_per_degree * Earth_radius_m)) dCol_3_sigma = NINT(three_sigma_m / (GSRM2_d_lon * COS(lat1 * radians_per_degree) * & & radians_per_degree * Earth_radius_m)) DO ii = -dRow_3_sigma, dRow_3_sigma ! both limits are integers iTarget = i + ii IF ((iTarget >= 1).AND.(iTarget <= nRows)) THEN DO jj = -dCol_3_sigma, dCol_3_sigma ! both limits are integers jTarget = j + jj IF (jTarget < 1) jTarget = jTarget + nColumns IF (jTarget > nColumns) jTarget = jTarget - nColumns IF ((iTarget == i).AND.(jTarget == j)) THEN distance_m = 0.0 ELSE lat2 = GSRM2_lat_max - (iTarget-1) * GSRM2_d_lat lon2 = GSRM2_lon_min + (jTarget-1) * GSRM2_d_lon CALL LonLat_2_Uvec(lon2, lat2, uvec2) distance_m = Arc(uvec1, uvec2) * Earth_radius_m END IF IF (distance_m <= three_sigma_m) THEN coefficient = (cell_area_m2(i) / (2.0 * Pi * sigma_m**2)) * & & EXP(-0.5 * distance_m**2 / sigma_m**2) eligable_receiver = underwater_TF_GRD(iTarget, jTarget) !Special exemption from the no-smoothing-on-land rule, for the New Britain area, because: ! (a) with only 2 GPS sites on island, rigidity cannot be demonstrated; and ! (b) if treated as "intraplate" it collects ~40 "intraplate" EQs in 1977-2004 and seriously biases the global rate! IF ((lon > 143.).AND.(lon < 156.).AND.(lat > -10.).AND.(lat < -1.)) THEN eligable_receiver = .TRUE. END IF IF (eligable_receiver) THEN ! add a bit of strain-rate to target: smoothed_GSRM2_grid(iTarget, jTarget, 1:3) = smoothed_GSRM2_grid(iTarget, jTarget, 1:3) + & & coefficient * GSRM2_grid(i, j, 1:3) ELSE ! target is illegal; do not transfer; return this bit to (the smoothed copy of) its donor cell. smoothed_GSRM2_grid(i , j , 1:3) = smoothed_GSRM2_grid(i , j , 1:3) + & & coefficient * GSRM2_grid(i, j, 1:3) END IF ! target cell is underwater, or not END IF END DO END IF END DO ELSE ! just add it, without modification: smoothed_GSRM2_grid(i, j, 1:3) = smoothed_GSRM2_grid(i, j, 1:3) + & & GSRM2_grid(i, j, 1:3) END IF END IF END DO END DO WRITE (*, "(' Isotropic smoothing of offshore non-SUB plate boundaries completed.')") WRITE (21, "('Isotropic smoothing of offshore non-SUB plate boundaries completed.')") IF (iVersion >= 3) THEN ! add directed smoothing of SUB boundaries toward outer rises: WRITE (*, *) WRITE (*, "(' Smoothing SUB rates toward the outer rise regions...')") WRITE (21, *) WRITE (21, "('Smoothing SUB rates toward the outer rise regions...')") !----Initialize IMSL version of random numbers (below): !CALL RNSET(234579) ! seeding the random number generator !--------------------------------------------------------- !----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 old_percent_done = 0 DO iteration = 1, 440 ! See MATLAB code ToSea.m for "toy" version in 1-D, used to calibrate. DO i = 2, (nRows-1) ! can't erode N-most or S-most row because I can't compare to their neighbors. lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat DO j = 1, nColumns IF (GSRM2_raw_SUB(i, j).AND.underwater_TF_GRD(i, j)) THEN ! Source: a previously-unsmoothed SUB region which is underwater: !decide transport distance (in m) first: !using IMSL method: !travel_m = 66000. * (0.7 * ANORIN(RNUNF()) + 1.0) ! 66 km * [Gaussian with mean 1 and s.d. of 0.6] !using MKL-VSL method: 66 km * [Gaussian with mean 1 and s.d. of 0.6] status = vsRngGaussian( method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, stream = my_stream, n = 1, r = MKL_VSL_random, a = 1.0, sigma = 0.6 ) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: status = vsRngGaussian() returned status = ',I12)") status CALL Pause(); STOP ELSE travel_m = 66000. * MKL_VSL_random(1) END IF travel_m = MAX(travel_m, 0.0) nAzim = GSRM2_closest_PB2002_seaward_azimuth(i, j) azimuth_radians = nAzim * radians_per_degree NS_fraction = ABS(COS(azimuth_radians)) ! raw, not yet normalized ... EW_fraction = ABS(SIN(azimuth_radians)) NS_fraction = NS_fraction / (NS_fraction + EW_fraction) ! normalize EW_fraction = EW_fraction / (NS_fraction + EW_fraction) !Decide which type of transport is active in this iteration: !IMSL method: !NS_this_time = (RNUNF() < NS_fraction) ! T or F; if F, then EW_this_time is implied !MKL-VSL method: 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 WRITE (*, "(' ERROR: status = vsRngUniform() returned status = ',I12)") status CALL Pause(); STOP ELSE NS_this_time = (MKL_VSL_random(1) < NS_fraction) END IF IF (NS_this_time) THEN ! decide based on rows to N & S: this_Sup = MAX(ABS(smoothed_GSRM2_grid(i,j,1)),ABS(smoothed_GSRM2_grid(i,j,2)),ABS(smoothed_GSRM2_grid(i,j,3))) N_Sup = MAX(ABS(smoothed_GSRM2_grid(i-1,j,1)),ABS(smoothed_GSRM2_grid(i-1,j,2)),ABS(smoothed_GSRM2_grid(i-1,j,3))) S_Sup = MAX(ABS(smoothed_GSRM2_grid(i+1,j,1)),ABS(smoothed_GSRM2_grid(i+1,j,2)),ABS(smoothed_GSRM2_grid(i+1,j,3))) IF (this_Sup >= 1.001 * (N_Sup + S_Sup)/2) THEN ! move 1% of strain-rates: bins = NINT(travel_m / (GSRM2_d_lat * radians_per_degree * Earth_radius_m)) IF (bins /= 0) THEN !decide direction: IF (COS(azimuth_radians) > 0.0) THEN ! move to N (even if off-the-grid): iTarget = i - bins ELSE ! move to S (even if off-the-grid): iTarget = i + bins END IF IF ((iTarget >= 1).AND.(iTarget <= nRows)) THEN !Target column may be deflected E or W by 1 column: !IMSL method: !tester = RNUNF() ! random number in [0., 1.] !MKL-VSL method: 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 WRITE (*, "(' ERROR: status = vsRngUniform() returned status = ',I12)") status CALL Pause(); STOP ELSE tester = MKL_VSL_random(1) END IF IF (tester < 0.1) THEN ! deflect to W: jTarget = j - 1 IF (jTarget < 1) jTarget = jTarget + nColumns ELSE IF (tester > 0.9) THEN ! deflect to E: jTarget = j + 1 IF (jTarget > nColumns) jTarget = jTarget - nColumns ELSE ! no deflection about 80% of the time: jTarget = j END IF !transfer only occurs if deposition point is underwater, too: IF (underwater_TF_GRD(iTarget, jTarget)) THEN smoothed_GSRM2_grid(iTarget, jTarget, 1:3) = smoothed_GSRM2_grid(iTarget, jTarget, 1:3) + & & 0.001 * smoothed_GSRM2_grid(i, j, 1:3) !remove material from source bin: smoothed_GSRM2_grid(i, j, 1:3) = 0.999 * smoothed_GSRM2_grid(i, j, 1:3) END IF END IF END IF ! bins /= 0 END IF ELSE ! EW_this_time; decide based on columns to W & E: this_Sup = MAX(ABS(smoothed_GSRM2_grid(i,j,1)),ABS(smoothed_GSRM2_grid(i,j,2)),ABS(smoothed_GSRM2_grid(i,j,3))) jBeside = j - 1; IF (jBeside == 0) jBeside = nColumns W_Sup = MAX(ABS(smoothed_GSRM2_grid(i,jBeside,1)),ABS(smoothed_GSRM2_grid(i,jBeside,2)),ABS(smoothed_GSRM2_grid(i,jBeside,3))) jBeside = j + 1; IF (jBeside > nColumns) jBeside = 1 E_Sup = MAX(ABS(smoothed_GSRM2_grid(i,jBeside,1)),ABS(smoothed_GSRM2_grid(i,jBeside,2)),ABS(smoothed_GSRM2_grid(i,jBeside,3))) IF (this_Sup >= 1.001 * (W_Sup + E_Sup)/2) THEN ! move 1% of strain-rates: bins = NINT(travel_m / (GSRM2_d_lon * COS(lat * radians_per_degree) * radians_per_degree * Earth_radius_m)) IF (bins /= 0) THEN !decide direction: IF (SIN(azimuth_radians) > 0.0) THEN ! move to E: jTarget = j + bins IF (jTarget > nColumns) jTarget = jTarget - nColumns ELSE ! move to W: jTarget = j - bins IF (jTarget < 1) jTarget = jTarget + nColumns END IF !Target row may be deflected N or S by 1 row: !tester = RNUNF() ! random number in [0., 1.] by IMSL method !now using MKL-VSL method instead: 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 WRITE (*, "(' ERROR: status = vsRngUniform(...) returns status = ',I12)") status CALL Pause(); STOP ELSE tester = MKL_VSL_random(1) END IF IF (tester < 0.1) THEN ! deflect to S: iTarget = i + 1 iTarget = MIN(iTarget, nRows) ELSE IF (tester > 0.9) THEN ! deflect to N: iTarget = i - 1 iTarget = MAX(iTarget, 1) ELSE ! no deflection about 80% of the time: iTarget = i END IF !transfer only occurs if deposition point is underwater, too: IF (underwater_TF_GRD(iTarget, jTarget)) THEN smoothed_GSRM2_grid(iTarget, jTarget, 1:3) = smoothed_GSRM2_grid(iTarget, jTarget, 1:3) + & & 0.001 * smoothed_GSRM2_grid(i, j, 1:3) !remove material from source bin: smoothed_GSRM2_grid(i, j, 1:3) = 0.999 * smoothed_GSRM2_grid(i, j, 1:3) END IF END IF ! bins /= 0 END IF END IF ! NS_this_time, or EW_this_time END IF ! a SUB point (in smoothed grid) END DO ! j = 1, nColumns END DO ! i = 1, nRows percent_done = (100. * iteration) / 440 IF (percent_done > old_percent_done) THEN WRITE (*, "(' ',I3,'% done...')") percent_done old_percent_done = percent_done END IF END DO ! 440 iterations WRITE (*, "(' Smoothing of SUB rates toward outer rises completed.')") WRITE (21, "('Smoothing of SUB rates toward outer rises completed.')") END IF ! iVersion >= 3 (directed smoothing of SUBs) IF (iVersion < 4) THEN ! N.B. In present context, we known iVersion >= 2, so the possibilities are 2 or 3. ! Code in this block is for models with spatial smoothing, but only 2 DOF ! (plate-boundary vs. intraplate constants). !In cases of iVersion == 2 or 3, !Re-compute total areas of straining plate boundaries, and of intraplate regions: plate_boundary_area_m2 = 0.0D0 intraplate_area_m2 = 0.0D0 DO i = 1, nRows DO j = 1, nColumns IF ((smoothed_GSRM2_grid(i,j,1)/=0.0).OR.(smoothed_GSRM2_grid(i,j,2)/=0.0).OR.(smoothed_GSRM2_grid(i,j,3)/=0.0)) THEN plate_boundary_area_m2 = plate_boundary_area_m2 + cell_area_m2(i) ELSE intraplate_area_m2 = intraplate_area_m2 + cell_area_m2(i) END IF END DO END DO percent_boundary_area = 100. * plate_boundary_area_m2 / sphere_area_m2 percent_intraplate_area = 100. * intraplate_area_m2 / sphere_area_m2 !Count shallow earthquakes in calibration period (both boundary & intraplate), !at or above the catalog (not forecast) threshold magnitude, !and prepare to save catalog earthquakes designated "intraplate" to a new .EQC file: boundary_calibration_count = 0 intraplate_calibration_count = 0 OPEN (UNIT = 2, FILE = TRIM(EQC_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios) new_EQC_file = "SHIFT_GSRM2"//alphabet(iVersion)//"_intraplate.eqc" OPEN (UNIT = 20, FILE = new_EQC_file) reticking: DO READ (2, 100, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data ! Note: PAD = "YES" permits READ with no appended data IF (ios /= 0) EXIT reticking IF (eq_mag >= catalog_threshold_magnitude) THEN i = 1 + NINT((GSRM2_lat_max - eq_Nlat) / GSRM2_d_lat) IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 j = 1 + NINT((eq_Elon - GSRM2_lon_min) / GSRM2_d_lon) IF ((smoothed_GSRM2_grid(i,j,1)/=0.0).OR.(smoothed_GSRM2_grid(i,j,2)/=0.0).OR.(smoothed_GSRM2_grid(i,j,3)/=0.0)) THEN boundary_calibration_count = boundary_calibration_count + 1 ELSE intraplate_calibration_count = intraplate_calibration_count + 1 WRITE (20, 100) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data END IF END IF END DO reticking CLOSE (2) ! EQC_file, for calibration CLOSE (20) ! new EQC file with events (above catalog threshold) selected as "intraplate" percent_boundary_eqs = 100. * boundary_calibration_count / & & (boundary_calibration_count + intraplate_calibration_count) percent_intraplate_eqs = 100. * intraplate_calibration_count / & & (boundary_calibration_count + intraplate_calibration_count) WRITE (*, *) WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate boundaries: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, boundary_calibration_count, percent_boundary_eqs WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs intraplate_seismicity_SI = intraplate_calibration_count / & & (intraplate_area_m2 * catalog_years * s_per_year) WRITE (*, "(' Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI WRITE (21, *) WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate boundaries: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, boundary_calibration_count, percent_boundary_eqs WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs WRITE (21, "('Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI !First-pass analysis of boundaries: Treat all as subduction zones !(N.B. This will be checked and corrected by a common factor, later on.) tentative_boundary_eqs_per_s = 0.0D0 DO i = 1, nRows DO j = 1, nColumns exx = smoothed_GSRM2_grid(i, j, 1) eyy = smoothed_GSRM2_grid(i, j, 2) exy = smoothed_GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value T5_column = 10 ! SUB CALL SHIFT_continuum_seismicity (catalog_threshold_moment_Nm, & ! INTENT(IN) & e1h, e2h, err, T5_column, & ! INTENT(IN) & seismicity) ! INTENT(OUT) tentative_boundary_eqs_per_s = tentative_boundary_eqs_per_s + seismicity * cell_area_m2(i) !Note that this summation uses the catalog (not forecast) threshold magnitude. !Also, it does not include the underlying "floor" of intraplate_seismicity_SI which is everywhere. END DO END DO tentative_boundary_eqs_per_catalog = tentative_boundary_eqs_per_s * catalog_years * s_per_year boundary_scaling_factor = ( boundary_calibration_count - & & (plate_boundary_area_m2 * intraplate_seismicity_SI * catalog_years * s_per_year) ) / & & tentative_boundary_eqs_per_catalog WRITE (*, *) WRITE (*, "(' Global plate-boundary seismicity scaling = ',F7.5,' (relative to SUB).')") boundary_scaling_factor WRITE (21, *) WRITE (21, "('Global plate-boundary seismicity scaling = ',F7.5,' (relative to SUB).')") boundary_scaling_factor !Scale the plate-boundary seismicity, and impose intraplate value as floor (both at forecast threshold): !First, determine how much intraplate seismicity shifts at forecast threshold (relative to catalog): !Use intraplate corner magnitude of 9.0 and beta = 0.63 inferred by Bird et al. [2010, SRL]. intraplate_forecast_over_catalog_factor = (( forecast_threshold_moment_Nm / catalog_threshold_moment_Nm )**(-0.63)) * & & EXP(( catalog_threshold_moment_Nm - forecast_threshold_moment_Nm ) / & & Moment(9.0) ) !using equation (7) in Bird et al. [2010, SRL]. intraplate_seismicity_SI = intraplate_seismicity_SI * intraplate_forecast_over_catalog_factor DO i = 1, nRows DO j = 1, nColumns exx = smoothed_GSRM2_grid(i, j, 1) eyy = smoothed_GSRM2_grid(i, j, 2) exy = smoothed_GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value T5_column = 10 ! SUB CALL SHIFT_continuum_seismicity (forecast_threshold_moment_Nm, & ! INTENT(IN) & e1h, e2h, err, T5_column, & ! INTENT(IN) & seismicity) ! INTENT(OUT) epicenter_rate_density(i, j) = boundary_scaling_factor * seismicity + intraplate_seismicity_SI IF ((exx /= 0.0).OR.(eyy /= 0.0).OR.(exy /= 0.0)) THEN effective_cz_meters(i, j) = NINT(boundary_scaling_factor * T5_coupledThickness_m(T5_column)) mean_G(i, j) = T5_assumed_mu_Pa(T5_column) END IF END DO END DO ELSE ! iVersion >= 4, so use tectonic zones [5 DOF] ------------------------------------------- !Get tectonic zones grid from file of Bird et al. [2010, Pure Appl. Geoph.]" WRITE (21, *) WRITE (21, "('File PB2002_tectonic_zones.grd must be available. Please check this.')") OPEN (UNIT = 1, FILE = "PB2002_tectonic_zones.grd", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Copy file PB2002_tectonic_zones.grd was not found in active folder.')") CALL Pause() STOP END IF READ (1, *) TZ_lon_min, TZ_d_lon, TZ_lon_max ! 0.0 0.1 360.0 (?) READ (1, *) TZ_lat_min, TZ_d_lat, TZ_lat_max ! -90.0 0.1 90.0 (?) CLOSE (1) nTZ_rows = NINT(1 + (TZ_lat_max - TZ_lat_min) / TZ_d_lat) nTZ_columns = NINT(1 + (TZ_lon_max - TZ_lon_min) / TZ_d_lon) ALLOCATE ( tectonic_zones_I1(nTZ_rows, nTZ_columns) ) OPEN (UNIT = 1, FILE = "PB2002_tectonic_zones.grd", STATUS = "OLD", IOSTAT = ios) READ (1, *) READ (1, *) READ (1, *)((tectonic_zones_I1(i, j), j = 1, nTZ_columns), i = 1, NTZ_rows) CLOSE (1) WRITE (*, "(' Tectonic zones grid has ',I6,' rows of ',I6' columns.')") nTZ_rows, nTZ_columns WRITE (*, "(' and was stored as INTEGER*1 to keep memory down to ~6 MB.')") WRITE (21, "('Tectonic zones grid has ',I6,' rows of ',I6' columns.')") nTZ_rows, nTZ_columns WRITE (21, "('and was stored as INTEGER*1 to keep memory down to ~6 MB.')") !Transfer to smoothed_GSRM2-sized grid, using zone point closest to cell center: ALLOCATE ( smoothed_GSRM2_tectonic_zone_I1(nRows, nColumns) ) DO i = 1, nRows lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat iTarget = NINT(1 + (TZ_lat_max - lat) / TZ_d_lat) DO j = 1, nColumns lon = GSRM2_lon_min + (j-1) * GSRM2_d_lon getLon = lon; IF (getLon < TZ_lon_min) getLon = getLon + 360.0 jTarget = NINT(1 + (getLon - TZ_lon_min) / TZ_d_lon) smoothed_GSRM2_tectonic_zone_I1(i, j) = tectonic_zones_I1(iTarget, jTarget) END DO END DO WRITE (*, "(' Finished selecting tectonic zones for GSRM2 grid from finer input grid.')") WRITE (21, "('Finished selecting tectonic zones for GSRM2 grid from finer input grid.')") !Assign tectonic zones to any active cells currently in zone 0, by referring to neighbors: WRITE (*, *) WRITE (*, "(' Inferring zone #s (> 0) for active cells in zone 0, by referring to neighbors.')") WRITE (21, *) WRITE (21, "('Inferring zone #s (> 0) for active cells in zone 0, by referring to neighbors.')") ALLOCATE ( snapshot_of_zones(nRows, nColumns) ) fiddling: DO iteration = 1, 1000 snapshot_of_zones = smoothed_GSRM2_tectonic_zone_I1 ! copy whole global matrix matrix nProblems = 0 ! just initializing new_assigns = 0 ! just initializing DO i = 1, nRows DO j = 1, nColumns IF (smoothed_GSRM2_tectonic_zone_I1(i, j) == 0) THEN IF ((smoothed_GSRM2_grid(i, j, 1) /= 0.).OR. & & (smoothed_GSRM2_grid(i, j, 2) /= 0.).OR. & & (smoothed_GSRM2_grid(i, j, 3) /= 0.)) THEN nProblems = nProblems + 1 neighbors_in_zone = 0 ! initializing whole array (4-vector, for the positive zone codes, 1:4) best_support = 0 ! initializing whole array (4-vector, for the positive zone codes, 1:4) !check neighbor to North: IF (i > 1) THEN it = snapshot_of_zones(i-1, j) IF (it > 0) THEN neighbors_in_zone(it) = neighbors_in_zone(it) + 1 best_support(it) = MAX(best_support(it), Get_Support(snapshot_of_zones, i-1, j, nRows, nColumns)) END IF END IF !check neighbor to South: IF (i < nColumns) THEN it = snapshot_of_zones(i+1, j) IF (it > 0) THEN neighbors_in_zone(it) = neighbors_in_zone(it) + 1 best_support(it) = MAX(best_support(it), Get_Support(snapshot_of_zones, i+1, j, nRows, nColumns)) END IF END IF !check neighbor to West: jTarget = j - 1; IF (jTarget < 1) jTarget = jTarget + nColumns it = snapshot_of_zones(i, jTarget) IF (it > 0) THEN neighbors_in_zone(it) = neighbors_in_zone(it) + 1 best_support(it) = MAX(best_support(it), Get_Support(snapshot_of_zones, i, jTarget, nRows, nColumns)) END IF !check neighbor to East: jTarget = j + 1; IF (jTarget > nColumns) jTarget = jTarget - nColumns it = snapshot_of_zones(i, jTarget) IF (it > 0) THEN neighbors_in_zone(it) = neighbors_in_zone(it) + 1 best_support(it) = MAX(best_support(it), Get_Support(snapshot_of_zones, i, jTarget, nRows, nColumns)) END IF fixing: DO neighbor_count = 4, 1, -1 ! preference is strongest for large # of (identical) neighbors DO support = 8, 0, -1 ! preferring potential parents with greater local support top_zone = 0 ! just initializing, as search loop "short" below will often fail short: DO k = 4, 1, -1 ! trial zone # !Note that search order implies a very subtle preference for parents at the active 4/3 end, not the quiet 1/2 end. !This is because propagation of zone 1 (active continent) away from tiny island spots is problematic, if encouraged. IF (best_support(k) == support) THEN top_zone = k EXIT short END IF END DO short IF (top_zone > 0) THEN ! some zone had the currently-expected support level IF (neighbors_in_zone(top_zone) == neighbor_count) THEN ! got an assignment! new_assigns = new_assigns + 1 smoothed_GSRM2_tectonic_zone_I1(i, j) = top_zone EXIT fixing ! and, this inner loop as well! END IF END IF END DO END DO fixing END IF ! problem cell; active! END IF ! cell currently assigned to zone 0; possible problem? END DO ! loop on columns END DO ! loop on rows WRITE (*, "(' Iteration ',I3,': ',I6,' problems; ',I6,' assigned (through common edges).')") iteration, nProblems, new_assigns WRITE (21, "('Iteration ',I3,': ',I6,' problems; ',I6,' assigned (through common edges).')") iteration, nProblems, new_assigns IF (nProblems == 0) EXIT fiddling IF (new_assigns == 0) EXIT fiddling END DO fiddling ! iteration = 1, 1000 DEALLOCATE ( snapshot_of_zones ) !GPBrezoning: !Apply patches to region of western North America where auto-zoning gives BAD results! DO i = 1, nRows lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat DO j = 1, nColumns lon = GSRM2_lon_min + (j-1) * GSRM2_d_lon ! -180 to +180 !Borderland of California should be "active continent", not "fast ridge": IF ((lat > 28.).AND.(lat < 40.).AND.(lon > -131.).AND.(lon < -115.)) THEN IF (smoothed_GSRM2_tectonic_zone_I1(i, j) /= 0) smoothed_GSRM2_tectonic_zone_I1(i, j) = 1 END IF !"Subduction zone" should include Mount Lassen and Mount Shasta in CA: IF ((lat >= 40.).AND.(lat <= 46.).AND.(lon > -126.).AND.(lon < -121.)) THEN IF (smoothed_GSRM2_tectonic_zone_I1(i, j) /= 0) smoothed_GSRM2_tectonic_zone_I1(i, j) = 4 END IF !"Subduction zone" should NOT extend further inland than the Cascade volcanic arc: IF ((lat >= 40.).AND.(lat <= 50.).AND.(lon > -121.).AND.(lon < -115.)) THEN IF (smoothed_GSRM2_tectonic_zone_I1(i, j) /= 0) smoothed_GSRM2_tectonic_zone_I1(i, j) = 1 END IF IF ((lat >= 44.).AND.(lat <= 54.).AND.(lon > (-122. - 1.5 * (lat - 50.))).AND.(lon < -115.)) THEN IF (smoothed_GSRM2_tectonic_zone_I1(i, j) /= 0) smoothed_GSRM2_tectonic_zone_I1(i, j) = 1 END IF END DO END DO !Set iZone to 0 for cells with zero strain-rate (wherever PB2002_tectonic_zones were broader than GSRM2, e.g., NE of Baikal). DO i = 1, nRows DO j = 1, nColumns IF ((smoothed_GSRM2_grid(i, j, 1) == 0.0).AND.(smoothed_GSRM2_grid(i, j, 2) == 0.0).AND.(smoothed_GSRM2_grid(i, j, 3) == 0.0)) THEN smoothed_GSRM2_tectonic_zone_I1(i, j) = 0 ! signifying "intraplate" END IF END DO END DO !Eliminate orphan cells (spatially isolated, probably due to sandblasting of raw SUBs). IF (nProblems > 0) THEN WRITE (*, "(' Now giving up and zeroing strain rates in ',I6,' orphan cells.')") nProblems WRITE (21, "('Now giving up and zeroing strain rates in ',I6,' orphan cells.')") nProblems biggest_eRate = 0.0 DO i = 1, nRows DO j = 1, nColumns IF (smoothed_GSRM2_tectonic_zone_I1(i, j) == 0) THEN biggest_eRate = MAX(biggest_eRate, ABS(smoothed_GSRM2_grid(i, j, 1)), & & ABS(smoothed_GSRM2_grid(i, j, 2)), & & ABS(smoothed_GSRM2_grid(i, j, 3))) smoothed_GSRM2_grid(i, j, 1:3) = 0.0 END IF END DO END DO biggest_n_per_year = biggest_eRate * 1.0E9 * 3.1536E7 WRITE (*, "(' Largest strain rate zeroed out was ',ES10.2,'/s.(',F6.1,' nanostrain/year).')") biggest_eRate, biggest_n_per_year WRITE (21, "('Largest strain rate zeroed out was ',ES10.2,'/s (',F6.1,' nanostrain/year).')") biggest_eRate, biggest_n_per_year END IF ! orphan cells need to be zeroed !OPEN (UNIT = 10, FILE = "SHIFT_GSRM2d_tectonic_zones.grd") !WRITE (10, "(3F10.2)") GSRM2_lon_min, GSRM2_d_lon, GSRM2_lon_max !WRITE (10, "(3F10.2)") GSRM2_lat_min, GSRM2_d_lat, GSRM2_lat_max !WRITE (10, "(30I2)") ((smoothed_GSRM2_tectonic_zone_I1(i, j), j = 1, nColumns), i = 1, nRows) !CLOSE (10) !WRITE (*, *) !WRITE (*, "(' GSRM2d_tectonic_zones.grd file was written.')") ! !DO i = 1, nRows ! DO j = 1, nColumns ! epicenter_rate_density(i, j) = 10.0**(smoothed_GSRM2_tectonic_zone_I1(i, J) - 20.0) ! END DO !END DO !CALL Pause() IF (iVersion >= 5) THEN !Find PB2002 plate boundary nearest to each newly-active (due to spreading) smoothed_GSRM2 grid cell, !but store it in the unused cells of the already-existing arrays: WRITE (*, *) WRITE (*, "(' Finding PB2002 plate boundary nearest to each newly-active GSRM2 grid cell...')") WRITE (21, *) WRITE (21, "('Finding PB2002 plate boundary nearest to each newly-active GSRM2 grid cell...')") old_percent_done = 0 DO i = 1, nRows lat = GSRM2_lat_max - (i-1) * GSRM2_d_lat DO j = 1, nColumns IF ((smoothed_GSRM2_tectonic_zone_I1(i, j) > 0).AND.(GSRM2_closest_PB2002_code(i, j) == 0)) THEN !conduct a search: lon = GSRM2_lon_min + (j-1) * GSRM2_d_lon CALL LonLat_2_Uvec(lon, lat, uvec1) standard = 1.0 ! not close at all! DO k = 1, PB2002_list_highest uvec2(1:3) = PB2002_uvecs(1:3, k) arc_radian2 = (uvec1(1) - uvec2(1))**2 + & & (uvec1(2) - uvec2(2))**2 + & & (uvec1(3) - uvec2(3))**2 IF (arc_radian2 < standard) THEN standard = arc_radian2 best_code = PB2002_code(k) best_elevation = PB2002_elevation(k) best_azimuth = PB2002_seaward_azimuth(k) best_velocity_mmpa = PB2002_velocity_mmpa(k) END IF END DO GSRM2_closest_PB2002_radian2(i, j) = standard GSRM2_closest_PB2002_code(i, j) = best_code GSRM2_closest_PB2002_elevation(i, j) = best_elevation GSRM2_closest_PB2002_velocity_mmpa(i, j) = best_velocity_mmpa GSRM2_closest_PB2002_seaward_azimuth(i, j) = best_azimuth ! INTEGER*2, and 0 for non-SUB GSRM2_raw_SUB(i, j) = .FALSE. END IF END DO percent_done = (100. * i) / nRows IF (percent_done > old_percent_done) THEN WRITE (*, "(' ',I3,'% done...')") percent_done old_percent_done = percent_done END IF END DO WRITE (*, "(' Search completed.')") WRITE (21, "('Search completed.')") END IF ! iVersion >= 5, so PB2002 data now needed for additional cells made active by spread strain rates. !In cases of iVersion >= 4 (4, 5, 6), !Re-compute total areas of straining plate boundaries (in each zone), and of intraplate regions: plate_boundary_area_m2_in_zone = 0.0D0 ! whole 4-vector intraplate_area_m2 = 0.0D0 DO i = 1, nRows DO j = 1, nColumns IF ((smoothed_GSRM2_grid(i,j,1)/=0.0).OR.(smoothed_GSRM2_grid(i,j,2)/=0.0).OR.(smoothed_GSRM2_grid(i,j,3)/=0.0)) THEN iZone = smoothed_GSRM2_tectonic_zone_I1(i, j) plate_boundary_area_m2_in_zone(iZone) = plate_boundary_area_m2_in_zone(iZone) + cell_area_m2(i) ELSE intraplate_area_m2 = intraplate_area_m2 + cell_area_m2(i) END IF END DO END DO WRITE (*, *) WRITE (21, *) DO k = 1, 4 percent_boundary_area = 100. * plate_boundary_area_m2_in_zone(k) / sphere_area_m2 WRITE (*, "(' Tectonic zone ',I2,': area ',ES12.4,' m^2 (',F6.2,'% of Earth)')") & & k, plate_boundary_area_m2_in_zone(k), percent_boundary_area WRITE (21, "('Tectonic zone ',I2,': area ',ES12.4,' m^2 (',F6.2,'% of Earth)')") & & k, plate_boundary_area_m2_in_zone(k), percent_boundary_area END DO percent_intraplate_area = 100. * intraplate_area_m2 / sphere_area_m2 WRITE (*, "(' Intraplate : area 'ES12.4,' m^2 (',F6.2,'% of Earth)')") & & intraplate_area_m2, percent_intraplate_area WRITE (21, "('Intraplate : area 'ES12.4,' m^2 (',F6.2,'% of Earth)')") & & intraplate_area_m2, percent_intraplate_area !Count shallow earthquakes in calibration period (both boundary/per-zone, & intraplate), !at or above the catalog (not forecast) threshold magnitude, !and prepare to save catalog earthquakes designated "intraplate" to a new .EQC file: boundary_calibration_count_in_zone = 0 ! whole 4-vector intraplate_calibration_count = 0 OPEN (UNIT = 2, FILE = TRIM(EQC_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios) new_EQC_file = "SHIFT_GSRM2"//alphabet(iVersion)//"_intraplate.eqc" OPEN (UNIT = 20, FILE = new_EQC_file) total_EQ_count_above_catalog_threshold = 0 ! just initializing before loop moreticking: DO READ (2, 100, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data ! Note: PAD = "YES" permits READ with no appended data IF (ios /= 0) EXIT moreticking IF (eq_mag >= catalog_threshold_magnitude) THEN total_EQ_count_above_catalog_threshold = total_EQ_count_above_catalog_threshold + 1 i = 1 + NINT((GSRM2_lat_max - eq_Nlat) / GSRM2_d_lat) IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((eq_Elon - GSRM2_central_lon) >= 360.0) eq_Elon = eq_Elon - 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 IF ((GSRM2_central_lon - eq_Elon) >= 360.0) eq_Elon = eq_Elon + 360.0 j = 1 + NINT((eq_Elon - GSRM2_lon_min) / GSRM2_d_lon) iZone = smoothed_GSRM2_tectonic_zone_I1(i, j) IF ((smoothed_GSRM2_grid(i,j,1)/=0.0).OR.(smoothed_GSRM2_grid(i,j,2)/=0.0).OR.(smoothed_GSRM2_grid(i,j,3)/=0.0)) THEN IF (iZone > 0) THEN ! expected boundary_calibration_count_in_zone(iZone) = boundary_calibration_count_in_zone(iZone) + 1 ELSE ! problem WRITE (*, "(' ERROR: Calibration earthquake at ',F9.3,' E, ',F7.3,' N')") eq_Elon, eq_Nlat WRITE (*, "(' falls into cell (',I6,', ',I6,')')") i, j WRITE (*, "(' which has non-zero strain-rate of: ',1P,3ES10.2)") smoothed_GSRM2_grid(i, j, 1:3) WRITE (*, "(' but has assigned iZone == 0.')") WRITE (21, "('ERROR: Calibration earthquake at ',F9.3,' E, ',F7.3,' N')") eq_Elon, eq_Nlat WRITE (21, "('falls into cell (',I6,', ',I6,')')") i, j WRITE (21, "('which has non-zero strain-rate of: ',1P,3ES10.2)") smoothed_GSRM2_grid(i, j, 1:3) WRITE (21, "('but has assigned iZone == 0.')") CALL Pause() STOP END IF ELSE intraplate_calibration_count = intraplate_calibration_count + 1 WRITE (20, 100) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data END IF END IF END DO moreticking CLOSE (2) ! EQC_file, for calibration CLOSE (20) ! new EQC file with events (above catalog threshold) selected as "intraplate" WRITE (*, *) WRITE (21, *) DO k = 1, 4 percent_boundary_eqs = 100. * boundary_calibration_count_in_zone(k) / total_EQ_count_above_catalog_threshold WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate boundary zone ',I2,': ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, k, boundary_calibration_count_in_zone(k), percent_boundary_eqs WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate boundary zone ',I2,': ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, k, boundary_calibration_count_in_zone(k), percent_boundary_eqs END DO percent_intraplate_eqs = 100. * intraplate_calibration_count / total_EQ_count_above_catalog_threshold WRITE (*, "(' Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs WRITE (21, "('Earthquakes of magnitude >= ',F6.3,' in plate interiors: ',I6,' (',F6.2,'%)')") & & catalog_threshold_magnitude, intraplate_calibration_count, percent_intraplate_eqs intraplate_seismicity_SI = intraplate_calibration_count / & & (intraplate_area_m2 * catalog_years * s_per_year) WRITE (*, "(' Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI WRITE (21, "('Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_SI !First-pass analyses of boundaries: Treat each zone as its dominant plate-boundary type: !(N.B. These will be checked and corrected by zone-specific correction factors, later on.) tentative_boundary_eqs_per_s_in_zone = 0.0D0 ! whole 4-vector DO i = 1, nRows DO j = 1, nColumns exx = smoothed_GSRM2_grid(i, j, 1) eyy = smoothed_GSRM2_grid(i, j, 2) exy = smoothed_GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive (negative) horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional (positive) horizontal principal value iZone = smoothed_GSRM2_tectonic_zone_I1(i, j) IF (iZone > 0) THEN !VVVVVV pass 1 VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV velocity_factor(1) = 1.0 ! and remains 1.0, except for SUB & CCB when iVersion >= 6 velocity_factor(2) = 1.0 ! and remains 1.0, except for SUB & CCB when iVersion >= 6 !Note: In present context, iVersion is >= 4. IF (iVersion == 4) THEN nTensorParts = 1 ! no partitioning if iVersion == 4 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! "Active continent" T5_column_list(1) = 2 ! CTF (intermediate between CRB and CCB) ELSE IF (iZone == 2) THEN ! "Slow ridge" T5_column_list(1) = 6 ! slow OTF (most of seismicity) ELSE IF (iZone == 3) THEN ! "Fast ridge" T5_column_list(1) = 8 ! fast OTF (most of seismicity) ELSE IF (iZone == 4) THEN ! "Trench" T5_column_list(1) = 10 ! SUB END IF ELSE ! iVersion >= 5; add classification & partitioning of strain rate tensors IF (iZone == 4) THEN ! no partitioning; only issue is: Whether SUB (&, fast or slow?), or OCB? nTensorParts = 1 ! no partitioning e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (GSRM2_closest_PB2002_code(i, j) == 7) THEN ! SUB is closest PB2002 boundary T5_column_list(1) = 10 ! SUB IF (iVersion >= 6) THEN ! add velocity-dependence of SUB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 66.0) THEN velocity_factor(1) = 1.206 ELSE velocity_factor(1) = 0.589 END IF END IF ELSE ! assume OCB: T5_column_list(1) = 9 ! OCB END IF ELSE ! iZone == 1, 2, 3; go ahead and partition IF (err > 0.0) THEN ! there is some thrust-faulting IF ((e1h * e2h) < 0.0) THEN ! there is also some strike-slip faulting: transpression nTensorParts = 2 ! partitioning !first part is the thrust-faulting component err_list(1) = err ! positive; all vertical strain is due to thrust-faulting e1h_list(1) = -err e2h_list(1) = 0.0 !second part is the strike-slip component err_list(2) = 0.0 e2h_list(2) = e2h e1h_list(2) = -e2h IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 3 ! CCB IF (iVersion >= 6) THEN ! add velocity-dependence of CCB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 24.2) THEN velocity_factor(1) = 1.383 ELSE velocity_factor(1) = 0.606 END IF END IF T5_column_list(2) = 2 ! CTF ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 9 ! OCB T5_column_list(2) = 6 ! OTF slow ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 9 ! OCB T5_column_list(2) = 8 ! OTF fast END IF ! iZone == 1, 2, or 3 ELSE ! e1h & e12 have same (negative) sign; only thrust faulting, in 2 conjugate sets. nTensorParts = 1 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 3 ! CCB IF (iVersion >= 6) THEN ! add velocity-dependence of CCB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 24.2) THEN velocity_factor(1:2) = 1.383 ELSE velocity_factor(1:2) = 0.606 END IF END IF ELSE IF (iZone == 2) THEN ! Slow ridge (including its transpressive transforms) T5_column_list(1) = 9 ! OCB ELSE IF (iZone == 3) THEN ! Fast ridge (including its transpressive transforms) T5_column_list(1) = 9 ! OCB END IF ! iZone == 1, 2, or 3 END IF ! any strike-slip faulting component, or not ELSE ! err <= 0.0; there is (usually) some normal-faulting IF ((e1h * e2h) < 0.0) THEN ! there is also some strike-slip faulting: transtension nTensorParts = 2 ! partitioning !first part is the normal-faulting component err_list(1) = err ! negative; all vertical strain is due to normal-faulting e1h_list(1) = 0.0 e2h_list(1) = -err ! positive !second part is the strike-slip component err_list(2) = 0.0 e2h_list(2) = -e1h ! positive e1h_list(2) = e1h ! negative IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 1 ! CRB T5_column_list(2) = 2 ! CTF ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 4 ! OSR normal T5_column_list(2) = 6 ! OTF slow ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 4 ! OSR normal T5_column_list(2) = 8 ! OTF fast END IF ! iZone == 1, 2, or 3 ELSE ! e1h & e12 have same (positive) sign; only normal faulting, in 2 conjugate sets. nTensorParts = 1 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 1 ! CRB ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 4 ! OSR normal ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 4 ! OSR normal END IF ! iZone == 1, 2, or 3 END IF ! any strike-slip faulting component, or not END IF ! err > 0.0 (thrust-faulting) or <= 0.0 (normal-faulting) END IF ! iZone == 4, OR {1, 2, 3}; no-partition? or partition? END IF ! iVersion = 4 (no classification or partitioning) or iVersion >= 5 (with classification and partitioning). !^^^^^^ pass 1 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ total_cell_seismicity = 0.0 ! initialize sum (of 1 or 2 terms) DO k = 1, nTensorParts CALL SHIFT_continuum_seismicity (catalog_threshold_moment_Nm, & ! INTENT(IN) & e1h_list(k), e2h_list(k), err_list(k), T5_column_list(k), & ! INTENT(IN) & seismicity) ! INTENT(OUT) total_cell_seismicity = total_cell_seismicity + seismicity * velocity_factor(k) END DO ! k = 1, nTensorParts tentative_boundary_eqs_per_s_in_zone(iZone) = tentative_boundary_eqs_per_s_in_zone(iZone) + & & total_cell_seismicity * cell_area_m2(i) !Note that this summation uses the catalog (not forecast) threshold magnitude. !Also, it does not include the underlying "floor" of intraplate_seismicity_SI which is everywhere. END IF ! iZone > 0 END DO ! loop j = 1, nColumns END DO ! loop i = 1, nRows tentative_boundary_eqs_per_catalog_in_zone = tentative_boundary_eqs_per_s_in_zone * catalog_years * s_per_year ! whole 4-vector WRITE (*, *) WRITE (21, *) DO iZone = 1, 4 boundary_scaling_factor_in_zone(iZone) = ( boundary_calibration_count_in_zone(iZone) - & & (plate_boundary_area_m2_in_zone(iZone) * intraplate_seismicity_SI * catalog_years * s_per_year) ) / & & tentative_boundary_eqs_per_catalog_in_zone(iZone) WRITE (*, "(' Global plate-boundary seismicity scaling factor in zone ',I2,' = ',F7.5,'.')") iZone, boundary_scaling_factor_in_zone(iZone) WRITE (21, "('Global plate-boundary seismicity scaling factor in zone ',I2,' = ',F7.5,'.')") iZone, boundary_scaling_factor_in_zone(iZone) END DO !Scale the plate-boundary seismicities (per zone), and impose intraplate value as floor (both at forecast threshold): !First, determine how much intraplate seismicity shifts at forecast threshold (relative to catalog): !Use intraplate corner magnitude of 9.0 and beta = 0.63 inferred by Bird et al. [2010, SRL]. intraplate_forecast_over_catalog_factor = (( forecast_threshold_moment_Nm / catalog_threshold_moment_Nm )**(-0.63)) * & & EXP(( catalog_threshold_moment_Nm - forecast_threshold_moment_Nm ) / & & Moment(9.0) ) !using equation (7) in Bird et al. [2010, SRL]. intraplate_seismicity_SI = intraplate_seismicity_SI * intraplate_forecast_over_catalog_factor epicenter_rate_density = intraplate_seismicity_SI ! whole array; only iZone>0 points will be revisited: DO i = 1, nRows DO j = 1, nColumns exx = smoothed_GSRM2_grid(i, j, 1) eyy = smoothed_GSRM2_grid(i, j, 2) exy = smoothed_GSRM2_grid(i, j, 3) !Characterize the strain-rate tensor of this cell: second_invariant = SQRT(exx**2 + eyy**2 + 2.0 *(exy**2)) ! (per Kreemer et al. [2002]; used in test plots) dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value iZone = smoothed_GSRM2_tectonic_zone_I1(i, j) IF (iZone > 0) THEN !VVVVVV pass 2 VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV velocity_factor(1) = 1.0 ! and remains 1.0, except for SUB & CCB when iVersion >= 6 velocity_factor(2) = 1.0 ! and remains 1.0, except for SUB & CCB when iVersion >= 6 !Note: In present context, iVersion is >= 4. IF (iVersion == 4) THEN nTensorParts = 1 ! no partitioning if iVersion == 4 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! "Active continent" T5_column_list(1) = 2 ! CTF (intermediate between CRB and CCB) ELSE IF (iZone == 2) THEN ! "Slow ridge" T5_column_list(1) = 6 ! slow OTF (most of seismicity) ELSE IF (iZone == 3) THEN ! "Fast ridge" T5_column_list(1) = 8 ! fast OTF (most of seismicity) ELSE IF (iZone == 4) THEN ! "Trench" T5_column_list(1) = 10 ! SUB END IF ELSE ! iVersion >= 5; add classification & partitioning of strain rate tensors IF (iZone == 4) THEN ! no partitioning; only issue is: Whether SUB (&, fast or slow?), or OCB? nTensorParts = 1 ! no partitioning e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (GSRM2_closest_PB2002_code(i, j) == 7) THEN ! SUB is closest PB2002 boundary T5_column_list(1) = 10 ! SUB IF (iVersion >= 6) THEN ! add velocity-dependence of SUB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 66.0) THEN velocity_factor(1) = 1.206 ELSE velocity_factor(1) = 0.589 END IF END IF ELSE ! assume OCB: T5_column_list(1) = 9 ! OCB END IF ELSE ! iZone == 1, 2, 3; go ahead and partition IF (err > 0.0) THEN ! there is some thrust-faulting IF ((e1h * e2h) < 0.0) THEN ! there is also some strike-slip faulting: transpression nTensorParts = 2 ! partitioning !first part is the thrust-faulting component err_list(1) = err ! positive; all vertical strain is due to thrust-faulting e1h_list(1) = -err e2h_list(1) = 0.0 !second part is the strike-slip component err_list(2) = 0.0 e2h_list(2) = e2h e1h_list(2) = -e2h IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 3 ! CCB IF (iVersion >= 6) THEN ! add velocity-dependence of CCB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 24.2) THEN velocity_factor(1) = 1.383 ELSE velocity_factor(1) = 0.606 END IF END IF T5_column_list(2) = 2 ! CTF ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 9 ! OCB T5_column_list(2) = 6 ! OTF slow ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 9 ! OCB T5_column_list(2) = 8 ! OTF fast END IF ! iZone == 1, 2, or 3 ELSE ! e1h & e12 have same (negative) sign; only thrust faulting, in 2 conjugate sets. nTensorParts = 1 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 3 ! CCB IF (iVersion >= 6) THEN ! add velocity-dependence of CCB coupling [Bird et al., 2009, BSSA] IF (GSRM2_closest_PB2002_velocity_mmpa(i, j) > 24.2) THEN velocity_factor(1:2) = 1.383 ELSE velocity_factor(1:2) = 0.606 END IF END IF ELSE IF (iZone == 2) THEN ! Slow ridge (including its transpressive transforms) T5_column_list(1) = 9 ! OCB ELSE IF (iZone == 3) THEN ! Fast ridge (including its transpressive transforms) T5_column_list(1) = 9 ! OCB END IF ! iZone == 1, 2, or 3 END IF ! any strike-slip faulting component, or not ELSE ! err <= 0.0; there is (usually) some normal-faulting IF ((e1h * e2h) < 0.0) THEN ! there is also some strike-slip faulting: transtension nTensorParts = 2 ! partitioning !first part is the normal-faulting component err_list(1) = err ! negative; all vertical strain is due to normal-faulting e1h_list(1) = 0.0 e2h_list(1) = -err ! positive !second part is the strike-slip component err_list(2) = 0.0 e2h_list(2) = -e1h ! positive e1h_list(2) = e1h ! negative IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 1 ! CRB T5_column_list(2) = 2 ! CTF ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 4 ! OSR normal T5_column_list(2) = 6 ! OTF slow ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 4 ! OSR normal T5_column_list(2) = 8 ! OTF fast END IF ! iZone == 1, 2, or 3 ELSE ! e1h & e12 have same (positive) sign; only normal faulting, in 2 conjugate sets. nTensorParts = 1 e1h_list(1) = e1h e2h_list(1) = e2h err_list(1) = err IF (iZone == 1) THEN ! Active continent T5_column_list(1) = 1 ! CRB ELSE IF (iZone == 2) THEN ! Slow ridge T5_column_list(1) = 4 ! OSR normal ELSE IF (iZone == 3) THEN ! Fast ridge T5_column_list(1) = 4 ! OSR normal END IF ! iZone == 1, 2, or 3 END IF ! any strike-slip faulting component, or not END IF ! err > 0.0 (thrust-faulting) or <= 0.0 (normal-faulting) END IF ! iZone == 4, OR {1, 2, 3}; no-partition? or partition? END IF ! iVersion = 4 (no classification or partitioning) or iVersion >= 5 (with classification and partitioning). !^^^^^^ pass 2 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ total_cell_seismicity = 0.0 ! initialize, before sum of 1 or 2 terms effective_cz_numerator = 0.0 ! initialize, before sum of 1 or 2 terms effective_cz_denominator = 0.0 ! initialize, before sum of 1 or 2 terms mean_G_numerator = 0.0 ! initialize, before sum of 1 or 2 terms mean_G_denominator = 0.0 ! initialize, before sum of 1 or 2 terms DO k = 1, nTensorParts CALL SHIFT_continuum_seismicity (forecast_threshold_moment_Nm, & ! INTENT(IN) & e1h_list(k), e2h_list(k), err_list(k), T5_column_list(k), & ! INTENT(IN) & seismicity) ! INTENT(OUT) total_cell_seismicity = total_cell_seismicity + seismicity * velocity_factor(k) e1_rate = MIN(e1h_list(k), e2h_list(k), err_list(k)) ! always negative e3_rate = MAX(e1h_list(k), e2h_list(k), err_list(k)) ! always positive e2_rate = 0.0 - e1_rate - e3_rate ! may have either sign IF (e2_rate < 0.0) THEN ! both e1_rate and e2_rate are negative; positive e3_rate is the unique axis: scalar_strainrate = e3_rate ELSE ! e2_rate >= 0.0 ; so both e2_rate and e3_rate are positive; negative e1_rate is the unique axis: scalar_strainrate = (-e1_rate) END IF effective_cz_numerator = effective_cz_numerator + T5_coupledThickness_m(T5_column_list(k)) * & & boundary_scaling_factor_in_zone(iZone) * & & velocity_factor(k) * & & scalar_strainrate effective_cz_denominator = effective_cz_denominator + scalar_strainrate mean_G_numerator = mean_G_numerator + ( T5_assumed_mu_Pa(T5_column_list(k)) * scalar_strainrate ) mean_G_denominator = mean_G_denominator + scalar_strainrate !save #s of columns used in this partitioning: jT5column(k, i, j) = T5_column_list(k) !save moment-rate/area from each of (1 or 2) partitioned components: (To be normalized later.) fractionT5column(k, i, j) = T5_coupledThickness_m(T5_column_list(k)) * & & T5_assumed_mu_Pa(T5_column_list(k)) * & & 2.0 * scalar_strainrate * & ! (N.B. Factor "* 2.0" could be omitted without changing output.) & boundary_scaling_factor_in_zone(iZone) * & ! (N.B. This factor could be omitted without changing output.) & velocity_factor(k) END DO ! k = 1, nTensorParts epicenter_rate_density(i, j) = boundary_scaling_factor_in_zone(iZone) * total_cell_seismicity + & & intraplate_seismicity_SI IF ((exx /= 0.0).OR.(eyy /= 0.0).OR.(exy /= 0.0)) THEN IF (effective_cz_denominator > 0.0) THEN effective_cz_meters(i, j) = NINT(effective_cz_numerator / effective_cz_denominator) END IF IF (mean_G_denominator > 0.0) THEN mean_G(i, j) = mean_G_numerator / mean_G_denominator END IF END IF !normalize components of partitioned moment rate and save as fractions: M_rate_per_m2 = fractionT5column(1, i, j) + fractionT5column(2, i, j) IF (M_rate_per_m2 > 0.0) THEN fractionT5column(1, i, j) = fractionT5column(1, i, j) / M_rate_per_m2 ! now, 0.0 ~ 1.0 fractionT5column(2, i, j) = fractionT5column(2, i, j) / M_rate_per_m2 ! now, 0.0 ~ 1.0 END IF END IF ! iZone > 0 END DO ! loop j = 1, nColumns END DO ! loop i = 1, nRows END IF ! iVersion == (2,3) or iVersion >= 4; without/with tectonic zones [2 DOF, or 5 DOF] ---- END IF ! iVersion == 1, or iVersion >= 2; without/with spatial smoothing -------------------------- !Output the seismicity (EQs/m^2/s) .grd file: WRITE (*, *) WRITE (*, "(' Writing seismicity (EQs/m^2/s) .grd file and computing area-integral...')") WRITE (21, *) WRITE (21, "('Writing seismicity (EQs/m^2/s) .grd file and computing area-integral...')") OPEN (UNIT = 3, FILE = TRIM(GRD_file), IOSTAT = ios) WRITE (3, "(3F11.5,' = lon_min, d_lon, lon_max Gridded long-term seismicity above magnitude ',F6.3,' in events/square-meter/s')") & & GSRM2_lon_min, GSRM2_d_lon, GSRM2_lon_max, forecast_threshold_magnitude WRITE (3, "(3F11.5,' = lat_min, d_lat, lat_max from program SHIFT_GSRM2x, P. Bird, UCLA, 2014')") & & GSRM2_lat_min, GSRM2_d_lat, GSRM2_lat_max WRITE (3, "(1P,10D10.2)") ((epicenter_rate_density(i, j), j = 1, nColumns), i = 1, nRows) !Note: Line breaks in output file have no relation to the logical shape of this array. ! Programs that read .grd files must take this into account, by reading the whole ! array with a single READ (as it was written), and not row-by-row. !----------------------------------------------------------------------------------- !Integrate global shallow seismicity: !Integrate seismicity over grid area: integral = 0.0D0 ! (just initializing) DO i = 1, nRows DO j = 1, nColumns integral = integral + epicenter_rate_density(i, j) * cell_area_m2(i) END DO END DO integral_converted = integral * 100.0 * s_per_year ! now, earthquakes per century WRITE (*, "(/' Area integral of seismicity =' & & /' ',ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted WRITE (3, "( 'Area integral of seismicity =' & & /ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted WRITE (21, "(/'Area integral of seismicity =' & & /ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted CLOSE (UNIT = 3, DISP = "KEEP") ! seismicity.grd file !Output the effective_cz_meters (= of Bird & Kagan [2004] * retrospective adjustment for zone) .grd file: WRITE (*, *) WRITE (*, "(' Writing effective_cz_meters .grd file...')") WRITE (21, *) WRITE (21, "('Writing effective_cz_meters .grd file...')") GRD_file = "SHIFT_GSRM2" // alphabet(iVersion) // "_effective_cz_meters.grd" OPEN (UNIT = 3, FILE = TRIM(GRD_file), IOSTAT = ios) WRITE (3, "(3F11.5,' = lon_min, d_lon, lon_max Effective , = Bird & Kagan [2004] * zonal retrospective adjustment')") & & GSRM2_lon_min, GSRM2_d_lon, GSRM2_lon_max WRITE (3, "(3F11.5,' = lat_min, d_lat, lat_max from program SHIFT_GSRM2x, P. Bird, UCLA, 2014')") & & GSRM2_lat_min, GSRM2_d_lat, GSRM2_lat_max WRITE (3, "(12I6)") ((effective_cz_meters(i, j), j = 1, nColumns), i = 1, nRows) !Note: Line breaks in output file have no relation to the logical shape of this array. ! Programs that read .grd files must take this into account, by reading the whole ! array with a single READ (as it was written), and not row-by-row. !----------------------------------------------------------------------------------- CLOSE (UNIT = 3, DISP = "KEEP") ! effective_cz_meters .grd file (for plotting with NeoKinema) !GPBhere !Output non-zero entries from effective_cz_meters in .dat format (for Corne Kreemer), including smoothed strain-rates & elastic shear moduli: WRITE (*, *) WRITE (*, "(' Writing moment_factors_in_SI.dat file...')") WRITE (21, *) WRITE (21, "('Writing moment_factors_in_SI.dat file...')") DAT_file = "SHIFT_GSRM2" // alphabet(iVersion) // "_moment_factors_in_SI.dat" OPEN (UNIT = 3, FILE = TRIM(DAT_file), IOSTAT = ios) total_moment_rate_SI = 0.0 ! initializing before sum DO i = nRows, 1, -1 lat = GSRM2_lat_max - (i - 1) * GSRM2_d_lat DO j = 1, nColumns lon = GSRM2_lon_min + (j - 1) * GSRM2_d_lon IF (effective_cz_meters(i, j) /= 0) THEN IF (iVersion == 1) THEN ! use original values exx = GSRM2_grid(i, j, 1) / strain_rate_units eyy = GSRM2_grid(i, j, 2) / strain_rate_units exy = GSRM2_grid(i, j, 3) / strain_rate_units ELSE ! use smoothed values exx = smoothed_GSRM2_grid(i, j, 1) / strain_rate_units eyy = smoothed_GSRM2_grid(i, j, 2) / strain_rate_units exy = smoothed_GSRM2_grid(i, j, 3) / strain_rate_units END IF !-------------------------------------------------------------------------------------------------------------------------------------- WRITE (3, "(F7.3, F9.3, 3F8.1, I6, ES10.3, I3, F5.2, I3, F5.2)") lat, lon, exx, eyy, exy, effective_cz_meters(i, j), mean_G(i, j), & & jT5column(1, i, j), fractionT5column(1, i, j), & ! e.g., " 10 1.00" & jT5column(2, i, j), fractionT5column(2, i, j) ! <-- may be " 0 0.00" !-------------------------------------------------------------------------------------------------------------------------------------- !repeat computation of scalar_strainrate, as above: exx = exx * strain_rate_units ! returning them to SI, /s units eyy = eyy * strain_rate_units exy = exy * strain_rate_units dilatation_rate = exx + eyy err = -dilatation_rate ! based on assumption of volume conservation in permanent, non-elastic strain !N.B. As strain-rates of GSRM2 may include large elastic parts, this is an APPROXIMATION when used here! center_normal_rate = (exx + eyy)/2.0 radius_rate = SQRT((exx - center_normal_rate)**2 + exy**2) e1h = center_normal_rate - radius_rate ! most-compressive horizontal principal value e2h = center_normal_rate + radius_rate ! most-extensional horizontal principal value e1_rate = MIN(e1h, e2h, err) ! always negative e3_rate = MAX(e1h, e2h, err) ! always positive e2_rate = 0.0 - e1_rate - e3_rate ! may have either sign IF (e2_rate < 0.0) THEN ! both e1_rate and e2_rate are negative; positive e3_rate is the unique axis: scalar_strainrate = e3_rate ELSE ! e2_rate >= 0.0 ; so both e2_rate and e3_rate are positive; negative e1_rate is the unique axis: scalar_strainrate = (-e1_rate) END IF total_moment_rate_SI = total_moment_rate_SI + ( cell_area_m2(i) * effective_cz_meters(i, j) * mean_G(i, j) * 2.0 * scalar_strainrate ) END IF END DO END DO CLOSE (UNIT = 3, DISP = "KEEP") ! moment_factors_in_SI.dat file (for Corne Kreemer) WRITE (*, "(' Total seismic moment rate in SI = ', ES10.4,' N m /s.')") total_moment_rate_SI WRITE (*, "(' Job completed.')") WRITE (21, "('Total seismic moment rate in SI = ', ES10.4,' N m /s.')") total_moment_rate_SI WRITE (21, "('Job completed.')") CALL Pause() CLOSE (21) ! log file CONTAINS INTEGER FUNCTION Get_Support(array, i, j, nRows, nColumns) ! Returns number of neighbors (0~8) of cell (i, j) which have same value. ! The array is a global grid on a sphere, with wrap-around to E and W, but not to N and S. IMPLICIT NONE INTEGER*1, DIMENSION(:, :), INTENT(IN) :: array INTEGER, INTENT(IN) :: i, j, nRows, nColumns INTEGER :: count, jTarget, seeking seeking = array(i, j) count = 0 IF (i > 1) THEN !look in next row to N: jTarget = j-1; IF (jTarget < 1) jTarget = jTarget + nColumns IF (array(i-1, jTarget) == seeking) count = count + 1 IF (array(i-1, j) == seeking) count = count + 1 jTarget = j+1; IF (jTarget > nColumns) jTarget = jTarget - nColumns IF (array(i-1, jTarget) == seeking) count = count + 1 END IF IF (i < nColumns) THEN !look in next row to S: jTarget = j-1; IF (jTarget < 1) jTarget = jTarget + nColumns IF (array(i+1, jTarget) == seeking) count = count + 1 IF (array(i+1, j) == seeking) count = count + 1 jTarget = j+1; IF (jTarget > nColumns) jTarget = jTarget - nColumns IF (array(i+1, jTarget) == seeking) count = count + 1 END IF !look to West: jTarget = j-1; IF (jTarget < 1) jTarget = jTarget + nColumns IF (array(i, jTarget) == seeking) count = count + 1 !look to East: jTarget = j+1; IF (jTarget > nColumns) jTarget = jTarget - nColumns IF (array(i, jTarget) == seeking) count = count + 1 Get_Support = count END FUNCTION Get_Support REAL FUNCTION Moment(magnitude) ! Returns scalar seismic moment in N m, per ! equation (1) of Bird & Kagan [2004; Bull. Seism. Soc. Am.]; ! originally from Hanks & Kanamori [1979; J. Geophys. Res., v. 84, p. 2348-2350]. IMPLICIT NONE REAL, INTENT(IN) :: magnitude moment = 10.**((1.5 * magnitude) + 9.05) END FUNCTION Moment SUBROUTINE Pause () IMPLICIT NONE CHARACTER*1 c1 WRITE (*,"(' Press [Enter] to continue...'\)") READ (*,"(A)") c1 END SUBROUTINE Pause SUBROUTINE SHIFT_continuum_seismicity (desired_threshold_moment_Nm, & & e1h, e2h, err, T5_column, & & seismicity) !Computes "seismicity" (# of earthquakes/m^2/s, including aftershocks) !above "desired_threshold_moment_Nm", !based on principal long-term strain rates !(e1h = most-compressive horizontal; ! e2h = least-compressive horizontal; ! err = vertical) !and choice of most comparable plate boundary expressed by !selection of column "T5_column" in the T5_{whatever} tables in !the main program above. !The only reason that this code is isolated in a subprogram !is that it is called more than once, !and redundant code would be needed if this were placed in-line. !The underlying assumptions and equations are those of the SHIFT !(Seismic Hazard Inferred From Tectonics) !method of Bird & Liu [2007, Seismol. Res. Lett.]. IMPLICIT NONE !But note that arrays T5_{whatever} of main program are read: INTENT(IN). DOUBLE PRECISION, INTENT(IN) :: desired_threshold_moment_Nm REAL, INTENT(IN) :: e1h, e2h, err INTEGER, INTENT(IN) :: T5_column REAL, INTENT(OUT) :: seismicity REAL :: CMT_threshold_moment, corner_moment, e1_rate, e2_rate, e3_rate, & & M_persec_per_m2, M_persec_per_m3, & & seismicity_at_CMT_threshold, tGR_beta DOUBLE PRECISION :: argument, G_function !Convert principal strain rates in continuum to seismic moment rate per unit volume, !per equation (7b) of Bird & Liu [2007]): e1_rate = MIN(e1h, e2h, err) ! always negative e3_rate = MAX(e1h, e2h, err) ! always positive e2_rate = 0.0 - e1_rate - e3_rate ! may have either sign IF (e2_rate < 0.0) THEN ! both e1_rate and e2_rate are negative; positive e3_rate is the unique axis: M_persec_per_m3 = T5_assumed_mu_Pa(T5_column) * 2.0 * e3_rate ELSE ! e2_rate >= 0.0 ; so both e2_rate and e3_rate are positive; negative e1_rate is the unique axis: M_persec_per_m3 = T5_assumed_mu_Pa(T5_column) * 2.0 * (-e1_rate) END IF !Convert moment rate per unit volume to moment rate per unit area: M_persec_per_m2 = M_persec_per_m3 * T5_coupledThickness_m(T5_column) !Determine subsampling-point value of seismicity (earthquakes/m**2/s), !per equation (4) of Bird & Liu [2007]: seismicity_at_CMT_threshold = T5_CMT_pure_event_rate(T5_column) * 1.0D0 * & & M_persec_per_m2 / T5_tGR_model_moment_rate(T5_column) !Caution: Must avoid underflow, as terms are likely to ! be roughly: 2E-7 * 3E-6 / 4E12 = 1E-25. ! Therefore, 1.0D0 is introduced into product ! to force double-precision computation. !Correct seismicity to desired_threshold_moment_Nm, !per equation (5) of Bird & Liu [2007]: !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 above: tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) CMT_threshold_moment = T5_CMT_threshold_moment(T5_column) !G = ((desired_threshold_moment_Nm / CMT_threshold_moment)**(-beta)) * & ! & EXP((CMT_threshold_moment - desired_threshold_moment_Nm) / corner_moment) argument = (CMT_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / CMT_threshold_moment)**(-tGR_beta)) * & & EXP(argument) END IF seismicity = G_function * seismicity_at_CMT_threshold END SUBROUTINE SHIFT_continuum_seismicity END PROGRAM SHIFT_GSRM2x