PROGRAM SubDigitise ! Reads a .dig file in (lon,lat) format. ! Subdivides any line segment longer than a user-supplied limit: ! * if endpoints are "same" longitude, result is a meridian; ! * if endpoints are "same" latitude, result is a parallel; ! * in other cases, result is an arc of a great circle. ! by Peter Bird, UCLA, October 1997; updated July 2020. ! (c) Copyright 1997, 2020 by Peter Bird and the Regents of the University of California. USE DSphere ! provided by P. Bird as file DSphere.f90 CHARACTER*80 :: in_file_name, line, out_file_name INTEGER :: i, ios, n_added, n_inside, n_points LOGICAL :: got_point, got_previous, includes_TF REAL*8 :: arc_km, azimuth_radians, bottom_km, frac, & & lat, lat1, lat2, limit_km, lon, lon1, lon2, & & radius_km REAL*8, PARAMETER :: km_per_degree = 111.1D0 REAL*8, DIMENSION(3) :: uvec, uvec1, uvec2 !------------------------------------------------------ WRITE (*,"(/' PROGRAM SubDigitise' & & //' Reads a .dig file in (lon,lat) format.' & & /' Subdivides any line segment longer than a user-supplied limit:' & & /' * if endpoints are same longitude, result is a meridian;' & & /' * if endpoints are same latitude, result is a parallel;' & & /' * in other cases, result is an arc of a great circle.' & & /' by Peter Bird, UCLA, October 1997')") ! - - - - - - - - - - - - -- - - - - - - 10 WRITE (*,"(/' Enter [drive:\][path\]name of a .dig file in (lon,lat) format:')") READ (*,"(A)") in_file_name OPEN (UNIT = 1, FILE = in_file_name, STATUS = 'OLD', PAD = 'YES', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' File not found in current directory; add path?')") CLOSE (1, IOSTAT = ios) GOTO 10 END IF n_points = 0 DO READ (1, "(A)", IOSTAT = ios) line IF (ios == -1) EXIT ! EOF READ (line, *, IOSTAT = ios) lon1, lat1 IF (ios == 0) n_points = n_points + 1 ! read a number END DO CLOSE (1) OPEN (UNIT = 1, FILE = in_file_name, STATUS = 'OLD', PAD = 'YES', IOSTAT = ios) WRITE (*, "(' ',I10,' points.')") n_points IF (n_points == 0) GOTO 10 ! - - - - - - - - - - - - - - - - WRITE (*,"(/' Enter radius of the planet, in km: '\)") READ (*, *) radius_km 20 WRITE (*,"(/' Enter longest permissable segment, in km: '\)") READ (*, *) limit_km ! - - - - - - - - - - - - - - - -- - 30 WRITE (*, "(/' Enter [drive:\][path\]name for output file: ')") READ (*, "(A)") out_file_name OPEN (UNIT = 2, FILE = out_file_name, STATUS = 'NEW', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' File already exists. Choose another name.')") CLOSE (2, IOSTAT = ios) GOTO 30 END IF ! - - - - - - - - - - - - - - - - - - n_added = 0 got_previous = .FALSE. DO ! main loop of program; reading lines READ (1, "(A)", IOSTAT = ios) line IF (ios == -1) EXIT ! EOF !Test for 'T', 't', 'F', 'f' in the line; !when one of these begins a word, list-directed input !processing treats it as .TRUE. or .FALSE. and converts !this value to 1.0 or 0.0, causing mis-alignment of !data and/or ficticious data points. For example, !the title line "TX" was read as .TRUE., converted to REAL !1.0, and assigned to the longitude, and then the !first number of the next line was read for the latitude! !Another time, the title line " 129 Tonga-1" was !interpreted as (+129.00E, +1.00N), causing an unwanted !great-circle arc to that point! includes_TF = (SCAN(line,'T') > 0).OR.(SCAN(line,'t') > 0).OR. & & (SCAN(line,'F') > 0).OR.(SCAN(line,'f') > 0) IF (includes_TF) THEN got_previous = .FALSE. WRITE (2, "(A)") TRIM(line) ELSE ! no T or F; safe to try reading for two numbers READ (line, * , IOSTAT = ios) lon2, lat2 got_point = (ios == 0) IF (got_point) THEN CALL DLonLat_2_Uvec (lon2, lat2, uvec2) IF (got_previous) THEN ! got a segment to check! !==================HEART OF PROGRAM======================== arc_km = DArc (uvec1, uvec2) * radius_km IF (arc_km > limit_km) THEN ! arc must be sub-divided n_inside = MAX(1,(NINT(arc_km / limit_km) - 1)) DO i = 1, n_inside n_added = n_added + 1 frac = (1.0D0 * i) / (n_inside + 1.0D0) !GPB IF (ABS(lon2 - lon1) < tolerance_degrees) THEN ! create a meridian at lon1 lon = (lon1 + lon2) / 2.0D0 lat = lat1 + frac * (lat2 - lat1) ELSE IF (ABS(lat2 - lat1) < tolerance_degrees) THEN ! create a parallel at lat1 lon = lon1 + frac * (lon2 - lon1) lat = (lat1 + lat2) / 2.0D0 ELSE ! create a great-circle arc CALL GreatCircle_Point (uvec1, uvec2, frac, & ! inputs & uvec, azimuth_radians) ! outputs CALL DUvec_2_LonLat (uvec, lon, lat) END IF ! meridian, parallel, or arc WRITE (2, "(1X, SP, ES12.5, ',', ES12.5)") lon, lat END DO ! inside points END IF ! arc must be subdivided (don't write end-point!) !========================================================== END IF ! got_previous WRITE (2, "(1X, SP, ES12.5, ',', ES12.5)") lon2, lat2 ! prepare to loop got_previous = .TRUE. lon1 = lon2 lat1 = lat2 uvec1(1:3) = uvec2(1:3) ELSE ! a termination or title line got_previous = .FALSE. WRITE (2, "(A)") TRIM(line) END IF ! numbers or other type of line END IF ! includes_TF, or not END DO ! - - - - - - - - - - - - - - - - - - WRITE (*,"(/' Job done. ', I10, ' points added.')") n_added CLOSE (1) CLOSE (2) CONTAINS SUBROUTINE GreatCircle_Point (from_uvec, to_uvec, s, & ! inputs & point_uvec, azimuth_radians) ! outputs ! Finds a point within a lesser arc of a great circle from ! "from_uvec" to "to_uvec" identified by dimensionless variable ! "s", which is 0.00 at "from_uvec" and 1.00 at "to_uvec. ! Also returns "azimuth_radians" which is the azimuth of travel ! (with increasing s) at that point. ! Useful for placing tick marks along a great-circle fault, for example. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec REAL*8, INTENT(IN) :: s REAL*8, DIMENSION(3), INTENT(OUT) :: point_uvec REAL*8, INTENT(OUT) :: azimuth_radians REAL*8 :: radians_to_s, radians_to_1, start_azimuth REAL*8, DIMENSION(3) :: omega_uvec, t_uvec start_azimuth = DRelative_Compass(from_uvec, to_uvec) radians_to_1 = DArc(from_uvec, to_uvec) radians_to_s = s * radians_to_1 !save to_uvec, in case CALLer modifies it by repeating actual argument t_uvec = to_uvec CALL DTurn_To (start_azimuth, from_uvec, radians_to_s, & ! inputs & omega_uvec, point_uvec) ! point_uvec might replace to_uvec! IF (ABS(1.00D0 - s) >= 0.5D0) THEN azimuth_radians = DRelative_Compass(point_uvec, t_uvec) ELSE azimuth_radians = Pi + DRelative_Compass(point_uvec, from_uvec) END IF END SUBROUTINE GreatCircle_Point SUBROUTINE Traceback () ! The sole function of this unit is to cause a traceable error, ! so that the route into the unit that called it is also traced. ! This unit is a good place to put a breakpoint while debugging! ! The intentional error must NOT be detected during compilation, ! but MUST cause a traceable error at run-time. ! If this routine has any error detected during compilation, ! then change its code to cause a different intentional error. IMPLICIT NONE INTEGER :: i REAL, DIMENSION(3) :: y WRITE (*,"(' -----------------------------------------------------')") WRITE (*,"(' Traceback was called to execute an intentional error:')") WRITE (*,"(' An array subscript will be intentionally out-of-range.')") DO i = 1, 4 y(i) = 1. * i END DO STOP ' ' END SUBROUTINE Traceback END PROGRAM SubDigitise