PROGRAM SETFEM9
C=======================================================================
C     DATA GENERATING PROGRAM FOR FEM4Q.FOR
C     DOMAIN SIZE: TLX BY TLY
C     BOUNDARY CONDITIONS
C        ON FACES OF -X AND -Y: T(X,Y)=0.
C        ON FACE OF +X: DT/DN = 0
C        ON FACE OF +Y: T(X,Y) = SIN ( 0.25*PI*X )
C     EIJI FUKUMORI NOVEMBER 28, 1995 
C                                  NUMBERING OF NODEAL POINTS
C  FOR CASE OF NEX=2 & NEY=2     21-----22-----23-----24-----25
C                                 |             |             |
C                                16     17     18     19     20
C                                 |             |             |
C                                11-----12-----13-----14-----15
C                                 |             |             |
C                                 6      7      8      9     10
C                                 |             |             |
C                                 1------2------3------4------5
C NOEDX OF 1ST ELEMENT = 1, 2, 3, 8, 13, 12, 11, 6, 7
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=9, MXE=130, MXN=441, MXB=100 )
      PARAMETER (   TLX = 2.,   TLY = 2. )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
     *          BVT(MXB), ITYPET(MXB)
      CHARACTER PROJECT*12,EXFILE*3
      LOGICAL YES
C=======================================================================
      DATA PROJECT / 'FEM9DATA.QQQ' /
C=======================================================================
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C=======================================================================
      WRITE (*,210)
  210 FORMAT( ' TYPE IN NUMBER OF ELEMENTS IN X-COORD (NEX) = ' $ )
      READ (*,*) NEX
      IF ( NEX .LT. 1 ) STOP 'NOT ACCEPTABLE NEX'
      WRITE (*,220)
  220 FORMAT( ' TYPE IN NUMBER OF ELEMENTS IN Y-COORD (NEY) = '  $ )
      READ (*,*) NEY
      IF ( NEY .LT. 1 ) STOP 'NOT ACCEPTABLE NEY'
      WRITE(*,*)' NEX =',NEX, '     NEY =', NEY
C
      DX = TLX / NEX
      DY = TLY / NEY
       EXX=1.
       EXY = 0.
       EYY=1.
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEY
      DO J = 1 , NEX
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = (I-1)*(NEX*4+2) + (j-1)*2 + 1
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,3) = NODEX(NE,2) + 1
      NODEX(NE,4) = NODEX(NE,3) + NEX*2 + 1
      NODEX(NE,5) = NODEX(NE,4) + NEX*2 + 1
      NODEX(NE,6) = NODEX(NE,5) - 1
      NODEX(NE,7) = NODEX(NE,6) - 1
      NODEX(NE,8) = NODEX(NE,7) - (NEX*2+1)
      NODEX(NE,9) = NODEX(NE,8) + 1
      END DO
      END DO
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NEY*2+1
      DO J = 1 , NEX*2+1
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) =  DX/2.*(J-1)
      YCOORD(NNODE) =  DY/2.*(I-1)
      END DO
      END DO
C=======================================================================
      PI = 4.D0* ATAN (1.D0)
C                    BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS AND HEAT EQUATION
      NBT  = 0
C--------- FACE OF +Y AND -Y
      DO J = 1 , NEX*2+1
C--------- -Y
      NBT = NBT + 1
      IBNDT(NBT) = J
      ITYPET(NBT) = 1
      BVT(NBT) = 0.
C--------- +Y
      NBT = NBT + 1
      IBNDT(NBT) = NEY*(NEX*4 + 2) + J
      ITYPET(NBT) = 1
      BVT(NBT) = DSIN ( 0.25D0*PI*XCOORD(IBNDT(NBT)))
      END DO
C--------- FACE OF +X AND -X
C           -X
      DO I = 1 , NEY
      NBT  = NBT  + 1
      IBNDT (NBT) = (I-1)*(NEX*4+2) + NEX*2 + 1 + 1
      ITYPET(NBT) = 1
      BVT(NBT ) = 0.
      NBT  = NBT  + 1
      IBNDT (NBT) = (I-1)*(NEX*4+2) + 2*(NEX*2 + 1) + 1
      ITYPET(NBT) = 1
      BVT(NBT ) = 0.
      END DO
      NBT = NBT - 1
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
      IR = 1
      OPEN ( IR, FILE=PROJECT, STATUS = EXFILE )
      WRITE (IR,*) EXX, EXY, EYY
      WRITE (IR,*) NE
      DO IEL = 1 , NE
        WRITE (IR,'(10I5)') IEL,(NODEX(IEL,J),J=1,ND)
      END DO
      WRITE (IR,*) NNODE
      DO NODE = 1 , NNODE
        WRITE (IR,*) NODE,XCOORD(NODE),YCOORD(NODE)
      END DO
      WRITE (IR,*) NBT
      DO I = 1 , NBT
        WRITE (IR,*) IBNDT(I), ITYPET(I), BVT(I)
      END DO
      CLOSE (1)
      STOP
      END