!
!
! Subroutine Read_Grid
!
!
SUBROUTINE Read_Grid
! 
! READs in an .FEG grid file. If errors occur, display message-box & close file.
! If READ succeeds, Global LOGICAL :: gridLoaded will be set to .TRUE. Otherwise, it will remain .FALSE. .
!
USE dflib
USE Global

IMPLICIT NONE
CHARACTER(len=10) ierr_c10, i_c10, j_c10
CHARACTER(len=80) longer_line, shorter_line
CHARACTER(len=132) input_record
CHARACTER(len=512) msg0, msg1
INTEGER(4) iRet, ierr, ierr2
INTEGER(4) i, j, LRi, n1, n2, n3, n4, no
REAL(4) lon, lat, elev, q, hc,hm, dips1,dips2, offst, density_anomaly, cooling_curvature
REAL(8) aSum, bSum, cSum, rt1, rt2, rt3

msg0 = ' 'C
msg1 = ' 'C
OPEN(UNIT = 100, FILE = feg_inp, STATUS = "OLD", ACTION = "READ", IOSTAT = ierr, PAD = "YES")
!NOTE: PAD = "YES" allows reading either old OrbData/Shells .feg, or newer OrbData5/Shells .feg
!      It also facilitates reading nodal lines which MAY have 4 data, or none.
!      It also facilitates reading element lines, which MAY have 3 data, or none.
IF (ierr /= 0) then
    WRITE (ierr_c10, "(I10)") ierr
    ierr_c10 = ADJUSTL(ierr_c10)
	msg1 = ' During OPEN, IOSTAT = ' // TRIM(ierr_c10) // ' 'C
    CALL error1(msg1)
	RETURN
END IF
READ(100, '(A80)') old_FEG_title_line
new_FEG_title_line = old_FEG_title_line ! (unless edited later)

! process node-count(s):

READ (100, *, IOSTAT = ierr) numNod, nRealN, nFakeN, n1000, brief
IF (ierr /= 0) THEN
    CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
    WRITE (ierr_c10, "(I10)") ierr
    ierr_c10 = ADJUSTL(ierr_c10)
	msg1 = ' At node-header, IOSTAT = ' // TRIM(ierr_c10) // ' 'C
    CALL error1(msg1)
	RETURN
END IF
!
! allocate nodal memory (if necessary) & initialize nodal memory
!
IF (.NOT.ALLOCATED(nodeABG)) THEN
    ALLOCATE( nodeABG(3, mxnode) )
    ALLOCATE( eqcm(6, mxnode) )
    ALLOCATE( NMemo(mxnode) )
END IF
nodeABG = 0.0 ! whole array
eqcm = 0.0    ! whole array
NMemo = 0     ! whole array
!
! READ nodes & their data(?)
!
aSum = 0.0D0; bSum = 0.0D0; cSum = 0.0D0 ! initialize sum of all node uvec's
OrbData5 = .FALSE. ! global switch initialized as FALSE, unless test below succeeds:
nodeInput: DO i = 1, numNod
    READ (100, "(A)", IOSTAT = ierr) input_record
    IF (ierr /= 0) THEN
        CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
        WRITE (i_c10, "(I10)") i
        i_c10 = ADJUSTL(i_c10)
        msg1 = " While trying to read node #" // TRIM(i_c10) // ", hit end-of-file!"C
        CALL error1(msg1)
	    RETURN
    END IF
    !Note: This indirect reading method prevents reading TWO records in the search for eqcm, and for density_anomaly!
    READ (input_record, *, IOSTAT = ierr) j, lon, lat, elev, q, hc, hm, density_anomaly, cooling_curvature
    IF (ierr /= 0) THEN ! perhaps the last 2 columns are not there?
        density_anomaly = 0.0
        cooling_curvature = 0.0
        READ (input_record, *, IOSTAT = ierr2) j, lon, lat, elev, q, hc, hm
        IF (ierr2 /= 0) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            i_c10 = ADJUSTL(i_c10)
            msg1 = " Not enough data for node #" // TRIM(i_c10) // ' 'C
            CALL error1(msg1)
            RETURN
        END IF
        IF (j /= i) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            WRITE (j_c10, "(I10)") j
            i_c10 = ADJUSTL(i_c10)
            j_c10 = ADJUSTL(j_c10)
            msg1 = " When expecting node #" // TRIM(i_c10) // ", encountered #" // TRIM(j_c10) // " instead."C
            CALL error1(msg1)
            RETURN
        END IF
    END IF
    lon = lon * .017453293 ! from degrees to radians
    lat = lat * .017453293
    nodeABG(1,i) = COS(lat) * COS(lon)
    nodeABG(2,i) = COS(lat) * SIN(lon)
    nodeABG(3,i) = SIN(lat)
    aSum = aSum + nodeABG(1,i)
    bSum = bSum + nodeABG(2,i)
    cSum = cSum + nodeABG(3,i)
    eqcm(1,i) = elev
    eqcm(2,i) = q
    eqcm(3,i) = hc
    eqcm(4,i) = hm
    eqcm(5,i) = density_anomaly   ! may be 0.0 in older files
    eqcm(6,i) = cooling_curvature ! may be 0.0 in older files
    OrbData5 = OrbData5 .OR. (density_anomaly /= 0.0) .OR. (cooling_curvature /= 0.0)
