PROGRAM GPS_Covariance ! Creates a GPS velocity covariance matrix (or reads an existing one), ! and adds any number of velocity-perturbation vectors, to represent ! any transient behaviour which is not part of the secular accumulation ! of tectonic strain. ! ! Also has an option to add a group of 3 velocity-perturbation vectors ! associated with uncertainty in the velocity reference frame. ! ! Input format is Peter Bird's .GPS format for the GPS data file ! (required) and Peter Bird's .GP2 format for the covariance file ! (optional input; always produced as output). ! Each is in the format documented in program NeoKinema ! (and its associated on-line tutorial pages). ! ! Data on active magma chambers can be read from a simple ASCII ! flat-file, with format specific to this program; ! see source code for examples. ! ! Data on opening/closing (mode-1) cracks can be read from a simple ASCII ! flat-file, with format specific to this program; ! see source code for examples. ! ! Data on fault-creep (shear dislocations) can be read from a simple ASCII ! flat-file, with format specific to this program; ! see source code for examples. ! ! Any other spatial pattern of transient or non-tectonic (noise) ! velocities can be input in .GPS format, using one file per noise ! source. The overall velocity amplitudes in such files are significant. ! ! If such local sources of covariance are added, then custom .GPS ! and .DIG files permit plotting of velocity perturbations and ! associated dislocation patches (if any), respectively. ! ! Versions 2+ (November 2017+) also produce an additional output file, ! GPS_Covariance_log.txt, which should be *SAVED* for eventual use ! with new utility program GPS_Postprocessor2. ! ! This is version 3 (September 2020) which adds mode-1 opening cracks. ! ! By Peter Bird, ! Dept. of Earth, Planetary, and Space Sciences ! UCLA ! 2014.10.19, 2015.01.03, 2015.02.12, 2016.03.02, 2017.11.06, 2020.09.27, ! with comments and printed texts further revised 2023.03.13. USE DSphere ! my code for spherical geometry, using unit vectors in ! a Cartesian 3-space, with DOUBLE PRECISION accuracy. ! From Peter Bird's file DSphere.f90. USE DDislocation ! my code to correct observed benchmark velocities by ! adding estimated long-term-average coseismic velocities. ! From Peter Bird's file DDislocation.f90. IMPLICIT NONE CHARACTER(3) :: file_index_c3 CHARACTER(5) :: zone ! time zone CHARACTER(8) :: date CHARACTER(10) :: clock_time ! wall clock time CHARACTER(15) :: velocity_perturbation_CLASS ! for HEADER in GPS_Covariance_log.txt CHARACTER(80) :: alleged_benchmark_name, c80, filename CHARACTER(134) :: GPS_pathfile, old_GP2_pathfile, new_GP2_pathfile, & & magma_chamber_pathfile, input_transient_GPS_file, afterslip_dislocation_pathfile, & & opening_crack_patches_pathfile CHARACTER(134) :: gps_title ! title (first) line of .gps file CHARACTER(134) :: gps_format ! format (second) line of .gps file CHARACTER(134) :: gps_header ! header (third) line of .gps file CHARACTER(134) :: crack_title ! title (first) line of opening/closing (mode-1) crack data file CHARACTER(134) :: crack_header ! format (second) line of opening/closing (mode-1) crack data file CHARACTER(134) :: crack_units ! header (third) line of opening/closing (mode-1) crack data file file CHARACTER(134) :: dislocation_title ! title (first) line of afterslipping-dislocations file CHARACTER(134) :: dislocation_header ! format (second) line of afterslipping-dislocations file CHARACTER(134) :: dislocation_units ! header (third) line of afterslipping-dislocations file CHARACTER(134) :: line, volcano_header CHARACTER(80), DIMENSION(:), ALLOCATABLE :: benchmark_name ! any site name or reference to source that follows the velocity reference frame in the .gps file; ! saved to be output into the g*.nko file. ! (1:external_benchmarks = external benchmark index); LATER: ! (1:internal_benchmarks = internal benchmark index) after compaction INTEGER :: count_affected, crack_count, & & external_benchmarks, extra_files, & & file_index, & & geodetic_NDOF, & & i, iE, iN, ios, irow, & & j, jcolumn, jE, jN, & & k, & & last_slash, lines_in_GP2, & & n, & & patch_count INTEGER, DIMENSION(8) :: datetimenumber ! for output from DATE_AND_TIME INTEGER, DIMENSION(:), ALLOCATABLE :: active_benchmarks LOGICAL :: afterslipping, extras, openingCracks, read_GP2, reframe, swelling REAL*8, PARAMETER :: Earth_radius_km = 6371.D0 REAL*8, PARAMETER :: s_per_year = 3.15576D7 REAL*8 :: argume, azimuth, azimuth_radians, btheta, bphi, correlation, & & del_V_East, del_V_North, depth_in_km, dip, dipf, duthet, duphi, & & Earth_radius_m, East_longitude, Elon, expansion, & & ftheta, fphi, frame_sigma_mmpa, & & highest_velocity_mmpa, & & influence_radius_km, influence_radius_radians, & & lat, length, lon, loosening_degpMa, loosening_radiansPerYear, & & Nlat, North_latitude, opening, peak_mm_per_year, & & r_km, reduction_factor, & & sine, sine_influence_radius_radians, sinistral, & & vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, velocity_mmpa, vertical, & & wedge, zBot, zTop REAL*8, DIMENSION(3) :: Euler, phi, slip, test_vector, theta, uvec, vec1, vec2, volcano_uvec REAL*8, DIMENSION(:), ALLOCATABLE :: benchmark_lat, benchmark_lon, benchmark_vE_sigma, benchmark_vN_sigma, benchmark_correlation, looseness REAL*8, DIMENSION(:,:), ALLOCATABLE :: benchmark_motions, benchmark_uvec ! location of geodetic benchmark; ! Cartesian unit vector from Earth center: ! (1:3 = x,y,z; 1:external_benchmarks = external benchmark index); LATER: ! (1:3 = x,y,z; 1:internal_benchmarks = internal benchmark index) after compaction DOUBLE PRECISION :: azimuth_radians_dp, far_radians_dp, lat_dp, lon_dp, lf, RHS_azimuth_radians_dp, variance DOUBLE PRECISION, DIMENSION(3) :: center_uvec_dp, omega_uvec_dp, & & uvec1_dp, uvec2_dp, uvec3_dp, uvec4_dp, uvec5_dp, uvec6_dp ! for outlining each dislocation patch, & projected surface trace DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: benchmark_covariance ! (2, 2, i) DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: covariance_mmpa2 ! (geodetic_NDOF, geodetic_NDOF) ! geodetic covariance matrix for horizontal velocity components, ! in units of (mm/a)**2. ! Ordering of benchmarks must be the same as in .gps file! ! vEast is numbered before vNorth at each benchmark. ! Matrix is symmetric, so half of the information is redundant. !(1:geodetic_nDOF, 1:geodetic_nDOF) Earth_radius_m = 1000.0D0 * Earth_radius_km WRITE (*, "(' PROGRAM GPS_covariance:')") WRITE (*, "(' ')") WRITE (*, "(' Creates a GPS velocity covariance matrix (or reads an existing one),')") WRITE (*, "(' and adds any number of velocity-perturbation vectors, to represent')") WRITE (*, "(' any transient behaviour which is not part of the secular accumulation')") WRITE (*, "(' of tectonic strain.')") WRITE (*, "(' ')") WRITE (*, "(' Also has an option to add a group of 3 velocity-perturbations vectors')") WRITE (*, "(' associated with uncertainty in the velocity reference frame.')") WRITE (*, "(' ')") WRITE (*, "(' Input format is Peter Bird''s .GPS format for the GPS data file')") WRITE (*, "(' (required) and Peter Bird''s .GP2 format for the covariance file')") WRITE (*, "(' (optional input; always produced as output).')") WRITE (*, "(' Each is in the format documented in program NeoKinema.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' Data on active magma chambers can be read from a simple ASCII')") WRITE (*, "(' flat-file, with format specific to this program;')") WRITE (*, "(' see source code for examples.')") WRITE (*, "(' ')") WRITE (*, "(' Data on opening/closing (mode-1) cracks can be read from a simple ASCII')") WRITE (*, "(' flat-file, with format specific to this program;')") WRITE (*, "(' see source code for examples.')") WRITE (*, "(' ')") WRITE (*, "(' Data on fault-creep (shear dislocations) can be read from a simple ASCII')") WRITE (*, "(' flat-file, with format specific to this program;')") WRITE (*, "(' see source code for examples.')") WRITE (*, "(' ')") WRITE (*, "(' Any other spatial pattern of transient or non-tectonic (noise)')") WRITE (*, "(' velocities can be input in .GPS format, using one file per noise')") WRITE (*, "(' source. The overall velocity amplitudes in such files are significant.')") WRITE (*, "(' ')") WRITE (*, "(' If such local sources of covariance are added, then custom .GPS')") WRITE (*, "(' and .DIG files permit plotting of velocity-perturbations and')") WRITE (*, "(' associated dislocation patches (if any), respectively.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' Versions 2+ (November 2017+) also produce an additional output file,')") WRITE (*, "(' GPS_Covariance_log.txt, which should be *SAVED* for eventual use')") WRITE (*, "(' with new utility program GPS_Postprocessor2.')") WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird,')") WRITE (*, "(' Dept. of Earth, Planetary, and Space Sciences')") WRITE (*, "(' University of California Los Angeles')") WRITE (*, "(' 2014.10.19, 2015.01.03, 2015.02.12, 2015.03.26,')") WRITE (*, "(' 2016.03.03, 2017.11.06, 2020.09.27,')") WRITE (*, "(' with comments and printed texts further revised 2023.03.13.')") WRITE (*, "(' ')") CALL Pause() !=============================================================================================================== !Open the log-file, to record GPS_Covariance operations (for later use by GPS_Postprocessor2): OPEN (UNIT = 51, FILE = "GPS_Covariance_log.txt") ! Absolute OPEN; overwrites any existing log. CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) WRITE (51,"('GPS_Covariance_log.txt, produced on ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & datetimenumber(1), datetimenumber(2), datetimenumber(3), & datetimenumber(5), datetimenumber(6), datetimenumber(7) !Read existing .GPS file (required): ! Each .gps file will have 3 lines of headers: ! ------------------------------------------------------ ! File name and source(s) ! FORTRAN format for reading the data lines that follow the headers ! [ column header labels, in a standard order: ] ! E_lon_deg N_lat_deg v_E_mmpa v_N_mmpa v_E_sigma v_N_sigma correlation frame identifier(s) ! ------------------------------------------------------------------------------------------ ! with the (obvious) meanings: ! E_lon_deg = longitude, in degrees from Greenwich meridian, with East positive ! N_lat_deg = latitude, in degrees from equator, with North positive ! v_E_mmpa = velocity component to East, in millimeters per year ! v_N_mmpa = velocity component to North, in millimeters per year ! v_E_sigma = standard deviation (1-sigma) of v_E_mmpa, also in mm/a ! v_N_sigma = standard deviation (1-sigma) of v_N_mmpa, also in mm/a ! correlation = coefficient of correlation between v_E_mmpa and v_N_mmpa ! reference_frame = reference frame for velocity, left-justified, limited to 15 bytes ! identifier(s) = optional station name and/or source reference, if a compilation ! Following these headers there is one line of data per benchmark. 10 WRITE (*, *) CALL DPrompt_for_String('[Drive:][\path\]filename of input .GPS file:', "*.GPS", GPS_pathfile) OPEN (UNIT = 1, FILE = TRIM(GPS_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This .GPS file not found (in current folder.')") CALL Pause() GO TO 10 END IF WRITE (51, "('.GPS input file = ', A)") TRIM(GPS_pathfile) PRINT "(' ',4X,'Counting geodetic benchmarks in ',A)", TRIM(GPS_pathfile) ! count benchmarks READ (1, "(A)") gps_title READ (1, "(A)") gps_format READ (1, "(A)") gps_header external_benchmarks = 0 benchmarking: DO READ (1, gps_format, IOSTAT = ios) lon, lat IF (ios == 0) THEN external_benchmarks = external_benchmarks + 1 ELSE EXIT benchmarking END IF END DO benchmarking CLOSE (1) WRITE (*, "(' Found ',I8,' benchmarks in this file.')") external_benchmarks ALLOCATE ( benchmark_uvec(3, external_benchmarks) ) ALLOCATE ( benchmark_covariance(2, 2, external_benchmarks) ) ALLOCATE ( benchmark_name(external_benchmarks) ) ALLOCATE ( benchmark_lon(external_benchmarks) ) ALLOCATE ( benchmark_lat(external_benchmarks) ) ALLOCATE ( benchmark_vE_sigma(external_benchmarks) ) ALLOCATE ( benchmark_vN_sigma(external_benchmarks) ) ALLOCATE ( benchmark_correlation(external_benchmarks) ) !read file again, to memorize contents: OPEN(UNIT = 1, FILE = TRIM(GPS_pathfile), STATUS = "OLD", ACTION = "READ", PAD = "YES") READ (1, "(A)") gps_title READ (1, "(A)") gps_format READ (1, "(A)") gps_header !second time through, data is remembered: DO i = 1, external_benchmarks READ (1, gps_format, IOSTAT = ios) lon, lat, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation, c80, benchmark_name(i) benchmark_lon(i) = lon benchmark_lat(i) = lat benchmark_vE_sigma(i) = vE_sigma benchmark_vN_sigma(i) = vN_sigma benchmark_correlation(i) = correlation CALL DXyz_from_lonlat(lon, lat, uvec) benchmark_uvec(1:3, i) = uvec(1:3) IF (vE_sigma <= 0.0D0) THEN ! prevent non-positive entries on diagonal PRINT "(' ',8X,'ERROR: vE_sigma is non-positive at benchmark ',A)", TRIM(benchmark_name(i)) STOP END IF IF (vN_sigma <= 0.0D0) THEN ! prevent non-positive entries on diagonal PRINT "(' ',8X,'ERROR: vN_sigma is non-positive at benchmark ',A)", TRIM(benchmark_name(i)) STOP END IF benchmark_covariance(1, 1, i) = vE_sigma**2 ! Note: Keeping external conventions: mm/year, and East-before-North benchmark_covariance(2, 2, i) = vN_sigma**2 benchmark_covariance(1, 2, i) = correlation * vN_sigma * vE_sigma benchmark_covariance(2, 1, i) = benchmark_covariance(1, 2, i) ! symmetric END DO WRITE (*, "(' Finished memorizing benchmark positions and 2x2 covariances.')") CALL Pause() !=============================================================================================================== !Create covariance matrix, in case none is available (or in case it is incomplete): WRITE (*, *) WRITE (*, "(' Building simple block-diagonal covariance matrix...')") geodetic_NDOF = 2 * external_benchmarks ALLOCATE ( covariance_mmpa2(geodetic_NDOF, geodetic_NDOF) ) covariance_mmpa2 = 0.0D0 ! whole matrix; just initializing DO i = 1, external_benchmarks iN = 2 * i iE = iN - 1 covariance_mmpa2(iE, iE) = benchmark_covariance(1, 1, i) covariance_mmpa2(iE, iN) = benchmark_covariance(1, 2, i) covariance_mmpa2(iN, iE) = benchmark_covariance(2, 1, i) covariance_mmpa2(iN, iN) = benchmark_covariance(2, 2, i) END DO WRITE (*, "(' done.')") !=============================================================================================================== !check for existing .GP2 file that can be used as a basis for expansion... WRITE (*, *) WRITE (*, "(' Do you wish to read a MATCHED, CORRESPONDING .GP2 file that')") CALL DPrompt_for_Logical('corrects and/or expands this covariance?', .TRUE., read_GP2) ! *.gp2 Covariance matrix for velocity components listed in the ! preceding .gps file. Order of degrees of freedom is: ! East component before North component for each benchmark, ! and benchmarks in same order as in the .gps file. ! Note that the .gp2 file is symmetrical, so it is sufficient ! to enter the diagonal and upper (or lower) triangle. ! Zero values do not need to be entered. ! Values can be entered in any order. ! Each line in the .gp2 file has: ! irow jcolumn variance ! where: irow is an integer ! jcolumn is an integer ! variance is a real number in units of (mm/year)**2. ! Note that variance is required to be positive for all diagonal entries; ! zero is not acceptable on the diagonal. ! If "variance" does not agree with information from the .gps file, ! then information from the .gps file is over-written and replaced. IF (read_GP2) THEN 20 CALL DPrompt_for_String('[Drive:][\path\]filename of MATCHED, CORRESPONDING .GP2 file:', "*.GP2", old_GP2_pathfile) OPEN (UNIT = 2, FILE = TRIM(old_GP2_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: This .GP2 file not found (in current folder.')") CALL Pause() GO TO 20 END IF WRITE (51, "('.GP2 input file = ', A)") TRIM(old_GP2_pathfile) ! record this choice in the log-file lines_in_GP2 = 0 base_covariance: DO READ (2, *, IOSTAT = ios) irow, jcolumn, variance IF (ios /= 0) EXIT base_covariance lines_in_GP2 = lines_in_GP2 + 1 IF ((irow < 1).OR.(irow > geodetic_NDOF)) THEN WRITE (*, "(' ERROR: Illegal irow = ',I12,' in line ',I12,' of .GP2 file.')") irow, lines_in_GP2 CALL Pause() STOP END IF IF ((jcolumn < 1).OR.(jcolumn > geodetic_NDOF)) THEN WRITE (*, "(' ERROR: Illegal jcolumn = ',I12,' in line ',I12,' of .GP2 file.')") irow, lines_in_GP2 CALL Pause() STOP END IF covariance_mmpa2(irow, jcolumn) = variance covariance_mmpa2(jcolumn, irow) = variance ! necessary for symmetry (although redundant along the diagonal). END DO base_covariance WRITE (*, "(' ',I8,' lines were read from this .GP2 file.')") lines_in_GP2 CALL Pause() ELSE ! no input .GP2 file; however, this program will create one. WRITE (51, "('.GP2 input file = ','[none]')") ! record this choice in the log-file END IF WRITE (51, "('===============================================================================')") ! end of header !=============================================================================================================== 100 WRITE (*, *) CALL DPrompt_for_Logical('Do you wish to add covariance representing the uncertainty in the definition of the velocity reference frame?', .TRUE., reframe) IF (reframe) THEN frame_sigma_mmpa = 0.0D0 CALL DPrompt_for_Real('Enter standard deviation of all GPS velocities due to reference frame issues (mm/year):', frame_sigma_mmpa, frame_sigma_mmpa) IF (frame_sigma_mmpa <= 0.0D0) THEN WRITE (*, "(' Error: frame_sigma_mmpa must be positive. Suggestion: 0.3 to 1.0 mm/a (?)')") CALL Pause() GO TO 100 END IF loosening_radiansPerYear = frame_sigma_mmpa * 1.0D-6 / Earth_radius_km loosening_degpMa = degrees_per_radian * loosening_radiansPerYear * 1.0D6 ! == degrees_per_radian * frame_sigma_mmpa / Earth_radius_km, but THAT formula looks erroneous! PRINT "(' ','Adding reference-frame-loosening rotations of ',F12.4,' degree/Ma')", loosening_degpMa ALLOCATE ( looseness(geodetic_nDOF) ) velocity_perturbation_CLASS = "REFERENCE-FRAME" ! for HEADER in GPS_Covariance_log.txt DO k = 1, 3 ! 3 rotation axes Euler = 0.0D0 ! whole 3-vector is zeroed. Euler(k) = loosening_degpMa / (degrees_per_radian * 1.0D6 * s_per_year) ! just one component is now non-zero. !now we have Euler in radians per second !Create memo in log-file, for eventual use by GPS_Postprocessor WRITE (51, "(A15)") velocity_perturbation_CLASS WRITE (51, "('User entered frame_sigma_mmpa = ',F8.3,' mm/a.')") frame_sigma_mmpa ! source of information WRITE (51, "('Euler vector = ('ES12.4,',',ES12.4,',',ES12.4,') radians/s')") Euler(1), Euler(2), Euler(3) ! specific perturbation vector WRITE (51, "('East_lon N_lat dV_East dV_North SITE')") ! header for columns to follow(?) DO i = 1, external_benchmarks uvec(1:3) = benchmark_uvec(1:3, i) CALL DCross (Euler, uvec, vec1) ! vec1 is in radians/s vec1 = vec1 * Earth_radius_km * 1.0D3 ! now, vec1 is in m/s vec1 = vec1 * 1.0D3 * s_per_year ! now, vec1 is in mm/a CALL DLocal_Theta(uvec, vec2) looseness(2*i-1) = DDot(vec1, vec2) ! vTheta produced by rotation, in mm/a CALL DLocal_Phi (uvec, vec2) looseness(2*i) = DDot(vec1, vec2) ! vPhi produced by rotation, in mm/a !Add memo to the log-file: del_V_East = +looseness(2*i) ! in mm/a del_V_North = -looseness(2*i-1) IF ((ABS(del_V_East) >= 0.001D0).OR.(ABS(del_V_North) >= 0.001D0)) THEN WRITE (51, "(F8.3, ' ', F7.3, ' ', 2F10.3, ' ', A)") benchmark_lon(i), benchmark_lat(i), del_V_East, del_V_North, TRIM(benchmark_name(i)) END IF END DO DO i = 1, geodetic_nDOF DO j = 1, geodetic_nDOF covariance_mmpa2(i, j) = covariance_mmpa2(i, j) + (looseness(i) * looseness(j)) END DO END DO WRITE (51, "('-------------------------------------------------------------------------------')") ! division between vectors END DO ! 3 rotations DEALLOCATE ( looseness ) END IF ! reframe !=============================================================================================================== 110 WRITE (*, *) CALL DPrompt_for_Logical('Do you have a file of active magma chambers?', .TRUE., swelling) IF (swelling) THEN velocity_perturbation_CLASS = "VOLCANO " ! for HEADER in GPS_Covariance_log.txt file_index = 0 ! but will be incremented for each local-transient V00n .GPS file written WRITE (*, *) WRITE (*, "(' Required contents of the magma-chamber data table (with only 1 case):')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, "(' {title line; perhaps file-name, and comments?}')") WRITE (*, "(' East_longitude North_latitude depth_in_km peak_mm_per_year')") WRITE (*, "(' -75.366 41.012 6.0 2.5')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, *) CALL DPrompt_for_String('[Drive:][\path\]filename for magma-chamber data file:', ' ', magma_chamber_pathfile) OPEN (UNIT = 11, FILE = TRIM(magma_chamber_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Magma-chamber file not found (in this folder).')") CALL Pause() GO TO 110 END IF WRITE (*, *) READ (11, "(A)") volcano_header WRITE (*, "(' ', A)") TRIM(volcano_header) READ (11, "(A)") volcano_header WRITE (*, "(' ', A)") TRIM(volcano_header) chambers: DO READ (11, "(A)" , IOSTAT = ios) line IF (ios /= 0) EXIT chambers WRITE (*, "(' ', A)") TRIM(line(1:78)) READ (line, *, IOSTAT = ios) East_longitude, North_latitude, depth_in_km, peak_mm_per_year ! , name? !header in log-file WRITE (51, "(A15)") velocity_perturbation_CLASS WRITE (51, "('magma chamber file = ', A)") TRIM(magma_chamber_pathfile) ! source of information WRITE (51, "('East_longitude = ',F8.3,', North_latitude = ',F7.3,', depth_in_km = 'F6.2,', peak_mm_per_year = ',F6.2)") East_longitude, North_latitude, depth_in_km, peak_mm_per_year ! specific perturbation WRITE (51, "('East_lon N_lat dV_East dV_North SITE')") ! header for columns to follow(?) IF (ios /= 0) EXIT chambers IF ((North_latitude < -90.D0).OR.(North_latitude > 90.D0)) THEN WRITE (*, "(' ERROR: Illegal North_latitude of ',F10.5)") North_latitude CALL Pause() STOP END IF IF (depth_in_km <= 0.0D0) THEN WRITE (*, "(' ERROR: Illegal (non-positive) depth of ',F10.5)") depth_in_km CALL Pause() STOP END IF peak_mm_per_year = ABS(peak_mm_per_year) ! just to be safe; when squared, this will become the eigenvalue! CALL DXyz_from_lonlat(East_longitude, North_latitude, volcano_uvec) !Limiting influence_radius is at perturbation of 0.01 mm/year, per Mogi solution: IF (peak_mm_per_year > 0.01D0) THEN reduction_factor = 0.01D0 / peak_mm_per_year influence_radius_km = (depth_in_km * 1.414D0) / DSQRT(reduction_factor) ! approximate radius to velocity perturbation of 0.01 mm/a influence_radius_radians = influence_radius_km / Earth_radius_km sine_influence_radius_radians = DSIN(influence_radius_radians) count_affected = 0 !On first pass, just count affected benchmarks... highest_velocity_mmpa = 0.0D0 DO i = 1, external_benchmarks uvec(1:3) = benchmark_uvec(1:3, i) CALL DCross(volcano_uvec, uvec, test_vector) sine = DMagnitude(test_vector) IF (sine <= sine_influence_radius_radians) THEN count_affected = count_affected + 1 r_km = DASIN(sine) * Earth_radius_km velocity_mmpa = peak_mm_per_year * ( r_km / ((r_km**2 + depth_in_km**2)**1.5D0) ) / & & ( (depth_in_km/1.414D0) / (((depth_in_km/1.414D0)**2 + depth_in_km**2)**1.5D0) ) highest_velocity_mmpa = MAX(highest_velocity_mmpa, velocity_mmpa) END IF END DO WRITE (*, "(' affects ',I4,' benchmarks; highest velocity = ',F10.4,' mm/a.')") count_affected, highest_velocity_mmpa !On second pass, record active benchmarks and their radial velocity components: ALLOCATE ( active_benchmarks(count_affected) ) ALLOCATE ( benchmark_motions(2, count_affected) ) n = 0 ! counter, 1 ... count_affected DO i = 1, external_benchmarks uvec(1:3) = benchmark_uvec(1:3, i) CALL DCross(volcano_uvec, uvec, test_vector) sine = DMagnitude(test_vector) IF (sine <= sine_influence_radius_radians) THEN n = n + 1 active_benchmarks(n) = i ! remember benchmark #s r_km = DASIN(sine) * Earth_radius_km velocity_mmpa = peak_mm_per_year * ( r_km / ((r_km**2 + depth_in_km**2)**1.5D0) ) / & & ( (depth_in_km/1.414D0) / (((depth_in_km/1.414D0)**2 + depth_in_km**2)**1.5D0) ) CALL DLocal_Theta(uvec, theta) ! unit vector pointing South CALL DLocal_Phi(uvec, phi) ! unit vector pointing East test_vector(1:3) = uvec(1:3) - volcano_uvec(1:3) ! quasi-horizontal vector, magnitude <<< 1.0 azimuth_radians = DATAN2(DDot(phi, test_vector), -DDot(theta, test_vector)) benchmark_motions(1, n) = velocity_mmpa * DSIN(azimuth_radians) ! East component in mm/year benchmark_motions(2, n) = velocity_mmpa * DCOS(azimuth_radians) ! North component in mm/year !Add memo to the log-file: del_V_East = benchmark_motions(1, n) ! in mm/a del_V_North = benchmark_motions(2, n) IF ((ABS(del_V_East) >= 0.001D0).OR.(ABS(del_V_North) >= 0.001D0)) THEN WRITE (51, "(F8.3, ' ', F7.3, ' ', 2F10.3, ' ', A)") benchmark_lon(i), benchmark_lat(i), del_V_East, del_V_North, TRIM(benchmark_name(i)) END IF END IF END DO !Write local transient gps file (for plotting): file_index = file_index + 1 WRITE (file_index_c3, "(I3)") file_index IF (file_index_c3(1:1) == ' ') file_index_c3(1:1) = '0' IF (file_index_c3(2:2) == ' ') file_index_c3(2:2) = '0' filename = "transient_V" // file_index_c3 // ".gps" WRITE (*, "(' Writing file ',A,' for plotting...')") TRIM(filename) CALL Write_transient_GPS(filename) WRITE (*, "()") !Actually increment the covariance matrix! DO i = 1, count_affected ! in short list of active benchmarks iN = 2 * active_benchmarks(i) iE = iN - 1 DO j = 1, count_affected ! in short list of active benchmarks jN = 2 * active_benchmarks(j) jE = jN - 1 covariance_mmpa2(iE, jE) = covariance_mmpa2(iE, jE) + benchmark_motions(1, i) * benchmark_motions(1, j) covariance_mmpa2(iE, jN) = covariance_mmpa2(iE, jN) + benchmark_motions(1, i) * benchmark_motions(2, j) covariance_mmpa2(iN, jE) = covariance_mmpa2(iN, jE) + benchmark_motions(2, i) * benchmark_motions(1, j) covariance_mmpa2(iN, jN) = covariance_mmpa2(iN, jN) + benchmark_motions(2, i) * benchmark_motions(2, j) END DO END DO DEALLOCATE ( active_benchmarks ) DEALLOCATE ( benchmark_motions ) END IF ! peak_mm_per_year > 0.01D0 WRITE (51, "('-------------------------------------------------------------------------------')") ! division between vectors END DO chambers END IF !=============================================================================================================== 120 WRITE (*, *) WRITE (*, "(' CAUTION: We have recently discovered some cases in which our method of')") WRITE (*, "(' augmented covariance leads to over-correction of geodetic transients by ~50%,')") WRITE (*, "(' causing errors in the long-term tectonic model which are opposite in sign')") WRITE (*, "(' (but smaller in magnitude) compared to the errors that would have occurred')") WRITE (*, "(' without augmentation. This problem can occur where: (a) the tectonic model')") WRITE (*, "(' is dependent on geodetic (as opposed to geologic) data to indicate the proper')") WRITE (*, "(' map-pattern of long-term deformation; and (b) the spatial pattern of transient')") WRITE (*, "(' velocities resembles a temporary amplification of long-term tectonic')") WRITE (*, "(' velocities. Examples of case (b) include post-seismic after-slip and/or')") WRITE (*, "(' viscoelastic relaxation. Therefore, we no longer recommend our method of')") WRITE (*, "(' augmented covariance as the optimal method for dealing with post-seismic')") WRITE (*, "(' after-slip and/or viscoelastic relaxation. Where these are suspected, it may')") WRITE (*, "(' be preferable to either: (1) create a quantitative model of the transient')") WRITE (*, "(' [e.g., E. Hearn, 2022, Seismological Research Letters, 93, 2973-2989,')") WRITE (*, "(' doi: 10.1785/0220220156] and subtract it from the geodetic velocities to be')") WRITE (*, "(' used for long-term tectonic modeling; or (2) simply edit the affected')") WRITE (*, "(' benchmarks out of the geodetic dataset to be used for long-term tectonic')") WRITE (*, "(' modeling. PB 2023.03.15')") WRITE (*, "(' ')") CALL Pause() CALL DPrompt_for_Logical('Do you have a file of post-seismic slip dislocations?', .TRUE., afterslipping) IF (afterslipping) THEN velocity_perturbation_CLASS = "POST-SEISMIC " ! for HEADER in GPS_Covariance_log.txt OPEN (UNIT = 44, FILE = "active_dislocation_patches.dig") ! can be plotted with NeoKineMap file_index = 0 ! but will be incremented for each local-transient F00n .GPS file written WRITE (*, *) WRITE (*, "(' Required contents of the afterslip data table (with only 1 case):')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, "(' {title line; perhaps file-name, and comments?}')") WRITE (*, "(' Elon Nlat Azimuth Dip Length Ztop Zbot Sinistral Opening Vertical')") WRITE (*, "(' (deg.) (deg.) (deg.) (deg.) (m) (m) (m) (mm/a) (mm/a) (mm/a)')") WRITE (*, "(' 13.350 42.290 318.0 130.0 28000. 0. 14000. 0.0 0.6427 -0.766')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, *) WRITE (*, "(' Where it is critical to observe the following conventions:')") WRITE (*, "(' -Coordinates are for center of fault trace, projected to surface.')") WRITE (*, "(' -East longitude is positive; West longitude is negative.')") WRITE (*, "(' -North latitude is positive; South latitude is negative.')") WRITE (*, "(' -Azimuth (of strike of patch) is clockwise from North.')") WRITE (*, "(' -Dip is measured from horizontal on right side of azimuth vector;')") WRITE (*, "(' therefore it may be greater than 90. if fault dips to left.')") WRITE (*, "(' -Length is the along-strike length of the slipping dislocation patch, in m.')") WRITE (*, "(' -Ztop (smaller) and Zbot (larger) are positive depths of the limits')") WRITE (*, "(' of the slipping dislocation patch, in m.')") WRITE (*, "(' -The 3 components of slip are all in mm/a.')") WRITE (*, "(' -Sinistral is strike-parallel slip, with left-lateral positive.')") WRITE (*, "(' -Opening is perpendicular horizontal displacement rate, extension positive.')") WRITE (*, "(' -The Vertical displacement rate is positive when right side of fault goes down.')") WRITE (*, *) CALL Pause() CALL DPrompt_for_String('[Drive:][\path\]filename for afterslip dislocation file:', 'afterslipping_dislocations.dat', afterslip_dislocation_pathfile) OPEN (UNIT = 12, FILE = TRIM(afterslip_dislocation_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Afterslipping-dislocation file not found (in this folder).')") CALL Pause() GO TO 120 END IF WRITE (*, *) READ (12, "(A)") dislocation_title WRITE(*, "(' ', A)") TRIM(dislocation_title(1:78)) READ (12, "(A)") dislocation_header WRITE(*, "(' ', A)") TRIM(dislocation_header(1:78)) READ (12, "(A)") dislocation_units WRITE(*, "(' ', A)") TRIM(dislocation_units(1:78)) patch_count = 0 patches: DO ! indefinite loop; terminated by reading EOF in input file READ (12, "(A)" , IOSTAT = ios) line IF (ios /= 0) THEN IF (ios < 0) THEN ! encountered EOF EXIT patches ELSE ! provide informative message WRITE (*, *) WRITE (*, "(' ERROR: While reading data regarding fault number ',I6,' in line ',I6,':')") (patch_count + 1), (patch_count + 4) WRITE (*, "(' ',A)") TRIM(line) WRITE (*, "(' the error condition IOSTAT = ',I6,' was returned.')") ios WRITE (*, "(' Please check the formats in your input file, after the 3 header lines.')") CALL Pause() STOP END IF ! ios > 0 END IF ! ios /= 0 WRITE (*, "(' ',A)") TRIM(line(1:78)) READ (line, *, IOSTAT = ios) Elon, Nlat, azimuth, dip, length, zTop, zBot, sinistral, opening, vertical IF (ios /= 0) THEN WRITE (*, *) WRITE (*, "(' ERROR: While reading data regarding fault number ',I6,' in line ',I6,':')") (patch_count + 1), (patch_count + 4) WRITE (*, "(' ',A)") TRIM(line) WRITE (*, "(' the error condition IOSTAT = ',I6,' was returned.')") ios WRITE (*, "(' Please check the formats in your input file, after the 3 header lines.')") CALL Pause() STOP END IF IF ((Nlat < -90.D0).OR.(Nlat > 90.D0)) THEN WRITE (*, "(' ERROR: Illegal North_latitude of ',F10.5)") Nlat CALL Pause() STOP END IF IF ((dip <= 0.D0).OR.(dip >= 180.D0)) THEN WRITE (*, "(' ERROR: Illegal dip of ',F10.5)") dip CALL Pause() STOP END IF IF (length <= 0.0D0) THEN WRITE (*, "(' ERROR: Illegal length of ',F10.5)") length CALL Pause() STOP END IF IF (zTop >= zBot) THEN WRITE (*, "(' ERROR: Illegal Ztop, Zbot of ',2F10.5)") zTop, zBot CALL Pause() STOP END IF patch_count = patch_count + 1 ! (only incremented after fault dislocation passes all error tests) !header in log-file WRITE (51, "(A15)") velocity_perturbation_CLASS WRITE (51, "('afterslip dislocation file = ', A)") TRIM(afterslip_dislocation_pathfile) ! source of information WRITE (51, "('Elon = ',F8.3,', Nlat = ',F7.3,', azimuth = ',F6.0,', dip = ',F5.1,', length = ',F7.0,', zTop = ',F6.0,', zBot = ',F6.0,' sinistral = ',F5.1,', opening = ',F5.1,', vertical = ',F5.1)") & & Elon, Nlat, azimuth, dip, length, zTop, zBot, sinistral, opening, vertical ! specific perturbation WRITE (51, "('East_lon N_lat dV_East dV_North SITE')") ! header for columns to follow(?) !Sample data: L'Aquila earthquake afterslip (assuming it is 1 mm/a), described approximately by: !Harv. CMT 2009.04.06 01:32:49.2 13.350 42.290 12 6.34 70 335 18 134 7 226 !corresponds to input-file data lines of: ! Elon Nlat Azimuth Dip Length Ztop Zbot Sinistral Opening Vertical Index ! (deg.) (deg.) (deg.) (deg.) (m) (m) (m) (mm/a) (mm/a) (mm/a) ! 13.350 42.290 318.0 130.0 28000. 0. 14000. 0.0 0.6427 -0.766 F001 argume = (180.0D0 - azimuth) / degrees_per_radian ! angle in radians, counterclockwise, from South to fault strike dipf = dip / degrees_per_radian ! measured from horizontal on R (e.g., 318 + 90 = 408 = 48 = NE) side of argume (strike) direction. ftheta = (90.0D0 - Nlat) / degrees_per_radian fphi = Elon / degrees_per_radian lf = 0.5D0 * length ! half-length of dislocation, measured along strike, in SI (m) like Earth_radius_m slip(1) = sinistral ! strike-slip, sinistral positive slip(2) = opening ! horizontal opening, perpendicular to strike slip(3) = vertical ! vertical motion (+ = up) on right (NE) side of fault. !Compute corners of dislocation patch, and output to the .DIG file on unit 44: !user-specified center-point (of projection of active dislocation plane to the surface): lon_dp = Elon ! convert to DOUBLE PRECISION lat_dp = Nlat CALL DLonLat_2_Uvec(lon_dp, lat_dp, center_uvec_dp) !uvec1_dp describes start-point of trace of dislocation patch plane, projected to surface: azimuth_radians_dp = azimuth / degrees_per_radian ! pointing forward from uvec1/3/5 toward uvec2/4/6 far_radians_dp = lf / Earth_radius_m CALL DTurn_To (azimuth_radians = (azimuth_radians_dp + Pi), & & base_uvec = center_uvec_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec1_dp) !uvec2_dp describes end-point of trace of dislocation patch plane, projected to surface: CALL DTurn_To (azimuth_radians = azimuth_radians_dp, & & base_uvec = center_uvec_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec2_dp) !project points to shallow depth along top of active dislocation patch: RHS_azimuth_radians_dp = azimuth_radians_dp + Pi_over_2 IF (ABS(DCOS(dip * radians_per_degree)) < 0.001D0) THEN ! avoid illegal function call DTAN(Pi/2) when fault is vertical. far_radians_dp = 0.0D0 ELSE far_radians_dp = (zTop / DTAN(dip * radians_per_degree)) / Earth_radius_m ! may be negative, if dip > 90. END IF !uvec3_dp describes start-point of trace of dislocation patch plane, at depth zTop: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec1_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec3_dp) !uvec4_dp describes end-point of trace of dislocation patch plane, at depth zTop: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec2_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec4_dp) !project points to deepest depth along base of active dislocation patch: IF (ABS(DCOS(dip * radians_per_degree)) < 0.001D0) THEN ! avoid illegal function call DTAN(Pi/2) when fault is vertical. far_radians_dp = 0.0D0 ELSE far_radians_dp = (zBot / DTAN(dip * radians_per_degree)) / Earth_radius_m ! may be negative, if dip > 90. END IF !uvec5_dp describes start-point of trace of dislocation patch plane, at depth zBot: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec1_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec5_dp) !uvec6_dp describes end-point of trace of dislocation patch plane, at depth zBot: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec2_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec6_dp) !output (projected) surface trace of dislocation patch, for plotting by NeoKineMap: CALL DUvec_2_LonLat(uvec1_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec2_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp WRITE (44, "('*** end of line segment ***')") !output outline of active dislocation patch, for plotting by NeoKineMap: CALL DUvec_2_LonLat(uvec3_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec5_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec6_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec4_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec3_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp WRITE (44, "('*** end of line segment ***')") !On first pass, just count affected benchmarks... !Locate and consider effects at nearby benchmarks: count_affected = 0 highest_velocity_mmpa = 0.0D0 wedge = 5.0D0 / degrees_per_radian ! faults of dip 85~95 will be treated as exactly vertical. DO i = 1, external_benchmarks btheta = (90.0D0 - benchmark_lat(i)) / degrees_per_radian bphi = benchmark_lon(i) / degrees_per_radian CALL DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & Earth_radius_m, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output ! This routine is a "driver" or "wrapper" for the ! Mansinha & Smylie routine Halo (vertical dislocation patch) ! or Aura (dipping dislocation patch) ! which permits them to be called from spherical-planet ! programs because it converts coordinates from ! (theta, phi) = (colatitude, longitude) in radians ! to the (z1,z2,z3) Cartesian units of a locally flat planet. ! Also, output (horizontal) vector (DUTHET, DUPHI) is expressed ! in the spherical-planet (theta,phi) system, so components are ! (South, East). Units are the same as input vector SLIP, ! which may be either a relative displacement or ! a relative velocity of the two sides of the fault, ! FTHETA and FPHI should give (theta, phi) coordinates of the ! midpoint of the trace of the plane containing the ! rectangular dislocation patch. ! BTHETA and BPHI should give (theta, phi) coordinates of the ! benchmark or test point at which the displacement (-rate) ! is desired. ! LF is the half-length (from center to end) of the dislocation ! patch, measured along a horizontal strike line. ! Units are the same as RADIUS, ZTOP, ZBOT (below). ! ARGUME is the argument of the trace of the dislocation patch, ! measured in radians counterclockwise from +theta (from South). ! DIPF is fault dip in radians measured clockwise ! (initially, down) from horizontal on the right side of the ! fault (when looking in direction ARGUME). ! SLIP must be already rotated into fault coordinates: ! SLIP(1) is the component parallel to the fault, ! and it is positive for sinistral sense of slip. ! SLIP(2) is the horizontal component in the direction ! perpendicular to the trace of the fault, ! and it is positive for divergence across the trace. ! SLIP(3) is the relative vertical component, and it ! is positive when the right side of the fault ! (when looking along direction ARGUME) is down. ! (Note that SLIP(1:3) is the only vector argument.) ! WEDGE is a tolerance (in radians) for dip; if DIPF is ! within WEDGE of Pi/2, then fault is considered vertical ! and routine Halo is called; otherwise, Aura is called. ! ZTOP and ZBOT are (positive) depths to top and bottom of ! the dislocation patch, respectively. Units are the same ! as RADIUS, which gives the size of the planet. ! NOTE: If distance from "F" point to "B" point exceeds ! (40 * ZBOT), result is approximated as zero and ! Halo/Aura are never called. velocity_mmpa = DSQRT(duthet**2 + duphi**2) IF (velocity_mmpa >= 0.001D0) count_affected = count_affected + 1 highest_velocity_mmpa = MAX(highest_velocity_mmpa, velocity_mmpa) END DO WRITE (*, "(' affects ',I4,' benchmarks; highest velocity = ',F10.4,' mm/a.')") count_affected, highest_velocity_mmpa !On second pass, record active benchmarks and their velocity components: ALLOCATE ( active_benchmarks(count_affected) ) ALLOCATE ( benchmark_motions(2, count_affected) ) n = 0 ! counter, 1 ... count_affected DO i = 1, external_benchmarks btheta = (90.0D0 - benchmark_lat(i)) / degrees_per_radian ! 270 = AQRA in devoti_2011.gps 13.374 42.366; on NE side bphi = benchmark_lon(i) / degrees_per_radian ! 270 = AQRA in devoti_2011.gps 13.374 42.366; on NE side CALL DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & Earth_radius_m, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output velocity_mmpa = DSQRT(duthet**2 + duphi**2) IF (velocity_mmpa >= 0.001D0) THEN n = n + 1 active_benchmarks(n) = i ! remember benchmark #s azimuth_radians = DATAN2(duphi, -duthet) ! clockwise, from North benchmark_motions(1, n) = velocity_mmpa * DSIN(azimuth_radians) ! East component in mm/year benchmark_motions(2, n) = velocity_mmpa * DCOS(azimuth_radians) ! North component in mm/year !Add a row to the log-file: del_V_East = benchmark_motions(1, n) ! in mm/a del_V_North = benchmark_motions(2, n) IF ((ABS(del_V_East) >= 0.001D0).OR.(ABS(del_V_North) >= 0.001D0)) THEN WRITE (51, "(F8.3, ' ', F7.3, ' ', 2F10.3, ' ', A)") benchmark_lon(i), benchmark_lat(i), del_V_East, del_V_North, TRIM(benchmark_name(i)) END IF END IF END DO ! i = 1, external_benchmarks WRITE (51, "('-------------------------------------------------------------------------------')") ! division between vectors !Actually increment the covariance matrix! DO i = 1, count_affected ! in short list of active benchmarks iN = 2 * active_benchmarks(i) iE = iN - 1 DO j = 1, count_affected ! in short list of active benchmarks jN = 2 * active_benchmarks(j) jE = jN - 1 covariance_mmpa2(iE, jE) = covariance_mmpa2(iE, jE) + benchmark_motions(1, i) * benchmark_motions(1, j) covariance_mmpa2(iE, jN) = covariance_mmpa2(iE, jN) + benchmark_motions(1, i) * benchmark_motions(2, j) covariance_mmpa2(iN, jE) = covariance_mmpa2(iN, jE) + benchmark_motions(2, i) * benchmark_motions(1, j) covariance_mmpa2(iN, jN) = covariance_mmpa2(iN, jN) + benchmark_motions(2, i) * benchmark_motions(2, j) END DO END DO !Write local transient gps file (for plotting): file_index = file_index + 1 WRITE (file_index_c3, "(I3)") file_index IF (file_index_c3(1:1) == ' ') file_index_c3(1:1) = '0' IF (file_index_c3(2:2) == ' ') file_index_c3(2:2) = '0' filename = "transient_F" // file_index_c3 // ".gps" WRITE (*, "(' Writing file ',A,' for plotting...')") TRIM(filename) CALL Write_transient_GPS(filename) WRITE (*, "()") DEALLOCATE ( active_benchmarks ) DEALLOCATE ( benchmark_motions ) END DO patches CLOSE (UNIT = 44) ! active_dislocation_patches.dig END IF ! afterslipping !=============================================================================================================== !GPBstart 125 WRITE (*, *) CALL DPrompt_for_Logical('Do you have a file of opening/closing (mode-1) cracks?', .TRUE., openingCracks) IF (openingCracks) THEN velocity_perturbation_CLASS = "OPENING CRACK " ! for HEADER in GPS_Covariance_log.txt OPEN (UNIT = 44, FILE = "opening_crack_patches.dig") ! can be plotted with NeoKineMap file_index = 0 ! but will be incremented for each local-transient C00n .GPS file written WRITE (*, *) WRITE (*, "(' Required contents of the mode-1 crack data table (with only 1 case):')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, "(' {title line; perhaps file-name, and comments?}')") WRITE (*, "(' Elon Nlat Azimuth Dip Length Ztop Zbot Expansion')") WRITE (*, "(' (deg.) (deg.) (deg.) (deg.) (m) (m) (m) (mm/a)')") WRITE (*, "(' 13.050 42.570 325.0 90.0 60000. 0. 8000. 20.0')") WRITE (*, "(' ----------------------------------------------------------------')") WRITE (*, *) WRITE (*, "(' Where it is critical to observe the following conventions:')") WRITE (*, "(' -Coordinates are for center of crack trace, after plane is projected to surface.')") WRITE (*, "(' -East longitude is positive; West longitude is negative.')") WRITE (*, "(' -North latitude is positive; South latitude is negative.')") WRITE (*, "(' -Azimuth (of strike of crack) is clockwise from North.')") WRITE (*, "(' -Dip is measured from horizontal on right side of azimuth vector;')") WRITE (*, "(' therefore it may be greater than 90. if crack dips to left.')") WRITE (*, "(' -Length is the along-strike length of the crack at depth, in m.')") WRITE (*, "(' -Ztop (smaller) and Zbot (larger) are positive depths of the limits')") WRITE (*, "(' of the opening/closing crack patch, in m.')") WRITE (*, "(' -The crack Expansion rate is in mm/a. This can be negative (for closing).')") WRITE (*, *) CALL Pause() CALL DPrompt_for_String('[Drive:][\path\]filename for mode-1 crack file:', 'opening_crack_patches.dat', opening_crack_patches_pathfile) OPEN (UNIT = 12, FILE = TRIM(opening_crack_patches_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Opening-crack patches data file not found (in this folder).')") CALL Pause() GO TO 125 END IF WRITE (*, *) READ (12, "(A)") crack_title WRITE(*, "(' ', A)") TRIM(crack_title(1:78)) READ (12, "(A)") crack_header WRITE(*, "(' ', A)") TRIM(crack_header(1:78)) READ (12, "(A)") crack_units WRITE(*, "(' ', A)") TRIM(crack_units(1:78)) crack_count = 0 cracks: DO ! indefinite loop; terminated by reading EOF in input file READ (12, "(A)" , IOSTAT = ios) line IF (ios /= 0) THEN IF (ios < 0) THEN ! encountered EOF EXIT cracks ELSE ! provide informative message WRITE (*, *) WRITE (*, "(' ERROR: While reading data regarding crack number ',I6,' in line ',I6,':')") (crack_count + 1), (crack_count + 4) WRITE (*, "(' ',A)") TRIM(line) WRITE (*, "(' the error condition IOSTAT = ',I6,' was returned.')") ios WRITE (*, "(' Please check the formats in your input file, after the 3 header lines.')") CALL Pause() STOP END IF ! ios > 0 END IF ! ios /= 0 WRITE (*, "(' ',A)") TRIM(line(1:78)) READ (line, *, IOSTAT = ios) Elon, Nlat, azimuth, dip, length, zTop, zBot, expansion IF (ios /= 0) THEN WRITE (*, *) WRITE (*, "(' ERROR: While reading data regarding crack number ',I6,' in line ',I6,':')") (crack_count + 1), (crack_count + 4) WRITE (*, "(' ',A)") TRIM(line) WRITE (*, "(' the error condition IOSTAT = ',I6,' was returned.')") ios WRITE (*, "(' Please check the formats in your input file, after the 3 header lines.')") CALL Pause() STOP END IF IF ((Nlat < -90.D0).OR.(Nlat > 90.D0)) THEN WRITE (*, "(' ERROR: Illegal North_latitude of ',F10.5)") Nlat CALL Pause() STOP END IF IF ((dip <= 0.D0).OR.(dip >= 180.D0)) THEN WRITE (*, "(' ERROR: Illegal dip of ',F10.5)") dip CALL Pause() STOP END IF IF (length <= 0.0D0) THEN WRITE (*, "(' ERROR: Illegal length of ',F10.5)") length CALL Pause() STOP END IF IF (zTop >= zBot) THEN WRITE (*, "(' ERROR: Illegal Ztop, Zbot of ',2F10.5)") zTop, zBot CALL Pause() STOP END IF crack_count = crack_count + 1 ! (only incremented after crack patch passes all error tests) !header in log-file WRITE (51, "(A15)") velocity_perturbation_CLASS ! == "OPENING CRACK " WRITE (51, "('opening crack patches file = ', A)") TRIM(opening_crack_patches_pathfile) ! source of information WRITE (51, "('Elon = ',F8.3,', Nlat = ',F7.3,', azimuth = ',F6.0,', dip = ',F5.1,', length = ',F7.0,', zTop = ',F6.0,', zBot = ',F6.0,' expansion = ',F5.1)") & & Elon, Nlat, azimuth, dip, length, zTop, zBot, expansion ! specific perturbation WRITE (51, "('East_lon N_lat dV_East dV_North SITE')") ! header for columns to follow(?) argume = (180.0D0 - azimuth) / degrees_per_radian ! angle in radians, counterclockwise, from South to fault strike dipf = dip / degrees_per_radian ! measured from horizontal on R side of argume (strike) direction. ftheta = (90.0D0 - Nlat) / degrees_per_radian fphi = Elon / degrees_per_radian lf = 0.5D0 * length ! half-length of dislocation, measured along strike, in SI (m) like Earth_radius_m slip(1) = 0.0D0 ! = sinistral strike-slip (none for an opening crack) slip(2) = +expansion * SIN(dipf) ! = horizontal opening, perpendicular to strike slip(3) = -expansion * COS(dipf) ! = vertical motion (+ = down) on right side of crack. (Usually negative for positive expansion.) !Compute corners of opening crack patch, and output to the .DIG file on unit 44: !user-specified center-point (of projection of active dislocation plane to the surface): lon_dp = Elon ! convert to DOUBLE PRECISION lat_dp = Nlat CALL DLonLat_2_Uvec(lon_dp, lat_dp, center_uvec_dp) !uvec1_dp describes start-point of trace of opening patch plane, projected to surface: azimuth_radians_dp = azimuth / degrees_per_radian ! pointing forward from uvec1/3/5 toward uvec2/4/6 far_radians_dp = lf / Earth_radius_m CALL DTurn_To (azimuth_radians = (azimuth_radians_dp + Pi), & & base_uvec = center_uvec_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec1_dp) !uvec2_dp describes end-point of trace of opening patch plane, projected to surface: CALL DTurn_To (azimuth_radians = azimuth_radians_dp, & & base_uvec = center_uvec_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec2_dp) !project points to shallow depth along top of opening patch: RHS_azimuth_radians_dp = azimuth_radians_dp + Pi_over_2 IF (ABS(DCOS(dip * radians_per_degree)) < 0.001D0) THEN ! avoid illegal function call DTAN(Pi/2) when crack is vertical. far_radians_dp = 0.0D0 ELSE far_radians_dp = (zTop / DTAN(dip * radians_per_degree)) / Earth_radius_m ! may be negative, if dip > 90. END IF !uvec3_dp describes start-point of trace of opening patch plane, at depth zTop: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec1_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec3_dp) !uvec4_dp describes end-point of trace of opening patch plane, at depth zTop: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec2_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec4_dp) !project points to deepest depth along base of active opening patch: IF (ABS(DCOS(dip * radians_per_degree)) < 0.001D0) THEN ! avoid illegal function call DTAN(Pi/2) when crack is vertical. far_radians_dp = 0.0D0 ELSE far_radians_dp = (zBot / DTAN(dip * radians_per_degree)) / Earth_radius_m ! may be negative, if dip > 90. END IF !uvec5_dp describes start-point of trace of opening patch plane, at depth zBot: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec1_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec5_dp) !uvec6_dp describes end-point of trace of opening patch plane, at depth zBot: CALL DTurn_To (azimuth_radians = RHS_azimuth_radians_dp, & & base_uvec = uvec2_dp, & & far_radians = far_radians_dp, & ! inputs & omega_uvec = omega_uvec_dp, result_uvec = uvec6_dp) !output (projected) surface trace of opening patch, for plotting by NeoKineMap: CALL DUvec_2_LonLat(uvec1_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec2_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp WRITE (44, "('*** end of line segment ***')") !output outline of opening/closing crack patch, for plotting by NeoKineMap: CALL DUvec_2_LonLat(uvec3_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec5_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec6_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec4_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp CALL DUvec_2_LonLat(uvec3_dp, lon_dp, lat_dp) WRITE (44, "(1X,SP,ES12.5,',',ES12.5)") lon_dp, lat_dp WRITE (44, "('*** end of line segment ***')") !On first pass, just count affected benchmarks... !Locate and consider effects at nearby benchmarks: count_affected = 0 highest_velocity_mmpa = 0.0D0 wedge = 5.0D0 / degrees_per_radian ! cracks of dip 85~95 degrees will be treated as exactly vertical (to avoid singularities in M&S formulas). DO i = 1, external_benchmarks btheta = (90.0D0 - benchmark_lat(i)) / degrees_per_radian bphi = benchmark_lon(i) / degrees_per_radian CALL DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & Earth_radius_m, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output ! This routine is a "driver" or "wrapper" for the ! Mansinha & Smylie routine Halo (vertical dislocation patch) ! or Aura (dipping dislocation patch) ! which permits them to be called from spherical-planet ! programs because it converts coordinates from ! (theta, phi) = (colatitude, longitude) in radians ! to the (z1,z2,z3) Cartesian units of a locally flat planet. ! Also, output (horizontal) vector (DUTHET, DUPHI) is expressed ! in the spherical-planet (theta,phi) system, so components are ! (South, East). Units are the same as input vector SLIP, ! which may be either a relative displacement or ! a relative velocity of the two sides of the fault, ! FTHETA and FPHI should give (theta, phi) coordinates of the ! midpoint of the trace of the plane containing the ! rectangular dislocation patch. ! BTHETA and BPHI should give (theta, phi) coordinates of the ! benchmark or test point at which the displacement (-rate) ! is desired. ! LF is the half-length (from center to end) of the dislocation ! patch, measured along a horizontal strike line. ! Units are the same as RADIUS, ZTOP, ZBOT (below). ! ARGUME is the argument of the trace of the dislocation patch, ! measured in radians counterclockwise from +theta (from South). ! DIPF is fault dip in radians measured clockwise ! (initially, down) from horizontal on the right side of the ! fault (when looking in direction ARGUME). ! SLIP must be already rotated into fault coordinates: ! SLIP(1) is the component parallel to the fault, ! and it is positive for sinistral sense of slip. ! SLIP(2) is the horizontal component in the direction ! perpendicular to the trace of the fault, ! and it is positive for divergence across the trace. ! SLIP(3) is the relative vertical component, and it ! is positive when the right side of the fault ! (when looking along direction ARGUME) is down. ! (Note that SLIP(1:3) is the only vector argument.) ! WEDGE is a tolerance (in radians) for dip; if DIPF is ! within WEDGE of Pi/2, then fault is considered vertical ! and routine Halo is called; otherwise, Aura is called. ! ZTOP and ZBOT are (positive) depths to top and bottom of ! the dislocation patch, respectively. Units are the same ! as RADIUS, which gives the size of the planet. ! NOTE: If distance from "F" point to "B" point exceeds ! (40 * ZBOT), result is approximated as zero and ! Halo/Aura are never called. velocity_mmpa = DSQRT(duthet**2 + duphi**2) IF (velocity_mmpa >= 0.001D0) count_affected = count_affected + 1 highest_velocity_mmpa = MAX(highest_velocity_mmpa, velocity_mmpa) END DO WRITE (*, "(' affects ',I4,' benchmarks; highest velocity = ',F10.4,' mm/a.')") count_affected, highest_velocity_mmpa !On second pass, record active benchmarks and their velocity components: ALLOCATE ( active_benchmarks(count_affected) ) ALLOCATE ( benchmark_motions(2, count_affected) ) n = 0 ! counter, 1 ... count_affected DO i = 1, external_benchmarks btheta = (90.0D0 - benchmark_lat(i)) / degrees_per_radian ! 270 = AQRA in devoti_2011.gps 13.374 42.366; on NE side bphi = benchmark_lon(i) / degrees_per_radian ! 270 = AQRA in devoti_2011.gps 13.374 42.366; on NE side CALL DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & Earth_radius_m, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output velocity_mmpa = DSQRT(duthet**2 + duphi**2) IF (velocity_mmpa >= 0.001D0) THEN n = n + 1 active_benchmarks(n) = i ! remember benchmark #s azimuth_radians = DATAN2(duphi, -duthet) ! clockwise, from North benchmark_motions(1, n) = velocity_mmpa * DSIN(azimuth_radians) ! East component in mm/year benchmark_motions(2, n) = velocity_mmpa * DCOS(azimuth_radians) ! North component in mm/year !Add a row to the log-file: del_V_East = benchmark_motions(1, n) ! in mm/a del_V_North = benchmark_motions(2, n) IF ((ABS(del_V_East) >= 0.001D0).OR.(ABS(del_V_North) >= 0.001D0)) THEN WRITE (51, "(F8.3, ' ', F7.3, ' ', 2F10.3, ' ', A)") benchmark_lon(i), benchmark_lat(i), del_V_East, del_V_North, TRIM(benchmark_name(i)) END IF END IF END DO ! i = 1, external_benchmarks WRITE (51, "('-------------------------------------------------------------------------------')") ! division between vectors !Actually increment the covariance matrix! DO i = 1, count_affected ! in short list of active benchmarks iN = 2 * active_benchmarks(i) iE = iN - 1 DO j = 1, count_affected ! in short list of active benchmarks jN = 2 * active_benchmarks(j) jE = jN - 1 covariance_mmpa2(iE, jE) = covariance_mmpa2(iE, jE) + benchmark_motions(1, i) * benchmark_motions(1, j) covariance_mmpa2(iE, jN) = covariance_mmpa2(iE, jN) + benchmark_motions(1, i) * benchmark_motions(2, j) covariance_mmpa2(iN, jE) = covariance_mmpa2(iN, jE) + benchmark_motions(2, i) * benchmark_motions(1, j) covariance_mmpa2(iN, jN) = covariance_mmpa2(iN, jN) + benchmark_motions(2, i) * benchmark_motions(2, j) END DO END DO !Write local transient gps file (for plotting): file_index = file_index + 1 WRITE (file_index_c3, "(I3)") file_index IF (file_index_c3(1:1) == ' ') file_index_c3(1:1) = '0' IF (file_index_c3(2:2) == ' ') file_index_c3(2:2) = '0' filename = "transient_C" // file_index_c3 // ".gps" WRITE (*, "(' Writing file ', A, ' for plotting...')") TRIM(filename) CALL Write_transient_GPS(filename) WRITE (*, "()") DEALLOCATE ( active_benchmarks ) DEALLOCATE ( benchmark_motions ) END DO cracks CLOSE (UNIT = 44) ! opening_crack_patches.dig END IF ! openingCracks !GPBend !=============================================================================================================== 130 WRITE (*, *) CALL DPrompt_for_Logical('Do you have any EXTRA transient motions to add?', .TRUE., extras) IF (extras) THEN velocity_perturbation_CLASS = "EXTRA " ! for HEADER in GPS_Covariance_log.txt WRITE (*, *) WRITE (*, "(' Requirements for each ""Extra"" transient-GPS file:')") WRITE (*, "(' *Each contains benchmark velocities due to only ONE scalar source.')") WRITE (*, "(' *Does not contain any secular tectonic strain accumulation.')") WRITE (*, *) WRITE (*, "(' Convenient conventions for ""Extra"" transient-GPS files:')") WRITE (*, "(' *The sign (sense) of motions is not important,')") WRITE (*, "(' as long as all simultaneous motions have consistent sense.')") WRITE (*, "(' *Uncertainties of transient motions are not used, and may')") WRITE (*, "(' all be set to "" 0.0 0.0 0.0"", if you like.')") WRITE (*, *) extra_files = 1 ! (?) CALL DPrompt_for_Integer('How many transient-GPS files do you want to read?', extra_files, extra_files) IF (extra_files > 0) THEN additions : DO file_index = 1, extra_files ! process requested number of Extra transient-GPS input files: WRITE (file_index_c3, "(I3)") file_index IF (file_index_c3(1:1) == ' ') file_index_c3(1:1) = '0' IF (file_index_c3(2:2) == ' ') file_index_c3(2:2) = '0' filename = "transient_E" // file_index_c3 // ".gps" WRITE (*, *) CALL DPrompt_for_String('[Drive:][\path\]filename for transient-GPS data file:', filename, input_transient_GPS_file) OPEN (UNIT = 13, FILE = TRIM(input_transient_GPS_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Extra transient-GPS file not found (in this folder).')") CALL Pause() GO TO 130 END IF !First reading is just to count the number of active benchmarks: READ (13, "(A)") gps_title WRITE (*, "(' ',A)") TRIM(gps_title) READ (13, "(A)") gps_format READ (13, "(A)") gps_header count_affected = 0 perusing: DO READ (13, gps_format, IOSTAT = ios) lon, lat, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation, c80, alleged_benchmark_name IF (ios /= 0) EXIT perusing count_affected = count_affected + 1 END DO perusing CLOSE (13) WRITE (*, "(' ',I6,' benchmarks found in this file.')") count_affected ALLOCATE ( active_benchmarks(count_affected) ) ALLOCATE ( benchmark_motions(2, count_affected) ) !Second pass memorizes critical parts of file contents: OPEN (UNIT = 13, FILE = TRIM(input_transient_GPS_file), STATUS = "OLD", IOSTAT = ios) READ (13, "(A)") gps_title READ (13, "(A)") gps_format READ (13, "(A)") gps_header i = 0 reusing: DO READ (13, gps_format, IOSTAT = ios) lon, lat, vE_mmpa, vN_mmpa, vE_sigma, vN_sigma, correlation, c80, alleged_benchmark_name IF (ios /= 0) EXIT reusing i = i + 1 !Identify which benchmark this is, in the global lists: n = 0 matching: DO j = 1, external_benchmarks IF (TRIM(alleged_benchmark_name) == TRIM(benchmark_name(j))) THEN IF (lon == benchmark_lon(j)) THEN IF (lat == benchmark_lat(j)) THEN n = j EXIT matching END IF ! lat matches END IF ! lon matches END IF ! name matches END DO matching IF (n == 0) THEN ! should not happen; matching failed! WRITE (*, "(' ERROR: Benchmark named ',A,' at (',F10.3,', ',F10.3,')')") TRIM(alleged_benchmark_name), lon, lat WRITE (*, "(' was not found in the current master list; OR, this name was found, but its')") WRITE (*, "(' (lon, lat) coordinates were not exactly the same in the master list.')") CALL Pause() STOP END IF ! matching failed !If we get to here, then matching succeeded... active_benchmarks(i) = n benchmark_motions(1, i) = vE_mmpa benchmark_motions(2, i) = vN_mmpa END DO reusing !header in log-file WRITE (51, "(A15)") velocity_perturbation_CLASS WRITE (51, "('input transient GPS file = ', A)") TRIM(input_transient_GPS_file) ! source of information WRITE (51, "('This file describes coordinated transient motions at ',I5,' benchmarks. The fast motions are:')") count_affected WRITE (51, "('East_lon N_lat dV_East dV_North SITE')") ! header for columns to follow(?) !Add as many rows as needed to the log-file: DO i = 1, count_affected del_V_East = benchmark_motions(1, i) ! in mm/a del_V_North = benchmark_motions(2, i) n = active_benchmarks(i) ! sequence # of the moving benchmark in the larger list of length external_benchmarks IF ((ABS(del_V_East) >= 0.001D0).OR.(ABS(del_V_North) >= 0.001D0)) THEN WRITE (51, "(F8.3, ' ', F7.3, ' ', 2F10.3, ' ', A)") benchmark_lon(n), benchmark_lat(n), del_V_East, del_V_North, TRIM(benchmark_name(n)) END IF END DO WRITE (51, "('-------------------------------------------------------------------------------')") ! division between vectors !Actually increment the covariance matrix! DO i = 1, count_affected ! in short list of active benchmarks iN = 2 * active_benchmarks(i) iE = iN - 1 DO j = 1, count_affected ! in short list of active benchmarks jN = 2 * active_benchmarks(j) jE = jN - 1 covariance_mmpa2(iE, jE) = covariance_mmpa2(iE, jE) + benchmark_motions(1, i) * benchmark_motions(1, j) covariance_mmpa2(iE, jN) = covariance_mmpa2(iE, jN) + benchmark_motions(1, i) * benchmark_motions(2, j) covariance_mmpa2(iN, jE) = covariance_mmpa2(iN, jE) + benchmark_motions(2, i) * benchmark_motions(1, j) covariance_mmpa2(iN, jN) = covariance_mmpa2(iN, jN) + benchmark_motions(2, i) * benchmark_motions(2, j) END DO END DO DEALLOCATE ( active_benchmarks ) DEALLOCATE ( benchmark_motions ) END DO additions END IF ! extra_files > 0 END IF ! extras !=============================================================================================================== !FINISHED adding new types of covariance; !Output the new, improved .GP2 covariance file: WRITE (*, *) WRITE (*, "(' Ready to output the new, improved .GP2 covariance file...')") 30 IF (read_GP2) THEN ! old_GP2_pathfile is defined; expand it: last_slash = INDEX(old_GP2_pathfile, '\', .TRUE.) IF (last_slash == 0) THEN new_GP2_pathfile = "expanded_" // TRIM(old_GP2_pathfile) ELSE new_GP2_pathfile = old_GP2_pathfile(1:last_slash) // "expanded_" // TRIM(old_GP2_pathfile((last_slash+1):134)) END IF ELSE new_GP2_pathfile = "expanded.GP2" END IF CALL DPrompt_for_String('[Drive:][\path\]filename for new, expanded .GP2 file:', new_GP2_pathfile, new_GP2_pathfile) IF (TRIM(new_GP2_pathfile) == TRIM(old_GP2_pathfile)) THEN WRITE (*, "(' ERROR: New .GP2 filename must be different from old .GP2 filename.')") CALL Pause() GO TO 30 END IF OPEN (UNIT = 3, FILE = TRIM(new_GP2_pathfile), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Cannot create file with this name. (Perhaps it already exists?)')") CALL Pause() GO TO 30 END IF !Write out upper-triangular entries ONLY to save filesize... DO i = 1, geodetic_NDOF DO j = i, geodetic_NDOF IF (covariance_mmpa2(i, j) /= 0.0D0) WRITE (3, "(I8, I8, ES12.4)") i, j, covariance_mmpa2(i, j) END DO END DO CLOSE (3) WRITE (*, "(' Writing of new, expanded .GP2 file was completed.')") WRITE (51, "('===============================================================================')") ! end of data memos; beginning of footer WRITE (51, "('.GP2 output file = ', A)") TRIM(new_GP2_pathfile) CALL Pause() !=============================================================================================================== WRITE (*, *) WRITE (*, "(' Job completed.')") CLOSE (51) ! GPS_Covariance_log.txt CALL Pause() CONTAINS CHARACTER*10 FUNCTION DASCII10(x) ! Returns a right-justified 10-byte (or shorter) version of a ! floating-point number, in "human" format, with no more ! than 4 significant digits. IMPLICIT NONE REAL*8, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL*8 :: x_log DOUBLE PRECISION :: y IF (x == 0.0D0) THEN DASCII10=' 0' RETURN ELSE IF (x > 0.0D0) THEN punt = (x >= 999999999.5D0).OR.(x < 0.0000100D0) ELSE ! x < 0.0 punt = (x <= -99999999.5D0).OR.(x > -0.000100D0) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,E10.3)') x !consider possible improvements, from left to right: IF (temp10(3:6) == '.000') THEN ! right-shift 4 spaces over it temp20(7:10) = temp10(7:10) temp20(5:6) = temp10(1:2) temp20(1:4) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(5:6) == '00') THEN ! right-shift 2 spaces over it temp20(7:10) = temp10(7:10) temp20(3:6) = temp10(1:4) temp20(1:2) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(6:6) == '0') THEN ! right-shift 1 space over it temp20(7:10) = temp10(7:10) temp20(2:6) = temp10(1:5) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(8:8) == '+') THEN ! right-shift over + sign in exponent temp20(9:10) = temp10(9:10) temp20(2:8) = temp10(1:7) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(9:9) == '0') THEN ! right-shift over leading 0 in exponent temp20(10:10) = temp10(10:10) temp20(2:9) = temp10(1:8) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF DASCII10 = temp10 ELSE ! can represent without exponential notation x_log = DLOG10(ABS(x)) zeros = DInt_Below(x_log) - 3 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0D0) y = -y WRITE (temp20,"(F20.9)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp20(19:20) == '01') temp20(19:20) = '00' !Find first important byte from right; change 0 -> ' ' k10 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 20, 12, -1 IF (temp20(j:j) == '0') THEN temp20(j:j) = ' ' ELSE k10 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0) THEN IF (temp20(10:11) == ' .') temp20(10:11) = '0.' ELSE ! x < 0.0 IF (k10 <= 18) THEN IF (temp20(9:11) == ' -.') temp20(9:11) = '-0.' END IF END IF k1 = k10 - 9 DASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION DASCII10 INTEGER FUNCTION DInt_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER :: i REAL*8 :: y i = INT(x) IF (x >= 0.0D0) THEN DInt_Below = i ELSE ! x < 0.0D0 y = 1.0D0*i IF (y <= x) THEN DInt_Below = i ELSE ! most commonly DInt_Below = i - 1 END IF END IF END FUNCTION DInt_Below SUBROUTINE 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 SUBROUTINE DPrompt_for_Integer (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with an integer value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text INTEGER, INTENT(IN) :: default INTEGER, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, trial, written LOGICAL :: finished WRITE (suggested,"(I11)") default suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an integer of 9 digits or less.')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Integer SUBROUTINE DPrompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE DPrompt_for_Logical SUBROUTINE DPrompt_for_Real (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a real value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text REAL*8, INTENT(IN) :: default REAL*8, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL*8 :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use DASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = DASCII10(default) !(2) Be sure that the number contains some sign that it is floating-point, not integer: IF (INDEX(suggested, '.') == 0) THEN IF ((INDEX(suggested, 'E') == 0).AND.(INDEX(suggested, 'e') == 0).AND. & & (INDEX(suggested, 'D') == 0).AND.(INDEX(suggested, 'd') == 0)) THEN suggested = ADJUSTL(suggested) point = LEN_TRIM(suggested) + 1 suggested(point:point) = '.' END IF END IF !------------------------------------------------------------------------------------ suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an real number using 11 characters (or less).')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Real SUBROUTINE DPrompt_for_String (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a character-string value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text CHARACTER*(*), INTENT(IN) :: default CHARACTER*(*), INTENT(OUT) :: answer CHARACTER*80 trial INTEGER :: blank_at, default_bytes, leftover, & & prompt_bytes, written prompt_bytes = LEN_TRIM(prompt_text) default_bytes = LEN_TRIM(default) written = 0 leftover = 79 - prompt_bytes - 4 ! unless changed below DO WHILE ((prompt_bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) leftover = 79 - (blank_at - (written+1) + 1) - 4 written = blank_at END DO IF (leftover >= default_bytes) THEN WRITE (*,"(' ',A,' [',A,']')") prompt_text(written+1:prompt_bytes), TRIM(default) ELSE WRITE (*,"(' ',A)") prompt_text(written+1:prompt_bytes) WRITE (*,"(' [',A,']')") TRIM(default) END IF WRITE (*,"(' ?: '\)") READ (*,"(A)") trial IF (LEN_TRIM(trial) == 0) THEN answer = TRIM(default) ELSE answer = TRIM(trial) END IF END SUBROUTINE DPrompt_for_String SUBROUTINE Write_transient_GPS(filename) IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: filename INTEGER :: i !All other data is in global variables: ! gps_title, gps_format, gps_header, ! benchmark_lon, benchmark_lat, benchmark_vE_sigma, benchmark_vN_sigma, benchmark_correlation, ! count_affected, active_benchmarks, benchmark_motions OPEN (UNIT = 33, FILE = TRIM(filename)) gps_title = "Local-transient GPS file " // TRIM(filename) // " from PROGRAM GPS_covariance" WRITE (33, "(A)") TRIM(gps_title) WRITE (33, "(A)") TRIM(gps_format) WRITE (33, "(A)") TRIM(gps_header) c80 = ' ' DO i = 1, count_affected WRITE (33, gps_format) benchmark_lon(active_benchmarks(i)), benchmark_lat(active_benchmarks(i)), & & benchmark_motions(1, i), benchmark_motions(2, i), & & benchmark_vE_sigma(active_benchmarks(i)), benchmark_vN_sigma(active_benchmarks(i)), benchmark_correlation(active_benchmarks(i)), & & c80, benchmark_name(active_benchmarks(i)) END DO CLOSE (33) END SUBROUTINE Write_transient_GPS SUBROUTINE DXyz_from_LonLat (lon, lat, vector) REAL*8, INTENT(IN) :: lon, lat REAL*8, DIMENSION(3), INTENT(OUT) :: vector ! assuming longitude is East longitude in degrees, ! and latitude is North latitude in degrees, ! computes 3 components of Cartesian unit vector ! from center of planet to surface point (on unit sphere). REAL*8, PARAMETER :: deg_per_rad = 57.2958D0 REAL*8 :: theta_, phi_, equat theta_ = (90.D0 - lat) / deg_per_rad phi_ = lon / deg_per_rad equat = DSIN(theta_) vector(1) = equat * DCOS(phi_) vector(2) = equat * DSIN(phi_) vector(3) = DCOS(theta_) END SUBROUTINE DXyz_from_LonLat END PROGRAM GPS_Covariance