PROGRAM SET4IMP1
C=======================================================================
C      DATA GENERATING PROGRAM FOR IMP4Q.FOR
C   PROJECT: CONFINED GROUNDWATER   PROJECT NAME: IMP4DATA.QQQ
C    DOMAIN: 
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     EIJI FUKUMORI DECEMBER 29, 1993
C                   FEBRUARY 05, 2013
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, MXE=20000, MXN=30000, MXB=3000 )
      PARAMETER ( HEIGHT=1.D0, TLX = HEIGHT,TLY = HEIGHT,
     * NEY = 20, NEX = 20, DX = TLX / NEX, DY = TLY / NEY,WAVEH=1.,
     * ALPHA=100.D0, NDX=NEX+1,NDY=NEY+1)
      PARAMETER ( INTEPT=2 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN),IBNDF(MXB),
     * BV(MXB), ITYPE(MXB), H(MXN)
      CHARACTER PROJECT*12, EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'IMP4Q.DAT' /
C=======================================================================
C      F(X) = X * X * ( 3.D0- 2.D0* X )
C      G(X) = X * ( (0.1D0-1.D0)* X + (2.D0- 0.1D0) )
      F(X) = X
      G(X) = X
C=======================================================================
      WRITE (*,*)' INPUT DATA FILE = ', PROJECT
C=======================================================================
C                    BASIC PARAMTERS
C------   EXX, EXY, EYY = TRANMISSIVITY
      EXX = 1.D0
      EXY = 0.D0
      EYY = 1.D0
C------ STORAGE = STRATIVITY
      STORAGE = 1.D0
      DT = 0.1D0
      MXNTIME = 100
      THETA = 0.5D0
C                    ELEMENT CREATION
      NE = 0
      DO J = 1 , NEX
      DO I = 1 , NEY
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = NDY*(J-1) + I
      NODEX(NE,2) = NODEX(NE,1) + NDY
      NODEX(NE,3) = NODEX(NE,2) + 1
      NODEX(NE,4) = NODEX(NE,1) + 1
      END DO
      END DO
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NDX
      DO J = 1 , NDY
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) = TLX * F( DX*(I-1)/TLX )
      YCOORD(NNODE) = TLY * G( DY*(J-1)/TLY )
      END DO
      END DO
C=======================================================================
C                    BOUNDARY CONDITIONS
C--------- FACE OF -X
      NBF = 0
      NDYPRT1 = NDY/3
CCCCCCC      DO J = 2 , NDY-1
      DO J = 1, NDYPRT1
      NBF = NBF + 1
      IBNDF(NBF) = J
      ITYPE(NBF) = 1
      BV(NBF) = WAVEH
      END DO
C      DO J = NDYPRT1+1,NDY-1
C      NBF = NBF + 1
C      IBNDF(NBF) = J
C      ITYPE(NBF) = 1
C      BV(NBF) = 0.D0
C      END DO
C--------- FACE OF -Y
C      DO I = 1 , NDX
C      NBF = NBF + 1
C      IBNDF(NBF) = 1 + NDY*(I-1)
C      ITYPE(NBF) = 1
C      BV(NBF) = 0.D0
C      END DO
C--------- FACE OF +Y
      DO I = 1 , NDX
      NBF = NBF + 1
      IBNDF(NBF) = I*NDY
      ITYPE(NBF) = 1
      BV(NBF) = 0.D0
      END DO
C--------- FACE OF +X
      DO I = 2 , NEY
      NBF = NBF + 1
      IBNDF(NBF) = NEX*NDY + I
      ITYPE(NBF) = 1
      BV(NBF) = 0.D0
      END DO
C=======================================================================
C                      INITIAL CONDITION
      DO I = 1 , NNODE
      H(I) = 0.D0
      END DO
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = PROJECT, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
C------ CONTROL PARAMETERS
      WRITE (1,*) EXX, EXY, EYY
      WRITE (1,*) STORAGE
      WRITE (1,*) DT, MXNTIME
      WRITE (1,*) THETA
C------ ELEMENTS
      WRITE(1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
      END DO
C------ NODAL COORDINATES
      WRITE(1,*) NNODE
      DO I = 1 , NNODE
      WRITE(1,*) I,XCOORD(I), YCOORD(I)
      END DO
C------ BOUNDARY CONDITIONS
      WRITE(1,*) NBF
      DO I = 1 , NBF
      WRITE (1,*) IBNDF(I), ITYPE(I), BV(I)
      END DO
      CLOSE (1)
C------ INITIAL CONDITION
      OPEN (1,FILE='SOLUTION.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE(1) ( H(I) , I=1,NNODE )
      CLOSE (1)
C=======================================================================
      CALL BANDWID ( ND, MXE, NE, NODEX , NBW  )
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
      OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) '      PARAMETER ( ND=',ND,')'
      WRITE (1,*) '      PARAMETER ( INTEPT=',INTEPT,')'
      WRITE (1,*) '      PARAMETER ( MXE=',NE,')'
      WRITE (1,*) '      PARAMETER ( MXN=',NNODE,')'
      WRITE (1,*) '      PARAMETER ( MXW=',NBW,')'
      WRITE (1,*) '      PARAMETER ( MXB=',NBF,')'
      CLOSE (1)
C=======================================================================
      STOP "NORMAL TERMINATION"
      END
C
C
      SUBROUTINE BANDWID ( ND, MXE, NE, NODEX , NBW  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND)
      NBW = 0
      DO  I = 1 , NE
      DO  J = 1 , ND-1
      DO  K = J+1 , ND
      NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
      END DO
      END DO
      END DO
      NBW = NBW + 1
      WRITE(*,*) 'HALF BANDWIDTH + 1 =', NBW
      RETURN
      END