PROGRAM SHIFT_GSRM2f_for_CSEP ! Applies the SHIFT (Seismic Hazard Inferred From Tectonics) ! hypotheses [Bird & Liu, 2007, Seismol. Res. Lett.] ! to kinematic model GSRM2.1 (Global Strain Rate Map, version 2.1) ! of Kreemer et al. [2014, Geochem. Geophys. Geosys.], ! to produce a long-term (or indefinite-term, or open-window) ! global forecast of shallow seismicity. ! ! Converting each 2-D strain-rate tensor to seimicity involves ! deciding which of the types of plate boundary studied by ! Bird & Kagan [2004, Bull. Seismol. Soc. Am.] is "most comparable." ! This is done with new logic, based on the global map of "tectonic zones" ! presented by Kagan et al. [2010, PAGEOPH]. ! (Tables 1 & 2 of Bird & Liu [2007] cannot be used, ! as they implicitly assumed that all subduction zones and spreading ridges ! and oceanic transforms would be represented by discrete faults.) ! ! An important approximation necessary to the conversion is the assumption ! that the (largely elastic) strain-rates measured by GPS geodesy and represented ! by GSRM2.1 are not very different in map pattern from long-term permanent ! strain rates. We consider the approximation adequate at the scale of this model. ! ! Program SHIFT_GSRM2x was written by Peter Bird, UCLA, in collaboration with ! Corne Kreemer. It was finalized in May 2014. ! This was the basis for the publication by Bird & Kreemer [2015, BSSA]. ! ! This program SHIFT_GSRM2f_for_CSEP applies the same logic to the same data, ! but produces output in the format of the Global Testing Region experiment ! at the Collaboratory for the Study of Earthquake Predictability (CSEP). ! Specifically, it increases nominal resolution to 0.1 x 0.1 degrees, ! and divides the seismicity of each cell among a number of magnitude bins, ! and it uses the XML format developed at CSEP. ! This version was written 2014.05.08 by Peter Bird at UCLA. !========================= 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 !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. !----------------------------------------------------------------------------------------- !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*6 :: lat_c6 CHARACTER*7 :: lon_c7 CHARACTER*12 :: duration_days_c12 CHARACTER*21 :: appended_data CHARACTER*132 :: first_line, last_line, line CHARACTER*160 :: averageStrainDatPath, calibrationCatalogDatPath, forecastDurationYears, & & forecastEndDate, forecastStartDate, forecastXmlPath, minTargetMagnitude, & & parameter_pathfile, plateBoundariesDatPath, tectonicZonesDatPath, text_line, & & underwaterDatPath INTEGER, PARAMETER :: text_line_length = 160 INTEGER :: best_azimuth, best_code, bins, & & boundary_calibration_count, & & code_now, CSEP_grid_columns, CSEP_grid_rows, & & dCol_3_sigma, dRow_3_sigma, & & eq_year, eq_depth_int, & & i, iCSEP, ii, iN, internal_points, intraplate_calibration_count, & & ios, iS, isolation_count, it, iTarget, iteration, iTest, iZone, & & j, jBeside, jCSEP, jE, jj, jTarget, jTargetE, jTargetW, jW, & & k, & & line_count, & & m, magnitude_bins, & & n, 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 :: pausing = .FALSE. ! Controls execution of all IF (pausing) CALL Pause() statements that allow user to read messages. 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, COS_latitude_radians, & & CSEP_grid_d_lat, CSEP_grid_d_lon, CSEP_grid_lat_max, CSEP_grid_lat_min, CSEP_grid_lon_max, CSEP_grid_lon_min, & & distance_m, duration_in_REAL_years, dx, dy, & & E_Sup, eq_Elon, eq_mag, eq_Nlat, & & e1_rate, e1h, e2_rate, e2h, e3_rate, err, EW_fraction, exx, exy, eyy, & & forecast_threshold_magnitude, fraction, & & G_function, getLon, GSRM2_Intraplate_corner_Nm, & & integral_converted, & & km_long, & & lat, lat1, lat2, latitude, lon, lon1, lon2, longitude, & & 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, single_bin_center_magnitude, & & sphere_area_m2, standard, strip_area_steradian, & & test_distance_m, test_radians2, tester, & & this_Sup, three_sigma_m, threshold_magnitude, top_lat, total_cell_seismicity, 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 :: magnitude_bin_threshold 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 :: argument, & & catalog_threshold_moment_NM, center_normal_rate, & & dilatation_rate, & & epicenters_per_m2_per_s_inMagBin, & & integral, intraplate_area_m2, intraplate_corner_Nm, intraplate_seismicity_GCMT, & & plate_boundary_area_m2, predicted_EQ_count, & & radius_rate, & & second_invariant, & & tentative_boundary_eqs_per_catalog, tentative_boundary_eqs_per_s, & & total_area_m2 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: intraplate_seismicity_SI DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: moment_bin_threshold 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 DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: epicenter_rate_density !**************************************************************************************** !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_GSRM2f_for_CSEP:')") WRITE (*, *) WRITE (*, "(' Applies the SHIFT (Seismic Hazard Inferred From Tectonics)')") WRITE (*, "(' hypotheses [Bird & Liu, 2007, Seismol. Res. Lett.]')") WRITE (*, "(' to kinematic model GSRM2.1 (Global Strain Rate Map, version 2.1)')") WRITE (*, "(' of Kreemer et al. [2014, Geochem. Geophys. Geosys.],')") WRITE (*, "(' to produce a long-term (or indefinite-term, or open-window)')") WRITE (*, "(' global forecast of shallow seismicity.')") WRITE (*, *) WRITE (*, "(' Converting each 2-D strain-rate tensor to seimicity involves')") WRITE (*, "(' deciding which of the types of plate boundary studied by')") WRITE (*, "(' Bird & Kagan [2004, Bull. Seismol. Soc. Am.] is ""most comparable.""')") WRITE (*, "(' This is done with new logic, based on the global map of ""tectonic zones""')") WRITE (*, "(' presented by Kagan et al. [2010, PAGEOPH].')") WRITE (*, "(' (Tables 1 & 2 of Bird & Liu [2007] cannot be used,')") WRITE (*, "(' as they implicitly assumed that all subduction zones and spreading ridges')") WRITE (*, "(' and oceanic transforms would be represented by discrete faults.)')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' An important approximation necessary to the conversion is the assumption')") WRITE (*, "(' that the (largely elastic) strain-rates measured by GPS geodesy and represented')") WRITE (*, "(' by GSRM2.1 are not very different in map pattern from long-term permanent')") WRITE (*, "(' strain rates. We consider the approximation adequate at the scale of this model.')") WRITE (*, *) WRITE (*, "(' Program SHIFT_GSRM2x was written by Peter Bird, UCLA, in collaboration with')") WRITE (*, "(' Corne Kreemer. It was finalized in May 2014.')") WRITE (*, "(' This was the basis for the publication by Bird & Kreemer [2015, BSSA].')") WRITE (*, *) WRITE (*, "(' This program SHIFT_GSRM2f_for_CSEP applies the same logic to the same data,')") WRITE (*, "(' but produces output in the format of the Global Testing Region experiment')") WRITE (*, "(' at the Collaboratory for the Study of Earthquake Predictability (CSEP).')") WRITE (*, "(' Specifically, it increases nominal resolution to 0.1 x 0.1 degrees,')") WRITE (*, "(' and divides the seismicity of each cell among a number of magnitude bins,')") WRITE (*, "(' and it uses the XML format developed at CSEP.')") WRITE (*, "(' This version was written 2014.05.08 by Peter Bird at UCLA.')") WRITE (*, *) OPEN (UNIT = 21, FILE = "SHIFT_GSRM2f_for_CSEP_log.txt") ! absolute OPEN; overwrites any existing WRITE (21, "('PROGRAM SHIFT_GSRM2f_for_CSEP:')") WRITE (21, *) WRITE (21, "('Applies the SHIFT (Seismic Hazard Inferred From Tectonics)')") WRITE (21, "('hypotheses [Bird & Liu, 2007, Seismol. Res. Lett.]')") WRITE (21, "('to kinematic model GSRM2.1 (Global Strain Rate Map, version 2.1)')") WRITE (21, "('of Kreemer et al. [2014, Geochem. Geophys. Geosys.],')") WRITE (21, "('to produce a long-term (or indefinite-term, or open-window)')") WRITE (21, "('global forecast of shallow seismicity.')") WRITE (21, *) WRITE (21, "('Converting each 2-D strain-rate tensor to seimicity involves')") WRITE (21, "('deciding which of the types of plate boundary studied by')") WRITE (21, "('Bird & Kagan [2004, Bull. Seismol. Soc. Am.] is ""most comparable.""')") WRITE (21, "('This is done with new logic, based on the global map of ""tectonic zones""')") WRITE (21, "('presented by Kagan et al. [2010, PAGEOPH].')") WRITE (21, "('(Tables 1 & 2 of Bird & Liu [2007] cannot be used,')") WRITE (21, "('as they implicitly assumed that all subduction zones and spreading ridges')") WRITE (21, "('and oceanic transforms would be represented by discrete faults.)')") WRITE (21, *) WRITE (21, "('An important approximation necessary to the conversion is the assumption')") WRITE (21, "('that the (largely elastic) strain-rates measured by GPS geodesy and represented')") WRITE (21, "('by GSRM2.1 are not very different in map pattern from long-term permanent')") WRITE (21, "('strain rates. We consider the approximation adequate at the scale of this model.')") WRITE (21, *) WRITE (21, "('Program SHIFT_GSRM2x was written by Peter Bird, UCLA, in collaboration with')") WRITE (21, "('Corne Kreemer. It was finalized in May 2014.')") WRITE (21, "('This was the basis for the publication by Bird & Kreemer [2015, BSSA].')") WRITE (21, *) WRITE (21, "('This program SHIFT_GSRM2f_for_CSEP applies the same logic to the same data,')") WRITE (21, "('but produces output in the format of the Global Testing Region experiment')") WRITE (21, "('at the Collaboratory for the Study of Earthquake Predictability (CSEP).')") WRITE (21, "('Specifically, it increases nominal resolution to 0.1 x 0.1 degrees,')") WRITE (21, "('and divides the seismicity of each cell among a number of magnitude bins,')") WRITE (21, "('and it uses the XML format developed at CSEP.')") WRITE (21, "('This version was written 2014.05.08 by Peter Bird at UCLA.')") WRITE (21, *) !------------------------------------------------------------------------------------- !Find, read, and memorize contents of the parameter file: ! ! SHIFT_GSRM2f_parameters.dat: ! ---------------------------------------------------- ! plateBoundariesDatPath=PB2002_steps.dat ! tectonicZonesDatPath=PB2002_tectonic_zones.grd ! averageStrainDatPath=GSRM_average_strain_v2.1.txt ! calibrationCatalogDatPath=GCMT_shallow_5p767_1977-2013.eqc ! forecastStartDate=2014-01-01T00:00:00Z ! forecastEndDate=2015-01-01T00:00:00Z ! forecastDurationYears=1.0 ! minTargetMagnitude=5.95 ! forecastXmlPath=SHIFT_GSRM2f.global.forecast.xml ! ---------------------------------------------------- !(1) Find out if there is a command-line argument which might be the name of the parameter file: n = COMMAND_ARGUMENT_COUNT() ! standard in Fortran 2003+ !n = 1 ! kludge for Fortran 95 testing !(2) Test the value of n: IF (n < 1) THEN WRITE (*, "(' ERROR: [/path/]filename of parameter file must follow in command line!')") WRITE (*, "(' SHIFT_GSRM2f_parameters.dat was supplied to CSEP.')") WRITE (21, "('ERROR: [/path/]filename of parameter file must follow in command line!')") WRITE (21, "('SHIFT_GSRM2f_parameters.dat was supplied to CSEP.')") IF (pausing) CALL Pause() STOP END IF !(3) Get the text of the first command-line argument: CALL GET_COMMAND_ARGUMENT (1, parameter_pathfile) ! standard in Fortran 2003 !parameter_pathfile = "SHIFT_GSRM2f_parameters.dat" ! kludge for Fortran 95 testing !(4) Try to open this file, and test for success: OPEN (UNIT = 11, FILE = TRIM(parameter_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: [/path/]filename of parameter file in command line is invalid!')") WRITE (21, "('ERROR: [/path/]filename of parameter file in command line is invalid!')") IF (pausing) CALL Pause() STOP END IF !(5) Read and remember parameter values: parameter_scanning: DO READ (11, "(A)", IOSTAT = ios) text_line IF (ios < 0) EXIT parameter_scanning IF (text_line(1:18) == "underwaterDatPath=") THEN underwaterDatPath = TRIM(text_line(19:text_line_length)) ELSE IF (text_line(1:23) == "plateBoundariesDatPath=") THEN plateBoundariesDatPath = TRIM(text_line(24:text_line_length)) ELSE IF (text_line(1:21) == "tectonicZonesDatPath=") THEN tectonicZonesDatPath = TRIM(text_line(22:text_line_length)) ELSE IF (text_line(1:21) == "averageStrainDatPath=") THEN averageStrainDatPath = TRIM(text_line(22:text_line_length)) ELSE IF (text_line(1:26) == "calibrationCatalogDatPath=") THEN calibrationCatalogDatPath = TRIM(text_line(27:text_line_length)) ELSE IF (text_line(1:18) == "forecastStartDate=") THEN forecastStartDate = TRIM(text_line(19:text_line_length)) ELSE IF (text_line(1:16) == "forecastEndDate=") THEN forecastEndDate = TRIM(text_line(17:text_line_length)) ELSE IF (text_line(1:22) == "forecastDurationYears=") THEN forecastDurationYears = TRIM(text_line(23:text_line_length)) READ (forecastDurationYears, * , IOSTAT = ios) duration_in_REAL_years ELSE IF (text_line(1:19) == "minTargetMagnitude=") THEN minTargetMagnitude = TRIM(text_line(20:text_line_length)) READ (minTargetMagnitude, *, IOSTAT = ios) threshold_magnitude single_bin_center_magnitude = threshold_magnitude + 0.05 ! assuming bins 0.1 magnitude units wide ELSE IF (text_line(1:16) == "forecastXmlPath=") THEN forecastXmlPath = TRIM(text_line(17:text_line_length)) END IF END DO parameter_scanning CLOSE (UNIT = 11, DISP = "KEEP") !------------------------------------------------------------------------------------- !Note that internal seismicity grid is GSRM2-sized (0.25 x 0.20 degree cells), !even though eventual output resolution will be nominally 0.1 x 0.1 degrees. nColumns = NINT(360.0 / GSRM2_d_lon) nRows = 1 + NINT((GSRM2_lat_max - GSRM2_lat_min) / GSRM2_d_lat) nGridPoints = nRows * nColumns ALLOCATE ( GSRM2_grid(nRows, nColumns, 3 ) ) GSRM2_grid = 0.0 !Describe CSEP's magnitude-bin system, and provide corresponding storage dimension: magnitude_bins = 31 ! widths of 0.10, and centers at 6.00, 6.10, ..., 8.90, 9.00. Last bin to be open on high-side. ALLOCATE ( magnitude_bin_threshold(magnitude_bins) ) ALLOCATE ( moment_bin_threshold(magnitude_bins) ) ALLOCATE ( intraplate_seismicity_SI(magnitude_bins) ) DO m = 1, magnitude_bins magnitude_bin_threshold(m) = 5.95 + (m-1) * 0.10 moment_bin_threshold(m) = Moment(magnitude_bin_threshold(m)) END DO ALLOCATE ( epicenter_rate_density(magnitude_bins, nRows, nColumns) ) epicenter_rate_density = 0.0D0 ! just to aid debugging; will be replaced with positive value WRITE (*, *) WRITE (*, "(' Note that internal calculation occurs on a')") WRITE (*, "(' GSRM2.1-type grid of 0.25 x 0.20 degree cells.')") WRITE (*, "(' Creating global grids of ',I6,' rows * ',I6,' columns = ',I12,' points')") nRows, nColumns, nGridPoints WRITE (*, "(' for input strain rate tensors, epicentral rate densities, & effective cz.')") WRITE (21, *) WRITE (21, "('Note that internal calculation occurs on a')") WRITE (21, "(' GSRM2.1-type grid of 0.25 x 0.20 degree cells.')") WRITE (21, "('Creating global grids of ',I6,' rows * ',I6,' columns = ',I12,' points')") nRows, nColumns, nGridPoints WRITE (21, "('for input strain rate tensors, 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 WRITE (*, *) IF (pausing) CALL Pause() !Read in Corne Kreemer's GSRM2 strain rates: WRITE (*, *) WRITE (*, "(' FILE ',A,' must be available to this program.')") TRIM(averageStrainDatPath) WRITE (21, *) WRITE (21, "('FILE ',A,' must be available to this program.')") TRIM(averageStrainDatPath) OPEN (UNIT = 1, FILE = TRIM(averageStrainDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ',A,' not found. Fix, and try again...')") TRIM(averageStrainDatPath) WRITE (21, "('ERROR: ',A,' not found. Fix, and try again...')") TRIM(averageStrainDatPath) IF (pausing) 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, !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 (*, *) WRITE (*, "(' FILE ',A,' (seismic catalog for calibration)')") TRIM(calibrationCatalogDatPath) WRITE (*, "(' is required to be available.')") WRITE (21, *) WRITE (21, "('FILE ',A,' (seismic catalog for calibration)')") TRIM(calibrationCatalogDatPath) WRITE (21, "(' is required to be available.')") OPEN (UNIT = 2, FILE = TRIM(calibrationCatalogDatPath), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ',A,' not found. Fix, and try again...')") TRIM(calibrationCatalogDatPath) IF (pausing) CALL Pause() STOP END IF 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 catalog_years = 37.0 ! <=== user choice removed; THIS LINE IS SPECIFIC TO THE CALIBRATION CATALOG YEARS 1977-2013 !!! WRITE (*, *) 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 catalog_threshold_magnitude = 5.767 ! <=== user choice removed WRITE (*, *) 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 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 !-------------------------------------------------------------------------------- ! Strain-rates must be smoothed in the offshore regions: !locate, read, and memorize underwater.grd: WRITE (21, *) WRITE (21, "('FILE ',A,' must be available. Please check this.')") TRIM(underwaterDatPath) 200 WRITE (*, *) WRITE (*, "(' FILE ',A,' must be available. Please check this.')") TRIM(underwaterDatPath) IF (pausing) CALL Pause() OPEN (UNIT = 1, FILE = TRIM(underwaterDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: FILE ', A,' must be available. Please check this.')") TRIM(underwaterDatPath) IF (pausing) 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 ',A,' must be available. Please check this.')") TRIM(plateBoundariesDatPath) WRITE (21, "('FILE ',A,' must also be available. Please check this.')") TRIM(tectonicZonesDatPath) 210 WRITE (*, *) WRITE (*, "(' FILE ',A,' must be available. Please check this.')") TRIM(plateBoundariesDatPath) WRITE (*, "(' FILE ',A,' must also be available. Please check this.')") TRIM(tectonicZonesDatPath) IF (pausing) CALL Pause() OPEN (UNIT = 1, FILE = TRIM(plateBoundariesDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: FILE ',A,' was not found in active folder.')") TRIM(plateBoundariesDatPath) WRITE (21, "('ERROR: FILE ',A,' was not found in active folder.')") TRIM(plateBoundariesDatPath) IF (pausing) 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 IF (pausing) 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.')") IF (pausing) 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 (*, "(' ',A,' has been converted to a list of ',I6,' internal points.')") TRIM(plateBoundariesDatPath), PB2002_list_highest WRITE (21, "( A,' has been converted to a list of ',I6,' internal points.')") TRIM(plateBoundariesDatPath), 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 ({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. !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) IF (pausing) 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.')") ! 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)") IF (pausing) 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 IF (pausing) 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 IF (pausing) 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 IF (pausing) 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 IF (pausing) 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 of directed smoothing of SUBs) ! Use tectonic zones [5 DOF] ------------------------------------------- !Get tectonic zones grid from file of Bird et al. [2010, Pure Appl. Geoph.]" OPEN (UNIT = 1, FILE = TRIM(tectonicZonesDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ',A,' was not found in active folder.')") TRIM(tectonicZonesDatPath) WRITE (21, "('ERROR: ',A,' was not found in active folder.')") TRIM(tectonicZonesDatPath) IF (pausing) 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 = TRIM(tectonicZonesDatPath), 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 !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.')") !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(calibrationCatalogDatPath), STATUS = "OLD", PAD = "YES", IOSTAT = ios) !new_EQC_file = "SHIFT_GSRM2f_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 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 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.')") IF (pausing) 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_GCMT = intraplate_calibration_count / & & (intraplate_area_m2 * catalog_years * s_per_year) WRITE (*, "(' Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_GCMT WRITE (*, "(' at stated threshold magnitude of the calibration catalog.')") WRITE (21, "('Intraplate seismicity rate (per m^2, per s) = ',ES12.5)") intraplate_seismicity_GCMT WRITE (21, "(' at stated threshold magnitude of the calibration catalog.')") !------------------------------------------------------------------------------------- !Precompute intraplate_seismicity (to save duplication, and to use a lower-limit): !Use intraplate corner magnitude of 9.0 and beta = 0.63 inferred by Bird et al. [2010, SRL]. intraplate_corner_Nm = Moment(9.0) DO m = 1, magnitude_bins argument = (catalog_threshold_moment_Nm - moment_bin_threshold(m)) / intraplate_corner_Nm IF (argument < -100.0D0) THEN G_function = 0.0D0 WRITE (*, "(' WARNING: argument underflow in EXP(argument): intraplate seismicity will be zero.')") WRITE (21, "('WARNING: argument underflow in EXP(argument): intraplate seismicity will be zero.')") IF (pausing) CALL Pause() ELSE G_function = ((moment_bin_threshold(m) / catalog_threshold_moment_Nm)**(-0.63)) * & & EXP(argument) END IF intraplate_seismicity_SI(m) = G_function * intraplate_seismicity_GCMT END DO !N.B. This uniform value will be applied to all iZone == 0 grid points, and will be used as a lower limit ! for the seismicity of other grid points as well. !------------------------------------------------------------------------------------- !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 velocity_factor(2) = 1.0 ! and remains 1.0, except for SUB & CCB ! 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 ! 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 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 ! 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 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 ! 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 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 of 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_GCMT * 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: DO i = 1, nRows DO j = 1, nColumns DO m = 1, magnitude_bins epicenter_rate_density(m, i, j) = intraplate_seismicity_SI(m) ! whole array; only iZone>0 points will be revisited: END DO END DO END DO WRITE (*, *) WRITE (*, "(' Computing gridded seismicity for each of ',I2,' magnitude thresholds:')") magnitude_bins WRITE (21, *) WRITE (21, "('Computing gridded seismicity for each of ',I2,' magnitude thresholds:')") magnitude_bins DO m = 1, magnitude_bins 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 velocity_factor(2) = 1.0 ! and remains 1.0, except for SUB & CCB ! 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 ! 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 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 ! 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 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 ! 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 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 of classification and partitioning) !^^^^^^ pass 2 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ total_cell_seismicity = 0.0 ! initialize, before sum of 1 or 2 terms DO k = 1, nTensorParts CALL SHIFT_continuum_seismicity (moment_bin_threshold(m), & ! 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 epicenter_rate_density(m, i, j) = boundary_scaling_factor_in_zone(iZone) * total_cell_seismicity + & & intraplate_seismicity_SI(m) END IF ! iZone > 0 END DO ! loop j = 1, nColumns END DO ! loop i = 1, nRows WRITE (*, "(' Completed magnitude-bin ',I2,' with lower threshold of ',F5.2)") m, magnitude_bin_threshold(m) WRITE (21, "('Completed magnitude-bin ',I2,' with lower threshold of ',F5.2)") m, magnitude_bin_threshold(m) END DO ! loop on m = 1, magnitude_bins !----------------------------------------------------------------------------------- ! Convert epicenter_rate_density from all-above-threshold to incremental only-those-within-bin format ! (except for last bin, m = magnitude_bins, which stays open): WRITE (*, *) WRITE (*, "(' Converting seismicity grid from rates-above-threshold to rates-in-magnitude-bin.')") WRITE (21, *) WRITE (21, "('Converting seismicity grid from rates-above-threshold to rates-in-magnitude-bin.')") DO m = 1, (magnitude_bins-1) DO i = 1, nRows DO j = 1, nColumns epicenter_rate_density(m, i, j) = epicenter_rate_density(m, i, j) - epicenter_rate_density(m+1, i, j) END DO END DO intraplate_seismicity_SI(m) = intraplate_seismicity_SI(m) - intraplate_seismicity_SI(m+1) END DO !----------------------------------------------------------------------------------- !Integrate seismicity over (internal, GSRM2-shaped) grid area: WRITE (*, *) WRITE (*, "(' Computing area-integral of seismicity...')") WRITE (21, *) WRITE (21, "('Computing area-integral of seismicity...')") integral = 0.0D0 ! (just initializing) DO i = 1, nRows DO j = 1, nColumns DO m = 1, magnitude_bins integral = integral + epicenter_rate_density(m, i, j) * cell_area_m2(i) END DO 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 (*, "(' above the CSEP lowest-magnitude-bin lower-threshold of 5.95.')") WRITE (21, "(/'Area integral of seismicity =' & & /ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century')") integral, integral_converted WRITE (21, "(' above the CSEP lowest-magnitude-bin lower-threshold of 5.95.')") !Provide global-mean frequency/magnitude curve, in terms of epicentral rate densities: WRITE (*, *) WRITE (*, "(' Computing global-mean frequency-density/magnitude curve, in SI...')") WRITE (21, *) WRITE (21, "('Computing global-mean frequency-density/magnitude curve, in SI...')") DO m = 1, magnitude_bins ! where m identifies the 1st bin used in the sum integral = 0.0D0 ! (just initializing) DO i = 1, nRows DO j = 1, nColumns DO n = m, magnitude_bins integral = integral + epicenter_rate_density(n, i, j) * cell_area_m2(i) END DO END DO END DO integral_converted = integral / sphere_area_m2 WRITE (*, "(' Above m = ',F5.2,' global-mean seismicity == ',ES12.4,' EQs/m^2/s.')") magnitude_bin_threshold(m), integral_converted WRITE (21, "('Above m = ',F5.2,' global-mean seismicity == ',ES12.4,' EQs/m^2/s.')") magnitude_bin_threshold(m), integral_converted END DO !----------------------------------------------------------------------------------- !Output the forecast .xml: WRITE (*, *) WRITE (*, "(' Writing .xml file...')") WRITE (21, *) WRITE (21, "('Writing .xml file...')") OPEN (UNIT = 3, FILE = TRIM(forecastXmlPath)) ! unconditional OPEN; overwrites any existing file WRITE (3, '("")') WRITE (3, '("")') WRITE (3, '(" ")') WRITE (3, '(" SHIFT_GSRM2f")') WRITE (3, '(" 1.0")') ! <=== v1.0=2014.03.21 WRITE (3, '(" Bird and Kreemer [2015, BSSA]")') WRITE (3, '(" smi://org.scec/persons/Bird_Peter")') WRITE (3, '(" 2014-05-08T23:59:59Z")') WRITE (3, '(" ",A,"")') TRIM(forecastStartDate) WRITE (3, '(" ",A,"")') TRIM(forecastEndDate) WRITE (duration_days_c12, '(F12.2)') 365.25*duration_in_REAL_years duration_days_c12 = ADJUSTL(duration_days_c12) WRITE (3, '(" ",A,"")') TRIM(duration_days_c12) WRITE (3, '(" ")') WRITE (3, '(" 0.1")') WRITE (3, '(" 5.95")') WRITE (3, '(" 1")') WRITE (3, '(" ")') !Define the grid of cell centers for which seismicity is to be computed and written out: CSEP_grid_lon_min = -179.95 ! degrees East CSEP_grid_lon_max = +179.95 ! degrees East CSEP_grid_lat_min = -89.95 ! degrees North CSEP_grid_lat_max = 89.95 ! degrees North !Note that this grid refers to the centers of the CSEP cells, not their outer edges. CSEP_grid_d_lon = 0.1 ! degrees. CSEP_grid_d_lat = 0.1 ! degrees CSEP_grid_rows = 1 + NINT((CSEP_grid_lat_max - CSEP_grid_lat_min) / CSEP_grid_d_lat) CSEP_grid_columns = 1 + NINT((CSEP_grid_lon_max - CSEP_grid_lon_min) / CSEP_grid_d_lon) !GPBhere DO jCSEP = 1, CSEP_grid_columns ! CSEP wants W to E, -179.95 to +179.95 longitude = CSEP_grid_lon_min + (jCSEP - 1) * CSEP_grid_d_lon !----------------------------------------------------------------------------------------- !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. !Make sure that longitude stays in range -179.999... to +179.999.... that matches GSRM2 and epicentral_rate_density grids: IF (longitude > +180.0) longitude = longitude - 360.0 IF (longitude < -180.0) longitude = longitude + 360.0 DO iCSEP = CSEP_grid_rows, 1, -1 ! (i = 1, grid_rows would be N to S, but CSEP wants S to N) latitude = CSEP_grid_lat_max - (iCSEP - 1) * CSEP_grid_d_lat WRITE (lat_c6, "(F6.2)") latitude lat_c6 = ADJUSTL(lat_c6) WRITE (lon_c7, "(F7.2)") longitude lon_c7 = ADJUSTL(lon_c7) WRITE (3, '(" ")') TRIM(lon_c7), TRIM(lat_c6) !spatial area of CSEP cell, in m**2: COS_latitude_radians = COS(latitude * radians_per_degree) dx = CSEP_grid_d_lon * COS_latitude_radians * radians_per_degree * Earth_radius_m dy = CSEP_grid_d_lat * radians_per_degree * Earth_radius_m !determine which cell (iTarget, jTarget) of internal grid will have the required epicentral_rate_density? iTarget = 1 + NINT((GSRM2_lat_max - latitude) / GSRM2_d_lat) ! caution: will be out-of-range near the poles iTarget = MAX(iTarget, 1) ! so values at lat = +89.95 are taken from the top-row cell iTarget = MIN(iTarget, nRows) ! so values at lat = -89.95 are taken from the last-row cells DO m = 1, magnitude_bins IF ((1 + MOD((jCSEP - 1), 5)) == 3) THEN ! CSEP cell is centered on the joint between two GSRM2 cells, to W and E jTargetW = 1 + NINT((longitude - (0.25 * CSEP_grid_d_lon) - GSRM2_lon_min) / GSRM2_d_lon) jTargetE = jTargetW + 1 epicenters_per_m2_per_s_inMagBin = 0.5 * epicenter_rate_density(m, iTarget, jTargetW) + & & 0.5 * epicenter_rate_density(m, iTarget, jTargetE) ELSE ! CSEP cell is entirely within one GSRM2 cell jTarget = 1 + NINT((longitude - GSRM2_lon_min) / GSRM2_d_lon) epicenters_per_m2_per_s_inMagBin = epicenter_rate_density(m, iTarget, jTarget) END IF ! hard or easy case predicted_EQ_count = epicenters_per_m2_per_s_inMagBin * dx * dy * duration_in_REAL_years * s_per_year WRITE (3, '(" ",ES10.3)') predicted_EQ_count END DO WRITE (3, '(" ")') END DO ! j = 1, grid_columns (W to E) END DO ! i = 1, grid_rows (N to S) WRITE (3, '(" ")') WRITE (3, '(" ")') WRITE (3, '("")') CLOSE (UNIT = 3, DISP = "KEEP") ! SHIFT_GSRM2f.global.forecast.xml file !----------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' Job completed.')") WRITE (21, *) WRITE (21, "('Job completed.')") IF (pausing) 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_GSRM2f_for_CSEP