P R O G R A M     SET3DSOUND1
C********************************************************************
C------------------- 8-NODED HEXA LINEAR ELEMENT --------------------
C********************************************************************
      PARAMETER( ND=8, INTEPT=2 )
      PARAMETER( IW=10,MXNMT=10,MXE=1000000,MXN=2000000 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN), ZCOORD(MXN)
      CHARACTER INFNAME*12
C
C--------------------------------------------------------------------
      INFNAME = 'EIGN3D8.DAT'
C====================================================================
C------------------------ SPEED OF SOUND ----------------------------
      WSPD = 340.D0
C-------------------------------------------------------
C TLX=LENGTH, TLY=HEIGHT, TLZ=DEPTH
      TLX = 4.D0*DATAN(1.D0)
      TLY = TLX/2.D0
      TLZ = TLX/3.D0
C----------------------------------------------------------
      NELEMENT = 4
      ALMIN = DMIN1 ( TLX, TLY, TLZ )
      NEX = NELEMENT * (TLX/ALMIN + 0.001D0)
      NEY = NELEMENT * (TLY/ALMIN + 0.001D0)
      NEZ = NELEMENT * (TLZ/ALMIN + 0.001D0)
C==========================================================
C--- APPROPRIATE VALUES:DELTA=0.1, MXNEIGEN=40
      AL = DMAX1 ( TLX, TLY, TLZ )
      DELTA = (4.D0*DATAN(1.D0)/AL)**2 * 0.9
C----------------------------------------------------------
      DX = TLX / NEX 
      DY = TLY / NEY
      DZ = TLZ / NEZ
C
      NDX = NEX + 1
      NDY = NEY + 1
      NDZ = NEZ + 1
C
      DO J = 1 , NEY
      DO K = 1 , NEZ
      DO I = 1 , NEX
      NE = NEX*(K-1)+I+(J-1)*NEZ*NEX
      NODEX(NE,1) = I + NDX*(K-1) + (J-1)*NDZ*NDX
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,5) = NODEX(NE,1) + NDX
      NODEX(NE,6) = NODEX(NE,5) + 1
      NODEX(NE,4) = NODEX(NE,1) + NDZ*NDX
      NODEX(NE,3) = NODEX(NE,4) + 1
      NODEX(NE,8) = NODEX(NE,4) + NDX
      NODEX(NE,7) = NODEX(NE,8) + 1
      END DO
      END DO
      END DO
C--------------------------------------------------------
C
C           COORDINATES OF NODE POINTS
      DO J = 1 , NDY
      DO K = 1 , NDZ
      DO I = 1 , NDX
      NNODE = I + (K-1) * NDX + (J-1) * NDZ * NDX
      XCOORD(NNODE) = DX * (I-1)
      YCOORD(NNODE) = DY * (J-1)
      ZCOORD(NNODE) = DZ * (K-1)
      END DO
      END DO
      END DO
C-----------------------------------------------------
      MXNEIGEN=MIN0 ( 20, NNODE/10 )
      IF ( MXNEIGEN .LT. 9 ) MXNEIGEN = 9
C----------------------------------------------------------
      DXMIN = DMIN1 ( DX, DY, DZ )
      EIGENMX = ( (4.D0*DATAN(1.D0))/(2.D0*DXMIN) ) ** 2
      TOTALVOL = TLX * TLY * TLZ
C----------------------------------------------------------
      IWW = 2
      OPEN  (IWW, FILE='DOMAIN.SUM', STATUS='UNKNOWN')
      WRITE (IWW,*) 'SPEED OF SOUND =', WSPD
      WRITE (IWW,*) 'TLX, TLY, TLZ =', TLX, TLY, TLZ
      WRITE (IWW,*) 'TOTAL VOLUME =', TOTALVOL
      WRITE (IWW,*) 'NEX, NEY, NEZ =', NEX, NEY, NEZ
      WRITE (IWW,*) 'MAXIMUM SURFACE LENGTH =', AL
      WRITE (IWW,*) 'SHIFT PARAMETER (DELTA) =', DELTA
      WRITE (IWW,*) 'MAX EIGEN (ALPHA SQ) =', EIGENMX
      WRITE (IWW,*) 'NE =', NE
      WRITE (IWW,*) 'NNODE =', NNODE
      WRITE (IWW,*) 'MXNEIGEN =', MXNEIGEN
      CLOSE (IWW)
C-----------------------------------------------------
      OPEN (IW, FILE=INFNAME, STATUS='UNKNOWN')
      WRITE (IW,*) DELTA, MXNEIGEN, WSPD
      WRITE (IW,*) NE
      DO I = 1, NE
      WRITE (IW,*) I,(NODEX(I,J),J=1,ND)
      END DO
C
      WRITE (IW,*) NNODE
      DO I = 1, NNODE
      WRITE (IW,*) I, XCOORD(I), YCOORD(I), ZCOORD(I)
      END DO
C
      CLOSE (IW)
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
      CALL BANDWID ( ND, MXE, NE, NODEX , NBW  )
      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 ( MXENGN=',MXNEIGEN,')'
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE BANDWID ( ND, MXE, NE, NODEX , NBW  )
      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