PROGRAM Extract_Rectangular_FEG

    !Reads a finite-element grid (.FEG) file in Shells or NeoKinema format,
    !and trims the grid near a user-specified latitude/longitude rectangle,
    !producing a smaller .FEG grid file.
    !All nodes and elements within the rectangle are kept, and renumbered.
    !Elements partially within the rectangle and their associated nodes
    !are also kept.  This produces an irregular grid outline, but the
    !irregular outline lies outside the specified rectangle.
    !Boundary nodes are renumbered but not moved; thus one can use
    !utility program complete_type2_bNKI to select boundary velocities
    !for these nodes based on a previous regional NeoKinema model.
    !By Peter Bird, UCLA, 2010.02.13.

    IMPLICIT NONE
    CHARACTER*80  :: old_FEG_pathfile = ' ', old_title_line = ' ', new_FEG_pathfile = ' ', new_title_line = ' '
    CHARACTER*132 :: old_FEG_file_line = ' '
    INTEGER :: blocks, count, i, inferred_fault_data, inferred_nodal_data, ios, j, k, &
             & n1, n2, n3, n4, new_numnod, new_numel, new_nfl, old_nfl, old_numel, old_numnod
    INTEGER, DIMENSION(:),   ALLOCATABLE :: new_node_index
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: old_elements, old_faults, new_elements, new_faults
    LOGICAL, DIMENSION(:),   ALLOCATABLE :: keep_this_element, keep_this_fault, keep_this_node, node_is_inside
    REAL*8 :: new_max_lat, new_max_lon, new_min_lat, new_min_lon
    REAL*8 :: old_max_lat, old_max_lon, old_min_lat, old_min_lon
    REAL*8,    DIMENSION(:),   ALLOCATABLE :: node_lat, node_lon, new_node_lat, new_node_lon
    REAL*8,    DIMENSION(:,:), ALLOCATABLE :: fault_data, nodal_data, new_nodal_data, new_fault_data

    WRITE (*, "(' PROGRAM Extract_Rectangular_FEG:')")
    WRITE (*, "(' ')")
    WRITE (*, "(' Reads a finite-element grid (.FEG) file in Shells or NeoKinema format,')")
    WRITE (*, "(' and trims the grid near a user-specified latitude/longitude rectangle,')")
    WRITE (*, "(' producing a smaller .FEG grid file.')")
    WRITE (*, "(' All nodes and elements within the rectangle are kept, and renumbered.')")
    WRITE (*, "(' Elements partially within the rectangle and their associated nodes')")
    WRITE (*, "(' are also kept.  This produces an irregular grid outline, but the')")
    WRITE (*, "(' irregular outline lies outside the specified rectangle.')")
    WRITE (*, "(' Boundary nodes are renumbered but not moved; thus one can use')")
    WRITE (*, "(' utility program complete_type2_bNKI to select boundary velocities')")
    WRITE (*, "(' for these nodes based on a previous regional NeoKinema model.')")
    WRITE (*, "(' By Peter Bird, UCLA, 2010.02.13.')")
    WRITE (*, "(' ')")

    !read and memorize existing .FEG file:
  1 WRITE (*, "(' Enter [Drive:][\path\]filename of existing large .FEG file: ')")
    READ (*, "(A)") old_FEG_pathfile
    OPEN (UNIT = 1, FILE = old_FEG_pathfile, STATUS = "OLD", PAD = "YES", IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: This .FEG file not found; please try again...')")
        CALL Pause()
        GO TO 1
    END IF
    READ (1, "(A)") old_title_line
    READ (1, *) old_numnod
    ALLOCATE ( node_lon(old_numnod) ); node_lon = 0.0D0
    ALLOCATE ( node_lat(old_numnod) ); node_lat = 0.0D0
    ALLOCATE ( nodal_data(6, old_numnod) ); nodal_data = 0.0D0 ! large enough for Shells-format
    DO i = 1, old_numnod
        READ (1, "(A)") old_FEG_file_line
            blocks = Scan_Line(old_FEG_file_line)
            inferred_nodal_data = blocks - 3
        IF (inferred_nodal_data > 0) THEN
            READ (old_FEG_file_line, *) j, node_lon(j), node_lat(j), (nodal_data(k, j), k = 1, inferred_nodal_data) 
        ELSE 
            READ (old_FEG_file_line, *) j, node_lon(j), node_lat(j)
        END IF
        ! PAD = "YES" means "0.0" will be inferred in shorter lines of NeoKinema-type .FEG files
    END DO ! i = 1, old_numnod; reading nodal data
    READ (1, *) old_numel
    ALLOCATE ( old_elements(3, old_numel) ); old_elements = 0
    DO i = 1, old_numel
        READ (1, *) j, (old_elements(k, j), k = 1, 3)
    END DO ! i = 1, old_numel; reading in element definitions
    READ (1, *) old_nfl
    IF (old_nfl > 0) THEN ! this would be a Shells-type .FEG file
        ALLOCATE ( old_faults(4, old_nfl) ); old_faults = 0
        ALLOCATE ( fault_data(3, old_nfl) ); fault_data = 0.0D0
        DO i = 1, old_nfl
            READ (1, "(A)") old_FEG_file_line
            blocks = Scan_Line(old_FEG_file_line)
            inferred_fault_data = blocks - 5
            IF (inferred_fault_data > 0) THEN
                READ (old_FEG_file_line, *) j, (old_faults(k, j), k = 1, 4), (fault_data(k, j), k = 1, inferred_fault_data)
            ELSE
                READ (old_FEG_file_line, *) j, (old_faults(k, j), k = 1, 4)
            END IF
        END DO
    END IF ! old_nfl > 0
    WRITE (*, "(' Existing .FEG file has been read and memorized.')")
    WRITE (*, "(' Inferred number of nodal data (after lon, lat) was ',I1)") inferred_nodal_data
    IF (old_nfl > 0) THEN
        WRITE (*, "(' Inferred number of fault data (after node list) was ',I1)") inferred_fault_data
    END IF
    old_max_lat = node_lat(1)
    old_min_lat = node_lat(1)
    old_max_lon = node_lon(1)
    old_min_lon = node_lon(1)
    DO i = 2, old_numnod
        old_max_lat = MAX(old_max_lat, node_lat(i))
        old_min_lat = MIN(old_min_lat, node_lat(i))
        old_max_lon = MAX(old_max_lon, node_lon(i))
        old_min_lon = MIN(old_min_lon, node_lon(i))
    END DO
    WRITE (*, "(' Node longitudes range from ',F8.3,' up to ',F8.3)") old_min_lon, old_max_lon
    WRITE (*, "(' Node latitudes  range from ',F8.3,' up to ',F8.3)") old_min_lat, old_max_lat

    !define the latitude/longitude rectangle:
    WRITE (*, *)
 10 WRITE (*, "(' Enter minimum (Westernmost) longitude for selection rectangle: ')")
    READ (*, *) new_min_lon
    WRITE (*, "(' Enter maximum (Easternmost) longitude for selection rectangle: ')")
    READ (*, *) new_max_lon
    IF (new_max_lon <= new_min_lon) THEN
        WRITE (*, "(' ERROR: Illogical.  Try again...')")
        CALL Pause()
        GO TO 10
    END IF
 11 WRITE (*, "(' Enter minimum (Southernmost) latitude for selection rectangle: ')")
    READ (*, *) new_min_lat
    WRITE (*, "(' Enter maximum (Northernmost) latitude for selection rectangle: ')")
    READ (*, *) new_max_lat
    IF (new_max_lat <= new_min_lat) THEN
        WRITE (*, "(' ERROR: Illogical.  Try again...')")
        CALL Pause()
        GO TO 11
    END IF

    !perform selection operations
    ALLOCATE ( node_is_inside(old_numnod) ); node_is_inside = .FALSE.
    count = 0
    DO i = 1, old_numnod
        IF ((node_lon(i) >= new_min_lon).AND.(node_lon(i) <= new_max_lon).AND. &
          & (node_lat(i) >= new_min_lat).AND.(node_lat(i) <= new_max_lat)) THEN
            node_is_inside(i) = .TRUE.
            count = count + 1
        END IF
    END DO
    WRITE (*, *)
    WRITE (*, "(' ',I8,' nodes are inside the selection rectangle.')") count
    IF (count == 0) THEN
        WRITE (*, "(' No output .FEG file will be produced.')")
        CALL Pause()
        STOP
    END IF
    ALLOCATE ( keep_this_element(old_numel) ); keep_this_element = .FALSE.
    new_numel = 0
    DO i = 1, old_numel
        n1 = old_elements(1, i)
        n2 = old_elements(2, i)
        n3 = old_elements(3, i)
        IF (node_is_inside(n1).OR.node_is_inside(n2).OR.node_is_inside(n3)) THEN
            keep_this_element(i) = .TRUE.
            new_numel = new_numel + 1
        END IF
    END DO
    WRITE (*, "(' ',I8,' continuum elements will be retained.')") new_numel
    IF (old_nfl > 0) THEN
        ALLOCATE ( keep_this_fault(old_nfl) ); keep_this_fault = .FALSE.
        new_nfl = 0
        DO i = 1, old_nfl
            n1 = old_faults(1, i)
            n2 = old_faults(2, i)
            n3 = old_faults(3, i)
            n4 = old_faults(4, i)
            IF (node_is_inside(n1).OR.node_is_inside(n2).OR.node_is_inside(n3).OR.node_is_inside(n4)) THEN
                keep_this_fault(i) = .TRUE.
                new_nfl = new_nfl + 1
            END IF
        END DO
        WRITE (*, "(' ',I8,' fault     elements will be retained.')") count
    ELSE
        new_nfl = 0
    END IF
    ALLOCATE ( keep_this_node(old_numnod) ); keep_this_node = .FALSE.
    DO i = 1, old_numel
        IF (keep_this_element(i)) THEN
            DO j = 1, 3
                keep_this_node(old_elements(j, i)) = .TRUE.
            END DO
        END IF
    END DO
    IF (new_nfl > 0) THEN
        DO i = 1, old_nfl
            IF (keep_this_fault(i)) THEN
                DO j = 1, 4
                    keep_this_node(old_faults(j, i)) = .TRUE.
                END DO
            END IF
        END DO
    END IF
    new_numnod = 0
    DO i = 1, old_numnod
        IF (keep_this_node(i)) new_numnod = new_numnod + 1
    END DO
    WRITE (*, "(' ',I8,' nodes (in & out)   will be retained.')") new_numnod
     
    !build (shorter) list of nodes in new grid, and keep track of # translations:
    ALLOCATE ( new_node_index(old_numnod) ); new_node_index = 0
    ALLOCATE ( new_node_lon(new_numnod) ); new_node_lon = 0.0D0
    ALLOCATE ( new_node_lat(new_numnod) ); new_node_lat = 0.0D0
    IF (inferred_nodal_data > 0) THEN
        ALLOCATE ( new_nodal_data(inferred_nodal_data, new_numnod) ); new_nodal_data = 0.0D0
    END IF
    j = 0
    DO i = 1, old_numnod
        IF (keep_this_node(i)) THEN
            j = j + 1
            new_node_lon(j) = node_lon(i)
            new_node_lat(j) = node_lat(i)
            IF (inferred_nodal_data > 0) THEN
                DO k = 1, inferred_nodal_data
                    new_nodal_data(k, j) = nodal_data(k, i)
                END DO
            END IF
            new_node_index(i) = j
        END IF
    END DO

    !build (shorter) list of retained elements, translating the node numbers
    ALLOCATE ( new_elements(3, new_numel) ); new_elements = 0
    j = 0
    DO i = 1, old_numel
        IF (keep_this_element(i)) THEN
            j = j + 1
            DO k = 1, 3
                new_elements(k, j) = new_node_index(old_elements(k, i))
            END DO
        END IF
    END DO

    !build (shorter) list of retained faults, translating the node numbers
    IF (new_nfl > 0) THEN
        ALLOCATE ( new_faults(4, new_nfl) ); new_faults = 0
        ALLOCATE ( new_fault_data(inferred_fault_data, new_nfl) ); new_fault_data = 0.0D0
        j = 0
        DO i = 1, old_nfl
            IF (keep_this_fault(i)) THEN
                j = j + 1
                DO k = 1, 4
                    new_faults(k, j) = new_node_index(old_faults(k, i))
                END DO
                DO k = 1, inferred_fault_data
                    new_fault_data(k, j) = fault_data(k, i)
                END DO
            END IF
        END DO
    END IF

    !output the new, smaller .FEG file
    WRITE (*, *)
 91 WRITE (*, "(' Enter name for new, smaller .FEG file: ')")
    READ (*, "(A)") new_FEG_pathfile
    IF (TRIM(new_FEG_pathfile) == TRIM(old_FEG_pathfile)) THEN
        WRITE (*, "(' ERROR: You must specify a NEW file name.  Try again...')") 
        CALL Pause()
        GO TO 91
    END IF
    OPEN (UNIT = 2, FILE = TRIM(new_FEG_pathfile)) ! unconditional OPEN; overwrites any existing file!
    WRITE (*, "(' Old title line was:'/' ',A)") TRIM(old_title_line)
    WRITE (*, "(' Enter desired NEW title line: ')")
    READ (*, "(A)") new_title_line
    WRITE (2, "(A)") TRIM(new_title_line)
    WRITE (2, "(I8,I8,'       0 1000000       F (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)')") new_numnod, new_numnod
    DO i = 1, new_numnod
        IF (inferred_nodal_data > 0) THEN
            WRITE (2, "(I6,F11.5,F11.5,6ES10.2)") i, new_node_lon(i), new_node_lat(i), (new_nodal_data(k, i), k = 1, inferred_nodal_data)
        ELSE
            WRITE (2, "(I6,F11.5,F11.5)") i, new_node_lon(i), new_node_lat(i)
        END IF
    END DO    
    WRITE (2, "(I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)')") new_numel
    DO i = 1, new_numel
        WRITE (2, "(I6,3I6)") i, (new_elements(k, i), k = 1, 3)
    END DO
    WRITE (2, "(I10,' (NFL =  NUMBER OF CURVILINEAR FAULT ELEMENTS)')") new_nfl
    IF (new_nfl > 0) THEN
        DO i = 1, new_nfl
            WRITE (2, "(I5,4I6,2F6.1,ES10.2)") i, (new_faults(k, i), k = 1, 4), (new_fault_data(k, i), k = 1, inferred_fault_data)
        END DO
    END IF
    CLOSE(2)

    WRITE (*, *)
    WRITE (*, "(' Job completed.')")
    WRITE (*, "(' Don''t forget to use OrbNumber to renumber this new .FEG file.')")
    CALL Pause()

CONTAINS

    SUBROUTINE Pause()
       IMPLICIT NONE
       WRITE (*,"(' Press [Enter]...'\)")
       READ  (*,*)
    END SUBROUTINE Pause

    INTEGER FUNCTION Scan_Line(line132)
        !Determine how many non-empty blocks of contiguous characters there are in line132.
        CHARACTER*132, INTENT(IN) :: line132
        INTEGER :: count, i
        LOGICAL :: was_blank, now_blank
        count = 0
        was_blank = .TRUE.
        DO i = 1, 132
            now_blank = (line132(i:i) == ' ')
            IF (was_blank.AND.(.NOT.now_blank)) count = count + 1
            was_blank = now_blank ! initializing for next time through this loop
        END DO
        Scan_Line = count
    END FUNCTION Scan_Line

END PROGRAM Extract_Rectangular_FEG