PROGRAM EQ_classification_II ! ! Reads: (1) PB2002_poles.dat (derived from PB2002_poles.xls by dropping headers, names, and areas) ! and: (2) PB2002_steps.dig (see Table 2 from Bird [2002]: step characteristics from model PB2002) ! and: (3) PB2002_orogens.dig (CCW outlines of orogens; events and steps here are *-ed) ! and: (4) XXXX.eqc (catalog of shallow earthquakes: e.g., CMT, or Pacheco and Sykes [1992]) ! ! Creates 2 * 8 = 16 .eqc files: (Note: "Pure" .eqc's don't include *-ed EQs or steps.) ! class 0: XXXX_PB2002_INT_pure.eqc and XXXX_PB2002_INT_all.eqc INTraplate (plate-INTerior) ! class 1: XXXX_PB2002_CCB_pure.eqc and XXXX_PB2002_CCB_all.eqc Continental Convergent Boundary ! class 2: XXXX_PB2002_CTF_pure.eqc and XXXX_PB2002_CTF_all.eqc Continental Transform Fault ! class 3: XXXX_PB2002_CRB_pure.eqc and XXXX_PB2002_CRB_all.eqc Continental Rift Boundary ! class 4: XXXX_PB2002_OSR_pure.eqc and XXXX_PB2002_OSR_all.eqc Oceanic Spreading Ridge ! class 5: XXXX_PB2002_OTF_pure.eqc and XXXX_PB2002_OTF_all.eqc Oceanic Transform Fault ! class 6: XXXX_PB2002_OCB_pure.eqc and XXXX_PB2002_OCB_all.eqc Oceanic Convergent Boundary ! class 7: XXXX_PB2002_SUB_pure.eqc and XXXX_PB2002_SUB_all.eqc SUBduction zone ! ! by Peter Bird, UCLA, 2003.02.14-2003.04.07 ! This is the second attempt at a classification scheme for shallow earthquakes. ! This method is based on the assessment of the relative probabilities that ! a given actual earthquake was caused by each of the plate boundary steps in PB2002. ! Factors considered include: ! A: relative seismicity of each step, based on its length, velocity, and class. ! B: spatial PDF around each boundary step based on known apparent half-widths, ! and earthquake and step location errors around its ends; ! And for each of a set of (up to) 5 model earthquakes for each boundary class: ! C: relative probability of this particular model earthquake; ! D: depth PDF for this model earthquake (if depth is known); ! E: rotation-angle PDF to test match between actual and model moment tensor; ! for Pacheco & Sykes [1992] catalog this factor must be simplified or dropped. ! ! Addition by Peter Bird, UCLA, in 2003.06: ! A set of 7 files "PB2002_XXX_pure.dat" (where XXX is CCB, CTF, CRB, OSR, OTF, OCB, SUB) ! are also created to record earthquake numbers and total moments associated with each ! pure plate boundary step of type XXX. ! A threshhold magnitude "class_threshhold magnitude" is imposed to make the sums consistent worldwide. ! After running this, the file names should probably be modified to show the ! origin of the catalog that was used (e.g., add "_CMT" to the file names, just before ".dat"). ! (It is understood that when creating these files, parameter Monte_Carlo_method should be .FALSE.) ! ! **** Read the REAL, PARAMETER section below and adjust the values! ***************** USE Sphere ! Peter Bird's personal library of spherical-trigonometry and Cartesian-vector routines USE Numerical_Libraries ! IMSL Library, used for ANORDF = distribution function of a standard normal (Gaussian) variable; ! and also for RNOPT/RNSET/RNUN if Monte_Carlo_method = .TRUE. USE Quaternion ! for DCROT function, from Kagan [1991], Geophys. J. Int., v. 106, p. 709-716. IMPLICIT NONE INTEGER, PARAMETER :: nPlates = 52 ! referring to PB2002 model CHARACTER*1 :: eq_tenths, pure_epicenter_c1, pure_step_c1 CHARACTER*2 :: eq_month, eq_day, eq_hour, eq_minute, eq_second CHARACTER*2, DIMENSION(nPlates) :: ID CHARACTER*3 :: c3 CHARACTER*5 :: zone ! time zone CHARACTER*8 :: date ! date program is run CHARACTER*9 :: catalog_name, prefix CHARACTER*10 :: clock_time ! wall clock time CHARACTER*27 :: c27 CHARACTER*80 :: appended_data, catalog_file_name, out_file_name CHARACTER*3, DIMENSION(0:7) :: class_c3 ! = INT, plus boundary classes 1-7 CHARACTER*1, DIMENSION(:), ALLOCATABLE :: step_continuity_c1, class_continuity_c1, star_c1 CHARACTER*3, DIMENSION(:), ALLOCATABLE :: class CHARACTER*5, DIMENSION(:), ALLOCATABLE :: c5 INTEGER, PARAMETER :: nOrogens = 13, mostInOneOrogen = 500 ! referring to PB2002_orogens.dig; really, most is 315 INTEGER :: i ! = indentifier for the PB2002 plate boundary step being considered (up to 5819 = nSteps) INTEGER :: j ! = indentifier for the model EQ being considered (1 to 4) INTEGER :: k ! = indentifier of the class of plate boundary step i; range 1..7 = CCB, CTF, CRB, OSR, OTF, OCB, SUB INTEGER :: best_i, best_k, catalog_ID, crossstrike_dip_degrees, & & eq_year, eq_depth_int, e1_plunge, e1_azimuth, & & e2_plunge, e2_azimuth, e3_plunge, e3_azimuth, & & i0, iEvent, iLeftPlate, iRightPlate, ios, iseed, & & n, nEvents, nSteps, & & slab_velocity_azimuth_integer, slab_downdip_azimuth_integer, & & P_azimuth_degrees, P_plunge_degrees, T_azimuth_degrees, T_plunge_degrees, & & velocity_dip_degrees INTEGER*2, DIMENSION(4) :: actual_TpaPpa_degrees, model_TpaPpa_degrees ! TYPE *must* agree with MODULE Quaternion! INTEGER, DIMENSION(7) :: best_i_in_class INTEGER, DIMENSION(7) :: N_over_mt ! corresponds to N_k in paper INTEGER, DIMENSION(7) :: relative_class_probability_int ! relative % likelihood of 7 classes, output for each EQ except INT's INTEGER, DIMENSION(0:7) :: count_all, count_pure, percentages ! note: 0 subscript corresponds to INT subcatalogs INTEGER, DIMENSION(8) :: datetimenumber ! for output from DATE_AND_TIME INTEGER, DIMENSION(nOrogens) :: nInEachOrogen INTEGER, DIMENSION(:), ALLOCATABLE :: azimuth_integer, age_Ma, & & elevation, eq_count_over_threshhold, & & k_of, & & sequence_number, velocity_azimuth_integer LOGICAL :: dip_is_to_right, eq_depth_known, inside, mechanism_known, normal, & & proximal, pure_epicenter, pure_step, relative_velocity_is_divergent, & & strikeslip, thrust, unknown LOGICAL, DIMENSION(7) :: in_competition REAL, PARAMETER :: radius_km = 6371.0 ! of the spherical model of the Earth !********************** USER SELECTIONS: REVIEW THESE CAREFULLY! ******************************************** !All boundaries: LOGICAL, PARAMETER :: Monte_Carlo_method = .FALSE. ! (alternative is assignment to maximum-probability step) REAL, PARAMETER :: alongStep_sigma_km = 32.5 ! so combined EQ mislocation and step-mislocation unlikely to exceed twice this value REAL, PARAMETER :: earthquake_depth_sigma_km = 15.0 ! always *some* probability down to 70 km when class123456_cutoff_depth_km = 25.0 REAL, PARAMETER :: Phi_max_degrees = 60.0 ! cut-off rotation for matching model and actual moment tensors REAL, PARAMETER :: class123456_cutoff_depth_km = 25.0 INTEGER, PARAMETER :: CCB_pure_totalEvents = 274 ! N_k values for m_t = 5.66, Harvard CMT 1977.01.01-2002.09.30 INTEGER, PARAMETER :: CTF_pure_totalEvents = 198 ! These values last adjusted 2003.04.08, 12:45, based on INTEGER, PARAMETER :: CRB_pure_totalEvents = 134 ! fifth iteration of boot-strap method! INTEGER, PARAMETER :: OSR_pure_totalEvents = 142 ! This iteration was run with Monte_Carlo_method = .FALSE. INTEGER, PARAMETER :: OTF_pure_totalEvents = 838 ! These are suggested for processing of the full, incomplete CMT catalog. INTEGER, PARAMETER :: OCB_pure_totalEvents = 105 ! INTEGER, PARAMETER :: SUB_pure_totalEvents = 2049 ! !INTEGER, PARAMETER :: CCB_pure_totalEvents = 10 ! N_k values for m_t = 6.99, Harvard CMT 1977.01.01-2002.09.30 !INTEGER, PARAMETER :: CTF_pure_totalEvents = 9 ! and Ekstrom and Nettles [1997] 1976 events. !INTEGER, PARAMETER :: CRB_pure_totalEvents = 4 ! These values last adjusted 2003.04.08, 14:18, based on !INTEGER, PARAMETER :: OSR_pure_totalEvents = 1 ! fifth iteration of boot-strap method !INTEGER, PARAMETER :: OTF_pure_totalEvents = 6 ! This iteration was run with Monte_Carlo_method = .FALSE. !INTEGER, PARAMETER :: OCB_pure_totalEvents = 10 ! These are suggested for max-prob processing of Pacheco & Sykes, 1900-1975. !INTEGER, PARAMETER :: SUB_pure_totalEvents = 111 ! ! V !INTEGER, PARAMETER :: CCB_pure_totalEvents = 237 ! 249 247 245 229 237 ! Record of 5 consecutive (post-5-for-convergence) iterations !INTEGER, PARAMETER :: CTF_pure_totalEvents = 202 ! 209 195 204 205 202 ! of the Monte Carlo method applied to CMT_shallow_complete.eqc !INTEGER, PARAMETER :: CRB_pure_totalEvents = 150 ! 140 138 136 144 150 ! (m_t = 5.66; complete). Computed 2003.04.08, 15:35-15:40. !INTEGER, PARAMETER :: OSR_pure_totalEvents = 153 ! 157 158 153 150 153 ! These can be recycled in classifications of the full !INTEGER, PARAMETER :: OTF_pure_totalEvents = 809 ! 819 811 827 816 809 ! CMT_shallow.eqc, to capture the effects of the !INTEGER, PARAMETER :: OCB_pure_totalEvents = 141 ! 128 124 122 132 141 ! uncertainties in the N_k under the Monte Carlo regime. !INTEGER, PARAMETER :: SUB_pure_totalEvents = 2046 ! 2037 2065 2051 2061 2046 ! ! V !INTEGER, PARAMETER :: CCB_pure_totalEvents = 10 ! 10 12 10 10 10 ! Record of 5 consecutive (post-5-for-convergence) iterations !INTEGER, PARAMETER :: CTF_pure_totalEvents = 6 ! 11 11 6 10 6 ! of the Monte Carlo method applied to CMT_shallow_m7_1976-2002.eqc !INTEGER, PARAMETER :: CRB_pure_totalEvents = 5 ! 4 3 5 4 5 ! (m_t = 6.99; Ekstrom & Nettles [1997] for 1976, + CMT)). !INTEGER, PARAMETER :: OSR_pure_totalEvents = 1 ! 1 1 1 1 1 ! These can be recycled in classifications of the Pacheco !INTEGER, PARAMETER :: OTF_pure_totalEvents = 6 ! 6 6 8 6 6 ! & Sykes catalog, to capture the effects of the !INTEGER, PARAMETER :: OCB_pure_totalEvents = 11 ! 20 8 11 10 11 ! uncertainties in the N_k under the Monte Carlo regime. !INTEGER, PARAMETER :: SUB_pure_totalEvents = 112 ! 109 110 110 110 112 ! REAL, PARAMETER :: CCB_pure_totalLength_km = 12516.0 REAL, PARAMETER :: CTF_pure_totalLength_km = 19375.0 REAL, PARAMETER :: CRB_pure_totalLength_km = 18126.0 REAL, PARAMETER :: OSR_pure_totalLength_km = 61807.0 ! these values last adjusted 2003.02.14, afternoon, REAL, PARAMETER :: OTF_pure_totalLength_km = 43902.0 ! based on PB2002_XXX_pure_lineIntegral.xls, same date REAL, PARAMETER :: OCB_pure_totalLength_km = 13236.0 REAL, PARAMETER :: SUB_pure_totalLength_km = 38990.0 REAL, PARAMETER :: CCB_pure_meanVelocity_mmpa = 18.2 REAL, PARAMETER :: CTF_pure_meanVelocity_mmpa = 21.5 REAL, PARAMETER :: CRB_pure_meanVelocity_mmpa = 19.0 REAL, PARAMETER :: OSR_pure_meanVelocity_mmpa = 48.3 ! these values last adjusted 2003.02.14, afternoon, REAL, PARAMETER :: OTF_pure_meanVelocity_mmpa = 40.4 ! based on PB2002_XXX_pure_lineIntegral.xls, same date REAL, PARAMETER :: OCB_pure_meanVelocity_mmpa = 19.2 REAL, PARAMETER :: SUB_pure_meanVelocity_mmpa = 61.5 !non-SUB boundaries (symmetrical): REAL, PARAMETER :: CCB_halfwidth_km = 189.0 REAL, PARAMETER :: CTF_halfwidth_km = 257.0 REAL, PARAMETER :: CRB_halfwidth_km = 154.0 REAL, PARAMETER :: OSR_halfwidth_km = 132.0 REAL, PARAMETER :: OTF_halfwidth_km = 128.0 REAL, PARAMETER :: OCB_halfwidth_km = 186.0 ! these values last adjusted 2003.02.12, 4:22 PM. !SUB boundaries (asymmetrical): REAL, PARAMETER :: SUB_landward_halfwidth_km = 220.0 ! landward width of subduction lanes, in km from trench REAL, PARAMETER :: SUB_seaward_halfwidth_km = 135.0 ! seaward half-width of subduction lanes, in km from trench REAL, PARAMETER :: SUB_peakSeismicity_km = 55.0 ! location of seismicity peak, relative to trench !The following lines create a standard quadratic model of depth to the slab top, !which is used (only) in depth PDF "D" in connection with model EQs on SUB steps: REAL, PARAMETER :: trench_depth_km = 8.25 ! p.221 of Le Pichon et al. [1976]; measured from sea level REAL, PARAMETER :: slope_at_trench = 0.1400 ! = 8 degrees; ibid (not a great source!) REAL, PARAMETER :: slab_depth_under_arc_km = 100.0 ! reference? Measured from sea level. REAL, PARAMETER :: arc_distance_km = 200.0 ! used in conjunction with above line to pin quadratic model of slab-top depth REAL, PARAMETER :: interplate_thrust_depth_sigma_fraction = 0.10 ! uncertainty in standard quadratic model of slab-top depth !Note: Use SUB_quadratic_profile.xls to experiment with these parameters. !Threshhold magnitudes by class: ! (Note: These values are used only in creating the 7 files PB2002_XXX_pure.dat, ! but when creating the subcatalog (.eqc) files, no threshhold is imposed. ! These following lines were added in the 2003.06 extension of this program.) REAL, PARAMETER :: CCB_threshhold_magnitude = 5.66 ! from Table 5 of Bird & Kagan, BSSA, 2003(?) REAL, PARAMETER :: CTF_threshhold_magnitude = 5.66 REAL, PARAMETER :: CRB_threshhold_magnitude = 5.33 REAL, PARAMETER :: OSR_threshhold_magnitude = 5.33 REAL, PARAMETER :: OTF_threshhold_magnitude = 5.50 REAL, PARAMETER :: OCB_threshhold_magnitude = 5.66 REAL, PARAMETER :: SUB_threshhold_magnitude = 5.66 ! ************************************************************************************************************ REAL :: A ! relative number (or rate) of EQs (over m_t) produced by the current plate boundary step REAL :: B ! PDF (per square km) that an EQ produced by this step is at the actual (longitude, latitude) REAL :: C ! relative probability of this model earthquake, compared to other (up to 3) models REAL :: D ! PDF (per km) that an EQ of this model occured at the actual EQ depth (if known; otherwise punt). REAL :: E ! PDF (per steradian) that this actual EQ mechanism (if known) matches the model EQ REAL :: sum_CDE ! SUM(j=1,4) {C * D * E} for all model earthquakes on this step REAL :: P_prime ! relative probability that this step produced the EQ in question REAL :: s_km ! local coordinate, measured along direction of plate boundary step, from 0 at initial point. REAL :: X ! cross-trace component of spatial PDF B; a function of cross-step offset (from peak of seismicity) REAL :: Y ! along-trace component of spatial PDF B; a function of along-step distance s REAL :: angle, angle_sum, argument, azimuth_radians, & & basis_length_km, & & closing_angle_radians, crossstrike_dip_radians, & & distance_km, distance_at_peak_km, ds, & & eq_Elon, eq_Nlat, eq_mag, eq_moment_Nm, & & fault_dip_degrees, fault_strike_degrees, fraction_toward_uvec1, from_trench_azimuth_radians, & & highest_P_prime, integral, & & interplate_thrust_depth_km, interplate_thrust_depth_sigma_km, & & lat, lane_width_km, lon, & & midpoint_compass, minimum_depth_km, & & offset_km, offside_radians, opening_angle_radians, & & partial_sum, & & quadratic, & & random, relative_velocity_azimuth_degrees, rotation_degpMa, & & s_radians, s1, s2, scalar_product, & & SUB_landward_radians, SUB_landward_sigma_km, SUB_seaward_radians, SUB_seaward_sigma_km, & & sum_P_prime_all_classes, & & to_uvec1_radians, to_uvec2_radians, t_latitude, t_longitude, t1, t2, & & velocity_dip_radians, velocity_East, velocity_North REAL*8 :: minimum_angle_degrees ! TYPE *must* agree with MODULE Quaternion! REAL, DIMENSION(3) :: c1_uvec, c2_uvec, c3_uvec, c4_uvec, & & eq_uvec, inward_uvec, omega_uvec, pole_uvec, tEuler, trench_uvec, tvec, & & uvec, uvec1, uvec2, velocity REAL, DIMENSION(4) :: dot_side ! sine of angle between EQ and sides of a subduction lane; positive inward REAL, DIMENSION(7) :: class_threshhold_magnitude REAL, DIMENSION(7) :: highest_P_prime_in_class REAL, DIMENSION(7) :: L_km ! corresponds to L_k in paper. REAL, DIMENSION(7) :: sum_P_prime_in_class ! sums of P_prime for all steps in one class REAL, DIMENSION(7) :: upper_fraction REAL, DIMENSION(7) :: V_mmpa ! corresponds to V_k in paper REAL, DIMENSION(6) :: h_km ! corresponds to h_k in paper; no entry for SUB because that is asymmetrical REAL, DIMENSION(3, nPlates) :: Euler REAL, DIMENSION(3, mostInOneOrogen, nOrogens) :: orogen_uvecs REAL, DIMENSION(:), ALLOCATABLE :: arc_km_from_step, & & dextral_mmpa, distance_km_from_step, divergence_mmpa, & & latitude, length_km, longitude, & & moment_sum_over_threshhold, & & old_longitude, old_latitude, & & subduction_azimuth_1, subduction_azimuth_2, & & velocity_mmpa, Y_factor_of_step REAL, DIMENSION(:,:), ALLOCATABLE :: begin_uvec, end_uvec REAL, DIMENSION(:,:,:), ALLOCATABLE :: lane_uvecs !======================================================================================================== WRITE (*, "(' Running EQ_classification_II.exe, from EQ_classification_II.f90....')") 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" N_over_mt(1) = CCB_pure_totalEvents N_over_mt(2) = CTF_pure_totalEvents N_over_mt(3) = CRB_pure_totalEvents N_over_mt(4) = OSR_pure_totalEvents N_over_mt(5) = OTF_pure_totalEvents N_over_mt(6) = OCB_pure_totalEvents N_over_mt(7) = SUB_pure_totalEvents L_km(1) = CCB_pure_totalLength_km L_km(2) = CTF_pure_totalLength_km L_km(3) = CRB_pure_totalLength_km L_km(4) = OSR_pure_totalLength_km L_km(5) = OTF_pure_totalLength_km L_km(6) = OCB_pure_totalLength_km L_km(7) = SUB_pure_totalLength_km V_mmpa(1) = CCB_pure_meanVelocity_mmpa V_mmpa(2) = CTF_pure_meanVelocity_mmpa V_mmpa(3) = CRB_pure_meanVelocity_mmpa V_mmpa(4) = OSR_pure_meanVelocity_mmpa V_mmpa(5) = OTF_pure_meanVelocity_mmpa V_mmpa(6) = OCB_pure_meanVelocity_mmpa V_mmpa(7) = SUB_pure_meanVelocity_mmpa h_km(1) = CCB_halfwidth_km h_km(2) = CTF_halfwidth_km h_km(3) = CRB_halfwidth_km h_km(4) = OSR_halfwidth_km h_km(5) = OTF_halfwidth_km h_km(6) = OCB_halfwidth_km ! no entry for SUB because it is asymmetrical SUB_landward_radians = SUB_landward_halfwidth_km / radius_km ! landward width of subduction lanes, in radians, from trench SUB_seaward_radians = SUB_seaward_halfwidth_km / radius_km ! seaward half-width of subduction lanes, in radians, from trench SUB_landward_sigma_km = 0.5 * (SUB_landward_halfwidth_km - SUB_peakSeismicity_km) ! sigmas for Gaussian approximation of B SUB_seaward_sigma_km = 0.5 * (SUB_seaward_halfwidth_km + SUB_peakSeismicity_km) class_threshhold_magnitude(1) = CCB_threshhold_magnitude class_threshhold_magnitude(2) = CTF_threshhold_magnitude class_threshhold_magnitude(3) = CRB_threshhold_magnitude class_threshhold_magnitude(4) = OSR_threshhold_magnitude class_threshhold_magnitude(5) = OTF_threshhold_magnitude class_threshhold_magnitude(6) = OCB_threshhold_magnitude class_threshhold_magnitude(7) = SUB_threshhold_magnitude !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Read 1: PB2002_poles.dat (derived from PB2002_poles.xls by dropping headers, names, and areas) OPEN (UNIT = 1, FILE = 'PB2002_poles.dat', STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: input file PB2002_poles.dat could not be opened.')") CALL Pause() STOP END IF WRITE(*, "(' Reading PB2002_poles.dat...')") DO i = 1, nPlates READ (1, "(A2,F9.3,F10.3,F9.4)", IOSTAT = ios) ID(i), lat, lon, rotation_degpMa IF (ios /= 0) THEN WRITE (*,"(' ERROR: could not read Euler pole for plate', I3)") i CALL Traceback() STOP END IF CALL LonLat_2_Uvec (lon, lat, uvec) uvec(1:3) = uvec(1:3) * rotation_degpMa Euler(1:3, i) = uvec(1:3) END DO ! i = 1, nPlates CLOSE (1) ! read (2) PB2002_steps.dat ----------------------------------------------------------------- ! first reading is just to count steps OPEN (UNIT = 2, FILE = "PB2002_steps.dat", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: PB2002_steps.dat not available to read in.')") CALL Pause() END IF nSteps = 0 count_steps: DO READ (2, *, IOSTAT = ios) i ! the sequence number; might not be in order(?) IF (ios /= 0) EXIT count_steps nSteps = nSteps + 1 END DO count_steps CLOSE (2) ALLOCATE ( sequence_number(nSteps) ) ALLOCATE ( step_continuity_c1(nSteps) ) ALLOCATE ( c5(nSteps) ) ALLOCATE ( old_longitude(nSteps) ) ALLOCATE ( old_latitude(nSteps) ) ALLOCATE ( longitude(nSteps) ) ALLOCATE ( latitude(nSteps) ) ALLOCATE ( length_km(nSteps) ) ALLOCATE ( azimuth_integer(nSteps) ) ALLOCATE ( velocity_mmpa(nSteps) ) ALLOCATE ( velocity_azimuth_integer(nSteps) ) ALLOCATE ( divergence_mmpa(nSteps) ) ALLOCATE ( dextral_mmpa(nSteps) ) ALLOCATE ( elevation(nSteps) ) ALLOCATE ( age_Ma(nSteps) ) ALLOCATE ( class_continuity_c1(nSteps) ) ALLOCATE ( class(nSteps) ) ALLOCATE ( star_c1(nSteps) ) ALLOCATE ( begin_uvec(3, nSteps) ) ALLOCATE ( end_uvec (3, nSteps) ) ALLOCATE ( k_of(nSteps) ) ALLOCATE ( lane_uvecs(3, 5, nSteps) ) ALLOCATE ( subduction_azimuth_1(nSteps) ) ALLOCATE ( subduction_azimuth_2(nSteps) ) ALLOCATE ( Y_factor_of_step(nSteps) ) ALLOCATE ( distance_km_from_step(nSteps) ) ! used only to output distance from classes 1..7 ALLOCATE ( arc_km_from_step(nSteps) ) ! used only to output great-circle-arc-distance from INT event to nearest plate boundary !The following arrays added 2003.06 to support creation of "PB2002_XXX_pure.dat" files: ALLOCATE ( eq_count_over_threshhold(nSteps) ) ALLOCATE ( moment_sum_over_threshhold(nSteps) ) ! second reading is to record the contents: OPEN (UNIT = 2, FILE = "PB2002_steps.dat", STATUS = "OLD", IOSTAT = ios) DO i = 1, nSteps READ (2, 1001) sequence_number(i), step_continuity_c1(i), c5(i), & & old_longitude(i), old_latitude(i), longitude(i), latitude(i), & & length_km(i), azimuth_integer(i), & & velocity_mmpa(i), velocity_azimuth_integer(i), & & divergence_mmpa(i), dextral_mmpa(i), & & elevation(i), age_Ma(i), & & class_continuity_c1(i), class(i), star_c1(i) 1001 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) CALL LonLat_2_Uvec(old_longitude(i), old_latitude(i), uvec) begin_uvec(1:3, i) = uvec(1:3) CALL LonLat_2_Uvec(longitude(i), latitude(i), uvec) end_uvec(1:3, i) = uvec(1:3) !determine class k (= 1..7) for this step: k_of(i) = 0 look_up: DO n = 1, 7 IF (class(i) == class_c3(n)) THEN k_of(i) = n EXIT look_up END IF END DO look_up IF (k_of(i) == 0) THEN WRITE (*, "(' ERROR: Illegal class ',A,' for step ',I6)") class(i), i CALL Pause() STOP END IF IF (class(i) == "SUB") THEN ! determine corners of this lane !define lane (parallelogram of great circle arcs) in counterclockwise direction: azimuth_radians = azimuth_integer(i) * radians_per_degree IF (c5(i)(3:3) == '/') THEN uvec1(1:3) = begin_uvec(1:3, i) uvec2(1:3) = end_uvec(1:3, i) ELSE ! (c5(i)(3:3) == '\') uvec1(1:3) = end_uvec(1:3, i) uvec2(1:3) = begin_uvec(1:3, i) END IF !compute subduction_azimuth_1(i) at location uvec1: iLeftPlate = iPlate(c5(i)(1:2)) iRightPlate = iPlate(c5(i)(4:5)) IF ((iLeftPlate == 0).OR.(iRightPlate == 0)) THEN WRITE (*, "(' Error: Segment title ',A,' includes invalid plate identifier.')") c5(i) CALL Pause() STOP END IF !Euler vector for this step, converting "degrees per Ma" to "km per Ma" or "mm/a" on a unit sphere: tEuler(1:3) = (radius_km / 57.29577) * (Euler(1:3, iLeftPlate) - Euler(1:3, iRightPlate)) !velocity vector of left plate wrt right plate, in mm/a (but in Cartesian coordinates): CALL Cross (tEuler, uvec1, velocity) CALL Local_Theta(uvec1, tvec) ! returns S-pointing vector in tvec velocity_North = -Dot(velocity, tvec) CALL Local_Phi(uvec1, tvec) ! returns E-pointing vector in tvec velocity_East = Dot(velocity, tvec) IF (c5(i)(3:3) == '/') THEN ! left plate is overriding subduction_azimuth_1(i) = Pi + ATAN2(velocity_East, velocity_North) ELSE ! (c5(i)(3:3) == '\'); left plate is subducting subduction_azimuth_1(i) = ATAN2(velocity_East, velocity_North) END IF !compute subduction_azimuth_2(i) at location uvec2 [re-using tEuler]: !velocity vector of left plate wrt right plate, in mm/a (but in Cartesian coordinates): CALL Cross (tEuler, uvec2, velocity) CALL Local_Theta(uvec2, tvec) ! returns S-pointing vector in tvec velocity_North = -Dot(velocity, tvec) CALL Local_Phi(uvec2, tvec) ! returns E-pointing vector in tvec velocity_East = Dot(velocity, tvec) IF (c5(i)(3:3) == '/') THEN ! left plate is overriding subduction_azimuth_2(i) = Pi + ATAN2(velocity_East, velocity_North) ELSE ! (c5(i)(3:3) == '\'); left plate is subducting subduction_azimuth_2(i) = ATAN2(velocity_East, velocity_North) END IF !find corner points of this lane: CALL Turn_To (azimuth_radians = subduction_azimuth_1(i) + Pi, & & base_uvec = uvec1, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c1_uvec) lane_uvecs(1:3, 1, i) = c1_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2(i) + Pi, & & base_uvec = uvec2, & & far_radians = SUB_seaward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c2_uvec) lane_uvecs(1:3, 2, i) = c2_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_2(i), & & base_uvec = uvec2, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c3_uvec) lane_uvecs(1:3, 3, i) = c3_uvec(1:3) CALL Turn_To (azimuth_radians = subduction_azimuth_1(i), & & base_uvec = uvec1, & & far_radians = SUB_landward_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = c4_uvec) lane_uvecs(1:3, 4, i) = c4_uvec(1:3) lane_uvecs(1:3, 5, i) = lane_uvecs(1:3, 1, i) lane_width_km = length_km(i) * ABS(SIN(subduction_azimuth_1(i) - (azimuth_integer(i) * radians_per_degree))) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" basis_length_km = lane_width_km ELSE basis_length_km = length_km(i) lane_uvecs(1:3, 1:5, i) = 0.0 ! just to be tidy; should never be used! END IF ! class(i) == "SUB", or not !Precompute Y_factor_of_step, to save time: integral = 0.0 s1 = 0.0 - 3.0 * alongStep_sigma_km s2 = basis_length_km + 3.0 * alongStep_sigma_km ds = (s2 - s1) / 100 DO n = 0, 100 s_km = s1 + (n / 100.0) * (s2 - s1) integral = integral + ds * ANORDF(s_km / alongStep_sigma_km) * ANORDF((basis_length_km - s_km) / alongStep_sigma_km) END DO Y_factor_of_step(i) = 1 / integral END DO ! i = 1, nSteps CLOSE (2) WRITE (*, "(' Finished reading and memorizing PB2002_steps.dat.')") ! Read (3): PB2002_orogens.dig (and store in memory) WRITE (*, "(' Begin reading and memorizing PB2002_orogens.dat...')") OPEN (UNIT = 3, FILE = 'PB2002_orogens.dig', STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: input file PB2002_orogens could not be opened.')") CALL Pause() END IF DO j = 1, nOrogens READ (3,"(A)", IOSTAT = ios) c27 IF (ios /= 0) THEN WRITE (*,"(' ERROR: could not read (all?) of PB2002_orogens.dig')") CALL Traceback() STOP END IF WRITE (*, "(' ',A)") TRIM(c27) points: DO i = 1, mostInOneOrogen + 1 READ (3, *, IOSTAT = ios) t_longitude, t_latitude IF (ios /= 0) THEN ! hit "*** end of line segment ***" nInEachOrogen(j) = i - 1 EXIT points END IF CALL LonLat_2_Uvec(t_longitude, t_latitude, uvec) orogen_uvecs(1:3, i, j) = uvec(1:3) END DO points END DO ! j = 1, norogens CLOSE (3) WRITE (*, "(' Finished reading and memorizing PB2002_orogens.dat.')") ! open the 16 .eqc files for the output that will soon follow ------------------------- 1100 WRITE (*, "(' Which type of seismic catalog will be processed in this run?')") WRITE (*, "(' (must always be a .eqc file, but formats vary for focal mechanism)')") WRITE (*, "(' (1) Harvard CMT')") WRITE (*, "(' (2) Pacheco and Sykes [1992]')") WRITE (*, "(' Enter integer to select: '\)") READ (*, *) catalog_ID IF ((catalog_ID < 1).OR.(catalog_ID > 2)) THEN WRITE (*, "(' ERROR: Illegal integer value; try again...')") GO TO 1100 END IF IF (catalog_ID == 1) THEN prefix = "CMT" ELSE IF (catalog_ID == 2) THEN prefix = "Pacheco" ELSE prefix = "other" END IF WRITE (*, "(' Enter name of seismic catalog (.eqc) file: '\)") READ (*, "(A)") catalog_file_name WRITE (*, "(' Opening 16 new .eqc files to be gradually filled with EQs...')") !class 1: XXXX_PB2002_CCB_pure.eqc and XXXX_PB2002_CCB_all.eqc Continental Convergent Boundary out_file_name = TRIM(prefix) // "_PB2002_CCB_pure.eqc" OPEN (UNIT = 11, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_CCB_all.eqc" OPEN (UNIT = 12, FILE = out_file_name) !class 2: XXXX_PB2002_CTF_pure.eqc and XXXX_PB2002_CTF_all.eqc Continental Transform Fault out_file_name = TRIM(prefix) // "_PB2002_CTF_pure.eqc" OPEN (UNIT = 21, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_CTF_all.eqc" OPEN (UNIT = 22, FILE = out_file_name) !class 3: XXXX_PB2002_CRB_pure.eqc and XXXX_PB2002_CRB_all.eqc Continental Rift Boundary out_file_name = TRIM(prefix) // "_PB2002_CRB_pure.eqc" OPEN (UNIT = 31, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_CRB_all.eqc" OPEN (UNIT = 32, FILE = out_file_name) !class 4: XXXX_PB2002_OSR_pure.eqc and XXXX_PB2002_OSR_all.eqc Oceanic Spreading Ridge out_file_name = TRIM(prefix) // "_PB2002_OSR_pure.eqc" OPEN (UNIT = 41, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_OSR_all.eqc" OPEN (UNIT = 42, FILE = out_file_name) !class 5: XXXX_PB2002_OTF_pure.eqc and XXXX_PB2002_OTF_all.eqc Oceanic Transform Fault out_file_name = TRIM(prefix) // "_PB2002_OTF_pure.eqc" OPEN (UNIT = 51, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_OTF_all.eqc" OPEN (UNIT = 52, FILE = out_file_name) !class 6: XXXX_PB2002_OCB_pure.eqc and XXXX_PB2002_OCB_all.eqc Oceanic Convergent Boundary out_file_name = TRIM(prefix) // "_PB2002_OCB_pure.eqc" OPEN (UNIT = 61, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_OCB_all.eqc" OPEN (UNIT = 62, FILE = out_file_name) !class 7: XXXX_PB2002_SUB_pure.eqc and XXXX_PB2002_SUB_all.eqc SUBduction zone out_file_name = TRIM(prefix) // "_PB2002_SUB_pure.eqc" OPEN (UNIT = 71, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_SUB_all.eqc" OPEN (UNIT = 72, FILE = out_file_name) !class 0 (or "10" for I/O purposes): XXXX_PB2002_INT_pure.eqc and XXXX_PB2002_INT_all.eqc INTraplate (plate-INTerior) out_file_name = TRIM(prefix) // "_PB2002_INT_pure.eqc" OPEN (UNIT = 101, FILE = out_file_name) out_file_name = TRIM(prefix) // "_PB2002_INT_all.eqc" OPEN (UNIT = 102, FILE = out_file_name) !common format for subcatalogs: 3002 FORMAT (A, & & I5,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & A21, & ! up to here, a normal .eqc format (with or w/o FPS) [ WRITE (...)..., TRIM(appended_data)! ] & 1X,A1,A3,A1,1X,I5, & ! additional information: " *SUB* 4392" & 7(1X,I3), & ! additional information: " 1 2 3 4 5 6 79" & F8.1) ! additional information: distance_km from most relevant step ! read (4) XXXX_shallow.eqc and process events (don't store them) ------------------------- count_all = 0 ! each of the (0:7) values, for "all" subcatalogs count_pure = 0 ! each of the (0:7) values, for "pure" subcatalogs eq_count_over_threshhold = 0 ! initialize count for each plate boundary step moment_sum_over_threshhold = 0 ! initialize sum for each plate boundary step WRITE (*, "(' Reading and classifying ', A, ' catalog...')") TRIM(prefix) iEvent = 0 ! just for displaying progress OPEN (UNIT = 4, FILE = catalog_file_name, STATUS = "OLD", IOSTAT = ios, PAD = "YES") ! PAD needed to read appended_data IF (ios /= 0) THEN WRITE (*,"(' ERROR: ', A, ' not available to read in.')") TRIM(catalog_file_name) CALL Pause() END IF WRITE (*, *) CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) IF (Monte_Carlo_method) THEN CALL RNOPT(1) ! Select type of random number generator: 1 = congruential, multiplier 16807, no shuffling [default]. iseed = (datetimenumber(1) - 2001) * 31557600 + & ! not precise; ignores leap years (datetimenumber(2) - 1) * 2629800 + & ! not precise; all months are not really the same length! (datetimenumber(3) - 1) * 86400 + & datetimenumber(5) * 3600 + & datetimenumber(6) * 60 + & datetimenumber(7) iseed = MOD(iseed, 2147483646) CALL RNSET(iseed) ! Seed the random number generator; iseed must be in the range (0, 2147483646) ! which is equivalent to the number of seconds in a span of 68 years or less. END IF next_event: DO READ (4, 3001, IOSTAT = ios) catalog_name, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data 3001 FORMAT ( A9, & & I5,'.',A2,'.',A2, 1X, & ! using I5 for negative years (B.C.) & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & A ) IF (ios /= 0) EXIT next_event iEvent = iEvent + 1 WRITE (*, "('+ processing event ',I10)") iEvent ! ===== heart of this program ================================================ !test whether earthquake is in any orogen: CALL LonLat_2_Uvec(eq_Elon, eq_Nlat, eq_uvec) pure_epicenter = .TRUE. ! initialization, before testing check_orogens: DO j = 1, nOrogens angle_sum = 0.0 ! initialize sum of angles subtended by orogen steps, as seen from test point: DO i = 2, nInEachOrogen(j) uvec1(1:3) = orogen_uvecs(1:3, i-1, j) uvec2(1:3) = orogen_uvecs(1:3, i , j) t1 = Relative_Compass(from_uvec = eq_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = eq_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 ! i = 2, nInEachOrogen(j) inside = ((angle_sum < -6.0).AND.(angle_sum > -6.6)) ! 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 pure_epicenter = .FALSE. EXIT check_orogens END IF END DO check_orogens ! j = 1, nOrogens IF (pure_epicenter) THEN pure_epicenter_c1 = ' ' ELSE ! in some orogen pure_epicenter_c1 = '*' END IF highest_P_prime = 0.0 ! initialization; if still true after loop, then EQ is an INT best_i = 0 ! will remember i index of step associated with highest_P_prime best_k = 0 ! will remember k index of class of the step associated with highest_P_prime, if any; ! if no association found, best_k remains 0, which is a signal for an INT EQ. highest_P_prime_in_class = 0.0 ! for all 7 classes best_i_in_class = 0 ! for all 7 classes sum_P_prime_in_class = 0.0 ! for all 7 classes sum_P_prime_all_classes = 0.0 distance_km_from_step = Pi * radius_km ! whole array (1:nSteps) arc_km_from_step = Pi * radius_km ! whole array (1:nSteps) DO i = 1, nSteps k = k_of(i) !------------------------------------------------------------------------------ ! A factor: relative number (or rate) of EQs (over m_t) produced by this step: A = N_over_mt(k) * (length_km(i) / L_km(k)) * (velocity_mmpa(i) / V_mmpa(k)) !------------------------------------------------------------------------------ !distances to end points: uvec1(1:3) = begin_uvec(1:3, i) to_uvec1_radians = Arc(eq_uvec, uvec1) uvec2(1:3) = end_uvec(1:3, i) to_uvec2_radians = Arc(eq_uvec, uvec2) !great-circle-arc distance computed ONLY to output for earthquakes that end up as INT's: arc_km_from_step(i) = radius_km * MIN(to_uvec1_radians, to_uvec2_radians) !Note that this may be reduced below to distance_km, if within end-points. IF (k < 7) THEN ! step is not a SUB; use normal distance/offset logic: !distance and offset measured via perpendicular great circle: CALL Cross(uvec1, uvec2, tvec) CALL Make_Uvec(tvec, pole_uvec) ! pole of the great-circle for this step offside_radians = ABS(Arc(pole_uvec, eq_uvec) - Pi_over_2) distance_km = radius_km * offside_radians ! always positive offset_km = distance_km ! no distinction for classes other than SUB !s_km measured along length of step (using longitudes relative to pole of step's great circle): t1 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec1) t2 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec2) !note: normally, t2 is slightly less than t1 IF (t2 > t1) t2 = t2 - Two_Pi angle = Relative_Compass(from_uvec = pole_uvec, to_uvec = eq_uvec) ! to the earthquake !check that angle is on the same cycle as t1, t2: midpoint_compass = (t1 + t2) / 2 IF ((angle - midpoint_compass) > Pi) angle = angle - Two_Pi IF ((midpoint_compass - angle) > Pi) angle = angle + Two_Pi s_radians = t1 - angle s_km = s_radians * radius_km !-------------------------------------------------------------------------------------------------------- ! X factor = cross-strike component of spatial PDF B: IF (offset_km > h_km(k)) THEN X = 0.0 ELSE X = (1.0 / 0.95) *(0.797885 / h_km(k)) * EXP(-2 * (offset_km / h_km(k))**2) END IF !-------------------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------------------------- ! Y factor causes PDF to fall off to ~0.5 as ends of step are approached (and even more beyond the ends) IF ((s_km > (-3.0 * alongStep_sigma_km)).AND.(s_km < (length_km(i) + 3.0 * alongStep_sigma_km))) THEN ! use complex formula: Y = Y_factor_of_step(i) * ANORDF(s_km / alongStep_sigma_km) * ANORDF((length_km(i) - s_km) / alongStep_sigma_km) distance_km_from_step(i) = distance_km ! save for possible output, if this step is selected ELSE ! too far away; don't risk numerical problems within ANORDF: Y = 0.0 END IF !------------------------------------------------------------------------------------------------------- IF ((s_km > 0.0).AND.(s_km < length_km(i))) THEN ! use perpendicular distance as shortest great-circle-arc: arc_km_from_step(i) = distance_km !Notes: Used only for output in the case that this earthquake ends up as an INT. ! Also, arc_km_from_step(i) was previously set to the lesser of the distance to the 2 endpoints of this step. END IF ELSE ! this is a SUB (k = 7) step; special offset logic is used, and - distances/offsets mean seaward of the trench! proximal = .TRUE. ! unless (usually) negated below: DO n = 1, 4 ! sides uvec1(1:3) = lane_uvecs(1:3, n, i) uvec2(1:3) = lane_uvecs(1:3, n+1, i) CALL Cross (uvec1, uvec2, tvec) ! direction toward pole is direction into the lane CALL Make_Uvec(tvec, inward_uvec) scalar_product = Dot(eq_uvec, inward_uvec) dot_side(n) = scalar_product proximal = proximal .AND. (scalar_product > -0.01) ! 63.7 km tolerance END DO ! n = 1, 4 sides ! Note: EQ is inside the lane if all 4 values in dot_side(1:4) are postive; ! these values are sines of angles from the sides, measured inwards. IF (proximal) THEN ! compute + or - distance/offset from trench, in km (even if slightly to one side of the lane): IF (c5(i)(3:3) == '/') THEN uvec1(1:3) = begin_uvec(1:3, i) uvec2(1:3) = end_uvec(1:3, i) ELSE ! (c5(i)(3:3) == '\') uvec1(1:3) = end_uvec(1:3, i) uvec2(1:3) = begin_uvec(1:3, i) END IF fraction_toward_uvec1 = dot_side(2) / (dot_side(2) + dot_side(4)) ! where all dot's are positive if inside; ! note that denominator is positive even if EQ (slightly) outside lane. fraction_toward_uvec1 = MAX(0.0, MIN(1.0, fraction_toward_uvec1)) ! This kludge is necessary to prevent pathalogical ! behavior when SUB steps are very short and oblique. tvec(1:3) = fraction_toward_uvec1 * uvec1(1:3) + & & (1.0 - fraction_toward_uvec1) * uvec2(1:3) CALL Make_Uvec(tvec, trench_uvec) ! at same relative left-right position in this lane distance_km = radius_km * Arc(eq_uvec, trench_uvec) ! always positive; must apply sign now: from_trench_azimuth_radians = Relative_Compass (from_uvec = trench_uvec, to_uvec = eq_uvec) IF (COS(from_trench_azimuth_radians - subduction_azimuth_1(i)) < 0.0) THEN ! eq_is on seaward side of trench; convert distance to negative: distance_km = -distance_km END IF !Note that "distance" is positive on landward side of trench, but negative on seaward side. distance_km_from_step(i) = distance_km ! save for possible output offset_km = distance_km - SUB_peakSeismicity_km ! preserving possibility of negative values; just shifting origin !-------------------------------------------------------------------------------------------------------- ! X factor = along-lane component of spatial PDF function B: IF (offset_km >= 0.0) THEN ! landward side of seismicity peak: IF (distance_km <= SUB_landward_halfwidth_km) THEN X = (1.0 / 0.95) * & & (2 * SUB_landward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_landward_sigma_km)) * & & EXP(-0.5 * (offset_km / SUB_landward_sigma_km)**2) ELSE X = 0.0 END IF ELSE ! seaward side of seismicity peak (but may not be seaward of the trench!) IF (distance_km >= -SUB_seaward_halfwidth_km) THEN X = (1.0 / 0.95) * & & (2 * SUB_seaward_sigma_km / (SUB_landward_sigma_km + SUB_seaward_sigma_km)) * & & (0.797885 / (2.0 * SUB_seaward_sigma_km )) * & & EXP(-0.5 * (offset_km / SUB_seaward_sigma_km)**2) ELSE X = 0.0 END IF END IF ! which side of seismicity peak? !-------------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------------- ! Y factor = cross-lane (orthogonal to X factor) component of spatial PDF function B: ! To evaluate Y(s_km), note that length_km(i) is not the basis length that we want to use, ! but rather the lane width: lane_width_km = radius_km * (dot_side(2) + dot_side(4)) lane_width_km = MAX(lane_width_km, 1.0) ! protect against singularity in case of extremely oblique "subduction" s_km = radius_km * dot_side(4) ! distance (if +) into lane, starting on the uvec1 end IF ((s_km > (-3.0 * alongStep_sigma_km)).AND.(s_km < (lane_width_km + 3.0 * alongStep_sigma_km))) THEN ! use complex formula: Y = Y_factor_of_step(i) * ANORDF(s_km / alongStep_sigma_km) * ANORDF((lane_width_km - s_km) / alongStep_sigma_km) ELSE ! too far away; don't risk numerical problems within ANORDF: Y = 0.0 END IF !-------------------------------------------------------------------------------------------------------- ELSE ! not proximal; punt X = 0.0 Y = 0.0 !Note: distance_km_from_step(i) was initialized as Pi * radius_km END IF END IF ! "normal" offset/distance logic (k = 1-6), or special logic for SUB steps (k = 7) B = X * Y ! spatial PDF is a product of cross-trace function and along-trace function IF ((A * B) > 0.0) THEN ! continue investigating this step; P_prime should be positive !------------------------------------------------------------------------------------ !Compute SUM (j=1,5) {C * D * E} for this step !extract depth reliability and mechanism classification, depending on catalog: IF (catalog_ID == 1) THEN ! CMT, has orientations of principal strain axes eq_depth_known = (eq_depth_int /= 33) READ (appended_data, "(3(1X,I2,1X,I3))", IOSTAT = ios) e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth IF (ios == 0) THEN mechanism_known = .TRUE. actual_TpaPpa_degrees(1) = e3_plunge ! plunge of T axis, in degrees actual_TpaPpa_degrees(2) = e3_azimuth ! azimuth of T axix, in degrees actual_TpaPpa_degrees(3) = e1_plunge ! plunge of P axis, in degrees actual_TpaPpa_degrees(4) = e1_azimuth ! azimuth of P axis, in degrees ELSE mechanism_known = .FALSE. WRITE (*, "(' Bad mechanism found in CMT-format file: ',A)") TRIM(appended_data) CALL Pause() STOP END IF ELSE IF (catalog_ID == 2) THEN ! Pacheco and Sykes [1992], has (some) type indicators eq_depth_known = (eq_depth_int /= 0) c3 = appended_data ! left-most 3 bytes, such as: " n", " u", " ts" IF (INDEX(c3, 'n') > 0) THEN !this is a NORMAL faulting event mechanism_known = .TRUE. unknown = .FALSE.; normal = .TRUE.; strikeslip = .FALSE.; thrust = .FALSE. ELSE IF ((INDEX(c3, 't') > 0).OR.(INDEX(c3, 'r') > 0).OR.(INDEX(c3, 'c') > 0)) THEN !this is a THRUST or REVERSE or COMPRESSIONAL faulting event mechanism_known = .TRUE. unknown = .FALSE.; normal = .FALSE.; strikeslip = .FALSE.; thrust = .TRUE. ELSE IF (INDEX(c3, 's') > 0) THEN !this is a STRIKE-SLIP faulting event mechanism_known = .TRUE. unknown = .FALSE.; normal = .FALSE.; strikeslip = .TRUE.; thrust = .FALSE. ELSE mechanism_known = .FALSE. unknown = .TRUE.; normal = .FALSE.; strikeslip = .FALSE.; thrust = .FALSE. END IF ELSE WRITE (*, "(' Rules not known for interpreting earthquake type in this catalog.')") CALL Traceback() END IF ! CMT, Pacheco and Sykes [1992], or other type of catalog sum_CDE = 0.0 ! initializing DO j = 1, 5 ! model earthquakes for this step (to skip unneeded ones, set C = 0.0) ! Assignment of C, the relative probability of this model EQ: SELECT CASE (k) CASE (1, 6) ! CCB and OCB: closing_angle_radians = ATAN2(-divergence_mmpa(i), dextral_mmpa(i)) ! syntax: ATAN2(y, x) !which will be close to +Pi/2 if convergence_mmpa >> ABS(dextral_mmpa) SELECT CASE (j) CASE (1) ! oblique thrust on fault dipping to left C = 0.333 * (1.0 - ABS(COS(closing_angle_radians))) CASE (2) ! oblique thrust on fault dipping to right C = 0.333 * (1.0 - ABS(COS(closing_angle_radians))) CASE (3) ! pure thrust with shortening perpendicular to plate boundary: C = 0.333 * (1.0 - ABS(COS(closing_angle_radians))) !Note: Distinction between the j=1 and j=2 models is academic for near-othogonal convergence. CASE (4) ! pure strike-slip on plane striking parallel to plate boundary: C = ABS(COS(closing_angle_radians)) CASE (5) ! (not used) C = 0.0 END SELECT ! on j CASE (2, 5) ! CTF and OTF: SELECT CASE (j) CASE (1) ! strike-slip on vertical plane of plate boundary: C = 0.830 ! = 1.00 - C_22 - C_23 CASE (2) ! normal on plane striking parallel to plate boundary: C = 0.085 ! limited, due to definition that CTF, OTF boundaries are within 20 deg. of plate velocity CASE (3) ! thrust on plane striking parallel to plate boundary: C = 0.085 ! limited, due to definition that CTF, OTF boundaries are within 20 deg. of plate velocity CASE (4) ! (not used) C = 0.0 CASE (5) ! (not used) C = 0.0 END SELECT ! on j CASE (3, 4) ! CRB and OSR: opening_angle_radians = ATAN2(divergence_mmpa(i), dextral_mmpa(i)) ! syntax: ATAN2(y, x) !which will be close to +Pi/2 if divergence_mmpa >> ABS(dextral_mmpa) SELECT CASE (j) CASE (1) ! oblique normal faulting, on plane to left of trace C = 0.333 * (1.0 - ABS(COS(opening_angle_radians))) CASE (2) ! oblique normal faulting, on plane to right of trace C = 0.333 * (1.0 - ABS(COS(opening_angle_radians))) !Note: Distinction between the j=1 and j=2 models is academic for near-othogonal spreading. CASE (3) ! pure normal faulting C = 0.333 * (1.0 - ABS(COS(opening_angle_radians))) CASE (4) ! pure strike-slip C = ABS(COS(opening_angle_radians)) CASE (5) ! (not used) C = 0.0 END SELECT ! on j CASE (7) ! SUB: SELECT CASE (j) CASE (1) ! oblique slip on main inter-plate thrust fault: C = 1689./6274. ! = 1.0 - C_72 - C_73 - C_74 CASE (2) ! pure dip-slip on main inter-plate thrust fault: C = 1689./6274. ! = 1.0 - C_72 - C_73 - C_74 CASE (3) ! strike-slip parallel to arc: C = 690./6274. ! based on CMT_PB2002_pure_SUB_breakdown.xls (from old EQ_classification.f90) CASE (4) ! down-dip extension due to bending and/or stress-guide: C = 1103./6274. ! based on CMT_PB2002_pure_SUB_breakdown.xls (from old EQ_classification.f90) CASE (5) ! down-dip compression due to bending and/or stress-guide: C = 1103./6274. ! arbitrarily assumed = to C_74, since hard to distinguish from C_71 and C_72 END SELECT ! on j END SELECT ! on k; END of C (relative probabilities of models) section IF (C > 0.0) THEN ! this model has a chance of happening ! Compute D, the depth PDF of this model EQ evaluated at the actual EQ depth: IF (eq_depth_known) THEN SELECT CASE (k) CASE (1:6) ! all EXCEPT SUB are the same, regardless of model (j): argument = (class123456_cutoff_depth_km - eq_depth_int) / earthquake_depth_sigma_km IF (argument > -3.0) THEN ! safe to call ANORDF: D = (1. / class123456_cutoff_depth_km) * ANORDF( argument ) ELSE ! ANORDF might underflow D = 0.0 END IF CASE (7) ! SUB is different quadratic = (slab_depth_under_arc_km - trench_depth_km - (slope_at_trench*arc_distance_km)) / arc_distance_km**2 distance_at_peak_km = -slope_at_trench / (2.0 * quadratic) minimum_depth_km = trench_depth_km + distance_at_peak_km * slope_at_trench + quadratic * distance_at_peak_km**2 IF (distance_km > distance_at_peak_km) THEN interplate_thrust_depth_km = trench_depth_km + distance_km * slope_at_trench + quadratic * distance_km**2 ELSE interplate_thrust_depth_km = minimum_depth_km END IF interplate_thrust_depth_sigma_km = interplate_thrust_depth_sigma_fraction * interplate_thrust_depth_km !Note: Use SUB_quadratic_profile.xls to experiment with these functions and parameters. SELECT CASE (j) CASE (1:2) ! main inter-plate thrust fault D = (1.0 / SQRT( Two_Pi * (earthquake_depth_sigma_km**2 + interplate_thrust_depth_sigma_km**2) )) * & & EXP( -(eq_depth_int - interplate_thrust_depth_km)**2 / & & (2 * (earthquake_depth_sigma_km**2 + interplate_thrust_depth_sigma_km**2)) ) CASE (3) ! strike-slip parallel to arc (same as classes 1-6 above) argument = (class123456_cutoff_depth_km - eq_depth_int) / earthquake_depth_sigma_km IF (argument > -3.0) THEN ! safe to call ANORDF: D = (1. / class123456_cutoff_depth_km) * ANORDF( argument ) ELSE ! ANORDF might underflow D = 0.0 END IF CASE (4:5) ! down-dip extension or compression due to bending and/or stress-guid argument = (eq_depth_int - interplate_thrust_depth_km) / earthquake_depth_sigma_km IF (argument > -3.0) THEN ! safe to call ANORDF: D = (1. / 60.0) * ANORDF( argument ) !where 60 km is the typical thickness of the slab stress guide ELSE ! ANORDF might underflow D = 0.0 END IF END SELECT ! on j END SELECT ! on k ELSE ! EQ depth is NOT known: D = 1./70. ! in units of /km END IF ! end of D (depth PDF) section IF (D > 0.0) THEN ! this EQ depth could match this model ! Compute E, the angular PDF of this model EQ evaluated at the actual EQ mechanism orientation: IF (mechanism_known) THEN ! use it to evaluate E: SELECT CASE (catalog_ID) CASE (1) ! CMT: SELECT CASE(k) ! boundary class CASE (1, 6) ! CCB or OCB: SELECT CASE (j) CASE (1) ! oblique thrust dipping 20 degrees to left of trace: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = 20.0 dip_is_to_right = .FALSE. relative_velocity_azimuth_degrees = velocity_azimuth_integer(i) relative_velocity_is_divergent = .FALSE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output model_TpaPpa_degrees(1) = e3_plunge model_TpaPpa_degrees(2) = e3_azimuth model_TpaPpa_degrees(3) = e1_plunge model_TpaPpa_degrees(4) = e1_azimuth CASE (2) ! oblique thrust dipping 20 degrees to right of trace: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = 20.0 dip_is_to_right = .TRUE. relative_velocity_azimuth_degrees = velocity_azimuth_integer(i) relative_velocity_is_divergent = .FALSE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output model_TpaPpa_degrees(1) = e3_plunge model_TpaPpa_degrees(2) = e3_azimuth model_TpaPpa_degrees(3) = e1_plunge model_TpaPpa_degrees(4) = e1_azimuth CASE (3) ! pure thrust with shortening perpendicular to plate boundary: model_TpaPpa_degrees(1) = 90 model_TpaPpa_degrees(2) = 0 model_TpaPpa_degrees(3) = 0 model_TpaPpa_degrees(4) = azimuth_integer(i) + 90 CASE (4) ! pure strike-slip on plane striking parallel to plate boundary: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(3) = 0 IF (dextral_mmpa(i) > 0.0) THEN ! right-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 135 model_TpaPpa_degrees(4) = azimuth_integer(i) + 45 ELSE ! left-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 45 model_TpaPpa_degrees(4) = azimuth_integer(i) + 135 END IF END SELECT ! on j CASE (2, 5) ! CTF or OTF: SELECT CASE (j) CASE (1) ! strike-slip on vertical plane of plate boundary: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(3) = 0 IF (dextral_mmpa(i) > 0.0) THEN ! right-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 135 model_TpaPpa_degrees(4) = azimuth_integer(i) + 45 ELSE ! left-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 45 model_TpaPpa_degrees(4) = azimuth_integer(i) + 135 END IF CASE (2) ! normal on plane striking parallel to plate boundary: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(2) = azimuth_integer(i) + 90 model_TpaPpa_degrees(3) = 90 model_TpaPpa_degrees(4) = 0 CASE (3) ! thrust on plane striking parallel to plate boundary: model_TpaPpa_degrees(1) = 90 model_TpaPpa_degrees(2) = 0 model_TpaPpa_degrees(3) = 0 model_TpaPpa_degrees(4) = azimuth_integer(i) + 90 END SELECT ! on j CASE (3, 4) ! CRB or OSR: SELECT CASE (j) CASE (1) ! oblique normal faulting, with dip 55 degrees to left of trace: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = 55.0 dip_is_to_right = .FALSE. relative_velocity_azimuth_degrees = velocity_azimuth_integer(i) relative_velocity_is_divergent = .TRUE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output model_TpaPpa_degrees(1) = e3_plunge model_TpaPpa_degrees(2) = e3_azimuth model_TpaPpa_degrees(3) = e1_plunge model_TpaPpa_degrees(4) = e1_azimuth CASE (2) ! oblique normal faulting, with dip 55 degrees to right of trace: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = 55.0 dip_is_to_right = .TRUE. relative_velocity_azimuth_degrees = velocity_azimuth_integer(i) relative_velocity_is_divergent = .TRUE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output model_TpaPpa_degrees(1) = e3_plunge model_TpaPpa_degrees(2) = e3_azimuth model_TpaPpa_degrees(3) = e1_plunge model_TpaPpa_degrees(4) = e1_azimuth CASE (3) ! pure normal faulting, with extension perpendicular to plate boundary: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(2) = azimuth_integer(i) + 90 model_TpaPpa_degrees(3) = 90 model_TpaPpa_degrees(4) = 0 CASE (4) ! partitioned strike-slip on plane striking parallel to plate boundary: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(3) = 0 IF (dextral_mmpa(i) > 0.0) THEN ! right-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 135 model_TpaPpa_degrees(4) = azimuth_integer(i) + 45 ELSE ! left-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 45 model_TpaPpa_degrees(4) = azimuth_integer(i) + 135 END IF END SELECT ! on j CASE (7) ! SUB: quadratic = (slab_depth_under_arc_km - trench_depth_km - (slope_at_trench*arc_distance_km)) / arc_distance_km**2 distance_at_peak_km = -slope_at_trench / (2.0 * quadratic) IF (distance_km > distance_at_peak_km) THEN velocity_dip_radians = ATAN(slope_at_trench + 2.0 * quadratic * distance_km) ELSE velocity_dip_radians = 0.0 END IF velocity_dip_degrees = NINT(velocity_dip_radians * degrees_per_radian) IF ((velocity_dip_degrees < 0).OR.(velocity_dip_degrees >= 90)) THEN WRITE (*, "(' Bad velocity_dip_degrees: ', I6)") velocity_dip_degrees CALL Traceback() STOP END IF !velocity_dip is measured along the velocity vector (like the quadratic profile); !crossStrike_dip is steeper in the case of oblique subduction, equal for orthogonal subduction: crossStrike_dip_radians = ATAN(TAN(velocity_dip_radians) / & & ABS(SIN((velocity_azimuth_integer(i)-azimuth_integer(i))*radians_per_degree))) crossStrike_dip_degrees = NINT(crossStrike_dip_radians * degrees_per_radian) IF ((crossStrike_dip_degrees < 0).OR.(crossStrike_dip_degrees > 90)) THEN WRITE (*, "(' Bad crossStrike_dip_degrees: ', I6)") crossStrike_dip_degrees CALL Traceback() STOP END IF crossStrike_dip_degrees = MIN(crossStrike_dip_degrees, 89) IF (c5(i)(3:3) == '/') THEN ! velocity_azimuth_integer(i) refers to the LEFT plate, going along the boundary: dip_is_to_right = .FALSE. slab_velocity_azimuth_integer = velocity_azimuth_integer(i) + 180 slab_downdip_azimuth_integer = azimuth_integer(i) - 90 ELSE ! c5(3:3) == '\'; left plate IS the slab: dip_is_to_right = .TRUE. slab_velocity_azimuth_integer = velocity_azimuth_integer(i) slab_downdip_azimuth_integer = azimuth_integer(i) + 90 END IF SELECT CASE (j) CASE (1) ! oblique slip on main inter-plate thrust fault: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = MAX(crossStrike_dip_degrees, 5) ! in case EQ is mislocated in outer rise relative_velocity_azimuth_degrees = velocity_azimuth_integer(i) relative_velocity_is_divergent = .FALSE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output T_plunge_degrees = e3_plunge T_azimuth_degrees = e3_azimuth P_plunge_degrees = e1_plunge P_azimuth_degrees = e1_azimuth IF (T_plunge_degrees > 90) THEN T_plunge_degrees = 180 - T_plunge_degrees T_azimuth_degrees = T_azimuth_degrees + 180 P_plunge_degrees = -P_plunge_degrees P_azimuth_degrees = P_azimuth_degrees - 180 END IF IF ((T_plunge_degrees < 0).OR.(T_plunge_degrees > 90).OR. & (P_plunge_degrees < 0).OR.(P_plunge_degrees > 90)) THEN WRITE (*, "(' Bad plunge')") CALL Traceback() STOP END IF !only truncate to INTEGER*2 (losing sign) *after* testing! model_TpaPpa_degrees(1) = T_plunge_degrees model_TpaPpa_degrees(2) = T_azimuth_degrees model_TpaPpa_degrees(3) = P_plunge_degrees model_TpaPpa_degrees(4) = P_azimuth_degrees CASE (2) ! pure thrust on main inter-plate thrust fault: fault_strike_degrees = azimuth_integer(i) fault_dip_degrees = MAX(crossStrike_dip_degrees, 5) ! in case EQ is mislocated in outer rise relative_velocity_azimuth_degrees = fault_strike_degrees + 90 relative_velocity_is_divergent = .FALSE. CALL Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth, e1_plunge, & & e2_azimuth, e2_plunge, & & e3_azimuth, e3_plunge) ! output T_plunge_degrees = e3_plunge T_azimuth_degrees = e3_azimuth P_plunge_degrees = e1_plunge P_azimuth_degrees = e1_azimuth IF (T_plunge_degrees > 90) THEN T_plunge_degrees = 180 - T_plunge_degrees T_azimuth_degrees = T_azimuth_degrees + 180 P_plunge_degrees = -P_plunge_degrees P_azimuth_degrees = P_azimuth_degrees - 180 END IF IF ((T_plunge_degrees < 0).OR.(T_plunge_degrees > 90).OR. & (P_plunge_degrees < 0).OR.(P_plunge_degrees > 90)) THEN WRITE (*, "(' Bad plunge')") CALL Traceback() STOP END IF !only truncate to INTEGER*2 (losing sign) *after* testing! model_TpaPpa_degrees(1) = T_plunge_degrees model_TpaPpa_degrees(2) = T_azimuth_degrees model_TpaPpa_degrees(3) = P_plunge_degrees model_TpaPpa_degrees(4) = P_azimuth_degrees CASE (3) ! strike-slip parallel to arc: model_TpaPpa_degrees(1) = 0 model_TpaPpa_degrees(3) = 0 IF (dextral_mmpa(i) > 0.0) THEN ! right-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 135 model_TpaPpa_degrees(4) = azimuth_integer(i) + 45 ELSE ! left-lateral: model_TpaPpa_degrees(2) = azimuth_integer(i) + 45 model_TpaPpa_degrees(4) = azimuth_integer(i) + 135 END IF CASE (4) ! down-dip extension due to bending and/or stress-guide: T_plunge_degrees = crossstrike_dip_degrees T_azimuth_degrees = slab_downdip_azimuth_integer P_plunge_degrees = 90 - crossstrike_dip_degrees P_azimuth_degrees = slab_downdip_azimuth_integer + 180 IF ((T_plunge_degrees < 0).OR.(T_plunge_degrees > 90).OR. & (P_plunge_degrees < 0).OR.(P_plunge_degrees > 90)) THEN WRITE (*, "(' Bad plunge')") CALL Traceback() STOP END IF !only truncate to INTEGER*2 (losing sign) *after* testing! model_TpaPpa_degrees(1) = T_plunge_degrees model_TpaPpa_degrees(2) = T_azimuth_degrees model_TpaPpa_degrees(3) = P_plunge_degrees model_TpaPpa_degrees(4) = P_azimuth_degrees CASE (5) ! down-dip compression due to bending and/or stress-guide: T_plunge_degrees = 90 - crossstrike_dip_degrees T_azimuth_degrees = slab_downdip_azimuth_integer + 180 P_plunge_degrees = crossstrike_dip_degrees P_azimuth_degrees = slab_downdip_azimuth_integer IF ((T_plunge_degrees < 0).OR.(T_plunge_degrees > 90).OR. & (P_plunge_degrees < 0).OR.(P_plunge_degrees > 90)) THEN WRITE (*, "(' Bad plunge')") CALL Traceback() STOP END IF !only truncate to INTEGER*2 (losing sign) *after* testing! model_TpaPpa_degrees(1) = T_plunge_degrees model_TpaPpa_degrees(2) = T_azimuth_degrees model_TpaPpa_degrees(3) = P_plunge_degrees model_TpaPpa_degrees(4) = P_azimuth_degrees END SELECT ! on j END SELECT ! on k; boundary class minimum_angle_degrees = DCROT ( actual_TpaPpa_degrees, model_TpaPpa_degrees ) E = (3 / (2 * Phi_max_degrees * radians_per_degree)) * & & MAX( 0.0, (1.0 - (minimum_angle_degrees / Phi_max_degrees)**2) ) CASE (2) ! Pacheco & Sykes [1992]: SELECT CASE(k) ! boundary class CASE (1, 6) ! CCB or OCB: SELECT CASE (j) CASE (1) ! thrust with shortening parallel to relative plate velocity: IF (thrust) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (2) ! thrust with partitioned shortening perpendicular to plate boundary: IF (thrust) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (3) ! partitioned strike-slip on plane striking parallel to plate boundary: IF (strikeslip) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF END SELECT ! on j CASE (2, 5) ! CTF or OTF: SELECT CASE (j) CASE (1) ! strike-slip on vertical plane of plate boundary: IF (strikeslip) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (2) ! normal on plane striking parallel to plate boundary: IF (normal) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (3) ! thrust on plane striking parallel to plate boundary: IF (thrust) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF END SELECT ! on j CASE (3, 4) ! CRB or OSR: SELECT CASE (j) CASE (1) ! normal faulting, with extension in direction of relative plate velocity: IF (normal) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (2) ! normal faulting, with partitioned extension perpendicular to plate boundary: IF (normal) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (3) ! partitioned strike-slip on plane striking parallel to plate boundary: IF (strikeslip) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF END SELECT ! on j CASE (7) ! SUB: SELECT CASE (j) CASE (1) ! main inter-plate thrust fault: IF (thrust) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (2) ! strike-slip parallel to arc: IF (strikeslip) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (3) ! down-dip extension due to bending and/or stress-guide: IF (normal) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF CASE (4) ! down-dip compression due to bending and/or stress-guide: IF (thrust) THEN ; E = 1.0 ; ELSE ; E = 0.0 ; END IF END SELECT ! on j END SELECT ! on k; boundary class END SELECT ! on catalog_ID (CMT, Pacheco & Sykes, ...) ELSE ! no mechanism information; punt: E = 1.0 END IF ! end of E (focal mechanism) section IF (E > 0.0) THEN ! this EQ mechanism could match this model sum_CDE = sum_CDE + C * D * E END IF ! E > 0.0; this EQ mechanism could match this model END IF ! D > 0.0; this EQ depth could match this model END IF ! this model has a chance of happening END DO ! j = 1, 4 ! (up to) 4 model earthquakes for this step !------------------------------------------------------------------------------------ P_prime = A * B * sum_CDE sum_P_prime_in_class(k) = sum_P_prime_in_class(k) + P_prime sum_P_prime_all_classes = sum_P_prime_all_classes + P_prime !Check whether a new overall maximum has been achieved: IF (P_prime > highest_P_prime) THEN highest_P_prime = P_prime best_i = i best_k = k END IF !Check whether a new maximum has been achieved for this class: IF (P_prime > highest_P_prime_in_class(k)) THEN highest_P_prime_in_class(k) = P_prime best_i_in_class(k) = i END IF END IF ! (A * B) > 0.0 (meaning this step has possibilities) END DO ! i = 1, nSteps; testing all steps as possible producers of this EQ IF (Monte_Carlo_method.AND.(best_k > 0)) THEN !Re-determine best_k, and therefore best_i (which often will turn out to be unchanged), !by comparing a random number to a set of adjacent bins whose widths are proportional !to sum_P_prime_in_class(1:7). !Set up the bins for classes which are in competition: partial_sum = 0.0 ! initialization before sum upper_fraction = 0.0 ! whole array; to simplify debugging DO k = 1, 7 in_competition(k) = (sum_P_prime_in_class(k) > 0.0) IF (in_competition(k)) THEN partial_sum = partial_sum + sum_P_prime_in_class(k) / sum_P_prime_all_classes upper_fraction(k) = partial_sum END IF END DO CALL RNUN(1, random) ! Retrieve next random number, in range (0, 1). !Compare random number to upper limit of each active bin. !Note: Order of comparison eliminates the need to test the lower limit. comparing: DO k = 1, 7 IF (in_competition(k)) THEN IF (random <= upper_fraction(k)) THEN best_k = k EXIT comparing END IF END IF END DO comparing best_i = best_i_in_class(best_k) END IF ! =========== end heart of this program ====================================================== ! At this point, critical information is contained in: ! pure_epicenter = .TRUE. if this earthquake was NOT in any orogen ! pure_epicenter_c1 = ' ' (not in any orogen) or '*' (orogen) ! best_k = 1-7 = index of class that this earthquake should be assigned to. ! best_i = 1-nStep = index of the plate boundary step this earthquake should be assigned to. ! sum_P_prime_in_class(1:7) = values of P_prime for steps, aggregated and summed by step class. ! sum_P_prime_all_classes = SUM(sum_P_prime_in_class) ! =========== end heart of this program ====================================================== ! output this event to appropriate catalog(s): SELECT CASE (best_k) CASE (0) ! INT count_all(0) = count_all(0) + 1 i0 = 0 ! place-holder, for name of associated step pure_step_c1 = ' ' ! because it is not associated with any step, step cannot make it impure relative_class_probability_int = 0 ! whole array distance_km = arc_km_from_step(1) DO i = 2, nSteps distance_km = MIN(distance_km, arc_km_from_step(i)) END DO WRITE (102, 3002) catalog_name, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data), & & pure_epicenter_c1, class_c3(0), pure_step_c1, & & i0, & ! additional information: step number is 0 & (relative_class_probability_int(n), n = 1, 7), & ! additional info: 7 0's (probabilities) & distance_km ! additional info: distance to closest step IF (pure_epicenter) THEN count_pure(0) = count_pure(0) + 1 WRITE (101, 3002) catalog_name, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data), & & pure_epicenter_c1, class_c3(0), pure_step_c1, & & i0, & ! additional information: step number is 0 & (relative_class_probability_int(n), n = 1, 7), & ! additional info: 7 0's (probabilities) & distance_km ! additional info: distance to closest step END IF CASE (1:7) ! CCB:SUB count_all(best_k) = count_all(best_k) + 1 pure_step_c1 = star_c1(best_i) pure_step = (pure_step_c1 /= '*') relative_class_probability_int = NINT(100 * sum_P_prime_in_class / sum_P_prime_all_classes) ! array of 7 %'s distance_km = distance_km_from_step(best_i) WRITE (10 * best_k + 2, 3002) catalog_name, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data), & & pure_epicenter_c1, class(best_i), pure_step_c1, best_i, & ! additional info. & (relative_class_probability_int(n), n = 1, 7), & ! additional info. & distance_km ! additional info: distance to closest step !Note: Unless there has been some error, we expect: ! class(best_i) == class_c3(best_k) !Either one could be output. IF (pure_epicenter.AND.pure_step) THEN count_pure(best_k) = count_pure(best_k) + 1 WRITE (10 * best_k + 1, 3002) catalog_name, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data), & & pure_epicenter_c1, class(best_i), pure_step_c1, best_i, & ! additional info. & (relative_class_probability_int(n), n = 1, 7), & ! additional info. & distance_km ! additional info: distance to closest step !Note: Unless there has been some error, we expect: ! class(best_i) == class_c3(best_k) !Either one could be output. !addition of 2003.06: IF (eq_mag > class_threshhold_magnitude(best_k)) THEN !GPBhere eq_count_over_threshhold(best_i) = eq_count_over_threshhold(best_i) + 1 eq_moment_Nm = 10.0**(9.05 + 1.5 * eq_mag) ! Hanks & Kanamori [1979] moment_sum_over_threshhold(best_i) = moment_sum_over_threshhold(best_i) + eq_moment_Nm END IF END IF END SELECT END DO next_event CLOSE (4) CLOSE (11) CLOSE (12) CLOSE (21) CLOSE (22) CLOSE (31) CLOSE (32) CLOSE (41) CLOSE (42) CLOSE (51) CLOSE (52) CLOSE (61) CLOSE (62) CLOSE (71) CLOSE (72) CLOSE (101) CLOSE (102) !GPBhere !Addition of 2003.06: output the 7 "PB2002_XXX_pure.dat" files with eq counts and moment sums: DO k = 1, 7 ! classes 1-7 = CCB...SUB out_file_name = "PB2002_" // class_c3(k) // "_pure.dat" OPEN (UNIT = 201, FILE = out_file_name) ! unconditional OPEN; overwrites any existing file DO i = 1, nSteps IF (class(i) == class_c3(k)) THEN IF (star_c1(i) == ' ') THEN ! pure step WRITE (201, 9001) sequence_number(i), step_continuity_c1(i), c5(i), & & old_longitude(i), old_latitude(i), longitude(i), latitude(i), & & length_km(i), azimuth_integer(i), & & velocity_mmpa(i), velocity_azimuth_integer(i), & & divergence_mmpa(i), dextral_mmpa(i), & & elevation(i), age_Ma(i), & & class_continuity_c1(i), class(i), star_c1(i), & ! so far, same as READ & eq_count_over_threshhold(i), moment_sum_over_threshhold(i) 9001 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, & ! so far, same as 1001 used in READ & 1X, I6, 1X, E10.3) END IF ! a pure step END IF ! correct plate boundary class for this file END DO ! i = 1, nSteps CLOSE (201) END DO ! k = 1, 7 !------------- More output to screen ----------------------------------------------- WRITE (*, "(' Finished reading and classifying ', A, ':')") TRIM(catalog_file_name) IF (Monte_Carlo_method) THEN WRITE ( *, "(' using Monte-Carlo randomization of class assignments')") ELSE WRITE ( *, "(' using maximum-probability step assignment to decide class:')") END IF WRITE ( *, *) WRITE ( *, "(' ',8(4X,A3))") (class_c3(n), n = 0, 7) WRITE ( *, "(' pure:',8(1X,I6))") (count_pure(n), n = 0, 7) nEvents = SUM(count_pure) percentages = NINT(100.0 * count_pure / nEvents) WRITE ( *, "(' pure:',8(3X,I3,'%'))") (percentages(n), n = 0, 7) WRITE ( *, *) WRITE ( *, "(' all:',8(1X,I6))") (count_all(n), n = 0, 7) nEvents = SUM(count_all) percentages = NINT(100.0 * count_all / nEvents) WRITE ( *, "(' all:',8(3X,I3,'%'))") (percentages(n), n = 0, 7) WRITE ( *, *) !------------------------------------------------------------------------------------ !------------- Output to log file ------------------------------------------------- out_file_name = TRIM(prefix) // "_PB2002_EQ_classification.txt" OPEN (UNIT = 99, FILE = out_file_name) WRITE (99, "('Running EQ_classification_II.exe, from EQ_classification_II.f90....')") WRITE (99, *) WRITE (99,"('Run on ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & datetimenumber(1), datetimenumber(2), datetimenumber(3), & datetimenumber(5), datetimenumber(6), datetimenumber(7) WRITE (99, *) WRITE (99, "('Catalog file = ',A)") TRIM(catalog_file_name) WRITE (99, *) WRITE (99, "('Parameter values (compiled-in):')") WRITE (99, "('-------------------------------')") WRITE (99, "('All boundaries:')") WRITE (99, "(L10, 2X,'Monte_Carlo_method')") Monte_Carlo_method IF (Monte_Carlo_method) THEN WRITE (99, "(I10, 2X,'iseed')") iseed END IF WRITE (99, "(F10.1,2X,'alongStep_sigma_km')") alongStep_sigma_km WRITE (99, "(F10.1,2X,'earthquake_depth_sigma_km')") earthquake_depth_sigma_km WRITE (99, "(F10.1,2X,'Phi_max_degrees')") Phi_max_degrees WRITE (99, "('Non-SUB boundaries:')") WRITE (99, "(F10.1,2X,'class123456_cutoff_depth_km')") class123456_cutoff_depth_km WRITE (99, "(' - - - - - - - - - - - - - - - - -')") WRITE (99, "(I10, 2X,'CCB_pure_totalEvents')") CCB_pure_totalEvents WRITE (99, "(I10, 2X,'CTF_pure_totalEvents')") CTF_pure_totalEvents WRITE (99, "(I10, 2X,'CRB_pure_totalEvents')") CRB_pure_totalEvents WRITE (99, "(I10, 2X,'OSR_pure_totalEvents')") OSR_pure_totalEvents WRITE (99, "(I10, 2X,'OTF_pure_totalEvents')") OTF_pure_totalEvents WRITE (99, "(I10, 2X,'OCB_pure_totalEvents')") OCB_pure_totalEvents WRITE (99, "(I10, 2X,'SUB_pure_totalEvents')") SUB_pure_totalEvents WRITE (99, "(' - - - - - - - - - - - - - - - - -')") WRITE (99, "(F10.1,2X,'CCB_pure_totalLength_km')") CCB_pure_totalLength_km WRITE (99, "(F10.1,2X,'CTF_pure_totalLength_km')") CTF_pure_totalLength_km WRITE (99, "(F10.1,2X,'CRB_pure_totalLength_km')") CRB_pure_totalLength_km WRITE (99, "(F10.1,2X,'OSR_pure_totalLength_km')") OSR_pure_totalLength_km WRITE (99, "(F10.1,2X,'OTF_pure_totalLength_km')") OTF_pure_totalLength_km WRITE (99, "(F10.1,2X,'OCB_pure_totalLength_km')") OCB_pure_totalLength_km WRITE (99, "(F10.1,2X,'SUB_pure_totalLength_km')") SUB_pure_totalLength_km WRITE (99, "(' - - - - - - - - - - - - - - - - -')") WRITE (99, "(F10.1,2X,'CCB_pure_meanVelocity_mmpa')") CCB_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'CTF_pure_meanVelocity_mmpa')") CTF_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'CRB_pure_meanVelocity_mmpa')") CRB_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'OSR_pure_meanVelocity_mmpa')") OSR_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'OTF_pure_meanVelocity_mmpa')") OTF_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'OCB_pure_meanVelocity_mmpa')") OCB_pure_meanVelocity_mmpa WRITE (99, "(F10.1,2X,'SUB_pure_meanVelocity_mmpa')") SUB_pure_meanVelocity_mmpa WRITE (99, "(' - - - - - - - - - - - - - - - - -')") WRITE (99, "('non-SUB boundaries (symmetrical):')") WRITE (99, "(F10.1,2X,'CCB_halfwidth_km')") CCB_halfwidth_km WRITE (99, "(F10.1,2X,'CTF_halfwidth_km')") CTF_halfwidth_km WRITE (99, "(F10.1,2X,'CRB_halfwidth_km')") CRB_halfwidth_km WRITE (99, "(F10.1,2X,'OSR_halfwidth_km')") OSR_halfwidth_km WRITE (99, "(F10.1,2X,'OTF_halfwidth_km')") OTF_halfwidth_km WRITE (99, "(F10.1,2X,'OCB_halfwidth_km')") OCB_halfwidth_km WRITE (99, "(' - - - - - - - - - - - - - - - - -')") WRITE (99, "('SUB boundaries (asymmetric):')") WRITE (99, "(F10.1,2X,'SUB_landward_halfwidth_km')") SUB_landward_halfwidth_km WRITE (99, "(F10.1,2X,'SUB_seaward_halfwidth_km')") SUB_seaward_halfwidth_km WRITE (99, "(F10.1,2X,'SUB_peakSeismicity_km')") SUB_peakSeismicity_km WRITE (99, "('quadratic model of depth to slab in SUBs:')") WRITE (99, "(F10.2,2X,'trench_depth_km')") trench_depth_km WRITE (99, "(F10.4,2X,'slope_at_trench')") slope_at_trench WRITE (99, "(F10.1,2X,'slab_depth_under_arc_km')") slab_depth_under_arc_km WRITE (99, "(F10.1,2X,'arc_distance_km')") arc_distance_km WRITE (99, "(F10.2,2X,'interplate_thrust_depth_sigma_fraction')") interplate_thrust_depth_sigma_fraction WRITE (99, "('-------------------------------------------------------------------')") WRITE (99, *) IF (Monte_Carlo_method) THEN WRITE (99, "('Using Monte-Carlo randomization of class assignments:')") ELSE WRITE (99, "('Using maximum-probability step assignment to decide class:')") END IF WRITE (99, *) WRITE (99, "(' ',8(4X,A3))") (class_c3(n), n = 0, 7) WRITE (99, "(' pure:',8(1X,I6))") (count_pure(n), n = 0, 7) nEvents = SUM(count_pure) percentages = NINT(100.0 * count_pure / nEvents) WRITE (99, "(' pure:',8(3X,I3,'%'))") (percentages(n), n = 0, 7) WRITE (99, *) WRITE (99, "(' all:',8(1X,I6))") (count_all(n), n = 0, 7) nEvents = SUM(count_all) percentages = NINT(100.0 * count_all / nEvents) WRITE (99, "(' all:',8(3X,I3,'%'))") (percentages(n), n = 0, 7) !------------------------------------------------------------------------------------ CLOSE(99) CALL Pause() CONTAINS INTEGER FUNCTION iPlate(c2) IMPLICIT NONE CHARACTER*2, INTENT(IN) :: c2 INTEGER :: i iPlate = 0 DO i = 1, nPlates IF (ID(i) == c2) THEN iPlate = i RETURN END IF END DO END FUNCTION iPlate SUBROUTINE Oblique_FPS (fault_strike_degrees, fault_dip_degrees, dip_is_to_right, & & relative_velocity_azimuth_degrees, relative_velocity_is_divergent, & ! input & e1_azimuth_degrees, e1_plunge_degrees, & & e2_azimuth_degrees, e2_plunge_degrees, & & e3_azimuth_degrees, e3_plunge_degrees) ! output !Computes principal strain axes (e1 most compressive, e2 neutral, e3 most extensional) !of the double-couple moment tensor produced by oblique slip on a dipping fault plane. !Warning: Do NOT call this routine with fault_dip_degrees = 90.0! Dip should be in range: 0.0 < dip < 90.0 ! Also, do NOT call this routine for events exactly at the N (or S) pole! USE Sphere ! Sphere.f90 by Peter Bird, UCLA IMPLICIT NONE LOGICAL, INTENT(IN) :: dip_is_to_right, relative_velocity_is_divergent REAL, INTENT(IN) :: fault_strike_degrees, fault_dip_degrees, relative_velocity_azimuth_degrees INTEGER, INTENT(OUT) :: e1_azimuth_degrees, e1_plunge_degrees, & & e2_azimuth_degrees, e2_plunge_degrees, & & e3_azimuth_degrees, e3_plunge_degrees INTEGER :: t_azimuth_degrees, t_plunge_degrees REAL :: crossStrike_component, downDip_component, strike_component, vertical_component REAL, DIMENSION(3) :: backSlip_uvec, & & crossStrike_uvec, & & downDip_uvec, & & e1_uvec, e2_uvec, e3_uvec, eX_uvec, & & faultNormal_uvec, & & heave_uvec, & & phi_uvec, & & slip_uvec, strike_uvec, & & tvec, theta_uvec, & & up_uvec !Validate input values: IF ((fault_dip_degrees <= 0.0).OR.(fault_dip_degrees >= 90.0)) THEN WRITE (*, "(' ERROR: Oblique_FPS requires 0.0 < fault_dip_degrees < 90.0')") CALL Traceback() END IF !Note: Without loss of generality, earthquake is assumed to be at (0 E, 0 N): up_uvec = (/ 1.0, 0.0, 0.0 /) !Rest of code should be correct, regardless of earthquake location. !Finish local coordinate system: CALL Local_Theta(up_uvec, theta_uvec) ! South CALL Local_Phi (up_uvec, phi_uvec) ! East !Describe orientation of fault: strike_uvec(1:3) = -COS(fault_strike_degrees * radians_per_degree) * theta_uvec(1:3) + & & SIN(fault_strike_degrees * radians_per_degree) * phi_uvec(1:3) IF (dip_is_to_right) THEN crossStrike_uvec(1:3) = -COS((fault_strike_degrees + 90.0) * radians_per_degree) * theta_uvec(1:3) + & & SIN((fault_strike_degrees + 90.0) * radians_per_degree) * phi_uvec(1:3) ELSE crossStrike_uvec(1:3) = -COS((fault_strike_degrees - 90.0) * radians_per_degree) * theta_uvec(1:3) + & & SIN((fault_strike_degrees - 90.0) * radians_per_degree) * phi_uvec(1:3) END IF ! Note: crossStrike_uvec is horizontal, and lies above the dipping fault. downDip_uvec(1:3) = COS(fault_dip_degrees * radians_per_degree) * crossStrike_uvec(1:3) - & & SIN(fault_dip_degrees * radians_per_degree) * up_uvec(1:3) IF (dip_is_to_right) THEN CALL Cross(downDip_uvec, strike_uvec, faultNormal_uvec) ELSE CALL Cross(strike_uvec, downDip_uvec, faultNormal_uvec) END IF !Note: faultNormal_uvec always has a positive (upward) vertical component. !Describe slip vector: heave_uvec(1:3) = -COS(relative_velocity_azimuth_degrees * radians_per_degree) * theta_uvec(1:3) + & & SIN(relative_velocity_azimuth_degrees * radians_per_degree) * phi_uvec(1:3) !Note: At this point, heave_uvec has same sense as input relative_velocity_azimuth_degrees IF (relative_velocity_is_divergent) THEN IF (Dot(heave_uvec, crossStrike_uvec) < 0.0) heave_uvec = -heave_uvec ELSE ! relative velocity is convergent IF (Dot(heave_uvec, crossStrike_uvec) > 0.0) heave_uvec = -heave_uvec END IF !Note: At this point, heave_uvec refers to motion of the hanging wall, relative to the footwall. strike_component = Dot(heave_uvec, strike_uvec) ! may be negative crossStrike_component = Dot(heave_uvec, crossStrike_uvec) ! positive for extension; negative for thrusting vertical_component = -crossStrike_component * TAN(fault_dip_degrees * radians_per_degree) ! negative for extension; positive for thrusting !Note: vertical_component magnitude is relative to ||heave_uvec|| = 1.0 downDip_component = crossStrike_component / COS(fault_dip_degrees * radians_per_degree) ! positive for extension; negative for thrusting !Note: downDip_component magnitude is relative to ||heave_uvec|| = 1.0 tvec(1:3) = strike_component * strike_uvec(1:3) + downDip_component * downDip_uvec(1:3) CALL Make_Uvec(tvec, slip_uvec) ! refers to motion of hanging wall relative to footwall !Describe moment tensor by its principal axes, getting e2 first: CALL Cross(faultNormal_uvec, slip_uvec, e2_uvec) !Find a vector normal to e2 (composed from slip_uvec and faultNormal_uvec) which is nearly horizontal: IF (Dot(slip_uvec, up_uvec) > 0.0) THEN backSlip_uvec = -slip_uvec ELSE backSlip_uvec = slip_uvec END IF !Note: backSlip_uvec points (obliquely) down in the fault plane; it might be parallel, or opposite, to slip_uvec. tvec(1:3) = (backSlip_uvec(1:3) + faultNormal_uvec(1:3)) / 2.0 ! averaging one pointing down with one pointing up. CALL Make_Uvec(tvec, eX_uvec) !Note: eX_uvec is normal to e2_uvec, and 45 degrees from both slip_uvec and faultNormal_uvec; ! therefore it should be either e1_uvec or e3_uvec. ! Choose the name of the principal axis which is closest to horizontal: IF (relative_velocity_is_divergent) THEN e3_uvec = eX_uvec CALL Cross (e3_uvec, e2_uvec, e1_uvec) ELSE ! relative velocity is convergent; a thrust e1_uvec = eX_uvec CALL Cross (e1_uvec, e2_uvec, e3_uvec) END IF !Be sure all principal axes have positive plunge, then convert to azimuth and plunge: IF (Dot(e1_uvec, up_uvec) > 0.0) e1_uvec = -e1_uvec IF (Dot(e2_uvec, up_uvec) > 0.0) e2_uvec = -e2_uvec IF (Dot(e3_uvec, up_uvec) > 0.0) e3_uvec = -e3_uvec CALL Uvec_2_PlungeAzimuth(up_uvec, e1_uvec, e1_plunge_degrees, e1_azimuth_degrees) CALL Uvec_2_PlungeAzimuth(up_uvec, e2_uvec, e2_plunge_degrees, e2_azimuth_degrees) CALL Uvec_2_PlungeAzimuth(up_uvec, e3_uvec, e3_plunge_degrees, e3_azimuth_degrees) END SUBROUTINE Oblique_FPS SUBROUTINE Pause () IMPLICIT NONE WRITE (*, "(/' Press [Enter] to continue...')") READ (*, *) RETURN END SUBROUTINE Pause END PROGRAM EQ_classification_II