PROGRAM SETMBC1
C=======================================================================
C               DATA GENERATING PROGRAM FOR FEM4Q.FOR
C                   DOMAIN : TLX BY TLY
C                  BOUNDARY CONDITIONS
C                       ON FACES OF -X : H1
C                       ON FACES OF -Y: Q=0.
C                       ON FACE OF +X: H2 AND SEEPAGE FACE
C                       ON FACE OF +Y: MOVING WITH Q=0.
C                   EIJI FUKUMORI NOVEMBER 28, 1995
C                         FEBRUARY 06, 2013
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER (ND=4,MXE=4000, MXN=5000,MXB=1000,MXG=100)
      PARAMETER (   TLX = 5.,   TLY = 10. )
      PARAMETER ( INTEPT=2, MXMVG=100,MXGEACH=100 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
     *          BVT(MXB), ITYPET(MXB),NMOVE(MXG), MOVEND(MXG,MXGEACH)
      CHARACTER PROJECT*12,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'MOVINGBC.DAT' /
C=======================================================================
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C=======================================================================
      MAXIT=200
      ERALLOW=0.0005D0
C
      NEX = 20
      NEY = 20
C
      DX = TLX / NEX
      DY = TLY / NEY
      NDX=NEX+1
      NDY=NEY+1
       EXX = 1.D0
       EXY = 0.D0
       EYY = 1.D0
C
C---------------NDMAXH2 = # OF NODES IN Y DOWN-STREAM HEAD IS ASSIGNED
      H1 = TLY
      NDMAXH2 = 0.25D0*NEY+1
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEX
      DO J = 1 , NEY
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = NDY*(I-1) + J
      NODEX(NE,2) = NODEX(NE,1) + NDY
      NODEX(NE,3) = NODEX(NE,2) + 1
      NODEX(NE,4) = NODEX(NE,3) - NDY
      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'
      NODE = NDY*(I-1) + J
      XCOORD(NODE) =  DX*(I-1)
      YCOORD(NODE) =  DY*(J-1)
      END DO
      END DO
C=======================================================================
      NBT  = 0
C--------- FACE OF -X
      DO J = 1 , NDY
      NBT = NBT + 1
      IBNDT(NBT) = J
      ITYPET(NBT) = 1
      BVT(NBT) = H1          
      END DO
C--------- FACE OF -Y
C   NO FLOW; Q=0
C--------- FACE OF +X
C== HYDRAULIC HEAD = CONSTANT NODE
      H2 = YCOORD(NDY*NEX+NDMAXH2)
      DO J = 1 , NDMAXH2
      NBT  = NBT  + 1
      IBNDT (NBT) = NDY*NEX + J
      ITYPET(NBT) = 1
      BVT(NBT ) = H2
      END DO
C== SEEPAGE FACE <---------------- ITYPET(NBT) = 3 -----------------
      DO J = NDMAXH2+1, NEY
      NBT  = NBT  + 1
      IBNDT (NBT) = NDY*NEX + J
      ITYPET(NBT) = 3
      BVT(NBT ) = YCOORD (IBNDT(NBT))
      END DO
C--------- FACE OF +Y
C   NO FLOW; Q=0 AND FREE SURFACE
C
C---------- MOVING NODES (GROUPING)
C*******************************************************************
      NGROUP = 0
      DO I = 2 , NDX
      II = I-1
      NGROUP = NGROUP + 1
      NMOVEI = 0
      ISTART = 1
      IF ( I .EQ. NDX ) ISTART = NDMAXH2
      DO J = ISTART, NDY
      NMOVEI = NMOVEI + 1
      JJ = J - ISTART+1
      MOVEND(II,JJ) = NDY*(I-1)+J
      END DO
      NMOVE(II) = NMOVEI
      END DO
C=======================================================================
      WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
      WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = PROJECT, EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
C---------- 'PROJECT'.JNK
      OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
C---------- MATERIAL PARAMETER
      WRITE (1,*) EXX, EXY, EYY
C---------- CONVERGENCE PARAMETER
      WRITE (1,*) MAXIT,ERALLOW
C---------- ELEMENT INFORMATION
      WRITE (1,*) NE
      DO IEL = 1 , NE
        WRITE (1,*) IEL,(NODEX(IEL,J),J=1,ND)
      END DO
C---------- NODAL COORDINATE INFOMATION
      WRITE (1,*) NNODE
      DO NODE = 1 , NNODE
        WRITE (1,*) NODE,XCOORD(NODE),YCOORD(NODE)
      END DO
C---------- BOUNDARY CONDITIONS
      WRITE (1,*) NBT
      DO I = 1 , NBT
        WRITE (1,*) IBNDT(I), ITYPET(I), BVT(I)
      END DO
C---------- MOVING NODE INFORMATION
      WRITE (1,*) NGROUP
      DO I = 1 , NGROUP
      NNMOVE = NMOVE(I)
      WRITE (1,*) NNMOVE, (MOVEND(I,J), J=1,NNMOVE)
      END DO
      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=',NBT,')'
      WRITE (1,*) '      PARAMETER ( MXMVG=',NGROUP,')'
      NMOVEMX = 0
      DO I = 1 , NGROUP
      NMOVEMX = MAX0 ( NMOVEMX , NMOVE(I) )
      END DO
      WRITE (1,*) '      PARAMETER ( MXMVEACH=',NMOVEMX,')'
      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