END DO nodeInput
!
!choose WindowPosition to display this .feg:
!
aSum = aSum / DBLE(numNod)
bSum = bSum / DBLE(numNod)
cSum = cSum / DBLE(numNod)
call Unitize(aSum, bSum, cSum)
nettempvec(1) = aSum
nettempvec(2) = bSum
nettempvec(3) = cSum
!
! process element-count:
!
READ (100, *, IOSTAT = ierr) numel
IF (ierr /= 0) THEN
    CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
    WRITE (ierr_c10, "(I10)") ierr
    ierr_c10 = ADJUSTL(ierr_c10)
	msg1 = ' At element-header, IOSTAT = ' // TRIM(ierr_c10) // ' 'C
    CALL error1(msg1)
	RETURN
END IF
!
! allocate element memory (if necessary) & initialize element memory:
!
IF (.NOT.ALLOCATED(nodes))        ALLOCATE ( nodes(3, mxel) )
IF (.NOT.ALLOCATED(EMemo))        ALLOCATE ( EMemo(mxel) )
IF (.NOT.ALLOCATED(element_data)) ALLOCATE ( element_data(3, mxel) )
nodes = 0 ! whole array
EMemo = 0
element_data = 0.0
IF (AppType == 1) THEN ! Shells mode; support v5.0+
    IF (.NOT. ALLOCATED(continuum_LRi)) ALLOCATE ( continuum_LRi(mxel) )
    continuum_LRi = 0 ! whole array
END IF
!
! READ elements & element data:
!
IF (appType == 3) THEN ! initializing of mu-memory required:
    lowest_mu = 1.0D0
    highest_mu = 0.0D0 ! (expecting both of these to be replaced)
