PROGRAM Fault_Corridors ! The process of adding narrow fault corridors to .FEG (Finite Element Grid) ! files intended for use with either NeoKinema or Restore involves 4 steps: ! (1) Create a regular regional grid in OrbWin, using the Tool ! TileRegionalGrid, and then tidy-up the edges. ! (2) Compare elements to digitized fault traces (f_.DIG), ! and delete any elements that are cut by any fault(s). ! (3) Add nodes and elements to create narrow fault-corridor ! strips of elements around each fault trace. ! (4) Hand-edit in OrbWin to stitch these separate grid regions ! together into one topologically-correct whole. ! This utility program automates steps (2) and (3) above. ! During the final hand-editing step (4), you can also discover and fix ! any element-overlaps or other problems that this program created. ! ! By Peter Bird, Department of Earth, Planetary, and Space Sciences, ! University of California, Los Angeles, Nov. 2017-Jan. 2018, & May 2021. USE DSphere ! Spherical-geometry utilities, from P. Bird, in file DSphere.f90 IMPLICIT NONE TYPE :: is123 ! element & internal coordinates of any point on surface INTEGER :: element REAL*8, DIMENSION(3) :: s END TYPE is123 CHARACTER*1 :: c, c1 CHARACTER*4 :: c4 CHARACTER*5 :: memo ! temporary memo to set translation_method = 1 at end of trace CHARACTER*6 :: fault_c6 CHARACTER*50 :: c50, fault_label CHARACTER*80 :: c80, input_fDIG_filename, input_faultOffset_filename, & & faultOffset_FORMAT, faultOffset_header, input_FEG_filename, output_FEG_filename CHARACTER*132 :: new_title_line, title_line CHARACTER*1, DIMENSION(:), ALLOCATABLE :: sense ! T = thrust, N = normal, D = detachment, R = dextral, L = sinistral ! (1:f_rst_count = offset datum index) CHARACTER*1, DIMENSION(:), ALLOCATABLE :: trace_type ! records 1-character *predominant* fault sense from f_dig (NOT from f_dat.nki). !(1:f_highest = trace index) CHARACTER*50, DIMENSION(10) :: master_faults CHARACTER*50, DIMENSION(:), ALLOCATABLE :: fault_name ! (1:f_rst_count = offset datum index) CHARACTER*80, DIMENSION(:), ALLOCATABLE :: f_dig_faultName_lines ! memorized so that they can be copied to palinspastic f_.dig files. !(1:f_highest = trace index) CHARACTER*80, DIMENSION(:), ALLOCATABLE :: f_dig_faultData_lines ! may include useful information such as dip, or translation method !(1:f_highest = trace index) INTEGER :: closest_k, closest_trace, f_DIG_count, f_highest, & & i, i_end, ios, iProgram, j, j_end, j1, j2, j3, k, k1, k2, k_end, & & line, look_at_element, master_fault_count, & & n, n1, n11, n12, n1future, n1origin, n1past, n1target, n2, n21, n22, n2future, n2origin, n2past, n2target, n3, n_faults, n_free_nodes, n_steps, & & n_lines_read, n_offsets, n_to_delete, n_type_1, n_type_2, n_type_3, n_type_4, n_type_5, & & nFl, number, numEl, numEl_limit, numNod, numNod_limit, & & read_status, seg_count, temp_int_1, temp_int_2 INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: end_nodes !New nodes that bracket the uniform-width parts of fault corridors, !between the complexities that arise near each termination/connection. !(1:2 = left/right (per digitizing direction); 1:2 = beginning/end of trace; 1:f_highest) INTEGER*2, DIMENSION(:, :), ALLOCATABLE :: end_type !Values 1, 2, 3, 4 indicate topology of each end of the fault trace: !(1:2 = initial/final; 1:f_highest) INTEGER*2, DIMENSION(:, :, :), ALLOCATABLE :: link_end_list !(1:10 neighbors; 1:2 ends; 1:f_highest Fnnnn #) INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: link_to_list !(1:10 neighbors; 1:2 ends; 1:f_highest Fnnnn #) INTEGER, DIMENSION(:, :), ALLOCATABLE :: link_list_length !records number of (non-zero) Fnnnn #s in link_to_list !(1:2 ends; 1:f_highest Fnnnn #s) INTEGER, DIMENSION(:, :), ALLOCATABLE :: neighbor ! list of neighboring finite elements !(1:3 = side crossed; 1:numEl = element index l_) INTEGER, DIMENSION(:), ALLOCATABLE :: new_nodes INTEGER, DIMENSION(:, :), ALLOCATABLE :: nodes INTEGER, DIMENSION(:), ALLOCATABLE :: onramp_count !count of other faults merging with this one, in type-4 terminations !(1:f_highest Fnnnn #) INTEGER, DIMENSION(:, :, :, :), ALLOCATABLE :: onramp_nodes !nodes of small elements created at type-4 "onramp" fault mergers: !(1:2 = left/right side of trace (in dig. direction); ! 1:2 = earlier/later pair of nodes (in dig. direction); ! 1:10 = list of onramps (see onramp_count); ! 1:f_highest = Fnnnn#) INTEGER, DIMENSION(:,:), ALLOCATABLE :: seg_def ! defines a fault segment by its fault trace # and element # !(1:2 = trace, element; 1:seg_count = segment index) INTEGER, DIMENSION(:), ALLOCATABLE :: seg_u_ ! u_ = 1, 2, or 3 to identify isolated node of segment !(1:seg_count = segment index) INTEGER, DIMENSION(:,:), ALLOCATABLE :: trace_loc ! gives locations within "trace" where one fault trace is found, ! both in terms of digitised points (in traces) and in terms of ! fault segments (in seg_def, seg_end, seg_end_is); ! (1:4 = first, last entry in "trace"; first, last segment; ! 0:f_highest = trace index (ties f_rst to f_dig)) INTEGER*2, DIMENSION(:), ALLOCATABLE :: translation_method ! Values may include: ! 0: default method (trace moves with surrounding element, ! except as needed to prevent hooking ends at fault intersections); ! 1: symmetric_spreading_system method (mean of velocities of two bounding plates). !(1:f_dig_count = in order read) INTEGER, DIMENSION(:), ALLOCATABLE :: which_trace ! which trace index goes with this fault offset datum? ! use value to read actual trace location from trace_loc ! (1:f_rst_count = fault offset datum index) LOGICAL :: alongside, consider_that, consider_this, dextral, end1_OK, end2_OK, & & folding, found_mates, got_index, got_point, i_end_written, in_trace, & & left_side, left_turn, named_that, named_this, overlap, & & Restore, right_side, RT_corner, spreading_on_i, spreading_on_j LOGICAL*1, DIMENSION(:), ALLOCATABLE :: delete_this LOGICAL*1, DIMENSION(:), ALLOCATABLE :: f_2_in ! Does this fault have at least 2 points inside the current grid? !(In that case, it can be transported through space and time.) !(1:f_highest = trace index) LOGICAL*1, DIMENSION(:), ALLOCATABLE :: f_relevant ! Did this fault trace result in ANY segments, ! when first processed by Def_seg_v2? If not, mark it as useless, ! so that we don't waste more time on it. !(Most of these distant traces are composed of tedious '***HARD*** steps.') !(1:f_highest = trace index) LOGICAL*1, DIMENSION(:), ALLOCATABLE :: free_node LOGICAL*1, DIMENSION(:, :), ALLOCATABLE :: ends_done !Records whether end-nodes of a fault corridor have already been created and recorded? !(1:2 = start/end, 1:f_highest = Fnnnn fault trace #) REAL*8, PARAMETER :: min_corridor_width_km = 4.0D0 ! about 13% (1/8th) of the size of a level-8 element. REAL*8, PARAMETER :: max_corridor_width_km = 50.0D0 ! should allow at least 1 m.y. of finite strain before manual re-gridding in OrbWin REAL*8, PARAMETER :: R_Earth_km = 6371.0D0 ! (to be compared to the values above, allowing conversions to radians) REAL*8, PARAMETER :: expansion_factor = 10.0D0 ! allowance of extra memory, to let the ALLOCATED arrays contain new nodes and elements REAL*8, PARAMETER :: cot_thrust_dip = 2.14451D0 ! 1./TAN(25.); should agree with same statement in Restore REAL*8, PARAMETER :: cot_normal_dip = 0.46631D0 ! 1./TAN(65.); should agree with same statement in Restore REAL*8 :: allowance, arc_length, away_radians, azimuth, azimuth1, azimuth2, azimuth3, azimuth_ahead, azimuth_along_i, azimuth_along_j, & & azimuth_back, azimuth_into_i, azimuth_limit, azimuth_out, azimuth_to_loose_end, & & back_azimuth_along_j, closest_radians, corridor_length_radians, corridor_width_radians, distance, half_R2, & & large_corridor_width, max_corridor_width_radians, min_corridor_width_radians, & & offset_km, old_Elon, old_Nlat, R_Earth_m, radius, & & s_temp, small_corridor_width, spreading_km, t1, t12, t13, t2, t3, t31, t32, temp_corridor_length_radians, temp_R8, test REAL*8, DIMENSION(3) :: correction, end_uvec, far_uvec, inner_uvec, middle_uvec, & & omega_uvec, point1_uvec, point2_uvec, pole_uvec, pole_1_uvec, pole_2_uvec, & & tv, tvec, tvec1, tvec2, tvec3, tvo, uvec, uvec1, uvec2, uvec3, uvec4, uvec5, uvec6 REAL*8, DIMENSION(:), ALLOCATABLE :: a_ ! area of plane triangle element beneath surface, m**2 !(1:numEl = element index l_) REAL*8, DIMENSION(:), ALLOCATABLE :: azimuth_radians_list REAL*8, DIMENSION(:, :), ALLOCATABLE :: center ! Cartesian unit vectors at center of finite elements !(1:3 = x,y,z; 1:numEl = element index l_) REAL*8, DIMENSION(:), ALLOCATABLE :: node_Elon, node_Nlat REAL*8, DIMENSION(:, :), ALLOCATABLE :: node_uvec REAL*8, DIMENSION(:, :), ALLOCATABLE :: onramp_s !internal coordinate (within trace of fault) of an "on-ramp" !(type-4 merger of fault traces). !(1:10 = list of onramps for this fault; 1:f_highest = Fnnnn#) REAL*8, DIMENSION(:), ALLOCATABLE :: opening_heave_km, right_heave_km !(1:f_highest = trace index) for both vectors above. REAL*8, DIMENSION(:, :, :), ALLOCATABLE :: seg_end !Cartesian unit vectors at each end of fault segment !(a fault segment is the intersection of an element with a trace) !(1:3 = xyz, 1:2 = beginning/end, 1:seg_count = segment index) REAL*8, DIMENSION(:), ALLOCATABLE :: seg_eta_ ! eta_ (+1 or -1) for each fault segment, depending on whether ! isolated node u_ is to right or to left. !(1:seg_count = segment index) REAL*8, DIMENSION(:), ALLOCATABLE :: seg_kappa_ ! kappa_ (relative length) of each fault segment !(1:seg_count = segment index) REAL*8, DIMENSION(:, :), ALLOCATABLE :: trace ! uvecs of 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) TYPE(is123), DIMENSION(:, :), ALLOCATABLE :: seg_end_is ! element number and internal coordinates at each end of fault segment !(a fault segment is the intersection of an element with a trace) !(1:2 = beginning/end, 1:seg_count = segment index) TYPE(is123), DIMENSION(:), ALLOCATABLE :: trace_is ! locations of fault traces in internal coordinates !(1:f_dig_count = in order read) !------------------------------------------------------------------------------------ R_Earth_m = R_Earth_km * 1.0D3 half_R2 = (R_Earth_m**2) / 2.0D0 max_corridor_width_radians = max_corridor_width_km * 1.0D3 / R_Earth_m min_corridor_width_radians = min_corridor_width_km * 1.0D3 / R_Earth_m WRITE (*, "(' PROGRAM Fault_Corridors:')") WRITE (*, "(' ')") WRITE (*, "(' The process of adding narrow fault corridors to .FEG (Finite Element Grid)')") WRITE (*, "(' files intended for use with either NeoKinema or Restore involves 4 steps:')") WRITE (*, "(' (1) Create a regular regional grid in OrbWin, using the Tool')") WRITE (*, "(' TileRegionalGrid, and then tidy-up the edges.')") WRITE (*, "(' (2) Compare elements to digitized fault traces (f_.DIG),')") WRITE (*, "(' and delete any elements that are cut by any fault(s).')") WRITE (*, "(' (3) Add nodes and elements to create narrow fault-corridor')") WRITE (*, "(' strips of elements around each fault trace.')") WRITE (*, "(' (4) Hand-edit in OrbWin to stitch these separate grid regions')") WRITE (*, "(' together into one topologically-correct whole.')") WRITE (*, "(' This utility program automates steps (2) and (3) above.')") WRITE (*, "(' During the final hand-editing step (4), you can also discover and fix')") WRITE (*, "(' any element-overlaps or other problems that this program created.')") WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird, Department of Earth, Planetary, and Space Sciences,')") WRITE (*, "(' University of California, Los Angeles, Nov. 2017-Jan. 2018,')") WRITE (*, "(' and May 2021.')") WRITE (*, "(' ')") CALL Pause() !------------------------------------------------------------------------------------ !Read existing .FEG with regular grid: 10 WRITE (*, *) WRITE (*, "(' Enter file-name of an existing .FEG grid with regular elements:')") READ (*, "(A)") input_FEG_filename OPEN (UNIT = 1, FILE = TRIM(input_FEG_filename), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This file not found (in current folder). Try again...')") CALL Pause() GO TO 10 END IF READ (1, "(A)") title_line WRITE (*, *) WRITE (*, "(' Title line of this .FEG file =')") WRITE (*, "(' ', A)") TRIM(title_line(1:79)) READ (1, *) numNod numNod_limit = NINT(numNod * expansion_factor) ALLOCATE ( node_Elon(numNod_limit) ) node_Elon = 0.0D0 ! to simplify debugging ALLOCATE ( node_Nlat(numNod_limit) ) node_Nlat = 0.0D0 ! to simplify debugging ALLOCATE ( node_uvec(3, numNod_limit) ) node_uvec = 0.0D0 ! to simplify debugging ALLOCATE ( free_node(numNod_limit) ) free_node = .FALSE. DO i = 1, numNod READ (1, *, IOSTAT = ios) j, node_Elon(i), node_Nlat(i) IF (j /= i) THEN WRITE (*, "(' ERROR: When expecting node #', I6, ', found node #', I6)") i, j CALL Pause() STOP END IF IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not read: node# longitude latitude of node#', I6)") i CALL Pause() END IF CALL DLonLat_2_Uvec(node_Elon(i), node_Nlat(i), uvec) node_uvec(1:3, i) = uvec(1:3) END DO ! i = 1, numNod READ (1, *) numEl numEl_limit = NINT(numEl * expansion_factor) ALLOCATE ( nodes(3, numEl_limit) ) nodes = 0 ! whole array; to simplify debugging ALLOCATE ( center(3, numEl_limit) ) center = 0.0D0 ALLOCATE ( neighbor(3, numEl_limit) ) neighbor = 0 ALLOCATE ( a_(numEl_limit) ) a_ = 0.0D0 ALLOCATE ( delete_this(numEl_limit) ) delete_this = .FALSE. DO i = 1, numEl READ (1, *, IOSTAT = ios) j, nodes(1, i), nodes(2, i), nodes(3, i) IF (j /= i) THEN WRITE (*, "(' ERROR: When expecting element #', I6, ', found element #', I6)") i, j CALL Pause() STOP END IF IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not read: element# node1 node2 node3 of node#', I6)") i CALL Pause() END IF END DO ! i = 1, numEl WRITE (*, "(' ', I6, ' nodes and ', I6, ' elements were read.')") numNod, numEl CALL Pause() READ (1, *) nFl ! Should be == 0 IF (nFl > 0) CALL No_Fault_Elements_Allowed() CLOSE (1) ! input .FEG file CALL Plane_area (folding, look_at_element) IF (folding) THEN WRITE (*, *) WRITE (*, "(' ERROR in grid topology; element #', I6, ' has non-positive area.')") look_at_element CALL PAUSE STOP END IF ! folding !------------------------------------------------------------------------------------ !Read digitized fault traces (f_.DIG): [Note that method is similar to code in Restore_v4; !However, it is slightly improved in the sense that this version IGNORES any duplicate points.] 20 WRITE (*, *) WRITE (*, "(' Enter file-name of an existing f_.DIG fault-traces file:')") READ (*, "(A)") input_fDIG_filename OPEN (UNIT = 3, FILE = TRIM(input_fDIG_filename), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This file not found (in current folder). Try again...')") CALL Pause() GO TO 20 END IF WRITE (*, "(' Reading fault traces from ', A)") TRIM(input_fDIG_filename) line = 0 ! Skim file and count number of data points f_dig_count = 0 loop_thru: DO READ (3, "(A)", IOSTAT = read_status) c50; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN f_dig_count = f_dig_count + 1 ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN READ (c50,"(1X,I4)", IOSTAT = read_status) i IF (read_status == 0) THEN IF (i < 0) CALL Prevent ('negative fault number', line, input_fDIG_filename) f_highest = MAX (f_highest, i) END IF END IF ELSE; EXIT loop_thru; END IF END DO loop_thru CLOSE (UNIT = 3) ! (will be re-read soon) ! allocate arrays ALLOCATE ( trace (3, f_dig_count) ) trace = 0.0D0 ! whole array ! ALLOCATE ( trace_is(f_dig_count) ) ALLOCATE ( translation_method(f_dig_count) ) translation_method = 0 ! default method; some may be set to 1 based on per-fault data. ALLOCATE ( trace_loc (4, 0:f_highest) ) trace_loc = 0 ! whole array! ALLOCATE ( trace_type(f_highest) ) trace_type = ' ' ALLOCATE ( f_2_in(f_highest) ) f_2_in = .TRUE. ! unless negated later ALLOCATE ( f_relevant(f_highest) ) f_relevant = .TRUE. ! unless negated later, inside Def_seg_v2 ALLOCATE ( f_dig_faultName_lines(f_highest) ) f_dig_faultName_lines = ' ' ! whole array set blank/empty. ALLOCATE ( f_dig_faultData_lines(f_highest) ) f_dig_faultData_lines = ' ' ! whole array set blank/empty. ALLOCATE ( onramp_count(f_highest) ) onramp_count = 0 ! whole list ALLOCATE ( onramp_nodes(2, 2, 10, f_highest) ) onramp_nodes = 0 ! whole array ALLOCATE ( onramp_s(10, f_highest) ) onramp_s = 0.0D0 ! whole array ! fill arrays on this pass: OPEN (UNIT = 3, FILE = input_fDIG_filename, STATUS = "OLD", ACTION = "READ", PAD = "YES") line = 0 in_trace = .FALSE. i = 0 read_dig: DO READ (3, "(A)", IOSTAT = read_status) c80; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c80(1:2) == ' +') .OR. (c80(1:2) == ' -')) THEN got_point = .TRUE. got_index = .FALSE. READ (c80,*) t1, t2 ! E longitude, N latitude IF (ABS(t2) > 90.001D0) THEN WRITE (*, "(' Bad latitude ',F10.2,' in line ',I10,' of ',A)") t2, line, TRIM(input_fDIG_filename) CALL Pause() STOP END IF ELSE IF ((c80(1:1) == 'F') .OR. (c80(1:1) == 'f')) THEN got_point = .FALSE. got_index = .TRUE. READ (c80,"(1X,I4,A)") j2, c ! new trace number and (primary) sense f_dig_faultName_lines(j2) = c80 ! memorize header IF (c == 't') c = 'T' IF (c == 'p') c = 'P' IF (c == 'n') c = 'N' IF (c == 'd') c = 'D' IF (c == 'r') c = 'R' IF (c == 'l') c = 'L' IF (.NOT.((c == 'T') .OR. (c == 'N') .OR. (c == 'R') .OR. (c == 'L') & .OR. (c == 'D') .OR. (c == 'P'))) THEN WRITE (*, "(' Illegal slip sense: ', A1, ' for digitized trace ', I4)") c, j2 CALL Pause() STOP END IF trace_type(j2) = c memo = "none" ! unless a data line with "symmetric_spreading_system" appears later ELSE ! EITHER '*** end of line segment ***' OR 'dip_degrees 22', etc. got_point = .FALSE. got_index = .FALSE. IF (c80(1:3) /= "***") THEN f_dig_faultData_lines(j2) = c80 ! memorize !Interpret data, to see if any special action is needed: IF (INDEX(c80, "symmetric_spreading_system") > 0) THEN IF (VERIFY(c, "DdLlRr") > 0) THEN ! value of c is not one of these 6 allowed characters... WRITE (*, "(' ERROR: Key-phrase ""symmetric_spreading_system"" is only allowed with sense D, L, R.')") WRITE (*, "(' The following fault trace violates this rule:')") WRITE (*, "(' ',A)") TRIM(f_dig_faultName_lines(j2)(1:79)) WRITE (*, "(' ',A)") TRIM(f_dig_faultData_lines(j2)(1:79)) CALL Pause() STOP END IF IF (INDEX(c80, "first") > 0) THEN memo = "first" ELSE IF (INDEX(c80, "last") > 0) THEN memo = "last" ELSE IF (INDEX(c80, "both") > 0) THEN memo = "both" ELSE ! input_syntax problem! WRITE (*, "(' ERROR: Key-phrase ""symmetric_spreading_system"" was used incorrectly.')") WRITE (*, "(' In the header lines for the following fault trace:')") WRITE (*, "(' ',A)") TRIM(f_dig_faultName_lines(j2)(1:79)) WRITE (*, "(' ',A)") TRIM(f_dig_faultData_lines(j2)(1:79)) WRITE (*, "(' the keyword must be followed by one of these specifiers:')") WRITE (*, "(' first')") WRITE (*, "(' last')") WRITE (*, "(' both')") CALL Pause() STOP END IF ! input syntax problem ELSE memo = "none" END IF ! contains "symmetric_spreading_system", or NOT END IF ! not a *** line; therefore probably a faultData line END IF ! (x, y)-line, F0001R line, or something else ELSE; EXIT read_dig; END IF IF (in_trace) THEN IF (got_point) THEN IF ((t1 /= old_Elon).OR.(t2 /= old_Nlat)) THEN ! ignore any duplicate points i = i + 1 CALL DLonLat_2_Uvec(t1, t2, tvo) trace(1:3, i) = tvo(1:3) trace_loc(2, j1) = i ! continually overwriting; thus, "last point" address is always current old_Elon = t1 old_Nlat = t2 END IF ! this point is different from previous one ELSE ! "*** end of line segment ***" OR "dip_degrees 22", ... IF (c80(1:3) == "***") THEN ! finish up this trace in_trace = .FALSE. !Apply alternate translation method #1 to this trace, if "memo" commands: IF (memo == "first") THEN translation_method(trace_loc(1, j1)) = 1 memo = "none" ! done, and crossed-off ELSE IF (memo == "last") THEN translation_method(trace_loc(2, j1)) = 1 memo = "none" ! done, and crossed-off ELSE IF (memo == "both") THEN ! usually only 2 points, but let's be careful in case there are more: k1 = trace_loc(1, j1) k2 = trace_loc(2, j1) DO k = k1, k2 translation_method(k) = 1 END DO memo = "none" ! done, and crossed-off END IF ! N.B. If memo == "none", no action is required, because translation_method was initialized as 0. END IF 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 DLonLat_2_Uvec(t1, t2, tv) trace(1:3, i) = tv(1:3) trace_loc(1, j1) = i ! address of first point in a trace in_trace = .TRUE. old_Elon = t1 ! needed to test later points for novelty old_Nlat = t2 ELSE ! "*** end of line segment ***" OR "dip_degrees 22", ... END IF END IF ! (in_trace) END DO read_dig CLOSE (UNIT = 3) ! close f_dig WRITE (*, "(' ', I6, ' fault-trace points were read')") f_dig_count !------------------------------------------------------------------------------------ 30 WRITE (*, *) WRITE (*, "(' WHICH FINITE-ELEMENT PROGRAM shall we prepare a grid for?')") WRITE (*, "(' -----------------------------')") WRITE (*, "(' 1 = NeoKinema (neotectonics)')") WRITE (*, "(' 2 = Restore (paleotectonics)')") WRITE (*, "(' -----------------------------')") WRITE (*, "(' Please enter 1 or 2: '\)") READ (*, *) iProgram IF ((iProgram < 1).OR.(iProgram > 2)) THEN WRITE (*, "(' Please enter an integer from this table!')") CALL Pause() GO TO 30 END IF Restore = (iProgram == 2) IF (Restore) THEN ! additional input data will be needed (finite fault offsets): ALLOCATE ( opening_heave_km(f_highest) ) ALLOCATE ( right_heave_km(f_highest) ) opening_heave_km = 0.0D0 ! whole vector; to be incremented from data file right_heave_km = 0.0D0 ! ditto !Now, locate, OPEN, and READ the fault-offset data file to fill these two vectors: 40 WRITE (*, *) WRITE (*, "(' Enter filename of table of finite fault offsets--')") WRITE (*, "(' the same file that you will supply to Restore:')") READ (*, "(A)") input_faultOffset_filename OPEN (UNIT = 3, FILE = TRIM(input_faultOffset_filename), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This file was not found (in current folder).')") WRITE (*, "(' Please correct filename (or move file) and try again...')") CALL Pause() GO TO 40 END IF READ (3, "(A)") faultOffset_FORMAT READ (3, "(A)") faultOffset_header n_lines_read = 2 ! just initializing this count n_offsets = 0 ! just initializing this count learning: DO ! indefinite loop, until EOF READ (3, faultOffset_FORMAT, IOSTAT = ios) fault_c6, fault_label, offset_km IF (ios == -1) EXIT learning n_lines_read = n_lines_read + 1 ! useful for error messages? n_offsets = n_offsets + 1 c1 = fault_c6(1:1) IF (.NOT.((c1 == 'F').OR.(c1 == 'f'))) THEN WRITE (*, "(' ERROR: First bytes of line #',I6,' of file must be F or f.')") n_lines_read CALL Pause() STOP END IF c4 = fault_c6(2:5) READ (c4, *, IOSTAT = ios) i ! fault-trace # IF (ios /= 0) THEN WRITE (*, "(' ERROR: Bytes #2-5 of line #',I6,' were ',A,'. Integer needed!')") n_lines_read, c4 CALL Pause() STOP END IF c1 = fault_c6(6:6) ! offset-type indicator IF ((c1 == 'R').OR.(c1 == 'r')) THEN right_heave_km(i) = right_heave_km(i) + offset_km ELSE IF ((c1 == 'L').OR.(c1 == 'l')) THEN right_heave_km(i) = right_heave_km(i) - offset_km ELSE IF ((c1 == 'D').OR.(c1 == 'd')) THEN opening_heave_km(i) = opening_heave_km(i) + offset_km ELSE IF ((c1 == 'P').OR.(c1 == 'p')) THEN opening_heave_km(i) = opening_heave_km(i) - offset_km ELSE IF ((c1 == 'N').OR.(c1 == 'n')) THEN opening_heave_km(i) = opening_heave_km(i) + (offset_km * cot_normal_dip) ELSE IF ((c1 == 'T').OR.(c1 == 't')) THEN opening_heave_km(i) = opening_heave_km(i) - (offset_km * cot_thrust_dip) ELSE WRITE (*, "(' ERROR: F',I4,' in line #',I6,' of file had illegal offset letter:',A)") i, n_lines_read, c1 CALL Pause() STOP END IF END DO learning CLOSE (3) ! input_faultOffset_filename WRITE (*, "(' ',I6,' fault offsets were read, classified, and summed.')") n_offsets CALL Pause() END IF !------------------------------------------------------------------------------------ WRITE (*, *) WRITE (*, "(' Computing internal coordinates of fault-trace points ...')") CALL Find_s1s2s3() !------------------------------------------------------------------------------------ !Survey (count), and then re-survey and store fault segments. !Notes: Uses arrays "center", "neighbor", ! and "trace_is" defined in Find_s1s2s3. WRITE (*, *) WRITE (*, "(' Counting, and then recording, the fault segments ...')") CALL Def_seg_v2 (seg_count = seg_count, savem = .FALSE.) IF (seg_count > 0) THEN WRITE (*, "(' ', I6, ' fault segments were counted.')") seg_count WRITE (*, "(' (A fault segment is the intersection of one fault trace with one element.)')") ALLOCATE( seg_def(2, seg_count) ) ALLOCATE( seg_end(1:3, 2, seg_count) ) ALLOCATE( seg_end_is(2, seg_count) ) ALLOCATE( seg_eta_(seg_count) ) ALLOCATE( seg_kappa_(seg_count) ) ALLOCATE( seg_u_(seg_count) ) !OPEN (UNIT = 21, FILE = "fault_segments.dig") ! graphical-test output from Def_seg_v2 CALL Def_seg_v2 (seg_count = seg_count, savem = .TRUE.) !CLOSE (21) ELSE WRITE (*, "(' No fault segments were identified.')") CALL Pause() STOP END IF ! seg_count > 0 !------------------------------------------------------------------------------------ !Mark faulted elements for deletion: WRITE (*, *) WRITE (*, "(' Marking all faulted elements for deletion ...')") n_to_delete = 0 DO i = 1, seg_count delete_this(seg_def(2, i)) = .TRUE. END DO ! i = 1, seg_count DO i = 1, numEl IF (delete_this(i)) n_to_delete = n_to_delete + 1 END DO ! i = 1, numEl WRITE (*, "(' ', I6, ' of the regular elements are faulted and will be deleted.')") n_to_delete CALL Pause() !------------------------------------------------------------------------------------ !Actually delete any unwanted regular elements: WRITE (*, *) WRITE (*, "(' Now deleting all unwanted elements ...')") i = 1 ! which will SOMETIMES be incremented on passes through the loop below: making_space: DO IF (delete_this(i)) THEN !Shift contents of all arrays with a per-element subscript: IF (i < numEl) THEN DO j = i, numEl nodes(1:3, j) = nodes(1:3, j+1) center(1:3, j) = center(1:3, j+1) neighbor(1:3, j) = neighbor(1:3, j+1) a_(j) = a_(j+1) delete_this(j) = delete_this(j+1) !COMMENT: I am NOT updating trace_is here, because later in this code it will be convenient ! to test the logical quantity (trace_is()%element > 0) to check whether a given ! fault-trace point is within the area of the original, regular .FEG before any edits. END DO ! j = i, numEl END IF ! i < numEl !Finally, note that there is now one less element: numEl = numEl - 1 !Note that we do NOT increment i along this branch, because i now points to what was previously i+1, and this needs to be checked. ELSE ! no action, except to prepare to check the next element (or stop): i = i + 1 END IF IF (i > numEl) EXIT making_space END DO making_space WRITE (*, "(' ... Done.')") WRITE (*, *) WRITE (*, "(' There are now ', I6, ' regular elements remaining.')") numEl CALL Pause() !------------------------------------------------------------------------------------ WRITE (*, *) WRITE (*, "(' Now, checking for un-connected nodes ...')") free_node = .TRUE. ! whole list; unless negated in the loop below... DO i = 1, numEl free_node(nodes(1, i)) = .FALSE. free_node(nodes(2, i)) = .FALSE. free_node(nodes(3, i)) = .FALSE. END DO ! i = 1, numEl n_free_nodes = 0 DO i = 1, numNod IF (free_node(i)) n_free_nodes = n_free_nodes + 1 END DO ! i = 1, numNod WRITE (*, "(' ', I6, ' un-connected nodes can be deleted.')") n_free_nodes CALL Pause() !------------------------------------------------------------------------------------ !Actually delete these un-connected nodes: WRITE (*, *) WRITE (*, "(' Deleting all un-connected nodes ...')") i = 1 ! which will SOMETIMES be incremented in the loop below: clearing_sky: DO IF (free_node(i)) THEN !Shift contents of all arrays with a per-node subscript: DO j = i, numNod ! (even if i == numNod; just shift in a "clean, empty" array element) node_Elon(j) = node_Elon(j+1) node_Nlat(j) = node_Nlat(j+1) node_uvec(1:3, j) = node_uvec(1:3, j+1) free_node(j) = free_node(j+1) END DO ! j = i, numNod DO j = 1, numEl DO k = 1, 3 IF (nodes(k, j) > i) nodes(k, j) = nodes(k, j) - 1 END DO ! k = 1, 3 END DO ! j = 1, numEl !Finally, reduce the number of nodes by 1: numNod = numNod - 1 !Note that we do NOT increment i along this branch; the node that was formerly #i+1 is now #i, and must be checked. ELSE ! no action, except to prepare for next loop (or exit): i = i + 1 END IF IF (i > numNod) EXIT clearing_sky END DO clearing_sky WRITE (*, "(' ... Done.')") WRITE (*, *) WRITE (*, "(' There are now ', I6, ' connected nodes remaining.')") numNod !------------------------------------------------------------------------------------ ! ! Decide, in advance, what types of ends each fault has? Codes are: ! 1 = fault leaves domain of the .FEG at its edge; corridor cut-off where last segment ends. ! 2 = close-coupling to another (designated "master" or "symmetric_spreading_system") fault; ! 3 = tectonic knot (non-close duple junction, 3-J, or N-J); ! 4 = merging, where current fault ends anywhere on the trace of another fault. ! 5 = fault termination inside .FEG area; corridor will taper to one node. !(When debugging, note that all indexes start as 0, and can only be promoted from 0. ! That is, all decisions are final and not subject to reconsideration.) OPEN (UNIT = 29, FILE = "fault_connections.txt") WRITE (29, "('Analyzing fault topology to classify fault end-type:')") WRITE (*, *) WRITE (*, "(' Analyzing fault topology to classify fault end-type:')") ALLOCATE ( end_type(2, f_highest) ) ! INTEGER*2, value 0...5 (1:2 = initial end {as digitized & segmented}, final end; 1:f_highest) end_type = 0 ! which is a flag meaning "undefined". Unused Fnnnn's will never change. ALLOCATE ( link_to_list(10, 2, f_highest) ) ! up to 10 other faults may be involved at each end link_to_list = 0 ! initially... ALLOCATE ( link_end_list(10, 2, f_highest) ) ! INTEGER*2; value 1 or 2 (which end of other fault?) link_end_list = 0 ! initially... ALLOCATE ( link_list_length(2, f_highest) ) ! counts entries (0...10) in link_to_list. link_list_length = 0 ! initially... ALLOCATE ( azimuth_radians_list(0:10) ) ! used as temporary storage while constructing a duple 2-J, triple 3-J, or n-ple n-J of type-3. ALLOCATE ( new_nodes(0:10) ) ! (see above) !Scan for type-1 first, since these are topologically and tectonically the most-important: ! 1 = fault leaves domain of the .FEG at its edge; corridor cut-off where last segment ends. n_type_1 = 0 DO i = 1, f_highest IF (f_relevant(i)) THEN ! (Don't do anything about faults that have 0 segments.) j1 = trace_loc(1, i) j2 = trace_loc(2, i) IF (trace_is(j1)%element == 0) THEN end_type(1, i) = 1 n_type_1 = n_type_1 + 1 END IF IF (trace_is(j2)%element == 0) THEN end_type(2, i) = 1 n_type_1 = n_type_1 + 1 END IF END IF ! f_relevant(i) END DO ! i = 1, f_highest WRITE (*, *) WRITE (*, "(' ', I6, ' fault-ends extend outside the .FEG area:')") n_type_1 WRITE (29, *) WRITE (29, "(I6, ' fault-ends extend outside the .FEG area:')") n_type_1 CALL Pause() IF (n_type_1 > 0) THEN DO i = 1, f_highest IF ((end_type(1, i) == 1).AND.(end_type(2, i) == 1)) THEN WRITE (*, "(' Both ends of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('Both ends of: ', A)") TRIM(f_dig_faultName_lines(i)) ELSE IF (end_type(1, i) == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) ELSE IF (end_type(2, i) == 1) THEN WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(i)) END IF END DO ! i = 1, f_highest CALL Pause() END IF ! n_type_1 > 0 !Establish high-priority close-coupling (type-2) for pieces of master fault systems: ! 2 = close-coupling to another (designated "master" or "symmetric_spreading_system") fault; master_faults = ' ' ! whole list master_fault_count = 0 ! initially... WRITE (*, *) WRITE (*, "(' Faults with the ""symmetric_spreading_system"" attribute')") WRITE (*, "(' in a trace-data line, within the input f_.DIG file already read')") WRITE (*, "(' are always considered closely-coupled, due to their large offsets.')") WRITE (*, "(' You may now designate other fault names, such as ""San Andreas"",')") WRITE (*, "(' to receive the same treatment, if you wish.')") WRITE (29, *) WRITE (29, "('Faults with the ""symmetric_spreading_system"" attribute')") WRITE (29, "('in a trace-data line, within the input f_.DIG file already read')") WRITE (29, "('are always considered closely-coupled, due to their large offsets.')") WRITE (29, "('You may now designate other fault names, such as ""San Andreas"",')") WRITE (29, "('to receive the same treatment, if you wish.')") naming_names: DO i = 1, 10 WRITE (*, "(' Enter a short, simple fault name (or, ""end""): '\)") READ (*, "(A)") master_faults(i) IF ((master_faults(i) == "end").OR.(master_faults(i) == "end").OR.(master_faults(i) == "End")) THEN EXIT naming_names ELSE master_fault_count = master_fault_count + 1 WRITE (29, "('Enter a short, simple fault name (or, ""end""): ', A)") TRIM(master_faults(i)) END IF END DO naming_names ! i = 1, 10 WRITE (*, *) WRITE (*, "(' Below are listed any fault-pairs considered closely-coupled:')") WRITE (*, "(' ------------------------------------------------------------')") WRITE (29, *) WRITE (29, "('Below are listed any fault-pairs considered closely-coupled:')") WRITE (29, "('------------------------------------------------------------')") n_type_2 = 0 ! initializing, before search... DO i = 1, (f_highest - 1) named_this = .FALSE. ! unless... DO k = 1, master_fault_count IF (INDEX(f_dig_faultName_lines(i), TRIM(master_faults(k))) > 0) named_this = .TRUE. END DO ! k = 1, master_fault_count DO i_end = 1, 2 IF (f_relevant(i)) THEN ! It's OK to use the (non-zero) values in trace_loc... consider_this = (end_type(i_end, i) == 0).AND.(trace_is(trace_loc(i_end, i))%element > 0).AND. & & ((translation_method(trace_loc(i_end, i)) == 1).OR.(named_this)) !This fault-end not already assigned; this fault-end within the grid area; and !this fault is either a symmetric_spreading_system or master fault. ELSE consider_this = .FALSE. END IF IF (consider_this) THEN mating: DO j = (i + 1), f_highest named_that = .FALSE. ! unless... DO k = 1, master_fault_count IF (INDEX(f_dig_faultName_lines(j), TRIM(master_faults(k))) > 0) named_that = .TRUE. END DO ! k = 1, master_fault_count DO j_end = 1, 2 IF (f_relevant(j)) THEN ! It's OK to use the (non-zero) values in trace_loc.... consider_that = (end_type(j_end, j) == 0).AND.(trace_is(trace_loc(j_end, j))%element > 0).AND. & & ((translation_method(trace_loc(j_end, j)) == 1).OR.(named_that)) !That fault-end not already assigned; that fault-end within the grid area; and !that fault is either a symmetric_spreading_system or master fault. ELSE consider_that = .FALSE. END IF IF (consider_that) THEN ! must STILL pass test of geographic contiguity: uvec1(1:3) = trace(1:3, trace_loc(i_end, i)) uvec2(1:3) = trace(1:3, trace_loc(j_end, j)) distance = DArc(uvec1, uvec2) IF (distance <= min_corridor_width_radians) THEN ! GOT A MATCH! end_type(i_end, i) = 2 ! (previously determined to be undefined) link_list_length(i_end, i) = link_list_length(i_end, i) + 1 link_to_list(link_list_length(i_end, i), i_end, i) = j link_end_list(link_list_length(i_end, i), i_end, i) = j_end end_type(j_end, j) = 2 ! (previously determined to be undefined) link_list_length(j_end, j) = link_list_length(j_end, j) + 1 link_to_list(link_list_length(j_end, j), j_end, j) = i link_end_list(link_list_length(j_end, j), j_end, j) = i_end IF (i_end == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) ELSE WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(i)) END IF IF (j_end == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(j)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(j)) ELSE WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(j)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(j)) END IF WRITE (*, *) WRITE (29, *) n_type_2 = n_type_2 + 1 EXIT mating ! because only one spouse is allowed (on each end)! END IF ! GOT A MATCH! END IF ! consider_that END DO ! j_end = 1, 2 END DO mating ! j = (i + 1), f_highest END IF ! consider_this END DO ! i_end = 1, 2 END DO ! i = 1, (f_highest - 1) WRITE (*, "(' ------------------------------------------------------------')") WRITE (*, "(' ', I6, ' close connections were found.')") n_type_2 WRITE (29, "('------------------------------------------------------------')") WRITE (29, "(I6, ' close connections were found.')") n_type_2 CALL Pause() !Search for junctions of type_3: ! 3 = tectonic knot (angled duple-junction, 3-J, or N-J); WRITE (*, *) WRITE (*, "(' Below are listed any fault-groups involved in 2-Js, 3-Js or N-Js:')") WRITE (*, "(' -----------------------------------------------------------------')") WRITE (29, *) WRITE (29, "('Below are listed any fault-groups involved in 2-Js, 3-Js or N-Js:')") WRITE (29, "('-----------------------------------------------------------------')") n_type_3 = 0 ! initializing, before search... DO i = 1, (f_highest - 1) DO i_end = 1, 2 i_end_written = .FALSE. ! until this fault name gets printed in table (ONCE per end, max). found_mates = .FALSE. IF (f_relevant(i)) THEN ! It is safe to inquire about trace_is(trace_loc(i_end, i))%element, because trace_loc(i_end, i) > 0. consider_this = (trace_is(trace_loc(i_end, i))%element > 0).AND.(end_type(i_end, i) == 0) ! Once we get past this test, loops below can add more than one link! ELSE ! no way! consider_this = .FALSE. END IF IF (consider_this) THEN DO j = (i + 1), f_highest DO j_end = 1, 2 consider_that = f_relevant(j).AND.(end_type(j_end, j) == 0) IF (consider_that) THEN ! must STILL pass test of geographic contiguity: uvec1(1:3) = trace(1:3, trace_loc(i_end, i)) uvec2(1:3) = trace(1:3, trace_loc(j_end, j)) distance = DArc(uvec1, uvec2) IF (distance <= min_corridor_width_radians) THEN ! GOT A MATCH! end_type(i_end, i) = 3 ! (previously determined to be undefined) link_list_length(i_end, i) = link_list_length(i_end, i) + 1 link_to_list(link_list_length(i_end, i), i_end, i) = j link_end_list(link_list_length(i_end, i), i_end, i) = j_end end_type(j_end, j) = 3 ! (previously determined to be undefined) link_list_length(j_end, j) = link_list_length(j_end, j) + 1 link_to_list(link_list_length(j_end, j), j_end, j) = i link_end_list(link_list_length(j_end, j), j_end, j) = i_end IF (.NOT.i_end_written) THEN IF (i_end == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) ELSE WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(i)) END IF i_end_written = .TRUE. ! So, don't print it again... n_type_3 = n_type_3 + 1 ! Count this junction only once, even if it grows to an N-J. END IF ! .NOT.i_end_written IF (j_end == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(j)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(j)) ELSE WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(j)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(j)) END IF found_mates = .TRUE. END IF ! GOT A MATCH! END IF ! consider_that END DO ! j_end = 1, 2 END DO ! j = (i + 1), f_highest !----Beginning of special kludge: Reject any duple-junctions with very acute angle at intersection: --------------- IF (link_list_length(i_end, i) == 1) THEN ! duple-junction IF (i_end == 1) THEN uvec1(1:3) = seg_end(1:3, 1, trace_loc(3, i)) uvec2(1:3) = seg_end(1:3, 2, trace_loc(3, i)) azimuth1 = DCompass(uvec1, uvec2) ! azimuth of fault #i, in direction that is away from the junction ELSE ! i_end == 2 uvec1(1:3) = seg_end(1:3, 2, trace_loc(4, i)) uvec2(1:3) = seg_end(1:3, 1, trace_loc(4, i)) azimuth1 = DCompass(uvec1, uvec2) ! azimuth of fault #i, in direction that is away from the junction END IF j = link_to_list(1, i_end, i) j_end = link_end_list(1, i_end, i) IF (j_end == 1) THEN uvec1(1:3) = seg_end(1:3, 1, trace_loc(3, j)) uvec2(1:3) = seg_end(1:3, 2, trace_loc(3, j)) azimuth2 = DCompass(uvec1, uvec2) ! azimuth of fault #j, in direction that is away from the junction ELSE ! j_end == 2 uvec1(1:3) = seg_end(1:3, 2, trace_loc(4, j)) uvec2(1:3) = seg_end(1:3, 1, trace_loc(4, j)) azimuth2 = DCompass(uvec1, uvec2) ! azimuth of fault #j, in direction that is away from the junction END IF test = COS(azimuth2 - azimuth1) ! should NOT be too close to 1.0 ! IF (test > 0.707106D0) THEN ! reject this acute-angle "duple-junction": WRITE (*, "(' but this ""duple-junction"" was rejected for its very acute angle.')") WRITE (29, "('but this ""duple-junction"" was rejected for its very acute angle.')") n_type_3 = n_type_3 - 1 end_type(i_end, i) = 0 ! for potential later assignment to type 4 or 5... end_type(j_end, j) = 0 link_to_list(1, i_end, i) = 0 link_to_list(1, j_end, j) = 0 link_end_list(1, i_end, i) = 0 link_end_list(1, j_end, j) = 0 link_list_length(i_end, i) = 0 link_list_length(j_end, j) = 0 END IF ! reject this acute-angle "duple-junction" END IF ! duple-junction !----End of special kludge: Reject any duple-junctions with very acute angle at intersection: --------------- IF (found_mates) THEN WRITE (*, *) ! terminating a group of 2+ consecutive fault names WRITE (29, *) END IF END IF ! consider_this (end #i_end of fault #i) END DO ! i_end = 1, 2 END DO ! i = 1, (f_highest - 1) WRITE (*, "(' -------------------------------------------------------------------------')") WRITE (*, "(' ', I6, ' tectonic knots (angled duple-junctions, 3-Js, or N-Js) were found.')") n_type_3 WRITE (29, "('-------------------------------------------------------------------------')") WRITE (29, "(I6, ' tectonic knots (angled duple-junctions, 3-Js, or N-Js) were found.')") n_type_3 CALL Pause() !Classify remaining fault-ends as either mergers (type-4)) or truly "loose ends" (type-5): ! 4 = merging, where current fault ends anywhere on the trace of another fault. ! 5 = fault termination inside .FEG area; corridor will taper to one node. WRITE (*, *) n_type_4 = 0 ! initializing, before search... n_type_5 = 0 ! initializing, before search... allowance = 0.5D0 * min_corridor_width_radians WRITE (*, *) WRITE (*, "(' Below are listed all faults that end by merging:')") WRITE (*, "(' -----------------------------------------------------------')") WRITE (29, *) WRITE (29, "('Below are listed all faults that end by merging:')") WRITE (29, "('-----------------------------------------------------------')") DO i = 1, f_highest DO i_end = 1, 2 IF (f_relevant(i)) THEN ! It is safe to inquire about trace_is(trace_loc(i_end, i))%element, because trace_loc(i_end, i) > 0. consider_this = (trace_is(trace_loc(i_end, i))%element > 0).AND.(end_type(i_end, i) == 0) ELSE ! no way! consider_this = .FALSE. END IF IF (consider_this) THEN closest_radians = 9.99D0 ! initializing as ridiculously large... closest_trace = 0 closest_k = 0 !Search for mid-trace points on fault #j that are close to fault-end (i_end, i): uvec3(1:3) = trace(1:3, trace_loc(i_end, i)) ! test point DO j = 1, f_highest IF (j /= i) THEN j1 = trace_loc(1, j) j2 = trace_loc(2, j) consider_that = (j2 > j1).AND.f_relevant(j) IF (consider_that) THEN uvec1(1:3) = trace(1:3, j1) ! initializing before first digitization step DO k = (j1 + 1), j2 ! which will execute at least once uvec2(1:3) = trace(1:3, k) CALL DCross(uvec1, uvec2, tvec) IF (DMagnitude(tvec) > 0.0D0) THEN CALL DMake_uvec(tvec, pole_uvec) away_radians = ABS(DDot(pole_uvec, uvec3)) IF (away_radians < closest_radians) THEN ! close enough to be considered more seriously, below: azimuth_along_j = DCompass(from_uvec = uvec1, to_uvec = uvec2) back_azimuth_along_j = DCompass(from_uvec = uvec2, to_uvec = uvec1) !Consider azimuths to uvec3 from each end-point [uvec1, uvec2]; are they acceptable? tvec(1:3) = uvec3(1:3) - uvec1(1:3) ! although I could also compare uvec3 to uvec2; same result would follow. left_side = (DDot(tvec, pole_uvec) > 0.0D0) right_side = (DDot(tvec, pole_uvec) < 0.0D0) ! If BOTH are false, then end-point is ON the trace of fault #f. IF ((DArc(uvec1, uvec3) == 0.0D0).OR.(DArc(uvec2, uvec3) == 0.0D0)) THEN ! uvec3 matches uvec1 or uvec2; ON the line! left_side = .FALSE. right_side = .FALSE. END IF IF (left_side) THEN ! termination-point uvec3 is on LEFT side of uvec1-->uvec2 digitization step !Decide end1_OK? azimuth_to_loose_end = DCompass(from_uvec = uvec1, to_uvec = uvec3) IF (azimuth_to_loose_end > azimuth_along_j) azimuth_to_loose_end = azimuth_to_loose_end - two_Pi IF (k == (j1+1)) THEN ! this is the first digitization step on trace #j; there is no previous step azimuth_limit = azimuth_along_j - Pi_over_2 ! somewhat arbitrary; could be loosened. ELSE ! there was a previous digitization step; get its starting-point uvec4(1:3) = trace(1:3, (k-2)) ! only possible if k > (j1+1), of course azimuth_back = DCompass(from_uvec = uvec1, to_uvec = uvec4) IF (azimuth_back > azimuth_along_j) azimuth_back = azimuth_back - two_Pi azimuth_limit = azimuth_along_j + 0.5D0 * (azimuth_back - azimuth_along_j) END IF ! first digitization step? Or, not. end1_OK = (azimuth_to_loose_end >= azimuth_limit) !--------------------------------------------------- !Decide end2_OK? azimuth_to_loose_end = DCompass(from_uvec = uvec2, to_uvec = uvec3) ! now, measuring from other end IF (azimuth_to_loose_end < back_azimuth_along_j) azimuth_to_loose_end = azimuth_to_loose_end + two_Pi IF (k == j2) THEN ! this is the last digitization step on trace #j; there is no further step azimuth_limit = back_azimuth_along_j + Pi_over_2 ! somewhat arbitrary; could be loosened. ELSE ! there is an additional digitization step; get its ending-point uvec4(1:3) = trace(1:3, (k+1)) ! only possible if k < j2, of course azimuth_ahead = DCompass(from_uvec = uvec2, to_uvec = uvec4) IF (azimuth_ahead < back_azimuth_along_j) azimuth_ahead = azimuth_ahead + two_Pi azimuth_limit = back_azimuth_along_j + 0.5D0 * (azimuth_ahead - back_azimuth_along_j) END IF ! first digitization step? Or, not. end2_OK = (azimuth_to_loose_end <= azimuth_limit) ELSE IF (right_side) THEN ! termination-point uvec3 is on the RIGHT side of uvec1-->uvec2 digitization step !Decide end1_OK? azimuth_to_loose_end = DCompass(from_uvec = uvec1, to_uvec = uvec3) IF (azimuth_to_loose_end < azimuth_along_j) azimuth_to_loose_end = azimuth_to_loose_end + two_Pi IF (k == (j1+1)) THEN ! this is the first digitization step on trace #j; there is no previous step azimuth_limit = azimuth_along_j + Pi_over_2 ! somewhat arbitrary; could be loosened. ELSE ! there was a previous digitization step; get its starting-point uvec4(1:3) = trace(1:3, (k-2)) ! only possible of k > (j1+1), of course azimuth_back = DCompass(from_uvec = uvec1, to_uvec = uvec4) IF (azimuth_back < azimuth_along_j) azimuth_back = azimuth_back + two_Pi azimuth_limit = azimuth_along_j + 0.5D0 * (azimuth_back - azimuth_along_j) END IF ! first digitization step? Or, not. end1_OK = (azimuth_to_loose_end <= azimuth_limit) !--------------------------------------------------- !Decide end2_OK? azimuth_to_loose_end = DCompass(from_uvec = uvec2, to_uvec = uvec3) ! now, measuring from other end IF (azimuth_to_loose_end > back_azimuth_along_j) azimuth_to_loose_end = azimuth_to_loose_end - two_Pi IF (k == j2) THEN ! this is the last digitization step on trace #j; there is no further step azimuth_limit = back_azimuth_along_j - Pi_over_2 ! somewhat arbitrary; could be loosened. ELSE ! there is an additional digitization step; get its ending-point uvec4(1:3) = trace(1:3, (k+1)) ! only possible of k < j2, of course azimuth_ahead = DCompass(from_uvec = uvec2, to_uvec = uvec4) IF (azimuth_ahead > back_azimuth_along_j) azimuth_ahead = azimuth_ahead - two_Pi azimuth_limit = back_azimuth_along_j + 0.5D0 * (azimuth_ahead - back_azimuth_along_j) END IF ! first digitization step? Or, not. end2_OK = (azimuth_to_loose_end >= azimuth_limit) ELSE ! end-point is ON this digitization step (or its extrapolation; check for that case, too): t1 = DArc(uvec1, uvec3) t2 = DArc(uvec2, uvec3) t3 = DArc(uvec1, uvec2) end1_OK = (t2 <= t3) end2_OK = (t1 <= t3) END IF ! left_side, or right_side, or ON the trace alongside = end1_OK .AND. end2_OK IF (alongside) THEN ! This becomes the current best-match: closest_radians = away_radians closest_trace = j closest_k = k - j1 ! switching from absolute position in trace(1:3, 1:f_highest) to relative position within trace #j END IF ! alongside END IF ! away_radians < closest_radians (new closest-distance) END IF ! DMagnitude(tvec) > 0.0D0; digitization step has positive length uvec1 = uvec2 ! preparing to loop END DO ! k = j1+1, j2 (all digitization steps along trace #j) END IF ! consider_that; fault trace #j should be considered as possible superhighway to merge into. END IF ! j /= i END DO ! j = 1, f_highest IF (closest_radians < allowance) THEN end_type(i_end, i) = 4 link_list_length(i_end, i) = link_list_length(i_end, i) + 1 link_to_list(link_list_length(i_end, i), i_end, i) = closest_trace !Note that there is no reciprocal link from fault #j to fault-end (i_end, i). link_end_list(link_list_length(i_end, i), i_end, i) = closest_k ! recording # of digitized-trace segment on fault #j n_type_4 = n_type_4 + 1 IF (i_end == 1) THEN WRITE (*, "(' Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('Begining of: ', A)") TRIM(f_dig_faultName_lines(i)) ELSE WRITE (*, "(' End of: ', A)") TRIM(f_dig_faultName_lines(i)) WRITE (29, "('End of: ', A)") TRIM(f_dig_faultName_lines(i)) END IF WRITE (*, "(' Merges into: ', A)") TRIM(f_dig_faultName_lines(closest_trace)) WRITE (*, *) WRITE (29, "('Merges into: ', A)") TRIM(f_dig_faultName_lines(closest_trace)) WRITE (29, *) END IF IF ((end_type(i_end, i) == 0).AND.(trace_is(trace_loc(i_end, i))%element > 0)) THEN end_type(i_end, i) = 5 ! loose-end (termination in continuum crust, within .FEG area) n_type_5 = n_type_5 + 1 END IF END IF ! consider_this END DO ! i_end = 1, 2 END DO ! i = 1, f_highest WRITE (*, "(' -------------------------------------------------------------------------')") WRITE (*, "(' ', I6, ' fault mergers (of named fault-end with mid-trace) were found.')") n_type_4 WRITE (*, "(' ', I6, ' fault terminations (in unbroken crust) were found.')") n_type_5 WRITE (29, "('-------------------------------------------------------------------------')") WRITE (29, "(I6, ' fault mergers (of named fault-end with mid-trace) were found.')") n_type_4 WRITE (29, "(I6, ' fault terminations (in unbroken crust) were found.')") n_type_5 CLOSE (29) ! fault_connections.txt CALL Pause() !------------------------------------------------------------------------------------ !Install new nodes & elements needed at fault ends (taking care to re-use nodes when possible): WRITE (*, *) WRITE (*, "(' Installing new nodes & elements at ends/junctions of fault corridors...')") ALLOCATE ( ends_done(2, f_highest) ) ends_done = .FALSE. ! for now... ALLOCATE ( end_nodes(2, 2, f_highest) ) ! (1:2 = left/right (per digitizing direction); 1:2 = start/end of trace; 1:f_highest = FnnnnX) end_nodes = 0 ! for now... DO i = 1, f_highest IF (f_relevant(i)) THEN DO i_end = 1, 2 IF (.NOT.ends_done(i_end, i)) THEN ! we have not already encountered this junction/termination, so do it now. IF (end_type(i_end, i) == 1) THEN ! NOTE: No action if end_type == 0; such ends SHOULD be outside the .FEG area. ! 1 = fault leaves domain of the .FEG at its edge; corridor cut-off where last segment ends. IF (Restore) THEN corridor_width_radians = MAX(min_corridor_width_radians, MIN(max_corridor_width_radians, (opening_heave_km(i)/R_Earth_km))) ELSE ! NeoKinema corridor_width_radians = min_corridor_width_radians END IF radius = 0.5D0 * corridor_width_radians IF (i_end == 1) THEN uvec(1:3) = seg_end(1:3, 1, trace_loc(3, i)) ! outer end of last segment inner_uvec(1:3) = seg_end(1:3, 2, trace_loc(3, i)) ! inner end of same segment ELSE ! i_end == 2 uvec(1:3) = seg_end(1:3, 2, trace_loc(4, i)) ! outer end of last segment inner_uvec(1:3) = seg_end(1:3, 1, trace_loc(4, i)) ! inner end of same segment END IF !Note that the last element in this fault-corridor, at edge of grid, will have its outer edge perpendicular to fault trace; !this may be oblique to the overall trend of the model-area perimeter in that region. azimuth_out = DCompass(from_uvec = inner_uvec, to_uvec = uvec) CALL DTurn_To (azimuth_radians = azimuth_out - Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL Add_node(uvec1); n1 = numNod ! on the left as you look outward CALL DTurn_To (azimuth_radians = azimuth_out + Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) CALL Add_node(uvec2); n2 = numNod ! on the right as you look outward !Remember the pair of end-nodes that will bracket the constant-width corridor: IF (i_end == 1) THEN ! looking outwards is looking backwards relative to trace-digitization: end_nodes(1, 1, i) = n2 ! on left side of start of trace, looking inward end_nodes(2, 1, i) = n1 ! on right side of start of trace, looking inward ELSE ! i_end == 2; looking outward is looking along trace-digitization direction: end_nodes(1, 2, i) = n1 ! on left side of end of trace end_nodes(2, 2, i) = n2 ! on right side of end of trace END IF ELSE IF (end_type(i_end, i) == 2) THEN ! 2 = close-coupling to another (designated "master" or "symmetric_spreading_system") fault. ! First, decide whether we are dealing with a RT_corner angle, or a ~straight, ~constant-width corridor? j = link_to_list(1, i_end, i) ! find the (only; #1) other closely-linked fault Fnnnn# = j j_end = link_end_list(1, i_end, i) ! find which end of fault #j is involved c = f_dig_faultName_lines(i)(6:6) ! slip-sense character for fault #i spreading_on_i = (c == 'D').OR.(c == 'd') c = f_dig_faultName_lines(j)(6:6) ! slip-sense character for fault #j spreading_on_j = (c == 'D').OR.(c == 'd') RT_corner = (spreading_on_i.NEQV.spreading_on_j) ! when spreading is on one fault only IF (RT_corner) THEN IF (Restore) THEN ! create wide corridors around spreading centers, to allow at least ~1 m.y. of restoration per FEG IF (spreading_on_i) THEN spreading_km = opening_heave_km(i) ELSE ! spreading_on_j spreading_km = opening_heave_km(j) END IF large_corridor_width = MIN(MAX(min_corridor_width_radians, (spreading_km/R_Earth_km)), max_corridor_width_radians) small_corridor_width = min_corridor_width_radians ELSE ! F-E grid for use in NeoKinema; no need for wide corridors around spreading centers large_corridor_width = min_corridor_width_radians small_corridor_width = min_corridor_width_radians END IF !Create 3 nodes, all rather distant from the present corner (to allow for area-loss during reconstruction), and one element: IF (spreading_on_j) THEN ! The 3 new nodes lie close to fault #i or its projection. !For a ~reliable azimuth, use the terminal segment from fault #i (instead of just the last 2 digitized points): end_uvec(1:3) = seg_end(1:3, i_end, trace_loc((2+i_end), i)) inner_uvec(1:3) = seg_end(1:3, (3-i_end), trace_loc((2+i_end), i)) azimuth_along_i = DCompass(from_uvec = end_uvec, to_uvec = inner_uvec) ! pointing from the RT corner, into fault #i CALL DTurn_To (azimuth_radians = azimuth_along_i, base_uvec = end_uvec, far_radians = 0.5D0 * large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = middle_uvec) CALL DTurn_To (azimuth_radians = azimuth_along_i + Pi_over_2, base_uvec = middle_uvec, far_radians = 0.5D0 * small_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) CALL Add_node(uvec2); n2 = numNod CALL DTurn_To (azimuth_radians = azimuth_along_i - Pi_over_2, base_uvec = middle_uvec, far_radians = 0.5D0 * small_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) CALL Add_node(uvec3); n3 = numNod !Decide sense of transform #i: c = f_dig_faultName_lines(i)(6:6) ! slip-sense character for fault #i dextral = (c == 'R').OR.(c == 'r') IF (dextral) THEN CALL DTurn_To (azimuth_radians = azimuth_along_i + Pi, base_uvec = uvec2, far_radians = large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! projected from node2 ELSE ! fault #i is sinistral CALL DTurn_To (azimuth_radians = azimuth_along_i + Pi, base_uvec = uvec3, far_radians = large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! projected from node3 instead END IF ! fault #i is dextral or sinistral CALL Add_node(uvec1); n1 = numNod CALL Add_element(n1, n2, n3) !Remember end-nodes for the simple, ~straight sections of these 2 fault corridors: IF (i_end == 1) THEN end_nodes(1, i_end, i) = n3 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n2 ! on R side of end #i_end of fault #i, looking in digitization direction ELSE ! i_end = 2 end_nodes(1, i_end, i) = n2 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n3 ! on R side of end #i_end of fault #i, looking in digitization direction END IF ends_done(i_end, i) = .TRUE. IF (.NOT.ends_done(j_end, j)) THEN IF (j_end == 1) THEN IF (dextral) THEN ! end_nodes(1, j_end, j) = n1 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n3 ! on R side of end #j_end of fault #j, looking in digitization direction ELSE ! sinistral end_nodes(1, j_end, j) = n2 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n1 ! on R side of end #j_end of fault #j, looking in digitization direction END IF ELSE ! j_end == 2 IF (dextral) THEN ! end_nodes(1, j_end, j) = n3 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n1 ! on R side of end #j_end of fault #j, looking in digitization direction ELSE ! sinistral end_nodes(1, j_end, j) = n1 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n2 ! on R side of end #j_end of fault #j, looking in digitization direction END IF END IF ends_done(j_end, j) = .TRUE. END IF ELSE ! spreading_on_i; the 3 new nodes lie close to fault #j or its projection. !For a ~reliable azimuth, use the terminal segment from fault #j (instead of just the last 2 digitized points): end_uvec(1:3) = seg_end(1:3, j_end, trace_loc((2+j_end), j)) inner_uvec(1:3) = seg_end(1:3, (3-j_end), trace_loc((2+j_end), j)) azimuth_along_j = DCompass(from_uvec = end_uvec, to_uvec = inner_uvec) ! pointing from the RT corner, into fault #j CALL DTurn_To (azimuth_radians = azimuth_along_j, base_uvec = end_uvec, far_radians = 0.5D0 * large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = middle_uvec) CALL DTurn_To (azimuth_radians = azimuth_along_j + Pi_over_2, base_uvec = middle_uvec, far_radians = 0.5D0 * small_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) CALL Add_node(uvec2); n2 = numNod CALL DTurn_To (azimuth_radians = azimuth_along_j - Pi_over_2, base_uvec = middle_uvec, far_radians = 0.5D0 * small_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) CALL Add_node(uvec3); n3 = numNod !Decide sense of transform #j: c = f_dig_faultName_lines(j)(6:6) ! slip-sense character for fault #j dextral = (c == 'R').OR.(c == 'r') IF (dextral) THEN CALL DTurn_To (azimuth_radians = azimuth_along_j + Pi, base_uvec = uvec2, far_radians = large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! projected from node2 ELSE ! fault #i is sinistral CALL DTurn_To (azimuth_radians = azimuth_along_j + Pi, base_uvec = uvec3, far_radians = large_corridor_width, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! projected from node3 instead END IF ! fault #i is dextral or sinistral CALL Add_node(uvec1); n1 = numNod CALL Add_element(n1, n2, n3) !Remember end-nodes for the simple, ~straight sections of these 2 fault corridors: IF (i_end == 1) THEN IF (dextral) THEN ! end_nodes(1, i_end, i) = n1 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n3 ! on R side of end #i_end of fault #i, looking in digitization direction ELSE ! sinistral end_nodes(1, i_end, i) = n2 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n1 ! on R side of end #i_end of fault #i, looking in digitization direction END IF ELSE ! i_end == 2 IF (dextral) THEN ! end_nodes(1, i_end, i) = n3 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n1 ! on R side of end #i_end of fault #i, looking in digitization direction ELSE ! sinistral end_nodes(1, i_end, i) = n1 ! on L side of end #i_end of fault #i, looking in digitization direction end_nodes(2, i_end, i) = n2 ! on R side of end #i_end of fault #i, looking in digitization direction END IF END IF ends_done(i_end, i) = .TRUE. IF (.NOT.ends_done(j_end, j)) THEN IF (j_end == 1) THEN end_nodes(1, j_end, j) = n3 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n2 ! on R side of end #j_end of fault #j, looking in digitization direction ELSE ! j_end = 2 end_nodes(1, j_end, j) = n2 ! on L side of end #j_end of fault #j, looking in digitization direction end_nodes(2, j_end, j) = n3 ! on R side of end #j_end of fault #j, looking in digitization direction END IF ends_done(j_end, j) = .TRUE. END IF END IF ! spreading_on_j, or spreading_on_i [two blocks of exactly parallel code, except swapping 'j' for 'i', and vice versa] ELSE ! an ordinary close-coupled junction in a ~straight, ~constant-width corridor: IF (Restore) THEN corridor_width_radians = MAX(min_corridor_width_radians, MIN(max_corridor_width_radians, (opening_heave_km(i)/R_Earth_km))) ELSE ! NeoKinema corridor_width_radians = min_corridor_width_radians END IF radius = 0.5D0 * corridor_width_radians tvec(1:3) = 0.5D0 * trace(1:3, trace_loc(i_end, i)) ! first half of sum (over 2 closely-linked end-nodes) tvec(1:3) = tvec(1:3) + 0.5D0 * trace(1:3, trace_loc(j_end, j)) ! second half of the sum CALL DMake_uvec(tvec, uvec) ! this uvec is now the center of the joint !For a ~reliable azimuth, use the terminal segment from fault #i (instead of just 2 digitized points): uvec1(1:3) = seg_end(1:3, i_end, trace_loc((2+i_end), i)) uvec2(1:3) = seg_end(1:3, (3-i_end), trace_loc((2+i_end), i)) azimuth_along_i = DCompass(from_uvec = uvec1, to_uvec = uvec2) ! points inward from active end of fault #i; not NECESSARILY == digitization direction! IF (i_end == 1) THEN ! digitization direction == "along_i" direction CALL DTurn_To (azimuth_radians = azimuth_along_i - Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL DTurn_To (azimuth_radians = azimuth_along_i + Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) ELSE ! digitization direction is OPPOSITE to "along_i" direction CALL DTurn_To (azimuth_radians = azimuth_along_i + Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL DTurn_To (azimuth_radians = azimuth_along_i - Pi_over_2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) END IF CALL Add_node(uvec1); n1 = numNod ! on the left as you look along fault #i in the digitization direction CALL Add_node(uvec2); n2 = numNod ! on the right as you look along fault #i in the digitization direction !Remember the pair of end-nodes that will bracket the constant-width corridor: end_nodes(1, i_end, i) = n1 ! on left side of end #i_end of trace #i end_nodes(2, i_end, i) = n2 ! on right side of end #i_end of trace #i ends_done(i_end, i) = .TRUE. IF (.NOT.ends_done(j_end, j)) THEN ! do it now !Check sense of digitization of fault #j, relative to fault #i and perhaps swap these node #s? IF (j_end == i_end) THEN ! faults #i and #j were digitized in OPPOSITE directions end_nodes(1, j_end, j) = n2 ! on left side of end #j_end of trace #j end_nodes(2, j_end, j) = n1 ! on right side of end #j_end of trace #j ELSE ! they were digitized in the same direction end_nodes(1, j_end, j) = n1 ! on left side of end #j_end of trace #j end_nodes(2, j_end, j) = n2 ! on right side of end #j_end of trace #j END IF ends_done(j_end, j) = .TRUE. ! flag, to avoid creating joint nodes & elements more than once in a junction END IF ! .NOT. ends_done(j_end, j) END IF ! RT_corner, or not ELSE IF (end_type(i_end, i) == 3) THEN ! 3 = tectonic knot (non-close duple junction, 3-J, or N-J); IF (.NOT.ends_done(i_end, i)) THEN n_faults = 1 + link_list_length(i_end, i) !Locate reference point at average location of all fault-ends: tvec(1:3) = trace(1:3, trace_loc(i_end, i)) DO j = 1, link_list_length(i_end, i) k = link_to_list(j, i_end, i) ! Fnnnn# of nearby trace k_end = link_end_list(j, i_end, i) ! which end of this nearby trace? tvec(1:3) = tvec(1:3) + trace(1:3, trace_loc(k_end, k)) END DO CALL DMake_uvec(tvec, uvec) ! average position of all nearby fault-ends. [Note; no need to divide /n_faults] !Find minimum radius that will encompass all fault-ends (then, multiply by a factor-of-safety): uvec1(1:3) = trace(1:3, trace_loc(i_end, i)) radius = DArc(uvec, uvec1) ! from constructed central-point to end of trace #i DO j = 1, link_list_length(i_end, i) k = link_to_list(j, i_end, i) ! Fnnnn# of nearby trace k_end = link_end_list(j, i_end, i) ! which end of this nearby trace? uvec2(1:3) = trace(1:3, trace_loc(k_end, k)) radius = MAX(radius, DArc(uvec, uvec2)) END DO radius = radius * 2.0D0 ! apply factor-of-safety IF (Restore) THEN ! consider possibility of larger radius, based on opening_heave_km() for all faults? radius = MAX(radius, (0.6 * opening_heave_km(i) / R_Earth_km)) DO j = 1, link_list_length(i_end, i) k = link_to_list(j, i_end, i) ! Fnnnn# of nearby trace radius = MAX(radius, (0.6 * opening_heave_km(k) / R_Earth_km)) END DO END IF radius = MAX(radius, 0.6D0 * min_corridor_width_radians) ! apply pre-defined reasonableness constraints radius = MIN(radius, 0.6D0 * max_corridor_width_radians) !Make a list of the azimuths at which each fault departs from the constructed center: azimuth_radians_list = 0.0D0 ! whole list, to assist debugging uvec1(1:3) = seg_end(1:3, i_end, trace_loc((2+i_end), i)) uvec2(1:3) = seg_end(1:3, (3-i_end), trace_loc((2+i_end), i)) azimuth_radians_list(0) = DCompass(from_uvec = uvec1, to_uvec = uvec2) ! Note use of subscript (0). DO j = 1, link_list_length(i_end, i) k = link_to_list(j, i_end, i) ! Fnnnn# of nearby trace k_end = link_end_list(j, i_end, i) ! which end of this nearby trace? uvec1(1:3) = seg_end(1:3, k_end, trace_loc((2+k_end), k)) uvec2(1:3) = seg_end(1:3, (3-k_end), trace_loc((2+k_end), k)) temp_R8 = DCompass(from_uvec = uvec1, to_uvec = uvec2) IF (temp_R8 < azimuth_radians_list(0)) temp_R8 = temp_R8 + two_Pi ! Be sure that all azimuths are greater than the (0) azimuth! azimuth_radians_list(j) = temp_R8 END DO !Sort faults in this tectonic knot (IF n_faults > 2) so that azimuths are in ascending order (clockwise), beginning with fault #i: IF (n_faults > 2) THEN DO j = 1, (link_list_length(i_end, i) - 1) DO k = (j + 1), link_list_length(i_end, i) IF (azimuth_radians_list(k) < azimuth_radians_list(j)) THEN ! swap this pair: temp_R8 = azimuth_radians_list(j) temp_int_1 = link_to_list(j, i_end, i) temp_int_2 = link_end_list(j, i_end, i) azimuth_radians_list(j) = azimuth_radians_list(k) link_to_list(j, i_end, i) = link_to_list(k, i_end, i) link_end_list(j, i_end, i) = link_end_list(k, i_end, i) azimuth_radians_list(k) = temp_R8 link_to_list(k, i_end, i) = temp_int_1 link_end_list(k, i_end, i) = temp_int_2 END IF ! swap is needed END DO ! inner loop of sort END DO ! outer loop of sort END IF ! sorting of nearby faults (by azimuth) is necessary azimuth_radians_list(n_faults) = two_Pi + azimuth_radians_list(0) ! copying the azimuth of fault #i to the end-of-list (as cyclicity eases programming). !Create nodes that split these azimuths (halfway-between), at distance radius: new_nodes = 0 ! whole list, to aid debugging temp_R8 = (azimuth_radians_list(0) + azimuth_radians_list(1)) / 2.0D0 CALL DTurn_To (azimuth_radians = temp_R8, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL Add_node(uvec1); new_nodes(1) = numNod ! Note that subscript of new_nodes() goes from 1...n_faults, or (n_faults+1) == (1). DO j = 1, link_list_length(i_end, i) temp_R8 = (azimuth_radians_list(j) + azimuth_radians_list(j+1)) / 2.0D0 CALL DTurn_To (azimuth_radians = temp_R8, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL Add_node(uvec1); new_nodes(j+1) = numNod ! This fills in new_nodes(2...n_faults) END DO new_nodes(n_faults + 1) = new_nodes(1) ! because cyclicity eases programming. !If necessary, create elements at this 3-J or n-J: IF (n_faults == 2) THEN !No action is needed; no small element added in junction. ELSE IF (n_faults == 3) THEN !Add a single element: n1 = new_nodes(3) n2 = new_nodes(2) ! Note that I reverse the order: clockwise ==> counter-clockwise. n3 = new_nodes(1) CALL Add_element(n1, n2, n3) ELSE ! n_faults >= 4 (rare case); add multiple local elements CALL Add_node(uvec); new_nodes(0) = numNod ! one more node is needed, at center of this junction DO j = 1, n_faults n1 = new_nodes(0) n2 = new_nodes(j+1) n3 = new_nodes(j) CALL Add_element(n1, n2, n3) END DO END IF !Record end-nodes for fault #i (always). IF (i_end == 1) THEN ! away-from-junction == digitization direction end_nodes(1, i_end, i) = new_nodes(n_faults) ! on L side of end # i_end of fault #i, looking in the digitization direction end_nodes(2, i_end, i) = new_nodes(1) ! on R side of end # i_end of fault #i, looking in the digitization direction ELSE ! into-junction == digitization direction end_nodes(1, i_end, i) = new_nodes(1) ! on L side of end # i_end of fault #i, looking in the digitization direction end_nodes(2, i_end, i) = new_nodes(n_faults) ! on R side of end # i_end of fault #i, looking in the digitization direction END IF ends_done(i_end, i) = .TRUE. !Consider whether this junction should be recorded (if not done already) on each of the nearby faults? DO j = 1, link_list_length(i_end, i) k = link_to_list(j, i_end, i) ! Fnnnn# of nearby trace k_end = link_end_list(j, i_end, i) ! which end of this nearby trace? IF (.NOT.ends_done(k_end, k)) THEN IF (k_end == 1) THEN ! away-from-junction == digitization direction end_nodes(1, k_end, k) = new_nodes(j) ! on L side of end # k_end of fault #k, looking in the digitization direction end_nodes(2, k_end, k) = new_nodes(j+1) ! on R side of end # k_end of fault #k, looking in the digitization direction ELSE ! into-junction == digitization direction end_nodes(1, k_end, k) = new_nodes(j+1) ! on L side of end # k_end of fault #k, looking in the digitization direction end_nodes(2, k_end, k) = new_nodes(j) ! on R side of end # k_end of fault #k, looking in the digitization direction END IF ends_done(k_end, k) = .TRUE. END IF ! end # k_end of fault #k has not been recorded. END DO END IF ! .NOT.ends_done(i_end, i) ELSE IF (end_type(i_end, i) == 4) THEN ! 4 = merging, where current fault ends anywhere on the trace of another fault. j = link_to_list(1, i_end, i) ! # of fault which is merged-into. k = link_end_list(1, i_end, i) ! # of digitization step within fault #j !Build entries describing 3 nodes & 1 element at the on-ramp joint ! (without worrying about the ordering of on-ramps along fault #j): !(1) Describe great-circle containing this end of fault #i: uvec1(1:3) = trace(1:3, trace_loc(i_end, i)) ! outer end of trace of fault #i IF (i_end == 1) THEN uvec2(1:3) = seg_end(1:3, 2, trace_loc(3, i)) ! inner end of terminal segment on fault #i ELSE ! i_end == 2 uvec2(1:3) = seg_end(1:3, 1, trace_loc(4, i)) ! inner end of terminal segment on fault #i END IF azimuth_into_i = DCompass(uvec1, uvec2) CALL DCross(uvec1, uvec2, tvec) CALL DMake_uvec(tvec, pole_1_uvec) ! pole of great circle describing terminal segment of fault #i !(2) Describe great-circle containing the critical digitization step on #j: uvec3(1:3) = trace(1:3, (trace_loc(1, j) + k - 1)) uvec4(1:3) = trace(1:3, (trace_loc(1, j) + k)) azimuth_along_j = DCompass(uvec3, uvec4) CALL DCross(uvec3, uvec4, tvec) CALL DMake_uvec(tvec, pole_2_uvec) ! pole of great circle describing critical digitization step on fault #j !(3) Spread the end-points (along the digitization step #k of fault #j), to provide for oblique-approach junctions: uvec3(1:3) = uvec3(1:3) - 1.00D0 * (uvec4(1:3) - uvec3(1:3)) CALL DMake_uvec(uvec3, uvec3) ! (Note that next formula is different, as uvec3 has already moved.) uvec4(1:3) = uvec4(1:3) + 0.50D0 * (uvec4(1:3) - uvec3(1:3)) CALL DMake_uvec(uvec4, uvec4) !(4) Find point of intersection: CALL DCircles_Intersect (pole_a_uvec = pole_1_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = uvec1, & ! <--- deliberately extrapolating trace of fault #i into a full great circle & pole_b_uvec = pole_2_uvec, dot_b = 0.0D0, first_b_uvec = uvec3, last_b_uvec = uvec4, & ! <--- but, requiring intersection on THIS side of globe & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number < 1) THEN !This happens sometimes, when end of fault #i is ~parallel to digitization-step #k of fault #j. !In that case, modify the trace of fault #i by adding one more step, deflected by 45 degrees: !Decide whether end *i_end of trace #i needs to deflect to the left? tvec1(1:3) = uvec1(1:3) - uvec2(1:3) ! short tvec1 points outward from end #i_end of trace #i tvec2(1:3) = ((uvec3(1:3) + uvec4(1:3)) / 2.0D0) - uvec2(1:3) ! short tvec2 points from uvec2 to middle of dig. step #k of trace #j CALL DCross(tvec1, tvec2, tvec3) ! now, tvec3 should point up if the necessary turn is to the left left_turn = (DDot(tvec3, uvec1) > 0.0D0) IF (left_turn) THEN azimuth = DCompass(from_uvec = uvec2, to_uvec = uvec1) - 0.25D0 * Pi ! 45 degrees to the left azimuth_into_i = azimuth_into_i - 0.25D0 * Pi ELSE ! right turn azimuth = DCompass(from_uvec = uvec2, to_uvec = uvec1) + 0.25D0 * Pi ! 45 degrees to the right azimuth_into_i = azimuth_into_i + 0.25D0 * Pi END IF !Find a far_uvec which is about 100 km away along this new direction: CALL DTurn_To (azimuth_radians = azimuth, base_uvec = uvec1, far_radians = 0.0157D0, & ! inputs (distance 0.0157 radians = 100 km) & omega_uvec = omega_uvec, result_uvec = far_uvec) !Try again for an intersection. This time, extend the dig-step #k of trace #j to a great-circle. Make the other arc 100 km long: CALL DCross(uvec1, far_uvec, tvec) CALL DMake_uvec(tvec, pole_1_uvec) CALL DCircles_Intersect (pole_a_uvec = pole_1_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = far_uvec, & ! arc of 100 km & pole_b_uvec = pole_2_uvec, dot_b = 0.0D0, first_b_uvec = uvec3, last_b_uvec = uvec3, & ! great circle & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number < 1) THEN WRITE (*, "(' CAUTION: Could not find type-4 termination of end ', I2,' of fault ', I4,' against fault ',I4)") i_end, i, j WRITE (*, "(' So, just averaging together the 3 known points for an approximate location.')") tvec(1:3) = uvec1(1:3) + uvec3(1:3) + uvec4(1:3) CALL DMake_uvec(tvec, point1_uvec) END IF ! no intersection found on second attempt END IF ! no intersection found on first attempt uvec(1:3) = point1_uvec(1:3) ! best-estimate of intersection point !(5) Build a new triangular element around this point: IF (Restore) THEN !Choose corridor-width based on the larger of the two opening_heave_km(i OR j): corridor_width_radians = MAX(min_corridor_width_radians, MIN(max_corridor_width_radians, & & (MAX(opening_heave_km(i), opening_heave_km(j))/R_Earth_km) ) ) ELSE ! NeoKinema corridor_width_radians = min_corridor_width_radians END IF radius = corridor_width_radians * 0.57735D0 IF (SIN(azimuth_into_i - azimuth_along_j) < 0.0D0) THEN ! fault #i is counter-clockwise from fault #j; so, fault #j is clockwise from #i azimuth1 = azimuth_along_j + Pi_over_2 CALL DTurn_To (azimuth_radians = azimuth1, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! redefining uvec3 as third node location CALL Add_node(uvec1); n1 = numNod azimuth2 = ATAN2( ((SIN(azimuth_into_i) + SIN(azimuth_along_j)) / 2.0D0), & & ((COS(azimuth_into_i) + COS(azimuth_along_j)) / 2.0D0) ) ! This complex way of averaging is insensitive to cycle-shifts. CALL DTurn_To (azimuth_radians = azimuth2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) ! redefining uvec1 as first node location CALL Add_node(uvec2); n2 = numNod azimuth3 = ATAN2( ((SIN(azimuth_along_j + Pi) + SIN(azimuth_into_i)) / 2.0D0), & & ((COS(azimuth_along_j + Pi) + COS(azimuth_into_i)) / 2.0D0) ) ! This complex way of averaging is insensitive to cycle-shifts. CALL DTurn_To (azimuth_radians = azimuth3, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) ! redefining uvec2 as second node location CALL Add_node(uvec3); n3 = numNod IF (i_end == 1) THEN end_nodes(1, 1, i) = n3 ! on left side of start of trace end_nodes(2, 1, i) = n2 ! on right side of start of trace ELSE ! i_end == 2 end_nodes(1, 2, i) = n2 ! on left side of end of trace end_nodes(2, 2, i) = n3 ! on right side of end of trace END IF ELSE ! fault #i is clockwise from fault #j, so fault #j is counter-clockwise from #i azimuth1 = azimuth_along_j - Pi_over_2 CALL DTurn_To (azimuth_radians = azimuth1, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) ! redefining uvec2 as second node location CALL Add_node(uvec1); n1 = numNod azimuth2 = ATAN2( ((SIN(azimuth_along_j + Pi) + SIN(azimuth_into_i)) / 2.0D0), & & ((COS(azimuth_along_j + Pi) + COS(azimuth_into_i)) / 2.0D0) ) ! This complex way of averaging is insensitive to cycle-shifts. CALL DTurn_To (azimuth_radians = azimuth2, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) ! redefining uvec3 as third node location CALL Add_node(uvec2); n2 = numNod azimuth3 = ATAN2( ((SIN(azimuth_into_i) + SIN(azimuth_along_j)) / 2.0D0), & & ((COS(azimuth_into_i) + COS(azimuth_along_j)) / 2.0D0) ) ! This complex way of averaging is insensitive to cycle-shifts. CALL DTurn_To (azimuth_radians = azimuth3, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) ! redefining uvec1 as first node location CALL Add_node(uvec3); n3 = numNod IF (i_end == 1) THEN end_nodes(1, 1, i) = n3 ! on left side of start of trace end_nodes(2, 1, i) = n2 ! on right side of start of trace ELSE ! i_end == 2 end_nodes(1, 2, i) = n2 ! on left side of end of trace end_nodes(2, 2, i) = n3 ! on right side of end of trace END IF END IF CALL Add_element(n1, n2, n3) ends_done(i_end, i) = .TRUE. onramp_count(j) = onramp_count(j) + 1 ! add another on-ramp to fault #j !Determine whether fault #i is (locally) on the left-side, as one proceeds in digitization-direction along trace #j? !(a) Get uvec of first internal point near the active end of trace #i: IF (i_end == 1) THEN far_uvec(1:3) = trace(1:3, (trace_loc(1, i) + 1)) ELSE ! i_end == 2 far_uvec(1:3) = trace(1:3, (trace_loc(2, i) - 1)) END IF !(b) Re-establish uvec3, uvec4 as 2 consecutive points on the trace of fault #j around the intersection. uvec3(1:3) = trace(1:3, (trace_loc(1, j) + k - 1)) uvec4(1:3) = trace(1:3, (trace_loc(1, j) + k)) !(c) Create two ~horizontal, non-unit vectors pointing away from intersection tvec1(1:3) = uvec4(1:3) - uvec3(1:3) ! short vector points ~horizontally along trace #j, in digitization direction tvec2(1:3) = far_uvec(1:3) - uvec(1:3) ! short vector points ~horizontally from intersection point into trace #i CALL DCross(tvec1, tvec2, tvec3) ! New tvec3 points up if left-side: left_side = (DDot(tvec3, uvec3) > 0.0D0) IF (left_side) THEN onramp_nodes(1, 1, onramp_count(j), j) = n3 ! left side, nearer to start onramp_nodes(2, 1, onramp_count(j), j) = n1 ! right side, nearer to start onramp_nodes(1, 2, onramp_count(j), j) = n2 ! left side, further from start onramp_nodes(2, 2, onramp_count(j), j) = n1 ! right side, further from start ELSE ! fault #i is (generally, globally) on the right side, as one proceeds in digitization-direction along trace #j. onramp_nodes(1, 1, onramp_count(j), j) = n1 ! left side, nearer to start onramp_nodes(2, 1, onramp_count(j), j) = n2 ! right side, nearer to start onramp_nodes(1, 2, onramp_count(j), j) = n1 ! left side, further from start onramp_nodes(2, 2, onramp_count(j), j) = n3 ! right side, further from start END IF ! left_side, or not. Note that, in either case, n1 is the node-number which appears twice. ELSE IF (end_type(i_end, i) == 5) THEN ! 5 = fault termination inside .FEG area; corridor will taper to one node. IF (Restore) THEN corridor_width_radians = MAX(min_corridor_width_radians, MIN(max_corridor_width_radians, (opening_heave_km(i)/R_Earth_km))) ELSE corridor_width_radians = min_corridor_width_radians END IF radius = corridor_width_radians * 0.57735D0 uvec(1:3) = trace(1:3, trace_loc(i_end, i)) ! end of trace will be centroid of new element !Take "pointing" direction of end-cap triangles from overall trace of fault: far_uvec(1:3) = trace(1:3, trace_loc((3 - i_end), i)) ! other end of fault trace azimuth_out = DCompass(from_uvec = far_uvec, to_uvec = uvec) CALL DTurn_To (azimuth_radians = azimuth_out, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL Add_node(uvec1); n1 = numNod CALL DTurn_To (azimuth_radians = azimuth_out - 2.094395D0, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) CALL Add_node(uvec2); n2 = numNod CALL DTurn_To (azimuth_radians = azimuth_out + 2.094395D0, base_uvec = uvec, far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) CALL Add_node(uvec3); n3 = numNod CALL Add_element(n1, n2, n3) ! Note that numEl increases by 1. !Remember the pair of end-nodes that will bracket the constant-width corridor: IF (i_end == 1) THEN end_nodes(1, 1, i) = n3 ! on left side of start of trace end_nodes(2, 1, i) = n2 ! on right side of start of trace ELSE ! i_end == 2 end_nodes(1, 2, i) = n2 ! on left side of end of trace end_nodes(2, 2, i) = n3 ! on right side of end of trace END IF END IF ! end_type(i_end, i) == 1, 2, 3, 4, 5 [Note: No action if end_type() == 0; such ends should be outside the .FEG area. END IF ! .NOT.ends_done(i_end, i) ends_done(i_end, i) = .TRUE. ! flag, to avoid creating joint nodes & elements more than once in a junction. (This may be redundant, if also set within paragraphs above.) END DO ! i_end = 1, 2 END IF ! f_relevant(i) END DO ! i = 1, f_highest WRITE (*, "(' ... DONE!')") !------------------------------------------------------------------------------------ !Define internal coordinate 's' (in range 0.0 <= s <= 1.0) for each type-4 fault-merger: DO i = 1, f_highest n = onramp_count(i) IF (n > 0) THEN DO j = 1, n ! (j = onramp #, within local list for this fault) uvec1(1:3) = trace(1:3, trace_loc(1, i)) ! start of trace uvec2(1:3) = trace(1:3, trace_loc(2, i)) ! end of trace IF (onramp_nodes(1, 1, j, i) == onramp_nodes(1, 2, j, i)) THEN ! unique single node is on the left side of trace #1 uvec3(1:3) = node_uvec(1:3, onramp_nodes(1, 1, j, i)) ! or, onramp_nodes(1, 2, j, i) would be the same ELSE ! unique single node is on the right side uvec3(1:3) = node_uvec(1:3, onramp_nodes(2, 1, j, i)) ! or, onramp_nodes(2, 2, j, i) would be the same END IF ! unique, single node is on the left/right side of trace #i t13 = DArc(uvec1, uvec3) t12 = DArc(uvec1, uvec2) onramp_s(j, i) = t13 / t12 !onramp_s(j, i) = MIN(1.0D0, MAX(0.0D0, onramp_s(j, i))) ! <-- commented-out because slightly out-of-range values may be useful for sorting; ! while two values that were both limited to 1.0 (or 0.0) would be indistinguishable. END DO ! j = 1, n (j = onramp #, within local list for this fault) END IF ! n > 0 END DO ! i = 1, f_highest !------------------------------------------------------------------------------------ !Sort the lists of onramps (for each fault) according to onramp_s(): DO i = 1, f_highest n = onramp_count(i) IF (n > 1) THEN ! list of at least 2 entries; needs to be sorted DO j = 1, (n-1) ! lower-index onramps DO k = (j+1), n ! higher-index onramps IF (onramp_s(k, i) < onramp_s(j, i)) THEN ! a swap is needed, between #j and #k in local list s_temp = onramp_s(j, i) ! saving a copy of #j n11 = onramp_nodes(1, 1, j, i) ! saving a copy of #j n12 = onramp_nodes(1, 2, j, i) ! saving a copy of #j n21 = onramp_nodes(2, 1, j, i) ! saving a copy of #j n22 = onramp_nodes(2, 2, j, i) ! saving a copy of #j onramp_s(j, i) = onramp_s(k, i) ! overwriting #j with #k onramp_nodes(1, 1, j, i) = onramp_nodes(1, 1, k, i) ! overwriting #j with #k onramp_nodes(1, 2, j, i) = onramp_nodes(1, 2, k, i) ! overwriting #j with #k onramp_nodes(2, 1, j, i) = onramp_nodes(2, 1, k, i) ! overwriting #j with #k onramp_nodes(2, 2, j, i) = onramp_nodes(2, 2, k, i) ! overwriting #j with #k onramp_s(k, i) = s_temp ! applying saved value (originally from #j) onramp_nodes(1, 1, k, i) = n11 ! applying saved value (originally from #j) onramp_nodes(1, 2, k, i) = n12 ! applying saved value (originally from #j) onramp_nodes(2, 1, k, i) = n21 ! applying saved value (originally from #j) onramp_nodes(2, 2, k, i) = n22 ! applying saved value (originally from #j) END IF ! a swap is needed, between #j and #k in local list END DO ! k = (j+1), n END DO ! j = 1, (n-1) END IF ! n > 1 END DO ! i = 1, f_highest !------------------------------------------------------------------------------------ !Build corridors along each fault, using a mix of pre-established and new nodes. WRITE (*, *) WRITE (*, "(' Completing fault corridors, with elements between junctions/ends...')") DO i = 1, f_highest IF (f_relevant(i)) THEN n11 = end_nodes(1, 1, i) n21 = end_nodes(2, 1, i) n12 = end_nodes(1, 2, i) n22 = end_nodes(2, 2, i) IF ((n11 > 0).AND.(n22 > 0).AND.(n21 > 0).AND.(n22 > 0)) THEN ! checking that all 4 bounding nodes are pre-defined? IF (Restore) THEN ! variable corridor width and length-steps, depending on fault offsets: dextral = (right_heave_km(i) >= 0.0D0) ! for program Restore corridor_length_radians = MIN( MAX( ABS(right_heave_km(i) / R_Earth_km), (4.0D0 * min_corridor_width_radians) ), max_corridor_width_radians) ! initial suggestion; will be adjusted up or down to fit in existing gap ELSE ! F-E grid for use in NeoKinema dextral = .TRUE. ! NeoKinema convention: diagonals of corridor elements reflect dextral shear in the past corridor_length_radians = 4.0D0 * min_corridor_width_radians ! initial suggestion; will be adjusted up or down to fit in existing gap END IF IF (onramp_count(i) == 0) THEN ! simple case; no on-ramps uvec1(1:3) = node_uvec(1:3, n11) uvec2(1:3) = node_uvec(1:3, n12) arc_length = DArc(uvec1, uvec2) n_steps = NINT(arc_length / corridor_length_radians) n_steps = MAX(1, n_steps) IF ((end_type(1, i) == 5).AND.(end_type(2, i) == 5)) n_steps = MAX(2, n_steps) ! to allow a triangular profile of slip (for isolated fault) corridor_length_radians = arc_length / n_steps ! final answer. (Note: OK to change value because there is only one step for this fault.) n1past = n11 ! left starting node for first step n2past = n21 ! right starting node for first step DO k = 1, n_steps IF (k == n_steps) THEN ! final step, so targets are predefined: n1target = n12 n2target = n22 ELSE ! must create 2 new nodes, to be the targets !(a) Initial, rough, straight-line estimates of positions of new nodes: uvec3(1:3) = node_uvec(1:3, n11) + ((k * 1.0D0)/(n_steps * 1.0D0)) * (node_uvec(1:3, n12) - node_uvec(1:3, n11)) CALL DMake_uvec(uvec3, uvec3) ! a small correction uvec4(1:3) = node_uvec(1:3, n21) + ((k * 1.0D0)/(n_steps * 1.0D0)) * (node_uvec(1:3, n22) - node_uvec(1:3, n21)) CALL DMake_uvec(uvec4, uvec4) ! a small correction !(b) A great-circle through these 2 positions defines the path along which we search for intersection with fault trace: CALL DCross(uvec3, uvec4, tvec) CALL DMake_uvec(tvec, pole_1_uvec) !(c) Try each digitization-step of the fault trace to see if it intersects? j1 = trace_loc(1, i) j2 = trace_loc(2, i) plodding: DO j = j1, (j2-1) ! different starting points (with absolute, large-# indexing) uvec5(1:3) = trace(1:3, j) uvec6(1:3) = trace(1:3, (j+1)) CALL DCross(uvec5, uvec6, tvec) CALL DMake_uvec(tvec, pole_2_uvec) ! for arc-of-great-circle representing this digitization step CALL DCircles_Intersect (pole_a_uvec = pole_1_uvec, dot_a = 0.0D0, first_a_uvec = uvec3, last_a_uvec = uvec3, & ! complete great-circle & pole_b_uvec = pole_2_uvec, dot_b = 0.0D0, first_b_uvec = uvec5, last_b_uvec = uvec6, & ! tiny arc of another great-circle & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number == 1) EXIT plodding ! with answer in point1_uvec END DO plodding ! j = j1, (j2-1) ! different starting points (with absolute, large-# indexing) correction(1:3) = point1_uvec(1:3) - 0.5D0 * (uvec3(1:3) + uvec4(1:3)) uvec3(1:3) = uvec3(1:3) + correction(1:3) CALL DMake_uvec(uvec3, uvec3) uvec4(1:3) = uvec4(1:3) + correction(1:3) CALL DMake_uvec(uvec4, uvec4) CALL Add_node(uvec3); n1target = numNod CALL Add_node(uvec4); n2target = numNod END IF ! using pre-defined, or new, node-pair for targets IF (dextral) THEN ! element diagonals should reflect right-lateral shear in the past CALL Add_element(n1past, n2past, n1target) CALL Add_element(n2target, n1target, n2past) ELSE ! element diagonals should reflect left-lateral shear in the past CALL Add_element(n2past, n2target, n1past) CALL Add_element(n1target, n1past, n2target) END IF ! dextral? or sinistral shear of element diagonals? n1past = n1target ! The former target becomes the new starting point, n2past = n2target ! to prepare for the next time through this loop. END DO ! k = 1, n_steps ELSE ! complex case, with on-ramps; onramp_count(i) > 0 DO n = 0, onramp_count(i) !Define starting and ending nodes (origin --> future; on left side of trace) for this segment: IF (n == 0) THEN ! beginning a new fault n1origin = n11 ! left starting node for first step within current segment uvec1(1:3) = node_uvec(1:3, n1origin) n2origin = n21 ! right starting node for first step within current segment ELSE ! not the first segment for this fault; starting at an onramp: n1origin = onramp_nodes(1, 2, n, i) uvec1(1:3) = node_uvec(1:3, n1origin) n2origin = onramp_nodes(2, 2, n, i) END IF IF (n == onramp_count(i)) THEN ! last segment for this fault n1future = n12 uvec2(1:3) = node_uvec(1:3, n1future) n2future = n22 ELSE ! not the last segment for this fault; segment ends at an onramp n1future = onramp_nodes(1, 1, (n + 1), i) uvec2(1:3) = node_uvec(1:3, n1future) n2future = onramp_nodes(2, 1, (n + 1), i) END IF arc_length = DArc(uvec1, uvec2) temp_corridor_length_radians = corridor_length_radians ! The left value will get changed in each DO; the right value will not vary. n_steps = NINT(arc_length / temp_corridor_length_radians) n_steps = MAX(1, n_steps) temp_corridor_length_radians = arc_length / n_steps ! final answer, for THIS pass through the DO loop of steps on one fault n1past = n1origin n2past = n2origin DO k = 1, n_steps IF (k == n_steps) THEN ! final step, so targets are predefined: n1target = n1future n2target = n2future ELSE ! must create 2 new nodes, to be the targets !(a) Initial, rough, straight-line estimates of positions of new nodes: uvec3(1:3) = node_uvec(1:3, n1origin) + ((k * 1.0D0)/(n_steps * 1.0D0)) * (node_uvec(1:3, n1future) - node_uvec(1:3, n1origin)) CALL DMake_uvec(uvec3, uvec3) ! a small correction uvec4(1:3) = node_uvec(1:3, n2origin) + ((k * 1.0D0)/(n_steps * 1.0D0)) * (node_uvec(1:3, n2future) - node_uvec(1:3, n2origin)) CALL DMake_uvec(uvec4, uvec4) ! a small correction !(b) A great-circle through these 2 positions defines the path along which we search for intersection with fault trace: CALL DCross(uvec3, uvec4, tvec) CALL DMake_uvec(tvec, pole_1_uvec) !(c) Try each digitization-step of the fault trace to see if it intersects? j1 = trace_loc(1, i) j2 = trace_loc(2, i) replodding: DO j = j1, (j2-1) ! different starting points (with absolute, large-# indexing) uvec5(1:3) = trace(1:3, j) uvec6(1:3) = trace(1:3, (j+1)) CALL DCross(uvec5, uvec6, tvec) CALL DMake_uvec(tvec, pole_2_uvec) ! for arc-of-great-circle representing this digitization step CALL DCircles_Intersect (pole_a_uvec = pole_1_uvec, dot_a = 0.0D0, first_a_uvec = uvec3, last_a_uvec = uvec3, & ! complete great-circle & pole_b_uvec = pole_2_uvec, dot_b = 0.0D0, first_b_uvec = uvec5, last_b_uvec = uvec6, & ! tiny arc of another great-circle & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number == 1) EXIT replodding ! with answer in point1_uvec END DO replodding ! j = j1, (j2-1) ! different starting points (with absolute, large-# indexing) correction(1:3) = point1_uvec(1:3) - 0.5D0 * (uvec3(1:3) + uvec4(1:3)) uvec3(1:3) = uvec3(1:3) + correction(1:3) CALL DMake_uvec(uvec3, uvec3) uvec4(1:3) = uvec4(1:3) + correction(1:3) CALL DMake_uvec(uvec4, uvec4) CALL Add_node(uvec3); n1target = numNod CALL Add_node(uvec4); n2target = numNod END IF ! using pre-defined, or new, node-pair for targets IF (dextral) THEN ! element diagonals should reflect right-lateral shear in the past CALL Add_element(n1past, n2past, n1target) CALL Add_element(n2target, n1target, n2past) ELSE ! element diagonals should reflect left-lateral shear in the past CALL Add_element(n2past, n2target, n1past) CALL Add_element(n1target, n1past, n2target) END IF ! dextral? or sinistral shear of element diagonals? n1past = n1target ! The former target becomes the new starting point, n2past = n2target ! to prepare for the next time through this loop. END DO ! k = 1, n_steps END DO ! n = 0, onramp_count(i) END IF ! simple or complex case (without, or with, on-ramp interruptions) ELSE ! previous code for defining end-nodes is incomplete or buggy! WRITE (*, "(' ERROR: Fault ',I6, ' has end-nodes n11=', I6, ', n12=', I6, ', n21=', I6, ', n22=',I6)") i, n11, n12, n21, n22 CALL Pause() STOP END IF ! all 4 bounding nodes are pre-defined END IF ! f_relevant(i) END DO ! i = 1, f_highest WRITE (*, "(' ... DONE!')") !------------------------------------------------------------------------------------ !Output the modified .FEG file: WRITE (*, *) WRITE (*, "(' Now, finally ready to write the modified .FEG file!')") 100 WRITE (*, "(' Enter NEW filename for the modified grid:')") READ (*, "(A)") output_FEG_filename OPEN (UNIT = 9, FILE = TRIM(output_FEG_filename), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This filename may already exist. Try again ...')") CALL Pause() GO TO 100 END IF new_title_line = TRIM(output_FEG_filename) // " = Fault_Corridors applied to: " // TRIM(input_FEG_filename) WRITE (9, "(A)") TRIM(new_title_line) WRITE (9, "(I8, I8, ' 0 1000000 T (numNod, nRealN, nFakeN, n1000, brief)')") numNod, numNod DO i = 1, numNod WRITE (9, "(I6, F11.5, F11.5, ' 0.00E+00 0.00E+00 0.00E+00 0.00E+00')") i, node_Elon(i), node_Nlat(i) END DO ! i = 1, numNod WRITE (9, "(I10, ' (numEl = number of triangular continuum elements)')") numEl DO i = 1, numEl WRITE (9, "(I9, 3I8)") i, nodes(1:3, i) END DO ! i = 1, numEl WRITE (9, "(' 0 (nFl = number of linear fault elements)')") CLOSE(9) !------------------------------------------------------------------------------------ WRITE (*, *) WRITE (*, "(' Job completed. Now, edit the new .FEG manually with OrbWin,')") WRITE (*, "(' until it passes the Perimeter/Area Tests and also the')") WRITE (*, "(' View Gaps/Overlaps test. Lastly, use OrbNumber.')") CALL Pause() CONTAINS SUBROUTINE Add_element(n1, n2, n3) ! Results in a new value of global variable numEl. IMPLICIT NONE INTEGER, INTENT(IN) :: n1, n2, n3 IF (numEl >= numEl_limit) THEN ! not possible! WRITE (*, "(' ERROR: Cannot add another element, due to current numEl_limit = ',I6)") numEl_limit CALL Pause() STOP END IF numEl = numEl + 1 nodes(1, numEl) = n1 nodes(2, numEl) = n2 nodes(3, numEl) = n3 free_node(n1) = .FALSE. free_node(n2) = .FALSE. free_node(n3) = .FALSE. !NOTE that: center(), neighbor(), a_(), delete_this() are not set here, ! as they are unlikely to be needed ! (in THIS particular application: Fault_Corridors). END SUBROUTINE Add_element SUBROUTINE Add_node(uvec) ! Results in a new value of global variable numNod. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8 :: Elon, Nlat IF (numNod >= numNod_limit) THEN ! not possible! WRITE (*, "(' ERROR: Cannot add another node, due to current numNod_limit = ',I6)") numNod_limit CALL Pause() STOP END IF numNod = numNod + 1 node_uvec(1:3, numNod) = uvec(1:3) CALL DUvec_2_LonLat(uvec, Elon, Nlat) node_Elon(numNod) = Elon node_Nlat(numNod) = Nlat free_node(numNod) = .TRUE. ! for the moment... END SUBROUTINE Add_node SUBROUTINE Def_seg_v2 (seg_count, savem) ! Defines fault segments. !(A fault segment is the intersection of one fault trace with one element.) ! IF (savem) THEN segments are recorded; otherwise just counted. ! This version 2 ("_v2") was created 24-27 November 2017, ! after it was discovered that the old Def_seg sometimes ! missed segments near the edges of the grid, in cases where ! adjacent digitized fault-trace points were far apart. !(Such cases were formerly rare, but now occur often because ! faults with the "symmetric_spreading_system" attribute ! are usually digitized using only two end-points.) IMPLICIT NONE CHARACTER(61) :: bar_graph INTEGER, INTENT(OUT) :: seg_count LOGICAL, INTENT(IN) :: savem INTEGER, PARAMETER :: library_size = 1000 ! Increase this if you get error messages. INTEGER :: a, b, e1, e2, ele, i, j, j1, j2, j_bar_graph, jin, jold, jout, k, & & n_hits, n_final_in_library, n_new_in_library, n_old_in_library, n_store, n_test_1, n_test_2, n_zero_length, number, & & temp_int, test_1st, test_2nd, which_neighbor, which_side INTEGER, DIMENSION(library_size) :: library_int ! Contains element# for each segment or sub-segment. LOGICAL :: neighbors, no_cut, overlap REAL*8 :: azimuth1, azimuth2, Elon, length, Nlat, s1, s2, s3 REAL*8, DIMENSION(3) :: omega_uvec, point1_uvec, point2_uvec, pole_a_uvec, pole_b_uvec, tuvec1, tuvec2, tvec, uvec1, uvec2, uvec3, vector REAL*8, DIMENSION(3) :: hit_s REAL*8, DIMENSION(3, 3) :: hit_list REAL*8, DIMENSION(7) :: temp_R8 REAL*8, DIMENSION(7, library_size) :: library_R8 ! 1:3 = starting uvec; 4 = s (of midpoint) within digitization step; 5:7 = final uvec REAL*8, DIMENSION(:, :, :), ALLOCATABLE :: pole_cloud ! (1:3 = uvec, 1:3 = side, 1...numEl) !----------------------------------------------------------------------------------------------- !To save time (in the ***HARD*** case below), create a dataset of poles to all element sides: ALLOCATE ( pole_cloud(3, 3, numEl) ) DO i = 1, numEl DO j = 1, 3 k = 1 + MOD (j, 3) ! has values 2, 3, 1 when j = 1, 2, 3 a = nodes(k, i) ! 1st node along side #j of element# i b = nodes(1 + MOD (k, 3), i) ! 2nd node along side #j of element# i uvec1(1:3) = node_uvec(1:3, a) ! 1st node of side# which_neighbor of element# e1 uvec2(1:3) = node_uvec(1:3, b) ! 2nd node of side# which_neighbor of element# e1 CALL DCross(uvec1, uvec2, tvec) CALL DMake_uvec(tvec, pole_a_uvec) pole_cloud(1:3, j, i) = pole_a_uvec(1:3) END DO ! j = 1, 3 END DO ! i = 1, numEl !----------------------------------------------------------------------------------------------- IF (savem) THEN bar_graph(1:41) = ' Recording fault segments ' ELSE bar_graph(1:41) = ' Counting fault segments ' END IF DO i = 42, 61 bar_graph(i:i) = CHAR(176) END DO PRINT "(' ',A)", bar_graph PRINT "('+',A)", bar_graph(1:41) jold = 0 ! memory of bar graph extent !----------------------------------------------------------------------------------------------- seg_count = 0 ! Initializing the global count of segments (including all faults). n_zero_length = 0 ! initializing this debugging variable, too. !----------------------------------------------------------------------------------------------- DO i = 1, f_highest ! Consider all possible fault trace #s (F0000, F0001, F0002, ... F9999?). f_highest is global trace_loc(3, i) = 0 ! In case no segments trace_loc(4, i) = 0 ! are found in grid. j1 = trace_loc(1, i) ! trace_loc is global. j2 = trace_loc(2, i) IF ((j1 > 0).AND.(j2 > j1).AND.f_relevant(i)) THEN ! This trace represented by at least 2 points in trace dataset, and did (or MIGHT) produce segment(s). !WRITE (*, "(' F',I4)") i !----------------------------------------------------------------------------------------------- !(Re-)initialize the library (which got dynamically re-ALLOCATED on each entry to Def_seg_v2) for this fault: n_old_in_library = 0 ! # of (sub-)segments in library, for this fault trace, but due to previous great-circle arcs n_new_in_library = 0 ! # of (sub-)segments in library, for this fault trace, and due to current great-circle arc !Note that we will NOT zero the library memory, as that wastes too much execution time! !----------------------------------------------------------------------------------------------- DO j = j1, (j2-1) ! consider each digitized point (except the last) as beginning of a great-circle arc e1 = trace_is(j)%element ! Element containing start-point (or 0, for outside-the-grid). e2 = trace_is(j+1)%element ! Element containing end-point (or 0, for outside-the-grid). neighbors = (e1 > 0).AND.(e2 > 0).AND.(e1 /= e2) ! As a minimum condition.... IF (neighbors) THEN ! Now, check more carefully... which_neighbor = 0 ! just initializing before tests... IF (e2 == neighbor(1, e1)) which_neighbor = 1 IF (e2 == neighbor(2, e1)) which_neighbor = 2 IF (e2 == neighbor(3, e1)) which_neighbor = 3 neighbors = (which_neighbor > 0) ELSE ! not neighbors which_neighbor = 0 END IF !----------------------------------------------------------------------------------------------------------- IF ((e1 == e2).AND.(e1 > 0)) THEN ! Easiest case; both start- and end-point are located in the same element. n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = e1 ! == e2 library_R8(1:3, n_store) = trace(1:3, j) ! Start-point of great-circle arc library_R8(4, n_store) = 0.5D0 ! (which will not be used; just for clarity) library_R8(5:7, n_store) = trace(1:3, j+1) ! End-point of great-circle arc !----------------------------------------------------------------------------------------------------------- ELSE IF (neighbors) THEN ! Start- and end-point are in different elements, but they are in neighboring elements. !Find crossing-point where great-circle arc hits side #which_neighbor of element #e1: !"Small" circle 'a' (actually a great-circle) describes the side of element #e1: k = 1 + MOD (which_neighbor, 3) ! has values 2, 3, 1 when which_neighbor = 1, 2, 3 a = nodes(k, e1) ! 1st node along side #which_neighbor of element# e1 b = nodes(1 + MOD (k, 3), e1) ! 2nd node along side #which_neighbor of element# e1 uvec1(1:3) = node_uvec(1:3, a) ! 1st node of side# which_neighbor of element# e1 uvec2(1:3) = node_uvec(1:3, b) ! 2nd node of side# which_neighbor of element# e1 CALL DCross(uvec1, uvec2, tvec) CALL DMake_uvec(tvec, pole_a_uvec) !"Small" circle 'b' (actually a great-circle) describes the digitization step along the fault trace: tuvec1(1:3) = trace(1:3, j) ! Start-point of great-circle arc (one dig-step in fault trace) tuvec2(1:3) = trace(1:3, j+1) ! End-point of great-circle arc (one dig-step in fault trace) CALL DCross(tuvec1, tuvec2, tvec) CALL DMake_uvec(tvec, pole_b_uvec) CALL DCircles_Intersect (pole_a_uvec = pole_a_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = uvec2, & & pole_b_uvec = pole_b_uvec, dot_b = 0.0D0, first_b_uvec = tuvec1, last_b_uvec = tuvec2, & ! input & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number == 1) THEN ! We hope so! !Store the (sub-)segment from the start-point to the crossing point: n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = e1 library_R8(1:3, n_store) = tuvec1(1:3) ! Start-point of great-circle arc (dig step) library_R8(4, n_store) = 0.5D0 * ( DArc(tuvec1, point1_uvec) / DArc(tuvec1, tuvec2) ) ! s value of mid-point (with dig step) library_R8(5:7, n_store) = point1_uvec(1:3) ! Crossing-point !Store the (sub-)segment from the crossing point to the end-point: n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = e2 library_R8(1:3, n_store) = point1_uvec(1:3) ! Crossing-point library_R8(4, n_store) = 0.5D0 * ( (DArc(tuvec1, point1_uvec) / DArc(tuvec1, tuvec2) ) + 1.0D0) ! s value of mid-point (with dig step) library_R8(5:7, n_store) = tuvec2(1:3) ! End-point of great-circle arc (dig step) ELSE ! number /= 1 (This should never happen.) WRITE (*, "(' ERROR: Topological contradiction in Def_seg_v2. Fix program logic.')") CALL Pause() STOP END IF !----------------------------------------------------------------------------------------------------------- ELSE ! ***HARD*** case; start- and end-points are far apart (and/or, either end is outside the grid); ! we must use a SLOW but fail-proof algorithm to find all possible segments: !WRITE (*, "(' includes a ***HARD*** step.')") !"Small" circle 'b' (actually a great-circle) describes the digitization step along the fault trace: tuvec1(1:3) = trace(1:3, j) ! Start-point of great-circle arc (one dig-step in fault trace) tuvec2(1:3) = trace(1:3, j+1) ! End-point of great-circle arc (one dig-step in fault trace) CALL DCross(tuvec1, tuvec2, tvec) CALL DMake_uvec(tvec, pole_b_uvec) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Check for an initial (sub-)segment starting inside element e1: IF (e1 > 0) THEN around_e1: DO which_side = 1, 3 !"Small" circle 'a' (actually a great-circle) describes the side of element #ele: k = 1 + MOD (which_side, 3) ! has values 2, 3, 1 when which_side = 1, 2, 3 a = nodes(k, e1) ! 1st node along side #which_side of element# e1 b = nodes(1 + MOD (k, 3), e1) ! 2nd node along side #which_side of element# e1 uvec1(1:3) = node_uvec(1:3, a) ! 1st node of side# which_neighbor of element# e1 uvec2(1:3) = node_uvec(1:3, b) ! 2nd node of side# which_neighbor of element# e1 pole_a_uvec(1:3) = pole_cloud(1:3, which_side, e1) ! pre-computed at top of this subprogram, for slightly better speed CALL DCircles_Intersect (pole_a_uvec = pole_a_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = uvec2, & & pole_b_uvec = pole_b_uvec, dot_b = 0.0D0, first_b_uvec = tuvec1, last_b_uvec = tuvec2, & ! input & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number > 0) THEN n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = e1 library_R8(1:3, n_store) = tuvec1(1:3) library_R8(4, n_store) = 0.5D0 * DArc(tuvec1, point1_uvec) / DArc(tuvec1, tuvec2) library_R8(5:7, n_store) = point1_uvec(1:3) EXIT around_E1 END IF ! number > 0 END DO around_e1 ! which_side = 1, 3 END IF ! e1 > 0; (sub-)segment starting inside element e1 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Check for any (sub-)segments that entirely cross an element? DO ele = 1, numEl n_hits = 0 ! before we have checked any of the 3 sides... DO which_side = 1, 3 !"Small" circle 'a' (actually a great-circle) describes the side of element #ele: k = 1 + MOD (which_side, 3) ! has values 2, 3, 1 when which_side = 1, 2, 3 a = nodes(k, ele) ! 1st node along side #which_side of element# ele b = nodes(1 + MOD (k, 3), ele) ! 2nd node along side #which_side of element# ele uvec1(1:3) = node_uvec(1:3, a) ! 1st node of side# which_neighbor of element# ele uvec2(1:3) = node_uvec(1:3, b) ! 2nd node of side# which_neighbor of element# ele pole_a_uvec(1:3) = pole_cloud(1:3, which_side, ele) ! pre-computed at top of this subprogram, for slightly better speed CALL DCircles_Intersect (pole_a_uvec = pole_a_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = uvec2, & & pole_b_uvec = pole_b_uvec, dot_b = 0.0D0, first_b_uvec = tuvec1, last_b_uvec = tuvec2, & ! input & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number > 0) THEN n_hits = n_hits + 1 ! probably will not exceed 2; definately no more than 3 hit_list(1:3, n_hits) = point1_uvec(1:3) hit_s(n_hits) = DArc(tuvec1, point1_uvec) / DArc(tuvec1, tuvec2) END IF ! number > 0; this dig-step intersected this element side! END DO ! which_side = 1, 3 !We can create a (sub-)segment for the library if we have 2 hits: IF (n_hits >= 2) THEN ! We hope it is not == 3; that should not happen, since no node can lie on a fault. !Check the ordering of the hit_s values of the (sub-)segment ends, and reverse if needed: IF (hit_s(2) < hit_s(1)) THEN ! reverse, using (3) as temporary storage: hit_list(1:3, 3) = hit_list(1:3, 1) ! #1 copied to temp storage hit_s(3) = hit_s(1) ! #1 copied to temp storage hit_list(1:3, 1) = hit_list(1:3, 2) ! #2 overwrites #1 hit_s(1) = hit_s(2) ! #2 overwrites #1 hit_list(1:3, 2) = hit_list(1:3, 3) ! old #1 retrieved from temp storage, placed in #2 hit_s(2) = hit_s(3) ! old #1 retrieved from temp storage, placed in #2 END IF ! reversing the order of the end-points !Now, record this new (sub-)segment in the library: n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = ele library_R8(1:3, n_store) = hit_list(1:3, 1) library_R8(4, n_store) = 0.5D0 * (hit_s(1) + hit_s(2)) library_R8(5:7, n_store) = hit_list(1:3, 2) END IF END DO ! ele = 1, numEl !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Lastly, check for a possible (sub-)segment ending inside element e2: IF (e2 > 0) THEN around_e2: DO which_side = 1, 3 !"Small" circle 'a' (actually a great-circle) describes the side of element #ele: k = 1 + MOD (which_side, 3) ! has values 2, 3, 1 when which_side = 1, 2, 3 a = nodes(k, e2) ! 1st node along side #which_side of element# e2 b = nodes(1 + MOD (k, 3), e2) ! 2nd node along side #which_side of element# e2 uvec1(1:3) = node_uvec(1:3, a) ! 1st node of side# which_neighbor of element# e2 uvec2(1:3) = node_uvec(1:3, b) ! 2nd node of side# which_neighbor of element# e2 pole_a_uvec(1:3) = pole_cloud(1:3, which_side, e2) ! pre-computed at top of this subprogram, for slightly better speed CALL DCircles_Intersect (pole_a_uvec = pole_a_uvec, dot_a = 0.0D0, first_a_uvec = uvec1, last_a_uvec = uvec2, & & pole_b_uvec = pole_b_uvec, dot_b = 0.0D0, first_b_uvec = tuvec1, last_b_uvec = tuvec2, & ! input & overlap = overlap, number = number, point1_uvec = point1_uvec, point2_uvec = point2_uvec) ! output IF (number > 0) THEN n_new_in_library = n_new_in_library + 1 n_store = n_old_in_library + n_new_in_library IF (n_store > library_size) THEN WRITE (*, "(' ERROR: Increase PARAMETER library_size in Def_seg_v2.')") CALL Pause(); STOP END IF library_int(n_store) = e2 library_R8(1:3, n_store) = point1_uvec(1:3) library_R8(4, n_store) = 0.5D0 * ( ( DArc(tuvec1, point1_uvec) / DArc(tuvec1, tuvec2) ) + 1.0D0 ) library_R8(5:7, n_store) = tuvec2(1:3) EXIT around_e2 END IF ! number > 0 END DO around_e2 ! which_side = 1, 3 END IF ! e2 > 0; (sub-)segment ending inside element e2 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Now, sort the new (sub-)segments in the library, if needed, according to internal coordinate s within the dig-step: IF (n_new_in_library > 1) THEN ! It is possible (likely!) that some 's' values are out of order. DO n_test_1 = (n_old_in_library + 1), (n_old_in_library + n_new_in_library - 1) DO n_test_2 = (n_test_1 + 1), (n_old_in_library + n_new_in_library) IF (library_R8(4, n_test_1) > library_R8(4, n_test_2)) THEN ! This pair needs to be reversed. temp_int = library_int(n_test_1) ! save a copy of n_test_1 value temp_R8(1:7) = library_R8(1:7, n_test_1) ! save copies of n_test_1 values library_int(n_test_1) = library_int(n_test_2) ! n_test_2 value overwrites n_test_1 value library_R8(1:7, n_test_1) = library_R8(1:7, n_test_2) ! n_test_2 values overwrite n_test_1 values library_int(n_test_2) = temp_int ! saved copy moved to n_test_2 library_R8(1:7, n_test_2) = temp_R8(1:7) ! saved copes moved to n_test_2 END IF ! This pair needed to be reversed. END DO ! n_test_2 = (n_test_1 + 1), (n_old_in_library + n_new_in_library) END DO ! n_test_1 = n_test_1 = (n_old_in_library + 1), (n_old_in_library + n_new_in_library - 1) END IF ! n_new_in_library > 1; sorting may be required END IF !end of consideration of this particular digitization-step on one fault n_old_in_library = n_old_in_library + n_new_in_library ! What was "new" is about to be "old", on the next pass through loop. n_new_in_library = 0 ! Re-initializing (this variable, but not the whole library). END DO ! j = j1, (j2-1); considering each digitized point (except the last) as beginning of a great-circle arc !At this point, we have considered the whole fault trace, and filled up the (temporary) library. !----------------------------------------------------------------------------------------------------------- !Compress the library by combining any sub-segments that are in the same element: n_final_in_library = n_old_in_library ! just initializing, before any compression: IF (n_old_in_library > 1) THEN ! compression is theoretically possible test_1st = 1 ! initializing (jagged) loop variable compressing: DO test_2nd = test_1st + 1 e1 = library_int(test_1st) e2 = library_int(test_2nd) IF (e1 == e2) THEN ! compression is required !Expand the lower-numbered sub-segment by incorporating the higher-numbered one: library_R8(5:7, test_1st) = library_R8(5:7, test_2nd) ! moving the end-point of this (sub-)segment. !Note that library_R8(4, ...) has served its purpose; it will no longer be maintained or used. !Decrement the size of the library: n_final_in_library = n_final_in_library - 1 !Shift all library contents that remain: IF (test_2nd <= n_final_in_library) THEN DO j = test_2nd, n_final_in_library library_int(j) = library_int(j+1) library_R8(1:7, j) = library_R8(1:7, j+1) END DO ! j = test_2nd, n_final_in_library END IF ! test_2nd <= n_final_in_library !Note that we do NOT increment test_1st++ in this branch, because it needs to be re-checked! This segment could grow again! ELSE ! no need for compression; but OK to increment jagged loop variable (or exit) test_1st = test_1st + 1 END IF ! compression? or not? IF (test_1st >= n_final_in_library) EXIT compressing ! Note the moving target! END DO compressing END IF ! more than one (sub-)segment in library; compression possible !----------------------------------------------------------------------------------------------------------- !Record(?) these newly-found segments (and finish them?) IF (savem) THEN DO j = 1, n_final_in_library seg_count = seg_count + 1 ele = library_int(j) seg_def(1, seg_count) = i ! fault-trace index (from outermost loop) seg_def(2, seg_count) = ele ! element# IF (trace_loc(3, i) == 0) trace_loc(3, i) = seg_count ! Replace any 0 with index of first segment (but only ONCE). trace_loc(4, i) = seg_count ! This may, or MAY NOT, be the final answer, as loop progresses... seg_end(1:3, 1, seg_count) = library_R8(1:3, j) seg_end(1:3, 2, seg_count) = library_R8(5:7, j) vector(1:3) = library_R8(1:3, j) ! start seg_end_is(1, seg_count)%element = ele CALL Dumb_s123 (node_uvec, ele, vector, s1, s2, s3) seg_end_is(1, seg_count)%s(1) = s1 seg_end_is(1, seg_count)%s(2) = s2 seg_end_is(1, seg_count)%s(3) = s3 vector(1:3) = library_R8(5:7, j) ! end seg_end_is(2, seg_count)%element = ele CALL Dumb_s123 (node_uvec, ele, vector, s1, s2, s3) seg_end_is(2, seg_count)%s(1) = s1 seg_end_is(2, seg_count)%s(2) = s2 seg_end_is(2, seg_count)%s(3) = s3 CALL Finish_seg (seg_count, .FALSE., jin, .FALSE., jout, no_cut) !The above computes jin, jout, and no_cut; !but more important are the side-effects, of computing seg_eta_, seg_kappa_, and seg_u_. !!Write this segment to fault_segments.dig, for graphical-test: !uvec1(1:3) = seg_end(1:3, 1, seg_count) !uvec2(1:3) = seg_end(1:3, 2, seg_count) !length = DArc(uvec1, uvec2) !IF (length > 0.0D0) THEN ! !creating a half-arrowhead to show sense (and end-point): ! azimuth1 = DCompass(from_uvec = uvec1, to_uvec = uvec2) ! in radians, clockwise from North ! azimuth2 = azimuth1 + (135.0D0 * radians_per_degree) ! for half-arrowhead, pointing back at to the right. ! length = 0.25D0 * length ! for half-arrowhead ! CALL DTurn_to(azimuth_radians = azimuth2, base_uvec = uvec2, far_radians = length, & ! inputs ! & omega_uvec = omega_uvec, result_uvec = uvec3) ! WRITE (21, "('F', I4, ', segment ',I6)") i, seg_count ! CALL DUvec_2_LonLat(uvec1, Elon, Nlat) ! WRITE (21, "(1X,SP,ES12.5,',',ES12.5)") Elon, Nlat ! CALL DUvec_2_LonLat(uvec2, Elon, Nlat) ! WRITE (21, "(1X,SP,ES12.5,',',ES12.5)") Elon, Nlat ! CALL DUvec_2_LonLat(uvec3, Elon, Nlat) ! WRITE (21, "(1X,SP,ES12.5,',',ES12.5)") Elon, Nlat ! WRITE (21, "('*** end of line segment ***')") !ELSE ! n_zero_length = n_zero_length + 1 !END IF END DO ! j = 1, n_final_in_library ELSE ! just count the segments, and also mark any faults that failed to produce any segments: IF (n_final_in_library > 0) THEN seg_count = seg_count + n_final_in_library ELSE ! mark this fault, to avoid analyzing it again! f_relevant(i) = .FALSE. END IF END IF ! savem !----------------------------------------------------------------------------------------------------------- ELSE ! This trace was NOT represented by at least 2 points in the external trace dataset: f_relevant(i) = .FALSE. ! (not necessary for Def_seg_v2 logic, but useful later in Fault_Corridors.) END IF ! This trace represented by at least 2 points in external trace dataset, .AND.f_relevant(i). !--------------------------------------------------------------------------------------------------------------- j_bar_graph = (20 * i) / f_highest IF (j_bar_graph > jold) THEN bar_graph(j_bar_graph+41:j_bar_graph+41) = CHAR(219) PRINT "('+',A)", bar_graph(1:j_bar_graph+41) jold = j_bar_graph END IF END DO ! i = 1, f_highest (F0000, F0001, F0002, ... F9999?) DEALLOCATE ( pole_cloud ) ! which would happen on RETURN anyway !IF (n_zero_length > 0) THEN ! WRITE (*, *) ! WRITE (*, "(' Caution: ', I6, ' segments had lengths of zero.')") n_zero_length ! CALL Pause() !END IF END SUBROUTINE Def_seg_v2 SUBROUTINE Dumb_s123 (node_uvec, element, vector, s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). IMPLICIT NONE REAL*8, DIMENSION(:, :), INTENT(IN) :: node_uvec INTEGER, INTENT(IN) :: element REAL*8, DIMENSION(3), INTENT(IN) :: vector REAL*8, INTENT(OUT) :: s1, s2, s3 INTEGER :: i1, i2, i3 REAL*8, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL*8 :: d1, dc, t IF (element == 0) THEN CALL Prevent('element = 0', 1, 'Dumb_s123') END IF i1 = nodes(1, element) i2 = nodes(2, element) i3 = nodes(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = DDot(vector, tv1) IF (dc <= 0.0D0) CALL Prevent('"Internal" vector >= 90 deg. from element', 1, 'Dumb_s123') tv2 = node_uvec(1:3, i1) d1 = DDot(tv2, tv1) t = d1 / dc v1 = t * vector tvi = node_uvec(1:3,i3) - node_uvec(1:3,i2) tvo = v1(1:3) - node_uvec(1:3,i3) CALL DCross(tvi, tvo, tv) s1 = DDot(tv1, tv) * half_R2 / a_(element) tvi = node_uvec(1:3,i1) - node_uvec(1:3,i3) tvo = v1(1:3) - node_uvec(1:3,i1) CALL DCross(tvi, tvo, tv) s2 = DDot(tv1, tv) * half_R2 / a_(element) s3 = 1.00D0 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE Dump_seg(limit, savem, punt) ! Called for debugging information: dump of segments list, ! when number becomes excessive. IMPLICIT NONE INTEGER, INTENT(IN) :: limit LOGICAL, INTENT(IN) :: savem LOGICAL, INTENT(OUT) :: punt INTEGER :: i IF (savem) THEN PRINT "(' Error: Number of fault segments exceeds plausible limit of ',I8)", limit PRINT "(' Probable infinite loop in Def_seg.')" PRINT "(' See file REPORT.txt for dump of segments.')" WRITE (21,"('Error: Number of fault segments exceeds plausible limit of ',I8)") limit WRITE (21,"('(2 * number of elements * number of traces)')") WRITE (21,"('Probable infinite loop in Def_seg.')") WRITE (21,"('Use ORBWEAVER to check for gaps/overlaps in area!')") WRITE (21,"('Dump of segments follows')") WRITE (21,"(/'Trace Element I1 S1 S2 S3 I2 S1 S2 S3')") DO i = 1, limit WRITE (21, "(I5,I8,I5,3F7.4,I5,3F7.4)") & & seg_def(1, i), seg_def(2,i), & & seg_end_is(1,i)%element, seg_end_is(1,i)%s(1:3), & & seg_end_is(2,i)%element, seg_end_is(2,i)%s(1:3) END DO !mmnn = '_BAD' !CALL Write_x_feg(shift = 0) !CALL Write_f_dig(shift = 0) STOP ELSE punt = .TRUE. WRITE (*, "(' Aborting search for fault segments at count of ',I8)") limit ENDIF END SUBROUTINE Dump_seg SUBROUTINE Dump_trace (i) ! i is like "460" in "F0460R" IMPLICIT NONE INTEGER, INTENT(IN) :: i INTEGER :: j, j1, j2 REAL*8 :: lat, lon REAL*8, DIMENSION(3) :: tv WRITE (21, "('Trace of fault ',I6,':')") i j1 = trace_loc(1, i) j2 = trace_loc(2, i) DO j = j1, j2 tv(1:3) = trace(1:3, j) CALL DUvec_2_LonLat(tv, lon, lat) WRITE (21,"(I4,F10.4,F9.3,I6,3F10.5)") i, & & lon, lat, & & trace_is(j)%element, & & trace_is(j)%s(1), trace_is(j)%s(2), trace_is(j)%s(3) END DO WRITE (21,*) END SUBROUTINE Dump_trace SUBROUTINE Find_out (element, inside, outside, & ! inputs & border, coords) ! outputs ! Given a directed arc of a great circle from point "inside" ! (which is inside, or on the border of "element") to point ! "outside" (both positions Cartesian unit vectors), finds ! the last point in contact with "element" and reports its ! Cartesian vector as "border" and its internal coordinates as "coords". IMPLICIT NONE INTEGER, INTENT(IN) :: element REAL*8, DIMENSION(3), INTENT(IN) :: inside, outside REAL*8, DIMENSION(3), INTENT(OUT):: border TYPE(is123), INTENT(OUT) :: coords INTEGER :: i, side, sidea, sideb INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL*8, DIMENSION(3) :: s, far, frac ! NOT vectors REAL*8 :: distance, slope coords%element = element CALL Dumb_s123 (node_uvec, element, inside, s(1), s(2), s(3)) CALL Pull_in(s) CALL Dumb_s123 (node_uvec, element, outside, far(1), far(2), far(3)) DO i = 1, 3 slope = far(i) - s(i) IF (slope < 0.0D0) THEN frac(i) = -s(i) / slope ELSE frac(i) = 9.99D+37 ! "huge"; will not be selected as minimum ENDIF END DO distance = MINVAL(frac) distance = MAX(0.0D0, MIN(1.0D0, distance)) array = MINLOC(frac) side = array(1) sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(side) = 0.0D0 s(sidea) = s(sidea) + distance * (far(sidea) - s(sidea)) s(sideb) = s(sideb) + distance * (far(sideb) - s(sideb)) CALL Pull_in(s) coords%s(1:3) = s CALL Interpolate(coords, border) END SUBROUTINE Find_out SUBROUTINE Find_s1s2s3 ! Determines internal coordinates of all fault trace points. ! Also, computes some speed-&-convenience arrays for general use. IMPLICIT NONE REAL*8, DIMENSION(3) :: tv, tvn, v1 INTEGER :: a, b, back1, back2, element, lastel INTEGER :: i, j, jold, j1, j2, jt, k, l_, m, n REAL*8 :: s1, s2, s3 CHARACTER(61) :: bar_graph LOGICAL :: debug = .FALSE., easy_match ! find element center points DO l_ = 1, numEl v1 = node_uvec(1:3, nodes(1,l_)) + node_uvec(1:3, nodes(2,l_)) + node_uvec(1:3, nodes(3,l_)) CALL DMake_Uvec(v1, tv) center(1:3, l_) = tv END DO ! find neighboring elements on all 3 sides of each neighbor = 0 ! whole array homes: DO l_ = 1, numEl sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) a = nodes(k, l_) ! 1st node along side b = nodes(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO m = 1, numEl IF ((nodes(1, m) == b) .AND. (nodes(2, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((nodes(2, m) == b) .AND. (nodes(3, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((nodes(3, m) == b) .AND. (nodes(1, m) == a)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides END DO homes IF (f_dig_count > 0) THEN bar_graph(1:41) = ' Locating fault traces ' IF (debug) WRITE (21, "('trace_is before correction:')") IF (debug) WRITE (21, "(' point element s1 s2 s3')") DO i = 42, 61 bar_graph(i:i) = CHAR(176) END DO PRINT "(' ',A)", bar_graph PRINT "('+',A)", bar_graph(1:41) jold = 0 DO i = 1, f_dig_count tv = trace(1:3, i) CALL Internal(tv, l_, s1, s2, s3) trace_is(i)%element = l_ trace_is(i)%s(1) = s1 trace_is(i)%s(2) = s2 trace_is(i)%s(3) = s3 IF (debug) WRITE (21, "(I6,I8,3F6.2)") i, trace_is(i)%element, & & (trace_is(i)%s(j3),j3=1,3) j = (20 * i) / f_dig_count IF (j > jold) THEN bar_graph(j+41:j+41) = CHAR(219) PRINT "('+',A)", bar_graph(1:j+41) jold = j END IF END DO ! Scan for cases where a trace wanders from element a briefly ! into element b and then back into a. Change element assignments ! of b -> a. This does not move the trace, but is does prevent ! serious problems with fault segmentation later. IF (debug) WRITE (21, "('Corrections to trace_is:')") IF (debug) WRITE (21, "('trace point old_ele old_s1 old_s2 old_s3 new_ele new_s1 new_s2 new_s3')") DO i = 1, f_highest CALL Unloop_Trace(i, debug) END DO ! i = 1, f_highest END IF ! f_dig_count > 0 IF (f_highest > 0) THEN DO i = 1, f_highest a = trace_loc(1, i) b = trace_loc(2, i) IF ((a > 0).AND.(b > a)) THEN ! We have at least 2 points to work with. n = 0 ! count points of trace inside the grid DO j = a, b IF (trace_is(j)%element > 0) n = n + 1 END DO IF (n < 2) f_2_in(i) = .FALSE. ELSE f_2_in(i) = .FALSE. END IF END DO END IF END SUBROUTINE Find_s1s2s3 SUBROUTINE Finish_seg (segment, cutin, jin, cutout, jout, no_cut) ! completes segment description by deciding u_, eta_, kappa_ ! IF (cutin), uses jin; else, computes it. ! IF (cutout), uses jout; else, computes it. ! IF (jin == jout) sets flag no_cut = T; segment should be junked. ! Refers to global arrays seg_... IMPLICIT NONE LOGICAL, INTENT(IN) :: cutin, cutout INTEGER, INTENT(INOUT) :: jin, jout INTEGER, INTENT(IN) :: segment LOGICAL, INTENT(OUT) :: no_cut INTEGER :: element, n REAL*8, DIMENSION(3) :: c, in, out, tv, tv1, tv2, vn, vs TYPE(is123) :: ist LOGICAL :: debug = .FALSE. ! IF (debug), table of segments goes to REPORT.txt. REAL*8 :: lat1, lat2, lon1, lon2 element = seg_def(2, segment) IF (cutin) THEN in = seg_end(1:3, 1, segment) ELSE ! project backward tv2 = seg_end(1:3, 2, segment) tv1 = seg_end(1:3, 1, segment) CALL Find_out (element = element, & & inside = tv2, & & outside = tv1, & & border = in, & & coords = ist) jin = Which_zero(ist) END IF IF (cutout) THEN out = seg_end(1:3, 2, segment) ELSE ! project forward tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) CALL Find_out (element = element, & & inside = tv1, & & outside = tv2, & & border = out, & & coords = ist) jout = Which_zero(ist) END IF IF (jin == jout) THEN no_cut = .TRUE. ELSE no_cut = .FALSE. n = 6 - jin - jout ! using local node numbering 1-3 IF ((n >=1) .AND. (n <= 3)) THEN seg_u_(segment) = n ELSE PRINT "(' Error: In Finish_seg, jin =',I2,', jout =',I2)", jin, jout PRINT "(' cutin = ',L1,', cutout = ',L1)", cutin, cutout PRINT "(' element =',I5,', segment =',I6)", element, segment STOP ENDIF IF (cutin .AND. cutout) THEN seg_kappa_(segment) = 1.0D0 ELSE tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) seg_kappa_(segment) = DArc(tv1, tv2) / & & DArc(in, out) ENDIF n = nodes(n, element) ! changing n to absolute node number vn = node_uvec(1:3, n) - seg_end(1:3, 1, segment) vs = seg_end(1:3, 2, segment) - seg_end(1:3, 1, segment) CALL DCross(vn, vs, c) tv = center(1:3, element) IF (DDot(c, tv) > 0.0D0) THEN seg_eta_(segment) = +1.0D0 ELSE seg_eta_(segment) = -1.0D0 END IF END IF IF (debug) THEN ! A table is printed in file REPORT.txt, listing all segments. tv = seg_end(1:3, 1, segment) CALL DUvec_2_LonLat (tv, lon1, lat1) tv = seg_end(1:3, 2, segment) CALL DUvec_2_LonLat (tv, lon2, lat2) WRITE (21, & &"(' seg& & def(1) def(2)& & tr_lo(3)& & tr_lo(4)& & seg_end_is(1)& & seg_end_is(2)& & seg_end(1)& & seg_end(2)& & cin jin cout jout no_cut& & u_ eta_ kappa_')") WRITE (21, 9921) & &segment,& &seg_def(1,segment),seg_def(2,segment),& &trace_loc(3,seg_def(1,segment)),& &trace_loc(4,seg_def(1,segment)),& &seg_end_is(1,segment)%element,(seg_end_is(1,segment)%s(j3),j3=1,3),& &seg_end_is(2,segment)%element,(seg_end_is(2,segment)%s(j3),j3=1,3),& &lon1,lat1,& &lon2,lat2,& &cutin,jin,cutout,jout,no_cut,& &seg_u_(segment),seg_eta_(segment),seg_kappa_(segment) 9921 FORMAT (I4,& &I7,I7,& &I9,& &I9,& &I4,3F5.2,& &I4,3F5.2,& &F7.2,F6.2,& &F7.2,F6.2,& &L4,I4,L5,I5,L7,& &I4,F5.0,F7.2) END IF ! debug END SUBROUTINE Finish_seg SUBROUTINE Internal (b_, iele, s1, s2, s3) IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: b_ INTEGER, INTENT(OUT) :: iele REAL*8, INTENT(OUT) :: s1, s2, s3 ! Input is a Cartesian unit vector in the unit sphere. ! Output is element number and s1, s2, s3 internal coordinates of the ! plane triangle finite element. INTEGER, PARAMETER :: memory = 18 ! rather arbitrary; should be big enough(?) INTEGER :: attempts, best_iet_so_far, i, iet, l_ INTEGER, DIMENSION(memory) :: iet_history LOGICAL :: in_loop REAL*8 :: best_minimum_so_far, lon, lat, lowest, & & r2, r2min, s1t, s2t, s3t, s_min_back1, s_min_back2, worst REAL*8, DIMENSION(3):: s_temp, tv ! establish defaults (not found) in case of quick exit s1 = 0.0D0; s2 = 0.0D0; s3 = 0.0D0 !find closest element center to initialize search r2min = 4.01D0 DO l_ = 1, numEl r2 = (b_(1) - center(1, l_))**2 +(b_(2) - center(2, l_))**2 +(b_(3) - center(3, l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 1 radian away, give up. tv = center(1:3, iet) IF (DDot(b_, tv) < 0.540D0) THEN iele = 0 RETURN END IF ! initialize search memory (with impossible numbers) iet_history = 0 ! whole list attempts = 0 best_iet_so_far = 0 ! (intended to be replaced right away) best_minimum_so_far = -999.9 ! (intended to be replaced right away) is_it_here: DO attempts = attempts + 1 IF (attempts > memory) THEN WRITE (*, "(' ERROR: Parameter memory in SUBROUTINE Internal is not large enough.')") WRITE (21, "('ERROR: Parameter memory in SUBROUTINE Internal is not large enough.')") STOP END IF iet_history(attempts) = iet ! first, check for infinite loop! in_loop = .FALSE. ! but that may change, below ... IF (attempts >= 3) THEN DO j = 1, (attempts - 1) IF (iet_history(j) == iet_history(attempts)) in_loop = .TRUE. END DO END IF IF (in_loop) THEN ! in loop; force location to be in the best element so far iet = best_iet_so_far CALL Dumb_s123 (node_uvec, iet, b_, s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL Dumb_s123 (node_uvec, iet, b_, s1t, s2t, s3t) ! - - - - (except, maintain memory of the process) - - - worst = MIN(s1t, s2t, s3t) IF (worst > best_minimum_so_far) THEN ! remember the new best... best_minimum_so_far = worst best_iet_so_far = iet END IF ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE iele = 0 RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE iele = 0 RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE iele = 0 RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iele = iet s1 = s1t s2 = s2t s3 = s3t END SUBROUTINE Internal SUBROUTINE Interpolate (coordinates, v) ! Input is an internal coordinate set. ! Output is v = unit Cartesian (xyz) location vector IMPLICIT NONE TYPE(is123) :: coordinates REAL*8, DIMENSION(3), INTENT(OUT) :: v INTEGER :: i1, i2, i3, iele REAL*8 :: s1, s2, s3 REAL*8, DIMENSION(3) :: vt iele = coordinates%element i1 = nodes(1, iele) i2 = nodes(2, iele) i3 = nodes(3, iele) s1 = coordinates%s(1) s2 = coordinates%s(2) s3 = coordinates%s(3) vt = s1 * node_uvec(1:3, i1) + s2 * node_uvec(1:3, i2) + s3 * node_uvec(1:3, i3) CALL DMake_Uvec(vt, v) END SUBROUTINE Interpolate SUBROUTINE No_Fault_Elements_Allowed() ! Called if nFl > 0 in the .feg file: ! Prints explanatory messages and stops execution. IMPLICIT NONE 101 FORMAT (& &/' This .feg (finite element grid) file contains fault elements!'& &/' Fault elements are not allowed the input .FEG grid, because:'& &/' I. NeoKinema & Restore do not use fault elements.'& &/' 1. Each program has logic to add the compliance of any number of'& &/' faults to the continuum (triangle) elements that contain them.'& &/' 2. Each program has logic to infer the heave-rate and slip-rate of'& &/' such implied fault(s) from the computed strain-rate of the'& &/' triangular continuum element(s).'& &/' 3. Graphics programs* have logic to plot the heave-rates'& &/' of these faults, and also velocity fields with fault '& &/' disontinuities, without the use of fault elements.' & &/' ' & &/' [ *NeoKineMap for NeoKinema; RetroMap for Restore ]') 102 FORMAT (& &/' II. Fault elements cause bad grid topology.'& &/' 1. Fault elements are not read or stored.'& &/' 2. With fault elements ignored, the grids on the two sides of'& &/' each fault are not connected in any way.'& &/' 3. The solution process may fail due to a singular stiffness'& &/' matrix during solution of the linear system.'& &/' 4. Even if the solution does not fail, its physical interpretation'& &/' will be problematical.') 103 FORMAT (& &/' III. Fault elements should be eliminated from the .feg file.'& &/' 1. Re-load this .feg file into OrbWin, select command Faults,'& &/' and use the right mouse button to Heal the fault cuts.'& &/' 2. Use command Adjust to move all nodes off of fault traces.'& &/' 3. Save the edited grid and re-number it with OrbNumber.'& &/' 4. Alternatively, use command Tile in OrbWin to create a'& &/' new .FEG grid with no fault elements.') WRITE (*, 101) CALL Pause() WRITE (*, 102) CALL Pause() WRITE (*, 103) CALL Pause() STOP END SUBROUTINE No_Fault_Elements_Allowed SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause SUBROUTINE Plane_area (folding, look_at_element) IMPLICIT NONE LOGICAL, INTENT(OUT) :: folding !puts areas of plane triangles (below surface) into array a_; !if any is zero or negative, reports folding INTEGER, INTENT(OUT) :: look_at_element !reports element # in which folding occurred (first); !if none, then reports zero. INTEGER :: i1, i2, i3, l_ REAL*8, DIMENSION(3) :: a, b, c, t folding = .FALSE. look_at_element = 0 DO l_ = 1, numEl ! global, like most variables i1 = nodes(1,l_) i2 = nodes(2,l_) i3 = nodes(3,l_) a = node_uvec(1:3,i2) - node_uvec(1:3,i1) b = node_uvec(1:3,i3) - node_uvec(1:3,i2) CALL DCross (a, b, c) t = node_uvec(1:3,i1) + node_uvec(1:3,i2) + node_uvec(1:3,i3) IF (DDot(t, c) > 0.0D0) THEN a_(l_) = DMagnitude(c) * half_R2 ELSE folding = .TRUE. look_at_element = l_ RETURN END IF END DO END SUBROUTINE Plane_area SUBROUTINE Prevent (bad_thing, line, filename) IMPLICIT NONE INTEGER, INTENT(IN) :: line CHARACTER(*), INTENT(IN) :: bad_thing, filename WRITE (*, "(' Error: ', A, ' is illegal in line ', I6 / ' of ', A)") TRIM(bad_thing), line, TRIM(filename) !CALL Traceback() CALL Pause() STOP END SUBROUTINE Prevent SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL*8 :: factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.0D0) THEN highest = MAXVAL(s) medium = 1.000D0 - lowest - highest IF (medium > 0.0D0) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0.0D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.000D0 / (1.000D0 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.000D0 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.000D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0.0D0 s(sideb) = 0.0D0 END IF END IF END SUBROUTINE Pull_in SUBROUTINE Unloop_Trace (i, debug) ! Scan for cases where a trace wanders from element "a" briefly ! into element "b" and then back into "a". Pull offending points ! into element "a", changing both external and internal ! coordinates of these points. This is necessary to prevent ! serious problems with fault segmentation later !(such as a segment which begins and ends on the same element ! side, thus failing to isolate any of the three corner nodes). IMPLICIT NONE INTEGER, INTENT(IN) :: i ! trace index LOGICAL, INTENT(IN) :: debug ! controls WRITE's INTEGER :: back1, back2, element, j, jt, j1, j2, lastel REAL*8 :: s1, s2, s3 REAL*8, DIMENSION(3) :: s, tv j1 = trace_loc(1, i) j2 = trace_loc(2, i) IF ((j1 > 0) .AND. (j2 > j1+1)) THEN ! trace has > 2 points lastel = -1 ! deliberately invalid initial values back1 = -1 back2 = -1 DO j = j1, j2 element = trace_is(j)%element IF (element /= lastel) THEN back2 = back1 ! cascade the memory (losing back2) back1 = lastel IF (element == back2) THEN ! MAY be crossing same side twice IF ((element > 0).AND.(back1 > 0)) THEN !Case where both are real elements; definately needs fixing; !move points from back1 to back2 == element. jt = j - 1 fix_back1: DO IF (trace_is(jt)%element == back1) THEN tv(1:3) = trace(1:3, jt) ! uncorrected location trace_is(jt)%element = element ! assign new element !find (s1, s2, s3) in new element (will not be in bounds): CALL Dumb_s123 (node_uvec = node_uvec, element = element, vector = tv, s1 = s(1), s2 = s(2), s3 = s(3)) CALL Pull_in (s) ! correct s1, s2, s3 to fall on bounds trace_is(jt)%s(1) = s(1) ! record corrected s1, s2, s3 trace_is(jt)%s(2) = s(2) trace_is(jt)%s(3) = s(3) CALL Interpolate (trace_is(jt), tv) ! correct the position trace(1:3, jt) = tv(1:3) ! record corrected position jt = jt - 1 ! prepare to loop ELSE ! we are past the defective loop EXIT fix_back1 END IF END DO fix_back1 back1 = -1 ! reset to a state of innocence back2 = -1 lastel = -1 ELSE IF (element > 0) THEN !Element is real, but back1 was 0 == outside grid. !I think that no action is required; Def_seg can handle this. ELSE !Element = 0 (outside), but back1 was inside. !I think that no action is required; Def_seg can handle this. END IF ! inside/inside, outside/inside, or inside/outside END IF ! element == back2 END IF ! element /= lastel ! prepare to loop lastel = element END DO ! all points in trace END IF ! length of trace exceeds two points END SUBROUTINE Unloop_Trace INTEGER FUNCTION Which_zero(is) ! Returns 1, 2, or 3 to identify which s_i is closest to zero. IMPLICIT NONE TYPE(is123), INTENT(IN) :: is INTEGER, DIMENSION(1) :: j REAL*8, DIMENSION(3) :: sabs sabs(1:3) = ABS(is%s(1:3)) j = MINLOC(sabs) Which_zero = j(1) END FUNCTION Which_zero END PROGRAM Fault_Corridors