PROGRAM GEAR1_for_CSEP ! Computes the Global Earthquake Activity Rate (GEAR) model, version 1, ! of long-term rates of shallow earthquakes on a global grid, as functions ! of longitude, latitude, and magnitude, as described by ! P. Bird, D.D. Jackson, Y.Y. Kagan, C. Kreemer, and R.S. Stein ! [2015, Seismol. Soc. Am. Bull.] ! This model is an optimal hybrid of the smoothed-Seismicity forecast ! of Kagan and Jackson [1994, 2010, 2011] with the Tectonics forecast of ! Bird and Kreemer (SHIFT_GSRM2f, Bird and Kreemer [2015, BSSA]). ! Specifically, it is the preferred hybrid H* described in our paper, ! with input datasets updated to end-2013 (or perhaps later?), ! and with an improved method for extrapolating the smoothed-seismicity ! parent forecast to high magnitudes. ! ! The primary computational steps are (in order): ! (1) Read the Kagan and Jackson smoothed-Seismicity long-term forecast ! grid and then scale it to 31 different threshold magnitudes; ! (2) Recreate the Tectonics forecast (SHIFT_GSRM2f) of Bird and Kreemer, ! exactly as in program SHIFT_GSRM2f_for_CSEP.f90, ! with forecast spatial grids for each of 31 threshold magnitudes; ! (3) Merge these 2 parent forecasts into hybrid forecast H*, ! with spatial grids for each of 31 threshold magnitudes; ! (4) Difference these grids to get seismicity in distinct magnitude bins. ! (5) Output the forecast in XML format. ! Run time is about ~5 hours if 16-fold parallel computing is available ! during step (1); otherwise, a few days may be needed. ! It is expected that 64-bit computing will be used because of the ! large memory requirements (4 GB) of this code. ! ! Note that random-number generators used in step (2) can come from ! either the Intel Math Kernel Library (MKL), or from the ! International Mathematics Subroutine Library (IMSL). ! A change-over would be effected by switching which lines of ! source code are commented-out, versus active. ! Search on "MKL" and also on "IMSL" to find these lines. ! The method of initialization of these generators is not important. ! ! This program produces output in the format of the Global Testing Region ! at the Collaboratory for the Study of Earthquake Predictability (CSEP). ! Specifically, it yields resolution of 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. ! The start-date and end-date of the forecast are specified by the user ! in the GEAR1_parameter.dat file (which is read from the command line); ! however, conceptually the forecast is open-ended toward the future. ! ! This source code (with 7 input datasets) is the primary means of ! distributing the forecast, for several reasons: ! (a) CSEP protocol requires submission of source code; ! (b) We support maximum transparency and open-source ethics; and ! (c) The output file is much larger (3.7 GB) than the sum of the ! source code and input files (about 470 MB, or 0.46 GB). ! ! This version was written 2014.06.04+ by Peter Bird at UCLA, ! using large chunks of code from the previous programs ! S_method2E_demo.f90, SHIFT_GSRM2f_for_CSEP.f90, ! and hybrid_forecasts_v3.f90. !========================= 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 SHIFT_GSRM2x model versions programmed to date; ! see user-prompt menu below at "SELECT DESIRED MODEL" ! for the meaning of each version. Here we use model "f". !---Parameters needed for statistical package Intel MKL/VSL: INTEGER :: status TYPE(VSL_STREAM_STATE) :: my_stream !----------------------------------------------------------------------------------------- !---Critical parameter from S_method2E_demo.f90, used in Step #1: ------------------------ REAL, PARAMETER :: smoothing_distance_m = 200000.0 !----------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------- !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 T_parent_coarse_3D 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, smoothing_distance_km_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, smoothedSeismicityDatPath, & & tectonicZonesDatPath, text_line, & & underwaterDatPath CHARACTER*79, DIMENSION(16) :: lines 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, & & fine_nColumns, fine_nRows, & ! <=== Kagan and Jackson smoothed-seismicity grid; also, GEAR1 output grid! & 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, low_argument_count, & & m, magnitude_bins, & & n, n_TINY, 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, & & smoothing_distance_km_INT, 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 the user some time 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, beta, biggest_eRate, biggest_n_per_year, & & bottom_lat, boundary_scaling_factor, & & catalog_threshold_magnitude, catalog_years, coefficient, corner_magnitude, 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, & & d, distance_m, duration_in_REAL_years, & & E_Sup, eq_Elon, eq_mag, eq_Nlat, & & e1_rate, e1h, e2_rate, e2h, e3_rate, err, EW_fraction, exx, exy, eyy, & & F_threshold_magnitude, F_threshold_moment, & ! Step #1 & fine_d_lat, fine_d_lon, fine_lat_max, fine_lat_min, fine_lon_max, fine_lon_min, & ! <== Kagan and Jackson smoothed-Seismicity; also GEAR1 output grid! & forecast_threshold_magnitude, fraction, & & getLon, GSRM2_Intraplate_corner_Nm, & & highest_beta0, highest_beta1, highest_beta2, highest_mc0, highest_mc1, highest_mc2, & ! Step #1 & integral_converted, & & km_long, & & lat, lat1, lat2, latitude, lon, lon1, lon2, longitude, & & lowest_beta0, lowest_beta1, lowest_beta2, lowest_mc0, lowest_mc1, lowest_mc2, & ! Step #1 & mean_beta0, mean_beta1, mean_beta2, mean_mc0, mean_mc1, mean_mc2, & ! Step #1 & 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, S_threshold_magnitude_Bird, S_threshold_magnitude_Kagan, & ! Step #1 & scalar_strainrate, seismicity, sigma_m, single_bin_center_magnitude, & & sphere_area_m2, standard, strip_area_steradian, & & standard_deviation_beta0, standard_deviation_beta1, standard_deviation_beta2, standard_deviation_mc0, & ! Step #1 & standard_deviation_mc1, standard_deviation_mc2, stretching_factor, & ! Step #1 & 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(0:4) :: beta_of_zone, corner_magnitude_of_zone ! Table 1 of Kagan et al. [2010]; used in Step #1 REAL, DIMENSION(:), ALLOCATABLE :: coarse_cell_area_m2, fine_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 REAL, DIMENSION(:, :), ALLOCATABLE :: beta_fine, smoothed_beta_fine, corner_magnitude_fine, smoothed_corner_magnitude_fine ! Step #1 (Seismicity) DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 DOUBLE PRECISION :: argument, & & catalog_threshold_moment_NM, center_normal_rate, coarse_grid_area_m2, correction, & & dilatation_rate, & & epicenters_per_m2_per_s_inMagBin, & & fine_grid_area_m2, floor_integral, floor_integral_per_year, & & G_function, & & H_floor, H_integral, H_integral_per_year, & & integral, intraplate_area_m2, intraplate_corner_Nm, intraplate_seismicity_GCMT, & & plate_boundary_area_m2, predicted_EQ_count, & & radius_rate, & & S_floor, S_integral, S_integral_per_year, S_lower_limit_at_m5p8, S_lower_limit_at_high_m, & & second_invariant, & & T_floor, T_integral, T_integral_per_year, target_N_per_year, tentative_boundary_eqs_per_catalog, & & tentative_boundary_eqs_per_s, total_area_m2, & & value 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 !HUGE ARRAYS (global grids in 2-D, which become 3-D when there are many magnitude thresholds: DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: S_parent_fine ! used by Step #1 (Seismicity input); then, re-used as temporary storage in Step #3 (hybridization) DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: S_parent_fine_3D ! created in Step #1 from contents of S_parent_fine; used in Step #3 DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: T_parent_coarse_3D ! created in Step #2 (Tectonics); used in Step #3 DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: T_parent_fine ! used temporarily in Step #3 (hybridization) to convert each T_parent map to 0.1 x 0.1 degree DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: H_star_fine ! used temporarily in Step #3 (hybridization) to create a map that will be saved in H_star_fine_3D DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: H_star_fine_3D ! created in Step #3 (hybridization); transformed in Step #4; output in Step #5 DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: floor_seismicity ! temporary array used in step #3. !----------------------------------------------------------------------------------------- !Definitions specific to parameters imported from Table 5 of Bird and 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 !****************** END SPECIFICATION SECTION; BEGIN EXECUTABLE CODE ************************************************* !GPBT5 ! Empirical coefficients from selected rows of Table 5 of Bird and 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 and 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 !----------------------------------------------------------------------------------------- R2 = Earth_radius_m**2 !----------------------------------------------------------------------------------------- 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' /) !----------------------------------------------------------------------------------------- !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 !----------------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' PROGRAM GEAR1_for_CSEP:')") WRITE (*, *) WRITE (*, "(' Computes the Global Earthquake Activity Rate (GEAR) model, version 1,')") WRITE (*, "(' of long-term rates of shallow earthquakes on a global grid, as functions')") WRITE (*, "(' of longitude, latitude, and magnitude, as described by')") WRITE (*, "(' P. Bird, D.D. Jackson, Y.Y. Kagan, C. Kreemer, and R.S. Stein')") WRITE (*, "(' [2015, Seismol. Soc. Am. Bull.]')") WRITE (*, "(' This model is an optimal hybrid of the smoothed-Seismicity forecast')") WRITE (*, "(' of Kagan and Jackson [1994, 2010, 2011] with the Tectonics forecast of ')") WRITE (*, "(' Bird and Kreemer (SHIFT_GSRM2f, Bird and Kreemer [2015, BSSA]).')") WRITE (*, "(' Specifically, it is the preferred hybrid H* described in our paper,')") WRITE (*, "(' with input datasets updated to end-2013 (or perhaps later?),')") WRITE (*, "(' and with an improved method of extrapolating the smoothed-seismicity')") WRITE (*, "(' parent forecast to high magnitudes.')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' The primary computational steps are (in order):')") WRITE (*, "(' (1) Read the Kagan and Jackson smoothed-Seismicity long-term forecast')") WRITE (*, "(' grid and then scale it to 31 different threshold magnitudes;')") WRITE (*, "(' (2) Recreate the Tectonics forecast (SHIFT_GSRM2f) of Bird and Kreemer,')") WRITE (*, "(' exactly as in program SHIFT_GSRM2f_for_CSEP.f90,')") WRITE (*, "(' with forecast spatial grids for each of 31 threshold magnitudes;')") WRITE (*, "(' (3) Merge these 2 parent forecasts into hybrid forecast H*,')") WRITE (*, "(' with spatial grids for each of 31 threshold magnitudes;')") WRITE (*, "(' (4) Difference these grids to get seismicity in distinct magnitude bins.')") WRITE (*, "(' (5) Output the forecast in XML format.')") WRITE (*, "(' Run time is about ~5 hours if 16-fold parallel computing is available')") WRITE (*, "(' during step (1); otherwise, a few days may be needed.')") WRITE (*, "(' It is expected that 64-bit computing will be used because of the')") WRITE (*, "(' large memory requirements (4 GB) of this code.')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' Note that random-number generators used in step (2) can come from')") WRITE (*, "(' either the Intel Math Kernel Library (MKL), or from the ')") WRITE (*, "(' International Mathematics Subroutine Library (IMSL).')") WRITE (*, "(' A change-over would be effected by switching which lines of ')") WRITE (*, "(' source code are commented-out, versus active.')") WRITE (*, "(' Search on ""MKL"" and also on ""IMSL"" to find these lines.')") WRITE (*, "(' The method of initialization of these generators is not important.')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' This program produces output in the format of the Global Testing Region')") WRITE (*, "(' at the Collaboratory for the Study of Earthquake Predictability (CSEP).')") WRITE (*, "(' Specifically, it yields resolution of 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 (*, "(' The start-date and end-date of the forecast are specified by the user')") WRITE (*, "(' in the GEAR1_parameter.dat file (which is read from the command line);')") WRITE (*, "(' however, conceptually the forecast is open-ended toward the future.')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' This source code (with 7 input datasets) is the primary means of ')") WRITE (*, "(' distributing the forecast, for several reasons:')") WRITE (*, "(' (a) CSEP protocol requires submission of source code;')") WRITE (*, "(' (b) We support maximum transparency and open-source ethics; and')") WRITE (*, "(' (c) The output file is much larger (3.7 GB) than the sum of the ')") WRITE (*, "(' source code and input files (about 470 MB, or 0.46 GB).')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) WRITE (*, "(' This version was written 2014.06.04+ by Peter Bird at UCLA,')") WRITE (*, "(' using large chunks of code from the previous programs')") WRITE (*, "(' S_method2E_demo.f90, SHIFT_GSRM2f_for_CSEP.f90,')") WRITE (*, "(' and hybrid_forecasts_v3.f90.')") WRITE (*, *) IF (pausing) CALL Pause() WRITE (*, *) OPEN (UNIT = 21, FILE = "GEAR1_for_CSEP_log.txt") ! absolute OPEN; overwrites any existing WRITE (21, "('PROGRAM GEAR1_for_CSEP:')") WRITE (21, *) WRITE (21, "('Computes the Global Earthquake Activity Rate (GEAR) model, version 1,')") WRITE (21, "('of long-term rates of shallow earthquakes on a global grid, as functions')") WRITE (21, "('of longitude, latitude, and magnitude, as described by')") WRITE (21, "('P. Bird, D.D. Jackson, Y.Y. Kagan, C. Kreemer, and R.S. Stein')") WRITE (21, "('[2015, Seismol. Soc. Am. Bull.]')") WRITE (21, "('This model is an optimal hybrid of the smoothed-Seismicity forecast')") WRITE (21, "('of Kagan and Jackson [1994, 2010, 2011] with the Tectonics forecast of ')") WRITE (21, "('Bird and Kreemer (SHIFT_GSRM2f, Bird and Kreemer [2015, BSSA]).')") WRITE (21, "('Specifically, it is the preferred hybrid H* described in our paper,')") WRITE (21, "('with input datasets updated to end-2013 (or perhaps later?).')") WRITE (21, "('and with an improved method of extrapolating the smoothed-seismicity')") WRITE (21, "('parent forecast to high magnitudes.')") WRITE (21, *) WRITE (21, "('The primary computational steps are (in order):')") WRITE (21, "('(1) Read the Kagan and Jackson smoothed-Seismicity long-term forecast')") WRITE (21, "(' grid and then scale it to 31 different threshold magnitudes;')") WRITE (21, "('(2) Recreate the Tectonics forecast (SHIFT_GSRM2f) of Bird and Kreemer,')") WRITE (21, "(' exactly as in program SHIFT_GSRM2f_for_CSEP.f90,')") WRITE (21, "(' with forecast spatial grids for each of 31 threshold magnitudes;')") WRITE (21, "('(3) Merge these 2 parent forecasts into hybrid forecast H*,')") WRITE (21, "(' with spatial grids for each of 31 threshold magnitudes;')") WRITE (21, "('(4) Difference these grids to get seismicity in distinct magnitude bins.')") WRITE (21, "('(5) Output the forecast in XML format.')") WRITE (21, "('Run time is about ~5 hours if 16-fold parallel computing is available')") WRITE (21, "('during step (1); otherwise, a few days may be needed.')") WRITE (21, "('It is expected that 64-bit computing will be used because of the')") WRITE (21, "('large memory requirements (4 GB) of this code.')") WRITE (21, *) WRITE (21, "('Note that random-number generators used in step (2) can come from')") WRITE (21, "('either the Intel Math Kernel Library (MKL), or from the ')") WRITE (21, "('International Mathematics Subroutine Library (IMSL).')") WRITE (21, "('A change-over would be effected by switching which lines of ')") WRITE (21, "('source code are commented-out, versus active.')") WRITE (21, "('Search on ""MKL"" and also on ""IMSL"" to find these lines.')") WRITE (21, "('The method of initialization of these generators is not important.')") WRITE (21, *) WRITE (21, "('This program produces output in the format of the Global Testing Region')") WRITE (21, "('at the Collaboratory for the Study of Earthquake Predictability (CSEP).')") WRITE (21, "('Specifically, it yields resolution of 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, "('The start-date and end-date of the forecast are specified by the user')") WRITE (21, "('in the GEAR1_parameter.dat file (which is read from the command line);')") WRITE (21, "('however, conceptually the forecast is open-ended toward the future.')") WRITE (21, *) WRITE (21, "('This source code (with 7 input datasets) is the primary means of ')") WRITE (21, "('distributing the forecast, for several reasons:')") WRITE (21, "('(a) CSEP protocol requires submission of source code;')") WRITE (21, "('(b) We support maximum transparency and open-source ethics; and')") WRITE (21, "('(c) The output file is much larger (3.7 GB) than the sum of the ')") WRITE (21, "(' source code and input files (about 470 MB, or 0.46 GB).')") WRITE (21, *) WRITE (21, "('This version was written 2014.06.04+ by Peter Bird at UCLA,')") WRITE (21, "('using large chunks of code from the previous programs')") WRITE (21, "('S_method2E_demo.f90, SHIFT_GSRM2f_for_CSEP.f90,')") WRITE (21, "('and hybrid_forecasts_v3.f90.')") WRITE (21, *) !------------------------------------------------------------------------------------- !Find, read, and memorize contents of the parameter file: ! GEAR1_parameters.dat should have contents similar to this: ! ---------------------------------------------------------- ! smoothedSeismicityDatPath=GL_HAZTBLT_M5_B2_2013.TMP ! tectonicZonesDatPath=PB2002_tectonic_zones.grd ! plateBoundariesDatPath=PB2002_steps.dat ! 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=GEAR1.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 (*, "(' GEAR1_parameters.dat was supplied to CSEP.')") WRITE (21, "('ERROR: [/path/]filename of parameter file must follow in command line!')") WRITE (21, "('GEAR1_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 = "GEAR1_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: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios 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:26) == "smoothedSeismicityDatPath=") THEN smoothedSeismicityDatPath = TRIM(text_line(27:text_line_length)) ELSE IF (text_line(1:21) == "tectonicZonesDatPath=") THEN tectonicZonesDatPath = TRIM(text_line(22:text_line_length)) ELSE 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) == "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") ! GEAR1_parameters.dat !========================================================================================================= WRITE (*, *) WRITE (*, "(' =======================================================================')") WRITE (*, "(' === STEP #1: Read smoothed-Seismicity & extrapolate ================')") WRITE (*, "(' =======================================================================')") WRITE (*, *) WRITE (21, *) WRITE (21, "('=======================================================================')") WRITE (21, "('=== STEP #1: Read smoothed-Seismicity & extrapolate ================')") WRITE (21, "('=======================================================================')") WRITE (21, *) !========================================================================================================= !Date from Table 1 of Kagan et al. [2010, Pure & Appl. Geophys.]: corner_magnitude_of_zone(0) = 7.72 ! Was 8.15 in Table 1 of Kagan et al.; but, here zones 0 & 1 are merged. corner_magnitude_of_zone(1) = 7.72 ! Was 7.59 in Table 1 of Kagan et al.; but, here zones 0 & 1 are merged. corner_magnitude_of_zone(2) = 7.38 corner_magnitude_of_zone(3) = 6.79 corner_magnitude_of_zone(4) = 9.50 ! Was 8.75 in Table 1 of Kagan et al.; but 9.5 supported by Bird and Kagan [2004]; Kagan and Jackson [2013]. beta_of_zone(0) = 0.645 ! Was 0.639 in Table 1 of Kagan et al.; but, here zones 0 & 1 are merged. beta_of_zone(1) = 0.645 ! Was 0.647 in Table 1 of Kagan et al.; but, here zones 0 & 1 are merged. beta_of_zone(2) = 0.812 beta_of_zone(3) = 0.767 beta_of_zone(4) = 0.639 smoothing_distance_km_int = NINT(smoothing_distance_m / 1000.) ! refers to REAL, PARAMETER :: smoothing_distance_m at top; probably 200000. WRITE (smoothing_distance_km_c3, "(I3)") smoothing_distance_km_int IF (smoothing_distance_km_c3(1:1) == ' ') smoothing_distance_km_c3(1:1) = '_' IF (smoothing_distance_km_c3(2:2) == ' ') smoothing_distance_km_c3(2:2) = '_' OPEN (UNIT = 1, FILE = TRIM(smoothedSeismicityDatPath), STATUS = "OLD", DISP = "KEEP", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (*, "(' ERROR: File ',A,' was not found. Try again...')") TRIM(smoothedSeismicityDatPath) WRITE (21, "('ERROR: File ',A,' was not found. Try again...')") TRIM(smoothedSeismicityDatPath) IF (pausing) CALL Pause() STOP END IF !Define the fine grid of Kagan and Jackson (which will also be the GEAR1 output grid): fine_d_lat = 0.1 fine_d_lon = 0.1 fine_lon_min = -179.95 fine_lon_max = +179.95 fine_lat_min = -89.95 fine_lat_max = +89.95 !Decide size of fine grids, and allocate them: fine_nRows = 1 + NINT((fine_lat_max - fine_lat_min) / fine_d_lat) fine_nColumns = 1 + NINT((fine_lon_max - fine_lon_min) / fine_d_lon) ALLOCATE ( S_parent_fine( fine_nRows, fine_nColumns) ) S_parent_fine = 0.0D0 ! whole array ALLOCATE ( beta_fine( fine_nRows, fine_nColumns) ) ALLOCATE ( smoothed_beta_fine( fine_nRows, fine_nColumns) ) ALLOCATE ( corner_magnitude_fine( fine_nRows, fine_nColumns) ) ALLOCATE ( smoothed_corner_magnitude_fine( fine_nRows, fine_nColumns) ) !Compute areas of fine grid cells and check the total: ALLOCATE ( fine_cell_area_m2(fine_nRows) ) total_area_m2 = 0.0D0 DO i = 1, fine_nRows ! N to S top_lat = fine_lat_max + 0.5 * fine_d_lat - (i-1) * fine_d_lat bottom_lat = top_lat - fine_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)) fine_cell_area_m2(i) = strip_area_steradian * R2 / fine_nColumns total_area_m2 = total_area_m2 + fine_nColumns * fine_cell_area_m2(i) END DO WRITE (*, "(' ')") WRITE (*, "(' Total fine grid area = ',ES12.5,' square meters,')") total_area_m2 WRITE (21, *) WRITE (21, "('Total fine 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 WRITE (*, *) !Temporarily borrow one of these empty arrays, to recompute area of the fine grid: beta_fine = 1.0 ! whole grid CALL Integrate_REAL_Grid(beta_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & fine_grid_area_m2) ! output WRITE (*, "(' Fine grid area from Integrate_REAL_Grid = ',ES12.4,' m^2.')") fine_grid_area_m2 WRITE (21, "('Fine grid area from Integrate_REAL_Grid = ',ES12.4,' m^2.')") fine_grid_area_m2 WRITE (*, *) WRITE (*, "(' Expecting 16 lines of column headers, shown here (between horizontal lines): ')") WRITE (*, "(' -----------------------------------------------------------------------------')") WRITE (21, *) WRITE (21, "('Expecting 16 lines of column headers, shown here (between horizontal lines): ')") WRITE (21, "('-----------------------------------------------------------------------------')") DO i = 1, 16 READ (1, "(A)") lines(i) WRITE (*, "(' ',A)") TRIM(lines(i)) WRITE (21, "(A)") TRIM(lines(i)) END DO WRITE (*, "(' -----------------------------------------------------------------------------')") WRITE (21, "('-----------------------------------------------------------------------------')") READ (lines(4), "(26X,F3.1)") S_threshold_magnitude_Kagan S_threshold_magnitude_Bird = S_threshold_magnitude_Kagan - 0.033 ! Kagan uses constant = 9; Bird uses constant = 9.05 (for moment in N m; before 2/3 factor) ! Each of us has built this formula into so many of our programs ! that it is easier to make adjustments (if necessary) during publication. WRITE (*, *) WRITE (*, "(' Threshold magnitude is ',F3.1,' (Kagan formula), or ',F5.3,' (Bird formula).')") & & S_threshold_magnitude_Kagan, S_threshold_magnitude_Bird WRITE (21, *) WRITE (21, "('Threshold magnitude is ',F3.1,' (Kagan formula), or ',F5.3,' (Bird formula).')") & & S_threshold_magnitude_Kagan, S_threshold_magnitude_Bird WRITE (*, "(' Reading ',A,',')") TRIM(smoothedSeismicityDatPath) WRITE (*, "(' which may take some minutes...')") WRITE (21, "('Reading ',A,',')") TRIM(smoothedSeismicityDatPath) WRITE (21, "(' which may take some minutes...')") DO i = 1, fine_nRows latitude = fine_lat_max - (i-1) * fine_d_lat DO j = 1, fine_nColumns longitude = fine_lon_min + (j-1) * fine_d_lon READ (1, *, IOSTAT = ios) lon, lat, value IF (ios /= 0) THEN WRITE (*, "(' ERROR in reading rates at latitude ',F6.1,', longitude ',F6.1,'; IOSTAT = ',I6)") latitude, longitude, ios WRITE (21, "('ERROR in reading rates at latitude ',F6.1,', longitude ',F6.1,'; IOSTAT = ',I6)") latitude, longitude, ios IF (pausing) CALL Pause() STOP END IF IF (ABS(lat - latitude) > 0.001) THEN WRITE (*, "(' ERROR: Read lat ',F6.1,' when expecting latitude ',F6.1)") lat, latitude WRITE (21, "('ERROR: Read lat ',F6.1,' when expecting latitude ',F6.1)") lat, latitude IF (pausing) CALL Pause() STOP END IF IF (ABS(lon - longitude) > 0.001) THEN WRITE (*, "(' ERROR: Read lon ',F6.1,' when expecting longitude ',F6.1)") lon, longitude WRITE (21, "('ERROR: Read lon ',F6.1,' when expecting longitude ',F6.1)") lon, longitude IF (pausing) CALL Pause() STOP END IF S_parent_fine(i, j) = value END DO END DO CLOSE (1) WRITE (*, "(' File was read successfully.')") WRITE (21, "('File was read successfully.')") S_parent_fine = S_parent_fine / (1.0E6 * 24 * 60 * 60) !from per square kilometer to per square meter, and from per day to per second. CALL Integrate_DP_Grid(S_parent_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & S_integral) ! output; should be in /s S_integral_per_year = S_integral * s_per_year WRITE (*, "(' Integral is ',F8.2,' shallow EQs above threshold per year.')") S_integral_per_year WRITE (21, "('Integral is ',F8.2,' shallow EQs above threshold per year.')") S_integral_per_year S_lower_limit_at_m5p8 = S_parent_fine(1, 1) DO i = 1, fine_nRows DO j = 1, fine_nColumns S_lower_limit_at_m5p8 = MIN(S_lower_limit_at_m5p8, S_parent_fine(i, j)) END DO END DO WRITE (*, *) WRITE (*, "(' Epicentral rate density at input threshold has lower limit of ',ES10.3)") S_lower_limit_at_m5p8 WRITE (21, *) WRITE (21, "('Epicentral rate density at input threshold has lower limit of ',ES10.3)") S_lower_limit_at_m5p8 !------------------------------------------------------------------------------------------------------- !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: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (*, "(' ERROR: Copy of file ',A,' was not found in active folder.')") TRIM(tectonicZonesDatPath) WRITE (21, "('ERROR: Copy of file ',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 (*, *) 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, *) 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.')") !------------------------------------------------------------------------------------------------------------ !LOOK-UP loop on all cells, to define tapered G-R parameters (prior to smoothing): DO i = 1, fine_nRows lat = fine_lat_max - (i-1) * fine_d_lon DO j = 1, fine_nColumns lon = fine_lon_min + (j-1) * fine_d_lon !adjust lon to 0~360, to agree with PB2002_tectonic_zones.GRD: IF (lon < 0.) lon = lon + 360.0 !find tectonic zone: iTarget = 1 + NINT((TZ_lat_max - lat) / TZ_d_lat) iTarget = MAX(iTarget, 1) iTarget = MIN(iTarget, nTZ_rows) jTarget = 1 + NINT((lon - TZ_lon_min) / TZ_d_lon) jTarget = MAX(jTarget, 1) jTarget = MIN(jTarget, nTZ_columns) iZone = tectonic_zones_I1(iTarget, jTarget) !SAVE tapered G-R parameters for this cell: beta_fine(i, j) = beta_of_zone(iZone) corner_magnitude_fine(i, j) = corner_magnitude_of_zone(iZone) END DO END DO !SMOOTH & STRETCH the maps of corner magnitude & beta (in the same ways): WRITE (*, *) WRITE (*, "(' Before smoothing, corner_magnitude_fine grid has')") WRITE (21, *) WRITE (21, "('Before smoothing, corner_magnitude_fine grid has')") CALL Scan_REAL_Grid(corner_magnitude_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_mc0, mean_mc0, standard_deviation_mc0, highest_mc0) ! WRITE (*, "(' low =',ES12.4,', mean+-s.d. ='p,ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc0, mean_mc0, standard_deviation_mc0, highest_mc0 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc0, mean_mc0, standard_deviation_mc0, highest_mc0 ! low = 6.7900E+00, mean+-s.d. = 7.8056E+00+- 4.6290E-01, high = 9.5000E+00 !CALL WRITE_REAL_GRD(corner_magnitude_fine, "S_method2_corner_magnitude.grd") IF (smoothing_distance_m > 0.0) THEN CALL Smooth_REAL_on_Sphere(corner_magnitude_fine, & ! <=== unsmoothed input grid & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & smoothing_distance_m, Earth_radius_m, & ! smoothing distance & Earth radius & smoothed_corner_magnitude_fine) ! <=== product WRITE (*, "(' After smoothing, smoothed_corner_magnitude_fine grid has')") WRITE (21, "('After smoothing, smoothed_corner_magnitude_fine grid has')") CALL Scan_REAL_Grid(smoothed_corner_magnitude_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_mc1, mean_mc1, standard_deviation_mc1, highest_mc1) WRITE (*, "(' low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc1, mean_mc1, standard_deviation_mc1, highest_mc1 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc1, mean_mc1, standard_deviation_mc1, highest_mc1 ! 200 km: low = 7.1363E+00, mean+-s.d. = 7.8056E+00+- 3.2164E-01, high = 9.5000E+00 stretching_factor = standard_deviation_mc0 / standard_deviation_mc1 ! > 1.0 WRITE (*, "(' Applying stretching factor of ',F6.3,' about the mean.')") stretching_factor WRITE (21, "('Applying stretching factor of ',F6.3,' about the mean.')") stretching_factor ! 200 km: Applying stretching factor of 1.439 about the mean. DO i = 1, fine_nRows DO j = 1, fine_nColumns smoothed_corner_magnitude_fine(i, j) = mean_mc0 + (smoothed_corner_magnitude_fine(i, j) - mean_mc0) * stretching_factor END DO END DO WRITE (*, "(' After stretching, smoothed_corner_magnitude_fine grid has')") WRITE (21, "('After stretching, smoothed_corner_magnitude_fine grid has')") CALL Scan_REAL_Grid(smoothed_corner_magnitude_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_mc2, mean_mc2, standard_deviation_mc2, highest_mc2) WRITE (*, "(' low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc2, mean_mc2, standard_deviation_mc2, highest_mc2 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_mc2, mean_mc2, standard_deviation_mc2, highest_mc2 ! 200 km: low = 6.8424E+00, mean+-s.d. = 7.8056E+00+- 4.6290E-01, high = 1.0244E+01 !output_pathfile = "S_method2E" // smoothing_distance_km_c3 // "km_corner_magnitude.grd" !CALL WRITE_REAL_GRD(smoothed_corner_magnitude_fine, TRIM(output_pathfile)) ELSE smoothed_corner_magnitude_fine = corner_magnitude_fine END IF WRITE (*, *) WRITE (*, "(' Before smoothing, beta_fine grid has')") WRITE (21, *) WRITE (21, "('Before smoothing, beta_fine grid has')") CALL Scan_REAL_Grid(beta_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_beta0, mean_beta0, standard_deviation_beta0, highest_beta0) WRITE (*, "(' low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta0, mean_beta0, standard_deviation_beta0, highest_beta0 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta0, mean_beta0, standard_deviation_beta0, highest_beta0 ! low = 6.3900E-01, mean+-s.d. = 6.5149E-01+- 3.1591E-02, high = 8.1200E-01 !CALL WRITE_REAL_GRD(beta_fine, "S_method2_beta.grd") IF (smoothing_distance_m > 0.0) THEN CALL Smooth_REAL_on_Sphere(beta_fine, & ! <=== unsmoothed input grid & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & smoothing_distance_m, Earth_radius_m, & ! smoothing distance & Earth radius & smoothed_beta_fine) ! <=== product WRITE (*, "(' After smoothing, smoothed_beta_fine grid has')") WRITE (21, "('After smoothing, smoothed_beta_fine grid has')") CALL Scan_REAL_Grid(smoothed_beta_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_beta1, mean_beta1, standard_deviation_beta1, highest_beta1) WRITE (*, "(' low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta1, mean_beta1, standard_deviation_beta1, highest_beta1 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta1, mean_beta1, standard_deviation_beta1, highest_beta1 ! 200 km: low = 6.3900E-01, mean+-s.d. = 6.5149E-01+- 1.8578E-02, high = 7.9756E-01 ! 300 km: low = 6.3900E-01, mean+-s.d. = 6.5149E-01+- 1.5327E-02, high = 7.7035E-01 stretching_factor = standard_deviation_beta0 / standard_deviation_beta1 ! > 1.0 WRITE (*, "(' Applying stretching factor of ',F6.3,' about the mean.')") stretching_factor WRITE (21, "('Applying stretching factor of ',F6.3,' about the mean.')") stretching_factor ! 200 km: Applying stretching factor of 1.700 about the mean. ! 300 km: Applying stretching factor of 2.061 about the mean. DO i = 1, fine_nRows DO j = 1, fine_nColumns smoothed_beta_fine(i, j) = mean_beta0 + (smoothed_beta_fine(i, j) - mean_beta0) * stretching_factor END DO END DO WRITE (*, "(' After stretching, smoothed_beta_fine grid has')") WRITE (21, "('After stretching, smoothed_beta_fine grid has')") CALL Scan_REAL_Grid(smoothed_beta_fine, & ! <=== grid to be summarized & fine_lon_min, fine_d_lon, fine_lon_max, & ! grid description & fine_lat_min, fine_d_lat, fine_lat_max, & ! grid description & lowest_beta2, mean_beta2, standard_deviation_beta2, highest_beta2) WRITE (*, "(' low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta2, mean_beta2, standard_deviation_beta2, highest_beta2 WRITE (21, "('low =',ES12.4,', mean+-s.d. =',ES12.4,'+-',ES11.4,', high =',ES12.4)") lowest_beta2, mean_beta2, standard_deviation_beta2, highest_beta2 ! 200 km: low = 6.3025E-01, mean+-s.d. = 6.5149E-01+- 3.1591E-02, high = 8.9988E-01 ! 300 km: low = 6.2575E-01, mean+-s.d. = 6.5149E-01+- 3.1591E-02, high = 8.9649E-01 !output_pathfile = "S_method2E" // smoothing_distance_km_c3 // "km_beta.grd" !CALL WRITE_REAL_GRD(smoothed_beta_fine, TRIM(output_pathfile)) ELSE smoothed_beta_fine = beta_fine END IF !-------------------------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' About to allocate array S_parent_fine_3D (1.6 GB) ...')") WRITE (21, *) WRITE (21, "('About to allocate array S_parent_fine_3D (1.6 GB) ...')") ALLOCATE ( S_parent_fine_3D( magnitude_bins, fine_nRows, fine_nColumns) ) ! created here in Step #1; used in Step #3 !PRODUCTION loops on all thresholds: DO m = 1, magnitude_bins F_threshold_magnitude = magnitude_bin_threshold(m) ! referring to CSEP's standard lower-limits of magnitude bins F_threshold_moment = Moment(F_threshold_magnitude) WRITE (*, "(' Extrapolating S forecast to threshold magnitude of ',F5.2,'...')") F_threshold_magnitude WRITE (21, "('Extrapolating S forecast to threshold magnitude of ',F5.2,'...')") F_threshold_magnitude !Decide what lower limit on seismicity will apply: corner_magnitude = 8.15 ! tectonic zone 0 (plate interior) beta = 0.639 ! tectonic zone 0 (plate interior) argument = (Moment(5.767) - F_threshold_moment) / Moment(corner_magnitude) IF (argument < -100.0D0) THEN G_function = 0.0D0 WRITE (*, "(' ERROR: G_function = 0.0 during lower-limit computation.')") WRITE (21, "('ERROR: G_function = 0.0 during lower-limit computation.')") IF (pausing) CALL Pause() STOP ELSE G_function = ((F_threshold_moment / Moment(5.767))**(-beta)) * & & EXP(argument) END IF S_lower_limit_at_high_m = S_lower_limit_at_m5p8 * G_function !PRODUCTION LOOP on all cells: n_TINY = 0 low_argument_count = 0 DO i = 1, fine_nRows lat = fine_lat_max - (i-1) * fine_d_lon DO j = 1, fine_nColumns lon = fine_lon_min + (j-1) * fine_d_lon !adjust lon to 0~360, to agree with PB2002_tectonic_zones.GRD: IF (lon < 0.) lon = lon + 360.0 ! Correct seismicity to desired_threshold_moment_Nm, ! per equation (5) of Bird and 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 corner_magnitude = smoothed_corner_magnitude_fine(i, j) beta = smoothed_beta_fine(i, j) !First, extrapolate according to tectonic zone, without applying limit: argument = (Moment(5.767) - F_threshold_moment) / Moment(corner_magnitude) IF (argument < -100.0D0) THEN low_argument_count = low_argument_count + 1 G_function = 0.0D0 IF (low_argument_count < 10) THEN !WRITE (*, *) !WRITE (*, "(' CAUTION: G_function = 0.0 in row ',I6,' and column ',I6)") i, j !WRITE (*, "(' Presumably this zero result will be replaced by the lower limit.')") !WRITE (21, *) !WRITE (21, "('CAUTION: G_function = 0.0 in row ',I6,' and column ',I6)") i, j !WRITE (21, "('Presumably this zero result will be replaced by the lower limit.')") ELSE IF (low_argument_count == 10) THEN !WRITE (*, *) !WRITE (*, "(' No more CAUTION statements will be printed.')") !WRITE (21, *) !WRITE (21, "('No more CAUTION statements will be printed.')") END IF ELSE G_function = ((F_threshold_moment / Moment(5.767))**(-beta)) * & & EXP(argument) END IF !Second, apply the lower limit: S_parent_fine_3D(m, i, j) = MAX((S_parent_fine(i, j) * G_function), S_lower_limit_at_high_m) IF (S_parent_fine_3D(m, i, j) < 1.D-37) n_TINY = n_TINY + 1 END DO ! j = 1, fine_nColumns END DO ! i = 1, fine_nRows !IF (n_TINY > 0) THEN ! WRITE (*, *) ! WRITE (*, "(' CAUTION: ',I10,' seismicity densities were less than 1.0D-37;')") n_TINY ! WRITE (*, "(' and these will have to be rounded-up for display in NeoKineMap.')") ! !IF (pausing) CALL Pause() !END IF !------------------------------------------------------------------------------------------------------------ !Write S_8p00method2E200km_2014plus.GRD (or similar file-name): !WRITE (mc_c5, "(F5.2)") F_threshold_magnitude !mc_c5(3:3) = 'p' ! replacing '.' !IF(mc_c5(1:1) == ' ') mc_c5(1:1) = '_' !output_pathfile = "S" // mc_c5 // "method2E" // smoothing_distance_km_c3 // "km_2014plus.GRD" !WRITE (*, "(' Writing output GRD file: ',A)") TRIM(output_pathfile) !CALL WRITE_DP_GRD(F_seismicity, TRIM(output_pathfile)) !------------------------------------------------------------------------------------------------------------ END DO ! m = 1, magnitude_bins WRITE (*, *) WRITE (*, "(' STEP #1 completed.')") WRITE (21, *) WRITE (21, "('STEP #1 completed.')") !========================================================================================================= WRITE (*, *) WRITE (*, "(' =======================================================================')") WRITE (*, "(' === STEP #2: Compute Tectonics forecast SHIFT_GSRM2f ================')") WRITE (*, "(' =======================================================================')") WRITE (*, *) WRITE (21, *) WRITE (21, "('=======================================================================')") WRITE (21, "('=== STEP #2: Compute Tectonics forecast SHIFT_GSRM2f ================')") WRITE (21, "('=======================================================================')") WRITE (21, *) !========================================================================================================= !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 WRITE (*, *) WRITE (*, "(' About to allocate array T_parent_coarse_3D (0.4 GB) ...')") WRITE (21, *) WRITE (21, "('About to allocate array T_parent_coarse_3D (0.4 GB) ...')") ALLOCATE ( T_parent_coarse_3D(magnitude_bins, nRows, nColumns) ) T_parent_coarse_3D = 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 coarse grid cells and check the total: ALLOCATE ( coarse_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)) coarse_cell_area_m2(i) = strip_area_steradian * R2 / nColumns total_area_m2 = total_area_m2 + nColumns * coarse_cell_area_m2(i) END DO WRITE (*, "(' ')") WRITE (*, "(' Total coarse grid area = ',ES12.5,' square meters,')") total_area_m2 WRITE (21, *) WRITE (21, "('Total coarse 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 (*, *) !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: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios 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 + coarse_cell_area_m2(i) ELSE intraplate_area_m2 = intraplate_area_m2 + coarse_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: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios 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) OPEN (UNIT = 1, FILE = TRIM(underwaterDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios 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) 210 WRITE (*, *) WRITE (*, "(' FILE ',A,' must be available. Please check this.')") TRIM(plateBoundariesDatPath) OPEN (UNIT = 1, FILE = TRIM(plateBoundariesDatPath), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: During file OPEN, return code IOSTAT = ',I12)") ios WRITE (21, "('ERROR: During file OPEN, return code IOSTAT = ',I12)") ios 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 and 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...')") WRITE (*, *) 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 (*, *) 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 and 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 and 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 and 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 = (coarse_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, "('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): WRITE (*, "(' Random numbers generated with brng = VSL_BRNG_MCG31 = ',I12,', seed = 1')") VSL_BRNG_MCG31 WRITE (21, "('Random numbers generated with brng = VSL_BRNG_MCG31 = ',I12,', seed = 1')") VSL_BRNG_MCG31 WRITE (*, *) ! to prevent overwriting previous message by "% complete" lines to come, below ... WRITE (21, *) status = VslNewStream( stream = my_stream, brng = VSL_BRNG_MCG31, seed = 1) IF (status /= VSL_STATUS_OK) THEN WRITE (*, "(' ERROR: When initializing my_stream with Intel MKL/VSL, status = ',I12)") status WRITE (*, "(' but desired status is VSL_STATUS_OK = ',I12)") VSL_STATUS_OK WRITE (21, "('ERROR: When initializing my_stream with Intel MKL/VSL, status = ',I12)") status WRITE (21, "(' but desired status is VSL_STATUS_OK = ',I12)") VSL_STATUS_OK 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 (*, *) 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] ------------------------------------------- !NOTE: PB2002_tectonic_zones.grd of Bird et al. [2010, Pure Appl. Geoph.] was already read & stored in Step #1." !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...')") WRITE (*, *) 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 (*, *) 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) + coarse_cell_area_m2(i) ELSE intraplate_area_m2 = intraplate_area_m2 + coarse_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 * coarse_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 T_parent_coarse_3D(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 T_parent_coarse_3D(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 !========================================================================================================= WRITE (*, *) WRITE (*, "(' =======================================================================')") WRITE (*, "(' === STEP #3: Merge to create hybrid forecast H* ================')") WRITE (*, "(' =======================================================================')") WRITE (*, *) WRITE (21, *) WRITE (21, "('=======================================================================')") WRITE (21, "('=== STEP #3: Merge to create hybrid forecast H* ================')") WRITE (21, "('=======================================================================')") WRITE (21, *) !========================================================================================================= WRITE (*, *) WRITE (*, "(' About to allocate array H_star_fine_3D (1.6 GB) ...')") WRITE (21, *) WRITE (21, "('About to allocate array H_star_fine_3D (1.6 GB) ...')") ALLOCATE ( H_star_fine_3D(magnitude_bins, fine_nRows, fine_nColumns) ) ! created here in Step #3; transformed in Step #4; output in Step #5 !Allocate smaller 2D arrays (less likely to cause an abend): ALLOCATE ( T_parent_fine(fine_nRows, fine_nColumns) ) ! used to resample one map at a time from T_parent_coarse_3D ALLOCATE ( floor_seismicity(fine_nRows, fine_nColumns) ) ! used only in this step ALLOCATE ( H_star_fine(fine_nRows, fine_nColumns) ) ! temporary storage for one map, which will be later loaded into H_star_fine_3D. DO m = 1, magnitude_bins !Copy 2-D S_parent_fine map from S_parent_fine_3D: DO i = 1, fine_nRows DO j = 1, fine_nColumns S_parent_fine(i, j) = S_parent_fine_3D(m, i, j) END DO END DO CALL Integrate_DP_Grid(S_parent_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & S_integral) ! output S_integral_per_year = S_integral * s_per_year WRITE (*, *) WRITE (*, "(' Resampling from T_parent_coarse_3D(m >',F6.3,') into 2-D T_parent_fine ...')") magnitude_bin_threshold(m) WRITE (21, *) WRITE (21, "('Resampling from T_parent_coarse_3D(m >',F6.3,') into 2-D T_parent_fine ...')") magnitude_bin_threshold(m) T_parent_fine = 0.0D0 ! whole 2-D array (to simplify debugging) DO i = 1, fine_nRows latitude = fine_lat_max - (i-1) * fine_d_lat 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 j = 1, fine_nColumns longitude = fine_lon_min + (j-1) * fine_d_lon IF ((1 + MOD((j-1), 5)) == 3) THEN ! fine cell centered on the joint between 2 coarse cells, to W and E jTargetW = 1 + NINT((longitude - (0.25 * fine_d_lon) - GSRM2_lon_min) / GSRM2_d_lon) jTargetE = jTargetW + 1 T_parent_fine(i, j) = 0.5 * T_parent_coarse_3D(m, iTarget, jTargetW) + 0.5 * T_parent_coarse_3D(m, iTarget, jTargetE) ELSE ! fine cell is entirely within one coarse cell jTarget = 1 + NINT((longitude - GSRM2_lon_min) / GSRM2_d_lon) T_parent_fine(i, j) = T_parent_coarse_3D(m, iTarget, jTarget) END IF ! hard or easy case END DO ! j = 1, fine_nColumns END DO ! i = 1, fine_nRows CALL Integrate_DP_Grid(T_parent_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & T_integral) ! output T_integral_per_year = T_integral * s_per_year WRITE (*, "(' After resampling to fine grid,')") WRITE (*, "(' Integral of fine T grid is ',F9.3,' shallow EQs above threshold per year.')") T_integral_per_year WRITE (21, "('After resampling to fine grid,')") WRITE (21, "('Integral of fine T grid is ',F9.3,' shallow EQs above threshold per year.')") T_integral_per_year !******************************************** d = 0.6 !*********************************** !******************************************** !determine desired "floor" value: S_floor = S_parent_fine(1, 1) T_floor = T_parent_fine(1, 1) DO i = 1, fine_nRows DO j = 1, fine_nColumns S_floor = MIN(S_floor, S_parent_fine(i, j)) T_floor = MIN(T_floor, T_parent_fine(i, j)) END DO END DO H_floor = (S_floor**d) * (T_floor**(1.0-d)) ! IFF this is high enough... H_floor = MAX(H_floor, MIN(S_floor, T_floor)) !compute rough forecast, before normalization: DO i = 1, fine_nRows DO j = 1, fine_nColumns H_star_fine(i, j) = (S_parent_fine(i, j)**d) * (T_parent_fine(i, j)**(1.0 - d)) H_star_fine(i, j) = MAX(H_star_fine(i, j), H_floor) END DO END DO !------------------------------------------------------------------- ! Decide target global shallow earthquake rate R, for this kind of hybrid: target_N_per_year = (S_integral_per_year**d) * (T_integral_per_year**(1.0 - d)) WRITE (*, "(' Hybrid global rate R_H* set to ',F9.3,' shallow EQs above threshold /year.')") target_N_per_year WRITE (21, "('Hybrid global rate R_H* set to ',F9.3,' shallow EQs above threshold /year.')") target_N_per_year !------------------------------------------------------------------- !scale hybrid forecast to desired standard integral, while preserving desired "floor" value: CALL Integrate_DP_Grid(H_star_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & H_integral) ! output H_integral_per_year = H_integral * s_per_year WRITE (*, "(' Raw H'' integral is ',F9.3,' shallow EQs above threshold per year.')") H_integral_per_year WRITE (21, "('Raw H'' integral is ',F9.3,' shallow EQs above threshold per year.')") H_integral_per_year floor_seismicity = H_floor ! whole grid CALL Integrate_DP_Grid(floor_seismicity, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & floor_integral) ! output floor_integral_per_year = floor_integral * s_per_year correction = (target_N_per_year - floor_integral_per_year) / (H_integral_per_year - floor_integral_per_year) H_star_fine = H_floor + correction * (H_star_fine - H_floor) ! whole grid CALL Integrate_DP_Grid(H_star_fine, fine_nRows, fine_d_lat, fine_lat_max, & & fine_nColumns, fine_d_lon, & ! inputs & H_integral) ! output H_integral_per_year = H_integral * s_per_year WRITE (*, "(' Adjusted to ',F9.3,' shallow EQs above threshold per year.')") H_integral_per_year WRITE (21, "('Adjusted to ',F9.3,' shallow EQs above threshold per year.')") H_integral_per_year !Save results in 3D array: DO i = 1, fine_nRows DO j = 1, fine_nColumns H_star_fine_3D(m, i, j) = H_star_fine(i, j) END DO END DO END DO ! m = 1, magnitude_bins WRITE (*, *) WRITE (*, "(' STEP #3 completed.')") WRITE (21, *) WRITE (21, "('STEP #3 completed.')") !========================================================================================================= WRITE (*, *) WRITE (*, "(' =======================================================================')") WRITE (*, "(' === STEP #4: Difference grids to bin by magnitude ================')") WRITE (*, "(' =======================================================================')") WRITE (*, *) WRITE (21, *) WRITE (21, "('=======================================================================')") WRITE (21, "('=== STEP #4: Difference grids to bin by magnitude ================')") WRITE (21, "('=======================================================================')") WRITE (21, *) !========================================================================================================= !----------------------------------------------------------------------------------- ! Convert H_star_fine_3D from all-above-threshold to incremental only-those-within-bin format ! (except for last bin, m = magnitude_bins, which stays open): WRITE (*, *) WRITE (*, "(' Converting H* 3D grid from rates-above-threshold to rates-in-magnitude-bin.')") WRITE (21, *) WRITE (21, "('Converting H* 3D grid from rates-above-threshold to rates-in-magnitude-bin.')") DO m = 1, (magnitude_bins-1) DO i = 1, fine_nRows DO j = 1, fine_nColumns H_star_fine_3D(m, i, j) = H_star_fine_3D(m, i, j) - H_star_fine_3D(m+1, i, j) END DO END DO intraplate_seismicity_SI(m) = intraplate_seismicity_SI(m) - intraplate_seismicity_SI(m+1) END DO !----------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' STEP #4 completed.')") WRITE (21, *) WRITE (21, "('STEP #4 completed.')") !========================================================================================================= WRITE (*, *) WRITE (*, "(' =======================================================================')") WRITE (*, "(' === STEP #5: Output GEAR1 forecast in XML format ================')") WRITE (*, "(' =======================================================================')") WRITE (*, *) WRITE (21, *) WRITE (21, "('=======================================================================')") WRITE (21, "('=== STEP #5: Output GEAR1 forecast in XML format ================')") WRITE (21, "('=======================================================================')") WRITE (21, *) !========================================================================================================= !Integrate seismicity over grid area: WRITE (*, *) WRITE (*, "(' Computing area-integral of GEAR1 seismicity...')") WRITE (21, *) WRITE (21, "('Computing area-integral of GEAR1 seismicity...')") integral = 0.0D0 ! (just initializing) DO i = 1, fine_nRows DO j = 1, fine_nColumns DO m = 1, magnitude_bins integral = integral + H_star_fine_3D(m, i, j) * fine_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',F6.3,'.')") magnitude_bin_threshold(1) 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',F6.3,'.')") magnitude_bin_threshold(1) !Provide global-mean frequency/magnitude curve, in terms of epicentral rate densities: WRITE (*, *) WRITE (*, "(' Computing GEAR1 global-mean frequency-density/magnitude curve, in SI...')") WRITE (21, *) WRITE (21, "('Computing GEAR1 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, fine_nRows DO j = 1, fine_nColumns DO n = m, magnitude_bins integral = integral + H_star_fine_3D(n, i, j) * fine_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, '(" GEAR1")') WRITE (3, '(" 1.0")') ! <=== v1.0=2014.03.21 WRITE (3, '(" Bird, Jackson, Kagan, Kreemer, and Stein [2015, BSSA]")') WRITE (3, '(" smi://org.scec/persons/Bird_Peter")') WRITE (3, '(" 2014-06-04T23: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) 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 !Make sure that longitude stays in range -179.999... to +179.999.... that matches H_star_fine_3D grid: 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) !determine which cell (iTarget, jTarget) of internal grid will have the required epicentral_rate_density? iTarget = 1 + NINT((fine_lat_max - latitude) / fine_d_lat) iTarget = MAX(iTarget, 1) iTarget = MIN(iTarget, fine_nRows) !NOTE: At this time of writing this, I expect iTarget = iCSEP, because CSEP grid is parallel to H_star_fine_3D. DO m = 1, magnitude_bins jTarget = jCSEP ! Much simplified from logic in SHIFT_GSRM2f_for_CSEP.f90, which had to resample at this point! epicenters_per_m2_per_s_inMagBin = H_star_fine_3D(m, iTarget, jTarget) predicted_EQ_count = epicenters_per_m2_per_s_inMagBin * fine_cell_area_m2(iTarget) * 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") ! GEAR1.global.forecast.xml file !----------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' STEP #5 completed.')") WRITE (21, *) WRITE (21, "('STEP #5 completed.')") WRITE (*, *) WRITE (*, "(' PROGRAM GEAR1_for_CSEP: Job completed.')") WRITE (21, *) WRITE (21, "('PROGRAM GEAR1_for_CSEP: Job completed.')") 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 SUBROUTINE Integrate_DP_Grid(seismicity, nRows, d_lat, grid_lat_max, & & nColumns, d_lon, & ! inputs & S_integral) ! output !Integrates grid of seismicity rates; units of answer depend on !units of input. Result is approximately the mean input grid value !times the area of box (around the grid) in square meters. IMPLICIT NONE DOUBLE PRECISION, DIMENSION(:, :), INTENT(IN) :: seismicity INTEGER, INTENT(IN) :: nRows, nColumns REAL, INTENT(IN) :: d_lat, d_lon, grid_lat_max DOUBLE PRECISION, INTENT(OUT) :: S_integral INTEGER :: i, j, use_nColumns LOGICAL :: extra_column REAL, PARAMETER :: Pi = 3.1415927 REAL, PARAMETER :: radians_per_degree = 0.0174533 REAL, PARAMETER :: Earth_radius_m = 6371000.0 REAL :: bottom_lat, cell_area_square_m, strip_area_steradian, top_lat DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 extra_column = (nColumns > NINT(360.0 / d_lon)) IF (extra_column) THEN use_nColumns = nColumns - 1 ELSE use_nColumns = nColumns END IF !Compute cell areas, in m^2, and build the integral: S_integral = 0.0D0 DO i = 1, nRows ! N to S: top_lat = MIN((grid_lat_max + 0.5 * d_lat - (i - 1) * d_lat), 90.0) bottom_lat = MAX((top_lat - d_lat), -90.0) strip_area_steradian = 2 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) & & - 2 * Pi * (unity - DCOS((90.0D0 - top_lat ) * radians_per_degree)) cell_area_square_m = strip_area_steradian * (Earth_radius_m)**2 * d_lon / 360.0 DO j = 1, use_nColumns S_integral = S_integral + cell_area_square_m * seismicity(i, j) END DO END DO END SUBROUTINE Integrate_DP_Grid SUBROUTINE Integrate_REAL_Grid(grid, nRows, d_lat, grid_lat_max, & & nColumns, d_lon, & ! inputs & integral) ! output !Integrates grid of seismicity rates; units of answer depend on !units of input. Result is approximately the mean input grid value !times the area of box (around the grid) in square meters. IMPLICIT NONE REAL, DIMENSION(:, :), INTENT(IN) :: grid INTEGER, INTENT(IN) :: nRows, nColumns REAL, INTENT(IN) :: d_lat, d_lon, grid_lat_max DOUBLE PRECISION, INTENT(OUT) :: integral INTEGER :: i, j, use_nColumns LOGICAL :: extra_column REAL, PARAMETER :: Pi = 3.1415927 REAL, PARAMETER :: radians_per_degree = 0.0174533 REAL, PARAMETER :: Earth_radius_m = 6371000.0 REAL :: bottom_lat, cell_area_square_m, strip_area_steradian, top_lat DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 extra_column = (nColumns > NINT(360.0 / d_lon)) IF (extra_column) THEN use_nColumns = nColumns - 1 ELSE use_nColumns = nColumns END IF !Compute cell areas, in m^2, and build the integral: integral = 0.0D0 DO i = 1, nRows ! N to S: top_lat = MIN((grid_lat_max + 0.5 * d_lat - (i - 1) * d_lat), 90.0) bottom_lat = MAX((top_lat - d_lat), -90.0) strip_area_steradian = 2 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) & & - 2 * Pi * (unity - DCOS((90.0D0 - top_lat ) * radians_per_degree)) cell_area_square_m = strip_area_steradian * (Earth_radius_m)**2 * d_lon / 360.0 DO j = 1, use_nColumns integral = integral + cell_area_square_m * grid(i, j) END DO END DO END SUBROUTINE Integrate_REAL_Grid REAL FUNCTION Moment(magnitude) ! Returns scalar seismic moment in N m, per ! equation (1) of Bird and Kagan [2004; Bull. Seism. Soc. Am.]; ! originally from Hanks and 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 Scan_REAL_Grid(grid, & ! <=== grid to be summarized & lon_min, d_lon, lon_max, & ! grid description & lat_min, d_lat, lat_max, & ! grid description & lowest, mean, standard_deviation, highest) IMPLICIT NONE REAL, DIMENSION(:, :), INTENT(IN) :: grid REAL, INTENT(IN) :: lon_min, d_lon, lon_max, & & lat_min, d_lat, lat_max REAL, INTENT(OUT) :: lowest, mean, standard_deviation, highest INTEGER :: i, j, nRows, nColumns REAL, PARAMETER :: Earth_radius_m = 6371000.0 REAL :: variance REAL, DIMENSION(:, :), ALLOCATABLE :: power_grid DOUBLE PRECISION :: integral nRows = 1 + NINT((lat_max - lat_min) / d_lat) nColumns = 1 + NINT((lon_max - lon_min) / d_lon) ALLOCATE ( power_grid(nRows, nColumns) ) ! for variance/standard_deviation computation highest = grid(1, 1) lowest = grid(1, 1) DO i = 1, nRows DO j = 1, nColumns highest = MAX(highest, grid(i, j)) lowest = MIN(lowest, grid(i, j)) END DO END DO CALL Integrate_REAL_Grid(grid, nRows, d_lat, lat_max, & & nColumns, d_lon, & ! inputs & integral) ! output mean = integral / fine_grid_area_m2 ! <==== global value; must be precomputed DO i = 1, nRows DO j = 1, nColumns power_grid(i, j) = (grid(i, j) - mean)**2 END DO END DO CALL Integrate_REAL_Grid(power_grid, nRows, d_lat, lat_max, & & nColumns, d_lon, & ! inputs & integral) ! output variance = integral / fine_grid_area_m2 standard_deviation = SQRT(variance) END SUBROUTINE Scan_REAL_Grid 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 and 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 and 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 and 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 and 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 SUBROUTINE Smooth_REAL_on_Sphere(grid, & ! <=== unsmoothed input grid & lon_min, d_lon, lon_max, & ! grid description & lat_min, d_lat, lat_max, & ! grid description & smoothing_distance_m, Earth_radius_m, & ! smoothing distance & Earth radius & smoothed_grid) ! <=== product USE Sphere ! .f90, by Peter Bird, for subprogram LonLat_2_Uvec (or, insert same logic here?) IMPLICIT NONE REAL, DIMENSION(:, :), INTENT(IN) :: grid REAL, INTENT(IN) :: lon_min, d_lon, lon_max, & & lat_min, d_lat, lat_max, & & smoothing_distance_m, Earth_radius_m REAL, DIMENSION(:, :), INTENT(OUT) :: smoothed_grid INTEGER :: counter, extent, highest_ii, i, ii, j, jj, lowest_ii, nRows, nColumns, & & old_percent_done, percent_done, small_part REAL :: bottom_lat, G_dA, Gaussian, horizon, horizon_squared, lat, lon, & & radians_away, radians_away_squared, radian_scale, radian_scale_squared, & & source_lat, source_lon, strip_area_steradian, & & target_lat, target_lon, top_lat REAL, DIMENSION(3) :: source_uvec, target_uvec, uvec REAL, DIMENSION(:), ALLOCATABLE :: relative_area ! of each row of the grid; units unimportant REAL, DIMENSION(:), ALLOCATABLE :: numerator_vector, denominator_vector REAL, DIMENSION(:, :, :), ALLOCATABLE :: uvecs DOUBLE PRECISION, PARAMETER :: unity = 1.0D0 !DOUBLE PRECISION, PARAMETER :: Pi = 3.141592654D0 is defined in MODULE Sphere DOUBLE PRECISION :: numerator, denominator WRITE (*, "(' Smoothing w/ Gaussian length-scale',F8.0,' on planet of radius',F9.0,' m ...')") smoothing_distance_m, Earth_radius_m WRITE (21, "('Smoothing w/ Gaussian length-scale',F8.0,' on planet of radius',F9.0,' m ...')") smoothing_distance_m, Earth_radius_m WRITE (*, *) ! preparing for next WRITE (*, "('+ ... radian_scale = smoothing_distance_m / Earth_radius_m radian_scale_squared = radian_scale**2 horizon = 3.0 * radian_scale ! arbitrary, but keeps the operator "local" horizon_squared = horizon**2 extent = 1.0 + (horizon / (d_lat * radians_per_degree)) ! truncated to INTEGER nRows = 1 + NINT((lat_max - lat_min) / d_lat) nColumns = 1 + NINT((lon_max - lon_min) / d_lon) small_part = (nRows * nColumns) / 100 ! approximate definition of each 1% completed ALLOCATE ( relative_area(nRows) ) old_percent_done = -1 DO i = 1, nRows ! N to S: top_lat = MIN((lat_max + 0.5 * d_lat - (i - 1) * d_lat), 90.0) bottom_lat = MAX((top_lat - d_lat), -90.0) strip_area_steradian = 2 * Pi * (unity - DCOS((90.0D0 - bottom_lat) * radians_per_degree)) & & - 2 * Pi * (unity - DCOS((90.0D0 - top_lat ) * radians_per_degree)) relative_area(i) = strip_area_steradian END DO ALLOCATE ( uvecs(3, nRows, nColumns) ) DO i = 1, nRows lat = lat_max - (i-1) * d_lat DO j = 1, nColumns lon = lon_min + (j-1) * d_lon CALL LonLat_2_Uvec(lon, lat, uvec) uvecs(1:3, i, j) = uvec(1:3) END DO END DO ALLOCATE ( numerator_vector(nColumns) ) ALLOCATE ( denominator_vector(nColumns) ) !MAIN TIME-CONSUMING PART (4 nested DO loops!): counter = 0 DO i = 1, nRows ! target row lowest_ii = MAX(1, (i - extent)) highest_ii = MIN(nRows, (i + extent)) DO j = 1, nColumns ! target column numerator = 0.0D0 denominator = 0.0D0 DO ii = lowest_ii, highest_ii ! potential source rows numerator_vector = 0.0 ! whole REAL vector denominator_vector = 0.0 ! whole REAL vector !dir$ loop count min(16) DO jj = 1, nColumns ! potential source columns radians_away_squared = (uvecs(1, ii, jj) - uvecs(1, i, j))**2 + & & (uvecs(2, ii, jj) - uvecs(2, i, j))**2 + & & (uvecs(3, ii, jj) - uvecs(3, i, j))**2 IF (radians_away_squared <= horizon_squared) THEN Gaussian = EXP(-0.5 * radians_away_squared / radian_scale_squared) G_dA = Gaussian * relative_area(ii) numerator_vector(jj) = G_dA * grid(ii, jj) denominator_vector(jj) = G_dA END IF ! local END DO ! jj = 1, nColumns; potential source columns DO jj = 1, nColumns numerator = numerator + numerator_vector(jj) denominator = denominator + denominator_vector(jj) END DO END DO ! ii = 1, nRows; potential source rows IF (denominator > 0.0D0) THEN smoothed_grid(i, j) = numerator / denominator ELSE smoothed_grid(i, j) = 0.0 END IF counter = counter + 1 percent_done = counter / small_part ! rounded down to INT IF (percent_done > old_percent_done) THEN old_percent_done = percent_done WRITE (*, "('+ ',I3,'% completed')") percent_done END IF END DO ! target column j = 1, nColumns END DO ! target row i = 1, nRows END SUBROUTINE Smooth_REAL_on_Sphere END PROGRAM GEAR1_for_CSEP