PROGRAM Delete_Cracked_Benchmarks ! A short utility program that reads a .GPS file ! with geodetic benchmark locations, velocities, and error ellipses !(in Peter Bird's .GPS format, compatible with NeoKinema and NeoKineMap), ! and also reads a (.DIG) file of digitized traces of active faults, ! and may also read a file of geologic slip rates for those faults, ! so as to classifiy benchmarks as CRACKED if they are too close ! to the trace of any fast-slipping fault. ! Extended in 2008 to also decimate and renumber any associated .GP2 file. ! By Peter Bird, UCLA, 2003.09.23, 2006.11.03, & 2008.01.22 USE Sphere ! Fortran 90 MODULE of spherical-trigonometry code by P. Bird IMPLICIT NONE REAL, PARAMETER :: R = 6371000.0 ! radius of Earth, for converting degrees to m and km REAL, PARAMETER :: cot_normal_dip = 0.7002 ! 1./TAN(55.); assumes 10-degree mean flattening of normal faults ! by finite-strain after initiation at 65 degrees dip. REAL, PARAMETER :: cot_thrust_dip = 2.14451 ! 1./TAN(25.) CHARACTER*1 :: c1 CHARACTER*6 :: c6 CHARACTER*50 :: c50 CHARACTER*80 :: f_dat, f_dig, in_gps_file_name, good_out_gps_file_name, bad_out_gps_file_name, & & in_gp2_file_name, good_out_gp2_file_name CHARACTER*132 :: f_dat_format, f_dat_titles, title_line, gps_format, header_line CHARACTER*30, DIMENSION(:), ALLOCATABLE :: frame CHARACTER*50, DIMENSION(:), ALLOCATABLE :: site__source__sequence INTEGER :: f_dig_count, f_highest, f_traces, & & i, i_plus_3, ios, i1, & & j, j1, j2, k, kp, last_DOF_index, & & L_method, lines, number_cracked, safe_benchmarks INTEGER, DIMENSION(:), ALLOCATABLE :: new_DOF_index INTEGER, DIMENSION(:,:),ALLOCATABLE :: trace_loc ! gives locations where each trace is found, ! in terms of digitised points (in traces); ! (1:2 = first, last entry in traces; ! 0:f_highest = trace index (ties f_dat to f_dig)) ! Note: If trace_loc(1, i) == 0 then no fault with this index number was read in. LOGICAL :: beside, fast_enough, got_index, got_point, in_trace, is_there_a_gp2 LOGICAL, DIMENSION(:), ALLOCATABLE :: cracked REAL :: a, azimuth_radians_1, azimuth_radians_2, & & b, bench_azimuth_radians, & & covariance, critical_heave_rate_mmpa, & & distance_km, distance_radians, & & factor, & & lon, & & t1, t2, this_heave_rate_mmpa REAL, DIMENSION(3) :: pole_uvec, tv, tvec, tvo, tv1, tv2, uvec, uvec1, uvec2 REAL, DIMENSION(:), ALLOCATABLE :: E_lon_deg, heave_rate_mmpa, N_lat_deg, v_E_mmpa, v_N_mmpa, & & v_E_sigma, v_N_sigma, correlation REAL, DIMENSION(:,:),ALLOCATABLE :: benchmark_uvec, trace ! all digitized points on all fault traces in one long list; access through trace_loc. ! (1:3 = x,y,z components of unit vector; 1:f_dig_count = in order read) WRITE (*, "(' PROGRAM Delete_Cracked_Benchmarks')") WRITE (*, "(' A short utility program that reads a .GPS file')") WRITE (*, "(' with geodetic benchmark locations, velocities, and error ellipses')") WRITE (*, "(' (in Peter Bird''s .GPS format, compatible with NeoKinema and NeoKineMap),')") WRITE (*, "(' and also reads a (.DIG) file of digitized traces of active faults,')") WRITE (*, "(' and may also read a file of geologic slip rates for those faults,')") WRITE (*, "(' so as to classify benchmarks as CRACKED if they are too close')") WRITE (*, "(' to the trace of any fast-slipping fault.')") WRITE (*, "(' Also decimates and renumbers any associated .GP2 file.')") WRITE (*, "(' By Peter Bird, UCLA, 2003.09.25, 2006.11.03, & 2008.01.22')") WRITE (*, "(' ')") 10 WRITE (*, "(' Enter name of input .GPS file: ')") READ (*, "(A)") in_gps_file_name OPEN (UNIT = 1, FILE = in_gps_file_name, STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' File not found (in current directory).')") CALL Pause() GOTO 10 END IF READ (1, "(A)") title_line READ (1, "(A)") gps_format READ (1, "(A)") header_line !count datum lines in input file lines = 0 counting: DO READ (1, gps_format, IOSTAT = ios) lon IF (ios /= 0) EXIT counting lines = lines + 1 END DO counting CLOSE (1) WRITE (*, "(' ',I8,' datum lines in input file.')") lines WRITE (*, "(' Is there an associated .GP2 file with covariance matrix? (T/F): '\)") READ (*, *) is_there_a_gp2 IF (is_there_a_gp2) THEN WRITE (*, "(' Enter name of input .gp2 file: ')") READ (*, "(A)") in_gp2_file_name END IF !allocate storage arrays for benchmarks and their velocities ALLOCATE ( E_lon_deg(lines) ) ALLOCATE ( N_lat_deg(lines) ) ALLOCATE ( v_E_mmpa(lines) ) ALLOCATE ( v_N_mmpa(lines) ) ALLOCATE ( v_E_sigma(lines) ) ALLOCATE ( v_N_sigma(lines) ) ALLOCATE ( correlation(lines) ) ALLOCATE ( frame(lines) ) ALLOCATE ( site__source__sequence(lines) ) ALLOCATE ( benchmark_uvec(3, lines) ) ALLOCATE ( cracked(lines) ) IF (is_there_a_gp2) THEN ALLOCATE ( new_DOF_index(2 * lines) ) DO i = 1, (2 * lines) new_DOF_index(i) = i ! initializing as "old_DOF_index" END DO END IF cracked = .FALSE. ! to assist debugging !----------------------------------------------------------------------------- !re-read input .gps file and memorize contents: OPEN (UNIT = 1, FILE = in_gps_file_name, STATUS = "OLD", PAD = "YES", IOSTAT = ios) READ (1, "(A)") title_line READ (1, "(A)") gps_format READ (1, "(A)") header_line DO i = 1, lines READ (1, gps_format, IOSTAT = ios) E_lon_deg(i), N_lat_deg(i), & & v_E_mmpa(i), v_N_mmpa(i), & & v_E_sigma(i), v_N_sigma(i), & & correlation(i), & & frame(i), & & site__source__sequence(i) IF (ios /= 0) THEN i_plus_3 = i + 3 WRITE (*, "(' Error during reading of datum number ',I8,' (file line number ',I8,')')") i, i_plus_3 CALL Pause() STOP END IF CALL LonLat_2_Uvec(E_lon_deg(i), N_lat_deg(i), tv) benchmark_uvec(1:3, i) = tv(1:3) END DO ! i = 1, lines CLOSE (1) !----------------------------------------------------------------------------- !read f.dig file (of active fault traces) 20 WRITE (*, "(' ')") WRITE (*, "(' Enter name of input file with digitized (.DIG) fault traces: ')") READ (*, "(A)") f_dig OPEN (UNIT = 2, FILE = f_dig, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) ! Skim file and count number of data points f_dig_count = 0 ! just initializing f_traces = 0 f_highest = 0 loop_thru: DO READ (2, "(A)", IOSTAT = ios) c50 IF (ios == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN f_dig_count = f_dig_count + 1 ELSE IF (c50(2:3) == "**") THEN ! most likely, '*** end of line segment ***' ELSE ! title line for the trace; either F0123N format, or AU\SU format f_traces = f_traces + 1 READ (c50,"(1X,I4)", IOSTAT = ios) i IF (ios == 0) THEN f_highest = MAX (f_highest, i) ELSE f_highest = MAX (f_highest, f_traces) END IF END IF ELSE; EXIT loop_thru; END IF END DO loop_thru CLOSE (UNIT = 2) ! (will be re-read just below) ! allocate arrays ALLOCATE ( trace (3, f_dig_count) ) trace = 0. ! whole array ! ALLOCATE ( trace_loc (2, 0:f_highest) ) trace_loc = 0 ! whole array; many entries will be filled in below ! fill arrays OPEN (UNIT = 2, FILE = f_dig, STATUS = "OLD", ACTION = "READ", PAD = "YES") in_trace = .FALSE. i = 0 f_traces = 0 read_dig: DO READ (2, "(A)", IOSTAT = ios) c50 IF (ios == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN got_point = .TRUE. got_index = .FALSE. READ (c50,*) t1, t2 ! E longitude, N latitude ELSE IF (c50(2:3) == "**") THEN ! most likely, '*** end of line segment ***' got_point = .FALSE. got_index = .FALSE. ELSE ! title line; may be in F0123N format, or AU\SU format! got_point = .FALSE. f_traces = f_traces + 1 READ (c50,"(1X,I4)", IOSTAT = ios) j2 ! new trace number (NOTE: Any slip-sense bytes are ignored.) IF (ios /= 0) j2 = f_traces ! consecutive numbering applied to faults which lack index numbers. got_index = .TRUE. END IF ELSE; EXIT read_dig; END IF IF (in_trace) THEN IF (got_point) THEN i = i + 1 CALL LonLat_2_Uvec(t1, t2, tvo) trace(1:3,i) = tvo trace_loc(2,j1) = i ELSE ! *** end ... in_trace = .FALSE. ENDIF ELSE ! (not in_trace) IF (got_index) THEN j1 = j2 ! new index becomes current index ELSE IF (got_point) THEN i = i + 1 CALL LonLat_2_Uvec(t1, t2, tv) trace(1:3,i) = tv trace_loc(1,j1) = i in_trace = .TRUE. ELSE ! *** end ... END IF END IF ! (in_trace) END DO read_dig CLOSE (UNIT =2) ! close f_dig for the last time WRITE (*, "(' Read ',I8,' fault traces'/' with ',I8,' total digitized points'/ & & ' and highest trace index of ',I8)") f_traces, f_dig_count, f_highest !----------------------------------------------------------------------------- 25 WRITE (*, "(/' Enter integer to select method of fault selection(?):'& &/' 0 means to use all digitized fault traces;' & &/' 1 means to read an f.dat file and only use traces'& &/' that exceed some minimum offset rate (e.g., 1 mm/a).')") READ (*, *) L_method IF ((L_method < 0).OR.(L_method > 1)) THEN WRITE (*, "(' ERROR: Illegal integer value entered. Try again.')") GO TO 25 END IF !----------------------------------------------------------------------------- IF (L_method == 1) THEN ALLOCATE ( heave_rate_mmpa(f_highest) ) heave_rate_mmpa = 0.0 ! whole array; in case no datum is read ! read f.dat and remember ONLY the highest heave-rate for each fault trace: 30 WRITE (*, "(' ')") WRITE (*, "(' Enter name of (.NKI) file with fault slip rates: ')") READ (*, "(A)") f_dat OPEN (UNIT = 3, FILE = f_dat, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) READ (3, "(A)") f_dat_format READ (3, "(A)") f_dat_titles read_f_dat: DO READ (3, f_dat_format, IOSTAT = ios) c6, c50, t1 IF (ios /= 0) EXIT read_f_dat READ (c6, "(1X,I4,1X)") i1 IF (i1 > f_highest) THEN WRITE (*,"('Illegally high trace index: ',I6)") i1 CALL Pause() STOP END IF c1 = c6(6:6) IF (c1 == 't') c1 = 'T' IF (c1 == 'p') c1 = 'P' IF (c1 == 'n') c1 = 'N' IF (c1 == 'd') c1 = 'D' IF (c1 == 'r') c1 = 'R' IF (c1 == 'l') c1 = 'L' IF (c1 == 's') c1 = 'S' IF (.NOT.((c1 == 'T') .OR. (c1 == 'N') .OR. (c1 == 'R') .OR. (c1 == 'L') & & .OR. (c1 == 'D') .OR. (c1 == 'P').OR. (c1 == 'S'))) THEN WRITE (*,"('Illegal slip sense: ',A1)") c1 CALL Pause() STOP END IF IF (c1 == 'T') THEN factor = cot_thrust_dip ELSE IF (c1 == 'P') THEN factor = 1. ELSE IF (c1 == 'S') THEN factor = 1. ELSE IF (c1 == 'N') THEN factor = cot_normal_dip ELSE IF (c1 == 'D') THEN factor = 1. ELSE IF (c1 == 'R') THEN factor = 1. ELSE IF (c1 == 'L') THEN factor = 1. END IF this_heave_rate_mmpa = factor * t1 heave_rate_mmpa(i1) = MAX(heave_rate_mmpa(i1), this_heave_rate_mmpa) END DO read_f_dat CLOSE (UNIT = 3) ! close f_dat !----------------------------------------------------------------------------- !ask user to decide what is a "fast-slipping" fault? WRITE (*, "(' ')") WRITE (*, "(' Enter critical heave rate (in mm/a) for faults to be ""fast-slipping"": ')") READ (*, *) critical_heave_rate_mmpa !----------------------------------------------------------------------------- END IF ! L_method == 1 !ask user to decide what is "too close" to a fault? WRITE (*, "(' ')") WRITE (*, "(' Enter minimum acceptable distance (in km) between benchmark and fault(s): ')") READ (*, *) distance_km distance_radians = (distance_km * 1000.0) / R !======================== end of Input, begin Computation ======================= !check each benchmark; if too close to fault, mark cracked(i) = .TRUE. and increment number_cracked DO j = 1, f_highest ! all potential trace numbers IF ((trace_loc(1,j) > 0).AND.(trace_loc(2,j) > trace_loc(1,j))) THEN IF (L_method == 0) THEN fast_enough = .TRUE. ELSE IF (L_method == 1) THEN fast_enough = (heave_rate_mmpa(j) >= critical_heave_rate_mmpa) END IF IF (fast_enough) THEN DO k = trace_loc(1,j), trace_loc(2,j)-1 kp = k + 1 uvec1(1:3) = trace(1:3, k) ! beginning uvec for segment uvec2(1:3) = trace(1:3, kp) ! ending uvec for segment IF ((uvec1(1) /= uvec2(1)).OR.(uvec1(2) /= uvec2(2)).OR.(uvec1(3) /= uvec2(3))) THEN ! non-trivial segment DO i = 1, lines ! that is, for all benchmarks IF (.NOT.cracked(i)) THEN ! (no point in cracking a benchmark twice!) uvec(1:3) = benchmark_uvec(1:3, i) !check for distance from beginning uvec: a = Arc(uvec, uvec1) IF (a <= distance_radians) THEN cracked(i) = .TRUE. ! <=============== ELSE ! make other tests IF (kp == trace_loc(2,j)) THEN ! last point on trace !check distance from end uvec a = Arc(uvec, uvec2) IF (a <= distance_radians) THEN cracked(i) = .TRUE. ! <=============== END IF END IF ! this is last segment in the trace !test for sideways distance IFF point is "beside" this segment: CALL Cross (uvec1, uvec2, tvec) CALL Make_uvec(tvec, pole_uvec) azimuth_radians_1 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec1) azimuth_radians_2 = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec2) !normally az1 will be greater than az2 IF (azimuth_radians_2 > azimuth_radians_1) THEN ! fix cycle shift to MAKE az1 > az2: azimuth_radians_2 = azimuth_radians_2 - Two_Pi END IF bench_azimuth_radians = Relative_Compass(from_uvec = pole_uvec, to_uvec = uvec) beside = (bench_azimuth_radians < azimuth_radians_1).AND. & & (bench_azimuth_radians > azimuth_radians_2) IF (beside) THEN a = Arc(pole_uvec, uvec) ! always positive; is it close enough to Pi/2 ??? b = ABS(Pi_over_2 - a) ! distance from segment, in radians IF (b <= distance_radians) THEN cracked(i) = .TRUE. ! <=============== END IF END IF ! beside END IF END IF ! not previously cracked END DO ! i = 1, lines (benchmarks), checking for cracks END IF ! non-trivial segment END DO ! k = trace_loc(1, j), trace_loc(2,j)-1 [segments] END IF ! this trace is "fast-slipping" END IF ! trace is non-empty END DO ! j = 1, f_highest ! all potential trace numbers IF (is_there_a_gp2) THEN ! construct array "new_DOF_index" whose subscript is the "old_DOF_index". last_DOF_index = 0 DO i = 1, lines IF (cracked(i)) THEN new_DOF_index(2 * i - 1) = 0 new_DOF_index(2 * i) = 0 ELSE new_DOF_index(2 * i - 1) = last_DOF_index + 1 new_DOF_index(2 * i) = last_DOF_index + 2 last_DOF_index = last_DOF_index + 2 END IF END DO END IF !======================== end of Computation; begin Output ====================== number_cracked = 0 DO i = 1, lines IF (cracked(i)) number_cracked = number_cracked + 1 END DO safe_benchmarks = lines - number_cracked WRITE (*, "(' ')") WRITE (*, "(' ',I8,' benchmarks are safely distant from all fault traces.')") safe_benchmarks WRITE (*, "(' ',I8,' benchmarks are unsafe because they are too close to a fault.')") number_cracked IF (number_cracked > 0) THEN WRITE (*, "(' ')") WRITE (*, "(' Enter name for new output .gps file with safe benchmarks: ')") READ (*, "(A)") good_out_gps_file_name IF (is_there_a_gp2) THEN WRITE (*, "(' Enter name for new output .gp2 file with covariance of safe benchmarks: ')") READ (*, "(A)") good_out_gp2_file_name END IF OPEN (UNIT = 4, FILE = good_out_gps_file_name) ! absolute unconditional OPEN; overwrites any existing file !write modified gps-like file with two new leading columns: WRITE (4, "(A)") TRIM(title_line) WRITE (4, "(A)") TRIM(gps_format) WRITE (4, "(A)") TRIM(header_line) WRITE (*, "(' ')") WRITE (*, "(' Enter name for new output .gps file with unsafe, cracked benchmarks: ')") READ (*, "(A)") bad_out_gps_file_name OPEN (UNIT = 5, FILE = bad_out_gps_file_name) ! absolute unconditional OPEN; overwrites any existing file !write modified gps-like file with two new leading columns: WRITE (5, "(A)") TRIM(title_line) WRITE (5, "(A)") TRIM(gps_format) WRITE (5, "(A)") TRIM(header_line) DO i = 1, lines IF (cracked(i)) THEN WRITE (5, gps_format, IOSTAT = ios) E_lon_deg(i), N_lat_deg(i), & & v_E_mmpa(i), v_N_mmpa(i), & & v_E_sigma(i), v_N_sigma(i), & & correlation(i), & & frame(i), & & site__source__sequence(i) ELSE WRITE (4, gps_format, IOSTAT = ios) E_lon_deg(i), N_lat_deg(i), & & v_E_mmpa(i), v_N_mmpa(i), & & v_E_sigma(i), v_N_sigma(i), & & correlation(i), & & frame(i), & & site__source__sequence(i) END IF ! cracked(i), or not END DO ! i = 1, lines CLOSE (4) CLOSE (5) END IF ! number_cracked > 0 IF (is_there_a_gp2) THEN OPEN (UNIT = 11, FILE = TRIM(in_gp2_file_name), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not find or open ',A)") TRIM(in_gp2_file_name) CALL Pause() STOP END IF OPEN (UNIT = 12, FILE = TRIM(good_out_gp2_file_name), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not create new ',A)") TRIM(good_out_gp2_file_name) CALL Pause() STOP END IF run:DO READ (11, *, IOSTAT = ios) i, j, covariance IF (ios /= 0) EXIT run IF ((new_DOF_index(i) > 0).AND.(new_DOF_index(j) > 0)) THEN WRITE (12, "(I4,1X,I4,1X,ES15.8)") new_DOF_index(i), new_DOF_index(j), covariance END IF END DO run CLOSE (11, DISP = "KEEP") CLOSE (12, DISP = "KEEP") END IF ! is_there_a_gp2: if so; decimate and renumber it. WRITE (*, "(' ')") WRITE (*, "(' Job done.')") CALL Pause() CONTAINS SUBROUTINE Pause () IMPLICIT NONE WRITE (*, "(/' Press [Enter] to continue...')") READ (*, *) RETURN END SUBROUTINE Pause END PROGRAM Delete_Cracked_Benchmarks