PROGRAM SET3D27NSEQ1_NEWFORMAT
C=======================================================================
C      DATA GENERATING PROGRAM FOR 3D-NSEQ8XX.FOR
C   PROJECT: DRIVEN CAVITY FLOW   PROJECT NAME: DCFR999
C    DOMAIN: SQUARE ( 1 X 1 X 1 ) 
C    ELEMENT: 3D HEXA 27-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 JULY 11, 2010
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=27, MXE=22000, MXN=53000, MXB=20000, NF=9 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN),
     * U(MXN),V(MXN),W(MXN),
     * IBNDFX(MXB),IBNDFY(MXB),IBNDFZ(MXB),BVX(MXB),BVY(MXB),BVZ(MXB)
      CHARACTER FNAME(NF)*11,EXTEN(NF)*4,PROJECT*7,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'DCFR999' /
      DATA EXTEN / '.JNK', '.ELE', '.NOD', '.BND', '.INI', '.BIN',
     *        '.XXX', '.MAX', '.STM' /
C=======================================================================
      F(X) = X * X * ( 3.D0- 2.D0* X )
      G(X) = X * ( (0.1D0-1.D0)* X + (2.D0- 0.1D0) )
C=======================================================================
C----------> FLOW DOMAIN PARARETERS
      TLX = 1.D0
      TLY = 1.D0
      TLZ = 1.D0
C----------> ELEMENT INFORMATION
      NEX = 10
      NEY = 10
      NEZ = 10
C----------> ELEMENT INFORMATION CONTINUED
      DX = TLX/(2*NEX)
      DY = TLY/(2*NEY)
      DZ = TLZ/(2*NEZ)
      NDX = 2*NEX + 1
      NDY = 2*NEY + 1
      NDZ = 2*NEZ + 1
C=======================================================================
C----------> FLUID PARAMETERS
      REYNLD = 1000.D0
      VRATIO = 1.D7
      UAVE = 1.D0
      VISCO = TLZ*UAVE/REYNLD
      FLMDA = VISCO*VRATIO
C=======================================================================
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  TLZ: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN
C=======================================================================
      CALL FILEMG ( NF, FNAME, EXTEN, PROJECT )
      WRITE (*,*)'DRIVEN CAVITY FLOW'
      WRITE (*,*)'REYNOLDS NUMBER OF ',REYNLD
      WRITE (*,*)'VISCOSITY RATIO OF ',VRATIO
      WRITE (*,*)'VISCOSITY ',VISCO 
      WRITE (*,*)'2ND VISCOSITY ',FLMDA
C=======================================================================
      DT = 0.005D0
      TMAX = 200.D0
      MAXITE = 10
      TIME = 0.D0
      ERMAX = 0.001D0
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NDX
      DO K = 1 , NDZ
      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 * F( DY*(J-1)/TLY )
      ZCOORD(NNODE) = TLZ * G( DZ*(K-1)/TLZ )
      IF ( K .EQ. NDZ ) ZCOORD(NNODE) = TLZ
      END DO
      END DO
      END DO
      WRITE (*,*) 'NNODE =', NNODE
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEX
      DO J = 1 , NEZ
      DO K = 1 , NEY
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
C------ BOTTOM
      NODEX(NE,1) = (I-1)*NDY*NDZ*2 + (J-1)*NDY*2 + (K-1)*2 + 1
      NODEX(NE,2) = NODEX(NE,1) + NDY*NDZ
      NODEX(NE,3) = NODEX(NE,2) + NDY*NDZ
      NODEX(NE,4) = NODEX(NE,3) + 1
      NODEX(NE,5) = NODEX(NE,4) + 1
      NODEX(NE,8) = NODEX(NE,1) + 1
      NODEX(NE,7) = NODEX(NE,8) + 1
      NODEX(NE,9) = NODEX(NE,2) + 1
      NODEX(NE,6) = NODEX(NE,9) + 1
