PROGRAM SETR4NSEQ1
C=======================================================================
C      DATA GENERATING PROGRAM FOR NSEQ8XX.FOR
C   PROJECT: DRIVEN CAVITY FLOW   PROJECT NAME: SEE BELOW DATA PROJECT
C    DOMAIN: SQUARE ( 1 X 1 )
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     REYNOLDS NUMBER: SPECIFIED IN PARAMETER AS REYNLD
C     VISCOSITY RATIO [VRATIO]: 2ND VISCOSITY / 1ST VISCOSITY
C     EIJI FUKUMORI DECEMBER 29, 1993 
C                   JAN.     29, 2013
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=4, MXE=12000, MXN=13000, MXB=1610, NF=3 )
      PARAMETER ( REYNLD=400.D0, VRATIO=1.D7, HEIGHT=1.D0, UAVE=1.D0,
     * VISCO=HEIGHT*UAVE/REYNLD,    FLMDA=VISCO*VRATIO,
     * DT = 0.01D0, TMAX=200.D0,  MAXITE=10, ERMAX=0.001D0, C1=0.D0,
     * TLX = HEIGHT,   TLY = HEIGHT,   NEY = 20, NEX = 20,
     * DX = TLX / NEX, DY = TLY / NEY,
     * NDX=NEX+1, NDY=NEY+1,  TIME = 0.D0, DTMAX=0.05D0, ITEFIX=3 )
      PARAMETER ( INTEPT=2 )
C=======================================================================
      DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN), U(MXN), V(MXN),
     * IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB),IBNDS(MXB),BVS(MXB)
      CHARACTER FNAME(NF)*19,EXTEN(NF)*4,PROJECT*15,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'DCF4REYNLD-0400' /
      DATA EXTEN / '.DAT', '.BIN', '.STM' /
C=======================================================================
      F(X) = X * X * ( 3.- 2.* X )
      G(X) = X * ( (0.1-1.)* X + (2.- 0.1) )
C=======================================================================
      FNAME(1)=PROJECT//EXTEN(1)
      FNAME(2)=PROJECT//EXTEN(2)
      FNAME(3)=PROJECT//EXTEN(3)
      OPEN ( UNIT=2, FILE='NSDATA4.FLN', STATUS='UNKNOWN')
      WRITE (2,'(A15)') PROJECT
      WRITE (2,'(A19)') FNAME(1)
      WRITE (2,'(A19)') FNAME(2)
      WRITE (2,'(A19)') FNAME(3)
      CLOSE (2)
      WRITE (*,*) PROJECT
      WRITE (*,*) FNAME(1)
      WRITE (*,*) FNAME(2)
      WRITE (*,*) FNAME(3)
C=======================================================================
C    C1: ADDED VISCOSITY COMPUTATION PARAMETER
C    C1=0. FOR NO ADDED VISCOSITY, C1=1. FOR UPWIND (OPTIMUM) ADDED 
C    C1 OPTION IS AVALABLE IN NSEQ8DD.FOR SOLVER
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C  MAXITE: MAXIMUM NUMBER OF ITERATIONS AT EACH TIME INCREMENT
C ITERATION PROCEDURE IS REQUIRED BECAUSE OF NON-LINEAR CONVECTIVE TERM.
C  ERMAX: MAXIMUM DIFFERENCE BETWEEN GUESSED {U&V} AND COMPUTED {U&V}
C  HEIGHT: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN
C=======================================================================
      WRITE (*,*)' DRIVEN CAVITY FLOW: REYNOLDS NUMBER OF ',REYNLD
      WRITE (*,*)' VISCOSITY RATIO OF ',VRATIO
C=======================================================================
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
      WRITE (*,*) 'NE=', NE
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO J = 1 , NDX
      DO I = 1 , NDY
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) = TLX * F( DX*(J-1)/TLX )
      YCOORD(NNODE) = TLY * G( DY*(I-1)/TLY )
      END DO
      END DO
      WRITE (*,*) 'NNODE=', NNODE
