PROGRAM Rake_And_Sliprate_Per_Fault ! Utility program in the NeoKinema family, which reads ! three files associated with some (successful) NeoKinema model !(the parameter file p*.nki, the fault-traces file f*.dig, and the output file f*.nko) ! and uses this information to compute the mean rake and mean slip-rate for each fault. ! Attaches two additional columns showing the dip used, and whether the fault creeps. ! Compatible with NeoKinema versions 4+. ! By Peter Bird, UCLA, 2011-2012; last updated 2019.04.25. USE DFLIB, ARCQQ => ARC ! provided with Digital Visual Fortran: ! Using GETFILEINFOQQ, which provides names of files ! matching spec.s like "v*.out". Helps user select correct file. ! If no substitute is available on your system when you compile, ! just omit SUBROUTINE File_List (and any CALLs to it). ! Also, using BEEPQQ to sound PC speaker when each task is done; ! again, this can simply be omitted if there is no substitute. ! However, not using ARC, because I have my own Arc; so I am ! renaming their ARC to ARCQQ to avoid conflicts. IMPLICIT NONE CHARACTER(1), PARAMETER :: TAB = CHAR(9) ! special "HT" tab character in ASCII sequence CHARACTER*1 :: component1_c1, component2_c1, creeps_c1, dipSlip_sense_c1, line1_c1, line2_c1, & & sense_c1, strikeSlip_sense_c1, total_c1 CHARACTER*4 :: c4, fault_dip_c4 CHARACTER*6 :: Fcode_c6, test_Fcode_c6, rake_c6 CHARACTER*7 :: slipRate_c7 CHARACTER*50 :: fault_name = ' ', test_fault_name = ' ' CHARACTER*79 :: f_nko_file = ' ', f_nko_FORMAT = ' ', extended_FORMAT = ' ', f_nko_headers = ' ', output_file = ' ' CHARACTER*80 :: line, line2 INTEGER :: components, i, ios, ios1, ios2, f_highest, internal_ios, line1_count, line2_count, & & read_loc, string_loc, total_count LOGICAL :: creeps, got_dip_degrees, got_i, reached_EOF, read_strikeSlip, read_dipSlip, this_fault_creeps REAL*8, PARAMETER :: radians_per_degree = 0.01745329252D0 REAL*8, PARAMETER :: normal_dip_degrees = 55.0D0 ! defaults, for NKv2 models, or NKv3/4 models with no "dip_degrees" data. REAL*8, PARAMETER :: thrust_dip_degrees = 20.0D0 ! defaults, for NKv2 models, or NKv3/4 models with no "dip_degrees" data. REAL*8, PARAMETER :: subduction_dip_degrees = 14.0D0 ! defaults, for NKv2 models, or NKv3/4 models with no "dip_degrees" data. REAL*8 :: dextral_rate, dip_degrees, input, normal_downDip_rate, output, rake, sigma, slipRate, t REAL*8, DIMENSION(:), ALLOCATABLE :: f_dip_degrees !--- begin variables needed by, or to be defined by "Get_Parameters" ------ CHARACTER*2 :: reference_plate_c2 CHARACTER*132 :: f_dat = ' ', f_dig = ' ', & & gps_file = ' ', gp2_file = ' ', & ! << f_dig & token are critical & parameter_file = ' ', s_dat = ' ', token = ' ', & & x_feg = ' ', x_bcs = ' ' INTEGER :: n_refine, stress_interpolation_method LOGICAL :: conservative_geodetic_adjustment, faults_give_sigma_1h, floating_frame, got_parameters REAL, PARAMETER :: m_per_km = 1000.0D0 REAL*8 :: A0, geodesy_weight, L0, locking_depth_m_max, locking_depth_m_min, & & locking_depth_m_subduction_max, locking_depth_m_subduction_min, & & mu_, R, sigma_offnormal_degrees, xi_ !--- end variables needed by, or to be defined by "Get_Parameters" ------ !================================================================================== WRITE (*, "(' PROGRAM Rake_And_Sliprate_Per_Fault')") WRITE (*, "(' ')") WRITE (*, "(' Utility program in the NeoKinema family, which reads')") WRITE (*, "(' three files associated with some (successful) NeoKinema model')") WRITE (*, "(' (the parameter file p*.nki,')") WRITE (*, "(' the fault-traces file f*.dig, and ')") WRITE (*, "(' the output file f*.nko)')") WRITE (*, "(' and uses this information to compute the mean rake and')") WRITE (*, "(' mean slip-rate for each fault.')") WRITE (*, "(' Attaches two additional columns showing the fault dip used,')") WRITE (*, "(' and whether the fault creeps.')") WRITE (*, "(' Compatible with NeoKinema versions 4+.')") WRITE (*, "(' By Peter Bird, UCLA, 2011-2012; last update 2019.04.25.')") CALL Get_Parameters() ! prompts for p*.nki file, and finds names of parameter_file, f_dig, and token. !This method GUARANTEES that the f_dig and the f_nko_file are compatible with parameter file. !Note that token may include a leading '_' or not; that is the user's option. f_nko_file = "f"//TRIM(token)//".nko" !read f_dig once to determine highest integer in any F1234 header: OPEN (UNIT = 1, FILE = TRIM(f_dig), STATUS = "OLD", IOSTAT = ios) f_highest = 0 f_dig_scan: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT f_dig_scan IF (line(1:1) == 'F') THEN c4 = line(2:5) READ (c4, *) i f_highest = MAX(f_highest, i) END IF END DO f_dig_scan CLOSE (1) !prepare to remember any dip_degrees given in f_dig: ALLOCATE ( f_dip_degrees(f_highest) ) f_dip_degrees = 0.0D0 ! whole array; only changed if a number is provided. !read f_dig again to memorize dips: OPEN (UNIT = 1, FILE = TRIM(f_dig), STATUS = "OLD", IOSTAT = ios) got_i = .FALSE. got_dip_degrees = .FALSE. f_dig_rescan: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT f_dig_rescan IF (line(1:1) == 'F') THEN c4 = line(2:5) READ (c4, *, IOSTAT = internal_ios) i IF (internal_ios == 0) THEN ! typically, got_i = .TRUE. got_dip_degrees = .FALSE. ! for now, but that could change ELSE WRITE (*, "(' ERROR: Bad fault index number in line:'/' ',A)") TRIM(line) CALL Pause() STOP ! not expected to happen, as NeoKinema already vetted this file format END IF ELSE IF ((line(1:2) == " +").OR.(line(1:2) == " -")) THEN ! a (lon, lat) pair; no action ELSE IF (line(1:3) == "***") THEN ! end of this fault; decide whether to save dip? IF (got_i .AND. got_dip_degrees) THEN f_dip_degrees(i) = dip_degrees END IF got_i = .FALSE. got_dip_degrees = .FALSE. ! must start over again ELSE ! perhaps this is a line with "dip_degrees" in it? string_loc = INDEX(line, "dip_degrees") IF (string_loc > 0) THEN ! this string was found read_loc = string_loc + 11 line2 = line(read_loc:LEN_TRIM(line)) READ (line2, *, IOSTAT = internal_ios) t IF (internal_ios == 0) THEN dip_degrees = t got_dip_degrees = .TRUE. ELSE WRITE (*, "(' ERROR: Unreadable number in following line:'/' ',A)") TRIM(line) CALL Pause() STOP ! not expected, as NeoKinema already vetted the format of this file. END IF END IF END IF END DO f_dig_rescan CLOSE (1) OPEN (UNIT = 1, FILE = TRIM(f_nko_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Necessary file ',A/' was not found in the current folder.')") TRIM(f_nko_file) CALL Pause() STOP END IF READ (1, "(A)") f_nko_FORMAT extended_FORMAT = ' '//TRIM(f_nko_FORMAT) READ (1, "(A)") f_nko_headers output_file = TRIM(token)//"_rake_and_sliprate_per_fault.txt" IF (output_file(1:1) == '_') output_file = output_file(2:LEN_TRIM(output_file)) OPEN (UNIT = 2, FILE = TRIM(output_file)) ! unconditional OPEN; overwrites any existing WRITE (2, "(A)") "Fnnnn"//tab//"Fault name"//tab//"Rake"//tab//"Slip-rate"//tab//"Dip"//tab//"Creeps?" WRITE (*, "(' ')") WRITE (*, "(' ')") !work is done inside an indefinite loop (terminated by EOF) reached_EOF = .FALSE. scanning: DO !initialize (mostly, a state of ignorance) for each new fault: got_i = .FALSE. read_strikeSlip = .FALSE. read_dipSlip = .FALSE. components = 0 !get first (perhaps only?) line related to a new fault: READ (1, f_nko_FORMAT, IOSTAT = ios) Fcode_c6, fault_name, input, sigma, creeps, output IF (ios /= 0) THEN IF (ios == -1) THEN reached_EOF = .TRUE. EXIT scanning ELSE WRITE (*, "(' ERROR: Could not interpret input line which follows: ')") WRITE (*, extended_FORMAT) Fcode_c6, fault_name, input, sigma, creeps, output CALL Pause() STOP END IF END IF c4 = Fcode_c6(2:5) READ (c4, *, IOSTAT = internal_ios) i IF (internal_ios == 0) THEN got_i = .TRUE. ELSE WRITE (*, "(' ERROR: Bad fault index number in line:'/' ',A)") TRIM(line) CALL Pause() STOP ! not expected to happen, as NeoKinema already vetted this file format END IF dip_degrees = f_dip_degrees(i) ! BUT, may be zero, meaning UNKNOWN. sense_c1 = Fcode_c6(6:6) component1_c1 = sense_c1 ! saved for final dip-output check components = components + 1 IF ((sense_c1 == 'L').OR.(sense_c1 == 'R')) THEN read_strikeSlip = .TRUE. this_fault_creeps = creeps strikeSlip_sense_c1 = sense_c1 IF (strikeSlip_sense_c1 == 'R') THEN dextral_rate = output ELSE ! must be 'L' dextral_rate = -output END IF !Caution: Do NOT overwrite the fault dip if it was previously set by a dip-slip data line! IF (.NOT.read_dipSlip) THEN IF (dip_degrees == 0.0D0) dip_degrees = 90.0D0 ! default for pure strike-slip END IF ELSE IF ((sense_c1 == 'D').OR.(sense_c1 == 'N').OR.(sense_c1 == 'P').OR.(sense_c1 == 'S').OR.(sense_c1 == 'T')) THEN read_dipSlip = .TRUE. this_fault_creeps = creeps dipSlip_sense_c1 = sense_c1 IF (dipSlip_sense_c1 == 'D') THEN IF (dip_degrees == 0.0D0) dip_degrees = normal_dip_degrees ! default value normal_downDip_rate = +output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'N') THEN IF (dip_degrees == 0.0D0) dip_degrees = normal_dip_degrees ! default value normal_downDip_rate = +output / DSIN(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'P') THEN IF (dip_degrees == 0.0D0) dip_degrees = thrust_dip_degrees ! default value normal_downDip_rate = -output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'S') THEN IF (dip_degrees == 0.0D0) dip_degrees = subduction_dip_degrees ! default value normal_downDip_rate = -output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'T') THEN IF (dip_degrees == 0.0D0) dip_degrees = thrust_dip_degrees ! default value normal_downDip_rate = -output / DSIN(dip_degrees * radians_per_degree) END IF ELSE ! bad sense bit WRITE (*, "(' ERROR: Bad slip-sense bit in the following line: ')") WRITE (*, extended_FORMAT) Fcode_c6, fault_name, input, sigma, creeps, output CALL Pause() STOP END IF !try to read a second line relating to the same fault(?): READ (1, f_nko_FORMAT, IOSTAT = ios) test_Fcode_c6, test_fault_name, input, sigma, creeps, output IF (ios /= 0) THEN IF (ios == -1) THEN reached_EOF = .TRUE. ! EXIT scanning will not occur yet; we have to make use of data already read in previous line ELSE WRITE (*, "(' ERROR: Could not interpret input line which follows: ')") WRITE (*, extended_FORMAT) Fcode_c6, fault_name, input, sigma, creeps, output CALL Pause() STOP END IF ELSE ! ios == 0; second line was read successfully: IF (test_Fcode_c6(1:5) /= Fcode_c6(1:5)) THEN BACKSPACE (1) ! save this data for the next pass through "scanning" loop ELSE ! new line matches previous line!!! !Note that i and dip_degrees do not require any update. dip_degrees = f_dip_degrees(i) ! BUT, may be zero, meaning UNKNOWN. sense_c1 = test_Fcode_c6(6:6) component2_c1 = sense_c1 ! saved for final dip-output check components = components + 1 IF ((sense_c1 == 'L').OR.(sense_c1 == 'R')) THEN read_strikeSlip = .TRUE. strikeSlip_sense_c1 = sense_c1 IF (strikeSlip_sense_c1 == 'R') THEN dextral_rate = output ELSE ! must be 'L' dextral_rate = -output END IF !Caution: Do NOT overwrite the fault dip if it was previously set by a dip-slip data line! IF (.NOT.read_dipSlip) THEN IF (dip_degrees == 0.0D0) dip_degrees = 90.0D0 ! default for pure strike-slip END IF ELSE IF ((sense_c1 == 'D').OR.(sense_c1 == 'N').OR.(sense_c1 == 'P').OR.(sense_c1 == 'S').OR.(sense_c1 == 'T')) THEN read_dipSlip = .TRUE. dipSlip_sense_c1 = sense_c1 IF (dipSlip_sense_c1 == 'D') THEN IF (dip_degrees == 0.0D0) dip_degrees = normal_dip_degrees normal_downDip_rate = +output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'N') THEN IF (dip_degrees == 0.0D0) dip_degrees = normal_dip_degrees normal_downDip_rate = +output / DSIN(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'P') THEN IF (dip_degrees == 0.0D0) dip_degrees = thrust_dip_degrees normal_downDip_rate = -output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'S') THEN IF (dip_degrees == 0.0D0) dip_degrees = subduction_dip_degrees normal_downDip_rate = -output / DCOS(dip_degrees * radians_per_degree) ELSE IF (dipSlip_sense_c1 == 'T') THEN IF (dip_degrees == 0.0D0) dip_degrees = thrust_dip_degrees normal_downDip_rate = -output / DSIN(dip_degrees * radians_per_degree) END IF ELSE ! bad sense bit WRITE (*, "(' ERROR: Bad slip-sense bit in the following line: ')") WRITE (*, extended_FORMAT) Fcode_c6, fault_name, input, sigma, creeps, output CALL Pause() STOP END IF ! OK sense bit, or bad one END IF ! new line mismatches previous, or matches previous END IF ! second line was read successfully IF (read_dipSlip .AND. read_strikeSlip) THEN slipRate = DSQRT(dextral_rate**2 + normal_downDip_rate**2) rake = 57.2957795130D0 * DATan2(-normal_downDip_rate, -dextral_rate) !attempt to combine count of geologic offset-rate data from the two adjacent lines, if possible: IF ((fault_name(43:43) == '(').AND.(test_fault_name(43:43) == '(').AND.& & (fault_name(45:50) == " data)").AND.(test_fault_name(45:50) == " data)")) THEN line1_c1 = fault_name(44:44) line2_c1 = test_fault_name(44:44) READ (line1_c1, "(I1)", IOSTAT = ios1) line1_count IF (ios1 == 0) THEN READ (line2_c1, "(I1)", IOSTAT = ios2) line2_count IF (ios2 == 0) THEN total_count = line1_count + line2_count IF (total_count < 10) THEN WRITE (total_c1, "(I1)") total_count fault_name(44:44) = total_c1 ELSE fault_name(44:44) = 'n' END IF IF (total_count > 0) fault_name(42:42) = '>' ELSE fault_name(44:44) = 'n' END IF ELSE fault_name(44:44) = 'n' END IF END IF ELSE IF (read_dipSlip) THEN slipRate = ABS(normal_downDip_rate) IF (normal_downDip_rate >= 0.0D0) THEN rake = -90.0D0 ELSE rake = 90.0D0 END IF ELSE IF (read_strikeSlip) THEN slipRate = ABS(dextral_rate) IF (dextral_rate >= 0.0D0) THEN rake = 180.0D0 ELSE rake = 0.0D0 END IF ELSE ! not expected ever to happen WRITE (*, "(' ERROR: Logical flow in PROGRAM Rake_And_Sliprate_Per_Fault.')") CALL Pause() STOP END IF !If dip was not specified with "dip_degrees" in f*.dig, and if fault had 2 components of slip, !then output the built-in dip that was used in computing the dip-slip rate: IF (f_dip_degrees(i) == 0.0D0) THEN ! we must output a built-in default dip; which one? IF (components == 2) THEN ! correct non-vertical dip may have been overwritten (if in component 1; restore): IF ((component1_c1 == 'D').OR.(component2_c1 == 'D').OR. & & (component1_c1 == 'N').OR.(component2_c1 == 'N')) THEN dip_degrees = normal_dip_degrees ELSE IF ((component1_c1 == 'P').OR.(component2_c1 == 'P').OR. & & (component1_c1 == 'T').OR.(component2_c1 == 'T')) THEN dip_degrees = thrust_dip_degrees ELSE IF ((component1_c1 == 'S').OR.(component2_c1 == 'S')) THEN dip_degrees = subduction_dip_degrees ELSE ! should not happen; 2-component slip-vector cannot be all strike-slip! dip_degrees = 90.0D0 END IF END IF END IF WRITE (rake_c6, "(F6.1)") rake WRITE (slipRate_c7, "(F7.3)") slipRate WRITE (fault_dip_c4, "(F4.1)") dip_degrees WRITE (creeps_c1, "(L1)") this_fault_creeps WRITE (2, "(A)") Fcode_c6(1:5)//TAB//TRIM(fault_name)//TAB//TRIM(ADJUSTL(rake_c6))//TAB//& & TRIM(ADJUSTL(slipRate_c7))//TAB//fault_dip_c4//TAB//creeps_c1 IF (reached_EOF) EXIT scanning END DO scanning CLOSE (UNIT = 1, DISP = "KEEP") ! input f_nko file CLOSE (UNIT = 2, DISP = "KEEP") ! new output file WRITE (*, "(' ')") WRITE (*, "(' Job completed. See new ', A)") TRIM(output_file) WRITE (*, "(' This tab-delimited ASCII table can be read by any spreadsheet program,')") WRITE (*, "(' in which it can be formatted for a more elegant display.')") CALL Pause() CONTAINS SUBROUTINE Bad_Parameters() !prints admonition screen, calls Pause(), and STOP's. IMPLICIT NONE WRITE (*, *) WRITE (*, "(' ==========================================================================')") WRITE (*, "(' ERROR DURING READING OF PARAMETER FILE')") WRITE (*, "(' ')") WRITE (*, "(' Please review the proper format, which is listed in the introductory')") WRITE (*, "(' comment lines of NeoKinema.f90.')") WRITE (*, "(' ')") WRITE (*, "(' This program is NOT smart enough to read NeoKinema v1~v3 format.')") WRITE (*, "(' NeoKinema v4+ input parameter files are required.')") WRITE (*, "(' ')") WRITE (*, "(' After you have read this message, this program will stop.')") WRITE (*, "(' Please correct the input parameter file and re-start it.')") WRITE (*, "(' ==========================================================================')") Call Pause() STOP END SUBROUTINE Bad_Parameters SUBROUTINE Could_Not_Find_File(pathfilename) !prevents multiple duplications of this very simple code: IMPLICIT NONE CHARACTER*(*), INTENT(IN), OPTIONAL :: pathfilename IF (PRESENT(pathfilename)) THEN WRITE (*, "(' ERROR: Could not find a necessary input file:'/' ',A)") TRIM(pathfilename) ELSE WRITE (*, "(' ERROR: Could not find a necessary input file')") END IF CALL Pause() STOP END SUBROUTINE Could_Not_Find_File ! SUBROUTINE File_List( file_type, & ! & suggested_file) ! ! Reports a list (on default device) of filenames of the type requested. ! ! ! ! Usage of CHARACTER*(*), INTENT(INOUT) :: suggested_file ! ! depends on how many files (of specified type) are ! ! found in the current using_path directory: ! ! * If none are found, suggested_file is unchanged (it may ! ! be a correct file name in some other directory). ! ! * If one file is found, suggested_file is changed to its name. ! ! * If multiple files are found: ! ! -if suggested_file is one of them, it is unchanged. ! ! -if suggested_file is not one, it is changed to ' '. ! ! ! ! Uses GETFILEINFOQQ of module DFLIB.F90 ! ! (DIGITAL Visual Fortran 5.0). ! IMPLICIT NONE ! CHARACTER*(*), INTENT(IN) :: file_type ! CHARACTER*(*), INTENT(INOUT) :: suggested_file ! CHARACTER*1 :: first_letter ! CHARACTER*70 :: line = ' ', old_name ! CHARACTER*132 :: string0, string1, string2 ! CHARACTER*255 :: files ! INTEGER :: count, full_to, handle, old_result, result ! LOGICAL :: duplicate, matched !! TYPE file$info ! this type as defined in DFLIB.F90 !! INTEGER(4) creation !! INTEGER(4) lastwrite !! INTEGER(4) lastaccess !! INTEGER(4) length !! INTEGER(4) permit !! CHARACTER(255) name !! END TYPE file$info ! TYPE (FILE$INFO) info ! this type as defined in DFLIB.F90 ! !10 count = 0 ! matched = .FALSE. ! until we find a file == suggested_file ! IF (file_type == "*.*") THEN ! WRITE (*,"(/' Here are all the files in the input directory:')") ! files = '*.*' ! ! ELSE IF (file_type == "*.eqc") THEN ! WRITE (*,"(/' The following appear to be EarthQuake catalog (.eqc) files:')") ! files = '*.EQC' ! ! ELSE IF (file_type == "*.dig") THEN ! WRITE (*,"(/' The following appear to be basemap (.dig) files:')") ! files = '*.DIG' ! ! ELSE IF (file_type == "*.feg") THEN ! WRITE (*,"(/' The following appear to be FE grid (.feg) files:')") ! files = '*.FEG' ! ELSE IF (file_type == "*.gps") THEN ! WRITE (*,"(/' The following appear to be geodetic-velocity (.gps) files:')") ! files = '*.GPS' ! ELSE IF (file_type == "*.grd") THEN ! WRITE (*,"(/' The following appear to be gridded data (.grd) files:')") ! files = '*.GRD' ! ELSE IF (file_type == "e*.nko") THEN ! WRITE (*,"(/' The following appear to be continuum strain-rate output (e*.nko) files:')") ! files = '*.NKO' ! ELSE IF (file_type == "f*.dig") THEN ! WRITE (*,"(/' The following appear to be fault-trace (f*.dig) files:')") ! files = '*.DIG' ! ! ELSE IF (file_type == "f*.nki") THEN ! WRITE (*,"(/' The following appear to be Fault offset-rate input (f*.nki) files:')") ! files = '*.NKI' ! ELSE IF (file_type == "f*.nko") THEN ! WRITE (*,"(/' The following appear to be Fault offset-rate output (f*.nko) files:')") ! files = '*.NKO' ! ELSE IF (file_type == "p*.nki") THEN ! WRITE (*,"(/' The following appear to be Parameter input (p*.nki) files:')") ! files = '*.NKI' ! ELSE IF (file_type == "s*.nki") THEN ! WRITE (*,"(/' The following appear to be Stress-direction input (s*.nki) files:')") ! files = '*.NKI' ! ELSE IF (file_type == "h*.nko") THEN ! WRITE (*,"(/' The following appear to be model Heave-rate output (h*.nko) files:')") ! files = '*.NKO' ! ELSE IF (file_type == "s*.nko") THEN ! WRITE (*,"(/' The following appear to be interpolated Stress-direction output (s*.nko) files:')") ! files = '*.NKO' ! ELSE IF (file_type == "g*.nko") THEN ! WRITE (*,"(/' The following appear to be reframed geodetic-velocity output (g*.nko) files:')") ! files = '*.NKO' ! ELSE IF (file_type == "v*.out") THEN ! WRITE (*,"(/' The following appear to be velocity output (v*.out) files:')") ! files = '*.OUT' ! (must also filter below to exclude force "f*.out" and text "t*.out" files) ! ELSE ! WRITE (*, "(' ERROR: Unknown file_type (',A,') requested from FileList.')") TRIM(file_type) ! CALL DTraceback ! END IF ! full_to = 0 ! keeps track of use of line ! handle = FILE$FIRST ! flag constant, defined in DFLIB as -1 ! old_result = -999 ! old_name = 'undefined' ! all_files: DO ! result = GETFILEINFOQQ (TRIM(files), info, handle) ! !check for duplicate return of last file (a bug in GETFILEINFOQQ): ! IF (result >= 1) THEN ! duplicate = (result == old_result) .AND. (info.name(1:result) == TRIM(old_name)) ! old_name = info.name(1:result) ! ELSE ! duplicate = .FALSE. ! old_name = ' ' ! END IF ! old_result = result ! !- - - - - - - - - - - - - - - - - - - ! IF (handle == FILE$ERROR) RETURN ! defined in DFLIB as -3 ! IF ((result == 0).OR.duplicate) THEN ! no (new) matching files found ! IF (full_to > 0) THEN ! WRITE (*,"(' ',A)") TRIM(line) ! GO TO 100 ! ELSE IF (count == 0) THEN ! WRITE (*,"(' No such files in directory.')") ! CALL Pause() ! STOP ! END IF ! END IF ! first_letter = info.name(1:1) ! !If looking for f*.dig, reject .dig files that don't start with 'f' ! IF ((file_type == "f*.dig").AND.(.NOT.((first_letter == 'F').OR.(first_letter == 'f')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for f*.nki, reject "input" files that don't start with 'f' ! IF ((file_type == "f*.nki").AND.(.NOT.((first_letter == 'F').OR.(first_letter == 'f')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for p*.nki, reject "parameter" files that don't start with 'p' ! IF ((file_type == "p*.nki").AND.(.NOT.((first_letter == 'P').OR.(first_letter == 'p')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for s*.nki, reject "input" files that don't start with 's' ! IF ((file_type == "s*.nki").AND.(.NOT.((first_letter == 'S').OR.(first_letter == 's')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for e*.nko, reject "output" files that don't start with 'e' ! IF ((file_type == "e*.nko").AND.(.NOT.((first_letter == 'E').OR.(first_letter == 'e')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for f*.nko, reject "output" files that don't start with 'f' ! IF ((file_type == "f*.nko").AND.(.NOT.((first_letter == 'F').OR.(first_letter == 'f')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for h*.nko, reject "output" files that don't start with 'h' ! IF ((file_type == "h*.nko").AND.(.NOT.((first_letter == 'H').OR.(first_letter == 'h')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for s*.nko, reject "output" files that don't start with 's' ! IF ((file_type == "s*.nko").AND.(.NOT.((first_letter == 'S').OR.(first_letter == 's')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for g*.nko, reject "output" files that don't start with 'g' ! IF ((file_type == "g*.nko").AND.(.NOT.((first_letter == 'G').OR.(first_letter == 'g')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If looking for v*.out, reject "velocity" files that don't start with 'v' ! IF ((file_type == "v*.out").AND.(.NOT.((first_letter == 'V').OR.(first_letter == 'v')))) THEN ! IF (handle == FILE$LAST) THEN ! GO TO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If we've gotten this far, we have a qualified file! ! count = count + 1 ! string0 = TRIM(suggested_file) ! CALL DUpper_Case(string0) ! string1 = info.name(1:result) ! string2 = string1 ! CALL DUpper_Case(string2) ! matched = matched .OR. (string0 == string2) ! IF ((full_to + 2 + result) > 70) THEN ! line would overflow ! WRITE (*,"(' ',A)") TRIM(line) ! full_to = 0 ! line = ' ' ! line = info.name(1:result) ! full_to = result ! ELSE ! line can accept this name ! IF (full_to == 0) THEN ! no leading spaces ! line = info.name(1:result) ! full_to = result ! ELSE ! use 2 leading spaces ! line = TRIM(line) // ' ' // info.name(1:result) ! full_to = full_to + 2 + result ! END IF ! END IF ! IF (handle == FILE$LAST) THEN ! IF (full_to > 0) WRITE (*,"(' ',A)") TRIM(line) ! GO TO 100 ! END IF ! END DO all_files ! 100 IF (count == 1) THEN ! collector point, replacing "RETURN" ! ! so that we can adjust suggested_file(?) ! suggested_file = TRIM(string1) ! ELSE IF (count > 1) THEN ! IF (.NOT.matched) THEN ! suggested_file = ' ' ! END IF ! END IF ! END SUBROUTINE File_List SUBROUTINE Get_filename (unit, filename, error) ! Obtains a filename from the beginning of the next line; ! name may be padded with blanks on both sides, and ! may be followed by comments. Note that ! "unit" should have been opened with PAD = "YES". IMPLICIT NONE INTEGER, INTENT (IN) :: unit ! Fortran device number to READ CHARACTER*132, INTENT(OUT) :: filename ! main result of this SUBR LOGICAL, INTENT(OUT) :: error ! IF (ios /= 0) error = .TRUE. CHARACTER(132) :: buffer INTEGER :: i, ios LOGICAL :: past error = .FALSE. ! unless changed below READ (unit,"(A)", IOSTAT = ios) buffer IF (ios /= 0) THEN error = .TRUE. RETURN END IF buffer = ADJUSTL(buffer) ! left-justify past = .FALSE. ! will be T when past end of filename blank_right: DO i=2,132 IF ((buffer(i:i) == ' ') .OR. & (buffer(i:i) == ',') .OR. & (buffer(i:i) == '=') .OR. & (buffer(i:i) == ':')) past = .TRUE. if (past) buffer(i:i) = ' ' END DO blank_right IF (((buffer(1:1) == 'N') .OR. (buffer(1:1) == 'n')) .AND. & ((buffer(2:2) == 'O') .OR. (buffer(2:2) == 'o')) .AND. & ((buffer(3:3) == 'N') .OR. (buffer(3:3) == 'n')) .AND. & ((buffer(4:4) == 'E') .OR. (buffer(4:4) == 'e')) .AND. & (buffer(5:5) == ' ')) buffer = 'none' filename = buffer(1:80) END SUBROUTINE Get_filename SUBROUTINE Get_Parameters ! Reads input parameter file p*.nki ! and memorizes its contents. ! Values reside in global variables. ! This is just a SUBR to avoid repetition of code. ! ! Revised 2004.12.08 to read EITHER NeoKinema v1.x or NeoKinema v.2 ! input parameter files. IMPLICIT NONE INTEGER :: ios, line LOGICAL :: bad_parameter, any_bad_parameters REAL*8 :: t, t1, t2 !CALL File_List( file_type = "p*.nki", & ! & suggested_file = parameter_file) WRITE (*, *) CALL DPrompt_for_String('Which parameter (p*.nki) file should be used?', parameter_file, parameter_file) !------- Try reading file by assuming that it is in NeoKinema v2.x format: OPEN (UNIT = 1, FILE = TRIM(parameter_file), STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) CALL Could_Not_Find_File(parameter_file) any_bad_parameters = .FALSE. ! until an error occurs (below) CALL Get_Filename(1, token, bad_parameter) any_bad_parameters = any_bad_parameters .OR. bad_parameter ! delaying reaction line = 1 ! L0 READ (1, *, IOSTAT = ios) L0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! A0 READ (1, *, IOSTAT = ios) A0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! number of refinements READ (1, *, IOSTAT = ios) n_refine ; line = line + 1 IF (n_refine < 0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! strain-rate uncertainty for rigid blocks READ (1, *, IOSTAT = ios) mu_ ; line = line + 1 IF (mu_ <= 0.0D0) bad_parameter = .TRUE. IF (mu_ < DSQRT(1.1D0 * TINY(mu_))) bad_parameter = .TRUE. IF (mu_ > 1.D-10) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! small strain-rate increment (xi_) for imposing stress-directions READ (1, *, IOSTAT = ios) xi_ ; line = line + 1 IF (xi_ <= 0.0D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] READ (1, *, IOSTAT = ios) sigma_offnormal_degrees ; line = line + 1 IF (sigma_offnormal_degrees < 0.0D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! radius of planet READ (1, *, IOSTAT = ios) t ; line = line + 1 IF (t <= 0.) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter R = t * m_per_km ! minimum and maximum locking depths of intraplate faults, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_min = t1 * m_per_km locking_depth_m_max = t2 * m_per_km ! minimum and maximum locking depths of subduction zones, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_subduction_min = t1 * m_per_km locking_depth_m_subduction_max = t2 * m_per_km ! do new active faults count as sigma_1h data? READ (1, *, IOSTAT = ios) faults_give_sigma_1h ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! names of additional input files (or, "none ") CALL Get_filename (1, f_dat, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter IF (f_dat(1:5) == 'none ') THEN READ (1,*) ; line = line + 1 ! read and ignore f_dig = 'skipped' ELSE CALL Get_filename (1, f_dig, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter END IF CALL Get_filename (1, s_dat, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! stress interpolation method (INTEGER index): {Added in NeoKinema v.4} line = line + 1 READ (1, *, IOSTAT = ios) stress_interpolation_method bad_parameter = (ios /= 0) any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gps_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gp2_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! is velocity reference frame for geodetic data allowed to float? READ (1, *, IOSTAT =ios) floating_frame ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! conservative_geodetic_adjustment? (using geologic slip rates) READ (1, *, IOSTAT = ios) conservative_geodetic_adjustment ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! weight factor for all geodetic data: geodesy_weight = 1.0D0 ! by definition in NeoKinema v2.x and all later versions. CALL Get_filename (1, x_feg, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, x_bcs, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter READ (1,"(A)", IOSTAT = ios) reference_plate_c2; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) CLOSE (UNIT = 1) ! close parameter_file IF (.NOT.any_bad_parameters) THEN got_parameters = .TRUE. RETURN END IF CALL Bad_Parameters() !=========== Note: We only reach here if there was an error during reading ! of parameters under assumption of NeoKinema v2~3 format. END SUBROUTINE Get_Parameters SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause 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". ! The same happens IF (mt_flashby), without waiting for the user. ! 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 SUBROUTINE DTraceback () ! The sole function of this unit is to cause a traceable error, ! so that the route into the unit that called it is also traced. ! This unit is a good place to put a breakpoint while debugging! ! The intentional error must NOT be detected during compilation, ! but MUST cause a traceable error at run-time. ! If this routine has any error detected during compilation, ! then change its code to cause a different intentional error. IMPLICIT NONE CHARACTER*80 instring INTEGER :: i REAL*8, DIMENSION(3) :: y WRITE (*,"(' -----------------------------------------------------')") WRITE (*,"(' Traceback was called to execute an intentional error:')") WRITE (*,"(' An array subscript will be intentionally out-of-range.')") WRITE (*,"(/' After you read this notice, press [Enter]' & & /' to stop the program (no other option): '\)") READ (*,"(A)") instring DO i = 1, 4 y(i) = 1.0D0 * i END DO STOP ' ' END SUBROUTINE DTraceback SUBROUTINE DUpper_Case (string) !Modifies its sole argument so that a..z --> A..Z !Useful* as a filter applied to filenames, before testing ! for a match (*at least, on Windows systems!) IMPLICIT NONE CHARACTER*(*), INTENT(INOUT) :: string INTEGER :: i, j, length length = LEN_TRIM(string) DO i = 1, length j = IACHAR(string(i:i)) IF ((j >= 97).AND.(j <= 122)) THEN ! a..z string(i:i) = ACHAR(j - 32) ! A..Z = 65..90 END IF END DO END SUBROUTINE DUpper_Case END PROGRAM Rake_And_Sliprate_Per_Fault