PROGRAM Long_Term_Seismicity_v3 !========================================================= !This program exists in more than one version. !For version number, date, and relevant text see !the first WRITE statement block below, at "!VersionHere". !========================================================= ! Predicts long-term-average shallow seismicity of Earth, ! without any time-dependence. Provides output in the choice of: ! !(1) Peter Bird's .grd format [used by NeoKineMap & FiniteMap]: ! regular (longitude, latitude) grids of numbers ! which are seismicity rates (above a specified threshold magnitude) ! in SI units of [events]/square meter/second. ! Because these numbers are very small, they are expressed as ! double-precision numbers to avoid underflow. ! However, it should be understood that the accuracy of these ! forecasts is probably no better than 1 significant digit! ! ! Suggestion: Use my programs NeoKineMap and/or FiniteMap to ! display the output as a colored bitmap, with overlays. ! NeoKineMap includes an option to display the log ! of seismicity rate, which is usually superior ! to a linear display. ! !(2) RELM (Regional Earthquake Likelihood Models) format, for ! SCEC (Southern California Earthquake Center) standardized ! testing and comparison of models. Events expected are listed ! per grid rectangle per 5-years, within successive magnitude limits ! ranging from 4.95-5.05 (first bin) to 8.85-8.95 (penultimate bin) ! and lastly 8.95-10.00 (ultimate bin). The details of this format ! were specified by Danijell Schorlemmer in emails in 2005. ! Note that these forecasts INCLUDE all aftershocks, so they should ! never be quantitatively tested against "declustered" catalogs! ! by Peter Bird, ! Department of Earth and Space Sciences, ! University of California, Los Angeles, CA 90095-1567 ! [For program date(s), see first block of WRITE statements, ! at location !VersionHere below.] ! e-mail: pbird@ess.ucla.edu ! URL: http://element.ess.ucla.edu ! Financial support (at various times) has been provided by: ! *University of California, Los Angeles; ! *National Earthquake Hazards Reduction Program ! of the U.S. Geological Survey; ! *National Science Foundation; ! *Southern California Earthquake Center. ! Scientific collaboration has been provided (at various times) by: ! *David D. Jackson, UCLA ! *Yan Y. Kagan, UCLA ! *Zhen Liu, UCLA ! However, all computational methods and numerical results ! represent solely the opinions of the author, and do not represent ! the offical position of the United States government or the ! government of the state of California. USE Sphere ! Sphere.f90 = module of spherical-geometry tools by P. Bird, UCLA USE Numerical_Libraries ! IMSL Library, used for ANORDF = distribution function of a standard normal (Gaussian) variable; IMPLICIT NONE CHARACTER*1 :: c1, class_continuity_c1, star_c1, step_continuity_c1, tab CHARACTER*2 :: issue_date_day, issue_date_month, start_date_day, start_date_month CHARACTER*3 :: class CHARACTER*4 :: issue_date_year, start_date_year CHARACTER*5 :: c5, c5a CHARACTER*6 :: c6 CHARACTER*6, DIMENSION(12) :: T5_column_header CHARACTER*6, DIMENSION(:), ALLOCATABLE :: segment_c6 CHARACTER*80 :: global_feg_pathfile, global_v_pathfile, line, output_grd_pathfile, & & PB2002_steps_pathfile, template_pathfile CHARACTER*80, DIMENSION(:), ALLOCATABLE :: e_nko_pathfile, h_nko_pathfile, orogen_feg_pathfile, & & orogen_outline_dig_pathfile CHARACTER*3, DIMENSION(0:7) :: class_c3 ! = INT, plus boundary classes 1-7 DOUBLE PRECISION :: angle_sum, argument, average, & & B, & & CMT_threshold_moment, corner_moment, & & desired_threshold_moment_Nm, dI, difference, dot_a, dot_b, dx, dy, & & expectation, & & G_function, & & moment_rate_of_step_SI, & & point_value, & & seismicity_at_CMT_threshold, & & tGR_beta, & & weight DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: EQs_per_s ! (subscript is threshold_index = 1-41 IFF format_selector == 2 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: EQs_per_s_per_m2 ! (subscript is threshold_index = 1-41 IFF format_selector == 2 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: integral ! (subscript is threshold_index = 1-41 IFF format_selector == 2 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: integral_converted ! (subscript is threshold_index = 1-41 IFF format_selector == 2 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: single_column, single_row DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: seismicity ! (rows {N->S}, columns {W->E}, threshold {low->high}) DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: PB2002_seismicity ! parallel structure to that of seismicity INTEGER, PARAMETER :: subsamplings = 4 ! Number of subdivisions of each grid step, in each horizontal direction. !(So, the number of integration points for each grid point is subsamplings**2.) ! Subsamplings must be positive. 1 gives very fast, but *inaccurate*, results. ! I suggest using ~4 when your grid spacing is 0.1 degree !(allowing you to run global grids in 1-2 GB of memory, on a 32-bit computer), ! and higher subsamplings when the grid spacing is larger than 0.1 degree. INTEGER :: age_Ma, azimuth_int, azimuth_integer, & & column, columns, crossings_to_W, & & element, elevation, & & format_selector, & & i, iele, ii, ios, & & j, jj, j_East, j_West, & & k, k_class_index, & & l_, larger_depth_integer, line_number, & & m, m_time, max_in_outline, MB, & & n, nr1, nr2, n_step, number, number_of_zeros, numnod, numel, nfl, & & orogen, orogens, other_segment, & & row, rows, & & segment, segments, sequence_number, smaller_depth_integer, & & T5_column, threshold_index, timesteps_needed, top_threshold_index, & & velocity_azimuth_integer, velocity_azimuth_int, & & which_orogen_at_beginning, which_orogen_at_end, wrap_type INTEGER(KIND = 2), DIMENSION(-1080:1080, 0:4319) :: ETOPO5 ! Note that row/column 0 means lat/lon of 0 INTEGER(KIND = 2), DIMENSION( -720:900, 0:3600) :: age_1p5 ! Note that row/column 0 means lat/lon of 0 INTEGER, DIMENSION(:), ALLOCATABLE :: n_in_outline ! (orogen = 1, ..., orogens) INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor, nodes LOGICAL :: cold_start, creeping, decision_in, & & include_in_integral, inside, maybe, oceanic, overlap, & & parallel, problem, proximal, repeat_sweep, success LOGICAL*1, DIMENSION(:), ALLOCATABLE :: segment_creeping, segment_used LOGICAL*1, DIMENSION(:,:), ALLOCATABLE :: continental, template_mask REAL, PARAMETER :: Earth_radius_m = 6371000.0 !non-SUB boundaries (symmetrical): REAL, PARAMETER :: CCB_halfwidth_km = 189.0 REAL, PARAMETER :: CTF_halfwidth_km = 257.0 REAL, PARAMETER :: CRB_halfwidth_km = 154.0 REAL, PARAMETER :: OSR_halfwidth_km = 132.0 REAL, PARAMETER :: OTF_halfwidth_km = 128.0 REAL, PARAMETER :: OCB_halfwidth_km = 186.0 ! these values from Table 1 of Bird & Kagan [2004] !SUB boundaries (asymmetrical): based on EQ_classification_II.f90 of Bird & Kagan [2004; BSSA] REAL, PARAMETER :: SUB_landward_halfwidth_km = 220.0 ! landward width of subduction lanes, in km from trench REAL, PARAMETER :: SUB_seaward_halfwidth_km = 135.0 ! seaward half-width of subduction lanes, in km from trench REAL, PARAMETER :: SUB_peakSeismicity_km = 55.0 ! location of seismicity peak, relative to trench !and one new SUB parameter added for v2 of LTS, following Bird et al. [2009, BSSA]: REAL, PARAMETER :: SUB_bending_fraction = 0.38 ! relative seismicity level when main interplate thrust is greased REAL, PARAMETER :: alongStep_sigma_km = 32.5 ! so combined EQ mislocation and step-mislocation unlikely to exceed twice this value REAL :: angle, azimuth_radians, & & basis_length_km, best_separation, & & center_latitude, center_longitude, & & CMT_duration_s, COS_latitude_radians, cross_trace_km, coupledThickness_m, & & d_lat, d_lon, d_S_radians, d_E_radians, d_t, d0_m, dd_lat, dd_lon, delta_degrees, delta_lon_from_grid_center, dextral_mmpa, & & diffusion_scale, distance, distance_km, divergence_mmpa, dot_test, ds, dx_meters, dy_meters, & & edge_longitude, equat, err, e1_rate, e1h, e2_rate, e2h, e3_rate, & & factor, fault_dip_degrees, fault_dip_radians, & & fraction, fraction_inside, fraction_toward_uvec1, fraction_toward_uvec2, & & from_trench_azimuth_radians, & & high_template_latitude, high_template_longitude, highest_latitude, highest_longitude, highest_magnitude, & & intraplate_rate, & & lane_width_km, largest_ei_persec, lat, lat_max, lat_min, lat_offset, lat1, lat2, latitude, & & length_km, length_kilometers, & & lithosphere_thickness_km, & & lon, lon_max, lon_min, lon_offset, lon1, lon2, longitude, & & low_template_latitude, low_template_longitude, lowest_latitude, lowest_longitude, lowest_magnitude, & & M_persec_per_m2, M_persec_per_m3, middle, middle_D, midpoint_compass, minimum_arc, model_heave_rate_mmpa, & & offset_km, offside_radians, old_latitude, old_longitude, & & parallelogram_area_m2, & & radius_km, & & s_km, s_per_year, s_radians, s1, s2, s3, safe_timestep, scalar_product, scales, separation, slip_rate_mmpa, & & SUB_landward_radians, SUB_landward_sigma_km, SUB_seaward_radians, SUB_seaward_sigma_km, & & subduction_azimuth_1, subduction_azimuth_2, & & t_lat, t_lon, t1, t2, temp_integral, theta_, threshold_magnitude, to_uvec1_radians, to_uvec2_radians, total_time, & & u1phi, u1theta, u2phi, u2theta, & & v_oblique_mmpa, v_oblique_mps, velocity_mmpa, vs_mmpa, & & X, & & Y, Y_factor_of_step REAL, DIMENSION(3) :: begin_uvec, & & c1_uvec, c2_uvec, c3_uvec, c4_uvec, center_uvec, & & end_uvec, eps_dot, & & first_a_uvec, first_b_uvec, & & inward_uvec, & & last_a_uvec, last_b_uvec, & & omega_uvec, & & point1_uvec, point2_uvec, pole_uvec, pole_a_uvec, pole_b_uvec, & & pole_uvec1, pole_uvec2, pole_uvec3, pole_uvec4, & & segment_vec, & & trench_uvec, tvec, & & uvec, uvec1, uvec2, uvec3, uvec4 REAL, DIMENSION(4) :: dot_side ! sine of angle between test point and sides of a subduction lane; positive inward REAL, DIMENSION(6) :: h_km ! corresponds to h_k in Bird & Kagan [2004]; no entry for SUB because that is asymmetrical REAL, DIMENSION(12) :: 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 REAL, DIMENSION(12) :: T5_supplement_seismicity_fraction_based_on_z ! == 1.0 if format_selector == 1; otherwise less (user input). REAL, DIMENSION(3,5) :: lane_uvecs REAL, DIMENSION(3,7) :: Gauss_point REAL, DIMENSION(3,2,2) :: G ! nodal functions of spherical triangle element REAL, DIMENSION(3,2,2,2):: dG ! derivitives of nodal functions G REAL, DIMENSION(2, 130) :: GCNocb ! digitized ocean-continent boundary offshore Gorda-Cal-Nev orogen REAL, DIMENSION(:), ALLOCATABLE :: a_ REAL, DIMENSION(:), ALLOCATABLE :: RELM_threshold_magnitude, RELM_threshold_moment_Nm REAL, DIMENSION(:), ALLOCATABLE :: segment_azimuth_radians, segment_slip_rate_mmpa, segment_heave_rate_mmpa REAL, DIMENSION(:), ALLOCATABLE :: vw REAL, DIMENSION(:,:), ALLOCATABLE :: center, element_principal_strainrates, node_uvec REAL, DIMENSION(:,:), ALLOCATABLE :: diffusivity ! (rows {N->S}, columns {W->E}) REAL, DIMENSION(:,:,:), ALLOCATABLE :: outline_uvecs ! (1 ~ 3, 1 ~ n_in_outline(orogen), 1 ~ orogens) REAL, DIMENSION(:,:,:), ALLOCATABLE :: segment_uvecs, segment_outer_uvecs REAL, DIMENSION(:,:,:,:,:), ALLOCATABLE :: subsample_uvec ! about 18.5 MB for 360 x 180 grid with subsamplings = 5, ! or 1.16 GB for global 0.10 degree grid, subsamplings = 4. DATA Gauss_point / & & 0.3333333333333333E0, 0.3333333333333333E0, 0.3333333333333333E0, & & 0.0597158733333333E0, 0.4701420633333333E0, 0.4701420633333333E0, & & 0.4701420633333333E0, 0.0597158733333333E0, 0.4701420633333333E0, & & 0.4701420633333333E0, 0.4701420633333333E0, 0.0597158733333333E0, & & 0.7974269866666667E0, 0.1012865066666667E0, 0.1012865066666667E0, & & 0.1012865066666667E0, 0.7974269866666667E0, 0.1012865066666667E0, & & 0.1012865066666667E0, 0.1012865066666667E0, 0.7974269866666667E0/ tab = CHAR(9) ! character "HT" for horizontal tab radius_km = Earth_radius_m / 1000.0 class_c3(0) = "INT" class_c3(1) = "CCB" class_c3(2) = "CTF" class_c3(3) = "CRB" class_c3(4) = "OSR" class_c3(5) = "OTF" class_c3(6) = "OCB" class_c3(7) = "SUB" h_km(1) = CCB_halfwidth_km h_km(2) = CTF_halfwidth_km h_km(3) = CRB_halfwidth_km h_km(4) = OSR_halfwidth_km h_km(5) = OTF_halfwidth_km h_km(6) = OCB_halfwidth_km ! no entry for SUB because it is asymmetrical SUB_landward_radians = SUB_landward_halfwidth_km / radius_km ! landward width of subduction lanes, in radians, from trench SUB_seaward_radians = SUB_seaward_halfwidth_km / radius_km ! seaward half-width of subduction lanes, in radians, from trench SUB_landward_sigma_km = 0.5 * (SUB_landward_halfwidth_km - SUB_peakSeismicity_km) ! sigmas for Gaussian approximation of B SUB_seaward_sigma_km = 0.5 * (SUB_seaward_halfwidth_km + SUB_peakSeismicity_km) ! **************************************************************************************** !GPBT5 ! Empirical coefficients from selected rows of Table 5 of Bird & Kagan [2004, BSSA] ! {qualification added for Version 2: ! with two changes: ! *distinction between slow and fast CCBs (old column 3 becomes new columns 3/ 4) ! *distinction between slow and fast SUBS (old column 10 becomes new columns 11/12) ! based on newer results in Bird et al. [2009, BSSA].} ! T5_column = 1 2 3 4 5 6 7 8 9 10 11 12 ! CRB CTF CCB CCB OSR OSR OTF OTF OTF OCB SUB SUB ! slow fast normal other slow medium fast slow fast T5_column_header = (/"CRB ","CTF ","CCBslo","CCBfas","OSRnor","OSRoth","OTFslo","OTFmed","OTFfas","OCB ","SUBslo","SUBfas"/) T5_CMT_pure_events = (/ 285.9, 198.5, 78.3, 181.1, 424.3, 77.0, 398.0, 406.9, 376.6, 117.7, 399.5, 1653.3/) T5_CMT_threshold_moment = (/ 1.13E17, 3.5E17, 3.5E17, 3.5E17, 1.13E17, 1.13E17, 2.0E17, 2.0E17, 2.0E17, 3.5E17, 3.5E17, 3.5E17/) T5_beta = (/ 0.65, 0.65, 0.62, 0.62, 0.92, 0.82, 0.64, 0.65, 0.73, 0.53, 0.64, 0.64/) T5_corner_magnitude = (/ 7.64, 8.01, 8.46, 8.46, 5.86, 7.39, 8.14, 6.55, 6.63, 8.04, 9.58, 9.58/) T5_tGR_model_moment_rate = (/ 1.67E12, 3.8E12, 3.20E12, 7.40E12, 6.7E11, 1.9E11, 6.7E12, 9.4E11, 9.0E11, 4.6E12, 5.55E13, 2.29E14/) T5_length_km = (/ 18126., 19375., 10461., 2054., 61807., 61807., 27220., 10351., 6331., 13236., 21491., 17499./) T5_mean_velocity_mmpa = (/ 18.95, 21.54, 10.8, 55.8, 46.40, 7.58, 20.68, 57.53, 97.11, 19.22, 37.00, 91.54/) T5_assumed_dip_degrees = (/ 55., 73., 20., 20., 55., 55., 73., 73., 73., 20., 14., 14./) T5_assumed_mu_GPa = (/ 27.7, 27.7, 27.7, 27.7, 25.7, 25.7, 25.7, 25.7, 25.7, 49., 49., 49./) T5_lineIntegral_Nps = (/ 5.5E8, 4.4E8, 2.9E8, 3.0E8, 5.0E9, 4.7E8, 5.2E8, 5.3E8, 5.5E8, 1.2E9, 5.24E9, 1.06E10/) T5_coupledThickness_km = (/ 3.0, 8.6, 10.9, 24.9, 0.13, 0.40, 13., 1.8, 1.6, 3.8, 10.6, 21.7/) T5_assumed_lithosphere_km =(/ 6., 12., 13., 24.9, 8., 8., 14., 14., 14., 14., 26., 26./) T5_coupling = (/ 0.50, 0.72, 0.84, 1.00, 0.016, 0.05, 0.93, 0.13, 0.11, 0.27, 0.41, 0.84/) ! and, some readily-derived quantitites: s_per_year = 365.25 * 24 * 60 * 60 CMT_duration_s = 25.7474 * s_per_year ! 1977.01.01 through 2002.09.30 <== Do NOT update! Must match Bird & Kagan [2004] calibration study. DO i = 1, 12 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 ! **************************************************************************************** WRITE (*, "(/' PROGRAM Long_Term_Seismicity_v3')") WRITE (*, "(/' Predicts long-term-average shallow seismicity of Earth,')") WRITE (*, "( ' without any time-dependence. This is a Poissonian forecast,')") WRITE (*, "( ' and it includes aftershocks as well as mainshocks, without distinction.')") !VersionHere WRITE (*, "(/' by Peter Bird, UCLA; version 3 of 29 April 2009.')") !----------------------------------------------------------------------------------- !Version 1 (2006.10.26) was based on assumptions, equations, and rules presented in: ! *Bird, P., and Z. Liu [2007] Seismic hazard inferred from tectonics: ! California, Seismol. Res. Lett., 78(1), 37-48. !and it also makes use of definitions and research findings in: ! *Bird, P. [2003] An updated digital model of plate boundaries, ! Geochemistry Geophysics Geosystems, 4(3), 1027, doi:10.1029/2001GC000252. ! *Bird, P., and Y. Y. Kagan [2004] Plate-tectonic analysis of shallow seismicity: ! Apparent boundary width, beta, corner magnitude, coupled lithosphere ! thickness, and coupling in seven tectonic settings, ! Bull. Seismol. Soc. Am., 94(6), 2380-2399, plus electronic supplement. ! (N.B. The most important results are in Table 5.) !Version 1 was used in the calculations that were reported in: ! *Bird, P., and Z. Liu [2007] Seismic hazard inferred from tectonics: ! California, Seismol. Res. Lett., 78(1), 37-48. ! *Liu, Z., & P. Bird [2008] Kinematic modelling of neotectonics in the Persia-Tibet-Burma ! orogen, Geophys. J. Int., 172(2), 779-797; ! doi:10.1111/j.1365-246X.2007.03640.x, + 3 digital file appendices. !----------------------------------------------------------------------------------- !Version 2 (2009.04.08) is almost the same, but incorporates 3 discoveries of: ! *Bird, P., Y. Y. Kagan, D. D. Jackson, F. P. Schoenberg, & M. J. Werner [2009], ! Linear and nonlinear relations between relative plate velocity and seismicity, ! Bull. Seismol. Soc. Am.: ! (1) Coupled thickness of Continental Convergent Boundaries is not uniform, but ! increases by a factor-of-2 above a critical relative plate velocity of ~24 mm/a. ! Instead of a global coupled thickness of ~18 km for CCBs, there is a coupled ! thickness of ~10.9 km below ~24 mm/a, but a coupled thickness of ~24.9 km ! at higher velocities. ! (2) Coupled thickness of SUBduction zones is not uniform, but increases by about ! a factor-of-2 above a critical relative plate velocity of about ~66 mm/a. ! Instead of a global coupled thickness of ~18 km for SUBs, there is a coupled ! thickness of ~10.6 km below ~66 mm/a, but a coupled thickness of ~21.7 km ! at higher velocities. ! (3) Subduction zones which have aseismic creep on the main interplate thrust fault ! (e.g., when lubricated by evaporites) should NOT be totally aseismic, but should ! still produce earthquakes due to plate-bending, and perhaps due to deformation ! in the upper plate (e.g., in partitioned oblique subduction). ! Such SUBs apparently have earthquake productivity which is ~38% of normal. ! This factor of "SUB_bending_fraction = 0.38" is now applied to F0001S ! faults (SUB faults) in any NeoKinema model of an orogen which have been ! designated as creeping in the f_{token}.nko file. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Version 2 also uses a different method to make the "seams" (outside orogens) ! between plate-interior seismicity from triangular continum F-Es of an EarthN grid ! and the Gaussian-smoothed seismicity of simple plate boundary steps from PB2002: ! Now I use the GREATER of these two contributions, instead of the SUM. ! (Either method avoids a discontinuity, but the old v1 method could potentially ! double-count some plate-boundary-adjacent deformation and seismicity.) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Finally, Version 2 allows the PB2002_steps.file to have a different name & length. !----------------------------------------------------------------------------------- !Version 3 adds diffusive smoothing of seismicity around orogen boundaries, ! both to eliminate unattractive discontinuities, and to acknowledge ! that my original choice of each orogen boundary was somewhat arbitrary ! and subject to uncertainty of its own. This smoothing is performed ! prior to the addition of narrow red ribbons representing seismicity of ! faults within the orogen, so their contributions are not smoothed. ! This smoothing conserves total forecast seismicity. ! (To be sure of this, the 4 quadrilateral bounds of the .GRD area are ! treated as impermeable to diffusion of seismicity.) ! Smoothing is performed by explicit F-D method, pretending that the ! the longitude/latitude grid is a Cartesian x/y grid. ! Because this method will no work around the poles, no smoothing is ! performed at latitudes higher than 75 degrees away from the equator. !----------------------------------------------------------------------------------- 1 WRITE (*, "(/' CHOICE OF OUTPUT FORMATS:')") WRITE (*, "( ' ')") WRITE (*, "( ' .grd format (code 1) RELM format (code 2)')") WRITE (*, "( ' ==================== ==========================================')") WRITE (*, "( ' Time window: 1 second 5 years')") WRITE (*, "( ' Area footprint: 1 square meter rectangle based on grid spacing')") WRITE (*, "( ' Depth range: 0 to 70 km user-controlled (e.g., 0 to 30 km)')") WRITE (*, "( ' Magnitudes: all ABOVE threshold binned: 4.95-5.05, 5.05-5.15, ..., 8.95-10')") WRITE (*, "( ' Format: Peter Bird''s .grd Danijell Scholemmer''s RELM')") WRITE (*, "(/' Enter (integer) code to select desired format: '\)") READ (*, *) format_selector IF ((format_selector < 1).OR.(format_selector > 2)) THEN WRITE (*, "(' ERROR: Try again, using integer 1 or 2:')") CALL Pause() GO TO 1 END IF IF (format_selector == 1) THEN WRITE (*, "(' Enter threshold (minimum) magnitude: '\)") READ (*, *) threshold_magnitude desired_threshold_moment_Nm = Moment(threshold_magnitude) WRITE (*, "(' Threshold (minimum) scalar moment =',ES11.3,' N m')") desired_threshold_moment_Nm ELSE IF (format_selector == 2) THEN !arrays of threshold magnitudes and moments will be created below. END IF WRITE (*, "(/' Describe region of desired output file:')") IF (format_selector == 2) WRITE (*, "(' RELM currently requires -125.35 for lowest longitude of a cell center.')") WRITE (*, "( ' Westernmost or left longitude limit (in degrees E of Greenwich): '\)") READ (*, *) lon_min IF (format_selector == 2) WRITE (*, "(' RELM currently requires -113.15 for highest longitude of a cell center.')") WRITE (*, "( ' Easternmost or right longitude limit (in degrees E of Greenwich): '\)") READ (*, *) lon_max IF (lon_max <= lon_min) THEN lon_max = lon_max + 360.0 WRITE (*, "(' NOTE: lon_max expressed as ',F8.3,' so as to be greater than lon_min.')") lon_max CALL Pause() END IF IF (format_selector == 2) WRITE (*, "(' RELM currently requires 31.55 for lowest latitude of a cell center.')") WRITE (*, "( ' Southernmost or bottom latitude limit (in degrees N of equator): '\)") READ (*, *) lat_min IF (format_selector == 2) WRITE (*, "(' RELM currently requires 42.95 for highest latitude of a cell center.')") WRITE (*, "( ' Northernmost or top latitude limit (in degrees N of equator): '\)") READ (*, *) lat_max IF (lat_max < lat_min) THEN WRITE (*, "( ' ERROR: Lat_max must be >= lat_min.')") CALL Pause() STOP END IF IF (format_selector == 1) THEN ! report all seismicity down to 70 km depth (from sea level), like Bird & Kagan [2004]. smaller_depth_integer = 0 ! seismicity forecast covers depth range starting at 0 km. larger_depth_integer = 70 ! seismicity forecast covers depth range down to 70 km (below sea level). T5_supplement_seismicity_fraction_based_on_z = 1.00 ! whole array ELSE IF (format_selector == 2) THEN ! request guidance on fractions of seismicity to report 2 WRITE (*, "(/' Specify smaller of the depth limits (as integer, in km) for this forecast): '\)") READ (*, *) smaller_depth_integer WRITE (*, "( ' Specify larger of the depth limits (as integer, in km) for this forecast): '\)") READ (*, *) larger_depth_integer WRITE (*, "(/' User must specify seismicity fractions based on this depth range of forecast:')") DO i = 1, 12 WRITE (*, "(' column ',I2,': ',A6,' seismicity fraction: '\)") i, T5_column_header(i) READ (*, *) T5_supplement_seismicity_fraction_based_on_z(i) IF ((T5_supplement_seismicity_fraction_based_on_z(i) < 0.0).OR. & & (T5_supplement_seismicity_fraction_based_on_z(i) > 1.0)) THEN WRITE (*, "(' ERROR! Each fraction must be in range 0.0 ~ 1.0. Start over!')") CALL Pause() GOTO 2 END IF END DO END IF IF (format_selector == 2) WRITE (*, "(/' RELM currently (2005.08) requires 0.1-degree bins.')") WRITE (*, "( ' Enter grid spacing in degrees: '\)") READ (*, *) delta_degrees IF (delta_degrees <= 0.0) THEN WRITE (*, "( ' ERROR: delta_degrees must be positive.')") CALL Pause() STOP END IF d_lat = delta_degrees d_lon = delta_degrees rows = 1 + NINT((lat_max - lat_min) / d_lat) columns = 1 + NINT((lon_max - lon_min) / d_lon) WRITE (*, "(' Output grid will have ',I6,' rows and ',I6,' columns.')") rows, columns IF (format_selector == 2) THEN ! read template file to set template_mask(rows, columns) = .TRUE. WRITE (*, "(/' Enter [Drive:][\path\]filename of template file from RELM (to define test area): ')") READ (*, "(A)") template_pathfile END IF !subsampling parameters (in neighborhood of each grid point): IF (subsamplings < 1) THEN WRITE (*, "(' ERROR: INTEGER, PARAMETER :: subsamplings must be positive.')") CALL Pause() STOP END IF dd_lat = d_lat / subsamplings dd_lon = d_lon / subsamplings middle = (1 + subsamplings) / 2 lat_offset = -dd_lat * middle lon_offset = -dd_lon * middle weight = 1.0D0 / (subsamplings**2) WRITE (*, "(/' Enter [Drive:][\path\]filename for output [.grd? or .dat?] file: ')") READ (*, *) output_grd_pathfile OPEN (UNIT = 1, FILE = output_grd_pathfile, STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Requested output file named' & & /' ',A, & & /' could not be created. Perhaps it already exists. IOSTAT = ',I6)") & & TRIM(output_grd_pathfile), ios CALL Pause() STOP END IF IF (format_selector == 1) THEN WRITE (1, "(3F11.5,' = lon_min, d_lon, lon_max Gridded long-term seismicity above magnitude ',F6.3,' in events/square-meter/s')") lon_min, d_lon, lon_max, threshold_magnitude WRITE (1, "(3F11.5,' = lat_min, d_lat, lat_max from program Long_Term_Seismicity_v3, by P. Bird, UCLA ')") lat_min, d_lat, lat_max ELSE IF (format_selector == 2) THEN WRITE (*, "(/' Next you will be prompted for header information, including:')") WRITE (*, "( ' modelname')") WRITE (*, "( ' version')") WRITE (*, "( ' author')") WRITE (*, "( ' (NOTE that the first 3 lines are free-format text.)')") WRITE (*, "( ' issue_date')") WRITE (*, "( ' forecast_start_date')") WRITE (*, "( ' forecast_duration=5,years')") WRITE (*, "( ' (NOTE that forecast_duration does not need to be entered.')") WRITE (*, "(' Enter model name (e.g., Long_Term_Seismicity_v3 based on NeoKinema):')") READ (*, "(A)") line WRITE (1, "(A)") "modelname=" // TRIM(line) WRITE (*, "(' Enter version number (e.g., 1.0?) on following line:')") READ (*, "(A)") line WRITE (1, "(A)") "version=" // TRIM(line) WRITE (*, "(' Enter author''s name(s) on following line:')") READ (*, "(A)") line WRITE (1, "(A)") "author=" // TRIM(line) WRITE (*, "(' Enter year of issue date (e.g., 2005): '\)") READ (*, *) issue_date_year WRITE (*, "(' Enter month of issue date (from 1-12): '\)") READ (*, *) issue_date_month WRITE (*, "(' Enter day of issue date (from 1-31): '\)") READ (*, *) issue_date_day WRITE (1, "(A)") "issue_date=" // & & TRIM(ADJUSTL(issue_date_year)) // ',' // & & TRIM(ADJUSTL(issue_date_month)) // ',' // & & TRIM(ADJUSTL(issue_date_day)) // ",0,0,0,0" WRITE (*, "(' Enter year of start date (e.g., 2006): '\)") READ (*, *) start_date_year WRITE (*, "(' Enter month of start date (from 1-12; e.g., 1): '\)") READ (*, *) start_date_month WRITE (*, "(' Enter day of start date (from 1-31; e.g., 1): '\)") READ (*, *) start_date_day WRITE (1, "(A)") "start_date=" // & & TRIM(ADJUSTL(start_date_year)) // ',' // & & TRIM(ADJUSTL(start_date_month)) // ',' // & & TRIM(ADJUSTL(start_date_day)) // ",0,0,0,0" WRITE (1, "('forecast_duration=5,years')") WRITE (1, "('begin_forecast')") END IF WRITE (*, "(/' Enter [Drive:][\path\]filename for the global .feg file' & & /' that was used with Shells to compute intraplate strain-rates: ')") READ (*, *) global_feg_pathfile WRITE (*, "(/' Enter [Drive:][\path\]filename for the global v_.out file' & & /' that came from Shells, using the .feg file above: ')") READ (*, *) global_v_pathfile WRITE (*, "(/' To decribe the global rigid-plate model PB2002 of Bird [2003], you will' & &/' also need a file like PB2002_steps.dat (from an on-line archive).' & &/' Enter [Drive:][\path\]filename for this steps file: ')") READ (*, *) PB2002_steps_pathfile WRITE (*, "(/' How many orogens have you modeled with NeoKinema (0, 1, 2, ... 13)? '\)") READ (*, *) orogens IF (orogens < 0) THEN WRITE (*, "(' ERROR: Orogen count may not be negative.')") CALL Pause() STOP END IF IF (orogens > 0) THEN ALLOCATE ( orogen_feg_pathfile(orogens) ) ALLOCATE ( orogen_outline_dig_pathfile(orogens) ) ALLOCATE ( e_nko_pathfile(orogens) ) ALLOCATE ( h_nko_pathfile(orogens) ) DO orogen = 1, orogens IF (orogen <= 9) THEN WRITE (*, "(/' Enter [Drive:][\path\]filename of _.feg file for orogen #',I1,':')") orogen READ (*, *) orogen_feg_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of outline.dig file for orogen #',I1,':')") orogen READ (*, *) orogen_outline_dig_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of e_.nko file for orogen #',I1,':')") orogen READ (*, *) e_nko_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of h_.nko file for orogen #',I1,':')") orogen READ (*, *) h_nko_pathfile(orogen) ELSE ! orogen >= 10 (two-digit): WRITE (*, "(/' Enter [Drive:][\path\]filename of _.feg file for orogen #',I2,':')") orogen READ (*, *) orogen_feg_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of outline.dig file for orogen #',I2,':')") orogen READ (*, *) orogen_outline_dig_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of e_.nko file for orogen #',I2,':')") orogen READ (*, *) e_nko_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of h_.nko file for orogen #',I2,':')") orogen READ (*, *) h_nko_pathfile(orogen) END IF END DO ! i = 1, orogens END IF ! orogens > 0 WRITE (*, "(/' *** INPUT FILES REQUIRED: ***')") WRITE (*, "( ' (Empirical coefficients from Bird & Kagan [2004, BSSA] are built-in.)')") WRITE (*, "( ' ---------------------------------------------------------------------')") IF (format_selector == 2) WRITE (*, "(' ',A)") TRIM(template_pathfile) WRITE (*, "(' ETOPO5.grd (topography of Earth with 5 min. grid: 4320 columns, 2161 rows)')") WRITE (*, "(' age_1p5.grd (seafloor ages, Mueller et al., 1997: 3601 columns, 1621 rows')") WRITE (*, "(' ',A)") TRIM(global_feg_pathfile) WRITE (*, "(' ',A)") TRIM(global_v_pathfile) WRITE (*, "(' ',A)") TRIM(PB2002_steps_pathfile) IF (orogens > 0) THEN DO orogen = 1, orogens WRITE (*, "(' ',A)") TRIM(orogen_feg_pathfile(orogen)) WRITE (*, "(' ',A)") TRIM(orogen_outline_dig_pathfile(orogen)) WRITE (*, "(' ',A)") TRIM(e_nko_pathfile(orogen)) WRITE (*, "(' ',A)") TRIM(h_nko_pathfile(orogen)) IF ((orogens > 1).AND.(orogen < orogens)) CALL Pause() END DO END IF WRITE (*, "( ' ---------------------------------------------------------------------')") WRITE (*, "( ' If any files must be moved to the current directory/folder, DO SO NOW!')") WRITE (*, "( ' (After this point, program execution may be left unattended.)')") CALL Pause() IF (format_selector == 2) CALL Test_Opening_of_Existing_File(TRIM(template_pathfile)) CALL Test_Opening_of_Existing_File("ETOPO5.grd") CALL Test_Opening_of_Existing_File("age_1p5.grd") CALL Test_Opening_of_Existing_File(TRIM(global_feg_pathfile)) CALL Test_Opening_of_Existing_File(TRIM(global_v_pathfile)) CALL Test_Opening_of_Existing_File(TRIM(PB2002_steps_pathfile)) IF (orogens > 0) THEN DO orogen = 1, orogens CALL Test_Opening_of_Existing_File(TRIM(orogen_feg_pathfile(orogen))) CALL Test_Opening_of_Existing_File(TRIM(orogen_outline_dig_pathfile(orogen))) CALL Test_Opening_of_Existing_File(TRIM(e_nko_pathfile(orogen))) CALL Test_Opening_of_Existing_File(TRIM(h_nko_pathfile(orogen))) END DO END IF IF (format_selector == 1) THEN top_threshold_index = 1 ALLOCATE ( EQs_per_s(1) ) ALLOCATE ( EQs_per_s_per_m2(1) ) ALLOCATE ( integral(1) ) ALLOCATE ( integral_converted(1) ) ALLOCATE ( seismicity(rows, columns, 1) ) ! NOTE: rows are numbered from 1 in the N to "rows" in the S; ! columns are numbered from 1 in the W to "columns" in the E. ALLOCATE ( PB2002_seismicity(rows, columns, 1) ) ELSE IF (format_selector == 2) THEN top_threshold_index = 41 ALLOCATE ( EQs_per_s(41) ) ALLOCATE ( EQs_per_s_per_m2(41) ) ALLOCATE ( integral(41) ) ALLOCATE ( integral_converted(41) ) ALLOCATE ( seismicity(rows, columns, 41) ) ! NOTE: rows are numbered from 1 in the N to "rows" in the S; ! columns are numbered from 1 in the W to "columns" in the E; ! thresholds increase from 4.95, 5.05, ..., 8.95 ALLOCATE ( PB2002_seismicity(rows, columns, 41) ) ALLOCATE ( template_mask(rows, columns) ) template_mask = .FALSE. ! whole array, until changed (below) be reading a template cell ALLOCATE ( RELM_threshold_magnitude(41) ) ALLOCATE ( RELM_threshold_moment_Nm(41) ) DO i = 1, 41 RELM_threshold_magnitude(i) = 4.85 + i * 0.10 RELM_threshold_moment_Nm(i) = Moment(RELM_threshold_magnitude(i)) END DO END IF seismicity = 0.0D0 ! whole area; all thresholds {if more than 1}; to be incremented below. PB2002_seismicity = 0.0D0 ! whole area; all thresholds. This will be used to accumulate the sum of ! seismicity due to PB2002 plate boundaries (only) with their Gaussian smoothing ! footprints. Then, existing seismicity() [from EarthN continuum strain, outside orogens] ! will be replaced with MAX(seismicity(), PB2002()) for a smooth seam. ! Note that LTS_v1 simply did the equivalent of: ! seismicity() = seismicity() + PB2002_seismicity() ! which could double-count some plate-boundary-adjacent deformation. IF (format_selector == 2) THEN ! read template file to set template_mask(rows, columns) = .TRUE. OPEN (UNIT = 11, file = TRIM(template_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (8, "(/' ERROR: File not found under this [Drive:][\path\]filename.')") CALL Pause() STOP END IF READ (11, *) ! modelname= READ (11, *) ! version= READ (11, *) ! author= READ (11, *) ! issue_date= READ (11, *) ! forecast_start_date= READ (11, *) ! forecast_duration= READ (11, *) ! begin_forecast indefinite: DO READ (11, "(A)") line IF (line(1:3) == "end") EXIT indefinite READ (line, *) low_template_longitude, high_template_longitude, low_template_latitude, high_template_latitude center_longitude = (low_template_longitude + high_template_longitude) / 2. center_latitude = (low_template_latitude + high_template_latitude ) / 2. row = 1 + NINT((lat_max - center_latitude) / d_lat) column = 1 + NINT((center_longitude - lon_min) / d_lon) IF ((row >= 1).AND.(row <= rows)) THEN IF ((column >= 1).AND.(column <= columns)) THEN template_mask(row, column) = .TRUE. END IF END IF !N.B. This algorithm sets the same mask element to .TRUE. multiple times (41x), but that does no harm. END DO indefinite CLOSE (UNIT = 11, DISP = "KEEP") ! template_pathfile END IF ALLOCATE ( continental(rows, columns) ) ! NOTE: rows are numbered from 1 in the N to "rows" in the S; ! columns are numbered from 1 in the W to "columns" in the E. continental = .FALSE. ! whole array; to be changed to .TRUE. in continental areas (below) !============================================================================================================= ! Decide which grid points are "continental," using criteria of ! Bird [2003, Geochemistry Geophysics Geosystems, 4(3), 1027, doi:10.1029/2001GC000252] !Read 2: ETOPO5.grd (whole planet with 5' grid: 4320 columns, 2161 rows, 17.8 MB as INTEGER(KIND=2)) !INTEGER(KIND = 2), DIMENSION(-1080:1080, 0:4319) :: ETOPO5 ! Note that row/column 0 means lat/lon of 0 OPEN (UNIT = 2, FILE = 'ETOPO5.grd', STATUS = 'OLD', DISP = 'KEEP', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: input file ETOPO5.grd could not be opened.')") CALL Traceback() STOP END IF WRITE(*, "(/' Reading ETOPO5.grd (and packing 63.9 MB into 17.8 MB)...')") READ (2, *) ! 0 8.333333E-02 359.9167 READ (2, *) ! -90 8.333333E-02 89.99999 READ (2, *, IOSTAT = ios) ((ETOPO5(i, j), j = 0, 4319), i = 1080, -1080, -1) IF (ios /= 0) THEN WRITE (*,"(' ERROR: could not read (all?) of ETOPO5.grd')") CALL Traceback() STOP END IF CLOSE (2) !Read 3: age_1p5.grd (from Mueller et al., 1997: 3601 columns, 1621 rows (+90 to -72N), 11.1 MB as INTEGER(KIND=2)) !INTEGER(KIND = 2), DIMENSION( -720:900, 0:3600) :: age1p5 ! Note that row/column 0 means lat/lon of 0 OPEN (UNIT = 3, FILE = 'age_1p5.grd', STATUS = 'OLD', DISP = 'KEEP', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: input file age_1p5.grd could not be opened.')") CALL Traceback() STOP END IF WRITE(*, "(' Reading age_1p5.grd (and packing 28.9 MB into 11.1 MB)...')") READ (3, *) ! 0 0.1 360.0 READ (3, *) ! -72 0.1 90 READ (3, *, IOSTAT = ios) ((age_1p5(i, j), j = 0,3600), i = 900, -720, -1) IF (ios /= 0) THEN WRITE (*,"(' ERROR: could not read (all?) of age_1p5.grd')") CALL Traceback() STOP END IF CLOSE (3) !KLUDGE: Rectangular region S of New Zealand is "unknown ocean floor" ! (220 Ma; = my special code; c.f., Mueller used 327.xx) ! in age_1p5.grd of Mueller et al., 1997 simply ! because of lack of data. Where AU-PA boundary passes over ! Macquarie Ridge and goes higher than -2000 m, this causes it ! to be called "continental" which is clearly wrong, because ! Macquarie Island is an ophiolite, and papers like Massell ! et al., 2000 and Mosher, 2002 clearly show that it is ocean ! floor formed by spreading more recently than anomaly 12 ! = 33 Ma. Therefore, this block of code resets "220" ages ! in that region to "30". DO i = -720, -500 ! -72 to -50 N = 72 to 50 S DO j = 1580, 1660 ! 158 to 166 E IF (age_1p5(i, j) == 220) age_1p5(i, j) = 30 END DO END DO !KLUDGE: Ocean-continent boundary offshore Gorda-California-Nevada orogen ! is seriously misplaced (too far East; along coast) in age_1p5.grd. ! The corrected boundary below follows Muehlberger [1992] and Crouch & Suppe [1993]. GCNocb(1:2, 1) = (/ -1.25968E+02,+4.03814E+01 /) GCNocb(1:2, 2) = (/ -1.25835E+02,+4.02391E+01 /) GCNocb(1:2, 3) = (/ -1.25675E+02,+4.01204E+01 /) GCNocb(1:2, 4) = (/ -1.25511E+02,+3.99907E+01 /) GCNocb(1:2, 5) = (/ -1.25354E+02,+3.98269E+01 /) GCNocb(1:2, 6) = (/ -1.25170E+02,+3.96724E+01 /) GCNocb(1:2, 7) = (/ -1.24941E+02,+3.95002E+01 /) GCNocb(1:2, 8) = (/ -1.24741E+02,+3.93567E+01 /) GCNocb(1:2, 9) = (/ -1.24598E+02,+3.92667E+01 /) GCNocb(1:2, 10) = (/ -1.24508E+02,+3.91928E+01 /) GCNocb(1:2, 11) = (/ -1.24397E+02,+3.90348E+01 /) GCNocb(1:2, 12) = (/ -1.24303E+02,+3.88736E+01 /) GCNocb(1:2, 13) = (/ -1.24197E+02,+3.86841E+01 /) GCNocb(1:2, 14) = (/ -1.24094E+02,+3.84923E+01 /) GCNocb(1:2, 15) = (/ -1.24016E+02,+3.83513E+01 /) GCNocb(1:2, 16) = (/ -1.23923E+02,+3.82036E+01 /) GCNocb(1:2, 17) = (/ -1.23824E+02,+3.80420E+01 /) GCNocb(1:2, 18) = (/ -1.23731E+02,+3.78889E+01 /) GCNocb(1:2, 19) = (/ -1.23623E+02,+3.77350E+01 /) GCNocb(1:2, 20) = (/ -1.23524E+02,+3.75699E+01 /) GCNocb(1:2, 21) = (/ -1.23451E+02,+3.74379E+01 /) GCNocb(1:2, 22) = (/ -1.23373E+02,+3.72974E+01 /) GCNocb(1:2, 23) = (/ -1.23323E+02,+3.71441E+01 /) GCNocb(1:2, 24) = (/ -1.23299E+02,+3.70539E+01 /) GCNocb(1:2, 25) = (/ -1.23282E+02,+3.69816E+01 /) GCNocb(1:2, 26) = (/ -1.23270E+02,+3.69137E+01 /) GCNocb(1:2, 27) = (/ -1.23171E+02,+3.68118E+01 /) GCNocb(1:2, 28) = (/ -1.23016E+02,+3.66792E+01 /) GCNocb(1:2, 29) = (/ -1.22893E+02,+3.65894E+01 /) GCNocb(1:2, 30) = (/ -1.22793E+02,+3.64707E+01 /) GCNocb(1:2, 31) = (/ -1.22690E+02,+3.63396E+01 /) GCNocb(1:2, 32) = (/ -1.22591E+02,+3.62227E+01 /) GCNocb(1:2, 33) = (/ -1.22501E+02,+3.60703E+01 /) GCNocb(1:2, 34) = (/ -1.22434E+02,+3.59626E+01 /) GCNocb(1:2, 35) = (/ -1.22359E+02,+3.58358E+01 /) GCNocb(1:2, 36) = (/ -1.22269E+02,+3.56649E+01 /) GCNocb(1:2, 37) = (/ -1.22221E+02,+3.55482E+01 /) GCNocb(1:2, 38) = (/ -1.22181E+02,+3.54432E+01 /) GCNocb(1:2, 39) = (/ -1.22150E+02,+3.53087E+01 /) GCNocb(1:2, 40) = (/ -1.22107E+02,+3.52175E+01 /) GCNocb(1:2, 41) = (/ -1.22067E+02,+3.51427E+01 /) GCNocb(1:2, 42) = (/ -1.22024E+02,+3.50798E+01 /) GCNocb(1:2, 43) = (/ -1.22034E+02,+3.50349E+01 /) GCNocb(1:2, 44) = (/ -1.22045E+02,+3.49768E+01 /) GCNocb(1:2, 45) = (/ -1.21965E+02,+3.48741E+01 /) GCNocb(1:2, 46) = (/ -1.21868E+02,+3.47430E+01 /) GCNocb(1:2, 47) = (/ -1.21781E+02,+3.46304E+01 /) GCNocb(1:2, 48) = (/ -1.21687E+02,+3.45081E+01 /) GCNocb(1:2, 49) = (/ -1.21584E+02,+3.43949E+01 /) GCNocb(1:2, 50) = (/ -1.21484E+02,+3.42887E+01 /) GCNocb(1:2, 51) = (/ -1.21417E+02,+3.41524E+01 /) GCNocb(1:2, 52) = (/ -1.21325E+02,+3.40060E+01 /) GCNocb(1:2, 53) = (/ -1.21277E+02,+3.39203E+01 /) GCNocb(1:2, 54) = (/ -1.21205E+02,+3.38648E+01 /) GCNocb(1:2, 55) = (/ -1.21118E+02,+3.38285E+01 /) GCNocb(1:2, 56) = (/ -1.21001E+02,+3.37476E+01 /) GCNocb(1:2, 57) = (/ -1.20894E+02,+3.36446E+01 /) GCNocb(1:2, 58) = (/ -1.20826E+02,+3.35294E+01 /) GCNocb(1:2, 59) = (/ -1.20732E+02,+3.33800E+01 /) GCNocb(1:2, 60) = (/ -1.20698E+02,+3.32664E+01 /) GCNocb(1:2, 61) = (/ -1.20672E+02,+3.31718E+01 /) GCNocb(1:2, 62) = (/ -1.20617E+02,+3.30730E+01 /) GCNocb(1:2, 63) = (/ -1.20548E+02,+3.29729E+01 /) GCNocb(1:2, 64) = (/ -1.20499E+02,+3.29209E+01 /) GCNocb(1:2, 65) = (/ -1.20458E+02,+3.27960E+01 /) GCNocb(1:2, 66) = (/ -1.20395E+02,+3.27209E+01 /) GCNocb(1:2, 67) = (/ -1.20311E+02,+3.26439E+01 /) GCNocb(1:2, 68) = (/ -1.20243E+02,+3.25257E+01 /) GCNocb(1:2, 69) = (/ -1.20216E+02,+3.24263E+01 /) GCNocb(1:2, 70) = (/ -1.20160E+02,+3.23445E+01 /) GCNocb(1:2, 71) = (/ -1.20068E+02,+3.22490E+01 /) GCNocb(1:2, 72) = (/ -1.19982E+02,+3.21905E+01 /) GCNocb(1:2, 73) = (/ -1.19898E+02,+3.21246E+01 /) GCNocb(1:2, 74) = (/ -1.19834E+02,+3.20625E+01 /) GCNocb(1:2, 75) = (/ -1.19770E+02,+3.20060E+01 /) GCNocb(1:2, 76) = (/ -1.19656E+02,+3.18905E+01 /) GCNocb(1:2, 77) = (/ -1.19461E+02,+3.17646E+01 /) GCNocb(1:2, 78) = (/ -1.19329E+02,+3.16775E+01 /) GCNocb(1:2, 79) = (/ -1.19170E+02,+3.15764E+01 /) GCNocb(1:2, 80) = (/ -1.19049E+02,+3.14727E+01 /) GCNocb(1:2, 81) = (/ -1.18919E+02,+3.13905E+01 /) GCNocb(1:2, 82) = (/ -1.18787E+02,+3.13157E+01 /) GCNocb(1:2, 83) = (/ -1.18696E+02,+3.12471E+01 /) GCNocb(1:2, 84) = (/ -1.18586E+02,+3.11740E+01 /) GCNocb(1:2, 85) = (/ -1.17693E+02,+3.15128E+01 /) GCNocb(1:2, 86) = (/ -1.17904E+02,+3.17194E+01 /) GCNocb(1:2, 87) = (/ -1.18147E+02,+3.20165E+01 /) GCNocb(1:2, 88) = (/ -1.18366E+02,+3.22861E+01 /) GCNocb(1:2, 89) = (/ -1.18580E+02,+3.26362E+01 /) GCNocb(1:2, 90) = (/ -1.18705E+02,+3.28368E+01 /) GCNocb(1:2, 91) = (/ -1.18780E+02,+3.30791E+01 /) GCNocb(1:2, 92) = (/ -1.18887E+02,+3.32703E+01 /) GCNocb(1:2, 93) = (/ -1.19061E+02,+3.34490E+01 /) GCNocb(1:2, 94) = (/ -1.18990E+02,+3.34972E+01 /) GCNocb(1:2, 95) = (/ -1.18789E+02,+3.32981E+01 /) GCNocb(1:2, 96) = (/ -1.18589E+02,+3.30846E+01 /) GCNocb(1:2, 97) = (/ -1.18373E+02,+3.28510E+01 /) GCNocb(1:2, 98) = (/ -1.18144E+02,+3.26726E+01 /) GCNocb(1:2, 99) = (/ -1.17951E+02,+3.25006E+01 /) GCNocb(1:2, 100) = (/ -1.17884E+02,+3.24360E+01 /) GCNocb(1:2, 101) = (/ -1.17612E+02,+3.22278E+01 /) GCNocb(1:2, 102) = (/ -1.17381E+02,+3.20332E+01 /) GCNocb(1:2, 103) = (/ -1.17206E+02,+3.18384E+01 /) GCNocb(1:2, 104) = (/ -1.17018E+02,+3.16450E+01 /) GCNocb(1:2, 105) = (/ -1.16863E+02,+3.14435E+01 /) GCNocb(1:2, 106) = (/ -1.16758E+02,+3.12846E+01 /) GCNocb(1:2, 107) = (/ -1.16693E+02,+3.10982E+01 /) GCNocb(1:2, 108) = (/ -1.16566E+02,+3.08563E+01 /) GCNocb(1:2, 109) = (/ -1.16465E+02,+3.06240E+01 /) GCNocb(1:2, 110) = (/ -1.16363E+02,+3.03666E+01 /) GCNocb(1:2, 111) = (/ -1.16240E+02,+3.01302E+01 /) GCNocb(1:2, 112) = (/ -1.16111E+02,+2.98577E+01 /) GCNocb(1:2, 113) = (/ -1.16594E+02,+2.96056E+01 /) GCNocb(1:2, 114) = (/ -1.17013E+02,+3.00922E+01 /) GCNocb(1:2, 115) = (/ -1.17089E+02,+3.01516E+01 /) GCNocb(1:2, 116) = (/ -1.17129E+02,+3.00398E+01 /) GCNocb(1:2, 117) = (/ -1.17133E+02,+2.99323E+01 /) GCNocb(1:2, 118) = (/ -1.17129E+02,+2.98031E+01 /) GCNocb(1:2, 119) = (/ -1.17037E+02,+2.96303E+01 /) GCNocb(1:2, 120) = (/ -1.16894E+02,+2.94435E+01 /) GCNocb(1:2, 121) = (/ -1.16761E+02,+2.93018E+01 /) GCNocb(1:2, 122) = (/ -1.16629E+02,+2.91350E+01 /) GCNocb(1:2, 123) = (/ -1.16476E+02,+2.89479E+01 /) GCNocb(1:2, 124) = (/ -1.16269E+02,+2.86913E+01 /) GCNocb(1:2, 125) = (/ -1.16107E+02,+2.84753E+01 /) GCNocb(1:2, 126) = (/ -1.15854E+02,+2.81582E+01 /) GCNocb(1:2, 127) = (/ -1.15616E+02,+2.78720E+01 /) GCNocb(1:2, 128) = (/ -1.15364E+02,+2.76446E+01 /) GCNocb(1:2, 129) = (/ -1.15194E+02,+2.74767E+01 /) GCNocb(1:2, 130) = (/ -1.15080E+02,+2.73858E+01 /) DO i = 274, 403 ! 27.4 to 40.3 N latitude = i * 0.10 DO j = 2340, 2450 ! 234.0 E to 245.0 E = -126.0 to -115.0 E = 126.0 to 115.0 W longitude = j * 0.10 - 360.0 crossings_to_W = 0 ! initialize sum DO k = 2, 130 ! all digitization steps in string above IF (((GCNocb(2, k-1) - latitude) * (GCNocb(2, k) - latitude)) < 0.0) THEN ! digitization step spans this latitude fraction = (latitude - GCNocb(2, k-1)) / (GCNocb(2, k) - GCNocb(2, k-1)) edge_longitude = GCNocb(1, k-1) + fraction * (GCNocb(1, k) - GCNocb(1, k-1)) IF (longitude > edge_longitude) crossings_to_W = crossings_to_W + 1 END IF END DO IF (MOD(crossings_to_W, 2) == 1) THEN ! "inside" (E or NE) of the line ! Set to "continental" code-age of 260 Ma: age_1p5(i, j) = 260 ! (But, note that this point may still be labelled as "oceanic" if elevation is less than -2000 m.) END IF END DO END DO DO i = 1, rows latitude = lat_max - (i - 1) * d_lat DO j = 1, columns longitude = lon_min + (j - 1) * d_lon !elevation at grid point: elevation = NINT(LookUp_ETOPO5(longitude, latitude)) !seafloor age at grid point: age_Ma = NINT(LookUp_age_1p5(longitude, latitude)) !decision: oceanic = (age_Ma < 180).OR.(elevation < -2000) continental (i, j) = .NOT. oceanic END DO END DO !============================================================================================================= ! Pre-compute unit vectors ("uvecs") for all subsampling points at all grid points: MB = 4 * 3 * subsamplings * subsamplings * rows * columns / 1048576 WRITE (*, "(' Memory required for subsample_uvec = ',I6,' MB.')") MB ALLOCATE ( subsample_uvec (3, subsamplings, subsamplings, rows, columns) ) WRITE (*,"(' Computing unit vectors at all subsampling points...')") DO i = 1, rows latitude = lat_max - (i - 1) * d_lat WRITE (*,"('+Computing unit vectors at all subsampling points...',I6,' rows out of ',I6)") i, rows DO j = 1, columns longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings t_lat = latitude + lat_offset + ii * dd_lat t_lat = MIN(90.0, MAX(-90.0, t_lat)) DO jj = 1, subsamplings t_lon = longitude + lon_offset + jj * dd_lon CALL LonLat_2_Uvec(t_lon, t_lat, uvec) subsample_uvec(1:3, ii, jj, i, j) = uvec(1:3) END DO END DO END DO END DO WRITE (*,"('+Computing unit vectors at all subsampling points...DONE ')") !============================================================================================================= ! INTRAPLATE SEISMICITY (based on strain-rates in continuum elements from ! a global dynamic model computed with Shells, e.g., ! Earth3 = Bird & Liu [1999, BSSA, 89(6), 1642-1647], or ! Earth5 = Bird et al.[2008, JGR, 113(B11), B11406, doi:10.1029/2007JB005460] ! or a later (2009) improvement of the Earth5 model in which heat-flow ! has been manually reduced in local regions to suppress landslides.) !Read the .feg file WRITE (*, "(' Reading ',A,'...')") TRIM(global_feg_pathfile) OPEN (UNIT = 4, FILE = global_feg_pathfile, STATUS = 'OLD', DISP = 'KEEP', IOSTAT = ios, PAD = 'YES') problem = (ios /= 0) READ (4, *, IOSTAT = ios) ! title line problem = problem .OR. (ios /= 0) READ (4, *, IOSTAT = ios) numnod problem = problem .OR. (ios /= 0) ALLOCATE ( node_uvec(3,numnod) ) DO i = 1, numnod READ (4, *, IOSTAT = ios) j, longitude, latitude problem = problem .OR. (ios /= 0) .OR. (j /= i) CALL LonLat_2_Uvec (longitude, latitude, uvec) node_uvec(1:3,i) = uvec(1:3) END DO ! i = 1, numnod READ (4, *, IOSTAT = ios) numel problem = problem .OR. (ios /= 0) ALLOCATE ( nodes(3,numel) ) DO i = 1, numel READ (4, *, IOSTAT = ios) j, nodes(1,i), nodes(2,i), nodes(3,i) problem = problem .OR. (ios /= 0) END DO ! i = 1, numel IF (problem) THEN WRITE (*,"(' ERROR: Necessary information absent or defective in file ',A)") TRIM(global_feg_pathfile) CLOSE (4) CALL Pause() STOP END IF !Note: Fault elements in this .feg file will not be used, so they are not read. ! Instead, the plate-boundary slip rates will be taken either from ! the PB2002 rigid-plate model of Bird [2003], or from local ! studies of individual orogens done with NeoKinema (version 2.0+). CLOSE (4) !Read the nodal velocities from the global solution: OPEN(UNIT = 5, FILE = global_v_pathfile, STATUS = 'OLD', DISP = 'KEEP', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File ',A,' not found.')") TRIM(global_v_pathfile) CLOSE (5) CALL Pause() STOP END IF DO i = 1, 3 READ (5,"(A)") line END DO ! first 3 lines of velocity file ALLOCATE ( vw(2*numnod) ) READ (5,*) (vw(i), i = 1, (2*numnod)) CLOSE (5) !Compute strain-rates in the horizontal plane at centers of continuum elements: ALLOCATE ( element_principal_strainrates(3, numel) ) DO l_ = 1, numel ! uvec1(1:3) = node_uvec(1:3, nodes(1, l_)) uvec2(1:3) = node_uvec(1:3, nodes(2, l_)) uvec3(1:3) = node_uvec(1:3, nodes(3, l_)) m = 1 ! (not using loop, since values for m =1,...,7 almost identical ! evaluate nodal function and derivitives uvec4(1:3) = Gauss_point(1,m) * node_uvec(1:3, nodes(1,l_)) + & & Gauss_point(2,m) * node_uvec(1:3, nodes(2,l_)) + & & Gauss_point(3,m) * node_uvec(1:3, nodes(3,l_)) CALL Make_Uvec (uvec4, uvec) ! center of element equat = SQRT(uvec(1)**2 + uvec(2)**2) IF (equat == 0.) THEN PRINT "(' Error: integration point ',I1,' of element ',I5,' is N or S pole.')", m, l_ CALL Pause() STOP END IF theta_ = ATAN2(equat, uvec(3)) CALL Gjxy (l_, uvec1, & & uvec2, & & uvec3, & ! INTENT(IN) & uvec, G) ! INTENT(OUT) CALL Del_Gjxy_del_thetaphi (l_, uvec1, & & uvec2, & & uvec3, & ! INTENT(IN) & uvec, dG) ! INTENT(OUT) CALL E_rate(Earth_radius_m, l_, nodes, G, dG, theta_, vw, & ! INTENT(IN) & eps_dot) ! INTENT(OUT) ! eps_dot(1) = epsilon dot sub theta theta = epsilon dot sub S S ! eps_dot(2) = epsilon dot sub theta phi = epsilon dot sub S E ! eps_dot(3) = epsilon dot sub phi phi = epsilon dot sub E E CALL Principal_Axes_22 (eps_dot(1),eps_dot(2),eps_dot(3), & & e1h, e2h, u1theta,u1phi, u2theta,u2phi) err = -(e1h + e2h) ! from conservation of mass & volume ! (remember, these are anelastic, permanent strain rates) element_principal_strainrates(1, l_) = e1h element_principal_strainrates(2, l_) = e2h element_principal_strainrates(3, l_) = err !(not using loop on m = 1,...,7 since values almost identical) END DO ! l_ = 1, numel, computing strainrates at element centers !Transfer the seismicity values to the .grd points (using subsampling): ALLOCATE ( a_(numel) ) ALLOCATE ( center(3, numel) ) ALLOCATE ( neighbor(3, numel) ) CALL Learn_Spherical_Triangles (numel, nodes, node_uvec, .TRUE., & & a_, center, neighbor) WRITE (*,"(' Computing seismicity rates in continuum...')") DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat WRITE (*,"('+Computing seismicity rates in continuum...',I6,' rows out of ',I6)") i, rows cold_start = .TRUE. DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat !t_lat = MIN(90.0, MAX(-90.0, t_lat)) DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) CALL Which_Spherical_Triangle (uvec, cold_start, & & numel, nodes, node_uvec, center, a_, neighbor, & & success, iele, s1, s2, s3) cold_start = .NOT.success ! for the next time through these loops IF (success) THEN !unpack the principal strain-rates: e1h = element_principal_strainrates(1, iele) e2h = element_principal_strainrates(2, iele) err = element_principal_strainrates(3, iele) !determine appropriate column in Table 5: IF (continental(i, j)) THEN IF (err == 0.0) THEN ! rare: ideal strike-slip T5_column = 2 ! CTF ELSE IF (err > 0.0) THEN ! thrusting (possibly minor?) IF (err > (0.364 * e2h)) THEN ! primarily thrusting T5_column = 3 ! slow CCB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for CTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 2 ! CTF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (0.364 * e1h)) THEN ! primarily normal-faulting T5_column = 1 ! CRB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for CTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 2 ! CTF END IF END IF ! thrusting or normal-faulting? ELSE ! oceanic IF (err == 0.0) THEN ! rare: ideal strike-slip T5_column = 7 ! slow OTF ELSE IF (err > 0.0) THEN ! thrusting (possibly minor?) IF (err > (0.364 * e2h)) THEN ! primarily thrusting T5_column = 10 ! OCB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for OTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 7 ! slow OTF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (0.364 * e1h)) THEN ! primarily normal-faulting T5_column = 10 ! OCB [Sic. Not OSR, because very slow-extending ! continuum regions will not have the thermal ! structure of a spreading ridge, but a cold structure ! more similar to that of an OCB.] ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for OTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 7 ! slow OTF END IF END IF ! thrusting or normal-faulting? END IF ! continental, or oceanic !convert principal strain rates in continuum to seismic moment rate per unit volume 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 postive; negative e1_rate is the unique axis: M_persec_per_m3 = T5_assumed_mu_Pa(T5_column) * 2.0 * (-e1_rate) END IF !N.B. {Note added 2006.10.26}: ! The formula above was changed on this date from that of equation (7a) of Bird & Liu [2007; SRL] ! to that of equation (7b) of the same paper. !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) * T5_supplement_seismicity_fraction_based_on_z(T5_column) !determine subsampling-point value of seismicity (earthquakes/m**2/s) and add to array: 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/magnitude, !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) DO threshold_index = 1, top_threshold_index ! N.B. Looping only occurs for format_selector == 2 IF (format_selector == 2) desired_threshold_moment_Nm = RELM_threshold_moment_Nm(threshold_index) 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 point_value = G_function * seismicity_at_CMT_threshold !GPBkludge: !point_value = 1E-25 seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + weight * point_value END DO END IF END DO END DO END DO END DO WRITE (*,"('+Computing seismicity rates in continuum...DONE ')") DEALLOCATE ( neighbor ) DEALLOCATE ( center ) DEALLOCATE ( a_ ) DEALLOCATE ( element_principal_strainrates ) DEALLOCATE ( vw ) DEALLOCATE ( nodes ) DEALLOCATE ( node_uvec ) !============================================================================================================= !ZERO-OUT GRID AREAS COVERED BY OROGENS that were modeled in detail with NeoKinema (version 2.0+): IF (orogens > 0) THEN ALLOCATE ( n_in_outline(orogens) ) max_in_outline = 0 DO orogen = 1, orogens ! read all outline.dig files once to determine maximum length OPEN (UNIT = 8, FILE = orogen_outline_dig_pathfile(orogen), STATUS = 'OLD', DISP = 'KEEP') READ (8, *) ! title line n_in_outline(orogen) = 0 count_out: DO READ (8, *, IOSTAT = ios) longitude, latitude IF (ios /= 0) EXIT count_out n_in_outline(orogen) = n_in_outline(orogen) + 1 END DO count_out CLOSE (8) max_in_outline = MAX(max_in_outline, n_in_outline(orogen)) END DO ALLOCATE ( outline_uvecs(3, max_in_outline, orogens) ) DO orogen = 1, orogens ! read again, and memorize !Note: Whether this outline.dig file was part of PB2002_orogens.dig, ! or was produced by NeoKineMap (along with overlay type 2), ! it should have one title line (not used), ! it should circle the orogen in the counterclockwise direction ! (as seen from outside the Earth), and ! its last digitized point should be the same as its first digitized point. OPEN (UNIT = 8, FILE = orogen_outline_dig_pathfile(orogen), STATUS = 'OLD', DISP = 'KEEP') READ (8, *) ! title line DO n = 1, n_in_outline(orogen) READ (8, *, IOSTAT = ios) longitude, latitude CALL LonLat_2_Uvec(longitude, latitude, uvec) outline_uvecs(1:3, n, orogen) = uvec(1:3) END DO CLOSE (8) END DO ! orogen = 1, orogens WRITE (*,"(' Zeroing seismicity in orogen(s)...')") DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat WRITE (*,"('+Zeroing seismicity in orogen(s)...',I6,' rows out of ',I6)") i, rows DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon fraction_inside = 0.0 ! but may be increased in inner loops below: DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) inside = .FALSE. ! initialization, before testing this one subsampling point zones: DO orogen = 1, orogens !decide whether subsampling point is "inside" this orogen: decision_in = .FALSE. !avoid abend by first testing for point ON the boundary: edging: DO k = 1, n_in_outline(orogen) uvec1(1:3) = outline_uvecs(1:3, k, orogen) IF (Same_Uvec(uvec, uvec1)) THEN inside = .TRUE. decision_in = .TRUE. EXIT edging END IF END DO edging IF (.NOT.decision_in) THEN angle_sum = 0.0D0 ! initialize sum of angles subtended by orogen steps, as seen from test point: DO k = 2, n_in_outline(orogen) uvec1(1:3) = outline_uvecs(1:3, k-1, orogen) uvec2(1:3) = outline_uvecs(1:3, k, orogen) t1 = Relative_Compass(from_uvec = uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = uvec, to_uvec = uvec2) IF ((t2 - t1) > Pi) t2 = t2 - Two_Pi IF ((t2 - t1) < -Pi) t2 = t2 + Two_Pi angle_sum = angle_sum + t2 - t1 END DO ! k = 2, n_in_outline inside = ((angle_sum < -6.0D0).AND.(angle_sum > -6.6D0)) ! inside a counterclockwise circuit !Note: Generalizing formula to allow for being inside a clockwise circuit will unfortunately ! cast virtual images of each orogen on the far side of the Earth! !Note: Formula gives unpredictable result in the special case where the test point ! is exactly on the boundary line. This is another reason why subsamplings > 1 ! is a desirable PARAMETER setting. END IF IF (inside) THEN fraction_inside = fraction_inside + weight EXIT zones END IF END DO zones END DO ! jj END DO ! ii DO threshold_index = 1, top_threshold_index ! N.B. Looping only occurs for format_selector == 2 seismicity (i, j, threshold_index) = seismicity(i, j, threshold_index) * (1. - fraction_inside) END DO END DO ! j END DO ! i WRITE (*,"('+Zeroing seismicity in orogen(s)...DONE ')") END IF ! orogens > 0 !============================================================================================================= ! ADD PLATE-BOUNDARY SEISMICITY based on gobal plate model PB2002 of Bird [2003, G**3], ! and seismicity calibration by Bird & Kagan [2004, BSSA]. ! (However, plate-boundary segments that lie within one of the orogens processed in ! this run are omitted, because they would be redundant with NeoKinema-based seismicity ! that will be added later. The '*' bytes in PB2002_steps.dat are NOT used to make ! this decision, because some orogen boundaries may have changed since Bird [2003].) OPEN (UNIT = 7, FILE = TRIM(PB2002_steps_pathfile), STATUS = "OLD", DISP = "KEEP", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(/' ERROR: ',A,' not found (in current folder/directory).')") TRIM(PB2002_steps_pathfile) CALL Pause() STOP END IF WRITE (*,"(' Stepping through ',A,'...')") TRIM(PB2002_steps_pathfile) n_step = 0 stepping: DO !read in one plate-boundary step: n_step = n_step + 1 WRITE (*,"('+Stepping through ',A,'...',I6,' steps completed')") TRIM(PB2002_steps_pathfile), n_step !Note: Following READ, variable names, and FORMAT were copied from program ! PB2002_steps.f90, which created the file. READ (7, 1101, IOSTAT = ios) & & sequence_number, step_continuity_c1, c5, & & old_longitude, old_latitude, longitude, latitude, & & length_km, azimuth_integer, & & velocity_mmpa, velocity_azimuth_integer, & & divergence_mmpa, dextral_mmpa, & & elevation, age_Ma, & & class_continuity_c1, class, star_c1 1101 FORMAT (I4, 1X, A1, A5, & & 1X, F8.3, 1X, F7.3, 1X, F8.3, 1X, F7.3, & & 1X, F5.1, 1X, I3, & & 1X, F5.1, 1X, I3, & & 1X, F6.1, 1X, F6.1, & & 1X, I6, 1X, I3, & & 1X, A1, A3, A1) IF (ios == -1) EXIT stepping ! EOF indicator IF (ios /= 0) THEN WRITE (*, "(/' ERROR: Line ',I6,' in ',A,' has incorrect format.')") n_step, TRIM(PB2002_steps_pathfile) CALL Pause() STOP END IF !compute unit vectors for this step: CALL LonLat_2_Uvec(old_longitude, old_latitude, begin_uvec) CALL LonLat_2_Uvec(longitude, latitude, end_uvec) !determine whether this step in inside (or partially inside) any of the orogens being !processed in this run. If so, omit (or shorten) this step. IF (orogens > 0) THEN !test begin_uvec against all orogens (unless result is already known): IF (step_continuity_c1 == ':') THEN which_orogen_at_beginning = which_orogen_at_end ! because begin_uvec of this step = end_uvec of previous step ELSE ! test to determine which_orogen_at_beginning which_orogen_at_beginning = 0 ! initialize before search inside = .FALSE. ! initialization, before testing this one uvec zones2: DO orogen = 1, orogens angle_sum = 0.0D0 ! initialize sum of angles subtended by orogen steps, as seen from test point: DO k = 2, n_in_outline(orogen) uvec1(1:3) = outline_uvecs(1:3, k-1, orogen) uvec2(1:3) = outline_uvecs(1:3, k, orogen) t1 = Relative_Compass(from_uvec = begin_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = begin_uvec, to_uvec = uvec2) IF ((t2 - t1) > Pi) t2 = t2 - Two_Pi IF ((t2 - t1) < -Pi) t2 = t2 + Two_Pi angle_sum = angle_sum + t2 - t1 END DO ! k = 2, n_in_outline inside = ((angle_sum < -6.0D0).AND.(angle_sum > -6.6D0)) ! inside a counterclockwise circuit !Note: Generalizing formula to allow for being inside a clockwise circuit will unfortunately ! cast virtual images of each orogen on the far side of the Earth! IF (inside) THEN which_orogen_at_beginning = orogen EXIT zones2 END IF END DO zones2 END IF ! recall, or compute, which_orogen_at_beginning IF (which_orogen_at_beginning > 0) THEN ! begin_uvec is in that orogen !test end_uvec in(?) same orogen: inside = .FALSE. ! initialization, before testing this one uvec which_orogen_at_end = 0 ! initialize before search angle_sum = 0.0D0 ! initialize sum of angles subtended by orogen steps, as seen from test point: DO k = 2, n_in_outline(which_orogen_at_beginning) uvec1(1:3) = outline_uvecs(1:3, k-1, which_orogen_at_beginning) uvec2(1:3) = outline_uvecs(1:3, k, which_orogen_at_beginning) t1 = Relative_Compass(from_uvec = end_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = end_uvec, to_uvec = uvec2) IF ((t2 - t1) > Pi) t2 = t2 - Two_Pi IF ((t2 - t1) < -Pi) t2 = t2 + Two_Pi angle_sum = angle_sum + t2 - t1 END DO ! k = 2, n_in_outline inside = ((angle_sum < -6.0D0).AND.(angle_sum > -6.6D0)) ! inside a counterclockwise circuit !Note: Generalizing formula to allow for being inside a clockwise circuit will unfortunately ! cast virtual images of each orogen on the far side of the Earth! !**************************************************************** IF (inside) THEN ! both ends are in the same orogen; which_orogen_at_end = which_orogen_at_beginning ! save for next pass through DO stepping CYCLE stepping ! stop processing and go get the next step !**************************************************************** ELSE ! end_uvec is NOT in the orogen, although begin_uvec was; shorten this step by replacing begin_uvec which_orogen_at_end = 0 first_a_uvec = begin_uvec last_a_uvec = end_uvec ! so, arc "a" is the plate boundary step CALL Cross(first_a_uvec, last_a_uvec, tvec) CALL Make_Uvec(tvec, pole_a_uvec) dot_a = Dot(first_a_uvec, pole_a_uvec) find_exit: DO k = 2, n_in_outline(which_orogen_at_beginning) first_b_uvec(1:3) = outline_uvecs(1:3, k-1, which_orogen_at_beginning) last_b_uvec(1:3) = outline_uvecs(1:3, k, which_orogen_at_beginning) !so, arc "b" is one piece of the outline of the orogen CALL Cross(first_b_uvec, last_b_uvec, tvec) CALL Make_Uvec(tvec, pole_b_uvec) dot_b = Dot(first_b_uvec, pole_b_uvec) CALL Circles_Intersect (pole_a_uvec, dot_a, first_a_uvec, last_a_uvec, & & pole_b_uvec, dot_b, first_b_uvec, last_b_uvec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (number > 0) THEN ! found the exit point begin_uvec = point1_uvec EXIT find_exit END IF END DO find_exit END IF ! end_uvec was in same orogen, or not ELSE ! begin_uvec was NOT in any orogen (but end_uvec might be) !test end_uvec against all orogens: inside = .FALSE. ! initialization, before testing this one uvec which_orogen_at_end = 0 ! initialize before search zones3: DO orogen = 1, orogens angle_sum = 0.0D0 ! initialize sum of angles subtended by orogen steps, as seen from test point: DO k = 2, n_in_outline(orogen) uvec1(1:3) = outline_uvecs(1:3, k-1, orogen) uvec2(1:3) = outline_uvecs(1:3, k, orogen) t1 = Relative_Compass(from_uvec = end_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = end_uvec, to_uvec = uvec2) IF ((t2 - t1) > Pi) t2 = t2 - Two_Pi IF ((t2 - t1) < -Pi) t2 = t2 + Two_Pi angle_sum = angle_sum + t2 - t1 END DO ! k = 2, n_in_outline inside = ((angle_sum < -6.0D0).AND.(angle_sum > -6.6D0)) ! inside a counterclockwise circuit !Note: Generalizing formula to allow for being inside a clockwise circuit will unfortunately ! cast virtual images of each orogen on the far side of the Earth! IF (inside) THEN which_orogen_at_end = orogen EXIT zones3 END IF END DO zones3 IF (which_orogen_at_end > 0) THEN ! end_uvec was in this orogen; shorten this step by replacing end_uvec first_a_uvec = begin_uvec last_a_uvec = end_uvec ! so, arc "a" is the plate boundary step CALL Cross(first_a_uvec, last_a_uvec, tvec) CALL Make_Uvec(tvec, pole_a_uvec) dot_a = Dot(first_a_uvec, pole_a_uvec) find_entry: DO k = 2, n_in_outline(which_orogen_at_end) first_b_uvec(1:3) = outline_uvecs(1:3, k-1, which_orogen_at_end) last_b_uvec(1:3) = outline_uvecs(1:3, k, which_orogen_at_end) !so, arc "b" is one piece of the outline of the orogen CALL Cross(first_b_uvec, last_b_uvec, tvec) CALL Make_Uvec(tvec, pole_b_uvec) dot_b = Dot(first_b_uvec, pole_b_uvec) CALL Circles_Intersect (pole_a_uvec, dot_a, first_a_uvec, last_a_uvec, & & pole_b_uvec, dot_b, first_b_uvec, last_b_uvec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (number > 0) THEN ! found the exit point end_uvec = point1_uvec EXIT find_entry END IF END DO find_entry END IF ! end_uvec was in any orogen, or not END IF ! begin_uvec was in any orogen, or not END IF ! orogens > 0 !Note: still inside DO loop on steps; which_orogen_at_end for this step will usually become ! which_orogen_at_beginning for the next step, IF (step_continuity_c1 == ':'). !determine midpoint of step tvec = 0.5 * (begin_uvec + end_uvec) Call Make_Uvec(tvec, center_uvec) !determine k_class_index (= 1..7) for this step (a subscript which relates to arrays of apparent half-width): k_class_index = 0 look_up: DO n = 1, 7 IF (class == class_c3(n)) THEN k_class_index = n EXIT look_up END IF END DO look_up IF (k_class_index == 0) THEN WRITE (*, "(' ERROR: Illegal class ',A,' for step ',I6)") class, n_step CALL Pause() STOP END IF IF (class == "SUB") THEN ! determine corners of this lane !define lane (parallelogram of great circle arcs) in counterclockwise direction: azimuth_radians = azimuth_integer * radians_per_degree IF (c5(3:3) == '/') THEN uvec1 = begin_uvec uvec2 = end_uvec ELSE ! (c5(i)(3:3) == '\') uvec1 = end_uvec uvec2 = begin_uvec END IF IF (c5(3:3) == '\') THEN ! left plate is subducting subduction_azimuth_1 = velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec1 subduction_azimuth_2 = velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec2 ELSE ! right plate is subducting subduction_azimuth_1 = Pi + velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec1 subduction_azimuth_2 = Pi + velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec2 END IF !find corner points of this lane: CALL Turn_To (azimuth_radians = subduction_azimuth_1 + Pi, & & base_uvec = uvec1, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c1_uvec) lane_uvecs(1:3, 1) = c1_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2 + Pi, & & base_uvec = uvec2, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c2_uvec) lane_uvecs(1:3, 2) = c2_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2, & & base_uvec = uvec2, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c3_uvec) lane_uvecs(1:3, 3) = c3_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_1, & & base_uvec = uvec1, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c4_uvec) lane_uvecs(1:3, 4) = c4_uvec(1:3) lane_uvecs(1:3, 5) = lane_uvecs(1:3, 1) lane_width_km = length_km * ABS(SIN(subduction_azimuth_1 - (azimuth_integer * radians_per_degree))) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" basis_length_km = lane_width_km ELSE ! not a SUB step; k_class_index = 1...6, so corners will not be needed basis_length_km = length_km lane_uvecs(1:3, 1:5) = 0.0 ! just to be tidy; should never be used! END IF ! class(i) == "SUB", or not !Precompute Y_factor_of_step, to save time: temp_integral = 0.0 s1 = 0.0 - 3.0 * alongStep_sigma_km s2 = basis_length_km + 3.0 * alongStep_sigma_km ds = (s2 - s1) / 100 DO n = 0, 100 s_km = s1 + (n / 100.0) * (s2 - s1) temp_integral = temp_integral + ds * ANORDF(s_km / alongStep_sigma_km) * ANORDF((basis_length_km - s_km) / alongStep_sigma_km) END DO Y_factor_of_step = 1 / temp_integral !relate "class" to "T5_column": IF (class == "CRB") THEN T5_column = 1 ELSE IF (class == "CTF") THEN T5_column = 2 ELSE IF (class == "CCB") THEN IF (velocity_mmpa <= 24.2) THEN T5_column = 3 ! slow ELSE T5_column = 4 ! fast END IF ELSE IF (class == "OSR") THEN T5_column = 5 ! normal !NB: However, coupled thickness will be based on equation of Bird et al. [2002]. !NB: Should I separate out the parallel velocity component following column 5 of Table 5 of Bird & Kagan [2004]? ELSE IF (class == "OTF") THEN IF (velocity_mmpa < 39.5) THEN T5_column = 7 ! slow ELSE IF (velocity_mmpa < 68.5) THEN T5_column = 8 ! medium ELSE T5_column = 9 ! fast END IF ELSE IF (class == "OCB") THEN T5_column = 10 ELSE IF (class == "SUB") THEN IF (velocity_mmpa <= 66.) THEN T5_column = 11 ! slow ELSE T5_column = 12 ! slow END IF END IF !determine coupled thickness of lithosphere (from Table 5 of Bird & Kagan [2004]; OSR includes later cz(v) model): IF (class == "OSR") THEN d0_m = 1480. ! from CMT_PB2002_OSR_pure_normal_seismicity.xls, based on Bird & Kagan [2004] vs_mmpa = 19. ! [ibid] coupledThickness_m = d0_m * EXP(-divergence_mmpa / vs_mmpa) ELSE ! all other classes coupledThickness_m = T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) END IF !seismicity rate of this step, in EQs per second, above the CMT threshold magnitude: v_oblique_mmpa = SQRT(dextral_mmpa**2 + (divergence_mmpa / COS(T5_assumed_dip_radians(T5_column)))**2) v_oblique_mps = v_oblique_mmpa * (0.001 / s_per_year) moment_rate_of_step_SI = coupledThickness_m * & ! coupled lithosphere thickness, in m & T5_assumed_mu_Pa(T5_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/SIN(T5_assumed_dip_radians(T5_column))) * & ! csc(theta) & (length_km * 1000.) ! along-strike length, in m seismicity_at_CMT_threshold = T5_CMT_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) !correct seismicity to desired threshold moment/magnitude, !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) DO threshold_index = 1, top_threshold_index ! N.B. Looping only occurs if format_selector == 2 IF (format_selector == 2) desired_threshold_moment_Nm = RELM_threshold_moment_Nm(threshold_index) 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 EQs_per_s(threshold_index) = G_function * seismicity_at_CMT_threshold END DO !compare B function of this step to all subsample_uvec's for possible overlap: DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat !t_lat = MIN(90.0, MAX(-90.0, t_lat)) DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) !reject points that are very far away: dot_test = uvec(1) * center_uvec(1) + uvec(2) * center_uvec(2) + uvec(3) * center_uvec(3) IF (dot_test > 0.987) THEN ! within ~1000 km; more precise distance will be computed !distances to end points of this step: uvec1 = begin_uvec to_uvec1_radians = Arc(uvec, uvec1) uvec2 = end_uvec to_uvec2_radians = Arc(uvec, uvec2) IF (k_class_index < 7) THEN ! step is not a SUB; use normal distance/offset logic: !distance and offset measured via perpendicular great circle: CALL Cross(uvec1, uvec2, tvec) CALL Make_Uvec(tvec, pole_uvec) ! pole of the great-circle for this step offside_radians = ABS(Arc(pole_uvec, uvec) - Pi_over_2) distance_km = radius_km * offside_radians ! always positive offset_km = distance_km ! no distinction for classes other than SUB !s_km measured along length of step (using longitudes relative to pole of step's great circle): t1 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec2) !note: normally, t2 is slightly less than t1 IF (t2 > t1) t2 = t2 - Two_Pi angle = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec) ! to the earthquake !check that angle is on the same cycle as t1, t2: midpoint_compass = (t1 + t2) / 2 IF ((angle - midpoint_compass) > Pi) angle = angle - Two_Pi IF ((midpoint_compass - angle) > Pi) angle = angle + Two_Pi s_radians = t1 - angle s_km = s_radians * radius_km !============================================================================== !|| NOTE important difference from EQ_classification_II.f90: || !|| When choosing earthquakes to associate with plate boundary steps, || !|| half-widths of Table 1 [Bird & Kagan, 2004, BSSA] were strict limits. || !|| However, when projecting long term seismicity (here), the tails of the || !|| distribution are restored. || !============================================================================== !-------------------------------------------------------------------------------------------------------- ! X factor = cross-strike component of spatial PDF B: argument = -2 * (offset_km / h_km(k_class_index))**2 IF (argument > -74.0) THEN X = (1.0 / 0.95) *(0.797885 / h_km(k_class_index)) * EXP(argument) ELSE ! prevent underflow in EXP or elsewhere X = 0.0 END IF !-------------------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------------------------- ! Y factor causes PDF to fall off to ~0.5 as ends of step are approached (and even more beyond the ends) IF ((s_km > (-3.0 * alongStep_sigma_km)).AND.(s_km < (length_km + 3.0 * alongStep_sigma_km))) THEN ! use complex formula: Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((length_km - s_km) / alongStep_sigma_km) ELSE ! too far away; don't risk numerical problems within ANORDF: Y = 0.0 END IF !------------------------------------------------------------------------------------------------------- ELSE ! this is a SUB (k_class_index = 7) step; special offset logic is used, and - distances/offsets mean seaward of the trench! !============================================================================== !|| NOTE important difference from EQ_classification_II.f90: || !|| When choosing earthquakes to associate with SUB zones, landward and || !|| seaward half-widths (of 220 and 135 km, respectively) were strict limits.|| !|| However, when projecting long term seismicity (here), the tails of the || !|| distribution are restored. This is done by omitting the part of the test|| !|| for the LOGICAL value "proximal" that checks whether the test point is || !|| within the ends of the lanes. || !============================================================================== proximal = .TRUE. ! unless (usually) negated below: DO n = 2, 4, 2 ! sides ! Note: In EQ_classification_II.f90, this was "DO n = 1, 4" uvec1(1:3) = lane_uvecs(1:3, n) uvec2(1:3) = lane_uvecs(1:3, n+1) CALL Cross (uvec1, uvec2, tvec) ! direction toward pole is direction into the lane CALL Make_Uvec(tvec, inward_uvec) scalar_product = Dot(uvec, inward_uvec) dot_side(n) = scalar_product proximal = proximal .AND. (scalar_product > -0.01) ! 63.7 km tolerance END DO ! n = 1, 4 sides ! Note: EQ is inside the lane if all 4 values in dot_side(1:4) are postive; ! these values are sines of angles from the sides, measured inwards. IF (proximal) THEN ! compute + or - distance/offset from trench, in km (even if slightly to one side of the lane): IF (c5(3:3) == '/') THEN uvec1 = begin_uvec uvec2 = end_uvec ELSE ! (c5(3:3) == '\') uvec1 = end_uvec uvec2 = begin_uvec END IF fraction_toward_uvec1 = dot_side(2) / (dot_side(2) + dot_side(4)) ! where all dot's are positive if inside; ! note that denominator is positive even if EQ (slightly) outside lane. fraction_toward_uvec1 = MAX(0.0, MIN(1.0, fraction_toward_uvec1)) ! This kludge is necessary to prevent pathalogical ! behavior when SUB steps are very short and oblique. tvec(1:3) = fraction_toward_uvec1 * uvec1(1:3) + & & (1.0 - fraction_toward_uvec1) * uvec2(1:3) CALL Make_Uvec(tvec, trench_uvec) ! at same relative left-right position in this lane distance_km = radius_km * Arc(uvec, trench_uvec) ! always zero or positive; if positive, must apply sign now: IF (distance_km > 0.0) THEN from_trench_azimuth_radians = Relative_Compass (from_uvec = trench_uvec, to_uvec = uvec) IF (COS(from_trench_azimuth_radians - subduction_azimuth_1) < 0.0) THEN ! test point is on seaward side of trench; convert distance to negative: distance_km = -distance_km END IF END IF !Note that "distance" is positive on landward side of trench, but negative on seaward side. offset_km = distance_km - SUB_peakSeismicity_km ! preserving possibility of negative values; just shifting origin !-------------------------------------------------------------------------------------------------------- ! X factor = along-lane component of spatial PDF function B: IF (offset_km >= 0.0) THEN ! landward side of seismicity peak: argument = -0.5 * (offset_km / SUB_landward_sigma_km)**2 IF (argument > -74.0) THEN X = (1.0 / 0.95) * & & (2 * SUB_landward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_landward_sigma_km)) * & & EXP(argument) ELSE X = 0.0 END IF ELSE ! seaward side of seismicity peak (but may not be seaward of the trench!) argument = -0.5 * (offset_km / SUB_seaward_sigma_km)**2 IF (argument > -74.0) THEN X = (1.0 / 0.95) * & & (2 * SUB_seaward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_seaward_sigma_km )) * & & EXP(argument) ELSE X = 0.0 END IF END IF ! which side of seismicity peak? !-------------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------------- ! Y factor = cross-lane (orthogonal to X factor) component of spatial PDF function B: ! To evaluate Y(s_km), note that length_km(i) is not the basis length that we want to use, ! but rather the lane width: lane_width_km = radius_km * (dot_side(2) + dot_side(4)) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" s_km = radius_km * dot_side(4) ! distance (if +) into lane, starting on the uvec1 end IF ((s_km > (-3.0 * alongStep_sigma_km)).AND.(s_km < (lane_width_km + 3.0 * alongStep_sigma_km))) THEN ! use complex formula: Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((lane_width_km - s_km) / alongStep_sigma_km) ELSE ! too far away; don't risk numerical problems within ANORDF: Y = 0.0 END IF !-------------------------------------------------------------------------------------------------------- ELSE ! not proximal; punt X = 0.0 Y = 0.0 END IF END IF ! "normal" offset/distance logic (k_class_index = 1-6), or special logic for SUB steps (k_class_index = 7) B = X * Y ! spatial PDF is a product of cross-trace function and along-trace function !At this point, B is in units of (km)**(-2). B = B * 1.0D-6 !Now, B is in units of m**(-2), which is SI. IF (B > 0.0D0) THEN DO threshold_index = 1, top_threshold_index ! N.B. Looping only occurs if format_selector == 2 PB2002_seismicity(i, j, threshold_index) = PB2002_seismicity(i, j, threshold_index) + weight * B * EQs_per_s(threshold_index) END DO END IF ! B > 0.0D0 END IF ! uvec is close enough to center_uvec for further testing END DO ! jj = 1, subsamplings END DO ! ii = 1, subsamplings END DO ! j = 1, columns END DO ! i = 1, rows END DO stepping ! through PB2002_steps_pathfile DO i = 1, rows DO j = 1, columns DO k = 1, top_threshold_index seismicity(i, j, k) = MAX(seismicity(i, j, k), PB2002_seismicity(i, j, k)) !Note innovation in v2 of LTS: !Gaussian-smoothed seismicity of simple PB2002 plate-boundary steps (outside orogens) !is now compared to continuum seismicity from triangular finite elements of an EarthN .feg, !and the greater is kept. This is better than the v1 method of simply adding them, !as his-res .feg models like Earth5 may display their own plate-edge deformation zones, !and we don't want to double-count this seismicity. END DO END DO END DO WRITE (*,"('+Stepping through ',A,'...DONE ')") TRIM(PB2002_steps_pathfile) CLOSE (7) ! PB2002_steps_pathfile !============================================================================================================= IF (orogens > 0) THEN ALLOCATE ( diffusivity(rows, columns) ) ! NOTE: rows are numbered from 1 in the N to "rows" in the S; ! columns are numbered from 1 in the W to "columns" in the E. ALLOCATE ( single_column(rows) ) ! temporary storage needed for ADE diffusion algorithm ALLOCATE ( single_row(columns) ) DO orogen = 1, orogens ! ADD SEISMICITY DUE TO CONTINUUM STRAIN in a NeoKinema model of an orogen !Read the .feg file WRITE (*, "(' Reading ',A,'...')") TRIM(orogen_feg_pathfile(orogen)) OPEN (UNIT = 4, FILE = orogen_feg_pathfile(orogen), STATUS = 'OLD', DISP = 'KEEP', IOSTAT = ios, PAD = 'YES') problem = (ios /= 0) READ (4, *, IOSTAT = ios) ! title line problem = problem .OR. (ios /= 0) READ (4, *, IOSTAT = ios) numnod problem = problem .OR. (ios /= 0) ALLOCATE ( node_uvec(3,numnod) ) DO i = 1, numnod READ (4, *, IOSTAT = ios) j, longitude, latitude problem = problem .OR. (ios /= 0) .OR. (j /= i) CALL LonLat_2_Uvec (longitude, latitude, uvec) node_uvec(1:3,i) = uvec(1:3) END DO ! i = 1, numnod READ (4, *, IOSTAT = ios) numel problem = problem .OR. (ios /= 0) ALLOCATE ( nodes(3,numel) ) DO i = 1, numel READ (4, *, IOSTAT = ios) j, nodes(1,i), nodes(2,i), nodes(3,i) problem = problem .OR. (ios /= 0) END DO ! i = 1, numel IF (problem) THEN WRITE (*,"(' ERROR: Necessary information absent or defective in this file.')") CLOSE (4) CALL Pause() STOP END IF READ (4, *) nfl IF (nfl > 0) CALL No_Fault_Elements_Allowed() CLOSE (4) !Read the e_.nko file (with continuum strain-rates) WRITE (*, "(' Reading ',A,'...')") TRIM(e_nko_pathfile(orogen)) OPEN(UNIT = 5, FILE = e_nko_pathfile(orogen), STATUS = 'OLD', DISP = 'KEEP', PAD = 'YES', IOSTAT = ios) ALLOCATE ( element_principal_strainrates(3, numel) ) DO l_ = 1, numel ! read continuum strainrates at element centers READ (5, *) i, maybe, eps_dot(1), eps_dot(2), eps_dot(3) !convert to scalar measure, for histogram CALL Principal_Axes_22 (eps_dot(1),eps_dot(2),eps_dot(3), & & e1h, e2h, u1theta,u1phi, u2theta,u2phi) err = -(e1h + e2h) ! from conservation of mass & volume ! (remember, these are anelastic, permanent strain rates) element_principal_strainrates(1, l_) = e1h element_principal_strainrates(2, l_) = e2h element_principal_strainrates(3, l_) = err END DO ! l_ = 1, numel, reading continuum strainrates and finding principal values CLOSE(5) !Transfer the seismicity values to the .grd points (using subsampling): ALLOCATE ( a_(numel) ) ALLOCATE ( center(3, numel) ) ALLOCATE ( neighbor(3, numel) ) CALL Learn_Spherical_Triangles (numel, nodes, node_uvec, .TRUE., & & a_, center, neighbor) WRITE (*,"(' Computing seismicity rates in continuum...')") DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat WRITE (*,"('+Computing seismicity rates in continuum...',I6,' rows out of ',I6)") i, rows cold_start = .TRUE. DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) CALL Which_Spherical_Triangle (uvec, cold_start, & & numel, nodes, node_uvec, center, a_, neighbor, & & success, iele, s1, s2, s3) cold_start = .NOT.success ! for the next time through these loops IF (success) THEN !unpack the principal strain-rates: e1h = element_principal_strainrates(1, iele) e2h = element_principal_strainrates(2, iele) err = element_principal_strainrates(3, iele) !determine appropriate column in Table 5: IF (continental(i, j)) THEN IF (err == 0.0) THEN ! rare: ideal strike-slip T5_column = 2 ! CTF ELSE IF (err > 0.0) THEN ! thrusting (possibly minor?) IF (err > (0.364 * e2h)) THEN ! primarily thrusting T5_column = 3 ! slow CCB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for CTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 2 ! CTF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (0.364 * e1h)) THEN ! primarily normal-faulting T5_column = 1 ! CRB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for CTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 2 ! CTF END IF END IF ! thrusting or normal-faulting? ELSE ! oceanic IF (err == 0.0) THEN ! rare: ideal strike-slip T5_column = 7 ! slow OTF ELSE IF (err > 0.0) THEN ! thrusting (possibly minor?) IF (err > (0.364 * e2h)) THEN ! primarily thrusting T5_column = 10 ! OCB ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for OTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 7 ! slow OTF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (0.364 * e1h)) THEN ! primarily normal-faulting T5_column = 10 ! OCB [Sic. Not OSR, because very slow-extending ! continuum regions will not have the thermal ! structure of a spreading ridge, but a cold structure ! more similar to that of an OCB.] ELSE ! vertical strain is minor; primarily strike-slip; ! falls within velocity-azimuth tolerance for OTF ! plate boundaries as defined in ! Bird [2003, Geochemistry Geophysics Geosystems] T5_column = 7 ! slow OTF END IF END IF ! thrusting or normal-faulting? END IF ! continental, or oceanic !convert principal strain rates in continuum to seismic moment rate per unit volume M_persec_per_m3 = 0.0 IF ((e1h * e2h) < 0.0) M_persec_per_m3 = M_persec_per_m3 + T5_assumed_mu_Pa(T5_column) * ABS(e1h - e2h) IF ((e1h * err) < 0.0) M_persec_per_m3 = M_persec_per_m3 + T5_assumed_mu_Pa(T5_column) * ABS(e1h - err) IF ((e2h * err) < 0.0) M_persec_per_m3 = M_persec_per_m3 + T5_assumed_mu_Pa(T5_column) * ABS(e2h - err) !N.B. Faulting only occurs between pairs of principal axes that have ! principal strain rates of opposite sign! (Kostrov got this wrong.) !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) * T5_supplement_seismicity_fraction_based_on_z(T5_column) !determine subsampling-point value of seismicity (earthquakes/m**2/s) and add to array: 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/magnitude, !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) DO threshold_index = 1, top_threshold_index ! N.B. Looping only occurs if format_selector == 2 IF (format_selector == 2) desired_threshold_moment_Nm = RELM_threshold_moment_Nm(threshold_index) 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 point_value = G_function * seismicity_at_CMT_threshold seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + weight * point_value END DO END IF END DO END DO END DO END DO WRITE (*,"('+Computing seismicity rates in continuum...DONE ')") DEALLOCATE ( neighbor ) DEALLOCATE ( center ) DEALLOCATE ( a_ ) DEALLOCATE ( element_principal_strainrates ) DEALLOCATE ( nodes ) DEALLOCATE ( node_uvec ) !============================================================================================================= !GPBhere ! SMOOTH SEISMICITY ACROSS OUTER BOUNDARY OF THE OROGEN: ! Note that outlines of all orogens are already stored in: outline_uvecs(1:3, k, orogen) ! and the upper limit on k is already stored in: n_in_outline(orogen). ! Storage for the "diffusion coefficient" was already allocated in diffusivity(rows, columns), ! but the values need to be replaced each time through this loop, for the new outline: diffusivity = 0.0 ! whole array (rows, columns). Some values will be increased below... WRITE (*,"(' Computing gridded diffusivities of seismicity around boundary of this orogen...')") diffusion_scale = 100./6371. ! (in km/km = radians). !This distance is used to set the half-width of the belt of positive diffusivity, !and then is used again to set the amount of virtual time that elapses during diffusion. !Note: Using only the actual .GRD points, not all the subsampling points (which are for numerical integration). DO k = 1, n_in_outline(orogen) uvec1(1:3) = outline_uvecs(1:3, k, orogen) DO i = 1, rows DO j = 1, columns uvec2(1:3) = subsample_uvec(1:3, middle, middle, i, j) distance = Arc(uvec1, uvec2) scales = distance / diffusion_scale ! both in radians IF (scales < 3.) THEN ! discard values less than 0.0001 diffusivity(i, j) = MAX(diffusivity(i, j), exp(-scales**2)) END IF END DO END DO END DO ! k = 1, n_in_outline(orogen) !Decide row and column limits of the region in which diffusion is permitted: IF (ABS((lon_max - lon_min) - 360.0) < (0.5 * d_lon)) THEN wrap_type = 2 ! first and last column are 360 degrees longitude apart, and overlie each other. ELSE IF (ABS((lon_max - lon_min) - (360.0 - d_lon)) < (0.5 * d_lon)) THEN wrap_type = 1 ! first and last column are (360 - d_lon) degrees longitude apart, and adjacent. ELSE wrap_type = 0 ! first and last columns are far apart and unrelated. END IF IF (lat_max <= 75.0) THEN nr1 = 1 ! diffuse all the way up to Northernmost row ELSE nr1 = 1 + NINT( (lat_max - 75.0) / d_lat ) nr1 = MIN(nr1, rows) END IF IF (lat_min >= -75.0) THEN nr2 = rows ELSE nr2 = 1 + NINT( (lat_max - (-75.0)) / d_lat ) nr2 = MIN(nr2, rows) END IF IF (nr2 > nr1) THEN ! sufficient grid between 75S and 75N limits for diffusion to occur !Use Alternating Direction Explicit method (simple, fast): minimum_arc = MIN( (d_lat * radians_per_degree), & & (d_lon * radians_per_degree * COS(MIN(lat_max, +75.0) * radians_per_degree)), & & (d_lon * radians_per_degree * COS(MAX(lat_min, -75.0) * radians_per_degree)) ) d_S_radians = d_lat * radians_per_degree safe_timestep = 0.1 * minimum_arc**2 ! because diffusivity is always 1 or less total_time = diffusion_scale**2 ! again, based on diffusivity of 1 timesteps_needed = NINT(total_time / safe_timestep) IF (timesteps_needed == 0) THEN safe_timestep = total_time timesteps_needed = 1 END IF WRITE (*, "(' Diffusing seismicity at boundary with ',I8,' virtual timesteps...')") timesteps_needed DO m_time = 1, timesteps_needed !North-South diffusion in one column at a time: factor = safe_timestep / (d_S_radians**2) DO j = 1, columns ! select column DO k = 1, top_threshold_index single_column(nr1:nr2) = seismicity(nr1:nr2, j, k) ! create snapshot of initial condition DO i = nr1, (nr2-1) middle_D = (diffusivity(i, j) + diffusivity(i+1, j)) * 0.5 IF (middle_D > 0.0) THEN difference = single_column(i+1) - single_column(i) seismicity(i, j, k) = seismicity(i, j, k) + factor * middle_D * difference seismicity(i+1, j, k) = seismicity(i+1, j, k) - factor * middle_D * difference END IF END DO ! i = nr1, (nr2-1) END DO ! k = 1, top_threshold_index END DO ! j = 1, columns !East-West diffusion one row at a time: DO i = nr1, nr2 ! select row latitude = lat_max - (i-1) * d_lat d_E_radians = d_lon * radians_per_degree * COS(latitude * radians_per_degree) factor = safe_timestep / (d_E_radians**2) DO k = 1, top_threshold_index single_row(1:columns) = seismicity(i, 1:columns, k) ! create snapshot of initial conditions DO j = 1, (columns-1) middle_D = (diffusivity(i, j) + diffusivity(i, j+1)) * 0.5 IF (middle_D > 0.0) THEN difference = single_row(j+1) - single_row(j) seismicity(i, j, k) = seismicity(i, j, k) + factor * middle_D * difference seismicity(i, j+1, k) = seismicity(i, j+1, k) - factor * middle_D * difference END IF END DO ! j = 1, (columns-1) IF (wrap_type == 1) THEN middle_D = (diffusivity(i, 1) + diffusivity(i, columns)) * 0.5 IF (middle_D > 0.0) THEN difference = single_row(1) - single_row(columns) seismicity(i, columns, k) = seismicity(i, columns, k) + factor * middle_D * difference seismicity(i, 1, k) = seismicity(i, 1, k) - factor * middle_D * difference END IF ELSE IF (wrap_type == 2) THEN average = (seismicity(i, 1, k) + seismicity(i, columns, k)) * 0.5D0 seismicity(i, 1, k) = average seismicity(i, columns, k) = average END IF ! wrapping END DO ! k = 1, top_threshold_index END DO ! i = nr1, nr2 END DO ! m_time = 1, timesteps_needed END IF ! diffusing !============================================================================================================= !ADD SEISMICITY DUE TO FAULT SEGMENTS in a NeoKinema (version 2.0+) models of orogen !GPBdebug: Use line below if outputting outlines of fault patches: !OPEN(UNIT = 72, FILE = 'patches.dig') WRITE (*,"(' Adding seismicity of NeoKinema fault segments...')") !skim h_nko_pathfile once just to count its lines: OPEN (UNIT = 9, FILE = h_nko_pathfile(orogen), STATUS = "OLD", DISP = "KEEP", PAD = "NO") line_number = 0 segmenting: DO ! a loop on the lines in the h_token_nko file: READ (9,"(A6, 13X, F12.3, 10X, F8.3, 1X, F7.3, 3X, F8.3, 1X, F7.3, 11X, I7, 12X, L1, 13X, F12.3)", IOSTAT = ios) & & c6, model_heave_rate_mmpa, lon1, lat1, lon2, lat2, element, creeping, slip_rate_mmpa IF (ios /= 0) THEN IF (line_number < 1) THEN ! bad format or file version WRITE (*, "(' ERROR while trying to read input file ',A,':')") TRIM(h_nko_pathfile(orogen)) WRITE (*, "(' The first line could not be interpreted and/or lacked essential data')") WRITE (*, "(' (LOGICAL creeping? flag, and REAL slip_rate_mmpa).')") WRITE (*, "(' Use only h_[token].nko files produced by NeoKinema versions 2.0+')") WRITE (*, "(' which include this necessary supplementary information.')") WRITE (*, "(' Note that slip_rate_mmpa should be 0.0 for all shadow strike-slip')") WRITE (*, "(' data that were added to increase flexibility on dip-slip faults;')") WRITE (*, "(' this component of the total slip rate is merged into the')") WRITE (*, "(' slip_rate_mmpa value of the previous line (F0000N/D/T/P/S).')") CALL Pause() STOP ELSE ! probably just reached end of file EXIT segmenting END IF END IF line_number = line_number + 1 END DO segmenting CLOSE (9) !allocate storage to memorize the segments in h_nko_pathfile: segments = line_number ALLOCATE ( segment_c6(segments) ) ALLOCATE ( segment_uvecs(3, 2, segments) ) ALLOCATE ( segment_creeping(segments) ) ALLOCATE ( segment_slip_rate_mmpa(segments) ) ALLOCATE ( segment_heave_rate_mmpa(segments) ) !read h_nko_pathfile the second time to memorize the data: OPEN (UNIT = 9, FILE = h_nko_pathfile(orogen), STATUS = "OLD", DISP = "KEEP", PAD = "NO") segment = 0 resegmenting: DO ! a loop on the lines in the h_token_nko file: READ (9,"(A6, 13X, F12.3, 10X, F8.3, 1X, F7.3, 3X, F8.3, 1X, F7.3, 11X, I7, 12X, L1, 13X, F12.3)", IOSTAT = ios) & & c6, model_heave_rate_mmpa, lon1, lat1, lon2, lat2, element, creeping, slip_rate_mmpa IF (ios /= 0) EXIT resegmenting segment = segment + 1 segment_c6(segment) = c6 CALL LonLat_2_Uvec(lon1, lat1, uvec1) segment_uvecs(1:3, 1, segment) = uvec1(1:3) CALL LonLat_2_Uvec(lon2, lat2, uvec2) segment_uvecs(1:3, 2, segment) = uvec2(1:3) segment_creeping(segment) = creeping segment_slip_rate_mmpa(segment) = slip_rate_mmpa segment_heave_rate_mmpa(segment) = model_heave_rate_mmpa ! needed for velocity-dependent cz of OSR's END DO resegmenting CLOSE (9) !find segment_outer_uvecs (to permit computing mean azimuths): ALLOCATE ( segment_used(segments) ) ALLOCATE ( segment_outer_uvecs(3, 2, segments) ) ALLOCATE ( segment_azimuth_radians(segments) ) DO segment = 1, segments !initialize; typically these uvecs will get replaced below, but sometimes not segment_outer_uvecs(1:3, 1:2, segment) = segment_uvecs(1:3, 1:2, segment) segment_used = .FALSE. ! whole array segment_used(segment) = .TRUE. !grow chain toward head end until no more links can be found: heading: DO repeat_sweep = .FALSE. ! unless reset to TRUE below check_others: DO other_segment = 1, segments IF (.NOT.segment_used(other_segment)) THEN IF (segment_c6(other_segment) == segment_c6(segment)) THEN ! other segment belongs to same fault; check for match at head end: uvec1(1:3) = segment_uvecs(1:3, 2, other_segment) ! tail of other segment uvec2(1:3) = segment_outer_uvecs(1:3, 1, segment) ! head of current chain IF (Same_Uvec(uvec1, uvec2)) THEN ! tail of other_segment matches head of current chain segment_outer_uvecs(1:3, 1, segment) = segment_uvecs(1:3, 1, other_segment) segment_used(other_segment) = .TRUE. repeat_sweep = .TRUE. EXIT check_others END IF ! tail of other segment matches head of current chain END IF ! other_segment belongs to same fault END IF END DO check_others ! other_segment = 1, segments IF (.NOT.repeat_sweep) EXIT heading END DO heading !grow chain toward tail end until no more links can be found: tailing: DO repeat_sweep = .FALSE. ! unless reset to TRUE below check_another: DO other_segment = 1, segments IF (.NOT.segment_used(other_segment)) THEN IF (segment_c6(other_segment) == segment_c6(segment)) THEN ! other segment belongs to same fault; check for match at tail end: uvec1(1:3) = segment_uvecs(1:3, 1, other_segment) ! head of other segment uvec2(1:3) = segment_outer_uvecs(1:3, 2, segment) ! tail of current chain IF (Same_Uvec(uvec1, uvec2)) THEN ! head of other_segment matches tail of current chain segment_outer_uvecs(1:3, 2, segment) = segment_uvecs(1:3, 2, other_segment) segment_used(other_segment) = .TRUE. repeat_sweep = .TRUE. EXIT check_another END IF ! head of other segment matches tail of current chain END IF ! other_segment belongs to same fault END IF END DO check_another ! other_segment = 1, segments IF (.NOT.repeat_sweep) EXIT tailing END DO tailing uvec1(1:3) = segment_outer_uvecs(1:3, 1, segment) uvec2(1:3) = segment_outer_uvecs(1:3, 2, segment) IF (Same_Uvec(uvec1, uvec2)) THEN !null, zero-length segment chain segment_azimuth_radians(segment) = 0.0 ELSE segment_azimuth_radians(segment) = Relative_Compass(from_uvec = uvec1, to_uvec = uvec2) END IF END DO ! segment = 1, segments; seeking segment_outer_uvecs !Add seismicity of (positive-rate, non-creeping OR subducting) segments to array, using different ! algorithms for "S" versus non-"S" segments: add_segment_seismicity: DO segment = 1, segments WRITE (*,"('+Adding seismicity of NeoKinema fault segment ',I6,' out of ',I6)") segment, segments IF ((segment_slip_rate_mmpa(segment) > 0.0).AND. & & ((.NOT.segment_creeping(segment)) .OR. (segment_c6(segment)(6:6) == 'S')) & & ) THEN !compute unit vectors for this step: begin_uvec(1:3) = segment_uvecs(1:3, 1, segment) end_uvec(1:3) = segment_uvecs(1:3, 2, segment) tvec = 0.5 * (begin_uvec + end_uvec) Call Make_Uvec(tvec, center_uvec) !compute geometric properties of NeoKinema segment (using same variable names as in SUB code above): length_km = radius_km * Arc(begin_uvec, end_uvec) IF (length_km <= 0.0) CYCLE add_segment_seismicity ! to prevent crash in Compass azimuth_integer = NINT(degrees_per_radian * Compass(from_uvec = begin_uvec, to_uvec = end_uvec)) IF (azimuth_integer < 0) azimuth_integer = azimuth_integer + 360 azimuth_radians = azimuth_integer * radians_per_degree segment_vec = end_uvec - begin_uvec ! NOT a unit vector! ! be CAREFUL not to overwrite these values when scanning PB2002_steps.dat!!! IF (segment_c6(segment)(6:6) == 'S') THEN ! subduction zone !subduction zones in NeoKinema have the same wide footprint !as subduction zones in PB2002 (both per Bird & Kagan [2004, BSSA]). k_class_index = 7 ! SUB IF (segment_slip_rate_mmpa(segment) <= 66.) THEN T5_column = 11 ! slow SUB ELSE T5_column = 12 ! fast SUB END IF !look up velocity_azimuth_integer in PB2002_steps_pathfile (using separation as criterion to choose): OPEN (UNIT = 7, FILE = TRIM(PB2002_steps_pathfile), STATUS = "OLD", DISP = "KEEP", PAD = "YES", IOSTAT = ios) best_separation = 4.0 ! in units of delta_uvec**2; initializing at antipodal value n_step = 0 stepping2: DO !read in one plate-boundary step: READ (7, 1101, IOSTAT = ios) & ! NOTE the name changes: & sequence_number, step_continuity_c1, c5a, & ! c5 --> c5a (so as not to overwrite) & old_longitude, old_latitude, longitude, latitude, & & length_kilometers, azimuth_int, & ! length_km --> length_kilometers; azimuth_integer --> azimuth_int (so as not to overwrite) & velocity_mmpa, velocity_azimuth_int, & ! velocity_azimuth_integer --> velocity_azimuth_int (ditto) & divergence_mmpa, dextral_mmpa, & & elevation, age_Ma, & & class_continuity_c1, class, star_c1 IF (ios == -1) EXIT stepping2 ! EOF indicator n_step = n_step + 1 IF ((c5a(3:3) == '/').OR.(c5a(3:3) == '\')) THEN ! a SUB step was read CALL LonLat_2_Uvec(old_longitude, old_latitude, uvec1) CALL LonLat_2_Uvec(longitude, latitude, uvec2) tvec = 0.5 * (uvec1 + uvec2) CALL Make_Uvec(tvec, uvec3) ! midpoint of step separation = (center_uvec(1) - uvec3(1))**2 + & & (center_uvec(2) - uvec3(2))**2 + & & (center_uvec(3) - uvec3(3))**2 IF (separation < best_separation) THEN best_separation = separation !check whether the segment and the step are parallel or antiparallel: tvec = uvec2 - uvec1 ! progress vector for SUB step just read parallel = (DOT(tvec, segment_vec) > 0.0) IF (parallel) THEN c5 = c5a ! segment value = step value velocity_azimuth_integer = velocity_azimuth_int ! segment value = step value ELSE ! antiparallel c5(1:2) = c5a(4:5) IF (c5a(3:3) == '\') THEN c5(3:3) = '/' ELSE ! c5a(3:3) == '/' c5(3:3) = '\' END IF c5(4:5) = c5a(1:2) velocity_azimuth_integer = velocity_azimuth_int - 180 IF (velocity_azimuth_integer < 0) velocity_azimuth_integer = velocity_azimuth_integer + 360 END IF ! parallel or antiparallel END IF ! new best separation was found END IF ! a SUB step was read END DO stepping2 CLOSE(7) ! PB2002_steps_pathfile !IMPORTANT COMMENT: ! Values defined by the block of code above ! (searching PB2002_steps_pathfile for the closest SUB step) ! include: ! *velocity_azimuth_integer ! *c5 ! and in order to run test cases with artificial SUB boundaries, ! it will be necessary to add kludgy statement overriding these values! ! For example, to use Pacific10.feg and f_SUB.nki in testing, ! so that an E-dipping trench runs S across the central Pacific, use: ! velocity_azimuth_integer = 270 ! for left plate w.r.t. right plate ! c5 = "CT/SF" ! so left plate is overriding, as you go S !NOTE: Following block of code is copied from PB2002_steps seismicity section far above: !define lane (parallelogram of great circle arcs) in counterclockwise direction: IF (c5(3:3) == '/') THEN uvec1 = begin_uvec uvec2 = end_uvec ELSE ! (c5(i)(3:3) == '\') uvec1 = end_uvec uvec2 = begin_uvec END IF IF (c5(3:3) == '\') THEN ! left plate is subducting subduction_azimuth_1 = velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec1 subduction_azimuth_2 = velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec2 ELSE ! right plate is subducting subduction_azimuth_1 = Pi + velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec1 subduction_azimuth_2 = Pi + velocity_azimuth_integer * radians_per_degree ! approximation of value at uvec2 END IF !find corner points of this lane: CALL Turn_To (azimuth_radians = subduction_azimuth_1 + Pi, & & base_uvec = uvec1, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c1_uvec) lane_uvecs(1:3, 1) = c1_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2 + Pi, & & base_uvec = uvec2, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c2_uvec) lane_uvecs(1:3, 2) = c2_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2, & & base_uvec = uvec2, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c3_uvec) lane_uvecs(1:3, 3) = c3_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_1, & & base_uvec = uvec1, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c4_uvec) lane_uvecs(1:3, 4) = c4_uvec(1:3) lane_uvecs(1:3, 5) = lane_uvecs(1:3, 1) lane_width_km = length_km * ABS(SIN(subduction_azimuth_1 - (azimuth_integer * radians_per_degree))) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" basis_length_km = lane_width_km !Compute Y_factor_of_step, based on length of NeoKinema segment: temp_integral = 0.0 s1 = 0.0 - 3.0 * alongStep_sigma_km s2 = basis_length_km + 3.0 * alongStep_sigma_km ds = (s2 - s1) / 100 DO n = 0, 100 s_km = s1 + (n / 100.0) * (s2 - s1) temp_integral = temp_integral + ds * ANORDF(s_km / alongStep_sigma_km) * ANORDF((basis_length_km - s_km) / alongStep_sigma_km) END DO Y_factor_of_step = 1 / temp_integral !seismicity rate of this step, in EQs per second, above the CMT threshold magnitude: v_oblique_mmpa = segment_slip_rate_mmpa(segment) v_oblique_mps = v_oblique_mmpa * (0.001 / s_per_year) moment_rate_of_step_SI = T5_coupledThickness_m(T5_column) * & ! coupled lithosphere thickness, in m & T5_supplement_seismicity_fraction_based_on_z(T5_column) * & & T5_assumed_mu_Pa(T5_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/SIN(T5_assumed_dip_radians(T5_column))) * & ! csc(theta) & (length_km * 1000.) ! along-strike length, in m seismicity_at_CMT_threshold = T5_CMT_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) !correct seismicity to desired threshold moment/magnitude, !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) DO threshold_index = 1, top_threshold_index IF (format_selector == 2) desired_threshold_moment_Nm = RELM_threshold_moment_Nm(threshold_index) 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 EQs_per_s(threshold_index) = G_function * seismicity_at_CMT_threshold END DO !compare B function of this step to all subsample_uvec's for possible overlap: DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat !t_lat = MIN(90.0, MAX(-90.0, t_lat)) DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) !reject points that are very far away: dot_test = uvec(1) * center_uvec(1) + uvec(2) * center_uvec(2) + uvec(3) * center_uvec(3) IF (dot_test > 0.987) THEN ! within ~1000 km; more precise distance will be computed !distances to end points of this step: uvec1 = begin_uvec to_uvec1_radians = Arc(uvec, uvec1) uvec2 = end_uvec to_uvec2_radians = Arc(uvec, uvec2) !============================================================================== !|| NOTE important difference from EQ_classification_II.f90: || !|| When choosing earthquakes to associate with SUB zones, landward and || !|| seaward half-widths (of 220 and 135 km, respectively) were strict limits.|| !|| However, when projecting long term seismicity (here), the tails of the || !|| distribution are restored. This is done by omitting the part of the test|| !|| for the LOGICAL value "proximal" that checks whether the test point is || !|| within the ends of the lanes. || !============================================================================== proximal = .TRUE. ! unless (usually) negated below: DO n = 2, 4, 2 ! sides ! Note: In EQ_classification_II.f90, this was "DO n = 1, 4" uvec1(1:3) = lane_uvecs(1:3, n) uvec2(1:3) = lane_uvecs(1:3, n+1) CALL Cross (uvec1, uvec2, tvec) ! direction toward pole is direction into the lane CALL Make_Uvec(tvec, inward_uvec) scalar_product = Dot(uvec, inward_uvec) dot_side(n) = scalar_product proximal = proximal .AND. (scalar_product > -0.01) ! 63.7 km tolerance END DO ! n = 1, 4 sides ! Note: EQ is inside the lane if all COMPUTED values in dot_side(1:4) are positive; ! these values are sines of angles from the sides, measured inwards. IF (proximal) THEN ! compute + or - distance/offset from trench, in km (even if slightly to one side of the lane): IF (c5(3:3) == '/') THEN uvec1 = begin_uvec uvec2 = end_uvec ELSE ! (c5(3:3) == '\') uvec1 = end_uvec uvec2 = begin_uvec END IF fraction_toward_uvec1 = dot_side(2) / (dot_side(2) + dot_side(4)) ! where all dot's are positive if inside; ! note that denominator is positive even if EQ (slightly) outside lane. fraction_toward_uvec1 = MAX(0.0, MIN(1.0, fraction_toward_uvec1)) ! This kludge is necessary to prevent pathalogical ! behavior when SUB steps are very short and oblique. tvec(1:3) = fraction_toward_uvec1 * uvec1(1:3) + & & (1.0 - fraction_toward_uvec1) * uvec2(1:3) CALL Make_Uvec(tvec, trench_uvec) ! at same relative left-right position in this lane distance_km = radius_km * Arc(uvec, trench_uvec) ! always zero or positive; if postive, must apply sign now: IF (distance_km > 0.0) THEN from_trench_azimuth_radians = Relative_Compass (from_uvec = trench_uvec, to_uvec = uvec) IF (COS(from_trench_azimuth_radians - subduction_azimuth_1) < 0.0) THEN ! test point is on seaward side of trench; convert distance to negative: distance_km = -distance_km END IF END IF !Note that "distance" is positive on landward side of trench, but negative on seaward side. offset_km = distance_km - SUB_peakSeismicity_km ! preserving possibility of negative values; just shifting origin !-------------------------------------------------------------------------------------------------------- ! X factor = along-lane component of spatial PDF function B: IF (offset_km >= 0.0) THEN ! landward side of seismicity peak: argument = -0.5 * (offset_km / SUB_landward_sigma_km)**2 IF (argument > -74.0) THEN X = (1.0 / 0.95) * & & (2 * SUB_landward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_landward_sigma_km)) * & & EXP(argument) ELSE X = 0.0 END IF ELSE ! seaward side of seismicity peak (but may not be seaward of the trench!) argument = -0.5 * (offset_km / SUB_seaward_sigma_km)**2 IF (argument > -74.0) THEN X = (1.0 / 0.95) * & & (2 * SUB_seaward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_seaward_sigma_km )) * & & EXP(argument) ELSE X = 0.0 END IF END IF ! which side of seismicity peak? !-------------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------------- ! Y factor = cross-lane (orthogonal to X factor) component of spatial PDF function B: ! To evaluate Y(s_km), note that length_km(i) is not the basis length that we want to use, ! but rather the lane width: lane_width_km = radius_km * (dot_side(2) + dot_side(4)) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" s_km = radius_km * dot_side(4) ! distance (if +) into lane, starting on the uvec1 end IF ((s_km > (-3.0 * alongStep_sigma_km)).AND.(s_km < (lane_width_km + 3.0 * alongStep_sigma_km))) THEN ! use complex formula: Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((lane_width_km - s_km) / alongStep_sigma_km) ELSE ! too far away; don't risk numerical problems within ANORDF: Y = 0.0 END IF !-------------------------------------------------------------------------------------------------------- ELSE ! not proximal; punt X = 0.0 Y = 0.0 END IF B = X * Y ! spatial PDF is a product of cross-trace function and along-trace function !At this point, B is in units of (km)**(-2). B = B * 1.0D-6 !Now, B is in units of m**(-2), which is SI. IF (B > 0.0D0) THEN DO threshold_index = 1, top_threshold_index IF (segment_creeping(segment)) THEN ! subduction seismicity is purely due to plate-bending: seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + & & SUB_bending_fraction * weight * B * EQs_per_s(threshold_index) ELSE ! normal, mostly-coupled subduction: seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + & & weight * B * EQs_per_s(threshold_index) END IF END DO END IF ! B > 0.0D0 END IF ! uvec is close enough to center_uvec for further testing END DO ! jj = 1, subsamplings END DO ! ii = 1, subsamplings END DO ! j = 1, columns END DO ! i = 1, rows ELSE ! not a subduction zone ('S') segment !N.B. begin_uvec, end_uvec, center_uvec, length_km, azimuth_integer, segment_vec, ! segment_azimuth_radians(segment), ! were already defined above for this segment. !Find closest grid point to segment and decide if it is oceanic or not: CALL Uvec_2_LonLat(center_uvec, lon, lat) ! but, lon may be on the wrong cycle! delta_lon_from_grid_center = lon - (lon_min + lon_max)/2.0 IF (delta_lon_from_grid_center < -180.) delta_lon_from_grid_center = delta_lon_from_grid_center + 360. IF (delta_lon_from_grid_center > +180.) delta_lon_from_grid_center = delta_lon_from_grid_center - 360. lon = delta_lon_from_grid_center + (lon_min + lon_max)/2.0 ! now on the correct cycle j = 1 + NINT((lon - lon_min) / d_lon) j = MAX(MIN(j, columns), 1) i = 1 + NINT((lat_max - lat) / d_lat) i = MAX(MIN(i, rows), 1) oceanic = .NOT.continental(i, j) !Choose which of the plate boundary types to calibrate to, !and which fault dip to assume (consistent with NeoKinema): c1 = segment_c6(segment)(6:6) IF ((c1 == 'L').OR.(c1 == 'R').OR.(c1 == 'l').OR.(c1 == 'r')) THEN IF (oceanic) THEN k_class_index = 5 ! OTF IF (segment_slip_rate_mmpa(segment) < 39.5) THEN T5_column = 7 ! slow OTF ELSE IF (segment_slip_rate_mmpa(segment) < 68.5) THEN T5_column = 8 ! medium OTF ELSE T5_column = 9 ! fast OTF END IF ELSE ! continental k_class_index = 2 ! CTF T5_column = 2 ! CTF END IF ELSE IF ((c1 == 'N').OR.(c1 == 'D').OR.(c1 == 'n').OR.(c1 == 'd')) THEN IF (oceanic) THEN k_class_index = 4 ! OSR T5_column = 5 ! OSR normal !NB: However, coupled thickness will be based on equation of Bird et al. [2002]. ELSE ! continental k_class_index = 3 ! CRB T5_column = 1 ! CRB END IF ELSE IF ((c1 == 'T').OR.(c1 == 'P').OR.(c1 == 't').OR.(c1 == 'p')) THEN IF (oceanic) THEN k_class_index = 6 ! OCB T5_column = 10 ! OCB ELSE ! continental k_class_index = 1 ! CCB IF (segment_slip_rate_mmpa(segment) <= 24.2) THEN T5_column = 3 ! slow CCB ELSE T5_column = 4 ! fast CCB END IF END IF ELSE WRITE (*, *) WRITE (*, "(' ERROR: Illegal byte ',A,' in fault descriptor ',A)") c1, segment_c6(segment) CALL Pause() STOP END IF fault_dip_degrees = T5_assumed_dip_degrees(T5_column) fault_dip_radians = T5_assumed_dip_radians(T5_column) lithosphere_thickness_km = T5_assumed_lithosphere_km(T5_column) ! used to compute width of parallelogram !Locate corner points of seismicity parallelogram (on planet's surface): cross_trace_km = lithosphere_thickness_km / TAN(fault_dip_radians) IF ((c1 == 'L').OR.(c1 == 'R').OR.(c1 == 'l').OR.(c1 == 'r')) THEN !dip sense is uncertain, so extend parallelograms to BOTH sides of the trace: CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) + Pi_over_2, & & base_uvec = begin_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec1) CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) + Pi_over_2, & & base_uvec = end_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec2) CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) - Pi_over_2, & & base_uvec = end_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec3) CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) - Pi_over_2, & & base_uvec = begin_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec4) parallelogram_area_m2 = 1.0E6 * length_km * 2.0 * cross_trace_km * & & ABS(COS(segment_azimuth_radians(segment) - (azimuth_integer * radians_per_degree))) !Note: segment_azimuth_radians(segment) is the mean azimuth of the whole fault; ! while azimuth_integer*radians_per_degree is the particular azimuth of this one segment. ELSE ! dip-slip fault; dip is to left of trace, by NeoKinema f.dig convention. uvec1 = begin_uvec uvec2 = end_uvec CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) - Pi_over_2, & & base_uvec = end_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec3) CALL Turn_To (azimuth_radians = segment_azimuth_radians(segment) - Pi_over_2, & & base_uvec = begin_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec4) parallelogram_area_m2 = 1.0E6 * length_km * cross_trace_km * & & ABS(COS(segment_azimuth_radians(segment) - (azimuth_integer * radians_per_degree))) !Note: segment_azimuth_radians(segment) is the mean azimuth of the whole fault; ! while azimuth_integer*radians_per_degree is the particular azimuth of this one segment. END IF !GPBdebug: Use the following 10 lines to output the outlines of the fault patches: !CALL Uvec_2_LonLat(uvec1, lon1, lat1) !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon1, lat1 !CALL Uvec_2_LonLat(uvec2, lon, lat) !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon, lat !CALL Uvec_2_LonLat(uvec3, lon, lat) !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon, lat !CALL Uvec_2_LonLat(uvec4, lon, lat) !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon, lat !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon1, lat1 !WRITE (72, "('*** end of line segment ***')") !determine coupled thickness of lithosphere (from Table 5 of Bird & Kagan [2004]; OSR includes later cz(v) model): IF (T5_column == 5) THEN ! OSR normal d0_m = 1480. ! from CMT_PB2002_OSR_pure_normal_seismicity.xls, based on Bird & Kagan [2004] vs_mmpa = 19. ! [ibid] coupledThickness_m = d0_m * EXP(-segment_heave_rate_mmpa(segment) / vs_mmpa) ELSE ! all other classes coupledThickness_m = T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) END IF !seismicity rate of this step, in EQs per second, above the CMT threshold magnitude: v_oblique_mmpa = segment_slip_rate_mmpa(segment) v_oblique_mps = v_oblique_mmpa * (0.001 / s_per_year) moment_rate_of_step_SI = coupledThickness_m * & ! coupled lithosphere thickness, in m & T5_assumed_mu_Pa(T5_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/SIN(T5_assumed_dip_radians(T5_column))) * & ! csc(theta) & (length_km * 1000.) ! along-strike length, in m seismicity_at_CMT_threshold = T5_CMT_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) !correct seismicity to desired threshold moment/magnitude, !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) DO threshold_index = 1, top_threshold_index IF (format_selector == 2) desired_threshold_moment_Nm = RELM_threshold_moment_Nm(threshold_index) 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 EQs_per_s(threshold_index) = G_function * seismicity_at_CMT_threshold !distribute seismicity evenly over the parallelogram: EQs_per_s_per_m2(threshold_index) = EQs_per_s(threshold_index) / parallelogram_area_m2 END DO !prepare 4 pole_uvec's (poles of short arcs from uvec1->uvec2, .... uvec4-->uvec1: CALL Cross(uvec1, uvec2, tvec) CALL Make_Uvec(tvec, pole_uvec1) CALL Cross(uvec2, uvec3, tvec) CALL Make_Uvec(tvec, pole_uvec2) CALL Cross(uvec3, uvec4, tvec) CALL Make_Uvec(tvec, pole_uvec3) CALL Cross(uvec4, uvec1, tvec) CALL Make_Uvec(tvec, pole_uvec4) !compare area of this parallelogram to all subsample_uvec's for possible overlap: DO i = 1, rows !latitude = lat_max - (i - 1) * d_lat DO j = 1, columns !longitude = lon_min + (j - 1) * d_lon DO ii = 1, subsamplings !t_lat = latitude + lat_offset + ii * dd_lat !t_lat = MIN(90.0, MAX(-90.0, t_lat)) DO jj = 1, subsamplings !t_lon = longitude + lon_offset + jj * dd_lon !CALL LonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) !reject points that are very far away, by testing with center_uvec (center of segment): dot_test = uvec(1) * center_uvec(1) + uvec(2) * center_uvec(2) + uvec(3) * center_uvec(3) IF (dot_test > 0.9967) THEN ! within ~500 km; more precise distance will be computed !check whether uvec is on correct side of all 4 sides of parallelogram: dot_test = uvec(1) * pole_uvec1(1) + uvec(2) * pole_uvec1(2) + uvec(3) * pole_uvec1(3) IF (dot_test > 0.0) THEN dot_test = uvec(1) * pole_uvec2(1) + uvec(2) * pole_uvec2(2) + uvec(3) * pole_uvec2(3) IF (dot_test > 0.0) THEN dot_test = uvec(1) * pole_uvec3(1) + uvec(2) * pole_uvec3(2) + uvec(3) * pole_uvec3(3) IF (dot_test > 0.0) THEN dot_test = uvec(1) * pole_uvec4(1) + uvec(2) * pole_uvec4(2) + uvec(3) * pole_uvec4(3) IF (dot_test >= 0.0) THEN inside = .TRUE. ELSE inside = .FALSE. END IF ELSE inside = .FALSE. END IF ELSE inside = .FALSE. END IF ELSE inside = .FALSE. END IF !==================================================================== IF (inside) THEN DO threshold_index = 1, top_threshold_index seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + & & weight * EQs_per_s_per_m2(threshold_index) END DO END IF ! B > 0.0D0 !==================================================================== END IF ! uvec is close enough to parallelogram for further testing END DO ! jj = 1, subsamplings END DO ! ii = 1, subsamplings END DO ! j = 1, columns END DO ! i = 1, rows END IF ! a subduction zone ('S') segment, or not END IF ! ((segment_slip_rate_mmpa(segment) > 0.0).AND.(.NOT.segment_creeping(segment))) END DO add_segment_seismicity ! loop on segments (parallel to lines in h_nko_pathfile, but memorized) WRITE (*,"('+Adding seismicity of NeoKinema fault segments...DONE ')") !GPBdebug: Use line below to close file with outlines of fault patches: !CLOSE(72) !============================================================================================================= DEALLOCATE ( segment_c6 ) DEALLOCATE ( segment_uvecs ) DEALLOCATE ( segment_creeping ) DEALLOCATE ( segment_slip_rate_mmpa ) DEALLOCATE ( segment_heave_rate_mmpa ) DEALLOCATE ( segment_used ) DEALLOCATE ( segment_outer_uvecs ) DEALLOCATE ( segment_azimuth_radians ) END DO ! orogen = 1, orogens END IF ! orogens > 0 DEALLOCATE ( subsample_uvec ) number_of_zeros = 0 ! initializing IF (format_selector == 1) THEN WRITE (1, "(1P,10D10.2)") ((seismicity(i, j, 1), j = 1, columns), i = 1, rows) !Note: Line breaks in output file have no relation to the logical shape of this array. ! Programs that read .grd files must take this into account, by reading the whole ! array with a single READ (as it was written), and not row-by-row. ELSE IF (format_selector == 2) THEN DO i = rows, 1, -1 ! working from S toward N center_latitude = lat_max - (i - 1) * d_lat lowest_latitude = center_latitude - 0.5 * d_lat highest_latitude = center_latitude + 0.5 * d_lat dy_meters = Earth_radius_m * (highest_latitude - lowest_latitude) / 57.29578 DO j = 1, columns ! working from W to E IF (template_mask(i, j)) THEN center_longitude = lon_min + (j - 1) * d_lon lowest_longitude = center_longitude - 0.5 * d_lon highest_longitude = center_longitude + 0.5 * d_lon dx_meters = Earth_radius_m * ((highest_longitude - lowest_longitude) / 57.29578) * & & COS (center_latitude / 57.29578) DO threshold_index = 1, 41 ! from low to high lowest_magnitude = RELM_threshold_magnitude(threshold_index) IF (threshold_index < 41) THEN highest_magnitude = RELM_threshold_magnitude(1 + threshold_index) EQs_per_s_per_m2(threshold_index) = seismicity(i, j, threshold_index) - & & seismicity(i, j, (1 + threshold_index)) ELSE ! effectively, no upper limit highest_magnitude = 10. EQs_per_s_per_m2(threshold_index) = seismicity(i, j, threshold_index) END IF EQs_per_s(threshold_index) = EQs_per_s_per_m2(threshold_index) * dx_meters * dy_meters expectation = EQs_per_s(threshold_index) * 5. * 3.15576E7 ! converting from 1 s to 5 years IF (expectation <= 0.0D0) number_of_zeros = number_of_zeros + 1 WRITE (1, 2999) lowest_longitude, tab, highest_longitude, tab, & & lowest_latitude, tab, highest_latitude, tab, & & smaller_depth_integer, tab, larger_depth_integer, tab, & & lowest_magnitude, tab, highest_magnitude, tab, & & expectation, tab ! finally, the mask "1" is supplied by the FORMAT 2999 FORMAT (F12.4,A,F12.4,A,F12.4,A,F12.4,A,I2,A,I2,A,F6.2,A,F6.2,A,ES10.2,A,'1') END DO ! threshold_index END IF ! template_mask(i, j) END DO ! columns END DO ! rows WRITE (1, "('end_forecast')") END IF IF (number_of_zeros > 0) THEN WRITE (*, "(/' WARNING! ',I10,' of the values printed were zeros!!!')") number_of_zeros CALL Pause() END IF !Integrate seismicity over grid area: integral = 0.0D0 ! whole array, if format_selector == 2 DO i = 1, rows latitude = lat_max - (i - 1) * d_lat COS_latitude_radians = COS(latitude * radians_per_degree) IF ((i == 1).OR.(i == rows)) THEN dy = 0.5D0 * d_lat * radians_per_degree * Earth_radius_m ELSE dy = d_lat * radians_per_degree * Earth_radius_m END IF DO j = 1, columns IF (format_selector == 2) THEN include_in_integral = template_mask(i, j) ELSE ! template_mask never allocated; don't refer to it! include_in_integral = .TRUE. END IF IF (include_in_integral) THEN IF ((j == 1).OR.(j == columns)) THEN dx = 0.5D0 * d_lon * COS_latitude_radians * radians_per_degree * Earth_radius_m ELSE dx = d_lon * COS_latitude_radians * radians_per_degree * Earth_radius_m END IF DO threshold_index = 1, top_threshold_index dI = seismicity(i, j, threshold_index) * dx * dy integral(threshold_index) = integral(threshold_index) + dI END DO END IF END DO END DO integral_converted = integral * 100.0 * 365.25 * 24 * 60 * 60 ! whole array, if format_selector == 2 IF (format_selector == 1) THEN WRITE (*, "(/' Area integral of seismicity =' & & /' ',ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted WRITE (1, "( 'Area integral of seismicity =' & & /ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") integral, integral_converted ELSE IF (format_selector == 2) THEN WRITE (*, "(/' Area integral (in template area) of seismicity (earthquakes/century)')") WRITE (*, "( ' ====================================================================')") WRITE (*, "( ' Threshold magnitude Earthquakes above this magnitude')") DO threshold_index = 1, top_threshold_index WRITE (*, "( ' ',F12.3,11X,ES12.3)") & & RELM_threshold_magnitude(threshold_index), integral_converted(threshold_index) END DO END IF CLOSE (1) ! output .grd or .dat file WRITE (*, "(/' All computations completed. See file ',A)") TRIM(output_grd_pathfile) CALL Pause() CONTAINS ! ************************************************************************ REAL FUNCTION ATAN2F(y, x) IMPLICIT NONE REAL, INTENT(IN) :: y, x IF ((y /= 0.0).OR.(x /= 0.0)) THEN ATAN2F = ATAN2(y, x) ELSE ATAN2F = 0.0 END IF END FUNCTION ATAN2F SUBROUTINE Del_Gjxy_del_thetaphi (l_, uvec1, uvec2, uvec3, r_, dG) INTEGER, INTENT(IN) :: l_ ! element number or code REAL, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & & r_ ! position vector REAL, DIMENSION (3,2,2,2), INTENT(OUT) :: dG ! Computes array of 2 derivitives of each of the 2 components of ! each of the 6 nodal functions for element l_, ! whose corners are uvec1, uvec2, and uvec3 (counterclockwise!), ! evaluated at position r_ (Cartesian unit vector). ! Results are in 1./radian (dimensionless), NOT 1./degree ! ! Also note that derivitives with respect to phi are algebraic, ! not spatial; that is, they have not been corrected for the ! changing length of a unit of phi with latitude, so these ! derivitives tend to become small near the poles, where ! a unit of phi does not cover very much distance. ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = 0 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2)? INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function? INTEGER :: m ! 1:2 = theta (S-ward) or phi (E-ward) derivitive? REAL, DIMENSION(3,2) :: del_r_ ! theta- and phi-derivitives of r_ (in 3-D) REAL, DIMENSION(3,2) :: local ! local Theta, Phi unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,2,2) :: del_local ! theta-, phi- derivitives of local REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, tv3, vfa, vfb ! temporary vector factors REAL :: cos_phi, cos_theta, phi, sin_phi, sin_theta INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ corner(1:3, 1) = uvec1(1:3) corner(1:3, 2) = uvec2(1:3) corner(1:3, 3) = uvec3(1:3) DO j = 1, 3 tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin calculations which depend on r_ CALL Local_Theta(r_, tv) local(1:3,1) = tv CALL Local_Phi(r_, tv) local(1:3,2) = tv ! Note: these functions will catch polar points; don't test again phi = ATAN2(r_(2), r_(1)) cos_phi = COS(phi) sin_phi = SIN(phi) cos_theta = r_(3) sin_theta = SQRT(r_(1)**2 + r_(2)**2) del_r_(1:3,1) = local(1:3,1) ! d.r_/d.theta = Theta del_r_(1:3,2) = local(1:3,2) * sin_theta ! d.r_/d.phi = Phi * SIN(theta) del_local(1:3,1,1) = - r_(1:3) ! d.Theta/d.theta = - r_ del_local(1:3,1,2) = local(1:3,2) * cos_theta ! d.Theta/d.phi = Phi * COS(theta) del_local(1:3,2,1) = (/ 0., 0., 0. /) ! d.Phi/d.theta = 0 del_local(1:3,2,2) = (/ -cos_phi, -sin_phi, 0. /) ! d.Phi/d.phi = (-COS(phi),-SIN(phi,0) DO j = 1, 3 ! 3 corner nodes of element i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vfa) vfb = vfa / Dot (tv1, vfa) DO x = 1, 2 ! unit velocity at node is S or E DO y = 1, 2 ! S- or E- component of nodal function tv1 = post(1:3, j, x) tvi = local(1:3, y) DO m = 1, 2 ! theta- or phi-derivitive tv = del_r_(1:3, m) tvo = del_local(1:3, y, m) dG(j, x, y, m) = & & (Dot(tv,vfb)*Dot(tv1,tvi)) + & & (Dot(r_,vfb)*Dot(tv1,tvo)) END DO END DO END DO END DO END SUBROUTINE Del_Gjxy_del_thetaphi SUBROUTINE Dumb_s123 (element, vector, node, xyz_nod, center, a_, & & s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). IMPLICIT NONE INTEGER, INTENT(IN) :: element ! element # REAL, DIMENSION(3), INTENT(IN) :: vector ! uvec to point INTEGER, DIMENSION(:,:), INTENT(IN) :: node ! element definitions REAL, DIMENSION(:,:), INTENT(IN) :: xyz_nod ! uvecs of nodes REAL, DIMENSION(:,:), INTENT(IN) :: center ! uvec of each element (uvec) REAL, DIMENSION(:), INTENT(IN) :: a_ ! element areas (plane; R == 1.0) REAL, INTENT(OUT) :: s1, s2, s3 ! results !- - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: i1, i2, i3 REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL :: d1, dc, t IF (element == 0) THEN WRITE (*,"(' ERROR: element = 0 in Dumb_s123')") CALL Traceback END IF i1 = node(1, element) i2 = node(2, element) i3 = node(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = DOT_PRODUCT(vector, tv1) IF (dc <= 0.) THEN WRITE (*,"(' ERROR: Internal vector >= 90 deg. from element in Dumb_s123')") CALL Traceback END IF tv2 = xyz_nod(1:3, i1) d1 = DOT_PRODUCT(tv2, tv1) t = d1 / dc v1 = t * vector tvi = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) tvo = v1(1:3) - xyz_nod(1:3,i3) CALL Cross(tvi, tvo, tv) s1 = 0.5 * DOT_PRODUCT(tv1, tv) / a_(element) tvi = xyz_nod(1:3,i1) - xyz_nod(1:3,i3) tvo = v1(1:3) - xyz_nod(1:3,i1) CALL Cross(tvi, tvo, tv) s2 = 0.5 * DOT_PRODUCT(tv1, tv) / a_(element) s3 = 1.00 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE E_rate(R, l_, nodes, G, dG, theta_, vw, eps_dot) ! evaluate strain-rate in spherical continuum element REAL, INTENT(IN) :: R ! radius of planet, in m INTEGER, INTENT(IN) :: l_ ! element number INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes REAL, DIMENSION(3,2,2) :: G ! nodal functions @ selected point REAL, DIMENSION(3,2,2,2):: dG ! derivitives of nodal functions REAL, INTENT(IN) :: theta_ ! colatitude, radians REAL, DIMENSION(:), INTENT(IN) :: vw REAL, DIMENSION(3), INTENT(OUT) :: eps_dot INTEGER :: iv, iw, j REAL :: cott, csct, prefix eps_dot = 0. ! (1..3) cott = 1. / TAN(theta_) csct = 1. / SIN(theta_) prefix = 1./R DO j = 1, 3 iv = 2 * nodes(j, l_) - 1 iw = iv + 1 ! epsilon_dot_sub_theta_theta eps_dot(1) = eps_dot(1) + & & vw(iv) * prefix * dG(j,1,1,1) + & & vw(iw) * prefix * dG(j,2,1,1) ! epsilon_dot_sub_theta_phi eps_dot(2) = eps_dot(2) + & & vw(iv) * prefix * 0.5 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) + & & vw(iw) * prefix * 0.5 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) ! epsilon_dot_sub_phi_phi eps_dot(3) = eps_dot(3) + & & vw(iv) * prefix * (csct * dG(j,1,2,2) + cott * G(j,1,1)) + & & vw(iw) * prefix * (csct * dG(j,2,2,2) + cott * G(j,2,1)) END DO ! 3 local nodes END SUBROUTINE E_rate SUBROUTINE Gjxy (l_, uvec1, uvec2, uvec3, r_, G) INTEGER, INTENT(IN) :: l_ ! element number or code REAL, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & ! corners & r_ ! position vector REAL, DIMENSION (3,2,2), INTENT(OUT) :: G ! Computes matrix of 6 vector nodal functions for element l_, ! whose corners are at uvec1, uvec2, and uvec3 (counterclockwise!), ! evaluated at position r_ (Cartesian unit vector). ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = -999 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2) INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function REAL, DIMENSION(3,2) :: local ! local unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tvi, tvo, tv1, tv2, tv3, vf ! temporary vector factor REAL :: f_sup_j ! as in Kong and Bird (1995) [ j == k ] INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ corner(1:3, 1) = uvec1(1:3) corner(1:3, 2) = uvec2(1:3) corner(1:3, 3) = uvec3(1:3) DO j = 1, 3 tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin computations which depend on r_ CALL Local_Theta(r_, tvo) local(1:3,1) = tvo CALL Local_Phi(r_, tvo) local(1:3,2) = tvo DO j = 1, 3 i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vf) f_sup_j = Dot(r_, vf) / Dot (tv1, vf) DO x = 1, 2 tv1 = post(1:3, j, x) DO y = 1, 2 tv2 = local(1:3, y) G(j, x, y) = f_sup_j * Dot(tv1, tv2) END DO END DO END DO END SUBROUTINE Gjxy INTEGER FUNCTION IBELOW (X) IMPLICIT NONE REAL, INTENT(IN) :: X IBELOW=INT(X) IF (X.LT.(1.*IBELOW)) IBELOW=IBELOW-1 RETURN END FUNCTION IBELOW SUBROUTINE Learn_Spherical_Triangles (numel, nodes, node_uvec, chatty, & & a_, center, neighbor) !Creates arrays needed by lookup subr. Which_Spherical_Triangle: ! a_ = area of plane triangle below element, when radius == 1.0 ! center = uvec pointing to center of element ! neighbor = neighboring spherical triangular elements on each ! side (or zero if none); note that algorithm ! depends on node-location match, not on node-number ! match, and therefore ignores intevening faults. !These arrays are only meaningful for finite element grids used ! with SHELLS and/or RESTORE. IMPLICIT NONE INTEGER, INTENT(IN) :: numel ! number of spherical triangle elements INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes ! element definitions REAL, DIMENSION(:,:), INTENT(IN) :: node_uvec ! uvecs of nodes LOGICAL, INTENT(IN) :: chatty REAL, DIMENSION(:), INTENT(OUT) :: a_ REAL, DIMENSION(:,:), INTENT(OUT) :: center INTEGER, DIMENSION(:,:), INTENT(OUT) :: neighbor INTEGER :: furthest, i, ia, ib, i1, i2, i3, j, j1, j2, j3, k, l_, m, step_aside REAL, DIMENSION(3) :: a, b, c, t, u IF (chatty) WRITE (*,"(' Learning the spherical triangles...')") furthest = (numel + 1) / 2 neighbor = 0 ! whole array, initialized to "no neighbor on this side" homes: DO l_ = 1, numel !first, a_ i1 = nodes(1,l_) i2 = nodes(2,l_) i3 = nodes(3,l_) a = node_uvec(1:3,i2) - node_uvec(1:3,i1) b = node_uvec(1:3,i3) - node_uvec(1:3,i2) CALL Cross (a, b, c) a_(l_) = 0.5 * Magnitude(c) !second, compute center t(1:3) = (node_uvec(1:3,i1)+node_uvec(1:3,i2)+node_uvec(1:3,i3))/3.0 CALL Make_Uvec(t, u) center(1:3, l_) = u(1:3) !third, find neighbor(?) for each side of element sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) ia = nodes(k, l_) ! 1st node along side ib = nodes(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO step_aside = 1, furthest m = l_ + step_aside ! I also try -step_aside, below m = 1 + MOD(m-1, numel) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (Same(j1, ib) .AND. Same(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j2, ib) .AND. Same(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j3, ib) .AND. Same(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF m = l_ - step_aside ! I also try +step_aside, above m = 1 + MOD(m-1+numel, numel) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (Same(j1, ib) .AND. Same(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j2, ib) .AND. Same(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j3, ib) .AND. Same(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides IF (chatty) WRITE (*,"('+Learning the spherical triangles...',I6)") l_ END DO homes IF (chatty) WRITE (*,"('+Learning the spherical triangles...DONE ')") END SUBROUTINE Learn_Spherical_Triangles REAL FUNCTION LookUp_ETOPO5 (longitude, latitude) IMPLICIT NONE REAL, INTENT(IN) :: longitude, latitude REAL :: lat, lon, & & rRow, rColumn, & & fracRow, fracColumn, & & lower_left, lower_right, & & upper_left, upper_right, base, top INTEGER :: iRow1, iRow2, jColumn1, jColumn2 lon = longitude IF (lon >= 360.0) lon = lon - 360.0 IF (lon < 0.0) lon = lon + 360.0 lat = MIN(90.0, MAX(-90.0, latitude)) rRow = 12.0 * lat rColumn = 12.0 * lon iRow1 = IBelow(rRow) iRow2 = iRow1 + 1 jColumn1 = IBelow(rColumn) jColumn2 = jColumn1 + 1 fracRow = rRow - iRow1 fracColumn = rColumn - jColumn1 IF (iRow1 < -1080) THEN iRow2 = iRow2 + (-1080 - iRow1) fracRow = fracRow - (-1080 - iRow1) iRow1 = -1080 END IF IF (iRow2 > +1080) THEN iRow1 = iRow1 - (iRow2 - 1080) fracRow = fracRow + (iRow2 - 1080) iRow2 = +1080 END IF IF (jColumn1 < 0) THEN jColumn2 = jColumn2 + (0 - jColumn1) fracColumn = fracColumn - (0 - jColumn1) jColumn1 = 0 END IF IF (jColumn1 == 4319) THEN jColumn2 = 0 END IF lower_left = ETOPO5(iRow1, jColumn1) lower_right = ETOPO5(iRow1, jColumn2) upper_left = ETOPO5(iRow2, jColumn1) upper_right = ETOPO5(iRow2, jColumn2) base = lower_left + fracColumn * (lower_right - lower_left) top = upper_left + fracColumn * (upper_right - upper_left) LookUp_ETOPO5 = base + fracRow * (top - base) END FUNCTION LookUp_ETOPO5 REAL FUNCTION LookUp_age_1p5 (longitude, latitude) IMPLICIT NONE REAL, INTENT(IN) :: longitude, latitude REAL :: lat, lon, & & rRow, rColumn, & & fracRow, fracColumn, & & lower_left, lower_right, & & upper_left, upper_right, base, top INTEGER :: iRow1, iRow2, jColumn1, jColumn2 lon = longitude IF (lon >= 360.0) lon = lon - 360.0 IF (lon < 0.0) lon = lon + 360.0 lat = MIN(90.0, MAX(-90.0, latitude)) rRow = 10.0 * lat rColumn = 10.0 * lon iRow1 = IBelow(rRow) iRow2 = iRow1 + 1 jColumn1 = IBelow(rColumn) jColumn2 = jColumn1 + 1 fracRow = rRow - iRow1 fracColumn = rColumn - jColumn1 IF (iRow1 < -720) THEN LookUp_age_1p5 = 220.0 ! unknown RETURN END IF IF (iRow2 > +900) THEN iRow1 = iRow1 - (iRow2 - 900) fracRow = fracRow + (iRow2 - 900) iRow2 = +900 END IF IF (jColumn1 < 0) THEN jColumn2 = jColumn2 + (0 - jColumn1) fracColumn = fracColumn - (0 - jColumn1) jColumn1 = 0 END IF IF (jColumn2 > 3600) THEN jColumn1 = jColumn1 - (jColumn2 - 3600) fracColumn = fracColumn + (jColumn2 - 3600) jColumn2 = 3600 END IF lower_left = age_1p5(iRow1, jColumn1) lower_right = age_1P5(iRow1, jColumn2) upper_left = age_1p5(iRow2, jColumn1) upper_right = age_1p5(iRow2, jColumn2) base = lower_left + fracColumn * (lower_right - lower_left) top = upper_left + fracColumn * (upper_right - upper_left) LookUp_age_1p5 = base + fracRow * (top - base) END FUNCTION LookUp_age_1p5 REAL FUNCTION Moment(magnitude) ! Returns scalar seismic moment in N m, per ! equation (1) of Bird & Kagan [2004; Bull. Seism. Soc. Am.]; ! originally from Hanks & Kanamori [1979; J. Geophys. Res., v. 84, p. 2348-2350]. IMPLICIT NONE REAL, INTENT(IN) :: magnitude moment = 10.**((1.5 * magnitude) + 9.05) END FUNCTION Moment SUBROUTINE No_Fault_Elements_Allowed() ! Called if nfl > 0 in the .feg file: ! Prints explanatory messages and stops execution. IMPLICIT NONE 101 FORMAT (& &/' This .feg (finite element grid) file contains fault elements!'& &/' Fault elements are not allowed in NeoKinema grids, because:'& &/' I. NeoKinema does not require fault elements.'& &/' 1. NeoKinema has logic to add the compliance of any number of'& &/' faults to the continuum (triangle) elements that contain them.'& &/' 2. NeoKinema has logic to infer the heave-rate and slip-rate of'& &/' such implied fault(s) from the computed strain-rate of the'& &/' triangular continuum element(s).'& &/' 3. Graphics program NeoKineMap has logic to plot the heave-rates'& &/' of these faults, and also velocity fields with fault '& &/' disontinuities, without the use of fault elements.') 102 FORMAT (& &/' II. Fault elements cause bad grid topology.'& &/' 1. Fault elements are not read or stored by NeoKineMap.'& &/' 2. With fault elements deleted, the grids on the two sides of'& &/' each fault are not connected in any way.'& &/' 3. The solution process may fail due to a singular stiffness'& &/' matrix during solution of the linear system.'& &/' 4. Even if the solution does not fail, its physical interpretation'& &/' will be problematical.') 103 FORMAT (& &/' III. Fault elements should be eliminated from the .feg file.'& &/' 1. Re-load this .feg file into OrbWeaver, select command Faults,'& &/' and use the right mouse button to heal the fault cuts.'& &/' 2. Use command Adjust to move all nodes off of fault traces.'& &/' 3. Save the edited grid and re-number it with OrbNumber.'& &/' 4. Alternatively, use command Tile in OrbWeaver to create a'& &/' new .feg grid with no fault elements.') WRITE (*, 101) CALL Pause() WRITE (*, 102) CALL Pause() WRITE (*, 103) CALL Pause() STOP END SUBROUTINE No_Fault_Elements_Allowed SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause SUBROUTINE Principal_Axes_22 (e11, e12, e22, & & e1, e2, u1x,u1y, u2x,u2y) ! Find principal values (e1,e2) of the symmetric 2x2 tensor ! e11 e12 ! e12 e22 ! and also the associated unit eigenvectors: ! #1 = (u1x, u1y); #2 = (u2x, u2y). ! The convention is that e1 <= e2. IMPLICIT NONE REAL, INTENT(IN) :: e11, e12, e22 REAL, INTENT(OUT) :: e1, e2, u1x, u1y, u2x, u2y REAL :: c, f1, f11, f12, f2, f22, r, scale, smallest, test, theta ! Smallest number that can be squared without underflowing: smallest = 1.1 * SQRT(TINY(f1)) ! First, check for trivial isotropic case: IF ((e11 == e22).AND.(e12 == 0.)) THEN ! In this case, directions are arbitrary: e1 = e11 u1x = 1. u1y = 0. e2 = e22 u2x = 0. u2y = 1. RETURN END IF ! Else, re-scale matrix to prevent under/overflows: scale = MAX(ABS(e11), ABS(e12), ABS(e22)) f11 = e11 / scale IF (ABS(f11) < smallest) f11 = 0. f12 = e12 / scale IF (ABS(f12) < smallest) f12 = 0. f22 = e22 / scale IF (ABS(f22) < smallest) f22 = 0. ! Find eigenvalues and eigenvectors of scaled matrix: r = SQRT(((f11 - f22)*0.5)**2 + f12**2) c = (f11 + f22)*0.5 f1 = c - r f2 = c + r test = 0.01 * MAX (ABS(f1), ABS(f2)) IF ((ABS(f12) > test).OR.(ABS(f11 - f1) > test)) THEN theta = Atan2F((f11 - f1), -f12) ELSE theta = Atan2F(f12, (f1 - f22)) END IF u1x = COS(theta) u1y = SIN(theta) u2x = u1y u2y = -u1x ! Undo the scaling e1 = scale * f1 e2 = scale * f2 END SUBROUTINE Principal_Axes_22 SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. IMPLICIT NONE REAL, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.) THEN highest = MAXVAL(s) medium = 1.00 - lowest - highest IF (medium > 0.) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0. sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00 / (1.00 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0. s(sideb) = 0. END IF END IF END SUBROUTINE Pull_in LOGICAL FUNCTION Same(i,j) ! Are node_uvec #i and #j the same vector? INTEGER, INTENT(IN) :: i, j !Refers to global array REAL, DIMENSION(3, numnod):: node_uvec !The logic is: !Same = (node_uvec(1,i) == node_uvec(1,j)).AND. & ! & (node_uvec(2,i) == node_uvec(2,j)).AND. & ! & (node_uvec(3,i) == node_uvec(3,j)) !But, it is written this way for speed: IF (node_uvec(1,i) == node_uvec(1,j)) THEN IF (node_uvec(2,i) == node_uvec(2,j)) THEN IF (node_uvec(3,i) == node_uvec(3,j)) THEN Same = .TRUE. ELSE Same = .FALSE. END IF ELSE Same = .FALSE. END IF ELSE Same = .FALSE. END IF END FUNCTION Same LOGICAL FUNCTION Same_Uvec(uvec1, uvec2) REAL, DIMENSION(3), INTENT(IN) :: uvec1, uvec2 IF (uvec1(1) == uvec2(1)) THEN IF (uvec1(2) == uvec2(2)) THEN IF (uvec1(3) == uvec2(3)) THEN Same_Uvec = .TRUE. ELSE Same_Uvec = .FALSE. END IF ELSE Same_Uvec = .FALSE. END IF ELSE Same_Uvec = .FALSE. END IF END FUNCTION Same_Uvec SUBROUTINE Test_Opening_of_Existing_File(pathfile) IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: pathfile INTEGER :: ios OPEN (UNIT = 123, FILE = pathfile, STATUS = "OLD", DISP = "KEEP", IOSTAT = ios) IF (ios /= 0) THEN WRITE(*) WRITE(*, "(' ERROR: Input file ',A,' could not be found. ios = ',I6)") TRIM(pathfile), ios CALL Pause() STOP END IF CLOSE (123) END SUBROUTINE Test_Opening_of_Existing_File SUBROUTINE Which_Spherical_Triangle (b_, cold_start, & & num_ele, node, xyz_node, center, a_, neighbor, & & success, iele, s1, s2, s3) !Locates a point (b_, a uvec) in element iele with internal !coordinates (s1, s2, s3) in a SHELLS or RESTORE .feg !and reports success. !If (cold_start), makes no use of input iele, s1, s2, s3 !If not, uses these values to initialize the search. ! !Note that Learn_Spherical_Triangles can be used to initialize !necessary arrays a_, center, and neighbor. ! !Beware of variable name changes (numel = num_ele, nodes = node, etc.). IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: b_ ! uvec of unknown point LOGICAL, INTENT(IN) :: cold_start ! mode switch INTEGER, INTENT(IN) :: num_ele ! count of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: node ! element definitions REAL, DIMENSION(:,:),INTENT(IN) :: xyz_node ! uvecs of nodes REAL, DIMENSION(:,:),INTENT(IN) :: center ! center uvecs of elements REAL, DIMENSION(:), INTENT(IN) :: a_ ! (plane) areas of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: neighbor ! neighbors of each element LOGICAL, INTENT(OUT) :: success ! OUTPUT INTEGER, INTENT(INOUT) :: iele ! OUTPUT REAL, INTENT(INOUT) :: s1, s2, s3 ! OUTPUT !- - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER, PARAMETER :: infinity = 1000000 ! used to test for complex infinite loops INTEGER :: back1, back2, back3, i, iet, l_, passes REAL :: lat, lon, r2, r2min, s1t, s2t, s3t REAL, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit success = .FALSE. IF (cold_start) THEN iele = 0 ! default s1 = 0.0; s2 = 0.0; s3 = 0.0 ! default !find closest element center to initialize search r2min = 4.01 ! radians DO l_ = 1, num_ele r2 = (b_(1) - center(1,l_))**2 +(b_(2) - center(2,l_))**2 +(b_(3) - center(3,l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 1 radian away, give up. tv = center(1:3, iet) IF (DOT_PRODUCT(b_, tv) < 0.540) RETURN END IF ! cold_start ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 back3 = -3 passes = 0 is_it_here: DO passes = passes + 1 ! do not permit infinite loop, no matter how complex: IF (passes >= infinity) THEN CALL Uvec_2_LonLat(b_, lon, lat) WRITE (*, *) WRITE (*, "(' ERROR: Infinite loop when seeking element containing (',F9.4,', ',F8.4,'):' & & /' iet = ',I6,', back1 = ',I6,', back2 = ',I6,', back3 = ',I6)") & & lon, lat, iet, back1, back2, back3 WRITE (*, *) RETURN END IF ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ! then, check for infinite loop between 3 elements! ELSE IF (iet == back3) THEN ! in loop; force location in one or the other! CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= -0.00001) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= -0.00001) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= -0.00001) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iele = iet s1 = s1t s2 = s2t s3 = s3t success = .TRUE. END SUBROUTINE Which_Spherical_Triangle END PROGRAM Long_Term_Seismicity_v3