C=======================================================================
C                    BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS
C--------- FACE OF -X
      NBFX = 0
      NBFY = 0
      U0 = 0.
      V0 = 0.
      DO J = 2 , NDY-1
      NBFX = NBFX + 1
      NBFY = NBFY + 1
      IBNDFX(NBFX) = J
      IBNDFY(NBFY) = J
      BVX(NBFX) = U0
      BVY(NBFY) = V0
      END DO
C--------- FACE OF -Y
      U0 = 0.
      V0 = 0.
      DO I = 1 , NDX
      NBFX = NBFX + 1
      IBNDFX(NBFX) = 1 + NDY*(I-1)
      NBFY = NBFY + 1
      IBNDFY(NBFY) = 1 + NDY*(I-1)
      BVX(NBFX) = U0
      BVY(NBFY) = V0
      END DO
C--------- FACE OF +Y
      U0 = UAVE
      V0 = 0.
      DO I = 1 , NDX
      NBFX = NBFX + 1
      IBNDFX(NBFX) = I*NDY
      NBFY = NBFY + 1
      IBNDFY(NBFY) = I*NDY
      BVX(NBFX) = U0
      BVY(NBFY) = V0
      END DO
C--------- FACE OF +X
      U0 = 0.
      V0 = 0.
      DO I = 2 , NEY
      NBFX = NBFX + 1
      IBNDFX(NBFX) = NEX*NDY + I
      NBFY = NBFY + 1
      IBNDFY(NBFY) = NEX*NDY + I
      BVX(NBFX) = U0
      BVY(NBFY) = V0
      END DO
      WRITE (*,*) 'NBFX=', NBFX
      WRITE (*,*) 'NBFY=', NBFY
C=======================================================================
C                       INITIAL CONDITIONS
      U0 = 0.
      V0 = 0.
      DO I = 1 , NNODE
      U(I) = U0
      V(I) = V0
      END DO
C=======================================================================
C                     B. C. FOR STREAM FUNCTION
      DDY = YCOORD(NDY) - YCOORD(NDY-1)
      Q = DDY * UAVE / 2.
      NBS = NBFX
      DO I = 1 , NBS
      IBNDS(I) = IBNDFX(I)
      BVS(I) = BVX(I)
      IF ( BVX(I) .EQ. UAVE ) BVS(I) = Q
      END DO
C=======================================================================
C                     MAKING DATA FILES
C---------- HEADER
      OPEN ( 1, FILE=FNAME(1), STATUS = 'UNKNOWN' )
      WRITE(1,*) VISCO,  FLMDA
      WRITE(1,*) TMAX ,  DTMAX 
      WRITE(1,*) MAXITE, ERMAX,   ITEFIX
C---------- ELEMENT
      WRITE(1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
      END DO
C---------- NODE
      WRITE(1,*) NNODE
      DO I = 1 , NNODE
      WRITE(1,*) I, XCOORD(I), YCOORD(I)
      END DO
C---------- BOUNDARY
      WRITE(1,*) NBFX
      DO I = 1 , NBFX
      WRITE (1,*) IBNDFX(I), BVX(I)
      END DO
      WRITE(1,*) NBFY
      DO I = 1 , NBFY
      WRITE (1,*) IBNDFY(I), BVY(I)
      END DO
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BIN
      OPEN ( 1, FILE=FNAME(2),STATUS ='UNKNOWN',FORM='UNFORMATTED' )
      WRITE(1) TIME, DT
      WRITE(1) ( U(I)    , I = 1 , NNODE )
      WRITE(1) ( V(I)    , I = 1 , NNODE )
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.STM
      OPEN ( 1, FILE=FNAME(3), STATUS = 'UNKNOWN' )
      WRITE(1,*) NBS
      WRITE (1,*) ( IBNDS(I) , BVS(I) , I = 1 , NBS )
      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 ( NF=',NF,'    )'
      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,')'
      NB = MAX0 (NBFX, NBFY, NBS )
      WRITE (1,*) '      PARAMETER ( MXB=',NB,'   )'
      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