PROGRAM NK_Interseismic_Vs_At_Benchmarks ! Utility program in the NeoKinema family; ! used to compute and report NeoKinema model interseismic velocities ! at GPS geodetic benchmarks. ! The method is to assume that the model errors are the same in ! the interseismic domain (of *.gps) and in the long-term domain ! (of g_[token].nko), and transfer these errors across files, ! for all benchmarks that were actually used in the final solution. ! Note that this is more accurate than a previous utility ! (now deprecated) which made the estimate by linear interpolation ! in space from v_interseismic_[token].out; the problem is that the ! fault locking/unlocking correction is very nonlinear in space ! near major faults. ! By Peter Bird, UCLA, 2012.01.19 IMPLICIT NONE CHARACTER*132 :: GPS_filename = ' ', GPS_titles = ' ', GPS_FORMAT = ' ', GPS_headers = ' ', & & g_token_NKO_filename = ' ', g_token_NKO_titles = ' ', g_token_NKO_FORMAT = ' ', g_token_NKO_headers = ' ', & & g_interseismic_token_NKO_filename = ' ', & & model_name_token = ' ' CHARACTER*20 :: reference_frame, site CHARACTER*20, DIMENSION(:), ALLOCATABLE :: GPS_site CHARACTER*20, DIMENSION(:), ALLOCATABLE :: g_token_NKO_site INTEGER :: GPS_count, g_token_NKO_count, i, ios, j, match REAL*8 :: E_lon_deg, N_lat_deg, v_E_mmpa, v_N_mmpa, v_E_sigma, v_N_sigma, correlation, E_error_mmpa, N_error_mmpa, & & interseismic_E_mmpa, interseismic_N_mmpa REAL*8, DIMENSION (:), ALLOCATABLE :: GPS_E_lon_deg, GPS_N_lat_deg REAL*8, DIMENSION (:), ALLOCATABLE :: GPS_v_E_mmpa, GPS_v_N_mmpa REAL*8, DIMENSION (:), ALLOCATABLE :: g_token_NKO_E_lon_deg, g_token_NKO_N_lat_deg REAL*8, DIMENSION (:), ALLOCATABLE :: g_token_NKO_E_error_mmpa, g_token_NKO_N_error_mmpa WRITE (*, "(' PROGRAM NK_Interseismic_Vs_At_Benchmarks')") WRITE (*, "(' ')") WRITE (*, "(' Utility program in the NeoKinema family;')") WRITE (*, "(' used to compute and report NeoKinema model interseismic velocities')") WRITE (*, "(' at GPS geodetic benchmarks.')") WRITE (*, "(' The method is to assume that the model errors are the same in')") WRITE (*, "(' the interseismic domain (of *.gps) and in the long-term domain')") WRITE (*, "(' (of g_[token].nko), and transfer these errors across files,')") WRITE (*, "(' for all benchmarks that were actually used in the final solution.')") WRITE (*, "(' Note that this is more accurate than a previous utility')") WRITE (*, "(' (now deprecated) which made the estimate by linear interpolation')") WRITE (*, "(' in space from v_interseismic_[token].out; the problem is that the')") WRITE (*, "(' fault locking/unlocking correction is very nonlinear in space')") WRITE (*, "(' near major faults.')") WRITE (*, "(' By Peter Bird, UCLA, 2012.01.19')") ! ---------------------------- read .GPS input file ----------------------------------- 10 WRITE (*, "(' ')") WRITE (*, "(' Enter name of geodetic input file (*.GPS): ')") READ (*, "(A)") GPS_filename OPEN (UNIT = 1, FILE = TRIM(GPS_filename), STATUS = "OLD", IOSTAT = ios, PAD = "YES") IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found (in current active folder).')") CALL Pause() GO TO 10 END IF READ (1, "(A)") GPS_titles WRITE (*, "(' ',A)") TRIM(GPS_titles) READ (1, "(A)") GPS_FORMAT WRITE (*, "(' ',A)") TRIM(GPS_FORMAT) READ (1, "(A)") GPS_headers WRITE (*, "(' ',A)") TRIM(GPS_headers) GPS_count = 0 counting: DO READ (1, *, IOSTAT = ios) E_lon_deg, N_lat_deg IF (ios == 0) THEN GPS_count = GPS_count + 1 ELSE EXIT counting END IF END DO counting WRITE (*, "(' contains ',I6,' benchmarks.')") GPS_count ALLOCATE ( GPS_E_lon_deg(GPS_count) ) ALLOCATE ( GPS_N_lat_deg(GPS_count) ) ALLOCATE ( GPS_v_E_mmpa(GPS_count) ) ! REAL ALLOCATE ( GPS_v_N_mmpa(GPS_count) ) ! REAL ALLOCATE ( GPS_site(GPS_count) ) ! C*20 CLOSE (1) OPEN (UNIT = 1, FILE = TRIM(GPS_filename), STATUS = "OLD", IOSTAT = ios, PAD = "YES") READ (1, "(A)") GPS_titles READ (1, "(A)") GPS_FORMAT READ (1, "(A)") GPS_titles DO i = 1, GPS_count READ (1, GPS_FORMAT) E_lon_deg, N_lat_deg, v_E_mmpa, v_N_mmpa, v_E_sigma, v_N_sigma, correlation, reference_frame, site GPS_E_lon_deg(i) = E_lon_deg GPS_N_lat_deg(i) = N_lat_deg GPS_v_E_mmpa(i) = v_E_mmpa GPS_v_N_mmpa(i) = v_N_mmpa GPS_site(i) = TRIM(site) END DO ! ---------------------------- read .g_[token].nko output file ----------------------------------- 20 WRITE (*, "(' ')") WRITE (*, "(' Enter name-token for NeoKinema model (with no _ or . included): ')") READ (* , "(A)") model_name_token g_token_NKO_filename = "g_" // TRIM(model_name_token) // ".nko" WRITE (*, "(' Looking for NeoKinema output file ',A, '...')") TRIM(g_token_NKO_filename) OPEN (UNIT = 2, FILE = TRIM(g_token_NKO_filename), STATUS = "OLD", IOSTAT = ios, PAD = "YES") IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found (in current active folder).')") CALL Pause() GO TO 20 END IF READ (2, "(A)") g_token_NKO_titles WRITE (*, "(' ',A)") TRIM(g_token_NKO_titles) READ (2, "(A)") g_token_NKO_FORMAT WRITE (*, "(' ',A)") TRIM(g_token_NKO_FORMAT) READ (2, "(A)") g_token_NKO_headers WRITE (*, "(' ',A)") TRIM(g_token_NKO_headers) g_token_NKO_count = 0 counting_again: DO READ (2, *, IOSTAT = ios) E_lon_deg, N_lat_deg IF (ios == 0) THEN g_token_NKO_count = g_token_NKO_count + 1 ELSE EXIT counting_again END IF END DO counting_again WRITE (*, "(' contains ',I6,' benchmarks.')") g_token_NKO_count ALLOCATE ( g_token_NKO_E_lon_deg(g_token_NKO_count) ) ALLOCATE ( g_token_NKO_N_lat_deg(g_token_NKO_count) ) ALLOCATE ( g_token_NKO_E_error_mmpa(g_token_NKO_count) ) ! REAL ALLOCATE ( g_token_NKO_N_error_mmpa(g_token_NKO_count) ) ! REAL ALLOCATE ( g_token_NKO_site(g_token_NKO_count) ) ! C*20 CLOSE (2) OPEN (UNIT = 2, FILE = TRIM(g_token_NKO_filename), STATUS = "OLD", IOSTAT = ios, PAD = "YES") READ (2, "(A)") g_token_NKO_titles READ (2, "(A)") g_token_NKO_FORMAT READ (2, "(A)") g_token_NKO_titles DO i = 1, g_token_NKO_count READ (2, g_token_NKO_FORMAT) E_lon_deg, N_lat_deg, v_E_mmpa, v_N_mmpa, v_E_sigma, v_N_sigma, correlation, reference_frame, site, & E_error_mmpa, N_error_mmpa g_token_NKO_E_lon_deg(i) = E_lon_deg g_token_NKO_N_lat_deg(i) = N_lat_deg g_token_NKO_E_error_mmpa(i) = E_error_mmpa g_token_NKO_N_error_mmpa(i) = N_error_mmpa g_token_NKO_site(i) = TRIM(site) END DO ! ---------------------------- combine files and write results ----------------------------------- g_interseismic_token_NKO_filename = "g_interseismic_" // TRIM(model_name_token) // ".nko" OPEN (UNIT = 3, FILE = TRIM(g_interseismic_token_NKO_filename)) ! unconditional OPEN; overwrites existing. WRITE (3, "('Interseismic benchmark velocities of NeoKinema model ',A)") TRIM(model_name_token) WRITE (3, "('from information in files ',A,' and ',A)") TRIM(GPS_filename), TRIM(g_token_NKO_filename) WRITE (3, "(' E_lon_deg N_lat_deg GPS_E_mm/a GPS_N_mm/a NK_E_mm/a NK_N_mm/a Error_E_mm/a Error_N_mm/a Site')") DO i = 1, g_token_NKO_count ! for every benchmark reported in g_[token].nko, match = 0 ! just initializing; should not be 0 at end! searching: DO j = 1, GPS_count ! try every benchmark in GPS_filename IF ((MOD(GPS_E_lon_deg(j) - g_token_NKO_E_lon_deg(i) + 720.0D0, 360.0D0) == 0.0D0).AND. & & (GPS_N_lat_deg(j) == g_token_NKO_N_lat_deg(i))) THEN match = j EXIT searching END IF END DO searching IF (match == 0) THEN ! failure! WRITE (*, "(' ERROR: Could not find benchmark #',I6,' at ',F10.2,', ',F10.2)") & & i, g_token_NKO_E_lon_deg(i), g_token_NKO_N_lat_deg(i) WRITE (*, "(' in the input GPS file ',A)") TRIM(GPS_filename) CALL Pause() END IF !Continue, assuming a match was found: interseismic_E_mmpa = GPS_v_E_mmpa(match) + g_token_NKO_E_error_mmpa(i) interseismic_N_mmpa = GPS_v_N_mmpa(match) + g_token_NKO_N_error_mmpa(i) WRITE (3, "(6F11.2,2F13.2,3X,A)") g_token_NKO_E_lon_deg(i), g_token_NKO_N_lat_deg(i), & & GPS_v_E_mmpa(match), GPS_v_N_mmpa(match), & & interseismic_E_mmpa, interseismic_N_mmpa, & & g_token_NKO_E_error_mmpa(i), g_token_NKO_N_error_mmpa(i), & & TRIM(GPS_site(match)) END DO WRITE (*, "(' ')") WRITE (*, "(' ')") WRITE (*, "(' ')") WRITE (*, "(' Job completed.')") CALL Pause() CLOSE (1) ! GPS_filename (OLD) CLOSE (2) ! g_token_NK0_filename (OLD) CLOSE (3) ! g_interseismic_token_NKO_filename (NEW) CONTAINS SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM NK_Interseismic_Vs_At_Benchmarks