C     ==================================================================
C     [][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
C     ==================================================================
      SUBROUTINE GRID(X1,X2,NEC,EEC)
      IMPLICIT NONE
      INTEGER NP,NE
      PARAMETER (NP=5433,NE=5262)
      INTEGER NEC(NE,4),NPC(NP,5),EEC(NE,4)
      REAL*8  X1(NP),X2(NP)
      REAL*8  XL_1,XL_2
      INTEGER I,I1,I2,I3,I4,I5,I6,L,L0,L1,M,M0,M1,M2,N,K,K0,K1,K2


C-----==================================================================
C-----READ MESH FILE (GAMBIT) ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      OPEN(10,FILE='bump.neu',FORM='FORMATTED')
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      READ(10,*)
      DO I=1,NP 
      READ(10,*) I1,X1(I),X2(I)
      END DO
      READ(10,*)
      READ(10,*)
      DO I=1,NE
      READ(10,*) I1,I2,I3,NEC(I,1),NEC(I,2),NEC(I,3),NEC(I,4)
      END DO
      CLOSE(10)

      
C-----==================================================================
C-----NODE ELEMENT CONNECTIVITY ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      NPC(:,:)=0
      DO     I=1,NE
      DO 100 L=1,4
      IF (NPC(NEC(I,L),1).EQ.0) THEN
      NPC(NEC(I,L),1)=I ;GOTO 100
      END IF
      IF (NPC(NEC(I,L),2).EQ.0) THEN
      NPC(NEC(I,L),2)=I ;GOTO 100
      END IF
      IF (NPC(NEC(I,L),3).EQ.0) THEN
      NPC(NEC(I,L),3)=I ;GOTO 100
      END IF
      IF (NPC(NEC(I,L),4).EQ.0) THEN
      NPC(NEC(I,L),4)=I ;GOTO 100
      END IF
      IF (NPC(NEC(I,L),5).EQ.0) THEN
      NPC(NEC(I,L),5)=I ;GOTO 100
      END IF
      WRITE(6,*) NEC(I,L),I,"STOP...MORE THAN FIVE CONNECTION."
  100 END DO
      END DO


C-----==================================================================
C-----REARANGE THEM IN CORRECT ORDER ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      DO 200 I=1,NP
C-----INITIAL K0 VALUE ON BOUNDARY ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      DO K=1,5
      IF (NPC(I,K).NE.0) THEN
      DO M=1,4
      IF (NEC(NPC(I,K),M).EQ.I) M0=M
      END DO
      IF (M0.EQ.4) M1=1
      IF (M0.NE.4) M1=M0+1
      K0=0
      DO N=1,5
      IF (NPC(I,N).NE.0.AND.N.NE.K) THEN
      DO M=1,4
      IF (NEC(NPC(I,K),M1).EQ.NEC(NPC(I,N),M)) K0=N
      END DO
      ENDIF
      END DO
      IF (K0.EQ.0) THEN
      I1=NPC(I,K)
      NPC(I,K)=NPC(I,1)
      NPC(I,1)=I1
      ENDIF
      ENDIF
      END DO

C-----KI VALUES ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      DO K=2,5
      IF (NPC(I,K).NE.0) THEN
      L0=0
      DO L=1,4
      IF (NEC(NPC(I,K-1),L).EQ.I) L0=L
      END DO
      IF (L0.EQ.0) PAUSE
      IF (L0-1.EQ.0) THEN 
      L1=4
      ELSE
      L1=L0-1
      ENDIF
      DO M=K,5
      IF (NPC(I,M).NE.0) THEN
      DO L=1,4
      IF (NEC(NPC(I,M),L).EQ.NEC(NPC(I,K-1),L1)) THEN
      I1=NPC(I,K)
      NPC(I,K)=NPC(I,M)
      NPC(I,M)=I1
      END IF
      END DO
      END IF
      END DO
      END IF
      END DO
C     WRITE(10,*) I,NPC(I,1),NPC(I,2),NPC(I,3),NPC(I,4),NPC(I,5)
  200 END DO


C-----==================================================================
C-----ELELEMT ELEMENT CONNECTION ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      EEC(:,:)=0
      DO I=1,NE
      DO 300 M1=1,4
      M2=M1+1
      IF (M1.EQ.4) M2=1
      DO K1=1,5
      DO K2=1,5
      IF (NPC(NEC(I,M1),K1).NE.0.AND.NPC(NEC(I,M2),K2).NE.0) THEN
      IF (NPC(NEC(I,M1),K1).NE.I.OR .NPC(NEC(I,M2),K2).NE.I) THEN
      IF (NPC(NEC(I,M1),K1).EQ.NPC(NEC(I,M2),K2)) THEN
      EEC(I,M1)=NPC(NEC(I,M1),K1)
      END IF
      END IF
      END IF
      END DO
      END DO
  300 END DO
      END DO


C-----==================================================================
C-----BOUNDARY CONDITION ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      DO I=1,NE
      DO M1=1,4
      M2=M1+1
      IF (M1.EQ.4) M2=1
      XL_1=0.5D00*(X1(NEC(I,M1))+X1(NEC(I,M2)))
      XL_2=0.5D00*(X2(NEC(I,M1))+X2(NEC(I,M2)))

C-----INFLOW BOUNDARY CONDITIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      IF (XL_1.LT.-1.4999999D00) THEN
      EEC(I,M1)=-1
      END IF

C-----OUTFLOW BOUNDARY CONDITIONS ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      IF (XL_1.GT. 1.4999999D00) THEN
      EEC(I,M1)=-2
      END IF
      END DO
      END DO


C-----==================================================================
C-----MESH OUTPUT FILE (TECHPLOT FORMAT) ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
C-----==================================================================
      OPEN(10,FILE='MESH.DAT',FORM='FORMATTED')
      WRITE(10,*) 'TITLE="CAVITY"'
      WRITE(10,*) 'VARIABLES="X","Y"'
      WRITE(10,*) 'ZONE  N=',NP,'E=',NE,'F=FEPOINT, ET=QUADRILATERAL'
      DO I=1,NP
      WRITE(10,*) X1(I),X2(I)
      END DO
      DO I=1,NE
      WRITE(10,*) NEC(I,1),NEC(I,2),NEC(I,3),NEC(I,4)
      END DO
      CLOSE(10)
      RETURN
      END