C------ MID
      NODEX(NE,10) = NODEX(NE,1) + NDY
      NODEX(NE,17) = NODEX(NE,10) + 1
      NODEX(NE,16) = NODEX(NE,17) + 1
      NODEX(NE,11) = NODEX(NE,10) + NDY*NDZ
      NODEX(NE,18) = NODEX(NE,11) + 1
      NODEX(NE,15) = NODEX(NE,18) + 1
      NODEX(NE,12) = NODEX(NE,11) + NDY*NDZ
      NODEX(NE,13) = NODEX(NE,12) + 1
      NODEX(NE,14) = NODEX(NE,13) + 1
C------ TOP
      NODEX(NE,19) = NODEX(NE,10) + NDY
      NODEX(NE,26) = NODEX(NE,19) + 1
      NODEX(NE,25) = NODEX(NE,26) + 1
      NODEX(NE,20) = NODEX(NE,19) + NDY*NDZ
      NODEX(NE,27) = NODEX(NE,20) + 1
      NODEX(NE,24) = NODEX(NE,27) + 1
      NODEX(NE,21) = NODEX(NE,20) + NDY*NDZ
      NODEX(NE,22) = NODEX(NE,21) + 1
      NODEX(NE,23) = NODEX(NE,22) + 1
      END DO
      END DO
      END DO
      WRITE (*,*) 'NE =', NE
C=======================================================================
      CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
C=======================================================================
C                    BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS
      NBFX = 0
      NBFY = 0
      NBFZ = 0
C--------- FACE OF -X
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -X
      DO K = 2 , NDZ-1
      DO J = 1 , NDY
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = J + (K-1)*NDY
      BVX(NBFX) = 0.D0
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = J + (K-1)*NDY
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = J + (K-1)*NDY
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +X
      DO K = 2 , NDZ-1
      DO J = 1 , NDY
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
      BVX(NBFX) = 0.D0
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -Y
      DO K = 2 , NDZ-1
      DO I = 2 , NDX-1
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
      BVX(NBFX) = 0.D0
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +Y
      DO K = 2 , NDZ-1
      DO I = 2 , NDX-1
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = K*NDY + (I-1)*NDY*NDZ
      BVX(NBFX) = 0.D0
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = K*NDY + (I-1)*NDY*NDZ
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = K*NDY + (I-1)*NDY*NDZ
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -Z
      DO I = 1 , NDX
      DO J = 1 , NDY
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = J + (I-1)*NDY*NDZ
      BVX(NBFX) = 0.D0
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = J + (I-1)*NDY*NDZ
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = J + (I-1)*NDY*NDZ
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +Z
      DO I = 1 , NDX
      DO J = 1 , NDY
C
      NBFX = NBFX + 1
      IBNDFX(NBFX) = I*NDY*NDZ - NDY + J
      BVX(NBFX) = UAVE
C
      NBFY = NBFY + 1
      IBNDFY(NBFY) = I*NDY*NDZ - NDY + J
      BVY(NBFY) = 0.D0
C
      NBFZ = NBFZ + 1
      IBNDFZ(NBFZ) = I*NDY*NDZ - NDY + J
      BVZ(NBFZ) = 0.D0
C
      END DO
      END DO

      WRITE (*,*) 'NBFX =', NBFX
      WRITE (*,*) 'NBFY =', NBFY
      WRITE (*,*) 'NBFZ =', NBFZ
C=======================================================================
C                       INITIAL CONDITIONS
      DO I = 1 , NNODE
      U(I) = 0.D0
      V(I) = 0.D0
      W(I) = 0.D0
      END DO
C=======================================================================
C                      DATA FILE INQUIRY
      EXFILE = 'NEW'
      INQUIRE ( FILE = FNAME(1), EXIST = YES )
      IF ( YES ) EXFILE='OLD'