END IF
eleInput: DO i = 1, numel

    READ (100, "(A)", IOSTAT = ierr) input_record
    IF (ierr /= 0) THEN
        CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
        WRITE (i_c10, "(I10)") i
        i_c10 = ADJUSTL(i_c10)
        msg1 = " While trying to read element #" // TRIM(i_c10) // ", hit end-of-file!"C
        CALL error1(msg1)
	    RETURN
    END IF
    IF (AppType == 1) THEN ! Shells mode; support v5.0+
        longer_line = input_record
        CALL Extract_LRi (longer_line, &     ! input
                        & LRi, shorter_line) ! output
        input_record = shorter_line
        READ (input_record, *, IOSTAT = ierr) j
        continuum_LRi(j) = LRi ! Save Lithospheric Rheology index (which will often be 0).
        !N.B. This block of code is intented to have no side-effects on older code that follows.
    END IF
    IF (appType == 3) then ! Restore3+ format, with per-element data ALLOWED (not mandatory).  Capture range of mu_element values.
        !First, try to capture element with 3 per-element data:
        READ (input_record, *, IOSTAT = ierr) j, n1, n2, n3, rt1, rt2, rt3
        IF (ierr == 0) THEN
           ! READ went well, and 3 per-element data were captured
            element_data(1, i) = rt1
            element_data(2, i) = rt2
            element_data(3, i) = rt3
            lowest_mu = MIN(lowest_mu, rt1, rt3)
            highest_mu = MAX(highest_mu, rt1, rt3)
        ELSE ! try to read a shorter record:
            !Try to read only 2 per-element data:
            READ (input_record, *, IOSTAT = ierr) j, n1, n2, n3, rt1, rt2
            IF (ierr == 0) THEN
                element_data(1, i) = rt1
                element_data(2, i) = rt2
                element_data(3, i) = 0.0
                lowest_mu = MIN(lowest_mu, rt1)
                highest_mu = MAX(highest_mu, rt1)
            ELSE ! try again to read a shorter record:
                !Try reading a record with only 1 per-element datum:
                READ (input_record, *, IOSTAT = ierr) j, n1, n2, n3, rt1
                IF (ierr == 0) THEN
                    element_data(1, i) = rt1
                    element_data(2, i) = 0.0
                    element_data(3, i) = 0.0
                    lowest_mu = MIN(lowest_mu, rt1)
                    highest_mu = MAX(highest_mu, rt1)
                ELSE ! try (for 3rd time) to read an even-shorter record, with only INTEGERs:
                    READ (input_record, *, IOSTAT = ierr) j, n1, n2, n3
                    IF (ierr == 0) THEN
                        element_data(1, i) = 0.0
                        element_data(2, i) = 0.0
                        element_data(3, i) = 0.0
                    ELSE
                        CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
                        WRITE (i_c10, "(I10)") i
                        i_c10 = ADJUSTL(i_c10)
                        msg1 = " Not enough data for element #" // TRIM(i_c10) // ' 'C
                        CALL error1(msg1)
	                    RETURN
                    END IF ! reading of just i, n1, n2, n3 good/bad
                END IF ! reading of i, n1, n2, n3, rt1 good/bad
            END IF ! reading of i, n1, n2, n3, rt1, rt2 good/bad
        END IF ! reading of i, n1, n2, n3, rt1, rt2, rt3 good/bad
    ELSE ! appType < 3; no per-element data allowed
        READ (input_record, *, IOSTAT = ierr) j, n1, n2, n3
        IF (ierr /= 0) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            i_c10 = ADJUSTL(i_c10)
            msg1 = " Not enough data for element #" // TRIM(i_c10) // ' 'C
            CALL error1(msg1)
	        RETURN
        END IF
    END IF
    IF (j /= i) THEN
        CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
        WRITE (i_c10, "(I10)") i
        WRITE (j_c10, "(I10)") j
        i_c10 = ADJUSTL(i_c10)
        j_c10 = ADJUSTL(j_c10)
        msg1 = " When expecting element #" // TRIM(i_c10) // ", encountered #" // TRIM(j_c10) // " instead."C
        CALL error1(msg1)
        RETURN
    END IF
	nodes(1, i) = n1
	nodes(2, i) = n2
	nodes(3, i) = n3
END DO eleInput ! i = 1, numel
IF (appType == 3) THEN ! check for lack-of-info case, and set mu-range to non-disastrous values (to set color-scale):
    IF ((lowest_mu == 1.0D0).AND.(highest_mu == 0.0D0)) THEN ! we never found ANY mu's in the input file; set defaults:
        lowest_mu = 0.0D-16
        highest_mu = 1.0D-14
    END IF
    highest_mu = MAX(highest_mu, (1.0D-15 + lowest_mu)) ! to guard against case where all mu's in input file were identical.
