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