PROGRAM Long_Term_Seismicity_v11 !========================================================= !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, ! Kagan_2009_GJI_I_scores, and pseudoCSEP]: ! regular (longitude, latitude) grid 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 typically 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 Danijel 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, Planetary, 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@epss.ucla.edu ! URL: http://peterbird.name ! 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 ! *Michele M. C. Carafa, INGV Roma & L'Aquila, Italia. ! 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 of the ! government of the state of California, or of the agencies listed. USE DSphere ! file DSphere.f90 = module of DOUBLE PRECISION Spherical-geometry tools by P. Bird, UCLA !========= IMSL option ==================================================================================================== !USE Numerical_Libraries ! IMSL Library, used for ANORDF = distribution function of a standard normal (Gaussian) variable !=============== OR ======================================================================================================= !========= Intel MKL-VSL option =========================================================================================== USE mkl_vsl ! for VSCdfNorm, VDCdfNorm = Single/Double-precision distribution function of a standard normal (Gaussian) variable. !========================================================================================================================== !Here are USAGE NOTES showing how to code each of these alternative functions: ! !--- IMSL version: ------------------------------------------------------------------------------- ! ! !ASS_quantile = ANORDF(sigmas) ! ! using standard Gaussian distribution function from IMSL ! ! !--- MKL version: -------------------------------------------------------------------------------- ! ! INTEGER(4), PARAMETER :: I4_one = 1 ! needed as input to some MKL subprograms, if default INTEGER ever changes to INTEGER(8). ! REAL, DIMENSION(1) :: VSCdfNorm_in, VSCdfNorm_out, VSCdfNormInv_in, VSCdfNormInv_out ! (and/or) ! DOUBLE PRECISION, DIMENSION(1) :: VDCdfNorm_in, VDCdfNorm_out, VDCdfNormInv_in, VDCdfNormInv_out ! ... ! VSCdfNorm_in(1) = sigmas ! CALL VSCdfNorm(I4_one, VSCdfNorm_in, VSCdfNorm_out) ! ASS_quantile = VSCdfNorm_out(1) ! (and/or) ! VDCdfNorm_in(1) = sigmas ! CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) ! ASS_quantile = VDCdfNorm_out(1) ! ! !------------------------------------------------------------------------------------------------- IMPLICIT NONE CHARACTER*1 :: c1, class_continuity_c1, first_byte, star_c1, step_continuity_c1, tab CHARACTER*2 :: c2, highlight_c2, issue_date_day, issue_date_month, start_date_day, start_date_month CHARACTER*3 :: class CHARACTER*4 :: c4, issue_date_year, start_date_year CHARACTER*5 :: c5, c5a CHARACTER*6 :: c6 CHARACTER*6, DIMENSION(12) :: T5_column_header, T5_catalog_name CHARACTER*6, DIMENSION(12) :: IT_column_header, IT_catalog_name CHARACTER*6, DIMENSION(:), ALLOCATABLE :: segment_c6 CHARACTER*80 :: fraction_GRD_filename, global_feg_pathfile, global_v_pathfile, question, output_grd_pathfile, & & PB2002_steps_pathfile, template_pathfile CHARACTER*80, DIMENSION(:), ALLOCATABLE :: e_nko_pathfile, h_nko_pathfile, orogen_f_dig_pathfile, orogen_feg_pathfile CHARACTER*132 :: line CHARACTER*3, DIMENSION(0:7) :: class_c3 ! = INT, plus boundary classes 1-7 DOUBLE PRECISION :: allowed_SS_impurity, angle_sum, argument, average, & & B, & & catalog_threshold_moment, continuum_fraction, continuum_percent, continuum_sum, corner_moment, & & desired_threshold_moment_Nm, dI, difference, dot_a, dot_b, dx, dy, & & expectation, & & fault_fraction, fault_sum, fault_percent, & & G_function, & & moment_rate_of_step_SI, & & point_value, & & seismicity_at_catalog_threshold, seismicity_increment, sum_total, & & 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 :: seismicity_bank ! parallel structure to that of seismicity DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: PB2002_seismicity ! parallel structure to that of seismicity DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: sum_of_fractions ! IF (subdivide), should equal array of 1.00 values. 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, & & beyond_loc, & & column, columns, crossings_to_W, & & dip_degrees, & & element, elevation, end_loc, & & format_selector, & & highest_fault_number, highlight, & & i, iele, ii, internal_ios, ios, IT_column, & & j, j_East, j_West, jj, jp, & & k, k_class_index, & & l, l_, larger_depth_integer, last_highlight, line_number, lp, & & m, m_time, ma, max_in_outline, MB, & ! Note that "MB" is also "mb" (with different, also temporary meaning) in a different block of code. & n, na, nb, nfl, nr1, nr2, n_step, nseg, number, number_of_zeros, numnod, numel, & & orogen, orogens, other_segment, & & region_index, row, rows, & & segment, segments, sequence_number, smaller_depth_integer, start_loc, & & T5_column, this_fault_number, threshold_index, timesteps_needed, top_threshold_index, trace_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(:,:), ALLOCATABLE :: age_1p5 ! MAY or MAY NOT be used, depending on region_index INTEGER(KIND = 2), DIMENSION(:,:), ALLOCATABLE :: continental_int2 INTEGER, DIMENSION(:), ALLOCATABLE :: n_in_outline ! (orogen = 1, ..., orogens) INTEGER, DIMENSION(:), ALLOCATABLE :: fault_dip_degrees, segment_dip_degrees INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor, nodes LOGICAL :: cold_start, creeping, decision_in, faults_used, got_dip_degrees, & & include_in_integral, inside, mated, maybe, oceanic, overlap, & & parallel, problem, proximal, repeat_sweep, skip_faults, strike_slip, subdivide, success LOGICAL*1, DIMENSION(:), ALLOCATABLE :: faults_used_in_orogen, segment_creeping, segment_used LOGICAL*1, DIMENSION(:,:), ALLOCATABLE :: continental, template_mask REAL*8, PARAMETER :: Earth_radius_m = 6371000.0D0 !non-SUB boundaries (symmetrical): REAL*8, PARAMETER :: CCB_halfwidth_km = 189.0D0 REAL*8, PARAMETER :: CTF_halfwidth_km = 257.0D0 REAL*8, PARAMETER :: CRB_halfwidth_km = 154.0D0 REAL*8, PARAMETER :: OSR_halfwidth_km = 132.0D0 REAL*8, PARAMETER :: OTF_halfwidth_km = 128.0D0 REAL*8, PARAMETER :: OCB_halfwidth_km = 186.0D0 ! these values from Table 1 of Bird & Kagan [2004] !SUB boundaries (asymmetrical): based on EQ_classification_II.f90 of Bird & Kagan [2004; BSSA] REAL*8, PARAMETER :: SUB_landward_halfwidth_km = 220.0D0 ! landward width of subduction lanes, in km from trench REAL*8, PARAMETER :: SUB_seaward_halfwidth_km = 135.0D0 ! seaward half-width of subduction lanes, in km from trench REAL*8, PARAMETER :: SUB_peakSeismicity_km = 55.0D0 ! location of seismicity peak, relative to trench !and one new SUB parameter added for v2 of LTS, following Bird et al. [2009, BSSA]: REAL*8, PARAMETER :: SUB_bending_fraction = 0.38D0 ! relative seismicity level when main interplate thrust is greased REAL*8, PARAMETER :: alongStep_sigma_km = 32.5D0 ! so combined EQ mislocation and step-mislocation unlikely to exceed twice this value REAL*8 :: angle, area_m2, azimuth_radians, & & basis_length_km, best_separation, & & center_latitude, center_longitude, & & 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, dip_radians, 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, & & fraction, fraction_inside, fraction_toward_uvec1, fraction_toward_uvec2, & & from_trench_azimuth_radians, & & half_R2, 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, other_azimuth_radians, & & parallelogram_area_m2, peak_diffusivity, & & quantile1, quantile2, & & radius_km, & & s_km, s_per_year, s_radians, s1, s2, s3, safe_timestep, scalar_product, scales, separation, SIGN, 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, triangle_area_m2, & & u1phi, u1theta, u2phi, u2theta, & & v_oblique_mmpa, v_oblique_mps, velocity_mmpa, vs_mmpa, & & X, & & Y, Y_factor_of_step REAL*8, 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, other_begin_uvec, other_end_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, & & v1Mv3, v2Mv1, v3Mv2 REAL*8, DIMENSION(4) :: dot_side ! sine of angle between test point and sides of a subduction lane; positive inward REAL*8, DIMENSION(6) :: h_km ! corresponds to h_k in Bird & Kagan [2004]; no entry for SUB because that is asymmetrical REAL*8, DIMENSION(12) :: T5_catalog_years, T5_catalog_s, & ! Table 5 of Bird & Kagan [2004, BSSA] & T5_catalog_area_m2, & & T5_catalog_pure_events, T5_catalog_pure_event_rate, & & T5_catalog_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*8, DIMENSION(12) :: IT_catalog_years, IT_catalog_s, & ! SUBSTITUTE "Italia" coefficients, added in LTS_v5+. & IT_catalog_area_m2, & & IT_catalog_assigned_events, IT_catalog_assigned_event_rate, & & IT_catalog_threshold_moment, & & IT_beta, & & IT_corner_magnitude, IT_corner_moment, & & IT_tGR_model_moment_rate, & & IT_assumed_dip_degrees, IT_assumed_dip_radians, & & IT_assumed_mu_GPa, IT_assumed_mu_Pa, & & IT_coupledThickness_km, IT_coupledThickness_m, & & IT_assumed_lithosphere_km, IT_assumed_lithosphere_m, & & IT_coupling REAL*8, DIMENSION(12) :: T5_supplement_seismicity_fraction_based_on_z ! == 1.0 if format_selector == 1; otherwise less (user input). REAL*8, DIMENSION(12) :: IT_supplement_seismicity_fraction_based_on_z ! == 1.0 if format_selector == 1; otherwise less (user input). REAL*8, DIMENSION(3,5) :: lane_uvecs REAL*8, DIMENSION(3,7) :: Gauss_point REAL*8, DIMENSION(3,2,2) :: G ! nodal functions of spherical triangle element REAL*8, DIMENSION(3,2,2,2):: dG ! derivitives of nodal functions G REAL*8, DIMENSION(2, 130) :: GCNocb ! digitized ocean-continent boundary offshore Gorda-Cal-Nev orogen REAL*8, DIMENSION(:), ALLOCATABLE :: a_ REAL*8, DIMENSION(:), ALLOCATABLE :: fault_azimuth_radians REAL*8, DIMENSION(:), ALLOCATABLE :: RELM_threshold_magnitude, RELM_threshold_moment_Nm REAL*8, DIMENSION(:), ALLOCATABLE :: segment_slip_rate_mmpa, segment_heave_rate_mmpa REAL*8, DIMENSION(:), ALLOCATABLE :: vw REAL*8, DIMENSION(:,:), ALLOCATABLE :: center, element_principal_strainrates, grid_cell_areas, node_uvec REAL*8, DIMENSION(:,:), ALLOCATABLE :: diffusivity ! (rows {N->S}, columns {W->E}) REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: outline_uvecs ! (1 ~ 3, 1 ~ n_in_outline(orogen), 1 ~ orogens) REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: outline_segments ! used in creation of .feg outlines (only) REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: segment_uvecs, segment_outer_uvecs REAL*8, 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. !-------------------------------------------------------------------------------------------------- !Dummy vectors (of 1 entry each) are needed to sync with MKL subprogramsL INTEGER(4), PARAMETER :: I4_one = 1 ! needed as input to some MKL subprograms, if default INTEGER ever changes to INTEGER(8). INTEGER(4) :: status ! specifying (4) in case the default INTEGER type ever changes due to compile-time switches(?) REAL, DIMENSION(1) :: VSCdfNorm_in, VSCdfNorm_out, VSCdfNormInv_in, VSCdfNormInv_out DOUBLE PRECISION, DIMENSION(1) :: VDCdfNorm_in, VDCdfNorm_out, VDCdfNormInv_in, VDCdfNormInv_out !-------------------------------------------------------------------------------------------------- !Angles assumed in conversion of (distributed permanent strain rates) to (seismic moment rate), !per Appendix of M.M.C. Carafa & P. Bird [2016?, J. Geophys. Res.-Solid Earth]: REAL*8 :: CF_dip_degrees, EF_dip_degrees, SS_alpha_degrees ! N.B. SS_alpha_degrees is NOT the dip; it is the angle between fault and horizontal epsilon^dot_great. REAL*8 :: alpha_1_degrees, alpha_2_degrees !These are the angles (in degrees) between the dominant (highest-rate) principal strain-rate !axis and the active fault planes in the more-active and less-active conjugate sets, respectively. !Note a choice of 45 degrees gives minimum forecast seismicity, and !that changing either angle away from 45 degrees will result in ! HIGHER seismicity rates due to distributed permanent deformation. !However, this may desirable; see the version-notes in this source code file for more discussion. !Also note that the formulas employing these two angles are symmetrical about 45 degrees, !so setting alpha_1_degrees to 25 degrees has the same effect as 65 degrees, and so on. REAL*8 :: alpha_1_radians, alpha_2_radians REAL*8 :: k_alpha_1, k_alpha_2 ! Each is equal to 2.0D0 when both alpha angles are 45 degrees; ! these coefficients will be larger for any other angle(s). REAL*8 :: epsilonDot_great, epsilonDot_medium, epsilonDot_small ! ordered absolute values of the principal strain-rates !-------------------------------------------------------------------------------------------------- DATA Gauss_point / & & 0.3333333333333333D0, 0.3333333333333333D0, 0.3333333333333333D0, & & 0.0597158733333333D0, 0.4701420633333333D0, 0.4701420633333333D0, & & 0.4701420633333333D0, 0.0597158733333333D0, 0.4701420633333333D0, & & 0.4701420633333333D0, 0.4701420633333333D0, 0.0597158733333333D0, & & 0.7974269866666667D0, 0.1012865066666667D0, 0.1012865066666667D0, & & 0.1012865066666667D0, 0.7974269866666667D0, 0.1012865066666667D0, & & 0.1012865066666667D0, 0.1012865066666667D0, 0.7974269866666667D0/ tab = CHAR(9) ! character "HT" for horizontal tab half_R2 = 0.5D0 * (Earth_radius_m**2) radius_km = Earth_radius_m / 1000.0D0 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.5D0 * (SUB_landward_halfwidth_km - SUB_peakSeismicity_km) ! sigmas for Gaussian approximation of B SUB_seaward_sigma_km = 0.5D0 * (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_catalog_name = (/"CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT "/) T5_catalog_years = (/ 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474/) T5_catalog_area_m2 = (/ 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14/) T5_catalog_threshold_moment = (/ 1.13D17, 3.5D17, 3.5D17, 3.5D17, 1.13D17, 1.13D17, 2.0D17, 2.0D17, 2.0D17, 3.5D17, 3.5D17, 3.5D17/) ! ---------------------------------------------------------------------------------------------------------------------------------------------- T5_catalog_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_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.67D12, 3.8D12, 3.20D12, 7.40D12, 6.7D11, 1.9D11, 6.7D12, 9.4D11, 9.0D11, 4.6D12, 5.55D13, 2.29D14/) 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.5D8, 4.4D8, 2.9D8, 3.0D8, 5.0D9, 4.7D8, 5.2D8, 5.3D8, 5.5D8, 1.2D9, 5.24D9, 1.06D10/) 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/) ! **************************************************************************************** ! Seismicity coefficients, ***WITH EDITS***, originally based on 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].} ! IT_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 IT_column_header = (/"CRB ","CTF ","CCBslo","CCBfas","OSRnor","OSRoth","OTFslo","OTFmed","OTFfas","OCB ","SUBslo","SUBfas"/) ! ---------------------------------------------------------------------------------------------------------------------------------------------- IT_catalog_name = (/"CPTI ","CPTI ","CPTI ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT ","CMT "/) IT_catalog_years = (/ 135.00, 135.00, 135.00, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474, 25.7474/) IT_catalog_area_m2 = (/ 1.41D12, 1.41D12, 1.41D12, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14, 4.72D14/) IT_catalog_threshold_moment = (/ 1.78D16, 1.78D16, 1.78D16, 3.5D17, 1.13D17, 1.13D17, 2.0D17, 2.0D17, 2.0D17, 3.5D17, 3.5D17, 3.5D17/) ! ---------------------------------------------------------------------------------------------------------------------------------------------- IT_catalog_assigned_events = (/ 154.3, 44.1, 147.7, 181.1, 424.3, 77.0, 398.0, 406.9, 376.6, 117.7, 399.5, 1653.3/) IT_beta = (/ 0.634, 0.714, 0.714, 0.62, 0.92, 0.82, 0.64, 0.65, 0.73, 0.53, 0.64, 0.64/) IT_corner_magnitude = (/ 7.170, 6.620, 6.620, 8.46, 5.86, 7.39, 8.14, 6.55, 6.63, 8.04, 9.58, 9.58/) IT_tGR_model_moment_rate = (/ 32.29D9, 3.83D9, 12.83D9, 7.40D12, 6.7D11, 1.9D11, 6.7D12, 9.4D11, 9.0D11, 4.6D12, 5.55D13, 2.29D14/) IT_assumed_dip_degrees = (/ 49., 76., 35., 20., 55., 55., 73., 73., 73., 20., 14., 14./) IT_assumed_mu_GPa = (/ 35.2, 35.2, 35.2, 27.7, 25.7, 25.7, 25.7, 25.7, 25.7, 49., 49., 49./) IT_coupledThickness_km = (/ 7.0, 5.0, 3.8, 24.9, 0.13, 0.40, 13., 1.8, 1.6, 3.8, 10.6, 21.7/) IT_assumed_lithosphere_km = (/ 12., 11., 11., 24.9, 8., 8., 14., 14., 14., 14., 26., 26./) IT_coupling = (/ 0.58, 0.45, 0.35, 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.25D0 * 24 * 60 * 60 DO i = 1, 12 T5_catalog_s(i) = T5_catalog_years(i) * s_per_year T5_catalog_pure_event_rate(i) = T5_catalog_pure_events(i) / T5_catalog_s(i) T5_corner_moment(i) = DMoment(T5_corner_magnitude(i)) T5_length_m(i) = T5_length_km(i) * 1000.0D0 T5_mean_velocity_mps(i) = T5_mean_velocity_mmpa(i) * 0.001D0 / 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.0D9 T5_coupledThickness_m(i) = T5_coupledThickness_km(i) * 1000.0D0 T5_assumed_lithosphere_m(i) = T5_assumed_lithosphere_km(i) * 1000.0D0 !-------------------------------------------------------------------------- IT_catalog_s(i) = IT_catalog_years(i) * s_per_year IT_catalog_assigned_event_rate(i) = IT_catalog_assigned_events(i) / IT_catalog_s(i) IT_corner_moment(i) = DMoment(IT_corner_magnitude(i)) IT_assumed_dip_radians(i) = IT_assumed_dip_degrees(i) * radians_per_degree IT_assumed_mu_Pa(i) = IT_assumed_mu_GPa(i) * 1.0D9 IT_coupledThickness_m(i) = IT_coupledThickness_km(i) * 1000.0D0 IT_assumed_lithosphere_m(i) = IT_assumed_lithosphere_km(i) * 1000.0D0 END DO ! **************************************************************************************** WRITE (*, "(/' PROGRAM Long_Term_Seismicity_v11')") WRITE (*, "(/' Predicts long-term-average shallow seismicity of Earth,')") WRITE (*, "( ' without any time-dependence. This is a stationary forecast,')") WRITE (*, "( ' and it includes aftershocks as well as mainshocks, without distinction.')") !VersionHere WRITE (*, "(/' by Peter Bird, UCLA; version 11 of 8 January 2018.')") !----------------------------------------------------------------------------------- !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 NOT work around the poles, no smoothing is ! performed at latitudes higher than 75 degrees away from the equator. !----------------------------------------------------------------------------------- !Version 3.1 adds reporting of 2 additional numbers: ! the fractions of seismicity coming from modeled faults, and the ! fraction coming from distributed permanent deformation. ! This improvement programmed 2012.01.31 for UCERF3. !----------------------------------------------------------------------------------- !Version 3.1_Mediterranean has 2 special changes: ! (1) age_1p5.grd is not used to decide continental/oceanic crust; instead, ! an elevation limit of -1500 m (from ETOPO5.grd) is used. ! (2) The smoothing distance around the edge of any NeoKinema model is ! reduced from 100 km to 30 km. !----------------------------------------------------------------------------------- !Version 4 has these new features: ! *The version 3.1_Mediterranean settings are now available by choosing region_index = 2 ! in the single, unified program Long_Term_Seismicity_v4. ! *All former REAL variables and calculations are now REAL*8, using DSphere and ! REAL*8 versions of the Shells/NeoKinema routines Gjxy, Del_Gjxy, & Erate. ! (This causes an important reduction in forecast seismicity of intraplate regions.) ! *Fault dips, where known, are now read from the f[token].dig file(s). ! *To speed up the diffusion step when the forecast grid does not include any ! orogen boundary/boundaries, a larger (but still safe) timestep is used. !----------------------------------------------------------------------------------- !Version 5a has these new features: ! *When "Mediterranean settings" are chosen (region_index == 2), the seismicity ! coefficients will NOT come from the T5_... variables that represent ! Table 5 of Bird & Kagan [2004, BSSA], plus the later findings of ! Bird et al. [2009, BSSA, ibid]. INSTEAD, they will come from the ! alternate variables ME_... which MAY have MANY identical values, but will ! also have SOME CHANGED VALUES (e.g., in tectonic zones CRB & slow-CCB?). ! This feature was added to support joint research with ! Dr. Michele M. C. Carafa of INGV, L'Aquila, Italy, in 2015, in order to ! permit more accurate simulations of seismicity in the Italian region. ! *Optional switch "subdivide", if set .TRUE. by user, will provide 12 additional ! output .GRD files, full of real numbers in range of [0., 1.] that express ! the fraction of total forecast seismicity that is attributable to coupling ! in one column of the T5_... (or IT_...) table. For example, extra ouput ! file column_01_fraction.grd shows the fraction of total forecast ! seismicity which is due to the coupled thickness (and other parameters) ! in the CRB (Continental Rift Boundary) column of the table. ! Choosing this option will increase execution time by x13 ! !Version 5b has these new features: ! *Conversion of distributed permanent strain-rates to seismic moment rates ! is according to the formula in the Appendix to M.C.C. Carafa & ! P. Bird [2016, J. Geophys. Res.-Solid Earth]. Note that the ! angles of theta_1_degrees = 45.0D0 and theta_2_degrees = 45.0D0 ! assigned when (region_index /= 2) will result in computation of ! MINIMUM seismic moment rates (and forecast seismicity). ! This will give results compatible with the formula used by ! Bird & Liu [2007, Seism. Res. Lett.] (their revised formula), ! Bird et al. [2010, SRL], and Bird & Kreemer [2015, Bull. Seismol. Soc. ! Am.]; HOWEVER, also note that in the latter 2 papers we also needed an ! empirical calibration factor, which ranged from 0.86 to 3.4; the need ! for this correction factor could be attributed to incorrect angles. ! Therefore, the user may prefer to set other angle values, and recompile? ! Currently, when (region_index == 2) the angles are taken from ! Table 4 in Carafa & Bird [2016, J.G.R.]. ! *This program now asks the user whether each NeoKinema model that is to be ! used has any active faults? If not, the user will not be asked to ! provide input file f_{token}.dig and h_{token}.nko, because they do ! not exist. !Version 5c has these changes: ! *Seismogenic thickness z, coupling coefficient c, and coupled thickness ! are changed to "Mediterranean" (Italian) values in the CRB and slow-CCB ! columns of the IT_ tables (relative to the old standard values in the ! same columns of the T5_ table). These values come from modeling of ! Italian neotectonics reported in Carafa & Bird [2016, J.G.R.]. ! Because of this change, LTSv5 will report different (higher) seismicity ! when region_index == 2 than when it has some other value. !----------------------------------------------------------------------------------- !Version 6a has these changes: ! *In a purely cosmetic change, the region_index == 2 options that were formerly ! known as "Mediterranean" (with arrays labelled "ME_...") are now known ! by the labels "Italia" and "IT_...". This better reflects the true ! origin and relevance of the parameters and choices made in joint research ! between Peter Bird (UCLA) and Michele Carafa (INGV). ! *In the computation of seismic moment rates from NeoKinema faults, the down-dip ! fault width is now computed with the user-supplied fault dip (if any) ! contained in the f.dig input file. Users who are unsure of fault dips ! should leave them unspecified; this program will use reasonable defaults ! (as NeoKinema does). If a fault dip is specified, and especially when it ! is a small angle, the user should be careful that it reflects the mean dip ! over the whole seismogenic depth range. Zero dip is not allowed. ! *For each NeoKinema F-E grid (.feg) file read, the user was formerly prompted ! to supply an outline.dig file (which can be obtained by running program ! NeoKineMap and asking it to plot the outline). However, this left room ! for errors in the outline format (e.g, a clockwise outline provided, when a ! counterclockwise outline was wanted), and for small inconsistencies arising ! when the .feg file was edited but the outline was not re-created. Therefore, ! the current Long_Term_Seismicity_v6 computes this outline for itself, and ! thus uses one less input file (per orogen modeled with NeoKinema). !----------------------------------------------------------------------------------- !Version 7 has this change: ! *New procedure/parameter option added to the REGIONAL SETTINGS list (below): ! 3: SyntheticExample ! This option is for processing simple, idealized NeoKinema models like that ! shown in Figures 2~6 of Bird & Carafa [2016, JGR]. ! No global datasets are read. The forecast region is assumed to be all- ! continental. No background Shells model is read, and no smoothing is applied ! around the edges of the NeoKinema model. Only a single NeoKinema model is input, ! and the forecast grid must be entirely within the domain of this model. !----------------------------------------------------------------------------------- !Version 8 has these changes: ! *Under the "Italia" options (region_index == 2) we now specify the dips of ! normal faults and thrust faults separately, instead of selecting one generic ! "alpha" angle between the principal strain-rate axis of greatest absolute ! magnitude and the two sets of conjugate faults. ! *Parameter "allowed_SS_impurity" is introduced, and set to 0.250 for the ! Italia region. (Elsewhere, and previously, it was 0.364.) This parameter ! governs how different the strain-rate may be from pure strike-slip before ! the strain-rate of an element is classified as having a compressional-boundary ! analog or an extensional-boundary analog instead of a transform-fault analog. ! *Values (especially coupled thickness of CTF, CRB, and CCB) in Table IT_ ! were updated (perhaps several times!). Check program date above. !----------------------------------------------------------------------------------- !Version 9 has no changes in algorithm or output. However, some parts of the ! source code were re-written to improve clarity and make program upkeep ! (and future extension) easier: ! *Because LTS uses both CMT catalog data (in table T5_) and CPTI catalog data ! (in certain columns of table IT_, used only when region_index == 2), ! it was not correct to assume that all seismicity calibrations in the code ! refer to CMT-derived quantities. Thus, the text-string "CMT_" has been ! replaced throughout by "catalog_". The identity of the catalog is now ! displayed in both table T5_ and table IT_ in the new text variables ! T5_catalog_name and IT_catalog_name in the 2nd row of each table. ! *Because the useful parts of different catalogs span different numbers of ! years, it was not correct to build the CMT catalog duration into program ! logic. Instead, the length of each catalog now appears (in both table T5_ ! and table IT_) in the new rows (#3) T5_catalog_years and IT_catalog_years. ! *To explicitly recognize that different catalogs cover different geographic ! areas, new rows (#4) T5_catalog_area_m2 and IT_catalog_area_m2 have been ! added to the tables. These are not currently used in calculations, but ! they help to further clarify the origins of the following data rows. ! *To complete the set of 4 adjacent rows (#2~#5) fully identifying each catalog ! (by name, duration, area, and threshold) the existing rows T5_catalog_threshold_moment ! and IT_catalog_threshold_moment have been moved UP, to come BEFORE the rows ! T5_catalog_pure_events and IT_catalog_assigned_events, respectively. ! *Because the "pure" (outside-of-orogen) distinction only applies to the ! study of the global CMT catalog by Bird & Kagan [2004], using the PB2002 ! plate model, the corresponding row name in the IT_ table has been changed ! to IT_assigned_events (meaning that they were assigned to this column, ! on the basis of their tectonic setting and focal-mechanism). ! *Certain rows of table T5_ cannot be supplied for a regional table like IT_ ! because they refer to rigid-plate models (outside of orogen regions) that ! were specific to plate model PB2002 and to the Bird & Kagan [2004] paper. ! To reduce confusion, the following table rows have been deleted: ! IT_length_km, IT_mean_velocity_mmpa, IT_lineIntegral_Nps. ! *To sum-up the usage of customized regional data presented in ! modified columns of table IT_ (and any future regional customizations): ! -->The computation of a local tectonic moment rate depends on the ! product of rows IT_coupledThickness_km * IT_assumed_mu_GPa ! (after both are converted to SI units, of course). ! -->The rows: IT_catalog_assigned_events and IT_tGR_model_moment_rate are ! used to convert this local model moment rate into an earthquake rate, above ! the threshold of the governing catalog. THEREFORE, it should not matter ! how many years, or how many square-meters, are covered by the governing ! catalog, as long as IT_tGR_model_moment_rate is computed for the same ! geographic area in which IT_catalog_assigned_events were observed. ! -->The 3 rows: IT_catalog_threshold_moment, IT_beta, and IT_corner_magnitude ! are all used in computing the tapered Gutenberg-Richter G-function that ! scales the earthquake rate from the catalog threshold magnitude to the ! desired magnitude theshold of the new seismicity model. ! -->The 2 rows: IT_assumed_dip_degrees and IT assumed_lithosphere_km are ! used together, to decide how far from the mapped trace the seismicity of a ! specific fault should extend (once the calculations above have been finished). ! (However, the measured dip of that specific fault will be substituted ! for the generic dip IT_assumed_dip_degrees, when it is known.) ! -->Row IT_coupling is currently not used in any computation, but serves to ! help clarify the difference between IT_coupledThickness_km and ! IT_assumed_lithosphere_km. When changing the data in IT_, be sure to ! enter values so that IT_coupling = IT_coupledThickness_km / ! IT_assumed_lithosphere_km. ! {PB, UCLA, 2016.07.11} !----------------------------------------------------------------------------------- !Version 10 adds a correction for missing (or duplicated) areas in the "gores" ! that occur down-dip of some significant bends in the traces of non-vertical ! NeoKinema faults, for more accurate depiction of the active fault areas and ! their implied seismicity maps. This change is consistent with the update ! from NeoKinema_v4.3 => NeoKinema_v5.0, although the actual changes to this code ! are not the same. ! -N.B. There is no need to add/subtract triangular areas from Subduction ! zones (either in PB2002 plate-boundary steps, or NeoKinema), because ! this code Long_Term_Seismicity has always used parallel ! "subduction lanes", which is another way of avoiding ! the problem of missing (or extra) triangular areas. ! However, Long_Term_Seismicity_v10 includes another improvement, ! so that Subduction faults modeled by NeoKinema are recognized as ! parts of the "same Subduction zone" (for defining seismicity lanes) ! EVEN IF the S-trace has been divided into different NeoKinema faults ! (F0001S, F0002S, F0003S) in order to reduce unwanted along-strike ! averaging of strike-slip components of slip rate (e.g., between an ! oblique portion and two adjacent orthogonally-subducting portions). ! {PB, INGV L'Aquila, 2016.11.07} !----------------------------------------------------------------------------------- !Version 11 fixes a bug caused by the new "gore" feature of Version 10: ! Where a negative (SIGN == -1.0D0) gore extends slightly beyond the two ! overlapped rectangular regions it was meant to correct, seismicity would be ! over-corrected and could become negative. We now prevent any correction ! so large that it would produce a negative seismicity. (Note that this is ! only a rough and partial first-order fix of a complex geometric problem.) ! {PB, UCLA, 2017.10.17} !----------------------------------------------------------------------------------- 1 WRITE (*, "(/' CHOICE OF REGIONAL SETTINGS:')") WRITE (*, "( ' ')") WRITE (*, "( ' 0: global, or default, settings')") WRITE (*, "( ' 1: western United States settings')") WRITE (*, "( ' 2: Italia settings')") WRITE (*, "( ' 3: SyntheticExample')") WRITE (*, "( ' ===================================')") WRITE (*, "(/' Enter (integer) code to select desired regional settings: '\)") READ (*, *) region_index IF ((region_index < 0).OR.(region_index > 3)) THEN WRITE (*, "(' ERROR: Try again, selecting an integer from the list:')") CALL Pause() GO TO 1 END IF !Determine how different the strain-rate of an element may be (in relative, dimensionless terms) !from pure strike-slip, before the element is classified as having a compressional-boundary !analog or an extensional-boundary analog, rather than a transform-fault analog: IF (region_index == 2) THEN ! Italia region allowed_SS_impurity = 0.250D0 ELSE ! any other region index allowed_SS_impurity = 0.364D0 END IF !Establish angles that will be used in conversion of distributed permanent strain-rate !to seismic moment rate (per Appendix of Carafa & Bird [2016, J.G.R.]): IF (region_index == 2) THEN ! Italia region; values from Carafa & Bird [2016?, JGR?]: CF_dip_degrees = 35.0D0 ! Based on mean dip in DISS fault database. EF_dip_degrees = 50.0D0 ! Based on mean dip in DISS fault database. SS_alpha_degrees = 37.5D0 ! SS_alpha_degrees = (35 + (90 - 50)) / 2 = 37.5 degrees, ! taking the average of the CF-type and EF-type behaviors. ELSE ! any other region index; angles per Bird & Liu [2007], Bird et al. [2010], Bird & Kreemer [2015]: CF_dip_degrees = 45.0D0 EF_dip_degrees = 45.0D0 SS_alpha_degrees = 45.0D0 ! but note that this choice gives the MINIMUM forecast seismicity, and may underestimate in some cases. END IF IF (region_index == 2) THEN WRITE (*, *) WRITE (*, "(' CAUTION: Seismicity coefficients of tectonic zones will NOT be taken')") WRITE (*, "(' from the T5_... variables representing Table 5 of Bird & Kagan [2004],')") WRITE (*, "(' plus the later improvements of Bird et al. [2009, BSSA];')") WRITE (*, "(' instead, they will be taken from the IT_... tables, with parallel')") WRITE (*, "(' structure and MANY equal values, but SOME altered values.')") WRITE (*, "(' Please inspect the source code for details (and attached comments).')") CALL Pause() END IF IF (region_index < 2) THEN ! we need to store the age_1p5 grid: ALLOCATE ( age_1p5( -720:900, 0:3600) ) ! Note that row/column 0 means lat/lon of 0 END IF 11 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 Danijel Schorlemmer''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 11 END IF IF (format_selector == 1) THEN WRITE (*, "(' Enter threshold (minimum) magnitude: '\)") READ (*, *) threshold_magnitude desired_threshold_moment_Nm = DMoment(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 (region_index == 3) WRITE (*, "(' Note that the forecast grid must be entirely INSIDE the NeoKinema model domain!')") 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). IF (region_index == 2) THEN IT_supplement_seismicity_fraction_based_on_z = 1.00 ! whole array ELSE T5_supplement_seismicity_fraction_based_on_z = 1.00 ! whole array END IF ELSE IF (format_selector == 2) THEN ! request guidance on fractions of seismicity to report 12 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 IF (region_index == 2) THEN WRITE (*, "(' column ',I2,': ',A6,' seismicity fraction: '\)") i, IT_column_header(i) READ (*, *) IT_supplement_seismicity_fraction_based_on_z(i) IF ((IT_supplement_seismicity_fraction_based_on_z(i) < 0.0D0).OR. & & (IT_supplement_seismicity_fraction_based_on_z(i) > 1.0D0)) THEN WRITE (*, "(' ERROR! Each fraction must be in range 0.0 ~ 1.0. Start over!')") CALL Pause() GOTO 12 END IF ELSE 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.0D0).OR. & & (T5_supplement_seismicity_fraction_based_on_z(i) > 1.0D0)) THEN WRITE (*, "(' ERROR! Each fraction must be in range 0.0 ~ 1.0. Start over!')") CALL Pause() GOTO 12 END IF 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.0D0) 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 !precompute grid cell areas: ALLOCATE ( grid_cell_areas(rows, columns) ) DO i = 1, rows latitude = lat_max - (i - 1) * d_lat area_m2 = Area_of_cell(d_lon, d_lat, latitude, Earth_radius_m) grid_cell_areas(i, 1:columns) = area_m2 END DO 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) ! Note: OPENed early, to prevent last-minute crashes! 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_v11, 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_v11 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 IF (region_index < 3) THEN 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 ELSE ! region_index == 3: SyntheticExample orogens = 1 WRITE (*, "(' Exactly ONE NeoKinema model must be supplied.')") END IF IF (orogens > 0) THEN ALLOCATE ( orogen_f_dig_pathfile(orogens) ) ALLOCATE ( orogen_feg_pathfile(orogens) ) ALLOCATE ( e_nko_pathfile(orogens) ) ALLOCATE ( h_nko_pathfile(orogens) ) ALLOCATE ( faults_used_in_orogen(orogens) ) DO orogen = 1, orogens IF (orogen <= 9) THEN WRITE (c1, "(I1)") orogen question = "Were faults used in the NeoKinema model of orogen #" // c1 // '?' CALL DPrompt_for_Logical(TRIM(question), .TRUE., faults_used) faults_used_in_orogen(orogen) = faults_used IF (faults_used) THEN WRITE (*, "(/' Enter [Drive:][\path\]filename of f_.dig file for orogen #',I1,':')") orogen READ (*, *) orogen_f_dig_pathfile(orogen) END IF WRITE (*, "(/' Enter [Drive:][\path\]filename of _.feg file for orogen #',I1,':')") orogen READ (*, *) orogen_feg_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of e_.nko file for orogen #',I1,':')") orogen READ (*, *) e_nko_pathfile(orogen) IF (faults_used) THEN WRITE (*, "(/' Enter [Drive:][\path\]filename of h_.nko file for orogen #',I1,':')") orogen READ (*, *) h_nko_pathfile(orogen) END IF ELSE ! orogen >= 10 (two-digit): WRITE (c2, "(I2)") orogen question = "Were faults used in the NeoKinema model of orogen #" // c2 // '?' CALL DPrompt_for_Logical(TRIM(question), .TRUE., faults_used) faults_used_in_orogen(orogen) = faults_used IF (faults_used) THEN WRITE (*, "(/' Enter [Drive:][\path\]filename of f_.dig file for orogen #',I2,':')") orogen READ (*, *) orogen_f_dig_pathfile(orogen) END IF WRITE (*, "(/' Enter [Drive:][\path\]filename of _.feg file for orogen #',I2,':')") orogen READ (*, *) orogen_feg_pathfile(orogen) WRITE (*, "(/' Enter [Drive:][\path\]filename of e_.nko file for orogen #',I2,':')") orogen READ (*, *) e_nko_pathfile(orogen) IF (faults_used) THEN WRITE (*, "(/' Enter [Drive:][\path\]filename of h_.nko file for orogen #',I2,':')") orogen READ (*, *) h_nko_pathfile(orogen) END IF 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) IF (region_index < 3) THEN WRITE (*, "(' ETOPO5.grd (topography of Earth with 5 min. grid: 4320 columns, 2161 rows)')") IF (region_index /= 2) 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) END IF ! region_index < 3 (N.B. No global datasets are needed for region_index == 3: SyntheticExample.) IF (orogens > 0) THEN DO orogen = 1, orogens IF (faults_used_in_orogen(orogen)) WRITE (*, "(' ',A)") TRIM(orogen_f_dig_pathfile(orogen)) WRITE (*, "(' ',A)") TRIM(orogen_feg_pathfile(orogen)) WRITE (*, "(' ',A)") TRIM(e_nko_pathfile(orogen)) IF (faults_used_in_orogen(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!')") IF (format_selector == 1) THEN WRITE (*, *) WRITE (*, "(' Optionally, you may choose to subdivide the results of LTS_v5+ by producing')") WRITE (*, "(' 12 additional output .GRD files, each of which displays the fraction of')") WRITE (*, "(' total forecast seismicity that is due to each column of the table of')") WRITE (*, "(' seismicity coupling factors (T5_... or IT_... table), present at the')") WRITE (*, "(' top of the LTS_v5+ source code.')") WRITE (*, "(' CAUTION: Selecting this option will increase execution time x13 !')") CALL Prompt_for_Logical('Subdivide forecast seismicity into 13 files?', .FALSE., subdivide) ELSE subdivide = .FALSE. ! this option only available together with choice of .GRD format. END IF WRITE (*, *) 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)) IF (region_index < 3) THEN CALL Test_Opening_of_Existing_File("ETOPO5.grd") IF (region_index /= 2) CALL Test_Opening_of_Existing_File("age_1p5.grd") ! not used under 2:Italia settings 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)) END IF ! region_index < 3 (N.B. No global datasets are needed for region_index == 3: SyntheticExample.) IF (orogens > 0) THEN DO orogen = 1, orogens IF (faults_used_in_orogen(orogen)) CALL Test_Opening_of_Existing_File(TRIM(orogen_f_dig_pathfile(orogen))) CALL Test_Opening_of_Existing_File(TRIM(orogen_feg_pathfile(orogen))) CALL Test_Opening_of_Existing_File(TRIM(e_nko_pathfile(orogen))) IF (faults_used_in_orogen(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) ) IF (subdivide) THEN ALLOCATE ( seismicity_bank(rows, columns, 1) ) ! used to save results of total forecast, as a reference during subdivision/highlighting. ALLOCATE ( sum_of_fractions(rows, columns) ) ! should evaluate to a matrix of 1.00 values, in the end sum_of_fractions = 0.0D0 END IF 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 IF (region_index < 3) ALLOCATE ( PB2002_seismicity(rows, columns, 41) ) ! array for seismicity of global plate boundaries (only). 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.85D0 + i * 0.10D0 RELM_threshold_moment_Nm(i) = DMoment(RELM_threshold_magnitude(i)) END DO END IF !===================================================================================================================== !If any orogens were modeled with NeoKinema, then pre-allocate and pre-compute the outlines of each .feg file... IF (orogens > 0) THEN WRITE (*, *) WRITE (*, "(' Computing & storing outlines of all ',I2,' NeoKinema .feg files...')") orogens !Start by scanning for greatest number of nodes in any NeoKinema .feg file. !This will be used as a (generous) proxy for the maximum number of points in any outline. max_in_outline = 0 ! just initializing; this will increase below... DO orogen = 1, orogens !Read the .feg file and check its node-count... WRITE (*, "(' Looking at headers of ',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 IF (problem) THEN WRITE (*,"(' ERROR: Necessary information (numnod) absent or defective in this file.')") CLOSE (4) CALL Pause() STOP END IF CLOSE(4) max_in_outline = MAX(max_in_outline, (numnod+1)) END DO ! orogen = 1, orogens (initial inspection, just to scan for numnod values) !Allocate (generous) space to hold all the outlines of all the NeoKinema .feg grids: WRITE (*, "(' ALLOCATE ( outline_uvecs(3,',I10,',',I3,') )')") max_in_outline, orogens ALLOCATE ( n_in_outline(orogens) ) ! will record how many outline uvecs are ACTUALLY USED in each orogen ALLOCATE ( outline_uvecs(3, max_in_outline, orogens) ) !Re-read (all of) each NeoKinema .feg file, and then compute and store its outline... DO orogen = 1, orogens !Read the .feg file and store it (temporarily) in memory... WRITE (*, "(' Re-reading ',A,' in full ...')") 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) IF (ALLOCATED(node_uvec)) DEALLOCATE ( node_uvec ) 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 DLonLat_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) IF (ALLOCATED(nodes)) DEALLOCATE ( nodes ) 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) WRITE (*, "(' Computing & storing the outline of this NeoKinema .feg file...')") IF (ALLOCATED( outline_segments )) DEALLOCATE ( outline_segments ) ALLOCATE ( outline_segments(3, 2, numnod) ) ! Build library of external edge outline_segments (unsorted) nseg = 0 DO i = 1, numel ! first, look for external edges of spherical triangles DO j = 1, 3 jp = 1 + MOD(j,3) na = nodes(j,i) nb = nodes(jp,i) mated = .FALSE. might_mate_1: DO k = 1, numel ! mating to another spherical triangle? DO l = 1, 3 lp = 1 + MOD(l,3) ma = nodes(l, k) mb = nodes(lp, k) IF ((na == mb).AND.(nb == ma)) THEN mated = .TRUE. EXIT might_mate_1 END IF ! mate was found END DO ! l = 1, 3 sides of trial spherical triangle END DO might_mate_1 ! k = 1, numel IF (.NOT.mated) THEN nseg = MIN(nseg + 1, numnod) ! no problem expected here outline_segments(1:3, 1, nseg) = node_uvec(1:3, na) outline_segments(1:3, 2, nseg) = node_uvec(1:3, nb) !note that outline_segments always proceed counterclockwise around grid END IF ! NOT mated END DO ! j = 1, 3 END DO ! i = 1, numel; looking for external edges of spherical triangles !OPEN (UNIT = 55, FILE = "outline_of_"//TRIM(orogen_feg_pathfile(orogen))//".dig") !WRITE (55, "('outline of ',A)") TRIM(orogen_feg_pathfile(orogen)) !link outline_segments to create outline j = 1 ! begin with first segment uvec1(1:3) = outline_segments(1:3, 1, j) outline_uvecs(1:3, 1, orogen) = outline_segments(1:3, 1, j) n_in_outline(orogen) = 1 !CALL DUvec_2_LonLat(uvec1, longitude, latitude) !WRITE (55, "(1X,SP,ES12.5,',',ES12.5)") longitude, latitude DO i = 1, nseg uvec2(1:3) = outline_segments(1:3, 2, j) n_in_outline(orogen) = n_in_outline(orogen) + 1 outline_uvecs(1:3, n_in_outline(orogen), orogen) = outline_segments(1:3, 2, j) !CALL DUvec_2_LonLat(uvec2, longitude, latitude) !WRITE (55, "(1X,SP,ES12.5,',',ES12.5)") longitude, latitude find_next: DO k = 2, nseg IF (uvec2(1) == outline_segments(1, 1, k)) THEN IF (uvec2(2) == outline_segments(2, 1, k)) THEN IF (uvec2(3) == outline_segments(3, 1, k)) THEN j = k EXIT find_next END IF END IF END IF END DO find_next !prepare to loop uvec1 = uvec2 END DO ! i = 1, nseg !WRITE (55, "('*** end of line segment ***')") !CLOSE(55) END DO ! orogen = 1, orogens (computing and storing outlines of each NeoKinema .feg file) END IF ! (orogens > 0), computing and storing outlines of each .feg file !===================================================================================================================== !DECIDE whether this program will have a single pass (or multiple passes) through the following code... IF (subdivide) THEN last_highlight = 12 ! == highest-numbered column in T5_... or IT_... table. ELSE last_highlight = 0 ! and following code will only be executed ONCE. END IF DO highlight = 0, last_highlight ! potentially repeating the entire calculation 13 times over! !===================================================================================================================== seismicity = 0.0D0 ! whole area; all thresholds {if more than 1}; to be incremented below. IF (region_index < 3) THEN 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. END IF fault_sum = 0.0D0 continuum_sum = 0.0D0 ! zeroing the scalar accumulators; to be incremented only for lowest threshold in the run 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.D0 center_latitude = (low_template_latitude + high_template_latitude ) / 2.D0 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 IF (highlight == 0) THEN ! first pass through (and perhaps the only pass?) !--------------------------------------------------------------------------------- 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. IF (region_index == 3) THEN continental = .TRUE. ! because SyntheticExample includes GPS data, which is rarely available in the ocean basins. ELSE ! global datasets will be consulted: continental = .FALSE. ! whole array; to be changed to .TRUE. in continental areas (below) !============================================================================================================= ! Decide which grid points are "continental," using criteria of elevation ALONE. !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 DTraceback() 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 DTraceback() STOP END IF CLOSE (2) IF (region_index /= 2) THEN ! note that age_1p5 is not read or used under 2:Italia. !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 DTraceback() 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 DTraceback() 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 END IF ! region_index /= 2, so age_1p5.grd should be read and memorized 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)) IF (region_index /= 2) THEN !seafloor age at grid point: age_Ma = NINT(LookUp_age_1p5(longitude, latitude)) !decision: oceanic = (age_Ma < 180).OR.(elevation < -2000) ELSE ! region_index == 2: Italia oceanic = (elevation < -1500) END IF continental(i, j) = .NOT. oceanic END DO END DO END IF ! region_index == 3 (SyntheticExample), or NOT !!output the continental_int2 array in .GRD format for map-plotting !ALLOCATE ( continental_int2(rows, columns) ) !DO i = 1, rows ! DO j = 1, columns ! IF (continental(i, j)) THEN ! continental_int2(i, j) = 1 ! ELSE ! continental_int2(i, j) = 0 ! END IF ! END DO !END DO !OPEN (UNIT = 21, FILE = "continental.grd") !WRITE (21, "(3F11.5,' = lon_min, d_lon, lon_max')") lon_min, d_lon, lon_max !WRITE (21, "(3F11.5,' = lat_min, d_lat, lat_max')") lat_min, d_lat, lat_max !WRITE (21, "(40I2)") ((continental_int2(i, j), 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. !CLOSE (21) ! output .grd file !============================================================================================================= ! Pre-compute unit vectors ("uvecs") for all subsampling points at all grid points: MB = 8 * 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.0D0, MAX(-90.0D0, t_lat)) DO jj = 1, subsamplings t_lon = longitude + lon_offset + jj * dd_lon CALL DLonLat_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 ')") !--------------------------------------------------------------------------------- END IF ! highlight == 0; first pass through IF (region_index < 3) THEN !============================================================================================================= ! 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) IF (ALLOCATED( node_uvec )) DEALLOCATE ( node_uvec ) 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 DLonLat_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) IF (ALLOCATED( nodes )) DEALLOCATE ( nodes ) 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 DMake_Uvec (uvec4, uvec) ! center of element equat = DSQRT(uvec(1)**2 + uvec(2)**2) IF (equat == 0.0D0) THEN PRINT "(' Error: integration point ',I1,' of element ',I5,' is N or S pole.')", m, l_ CALL Pause() STOP END IF theta_ = DATAN2(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 DPrincipal_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 DLearn_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 DLonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) CALL DWhich_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.0D0) THEN ! rare: ideal strike-slip IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF ELSE IF (err > 0.0D0) THEN ! thrusting (possibly minor?) IF (err > (allowed_SS_impurity * e2h)) THEN ! primarily thrusting IF (region_index == 2) THEN IT_column = 3 ! slow CCB ELSE T5_column = 3 ! slow CCB END IF 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] IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (allowed_SS_impurity * e1h)) THEN ! primarily normal-faulting IF (region_index == 2) THEN IT_column = 1 ! CRB ELSE T5_column = 1 ! CRB END IF 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] IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF END IF END IF ! thrusting or normal-faulting? ELSE ! oceanic IF (err == 0.0D0) THEN ! rare: ideal strike-slip IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF ELSE IF (err > 0.0D0) THEN ! thrusting (possibly minor?) IF (err > (allowed_SS_impurity * e2h)) THEN ! primarily thrusting IF (region_index == 2) THEN IT_column = 10 ! OCB ELSE T5_column = 10 ! OCB END IF 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] IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (allowed_SS_impurity * e1h)) THEN ! primarily normal-faulting IF (region_index == 2) THEN IT_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 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.] END IF 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] IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF 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.0D0 - e1_rate - e3_rate ! may have either sign epsilonDot_great = MAX(ABS(e1_rate), ABS(e2_rate), ABS(e3_rate)) ! per Appendix of Carafa & Bird [2016, J. Geophys. Res.] epsilonDot_small = MIN(ABS(e1_rate), ABS(e2_rate), ABS(e3_rate)) ! per Appendix of Carafa & Bird [2016, J. Geophys. Res.] epsilonDot_medium = epsilonDot_great - epsilonDot_small ! by definition (ibid) IF (epsilonDot_great == ABS(err)) THEN ! CF+CF or EF+EF regime (i.e., 2 conjugate sets of thrusts, or 2 conjugate sets of normal faults) IF (err > 0.0D0) THEN ! CF+CF regime (2 conjugate sets of thrusts) alpha_1_degrees = CF_dip_degrees alpha_2_degrees = CF_dip_degrees ELSE ! err < 0.0D0; EF+EF regime (2 conjugate sets of normal faults) alpha_1_degrees = (90.0D0 - EF_dip_degrees) alpha_2_degrees = (90.0D0 - EF_dip_degrees) END IF ELSE ! epsilonDot_great is horizontal, so we have either SS+CF or CF+SS or SS+EF or EF+SS (transpression or transtension). IF (err > 0.0D0) THEN ! case of SS+CF or CF+SS: IF (epsilonDot_medium == ABS(err)) THEN ! case of CF+SS: alpha_1_degrees = CF_dip_degrees alpha_2_degrees = SS_alpha_degrees ELSE ! case of SS+CF: alpha_1_degrees = SS_alpha_degrees alpha_2_degrees = CF_dip_degrees END IF ELSE ! err <= 0.0D0; case of SS+EF or EF+SS: IF (epsilonDot_medium == ABS(err)) THEN ! case of EF+SS: alpha_1_degrees = (90.0D0 - EF_dip_degrees) alpha_2_degrees = SS_alpha_degrees ELSE ! case of SS+EF: alpha_1_degrees = SS_alpha_degrees alpha_2_degrees = (90.0D0 - EF_dip_degrees) END IF END IF END IF alpha_1_radians = alpha_1_degrees * radians_per_degree alpha_2_radians = alpha_2_degrees * radians_per_degree k_alpha_1 = 1.0D0 / (DSIN(alpha_1_radians) * DCOS(alpha_1_radians)) ! dimensionless; greater than or equal to 2.0D0 k_alpha_2 = 1.0D0 / (DSIN(alpha_2_radians) * DCOS(alpha_2_radians)) ! dimensionless; greater than or equal to 2.0D0 IF (region_index == 2) THEN M_persec_per_m3 = IT_assumed_mu_Pa(IT_column) * ((k_alpha_1 * epsilonDot_medium) + (k_alpha_2 * epsilonDot_small)) ELSE M_persec_per_m3 = T5_assumed_mu_Pa(T5_column) * ((k_alpha_1 * epsilonDot_medium) + (k_alpha_2 * epsilonDot_small)) END IF !convert moment rate per unit volume to moment rate per unit area: IF (region_index == 2) THEN M_persec_per_m2 = M_persec_per_m3 * IT_coupledThickness_m(IT_column) * IT_supplement_seismicity_fraction_based_on_z(IT_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) M_persec_per_m2 = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE M_persec_per_m2 = M_persec_per_m3 * T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) M_persec_per_m2 = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF !determine subsampling-point value of seismicity (earthquakes/m**2/s) and add to array: IF (region_index == 2) THEN seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * 1.0D0 * & & M_persec_per_m2 / IT_tGR_model_moment_rate(IT_column) ELSE seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * 1.0D0 * & & M_persec_per_m2 / T5_tGR_model_moment_rate(T5_column) END IF !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: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF !G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & EXP(argument) END IF point_value = G_function * seismicity_at_catalog_threshold seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + weight * point_value IF (threshold_index == 1) continuum_sum = continuum_sum + weight * point_value * grid_cell_areas(i, j) 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 ) END IF ! region_index < 3 (The code above is not needed for SyntheticExample cases.) !============================================================================================================= !ZERO-OUT GRID AREAS COVERED BY OROGENS that were modeled in detail with NeoKinema (version 2.0+): IF (region_index < 3) THEN IF (orogens > 0) THEN 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.0D0 ! 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 DLonLat_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 (DSaIT_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 = DRelative_Compass(from_uvec = uvec, to_uvec = uvec1) t2 = DRelative_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 IF (threshold_index == 1) continuum_sum = continuum_sum - seismicity(i, j, threshold_index) * fraction_inside * grid_cell_areas(i, j) 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 END IF ! region_index < 3 IF (region_index < 3) THEN !============================================================================================================= ! ! ADD PLATE-BOUNDARY SEISMICITY based on gobal plate model PB2002 of Bird [2003, G**3], ! and seismicity calibrations { , m_c } by Bird & Kagan [2004, BSSA] ! and velocity-dependent adjustments to by Bird et al. [2009, 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], ! especially if other researchers are doing the NeoKinema modeling!) 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 DLonLat_2_Uvec(old_longitude, old_latitude, begin_uvec) CALL DLonLat_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 = DRelative_Compass(from_uvec = begin_uvec, to_uvec = uvec1) t2 = DRelative_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 = DRelative_Compass(from_uvec = end_uvec, to_uvec = uvec1) t2 = DRelative_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 this formula to allow for being inside a clockwise circuit would (unacceptably!) ! 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 DCross(first_a_uvec, last_a_uvec, tvec) CALL DMake_Uvec(tvec, pole_a_uvec) dot_a = DDot(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 DCross(first_b_uvec, last_b_uvec, tvec) CALL DMake_Uvec(tvec, pole_b_uvec) dot_b = DDot(first_b_uvec, pole_b_uvec) CALL DCircles_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 = DRelative_Compass(from_uvec = end_uvec, to_uvec = uvec1) t2 = DRelative_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 DCross(first_a_uvec, last_a_uvec, tvec) CALL DMake_Uvec(tvec, pole_a_uvec) dot_a = DDot(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 DCross(first_b_uvec, last_b_uvec, tvec) CALL DMake_Uvec(tvec, pole_b_uvec) dot_b = DDot(first_b_uvec, pole_b_uvec) CALL DCircles_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.5D0 * (begin_uvec + end_uvec) Call DMake_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 DTurn_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 DTurn_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 DTurn_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 DTurn_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(DSIN(subduction_azimuth_1 - (azimuth_integer * radians_per_degree))) lane_width_km = MAX(lane_width_km, 1.0D0) ! 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.0D0 ! 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.0D0 s1 = 0.0D0 - 3.0D0 * alongStep_sigma_km s2 = basis_length_km + 3.0D0 * alongStep_sigma_km ds = (s2 - s1) / 100.D0 DO n = 0, 100 s_km = s1 + (n / 100.0D0) * (s2 - s1) !========== IMSL version: ===================================================================================================== !temp_integral = temp_integral + ds * ANORDF(s_km / alongStep_sigma_km) * ANORDF((basis_length_km - s_km) / alongStep_sigma_km) !========== MKL version: ====================================================================================================== VDCdfNorm_in(1) = s_km / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile1 = VDCdfNorm_out(1) VDCdfNorm_in(1) = (basis_length_km - s_km) / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile2 = VDCdfNorm_out(1) temp_integral = temp_integral + ds * quantile1 * quantile2 !============================================================================================================================== END DO Y_factor_of_step = 1 / temp_integral !relate "class" to "T5_column" (or, "IT_column"): IF (class == "CRB") THEN IF (region_index == 2) THEN IT_column = 1 ELSE T5_column = 1 END IF ELSE IF (class == "CTF") THEN IF (region_index == 2) THEN IT_column = 2 ELSE T5_column = 2 END IF ELSE IF (class == "CCB") THEN IF (velocity_mmpa <= 24.2D0) THEN IF (region_index == 2) THEN IT_column = 3 ! slowCCB ELSE T5_column = 3 ! slowCCB END IF ELSE IF (region_index == 2) THEN IT_column = 4 ! fastCCB ELSE T5_column = 4 ! fastCCB END IF END IF ELSE IF (class == "OSR") THEN IF (region_index == 2) THEN IT_column = 5 ! normalOSR ELSE T5_column = 5 ! normalOSR END IF !N.B. However, coupled thickness will be based on equation of Bird et al. [2002]. !N.B. 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.5D0) THEN IF (region_index == 2) THEN IT_column = 7 ! slowOTF ELSE T5_column = 7 ! slowOTF END IF ELSE IF (velocity_mmpa < 68.5D0) THEN IF (region_index == 2) THEN IT_column = 8 ! mediumOTF ELSE T5_column = 8 ! mediumOTF END IF ELSE IF (region_index == 2) THEN IT_column = 9 ! fastOTF ELSE T5_column = 9 ! fastOTF END IF END IF ELSE IF (class == "OCB") THEN IF (region_index == 2) THEN IT_column = 10 ! OCB ELSE T5_column = 10 ! OCB END IF ELSE IF (class == "SUB") THEN IF (velocity_mmpa <= 66.D0) THEN IF (region_index == 2) THEN IT_column = 11 ! slowSUB ELSE T5_column = 11 ! slowSUB END IF ELSE IF (region_index == 2) THEN IT_column = 12 ! fastSUB ELSE T5_column = 12 ! fastSUB END IF END IF ELSE WRITE (*, "(' ERROR: Invalid plate-boundary class, ',A,' in PB2002_steps.dat.')") TRIM(class) CALL Pause() STOP 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.D0 ! from catalog_PB2002_OSR_pure_normal_seismicity.xls, based on Bird & Kagan [2004] vs_mmpa = 19.D0 ! [ibid] coupledThickness_m = d0_m * DEXP(-divergence_mmpa / vs_mmpa) IF (region_index == 2) THEN IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF ELSE ! all other classes IF (region_index == 2) THEN coupledThickness_m = IT_coupledThickness_m(IT_column) * IT_supplement_seismicity_fraction_based_on_z(IT_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE coupledThickness_m = T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF END IF !seismicity rate of this step, in EQs per second, above the CMT threshold magnitude: IF (region_index == 2) THEN v_oblique_mmpa = DSQRT(dextral_mmpa**2 + (divergence_mmpa / DCOS(IT_assumed_dip_radians(IT_column)))**2) v_oblique_mps = v_oblique_mmpa * (0.001D0 / s_per_year) moment_rate_of_step_SI = coupledThickness_m * & ! coupled lithosphere thickness, in m & IT_assumed_mu_Pa(IT_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/DSIN(IT_assumed_dip_radians(IT_column))) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * & & moment_rate_of_step_SI / IT_tGR_model_moment_rate(IT_column) ELSE v_oblique_mmpa = DSQRT(dextral_mmpa**2 + (divergence_mmpa / DCOS(T5_assumed_dip_radians(T5_column)))**2) v_oblique_mps = v_oblique_mmpa * (0.001D0 / 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/DSIN(T5_assumed_dip_radians(T5_column))) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) END IF !correct seismicity to desired threshold moment/magnitude, !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 (or, table IT_) above: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF !G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & EXP(argument) END IF EQs_per_s(threshold_index) = G_function * seismicity_at_catalog_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 DLonLat_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.987D0) THEN ! within ~1000 km; more precise distance will be computed !distances to end points of this step: uvec1 = begin_uvec to_uvec1_radians = DArc(uvec, uvec1) uvec2 = end_uvec to_uvec2_radians = DArc(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 DCross(uvec1, uvec2, tvec) CALL DMake_Uvec(tvec, pole_uvec) ! pole of the great-circle for this step offside_radians = ABS(DArc(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 = DRelative_Compass(from_uvec = pole_uvec, to_uvec = uvec1) t2 = DRelative_Compass(from_uvec = pole_uvec, to_uvec = uvec2) !note: normally, t2 is slightly less than t1 IF (t2 > t1) t2 = t2 - Two_Pi angle = DRelative_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.0D0 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.0D0 * (offset_km / h_km(k_class_index))**2 IF (argument > -74.0D0) THEN X = (1.0D0 / 0.95D0) *(0.797885D0 / h_km(k_class_index)) * DEXP(argument) ELSE ! prevent underflow in DEXP or elsewhere X = 0.0D0 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.0D0 * alongStep_sigma_km)).AND.(s_km < (length_km + 3.0D0 * alongStep_sigma_km))) THEN ! use complex formula: !======= IMSL version: ==================================================================================== !Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((length_km - s_km) / alongStep_sigma_km) !======= MKL version: ===================================================================================== VDCdfNorm_in(1) = s_km / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile1 = VDCdfNorm_out(1) VDCdfNorm_in(1) = (length_km - s_km) / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile2 = VDCdfNorm_out(1) Y = Y_factor_of_step * quantile1 * quantile2 !=========================================================================================================== ELSE ! too far away; don't risk numerical problems within ANORDF or VDCdfNorm: Y = 0.0D0 END IF !------------------------------------------------------------------------------------------------------- fault_fraction = 1.0D0 continuum_fraction = 0.0D0 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 DCross (uvec1, uvec2, tvec) ! direction toward pole is direction into the lane CALL DMake_Uvec(tvec, inward_uvec) scalar_product = DDot(uvec, inward_uvec) dot_side(n) = scalar_product proximal = proximal .AND. (scalar_product > -0.01D0) ! 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.0D0, 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 DMake_Uvec(tvec, trench_uvec) ! at same relative left-right position in this lane distance_km = radius_km * DArc(uvec, trench_uvec) ! always zero or positive; if positive, must apply sign now: IF (distance_km > 0.0D0) THEN from_trench_azimuth_radians = DRelative_Compass (from_uvec = trench_uvec, to_uvec = uvec) IF (DCOS(from_trench_azimuth_radians - subduction_azimuth_1) < 0.0D0) 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.0D0) THEN ! landward side of seismicity peak: argument = -0.5D0 * (offset_km / SUB_landward_sigma_km)**2 IF (argument > -74.0D0) THEN X = (1.0D0 / 0.95D0) * & & (2.0D0 * SUB_landward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885D0 / (2.0D0 * SUB_landward_sigma_km)) * & & DEXP(argument) ELSE X = 0.0D0 END IF ELSE ! seaward side of seismicity peak (but may not be seaward of the trench!) argument = -0.5D0 * (offset_km / SUB_seaward_sigma_km)**2 IF (argument > -74.0D0) THEN X = (1.0D0 / 0.95D0) * & & (2.0D0 * SUB_seaward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885D0 / (2.0D0 * SUB_seaward_sigma_km )) * & & DEXP(argument) ELSE X = 0.0D0 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.0D0) ! 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.0D0 * alongStep_sigma_km)).AND.(s_km < (lane_width_km + 3.0D0 * alongStep_sigma_km))) THEN ! use complex formula: !======= IMSL version: ==================================================================================== !Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((lane_width_km - s_km) / alongStep_sigma_km !======= MKL version: ===================================================================================== VDCdfNorm_in(1) = s_km / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile1 = VDCdfNorm_out(1) VDCdfNorm_in(1) = (lane_width_km - s_km) / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile2 = VDCdfNorm_out(1) Y = Y_factor_of_step * quantile1 * quantile2 !=========================================================================================================== ELSE ! too far away; don't risk numerical problems within ANORDF or VDCdfNorm: Y = 0.0D0 END IF !-------------------------------------------------------------------------------------------------------- ELSE ! not proximal; punt X = 0.0D0 Y = 0.0D0 END IF fault_fraction = (1.0D0 - SUB_bending_fraction) continuum_fraction = SUB_bending_fraction !See caption to Figure 13 in Bird et al. [2009, BSSA, 99(6), 3097-3113, doi:10.1785/0120090082] 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 IF (threshold_index == 1) THEN fault_sum = fault_sum + fault_fraction * weight * B * EQs_per_s(threshold_index) * grid_cell_areas(i, j) continuum_sum = continuum_sum + fault_fraction * weight * B * EQs_per_s(threshold_index) * grid_cell_areas(i, j) END IF 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 !--------------------------------------------------------------------------------------------------------- ! NONLINEAR STEP (which compromises the breakdown-by-columns attempted under the "subdivide = .TRUE" ! debugging option, leading to output sum_of_fractions > 1.00 outside the NeoKinema model(s), ! and also in marginal parts of the NeoKinema model(s) that are affected by lateral-diffusion): DO i = 1, rows DO j = 1, columns IF (seismicity(i, j, 1) >= PB2002_seismicity(i, j, 1)) THEN !Earth5-based continuum seismicity is greater than tails of PB2002 fault-based seismicity; !continuum seismicity will set the total level, and fault-based tail will be dropped. !Adjust fractional sums accordingly: fault_sum = fault_sum - PB2002_seismicity(i, j, 1) * grid_cell_areas(i, j) ELSE !PB2002-based fault seismicity is greater than Earth5-based continuum seismicity; !PB2002 will set the level, and continuum seismicity will be dropped. !Adjust fractional sums accordingly: continuum_sum = continuum_sum - seismicity(i, j, 1) * grid_cell_areas(i, j) END IF 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 END IF ! region_index < 3 (The code above is not needed for case 3: SyntheticExample.) !============================================================================================================= IF (orogens > 0) THEN IF (region_index < 3) THEN ! smoothing will be applied around orogen edges. We need smoothing arrays: IF (.NOT.ALLOCATED(diffusivity)) 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. IF (.NOT.ALLOCATED(single_column)) ALLOCATE ( single_column(rows) ) ! temporary storage needed for ADE diffusion algorithm IF (.NOT.ALLOCATED(single_row)) ALLOCATE ( single_row(columns) ) END IF ! region_index < 3 DO orogen = 1, orogens ! ADD SEISMICITY DUE TO DISTRIBUTED PERMANENT ("off-[modeled-]fault") STRAIN-RATES in a NeoKinema model of an orogen !Read the .feg file and store it (temporarily) in memory... 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) IF (ALLOCATED(node_uvec)) DEALLOCATE ( node_uvec ) 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 DLonLat_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) IF (ALLOCATED(nodes)) DEALLOCATE ( nodes ) 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) IF (.NOT.ALLOCATED(element_principal_strainrates)) ALLOCATE ( element_principal_strainrates(3, numel) ) element_principal_strainrates = 0.0D0 ! just in case some element never appears in the e_nko file(?) 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 DPrincipal_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, i) = e1h element_principal_strainrates(2, i) = e2h element_principal_strainrates(3, i) = err !NOTE that LOGICAL column "maybe" is not memorized or used--So, it need not be accurate. ! (Actually, this column is only used by plotting program NeoKineMap.) END DO ! l_ = 1, numel, reading continuum strainrates and finding principal values CLOSE(5) !Transfer the seismicity values to the .grd points (using subsampling): IF (.NOT.ALLOCATED(a_)) ALLOCATE ( a_(numel) ) IF (.NOT.ALLOCATED(center)) ALLOCATE ( center(3, numel) ) IF (.NOT.ALLOCATED(neighbor)) ALLOCATE ( neighbor(3, numel) ) CALL DLearn_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 DLonLat_2_Uvec(t_lon, t_lat, uvec) !recall from memory: uvec(1:3) = subsample_uvec(1:3, ii, jj, i, j) CALL DWhich_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 (or table IT_): IF (continental(i, j)) THEN IF (err == 0.0D0) THEN ! rare: ideal strike-slip IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF ENDIF ELSE IF (err > 0.0D0) THEN ! thrusting (possibly minor?) IF (err > (allowed_SS_impurity * e2h)) THEN ! primarily thrusting IF (region_index == 2) THEN IT_column = 3 ! slow CCB ELSE T5_column = 3 ! slow CCB END IF 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] IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (allowed_SS_impurity * e1h)) THEN ! primarily normal-faulting IF (region_index == 2) THEN IT_column = 1 ! CRB ELSE T5_column = 1 ! CRB END IF 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] IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF END IF END IF ! thrusting or normal-faulting? ELSE ! oceanic IF (err == 0.0D0) THEN ! rare: ideal strike-slip IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF ELSE IF (err > 0.0D0) THEN ! thrusting (possibly minor?) IF (err > (allowed_SS_impurity * e2h)) THEN ! primarily thrusting IF (region_index == 2) THEN IT_column = 10 ! OCB ELSE T5_column = 10 ! OCB END IF 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] IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF END IF ELSE ! err < 0.0; normal-faulting (possibly minor?) IF (err < (allowed_SS_impurity * e1h)) THEN ! primarily normal-faulting IF (region_index == 2) THEN IT_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 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.] END IF 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] IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF 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, !per Appendix of Carafa & Bird [2016, J. Geophys. Res.-Solid Earth]: epsilonDot_great = MAX(ABS(e1h), ABS(e2h), ABS(err)) epsilonDot_small = MIN(ABS(e1h), ABS(e2h), ABS(err)) epsilonDot_medium = epsilonDot_great - epsilonDot_small IF (epsilonDot_great == ABS(err)) THEN ! CF+CF or EF+EF regime (i.e., 2 conjugate sets of thrusts, or 2 conjugate sets of normal faults) IF (err > 0.0D0) THEN ! CF+CF regime (2 conjugate sets of thrusts) alpha_1_degrees = CF_dip_degrees alpha_2_degrees = CF_dip_degrees ELSE ! err < 0.0D0; EF+EF regime (2 conjugate sets of normal faults) alpha_1_degrees = (90.0D0 - EF_dip_degrees) alpha_2_degrees = (90.0D0 - EF_dip_degrees) END IF ELSE ! epsilonDot_great is horizontal, so we have either SS+CF or CF+SS or SS+EF or EF+SS (transpression or transtension). IF (err > 0.0D0) THEN ! case of SS+CF or CF+SS: IF (epsilonDot_medium == ABS(err)) THEN ! case of CF+SS: alpha_1_degrees = CF_dip_degrees alpha_2_degrees = SS_alpha_degrees ELSE ! case of SS+CF: alpha_1_degrees = SS_alpha_degrees alpha_2_degrees = CF_dip_degrees END IF ELSE ! err <= 0.0D0; case of SS+EF or EF+SS: IF (epsilonDot_medium == ABS(err)) THEN ! case of EF+SS: alpha_1_degrees = (90.0D0 - EF_dip_degrees) alpha_2_degrees = SS_alpha_degrees ELSE ! case of SS+EF: alpha_1_degrees = SS_alpha_degrees alpha_2_degrees = (90.0D0 - EF_dip_degrees) END IF END IF END IF alpha_1_radians = alpha_1_degrees * radians_per_degree alpha_2_radians = alpha_2_degrees * radians_per_degree k_alpha_1 = 1.0D0 / (DSIN(alpha_1_radians) * DCOS(alpha_1_radians)) ! dimensionless; greater than or equal to 2.0D0 k_alpha_2 = 1.0D0 / (DSIN(alpha_2_radians) * DCOS(alpha_2_radians)) ! dimensionless; greater than or equal to 2.0D0 IF (region_index == 2) THEN M_persec_per_m3 = IT_assumed_mu_Pa(IT_column) * ((k_alpha_1 * epsilonDot_medium) + (k_alpha_2 * epsilonDot_small)) ELSE M_persec_per_m3 = T5_assumed_mu_Pa(T5_column) * ((k_alpha_1 * epsilonDot_medium) + (k_alpha_2 * epsilonDot_small)) END IF !convert moment rate per unit volume to moment rate per unit area: IF (region_index == 2) THEN M_persec_per_m2 = M_persec_per_m3 * IT_coupledThickness_m(IT_column) * IT_supplement_seismicity_fraction_based_on_z(IT_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) M_persec_per_m2 = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE M_persec_per_m2 = M_persec_per_m3 * T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) M_persec_per_m2 = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF !determine subsampling-point value of seismicity (earthquakes/m**2/s) and add to array: IF (region_index == 2) THEN seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * & & M_persec_per_m2 / IT_tGR_model_moment_rate(IT_column) ELSE seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * & & M_persec_per_m2 / T5_tGR_model_moment_rate(T5_column) END IF !correct seismicity to desired threshold moment/magnitude, !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 above: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF ! G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & DEXP(argument) END IF point_value = G_function * seismicity_at_catalog_threshold seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + weight * point_value IF (threshold_index == 1) continuum_sum = continuum_sum + weight * point_value * grid_cell_areas(i, j) 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 ) !============================================================================================================= IF (region_index < 3) THEN ! lateral smoothing is used. (It is NOT used for 3: SyntheticExample.) ! 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.0D0 ! whole array (rows, columns). Some values will be increased below... peak_diffusivity = 0.0D0 ! just initializing; will be increased... WRITE (*,"(' Computing gridded diffusivities of seismicity around boundary of this orogen...')") IF (region_index <= 1) THEN ! Global/default, or Western US settings: diffusion_scale = 100.D0 / 6371.D0 ! (in km/km = radians). Note that this was 100 km in LTS_v3. ELSE IF (region_index == 2) THEN ! Italia settings: diffusion_scale = 30.D0 / 6371.D0 ! (in km/km = radians). Note that this was 100 km in LTS_v3. END IF !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 = DArc(uvec1, uvec2) scales = distance / diffusion_scale ! both in radians IF (scales < 3.D0) THEN ! discard values less than 0.0001 diffusivity(i, j) = MAX(diffusivity(i, j), DEXP(-scales**2)) peak_diffusivity = MAX(peak_diffusivity, diffusivity(i, j)) 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.0D0) < (0.5D0 * 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.0D0 - d_lon)) < (0.5D0 * 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.0D0) THEN nr1 = 1 ! diffuse all the way up to Northernmost row ELSE nr1 = 1 + NINT( (lat_max - 75.0D0) / d_lat ) nr1 = MIN(nr1, rows) END IF IF (lat_min >= -75.0D0) THEN nr2 = rows ELSE nr2 = 1 + NINT( (lat_max - (-75.0D0)) / 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 * DCOS(MIN(lat_max, +75.0D0) * radians_per_degree)), & & (d_lon * radians_per_degree * DCOS(MAX(lat_min, -75.0D0) * radians_per_degree)) ) d_S_radians = d_lat * radians_per_degree IF (peak_diffusivity > 0.0D0) THEN safe_timestep = 0.1D0 * minimum_arc**2 / peak_diffusivity ELSE ! safe timestep is irrelevant; there is no need for timesteps. safe_timestep = 1.0D0 ! (an arbitrary large number) END IF total_time = diffusion_scale**2 ! again, based on diffusivity of 1 along margins of this orogen 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.0D0) 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 * DCOS(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.0D0) 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.0D0) 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 END IF ! region_index < 3, so smoothing is applied to orogen edges. !============================================================================================================= IF (faults_used_in_orogen(orogen)) THEN !ADD SEISMICITY DUE TO FAULT SEGMENTS in a NeoKinema (version 2.0+) models of orogen !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 ! empty file, bad format, (e.g., wrong h_.nko 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 (*, *) WRITE (*, "(' Do you wish to SKIP this step, and compute WITHOUT FAULTS? [T/F]:')") READ (*, *) skip_faults IF (skip_faults) THEN EXIT segmenting ELSE 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 END IF ! skip_faults, or not 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 DLonLat_2_Uvec(lon1, lat1, uvec1) segment_uvecs(1:3, 1, segment) = uvec1(1:3) CALL DLonLat_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) !Decide on dip angles (in degrees, 1~90) to be used for each segment: IF (.NOT.ALLOCATED(segment_dip_degrees)) ALLOCATE ( segment_dip_degrees(segments) ) ! INTEGER vector !Set whole vector to zero (so I can see which values are actually known!): segment_dip_degrees = 0 !Set up storage for known dip angles (indexed by F0001R number, NOT by segment number): highest_fault_number = 0 ! just initializing, before search DO i = 1, segments c4 = segment_c6(i)(2:5) ! isolating "1234" from "F1234R", for example READ (c4, *) this_fault_number highest_fault_number = MAX(highest_fault_number, this_fault_number) END DO IF (.NOT.ALLOCATED(fault_dip_degrees)) ALLOCATE ( fault_dip_degrees(highest_fault_number) ) ! INTEGER vector fault_dip_degrees = 0 ! whole array, to clarify which have actually been read from f___.dig !Read f_____.dig file and look for "dip_degrees" and memorize them: 321 OPEN (UNIT = 11, FILE = TRIM(orogen_f_dig_pathfile(orogen)), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Required file ',A,' is now missing. Replace it!')") TRIM(orogen_f_dig_pathfile(orogen)) CALL Pause GO TO 321 END IF mining_for_dips: DO READ (11,"(A)", IOSTAT = ios) line IF (ios == -1) EXIT mining_for_dips ! EOF first_byte = line(1:1) IF ((first_byte == 'F').OR.(first_byte == 'f')) THEN READ (line,"(1X,I4)") trace_index got_dip_degrees = .FALSE. ! until... ELSE ! line was not a title; this leaves 3 possibilities: dip_degrees, (lon, lat), or *** end !try to read "dip_degrees"... BACKSPACE (11) ! instead of READ(line,*,...) which is BUGGY! READ (11, "(A)", IOSTAT = internal_ios) line start_loc = INDEX(line, "dip_degrees") IF (start_loc > 0) THEN ! found "dip_degrees" in line beyond_loc = start_loc + 11 ! first byte which is not part of "dip_degrees" end_loc = LEN_TRIM(line) line = line(beyond_loc:end_loc) READ (line, *, IOSTAT = internal_ios) dip_degrees IF (internal_ios == 0) THEN ! only use numbers that were read successfully! IF (trace_index <= highest_fault_number) THEN ! only record dips of F#s contained in our working array! fault_dip_degrees(trace_index) = dip_degrees END IF END IF END IF END IF END DO mining_for_dips CLOSE (11) ! orogen_f_dig_pathfile(orogen) !Transfer known dip angles from fault_dip_degrees to segment_dip_degrees: DO i = 1, segments c4 = segment_c6(i)(2:5) ! isolating "1234" from "F1234R", for example READ (c4, *) this_fault_number IF (fault_dip_degrees(this_fault_number) > 0) segment_dip_degrees(i) = fault_dip_degrees(this_fault_number) END DO !Note that unknown (generic) segment dips are still left as zeros, at this point. !Find segment_outer_uvecs (to permit computing mean azimuths). !NOTE that subduction zones that are contiguous will be recognized as "one fault" for this purpose ! even if the FnnnnS number is different for an oblique-subduction section ! (to give a more accurate NeoKinema solution, by limiting the along-strike averaging of strike-slip). IF (.NOT.ALLOCATED(segment_used)) ALLOCATE ( segment_used(segments) ) IF (.NOT.ALLOCATED(segment_outer_uvecs)) ALLOCATE ( segment_outer_uvecs(3, 2, segments) ) IF (.NOT.ALLOCATED(fault_azimuth_radians)) ALLOCATE ( fault_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)).OR. & & ((segment_c6(segment)(6:6) == 'S').AND.(segment_c6(other_segment)(6:6) == 'S'))) 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 (DSaIT_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)).OR. & & ((segment_c6(segment)(6:6) == 'S').AND.(segment_c6(other_segment)(6:6) == 'S'))) 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 (DSaIT_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 (DSaIT_Uvec(uvec1, uvec2)) THEN !null, zero-length segment chain fault_azimuth_radians(segment) = 0.0D0 ELSE fault_azimuth_radians(segment) = DRelative_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')) & ! N.B. "S" Subduction zones generate plate-bending EQs, even when creeping! & ) 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.5D0 * (begin_uvec + end_uvec) Call DMake_Uvec(tvec, center_uvec) !compute geometric properties of NeoKinema segment (using same variable names as in SUB code above): length_km = radius_km * DArc(begin_uvec, end_uvec) IF (length_km <= 0.0D0) CYCLE add_segment_seismicity ! to prevent crash in Compass azimuth_integer = NINT(degrees_per_radian * DCompass(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 seismicity footprint !as subduction zones in PB2002 (both per Bird & Kagan [2004, BSSA]). k_class_index = 7 ! SUB IF (segment_slip_rate_mmpa(segment) <= 66.D0) THEN IF (region_index == 2) THEN IT_column = 11 ! slow SUB ELSE T5_column = 11 ! slow SUB END IF ELSE IF (region_index == 2) THEN IT_column = 12 ! fast SUB ELSE T5_column = 12 ! fast SUB END IF END IF IF (region_index < 3) THEN ! real-Earth case; consult global plate model: !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.0D0 ! 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 DLonLat_2_Uvec(old_longitude, old_latitude, uvec1) CALL DLonLat_2_Uvec(longitude, latitude, uvec2) tvec = 0.5D0 * (uvec1 + uvec2) CALL DMake_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 = (DDot(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 moving W w.r.t. right plate ! c5 = "CT/SF" ! so left plate (as you go S) is overriding ELSE ! region_index == 3; Subduction zone is hypothetical/ideal; ! assume simple orthogonal subduction {i.e., orthogonal to overall trend of trench} to define lanes: !recall the outer limits of linked segments (that were found in code above): uvec1(1:3) = segment_outer_uvecs(1:3, 1, segment) uvec2(1:3) = segment_outer_uvecs(1:3, 2, segment) !and use the mean azimuth of the whole trench in the following formula. !Note that velocity_azimuth_integer refers to the relative velocity of the over-riding plate: velocity_azimuth_integer = NINT(90.0D0 + degrees_per_radian * DRelative_Compass(from_uvec = uvec1, to_uvec = uvec2)) c5 = "OP/SP" !Overriding-plate/Subducting-plate; right plate subducts (according to NeoKinema and Restore f_.dig conventions). END IF ! real-Earth or SyntheticExample Subduction zone !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 DTurn_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 DTurn_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 DTurn_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 DTurn_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(DSIN(subduction_azimuth_1 - (azimuth_integer * radians_per_degree))) lane_width_km = MAX(lane_width_km, 1.0D0) ! 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.0D0 s1 = 0.0D0 - 3.0D0 * alongStep_sigma_km s2 = basis_length_km + 3.0D0 * alongStep_sigma_km ds = (s2 - s1) / 100.D0 DO n = 0, 100 s_km = s1 + (n / 100.0D0) * (s2 - s1) !======= IMSL version: ==================================================================================== !temp_integral = temp_integral + ds * ANORDF(s_km / alongStep_sigma_km) * ANORDF((basis_length_km - s_km) / alongStep_sigma_km) !======= MKL version: ===================================================================================== VDCdfNorm_in(1) = s_km / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile1 = VDCdfNorm_out(1) VDCdfNorm_in(1) = (basis_length_km - s_km) / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile2 = VDCdfNorm_out(1) temp_integral = temp_integral + ds * quantile1 * quantile2 !=========================================================================================================== END DO Y_factor_of_step = 1.0D0 / 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.001D0 / s_per_year) IF (region_index == 2) THEN moment_rate_of_step_SI = IT_coupledThickness_m(IT_column) * & ! coupled lithosphere thickness, in m & IT_supplement_seismicity_fraction_based_on_z(IT_column) * & & IT_assumed_mu_Pa(IT_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/DSIN(IT_assumed_dip_radians(IT_column))) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) moment_rate_of_step_SI = 0.0 ! Set c = 0., preventing any M-rate or seismicity. seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * & & moment_rate_of_step_SI / IT_tGR_model_moment_rate(IT_column) ELSE 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/DSIN(T5_assumed_dip_radians(T5_column))) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) moment_rate_of_step_SI = 0.0 ! Set c = 0., preventing any M-rate or seismicity. seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) END IF !correct seismicity to desired threshold moment/magnitude, !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 above: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF ! G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & DEXP(argument) END IF EQs_per_s(threshold_index) = G_function * seismicity_at_catalog_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 DLonLat_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.987D0) THEN ! within ~1000 km; more precise distance will be computed !distances to end points of this step: uvec1 = begin_uvec to_uvec1_radians = DArc(uvec, uvec1) uvec2 = end_uvec to_uvec2_radians = DArc(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 DCross (uvec1, uvec2, tvec) ! direction toward pole is direction into the lane CALL DMake_Uvec(tvec, inward_uvec) scalar_product = DDot(uvec, inward_uvec) dot_side(n) = scalar_product proximal = proximal .AND. (scalar_product > -0.01D0) ! 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.0D0, MIN(1.0D0, 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 DMake_Uvec(tvec, trench_uvec) ! at same relative left-right position in this lane distance_km = radius_km * DArc(uvec, trench_uvec) ! always zero or positive; if postive, must apply sign now: IF (distance_km > 0.0D0) THEN from_trench_azimuth_radians = DRelative_Compass (from_uvec = trench_uvec, to_uvec = uvec) IF (DCOS(from_trench_azimuth_radians - subduction_azimuth_1) < 0.0D0) 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.0D0) THEN ! landward side of seismicity peak: argument = -0.5D0 * (offset_km / SUB_landward_sigma_km)**2 IF (argument > -74.0D0) THEN X = (1.0D0 / 0.95D0) * & & (2.D0 * SUB_landward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885D0 / (2.0D0 * SUB_landward_sigma_km)) * & & DEXP(argument) ELSE X = 0.0D0 END IF ELSE ! seaward side of seismicity peak (but may not be seaward of the trench!) argument = -0.5D0 * (offset_km / SUB_seaward_sigma_km)**2 IF (argument > -74.0D0) THEN X = (1.0D0 / 0.95D0) * & & (2.D0 * SUB_seaward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885D0 / (2.0D0 * SUB_seaward_sigma_km )) * & & DEXP(argument) ELSE X = 0.0D0 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.0D0) ! 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.0D0 * alongStep_sigma_km)).AND.(s_km < (lane_width_km + 3.0D0 * alongStep_sigma_km))) THEN ! use complex formula: !======= IMSL version: ==================================================================================== !Y = Y_factor_of_step * ANORDF(s_km / alongStep_sigma_km) * ANORDF((lane_width_km - s_km) / alongStep_sigma_km) !======= MKL version: ===================================================================================== VDCdfNorm_in(1) = s_km / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile1 = VDCdfNorm_out(1) VDCdfNorm_in(1) = (lane_width_km - s_km) / alongStep_sigma_km CALL VDCdfNorm(I4_one, VDCdfNorm_in, VDCdfNorm_out) quantile2 = VDCdfNorm_out(1) Y = Y_factor_of_step * quantile1 * quantile2 !=========================================================================================================== ELSE ! too far away; don't risk numerical problems within ANORDF or VDCdfNorm: Y = 0.0D0 END IF !-------------------------------------------------------------------------------------------------------- ELSE ! not proximal; punt X = 0.0D0 Y = 0.0D0 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) IF (threshold_index == 1) THEN continuum_sum = continuum_sum + SUB_bending_fraction * weight * B * EQs_per_s(threshold_index) * grid_cell_areas(i, j) END IF ELSE ! normal, mostly-coupled subduction: seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + & & weight * B * EQs_per_s(threshold_index) IF (threshold_index == 1) THEN fault_sum = fault_sum + (1.0D0 - SUB_bending_fraction) * weight * B * EQs_per_s(threshold_index) * grid_cell_areas(i, j) continuum_sum = continuum_sum + SUB_bending_fraction * weight * B * EQs_per_s(threshold_index) * grid_cell_areas(i, j) END IF 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 fault_fraction = 1.0D0 continuum_fraction = 0.0D0 !N.B. begin_uvec, end_uvec, center_uvec, length_km, azimuth_integer, azimuth_radians, segment_vec, ! & fault_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 DUvec_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.0D0 IF (delta_lon_from_grid_center < -180.D0) delta_lon_from_grid_center = delta_lon_from_grid_center + 360.D0 IF (delta_lon_from_grid_center > +180.D0) delta_lon_from_grid_center = delta_lon_from_grid_center - 360.D0 lon = delta_lon_from_grid_center + (lon_min + lon_max)/2.0D0 ! 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 generic/default fault dip to assume (consistent with NeoKinema, & Bird and Kagan [2004]): c1 = segment_c6(segment)(6:6) strike_slip = (c1 == 'L').OR.(c1 == 'R').OR.(c1 == 'l').OR.(c1 == 'r') IF (strike_slip) THEN IF (oceanic) THEN k_class_index = 5 ! OTF IF (segment_slip_rate_mmpa(segment) < 39.5D0) THEN IF (region_index == 2) THEN IT_column = 7 ! slow OTF ELSE T5_column = 7 ! slow OTF END IF ELSE IF (segment_slip_rate_mmpa(segment) < 68.5D0) THEN IF (region_index == 2) THEN IT_column = 8 ! medium OTF ELSE T5_column = 8 ! medium OTF END IF ELSE IF (region_index == 2) THEN IT_column = 9 ! fast OTF ELSE T5_column = 9 ! fast OTF END IF END IF ELSE ! continental k_class_index = 2 ! CTF IF (region_index == 2) THEN IT_column = 2 ! CTF ELSE T5_column = 2 ! CTF END IF END IF ELSE IF ((c1 == 'N').OR.(c1 == 'D').OR.(c1 == 'n').OR.(c1 == 'd')) THEN IF (oceanic) THEN k_class_index = 4 ! OSR IF (region_index == 2) THEN IT_column = 5 ! OSR normal ELSE T5_column = 5 ! OSR normal END IF !NB: However, coupled thickness will be based on equation of Bird et al. [2002]. ELSE ! continental k_class_index = 3 ! CRB IF (region_index == 2) THEN IT_column = 1 ! CRB ELSE T5_column = 1 ! CRB END IF END IF ELSE IF ((c1 == 'T').OR.(c1 == 'P').OR.(c1 == 't').OR.(c1 == 'p')) THEN IF (oceanic) THEN k_class_index = 6 ! OCB IF (region_index == 2) THEN IT_column = 10 ! OCB ELSE T5_column = 10 ! OCB END IF ELSE ! continental k_class_index = 1 ! CCB IF (segment_slip_rate_mmpa(segment) <= 24.2D0) THEN IF (region_index == 2) THEN IT_column = 3 ! slow CCB ELSE T5_column = 3 ! slow CCB END IF ELSE IF (region_index == 2) THEN IT_column = 4 ! fast CCB ELSE T5_column = 4 ! fast CCB END IF END IF END IF ELSE WRITE (*, *) WRITE (*, "(' ERROR: Illegal byte ',A,' in fault descriptor ',A)") c1, segment_c6(segment) CALL Pause() STOP END IF IF (region_index == 2) THEN dip_degrees = IT_assumed_dip_degrees(IT_column) ! default or generic dip, based on L/R/N/D/P/T ELSE dip_degrees = T5_assumed_dip_degrees(T5_column) ! default or generic dip, based on L/R/N/D/P/T END IF IF (segment_dip_degrees(segment) > 0) THEN ! geologic dip is known IF ((.NOT.strike_slip).OR.(segment_dip_degrees(segment) < dip_degrees)) THEN ! geologic dip is safe to use dip_degrees = segment_dip_degrees(segment) ! use actual value, if known & safe. END IF END IF dip_radians = dip_degrees * radians_per_degree IF (region_index == 2) THEN lithosphere_thickness_km = IT_assumed_lithosphere_km(IT_column) ! used to compute width of parallelogram ELSE lithosphere_thickness_km = T5_assumed_lithosphere_km(T5_column) ! used to compute width of parallelogram END IF !===== START OF: Create and apply dippling rectangular patch to the L of segment trace ====================== !Locate corner points of seismicity rectangle (on Earth's surface): cross_trace_km = lithosphere_thickness_km / TAN(dip_radians) ! where 0 < dip_radians < Pi_over_2, so TAN() > 0.0D0 IF (strike_slip) THEN !dip sense is uncertain, so extend parallelograms to BOTH sides of the trace, using mean fault azimuth: CALL DTurn_To (azimuth_radians = fault_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 DTurn_To (azimuth_radians = fault_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 DTurn_To (azimuth_radians = fault_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 DTurn_To (azimuth_radians = fault_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.0D6 * length_km * 2.0D0 * cross_trace_km * & & ABS(DCOS(fault_azimuth_radians(segment) - (azimuth_integer * radians_per_degree))) !on the surface of the planet. !Note: fault_azimuth_radians(segment) is the mean azimuth of the whole fault; ! while azimuth_radians 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 DTurn_To (azimuth_radians = azimuth_radians - Pi_over_2, & & base_uvec = end_uvec, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec3) CALL DTurn_To (azimuth_radians = azimuth_radians - 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.0D6 * length_km * cross_trace_km !on the surface of the planet. !Note: fault_azimuth_radians(segment) is the mean azimuth of the whole fault; ! while azimuth_radians is the particular azimuth of this one segment. END IF !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], or table IT_; OSR includes later cz(v) model): IF (region_index == 2) THEN IF (IT_column == 5) THEN ! OSR normal d0_m = 1480.D0 ! from catalog_PB2002_OSR_pure_normal_seismicity.xls, based on Bird & Kagan [2004] vs_mmpa = 19.D0 ! [ibid] coupledThickness_m = d0_m * DEXP(-segment_heave_rate_mmpa(segment) / vs_mmpa) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE ! all other classes coupledThickness_m = IT_coupledThickness_m(IT_column) * IT_supplement_seismicity_fraction_based_on_z(IT_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= IT_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF ELSE IF (T5_column == 5) THEN ! OSR normal d0_m = 1480.D0 ! from catalog_PB2002_OSR_pure_normal_seismicity.xls, based on Bird & Kagan [2004] vs_mmpa = 19.D0 ! [ibid] coupledThickness_m = d0_m * DEXP(-segment_heave_rate_mmpa(segment) / vs_mmpa) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. ELSE ! all other classes coupledThickness_m = T5_coupledThickness_m(T5_column) * T5_supplement_seismicity_fraction_based_on_z(T5_column) IF (subdivide .AND. (highlight > 0) .AND. (highlight /= T5_column)) coupledThickness_m = 0.0 ! Set c = 0., preventing any M-rate or seismicity. END IF 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.001D0 / s_per_year) IF (region_index == 2) THEN moment_rate_of_step_SI = coupledThickness_m * & ! coupled lithosphere thickness, in m & IT_assumed_mu_Pa(IT_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (1.D0/DSIN(dip_radians)) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * & & moment_rate_of_step_SI / IT_tGR_model_moment_rate(IT_column) ELSE 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/DSIN(dip_radians)) * & ! csc(theta) & (length_km * 1000.D0) ! along-strike length, in m seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) END IF !correct seismicity to desired threshold moment/magnitude, !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 above: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF ! G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & DEXP(argument) END IF EQs_per_s(threshold_index) = G_function * seismicity_at_catalog_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 DCross(uvec1, uvec2, tvec) CALL DMake_Uvec(tvec, pole_uvec1) CALL DCross(uvec2, uvec3, tvec) CALL DMake_Uvec(tvec, pole_uvec2) CALL DCross(uvec3, uvec4, tvec) CALL DMake_Uvec(tvec, pole_uvec3) CALL DCross(uvec4, uvec1, tvec) CALL DMake_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 DLonLat_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.9967D0) 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.0D0) THEN dot_test = uvec(1) * pole_uvec2(1) + uvec(2) * pole_uvec2(2) + uvec(3) * pole_uvec2(3) IF (dot_test > 0.0D0) THEN dot_test = uvec(1) * pole_uvec3(1) + uvec(2) * pole_uvec3(2) + uvec(3) * pole_uvec3(3) IF (dot_test > 0.0D0) THEN dot_test = uvec(1) * pole_uvec4(1) + uvec(2) * pole_uvec4(2) + uvec(3) * pole_uvec4(3) IF (dot_test >= 0.0D0) 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) IF (threshold_index == 1) THEN fault_sum = fault_sum + fault_fraction * weight * EQs_per_s_per_m2(threshold_index) * grid_cell_areas(i, j) continuum_sum = continuum_sum + continuum_fraction * weight * EQs_per_s_per_m2(threshold_index) * grid_cell_areas(i, j) END IF 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 OF: Create and apply dippling rectangular patch to the L of segment trace ====================== !========START OF: Create(?) and apply triangular gore of seismicity (before "next" linked segment) =========== !Determine whether there is any "next" segment, of same fault, linked at the end_uvec end? check_others2: DO other_segment = 1, segments IF (other_segment /= segment) THEN IF (segment_c6(other_segment) == segment_c6(segment)) THEN ! other segment belongs to same fault other_begin_uvec(1:3) = segment_uvecs(1:3, 1, other_segment) ! head of other_segment other_end_uvec(1:3) = segment_uvecs(1:3, 2, other_segment) ! end of other segment ! check that "other" segment has non-zero length (because the h_{token}.nko report has limited precision!): IF (.NOT.DSaIT_Uvec(other_begin_uvec, other_end_uvec)) THEN ! OK to proceed, investigating this "other" segment: ! check for match of present tail with other head? IF (DSaIT_Uvec(end_uvec, other_begin_uvec)) THEN !head of other_segment matches tail of current segment !check for significant angle between these 2 adjacent segments? other_azimuth_radians = DCompass(from_uvec = other_begin_uvec, to_uvec = other_end_uvec) IF (ABS(SIN(other_azimuth_radians - azimuth_radians)) >= 0.0349065D0) THEN ! a significant bend (>= 2 degrees) occurs at this joint between segments !Create and apply a triangular patch of seismicity using the Add/Subtract factor SIGN: !NOTE a fairly arbitrary convention: The triangular patch inherits its values of: ! oceanic, continental, k_class_index, T5_column OR IT_column, ! lithosphere_thickness_km, & coupledThickness_m, ! from the current "segment", not from any average with properties of the upcoming "other_segment". ! We also continue to use the same values of: ! dip_degrees, dip_radians, & cross_trace_km ! but this is not an approximation, because dip is constant along any NeoKinema fault. IF (SIN(other_azimuth_radians - azimuth_radians) > 0.0D0) THEN ! bend is to the right; gore is open; ADD a triangle of seismicity SIGN = +1.0D0 ! Dip is always to left of trace, by NeoKinema f.dig convention. ! The triangle is ALWAYS outlined in counterclockwise order (as seen from above): uvec1 = end_uvec ! of current "segment"; this is the shallow vertex of the triangle CALL DTurn_To (azimuth_radians = other_azimuth_radians - Pi_over_2, & & base_uvec = uvec1, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec2) CALL DTurn_To (azimuth_radians = azimuth_radians - Pi_over_2, & & base_uvec = uvec1, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec3) ELSE ! bend is the the left; gore has double-counted seismicity already; SUBTRACT a triangle of seismicity SIGN = -1.0D0 ! Dip is always to left of trace, by NeoKinema f.dig convention. ! The triangle is ALWAYS outlined in counterclockwise order (as seen from above): uvec1 = end_uvec ! of current "segment"; this is the shallow vertex of the triangle CALL DTurn_To (azimuth_radians = azimuth_radians - Pi_over_2, & & base_uvec = uvec1, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec2) CALL DTurn_To (azimuth_radians = other_azimuth_radians - Pi_over_2, & & base_uvec = uvec1, & & far_radians = cross_trace_km / radius_km, & ! inputs & omega_uvec = omega_uvec, & & result_uvec = uvec3) END IF ! significant bend is to the right, or the the left? v2Mv1(1:3) = uvec2(1:3) - uvec1(1:3) v3Mv2(1:3) = uvec3(1:3) - uvec2(1:3) CALL DCross(v2Mv1, v3Mv2, tvec) triangle_area_m2 = DMagnitude(tvec) * half_R2 !on the surface of the planet (projected upward). !Use the following lines to output the outlines of the gore triangles: !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 !WRITE (72, "(1X,SP,ES12.5,',',ES12.5)") lon1, lat1 !WRITE (72, "('*** end of line segment ***')") !seismicity rate of this triangle, in EQs per second, above the CMT threshold magnitude: v_oblique_mmpa = (segment_slip_rate_mmpa(segment) + segment_slip_rate_mmpa(other_segment)) / 2.0D0 v_oblique_mps = v_oblique_mmpa * (0.001D0 / s_per_year) !Note that rate will be always positive (for now); use of factor SIGN is reserved for later. length_km = DMagnitude(v3Mv2) * radius_km ! for the deep base of the triangle IF (region_index == 2) THEN moment_rate_of_step_SI = coupledThickness_m * & ! coupled lithosphere thickness, in m & IT_assumed_mu_Pa(IT_column) * & ! elastic modulus mu, in Pa & v_oblique_mps * & ! velocity in fault plane, m/s & (0.5D0/DSIN(dip_radians)) * & ! csc(theta) & (length_km * 1000.D0) ! along-base length, in m seismicity_at_catalog_threshold = IT_catalog_assigned_event_rate(IT_column) * & & moment_rate_of_step_SI / IT_tGR_model_moment_rate(IT_column) ELSE 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 & (0.5D0/DSIN(dip_radians)) * & ! csc(theta) & (length_km * 1000.D0) ! along-base length, in m seismicity_at_catalog_threshold = T5_catalog_pure_event_rate(T5_column) * & & moment_rate_of_step_SI / T5_tGR_model_moment_rate(T5_column) END IF !correct seismicity to desired threshold moment/magnitude, !using the tapered Gutenberg-Richter distribution !with coefficients from Table 5 above: IF (region_index == 2) THEN tGR_beta = IT_beta(IT_column) corner_moment = IT_corner_moment(IT_column) catalog_threshold_moment = IT_catalog_threshold_moment(IT_column) ELSE tGR_beta = T5_beta(T5_column) corner_moment = T5_corner_moment(T5_column) catalog_threshold_moment = T5_catalog_threshold_moment(T5_column) END IF ! G = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-beta)) * & ! & EXP((catalog_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 = (catalog_threshold_moment - desired_threshold_moment_Nm) / corner_moment IF (argument < -100.0D0) THEN G_function = 0.0D0 ELSE G_function = ((desired_threshold_moment_Nm / catalog_threshold_moment)**(-tGR_beta)) * & & DEXP(argument) END IF EQs_per_s(threshold_index) = G_function * seismicity_at_catalog_threshold !distribute seismicity evenly over the triangle (on the Earth surface): EQs_per_s_per_m2(threshold_index) = EQs_per_s(threshold_index) / triangle_area_m2 END DO !prepare 3 pole_uvec's (poles of short arcs from uvec1->uvec2, uvec2->uvec3, uvec3-->uvec1: CALL DCross(uvec1, uvec2, tvec) CALL DMake_Uvec(tvec, pole_uvec1) CALL DCross(uvec2, uvec3, tvec) CALL DMake_Uvec(tvec, pole_uvec2) CALL DCross(uvec3, uvec1, tvec) CALL DMake_Uvec(tvec, pole_uvec3) !compare the (surface) area of this triangle 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 DLonLat_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 current segment): dot_test = uvec(1) * center_uvec(1) + uvec(2) * center_uvec(2) + uvec(3) * center_uvec(3) IF (dot_test > 0.9967D0) THEN ! within ~500 km; more precise distance will be computed !check whether uvec is on correct side of all 3 sides of (surface projection of) gore triangle: dot_test = uvec(1) * pole_uvec1(1) + uvec(2) * pole_uvec1(2) + uvec(3) * pole_uvec1(3) IF (dot_test > 0.0D0) THEN dot_test = uvec(1) * pole_uvec2(1) + uvec(2) * pole_uvec2(2) + uvec(3) * pole_uvec2(3) IF (dot_test > 0.0D0) THEN dot_test = uvec(1) * pole_uvec3(1) + uvec(2) * pole_uvec3(2) + uvec(3) * pole_uvec3(3) IF (dot_test > 0.0D0) THEN inside = .TRUE. 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_increment = SIGN * EQs_per_s_per_m2(threshold_index) !Prevent negative seismicity forecasts, which otherwise could happen where !negative (folded, SIGN == -1.0D0) gores extend slightly beyond the two !rectangular regions that they are meant to moderate! seismicity_increment = MAX(seismicity_increment, -0.9D0 * seismicity(i, j, threshold_index)) seismicity(i, j, threshold_index) = seismicity(i, j, threshold_index) + seismicity_increment * weight IF (threshold_index == 1) THEN fault_sum = fault_sum + fault_fraction * seismicity_increment * weight * grid_cell_areas(i, j) continuum_sum = continuum_sum + continuum_fraction * seismicity_increment * weight * grid_cell_areas(i, j) !N.B. The line above has no effect, because in this block of code (non-S) continuum_fraction == 0.0D0. ! The line is just kept here for parallel structure with the S(ubduction-zone) code. END IF 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 ! there is a significant bend (azimuth change) between these 2 adjacent segments EXIT check_others2 ! because there cannot be more than one match END IF ! head of other segment matches tail of current segment END IF ! this "other" segment has non-zero length END IF ! other_segment belongs to same fault END IF ! other_segment /= segment END DO check_others2 ! other_segment = 1, segments !========= END OF: Create(?) and apply triangular gore of seismicity (before "next" linked segment) =========== 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 ')") !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 ( fault_azimuth_radians ) END IF ! faults_used_in_orogen(orogen) END DO ! orogen = 1, orogens END IF ! orogens > 0 IF (subdivide.AND.(highlight == 0)) seismicity_bank = seismicity ! save whole array for future reference IF (highlight == 0) THEN ! full reporting on first (perhaps only?) pass through this loop !Note that Fortran unit 1 was already OPENed, at the top of this code. !Also, the header lines for this file were already WRITtEn. 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.5D0 * d_lat highest_latitude = center_latitude + 0.5D0 * d_lat dy_meters = Earth_radius_m * (highest_latitude - lowest_latitude) / 57.29578D0 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.5D0 * d_lon highest_longitude = center_longitude + 0.5D0 * d_lon dx_meters = Earth_radius_m * ((highest_longitude - lowest_longitude) / 57.29578D0) * & & DCOS (center_latitude / 57.29578D0) 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.D0 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.D0 * 3.15576D7 ! 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 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 DO threshold_index = 1, top_threshold_index dI = seismicity(i, j, threshold_index) * grid_cell_areas(i, j) integral(threshold_index) = integral(threshold_index) + dI END DO END IF END DO END DO integral_converted = integral * 100.0D0 * 365.25D0 * 24 * 60 * 60 ! whole array, if format_selector == 2 sum_total = fault_sum + continuum_sum fault_fraction = fault_sum / sum_total continuum_fraction = continuum_sum / sum_total fault_percent = 100.D0 * fault_fraction continuum_percent = 100.D0 * continuum_fraction 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 WRITE (*, "(/' Breakdown of seismicity:')") WRITE (*, "(/' ',ES11.4,F11.1,'% = fault fraction')") fault_sum, fault_percent WRITE (*, "(/' ',ES11.4,F11.1,'% = continuum fraction')") continuum_sum, continuum_percent WRITE (1, "('Breakdown of seismicity:')") WRITE (1, "(ES11.4,F11.1,'% = fault fraction')") fault_sum, fault_percent WRITE (1, "(ES11.4,F11.1,'% = continuum fraction')") continuum_sum, continuum_percent ELSE IF (format_selector == 2) THEN WRITE (*, "(/' Breakdown of seismicity, at lowest threshold:')") WRITE (*, "(/' ',ES11.4,F11.1,'% = fault fraction')") fault_sum, fault_percent WRITE (*, "(/' ',ES11.4,F11.1,'% = continuum fraction')") continuum_sum, continuum_percent 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 ELSE ! highlight > 0 (which implies both subdivide = .TRUE. and format_specifier == 1); compute and output column_0n_fraction.GRD: DO i = 1, rows DO j = 1, columns IF (seismicity_bank(i, j, 1) > 0.0) THEN seismicity(i, j, 1) = seismicity(i, j, 1) / seismicity_bank(i, j, 1) ELSE seismicity(i, j, 1) = 0. ! should not happen; just for insurance. END IF END DO END DO !whole array now reduced to fractions-of-total; should range from [0., 1.] !add values to sum_of_fractions matrix (which was zeroed at creation): DO i = 1, rows DO j = 1, columns sum_of_fractions(i, j) = sum_of_fractions(i, j) + seismicity(i, j, 1) END DO END DO !write column_0n_fraction.grd: WRITE (highlight_c2, "(I2)") highlight IF (highlight_c2(1:1) == ' ') highlight_c2(1:1) = '0' fraction_GRD_filename = "column_" // highlight_c2 // "_fraction.grd" OPEN (UNIT = 1, FILE = TRIM(fraction_GRD_filename)) ! absolute OPEN; overwrites any previous file of same name WRITE (1, "(3F11.5,' = lon_min, d_lon, lon_max')") lon_min, d_lon, lon_max WRITE (1, "(3F11.5,' = lat_min, d_lat, lat_max')") lat_min, d_lat, lat_max WRITE (1, "(10F7.4)") ((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. CLOSE (1) ! output .grd file END IF ! highlight == 0, or greater? IF (subdivide) THEN WRITE (*, *) WRITE (*, "(' -------------------------------------------------------------------')") WRITE (*, "(' Completed highlight pass ',I2,' out of ',I2,'.')") highlight, last_highlight WRITE (*, "(' -------------------------------------------------------------------')") WRITE (*, *) END IF !===================================================================================================================== END DO ! highlight = 0, last_highlight ! potentially repeating the entire calculation 13 times over! !===================================================================================================================== DEALLOCATE ( subsample_uvec ) IF (subdivide) THEN ! write out internal-control array sum_of_fractions: OPEN (UNIT = 1, FILE = "sum_of_fractions.grd") ! absolute OPEN; overwrites any previous file of same name WRITE (1, "(3F11.5,' = lon_min, d_lon, lon_max')") lon_min, d_lon, lon_max WRITE (1, "(3F11.5,' = lat_min, d_lat, lat_max')") lat_min, d_lat, lat_max WRITE (1, "(10F7.3)") ((sum_of_fractions(i, j), 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. CLOSE (1) ! output .grd file WRITE (*, *) WRITE (*, "(' CAUTION: The sum_of_fractions.grd file, if plotted, should show value 1.00')") WRITE (*, "(' in the central parts of NeoKinema model(s).')") WRITE (*, "(' However, it will typically show sum_of_fractions > 1.00 in regions')") WRITE (*, "(' outside the NeoKinema orogen(s), and also in marginal areas of the')") WRITE (*, "(' NeoKinema orogen(s) that are affected by diffusive smoothing.')") WRITE (*, "(' This happens because source code in the NONLINEAR STEP paragraph above')") WRITE (*, "(' makes it impossible to subdivide seismicity by column in these regions.')") WRITE (*, "(' This known issue only (slightly) complicates the use of the debugging maps')") WRITE (*, "(' provided under the subdivide = .TRUE. option; it does not indicate')") WRITE (*, "(' any problem with the total-seismicity predictions of LTS.')") END IF WRITE (*, *) WRITE (*, "(/' All computations completed. See output file ',A)") TRIM(output_grd_pathfile) CALL Pause() CONTAINS ! ************************************************************************ REAL*8 FUNCTION Area_of_cell(d_lon, d_lat, center_lat, R_m) !Returns area (in m**2) of a grid cell with width d_lon (degrees) !and height of d_lat (degrees), with central-point latitude of center_lat (degrees), !on a planet of radius R_m (meters). !Uses the exact (calculus) solution, not the approximate ~d_lat*d_lan*COS(center_lat)... USE DSphere ! provided by Peter Bird in file DSphere.f90 IMPLICIT NONE REAL*8, INTENT(IN) :: d_lon, d_lat, center_lat, R_m REAL*8 :: theta1, theta2 theta1 = (90.0D0 - (center_lat + d_lat/2.0D0)) * radians_per_degree ! integration limits, measured theta2 = (90.0D0 - (center_lat - d_lat/2.0D0)) * radians_per_degree ! from N pole, in radians theta1 = MAX(0.0D0, MIN(Pi, theta1)) ! just for good luck theta2 = MAX(0.0D0, MIN(Pi, theta2)) Area_of_cell = Two_Pi * (R_m**2) * (d_lon / 360.0D0) * (COS(theta1) - COS(theta2)) ! sic; !last term is equal to [ -COS(theta2) - (-COS(theta1) ]. END FUNCTION Area_of_cell REAL FUNCTION DATAN2F(y, x) IMPLICIT NONE REAL*8, INTENT(IN) :: y, x IF ((y /= 0.0D0).OR.(x /= 0.0D0)) THEN DATAN2F = DATAN2(y, x) ELSE DATAN2F = 0.0D0 END IF END FUNCTION DATAN2F SUBROUTINE Del_Gjxy_del_thetaphi (l_, uvec1, uvec2, uvec3, r_, dG) INTEGER, INTENT(IN) :: l_ ! element number or code REAL*8, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & & r_ ! position vector DOUBLE PRECISION, 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_. INTEGER, SAVE :: 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? DOUBLE PRECISION, DIMENSION(3,2) :: del_r_ ! theta- and phi-derivitives of r_ (in 3-D) DOUBLE PRECISION, DIMENSION(3,2) :: local ! local Theta, Phi unit vectors at r_ (xyz, SE) DOUBLE PRECISION, DIMENSION(3,2,2) :: del_local ! theta-, phi- derivitives of local DOUBLE PRECISION, DIMENSION(3,3), SAVE :: corner ! positions vector of corner nodes (xyz, 123) DOUBLE PRECISION, DIMENSION(3,3,2), SAVE :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) DOUBLE PRECISION, DIMENSION(3) :: tr_, tv, tvi, tvo, tv1, tv2, tv3, vfa, vfb ! temporary vector factors DOUBLE PRECISION :: cos_phi, cos_theta, factor, 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_ tvi(1:3) = uvec1(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 1) = tvo(1:3) tvi(1:3) = uvec2(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 2) = tvo(1:3) tvi(1:3) = uvec3(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 3) = tvo(1:3) DO j = 1, 3 tvi(1:3) = corner(1:3, j) CALL DLocal_Theta(tvi, tvo) post(1:3, j, 1) = tvo(1:3) CALL DLocal_Phi (tvi, tvo) post(1:3, j, 2) = tvo(1:3) END DO END IF ! begin calculations which depend on r_ tvi(1:3) = r_(1:3) ! upgrade to REAL8 CALL DMake_Uvec(tvi, tr_) CALL DLocal_Theta(tr_, tv) local(1:3,1) = tv(1:3) CALL DLocal_Phi(tr_, tv) local(1:3,2) = tv(1:3) ! Note: these functions will catch polar points; don't test again phi = DATAN2(tr_(2), tr_(1)) cos_phi = DCOS(phi) sin_phi = DSIN(phi) cos_theta = tr_(3) sin_theta = DSQRT(tr_(1)**2 + tr_(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) = - tr_(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.0D0, 0.0D0, 0.0D0 /) ! d.Phi/d.theta = 0 del_local(1:3,2,2) = (/ -cos_phi, -sin_phi, 0.0D0 /) ! 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(1:3) = corner(1:3, i1) tv2(1:3) = corner(1:3, i2) tv3(1:3) = corner(1:3, i3) CALL DCross(tv2, tv3, vfa) factor = 1.0D0 / DDot (tv1, vfa) vfb(1:3) = vfa(1:3) * factor DO x = 1, 2 ! unit velocity at node is S or E DO y = 1, 2 ! S- or E- component of nodal function tv1(1:3) = post(1:3, j, x) tvi(1:3) = local(1:3, y) DO m = 1, 2 ! theta- or phi-derivitive tv(1:3) = del_r_(1:3, m) tvo(1:3) = del_local(1:3, y, m) dG(j, x, y, m) = & & (DDot(tv,vfb)*DDot(tv1,tvi)) + & & (DDot(tr_,vfb)*DDot(tv1,tvo)) END DO END DO END DO END DO END SUBROUTINE Del_Gjxy_del_thetaphi SUBROUTINE DDumb_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*8, DIMENSION(3), INTENT(IN) :: vector ! uvec to point INTEGER, DIMENSION(:,:), INTENT(IN) :: node ! element definitions REAL*8, DIMENSION(:,:), INTENT(IN) :: xyz_nod ! uvecs of nodes REAL*8, DIMENSION(:,:), INTENT(IN) :: center ! uvec of each element (uvec) REAL*8, DIMENSION(:), INTENT(IN) :: a_ ! element areas (plane; R == 1.0) REAL*8, INTENT(OUT) :: s1, s2, s3 ! results !- - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: i1, i2, i3 REAL*8, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL*8 :: d1, dc, t IF (element == 0) THEN WRITE (*,"(' ERROR: element = 0 in DDumb_s123')") CALL DTraceback 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.D0) THEN WRITE (*,"(' ERROR: Internal vector >= 90 deg. from element in DDumb_s123')") CALL DTraceback 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 DCross(tvi, tvo, tv) s1 = 0.5D0 * 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 DCross(tvi, tvo, tv) s2 = 0.5D0 * DOT_PRODUCT(tv1, tv) / a_(element) s3 = 1.00D0 - s1 - s2 END SUBROUTINE DDumb_s123 SUBROUTINE E_rate(R, l_, nodes, G, dG, theta_, vw, eps_dot) ! Evaluate strain-rate at one point in one spherical continuum element (# l_); ! the specific point is implied by the values of arrays G and dG supplied, ! but note that the value of theta_ must also be consistent. REAL*8, INTENT(IN) :: R ! radius of planet, in m INTEGER, INTENT(IN) :: l_ ! element number INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes DOUBLE PRECISION, DIMENSION(3,2,2), INTENT(IN) :: G ! nodal functions @ selected point DOUBLE PRECISION, DIMENSION(3,2,2,2), INTENT(IN):: dG ! derivitives of nodal functions @ selected point REAL*8, INTENT(IN) :: theta_ ! colatitude, radians DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: vw REAL*8, DIMENSION(3), INTENT(OUT) :: eps_dot DOUBLE PRECISION, DIMENSION(3) :: sums INTEGER :: iv, iw, j DOUBLE PRECISION :: cott, csct, prefix cott = 1.0D0 / DTAN(1.0D0 * theta_) csct = 1.0D0 / DSIN(1.0D0 * theta_) prefix = 1.0D0 / R sums(1:3) = 0.0D0 DO j = 1, 3 iv = 2 * nodes(j, l_) - 1 iw = iv + 1 ! epsilon_dot_sub_theta_theta sums(1) = sums(1) + & & vw(iv) * prefix * dG(j,1,1,1) + & & vw(iw) * prefix * dG(j,2,1,1) ! epsilon_dot_sub_theta_phi sums(2) = sums(2) + & & vw(iv) * prefix * 0.5D0 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) + & & vw(iw) * prefix * 0.5D0 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) ! epsilon_dot_sub_phi_phi sums(3) = sums(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 eps_dot(1:3) = sums(1:3) END SUBROUTINE E_rate SUBROUTINE Gjxy (l_, uvec1, uvec2, uvec3, r_, G) INTEGER, INTENT(IN) :: l_ ! element number or code REAL*8, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & ! corners & r_ ! position vector DOUBLE PRECISION, 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_. INTEGER, SAVE :: 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 DOUBLE PRECISION, DIMENSION(3,2) :: local ! local unit vectors at r_ (xyz, SE) DOUBLE PRECISION, DIMENSION(3,3), SAVE :: corner ! positions vector of corner nodes (xyz, 123) DOUBLE PRECISION, DIMENSION(3,3,2), SAVE :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) DOUBLE PRECISION, DIMENSION(3) :: tr_, tvi, tvo, tv1, tv2, tv3, vf ! temporary vector factors DOUBLE PRECISION :: 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_ tvi(1:3) = uvec1(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 1) = tvo(1:3) tvi(1:3) = uvec2(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 2) = tvo(1:3) tvi(1:3) = uvec3(1:3) CALL DMake_Uvec(tvi, tvo) corner(1:3, 3) = tvo(1:3) DO j = 1, 3 tvi(1:3) = corner(1:3, j) CALL DLocal_Theta(tvi, tvo) post(1:3, j, 1) = tvo(1:3) CALL DLocal_Phi (tvi, tvo) post(1:3, j, 2) = tvo(1:3) END DO END IF ! begin computations which depend on r_ tvi(1:3) = r_(1:3) CALL DMake_UVec(tvi, tr_) CALL DLocal_Theta(tr_, tvo) local(1:3,1) = tvo(1:3) CALL DLocal_Phi(tr_, tvo) local(1:3,2) = tvo(1:3) 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 DCross(tv2, tv3, vf) f_sup_j = DDot(tr_, vf) / DDot (tv1, vf) DO x = 1, 2 tv1(1:3) = post(1:3, j, x) DO y = 1, 2 tv2(1:3) = local(1:3, y) G(j, x, y) = f_sup_j * DDot(tv1, tv2) END DO END DO END DO END SUBROUTINE Gjxy INTEGER FUNCTION DIBelow (x) IMPLICIT NONE REAL*8, INTENT(IN) :: x DIBelow = INT(x) IF (x < (1.0D0 * DIBelow)) DIBelow = DIBelow - 1 RETURN END FUNCTION DIBelow SUBROUTINE DLearn_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*8, DIMENSION(:,:), INTENT(IN) :: node_uvec ! uvecs of nodes LOGICAL, INTENT(IN) :: chatty REAL*8, DIMENSION(:), INTENT(OUT) :: a_ REAL*8, 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*8, 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 DCross (a, b, c) a_(l_) = 0.5D0 * DMagnitude(c) !second, compute center t(1:3) = (node_uvec(1:3,i1)+node_uvec(1:3,i2)+node_uvec(1:3,i3))/3.0D0 CALL DMake_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 (SaIT_Node_Uvec(j1, ib) .AND. SaIT_Node_Uvec(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (SaIT_Node_Uvec(j2, ib) .AND. SaIT_Node_Uvec(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (SaIT_Node_Uvec(j3, ib) .AND. SaIT_Node_Uvec(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 (SaIT_Node_Uvec(j1, ib) .AND. SaIT_Node_Uvec(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (SaIT_Node_Uvec(j2, ib) .AND. SaIT_Node_Uvec(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (SaIT_Node_Uvec(j3, ib) .AND. SaIT_Node_Uvec(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 DLearn_Spherical_Triangles REAL FUNCTION LookUp_ETOPO5 (longitude, latitude) IMPLICIT NONE REAL*8, INTENT(IN) :: longitude, latitude REAL*8 :: 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.0D0) lon = lon - 360.0D0 IF (lon < 0.0D0) lon = lon + 360.0D0 lat = MIN(90.0D0, MAX(-90.0D0, latitude)) rRow = 12.0D0 * lat rColumn = 12.0D0 * lon iRow1 = DIBelow(rRow) iRow2 = iRow1 + 1 jColumn1 = DIBelow(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*8, INTENT(IN) :: longitude, latitude REAL*8 :: 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.0D0) lon = lon - 360.0D0 IF (lon < 0.0D0) lon = lon + 360.0D0 lat = MIN(90.0D0, MAX(-90.0D0, latitude)) rRow = 10.0D0 * lat rColumn = 10.0D0 * lon iRow1 = DIBelow(rRow) iRow2 = iRow1 + 1 jColumn1 = DIBelow(rColumn) jColumn2 = jColumn1 + 1 fracRow = rRow - iRow1 fracColumn = rColumn - jColumn1 IF (iRow1 < -720) THEN LookUp_age_1p5 = 220.0D0 ! 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 DMoment(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*8, INTENT(IN) :: magnitude DMoment = 10.D0**((1.5D0 * magnitude) + 9.05D0) END FUNCTION DMoment 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 DPrincipal_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*8, INTENT(IN) :: e11, e12, e22 REAL*8, INTENT(OUT) :: e1, e2, u1x, u1y, u2x, u2y REAL*8 :: c, f1, f11, f12, f2, f22, r, scale, smallest, test, theta ! Smallest number that can be squared without underflowing: smallest = 1.1D0 * DSQRT(TINY(f1)) ! First, check for trivial isotropic case: IF ((e11 == e22).AND.(e12 == 0.D0)) THEN ! In this case, directions are arbitrary: e1 = e11 u1x = 1.0D0 u1y = 0.0D0 e2 = e22 u2x = 0.0D0 u2y = 1.0D0 RETURN END IF ! Else, re-scale matrix to prevent under/overflows: scale = MAX(DABS(e11), DABS(e12), DABS(e22)) f11 = e11 / scale IF (DABS(f11) < smallest) f11 = 0.0D0 f12 = e12 / scale IF (DABS(f12) < smallest) f12 = 0.0D0 f22 = e22 / scale IF (DABS(f22) < smallest) f22 = 0.0D0 ! Find eigenvalues and eigenvectors of scaled matrix: r = DSQRT(((f11 - f22)*0.5D0)**2 + f12**2) c = (f11 + f22)*0.5D0 f1 = c - r f2 = c + r test = 0.01D0 * MAX (DABS(f1), DABS(f2)) IF ((DABS(f12) > test).OR.(DABS(f11 - f1) > test)) THEN theta = DAtan2F((f11 - f1), -f12) ELSE theta = DAtan2F(f12, (f1 - f22)) END IF u1x = DCOS(theta) u1y = DSIN(theta) u2x = u1y u2y = -u1x ! Undo the scaling e1 = scale * f1 e2 = scale * f2 END SUBROUTINE DPrincipal_Axes_22 SUBROUTINE DPrompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE DPrompt_for_Logical SUBROUTINE DPull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL*8 :: factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.D0) THEN highest = MAXVAL(s) medium = 1.00D0 - lowest - highest IF (medium > 0.D0) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0.D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00D0 / (1.00D0 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00D0 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0.D0 s(sideb) = 0.D0 END IF END IF END SUBROUTINE DPull_in LOGICAL FUNCTION SaIT_Node_Uvec(i,j) ! Are node_uvec #i and #j the same vector? INTEGER, INTENT(IN) :: i, j !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 SaIT_Node_Uvec = .TRUE. ELSE SaIT_Node_Uvec = .FALSE. END IF ELSE SaIT_Node_Uvec = .FALSE. END IF ELSE SaIT_Node_Uvec = .FALSE. END IF END FUNCTION SaIT_Node_Uvec LOGICAL FUNCTION DSaIT_Uvec(uvec1, uvec2) REAL*8, DIMENSION(3), INTENT(IN) :: uvec1, uvec2 IF (uvec1(1) == uvec2(1)) THEN IF (uvec1(2) == uvec2(2)) THEN IF (uvec1(3) == uvec2(3)) THEN DSaIT_Uvec = .TRUE. ELSE DSaIT_Uvec = .FALSE. END IF ELSE DSaIT_Uvec = .FALSE. END IF ELSE DSaIT_Uvec = .FALSE. END IF END FUNCTION DSaIT_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 DWhich_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*8, 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*8, DIMENSION(:,:),INTENT(IN) :: xyz_node ! uvecs of nodes REAL*8, DIMENSION(:,:),INTENT(IN) :: center ! center uvecs of elements REAL*8, 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*8, INTENT(INOUT) :: s1, s2, s3 ! OUTPUT !- - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: back1, back2, back3, i, iet, l_ REAL*8 :: r2, r2min, s1t, s2t, s3t REAL*8, 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.0D0; s2 = 0.0D0; s3 = 0.0D0 ! default !find closest element center to initialize search r2min = 4.01D0 ! 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.540D0) RETURN ELSE ! warm (re)start; re-use last element# iet = iele END IF ! cold_start ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 back3 = -3 is_it_here: DO ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! CALL DDumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL DPull_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 DDumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL DPull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL DDumb_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.D0) 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.D0) 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.D0) 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 DWhich_Spherical_Triangle SUBROUTINE Prompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE Prompt_for_Logical END PROGRAM Long_Term_Seismicity_v11