END IF
!
! Read faults and fault data (if any):
!
IF (appType == 1) THEN ! Shells-type .FEG, so the number of fault elements must be supplied as a header (even if 0).
    READ (100, *, IOSTAT = ierr) nfl
    IF (ierr /= 0) THEN
        CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
        WRITE (ierr_c10, "(I10)") ierr
        ierr_c10 = ADJUSTL(ierr_c10)
	    msg1 = ' At fault-header, IOSTAT = ' // TRIM(ierr_c10) // ' 'C
        CALL error1(msg1)
	    RETURN
    END IF
    IF (.NOT.ALLOCATED(nodef)) THEN
        ALLOCATE ( nodef(4, mxfel) )
        ALLOCATE ( fdip(2, mxfel) )
        ALLOCATE ( offset(mxfel) )
        ALLOCATE ( fault_LRi(mxfel) )
    END IF
    nodef = 0 ! whole arrays
    fdip = 0.0
    offset = 0.0    
    fault_LRi = 0
    !
    ! READ faults & fault data:
    !
    faultInput: DO i = 1, nfl
        READ (100, "(A)", IOSTAT = ierr) longer_line
        IF (ierr /= 0) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            i_c10 = ADJUSTL(i_c10)
            msg1 = " Not enough data for fault #" // TRIM(i_c10) // ' 'C
            CALL error1(msg1)
            RETURN
        END IF
        CALL Extract_LRi (longer_line, &     ! input
                        & LRi, shorter_line) ! ouput
        READ (shorter_line, *, IOSTAT = ierr) j, n1, n2, n3, n4, dips1, dips2, offst
        IF (ierr /= 0) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            i_c10 = ADJUSTL(i_c10)
            msg1 = " Not enough data for fault #" // TRIM(i_c10) // ' 'C
            CALL error1(msg1)
            RETURN
        END IF
        IF (j /= i) THEN
            CLOSE (UNIT=100) ! giving control back to the user, who will make corrections...
            WRITE (i_c10, "(I10)") i
            WRITE (j_c10, "(I10)") j
            i_c10 = ADJUSTL(i_c10)
            j_c10 = ADJUSTL(j_c10)
            msg1 = " When expecting fault #" // TRIM(i_c10) // ", encountered #" // TRIM(j_c10) // " instead."C
            CALL error1(msg1)
            RETURN
        END IF
        nodef(1, i) = n1
        nodef(2, i) = n2
        nodef(3, i) = n3
        nodef(4, i) = n4
        fdip(1, i) = dips1
        fdip(2, i) = dips2
        ! Dips are changed to [1 179] within OrbWin if dips < 0
        DO j = 1, 2
            IF (fdip(j, i) > 0.0) THEN
                !(no action needed)
            ELSE
	            fdip(j, i) = 180.0 + fdip(j, i)
            END IF
        END DO
        offset(i) = offst
        fault_LRi(i) = LRi
    END DO faultInput
ELSE ! appType > 1; not a Shells-type .FEG, so no faults allowed
    nfl = 0
END IF ! appType == 1, or not
!
CLOSE (UNIT=100)
!
! Set global control switches:
gridLoaded = .True.
plot_eleCenter_icons = (numel < 2000) ! usually .FALSE., since initial view of new grid always has windowheight = 2.2

CONTAINS

      SUBROUTINE Extract_LRi (longer_line, &     ! input
                            & LRi, shorter_line) ! output                     
      ! New routine added for Shells_v5.0+ to support multiple
      !"Lithospheric Rheology" (abbreviated as "LR") integer codes,
      ! in any line of the input .feg file which define an element
      !(either a triangular continuum element, or a
      !          linear fault element).
      ! CHARACTER*80, INTENT(IN) :: longer_line is the whole 
      !                             element-definition line from the .feg file.
      ! INTEGER, INTENT(OUT) :: LRi is the rheologic code
      !                        (or 0, if no such code was found).
      ! CHARACTER*80, INTENT(OUT) :: shorter_line has the " LRi" portion removed (if any),
      !                              so it can be interpreted by the same code as in Shells_v4.1-.
      IMPLICIT NONE
      CHARACTER*80, INTENT(IN) :: longer_line
      INTEGER, INTENT(OUT) :: LRi
      CHARACTER*80, INTENT(OUT) :: shorter_line
      CHARACTER*80 :: string
      INTEGER :: longer_length, LR_start_byte
      longer_length = LEN_TRIM(longer_line)
      LR_start_byte = INDEX(longer_line, "LR")
      IF (LR_start_byte > 0) THEN ! the "LR" flag was found
          IF (longer_length > (LR_start_byte + 1)) THEN ! some byte(s) follow the "LR"
              string = longer_line((LR_start_byte + 2):longer_length)
              READ (string, *) LRi
              shorter_line = longer_line(1:(LR_start_byte - 1))
          ELSE ! "LR" is present, but nothing follows it; infer 0.
              LRi = 0
              shorter_line = longer_line(1:(LR_start_byte - 1))
          END IF
      ELSE ! no "LR" flag is present
          LRi = 0
          shorter_line = longer_line
      END IF
      END SUBROUTINE Extract_LRi      
!
END SUBROUTINE Read_Grid