戻る

C ******************************************************************** C
C C
C F 3 D - F E M - R C
C C
C Fluid Analysis Software by Finite Element Method C
C C
C ( Verion : 3.8 ) C
C C
C ( Rectangular Parallelepiped Elements ) C
C C
C Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
PROGRAM MAIN
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
C ***** DATA ( 1 ) ****************************************************
C
** PARAMETER ( NEX = 2, NEY = 2, NEZ = 2,
C
PARAMETER ( NEX = 3, NEY = 3, NEZ = 3,
C
** PARAMETER ( NEX = 4, NEY = 4, NEZ = 4,
** PARAMETER ( NEX = 5, NEY = 5, NEZ = 5,
** PARAMETER ( NEX = 6, NEY = 6, NEZ = 6,
** PARAMETER ( NEX = 10, NEY = 10, NEZ = 10,
** PARAMETER ( NEX = 15, NEY = 15, NEZ = 15,
** PARAMETER ( NEX = 20, NEY = 20, NEZ = 20,
** PARAMETER ( NEX = 30, NEY = 30, NEZ = 30,
** PARAMETER ( NEX = 35, NEY = 35, NEZ = 35,
C
C ( Attention : Memory Check ) 40 x 40 x 40 ==> XXX
C
C *********************************************************************
C
1 LX = NEX + 1, LY = NEY + 1, LZ = NEZ + 1,
2 LZ2 = NEZ* 2 + 1, NNP = LX* LY* LZ,
3 NEL = NEX* NEY* NEZ, NP2 = LX* LY* LZ2,
C
4 NXY = NEX* NEY, NYZ = NEY* NEZ, NZX = NEZ* NEX,
C
5 NM = 13* NEX* NEY* NEZ + 4* ( NXY+NYZ+NZX ) + NEX+NEY+NEZ,
6 NNZ = NNP + NM, NZ2 = NNP + 2* NM )
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), P ( NNP),
1 DLP( NNP), SU( NNP), SV( NNP), SW( NNP),
2 SP ( NNP), U1( NP2), V1( NP2), W1( NNP),
3 X ( NNP), Y ( NNP), Z ( NNP), B ( NNP),
4 F ( NNP), FF( NNP,3)
C
DIMENSION XL( LX ), YL ( LY ), ZL ( LY ), HX( LX+1),
1 HY( LY+1), HZ ( LZ+1), GU ( NNP), GV( NNP),
2 GW( NNP), NC1( NNP), NC2( NNP,8)
C
DIMENSION DLX( NEL), DLY( NEL), DLZ( NEL), DE( NEL),
1 DI ( NNP), DD ( NNZ), PP ( NNZ), D ( NNZ),
2 A ( NNZ), NDE( NEL,8), KCC( NNP),
3 NCA( 2,NNZ), NCB( 2,NZ2)
C
DIMENSION NBP( NNP), NBU( NNP), NBV( NNP), NBW( NNP),
1 LXY( NXY,4), LYZ( NYZ,4), LZX( NZX,4)
C
C ----- Max. Var. of Velocity --- Loop < = 2000 ---
C
COMMON /CON1/ OMG, ERR, NTR
COMMON /CNST/ DT, COF, PHI
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
COMMON /SINF/ ITU, ITV, ITW, ITP, IU2, IV2, IW2, IP2
C
C ***** DATA ( 2 ) ***************************************************
C
** OPEN ( 10, FILE='E2_F3D_FEM_OUT' )
** OPEN ( 11, FILE='E2_F3D_FEM_FIG' )
C
** OPEN ( 10, FILE='E3_F3D_FEM_OUT' )
** OPEN ( 11, FILE='E3_F3D_FEM_FIG' )
C
** OPEN ( 10, FILE='E4_F3D_FEM_OUT' )
** OPEN ( 11, FILE='E4_F3D_FEM_FIG' )
C
** OPEN ( 10, FILE='E6_F3D_FEM_OUT' )
** OPEN ( 11, FILE='E6_F3D_FEM_FIG' )
C
OPEN ( 10, FILE='UE3_F3D_FEM_OUT' )
OPEN ( 11, FILE='UE3_F3D_FEM_FIG' )
C
*** OPEN ( 10, FILE='UE4_F3D_FEM_OUT' )
*** OPEN ( 11, FILE='UE4_F3D_FEM_FIG' )
C
*** OPEN ( 10, FILE='UE10_F3D_FEM_OUT' )
*** OPEN ( 11, FILE='UE10_F3D_FEM_FIG' )
C
*** OPEN ( 10, FILE='UE20_F3D_FEM_OUT' )
*** OPEN ( 11, FILE='UE20_F3D_FEM_FIG' )
C
*** OPEN ( 10, FILE='UE30_F3D_FEM_R1000_OUT' )
*** OPEN ( 11, FILE='UE30_F3D_FEM_R1000_FIG' )
C
C ********************************************************************
C
CALL GMSH ( XL, YL, ZL, HX, HY, HZ, LX, LY, LZ, MES )
C
CALL IPSB ( U, V, W, P, DLP, SP, X, Y, Z,
1 XL, YL, ZL, NDE, NBU, NBV, NBW, NBP, NCA,
2 NC1, NC2, NCB, DD, PP, KCC, LX, LY, LZ,
3 NNP, NEL, NNZ, NZ2 )
C
LOP = 0
C --------------------------------------------------------------------
999 LOP = LOP + 1
C
C ----- Corrected Diffusivity etc ---
C
CALL CKAK ( U, V, W, GU, GV, GW, XL, YL, ZL,
1 HX, HY, HZ, FF, LX, LY, LZ, NNP )
C
OMG = 1.1D0
C
C ----- Pseudo-Velocities : U, V, W ---
C
CALL CMPA ( U, V, W, U, X, Y, Z, NNP, NDE,
1 NEL, P, PP, B, A, D, F, 1, U1,
2 NNZ, NZ2, NBU, NCA, NCB, KCC, ITU )
C
CALL CMPA ( U, V, W, V, X, Y, Z, NNP, NDE,
1 NEL, P, PP, B, A, D, F, 2, V1,
2 NNZ, NZ2, NBV, NCA, NCB, KCC, ITV )
C
CALL CMPA ( U, V, W, W, X, Y, Z, NNP, NDE,
1 NEL, P, PP, B, A, D, F, 3, W1,
2 NNZ, NZ2, NBW, NCA, NCB, KCC, ITW )
C
OMG = 1.8D0
C ================
C
C ----- Guessed Pressure ---
C
CALL DPRS ( P, U1, V1, W1, X, Y, Z, NDE, NCB,
1 KCC, DD, B, NBP, NNP, NEL, NNZ, NZ2, ITP )
C
OMG = 0.8D0
C ================
C
C ----- Guessed Velocities : U2, V2, W2 ---
C
CALL CMPB ( U, V, W, GU, GV, GW, U, X, Y, Z,
1 NDE, P, PP, B, A, D, F, FF, 1, SU,
2 NNP, NEL, NNZ, NZ2, NBU, NCA, NCB, KCC, IU2 )
C
CALL CMPB ( U, V, W, GU, GV, GW, V, X, Y, Z,
1 NDE, P, PP, B, A, D, F, FF, 2, SV,
2 NNP, NEL, NNZ, NZ2, NBV, NCA, NCB, KCC, IV2 )
C
CALL CMPB ( U, V, W, GU, GV, GW, W, X, Y, Z,
1 NDE, P, PP, B, A, D, F, FF, 3, SW,
2 NNP, NEL, NNZ, NZ2, NBW, NCA, NCB, KCC, IW2 )
C
OMG = 1.6D0
C
C ----- Pressure Corrrection ---
C
CALL DPRS ( DLP, SU, SV, SW, X, Y, Z, NDE, NCB,
1 KCC, DD, B, NBP, NNP, NEL, NNZ, NZ2, IP2 )
C
C ----- Pressure & Velocities ---
C
CALL VLPR ( P, DLP, DE, DI, DLX, DLY, DLZ, NDE, NC2,
1 NC1, SU, SV, SW, SP, X, Y, Z, NBU,
2 NBV, NBW, NNP, NEL)
C
C ----- Max. Relative Variance for Velocities ---
C
CALL CVAR ( U, V, W, P, SU, SV, SW, SP, NNP )
C
C ----- Outout ---
C
CALL OUTP ( SU, SV, SW, P, X, Y, Z, XL, YL,
1 ZL, NDE, LX, LY, LZ, NNP, NEL, *999 )
C
C ----- Output for Plotting ---
C
CALL PLTO ( SU, SV, SW, P, X, Y, Z, XL, YL,
1 ZL, LXY, LYZ, LZX, LX, LY, LZ, NNP, NXY,
2 NYZ, NZX, MES )
C
C ----- Closing Output ---
C
CALL UVCK ( U, V, X, Y, XL, LX, LY, LZ, NNP )
C
CLOSE ( 10)
CLOSE ( 11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE GMSH ( XL, YL, ZL, HX, HY, HZ, LX, LY, LZ, MES )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XL( LX), YL ( LY), ZL( LZ), HX( LX+1), HY( LY+1),
1 HZ( LZ+1)
C
COMMON /ZDFC/ ZFC
C
C -----------------
C
CALL TTLE
C
C -----------------
C
C ----- Scaling Factor for Z-Direction ---
C
ZFC = 1.D0
C ================
C
C ------------------------------------------------
C
C MES = 1 : Meshes of equal length ----- EM
C
C MES = 2 : JSME-Paper Version --------- UEM
C
C ------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** MES = 1 ( Meshes of equal length ) or 2 ( Meshes o
1f unequal length ) ?'
*** READ(5,*) MES
MES = 2
WRITE(6,2000) MES
2000 FORMAT(/,' ** MES =',I3)
WRITE(6,*) ' '
C
IF ( MES.EQ.1) GO TO 333
IF ( MES.EQ.2) GO TO 555
C
333 CONTINUE
C -------------------------------------------
C
WRITE(6,*) '* Meshes of equal length '
C
C -------------------------------------------
DO 100 I = 1, LX
C
XL( I) = DFLOAT( I-1)/ DFLOAT( LX-1)
100 CONTINUE
C --------------------------------------------
DO 200 J = 1, LY
C
YL( J) = DFLOAT( J-1)/ DFLOAT( LY-1)
200 CONTINUE
C --------------------------------------------
DO 300 K = 1, LZ
C
ZL( K) = DFLOAT( K-1)/ DFLOAT( LZ-1)
300 CONTINUE
C --------------------------------------------
GO TO 999
C
555 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,*) '*** Meshes of unequal length ( UEM ) ***'
C
C ----------------------------------------------------------
C
C ----- JSME-Paper Ref.---
C
AA = 2.D0
C
DX = 2.D0/ DFLOAT( LX-1)
DY = 2.D0/ DFLOAT( LY-1)
DZ = 2.D0/ DFLOAT( LZ-1)
C
C ----- Important ---
C
MDX = IDNINT( DFLOAT( LX+1)/ 2.D0 )
MX1 = MDX - 1
C
MDY = IDNINT( DFLOAT( LY+1)/ 2.D0 )
MY1 = MDY - 1
C
MDZ = IDNINT( DFLOAT( LZ+1)/ 2.D0 )
MZ1 = MDZ - 1
C
C -------------------
C
XL( MDX) = 0.5D0
YL( MDY) = 0.5D0
ZL( MDZ) = 0.5D0
C
C ----------------------------------------------------
DO 400 I = 1, MX1
C
XL( I) = ( DEXP( AA* DFLOAT( I-1)* DX) - 1.D0 )
1 /( 2.D0* ( DEXP( AA) - 1.D0 ))
XL( LX+1-I ) = 1.D0 - XL ( I)
C
400 CONTINUE
C ----------------------------------------------------
DO 500J = 1, MY1
C
YL( J) = ( DEXP( AA* DFLOAT( J-1)* DY) - 1.D0 )
1 /( 2.D0* ( DEXP( AA) - 1.D0 ))
YL( LY+1-J) = 1.D0 - YL( J)
C
500 CONTINUE
C ----------------------------------------------------
DO 600 K = 1, MZ1
C
ZL( K) = ( DEXP( AA* DFLOAT( K-1)* DZ) - 1.D0 )
1 /( 2.D0* ( DEXP( AA) - 1.D0 ))
ZL( LZ+1-K) = 1.D0 - ZL( K)
C
600 CONTINUE
C -----------------------------------------------------
C
999 CONTINUE
C -----------------------------------------------------
DO 700 I = 1, LX-1
C
HX( I) = DABS( XL( I+1) - XL( I))
700 CONTINUE
C -----------------------------------------------------
DO 800 J = 1, LY-1
C
HY( J) = DABS( YL( J+1) - YL( J))
800 CONTINUE
C -----------------------------------------------------
DO 900 J = 1, LZ-1
C
HZ( J) = DABS( ZL( J+1) - ZL( J))
900 CONTINUE
C ------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMPB ( U, V, W, GU, GV, GW, C, X,
1 Y, Z, NDE, P, PP, B, A, D,
2 F, FF, ID, CW, NNP, NEL, NNZ, NZ2,
3 NBC, NCA, NCB, KCC, NTL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), GU( NNP),
1 GV( NNP), GW( NNP), C ( NNP), X ( NNP),
2 Y ( NNP), Z ( NNP), NDE( NEL,8), P ( NNP),
3 PP( NNZ), B ( NNP), A ( NNZ), D ( NNZ),
4 F ( NNP), FF( NNP,3), CW ( NNP)
C
DIMENSION XX ( 8), YY ( 8), ZZ ( 8), EXX( 8,8),
1 EYY( 8,8), EZZ( 8,8), EPP( 8,8), EFX( 8,8),
2 EFY( 8,8), EFZ( 8,8), PE ( 8)
C
DIMENSION KCC( NNP), NBC( NNP), NCA( 2, NNZ), NCB( 2,NZ2)
C
COMMON /CNST/ DT, COF, PHI
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
C
C -------------------------------------
DO 100 I = 1, NNP
C
F ( I) = 0.D0
B ( I) = 0.D0
CW( I) = C( I)
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, NNZ
C
D( I) = 0.D0
200 CONTINUE
C --------------------------------------
C
FSX = 0.D0
FSY = 0.D0
FSZ = 0.D0
FH = 0.D0
C
IF ( ID.EQ.1) FSX = 1.D0
IF ( ID.EQ.2) FSY = 1.D0
IF ( ID.EQ.2) FH = 1.D0
IF ( ID.EQ.3) FSZ = 1.D0
C
C -------------------------------------
DO 300 INEL = 1, NEL
C
UE = 0.D0
VE = 0.D0
WE = 0.D0
C
FF1 = 0.D0
FF2 = 0.D0
FF3 = 0.D0
C
GUE = 0.D0
GVE = 0.D0
GWE = 0.D0
C
C ----------------------------
DO 350 I = 1, 8
C
N = NDE( INEL,I)
XX( I) = X( N)
YY( I) = Y( N)
ZZ( I) = Z( N)
C
PE( I) = P( N)
C
UE = UE + U( N)/ 8.D0
VE = VE + V( N)/ 8.D0
WE = WE + W( N)/ 8.D0
C
FF1 = FF1 + FF( N,1)/ 8.D0
FF2 = FF2 + FF( N,2)/ 8.D0
FF3 = FF3 + FF( N,3)/ 8.D0
C
GUE = GUE + GU( N)/ 8.D0
GVE = GVE + GV( N)/ 8.D0
GWE = GWE + GW( N)/ 8.D0
C
350 CONTINUE
C ------------------------------------------------------------------
C
CALL JCPD ( XX, YY, ZZ, EXX, EYY, EZZ, EFX, EFY, EFZ, EPP )
C
C ------------------------------------------------------------------
DO 300 I = 1, 8
DO 300 J = 1, 8
C --------------------------------------
C
II = NDE( INEL,I)
JJ = NDE( INEL,J)
C
DS = COF* ( FF1* EXX( I,J) + FF2* EYY( I,J) + FF3* EZZ( I,J))
C
C ----------------------------------------------
C
CALL DCON ( DS, II, JJ, D, NCA, NNZ, L )
C
C ----------------------------------------------
F( II) = F( II)
1 + ( GUE* EFX( I,J) + GVE* EFY( I,J) + GWE* EFZ( I,J))* C( JJ)
C
2 + FSX* EFX( I,J)* PE( J)
3 + FSY* EFY( I,J)* PE( J)
4 + FSZ* EFZ( I,J)* PE( J)
C
300 CONTINUE
C --------------------------------------
JS = 1
IF ( ISL.EQ.0) GO TO 999
C
C ***** Implicit Method ***
C
C --------------------------------------
DO 400 I = 1, NNP
C
JL = KCC( I)
C
C ----------------------------
DO 450 J = JS, JL
C
KR = NCB( 1,J)
KL = NCB( 2,J)
C
A( KR) = PP( KR)/ DT + PHI* D( KR)
B( I) = B( I) + ( PP( KR)/ DT - ( 1.D0-PHI)* D( KR))* C( KL)
C
450 CONTINUE
C ----------------------------
C
B( I)= B( I) - F( I)
JS = JL + 1
C
400 CONTINUE
C --------------------------------------
C
CALL DSLVE ( NNP, NNZ, NZ2, B, CW, A, NCB, KCC, NBC, NTL )
C
C --------------------------------------
GO TO 990
C
C ***** Explicit Method ***
C
999 CONTINUE
C --------------------------------------
DO 500 I = 1, NNP
C
JL = KCC( I)
FW = 0.D0
C
C ---------------------------
DO 550 J = JS, JL
C
KR = NCB( 1,J)
KL = NCB( 2,J)
C
FW = FW + PP( KR)
B( I) = B( I) + ( PP( KR)/ DT - D( KR))* C( KL)
C
550 CONTINUE
C ----------------------------
IF ( NBC( I).NE.9999 ) CW( I) = ( B( I) - F( I))/ FW* DT
JS = JL + 1
C
500 CONTINUE
C --------------------------------------
C
990 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMPA ( U, V, W, C, X, Y, Z, NNP,
1 NDE, NEL, P, PP, B, A, D, F,
2 ID, CW, NNZ, NZ2, NBC, NCA, NCB, KCC,
3 NTL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION X ( NNP), Y( NNP), Z ( NNP), U( NNP), V ( NNP),
1 W ( NNP), P( NNP), F ( NNP), B( NNP), C ( NNP),
2 CW( NNP), A( NNZ), PP( NNZ), D( NNZ), NDE( NEL,8)
C
DIMENSION XX ( 8), YY ( 8), ZZ ( 8), EXX( 8,8),
1 EYY( 8,8), EZZ( 8,8), EPP( 8,8), EFX( 8,8),
2 EFY( 8,8), EFZ( 8,8), PE ( 8)
C
DIMENSION KCC( NNP), NBC( NNP), NCA( 2,NNZ), NCB( 2,NZ2)
C
COMMON /CNST/ DT, COF, PHI
C
C ----------------------------
DO 100 I = 1, NNP
C
F ( I) = 0.D0
B ( I) = 0.D0
CW( I) = C( I)
C
100 CONTINUE
C ----------------------------
DO 200 I = 1, NNZ
C
D( I) = 0.D0
200 CONTINUE
C ----------------------------
FH = 0.D0
C
IF ( ID.EQ.2) FH = 1.D0
C ----------------------------------------------------
DO 300 INEL = 1, NEL
C
UE = 0.D0
VE = 0.D0
WE = 0.D0
C
C ----------------------------
DO 350 I = 1, 8
C
N = NDE( INEL,I)
XX( I) = X( N)
YY( I) = Y( N)
ZZ( I) = Z( N)
PE( I) = P( N)
C
UE = UE + U( N)/ 8.D0
VE = VE + V( N)/ 8.D0
WE = WE + W( N)/ 8.D0
C
350 CONTINUE
C ---------------------------
C
CALL JCPD ( XX, YY, ZZ, EXX, EYY, EZZ, EFX, EFY, EFZ, EPP )
C
C ---------------------------------
DO 300 I = 1, 8
DO 300 J = 1, 8
C ---------------------------------
C
II = NDE( INEL,I)
JJ = NDE( INEL,J)
DS = COF* ( EXX( I,J) + EYY( I,J) + EZZ( I,J))
C
C ----------------------------------------------
C
CALL DCON ( DS, II, JJ, D, NCA, NNZ, L )
C
C ----------------------------------------------
C
F( II) = F( II)
1 + ( UE* EFX( I,J) + VE* EFY( I,J) + WE* EFZ( I,J))* C( JJ)
C
300 CONTINUE
C ---------------------------------
JS = 1
C
C ---------------------------------
DO 400 I = 1, NNP
C
JL = KCC( I)
C
C ------------------------
DO 450 J = JS, JL
C
KR = NCB( 1,J)
KL = NCB( 2,J)
C
A( KR) = PP( KR)/ DT
B( I) = B( I) + ( PP( KR)/ DT - D( KR))* C( KL)
C
450 CONTINUE
C ------------------------
C
B( I)= B( I) - F( I)
JS = JL + 1
C
400 CONTINUE
C ---------------------------------------------------------------
C
CALL DSLVE ( NNP, NNZ, NZ2, B, CW, A, NCB, KCC, NBC, NTL )
C
C ---------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CVAR ( U, V, W, P, SU, SV, SW, SP, NNP )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), P( NNP), SU( NNP),
1 SV( NNP), SW( NNP), SP( NNP)
C
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
C
EPS = 0.D0
UMM = 0.D0
C
EPX = 0.D0
UMX = 0.D0
C
ESP = 0.D0
PMX = 0.D0
C
C ------------------------------------------
DO 100 I = 1, NNP
C
UMM = DMAX1( UMM, DABS( SU( I)))
UMM = DMAX1( UMM, DABS( SV( I)))
UMM = DMAX1( UMM, DABS( SW( I)))
C
UMX = DMAX1( UMX, DSQRT( SU( I)** 2 + SV( I)** 2 + SW( I)** 2 ))
PMX = DMAX1( PMX, DABS ( SP( I)))
C
100 CONTINUE
C ------------------------------------------
C
IF ( LOP.EQ.1) SPX = PMX
C
C ------------------------------------------
DO 200 I = 1, NNP
C
UN = SU( I)
VN = SV( I)
WN = SW( I)
PN = SP( I)
C
IF ( DABS( UN).GT.UMM* 0.01D0 )
1 EPS = DMAX1( EPS, DABS(( U( I) - UN)/ UN))
C
IF ( DABS( VN).GT.UMM* 0.01D0 )
1 EPS = DMAX1( EPS, DABS(( V( I) - VN)/ VN))
C
IF ( DABS( WN).GT.UMM* 0.01D0 )
1 EPS = DMAX1( EPS, DABS(( W( I) - WN)/ WN))
C
C ---------------------------------------------
C
IF ( DABS( PN).GT.PMX* 0.01D0 )
1 ESP = DMAX1( ESP, DABS(( PMX - SPX)/ SPX))
C
UXN = DSQRT( SU( I)** 2 + SV( I)** 2 + SW( I)** 2 )
UX = DSQRT( U ( I)** 2 + V ( I)** 2 + W ( I)** 2 )
C
IF ( UXN.GT.UMX* 0.01D0)
1 EPX = DMAX1( EPX, DABS(( UX - UXN)/ UXN))
C
U( I) = UN
V( I) = VN
W( I) = WN
C
200 CONTINUE
C -------------------------------------------
C
C ----- (%) --- Computation for each LOP ---
C
VVR( LOP) = EPX* 100.D0
PVR( LOP) = ESP* 100.D0
C
SPX = PMX
C ==============
C
C ---------------------------
DO 300 I = 1, NNP
C
SP( I) = P( I)
C
300 CONTINUE
C ---------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE UVCK ( U, V, X, Y, XL, LX, LY, LZ, NNP )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /CNST/ DT, COF, PHI
C
DIMENSION U( NNP), V( NNP), X( NNP), Y( NNP), XL( LX)
C
L1 = LX* LY* ( LZ - 1)/ 2
L2 = ( LX - 1)/ 2 + 1
L3 = ( LY - 1)/ 2
C
U2MX = 0.D0
V2MX = 0.D0
C
C --------------------------------------
DO 100 I = 1, LY
C
LL = L1 + L2 + LX* ( I-1)
U2MX = DMAX1( U2MX, U( LL))
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, LX
C
LL = L1 + LX* L3 + I
V2MX = DMAX1( V2MX, V( LL))
C
200 CONTINUE
C --------------------------------------
C
UMX = 0.D0
VMX = 0.D0
C --------------------------------------
DO 300 I = 1, NNP
C
UMX = DMAX1( UMX, U( I))
VMX = DMAX1( VMX, V( I))
C
300 CONTINUE
C --------------------------------------
C
WRITE(6,2000) UMX, VMX
2000 FORMAT(/,' * Umax. =', F8.3,5X,' Vmax. =',F8.3)
WRITE(10,2000) UMX, VMX
C
C ----- Attention -----------
C
HHX = XL( 2) - XL( 1)
C ---------------------------
C
BXM = UMX* DT/ HHX
BYM = VMX* DT/ HHX
C
WRITE(6,2100) BXM, BYM
2100 FORMAT(' * bx_max.=', F8.3,5X,' by_max.=', F8.3)
WRITE(10,2100) BXM, BYM
WRITE(6, *) ' '
WRITE(6, *) '*** End of Computation ***'
WRITE(10,*) ' '
WRITE(10,*) '*** End of Computation ***'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCON ( VAL, I,J, RMT, NCA, NNZ, IC )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION RMT( NNZ), NCA( 2, NNZ)
C
IC = 0
C
IF ( I.GT.J) RETURN
C
C ----- Found I,J & store valkue with additional old value ---
C
N1 = 1
N2 = NNZ
99 IF (( N1 + N2)/ 2.EQ.IC) N2 = N2 + 1
IC = ( N1 + N2)/ 2
C
C ----------------------------
IF ( NCA( 1,IC).LT.I.OR.
1 ( NCA( 1,IC).EQ.I.AND.NCA( 2,IC).LT.J)) THEN
C
N1 = IC
GO TO 99
ELSE IF ( NCA( 1,IC).GT.I.OR.
1 ( NCA( 1,IC).EQ.I.AND.NCA( 2,IC).GT.J)) THEN
N2 = IC
GO TO 99
C
END IF
C ---------------------------
C
RMT( IC) = RMT( IC) + VAL
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCNL ( PDV, NCA, I, J, PD, K, NNZ )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION PD( NNZ), NCA( 2,NNZ)
C
IF ( I.GT.J) GO TO 99
IE = K
N1 = 1
IF ( K.EQ.0) N2 = K + 1
IF ( K.EQ.0) GO TO 10
N2 = K
10 N = N2
IF ( IE.EQ.0) IE = 1
C
C ----------------------------------------------------------
DO 100 II = 1, IE
C
IF ( K.EQ.0) GO TO 43
IF ( I.LT.NCA( 1,N).OR.( I.EQ.NCA( 1,N).AND.J.LT.NCA( 2,N)))
1 GO TO 21
IF ( I.EQ.NCA( 1,N).AND.J.EQ.NCA( 2,N)) GO TO 22
GO TO 23
C
21 CONTINUE
N2 = N
N = N2 - ( N2 - N1)/ 2
IF ( N.NE.N2) GO TO 100
IF ( I.LT.NCA( 1,N1).OR.( I.EQ.NCA( 1,N1)
1 .AND.J.LT.NCA( 2,N1))) GO TO 31
IF ( I.EQ.NCA( 1,N1).AND.J.EQ.NCA( 2,N1)) GO TO 32
GO TO 33
C
31 N = N1
33 CONTINUE
C -------------------------------------------------
C
CALL SUBS ( PDV, I, J, N, PD, NCA, K, NNZ )
C
C -------------------------------------------------
GO TO 99
C
32 N = N1
GO TO 22
C
42 N = N2
22 CONTINUE
PD( N) = PD( N) + PDV
GO TO 99
C
23 CONTINUE
N1 = N
N = N1 + ( N2 - N1)/ 2
IF ( N.NE.N1) GO TO 100
IF ( I.LT.NCA( 1,N2).OR.( I.EQ.NCA( 1,N2)
1 .AND. J.LT.NCA( 2,N2))) GO TO 41
IF ( I.EQ.NCA( 1,N2).AND. J.EQ.NCA( 2,N2)) GO TO 42
GO TO 43
C
41 N = N2
GO TO 33
C
43 K = K + 1
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' ** Error ** ''NNZ'' in ''SOR'' is too small.**')
PD ( K) = PDV
NCA( 1,K) = I
NCA( 2,K) = J
GO TO 99
C
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error --- DO Loop Strange, K =',I10)
STOP
C
99 CONTINUE
RETURN
END
C **********************************************************************
C
SUBROUTINE DCIJ ( KCC, NCA, NZ, NNP, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION KCC( NNP), NCA( 2,NZ2)
C
C NP, NZ : Current length of KCC, NCA
C NNP, NZ2 : Max.length of KCC, NCA
C
CHARACTER*3 ANNP, ANZ2
C
DATA ANNP /'NNP'/, ANZ2 /'NZ2'/
C
NP = 1
C
C --------------------------------------
DO 100 I = 1, NZ
C
IF ( NP.EQ.NCA( 1,I)) GO TO 50
C
C ----- To set KCC ( NP) ---
C
KCC( NP) = I-1
NP = NCA( 1,I)
IF ( NP.GT.NNP ) GO TO 30
50 NCA( 1,I) = I
C
100 CONTINUE
C --------------------------------------
C
C ----- To set Last row's counter ---
C
KCC( NP) = NZ
C
C ----- To set ( I,J) in NCA for I > J ---
C
C ----------------------------------------------------------
DO 200 II = 2, NP
C
IM1 = II-1
IS = KCC( IM1)
C
C ---------------------------------------------------
DO 220 J = 1, IM1
C
JS = 1
IF ( J.GT.1) JS = KCC( J-1) + 1
JL = KCC( J)
C
C ----------------------------------------
DO 240 K = JS, JL
C
IF ( NCA( 2,K).NE.II ) GO TO 240
K1 = K
GO TO 20
C
240 CONTINUE
C ----------------------------------------
C
C ----- ( I,J) Not exist ---
C
GO TO 220
C
20 CONTINUE
C
C ------ To shift NCA ---
C
NZIS = NZ - IS
C ----------------------------------------
DO 260 L = 1, NZIS
C
NCA( 1,NZ-L+2) = NCA( 1,NZ-L+1)
NCA( 2,NZ-L+2) = NCA( 2,NZ-L+1)
C
260 CONTINUE
C ----------------------------------------
C
NZ = NZ + 1
C *** Attention **
C
IF ( NZ.GT.NZ2) GO TO 40
C
C ----- To increase KCC ---
C
C ----------------------------
DO 280 L = II, NP
C
KCC( L) = KCC( L)+ 1
C
280 CONTINUE
C ----------------------------
C
C.----- To insert ( I,J) element ---
C
IS = IS + 1
NCA( 1,IS) = NCA( 1,K1)
NCA( 2,IS) = J
C
220 CONTINUE
C ---------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
RETURN
C
30 WRITE(6,2000) ANNP, NNP, ANNP, NP
STOP
C
40 WRITE(6,2000) ANZ2, NZ2, ANZ2, NZ
2000 FORMAT(' ** Error ( DCIJ) ',A3,'(',I5,') too small',A2,'=',I5 )
STOP
END
C **********************************************************************
C
SUBROUTINE DPRS ( DLP, U, V, W, X, Y, Z, NDE, NCB,
1 KCC, DD, B, NBP, NNP, NEL, NNZ, NZ2, NTL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION DLP( NNP), U ( NNP), V ( NNP), W ( NNP),
1 X ( NNP), Y ( NNP), Z ( NNP), NDE( NEL,8),
2 NCB( 2,NZ2), KCC( NNP), DD( NNZ), B ( NNP),
3 NBP( NNP)
C
DIMENSION XX( 8), YY( 8), ZZ( 8), UE( 8), VE( 8), WE( 8)
C
COMMON /CNST/ DT, COF, PHI
COMMON /CON1/ OMG, ERR, NTR
C
C ----- Generation for Pressure Poisson Equation ---
C
DO 100 I = 1, NNP
C
B( I) = 0.D0
C
100 CONTINUE
C
C ------------------------------------------------
DO 200 INEL = 1, NEL
C
C ----------------------------
DO 220 I = 1, 8
C
N = NDE( INEL,I)
C
XX( I) = X( N)
YY( I) = Y( N)
ZZ( I) = Z( N)
C
UE( I) = U( N)
VE( I) = V( N)
WE( I) = W( N)
C
220 CONTINUE
C ----------------------------
C
C --- Global matrix ---
C
C ----------------------------
DO 240 J = 1, 8
C
JJ = NDE ( INEL,J)
C
UX = DVLX( X, U, NDE, NNP, NEL, INEL )
VY = DVLY( Y, V, NDE, NNP, NEL, INEL )
WZ = DVLZ( Z, W, NDE, NNP, NEL, INEL )
C
VV = DABS((XX( 7)-XX( 1))* ( YY( 7)-YY( 1))* ( ZZ( 7)-ZZ( 1)))
C
B( JJ) = B( JJ) - VV/ 8.D0* ( UX + VY + WZ )/ DT
C
240 CONTINUE
C ---------------------------
C
200 CONTINUE
C ------------------------------------------------
C
C ------------------------------------------------------------------
C
CALL SOLVP ( NNP, NNZ, NZ2, B, DLP, DD, NCB, KCC, NBP, NTL )
C
C ------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION DVLX ( X, C, NDE, NNP, NEL, I )
C
C ----- Differentiated values ( 3D ) ---
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION X( NNP), C( NNP), NDE( NEL,8)
C
HX = DABS( X( NDE( I,7)) - X( NDE( I,1)))
C
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
C5 = C( NDE( I,5))
C6 = C( NDE( I,6))
C7 = C( NDE( I,7))
C8 = C( NDE( I,8))
C
DVLX = ( - C1 + C2 + C3 - C4 - C5 + C6 + C7 - C8 )/( 4.D0* HX )
C
RETURN
END
C **********************************************************************
C
FUNCTION DVLY ( Y, C, NDE, NNP, NEL, I )
C
C ----- Differentiated values ( 3D ) ---
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION Y( NNP), C( NNP), NDE( NEL,8)
C
HY = DABS( Y( NDE( I,7)) - Y( NDE( I,1)))
C
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
C5 = C( NDE( I,5))
C6 = C( NDE( I,6))
C7 = C( NDE( I,7))
C8 = C( NDE( I,8))
C
DVLY = ( - C1 - C2 + C3 + C4 - C5 - C6 + C7 + C8 )/( 4.D0* HY )
C
RETURN
END
C **********************************************************************
C
FUNCTION DVLZ ( Z, C, NDE, NNP, NEL, I )
C
C ----- Differentiated values ( 3D ) ---
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION Z( NNP), C( NNP), NDE( NEL,8)
C
HZ = DABS( Z( NDE( I,7)) - Z( NDE( I,1)))
C
C1 = C( NDE( I,1))
C2 = C( NDE( I,2))
C3 = C( NDE( I,3))
C4 = C( NDE( I,4))
C
C5 = C( NDE( I,5))
C6 = C( NDE( I,6))
C7 = C( NDE( I,7))
C8 = C( NDE( I,8))
C
DVLZ = ( - C1 - C2 - C3 - C4 + C5 + C6 + C7 + C8 )/( 4.D0* HZ )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C F 3 D - F E M - R C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( V.3.8 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 3D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Rectangular Parallelepiped Elements ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE IPSB ( U, V, W, P, DLP, SP, X, Y,
1 Z, XL, YL, ZL, NDE, NBU, NBV, NBW,
2 NBP, NCA, NC1, NC2, NCB, DD, PP, KCC,
3 LX, LY, LZ, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), P ( NNP),
1 DLP( NNP), SP ( NNP), X ( NNP), Y ( NNP),
2 Z ( NNP), NDE( NEL,8), NBU( NNP), NBV( NNP),
3 NBW( NNP), NBP( NNP), NCA( 2,NNZ), NC1( NNP),
4 NC2( NNP), NCB( 2,NZ2), DD ( NNZ), PP ( NNZ),
5 KCC( NNP)
C
DIMENSION XL( LX), YL( LY), ZL( LZ)
C
DIMENSION EPP( 8,8), EDD( 8,8), XX( 8), YY( 8), ZZ( 8)
C
COMMON /CON1/ OMG, ERR, NTR
COMMON /CNST/ DT, COF, PHI
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
C
CALL DINT ( P, DLP, SP, NBP, KCC, NCA, NCB, LX, LY, LZ,
1 NNP, NEL, NNZ, NZ2 )
C
C ----------------------------------------------------------------------
C
CALL IPGN ( X, Y, Z, U, V, W, NDE, NBU, NBV,
1 NBW, XL, YL, ZL, LX, LY, LZ, NEL, NNP,
2 NNZ, NZ2 )
C
C ----------------------------------------------------------------------
C
CALL STNC ( NCA, NNZ, NDE, NEL, 8, PP )
C
C ----------------------------------------------------------------------
DO 100 I = 1, NNZ
C
NCB( 1,I) = NCA( 1,I)
NCB( 2,I) = NCA( 2,I)
C
100 CONTINUE
C ----------------------------
C
NZ = NNZ
C
WRITE(6,*) ' ** In Processing **'
** WRITE(6,*) ' '
C ---------------------------------------------
C
CALL DCIJ ( KCC, NCB, NZ, NNP, NZ2 )
C
C ---------------------------------------------
WRITE(6,*) ' ** Computation started **'
WRITE(6,*) ' '
C
C ---------------------------
C
DO 200 I = 1, NNZ
C
DD( I) = 0.D0
PP( I) = 0.D0
C
200 CONTINUE
C
C -----------------------------------------------
DO 300 INEL = 1, NEL
C
C ----------------------------
DO 350 I = 1, 8
C
N = NDE( INEL,I)
XX( I) = X( N)
YY( I) = Y( N)
ZZ( I) = Z( N)
C
350 CONTINUE
C ----------------------------
C
CALL JACP ( XX, YY, ZZ, EPP, EDD )
C
C --------------------------------------
DO 370 I = 1, 8
DO 370 J = 1, 8
C --------------------------------------
C
II = NDE( INEL,I)
JJ = NDE( INEL,J)
C ----------------------------------------------------------
C
CALL DCON ( EDD( I,J), II, JJ, DD, NCA, NNZ, L )
C
C ----------------------------------------------------------
IF ( L.NE.0) PP( L) = PP( L) + EPP( I,J)
C
370 CONTINUE
C --------------------------------------
C
300 CONTINUE
C -------------------------------------------
C
CALL STNN ( NDE, NC1, NC2, NNP, NEL )
C
C -------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DINT ( P, DLP, SP, NBP, KCC, NCA, NCB, LX, LY,
1 LZ, NNP, NEL, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /ZDFC/ ZFC
COMMON /CNST/ DT, COF, PHI
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
C
DIMENSION P ( NNP), DLP( NNP), SP ( NNP), NBP( NNP),
1 KCC( NNP), NCA( 2,NNZ), NCB( 2,NZ2)
C
C --------------------------------------
WRITE(6,2000) ZFC
2000 FORMAT(/,' *** ZFC =',F6.3)
WRITE(10,2000) ZFC
WRITE(6,2050)
2050 FORMAT(/,' ** Re = ?')
*** READ(5,*) RE
C
RE = 100.D0
COF = 1.D0 / RE
C --------------------------------------
C
WRITE(6,2100)
2100 FORMAT(' ** Time Step ? ')
*** READ(5,*) DT
DT = 0.3D0
WRITE(6,2200)
2200 FORMAT(' ** How many Loops ? & NPR = ?')
C
*** READ(5,*) LPE, NPR
C
LPE = 20
NPR = 10
C --------------------------------------
C
WRITE(6,2300)
2300 FORMAT(' ** Implicit Method = 1 * Explicit Method = 0 ')
WRITE(6,*) ' '
*** READ(5,*) ISL
C
ISL = 1
C
EPL = 0.001D0
C
WRITE(6,2400) RE, DT, LPE, ISL, NPR
2400 FORMAT(' *** Re =',F7.1,' DT =',F7.4,' LPE =',I5,' *ISL =',
1 I2,' NPR =',I3 )
WRITE(10,2400) RE, DT, LPE, ISL, NPR
WRITE(6,2500) LX-1, LY-1, LZ-1, NNP, NEL, NNZ, NZ2
2500 FORMAT(' * NEX =',I3,' NEY =',I3,' NEZ =',I3,' NNP=',I5,' NEL=',
1 I5,' NNZ =',I7,' NZ2=',I7 )
WRITE(10,2500) LX-1, LY-1, LZ-1, NNP, NEL, NNZ, NZ2
WRITE(6,*) ' '
C
C ------ Initialization ----------------
DO 100 I = 1, NNP
C
KCC( I) = 0
DLP( I) = 0.D0
NBP( I) = 0
P ( I) = 0.D0
SP ( I) = 1.D0
C
100 CONTINUE
C -------------------------------------
DO 200 I = 1, 2
DO 200 J = 1, NNZ
C
NCA( I,J) = 0
C
200 CONTINUE
C -------------------------------------
DO 300 I = 1, 2
DO 300 J = 1, NZ2
C
NCB( I,J) = 0
C
300 CONTINUE
C -------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE IPGN ( X, Y, Z, U, V, W, NDE, NBU,
1 NBV, NBW, XL, YL, ZL, LX, LY, LZ,
2 NEL, NNP, NNZ, NZ2 )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), Z ( NNP), U ( NNP),
1 V ( NNP), W ( NNP), NDE( NEL,8), NBU( NNP),
2 NBV( NNP), NBW( NNP), XL ( LX), YL ( LY),
3 ZL ( LZ)
C
COMMON /CNST/ DT, COF, PHI
COMMON /CON1/ OMG, ERR, NTR
C
C ----- Attention ---
C
PHI = 1.D0
C ================
C
ERR = 0.0001D0
NTR = 200
C
C ----------------------------
DO 100 I = 1, NNP
C
U( I) = 0.D0
V( I) = 0.D0
W( I) = 0.D0
C
100 CONTINUE
C ----------------------------
C
J = 1
C ----------------------------
DO 200 IZ = 1, LZ
DO 200 IY = 1, LY
DO 200 IX = 1, LX
C
X( J) = XL( IX)
Y( J) = YL( IY)
Z( J) = ZL( IZ)
J = J + 1
C
200 CONTINUE
C ----------------------------
C
J = 1
C ----------------------------
DO 300 IZ = 1, LZ-1
DO 300 IY = 1, LY-1
DO 300 IX = 1, LX-1
C ----------------------------
C
ST= IX + ( IY-1)* LX + ( IZ-1)* LX* LY
C
NDE( J,1) = ST
NDE( J,2) = ST + 1
NDE( J,3) = ST + LX + 1
NDE( J,4) = ST + LX
C
NDE( J,5) = ST + LX* LY
NDE( J,6) = ST + LX* LY + 1
NDE( J,7) = ST + LX* ( LY+1) + 1
NDE( J,8) = ST + LX* ( LY+1)
C
J = J + 1
C
300 CONTINUE
C ----------------------------
J = 1
C
C ----- Boundary conditions -----------
C
DO 400 IZ = 1, LZ
DO 400 IY = 1, LY
DO 400 IX = 1, LX
C -------------------------------------
C
IF ( IY.EQ.LY ) THEN
C
C ----- Attention ----
C
U( J) = 1.D0
C ====================
NBU( J) = 9999
C --------------------
C
V ( J) = 0.D0
NBV( J) = 9999
C
W ( J) = 0.D0
NBW( J) = 9999
C
C --- U = 0, V = 0, W = 0 ( When IX= 0 or IY= 0 or IY= LY or IZ= 0 or IZ= LZ )
C
ELSE
1 IF ( IX.EQ.1.OR.IX.EQ.LX .OR. IY.EQ.1 .OR. IZ.EQ.1) THEN
U ( J) = 0.D0
NBU( J) = 9999
C
V ( J) = 0.D0
NBV( J) = 9999
C
W ( J) = 0.D0
NBW( J) = 9999
C
ELSE IF ( IZ.EQ.LZ) THEN
U ( J) = 0.D0
NBU( J) = 9999
C
V ( J) = 0.D0
NBV( J) = 9999
C
W ( J) = 0.D0
NBW( J) = 9999
C
END IF
C ---------------------------
C
J = J + 1
C
400 CONTINUE
C --------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE JACP ( XX, YY, ZZ, EPP, EDD )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION XX( 8), YY( 8), ZZ( 8), EPP( 8,8), EDD( 8,8),
1 P ( 64), DX( 64), DY( 64), DZ ( 64)
C
DATA P / 8.D0, 4.D0, 2.D0, 4.D0, 4.D0, 2.D0, 1.D0, 2.D0,
1 4.D0, 8.D0, 4.D0, 2.D0, 2.D0, 4.D0, 2.D0, 1.D0,
2 2.D0, 4.D0, 8.D0, 4.D0, 1.D0, 2.D0, 4.D0, 2.D0,
3 4.D0, 2.D0, 4.D0, 8.D0, 2.D0, 1.D0, 2.D0, 4.D0,
C
4 4.D0, 2.D0, 1.D0, 2.D0, 8.D0 ,4.D0, 2.D0, 4.D0,
5 2.D0, 4.D0, 2.D0, 1.D0, 4.D0, 8.D0, 4.D0, 2.D0,
6 1.D0, 2.D0, 4.D0, 2.D0, 2.D0, 4.D0, 8.D0, 4.D0,
7 2.D0, 1.D0, 2.D0, 4.D0, 4.D0, 2.D0, 4.D0, 8.D0 /
C
DATA DX / 4.D0, -4.D0, -2.D0, 2.D0, 2.D0, -2.D0, -1.D0, 1.D0,
1 -4.D0, 4.D0, 2.D0, -2.D0, -2.D0, 2.D0, 1.D0, -1.D0,
2 -2.D0, 2.D0, 4.D0, -4.D0, -1.D0, 1.D0, 2.D0, -2.D0,
3 2.D0, -2.D0, -4.D0, 4.D0, 1.D0, -1.D0, -2.D0, 2.D0,
C
4 2.D0, -2.D0, -1.D0, 1.D0, 4.D0, -4.D0, -2.D0, 2.D0,
5 -2.D0, 2.D0, 1.D0, -1.D0, -4.D0, 4.D0, 2.D0, -2.D0,
6 -1.D0, 1.D0, 2.D0, -2.D0, -2.D0, 2.D0, 4.D0, -4.D0,
7 1.D0, -1.D0, -2.D0, 2.D0, 2.D0, -2.D0, -4.D0, 4.D0 /
C
DATA DY / 4.D0, 2.D0, -2.D0, -4.D0, 2.D0, 1.D0, -1.D0, -2.D0,
1 2.D0, 4.D0, -4.D0, -2.D0, 1.D0, 2.D0, -2.D0, -1.D0,
2 -2.D0, -4.D0, 4.D0, 2.D0, -1.D0, -2.D0, 2.D0, 1.D0,
3 -4.D0, -2.D0, 2.D0, 4.D0, -2.D0, -1.D0, 1.D0, 2.D0,
C
4 2.D0, 1.D0, -1.D0, -2.D0, 4.D0, 2.D0, -2.D0, -4.D0,
5 1.D0, 2.D0, -2.D0, -1.D0, 2.D0, 4.D0, -4.D0, -2.D0,
6 -1.D0, -2.D0, 2.D0, 1.D0, -2.D0, -4.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, 1.D0, 2.D0, -4.D0, -2.D0, 2.D0, 4.D0 /
C
DATA DZ / 4.D0, 2.D0, 1.D0, 2.D0, -4.D0, -2.D0, -1.D0, -2.D0,
1 2.D0, 4.D0, 2.D0, 1.D0, -2.D0, -4.D0, -2.D0, -1.D0,
2 1.D0, 2.D0, 4.D0, 2.D0, -1.D0, -2.D0, -4.D0, -2.D0,
3 2.D0, 1.D0, 2.D0, 4.D0, -2.D0, -1.D0, -2.D0, -4.D0,
C
4 -4.D0, -2.D0, -1.D0, -2.D0, 4.D0, 2.D0, 1.D0, 2.D0,
5 -2.D0, -4.D0, -2.D0, -1.D0, 2.D0, 4.D0, 2.D0, 1.D0,
6 -1.D0, -2.D0, -4.D0, -2.D0, 1.D0, 2.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, -2.D0, -4.D0, 2.D0, 1.D0, 2.D0, 4.D0 /
C
C --------------------------------------------------------------------
C
A = ( DABS( XX( 2) - XX( 1)))/ 2.D0
B = ( DABS( YY( 4) - YY( 1)))/ 2.D0
C = ( DABS( ZZ( 5) - ZZ( 1)))/ 2.D0
C
RP = A* B* C / 27.D0 ! RP = 8.D0* A* B* C / 216.D0
C
RDX = ( B* C )/( 18.D0* A )
RDY = ( C* A )/( 18.D0* B )
RDZ = ( A* B )/( 18.D0* C )
C
K = 1
C ----------------------------------------------------------
DO 100 I = 1, 8
DO 100 J = 1, 8
C
EPP( I,J) = RP* P( K)
EDD( I,J) = RDX* DX( K) + RDY* DY( K) + RDZ* DZ( K)
K = K + 1
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE JCPD ( XX, YY, ZZ, EXX, EYY, EZZ, EFX, EFY,
1 EFZ, EPP )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION XX ( 8), YY ( 8), ZZ ( 8), EXX( 8,8),
1 EYY( 8,8), EZZ( 8,8), EFX( 8,8), EFY( 8,8),
2 EFZ( 8,8), EPP( 8,8)
C
DIMENSION DX( 64), DY( 64), DZ( 64), FX( 64),
1 FY( 64), FZ( 64), P ( 64)
C
DATA P / 8.D0, 4.D0, 2.D0, 4.D0, 4.D0, 2.D0, 1.D0, 2.D0,
1 4.D0, 8.D0, 4.D0, 2.D0, 2.D0, 4.D0, 2.D0, 1.D0,
2 2.D0, 4.D0, 8.D0, 4.D0, 1.D0, 2.D0, 4.D0, 2.D0,
3 4.D0, 2.D0, 4.D0, 8.D0, 2.D0, 1.D0, 2.D0, 4.D0,
C
4 4.D0, 2.D0, 1.D0, 2.D0, 8.D0, 4.D0, 2.D0, 4.D0,
5 2.D0, 4.D0, 2.D0, 1.D0, 4.D0, 8.D0, 4.D0, 2.D0,
6 1.D0, 2.D0, 4.D0, 2.D0, 2.D0, 4.D0, 8.D0, 4.D0,
7 2.D0, 1.D0, 2.D0, 4.D0, 4.D0, 2.D0, 4.D0, 8.D0 /
C
DATA DX / 4.D0, -4.D0, -2.D0, 2.D0, 2.D0, -2.D0, -1.D0, 1.D0,
1 -4.D0, 4.D0, 2.D0, -2.D0, -2.D0, 2.D0, 1.D0, -1.D0,
2 -2.D0, 2.D0, 4.D0, -4.D0, -1.D0, 1.D0, 2.D0, -2.D0,
3 2.D0, -2.D0, -4.D0, 4.D0, 1.D0, -1.D0, -2.D0, 2.D0,

4 2.D0, -2.D0, -1.D0, 1.D0, 4.D0, -4.D0, -2.D0, 2.D0,
5 -2.D0, 2.D0, 1.D0, -1.D0, -4.D0, 4.D0, 2.D0, -2.D0,
6 -1.D0, 1.D0, 2.D0, -2.D0, -2.D0, 2.D0, 4.D0, -4.D0,
7 1.D0, -1.D0, -2.D0, 2.D0, 2.D0, -2.D0, -4.D0, 4.D0 /
C
DATA DY / 4.D0, 2.D0, -2.D0, -4.D0, 2.D0, 1.D0, -1.D0, -2.D0,
1 2.D0, 4.D0, -4.D0, -2.D0, 1.D0, 2.D0, -2.D0, -1.D0,
2 -2.D0, -4.D0, 4.D0, 2.D0, -1.D0, -2.D0, 2.D0, 1.D0,
3 -4.D0, -2.D0, 2.D0, 4.D0, -2.D0, -1.D0, 1.D0, 2.D0,

4 2.D0, 1.D0, -1.D0, -2.D0, 4.D0, 2.D0, -2.D0, -4.D0,
5 1.D0, 2.D0, -2.D0, -1.D0, 2.D0, 4.D0, -4.D0, -2.D0,
6 -1.D0, -2.D0, 2.D0, 1.D0, -2.D0, -4.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, 1.D0, 2.D0, -4.D0, -2.D0, 2.D0, 4.D0 /
C
DATA DZ / 4.D0, 2.D0, 1.D0, 2.D0, -4.D0, -2.D0, -1.D0, -2.D0,
1 2.D0, 4.D0, 2.D0, 1.D0, -2.D0, -4.D0, -2.D0, -1.D0,
2 1.D0, 2.D0, 4.D0, 2.D0, -1.D0, -2.D0, -4.D0, -2.D0,
3 2.D0, 1.D0, 2.D0, 4.D0, -2.D0, -1.D0, -2.D0, -4.D0,

4 -4.D0, -2.D0, -1.D0, -2.D0, 4.D0, 2.D0, 1.D0, 2.D0,
5 -2.D0, -4.D0, -2.D0, -1.D0, 2.D0, 4.D0, 2.D0, 1.D0,
6 -1.D0, -2.D0, -4.D0, -2.D0, 1.D0, 2.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, -2.D0, -4.D0, 2.D0, 1.D0, 2.D0, 4.D0 /
C
DATA FX / -4.D0, 4.D0, 2.D0, -2.D0, -2.D0, 2.D0, 1.D0, -1.D0,
1 -4.D0, 4.D0, 2.D0, -2.D0, -2.D0, 2.D0, 1.D0, -1.D0,
2 -2.D0, 2.D0, 4.D0, -4.D0, -1.D0, 1.D0, 2.D0, -2.D0,
3 -2.D0, 2.D0, 4.D0, -4.D0, -1.D0, 1.D0, 2.D0, -2.D0,

4 -2.D0, 2.D0, 1.D0, -1.D0, -4.D0, 4.D0, 2.D0, -2.D0,
5 -2.D0, 2.D0, 1.D0, -1.D0, -4.D0, 4.D0, 2.D0, -2.D0,
6 -1.D0, 1.D0, 2.D0, -2.D0, -2.D0, 2.D0, 4.D0, -4.D0,
7 -1.D0, 1.D0, 2.D0, -2.D0, -2.D0, 2.D0, 4.D0, -4.D0 /
C
DATA FY / -4.D0, -2.D0, 2.D0, 4.D0, -2.D0, -1.D0, 1.D0, 2.D0,
1 -2.D0, -4.D0, 4.D0, 2.D0, -1.D0, -2.D0, 2.D0, 1.D0,
2 -2.D0, -4.D0, 4.D0, 2.D0, -1.D0, -2.D0, 2.D0, 1.D0,
3 -4.D0, -2.D0, 2.D0, 4.D0, -2.D0, -1.D0, 1.D0, 2.D0,

4 -2.D0, -1.D0, 1.D0, 2.D0, -4.D0, -2.D0, 2.D0, 4.D0,
5 -1.D0, -2.D0, 2.D0, 1.D0, -2.D0, -4.D0, 4.D0, 2.D0,
6 -1.D0, -2.D0, 2.D0, 1.D0, -2.D0, -4.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, 1.D0, 2.D0, -4.D0, -2.D0, 2.D0, 4.D0 /
C
DATA FZ / -4.D0, -2.D0, -1.D0, -2.D0, 4.D0, 2.D0, 1.D0, 2.D0,
1 -2.D0, -4.D0, -2.D0, -1.D0, 2.D0, 4.D0, 2.D0, 1.D0,
2 -1.D0, -2.D0, -4.D0, -2.D0, 1.D0, 2.D0, 4.D0, 2.D0,
3 -2.D0, -1.D0, -2.D0, -4.D0, 2.D0, 1.D0, 2.D0, 4.D0,
C
4 -4.D0, -2.D0, -1.D0, -2.D0, 4.D0, 2.D0, 1.D0, 2.D0,
5 -2.D0, -4.D0, -2.D0, -1.D0, 2.D0, 4.D0, 2.D0, 1.D0,
6 -1.D0, -2.D0, -4.D0, -2.D0, 1.D0, 2.D0, 4.D0, 2.D0,
7 -2.D0, -1.D0, -2.D0, -4.D0, 2.D0, 1.D0, 2.D0, 4.D0 /
C
A = ( DABS( XX( 7) - XX( 1)))/ 2.D0
B = ( DABS( YY( 7) - YY( 1)))/ 2.D0
C = ( DABS( ZZ( 7) - ZZ( 1)))/ 2.D0
C
RP = 8.D0* A* B* C / 216.D0
C
RDX = ( B* C)/( 18.D0* A)
RDY = ( C* A)/( 18.D0* B)
RDZ = ( A* B)/( 18.D0* C)
C
RFX = ( B* C)/ 18.D0
RFY = ( C* A)/ 18.D0
RFZ = ( A* B)/ 18.D0
C
K = 1
C --------------------------------------
DO 100 I = 1, 8
DO 100 J = 1, 8
C --------------------------------------
C
EXX( I,J) = RDX* DX( K)
EYY( I,J) = RDY* DY( K)
EZZ( I,J) = RDZ* DZ( K)
C
EFX( I,J) = RFX* FX( K)
EFY( I,J) = RFY* FY( K)
EFZ( I,J) = RFZ* FZ( K)
EPP( I,J) = RP * P ( K)
C
K = K + 1
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CKAK ( U, V, W, GU, GV, GW, XL, YL, ZL,
1 HX, HY, HZ, FF, LX, LY, LZ, NNP )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NNP), V ( NNP), W ( NNP), GU( NNP),
1 GV( NNP), GW( NNP), XL( LX ), YL( LX ),
2 ZL( LZ), HX( LX+1), HY( LY+1), HZ( LZ+1),
3 FF( NNP,3)
C
COMMON /CNST/ DT, COF, PHI
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
C
NX2 = LX + 1
NY2 = LY + 1
NZ2 = LZ + 1
C
C --------------------------------------
DO 100 I = 2, NX2-1
C
HX( I) = XL( I) - XL( I-1)
C
100 CONTINUE
C --------------------------------------
HX( 1) = HX( 2)
HX( NX2) = HX( NX2-1)
C
C --------------------------------------
DO 200 I = 2, NY2-1
C
HY( I) = YL( I) - YL( I-1)
C
200 CONTINUE
C -------------------------------------
HY( 1) = HY( 2)
HY( NY2) = HY( NY2-1)
C
C --------------------------------------
DO 300 I = 2, NZ2-1
C
HZ( I) = ZL( I) - ZL( I-1)
C
300 CONTINUE
C --------------------------------------
HZ( 1) = HZ( 2)
HZ( NZ2) = HZ( NZ2-1)
C
C ------------------------------------------------
L = 1
C ------------------------------------------------
DO 400 K = 1, LZ
DO 400 J = 1, LY
DO 400 I = 1, LX
C ------------------------------------------------
C
AX = HX( I+1)/ HX( I)
AY = HY( J+1)/ HY( J)
AZ = HZ( K+1)/ HZ( K)
C
BX = U( L)* DT/ HX( I)
BY = V( L)* DT/ HY( J)
BZ = W( L)* DT/ HZ( K)
C
RX = COF* DT/( HX( I)* HX( I))
RY = COF* DT/( HY( J)* HY( J))
RZ = COF* DT/( HZ( K)* HZ( K))
C
BXU = DT/ HX( I)
BYV = DT/ HY( J)
BZW = DT/ HZ( K)
C
IF ( ISL.EQ.0) GO TO 420
C
C ----- Corrected coefficient of Diffusivity --------------------
C
FF( L,1) = 1.D0
C
1 + 0.5D0* (( BX+BY+BZ )* ( BX+BY+BZ ))/( RX+RY+RZ )
C
2 + ( 1.D0/ 6.D0 )* (( AX-1.D0)* BX + ( AY-1.D0)* BY
3 + ( AZ-1.D0)* BZ )/( RX + RY + RZ )
C
C ----------------------------------------------------------------
C
FF( L,2) = FF( L,1)
FF( L,3) = FF( L,1)
C
GU( L) = U( L)
GV( L) = V( L)
GW( L) = W( L)
GO TO 440
C
420 CONTINUE
C
C ----- Corrected coefficient of Diffusivity ---------------------------
C
FF( L,1) = 1.D0
C
1 + 0.5D0* (( BX+BY+BZ )* ( BX+BY+BZ ))/( RX+RY+RZ )
2 - 0.5D0/( RX+RY+RZ )
C
C ----------------------------------------------------------------------
C
FF( L,2) = FF( L,1)
FF( L,3) = FF( L,1)
C
GU( L) = U( L)
GV( L) = V( L)
GW( L) = W( L)
C
440 CONTINUE
C
L = L + 1
C
400 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE OUTP ( U, V, W, P, X, Y, Z, XL, YL,
1 ZL, NDE, LX, LY, LZ, NNP, NEL, * )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U( NNP), V( NNP), W ( NNP), P ( NNP), X ( NNP),
1 Y( NNP), Z( NNP), XL( LX ), YL( LY ), ZL( LZ),
2 NDE( NEL,8)
C
COMMON /CNST/ DT, COF, PHI
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
COMMON /SINF/ ITU, ITV, ITW, ITP, IU2, IV2, IW2, IP2
C
SASV = 0.D0
C ------------------------------------------------
DO 100 I = 1, NEL
C
DVX = DVLX( X, U, NDE, NNP, NEL, I )
DVY = DVLY( Y, V, NDE, NNP, NEL, I )
DVZ = DVLZ( Z, W, NDE, NNP, NEL, I )
C
ASDV = DVX + DVY + DVZ
SASV = SASV + ASDV
C
100 CONTINUE
C ------------------------------------------------
C
AVRV = DABS( SASV/ DFLOAT( NEL))
C
C ------------------------------------------------
IF ( MOD( LOP,NPR).EQ.0 ) THEN
WRITE(6,2000) LOP, AVRV, ITU, ITV, ITW, ITP, IU2, IV2, IW2,
1 IP2, VVR( LOP), PVR( LOP)
C
WRITE(10,2000) LOP, AVRV, ITU, ITV, ITW, ITP, IU2, IV2, IW2,
1 IP2, VVR( LOP), PVR( LOP)
C
2000 FORMAT(' L=',I4,' Av_R=',D9.2,' It=', 8I3,' V_Vr=',F6.2,' %',
1 ' V_Pm=',F6.2,' %')
END IF
C ------------------------------------------------
C
IF ( EPX.LT.EPL) GO TO 150 ! < 0.1(%) ==> Var. of Max_Velocity
C ---------------------------------------------------------------------
C
** IF ( EPX.LT.EPL .AND. ESP.LT.EPL ) GO TO 150 ! < 0.1(%)
IF ( UMM.GT.1.1D0) GO TO 250
IF ( LOP.LT.LPE) RETURN 1
C
WRITE(6,*) ' *** Loop over ***'
WRITE(10,*) ' *** Loop over ***'
RETURN
C
250 WRITE(6,*) ' *** Umax. > 1.1 * STOP * Umax.=', UMM
WRITE(10,*) ' *** Umax. > 1.1 * STOP * Umax.=', UMM
STOP
C
150 CONTINUE
WRITE(6,2100) LOP
2100 FORMAT(/' *** Stationary Condition Reached : Loop =',I4)
WRITE(6,2200) AVRV, ITU, ITV, ITW, ITP, IU2, IV2,
1 IW2, IP2, VVR( LOP), PVR( LOP)
2200 FORMAT(' * Av_R=',D9.2,' It=', 8I3,' V_Vr=',F6.2,' %',
1 ' V_Pm=',F6.2,' %')
WRITE(10,2100) LOP
WRITE(10,2200) AVRV, ITU, ITV, ITW, ITP, IU2, IV2,
1 IW2, IP2, VVR( LOP), PVR( LOP)
C
RETURN
END
C *********************************************************************
C
SUBROUTINE PLTO ( SU, SV, SW, P, X, Y, Z, XL,
1 YL, ZL, LXY, LYZ, LZX, LX, LY, LZ,
2 NNP, NXY, NYZ, NZX, MES )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
C
DIMENSION SU( NNP), SV( NNP), SW( NNP), P ( NNP), X ( NNP),
1 Y ( NNP), Z ( NNP), XL( LX), YL( LY), ZL( LY)
C
DIMENSION LXY( NXY,4), LYZ( NYZ,4), LZX( NZX,4)
C
COMMON /ZDFC/ ZFC
COMMON /CNST/ DT, COF, PHI
COMMON /CON2/ RE, LPE, ISL, EPL, EPX, UMM, ESP, SPX
C
NX = LX - 1
NY = LY - 1
NZ = LZ - 1
C
NNP = LX* LY* LZ
J = 1
C --------------------------------------
DO 100 IZ = 1, LZ
DO 100 IY = 1, LY
DO 100 IX = 1, LX
C --------------------------------------
C
X( J) = XL( IX)
Y( J) = YL( IY)
Z( J) = ZL( IZ)
IF ( IZ.EQ.LZ) SW( J) = 0.D0
J = J + 1
C
100 CONTINUE
C ---------------------------------------
C
IZI = LZ
IXI = LZ
IYI = LZ
NPXY = LX* LY
NPYZ = LY* LZ
NPZX = LZ* LX
C
J = 1
C --------------------------------------
DO 200 IY = 1, LY-1
DO 200 IX = 1, LX-1
C --------------------------------------
C
IST = IX + ( IY-1)* LX + ( IZI-1)/ 2* LX* LY
LXY( J,1) = IST
LXY( J,2) = IST + 1
LXY( J,3) = IST + LX + 1
LXY( J,4) = IST + LX
J = J + 1
C
200 CONTINUE
C --------------------------------------
C
J = 1
C --------------------------------------
DO 300 IZ = 1, LZ-1
DO 300 IY = 1, LY-1
C --------------------------------------
C
IST = ( IXI-1)/2 + 1 + ( IY-1)* LX + ( IZ-1)* LX* LY
LYZ( J,1) = IST
LYZ( J,2) = IST + LX
LYZ( J,3) = IST + LX* ( LY+1)
LYZ( J,4) = IST + LX* LY
J = J + 1
C
300 CONTINUE
C -------------------------------------
C
J = 1
C -------------------------------------
DO 400 IX = 1, LX-1
DO 400 IZ = 1, LZ-1
C -------------------------------------
C
IST = IX + ( IYI-1)/2* LX + ( IZ-1)* LX* LY
LZX( J,1) = IST
LZX( J,2) = IST + LX* LY
LZX( J,3) = IST + LX* LY + 1
LZX( J,4) = IST + 1
J = J + 1
C
400 CONTINUE
C -------------------------------------
C
WRITE(11,2000) LX, LY, LZ, LX-1, MES
2000 FORMAT(5I3)
WRITE(11,2100) ( XL( I), I=1, LX ), ( YL( I), I=1, LY )
WRITE(11,2200) RE, DT, LOP, ZFC
2200 FORMAT(F7.1,F7.4,I5,D11.4)
WRITE(11,2100) ( VVR( I), I = 1, LOP)
WRITE(11,2100) ( PVR( I), I = 1, LOP)
WRITE(11,2300) NNP
2300 FORMAT(I5)
WRITE(11,2100) ( X ( I), Y ( I), Z ( I), I=1, NNP )
WRITE(11,2100) ( SU( I), SV( I), SW( I), I=1, NNP )
2100 FORMAT(9D14.7)
WRITE(11,2010) NXY, NYZ, NZX
2010 FORMAT(3I5)
WRITE(11,2100) ( P( I), I = 1, NNP)
WRITE(11,2400) (( LXY( I,J), I=1, NXY), J = 1, 4 )
WRITE(11,2400) (( LYZ( I,J), I=1, NYZ), J = 1, 4 )
WRITE(11,2400) (( LZX( I,J), I=1, NZX), J = 1, 4 )
2400 FORMAT(16I5)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STNC ( NCA, NNZ, NDE, NEL, NP, P )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION NDE( NEL,NP), P( NNZ), NCA( 2,NNZ)
C
K = 0
C -------------------------------
DO 100 INEL = 1, NEL
DO 100 I = 1, NP
C -------------------------------
C
II = NDE( INEL,I)
C
C -------------------------------
DO 100 J = 1, NP
C -------------------------------
C
JJ = NDE( INEL,J)
C
CALL DCNL ( 0.D0, NCA, II,JJ, P, K, NNZ )
C
100 CONTINUE
C -------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STNN ( NDE, NC1, NC2, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION NDE( NEL,8), NC1( NNP), NC2( NNP,8)
C
C --------------------------------------
DO 100 I = 1, NNP
C
NC1( I) = 0
100 CONTINUE
C --------------------------------------
DO 200 J = 1, 8
DO 200 I = 1, NNP
C
NC2( I,J) = 0
200 CONTINUE
C --------------------------------------
DO 300 J = 1, 8
DO 300 I = 1, NEL
C
N = NDE( I,J)
NC1( N) = NC1( N)+ 1
NC2( N,NC1( N)) = I
C
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DSLVE ( NNP, NNZ, NZ2, FVC, C, TSM, NCA, KCC, NBC,
1 NTL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION FVC( NNP), C( NNP), KCC( NNP), NBC( NNP), TSM( NNZ),
1 NCA( 2,NZ2)
C
COMMON /CON1/ OMG, ERR, NTR
COMMON /VRDT/ LOP, VVR( 2000), PVR( 2000), NPR
C
C ----- S.O.R. method ---
C
CMAX = 0.D0
C
C --------------------------------------
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMAX) CMAX = AC
C
100 CONTINUE
C --------------------------------------
C
EPS = CMAX* 0.01D0
LMP = 3* NNP
KTR = 2* NNP
C ----------------------------
C
NTL = 0
C
150 CONTINUE
NTL = NTL + 1
ERX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999 ) GO TO 200
C
IF ( I.EQ.1) GO TO 220
MS = KCC( I-1)+ 1
GO TO 240
C
220 MS = 1
240 ML = KCC( I)
IF ( MS.GT.ML) GO TO 290
DK = FVC( I)
C
C --------------------------------------
DO 260 M = MS, ML
C
M1 = NCA( 1,M)
M2 = NCA( 2,M)
IF ( M2.EQ.I) DSM = 1.D0 / TSM( M1)
IF ( M2.EQ.I) GO TO 260
C
DK = DK - TSM( M1)* C( M2 )
C
260 CONTINUE
C --------------------------------------
C
DC = - C( I) + DSM* DK
C( I) = C( I) + OMG* DC
TEMP = DABS( C( I))
C
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS( DC/ C( I))
C
IF ( TEMP.LE.ERX) GO TO 200
ERX = TEMP
NRX = I
C
200 CONTINUE
C ------------------------------------------------
C
IF ( ERX.LE.ERR) RETURN
IF ( NTL.GE.KTR) GO TO 320
IF ( NTL.GE.LMP) GO TO 340
GO TO 150
C
340 CONTINUE
WRITE(6,2000) EPS, NTL, NRX, ERX
2000 FORMAT(' * EPS=',E11.4,' CITAT=',I6,2X,' NDE=',I5,2X,' ERR=',
1 E11.4)
WRITE(6,2100) ( C( I),I = 1, NNP)
2100 FORMAT(' ',6E18.6)
GO TO 150
C
320 CONTINUE
WRITE(6,2200) NTL, LOP
2200 FORMAT(' *** No. of iterations exceeded ( DSLVE ), NTL =',I5,
1 ' at LOP=',I5)
WRITE(10,2200) NTL, LOP
RETURN
C
290 NER = 999
WRITE(6,2300) NER, I, MS, ML
2300 FORMAT(' ** Error',4I6)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SOLVP ( NNP, NNZ, NZ2, FVC, C, TSM, NCA, KCC, NBC,
1 NTL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION FVC( NNP), C( NNP), KCC( NNP), NBC( NNP), TSM( NNZ),
1 NCA ( 2,NZ2)
C
COMMON /CON1/ OMG, ERR, NTR
C
C ----- S.O.R. method ---
C
CMAX = 0.D0
C
C -----------------------------------------
DO 100 I = 1, NNP
C
AC = DABS( C( I))
IF ( AC.GE.CMAX) CMAX = AC
C
100 CONTINUE
C -----------------------------------------
C
EPS = CMAX* 0.01D0
LMP = 3* NNP
C
C ----- Attention ---
C
KTR = 2* NNP
NTL = 0
C
150 NTL = NTL + 1
ERX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999 ) GO TO 200
IF ( I.EQ.1) GO TO 220
MS = KCC( I-1)+ 1
GO TO 240
C
220 MS = 1
240 ML = KCC( I)
IF ( MS.GT.ML) GO TO 290
DK = FVC( I)
C
C --------------------------------------
DO 260 M = MS, ML
C
M1 = NCA( 1,M )
M2 = NCA( 2,M )
IF ( M2.EQ.I) DSM = 1.D0 / TSM( M1)
IF ( M2.EQ.I) GO TO 260
DK = DK - TSM( M1)* C( M2)
C
260 CONTINUE
C --------------------------------------
C
DC = - C( I) + DSM* DK
C( I) = C( I) + OMG* DC
TEMP = DABS( C( I))
IF ( TEMP.LE.EPS) GO TO 200
C
TEMP = DABS( DC/ C( I))
IF ( TEMP.LE.ERX) GO TO 200
ERX = TEMP
NRX = I
C
200 CONTINUE
C ------------------------------------------------
IF ( ERX.LE.ERR ) RETURN
IF ( NTL.GE.KTR ) GO TO 320
IF ( NTL.GE.LMP ) GO TO 340
GO TO 150
C
340 CONTINUE
WRITE(6,2000) EPS, NTL, NRX, ERX
2000 FORMAT(' EPS=',E11.4,' CITAT=',I6,2X,' NDE=',I5,2X,' ERR=',
1 E11.4)
WRITE(6,2100) ( C( I),I = 1, NNP)
2100 FORMAT(' ',6E18.6)
GO TO 150
C
320 CONTINUE
WRITE(6,2300) NTL
2300 FORMAT(/' *** No. of iterations exceeded ( SOLVP), NTL =',I5)
WRITE(10,2300) NTL
RETURN
C
290 NER = 999
WRITE(6,2400) NER, I, MS, ML
2400 FORMAT(' ***** Error ***', 4I6)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( PDV, I,J, N, PD, NCA, K, NNZ )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION PD( NNZ), NCA( 2,NNZ)
C
K = K + 1
C
IF ( K.GT.NNZ) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NNZ *')
C
JE = K - N
C --------------------------------------
DO 100 JJ = 1, JE
C
J1 = K - JJ + 1
PD ( J1) = PD ( J1-1)
NCA( 1,J1) = NCA( 1,J1-1)
NCA( 2,J1) = NCA( 2,J1-1)
C
100 CONTINUE
C --------------------------------------
C
PD ( N) = PDV
NCA( 1,N) = I
NCA( 2,N) = J
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLPR ( P, DLP, DE, DI, DLX, DLY, DLZ,
1 NDE, NC2, NC1, SU, SV, SW, SP,
2 X, Y, Z, NBU, NBV, NBW, NNP,
3 NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION P ( NNP), DLP( NNP), DE ( NEL), DI ( NNP),
1 DLX( NEL), DLY( NEL), DLZ( NEL), NDE( NEL,8),
2 NC2( NNP,8), NC1( NNP), SU ( NNP), SV ( NNP),
3 SW ( NNP), SP( NNP), X ( NNP), Y ( NNP),
4 Z ( NNP), NBU( NNP), NBV( NNP), NBW( NNP)
C
COMMON /CNST/ DT, COF, PHI
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
DLX( I) = DVLX( X, DLP, NDE, NNP, NEL, I )
DLY( I) = DVLY( Y, DLP, NDE, NNP, NEL, I )
DLZ( I) = DVLZ( Z, DLP, NDE, NNP, NEL, I )
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NEL
C
DE( I) = DLX( I)
C
200 CONTINUE
C ----------------------------------------------------------
C
CALL VLNDE ( DI, DE, NC2, NC1, NNP, NEL )
C
C ----------------------------------------------------------
DO 300 I = 1, NNP
C
IF ( NBU( I).NE.9999) SU( I) = SU( I) - DT* DI( I)
C
300 CONTINUE
C ----------------------------------------------------------
DO 400 I = 1, NEL
C
DE( I) = DLY( I)
C
400 CONTINUE
C ------------------------------------------------
C
CALL VLNDE ( DI, DE, NC2, NC1, NNP, NEL )
C
C ---------------------------------------------------------
DO 500 I = 1, NNP
C
IF ( NBV( I).NE.9999 ) SV( I) = SV( I) - DT* DI( I)
C
500 CONTINUE
C ---------------------------------------------------------
DO 600 I = 1, NEL
C
DE( I) = DLZ( I)
C
600 CONTINUE
C ------------------------------------------------
C
CALL VLNDE ( DI, DE, NC2, NC1, NNP, NEL )
C
C ------------------------------------------------
DO 700 I = 1, NNP
C
IF ( NBW( I).NE.9999 ) SW( I) = SW( I) - DT* DI( I)
C
700 CONTINUE
C
C ----- Pressure -----------------
C
DO 800 I = 1, NNP
C
P( I) = P( I) + DLP( I)
C
800 CONTINUE
C --------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLNDE ( CI, CE, NC2, NC1, NNP, NEL )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION CI( NNP), CE( NEL), NC2( NNP,8), NC1( NNP)
C
C --------------------------------------
DO 100 I = 1, NNP
C
NQ = NC1( I)
FNN = 1.D0/ DFLOAT( NQ)
CT = 0.D0
C
C ----------------------------
DO 150 J = 1, NQ
C
JJ = NC2( I,J)
CT = CT + CE( JJ)
C
150 CONTINUE
C ----------------------------
C
CI( I) = CT* FNN
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************