C=======================================================================
C                     MAKING DATA FILES
C---------- 'PROJECT'.JNK
      OPEN ( 1, FILE=FNAME(1), STATUS = EXFILE )
      WRITE(1,*) VISCO, FLMDA
      WRITE(1,*) DT, TMAX
      WRITE(1,*) MAXITE, ERMAX
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.ELE
      OPEN ( 1, FILE=FNAME(2), STATUS = EXFILE )
      WRITE(1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
      END DO
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.NOD
      OPEN ( 1, FILE=FNAME(3), STATUS = EXFILE )
      WRITE(1,*) NNODE
      DO  I = 1, NNODE
      WRITE(1,*) I,XCOORD(I), YCOORD(I), ZCOORD(I)
      END DO
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BND
      OPEN ( 1, FILE=FNAME(4), STATUS = EXFILE )
      WRITE(1,*) NBFX
      DO I = 1 , NBFX
      WRITE (1,*) IBNDFX(I), BVX(I)
      END DO
C
      WRITE(1,*) NBFY
      DO  I = 1 , NBFY
      WRITE (1,*) IBNDFY(I), BVY(I)
      END DO
C
      WRITE(1,*) NBFZ
      DO  I = 1 , NBFZ
      WRITE (1,*) IBNDFZ(I), BVZ(I)
      END DO
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BIN
      OPEN ( 1, FILE=FNAME(6), STATUS = EXFILE,FORM='UNFORMATTED' )
      WRITE(1) TIME
      WRITE(1) ( U(I)    , I = 1 , NNODE )
      WRITE(1) ( V(I)    , I = 1 , NNODE )
      WRITE(1) ( W(I)    , I = 1 , NNODE )
      CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.MAX
      OPEN ( 1, FILE=FNAME(8), STATUS = EXFILE )
      UMIN = 0.
      UMAX = 0.
      VMIN = 0.
      VMAX = 0.
      PMIN = 0.
      PMAX = 0.
      SMIN = 0.
      SMAX = 0.
      WRITE (1,*) UMIN , UMAX
      WRITE (1,*) VMIN , VMAX
      WRITE (1,*) PMIN , PMAX
      WRITE (1,*) SMIN , SMAX
      CLOSE (1)
C=======================================================================
      STOP
      END
C
C
      SUBROUTINE FILEMG ( NF, FNAME, EXTEN, PROJECT )
      CHARACTER FNAME(NF)*11, PROJECT*7, EXTEN(NF)*4, FNWOE(7)*1,
     *          EXFILE*3
      LOGICAL YES
C------- RETURN VALUE: FNAME
C
      EXFILE ='NEW'
      INQUIRE ( FILE='NSDATA.FLN', EXIST=YES )
      IF ( YES ) EXFILE='OLD'
      OPEN ( UNIT=2, FILE='NSDATA.FLN', STATUS=EXFILE )
C
      NCHAR = 7
      DO I = 1 , NCHAR
      FNWOE(I) = PROJECT (I:I)
      END DO
      CALL GENFNM ( NCHAR, FNWOE, NF, FNAME, EXTEN )
      WRITE (2,'(I5)') NCHAR
      WRITE (2,'(7A1)') ( FNWOE(I) , I = 1 , NCHAR )
      DO I = 1 , NF
      WRITE (*,'(1X,A11)') FNAME(I)
      WRITE (2,'(   A11)') FNAME(I)
      END DO
      CLOSE (2)
      RETURN
      END
C
C
      SUBROUTINE GENFNM ( NCHAR, FNWOE, NF, FNAME, EXTEN )
      CHARACTER*1  FNWOE(7)
      CHARACTER*4  EXTEN(NF)
      CHARACTER*11 FNAME(NF)
      DO I = 1 , NF
      DO J = 1 , NCHAR
      FNAME (I) (J:J) = FNWOE(J)
      END DO
      DO J = NCHAR+1,NCHAR+4
      FNAME (I) (J:J) = EXTEN(I) (J-NCHAR:J-NCHAR)
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
      DIMENSION NODEX(MXE,ND)
C------- RETURN VALUE: NBW
      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 (*,*) ' HALH BANDWIDTH =', NBW
      RETURN
      END