PROGRAM GPS_Postprocessor ! Utility program to classify residual geodetic velocity misfits ! in a relatively-successful neotectonic model computed with ! NeoKinema in which an augmented .GP2 matrix created by ! GPS_Covariance was used (see Bird & Carafa [2016, JGR]). ! ! It reloads the same .GPS and .GP2 data files used by NeoKinema. ! Then, it looks at benchmarks still present in the g_{token}.nko ! file output from NeoKinema, and shrinks the .GP2 matrix to ! eliminate any benchmarks that were left out of the solution ! (either because they fell outside the .FEG area, or because the ! user was told to remove them, using Delete_Cracked_Benchmarks). ! ! Then, uses MKL to perform an eigenvector/eigenvalue analysis ! on this shrunken (reduced-rank) .GP2 matrix. ! ! Reads the GPS_Covariance_log.txt file (saved from the original ! run of GPS_Covariance that created the augmented .GP2 file), ! and uses this info to assign each eigenvector to one of 3 groups: ! (1) Reference-frame loosening (usually 3 eigenvectors); ! (2) Local noise sources (volcano swelling, postseismic fault afterslip, ! hydrologic opening of mode-1 cracks, very slow landslides, ...); ! (3) Other velocity misfits that are not understood or wanted. ! ! Divides the geodetic velocity misfits reported in the g_{token}.nko ! file (from NeoKinema) into 3 vector components of misfit, at ! each benchmark, that sum to the total velocity misfit. ! ! Finally, uses methods similar to those of NeoKineMap to produce ! a set of 3 related Adobe Illustrator (.AI) files, showing these ! 3 components of the unfit residual geodetic velocities, ! relative to the relevant uncertainty ellipses. ! ! By Peter Bird, UCLA, 2017.11.07 & 2019.04.29 & 2020.09.29 !The following USE statements are needed to produce maps of the misfit components: USE DAdobe_Illustrator ! file DAdobe_Illustrator.f90, from Peter Bird USE DMap_Projections ! file DMap_Projections.f90, from Peter Bird USE DMap_Tools ! file DMap_Tools.f90, from Peter Bird !The following USE statements lead to the eigenvector/eigenvalue routines we need: USE MKL95_PRECISION ! Intel Math Kernel Library, Linear Algebra Package; available in file LAPACK.f90 USE MKL95_LAPACK ! Intel Math Kernel Library, Linear Algebra Package; available in file LAPACK.f90 ! ***REMEMBER*** in the Microsoft Visual Studio GUI provided with ! Intel Parallel Studio XE 2013, you need to set ! Project / Properties / Fortran / Libraries / ! Use Math Kernel Library {to some choice other than "No"}! USE DFLIB, ARCQQ => ARC ! provided with Digital Visual Fortran, and ALSO ! available from within Intel Parallel Studio XE 2013 !(and its GUI of Microsoft Visual Studio), but ONLY by selecting ! Project/Properties/Fortran/Compatibility/Use PowerStation ! Portability Library. {Caution: May only work in 32-bit realm.} ! I am using only GETFILEINFOQQ, which provides names of files ! matching specifications like "p*.nki". Helps user select input file. ! If no substitute is available on your system when you compile, ! just omit SUBROUTINE File_List (and any CALLs to it). ! Note that I am not using their ARC, because I have my own Arc; so I am ! renaming their ARC to ARCQQ [sic] to avoid conflicts. IMPLICIT NONE !GPBglobals CHARACTER*1 :: jobz, trans, uplo ! two one-byte control parameters in some calls to MKL CHARACTER*1 :: class_index_c1, cN, cE CHARACTER*4 :: c4 CHARACTER*5 :: zone ! time zone CHARACTER*6 :: c6 CHARACTER*8 :: date, number8 CHARACTER*9 :: c9 CHARACTER*10 :: clock_time ! wall clock time CHARACTER*12 :: grid_units CHARACTER*80 :: c80 CHARACTER*132 :: path_in = ' ' ! Not using remote folders; all files sit in one current folder. CHARACTER*132 :: path_out = ' ' ! Not using remote folders; all files sit in one current folder. CHARACTER*132 :: f_dig, feg_file, feg_pathfile, g_nko_format, g_nko_header, g_nko_title, g_token_nko_file, & & GPS_file, GPS_Covariance_log_file, GP2_file, GPS_format, GPS_header, GPS_title, & & line, lines_basemap_file, lines_basemap_pathfile, parameter_file, & & suggested_new_ai_filename, & & temp_path_in, temp_title, token, trimmed_token, traces_file, traces_pathfile, volcano_file, volcano_pathfile, x_feg CHARACTER*132 :: bottom_line = ' ', bottom_line_memo = ' ', f_nko_format = ' ', top_line = ' ', top_line_memo = ' ' CHARACTER*132, DIMENSION(20) :: titles CHARACTER*15, DIMENSION(:), ALLOCATABLE :: perturbation_type_c15 CHARACTER*35, DIMENSION(:), ALLOCATABLE :: benchmark_label CHARACTER*80, DIMENSION(:), ALLOCATABLE :: external_GPS_site, internal_GPS_site CHARACTER*132, DIMENSION(:), ALLOCATABLE :: perturbation_source, perturbation_title INTEGER :: benchmarks, choice, class_index, dig_title_method, & & external_benchmarks, external_GPS_DOF, & & gps_type, & & i, ii, info, internal_benchmarks, internal_GPS_DOF, ios, & & j, jBest, jj, jp, jp1, k, kilometers, kk, l, label_thinner, liwork, ll, lp, lwork, m, ma, mb, minutes, & & n, n_checked, na, nb, negative_eigenvalues, nfl, np1, nrhs, nseg, null_perturbations, numel, numnod, & & perturbations, positive_eigenvalues, title_choice, title_count, trace_choice, zero_eigenvalues INTEGER, DIMENSION(8) :: datetimenumber ! for output from DATE_AND_TIME INTEGER, DIMENSION(:), ALLOCATABLE :: eigen_class123, external_number, internal_number, iwork_space, match, perturbation_class123 INTEGER, DIMENSION(:,:), ALLOCATABLE :: nodes LOGICAL :: add_titles, any_titles, bottom, class1_exists, class2_exists, confirm, & & dig_is_lonlat, do_more_overlays, do_overlay, ellipses, first_pass, got_parameters, in_OK, & & mated, plot_dig_titles, polygons, right, verbose, virgin, visible, xy_defined, xy_mode LOGICAL, DIMENSION(:), ALLOCATABLE :: available, taken REAL*8, PARAMETER :: m_per_km = 1.0D3 REAL*8, PARAMETER :: s_per_year = 365.25D0 * 24.D0 * 60.D0 * 60.D0 REAL*8, PARAMETER :: bottomlegend_gap_points = 14.0D0 REAL*8, PARAMETER :: rightlegend_gap_points = 14.0D0 REAL*8 :: arc2, arc3, az1, az2, az3, aze2, aze3, & & benchmark_points, bottomlegend_used_points, & & c_EE, c_EN, c_NE, c_NN, correlation, covariance_11, covariance_12, covariance_22, & & Delta_Elon_deg, dl2, dl3, ds2, ds3, dV_East, dV_North, & & e1, e2, & & factor, & & lat, latitude, leg, lon, longitude, & & maximum, minimum, & & node_radius_points, node_radius_radians, & & R, r2, r2Best, rad, radians_per_point, & & rel_az2, rel_az3, rightlegend_used_points, & & size, start_azimuth, & & t, t1, tick_points, & & u1x, u1y, u2x, u2y, & & v_East_mps, v_South_mps, v_mma, v_mps, velocity_Ma, velocity_magnitude, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, volcano_Elon, volcano_Nlat, volcano_points, & & x0p, x1p, x1_points, x2_points, xp, xpt, y0p, y1p, y1_points, y2_points, yp, ypt REAL*8, DIMENSION(3) :: measure, omega_uvec, result_uvec, uvec, uvec1, uvec2, uvec3 REAL*8, DIMENSION(:), ALLOCATABLE :: amplitudes, benchmark_correlation, & & benchmark_E_sigma, benchmark_E_velocity, benchmark_hypotenuse, benchmark_N_sigma, benchmark_N_velocity, & & class1_vE_sigma, class1_vN_sigma, class1_correlation, & & class2_vE_sigma, class2_vN_sigma, class2_correlation, & & class3_vE_sigma, class3_vN_sigma, class3_correlation, & & copied_vector, & & E_error_mmpa, eigenvalues_list, external_GPS_Elon, external_GPS_Nlat, & & internal_GPS_Elon, internal_GPS_Nlat, & & N_error_mmpa, work, & & perturbation_sizes_mmpa, & & train, & & unfit_eigenamplitudes, unfit_EN, unfit_EN_class1, unfit_EN_class2, & & unfit_EN_class1and2, unfit_EN_class3, unfit_supervector REAL*8, DIMENSION(:,:), ALLOCATABLE :: benchmark_uvec, copied_matrix, & & external_GP2, internal_GP2, eigenvectors_block, & & node_uvec, & & perturbation_unit_vectors, perturbation_vectors, principal_axes REAL*8, DIMENSION(:,:,:),ALLOCATABLE :: segments !==================================================================================================== ! Conventions observed throughout this code: ! (1) GPS velocities and misfits are in units of mm/a; ! therefore, velocity covariance matrix is in units of (mm/a)**2. ! (2) Order of benchmarks in the covariance matrices is the same as in the corresponding .GPS file ! (for the original benchmark list, input to NeoKinema) or g_{token}.nko file ! (output from NeoKinema, and possibly including fewer benchmarks). ! (3) Ordering of velocity components in all covariance matrices is dV_East before dV_North. !==================================================================================================== !Begin recording all output in a log-file: OPEN (UNIT = 71, FILE = "GPS_Postprocessor_log.txt") ! absolute OPEN; overwrites any existing file with same name CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) WRITE (71,"('GPS_Postprocessor_log.txt, produced on ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & datetimenumber(1), datetimenumber(2), datetimenumber(3), & datetimenumber(5), datetimenumber(6), datetimenumber(7) WRITE (*, "(' ')") WRITE (*, "(' PROGRAM GPS_Postprocessor:')") WRITE (*, "(' ')") WRITE (*, "(' Utility program to classify residual geodetic velocity misfits')") WRITE (*, "(' in a relatively-successful neotectonic model computed with')") WRITE (*, "(' NeoKinema in which an augmented .GP2 matrix created by')") WRITE (*, "(' GPS_Covariance was used (see Bird & Carafa [2016, JGR]).')") WRITE (*, "(' ')") WRITE (*, "(' It reloads the same .GPS and .GP2 data files used by NeoKinema.')") WRITE (*, "(' Then, it looks at benchmarks still present in the g_{token}.nko')") WRITE (*, "(' file output from NeoKinema, and shrinks the .GP2 matrix to')") WRITE (*, "(' eliminate any benchmarks that were left out of the solution')") WRITE (*, "(' (either because they fell outside the .FEG area, or because the')") WRITE (*, "(' user was told to remove them, using Delete_Cracked_Benchmarks).')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' Then, uses MKL to perform an eigenvector/eigenvalue analysis')") WRITE (*, "(' on this shrunken (reduced-rank) .GP2 matrix.')") WRITE (*, "(' ')") WRITE (*, "(' Reads the GPS_Covariance_log.txt file (saved from the original')") WRITE (*, "(' run of GPS_Covariance that created the augmented .GP2 file),')") WRITE (*, "(' and uses this info to assign each eigenvector to one of 3 groups:')") WRITE (*, "(' (1) Reference-frame loosening (usually 3 eigenvectors);')") WRITE (*, "(' (2) Local noise sources (volcano swelling, postseismic fault afterslip,')") WRITE (*, "(' hydrologic opening of mode-1 cracks, very slow landslides, ...);')") WRITE (*, "(' (3) Other velocity misfits that are not understood or wanted.')") WRITE (*, "(' ')") WRITE (*, "(' Divides the geodetic velocity misfits reported in the g_{token}.nko')") WRITE (*, "(' file (from NeoKinema) into 3 vector components of misfit, at')") WRITE (*, "(' each benchmark, that sum to the total velocity misfit.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' Finally, uses methods similar to those of NeoKineMap to produce')") WRITE (*, "(' a set of 3 related Adobe Illustrator (.AI) files, showing these')") WRITE (*, "(' 3 components of the unfit residual geodetic velocities,')") WRITE (*, "(' relative to the relevant uncertainty ellipses.')") WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird, UCLA, 2017.11.07 & 2019.04.29 & 2020.09.29')") WRITE (*, *) CALL Pause() WRITE (71, *) WRITE (71, "('PROGRAM GPS_Postprocessor:')") WRITE (71, *) WRITE (71, "('Utility program to classify residual geodetic velocity misfits')") WRITE (71, "('in a relatively-successful neotectonic model computed with')") WRITE (71, "('NeoKinema in which an augmented .GP2 matrix created by')") WRITE (71, "('GPS_Covariance was used (see Bird & Carafa [2016, JGR]).')") WRITE (71, *) WRITE (71, "('It reloads the same .GPS and .GP2 data files used by NeoKinema.')") WRITE (71, "('Then, it looks at benchmarks still present in the g_{token}.nko')") WRITE (71, "('file output from NeoKinema, and shrinks the .GP2 matrix to')") WRITE (71, "('eliminate any benchmarks that were left out of the solution')") WRITE (71, "('(either because they fell outside the .FEG area, or because the')") WRITE (71, "('user was told to remove them, using Delete_Cracked_Benchmarks).')") WRITE (71, *) WRITE (71, "('Then, uses MKL to perform an eigenvector/eigenvalue analysis')") WRITE (71, "('on this shrunken (reduced-rank) .GP2 matrix.')") WRITE (71, *) WRITE (71, "('Reads the GPS_Covariance_log.txt file (saved from the original')") WRITE (71, "('run of GPS_Covariance that created the augmented .GP2 file),')") WRITE (71, "('and uses this info to assign each eigenvector to one of 3 groups:')") WRITE (71, "('(1) Reference-frame loosening (usually 3 eigenvectors);')") WRITE (71, "('(2) Local noise sources (volcano swelling, postseismic fault afterslip,')") WRITE (71, "(' hydrologic opening of mode-1 cracks, very slow landslides, ...);')") WRITE (71, "('(3) Other velocity misfits that are not understood or wanted.')") WRITE (71, *) WRITE (71, "('Divides the geodetic velocity misfits reported in the g_{token}.nko')") WRITE (71, "('file (from NeoKinema) into 3 vector components of misfit, at')") WRITE (71, "('each benchmark, that sum to the total velocity misfit.')") WRITE (71, *) WRITE (71, "('Finally, uses methods similar to those of NeoKineMap to produce')") WRITE (71, "('a set of 3 related Adobe Illustrator (.AI) files, showing these')") WRITE (71, "('3 components of the unfit residual geodetic velocities,')") WRITE (71, "('relative to the relevant uncertainty ellipses.')") WRITE (71, *) WRITE (71, "('By Peter Bird, UCLA, 2017.11.07 & 2019.04.29 & 2020.09.29')") WRITE (71, *) ! open the Parameter (p*.nki) file and get the name token from the first line: WRITE (*, *) WRITE (*, "(' Enter the name of a parameter file (p*.nki) used in the NeoKinema run:')") parameter_file = ' ' !CALL File_List( file_type = "p*.nki", & ! & suggested_file = parameter_file, & ! & using_path = path_in ) CALL DPrompt_for_String('Which parameter file was used?', parameter_file, parameter_file) WRITE (71, *) WRITE (71, "('Enter the name of a parameter file (p*.nki) used in the NeoKinema run:')") WRITE (71, "(A)") TRIM(parameter_file) OPEN (UNIT = 1, FILE = parameter_file, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This parameter file was not found (in current folder).')") CALL Pause() STOP END IF token = Get_Filename(unit = 1) ! will generate CALL_Bad_Parameters() if an error occurs IF (token(1:1) == '_') THEN trimmed_token = token(2:132) ELSE trimmed_token = token END IF WRITE (*, "(' Name-token for this NeoKinema run was: ', A)") TRIM(token) !Skip over the next 6 lines of the p_.nki file: DO i = 1, 6 READ (1, *) END DO !Read Earth radius, in km: READ (1, *) t R = t * m_per_km !Skip over the next 4 lines of the p_.nki file: DO i = 1, 4 READ (1, *) END DO f_dig = Get_Filename(unit = 1) ! name of the fault-traces.DIG file used by NeoKinema (for possible plotting) READ (1, *) ! s_.nki READ (1, *) ! stress interpolation method !Read names of .GPS file and .GP2 file: GPS_file = Get_Filename(unit = 1) GP2_file = Get_Filename(unit = 1) IF ((GP2_file == "none").OR.(GP2_file == "NONE").OR.(GP2_file == "None")) THEN WRITE (*, "(' ERROR: It is not appropriate to use GPS_Postprocessor to analyze this NeoKinema model,')") WRITE (*, "(' because no velocity covariance (.GP2) file was used in its creation.')") WRITE (71, "('ERROR: It is not appropriate to use GPS_Postprocessor to analyze this NeoKinema model,')") WRITE (71, "(' because no velocity covariance (.GP2) file was used in its creation.')") CALL Pause() STOP END IF READ (1, *) ! floating_frame? READ (1, *) ! conservative_geodetic_adjustment? x_feg = Get_filename(unit = 1) ! name of .FEG file used by NeoKinema CLOSE (1) ! p_{token}.nki file, no longer needed !Warn user that these .GPS and .GP2 files will be needed. WRITE (*, *) WRITE (*, "(' NOTE that the following two files, used by NeoKinema, need to be available:')") WRITE (*, "(' ', A)") TRIM(GPS_file) WRITE (*, "(' ', A)") TRIM(GP2_file) WRITE (*, "(' (If they are not in the current folder, please move them in now.)')") CALL Pause() WRITE (71, *) WRITE (71, "('NOTE that the following two files, used by NeoKinema, need to be available:')") WRITE (71, "(' ', A)") TRIM(GPS_file) WRITE (71, "(' ', A)") TRIM(GP2_file) WRITE (71, "('(If they are not in the current folder, please move them in now.)')") !Open the .GPS file and count the benchmarks: OPEN(UNIT = 10, FILE = TRIM(GPS_file), STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This .GPS file was not found (in current folder).')") WRITE (71, "('ERROR: This .GPS file was not found (in current folder).')") CALL Pause() STOP END IF 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 CLOSE (10) ALLOCATE ( external_GPS_Elon(external_benchmarks) ) ALLOCATE ( external_GPS_Nlat(external_benchmarks) ) ALLOCATE ( external_GPS_site(external_benchmarks) ) !OPEN the .GPS file a second time, and memorize benchmark locations and names: 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, site names and (Elon, Nlat) are remembered: DO i = 1, external_benchmarks READ (10, GPS_format, IOSTAT = ios) external_GPS_Elon(i), external_GPS_Nlat(i), & & vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation, c80, external_GPS_site(i) IF (ios /= 0) THEN WRITE (*, "(' ERROR while reading (Elon, Nlat) and SITE name of benchmark #',I6)") i WRITE (71, "('ERROR while reading (Elon, Nlat) and SITE name of benchmark #',I6)") i CALL Pause() STOP END IF END DO CLOSE (10) WRITE (*, *) WRITE (*, "(' ',I8,' GPS sites were read from ', A)") external_benchmarks, TRIM(GPS_file) WRITE (71, "(I8,' GPS sites were read from ', A)") external_benchmarks, TRIM(GPS_file) !Read the original (external) .GP2 file: external_GPS_DOF = 2 * external_benchmarks ALLOCATE ( external_GP2(external_GPS_DOF, external_GPS_DOF) ) external_GP2 = 0.0D0 ! whole matrix, because zero entries may be omitted from the input file. OPEN (UNIT = 11, FILE = TRIM(GP2_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This expected .GP2 file was not found (in current folder).')") WRITE (71, "('ERROR: This expected .GP2 file was not found (in current folder).')") CALL Pause() STOP END IF entering: DO READ (11, *, IOSTAT = ios) i, j, t IF (ios /= 0) EXIT entering IF ((i == j).AND.(t <= 0.0D0)) THEN ! prevent non-positive entries on diagonal WRITE (*, "(' ERROR: Diagonal element is non-positive: ', 2I6, E12.4)") i, j, t WRITE (71, "('ERROR: Diagonal element is non-positive: ', 2I6, E12.4)") i, j, t CALL Pause() STOP END IF external_GP2(i, j) = t external_GP2(j, i) = t ! Because only upper-triangle + diagonal (or lower-triangle + diagonal) is required in the .GP2 file. END DO entering CLOSE (11) WRITE (*, "(' .GP2 file ', A, ' has also been read.')") TRIM(GP2_file) WRITE (71, "('.GP2 file ', A, ' has also been read.')") TRIM(GP2_file) !Build g_token_nko_file name and warn user that it is needed: g_token_nko_file = 'g' // TRIM(token) // ".nko" WRITE (*, *) WRITE (*, "(' NOTE that the following output file, created by NeoKinema, must be available:')") WRITE (*, "(' ', A)") TRIM(g_token_nko_file) WRITE (*, "(' (If it is not in the current folder, please move it in now.)')") CALL Pause() WRITE (71, *) WRITE (71, "('NOTE that the following output file, created by NeoKinema, must be available:')") WRITE (71, "(' ', A)") TRIM(g_token_nko_file) WRITE (71, "('(If it is not in the current folder, please move it in now.)')") !Read the g_token_nko file once in order to count benchmarks: OPEN (UNIT = 13, FILE = TRIM(g_token_nko_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This expected file was not found (in current folder).')") WRITE (71, "('ERROR: This expected file was not found (in current folder).')") CALL Pause() STOP END IF READ (13, "(A)") g_nko_title READ (13, "(A)") g_nko_format READ (13, "(A)") g_nko_header internal_benchmarks = 0 ! just initializing interating: DO READ (13, g_nko_format, IOSTAT = ios) lon, lat IF (ios /= 0) EXIT interating internal_benchmarks = internal_benchmarks + 1 END DO interating CLOSE (13) WRITE (*, *) WRITE (*, "(' ', I6, ' benchmarks found in file ', A)") internal_benchmarks, TRIM(g_token_nko_file) WRITE (71, "(I6, ' benchmarks found in file ', A)") internal_benchmarks, TRIM(g_token_nko_file) !Memorize essential information from the g_token_nko file: ALLOCATE ( internal_GPS_Elon(internal_benchmarks) ) ALLOCATE ( internal_GPS_Nlat(internal_benchmarks) ) ALLOCATE ( internal_GPS_site(internal_benchmarks) ) ALLOCATE ( E_error_mmpa(internal_benchmarks) ) ALLOCATE ( N_error_mmpa(internal_benchmarks) ) ALLOCATE ( class3_vE_sigma(internal_benchmarks) ) ! to be remembered for plotting maps of Class-3 unfit residuals ALLOCATE ( class3_vN_sigma(internal_benchmarks) ) ! ditto ALLOCATE ( class3_correlation(internal_benchmarks) ) ! ditto OPEN (UNIT = 13, FILE = TRIM(g_token_nko_file), STATUS = "OLD") READ (13, "(A)") g_nko_title READ (13, "(A)") g_nko_format READ (13, "(A)") g_nko_header DO i = 1, internal_benchmarks READ (13, g_nko_format, IOSTAT = ios) internal_GPS_Elon(i), internal_GPS_Nlat(i), & & vE_mmpa, vN_mmpa, class3_vE_sigma(i), class3_vN_sigma(i), class3_correlation(i), & & c9, internal_GPS_site(i), & & E_error_mmpa(i), N_error_mmpa(i) IF (ios /= 0) THEN WRITE (*, "(' ERROR reading (Elon, Nlat), sitename, (err_E, err_N) of benchmark #', I6)") i WRITE (71, "('ERROR reading (Elon, Nlat), sitename, (err_E, err_N) of benchmark #', I6)") i CALL Pause() STOP END IF internal_GPS_site(i) = ADJUSTL(internal_GPS_site(i)) END DO CLOSE (13) !Establish cross-reference between internal (fewer?) and external (more?) numbering of benchmarks. ALLOCATE ( internal_number(external_benchmarks) ) ! longer list, but some entries will remain 0 internal_number = 0 ! initializing whole list ALLOCATE ( external_number(internal_benchmarks) ) ! shorter list, and no 0's should remain! external_number = 0 ! initializing whole list IF (internal_benchmarks == external_benchmarks) THEN DO i = 1, internal_benchmarks internal_number(i) = i external_number(i) = i END DO ELSE ! more common case, where internal_number < external_number ALLOCATE ( available(external_benchmarks) ) available = .TRUE. ! for all external sites, initially DO i = 1, internal_benchmarks r2Best = 999.9D9 ! just initializing as a "huge number" of degrees**2; will be replaced jBest = 0 ! just initializing; will be replaced DO j = 1, external_benchmarks IF (available(j)) THEN Delta_Elon_deg = (external_GPS_Elon(j) - internal_GPS_Elon(i)) ! BUT, this may contain a cycle-shift, or even two! IF (Delta_Elon_deg > 180.0D0) Delta_Elon_deg = Delta_Elon_deg - 360.0D0 IF (Delta_Elon_deg > 180.0D0) Delta_Elon_deg = Delta_Elon_deg - 360.0D0 IF (Delta_Elon_deg < -180.0D0) Delta_Elon_deg = Delta_Elon_deg + 360.0D0 IF (Delta_Elon_deg < -180.0D0) Delta_Elon_deg = Delta_Elon_deg + 360.0D0 r2 = Delta_Elon_deg**2 + (external_GPS_Nlat(j) - internal_GPS_Nlat(i))**2 IF (r2 < r2Best) THEN r2Best = r2 jBest = j available(j) = .FALSE. ! because it is already mated END IF END IF END DO ! survey over all external benchmarks external_number(i) = jBest internal_number(jBest) = i END DO ! for all internal benchmarks END IF ! internal_benchmarks == external_benchmarks, or NOT !Create internal_GP2 matrix: internal_GPS_DOF = 2 * internal_benchmarks ALLOCATE ( internal_GP2(internal_GPS_DOF, internal_GPS_DOF) ) internal_GP2 = 0.0D0 ! just initializing, to simplify debugging IF (internal_benchmarks == external_benchmarks) THEN ! There is no difference between these two matrices. internal_GP2 = external_GP2 ELSE ! more common case; internal_GP2 describes fewer benchmarks !Pick and choose which numbers will get preserved: DO i = 1, internal_benchmarks ! row number, in terms of benchmarks DO j = 1, internal_benchmarks ! column number, in terms of benchmarks !(dV_E, dV_E): ii = 2 * i - 1 jj = 2 * j - 1 kk = 2 * external_number(i) - 1 ll = 2 * external_number(j) - 1 internal_GP2(ii, jj) = external_GP2(kk, ll) internal_GP2(jj, ii) = external_GP2(kk, ll) !(dV_E, dV_N): ii = 2 * i - 1 jj = 2 * j kk = 2 * external_number(i) - 1 ll = 2 * external_number(j) internal_GP2(ii, jj) = external_GP2(kk, ll) internal_GP2(jj, ii) = external_GP2(kk, ll) !(dV_N, dV_E): ii = 2 * i jj = 2 * j - 1 kk = 2 * external_number(i) ll = 2 * external_number(j) - 1 internal_GP2(ii, jj) = external_GP2(kk, ll) internal_GP2(jj, ii) = external_GP2(kk, ll) !(dV_N, dV_N): ii = 2 * i jj = 2 * j kk = 2 * external_number(i) ll = 2 * external_number(j) internal_GP2(ii, jj) = external_GP2(kk, ll) internal_GP2(jj, ii) = external_GP2(kk, ll) END DO END DO END IF WRITE (*, *) WRITE (*, "(' Attempting eigenvector/eigenvalue analysis of the covariance matrix.')") WRITE (*, *) ! (to guarantee that any error messages will be readable) WRITE (71, *) WRITE (71, "('Attempting eigenvector/eigenvalue analysis of the covariance matrix.')") jobz = 'V' ! required when you want eigenvectors as well as eigenvalues uplo = 'U' ! although 'L' would also work lwork = 2 * internal_GPS_DOF**2 + 6 * internal_GPS_DOF + 1 ! see documentation for MKL routine "dsyevd" for value ALLOCATE ( work(lwork) ) ! "work" is REAL*8, DIMENSION(:) work = 1.0D-30 ! prevent any NaN, LOG(<=0.), or SQRT(<0.) AbEnds due to speculative execution (or other bug) in MKL. liwork = 5 * internal_GPS_DOF + 3 ! see documentation for MKL routine "dsyevd" for value ALLOCATE ( iwork_space(liwork) ) ! INTEGER, DIMENSION(:); local in this SUBR iwork_space = 1 ! prevent any subscript-out-of-range or memory-protextion AbEnds due to speculative execution (or other bug) in MKL. ALLOCATE ( eigenvalues_list(internal_GPS_DOF) ) eigenvalues_list = 1.0D-30 ! prevent any NaN, LOG(<=0.), or SQRT(<0.) AbEnds due to speculative execution (or other bug) in MKL. ALLOCATE ( eigenvectors_block(internal_GPS_DOF, internal_GPS_DOF) ) eigenvectors_block = 1.0D-30 ! prevent any NaN, LOG(<=0.), or SQRT(<0.) AbEnds due to speculative execution (or other bug) in MKL. eigenvectors_block = internal_GP2 ! copy input matrix into eventual output array CALL dsyevd(jobz, uplo, internal_GPS_DOF, eigenvectors_block, internal_GPS_DOF, eigenvalues_list, work, lwork, iwork_space, liwork, info) ! using Fortran77 call because F95 call is buggy. !WARNING: The CALL above sometimes leads to mysterious illegal-value AbEnds, unless all work arrays are predefined to safe values! IF (info /= 0) THEN WRITE (*, "(' ERROR: dsyevd of MKL could not solve eigenvalues/vectors of internal_GP2'/' info = ',I12)") info WRITE (71, "('ERROR: dsyevd of MKL could not solve eigenvalues/vectors of internal_GP2'/' info = ',I12)") info CALL Pause() STOP END IF DEALLOCATE ( work ) ! so that a new, different "work" vector can be created later for another CALL to MKL/LAPACK. !Check that all eigenvalues are positive(?) positive_eigenvalues = 0 ! just initializing, before sum negative_eigenvalues = 0 ! just initializing, before sum zero_eigenvalues = 0 ! just initializing, before sum DO i = 1, internal_GPS_DOF IF (eigenvalues_list(i) > 0.0D0) THEN positive_eigenvalues = positive_eigenvalues + 1 ELSE IF (eigenvalues_list(i) < 0.0D0) THEN negative_eigenvalues = negative_eigenvalues + 1 ELSE zero_eigenvalues = zero_eigenvalues + 1 END IF END DO IF (positive_eigenvalues /= internal_GPS_DOF) THEN WRITE (*, *) WRITE (*, "(' ERROR: Eigenanalysis of internal_GP2 returned non-positive eigenvalue(s):')") WRITE (*, "(' negative: ', I6, ', zero: ', I6,', positive: ', I6)") negative_eigenvalues, zero_eigenvalues, positive_eigenvalues WRITE (71, *) WRITE (71, "('ERROR: Eigenanalysis of internal_GP2 returned non-positive eigenvalue(s):')") WRITE (71, "('negative: ', I6, ', zero: ', I6,', positive: ', I6)") negative_eigenvalues, zero_eigenvalues, positive_eigenvalues CALL Pause() DO i = 1, internal_GPS_DOF WRITE (*, "(' ', I10, ES15.6)") i, eigenvalues_list(i) WRITE (71, "(I10, ES15.6)") i, eigenvalues_list(i) END DO WRITE (*, "(' First (most-negative) eigenvalue has this eigenvector:')") WRITE (71, "('First (most-negative) eigenvalue has this eigenvector:')") DO i = 1, internal_GPS_DOF WRITE (*, "(' ', I10, ES15.6)") i, eigenvectors_block(i, 1) WRITE (71, "(I10, ES15.6)") i, eigenvectors_block(i, 1) END DO STOP END IF WRITE (*, *) WRITE (*, "(' Eigenanalysis finished and succeeded.')") WRITE (71, *) WRITE (71, "('Eigenanalysis finished and succeeded.')") !Describe unfit velocity components as one supervector: ALLOCATE ( unfit_supervector(internal_GPS_DOF) ) unfit_supervector = 0.0D0 ! but all will be replaced immediately... DO i = 1, internal_benchmarks unfit_supervector(2 * i - 1) = -E_error_mmpa(i) unfit_supervector(2 * i) = -N_error_mmpa(i) !Note that the negative sign converts "model error" velocities to "unfit" velocities (i.e., problems in the data). END DO !Describe unfit velocity components as unfit_eigenamplitudes: ALLOCATE ( unfit_eigenamplitudes(internal_GPS_DOF) ) unfit_eigenamplitudes = 0.0D0 ! initializing, before sums below... !Notes: Each column of eigenvectors_block (eigenvectors_block(j, i), j = 1, internal_GPS_DOF) is one eigenvector, #i. ! The eigenvalues are arranged in ascending order in eigenvalues_list. DO i = 1, internal_GPS_DOF DO j = 1, internal_GPS_DOF unfit_eigenamplitudes(i) = unfit_eigenamplitudes(i) + unfit_supervector(j) * eigenvectors_block(j, i) END DO END DO !Describe results of eigenanalysis of internal_GP2 with a matrix of principal_axes (analogous to velocity perturbations): ALLOCATE ( principal_axes(internal_GPS_DOF, internal_GPS_DOF) ) principal_axes = 0.0D0 ! but these will be replaced, immediately... DO j = 1, internal_GPS_DOF ! working on one column at a time... factor = SQRT(eigenvalues_list(j)) DO i = 1, internal_GPS_DOF principal_axes(i, j) = factor * eigenvectors_block(i, j) END DO END DO !Open the GPS_Covariance_log.txt WRITE (*, *) WRITE (*, "(' Attempting to read GPS_Covariance_log.txt (or re-named equivalent)...')") WRITE (71, *) WRITE (71, "('Attempting to read GPS_Covariance_log.txt (or re-named equivalent)...')") CALL DPrompt_for_String("What is the name of your GPS_Covariance log-file?", "GPS_Covariance_log.txt", GPS_Covariance_log_file) OPEN (UNIT = 15, FILE = TRIM(GPS_Covariance_log_file), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ', A, ' not found (in current folder).')") TRIM(GPS_Covariance_log_file) WRITE (71, "('ERROR: ', A, ' not found (in current folder).')") TRIM(GPS_Covariance_log_file) CALL Pause() STOP END IF !Display header and footer from this file, and ask user to confirm... WRITE (*, *) WRITE (*, "(' Below are the headers and footers from this file:')") WRITE (*, "(' ===============================================================================')") WRITE (71, *) WRITE (71, "('Below are the headers and footers from this file:')") WRITE (71, "('===============================================================================')") verbose = .TRUE. perturbations = 0 ! just initializing this count class1_exists = .FALSE. ! unless changed below class2_exists = .FALSE. ! unless changed below browsing: DO READ (15, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT browsing IF (verbose) THEN WRITE (*, "(' ', A)") TRIM(line) WRITE (71, "(A)") TRIM(line) END IF IF (line(1:5) == "=====") verbose = .NOT.verbose IF ((line(1:15) == "REFERENCE-FRAME").OR.(line(1:7) == "VOLCANO").OR.(line(1:12) == "POST-SEISMIC").OR.(line(1:13) == "OPENING CRACK").OR.(line(1:5) == "EXTRA")) THEN perturbations = perturbations + 1 IF (line(1:15) == "REFERENCE-FRAME") class1_exists = .TRUE. IF ((line(1:7) == "VOLCANO").OR.(line(1:12) == "POST-SEISMIC").OR.(line(1:13) == "OPENING CRACK").OR.(line(1:5) == "EXTRA")) class2_exists = .TRUE. END IF ! a CAPITALIZED TITLE was read END DO browsing WRITE (*, "(' ===============================================================================')") WRITE (*, "(' ', I6, ' velocity-perturbations vectors are listed in this file.')") perturbations WRITE (71, "('===============================================================================')") WRITE (71, "(I6, ' velocity-perturbations vectors are listed in this file.')") perturbations CLOSE (15) ! GPS_Covariance_log.txt IF (perturbations == 0) THEN WRITE (*, "(' ERROR: NO velocity-perturbation vectors were found.')") WRITE (*, "(' If this IS the right log-file, then you cannot use GPS_Postprocessor.')") WRITE (71, "('ERROR: NO velocity-perturbation vectors were found.')") WRITE (71, "('If this IS the right log-file, then you cannot use GPS_Postprocessor.')") CALL Pause() END IF confirm = .TRUE. ! for default suggestion CALL DPrompt_for_Logical('Can you confirm that this is the RIGHT log file?', confirm, confirm) WRITE (71, "('Can you confirm that this is the RIGHT log file?')") WRITE (71, *) confirm IF (.NOT.confirm) THEN WRITE (*, "(' Please put the CORRECT GPS_Covariance_log.txt in the current folder, & restart.')") WRITE (71, "('Please put the CORRECT GPS_Covariance_log.txt in the current folder, & restart.')") CALL Pause() STOP END IF !Re-read GPS_Covariance_log.txt, and memorize essential information: ALLOCATE ( perturbation_class123(perturbations) ) perturbation_class123 = 0 ! just initializing; each to become 1 or 2 later on... ALLOCATE ( perturbation_vectors(internal_GPS_DOF, perturbations) ) perturbation_vectors = 0.0D0 ! but some will be replaced based on file contents. ALLOCATE ( perturbation_type_c15(perturbations) ) ! CHARACTER*15 perturbation_type_c15 = ' ' ALLOCATE ( perturbation_title(perturbations) ) ! CHARACTER*132 perturbation_title = ' ' ALLOCATE ( perturbation_source(perturbations) ) ! CHARACTER*132 perturbation_source = ' ' OPEN (UNIT = 15, FILE = TRIM(GPS_Covariance_log_file), STATUS = "OLD", PAD = "YES") DO i = 1, 4 ! get past header READ (15, *) END DO DO i = 1, perturbations READ (15, "(A)") perturbation_type_c15(i) IF (perturbation_type_c15(i) == "REFERENCE-FRAME") THEN perturbation_class123(i) = 1 ELSE ! loose group including VOLCANO, POST-SEISMIC, OPENING CRACK, & EXTRA perturbation_class123(i) = 2 END IF READ (15, "(A)") perturbation_source(i) READ (15, "(A)") perturbation_title(i) READ (15, *) ! column-headers line unpacking: DO ! indefinite loop, because number of moving benchmarks is unknown READ (15, "(A)") line IF (line(1:5) == "-----") EXIT unpacking ! ; ELSE... READ (line, *) lon, lat, dV_East, dV_North ! in mm/a !Find this benchmark in the larger external list: jBest = 0 r2Best = 999.9D9 DO j = 1, external_benchmarks Delta_Elon_deg = (lon - external_GPS_Elon(j)) IF (Delta_Elon_deg > 180.0D0) Delta_Elon_deg = Delta_Elon_deg - 360.0D0 IF (Delta_Elon_deg > 180.0D0) Delta_Elon_deg = Delta_Elon_deg - 360.0D0 IF (Delta_Elon_deg < -180.0D0) Delta_Elon_deg = Delta_Elon_deg + 360.0D0 IF (Delta_Elon_deg < -180.0D0) Delta_Elon_deg = Delta_Elon_deg + 360.0D0 r2 = Delta_Elon_deg**2 + (lat - external_GPS_Nlat(j))**2 IF (r2 < r2Best) THEN r2Best = r2 jBest = j END IF END DO ! jBest now points to this benchmark in the external list IF (internal_number(jBest) > 0) THEN ! this benchmark is also on the internal list jBest = internal_number(jBest) ! switch pointer to internal benchmark list perturbation_vectors(2 * jBest - 1, i) = dV_East perturbation_vectors(2 * jBest , i) = dV_North END IF ! this benchmark is also on the internal list END DO unpacking END DO ! i = 1, perturbations CLOSE (15) !Convert perturbation_vectors to perturbation_unit_vectors: ALLOCATE ( perturbation_unit_vectors(internal_GPS_DOF, perturbations) ) ALLOCATE ( perturbation_sizes_mmpa(perturbations) ) null_perturbations = 0 ! (and we hope it stays there...) perturbation_unit_vectors = 0.0D0 DO i = 1, perturbations r2 = 0.0D0 DO j = 1, internal_GPS_DOF r2 = r2 + perturbation_vectors(j, i)**2 END DO size = SQRT(r2) IF (size > 0.0D0) THEN ! usual case. BUT, can be 0.0D0 when one "noise" source is incredibly far from ALL benchmarks... (e.g., Yellowstone, re CA GPS). perturbation_unit_vectors(1:internal_GPS_DOF, i) = perturbation_vectors(1:internal_GPS_DOF, i) / size ELSE ! size == 0.0D0; bad perturbation! null_perturbations = null_perturbations + 1 perturbation_unit_vectors(1:1 , i) = 1.0D0 ! This is a BOGUS unit-vector, but we have no intention of actually USING it. perturbation_unit_vectors(2:internal_GPS_DOF, i) = 0.0D0 END IF perturbation_sizes_mmpa(i) = size ! save (including 0's) for later use, when assessing significance of amplitudes(1:perturbations) in sigma-type units. END DO IF (null_perturbations == 1) THEN WRITE (*, *) WRITE (71, *) WRITE (*, "(' Actually, ONE perturbation was null; after condensing this out,'/ ' perturbations = ', I6)") (perturbations - 1) WRITE (71, "('Actually, ONE perturbation was null; after condensing this out, perturbations = ', I6)") (perturbations - 1) CALL Pause() oneShift: DO i = 1, (perturbations-1) ! If null, then shift-down all higher-numbered entries: IF (perturbation_sizes_mmpa(i) == 0.0D0) THEN DO j = i, (perturbations-1) ! entries that need to be redefined, as equal to higher-#ed values: perturbation_sizes_mmpa(j) = perturbation_sizes_mmpa(j+1) perturbation_class123(j) = perturbation_class123(j+1) perturbation_vectors(1:internal_GPS_DOF, j) = perturbation_vectors(1:internal_GPS_DOF, j+1) perturbation_type_c15(j) = perturbation_type_c15(j+1) perturbation_title(j) = perturbation_title(j+1) perturbation_source(j) = perturbation_source(j+1) END DO EXIT oneShift END IF END DO oneShift ELSE IF (null_perturbations >= 2) THEN WRITE (*, *) WRITE (71, *) WRITE (*, "(' Actually, ', I6, ' perturbations were null; after condensing them out,'/ ' perturbations = ', I6)") (perturbations - null_perturbations) WRITE (71, "('Actually, ', I6, ' perturbations were null; after condensing them out, perturbations = ', I6)") (perturbations - null_perturbations) CALL Pause() !For multiple compactions, must use a loop with (possibly) non-incrementing subscript! n_checked = 0 i = 1 ! initializing (possibly non-incrementing) loop variable multiShift: DO IF (perturbation_sizes_mmpa(i) > 0.0D0) THEN ! good, lets move along to next... n_checked = n_checked + 1 IF (n_checked >= (perturbations - null_perturbations)) EXIT multiShift IF (i >= perturbations) EXIT multiShift i = i + 1 ! else, increment i++ and go again ELSE ! shift everything down by one (without checking whether replacement is ALSO null). DO j = i, (perturbations-1) ! entries that need to be redefined, as equal to higher-#ed values: perturbation_sizes_mmpa(j) = perturbation_sizes_mmpa(j+1) perturbation_class123(j) = perturbation_class123(j+1) perturbation_vectors(1:internal_GPS_DOF, j) = perturbation_vectors(1:internal_GPS_DOF, j+1) perturbation_type_c15(j) = perturbation_type_c15(j+1) perturbation_title(j) = perturbation_title(j+1) perturbation_source(j) = perturbation_source(j+1) END DO !Note that i is NOT incremented in this branch! So, it will be checked, in the next pass! END IF END DO multiShift END IF perturbations = perturbations - null_perturbations IF (perturbations == 0) THEN WRITE (*, "(' ERROR: NO velocity-perturbation vectors reamin.')") WRITE (*, "(' If this IS the right log-file, then you cannot use GPS_Postprocessor.')") WRITE (71, "('ERROR: NO velocity-perturbation vectors remain.')") WRITE (71, "('If this IS the right log-file, then you cannot use GPS_Postprocessor.')") CALL Pause() STOP END IF !Match each perturbation vector to the closest principal-axis (or negative of principal-axis): ALLOCATE ( taken(internal_GPS_DOF) ) ! LOGICAL taken = .FALSE. ! but when one principal-axis is chosen, we will cross it off as unavailable in future ALLOCATE ( match(perturbations) ) ! INTEGER match = 0 ! whole list DO i = 1, perturbations IF (perturbation_sizes_mmpa(i) > 0.0D0) THEN r2Best = 999.9D99 jBest = 0 DO j = 1, internal_GPS_DOF ! considering all the available principal-axes: IF (.NOT.taken(j)) THEN !Try matching to this principal-axis: r2 = 0.0D0 ! just initializing the trial sum-of-squared-component-differences DO k = 1, internal_GPS_DOF r2 = r2 + (perturbation_vectors(k, i) - principal_axes(k, j))**2 !r2 = r2 + (perturbation_unit_vectors(k, i) - eigenvectors_block(k, j))**2 ! <== Gives the same result in test case, but perhaps more dangerous? END DO IF (r2 < r2Best) THEN r2Best = r2 jBest = j END IF !ALSO, try matching to the other end (negative) of this principal-axis: r2 = 0.0D0 ! just initializing the trial sum-of-squared-component-differences DO k = 1, internal_GPS_DOF r2 = r2 + (perturbation_vectors(k, i) + principal_axes(k, j))**2 !r2 = r2 + (perturbation_unit_vectors(k, i) + eigenvectors_block(k, j))**2 ! <== Gives the same result in test case, but perhaps more dangerous? END DO IF (r2 < r2Best) THEN r2Best = r2 jBest = j END IF END IF ! this principal-axis #j is still available END DO ! j = 1, internal_GPS_DOF; considering all the available principal-axes match(i) = jBest ! record the closest match found taken(jBest) = .TRUE. ! so no other perturbation can claim it END IF ! size of this unit-vector was positive END DO ! i = 1, perturbations !Assign each eigenvector/eigenvalue of internal_GP2 to one of three classes: ! (1) Reference-frame loosening (usually 3 eigenvectors); ! (2) Local noise sources (volcanos, post-seismic, opening crack, landslides, ...); ! (3) Other velocity misfits that are not understood or wanted. ALLOCATE ( eigen_class123(internal_GPS_DOF) ) ! INTEGER list eigen_class123 = 3 ! just initializing; a FEW of these will later be changed to 1 or 2. DO i = 1, perturbations IF (perturbation_sizes_mmpa(i) > 0.0D0) THEN eigen_class123(match(i)) = perturbation_class123(i) ! either 1 (reference-frame) or 2 (local noise source) ELSE ! this "noise source" produced no velocity perturbations at all END IF END DO !Divide unfit velocity residuals into: class1and2 versus class3, and report measures: ALLOCATE ( unfit_EN_class1and2(internal_GPS_DOF) ) unfit_EN_class1and2 = 0.0D0 ALLOCATE ( unfit_EN_class3(internal_GPS_DOF) ) unfit_EN_class3 = 0.0D0 DO i = 1, internal_GPS_DOF IF ((eigen_class123(i) == 1).OR.(eigen_class123(i) == 2)) THEN ! Note that we do not attempt to distinguish between these (at this stage) because ! often they have recombined in complex ways to create "mixed class1and2" eigenvectors of internal_GP2. ! (There will be an attempt to separate them later.) DO j = 1, internal_GPS_DOF unfit_EN_class1and2(j) = unfit_EN_class1and2(j) + unfit_eigenamplitudes(i) * eigenvectors_block(j, i) END DO ELSE IF (eigen_class123(i) == 3) THEN DO j = 1, internal_GPS_DOF unfit_EN_class3(j) = unfit_EN_class3(j) + unfit_eigenamplitudes(i) * eigenvectors_block(j, i) END DO END IF END DO !General header for class1and2 + class3 = unfit_supervector WRITE (*, *) WRITE (*, "(' Amplitude measures of various classes of')") WRITE (*, "(' actual post-NeoKinema unfit-velocity residuals:')") WRITE (*, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, *) WRITE (71, "('Amplitude measures of various classes of')") WRITE (71, "('actual post-NeoKinema unfit-velocity residuals:')") WRITE (71, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' -------------------------------------------------------------------------------')") !Class1and2: WRITE (*, "(' Classes 1 & 2 = Frame-loosening residuals & local noise residuals:')") WRITE (71, "('Classes 1 & 2 = Frame-loosening residuals & local noise residuals:')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN_class1and2(2 * j - 1)**2 + unfit_EN_class1and2(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") !Class3: WRITE (*, "(' Class 3 = Residuals controlled by original .GPS uncertainties:')") WRITE (71, "('Class 3 = Residuals controlled by original .GPS uncertainties:')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN_class3(2 * j - 1)**2 + unfit_EN_class3(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") !unfit_supervector WRITE (*, "(' Total unfit velocity residuals (Classes 1, 2, & 3):')") WRITE (71, "('Total unfit velocity residuals (Classes 1, 2, & 3):')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_supervector(2 * j - 1)**2 + unfit_supervector(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") CALL Pause() ! Attempt to separate the components of unfit_EN_class1and2 ! (which are hopelessly mixed together in the eigenanalysis of internal_GP2) ! by using the GPS_Covariance perturbation velocity vectors reported by the ! GPS_Covariance_log.txt file, and the least-squares solution of an ! over-determined linear system: ! ! perturbation_unit_vectors(1:internal_GPS_DOF, 1:perturbations) * amplitudes(1:perturbations) =~= unfit_EN_class1and2(1:internal_GPS_DOF) ! ! to get the desired amplitudes(1:perturbations). WRITE (*, *) WRITE (*, "(' Attempting to solve over-determined linear system for amplitudes of ',I4,' perturbations.')") perturbations WRITE (71, *) WRITE (71, "('Attempting to solve over-determined linear system for amplitudes of ',I4,' perturbations.')") perturbations trans = 'N' nrhs = 1 ALLOCATE ( copied_matrix(internal_GPS_DOF, perturbations) ) copied_matrix = perturbation_unit_vectors ! because the CALL below will overwrite this information! ALLOCATE ( copied_vector(internal_GPS_DOF) ) ! NOTE: Needs to have length m, not shorter length n, because it holds forcing vector at input time. copied_vector = unfit_EN_class1and2 ! because the CALL below will overwrite this information with the solution! lwork = MIN(internal_GPS_DOF, perturbations) + MAX(1, internal_GPS_DOF, perturbations, nrhs) ! See documentation for MKL/LAPACK routine gels (or dgels). ALLOCATE ( work(lwork) ) work = 0.0D0 ! (probably not necessary) CALL dgels(trans, internal_GPS_DOF, perturbations, nrhs, copied_matrix, internal_GPS_DOF, copied_vector, & & internal_GPS_DOF, work, lwork, info) ! Math Kernel Library (MKL) routine, from the Linear Analysis Package (LAPACK) section. IF (info /= 0) THEN WRITE (*, "(' ERROR: dgels of MKL could not solve overdetermined amplitudes problem by least-squares.'/' info = ',I12)") info WRITE (71, "('ERROR: dgels of MKL could not solve overdetermined amplitudes problem by least-squares.'/' info = ',I12)") info CALL Pause() STOP END IF DEALLOCATE ( work ) ! so that a new, different "work" vector can be created later for another CALL to MKL/LAPACK. ALLOCATE ( amplitudes(perturbations) ) amplitudes(1:perturbations) = copied_vector(1:perturbations) ! unpack solution values from location where MKL left them. !----------------------------------------------------------------------------------------------- !REPORT section: Give information on the amplitude of each perturbation added by GPS_Covariance ! in actual post-NeoKinema unfit velocity residuals. ALLOCATE ( unfit_EN(internal_GPS_DOF) ) unfit_EN = 0.0D0 WRITE (*, *) WRITE (*, "(' Amplitude measures of each perturbation added by GPS_Covariance')") WRITE (*, "(' in actual post-NeoKinema unfit-velocity residuals:')") WRITE (*, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, *) WRITE (71, "('Amplitude measures of each perturbation added by GPS_Covariance')") WRITE (71, "('in actual post-NeoKinema unfit-velocity residuals:')") WRITE (71, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (71, "('-------------------------------------------------------------------------------')") DO i = 1, perturbations IF (perturbation_sizes_mmpa(i) > 0.0D0) THEN WRITE (*, "(' ', A)") perturbation_type_c15(i) IF (perturbation_type_c15(i)(1:5) == "EXTRA") THEN WRITE (*, "(' ', A)") TRIM(perturbation_source(i)) ! more distinctive than title, which was filled in (here) with generic ELSE WRITE (*, "(' ', A)") TRIM(perturbation_title(i)) END IF WRITE (71, "(A)") perturbation_type_c15(i) IF (perturbation_type_c15(i)(1:5) == "EXTRA") THEN WRITE (71, "(A)") TRIM(perturbation_source(i)) ! more distinctive than title, which was filled in (here) with generic ELSE WRITE (71, "(A)") TRIM(perturbation_title(i)) END IF DO k = 1, internal_GPS_DOF unfit_EN(k) = amplitudes(i) * perturbation_unit_vectors(k, i) END DO measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN(2 * j - 1)**2 + unfit_EN(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' Amplitude of this perturbation = ',F7.3,' mm/a')") amplitudes(i) WRITE (*, "(' Amplitude in units of sigma (dimensionless) = ',F7.3)") amplitudes(i) / perturbation_sizes_mmpa(i) !WRITE (*, "(8(' ',10F7.3/))") (perturbation_unit_vectors(k, i), k = 1, internal_GPS_DOF) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('Amplitude of this perturbation = ',F7.3,' mm/a')") amplitudes(i) WRITE (71, "('Amplitude in units of sigma (dimensionless) = ',F7.3)") amplitudes(i) / perturbation_sizes_mmpa(i) !WRITE (71, "(8(10F7.3/))") (perturbation_unit_vectors(k, i), k = 1, internal_GPS_DOF) WRITE (71, "('-------------------------------------------------------------------------------')") END IF END DO ! i = 1, perturbations CALL Pause() !----------------------------------------------------------------------------------------------- !Give further statistics on the breakdown of unfit_EN_class1and2 = unfit_EN_class1 + unfit_EN_class2: ALLOCATE ( unfit_EN_class1(internal_GPS_DOF) ) unfit_EN_class1 = 0.0D0 ! just initializing before sum ALLOCATE ( unfit_EN_class2(internal_GPS_DOF) ) unfit_EN_class2 = 0.0D0 ! just initialing before sum DO i = 1, perturbations IF (perturbation_class123(i) == 1) THEN DO j = 1, internal_GPS_DOF unfit_EN_class1(j) = unfit_EN_class1(j) + amplitudes(i) * perturbation_unit_vectors(j, i) END DO ELSE IF (perturbation_class123(i) == 2) THEN DO j = 1, internal_GPS_DOF unfit_EN_class2(j) = unfit_EN_class2(j) + amplitudes(i) * perturbation_unit_vectors(j, i) END DO END IF END DO !General header for Classes 1, 2, & 3 = unfit_supervector WRITE (*, *) WRITE (*, "(' Amplitude measures of various classes of')") WRITE (*, "(' actual post-NeoKinema unfit-velocity residuals:')") WRITE (*, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, *) WRITE (71, "('Amplitude measures of various classes of')") WRITE (71, "('actual post-NeoKinema unfit-velocity residuals:')") WRITE (71, "(' m1 = mean velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m2 = RMS velocity magnitude across benchmarks, mm/a')") WRITE (71, "(' m3 = peak velocity magnitude across benchmarks, mm/a')") WRITE (71, "('-------------------------------------------------------------------------------')") !Class 1: WRITE (*, "(' Class 1 = Unfit velocity residuals attributed to reference frame error:')") WRITE (71, "('Class 1 = Unfit velocity residuals attributed to reference frame error:')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN_class1(2 * j - 1)**2 + unfit_EN_class1(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") !Class 2 WRITE (*, "(' Class 2 = Unfit velocity residuals attributed to local noise sources:')") WRITE (71, "('Class 2 = Unfit velocity residuals attributed to local noise sources:')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN_class2(2 * j - 1)**2 + unfit_EN_class2(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") !Class3: WRITE (*, "(' Class 3 = Residuals controlled by original .GPS uncertainties:')") WRITE (71, "('Class 3 = Residuals controlled by original .GPS uncertainties:')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_EN_class3(2 * j - 1)**2 + unfit_EN_class3(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") !unfit_supervector WRITE (*, "(' Total unfit GPS velocity residuals (Classes 1, 2, & 3):')") WRITE (71, "('Total unfit GPS velocity residuals (Classes 1, 2, & 3):')") measure = 0.0D0 ! just initializing, before scan of benchmarks, below... DO j = 1, internal_benchmarks velocity_magnitude = SQRT(unfit_supervector(2 * j - 1)**2 + unfit_supervector(2 * j)**2) measure(1) = measure(1) + velocity_magnitude measure(2) = measure(2) + velocity_magnitude**2 measure(3) = MAX(measure(3), velocity_magnitude) END DO measure(1) = measure(1) / internal_benchmarks measure(2) = SQRT(measure(2) / internal_benchmarks) WRITE (*, "(' m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (71, "('m1 = ',F7.3,' mm/a, m2 = ',F7.3,' mm/a, m3 = ',F7.3,' mm/a')") & & measure(1), measure(2), measure(3) WRITE (71, "('-------------------------------------------------------------------------------')") CALL Pause() !Create uncertainty-ellipse parameters for the 3 maps to be produced below ! unfit GPS velocity residuals, divided into Class 1, Class 2, and Class 3): ALLOCATE ( class1_vE_sigma(internal_benchmarks) ) ALLOCATE ( class1_vN_sigma(internal_benchmarks) ) ALLOCATE ( class1_correlation(internal_benchmarks) ) ALLOCATE ( class2_vE_sigma(internal_benchmarks) ) ALLOCATE ( class2_vN_sigma(internal_benchmarks) ) ALLOCATE ( class2_correlation(internal_benchmarks) ) !NOTE that class3 arrays were allocated earlier, and their contents were read from the g_{token}.nko file. ! (They could also have been read from the original .GPS file, but that list may have had extra benchmarks, not used by NeoKinema.) !Class 1 (reference-frame loosening): DO i = 1, internal_benchmarks c_EE = 0.0D0; c_EN = 0.0D0; c_NE = 0.0D0; c_NN = 0.0D0 ! initializing local 2x2 velocity-covariance matrix; (E, N), and (mm/a)**2. DO j = 1, perturbations IF (perturbation_class123(j) == 1) THEN c_EE = c_EE + perturbation_vectors(2*i-1, j) * perturbation_vectors(2*i-1, j) c_EN = c_EN + perturbation_vectors(2*i-1, j) * perturbation_vectors(2*i , j) c_NE = c_NE + perturbation_vectors(2*i , j) * perturbation_vectors(2*i-1, j) ! Note: c_NE = c_EN, always. c_NN = c_NN + perturbation_vectors(2*i , j) * perturbation_vectors(2*i , j) END IF ! this perturbation is Class-1 END DO ! all perturbations IF ((c_EE > 0.0D0).AND.(c_NN > 0.0D0)) THEN class1_vE_sigma(i) = SQRT(c_EE) class1_vN_sigma(i) = SQRT(c_NN) class1_correlation(i) = c_EN / (class1_vE_sigma(i) * class1_vN_sigma(i)) ! which could also be computed using c_NE instead of c_EN. ELSE ! avoid division-by-zero class1_vE_sigma(i) = 0.0D0 class1_vN_sigma(i) = 0.0D0 class1_correlation(i) = 0.0D0 END IF END DO ! all internal benchmarks !Class 2 (local noise sources, such as volcanos, post-seismic, landslides, ...): DO i = 1, internal_benchmarks c_EE = 0.0D0; c_EN = 0.0D0; c_NE = 0.0D0; c_NN = 0.0D0 ! initializing local 2x2 velocity-covariance matrix; (E, N), and (mm/a)**2. DO j = 1, perturbations IF (perturbation_class123(j) == 2) THEN c_EE = c_EE + perturbation_vectors(2*i-1, j) * perturbation_vectors(2*i-1, j) c_EN = c_EN + perturbation_vectors(2*i-1, j) * perturbation_vectors(2*i , j) c_NE = c_NE + perturbation_vectors(2*i , j) * perturbation_vectors(2*i-1, j) ! Note: c_NE = c_EN, always. c_NN = c_NN + perturbation_vectors(2*i , j) * perturbation_vectors(2*i , j) END IF ! this perturbation is Class-2 END DO ! all perturbations IF ((c_EE > 0.0D0).AND.(c_NN > 0.0D0)) THEN class2_vE_sigma(i) = SQRT(c_EE) class2_vN_sigma(i) = SQRT(c_NN) class2_correlation(i) = c_EN / (class2_vE_sigma(i) * class2_vN_sigma(i)) ! which could also be computed using c_NE instead of c_EN. ELSE ! avoid division-by-zero class2_vE_sigma(i) = 0.0D0 class2_vN_sigma(i) = 0.0D0 class2_correlation(i) = 0.0D0 END IF END DO ! all internal benchmarks !Output a table of results: WRITE (*, *) WRITE (*, "(' Parameters of uncertainty ellipses for internal benchmarks, by misfit class:')") WRITE (*, "(' CLASS-1: CLASS-2: CLASS-3: ')") WRITE (*, "(' Elon Nlat vEsig vNsig corre vEsig vNsig corre vEsig vNsig corre SITE')") WRITE (71, *) WRITE (71, "('Parameters of uncertainty ellipses for internal benchmarks, by misfit class:')") WRITE (71, "(' CLASS-1: CLASS-2: CLASS-3: ')") WRITE (71, "(' Elon Nlat vEsig vNsig corre vEsig vNsig corre vEsig vNsig corre SITE')") DO i = 1, internal_benchmarks WRITE (*, "(' ',F8.3,F7.3,9F6.2,' ',A)") internal_GPS_Elon(i), internal_GPS_Nlat(i), & & class1_vE_sigma(i), class1_vN_sigma(i), class1_correlation(i), & & class2_vE_sigma(i), class2_vN_sigma(i), class2_correlation(i), & & class3_vE_sigma(i), class3_vN_sigma(i), class3_correlation(i), & & TRIM(internal_GPS_site(i)) WRITE (71, "(F8.3,F7.3,9F6.2,' ',A)") internal_GPS_Elon(i), internal_GPS_Nlat(i), & & class1_vE_sigma(i), class1_vN_sigma(i), class1_correlation(i), & & class2_vE_sigma(i), class2_vN_sigma(i), class2_correlation(i), & & class3_vE_sigma(i), class3_vN_sigma(i), class3_correlation(i), & & TRIM(internal_GPS_site(i)) END DO WRITE (71, *) WRITE (71, "('Linear-algebra section completed.')") CLOSE (71) ! GPS_Postprocessor_log.txt WRITE (*, *) WRITE (*, "(' Linear-algebra section completed.')") CALL Pause() !============================================== BEGIN GRAPHICS SECTION =========================================== WRITE (*, *) WRITE (*, "(' ========================== BEGIN GRAPHICS SECTION ==============================')") WRITE (*, "(' The last part of this program uses code from NeoKineMap to produce 3 maps of')") WRITE (*, "(' the 3 components of the residual GPS velocity misfit in the NeoKinema solution.')") WRITE (*, *) WRITE (*, "(' For success, you need to have 2 more files available in the current folder:')") WRITE (*, "(' * AI7frame.ai (Adobe Illustrator template file, from Peter Bird)')") WRITE (*, "(' * Map_Tools.ini, copied from your NeoKineMap folder;')") WRITE (*, "(' this will give you the same map-projection as in other plots.')") WRITE (*, "(' PLEASE copy these files into the current folder NOW if they are not present.')") CALL Pause() WRITE (*, *) WRITE (*, "(' NOTE that each plot will begin with the unfit velocity residual vectors')") WRITE (*, "(' (and their associated benchmarks, and uncertainty ellipses),')") WRITE (*, "(' but that you will have a chance to add additional overlays,')") WRITE (*, "(' such as lines from a .DIG file, faults, and/or the outline of the .FEG.')") CALL Pause() !Initialize critical variables (which, in NeoKineMap, would have been read from NeoKineMap.ini). plot_dig_titles = .FALSE. dig_title_method = 1 grid_units = 'm' feg_file = ' ' lines_basemap_file = ' ' tick_points = 6.0D0 node_radius_points = 0.0D0 traces_file = ' ' benchmark_points = 4.0D0 velocity_Ma = 10.0D0 volcano_file = "Volcanoes.dat" !N.B. This refers to volcanoes-of-the-Earth, a potential map overlay, NOT to swelling-volcanoes used as GPS noise sources. volcano_points = 7.0D0 label_thinner = 1 minutes = 60 kilometers = 100 top_line_memo = ' ' bottom_line_memo = ' ' !This main, outer loop of the graphics section will produce 3 related maps, with slightly different names: ! g_unfit_residuals-class1.ai g_unfit_residuals-class2.ai g_unfit_residuals-class3.ai class_pictures: DO class_index = 1, 3 IF ((class_index == 1).AND.(.NOT.class1_exists)) CYCLE class_pictures IF ((class_index == 2).AND.(.NOT.class2_exists)) CYCLE class_pictures !State variables specific to each map page: title_count = 0 bottomlegend_used_points = 0.D0 ! records filling of bottom legend, from left rightlegend_used_points = 0.D0 ! records filling of right legend, from top WRITE (class_index_c1, "(I1)") class_index suggested_new_ai_filename = "g_unfit_residuals-class" // class_index_c1 // ".ai" WRITE (*, *) WRITE (*, "(' -------------------------------------------------------------------------------')") WRITE (*, "(' In next question, CHANGE the output .AI filename to ', A)") TRIM(suggested_new_ai_filename) WRITE (*, "(' -------------------------------------------------------------------------------')") CALL Pause() CALL DPrompter (xy_mode = .TRUE., lonlat_mode = .TRUE., path_out = path_out, & ! inputs & xy_defined = xy_defined) ! output; reports whether user set (x,y) system !NOTE: Prompter opens AI7Frame.ai, begins new graphics file. ! At this stage, we are ready to write on the page! got_parameters = .TRUE. ! p_{token)_.nki was already read-in; we know all about this NeoKinema model. !-------------------------- OVERLAYS ------------------------------ !----- (symbols composed mostly of lines; mostly transparent) ----- first_pass = .TRUE. ! So, the first type of overlay (8, type 5) will be pre-selected. !However, if the user loops-back to line 2000, this switch will then be .FALSE. 2000 WRITE (*,"(//' ----------------------------------------------------------------------')") IF (first_pass) THEN do_overlay = .TRUE. choice = 8 ! "geodetic benchmarks with velocities" first_pass = .FALSE. ! for next time (if any) ELSE ! User is returning to add additional overlay(s). WRITE (*,"( /' LINE AND SYMBOL OVERLAY LAYERS AVAILABLE:')") !GPBoverlays WRITE (*,"( ' 1 :: digitized basemap (lines type)')") WRITE (*,"( ' 2 :: outline of finite-element grid')") WRITE (*,"( ' 4 :: fault traces')") WRITE (*,"( ' 15 :: volcanoes (Recent, subaerial) from file volcanoes.dat')") WRITE (*,"( ' ----------------------------------------------------------------------')") CALL DPrompt_for_Logical('Do you want one (or more) of these overlays?', .TRUE., do_overlay) IF (do_overlay) THEN choice = 1 ! (just a suggestion) CALL DPrompt_for_Integer('Which overlay type should be added?', choice, choice) IF (.NOT.((choice == 1).OR.(choice == 2).OR.(choice == 4).OR.(choice == 15))) THEN WRITE (*,"(/' ERROR: Please select an integer from the list above!')") CALL DPress_Enter mt_flashby = .FALSE. GO TO 2000 END IF ! illegal or legal choice END IF ! do_overlay END IF ! first_pass, or .NOT. IF (do_overlay) THEN SELECT CASE (choice) CASE (1) ! basemap (lines type) 2010 temp_path_in = path_in !CALL File_List( file_type = "*.dig", & ! & suggested_file = lines_basemap_file, & ! & using_path = temp_path_in) CALL DPrompt_for_String('Which file should be plotted?', lines_basemap_file, lines_basemap_file) lines_basemap_pathfile = TRIM(temp_path_in)//TRIM(lines_basemap_file) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) CALL DSet_Stroke_Color (color_name = 'foreground') CALL Dig_Type (lines_basemap_pathfile, 21, dig_is_lonlat, any_titles) CALL DPrompt_for_Logical('Is this basemap written in (lon,lat) coordinates?', dig_is_lonlat, dig_is_lonlat) IF ((.NOT.xy_defined).AND.(.NOT.dig_is_lonlat)) THEN WRITE (*,"(' ERROR: Start GPS_Postprocessor over again and define the (x,y) coordinates.')") CALL DPress_Enter STOP ' ' END IF IF (any_titles) THEN WRITE (*, "(/' Title lines were detected in this .dig file;')") CALL DPrompt_for_Logical('do you want to include these titles in the plot?', plot_dig_titles, plot_dig_titles) IF (plot_dig_titles) THEN WRITE (*,"(' -------------------------------------------------')") WRITE (*,"(' Choose Alignment Method for Titles:')") WRITE (*,"(' 1: upright, at geometric center of polyline')") WRITE (*,"(' 2: parallel to first segment of polyline')") WRITE (*,"(' -------------------------------------------------')") 2011 CALL DPrompt_for_Integer('Which alignment method?', dig_title_method, dig_title_method) IF ((dig_title_method < 1).OR.(dig_title_method > 2)) THEN WRITE (*,"(' ERROR: Illegal choice of method.')") mt_flashby = .FALSE. GO TO 2011 END IF ! bad choice END IF ! plot_dig_titles END IF ! any_titles WRITE (*,"(/' Working on basemap....')") polygons = .FALSE. IF (dig_is_lonlat) THEN CALL DPlot_Dig (7, lines_basemap_pathfile, polygons, 21, in_ok) ELSE CALL DPlot_Dig (3, lines_basemap_pathfile, polygons, 21, in_ok) END IF IF (.NOT.in_ok) THEN mt_flashby = .FALSE. GO TO 2010 END IF IF (any_titles .AND. plot_dig_titles) THEN CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') IF (dig_is_lonlat) THEN CALL DPlot_Dig (7, lines_basemap_pathfile, polygons, 21, in_ok, dig_title_method) ELSE CALL DPlot_Dig (3, lines_basemap_pathfile, polygons, 21, in_ok, dig_title_method) END IF END IF ! any_titles .AND. plot_dig_titles WRITE (*,"('+Working on basemap....DONE.')") ! possible plot titles: filename, and first 3 lines CALL Add_Title(lines_basemap_file) OPEN (UNIT = 22, FILE = lines_basemap_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) CALL Could_Not_Find_File(lines_basemap_pathfile) DO i = 1, 3 READ (22, *, IOSTAT = ios) lon, lat IF (ios /= 0) THEN ! got possible title BACKSPACE (21) READ (22,"(A)") line CALL Add_Title(line) END IF END DO CLOSE (22) CALL BEEPQQ (440, 250) ! end of basemap overlay CASE (2) ! finite-element grid outline !Note: Code added to produce output file "outline.dig" which is needed by Long_Term_Seismicity_v1/2/3/4/5; ! however, in Long_Term_Seismicity_v6+ this functionality is built-in (by copying the following code). CALL Add_Title('Outline of Finite Element Grid') 2020 temp_path_in = path_in 2030 IF (got_parameters) THEN feg_file = x_feg ELSE !CALL File_List( file_type = "*.feg", & ! & suggested_file = feg_file, & ! & using_path = temp_path_in) CALL DPrompt_for_String('Which .FEG file should be outlined?', feg_file, feg_file) END IF feg_pathfile = TRIM(temp_path_in)//TRIM(feg_file) OPEN(UNIT = 22, FILE = feg_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22) CALL DPress_Enter mt_flashby = .FALSE. got_parameters = .FALSE. GO TO 2020 END IF READ (22,"(A)") line CALL Add_Title(line) READ (22,*) numnod ALLOCATE ( node_uvec(3,numnod), segments(3,2,numnod) ) DO i = 1, numnod READ (22,*) j, lon, lat CALL DLonLat_2_Uvec (lon, lat, uvec1) node_uvec(1:3,i) = uvec1(1:3) END DO ! i = 1, numnod READ (22,*) numel ALLOCATE ( nodes(3,numel) ) DO i = 1, numel READ (22,*) j, nodes(1,i), nodes(2,i), nodes(3,i) END DO ! i = 1, numel READ (22, *) nfl IF (nfl > 0) CALL No_Fault_Elements_Allowed() CLOSE(22) IF (choice == 2) THEN ! plot outline only WRITE (*,"(/' Working on outline of grid....')") CALL DBegin_Group() ! This might seem unecessary for a simply-connected perimeter; ! however, after L2-->L1 windowing it may no longer be one simply-connected piece; ! it may be chopped into bits where it leaves and re-enters the map window. ! Build library of external edge segments (unsorted) nseg = 0 DO i = 1, numel ! first, look for external edges of spherical triangles DO j = 1, 3 jp = 1 + MOD(j,3) na = nodes(j,i) nb = nodes(jp,i) mated = .FALSE. might_mate_1: DO k = 1, numel ! mating to another spherical triangle? DO l = 1, 3 lp = 1 + MOD(l,3) ma = nodes(l,k) mb = nodes(lp,k) IF ((na == mb).AND.(nb == ma)) THEN mated = .TRUE. EXIT might_mate_1 END IF ! mate was found END DO ! l = 1, 3 sides of trial spherical triangle END DO might_mate_1 ! k = 1, numel IF (.NOT.mated) THEN nseg = MIN(nseg + 1, numnod) ! no problem expected segments(1:3,1,nseg) = node_uvec(1:3,na) segments(1:3,2,nseg) = node_uvec(1:3,nb) !note that segments always proceed counterclockwise around grid END IF ! NOT mated END DO ! j = 1, 3 END DO ! i = 1, numel; looking for external edges of spherical triangles !link segments to create outline CALL DSet_Line_Style (width_points = 4.0D0, dashed = .FALSE.) CALL DSet_Stroke_Color ('gray______') j = 1 ! begin with first segment uvec1(1:3) = segments(1:3,1,j) CALL DNew_L45_Path (5, uvec1) OPEN (UNIT = 55, FILE = TRIM(path_out)//"outline_of_"//TRIM(feg_file)//".dig") WRITE (55, "('outline of ',A)") TRIM(feg_file) CALL DUvec_2_LonLat(uvec1, longitude, latitude) WRITE (55, "(1X,SP,ES12.5,',',ES12.5)") longitude, latitude DO i = 1, nseg uvec2(1:3) = segments(1:3,2,j) CALL DGreat_to_L45 (uvec2) CALL DUvec_2_LonLat(uvec2, longitude, latitude) WRITE (55, "(1X,SP,ES12.5,',',ES12.5)") longitude, latitude find_next: DO k = 2, nseg IF (uvec2(1) == segments(1,1,k)) THEN IF (uvec2(2) == segments(2,1,k)) THEN IF (uvec2(3) == segments(3,1,k)) THEN j = k EXIT find_next END IF END IF END IF END DO find_next !prepare to loop uvec1 = uvec2 END DO ! i = 1, nseg CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) WRITE (55, "('*** end of line segment ***')") CLOSE(55) CALL DEnd_Group() WRITE (*,"('+Working on outline of grid....DONE.')") END IF ! choice = 2 (outline) or 3 (nodes and elements) DEALLOCATE ( node_uvec, segments ) DEALLOCATE ( nodes ) CALL BEEPQQ (440, 250) ! end of (2) outline or (3) full finite element grid (which has no fault elements) CASE (4) ! fault traces CALL Add_Title("Fault Traces") 2040 temp_path_in = path_in IF (got_parameters) THEN traces_file = f_dig ELSE !don't require use of parameter file; it may not have been created yet. !CALL File_List(file_type = "f*.dig", & ! & suggested_file = traces_file, & ! & using_path = temp_path_in) CALL DPrompt_for_String('Which file should be plotted?',traces_file,traces_file) END IF traces_pathfile = TRIM(temp_path_in)//TRIM(traces_file) OPEN(UNIT = 21, FILE = traces_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (21) WRITE (*,"(' Press any key to continue...'\)") READ (*,"(A)") line got_parameters = .FALSE. GOTO 2040 END IF CLOSE (21) CALL DPrompt_for_Real('How large are the dip ticks (in points)?', tick_points, tick_points) WRITE (*,"(/' Working on fault traces....')") CALL Fault_Traces (trace_choice = 0) ! N.B. This subprogram will read "dip_degrees" if present. WRITE (*,"('+Working on fault traces....DONE.')") CALL Add_Title(traces_file) CALL BEEPQQ (440, 250) CALL Chooser (bottom, right) IF (right) THEN CALL DBegin_Group CALL Fault_Key_Right (rightlegend_gap_points, rightlegend_used_points, tick_points) CALL DEnd_Group ELSE IF (bottom) THEN CALL DBegin_Group CALL Fault_Key_Bottom (bottomlegend_gap_points, bottomlegend_used_points, tick_points) CALL DEnd_Group END IF ! end of (4) fault traces overlay CASE (8) ! GPS benchmarks with velocities 2080 temp_path_in = path_in !WRITE (*,"(/' Which kind of geodetic velocities do you wish to plot?')") !WRITE (*,"( ' (5) unfit portion of GPS velocities (data - model residuals)')") gps_type = 5 ! No choice is given to the user. ALSO, all source code related to the other 4 options has been deleted! temp_title = "Unfit GPS Velocities (data-model residuals) in model " // TRIM(trimmed_token) CALL Add_Title(temp_title) IF (class_index == 1) THEN CALL Add_Title("Class 1: Inconsistency of velocity reference frames for GPS vs. VBCs") ELSE IF (class_index == 2) THEN CALL Add_Title("Class 2: Local sources of non-stationary or non-tectonic noise") ELSE IF (class_index == 3) THEN CALL Add_Title("Class 3: Remaining unexplained (and undesirable) residuals") END IF benchmarks = internal_benchmarks ! We will only plot those that were used in the NeoKinema solution. ALLOCATE ( benchmark_label (benchmarks) ) ALLOCATE ( benchmark_uvec (3,benchmarks) ) ALLOCATE ( benchmark_N_velocity (benchmarks) ) ALLOCATE ( benchmark_N_sigma (benchmarks) ) ALLOCATE ( benchmark_E_velocity (benchmarks) ) ALLOCATE ( benchmark_E_sigma (benchmarks) ) ALLOCATE ( benchmark_correlation (benchmarks) ) ALLOCATE ( benchmark_hypotenuse (benchmarks) ) ! Pythagorean sum of dV_East and dV_North components, for convenience DO i = 1, benchmarks benchmark_label(i) = internal_GPS_site(i) longitude = internal_GPS_Elon(i) latitude = internal_GPS_Nlat(i) CALL DLonLat_2_Uvec(longitude, latitude, uvec) benchmark_uvec(1:3, i) = uvec(1:3) IF (class_index == 1) THEN benchmark_E_velocity(i) = unfit_EN_class1(2*i-1) benchmark_N_velocity(i) = unfit_EN_class1(2*i ) benchmark_E_sigma(i) = class1_vE_sigma(i) benchmark_N_sigma(i) = class1_vN_sigma(i) benchmark_correlation(i) = class1_correlation(i) ELSE IF (class_index == 2) THEN benchmark_E_velocity(i) = unfit_EN_class2(2*i-1) benchmark_N_velocity(i) = unfit_EN_class2(2*i ) benchmark_E_sigma(i) = class2_vE_sigma(i) benchmark_N_sigma(i) = class2_vN_sigma(i) benchmark_correlation(i) = class2_correlation(i) ELSE IF (class_index == 3) THEN benchmark_E_velocity(i) = unfit_EN_class3(2*i-1) benchmark_N_velocity(i) = unfit_EN_class3(2*i ) benchmark_E_sigma(i) = class3_vE_sigma(i) benchmark_N_sigma(i) = class3_vN_sigma(i) benchmark_correlation(i) = class3_correlation(i) END IF !GPBhere benchmark_hypotenuse(i) = SQRT(benchmark_E_velocity(i)**2 + benchmark_N_velocity(i)**2) END DO ALLOCATE ( train (benchmarks) ) k = 0 DO i = 1, benchmarks uvec(1:3) = benchmark_uvec(1:3, i) visible = DL5_In_Window(uvec) IF (visible) THEN k = k + 1 train(k) = benchmark_hypotenuse(i) END IF END DO WRITE (*,"(/' Here is the distribution of visible misfit velocities (in mm/a):')") CALL Histogram (train, k, .FALSE., maximum, minimum) DEALLOCATE ( train ) CALL DPrompt_for_Real('For how many Ma should velocity misfits be projected?', velocity_Ma, velocity_Ma) CALL DPrompt_for_Real('How large (in points) should benchmark locations be plotted?', benchmark_points, benchmark_points) WRITE (*,"(/' Working on benchmark velocity vectors....')") !-------------------------------------------------------------------- !create group of error ellipses: ellipses = .FALSE. ! usually reversed by any finite ellipse, below CALL DSet_Join_to_Round CALL DSet_Line_Style (width_points = 0.5D0, dashed = .FALSE.) CALL DSet_Stroke_Color ('foreground') t = (velocity_Ma * 1.0D6) * 0.001D0 / mp_radius_meters ! arc-radians per (mm/a) IF (velocity_Ma /= 0.0D0) THEN CALL DBegin_Group DO i = 1, benchmarks IF ((benchmark_N_sigma(i) > 0.0D0).AND.(benchmark_E_sigma(i) > 0.0D0)) THEN ! Maybe I should write ".OR." ? ellipses = .TRUE. uvec(1:3) = benchmark_uvec(1:3,i) t1 = t * DConformal_Deflation (uvec) ! arc-radians per (mm/a) ! (Benchmark location will be the center of the ellipse.) ! Find rotations of principal axes of ellipse: covariance_11 = benchmark_E_sigma(i)**2 covariance_22 = benchmark_N_sigma(i)**2 covariance_12 = benchmark_N_sigma(i) * benchmark_E_sigma(i) * benchmark_correlation(i) CALL DPrincipal_Axes_22 (covariance_11, covariance_12, covariance_22, & & e1, e2, u1x,u1y, u2x,u2y) e1 = MAX(e1, 0.0D0) ! prevent SQRT(<0.0) {probably just very slightly, due to round-off errors} e2 = MAX(e2, 0.0D0) ! prevent SQRT(<0.0) {probably just very slightly, due to round-off errors} e1 = 1.96D0 * DSQRT(e1) e2 = 1.96D0 * DSQRT(e2) ! in units of mm/a, but now amplified by *1.96, to convert from 1-sigma to 95%-confidence start_azimuth = Pi_over_2 - DATan2F(u1y, u1x) ! smallest axis, in radians clockwise from North !find initial point at top of ellipse: CALL DTurn_To (azimuth_radians = start_azimuth, base_uvec = uvec, far_radians = t1 * e1, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL DNew_L45_Path (5, uvec1) ! beginning (e1 axis) of ellipse DO j = 1, 24 ! 24 15-degree sectors, counterclockwise from e1 axis: rel_az2 = -(j - 0.5D0) * 15.0D0 * radians_per_degree ! mid-point; relative to e1 axis rel_az3 = -j * 15.0D0 * radians_per_degree ! end-point; relative to e1 axis az2 = start_azimuth + rel_az2 ! mid-point, in radians clockwise from N az3 = start_azimuth + rel_az3 ! end-point, in radians clockwise from N ds2 = DCOS(rel_az2) * t1 * e1 ! arc-radians dl2 = DSIN(rel_az2) * t1 * e2 arc2 = DSQRT(ds2**2 + dl2**2) aze2 = start_azimuth + DATan2F(dl2,ds2) CALL DTurn_To (azimuth_radians = aze2, base_uvec = uvec, far_radians = arc2, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) ds3 = DCOS(rel_az3) * t1 * e1 ! arc-radians dl3 = DSIN(rel_az3) * t1 * e2 arc3 = DSQRT(ds3**2 + dl3**2) aze3 = start_azimuth + DATan2F(dl3,ds3) CALL DTurn_To (azimuth_radians = aze3, base_uvec = uvec, far_radians = arc3, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) CALL DSmall_Through_L45 (uvec2, uvec3) ! through uvec2 to uvec3 END DO ! j = 1, 12 ! 30-degree sectors forming a circle CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END IF ! ellipise has positive dimensions END DO ! i = 1, benchmarks CALL DEnd_Group ! of error ellipses END IF ! velocity_Ma /= 0.0D0 !-------------------------------------------------------------------- !create group of benchmarks: IF (benchmark_points > 0.0D0) THEN CALL DSet_Join_to_Mitre CALL DSet_Line_Style (width_points = 0.5D0, dashed = .FALSE.) CALL DSet_Stroke_Color ('foreground') t = 0.6667D0 * mp_meters_per_point * benchmark_points / mp_radius_meters CALL DBegin_Group ! of benchmark triangles DO i = 1, benchmarks uvec(1:3) = benchmark_uvec(1:3,i) t1 = t * DConformal_Deflation (uvec) CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, far_radians = t1, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec1) CALL DNew_L45_Path (5, uvec1) CALL DTurn_To (azimuth_radians = 4.188D0, base_uvec = uvec, far_radians = t1, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec2) CALL DGreat_To_L45 (uvec2) CALL DTurn_To (azimuth_radians = 2.094D0, base_uvec = uvec, far_radians = t1, & ! inputs & omega_uvec = omega_uvec, result_uvec = uvec3) CALL DGreat_To_L45 (uvec3) CALL DGreat_To_L45 (uvec1) CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END DO ! i = 1, benchmarks CALL DEnd_Group ! of benchmark triangles END IF ! benchmark_points > 0.0 !-------------------------------------------------------------------- !create group of benchmark labels: CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DBegin_Group ! of benchmark labels DO i = 1, benchmarks uvec(1:3) = benchmark_uvec(1:3,i) c4 = benchmark_label(i)(1:4) CALL DL5_Text (uvec = uvec, angle_radians = 0.0D0, from_east = .TRUE., & & font_points = 6, lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = c4) END DO ! i = 1, benchmarks CALL DEnd_Group ! of benchmark triangles !-------------------------------------------------------------------- !create group of velocity vectors: IF (velocity_Ma /= 0.0D0) THEN CALL DSet_Join_to_Mitre CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DBegin_Group CALL DSet_Stroke_Color ('foreground') DO i = 1, benchmarks uvec1(1:3) = benchmark_uvec(1:3,i) v_South_mps = -0.001D0 * benchmark_N_velocity(i) / s_per_year v_East_mps = +0.001D0 * benchmark_E_velocity(i) / s_per_year CALL DVelocity_Vector_on_Sphere (from_uvec = uvec1, & & v_theta_mps = v_South_mps, v_phi_mps = v_East_mps, & & dt_sec = velocity_Ma * 1.D6 * s_per_year, deflate = .TRUE.) END DO ! actually plotting benchmark velocity vectors CALL DEnd_Group END IF ! velocity_Ma /= 0.0D0 !-------------------------------------------------------------------- IF (ai_using_color) CALL DSet_Stroke_Color ('foreground') CALL Chooser(bottom, right) IF (right) THEN CALL DReport_RightLegend_Frame (x1_points, x2_points, y1_points, y2_points) y2_points = y2_points - rightlegend_used_points - rightlegend_gap_points CALL DBegin_Group CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') IF (gps_type == 1) THEN CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = 'geodetic') ELSE IF (gps_type == 2) THEN CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = 'adjustment to') ELSE IF (gps_type == 3) THEN CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = 'adjusted') ELSE IF (gps_type == 4) THEN CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = 'NeoKinema') END IF CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points - 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = 'velocity') number8 = ADJUSTL(DASCII8(velocity_Ma)) CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points - 24.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = '(x '//TRIM(number8)//' Ma):') x0p = 0.5D0 * (x1_points + x2_points) - 14.17D0 x1p = x0p + 2 * 14.17D0 ! 1-cm-long vector, expressed in points y0p = y2_points - 47.0D0 CALL DSet_Line_Style (width_points = 1.0D0, dashed = .FALSE.) CALL DSet_Stroke_Color ('foreground') IF (ellipses) THEN CALL DCircle_on_L12 (level = 1, x = x0p, y = y0p, radius = 6.0D0, stroke = .TRUE., fill = .FALSE.) END IF IF (benchmark_points > 0.0D0) THEN CALL DNew_L12_Path (1, x0p, y0p + 0.6667D0 * benchmark_points) CALL DLine_to_L12 (x0p - 0.577D0 * benchmark_points, y0p - 0.333D0 * benchmark_points) CALL DLine_to_L12 (x0p + 0.577D0 * benchmark_points, y0p - 0.333D0 * benchmark_points) CALL DLine_to_L12 (x0p, y0p + 0.6667D0 * benchmark_points) CALL DEnd_L12_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END IF CALL DSet_Line_Style (width_points = 1.5D0, dashed = .FALSE.) IF (ai_using_color) THEN IF (gps_type == 1) THEN ! original *.gps file CALL DSet_Stroke_Color ('foreground') ELSE IF (gps_type == 2) THEN ! adjustments to velocity CALL DSet_Stroke_Color ('gray______') ELSE IF (gps_type == 3) THEN ! modified g*.nko file CALL DSet_Stroke_Color ('dark_blue_') ELSE IF (gps_type == 4) THEN ! color will vary, green -> bronze -> red with error (in sigmas) CALL DSet_Stroke_Color ('red_______') END IF ELSE CALL DSet_Stroke_Color ('foreground') END IF CALL DVector_in_Plane (level = 1, & ! 1-cm-long vector & from_x = x0p, from_y = y0p, & & to_x = x1p, to_y = y0p) IF (ai_using_color) CALL DSet_Stroke_Color ('foreground') v_mps = 0.01D0 * mp_scale_denominator / (velocity_Ma * 1.D6 * s_per_year) v_mma = v_mps * 1000.D0 * s_per_year number8 = ADJUSTL(DASCII8(v_mma)) CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points - 48.D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = TRIM(number8)//' mm/a') IF (ellipses) THEN CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points - 60.D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = "(95%-c.") CALL DL12_Text (level = 1, & & x_points = 0.5D0*(x1_points + x2_points), & & y_points = y2_points - 72.D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 1.0D0, & & text = "ellipse)") END IF ! ellipses CALL DEnd_Group rightlegend_used_points = rightlegend_used_points + rightlegend_gap_points + 60.0D0 IF (ellipses) rightlegend_used_points = rightlegend_used_points + 24.0D0 ! for "(95%-c./ellipse)" ELSE IF (bottom) THEN CALL DReport_BottomLegend_Frame (x1_points, x2_points, y1_points, y2_points) x1_points = x1_points + bottomlegend_used_points + bottomlegend_gap_points CALL DBegin_Group CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') IF (gps_type == 1) THEN CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) + 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'geodetic') ELSE IF (gps_type == 2) THEN CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) + 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'adjustment to') ELSE IF (gps_type == 3) THEN CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) + 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'adjusted') ELSE IF (gps_type == 4) THEN CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) + 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'Neokinema') END IF CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points), & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'velocity') number8 = ADJUSTL(DASCII8(velocity_Ma)) CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) - 12.0D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = '(x '//TRIM(number8)//' Ma):') x0p = (x1_points + 29.0D0) - 14.17D0 x1p = x0p + 2 * 14.17D0 ! 1-cm-long vector, expressed in points y0p = 0.5D0 * (y1_points + y2_points) - 22.0D0 CALL DSet_Line_Style (width_points = 1.0D0, dashed = .FALSE.) CALL DSet_Stroke_Color ('foreground') IF (ellipses) THEN CALL DCircle_on_L12 (level = 1, x = x0p, y = y0p, radius = 6.0D0, stroke = .TRUE., fill = .FALSE.) CALL DL12_Text (level = 1, & & x_points = x1p + 9.D0, & & y_points = y0p, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.0D0, ud_fraction = 0.4D0, & & text = '95%-c.') END IF IF (benchmark_points > 0.0D0) THEN CALL DNew_L12_Path (1, x0p, y0p + 0.6667D0 * benchmark_points) CALL DLine_to_L12 (x0p - 0.577D0 * benchmark_points, y0p - 0.333D0 * benchmark_points) CALL DLine_to_L12 (x0p + 0.577D0 * benchmark_points, y0p - 0.333D0 * benchmark_points) CALL DLine_to_L12 (x0p, y0p + 0.6667D0 * benchmark_points) CALL DEnd_L12_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END IF CALL DSet_Line_Style (width_points = 1.5D0, dashed = .FALSE.) IF (ai_using_color) THEN IF (gps_type == 1) THEN ! original *.gps file CALL DSet_Stroke_Color ('foreground') ELSE IF (gps_type == 2) THEN ! adjustments to velocity CALL DSet_Stroke_Color ('magenta___') ELSE IF (gps_type == 3) THEN ! modified g*.nko file CALL DSet_Stroke_Color ('dark_blue_') ELSE IF (gps_type == 4) THEN ! color will vary, green -> bronze -> red with error (in sigmas) CALL DSet_Stroke_Color ('red_______') END IF ELSE CALL DSet_Stroke_Color ('foreground') END IF CALL DVector_in_Plane (level = 1, & ! 1-cm-long vector & from_x = x0p, from_y = y0p, & & to_x = x1p, to_y = y0p) IF (ai_using_color) CALL DSet_Stroke_Color ('foreground') v_mps = 0.01D0 * mp_scale_denominator / (velocity_Ma * 1.D6 * s_per_year) v_mma = v_mps * 1000.D0 * s_per_year number8 = ADJUSTL(DASCII8(v_mma)) CALL DL12_Text (level = 1, & & x_points = x1_points + 29.D0, & & y_points = 0.5D0*(y1_points + y2_points) - 36.D0, & & angle_radians = 0.D0, & & font_points = 12, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = TRIM(number8)//' mm/a') CALL DEnd_Group bottomlegend_used_points = bottomlegend_used_points + bottomlegend_gap_points + 58.D0 IF (ellipses) bottomlegend_used_points = bottomlegend_used_points + 36.D0 ! for "95%-c." END IF ! bottom or right legend WRITE (*,"('+Working on benchmark velocity vectors....DONE.')") DEALLOCATE ( benchmark_hypotenuse ) ! in LIFO order DEALLOCATE ( benchmark_correlation ) DEALLOCATE ( benchmark_E_sigma ) DEALLOCATE ( benchmark_E_velocity ) DEALLOCATE ( benchmark_N_sigma ) DEALLOCATE ( benchmark_N_velocity ) DEALLOCATE ( benchmark_uvec ) DEALLOCATE ( benchmark_label ) CALL BEEPQQ (440, 250) WRITE (*, *) WRITE (*, "(' -------------------------------------------------------------------------')") WRITE (*, "(' Class ', A, ' of unfit GPS velocity residuals have been plotted.')") class_index_c1 WRITE (*, "(' -------------------------------------------------------------------------')") ! end of 8: geodetic velocities of benchmarks !GPBhere CASE (15) ! volcanoes of the Earth (an overlay of cone-icons) 2150 temp_path_in = path_in !CALL File_List( file_type = "*.*", & ! & suggested_file = volcano_file, & ! & using_path = temp_path_in) 2151 WRITE (*,*) CALL DPrompt_for_String ('Which file has the volcano locations?',volcano_file,volcano_file) IF (LEN_TRIM(temp_path_in) > 0) THEN volcano_pathfile = TRIM(temp_path_in) // TRIM(volcano_file) ELSE volcano_pathfile = TRIM(volcano_file) END IF OPEN (UNIT = 22, FILE = volcano_pathfile, STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22, IOSTAT = ios) CALL DPress_Enter mt_flashby = .FALSE. GO TO 2151 END IF CALL Add_Title ('Volcanoes') IF (TRIM(volcano_file) == 'Volcanoes.dat') CALL Add_Title & & ("Smithsonian Institution, Global Volcanism Project") WRITE (*,*) CALL DPrompt_for_Real('How many points high shall symbols be?', volcano_points, volcano_points) ! Vents have white fill (snow) with black outlines (to separate them), ! unless map is b/w, in which case they are gray. IF (ai_using_color) THEN CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='background') ELSE ! b/w figure CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='gray______') END IF CALL DSet_Stroke_Color ('foreground') CALL DSet_Line_Style (0.6D0, .FALSE.) WRITE (*,"(/' Working on volcanoes....')") CALL DBegin_Group volcano_reading: DO READ (22, "(A)", IOSTAT = ios) line IF (ios == -1) EXIT volcano_reading READ (line, "(61X,F6.3,1X,A1,1X,F7.3,1X,A1)", IOSTAT = ios) & & volcano_Nlat, cN, volcano_Elon, cE IF (ios /= 0) THEN WRITE (*,"(/' ERROR: Bad line of data in file ',A,':')") TRIM(volcano_file) WRITE (*,"(' ',A)") TRIM(line) STOP END IF IF (cN == 'S') volcano_Nlat = -volcano_Nlat IF (cE == 'W') volcano_Elon = -volcano_Elon CALL DLonLat_2_Uvec (volcano_Elon, volcano_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) leg = 1.1547D0 * volcano_points * radians_per_point rad = 0.6666D0 * volcano_points * radians_per_point CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = rad, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) uvec(1:3) = result_uvec(1:3) CALL DNew_L45_Path (5, uvec) CALL DTurn_To (azimuth_radians = 3.665D0, base_uvec = uvec, & & far_radians = leg, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DGreat_to_L45 (result_uvec) CALL DTurn_To (azimuth_radians = 2.618D0, base_uvec = uvec, & & far_radians = leg, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DGreat_to_L45 (result_uvec) CALL DGreat_to_L45 (uvec) CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .TRUE.) END DO volcano_reading CALL DEnd_Group CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CLOSE(22) CALL Chooser (bottom, right) IF (bottom) THEN CALL DReport_BottomLegend_Frame (x1_points, x2_points, y1_points, y2_points) CALL DBegin_Group xp = x1_points + bottomlegend_used_points + bottomlegend_gap_points + 23.0D0 yp = (y1_points + y2_points)/2.0D0 + 12.0D0 + volcano_points IF (ai_using_color) THEN CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='background') ELSE ! b/w figure CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='gray______') END IF CALL DSet_Stroke_Color ('foreground') CALL DSet_Line_Style (0.6D0, .FALSE.) CALL DNew_L12_Path(1, xp, yp) xpt = xp - 0.57735D0 * volcano_points ypt = yp - volcano_points CALL DLine_To_L12 (xpt, ypt) xpt = xp + 0.57735D0 * volcano_points CALL DLine_To_L12 (xpt, ypt) CALL DLine_To_L12 (xp, yp) CALL DEnd_L12_Path (close = .TRUE., stroke = .TRUE., fill = .TRUE.) ypt = ypt - 12.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'Recent') ypt = ypt - 12.0D0 CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'volcano') CALL DEnd_Group bottomlegend_used_points = bottomlegend_used_points + bottomlegend_gap_points + 46.0D0 ELSE IF (right) THEN CALL DReport_RightLegend_Frame (x1_points, x2_points, y1_points, y2_points) CALL DBegin_Group xp = (x1_points + x2_points)/2.0D0 yp = y2_points - rightlegend_used_points + rightlegend_gap_points IF (ai_using_color) THEN CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='background') ELSE ! b/w figure CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='gray______') END IF CALL DSet_Stroke_Color ('foreground') CALL DSet_Line_Style (0.6D0, .FALSE.) CALL DNew_L12_Path(1, xp, yp) xpt = xp - 0.57735D0 * volcano_points ypt = yp - volcano_points CALL DLine_To_L12 (xpt, ypt) xpt = xp + 0.57735D0 * volcano_points CALL DLine_To_L12 (xpt, ypt) CALL DLine_To_L12 (xp, yp) CALL DEnd_L12_Path (close = .TRUE., stroke = .TRUE., fill = .TRUE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') ypt = ypt - 12.0D0 CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'Recent') ypt = ypt - 12.0D0 CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'volcano') CALL DEnd_Group rightlegend_used_points = rightlegend_used_points + rightlegend_gap_points + volcano_points + 24.0D0 END IF ! bottom or right legend WRITE (*,"('+Working on volcanoes....DONE.')") CALL BEEPQQ (440, 250) ! end of 15: volcanoes END SELECT ! (choice) = overlay type WRITE (*, *) CALL DPrompt_for_Logical('Do you want additional overlays?', .TRUE., do_more_overlays) IF (do_more_overlays) GO TO 2000 END IF ! do overlay !-------------------------------------------------------------------- !Graticule of parallels and meridians CALL DSet_Line_Style (width_points = 0.25D0, dashed = .FALSE.) CALL DSet_Stroke_Color (color_name = 'foreground') WRITE (*,"(' ')") IF (mp_projection_number == 0) THEN ! (x,y) axes desired 3010 CALL DPrompt_for_Integer('How many kilometers apart should fiducial lines& & of constant x and constant y be plotted?', kilometers, kilometers) IF (kilometers < 1) THEN WRITE (*, "(' ERROR: This value must be an integer >= 1')") CALL DPress_Enter mt_flashby = .FALSE. GO TO 3010 END IF CALL DWire_Mesh (kilometers) ELSE ! parallels and meridians desired 3020 CALL DPrompt_for_Integer('How many minutes apart should parallels& & and meridians be plotted?', minutes, minutes) IF (minutes < 1) THEN WRITE (*, "(' ERROR: This value must be an integer >= 1')") CALL DPress_Enter mt_flashby = .FALSE. GO TO 3020 END IF CALL DGraticule (minutes) END IF !numbered margin: IF (mp_projection_number == 0) THEN ! (x,y) axes desired CALL DKilometer_Frame (kilometers) ELSE ! parallels and meridians desired CALL DLonLat_Frame (minutes) END IF !Titles at top of map IF (ai_toptitles_reserved) THEN WRITE (*,"(' ')") mt_flashby = .FALSE. ! Do NOT flash by the prompts for titles, if there is space! CALL DPrompt_for_Logical('Do you want to add a title to this map?',.TRUE.,add_titles) IF (add_titles) THEN CALL Add_Title(top_line_memo) CALL Add_Title(bottom_line_memo) 900 WRITE (*,"(/' ----------------------------------------------------------------------')") WRITE (*,"(' SOME SUGGESTED TITLE OPTIONS')") WRITE (*,"(' (culled from files opened for this map)')") WRITE (*,"(/' 0 :: ANYTHING YOU CHOOSE TO TYPE!')") DO i = 1, title_count WRITE (*,"(' ',I2,' :: ',A)") i, TRIM(titles(i)) END DO ! i = 1, title_count WRITE (*,"(' ----------------------------------------------------------------------')") CALL DPrompt_for_Integer('Which option do you want for the upper line?',0,title_choice) IF ((title_choice < 0).OR.(title_choice > title_count)) THEN WRITE (*,"(' ERROR: Please choose a number from the list.')") mt_flashby = .FALSE. GO TO 900 END IF ! illegal choice IF (title_choice == 0) THEN CALL DPrompt_for_String('Enter top title (or one space for none)',' ',top_line) top_line_memo = top_line ELSE ! selection from list top_line = TRIM(titles(title_choice)) END IF CALL DPrompt_for_Integer('Which option do you want for the lower line?', 0, title_choice) IF ((title_choice < 0).OR.(title_choice > title_count)) THEN WRITE (*,"(' ERROR: Please choose a number from the list.')") mt_flashby = .FALSE. GO TO 900 END IF ! illegal choice IF (title_choice == 0) THEN CALL DPrompt_for_String('Enter sub-title (or one space for none)',' ',bottom_line) bottom_line_memo = bottom_line ELSE ! selection from list bottom_line = TRIM(titles(title_choice)) END IF CALL DTop_Titles (top_line, & & bottom_line) END IF ! add_titles END IF ! ai_toptitles_reserved CALL DEnd_Page END DO class_pictures ! class_index = 1, 3 (for 3 maps) WRITE (*, *) WRITE (*, "(' Graphics section completed.')") CALL Pause() CONTAINS !GPBSUBRs SUBROUTINE Add_Title(line) ! Adds "line" to global array "titles" and bumps global "title_count" ! if "line" is non-blank and also is novel. CHARACTER*(*), INTENT(IN) :: line CHARACTER*132 :: copy LOGICAL :: blank, novel INTEGER :: i blank = LEN_TRIM(line) <= 0 IF (.NOT.blank) THEN copy = ADJUSTL(TRIM(line)) novel = .TRUE. IF (title_count > 0) THEN DO i = 1, title_count IF (TRIM(copy) == TRIM(titles(i))) novel = .FALSE. END DO ! i = 1, title_count END IF ! have stored titles already IF (novel) THEN title_count = MIN(20, title_count + 1) titles(title_count) = TRIM(copy) END IF ! novel END IF ! not blank END SUBROUTINE Add_Title SUBROUTINE Bad_Parameters() !prints admonition screen, calls Pause(), and STOPs. 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 (*, "(' The number (and ordering) of parameters changed between NeoKinema v1.x')") WRITE (*, "(' and NeoKinema v2.x~3.x; then changed again with NeoKinema v4.x.')") WRITE (*, "(' ')") WRITE (*, "(' Errors may occur when input parameter files prepared for NeoKinema v1.x,')") WRITE (*, "(' 2.x, or 3.x are read by this program, NeoKinema v.4, which has different')") WRITE (*, "(' parameters and/or parameter ordering.')") WRITE (*, "(' ')") WRITE (*, "(' After you have read this message, GPS_Postprocessor will stop.')") WRITE (*, "(' Please correct the input parameter file and re-start it.')") WRITE (*, "(' ==========================================================================')") Call Pause() STOP END SUBROUTINE Bad_Parameters SUBROUTINE Chooser(bottom, right) ! Decides whether there is more margin space at "bottom" or "right". ! Will return both = F if NOT (ai_bottomlegend_reserved OR ai_rightlegend_reserved). ! Refers to NeoKineMap global variables: bottomlegend_used_points, rightlegend_used_points. LOGICAL, INTENT(out) :: bottom, right REAL*8 :: bottomlegend_free_points, rightlegend_free_points, & & x1_points, x2_points, y1_points, y2_points bottom = ai_bottomlegend_reserved right = ai_rightlegend_reserved IF (bottom.AND.right) THEN ! must choose one CALL DReport_BottomLegend_Frame (x1_points, x2_points, y1_points, y2_points) bottomlegend_free_points = x2_points - x1_points - bottomlegend_used_points CALL DReport_RightLegend_Frame (x1_points, x2_points, y1_points, y2_points) rightlegend_free_points = y2_points - y1_points - rightlegend_used_points IF (rightlegend_free_points >= bottomlegend_free_points) THEN right = .TRUE. bottom = .FALSE. ELSE right = .FALSE. bottom = .TRUE. END IF END IF ! choice is needed END SUBROUTINE Chooser 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 Dig_Type (dig_pathfile, free_unit, dig_is_lonlat, any_titles) ! Decide whether dig_pathfile is (lon,lat) or (x,y) based ! on the extreme range displayed in the y (or latitude) ! component. ! Also reports "any_titles" = T/F. ! Note that there can be trouble when a title like "TX" is ! interpreted by (*) format-free READ as x, and then y is ! taken from the start of the next line (a longitude!). ! So, we have to test the first two bytes to rule out titles. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: dig_pathfile ! points to .dig file INTEGER, INTENT(IN) :: free_unit ! Fortran device number LOGICAL, INTENT(OUT) :: dig_is_lonlat, any_titles ! Yes or No CHARACTER*2 :: c2 CHARACTER*26 :: line INTEGER :: ios LOGICAL :: first REAL*8 :: high_y, low_y, x, y OPEN (UNIT = free_unit, FILE = dig_pathfile, STATUS = 'OLD', & & PAD = 'YES', IOSTAT = ios) IF (ios /= 0) CALL Could_Not_Find_File(dig_pathfile) IF (ios /= 0) THEN WRITE (*,"(' ERROR in Dig_Type: Following file cannot be opened:' & & /' ',A)") TRIM(dig_pathfile) CALL DTraceback END IF any_titles = .FALSE. ! unless changed below... WRITE (*,"(/' Scanning through .dig file...')") first = .TRUE. scanning: DO READ (free_unit, "(A)", IOSTAT = ios) c2 IF (ios == -1) EXIT scanning ! EOF IF ((c2 == ' +').OR.(c2 == ' -')) THEN BACKSPACE (free_unit) READ (free_unit, *, IOSTAT = ios) x, y IF (ios == 0) THEN IF (first) THEN first = .FALSE. high_y = y low_y = y END IF high_y = MAX(high_y, y) low_y = MIN(low_y, y) END IF ! read was successful ELSE ! not a number line; either *** or a title BACKSPACE (free_unit) READ (free_unit, "(A)") line any_titles = any_titles .OR. (line(1:3) /= '***') END IF ! line has two numbers, or not END DO scanning CLOSE (UNIT = free_unit, IOSTAT = ios) WRITE (*,"('+Scanning through .dig file...DONE')") dig_is_lonlat = (low_y > -91.0D0).AND.(high_y < 91.0D0) END SUBROUTINE Dig_Type SUBROUTINE Fault_Key_Bottom (bottomlegend_gap_points, bottomlegend_used_points, tick_points) ! Plots 6 sample fault traces in bottom legend, ! in color (if ai_using_color) or b/w, with tick size tick_points; ! then adds to bottomlegend_used_points to record size of block. ! Note that these 6 samples are not grouped; you may wish to do ! this externally. IMPLICIT NONE REAL*8, INTENT(INOUT) :: bottomlegend_used_points REAL*8, INTENT(IN) :: bottomlegend_gap_points, tick_points INTEGER :: fp = 12 ! font points REAL*8 :: x1_points, x2_points, xc, xl, xr, & & y1, y1_points, y2, y2_points, y3 IF (.NOT.ai_bottomlegend_reserved) THEN WRITE (*,"(' Error: Fault_Key_Bottom requires ai_bottomlegend_reserved = T')") CALL DTraceback END IF CALL DReport_BottomLegend_Frame (x1_points, x2_points, y1_points, y2_points) !center of first block: xc = x1_points + bottomlegend_used_points + bottomlegend_gap_points + 31.0D0 !left and right limits of fault traces xl = xc - 22.0D0 xr = xc + 22.0D0 !baselines of 3 lines of text/fault trace: y3 = y1_points + MAX(2.0D0, tick_points) y2 = y3 + MAX(1.0D0*fp, 0.4D0*fp + tick_points) y1 = y2 + fp CALL DSet_Stroke_Color ('foreground') ! will remain, if .NOT. ai_using_color ! low-angle thrust (S or P or T with no dip_degrees; or, with dip_degrees<=45): dark_blue_, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DL12_Text (level = 1, x_points = xc, y_points = y1, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'Low-angle') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'thrust:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('dark_blue_') CALL DSet_Fill_or_Pattern (.FALSE., 'dark_blue_') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'P', size_points = tick_points, offset_points = 1.0D0) ! high-angle reverse fault (T or P with dip_degrees>45): mid_blue__, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') xc = xc + 62.0D0 xl = xc - 22.0D0 xr = xc + 22.0D0 CALL DL12_Text (level = 1, x_points = xc, y_points = y1, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'High-angle') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'reverse:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('mid_blue__') CALL DSet_Fill_or_Pattern (.FALSE., 'mid_blue__') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'T', size_points = tick_points, offset_points = 1.0D0) ! dextral (R): green_____, xc = xc + 62.0D0 xl = xc - 22.0D0 xr = xc + 22.0D0 CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'Dextral:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('green_____') CALL DSet_Fill_or_Pattern (.FALSE., 'green_____') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'R', size_points = tick_points, offset_points = 1.0D0) ! sinistral (L): brown_____, xc = xc + 62.0D0 xl = xc - 22.0D0 xr = xc + 22.0D0 CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'Sinistral:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('brown_____') CALL DSet_Fill_or_Pattern (.FALSE., 'brown_____') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'L', size_points = tick_points, offset_points = 1.0D0) ! high-angle normal (N or D with no dip_degrees; or, with dip_degrees>45): bronze____, xc = xc + 62.0D0 xl = xc - 22.0D0 xr = xc + 22.0D0 CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DL12_Text (level = 1, x_points = xc, y_points = y1, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'High-angle') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'normal:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('bronze____') CALL DSet_Fill_or_Pattern (.FALSE., 'bronze____') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'N', size_points = tick_points, offset_points = 1.0D0) ! low-angle detachment (N or D with dip_degrees<=45): red_______, xc = xc + 62.0D0 xl = xc - 22.0D0 xr = xc + 22.0D0 CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') CALL DL12_Text (level = 1, x_points = xc, y_points = y1, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'Low-angle') CALL DL12_Text (level = 1, x_points = xc, y_points = y2, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.0D0, & & text = 'detachment:') CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('red_______') CALL DSet_Fill_or_Pattern (.FALSE., 'red_______') END IF CALL DNew_L12_Path (1, xl, y3) CALL DLine_to_L12 (xr, y3) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = y3, dip_angle_radians = Pi/2.D0, & & style_byte = 'D', size_points = tick_points, offset_points = 1.0D0) ! wrap-up: adjust bottomlegend_used_points bottomlegend_used_points = bottomlegend_used_points + 365.0D0 END SUBROUTINE Fault_Key_Bottom SUBROUTINE Fault_Key_Right (rightlegend_gap_points, rightlegend_used_points, tick_points) ! Plots 6 sample fault traces in right legend, ! in color (if ai_using_color) or b/w, with tick size tick_points; ! then adds to rightlegend_used_points to record size of block. ! Note that these 6 samples are not grouped; you may wish to do ! this externally. IMPLICIT NONE REAL*8, INTENT(INOUT) :: rightlegend_used_points REAL*8, INTENT(IN) :: rightlegend_gap_points, tick_points INTEGER :: fp = 12 ! font points REAL*8 :: x1_points, x2_points, xc, xl, xr, y1_points, y2_points, yc IF (.NOT.ai_rightlegend_reserved) THEN WRITE (*,"(' Error: Fault_Key_Right requires ai_rightlegend_reserved = T')") CALL DTraceback END IF CALL DReport_RightLegend_Frame (x1_points, x2_points, y1_points, y2_points) yc = y2_points - rightlegend_used_points - rightlegend_gap_points xc = (x1_points + x2_points) / 2.0D0 xl = 0.8D0*x1_points + 0.2D0*x2_points xr = 0.2D0*x1_points + 0.8D0*x2_points CALL DSet_Stroke_Color ('foreground') ! will remain, if .NOT. ai_using_color ! low-angle thrust (S or P or T with no dip_degrees; or, with dip_degrees<=45): dark_blue_, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - 0.4D0*fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'Low-angle') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'thrust:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('dark_blue_') CALL DSet_Fill_or_Pattern (.FALSE., 'dark_blue_') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'P', size_points = tick_points, offset_points = 1.0D0) ! high-angle reverse (T or P and dip_degrees>45): mid_blue__, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'High-angle') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'reverse:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('mid_blue__') CALL DSet_Fill_or_Pattern (.FALSE., 'mid_blue__') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'T', size_points = tick_points, offset_points = 1.0D0) ! dextral (R): green_____, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'Dextral:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('green_____') CALL DSet_Fill_or_Pattern (.FALSE., 'green_____') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'R', size_points = tick_points, offset_points = 1.0D0) yc = yc - tick_points ! sinistral (L): brown_____, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'Sinistral:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('brown_____') CALL DSet_Fill_or_Pattern (.FALSE., 'brown_____') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'L', size_points = tick_points, offset_points = 1.0D0) yc = yc - tick_points ! high-angle normal (N or D with no dip_degrees; or, with dip_degrees>45): bronze____, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'High-angle') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'normal:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('bronze____') CALL DSet_Fill_or_Pattern (.FALSE., 'bronze____') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'N', size_points = tick_points, offset_points = 1.0D0) ! low-angle detachment (D or N with dip_degrees<=45): red_______, CALL DSet_Fill_or_Pattern (.FALSE., 'foreground') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'Low-angle') yc = yc - fp CALL DL12_Text (level = 1, x_points = xc, y_points = yc, & & angle_radians = 0.D0, font_points = fp, & & lr_fraction = 0.5D0, ud_fraction = 0.4D0, & & text = 'detachment:') yc = yc - MAX(1.0D0*fp, 0.4D0*fp + tick_points) CALL DSet_Line_Style (width_points = 2.0D0, dashed = .FALSE.) IF (ai_using_color) THEN CALL DSet_Stroke_Color ('red_______') CALL DSet_Fill_or_Pattern (.FALSE., 'red_______') END IF CALL DNew_L12_Path (1, xl, yc) CALL DLine_to_L12 (xr, yc) CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DDipTick_in_Plane (level = 1, x = xc, y = yc, dip_angle_radians = Pi/2.D0, & & style_byte = 'D', size_points = tick_points, offset_points = 1.0D0) ! wrap-up: adjust rightlegend_used_points rightlegend_used_points = y2_points - yc + 0.5D0 END SUBROUTINE Fault_Key_Right SUBROUTINE Fault_Traces (trace_choice, & ! NOTE: following parameters only needed for trace_choice >0: & width_array, dw_scale_amount, dw_scale_points, sense) ! optional parameters ! Draws fault traces on an open page, in 3 groups: ! (1) dip ticks (allowing for deletion; changes of color ! and/or pen width are not practical due to the mixture ! of stroked-but-unfilled with filled-but-unstroked kinds; ! (2) traces (allowing for width change if trace_choice = 1, ! or shadowing in either case); ! (3) text labels. ! Different variants are plotted according to "trace_choice" = ! 0 :: plot all traces with equal width, annotate with id number. ! In this case, slip sense is from byte(s) in f.dig. ! 1 :: plot traces with width proportional to "width_array"; ! scaled using parameters "dw_scale_amount" and "dw_scale_points"; ! annotate with value (rounded to two significant digits). ! In this case, slip sense is from array "sense". ! The following variables must be pre-declared, and values ! defined, in the calling program: ! CHARACTER*(*) :: traces_pathfile ! [path\]name of F.DIG file. ! REAL*8 :: tick_points ! desired size of dip ticks, in points ! LOGICAL :: ai_using_color ! from Adobe_Illustrator module data. ! REAL*8 :: mp_radius_meters, mp_scale_denominator ! from Map_Projections module data. IMPLICIT NONE INTEGER, INTENT(IN) :: trace_choice REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: width_array REAL*8, INTENT(IN), OPTIONAL :: dw_scale_points, dw_scale_amount CHARACTER*1, DIMENSION(:), INTENT(IN), OPTIONAL :: sense CHARACTER*1 :: another_byte, dip_byte, f_byte, one_byte, second_dip_byte, tick_byte CHARACTER*10 :: color_name, label CHARACTER*132 :: line INTEGER :: beyond_loc, count, end_loc, fault_number, & & i, i1, i2, internal_ios, ios, longest, nsteps, start_loc, which_byte LOGICAL :: got_dip_degrees, plotting REAL*8 :: angle, dip_azimuth_radians, dip_degrees, goal_radians, & & lat, lon, lr_fraction, nudge_radians, offset_points, & & r1, r2, step_radians, & & trace_points, trace_radians, ud_fraction, width, & & x1_meters, x2_meters, x1_points, x2_points, & & y1_meters, y2_meters, y1_points, y2_points REAL*8, DIMENSION(3) :: tvec, uvec1, uvec2, uvec3 REAL*8, DIMENSION(:,:), ALLOCATABLE :: one_trace !Read through and determine length of longest trace OPEN(UNIT = 21, FILE = traces_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) CALL Could_Not_Find_File(traces_pathfile) longest = 0 count = 0 DO READ (21, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, *, IOSTAT = internal_ios) lon, lat IF (internal_ios == 0) THEN ! read lon, lat OK count = count + 1 longest = MAX(longest,count) ELSE ! encountered ***END or title line or additional header count = 0 END IF END DO CLOSE (21) ALLOCATE ( one_trace(2,longest) ) ! First group will have dip ticks only. OPEN(UNIT = 21, FILE = traces_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) CALL Could_Not_Find_File(traces_pathfile) CALL DBegin_Group ! group of all dip ticks count = 0 DO ! reading lines of traces file READ (21, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A1,I4,A1,A1)", IOSTAT = internal_ios) f_byte, i, one_byte, another_byte IF ((internal_ios == 0).AND.((f_byte == 'F').OR.(f_byte == 'f'))) THEN fault_number = i IF (trace_choice == 0) THEN dip_byte = one_byte second_dip_byte = another_byte ELSE ! trace_choice /= 0; probably 1? dip_byte = sense(i) ! primary sense may have reversed! second_dip_byte = ' ' ! don't try to display secondary sense END IF SELECT CASE (trace_choice) CASE (0) plotting = .TRUE. CASE (1) width = width_array(fault_number) * dw_scale_points / dw_scale_amount plotting = (width >= 0.1) ! to avoid orphan dip-ticks after trace dissapears in slide/print! END SELECT got_dip_degrees = .FALSE. ! until... ELSE ! line was not a title; this leaves 3 possibilities: dip_degrees, (lon, lat), or *** end !try to read "dip_degrees"... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A)", IOSTAT = internal_ios) line start_loc = INDEX(line, "dip_degrees") IF (start_loc > 0) THEN ! found "dip_degrees" in line beyond_loc = start_loc + 11 ! first byte which is not part of "dip_degrees" end_loc = LEN_TRIM(line) line = line(beyond_loc:end_loc) READ (line, *, IOSTAT = internal_ios) dip_degrees IF (internal_ios == 0) got_dip_degrees = .TRUE. ELSE ! "dip_degrees" not there, so either (lon, lat), or *** end !continue, on assumption that next line is probably a (lon, lat) pair... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, *, IOSTAT = internal_ios) lon, lat IF (internal_ios == 0) THEN ! read lon, lat OK count = count + 1 one_trace(1,count) = lon one_trace(2,count) = lat ELSE ! encountered "*** end of line segment ***"; ! finished reading a trace IF (plotting.AND.(count > 1)) THEN ! got a trace to plot!!! SELECT CASE (trace_choice) CASE (0) ! constant width offset_points = 0.8D0 ! 40% of 2.0 points CASE (1) ! width prop. to trace_mma(fault_number) offset_points = 0.5D0 * width_array(fault_number) * dw_scale_points / dw_scale_amount END SELECT DO which_byte = 1, 2 IF (which_byte == 2) THEN IF (dip_byte /= ' ') dip_byte = second_dip_byte !Note that we don't plot the "second dip byte" if the first was blank, !because then the "second dip byte" might actually be first character of fault name! END IF IF (ai_using_color) THEN IF ((dip_byte == 'L').OR.(dip_byte == 'l')) THEN color_name = 'brown_____' tick_byte = dip_byte ELSE IF ((dip_byte == 'R').OR.(dip_byte == 'r')) THEN color_name = 'green_____' tick_byte = dip_byte ELSE IF ((dip_byte == 'T').OR.(dip_byte == 't')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' tick_byte = 'P' ELSE color_name = 'mid_blue__' tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' tick_byte = 'P' END IF ELSE IF ((dip_byte == 'P').OR.(dip_byte == 'p')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' tick_byte = 'P' ELSE color_name = 'mid_blue__' tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' tick_byte = 'P' END IF ELSE IF ((dip_byte == 'S').OR.(dip_byte == 's')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' tick_byte = 'P' ELSE color_name = 'mid_blue__' tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' tick_byte = 'P' END IF ELSE IF ((dip_byte == 'D').OR.(dip_byte == 'd')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'red_______' tick_byte = 'D' ELSE color_name = 'bronze____' tick_byte = 'N' END IF ELSE ! assume steep dip as the default: color_name = 'bronze____' tick_byte = 'N' END IF ELSE IF ((dip_byte == 'N').OR.(dip_byte == 'n')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'red_______' tick_byte = 'D' ELSE color_name = 'bronze____' tick_byte = 'N' END IF ELSE ! assume steep dip as the default: color_name = 'bronze____' tick_byte = 'N' END IF ELSE IF (dip_byte == ' ') THEN ! provide for case of blank second dip byte color_name = 'background' tick_byte = ' ' ELSE ! any other unexpected code color_name = 'foreground' tick_byte = ' ' END IF CALL DSet_Stroke_Color (color_name) ELSE ! .NOT.ai_using_color IF (dip_byte == ' ') THEN color_name = 'background' tick_byte = ' ' ELSE ! any non-blank byte color_name = 'foreground' IF ((dip_byte == 'L').OR.(dip_byte == 'l')) THEN tick_byte = dip_byte ELSE IF ((dip_byte == 'R').OR.(dip_byte == 'r')) THEN tick_byte = dip_byte ELSE IF ((dip_byte == 'T').OR.(dip_byte == 't')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN tick_byte = 'P' ELSE tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: tick_byte = 'P' END IF ELSE IF ((dip_byte == 'P').OR.(dip_byte == 'p')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN tick_byte = 'P' ELSE tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: tick_byte = 'P' END IF ELSE IF ((dip_byte == 'S').OR.(dip_byte == 's')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN tick_byte = 'P' ELSE tick_byte = 'T' END IF ELSE ! assume gentle dip as the default: tick_byte = 'P' END IF ELSE IF ((dip_byte == 'D').OR.(dip_byte == 'd')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN tick_byte = 'D' ELSE tick_byte = 'N' END IF ELSE ! assume steep dip as the default: tick_byte = 'N' END IF ELSE IF ((dip_byte == 'N').OR.(dip_byte == 'n')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN tick_byte = 'D' ELSE tick_byte = 'N' END IF ELSE ! assume steep dip as the default: tick_byte = 'N' END IF ELSE IF (dip_byte == ' ') THEN ! provide for case of blank second dip byte tick_byte = ' ' ELSE ! any other unexpected code tick_byte = ' ' END IF END IF CALL DSet_Stroke_Color (color_name) END IF ! ai_using_color or not trace_radians = 0.0D0 CALL DLonLat_2_Uvec(one_trace(1,1),one_trace(2,1),uvec1) DO i = 2, count CALL DLonLat_2_Uvec(one_trace(1,i),one_trace(2,i),uvec2) tvec = uvec2 - uvec1 trace_radians = trace_radians + DLength(tvec) uvec1 = uvec2 ! preparing to loop END DO trace_points = trace_radians * (mp_radius_meters/mp_scale_denominator) * 39.37D0 * 72.D0 IF ((trace_points >= (2.D0 * tick_points)).AND.(tick_points > 0.0D0)) THEN ! long enough trace to plot a dip tick on CALL DSet_Line_Style (width_points = 0.75D0, dashed = .FALSE.) CALL DSet_Fill_or_Pattern (use_pattern = .FALSE., color_name = color_name) nsteps = MAX(2,NINT(trace_points/(4.D0*tick_points))) step_radians = 1.001D0 * trace_radians / nsteps CALL DLonLat_2_Uvec(one_trace(1,1),one_trace(2,1),uvec1) trace_radians = 0.D0 goal_radians = step_radians DO i = 2, count CALL DLonLat_2_Uvec(one_trace(1,i),one_trace(2,i),uvec2) tvec = uvec2 - uvec1 nudge_radians = DLength(tvec) trace_radians = trace_radians + nudge_radians IF (trace_radians > goal_radians) THEN IF (which_byte == 1) THEN ! primary dip tick !place primary dip tick mid-segment to reduce chance of misalignment tvec = 0.5D0*(uvec1 + uvec2) ELSE ! secondary dip tick !place secondary dip tick near end of segment to reduce chance of collision with primary: tvec = 0.15D0 * uvec1 + 0.85D0 * uvec2 END IF CALL DMake_Uvec(tvec, uvec3) IF ((uvec1(1) /= uvec2(1)).OR.(uvec1(2) /= uvec2(2)).OR.(uvec1(3) /= uvec2(3))) THEN dip_azimuth_radians = DRelative_Compass(uvec1,uvec2) - Pi_over_2 CALL DDipTick_on_Sphere (uvec = uvec3, & & dip_azimuth_radians = dip_azimuth_radians, & & style_byte = tick_byte, & & size_points = tick_points, & & offset_points = offset_points) END IF goal_radians = goal_radians + step_radians END IF ! plotted a tick uvec1 = uvec2 ! preparing to loop END DO END IF ! long enough trace to plot on END DO ! which_byte = 1, 2 END IF ! got a trace to plot!!! count = 0 END IF ! hit end of segment at *** end END IF ! line did not contain "dip_degrees", so either (lon, lat) or *** end END IF ! line was not a title, so either dip_degrees, (lon, lat), or *** end END DO ! reading line of traces file (for 1st time, out of 3!) CALL DEnd_Group ! of all dip ticks CLOSE(21) ! Second group will have fault traces only. OPEN(UNIT = 21, FILE = traces_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) CALL Could_Not_Find_File(traces_pathfile) CALL DBegin_Group ! group of all fault traces count = 0 DO ! reading lines of traces file READ (21, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A1,I4,A1)", IOSTAT = internal_ios) f_byte, i, one_byte IF ((internal_ios == 0).AND.((f_byte == 'F').OR.(f_byte == 'f'))) THEN fault_number = i IF (trace_choice == 0) THEN dip_byte = one_byte ELSE dip_byte = sense(i) END IF SELECT CASE (trace_choice) CASE (0) plotting = .TRUE. CASE (1) plotting = (width_array(fault_number) > 0.0D0) END SELECT got_dip_degrees = .FALSE. ! until... ELSE ! line was not a title; this leaves 3 possibilities: dip_degrees, (lon, lat), or *** end !try to read "dip_degrees"... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A)", IOSTAT = internal_ios) line start_loc = INDEX(line, "dip_degrees") IF (start_loc > 0) THEN ! found "dip_degrees" in line beyond_loc = start_loc + 11 ! first byte which is not part of "dip_degrees" end_loc = LEN_TRIM(line) line = line(beyond_loc:end_loc) READ (line, *, IOSTAT = internal_ios) dip_degrees IF (internal_ios == 0) got_dip_degrees = .TRUE. ELSE ! "dip_degrees" not there, so either (lon, lat), or *** end !continue, on assumption that next line is probably a (lon, lat) pair... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, *, IOSTAT = internal_ios) lon, lat IF (internal_ios == 0) THEN ! read lon, lat OK count = count + 1 one_trace(1,count) = lon one_trace(2,count) = lat ELSE ! encountered *** end of segment ***; ! finished reading a trace IF (plotting.AND.(count > 1)) THEN ! got a trace to plot!!! SELECT CASE (trace_choice) CASE (0) ! constant width width = 2.0D0 ! points CASE (1) ! width prop. to trace_mma(fault_number) width = width_array(fault_number) * dw_scale_points / dw_scale_amount END SELECT IF (ai_using_color) THEN IF ((dip_byte == 'L').OR.(dip_byte == 'l')) THEN color_name = 'brown_____' ELSE IF ((dip_byte == 'R').OR.(dip_byte == 'r')) THEN color_name = 'green_____' ELSE IF ((dip_byte == 'T').OR.(dip_byte == 't')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' ELSE color_name = 'mid_blue__' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' END IF ELSE IF ((dip_byte == 'P').OR.(dip_byte == 'p')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' ELSE color_name = 'mid_blue__' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' END IF ELSE IF ((dip_byte == 'S').OR.(dip_byte == 's')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'dark_blue_' ELSE color_name = 'mid_blue__' END IF ELSE ! assume gentle dip as the default: color_name = 'dark_blue_' END IF ELSE IF ((dip_byte == 'D').OR.(dip_byte == 'd')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'red_______' ELSE color_name = 'bronze____' END IF ELSE ! assume steep dip as the default: color_name = 'bronze____' END IF ELSE IF ((dip_byte == 'N').OR.(dip_byte == 'n')) THEN IF (got_dip_degrees) THEN IF (dip_degrees <= 45.0D0) THEN color_name = 'red_______' ELSE color_name = 'bronze____' END IF ELSE ! assume steep dip as the default: color_name = 'bronze____' END IF ELSE ! any unexpected code color_name = 'foreground' END IF CALL DSet_Stroke_Color (color_name) ELSE CALL DSet_Stroke_Color (color_name = 'foreground') END IF ! ai_using_color or not CALL DSet_Line_Style (width_points = width, dashed = .FALSE.) CALL DNew_L67_Path (7,one_trace(1,1),one_trace(2,1)) DO i = 2, count CALL DGreat_to_L67(one_trace(1,i),one_trace(2,i)) END DO CALL DEnd_L67_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) END IF ! got a trace to plot!!! count = 0 END IF ! hit end of segment at *** end END IF ! line did not contain "dip_degrees", so either (lon, lat) or *** end END IF ! line was not a title, so either dip_degrees, (lon, lat), or *** end END DO ! reading line of traces file (for 2nd time) CALL DEnd_Group ! of fault traces CLOSE(21) ! Third group will have text labels. OPEN(UNIT = 21, FILE = traces_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) CALL Could_Not_Find_File(traces_pathfile) CALL DSet_Fill_or_Pattern (use_pattern = .FALSE., color_name = 'foreground') CALL DBegin_Group ! group of text labels count = 0 DO ! reading lines of traces file READ (21, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A1,I4,A1)", IOSTAT = internal_ios) f_byte, i, one_byte IF ((internal_ios == 0).AND.((f_byte == 'F').OR.(f_byte == 'f'))) THEN fault_number = i IF (trace_choice == 0) THEN dip_byte = one_byte ELSE dip_byte = sense(i) END IF SELECT CASE (trace_choice) CASE (0) plotting = .TRUE. CASE (1) width = width_array(fault_number) * dw_scale_points / dw_scale_amount plotting = (width >= 0.1D0) !to avoid orphan text after the trace dissapears in slide/print! END SELECT got_dip_degrees = .FALSE. ! until... ELSE ! line was not a title; this leaves 3 possibilities: dip_degrees, (lon, lat), or *** end !try to read "dip_degrees"... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, "(A)", IOSTAT = internal_ios) line start_loc = INDEX(line, "dip_degrees") IF (start_loc > 0) THEN ! found "dip_degrees" in line beyond_loc = start_loc + 11 ! first byte which is not part of "dip_degrees" end_loc = LEN_TRIM(line) line = line(beyond_loc:end_loc) READ (line, *, IOSTAT = internal_ios) dip_degrees IF (internal_ios == 0) got_dip_degrees = .TRUE. ELSE ! "dip_degrees" not there, so either (lon, lat), or *** end !continue, on assumption that next line is probably a (lon, lat) pair... BACKSPACE (21) ! instead of READ(line,*,...) which is BUGGY! READ (21, *, IOSTAT = internal_ios) lon, lat IF (internal_ios == 0) THEN ! read lon, lat OK count = count + 1 one_trace(1,count) = lon one_trace(2,count) = lat ELSE ! encountered ***END OF SEGMENT***; ! finished reading a trace IF (plotting.AND.(count > 1)) THEN ! got a trace to plot!!! CALL DLonLat_2_Uvec (one_trace(1,1),one_trace(2,1),uvec1) CALL DLonLat_2_Uvec (one_trace(1,count),one_trace(2,count),uvec2) tvec = uvec2 - uvec1 IF (DLength(tvec) > 0.0D0) THEN SELECT CASE (trace_choice) CASE (0) ! label with id number WRITE (label,"(I10)") fault_number label = ADJUSTL(label) r1 = one_trace(1,1) r2 = one_trace(2,1) angle = Pi_over_2 - DRelative_Compass (uvec1, uvec2) lr_fraction = 0.0D0 ud_fraction = 1.0D0 CALL DL67_Text (level = 7, & & r1 = r1, r2 = r2, & & angle_radians = angle, from_east = .TRUE., & & font_points = 8, & & lr_fraction = lr_fraction, ud_fraction = ud_fraction, & & text = TRIM(label)) CASE (1) ! label with mm/a label = ADJUSTL(DASCII8(width_array(fault_number))) offset_points = 0.5D0 * width_array(fault_number) * dw_scale_points / dw_scale_amount !because of the need for this offset, use level 2 text. i1 = DInt_Below(count/2.0D0) i2 = i1 + 1 CALL DLonLat_2_Uvec (one_trace(1,i1),one_trace(2,i1),uvec1) CALL DLonLat_2_Uvec (one_trace(1,i2),one_trace(2,i2),uvec2) CALL DProject (uvec = uvec1, x = x1_meters, y = y1_meters) CALL DProject (uvec = uvec2, x = x2_meters, y = y2_meters) CALL DMeters_2_Points (x1_meters,y1_meters, x1_points,y1_points) CALL DMeters_2_Points (x2_meters,y2_meters, x2_points,y2_points) ! angle is of text baseline, counterclockwise from right, radians angle = DAtan2f((y1_points - y2_points),(x1_points - x2_points)) r1 = 0.5D0 * (x1_points + x2_points) - offset_points * DSIN(angle) r2 = 0.5D0 * (y1_points + y2_points) + offset_points * DCOS(angle) lr_fraction = 0.5D0 ud_fraction = -0.2D0 CALL DL12_Text (level = 2, & & x_points = r1, y_points = r2, & & angle_radians = angle, & & font_points = 8, & & lr_fraction = lr_fraction, ud_fraction = ud_fraction, & & text = TRIM(label)) END SELECT END IF ! trace has different 1st, last points END IF ! got a trace to plot!!! count = 0 END IF ! hit end of segment at *** end END IF ! line did not contain "dip_degrees", so either (lon, lat) or *** end END IF ! line was not a title, so either dip_degrees, (lon, lat), or *** end END DO ! reading line of traces file (for 3rd time) CALL DEnd_Group ! of text labels CLOSE(21) DEALLOCATE ( one_trace ) END SUBROUTINE Fault_Traces ! SUBROUTINE File_List( file_type, & ! & suggested_file, & ! & using_path ) ! ! 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, using_path ! CHARACTER*1 :: first_letter ! CHARACTER*70 :: line = ' ', old_name ! CHARACTER*80 :: 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 ! NOTE: Quoted here for easy reference; ! ! INTEGER(4) length ! however, do NOT remove the leading '!' comment-out characters, ! ! INTEGER(4) permit ! or you will get a cryptic compiler error! ! ! 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 = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.*' ! ! ELSE IF (file_type == "*.dig") THEN ! WRITE (*,"(/' The following appear to be basemap (.dig) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.DIG' ! ! ELSE IF (file_type == "*.cat") THEN ! WRITE (*,"(/' The following appear to be seismic catalog (.cat) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.CAT' ! ! ELSE IF (file_type == "*.feg") THEN ! WRITE (*,"(/' The following appear to be FE grid (.feg) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.FEG' ! ELSE IF (file_type == ".grd") THEN ! WRITE (*,"(/' The following appear to be gridded data (.grd) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.GRD' ! ELSE IF (file_type == "p*.nki") THEN ! WRITE (*,"(/' The following appear to be Parameter input (p*.nki) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.NKI' ! ELSE IF (file_type == "s*.nki") THEN ! WRITE (*,"(/' The following appear to be Stress-direction input (s*.nki) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.NKI' ! ELSE IF (file_type == "s*.nko") THEN ! WRITE (*,"(/' The following appear to be interpolated Stress-direction output (s*.nko) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.NKO' ! ELSE IF (file_type == "v*.out") THEN ! WRITE (*,"(/' The following appear to be velocity output (v*.out) files:')") ! files = TRIM(using_path) // & ! defined in NeoKineMap above ! & '*.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 ',A,';')") TRIM(using_path) ! CALL DPrompt_for_String('Select new directory (for this file only)?',using_path,using_path) ! GO TO 10 ! ELSE ! count > 0, but line empty ! GO TO 100 ! END IF ! END IF ! first_letter = info.name(1:1) ! !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 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 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 CHARACTER(80) FUNCTION Get_filename (unit) ! 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". CHARACTER(132) :: buffer INTEGER :: i, ios LOGICAL :: past INTEGER, INTENT (IN) :: unit ! Fortran device number READ (unit,"(A)", IOSTAT = ios) buffer IF (ios /= 0) CALL Bad_Parameters() 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' Get_filename = buffer(1:80) END FUNCTION Get_filename SUBROUTINE Histogram (real_list, list_length, skip_zeros, maximum, minimum) ! Puts a printer-plot on the default device, no more than ! (n15 + 2) rows tall by n70 bytes wide, showing the range and ! distribution of values within real_list. IMPLICIT NONE REAL*8, DIMENSION(:), INTENT(IN) :: real_list INTEGER, INTENT(IN) :: list_length LOGICAL, INTENT(IN) :: skip_zeros REAL*8, INTENT(OUT) :: maximum, minimum INTEGER, PARAMETER :: n15 = 15, n70 = 70 REAL*8, PARAMETER :: Huge = 9.99D59 CHARACTER*10 :: number10 CHARACTER*(n70) :: line INTEGER :: highest, i, j, length INTEGER, DIMENSION(:), ALLOCATABLE :: counters REAL*8 :: dx, factor IF (list_length < 1) RETURN IF (skip_zeros) THEN maximum = -Huge minimum = +Huge DO i = 1, list_length IF (real_list(i) /= 0.0D0) THEN maximum = MAX(maximum, real_list(i)) minimum = MIN(minimum, real_list(i)) END IF END DO ELSE maximum = real_list(1) minimum = maximum DO i = 2, list_length maximum = MAX(maximum, real_list(i)) minimum = MIN(minimum, real_list(i)) END DO END IF dx = (maximum - minimum) / (n70 - 1) IF (dx == 0.0D0) dx = 1.00D0 ! avoid divide-by-zero ALLOCATE ( counters(n70) ) counters = 0 ! whole array DO i = 1, list_length IF (skip_zeros) THEN IF (real_list(i) /= 0.0D0) THEN j = 1 + DInt_Below((real_list(i) - minimum) / dx) counters(j) = counters(j) + 1 END IF ELSE j = 1 + DInt_Below((real_list(i) - minimum) / dx) counters(j) = counters(j) + 1 END IF END DO highest = 0 DO j = 1, n70 highest = MAX(highest,counters(j)) END DO factor = (1.D0 * n15) / (1.D0 * highest) DO i = n15, 1, -1 ! rows, from top to bottom line = ' ' DO j = 1, n70 ! columns !In bottom row only, put "." to show non-zero contents: IF (i == 1) THEN ! bottom row IF (counters(j) > 0) line(j:j) = '.' END IF IF (NINT(factor * counters(j)) >= i) line(j:j) = '*' END DO ! columns WRITE (*,"(' ',A)") TRIM(line) END DO ! rows line = REPEAT('-',n70) WRITE (*,"(' ', A)") line number10 = DASCII10(minimum) line = TRIM(ADJUSTL(number10)) number10 = DASCII10(maximum) number10 = ADJUSTL(number10) length = LEN_TRIM(number10) line((n70 - length + 1):n70) = TRIM(number10) WRITE (*,"(' ', A)") line DEALLOCATE ( counters ) END SUBROUTINE Histogram 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 in NeoKinema grids, because:'& &/' I. NeoKinema does not require fault elements.'& &/' 1. NeoKinema has logic to add the compliance of any number of'& &/' faults to the continuum (triangle) elements that contain them.'& &/' 2. NeoKinema 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 program NeoKineMap has logic to plot the heave-rates'& &/' of these faults, and also velocity fields with fault '& &/' disontinuities, without the use of fault elements.') 102 FORMAT (& &/' II. Fault elements cause bad grid topology.'& &/' 1. Fault elements are not read or stored by NeoKineMap.'& &/' 2. With fault elements deleted, 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 OrbWeaver, 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 OrbWeaver 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() !The purpose of this routine is to give the user time to read an error message !before the program stops (and the operating system closes its window). !The programmer should WRITE the error, CALL Pause(), and then STOP. IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM GPS_Postprocessor