PROGRAM Coseismic_Lurch_Rate ! ! A utility program in the NeoKinema family, ! whose purpose is to estimate the horizontal displacement rates (Lurch_Rates) ! at geodetic benchmarks due to occasional slip-events (Coseismic) on a limited set of faults, ! such as the Ionian subduction megathrust, that lie outside a particular NeoKinema F-E grid. ! Specifically, it estimates the (mean coseismic rate) in the following general equation ! regarding horizontal velocities at geodetic benchmarks: ! ! (long-term-average rate) = (interseismic rate) + (mean coseismic rate) ! ! where (interseismic rate) is what is usually available from GPS (or GNSS) data, ! and (long-term-average rate) is the adjusted velocity that should be compared to ! plate-tectonic velocity boundary conditions, and/or long-term geologic fault slip rates. ! ! If all the velocity components coming from this program are negated ! (i.e., flipped by 180 degrees in azimuth), then they will represent ! the perturbation of observed GPS/GNSS horizontal velocities due to ! temporary *interseismic locking* of important faults outside a NeoKinema model area ! (relative to the case where those faults slip continuously at constant rate, ! or relative to the case where those faults do not exist). ! ! Most code in this program is excerpted from NeoKinema_v5.4.f90. ! In fact, if the most active faults (e.g., subduction megathrusts) can be INCLUDED ! in a NeoKinema F-E grid, that will usually be the better way to obtain self-consistent ! estimates of all 3 terms in the equation above. That is because NeoKinema solutions ! can also take advantage of pre-processing of geology fault offset rates by Slippery, ! and also pre-processing of GPS (or GNSS) uncertainties by GPS_Covariance. ! Therefore, this program should only be used if it is impossible/impractical/impolitic ! to include an adjacent subduction zone megathrust in the NeoKinema F-E grid, OR ! if it is impossible/impractical to compute a NeoKinema solution. ! ! By Peter Bird, Professor Emeritus, ! Department of Earth, Planetary, and Space Sciences ! University of California, Los Angeles, California, USA ! February 2026 ! based on suggestions from Dr. Michele M. C. Carafa of INGV, L'Aquila, Italy. !--------------------------------------------------------------------------------- USE DSphere ! Peter Bird's spherical-geometry code in file DSphere.f90; needed for DCircles_Intersect, etc. USE DDislocation ! Peter Bird's code, in DDislocation.f90, to correct observed benchmark velocities by ! adding estimated long-term-average coseismic velocities of rectangular dislocation patches. USE DTriangular ! Peter Bird's code, in DTriangular.f90, to correct benchmark velocities by ! adding estimated long-term-average coseismic velocities of triangular dislocation patches. !NOTE that the initial 'D' in each MODULE/file name above indicates REAL*8 (i.e., DOUBLE PRECISION) ! floating-point variables are used throughout, and expected as arguments in CALLs. !------------------------------------------------------------------------- IMPLICIT NONE ! *SCALAR VARIABLES*: !GPBdeclare CHARACTER*1 :: c CHARACTER*3 :: c3 CHARACTER*4 :: c4 CHARACTER*6 :: c6 CHARACTER*50 :: c50 CHARACTER*60 :: f_dat ! f____.nki = existing input dataset with fault offset rates CHARACTER*60 :: f_dig ! f____.dig = existing input dataset with digitized fault traces (and dip_degrees?) CHARACTER*60 :: gps_file ! list of test points where coseismic_lurch_rate should be computed (in .GPS format) CHARACTER*60 :: output_GPS ="coseismic_lurch_rate.gps" ! name of (new) output file CHARACTER*80 :: c80 CHARACTER*120 :: f_dat_format CHARACTER*120 :: f_dat_titles CHARACTER*134 :: c134, gps_title, gps_format, gps_header INTEGER:: i, i1, ios, j, j1, j2, k INTEGER :: internal_ios, loc_in_c_1, loc_in_c_2, loc_in_c_3 INTEGER :: f_dig_count, f_highest, line, read_status, trace_count INTEGER :: f_dat_dimension, f_dat_count INTEGER :: external_benchmarks, internal_benchmarks INTEGER :: n, n1, n2 LOGICAL :: creeping, in_trace, got_point, got_index, got_dip_degrees LOGICAL :: any_shadow_pseudodata, warned_of_big_offset_sigma, ignore_warning REAL*8, PARAMETER :: s_per_year = 365.25D0 * 24.0D0 * 60.0D0 * 60.0D0 REAL*8, PARAMETER :: normal_dip_degrees = 55.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity REAL*8, PARAMETER :: thrust_dip_degrees = 20.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity REAL*8, PARAMETER :: subduction_dip_degrees = 14.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity REAL*8, PARAMETER :: small_rate_in_mps = 0.1D0 * 0.001D0 / s_per_year ! 0.1 mm/year REAL*8 :: half_R2, R, t, m_per_km = 1000.0D0 REAL*8 :: t1, t2, locking_depth_m_min, locking_depth_m_max REAL*8 :: locking_depth_m_subduction_min, locking_depth_m_subduction_max REAL*8 :: dip_degrees REAL*8 :: z1, z2 REAL*8 :: minimal_bracket_gap_mps, bracket_low, bracket_high REAL*8 :: sigma_offnormal_degrees REAL*8 :: lon, lat, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation REAL*8 :: dipf, rate_mps,fphi, ftheta, lf, argume, duthet, duphi, prebend_trace_azim, postbend_trace_azim REAL*8 :: P1_depth_m, P2_depth_m, P3_depth_m, SIGN, base_azimuth_radians, heave_rate_mps, heave_azimuth_radians REAL*8 :: zBot, zTop, Ds, Ss, Ts, Poisson, test_depth_m, uSouth, uEast, uRadial REAL*8 :: vS_mmpa, wE_mmpa ! *ARRAY VARIABLES*: CHARACTER*1, DIMENSION(:), ALLOCATABLE :: f_sense CHARACTER*50, DIMENSION(:), ALLOCATABLE :: fault_name CHARACTER*80, DIMENSION(:), ALLOCATABLE :: benchmark_name INTEGER, DIMENSION(:, :), ALLOCATABLE :: trace_loc INTEGER, DIMENSION(:), ALLOCATABLE :: which_trace INTEGER, DIMENSION(:), ALLOCATABLE :: external_benchmark_index, internal_benchmark_index LOGICAL(1), DIMENSION(:), ALLOCATABLE :: f_dat_shadow, f_creeping, trace_has_dipslip_rate, trace_has_strikeslip_rate REAL*8, DIMENSION(3) :: omega_uvec, P1_uvec, P2_uvec, P3_uvec, test_uvec, tv, tv1, tv2, tv3, tvo, uvec REAL*8, DIMENSION(3) :: sliprate_mps_x1x2x3 REAL*8, DIMENSION(:, :), ALLOCATABLE :: trace REAL*8, DIMENSION(:), ALLOCATABLE :: f_dip_degrees REAL*8, DIMENSION(:), ALLOCATABLE :: f_dat_dip_degrees, f_locking_depth_m_max, f_locking_depth_m_min REAL*8, DIMENSION(:), ALLOCATABLE :: f_offset_rate, f_offset_rate_bracketed, f_model_offset_rate REAL*8, DIMENSION(:), ALLOCATABLE :: f_offset_rate_sigma_, f_offset_rate_floor, f_offset_rate_ceiling REAL*8, DIMENSION(:, :), ALLOCATABLE :: f_divide REAL*8, DIMENSION(:), ALLOCATABLE :: benchmark_Elon, benchmark_Nlat REAL*8, DIMENSION(:), ALLOCATABLE :: benchmark_theta, benchmark_phi, benchmark_vw, benchmark_unlocked_vw REAL*8, DIMENSION(:, :), ALLOCATABLE :: benchmark_uvec !------------------------------------------------------------------------- ! Introduction to user (on-screen): WRITE (*, "(' PROGRAM Coseismic_Lurch_Rate')") WRITE (*, "(' ')") WRITE (*, "(' A utility program in the NeoKinema family,')") WRITE (*, "(' whose purpose is to estimate the horizontal displacement rates (Lurch_Rates)')") WRITE (*, "(' at geodetic benchmarks due to occasional slip-events (Coseismic) on a limited set of faults,')") WRITE (*, "(' such as the Ionian subduction megathrust, that lie outside a particular NeoKinema F-E grid.')") WRITE (*, "(' Specifically, it estimates the (mean coseismic rate) in the following (vector) equation')") WRITE (*, "(' regarding horizontal velocities at geodetic benchmarks:')") WRITE (*, "(' ')") WRITE (*, "(' (long-term-average rate) = (interseismic rate) + (mean coseismic rate)')") WRITE (*, "(' ')") WRITE (*, "(' where (interseismic rate) is what is usually available from GPS (or GNSS) data,')") WRITE (*, "(' and (long-term-average rate) is the adjusted velocity that should be compared to')") WRITE (*, "(' plate-tectonic velocity boundary conditions, and/or long-term geologic fault slip rates.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' If all the velocity components coming from this program are negated')") WRITE (*, "(' (i.e., flipped by 180 degrees in azimuth), then they will represent')") WRITE (*, "(' the perturbation of observed GPS/GNSS horizontal velocities due to')") WRITE (*, "(' temporary *interseismic locking* of important faults outside a NeoKinema model area')") WRITE (*, "(' (relative to the case where those faults slip continuously at constant rate,')") WRITE (*, "(' or relative to the case where those faults do not exist).')") WRITE (*, "(' ')") WRITE (*, "(' Most code in this program is excerpted from NeoKinema_v5.4.f90.')") WRITE (*, "(' In fact, if the most active faults (e.g., subduction megathrusts) can be INCLUDED')") WRITE (*, "(' in a NeoKinema F-E grid, that will usually be the better way to obtain self-consistent')") WRITE (*, "(' estimates of all 3 terms in the equation above. That is because NeoKinema solutions')") WRITE (*, "(' can also take advantage of pre-processing of geology fault offset rates by Slippery,')") WRITE (*, "(' and also pre-processing of GPS (or GNSS) uncertainties by GPS_Covariance.')") WRITE (*, "(' Therefore, this program should only be used if it is impossible/impractical/impolitic')") WRITE (*, "(' to include an adjacent subduction zone megathrust in the NeoKinema F-E grid, OR')") WRITE (*, "(' if it is impossible/impractical to compute a NeoKinema solution.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird, Professor Emeritus,')") WRITE (*, "(' Department of Earth, Planetary, and Space Sciences')") WRITE (*, "(' University of California, Los Angeles, California, USA')") WRITE (*, "(' February 2026')") WRITE (*, "(' based on suggestions from Dr. Michele M. C. Carafa of INGV, L''Aquila, Italy.')") WRITE (*, "(' ')") CALL Pause() !------------------------------------------------------------------------- ! Read essential scalar parameters: WRITE (*, "(' ------------------------------------------------------------')") WRITE (*, "(' INPUT of essential parameter values:')") ! radius of planet WRITE (*, "(' ')") CALL DPrompt_for_Real("Radius of planet (km)?", 6371.0D0, t) R = t * m_per_km half_R2 = (R**2) / 2.0D0 ! minimum and maximum default locking depths of intraplate faults, in km WRITE (*, "(' ')") CALL DPrompt_for_Real("Default minimum locking depth of intraplate faults (km)?", 0.0D0, t1) CALL DPrompt_for_Real("Default maximum locking depth of intraplate faults (km)?", 15.0D0, t2) locking_depth_m_min = t1 * 1000.0D0 locking_depth_m_max = t2 * 1000.0D0 ! minimum and maximum default locking depths of subduction zones, in km WRITE (*, "(' ')") CALL DPrompt_for_Real("Default minimum locking depth of subduction megathrusts (km)?", 15.0D0, t1) CALL DPrompt_for_Real("Default maximum locking depth of subduction megathrusts (km)?", 35.0D0, t2) locking_depth_m_subduction_min = t1 * 1000.0D0 locking_depth_m_subduction_max = t2 * 1000.0D0 ! sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] WRITE (*, "(' ')") CALL DPrompt_for_Real("Standard deviation (degrees) between [heave vectors of dip-slip faults] & [trace-normal direction]?", 20.0D0, sigma_offnormal_degrees) ! names of fault traces input file (or, "none") WRITE (*, "(' ')") CALL DPrompt_for_String("Filename containing digitized fault traces?", "f____.dig", f_dig) ! names of fault offset-rate input file: WRITE (*, "(' ')") CALL DPrompt_for_String("Fault offset-rate dataset filename?", "f____.nki" , f_dat) ! list of test points where coseismic_lurch_rate should be computed (in .GPS format): WRITE (*, "(' ')") CALL DPrompt_for_String("File with list of test points (in .GPS format)?", "test_points.gps", gps_file) !------------------------------------------------------------------------- ! Read (small) f___.dig file of digitized fault traces describing the megathrust(s): WRITE (*, "(' ')") WRITE (*, "(' ------------------------------------------------------------')") PRINT "(' ','Begin reading input data files')" PRINT "(' ',4X,'Reading fault traces from ',A)", TRIM(f_dig) OPEN (UNIT = 3, FILE = f_dig, STATUS = "OLD", ACTION = "READ",PAD = "YES"); line = 0 ! Skim file and count number of data points f_dig_count = 0 f_highest = 0 trace_count = 0 line = 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) THEN WRITE (*, "('Illegal negative fault number ', I8, ' in line ', I8)") i, line CALL Pause() STOP END IF f_highest = MAX (f_highest, i) ELSE WRITE (*, "(' ERROR: Cannot interpret fault# from the following input line in f*.dig:')") WRITE (*, "(' ',A)") TRIM(c50) CALL Pause() STOP END IF END IF !N.B. In this pass, any "dip_degrees" headers are just ignored. ELSE; EXIT loop_thru; END IF END DO loop_thru CLOSE (UNIT = 3) ! (will be re-read) ! allocate arrays ALLOCATE ( trace (3, f_dig_count) ) trace = 0.0D0 ! whole array ! ALLOCATE ( f_dip_degrees (0:f_highest) ) ! N.B. Inclusion of "0" protects against ! abend if "dip_degrees" precedes "F0001" in the f*.dig file. f_dip_degrees = 0.0D0 ! and will remain 0.0 unless "dip_degrees" is found for this fault trace. ALLOCATE ( trace_loc (4, f_highest) ) trace_loc = 0 ! whole array ! fill arrays OPEN (UNIT = 3, FILE = f_dig, STATUS = "OLD", ACTION = "READ", PAD = "YES") ; line = 0 in_trace = .FALSE. i = 0 line = 0 read_dig: 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 got_point = .TRUE. got_index = .FALSE. READ (c50,*) t1, t2 ! E longitude, N latitude IF (ABS(t2) > 90.001) THEN PRINT "(' Bad latitude ',F10.2,' in line ',I10,' of ',A)", t2, line, TRIM(f_dig) STOP END IF ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN got_point = .FALSE. READ (c50,"(1X,I4)") j2 ! new trace number (NOTE: Any slip-sense bytes are ignored.) trace_count = trace_count + 1 got_index = .TRUE. got_dip_degrees = .FALSE. ! (but, it may follow in the next line?) ELSE IF (c50(1:3) == "***") THEN ! '*** end of line segment ***' got_point = .FALSE. got_index = .FALSE. ELSE ! check for optional "dip_degrees" header" loc_in_c_1 = INDEX(c50, "dip_degrees") IF (loc_in_c_1 > 0) THEN ! found this flag loc_in_c_2 = loc_in_c_1 + 10 ! location of final 's' in "dip_degrees" loc_in_c_3 = loc_in_c_2 + 1 ! potential start of number c50 = c50(loc_in_c_3:50) ! strip out the text flag, leaving the number READ (c50, *, IOSTAT = internal_ios) t IF (internal_ios == 0) THEN dip_degrees = t got_dip_degrees = .TRUE. ELSE WRITE (c4, "(I4)") j2 IF (c4(1:1) == ' ') c4(1:1) = '0' IF (c4(2:2) == ' ') c4(2:2) = '0' IF (c4(3:3) == ' ') c4(3:3) = '0' WRITE (*, "(' ERROR: Number(?) ',A,' following label ""dip_degrees""')") TRIM(c50) WRITE (*, "(' under fault F',A4,' in file ',A)") c4, TRIM(f_dig) WRITE (*, "(' could not be interpreted.')") CALL Pause() STOP END IF END IF END IF ! (lon, lat) line, or F1234 line, or "*** end" line, or "dip_degrees" line ELSE; EXIT read_dig; END IF IF (in_trace) THEN IF (got_point) THEN i = i + 1 CALL DLonLat_2_Uvec(t1, t2, tvo) trace(1:3,i) = tvo trace_loc(2,j1) = i ELSE ! *** end ... in_trace = .FALSE. ENDIF ELSE ! (not in_trace) IF (got_index) THEN j1 = j2 ! new index becomes current index ELSE IF (got_point) THEN i = i + 1 CALL DLonLat_2_Uvec(t1, t2, tv) trace(1:3,i) = tv trace_loc(1,j1) = i in_trace = .TRUE. ELSE ! *** end ... END IF END IF ! (in_trace) IF (got_dip_degrees) THEN IF ((dip_degrees >= 5.0D0).AND.(dip_degrees <= 90.0D0)) THEN f_dip_degrees(j2) = dip_degrees ELSE PRINT "(' ','ERROR: dip_degrees of ',F10.2,' for fault trace ',I4)", dip_degrees, j2 PRINT "(' ','is outside the legal range of 5~90 degrees.')" CALL Pause() STOP END IF END IF END DO read_dig CLOSE (UNIT = 3) ! close f_dig PRINT "(' ',I8,' fault-trace points were read')", f_dig_count !Note: Typically this output will look like "4X,I4,..." but I want to allow for larger integers. !------------------------------------------------------------------------- ! Read (small) f___.nki file of fault offset-rates: ! read f.nki PRINT "(' ',4X,'Reading fault offset-rate data from ',A)", TRIM(f_dat) OPEN (UNIT = 2, FILE = f_dat, STATUS = "OLD", ACTION = "READ", PAD = "YES") ! PAD = YES necessary to interpret short FORMAT string on line 1, ! and to allow reading any lines that lack optional columns #8,9 on right. READ (2, "(A)") f_dat_format ! N.B. This FORMAT **MUST** include entries for optional columns #8,9 on right, ! even if those columns will be blank/missing in every row! READ (2, "(A)") f_dat_titles ! Skim file and count number of data lines, highest fault index; ! also, IF (sigma_offnormal_degrees > 0.0).OR.(S fault(s)), then allow for extra R/L offset-rate components. f_dat_dimension = 0 ! initializing; to be augmented below line = 0 get_real_offset_rate_count: DO READ (2, "(A)", IOSTAT = read_status) c134 IF (read_status == 0) THEN ! read was successful line = line + 1 IF (LEN_TRIM(c134) > 0) THEN ! ignore empty lines (typically, at the end) f_dat_dimension = f_dat_dimension + 1 ! one for each datum line int the input file READ (c134, "(1X,I4,A1)", IOSTAT = read_status) i, c IF (read_status /= 0) THEN WRITE (*, *) WRITE (*, "(' ERROR: Could not read integer F# from bytes #2~#5 of the following line:')") WRITE (*, "(' ',A)") TRIM(c134) CALL Pause() STOP END IF IF (i > f_highest) THEN WRITE (*, *) WRITE (*, "(' ERROR: The following input data record:'/' ',A)") TRIM(c134) WRITE (*, "(' refers to F#',I4,', but the highest F# in the .dig file is ',I4,'.')") i, f_highest CALL Pause() STOP END IF IF (c == 't') c = 'T' IF (c == 'p') c = 'P' IF (c == 'n') c = 'N' IF (c == 'd') c = 'D' IF (c == 's') c = 'S' IF (c == 'S') f_dat_dimension = f_dat_dimension + 1 ! allowing space in arrays for possible strike-slip shadow pseudo-datum IF (((c == 'T').OR.(c == 'P').OR.(c == 'N').OR.(c == 'D')) & &.AND.(sigma_offnormal_degrees > 0.0)) f_dat_dimension = f_dat_dimension + 1 ! allowing space in arrays for possible strike-slip shadow pseudo-datum END IF ELSE EXIT get_real_offset_rate_count END IF END DO get_real_offset_rate_count CLOSE (UNIT = 2) ! (will be re-read twice more below) ! allocate arrays ALLOCATE ( which_trace (f_dat_dimension) ) ALLOCATE ( f_sense (f_dat_dimension) ) ALLOCATE ( f_dat_dip_degrees(f_dat_dimension) ) ALLOCATE ( fault_name (f_dat_dimension) ) ALLOCATE ( f_offset_rate (f_dat_dimension) ) ALLOCATE ( f_offset_rate_bracketed (f_dat_dimension) ) ALLOCATE ( f_model_offset_rate (f_dat_dimension) ) ALLOCATE ( f_offset_rate_sigma_ (f_dat_dimension) ) ALLOCATE ( f_offset_rate_floor (f_dat_dimension) ) ALLOCATE ( f_offset_rate_ceiling (f_dat_dimension) ) ALLOCATE ( f_divide(2, f_dat_dimension) ) ALLOCATE ( f_dat_shadow(f_dat_dimension) ) ALLOCATE ( f_creeping(f_dat_dimension) ) ALLOCATE ( trace_has_dipslip_rate(f_highest) ) ALLOCATE ( trace_has_strikeslip_rate(f_highest) ) ALLOCATE ( f_locking_depth_m_max(f_dat_dimension) ) ALLOCATE ( f_locking_depth_m_min(f_dat_dimension) ) trace_has_dipslip_rate = .FALSE. ! initializing; some value to be changed below trace_has_strikeslip_rate = .FALSE. ! initializing; some value to be changed below OPEN (UNIT = 2, FILE = f_dat, STATUS = "OLD", ACTION = "READ", PAD = "YES") ; line = 0 READ (2,*) ; line = line + 1 READ (2,*) ; line = line + 1 scan_f_dat: DO READ (2, f_dat_format, IOSTAT = ios) c6, c50, t1, t2, creeping, z1, z2; line = line + 1 IF (ios /= 0) EXIT scan_f_dat IF (c6 == ' ') EXIT scan_f_dat ! typically, an accidental blank line at end of file READ (c6, "(1X,I4,1X)") i1 IF (i1 > f_highest) THEN PRINT "(' Illegally high trace index: ',I6)", i1 STOP END IF c = c6(6:6) IF (c == 'd') c = 'D' IF (c == 'l') c = 'L' IF (c == 'n') c = 'N' IF (c == 'p') c = 'P' IF (c == 'r') c = 'R' IF (c == 's') c = 'S' IF (c == 't') c = 'T' IF (.NOT.((c == 'T') .OR. (c == 'N') .OR. (c == 'R') .OR. (c == 'L') & .OR. (c == 'D') .OR. (c == 'P').OR. (c == 'S'))) THEN PRINT "(' Illegal slip sense: ',A1)", c STOP END IF !finally, the main point of this loop, and this reading of the file: IF ((c == 'D').OR.(c == 'N').OR.(c == 'P').OR.(c == 'S').OR.(c == 'T')) THEN IF (trace_has_dipslip_rate(i1)) THEN PRINT "(' ERROR: More than one dip-slip offset rate datum for trace F',I4)", i1 STOP ELSE trace_has_dipslip_rate(i1) = .TRUE. END IF END IF IF ((c == 'L').OR.(c == 'R')) THEN IF (trace_has_strikeslip_rate(i1)) THEN PRINT "(' ERROR: More than one strike-slip offset rate datum for trace F',I4)", i1 STOP ELSE trace_has_strikeslip_rate(i1) = .TRUE. END IF END IF END DO scan_f_dat CLOSE (UNIT = 2) ! close f_dat (to allow re-opening below) ! Second, re-read same file and record its data; add shadow strike-slip pseudo-data if appropriate. ! Also, scan for possible columns #8 and #9 with f_offset_rate_floor and f_offset_rate_ceiling. any_shadow_pseudodata = .FALSE. ! just initializing; may be changed below f_offset_rate_floor = -999.0D0 * 0.001D0 / s_per_year ! whole array; in case no column #8 found. f_offset_rate_ceiling = +999.0D0 * 0.001D0 / s_per_year ! whole array; in case no column #9 found. f_offset_rate_bracketed = .FALSE. ! whole array; may be switched (only to TRUE) during computation. minimal_bracket_gap_mps = 0.0008D0 * 0.001D0 / s_per_year ! leave at least 0.0008 mm/a between target and either bracket! OPEN (UNIT = 2, FILE = f_dat, STATUS = "OLD", ACTION = "READ", PAD = "NO") ; line = 0 ! Note: PAD = "NO" is critical for logic below which detects optional columns #8,9 READ (2,*) ; line = line + 1 ! FORMAT line READ (2,*) ; line = line + 1 ! Headers line f_dat_count = 0 ! cannot use "DO i" because dextral components may be inserted for dip-slip faults! warned_of_big_offset_sigma = .FALSE. read_f_dat: DO READ (2, f_dat_format, IOSTAT = ios) c6, c50, t1, t2, creeping, z1, z2; line = line + 1 IF (ios /= 0) EXIT read_f_dat f_dat_count = f_dat_count + 1 ! use this index for storage READ (c6, "(1X,I4,1X)") i1 IF (i1 > f_highest) THEN PRINT "(' Illegally high trace index: ',I6)", i1 STOP END IF c = c6(6:6) IF (c == 'd') c = 'D' IF (c == 'l') c = 'L' IF (c == 'n') c = 'N' IF (c == 'p') c = 'P' IF (c == 'r') c = 'R' IF (c == 's') c = 'S' IF (c == 't') c = 'T' IF (.NOT.((c == 'T') .OR. (c == 'N') .OR. (c == 'R') .OR. (c == 'L') & .OR. (c == 'D') .OR. (c == 'P').OR. (c == 'S'))) THEN PRINT "(' Illegal slip sense: ',A1)", c STOP END IF which_trace(f_dat_count) = i1 !Note: The present context is creating primary data; shadow "R" data will be added below: dip_degrees = f_dip_degrees(i1) IF (dip_degrees > 0.0D0) THEN ! use this value: f_dat_dip_degrees(f_dat_count) = dip_degrees !test for impossible cases (that would abend with divide-by-zero): IF ((dip_degrees == 90.0D0).AND. & & ((c == 'T').OR.(c == 'P').OR.(c == 'S').OR. & & (c == 'N').OR.(c == 'D'))) THEN WRITE (c4, "(I4)") i1 IF (c4(1:1) == ' ') c4(1:1) = '0' IF (c4(2:2) == ' ') c4(2:2) = '0' IF (c4(3:3) == ' ') c4(3:3) = '0' WRITE (*, "(' ERROR: Fault F',A4,' has 90-degree dip in f*.dig,')") c4 WRITE (*, "(' and has a T, P, S, N, or D component of offset in f*.nki.')") WRITE (*, "(' These choices are incompatible! Please reduce the dip significantly,')") WRITE (*, "(' or just omit the dip_degrees line from f*.dig.')") CALL Pause() STOP END IF ELSE ! use default dip values based on sense: IF (c == 'D') f_dat_dip_degrees(f_dat_count) = normal_dip_degrees IF (c == 'L') f_dat_dip_degrees(f_dat_count) = 90.0D0 IF (c == 'N') f_dat_dip_degrees(f_dat_count) = normal_dip_degrees IF (c == 'P') f_dat_dip_degrees(f_dat_count) = thrust_dip_degrees IF (c == 'R') f_dat_dip_degrees(f_dat_count) = 90.0D0 IF (c == 'S') f_dat_dip_degrees(f_dat_count) = subduction_dip_degrees IF (c == 'T') f_dat_dip_degrees(f_dat_count) = thrust_dip_degrees END IF f_sense(f_dat_count) = c fault_name(f_dat_count) = c50 f_offset_rate(f_dat_count) = t1 * 0.001D0 / s_per_year f_offset_rate_sigma_(f_dat_count) = t2 * 0.001D0 / s_per_year IF (t2 >= 50.01) THEN ! consider warning about excessive offset-rate sigma of more than 50 mm/a IF (.NOT.warned_of_big_offset_sigma) THEN WRITE (*, *) WRITE (*, "(' CAUTION: Some uncertainties (sigmas, standard errors) in offset rates')") WRITE (*, "(' are entered as >50 mm/a. Experience shows that this may lead to')") WRITE (*, "(' ill-conditioned linear systems and spurious volatility in fault')") WRITE (*, "(' slip rates predicted by NeoKinema. I suggest that you edit the')") WRITE (*, "(' f_*.nki file to limit offset-rate uncertainties to <50 mm/a.')") CALL DPrompt_for_Logical('Do you want to IGNORE this warning and continue?',.FALSE.,ignore_warning) IF (.NOT.ignore_warning) THEN CALL Pause() STOP END IF warned_of_big_offset_sigma = .TRUE. END IF END IF f_dat_shadow(f_dat_count) = .FALSE. f_creeping(f_dat_count) = creeping ! Trap cases of missing depth-range columns #6, 7, or reversed min/max locking depths, ! or illogical input of a "locked" fault with zero depth range for its locking !(which might result from missing entries, defaulting both min and max to zero). ! However, note that I do NOT trap illogical cases of a "creeping" fault with ! a finite depth range for locking; this illogic is a prominent feature of the ! WGCEP Fault Models 2.x. Instead, we just omit dislocation correction entirely ! for any fault with creeping == .TRUE.. IF ((.NOT.creeping).AND.(z1 >= 0.0D0).AND.(z2 >= 0.0D0)) THEN IF (z1 >= z2) THEN PRINT "(' Illogical or missing creeping/locking data for fault ',A6)", c6 STOP END IF END IF IF (z1 >= 0.0D0) THEN f_locking_depth_m_min(f_dat_count) = z1 * 1000.0D0 ELSE IF (c == 'S') THEN f_locking_depth_m_min(f_dat_count) = locking_depth_m_subduction_min ELSE f_locking_depth_m_min(f_dat_count) = locking_depth_m_min END IF END IF IF (z2 >= 0.0D0) THEN f_locking_depth_m_max(f_dat_count) = z2 * 1000.0D0 ELSE IF (c == 'S') THEN f_locking_depth_m_max(f_dat_count) = locking_depth_m_subduction_max ELSE f_locking_depth_m_max(f_dat_count) = locking_depth_m_max END IF END IF !See if f_offset_rate_floor and/or f_offset_rate_ceiling are contained in optional columns #8,9: BACKSPACE (2) READ (2, f_dat_format, IOSTAT = ios) c6, c50, t1, t2, creeping, z1, z2, bracket_low IF (ios == 0) THEN f_offset_rate_floor(f_dat_count) = bracket_low * 0.001D0 / s_per_year IF ((f_offset_rate(f_dat_count) - f_offset_rate_floor(f_dat_count)) < minimal_bracket_gap_mps) THEN WRITE (*, "(' ERROR in line ',I6,' of ',A)") line, TRIM(f_dat) WRITE (*, "(' Left bracket (lower limit) on offset rate must be less than rate.')") CALL Pause() STOP END IF END IF BACKSPACE (2) READ (2, f_dat_format, IOSTAT = ios) c6, c50, t1, t2, creeping, z1, z2, bracket_low, bracket_high IF (ios == 0) THEN f_offset_rate_ceiling(f_dat_count) = bracket_high * 0.001D0 / s_per_year IF ((f_offset_rate_ceiling(f_dat_count) - f_offset_rate(f_dat_count)) < minimal_bracket_gap_mps) THEN WRITE (*, "(' ERROR in line ',I6,' of ',A)") line, TRIM(f_dat) WRITE (*, "(' Right bracket (upper limit) on offset rate must be greater than rate.')") CALL Pause() STOP END IF END IF !consider creating a "shadow" strike-slip component datum for dip-slip faults: IF ((c == 'T').OR.(c == 'P').OR.(c == 'N').OR.(c == 'D').OR.(c == 'S')) THEN IF ((sigma_offnormal_degrees > 0.0D0).OR.(c == 'S')) THEN IF (.NOT.trace_has_strikeslip_rate(i1)) THEN ! no user-specified strike-slip rate; OK to add shadow component: any_shadow_pseudodata = .TRUE. ! there is now at least one f_dat_count = f_dat_count + 1 which_trace (f_dat_count) = i1 f_dat_dip_degrees (f_dat_count) = f_dat_dip_degrees(f_dat_count - 1) ! inherited from parent dip-slip component f_sense (f_dat_count) = 'R' ! generic for strike-slip; if it's really L, heave rate will go negative fault_name (f_dat_count) = "shadow datum, which adds strike-slip flexibility" f_offset_rate (f_dat_count) = 0.0D0 ! target rate is always zero for strike-slip component f_offset_rate_sigma_ (f_dat_count) = small_rate_in_mps ! start small; may be increased in proportion to dip-slip f_dat_shadow (f_dat_count) = .TRUE. f_creeping (f_dat_count) = creeping ! inherited from parent dip-slip component f_locking_depth_m_min(f_dat_count) = f_locking_depth_m_min(f_dat_count - 1) ! inherited f_locking_depth_m_max(f_dat_count) = f_locking_depth_m_max(f_dat_count - 1) ! inherited END IF ! this trace does not already have a strike-slip component entered by the user END IF END IF END DO read_f_dat IF (any_shadow_pseudodata) THEN PRINT "(' ',I8,' offset-rate data (& shadow strike-slip components) were registered')", f_dat_count !Note: Typically this output will look like "4X,I4,..." but I want to allow for larger integers. ELSE PRINT "(' ',8X,I4,' offset-rate data were read')", f_dat_count END IF CLOSE (UNIT = 2) ! close f_dat ! Check that all necessary traces where actually read. j = 0 ! number of critical traces missing from f_dig DO i = 1, f_dat_count IF (trace_loc(2,which_trace(i)) == 0) j = j + 1 END DO IF (j > 0) THEN PRINT "(' Error: The following fault traces were missing:')" DO i = 1, f_dat_count IF (trace_loc(2,which_trace(i)) == 0) THEN WRITE (c4,"(I4)") which_trace(i) !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. DO j=1,4 IF (c4(j:j) == ' ') c4(j:j) = '0' END DO PRINT "(' F',A4,A1,' ',A)", & c4, f_sense(i), TRIM(fault_name(i)) END IF END DO STOP END IF ! IF (j > 0) !------------------------------------------------------------------------- ! Read list of GPS (or GNSS) benchmark locations in .GPS file format. ! NOTE that only the (longitude, latitude) and (name) columns will actually be used! ! read .gps file (if any) PRINT "(' ',4X,'Reading geodetic-velocity sites from ',A)", TRIM(gps_file) ! open 1st time to get FORMAT, and count benchmarks OPEN(UNIT = 10, FILE = gps_file, STATUS = "OLD", ACTION = "READ", PAD = "YES") READ (10, "(A)") gps_title READ (10, "(A)") gps_format READ (10, "(A)") gps_header external_benchmarks = 0 benchmarking: DO READ (10, gps_format, IOSTAT = ios) lon, lat IF (ios == 0) THEN external_benchmarks = external_benchmarks + 1 ELSE EXIT benchmarking END IF END DO benchmarking internal_benchmarks = external_benchmarks ! different from NeoKinema, since no F-E grid will be used. CLOSE (10) ALLOCATE ( external_benchmark_index(external_benchmarks) ) ALLOCATE ( internal_benchmark_index(external_benchmarks) ) ALLOCATE ( benchmark_Elon(external_benchmarks) ) ALLOCATE ( benchmark_Nlat(external_benchmarks) ) ALLOCATE ( benchmark_theta(external_benchmarks) ) ALLOCATE ( benchmark_phi(external_benchmarks) ) ALLOCATE ( benchmark_uvec(3, external_benchmarks) ) ALLOCATE ( benchmark_unlocked_vw(2 * external_benchmarks) ) benchmark_unlocked_vw = 0.0D0 ! whole array (both components at each site) ALLOCATE ( benchmark_name(external_benchmarks) ) OPEN(UNIT = 10, FILE = gps_file, STATUS = "OLD", ACTION = "READ", PAD = "YES") READ (10, "(A)") gps_title READ (10, "(A)") gps_format READ (10, "(A)") gps_header !second time through, data is remembered (using internal NeoKinema conventions but external indices) DO i = 1, external_benchmarks READ (10, gps_format, IOSTAT = ios) lon, lat, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation, c80, benchmark_name(i) external_benchmark_index(i) = i ! typically changed when benchmarks outside .feg are dropped; ! but needed to access externally-provided .gp2 file, if any. internal_benchmark_index(i) = i ! To become different in a later code section. [But, NOT in this program!] benchmark_Elon(i) = lon benchmark_Nlat(i) = lat benchmark_theta(i) = (90.0D0 - lat) / degrees_per_radian benchmark_phi(i) = lon / degrees_per_radian CALL DLonLat_2_Uvec(lon, lat, uvec) benchmark_uvec(1:3, i) = uvec(1:3) END DO CLOSE (10) PRINT "(' ',I8,' geodetic-velocity sites were read')", external_benchmarks !Note: Typically this output will look like "4X,I4,..." but I want to allow for larger integers. !------------------------------------------------------------------------- ! MAIN COMPUTATION: ! Correct from observed benchmark_vw --> (estimated) benchmark_unlocked_vw IF (f_dat_count > 0) THEN benchmark_unlocked_vw = 0.0D0 ! whole array; initialize before summing small corrections DO i = 1, f_dat_count !depths of top and bottom of the locked patch: ztop = 0.0D0 ! Was: "ztop = f_locking_depth_m_min(i)" in NeoKinema versions before 2.3 ! Note that logic within subprogram DChange will substitute a small positive value (e.g., 0.01% of zbot). zbot = f_locking_depth_m_max(i) !Note that these were already converted from km-->m, and replaced with default values if negative, in the input section. IF ((.NOT.f_creeping(i)).AND.(zbot > ztop)) THEN ! compute correction only if locked area is positive!!! !Fault dip !Note: The forward direction is the digitizing direction, which is L-->R on the footwall side; ! thus dips measured from the RHS are all 90~175 degrees (but, converted to radians). dipf = Pi - (f_dat_dip_degrees(i) * radians_per_degree) !Slip-rate, in m/s, in the (x1, x2, x3) coordinate system of Aura, in which: ! x1 is along the trace, and the sliprate is positive if sinistral; ! x2 is horizontal and perpendicular to trace, sliprate is positive if divergent; ! x3 is down, and sliprate is positive when the RHS (looking in the forward direction) moves relatively down. !the fault slip rate used for coseismic adjustment is fixed at the geologic slip rate: rate_mps = f_offset_rate(i) ! in m/s; heave or throw depends on f_sense(i) IF (rate_mps /= 0.0D0) THEN ! contribution of this fault trace to geodetic adjustment is non-zero: IF (f_sense(i) == 'L') THEN sliprate_mps_x1x2x3(1) = rate_mps sliprate_mps_x1x2x3(2) = 0.0D0 sliprate_mps_x1x2x3(3) = 0.0D0 ELSE IF (f_sense(i) == 'R') THEN sliprate_mps_x1x2x3(1) = -rate_mps sliprate_mps_x1x2x3(2) = 0.0D0 sliprate_mps_x1x2x3(3) = 0.0D0 ELSE IF (f_sense(i) == 'D') THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = rate_mps sliprate_mps_x1x2x3(3) = rate_mps * DTAN(dipf) ELSE IF (f_sense(i) == 'N') THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = rate_mps / ABS(DTAN(dipf)) sliprate_mps_x1x2x3(3) = -rate_mps ELSE IF (f_sense(i) == 'T') THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = -rate_mps / ABS(DTAN(dipf)) sliprate_mps_x1x2x3(3) = rate_mps ELSE IF ((f_sense(i) == 'P').OR.(f_sense(i) == 'S')) THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = -rate_mps sliprate_mps_x1x2x3(3) = -rate_mps * DTAN(dipf) END IF k = which_trace(i) n1 = trace_loc(1, k) n2 = trace_loc(2, k) - 1 ! first and last digitizing segments associated with this trace; IDed by initial point DO n = n1, n2 ! segment index tv1(1:3) = trace(1:3, n) ! Cartesian unit vector of initial point tv2(1:3) = trace(1:3, n+1) ! and final point of this digitized segment tv(1:3) = (tv1(1:3) + tv2(1:3)) / 2.0D0 ! location of midpoint fphi = DATan2F(tv(2), tv(1)) ! phi, or longitude, in radians ftheta = DATan2F(DSQRT(tv(1)**2 + tv(2)**2), tv(3)) ! theta, or colatitude, in radians lf = 0.5D0 * R * DSQRT((tv2(1)-tv1(1))**2 + (tv2(2)-tv1(2))**2 + (tv2(3)-tv1(3))**2) ! half-length, in m argume = Pi - Get_Azimuth(tv1, tv2) ! radians counterclockwise from +theta (South) DO j = 1, internal_benchmarks CALL DChange (argume = argume, & & btheta = benchmark_theta(j), bphi = benchmark_phi(j), & & dipf = dipf, lf = lf, & & ftheta = ftheta, fphi = fphi, & & radius = R, & & slip = sliprate_mps_x1x2x3, & & wedge = 0.2618D0, & & ztop = ztop, zbot = zbot, & ! inputs & duthet = duthet, duphi = duphi) ! output ! For a typical year with no earthquakes, ! correct the observed benchmark velocities ! by adding the (estimated) missing coseismic part: benchmark_unlocked_vw(2 * j - 1) = benchmark_unlocked_vw(2 * j - 1) + duthet benchmark_unlocked_vw(2 * j ) = benchmark_unlocked_vw(2 * j ) + duphi END DO ! j = 1, internal_benchmarks !-------------------------------------------------------------------------------------------- ! NEW code for NeoKinema_v5.0: Add triangular dislocation patches in gores of non-vertical ! faults which have any significant bend(s) within their traces (i.e., bends between consecutive ! element-sized segments, not just between adjacent digitized points in the input f_[token].DIG file). IF (ABS(Pi_over_2 - dipf) > 0.2618D0) THEN ! fault is significantly non-vertical (small-angle dip is less than 75 degrees). IF (n < n2) THEN ! there will be another segment, so a bend going into the next segment is possible !tv1 and tv2 were already defined by code above; add tv3 = end of next segment: tv3(1:3) = trace(1:3, n+2) ! which is legal since n < n2; so (n+2) <= last digitized point in trace array. prebend_trace_azim = DCompass(from_uvec = tv1, to_uvec = tv2) postbend_trace_azim = DCompass(from_uvec = tv2, to_uvec = tv3) IF (ABS(SIN(postbend_trace_azim - prebend_trace_azim)) > 0.05236D0) THEN ! significant bend, > 3 degrees (guarding against cycle-shifts) P1_uvec = tv2 ! shallow vertex of the TD is at the point of the bend P1_depth_m = 1.0D0 ! burying the shallowest vertex by 1 m below the free surface, to avoid singularities and abends IF (SIN(postbend_trace_azim - prebend_trace_azim) > 0.0D0) THEN ! bend is to the right (seen from above, moving in digitization direction) SIGN = 1.0D0 ! positive sign; must add a triangular dislocation to fill an empty gore. CALL DTurn_To (azimuth_radians = postbend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P2_uvec) ! output CALL DTurn_To (azimuth_radians = prebend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P3_uvec) ! output ELSE ! bend is to the left (seen from above, moving in digitization direction) SIGN = -1.0D0 ! negative sign; must subtract a triangular dislocation to cancel part of one redundant rectangular patch CALL DTurn_To (azimuth_radians = prebend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P2_uvec) ! output CALL DTurn_To (azimuth_radians = postbend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P3_uvec) ! output END IF P2_depth_m = zBot P3_depth_m = zBot base_azimuth_radians = DCompass(from_uvec = P2_uvec, to_uvec = P3_uvec) ! determined in the counterclockwise direction; measured clockwise from N !NOTE that this deep "base" of the triangle is always horizontal, because there is always a uniform zBot locking depth along one NeoKinema fault. IF (f_sense(i) == 'L') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians ! for upper plate, above the triangular dislocation ELSE IF (f_sense(i) == 'R') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians + Pi ! for upper plate, above the triangular dislocation ELSE IF (f_sense(i) == 'D') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians + Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF (f_sense(i) == 'N') THEN heave_rate_mps = rate_mps / ABS(TAN(dipf)) heave_azimuth_radians = base_azimuth_radians + Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF (f_sense(i) == 'T') THEN heave_rate_mps = rate_mps / ABS(TAN(dipf)) heave_azimuth_radians = base_azimuth_radians - Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF ((f_sense(i) == 'P').OR.(f_sense(i) == 'S')) THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians - Pi_over_2 ! for upper plate, above the triangular dislocation END IF !NOTE that the heave(-rate) vector is always horizontal, by definition. !Thus, the angle between "base" and "heave" is always in the horizontal plane. Ds = SIGN * (heave_rate_mps/COS(Pi - dipf)) * SIN(-(heave_azimuth_radians - base_azimuth_radians)) ! where thrusting is defined as positive. Ss = SIGN * heave_rate_mps * COS(heave_azimuth_radians - base_azimuth_radians) ! where left-lateral is defined as positive Ts = 0.0D0 ! no sill or dike intrusion component Poisson = 0.25D0 DO j = 1, internal_benchmarks CALL DThetaPhi_2_Uvec(benchmark_theta(j), benchmark_phi(j), test_uvec) test_depth_m = 0.0D0 ! always compute the result on the free surface, which is were benchmarks are placed CALL Gore(R, P1_uvec, P2_uvec, P3_uvec, P1_depth_m, P2_depth_m, P3_depth_m, & ! inputs & Ds, Ss, Ts, Poisson, & ! inputs & test_uvec, test_depth_m, & ! inputs & uSouth, uEast, uRadial) ! outputs ! For a typical year with no earthquakes, ! correct the observed benchmark velocities ! by adding the (estimated) missing coseismic part: benchmark_unlocked_vw(2 * j - 1) = benchmark_unlocked_vw(2 * j - 1) + uSouth benchmark_unlocked_vw(2 * j ) = benchmark_unlocked_vw(2 * j ) + uEast END DO ! j = 1, internal_benchmarks (2nd such loop, for possible triangular dislocations) END IF ! bend at this point is significant (more than 5 degrees). END IF ! (n < n2); there will be another segment, so a bend is possible. END IF! this fault is non-vertical (small-angle dip is less than 75 degrees). !-------------------------------------------------------------------------------------------- END DO ! n = n1, n2; digitizing step index (IDed by initial point, in array trace_loc) END IF ! rate_mpa /= 0.0 for this trace END IF ! zbot > ztop, so locked area is positive END DO ! i = 1, f_dat_count ELSE ! no faults to correct for! benchmark_unlocked_vw = benchmark_vw ! whole array; no correction performed END IF ! f_dat_count > 0 (fault locking to correct for), or not !------------------------------------------------------------------------- ! Output Coseismic_Lurch_Rates in .GPS file format !GPBhere OPEN (UNIT = 11, FILE = TRIM(output_GPS)) ! unconditional OPEN; overwrites any older results WRITE (11, "(A)") TRIM(output_GPS) WRITE (11, "(A)") TRIM(gps_format) WRITE (11, "(A)") TRIM(gps_header) t = 0.0D0 c3 = "CLR" DO i = 1, external_benchmarks vS_mmpa = benchmark_unlocked_vw(2 * i - 1) * 1000.0D0 * s_per_year wE_mmpa = benchmark_unlocked_vw(2 * i ) * 1000.0D0 * s_per_year vN_mmpa = -vS_mmpa WRITE (11, gps_format) benchmark_Elon(i), benchmark_Nlat(i), wE_mmpa, vN_mmpa, t, t, t, c3, TRIM(benchmark_name(i)) END DO CLOSE (11) !------------------------------------------------------------------------- WRITE (*, "(' =============================================================================')") WRITE (*, "(' ')") WRITE (*, "(' Coseismic_Lurch_Rates has finished normally. See output file ', A)") TRIM(output_GPS) WRITE (*, "(' Suggestion: Plot these vectors with NeoKineMap, and then decide whether to add')") WRITE (*, "(' them to your observed interseismic GPS data.')") CALL Pause() CONTAINS SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause CHARACTER*10 FUNCTION DASCII10(x) ! Returns a right-justified 10-byte (or shorter) version of a ! REAL*8 number, in "human" format, with no more ! than 4 significant digits. IMPLICIT NONE REAL*8, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL*8 :: x_log DOUBLE PRECISION :: y IF (x == 0.0D0) THEN DASCII10=' 0' RETURN ELSE IF (x > 0.0D0) THEN punt = (x >= 999999999.5D0).OR.(x < 0.0000100D0) ELSE ! x < 0.0 punt = (x <= -99999999.5D0).OR.(x > -0.000100D0) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,E10.3)') x !consider possible improvements, from left to right: IF (temp10(3:6) == '.000') THEN ! right-shift 4 spaces over it temp20(7:10) = temp10(7:10) temp20(5:6) = temp10(1:2) temp20(1:4) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(5:6) == '00') THEN ! right-shift 2 spaces over it temp20(7:10) = temp10(7:10) temp20(3:6) = temp10(1:4) temp20(1:2) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(6:6) == '0') THEN ! right-shift 1 space over it temp20(7:10) = temp10(7:10) temp20(2:6) = temp10(1:5) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(8:8) == '+') THEN ! right-shift over + sign in exponent temp20(9:10) = temp10(9:10) temp20(2:8) = temp10(1:7) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(9:9) == '0') THEN ! right-shift over leading 0 in exponent temp20(10:10) = temp10(10:10) temp20(2:9) = temp10(1:8) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF DASCII10 = temp10 ELSE ! can represent without exponential notation x_log = DLOG10(DABS(x)) zeros = DInt_Below(x_log) - 3 y = (10.D0**zeros) * NINT(DABS(x) / (10.D0**zeros)) IF (x < 0.0D0) y = -y WRITE (temp20,"(F20.9)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp20(19:20) == '01') temp20(19:20) = '00' !Find first important byte from right; change 0 -> ' ' k10 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 20, 12, -1 IF (temp20(j:j) == '0') THEN temp20(j:j) = ' ' ELSE k10 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0D0) THEN IF (temp20(10:11) == ' .') temp20(10:11) = '0.' ELSE ! x < 0.0 IF (k10 <= 18) THEN IF (temp20(9:11) == ' -.') temp20(9:11) = '-0.' END IF END IF k1 = k10 - 9 DASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION DASCII10 INTEGER FUNCTION DInt_Below (x) ! Returns integer equal to, or less than, REAL*8 x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER :: i REAL*8 :: y i = DINT(x) IF (x >= 0.D0) THEN DInt_Below = i ELSE ! x < 0. y = 1.0D0 * i IF (y <= x) THEN DInt_Below = i ELSE ! most commonly DInt_Below = i - 1 END IF END IF END FUNCTION DInt_Below SUBROUTINE DPrompt_for_Integer (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with an integer value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text INTEGER, INTENT(IN) :: default INTEGER, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, trial, written LOGICAL :: finished WRITE (suggested,"(I11)") default suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an integer of 9 digits or less.')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Integer SUBROUTINE DPrompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE DPrompt_for_Logical SUBROUTINE DPrompt_for_Real (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a real value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text REAL*8, INTENT(IN) :: default REAL*8, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL*8 :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use ASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = DASCII10(default) !(2) Be sure that the number contains some sign that it is floating-point, not integer: IF (INDEX(suggested, '.') == 0) THEN IF ((INDEX(suggested, 'E') == 0).AND.(INDEX(suggested, 'e') == 0).AND. & & (INDEX(suggested, 'D') == 0).AND.(INDEX(suggested, 'd') == 0)) THEN suggested = ADJUSTL(suggested) point = LEN_TRIM(suggested) + 1 suggested(point:point) = '.' END IF END IF !------------------------------------------------------------------------------------ suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an real number using 11 characters (or less).')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Real SUBROUTINE DPrompt_for_String (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a character-string value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text CHARACTER*(*), INTENT(IN) :: default CHARACTER*(*), INTENT(OUT) :: answer CHARACTER*80 trial INTEGER :: blank_at, default_bytes, leftover, & & prompt_bytes, written prompt_bytes = LEN_TRIM(prompt_text) default_bytes = LEN_TRIM(default) written = 0 leftover = 79 - prompt_bytes - 4 ! unless changed below DO WHILE ((prompt_bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) leftover = 79 - (blank_at - (written+1) + 1) - 4 written = blank_at END DO IF (leftover >= default_bytes) THEN WRITE (*,"(' ',A,' [',A,']')") prompt_text(written+1:prompt_bytes), TRIM(default) ELSE WRITE (*,"(' ',A)") prompt_text(written+1:prompt_bytes) WRITE (*,"(' [',A,']')") TRIM(default) END IF WRITE (*,"(' ?: '\)") READ (*,"(A)") trial IF (LEN_TRIM(trial) == 0) THEN answer = TRIM(default) ELSE answer = TRIM(trial) END IF END SUBROUTINE DPrompt_for_String REAL*8 FUNCTION Get_azimuth (v1, v2) REAL*8, DIMENSION(3), INTENT(IN) :: v1, v2 REAL*8, DIMENSION(3) :: Phi, step, Theta REAL*8 :: del_phi_, del_theta_, lon, lat ! v1 and v2 are Cartesian unit vectors. ! Result is azimuth from v1 to v2, measured at v1, in radians, ! and measured clockwise from North. IF (v1(1) == v2(1)) THEN IF (v1(2) == v2(2)) THEN IF (v1(3) == v2(3)) THEN CALL DUvec_2_LonLat(v1, lon, lat) PRINT "(' Error: Get_azimuth of coincident points at ',F8.3,'E, ',F7.3,'N')", lon, lat WRITE (21,"('Error: Get_azimuth of coincident points at ',F8.3,'E, ',F7.3,'N')") lon, lat STOP END IF END IF END IF step = v2 - v1 CALL DLocal_Phi (v1, Phi) CALL DLocal_Theta(v1, Theta) del_phi_ = DDot(step, Phi) del_theta_ = DDot(step, Theta) Get_azimuth = DATAN2(del_phi_, -del_theta_) END FUNCTION Get_azimuth END PROGRAM Coseismic_Lurch_Rate