PROGRAM Extract_Regional_GRD !Utility program to extract a smaller regional GRD file !from a larger (possibly global) GRD file, !and write it as a new file. !If the original file has two lines after the data that !give the integrated seismicity, then the extracted file !will also have two similar lines, but with a smaller integral. !By Peter Bird, UCLA, 2014.04.02 IMPLICIT NONE CHARACTER*13 :: c13 CHARACTER*80 :: big_GRD_file, small_GRD_file, line1_suffix, line2_suffix CHARACTER*132 :: line1, line2 DOUBLE PRECISION :: EQs_per_s, EQs_per_century INTEGER :: choice, i, ios, iSource, j, jSource, nColsBIG, nColsSmall, nRowsBIG, nRowsSmall LOGICAL :: got_integral REAL, PARAMETER :: Earth_R_meters = 6371000. REAL, PARAMETER :: radians_per_degree = 0.0174533 REAL :: d_lat, d_lon, dx, dy, & & lat, lon, & & lat_max_BIG, lat_min_BIG, lon_max_BIG, lon_min_BIG, & & lat_max_small, lat_min_small, lon_max_small, lon_min_small REAL, DIMENSION(:, :), ALLOCATABLE :: BIG_GRD, small_GRD !Header lines WRITE (*, "(' PROGRAM Extract_Regional_GRD')") WRITE (*, *) WRITE (*, "(' Utility program to extract a smaller regional GRD file')") WRITE (*, "(' from a larger (possibly global) GRD file,')") WRITE (*, "(' and write it as a new file.')") WRITE (*, "(' If the original file has two lines after the data that')") WRITE (*, "(' give the integrated seismicity, then the extracted file')") WRITE (*, "(' will also have two similar lines, but with a smaller integral.')") WRITE (*, "(' By Peter Bird, UCLA, 2014.04.02')") WRITE (*, *) !Get name of large GRD file and OPEN. big_GRD_file = ' ' 10 CALL Prompt_for_String('[Drive:][\path\]filename of existing large GRD file:', big_GRD_file, big_GRD_file) OPEN (UNIT = 1, FILE = TRIM(big_GRD_file), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found in current folder (or along typed path).')") WRITE (*, "(' Please try again...')") CALL Pause() GO TO 10 END IF !Read and memorize large GRD file: READ (1, *) lon_min_BIG, d_lon, lon_max_BIG READ (1, *) lat_min_BIG, d_lat, lat_max_BIG nRowsBIG = 1 + NINT((lat_max_BIG - lat_min_BIG) / d_lat) nColsBIG = 1 + NINT((lon_max_BIG - lon_min_BIG) / d_lon) ALLOCATE ( BIG_GRD(nRowsBIG, nColsBIG) ) READ (1, *) ((BIG_GRD(i, j), j = 1, nColsBIG), i = 1, nRowsBIG) WRITE (*, "(' Large GRD file has been read and memorized.')") !Test for presence of 2 extra lines: READ (1, "(A)", IOSTAT = ios) c13 ! "Area integral"? got_integral = (ios == 0).AND.(c13 == "Area integral") IF (got_integral) THEN READ (1, *) EQs_per_s EQs_per_century = EQs_per_s * 3.155769E9 CLOSE (1) OPEN (UNIT = 1, FILE = TRIM(big_GRD_file), STATUS = "OLD", IOSTAT = ios) READ (1, "(A)") line1 READ (1, "(A)") line2 line1_suffix = line1(60:132) line2_suffix = line2(60:132) ELSE EQs_per_s = 0.0D0 EQs_per_century = 0.0D0 line1_suffix = ' From extract_regional_GRD.exe' line2_suffix = ' by Peter Bird, UCLA, 2014.04' END IF CLOSE (1) !Ask user to define rectangle, and ALLOCATE small GRD: WRITE (*, *) WRITE (*, "(' Extreme coordinate values in large GRD: ')") WRITE (*, "(' ',12X,F12.4)") lat_max_BIG WRITE (*, "(' ',F12.4,12X,F12.4)") lon_min_BIG, lon_max_BIG WRITE (*, "(' ',12X,F12.4)") lat_min_BIG WRITE (*, *) WRITE (*, "(' Implied outer-limits of cells (for a seismicity GRD): ')") WRITE (*, "(' ',12X,F12.4)") lat_max_BIG + 0.5 * d_lat WRITE (*, "(' ',F12.4,12X,F12.4)") lon_min_BIG - 0.5 * d_lon, lon_max_BIG + 0.5 * d_lon WRITE (*, "(' ',12X,F12.4)") lat_min_BIG - 0.5 * d_lat WRITE (*, *) WRITE (*, "(' How do you wish to define the limits of the new, smaller GRD?')") WRITE (*, "(' (1) By extreme coordinates of grid points?')") WRITE (*, "(' (2) By outer limits of implied seismicity cells?')") 20 choice = 1 CALL Prompt_for_Integer(' Enter 1 or 2:', choice, choice) IF ((choice < 1).OR.(choice > 2)) THEN WRITE (*, "(' Illegal choice. Please try again.')") CALL Pause() GO TO 20 END IF IF (choice == 1) THEN lon_min_small = lon_min_BIG CALL Prompt_for_Real('Minimum longitude (+ for E): ', lon_min_small, lon_min_small) lon_max_small = lon_max_BIG CALL Prompt_for_Real('Maximum longitude (+ for E): ', lon_max_small, lon_max_small) lat_min_small = lat_min_BIG CALL Prompt_for_Real('Minimum latitude (+ for N): ', lat_min_small, lat_min_small) lat_max_small = lat_max_BIG CALL Prompt_for_Real('Maximum latitude (+ for N): ', lat_max_small, lat_max_small) ELSE IF (choice == 2) THEN lon_min_small = lon_min_BIG - 0.5 * d_lon CALL Prompt_for_Real('Minimum longitude (+ for E): ', lon_min_small, lon_min_small) lon_min_small = lon_min_small + 0.5 * d_lon lon_max_small = lon_max_BIG + 0.5 * d_lon CALL Prompt_for_Real('Maximum longitude (+ for E): ', lon_max_small, lon_max_small) lon_max_small = lon_max_small - 0.5 * d_lon lat_min_small = lat_min_BIG - 0.5 * d_lat CALL Prompt_for_Real('Minimum latitude (+ for N): ', lat_min_small, lat_min_small) lat_min_small = lat_min_small + 0.5 * d_lat lat_max_small = lat_max_BIG + 0.5 * d_lat CALL Prompt_for_Real('Maximum latitude (+ for N): ', lat_max_small, lat_max_small) lat_max_small = lat_max_small - 0.5 * d_lat END IF nRowsSmall = 1 + NINT((lat_max_small - lat_min_small) / d_lat) nColsSmall = 1 + NINT((lon_max_small - lon_min_small) / d_lon) ALLOCATE ( small_GRD(nRowsSmall, nColsSmall) ) !fill up the small GRD array: DO i = 1, nRowsSmall lat = lat_max_small - (i - 1) * d_lat iSource = 1 + NINT((lat_max_BIG - lat) / d_lat) IF ((iSource < 1).OR.(iSource > nRowsBIG)) THEN WRITE (*, "(' ERROR: Small-grid latitude ',F10.2,' requires accessing non-existant row ',I6,' of large grid.')") lat, iSource CALL Pause() STOP END IF DO j = 1, nColsSmall lon = lon_min_small + (j - 1) * d_lon jSource = 1 + NINT((lon - lon_min_BIG) / d_lon) IF ((jSource < 1).OR.(jSource > nColsBIG)) THEN WRITE (*, "(' ERROR: Small-grid longitude ',F10.2,' requires accessing non-existant column ',I6,' of large grid.')") lon, jSource CALL Pause() STOP END IF !-------------------------------------------- small_GRD(i, j) = BIG_GRD(iSource, jSource) !-------------------------------------------- END DO ! j = 1, nColsSmall END DO ! i = 1, nRowsSmall WRITE (*, *) WRITE (*, "(' Grid values transferred to new, small GRD.')") !Get filename for output file, OPEN, and WRITE: small_GRD_file = "extract_regional_GRD.grd" CALL Prompt_for_String('[Drive:][\path\]filename for new, small GRD file:', small_GRD_file, small_GRD_file) OPEN (UNIT = 2, FILE = TRIM(small_GRD_file)) ! unconditional OPEN, overwrites any existing file WRITE (2, "(3F11.5,' = lon_min, d_lon, lon_max',A)") lon_min_small, d_lon, lon_max_small, TRIM(line1_suffix) WRITE (2, "(3F11.5,' = lat_min, d_lat, lat_max',A)") lat_min_small, d_lat, lat_max_small, TRIM(line2_suffix) WRITE (2, "(10ES10.2)") ((small_GRD(i, j), j = 1, nColsSmall), i = 1, nRowsSmall) ! normal format, for REAL and DP data !WRITE (2, "(50I2)") ((NINT(small_GRD(i, j)), j = 1, nColsSmall), i = 1, nRowsSmall) ! variant FORMAT, for tectonic zones IF (got_integral) THEN WRITE (2, "('Area integral of seismicity =')") EQs_per_s = 0.0D0 ! just initializing before integration dy = Earth_R_meters * d_lat * radians_per_degree DO i = 1, nRowsSmall lat = lat_max_small - (i - 1) * d_lat dx = Earth_R_meters * d_lon * radians_per_degree * COS(lat * radians_per_degree) DO j = 1, nColsSmall EQs_per_s = EQs_per_s + dx * dy * small_GRD(i, j) END DO END DO EQs_per_century = EQs_per_s * 3.15576E9 WRITE (2, "(ES11.4,' earthquakes/s = ',ES11.4,' earthquakes/century.')") EQs_per_s, EQs_per_century END IF CLOSE (2) WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS CHARACTER*10 FUNCTION ASCII10(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, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL :: x_log DOUBLE PRECISION :: y IF (x == 0.0) THEN ASCII10=' 0' RETURN ELSE IF (x > 0.0) THEN punt = (x >= 999999999.5).OR.(x < 0.0000100) ELSE ! x < 0.0 punt = (x <= -99999999.5).OR.(x > -0.000100) 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 ASCII10 = temp10 ELSE ! can represent without exponential notation x_log = LOG10(ABS(x)) zeros = Int_Below(x_log) - 3 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0) 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 ASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION ASCII10 INTEGER FUNCTION Int_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL, INTENT(IN) :: x INTEGER :: i REAL :: y i = INT(x) IF (x >= 0.) THEN Int_Below = i ELSE ! x < 0. y = 1.*i IF (y <= x) THEN Int_Below = i ELSE ! most commonly Int_Below = i - 1 END IF END IF END FUNCTION Int_Below SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause SUBROUTINE Prompt_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 Prompt_for_Integer SUBROUTINE Prompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE Prompt_for_Logical SUBROUTINE Prompt_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, INTENT(IN) :: default REAL, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use ASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = ASCII10(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 Prompt_for_Real SUBROUTINE Prompt_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 Prompt_for_String END PROGRAM Extract_Regional_GRD