戻る

C *******************************************************************C
C C
C H F 2 D - F E T C
C C
C Thermal Fluid Analysis Software by Finite Element Method C
C C
C ( Triangular-Linear Element ) C
C C
C ( V.4.3 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C *******************************************************************C
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, ITB, IEB, KCC, MTR, IBP, IBV, IBT, LBC,
1 NED, NDB, NCN, NTB, NC2, IGL, IRS, NGL, NRS,
2 ICN
C
C ----- Attention Data for ( 17 x 17 x 4 ) ---
C
C ( 113 nodes and 196 elements )
C
DIMENSION X ( 113), Y ( 113), NDE( 196,3), ITB( 113),
1 IEB( 196), AT ( 196), U ( 196), V ( 196),
2 UN ( 113), VN ( 113), MTR( 196), LBC( 28),
3 P1 ( 421), P2 ( 421), NC2( 2,729), KCC( 113),
4 NCW( 421), FVC( 113), BE ( 196,3), CE ( 196,3)
C
DIMENSION DP ( 421), D1 ( 421), D2 ( 421), NED( 28),
1 DSH( 28), SH ( 28), DX1( 196), DY1( 196),
2 DX2( 196), DY2( 196), PSI( 113), PSO( 113),
3 IBP( 113), VOR( 113), VRO( 113), IBV( 113),
4 TEP( 113), TPO( 113), IBT( 113), QWK( 196),
5 QE ( 196)
C
DIMENSION NTB( 113,10), NDB( 28,10), NCN( 113), NCC1( 729),
1 IGL( 56), IRS( 616)
C
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN

COMMON /FLW/ NVRC, IVC( 50), VRC( 50), NPC, IPC( 50),
1 PSC( 50), NTPC, ITC( 50), TPC( 50), COF(100)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /ITL/ IPS, ITP, IVR
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
C
C ----- Attention --- < = 2000 ---
C
COMMON /VVT/ VARV( 2000), VART( 2000), VARP( 2000)
C
C ----- Attention --- 17x17x4 ---
C
OPEN ( 1, FILE='SORTIN.DAT')
OPEN ( 2, FILE='SORTOUT.DAT')
OPEN ( 7, FILE='HFE_T_I7.DAT')
OPEN ( 8, FILE='HFE_T_OUT.DAT')
OPEN (10, FILE='HFE_T_FIG.DAT')
C
CALL TTLE
C
CALL INPT ( X, Y, NDE, ITB, IEB, AT, U, V, DSH, SH,
1 NDB, NED, MTR, LBC, NNZ, NZ2, ICN )
C
CALL BILD ( NDE, IGL, IRS, NEL, NGL, NRS )
C
LOP = 0
C
CALL INIT ( IBP, IBV, LBC, X, Y, NDE, BE, CE, PSI,
1 PSO, VOR, VRO, AT, U, V, NCN, NTB, DX1,
2 DY1, DX2, DY2, MTR, IBT, TEP, TPO, QE, ICN )
C
C ----- P1 - Matrix ---
C
CALL PMAT ( AT, NDE, P1, NCC1, NCW, KCC, NNZ, NZ2, NC2 )
C
C ----- DP & D1 Matrix ---
C
CALL DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, DP, NCW, NEL,
1 NNZ, 1 )
C
CALL DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, D1, NCW, NEL,
1 NNZ, 2 )
C
CALL DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, D2, NCW, NEL,
1 NNZ, 3 )
C
C ----- Vorticity on Walls ---
C
CALL VORW ( SH, DSH, PSI, VRO, NED, LBC, NDB, NNP, LUV, ICN )
C
5555 CONTINUE
C
LOP = LOP + 1
C
C ----- Stream function equation ---
C
CALL CPSI ( PSI, PSO, IBP, VRO, DP, AT, BE, CE, DX1,
1 DY1, DX2, DY2, FVC, NC2, KCC, NDE, U, V,
2 MTR, COF, NNZ, NZ2, ICN )
C
C ----- Vorticity on Walls ---
C
CALL VORW ( SH, DSH, PSI, VOR, NED, LBC, NDB, NNP, LUV, ICN )
C
C ----- To modify D1 & D2 ---
C
CALL MODD ( D1, P1, P2, AT, BE, CE, DX1, DY1, DX2, DY2,
1 NCW, NDE, NEL, NNZ, 2 )
C
IF ( ICN.EQ.0) GO TO 100
C
CALL CTEP ( TEP, TPO, IBT, D1, P2, AT, BE, CE, U, V,
1 QWK, FVC, QE, UN, VN, NC2, KCC, NDE, NTB, NCN,
2 LBC, LUV, NNP, NNZ, NZ2, NEL, ICN )
C
100 CONTINUE
C
C ----- Vorticity Equation ---
C
CALL MODD ( D2, P1, P2, AT, BE, CE, DX1, DY1, DX2,
1 DY2, NCW, NDE, NEL, NNZ, 3 )
C
CALL LPVR ( VOR, VRO, IBV, D2, P2, AT, BE, CE, U,
1 V, UN, VN, FVC, NC2, KCC, NDE, NTB, NCN,
2 LBC, NNZ, NZ2, TEP, QWK, ICN )
C
CALL PRTS ( PSI, VOR, TEP, U, V, ICN )
C
C ----- Stationrity check ---
C
CALL CHCK ( PSI, PSO, VOR, VRO, TEP, TPO, U, V, NNP, NEL,
1 ICK, ICN )
C
IF ( ICK.EQ.1) GO TO 999
GO TO 5555
C
999 CONTINUE
C
CALL PLTD ( PSI, VOR, TEP, UN, VN, X, Y, NDE, IRS, NRS, ICN )
C
CLOSE( 1)
CLOSE( 2)
CLOSE( 7)
CLOSE( 8)
CLOSE(10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE BCMAT ( X, Y, NDE, BE, CE, NNP, NEL)
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,3), BE( NEL,3), CE( NEL,3)
C
C --------------------------------------
DO 100 I = 1 ,NEL
C
N1 = NDE( I,1)
N2 = NDE( I,2)
N3 = NDE( I,3)
C
BE( I,1)= Y( N2) - Y( N3)
BE( I,2)= Y( N3) - Y( N1)
BE( I,3)= Y( N1) - Y( N2)
C
CE( I,1)= X( N3) - X( N2)
CE( I,2)= X( N1) - X( N3)
CE( I,3)= X( N2) - X( N1)
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BCSB ( DSH, SH, NDB, NED, LUV )
C
C ----- Boundary Condition ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDB, NED
C
DIMENSION NDB( LUV,10), NED( LUV), DSH( LUV), SH( LUV)
C
C ----- Attention : DATA for 113 / 196 -----
C
DO 100 I = 1, LUV
C
SH ( I)= 21.D0
DSH( I)= 147.D0
NED( I)= 4 ! 4 corners
C
100 CONTINUE
C ------------------------------------------
C
C ----- Attention : Vort.=0 at 4 corners ---
C
NED ( 1)= 0
NED ( 8)= 0
NED ( 15)= 0
NED ( 22)= 0
C
NDB ( 1,1)= 9
NDB ( 2,1)= 17
NDB ( 3,1)= 32
NDB ( 4,1)= 47
NDB ( 5,1)= 62
NDB ( 6,1)= 77
NDB ( 7,1)= 92
NDB ( 8,1)= 99
NDB ( 9,1)= 92
NDB (10,1)= 93
C
NDB (11,1)= 94
NDB (12,1)= 95
NDB (13,1)= 96
NDB (14,1)= 97
NDB (15,1)=105
NDB (16,1)= 97
NDB (17,1)= 82
NDB (18,1)= 67
NDB (19,1)= 52
NDB (20,1)= 37
C
NDB (21,1)= 22
NDB (22,1)= 15
NDB (23,1)= 22
NDB (24,1)= 21
NDB (25,1)= 20
NDB (26,1)= 19
NDB (27,1)= 18
NDB (28,1)= 17
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BILD ( NDE, IGL, IRS, NEL, NGL, NRS )
C
C ----- Outline & Mesh-Line ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, IGL, IRS, NGL
DIMENSION NDE( NEL,3), IGL( 56), IRS( 616), KZ5( 2,3)
C
DATA KZ5 / 1,2, 2,3, 3,1 /
C
C ----- SORT ----------------------
DO 100 I = 1, NEL
DO 100 J = 1, 3
C ---------------------------------
C
N1 = NDE( I, KZ5(1,J))
N2 = NDE( I, KZ5(2,J))
C
IF ( N1.LE.N2) GO TO 150
N3 = N1
N1 = N2
N2 = N3
150 CONTINUE
WRITE(1,2000) N1, N2
2000 FORMAT(2I5)
C
100 CONTINUE
C----------------------------------
C
CALL SORTM ( NEL)
C
C ---------------------------------
C
NGL = 0
NRS = 0
C
IAL = 0
ILY = 0
I = 0
C
C ----- IGL : Comp. element for outline
C IRS : Comp. element for Meshes ---
C
REWIND 2
C
190 READ(2,2000,END=240) ND1O, ND2O
I = I + 1
GO TO 200
C
240 IF(I.NE.0) GO TO 320
WRITE(6,2200)
2200 FORMAT(' *** No Data in Sortout File *** STOP (SUB.BILD)')
STOP
C
200 READ(2,2000,END=230) ND1N, ND2N
I = I + 1
GO TO 250
C
230 ND1N = 0
ND2N = 0
250 IF ( ND1O.EQ.ND1N .AND. ND2O.EQ.ND2N) GO TO 300
C
IAL = IAL + 1
ILY = ILY + 1
C
IRS( IAL* 2-1) = ND1O
IRS( IAL* 2 ) = ND2O
C
IGL( ILY* 2-1) = ND1O
IGL( ILY* 2 ) = ND2O
C
IF ( ND1N.LE.0) GO TO 320
ND1O = ND1N
ND2O = ND2N
GO TO 200
C
300 IAL = IAL + 1
IRS( IAL* 2-1) = ND1O
IRS( IAL* 2 ) = ND2O
GO TO 190
C
320 CONTINUE
NGL = ILY
NRS = IAL
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( PSI, PSO, VOR, VRO, TEP, TPO, U, V, NNP,
1 NEL, ICK, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICN
C
DIMENSION PSI( NNP), PSO( NNP), VOR( NNP), VRO( NNP),
1 TEP( NNP), TPO( NNP), U ( NEL), V ( NEL)
C
COMMON /ITL/ IPS, ITP, IVR
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /VVT/ VARV( 2000), VART( 2000), VARP( 2000)
C
ICK = 0
EPS = 0.1 ! 0.1(%)
C
C ------------------------------------------------
CMAX = -1.D10
DO 100 I = 1, NNP
C
CMAX = DMAX1( CMAX, DABS( PSI( I)))
100 CONTINUE
PMAX = CMAX* 0.001D0 ! 0.1(%) Up
C ------------------------------------------------
C
DMAX = -1.D10
DO 200 I = 1, NNP
C
IF ( DABS( PSI( I)).LE.PMAX ) GO TO 200
DMAX = DMAX1( DMAX, DABS(( PSI( I) - PSO( I))/PSI( I)))
C
200 CONTINUE
C ------------------------------------------------------------------
C
C ----- Attention : All PSI, PSO and VRO at t = 0 ---
C
C ==> Ww = 0 ( Not moving wall only : Natural convection )
C ==> Var_P ( at LOP = 2) ==> 100
C
IF ( LOP.EQ.1 ) DMAX = 1.D0
VARP( LOP) = DMAX* 100.D0
C
C -------------------------------------------
CMAX = -1.D10
DO 300 I = 1, NNP
C
CMAX = DMAX1( CMAX, DABS( VOR( I)) )
300 CONTINUE
VMAX = CMAX* 0.001D0 ! 0.1(%) Up
C -------------------------------------------
DMAX = -1.E10
DO 400 I = 1, NNP
C
IF ( DABS(VOR( I)).LE.VMAX) GO TO 400
DMAX = DMAX1( DMAX, DABS( ( VOR( I)-VRO( I))/ VOR( I)))
C
400 CONTINUE
VARV( LOP)= DMAX* 100.D0
C ------------------------------------------------
C
IF ( DMAX.GE.EPS .AND. ICN.EQ.0 ) GO TO 333
C
C ------------------------------------------------
CMAX = -1.D10
DO 500 I = 1, NNP
C
CMAX = DMAX1( CMAX, DABS( TEP( I)))
500 CONTINUE
TMAX = CMAX* 0.001D0 ! 0.1(%) Up
C ------------------------------------------------
DMAX = -1.D10
DO 600 I = 1, NNP
C
IF ( DABS(TEP( I)).LE.TMAX) GO TO 600
DMAX = DMAX1( DMAX, DABS( ( TEP( I) - TPO( I))/ TEP( I)))
C
600 CONTINUE
VART( LOP)= DMAX* 100.D0
C ----------------------------------------------------------
IF ( VART( LOP).LE.EPS.AND. VARV( LOP).LE.EPS ) THEN
C
C ------------------------------------------------------
C
CALL LPST ( PSI, VOR, TEP, U, V, NNP, NEL, ICN )
C
C ------------------------------------------------------
ICK = 1
WRITE(6,*) ' '
C
END IF
C ----------------------------------------------------------
C
333 CONTINUE
C -----------------------------------------------------
DO 700 I = 1, NEL
C
IF ( DABS( U( I)).GE.10.D0) GOTO 900
IF ( DABS( V( I)).GE.10.D0) GOTO 900
C
700 CONTINUE
C -----------------------------------------------------
C
IF ( ICN.EQ.1)
1WRITE(6,2000) LOP, VARP( LOP), VARV( LOP), VART( LOP)
2000 FORMAT(' * LOP =', I4,' Var_Psi=', F8.3,' (%)',' Var_V =',F8.3,
1 ' (%)',' Var_T =',F8.3,' (%)')
C -----------------------------------------------------
C
IF ( ICN.EQ.0)
1 WRITE(6,2100) LOP, VARP( LOP), VARV( LOP)
2100 FORMAT(' * LOP =', I4,' Var_P=', F8.3,' (%)',' Var_Vr=', F8.3,
1 ' (%)')
C
C -----------------------------------------
C
CALL MOVER ( PSI, 1.D0, PSO, NNP )
C
C -----------------------------------------
C
CALL MOVER ( VOR, 1.D0, VRO, NNP )
C
C -----------------------------------------
C
CALL MOVER ( TEP, 1.D0, TPO, NNP )
C
C ------------------------------------------------------
IF ( LOP.GE.LMT) THEN
C
CALL LPOV ( PSI, VOR, TEP, U, V, NNP, NEL, ICN )
C
ICK = 1
C
END IF
C ------------------------------------------------------
RETURN
C
900 WRITE(6,2200)
2200 FORMAT(/,' *** U or V .GE. 10 *** ')
C
STOP
END
C **********************************************************************
C
SUBROUTINE CLEAR8 ( X, V, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( N )
C
DO 100 I = 1, N
C
X( I) = V
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLEAR2 ( IX2, IV2, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IX2( N )
C
DO 100 I = 1, N
C
IX2 ( I) = IV2
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCNL ( AV, I, J, A, L, K, NT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( NT), L( NT)
C
C ----- Attention --- No,of Nodes.LE.32768 ---
C
IF ( I.GT.J ) RETURN
C
M = I* 200 + J
C
IE = K
N1 = 1
IF ( K.EQ.0) N2 = K+1
IF ( K.EQ.0) GOTO 10
C
N2 = K
10 N = N2
C
IF ( IE.EQ.0) IE = 1
C
C ------------------------------------------------
DO 100 II = 1, IE
C
IF ( K.EQ.0) GO TO 600
IF ( M - L( N)) 200, 300, 400
C
200 CONTINUE
N2 = N
N = N2 - ( N2-N1)/2
IF ( N.NE.N2) GO TO 100
IF ( M - L( N1)) 110, 115, 120
110 N = N1
C
120 CONTINUE
C
C -------------------------------------------
C
CALL SUBS ( AV, M, N, A, L, K, NT)
C
C -------------------------------------------
RETURN
C
115 N = N1
GO TO 300
C
550 N = N2
300 CONTINUE
A( N )= A( N ) + AV
RETURN
C
400 CONTINUE
N1 = N
N = N1 + ( N2-N1)/2
IF ( N.NE.N1) GO TO 100
IF ( M - L( N2)) 500, 550, 600
500 N = N2
GOTO 120
C
600 K = K + 1
IF ( K.GT.NT) WRITE(6,2000)
2000 FORMAT(' ** Error ** ''NT'' IN ''SOR'' IS TOO SMALL **')
C
A( K)= AV
L( K)= M
C
RETURN
C
100 CONTINUE
C ------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' ** Error --- DO Loop Strange, K =',I10)
STOP
C
END
C **********************************************************************
C
SUBROUTINE DCNIJ ( IA, JA, NI, NJ, NIX, NJX )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IA( NIX ), JA( 2,NJX )
C
CHARACTER*4 ANIX, ANJX
C
C ----- NI, NJ : Current length of IA, JA
C
C NIX, NJX : Max. length of IA, JA
C
DATA ANIX/'NIX '/ , ANJX/'NJX '/
C
C ----- To count every ROW's element ---
C
NI = 1
C
C ------------------------------------------------
DO 100 I = 1, NJ
C
IF ( NI.EQ.JA( 1,I)) GO TO 101
IA( NI)= I - 1
NI = JA(1,I)
IF ( NI.GT.NIX) GO TO 601
101 JA( 1,I)= I
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Last row's counter---
C
IA( NI )= NJ
C
C ----- To set ( I,J ) in JA for I > J ---
C
C ----------------------------------------------------------
DO 200 I = 2, NI
C
I1 = I - 1
IS = IA( I1)
C
C -----------------------------------------------
DO 300 J = 1, I1
C
JS = 1
IF ( J.GT.1) JS = IA(J-1) + 1
JL = IA(J)
C
C --------------------------------------
DO 400 K = JS, JL
C
IF ( JA(2,K).NE.I ) GO TO 400
K1 = K
GO TO 410
C
400 CONTINUE
C --------------------------------------
C
C ----- (I,J) not exist ---
C
GO TO 300
C
410 CONTINUE
C
C ----- To shift JA ---
C
NJIS = NJ - IS
C
C --------------------------------------
DO 500 L = 1, NJIS
C
JA( 1,NJ-L+2)= JA( 1,NJ-L+1)
JA( 2,NJ-L+2)= JA( 2,NJ-L+1)
C
500 CONTINUE
C --------------------------------------
C
NJ = NJ + 1
IF ( NJ.GT.NJX) GO TO 602
C
C ----- To increase IA ---
C
C --------------------------------------
DO 600 L = I, NI
C
IA( L)= IA( L) + 1
C
600 CONTINUE
C --------------------------------------
C
C ----- To insert ( I,J ) element ---
C
IS = IS + 1
JA( 1,IS)= JA( 1,K1)
JA( 2,IS)= J
C
300 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
C
601 WRITE(6,2000) ANIX, NIX, ANIX, NI
STOP
C
602 WRITE(6,2000) ANJX, NJX, ANJX, NJ
2000 FORMAT(' ** Error ( DCNIJ ) ',A3,'(',I5,') too small',A2,'=',I5)
C
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NNL/ NNP, NEL, LUV
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C H F 2 D - F E T C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version 4.3 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2D Thermal Fluid Analysis by FEM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Triangular Linear Element ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) ' '
C
READ(7,1000) NNP, NEL, LUV
1000 FORMAT(10X,I5,5X,I5,5X,I5)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, D,
1 NCW, NEL, NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION AT( NEL), D ( NNZ), NCW( NNZ), NDE( NEL,3),
1 BE( NEL,3), CE( NEL,3), DX1( NEL), DY1( NEL),
2 A1( 3,3), N ( 6), DX2( NEL), DY2( NEL)
C
CALL CLEAR8 ( D, 0.D0, NNZ )
C
NZK = NNZ
C
C ------------------------------------------------
DO 100 IEL = 1, NEL
C
B1 = BE( IEL,1)
B2 = BE( IEL,2)
B3 = BE( IEL,3)
C
C1 = CE( IEL,1)
C2 = CE( IEL,2)
C3 = CE( IEL,3)
C
N ( 1)= NDE( IEL,1)
N ( 2)= NDE( IEL,2)
N ( 3)= NDE( IEL,3)
C
IF ( NDM.EQ.1) THEN
C
W1 = 0.25D0/ AT( IEL)
W2 = W1
C
ELSE
C
C ----- Attention Dynamic Viscosity = Thermal Conductivity ---
C
IF ( NDM.EQ.2) THEN
W1 = DX1( IEL)* 0.25D0/ AT( IEL)
W2 = DY1( IEL)* 0.25D0/ AT( IEL)
C
ELSE
W1 = DX2( IEL)* 0.25D0/ AT( IEL)
W2 = DY2( IEL)* 0.25D0/ AT( IEL)
END IF
C
END IF
C
A1( 1,1)= W1* B1* B1 + W2* C1* C1
A1( 1,2)= W1* B1* B2 + W2* C1* C2
A1( 1,3)= W1* B1* B3 + W2* C1* C3
C
A1( 2,1)= W1* B2* B1 + W2* C2* C1
A1( 2,2)= W1* B2* B2 + W2* C2* C2
A1( 2,3)= W1* B2* B3 + W2* C2* C3
C
A1( 3,1)= W1* B3* B1 + W2* C3* C1
A1( 3,2)= W1* B3* B2 + W2* C3* C2
A1( 3,3)= W1* B3* B3 + W2* C3* C3
C
C --------------------------------------
DO 150 I = 1, 3
DO 150 J = 1, 3
C --------------------------------------
C
II = N ( I)
JJ = N ( J)
AV = A1( I,J)
C
CALL DCNL ( AV, II, JJ, D, NCW, NZK, NNZ )
C
150 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DSLVE ( FVC, C, TSM, NCC, KCC, NBC, NNP, NNZ,
1 NZ2, OMG, ERR, INL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, NBC, NCC
C
DIMENSION FVC( NNP), C( NNP), TSM( NNZ), KCC( NNP), NCC( 2,NZ2),
1 NBC( NNP)
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
KIT = 2* NNP
NITL = 0
900 NITL = NITL + 1
ERMX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GOTO 200
IF ( I.EQ.1) GOTO 210
MS = KCC(I-1) + 1
GOTO 220
C
210 MS = 1
220 ML = KCC( I)
IF ( MS.GT.ML) GO TO 1000
DK = FVC( I)
C
C --------------------------------------
DO 300 M = MS, ML
C
M1 = NCC( 1,M)
M2 = NCC( 2,M)
IF ( M2.EQ.I ) DSM = 1.D0/ TSM( M1)
IF ( M2.EQ.I ) GO TO 300
C
DK = DK - TSM( M1)* C( M2 )
C
300 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.ERMX) GO TO 200
ERMX = TEMP
C
200 CONTINUE
C ------------------------------------------------
C
IF ( ERMX.LE.ERR) THEN
INL = NITL
RETURN
END IF
C
IF ( NITL.GE.KIT) RETURN
IF ( NITL.GE.LMP) GO TO 900
GO TO 900
C
1000 WRITE(6,2000) I, MS, ML
2000 FORMAT ('1',' ** Error in "DSLVE" ** ',3I6)
C
RETURN
END
C **********************************************************************
C
FUNCTION DVAL1 ( U, BE, AT, NDE, NNP, NEL, I )
C
C ----- Differenciation ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE( NEL,3)
C
DIMENSION U( NNP), BE( NEL,3), AT( NEL)
C
U1 = U( NDE( I,1))
U2 = U( NDE( I,2))
U3 = U( NDE( I,3))
C
DVAL1 = ( BE( I,1)* U1 + BE( I,2 )* U2 + BE( I,3)* U3 )
1 /( 2.D0* AT( I))
C
RETURN
END
C **********************************************************************
C
SUBROUTINE POUT ( X, Y, AT, NDE, MTR, ITB, IEB, LBC,
1 NLUV, NNP, NEL, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR, ITB, IEB, LBC, ICN
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,3), AT( NEL),
1 MTR( NEL), X0 ( 3), Y0 ( 3), A ( 3),
2 ITB( NNP), IEB( NEL), LBC( NLUV)
C
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
COMMON /FLW/ NVRC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTPC, ITC( 50), TPC( 50), COF(100)
C
WRITE(8,*) ' '
WRITE(8,*) '*****************************************************'
WRITE(8,*) ' '
WRITE(8,*) ' ******* H F 2 D - F E T ( SV ) ( V.4.3) *******'
WRITE(8,*) ' '
WRITE(8,*) ' Copyright 2011 : Y. Matsuda'
WRITE(8,*) ' '
WRITE(8,*) '*****************************************************'
C
C ----AREA ---
WRITE(8,2000) NNP,NEL
2000 FORMAT(//,4X,'/----- AREA Information -----/'//
1 ' ',9X,'NNP ( NO. OF NDE )---------------',I5/
2 ' ',9X,'NEL ( NO. OF ELEMENT)--------------',I5)
C
C ----- RUN ---
C
WRITE(8,2100)
IF ( ICN.EQ.0) THEN
WRITE(8,2300)
ELSE
WRITE(8,2200)
END IF
C
WRITE(8,2400) ALP, RE
WRITE(8,2500) DT, LMT, LOP1
2100 FORMAT(//,4X,'/----- RUN Information -----/'/)
2400 FORMAT(9X,' Time Scheme Parameter =',F8.2/
1 ,9X,' Reynolds Number =',F8.2/)
2500 FORMAT(9X,' TIME Constants'/
1 ,9X,' DELT =',F8.2/
2 ,9X,' No. of LOP Limit =',I7/
3 ,9X,' No. of Print step LOP =',I7)
2300 FORMAT(9X,'FLW Computation'/)
2200 FORMAT(9X,'FLW & HEAT Computation'/)
C
C ----- CST ---
C
WRITE(8,2600)
C
C --------------------------------------
DO 100 I = 1, 3
C
GO TO ( 24, 25, 26),I
C
24 CONTINUE
WRITE(8,2700)
WRITE(8,2011) OMG ( I)
WRITE(8,2012) ERR( I)
GO TO 100
C
25 CONTINUE
WRITE(8,2010)
WRITE(8,2011) OMG ( I)
WRITE(8,2012) ERR( I)
WRITE(8,2013) DLT ( I)
GO TO 100
C
26 IF ( ICN.EQ.'0') GO TO 100
WRITE(8,2030)
WRITE(8,2011) OMG ( I)
WRITE(8,2012) ERR( I)
WRITE(8,2013) DLT ( I)
WRITE(8,2031) PRN, RA
C
100 CONTINUE
C --------------------------------------
C
2600 FORMAT(//,4X,'/----- CONSTANT Information -----/'/)
2700 FORMAT(9X,' Stream Function comp. ')
2010 FORMAT(9X,' Vorticity computation ')
2030 FORMAT(/9X,' Temperature computation ')
2011 FORMAT( 11X,'S.O.R. Acceleration factor =',F7.2)
2012 FORMAT( 11X,'S.O.R. Convergency CHCK const.=',F10.5)
2013 FORMAT(/11X,'Stationarity CHCK constant =',F8.3)
2031 FORMAT(/9X,' * Prandtle Number =',F10.4,4X,' * Rayleigh Number =',
1 F14.4)
C
C ----- NDE ---
C
WRITE(8,2014) NNP
WRITE(8,2015) (ITB( I), X( I), Y( I), I = 1, NNP )
C
2014 FORMAT(/,4X,'/----- Node Information -----/' //
1 ,9X,' No. of NDEs =',I5//
2 ,9X,3(' Node NO.',4X,'X ',4X,'Y ',5X)/)
2015 FORMAT(9X,I6,4X,2F10.3,5X,I6,4X,2F10.3,5X,I6,4X,2F10.3,5X)
C
C ----- ELEMENT ---
C
WRITE(8,2016) NEL
2016 FORMAT(/,4X,'/----- Element Information -----/'//
1 9X,'NO.OF ELEMENT =',I5 /)
WRITE(8,2017)
2017 FORMAT( 6X,'ELMT NO. NDE-I NDE-J NDE-K AREA',
1 6X,'MAT# Viscosity'/)
C
C --------------------------------------
DO 200 I = 1, NEL
C
C ----------------------------
DO 210 J = 1, 3
C
X0(J)= X( NDE(I,J))
Y0(J)= Y( NDE(I,J))
C
210 CONTINUE
C ----------------------------
C
C ----- To compute AREA ---
C
A( 1)= X0( 1)* ( Y0( 2) - Y0( 3))
A( 2)= X0( 2)* ( Y0( 3) - Y0( 1))
A( 3)= X0( 3)* ( Y0( 1) - Y0( 2))
C
D = ( A( 1) + A( 2 ) + A(3)) / 2.D00
C
K = MTR( I)
C
I1 = ITB( NDE( I,1))
I2 = ITB( NDE( I,2))
I3 = ITB( NDE( I,3))
C
WRITE(8,2018) IEB( I), I1, I2, I3, D, K, COF( K)
2018 FORMAT(8X,I5,3I9,4X,F8.4,4X,I4,F13.4)
C
AT( I)= D
C
IF ( D.LE.0.D0) WRITE (6,2021) I
IF ( D.LE.0.D0) STOP
C
200 CONTINUE
C --------------------------------------
C
C ----- WALL Node Number ---
C
WRITE(8,2019) NLUV
2019 FORMAT(/,4X,'/--- Boundary NDE Information & ',
1 'Stream Function Value ---/'//
2 9X,'No. Boundary NDEs =',I5//
3 3(' NO. NDE NO. PSI. ')/)
WRITE(8,2020) (I, ITB( LBC( I)), PSC( I), I = 1, NLUV )
2020 FORMAT(3(6X,I2,7X,I3,5X,F3.1,5X))
C
IF ( ICN.EQ.0) RETURN
C
WRITE(8,2050) NLUV
2050 FORMAT(/,4X,'/--- Boundary NDE Information & ',
1 'Temperature Value ---/'//
2 9X,'No. of Boundary NDEs =',I5//
3 3(' NO. NDE NO. TEP. ')/)
WRITE(8,2051) ( I, ITC( I), TPC( I),I = 1, NTPC )
2051 FORMAT(3(6X,I2,7X,I3,5X,F3.1,5X))
2021 FORMAT(' *** Area < 0. *** ELMT NO.=',I3,'*STOP*')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVCTN ( D, CO, F, NCC, KCC, NNP, NNZ, NZ2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NCC, KCC
C
DIMENSION D( NNZ), CO( NNP), F( NNP), NCC( 2,NZ2 ), KCC( NNP)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
JS = 1
JL = KCC( I)
C
IF ( I.GT.1) JS = KCC( I-1) + 1
FW = - F( I)
C
C --------------------------------------
DO 150 J = JS, JL
C
KPR = NCC( 1,J)
KCL = NCC( 2,J)
FW = FW + D( KPR)* CO( KCL)
C
150 CONTINUE
C --------------------------------------
C
F( I) = FW
C
100 CONTINUE
C -----------------------------------------------
C
RETURN
END
C **********************************************************************

SUBROUTINE FVOR ( VOR, NDE, AT, NEL, NNP, FVC )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION VOR( NNP), AT( NEL), FVC( NNP), NDE( NEL,3)
C
CALL CLEAR8 ( FVC, 0.D0, NNP)
C
DO 100 NE = 1, NEL
C
S = AT( NE)/ 12.D0
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
C
FVC( N1)= FVC( N1) + ( 2.D0* VOR( N1) + VOR( N2) + VOR( N3))* S
C
FVC( N2)= FVC( N2) + ( VOR( N1) + 2.D0* VOR( N2) + VOR( N3))* S
C
FVC( N3)= FVC( N3) + ( VOR( N1) + VOR( N2) + 2.D0* VOR( N3))* S
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( IBP, IBV, LBC, X, Y, NDE, BE, CE,
1 PSI, PSO, VOR, VRO, AT, U, V, NCN,
2 NTB, DX1, DY1, DX2, DY2, MTR, IBT, TEP,
4 TPO, QE, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR, IBP, IBV, IBT, LBC, NCN, NTB, ICN
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,3), AT ( NEL),
1 U ( NEL), V ( NEL), MTR( NEL), BE ( NEL,3),
2 CE ( NEL,3), DX1( NEL), DY1( NEL), PSI( NNP),
3 PSO( NNP), DX2( NEL), DY2( NEL), VOR( NNP),
4 VRO( NNP), LBC( LUV), IBP( NNP), TEP( NNP),
5 TPO( NNP), IBT( NNP), QE ( NEL), IBV( NNP),
6 NTB( NNP,10), NCN( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /ITL/ IPS, ITP, IVR
COMMON /FLW/ NVRC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTPC, ITC( 50), TPC( 50), COF(100)
C
IPS = 0
ITP = 0
IVR = 0
C
CALL CLEAR2 ( IBP, 1, NNP )
CALL CLEAR2 ( IBV, 1, NNP )
CALL CLEAR2 ( IBT, 1, NNP )
CALL CLEAR8 ( QE, 0.D0, NEL )
C
C ----- To construct element 'B' & 'C' matrices ---
C
CALL BCMAT ( X ,Y ,NDE ,BE ,CE, NNP, NEL )
C
C ----- To set initial value for PSI & VOR ---
C
CALL CLEAR8 ( PSI, 0.D0, NNP )
CALL CLEAR8 ( PSO, 0.D0, NNP )
C
CALL CLEAR8 ( VOR, 0.D0, NNP )
CALL CLEAR8 ( VRO, 0.D0, NNP )
C
CALL CLEAR8 ( TEP, 0.D0, NNP )
CALL CLEAR8 ( TPO, 0.D0, NNP )
C
C -----------------------------------------------------
DO 100 I = 1, NEL
C
U( I)= DVAL1( PSI, CE, AT, NDE, 1, NEL, I )
V( I)= - DVAL1( PSI, BE, AT, NDE, 1, NEL, I )
C
100 CONTINUE
C -----------------------------------------------------
C
IF ( NPC.LE.0) GO TO 110
C
C --------------------------------------
DO 200 I = 1, NPC
C
IBP( IPC( I))= 9999
PSI( IPC( I))= PSC( I)
C
200 CONTINUE
C --------------------------------------
C
110 IF ( NVRC.LE.0) GO TO 120
C
C --------------------------------------
DO 300 I = 1, NVRC
C
IBV( IVC( I))= 9999
VOR( IVC( I))= VRC( I)
C
300 CONTINUE
C --------------------------------------
C
120 IF ( NTPC.LE.0) GO TO 130
C
C --------------------------------------
DO 400 I = 1, NTPC
C
IBT( ITC( I))= 9999
TEP( ITC( I))= TPC( I)
C
400 CONTINUE
C --------------------------------------
C
130 IF ( LUV.LE.0) GO TO 150
C
C --------------------------------------
DO 500 I = 1 ,LUV
C
IBV( LBC( I))= 9999
C
500 CONTINUE
C --------------------------------------
C
150 CONTINUE
C
C ----- NO. OF MESH AT THE NDE ---
C
C --------------------------------------
DO 600 I = 1, NNP
C
NCN ( I)= 0
C
600 CONTINUE
C --------------------------------------
C
DO 700 I = 1 ,NEL
DO 700 J = 1 ,3
C --------------------------------------
C
NE = NDE( I,J )
NCN ( NE )= NCN( NE ) + 1
C
IF ( NCN( NE).GT.10) WRITE(6,2000) NCN( NE)
2000 FORMAT(' *** NO. OF ELEMEN AT NDE=',I5,' GT 10 *** STOP' )
IF ( NCN( NE).GT.10) STOP
NTB( NE,NCN( NE ))= I
C
700 CONTINUE
C --------------------------------------
DO 800 I = 1 ,NEL
C
DX1( I)= COF( MTR( I))
DY1( I)= COF( MTR( I))
C
800 CONTINUE
C --------------------------------------
C
CALL MOVER ( PSI, 1.D0 , PSO, NNP )
C
CALL MOVER ( VOR, 1.D0 , VRO, NNP )
C
CALL MOVER ( TEP, 1.D0 , TPO, NNP )
C
CALL PRTS ( PSI, VOR, TEP, U, V, ICN )
C
CALL DIFC ( DX1, DY1, DX2, DY2, U, V, MTR, COF, NEL, ICN )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, NDE, ITB, IEB, AT, U, V, DSH,
1 SH, NDB, NED, MTR, LBC, NNZ, NZ2, ICN )
C
C ----- INPT- DATA READ & PRINT ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, ITB, IEB, MTR, LBC, NDB, NED, ICN
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,3), AT ( NEL),
1 ITB( NNP), IEB( NEL), U ( NEL), V ( NEL),
2 MTR( NEL), LBC( LUV), DSH( LUV), SH ( LUV),
3 NDB( LUV,10), NED( LUV), IWK(10), WKX( 10)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
COMMON /FLW/ NVRC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTPC, ITC( 50), TPC( 50), COF(100)
C
NNPX = 0
NELX = 0
NVRC = 0
C
NPC = 0
NTPC = 0
NLUV = 0
C
CALL CLEAR8 ( COF, 1.D0, 100 )
C
CALL CLEAR8 ( U, 0.D0, NEL )
C
CALL CLEAR8 ( V, 0.D0, NEL )
C
C ----- Attention -----
DO 100 I = 1, NEL
C
MTR( I)= 1
C
100 CONTINUE
C ----------------------
C
RE = 0.D0
C
OMG( 1)= 1.6D0
OMG( 2)= 1.1D0
OMG( 3)= 1.1D0
C
ERR( 1)= 0.0001D0
ERR( 2)= 0.0001D0
ERR( 3)= 0.0001D0
C
DLT( 1)= 0.01D0
DLT( 2)= 0.01D0
DLT( 3)= 0.01D0
C
C ----- Time scheme parameter ---
C
ALP = 1.D0
C
C ---- Prandtle number---
C
PRN = 0.71D0
C ==================
C
C ----- NORS : Number of "Sembun" ---
C
NRS = 308
C ====== Attention ===
C
NNZ = NRS + NNP
NZ2 = NRS*2 + NNP
C
C ----- TIM ---
C
READ(7,1000) DT, LMT, LOP1
1000 FORMAT(10X,F5.0,5X,I5,5X,I5)
C
IF ( LOP1.EQ.0) LOP1 = LMT
C
C- ---- JTN ---
C
C ------------------------------------------
DO 200 I = 1, NNP
C
READ(7,1100) ITB( I), X( I), Y( I)
1100 FORMAT(10X,I5,5X,2F10.0)
C
200 CONTINUE
C -------------------------------------------
C
C ----- ELM ---
C
C ----------------------------------------------------------
DO 300 I = 1, NEL
C
READ(7,1200) IEB( I), NDE(I,1), NDE(I,2), NDE(I,3)
1200 FORMAT(10X,4I5)
C
300 CONTINUE
C ----------------------------------------------------------
C
C ----- PSI ---
C
C -----------------------------------------------------
DO 400 II = 1, 10
C
READ(7,1300) ( IWK( I), WKX( I), I = 1, 3 )
1300 FORMAT( 10X, 3(I5,5X,F10.0))
C
C -------------------------------------------
DO 450 I = 1, 3
C
IF ( IWK( I).EQ.0) GO TO 450
NPC = NPC + 1
IPC( NPC)= IWK( I)
PSC( NPC)= WKX( I)
C
450 CONTINUE
C -------------------------------------------
C
400 CONTINUE
C -----------------------------------------------------
C
C ----- TEP ---
C
C -----------------------------------------------------
DO 500 II = 1, 6
C ----- Attention -----
C
READ(7,1300) ( IWK( I), WKX( I),I = 1, 3)
C
C -------------------------------------------
DO 550 I = 1, 3
C
IF ( IWK( I).EQ.0) GO TO 550
NTPC = NTPC + 1
ITC( NTPC) = IWK( I)
TPC( NTPC) = WKX( I)
C
550 CONTINUE
C -------------------------------------------
C
500 CONTINUE
C -----------------------------------------------------
C
C ----- WAL ---
C
C ------------------------------------------------
DO 600 II = 1, 3
C
READ(7,1400) (IWK( I),I=1,10)
1400 FORMAT(10X,10I5)
C
C --------------------------------------
DO 650 I = 1, 10
C
IF ( IWK( I).EQ.0) GO TO 650
NLUV = NLUV + 1
LBC( NLUV)= IWK( I)
C
650 CONTINUE
C --------------------------------------
C
600 CONTINUE
C ------------------------------------------------
C
WRITE(6,2000)
2000 FORMAT(' *** Flow & Heat ( 1) or Flow only ( 0) ')
C
*** READ(5,*) ICN
ICN = 1
C
C ----- ATTENTION ---
C
IF ( ICN.EQ.0) THEN
C
WRITE(6,2100)
2100 FORMAT(' ',' ** Re = ?')
C
READ(5,*) RE
C
C --------------------------------------
DO 700 I = 1, NTPC
C
TPC( I)= 0.D0
C
700 CONTINUE
C --------------------------------------
GO TO 990
C
ELSE
C
WRITE(6,2200)
2200 FORMAT(' ',' ** Ra = ?')
C
*** READ(5,*) RA
RA = 710.D0
C
PRN = 0.71D0
C -----------------
C
END IF
C
990 CONTINUE
C
C ----- Attention ---
C
COF( 1)= 1.D0
C
WRITE(6,2300)
2300 FORMAT(' ',' ** DT ? & How many Loops ?')
C
** READ(5,*) DT, LMT
DT = 0.01D0
LMT = 100
C
C --------------------------------------------------------------------
C
CALL POUT ( X, Y, AT, NDE, MTR, ITB, IEB, LBC, NLUV,
1 NNP, NEL, ICN )
C
C --------------------------------------------------------------------
C
LUV = NLUV
C
C ---------------------------------------------
C
CALL BCSB ( DSH, SH, NDB, NED, LUV )
C
C ---------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DIFC ( DX1, DY1, DX2, DY2, U, V, MTR, COF, NEL, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 MTR, ICN
C
DIMENSION DX1( NEL), DY1( NEL), DX2( NEL), DY2( NEL),
1 U ( NEL), V ( NEL), MTR( NEL), COF( NEL)
C
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
C
C ----- ( Averaged ) Correction COFficient by Error Analysis ---
C
C for Triangular meshes with equal size ---
C
C --------------------------------------
DO 100 I = 1, NEL
C
C ----- Attention --- For Energy Equation --- ( COF: 1 ) ---
C
DX1( I)= COF( 1) + 0.25D0* DT* ( U( I) + V( I))** 2
DY1( I)= COF( 1) + 0.25D0* DT* ( U( I) + V( I))** 2
C
C ----- Attention --- For Vorticity Equation ---- ( = PRN ) ---
C
IF ( ICN.EQ.1 ) THEN
C
DX2( I)= PRN + 0.25D0* DT* ( U( I) + V( I))** 2
DY2( I)= PRN + 0.25D0* DT* ( U( I) + V( I))** 2
END IF
C
C ----- Attention --- For Vorticity Equation ---- ( = 1/ RE ) ---
C
IF ( ICN.EQ.0) THEN
C
DX2( I)= ( 1.D0/ RE ) + 0.25D0* DT* ( U( I) + V( I))** 2
DY2( I)= ( 1.D0/ RE ) + 0.25D0* DT* ( U( I) + V( I))** 2
END IF
C
100 CONTINUE
C ---------------------------------------
C
RETURN
C
END
C ***********************************************************************
C
SUBROUTINE CPSI ( PSI, PSO, IBP, VOR, DP, AT, BE, CE,
1 DX1, DY1, DX2, DY2, FVC, NCC, KCC, NDE,
2 U, V, MTR, COF, NNZ, NZ2, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, NDE, IBP, MTR, NCC, ICN
C
DIMENSION PSI( NNP), PSO( NNP), IBP( NNP),
1 VOR( NNP), DX1( NEL), DY1( NEL),
2 DX2( NEL), DY2( NEL), DP ( NNZ),
3 FVC( NNP), NCC( 2,NZ2), KCC( NNP),
4 U ( NEL), V ( NEL), AT ( NEL),
5 MTR( NEL), COF( NEL), NDE( NEL,3),
6 BE ( NEL,3), CE ( NEL,3)
C
COMMON /ITL/ IPS, ITP, IVR
COMMON /NNL/ NNP, NEL, LUV
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
C
C ----- To construct F & to solve equation ---
C
CALL FVOR ( VOR, NDE, AT, NEL, NNP, FVC )
C
CALL MOVER ( PSO, 1.D0, PSI, NNP )
C
CALL DSLVE( FVC, PSI, DP, NCC, KCC, IBP, NNP, NNZ,
1 NZ2, OMG(1), ERR( 1), IPS )
C
C ----- Velocity ---
C
DO 100 I = 1, NEL
C
U( I)= DVAL1( PSI, CE, AT, NDE, 1, NEL, I )
V( I)= - DVAL1( PSI, BE, AT, NDE, 1, NEL, I )
C
100 CONTINUE
C
C ----- To re-compute DX1 & DY1 ---
C
C ----------------------------------------------------------------
C
CALL DIFC ( DX1, DY1, DX2, DY2, U, V, MTR, COF, NEL, ICN )
C
C ----------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CTEP ( TEP, TPO, IBT, D, P2, AT, BE, CE,
1 U, V, QWK, FVC, QE, UN, VN, NCC,
2 KCC, NDE, NTB, NCN, LBC, LUV, NNP, NNZ,
3 NZ2, NEL, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, IBT, NDE, NTB( NNP,10), NCN, LBC, NCC, ICN
C
DIMENSION TEP( NNP), TPO( NNP), IBT( NNP), U ( NEL),
1 V ( NEL), D ( NNZ), FVC( NNP), NCC( 2,NZ2),
2 KCC( NNP), AT ( NEL), BE ( NEL,3), NDE( NEL,3),
3 CE ( NEL,3), QWK( NEL), P2 ( NNZ), UN ( NNP),
4 VN ( NNP), QE ( 1)
C
DIMENSION NCN( NNP), LBC( LUV)
C
COMMON /ITL/ IPS, ITP, IVR
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
C
CALL MOVER ( QE, 1.D0, QWK, NEL)
C
CALL MQV12 ( AT, TEP, BE, CE, FVC, UN, VN, U, V, LBC,
1 NDE, NTB, NCN, QWK, ICN )
C
CALL FVCTN ( D, TPO, FVC, NCC, KCC, NNP, NNZ, NZ2 )
C
IF ( LOP.GT.1) CALL MOVER ( TPO, 1.D0, TEP, NNP )
C
CALL DSLVE ( FVC, TEP, P2, NCC, KCC, IBT, NNP, NNZ,
1 NZ2, OMG(3), ERR(3), ITP )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPVR ( VOR, VRO, IBV, D, P2, AT, BE, CE,
1 U, V, UN, VN, FVC, NCC, KCC, NDE,
2 NTB, NCN, LBC, NNZ, NZ2, TEP, QWK, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, IBV, NDE, NTB, NCN, LBC, NCC, ICN
C
DIMENSION VOR( NNP), VRO( NNP), IBV( NNP), U ( NEL),
1 V ( NEL), D ( NNZ), FVC( NNP), NCC( 2,NZ2),
2 KCC( NNP), AT ( NEL), BE ( NEL,3), NDE( NEL,3),
3 CE ( NEL,3), P2 ( NNZ), UN ( NNP), VN ( NNP),
4 TEP( NNP), QWK( NEL), LBC( LUV)
C
DIMENSION NTB( NNP,10), NCN( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /ITL/ IPS, ITP, IVR
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
COMMON /FLW/ NVRC, IVC( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 NTPC, ITC( 50), TPC( 50), COF(100)
C
C ----- To ICSTruct F ---
C
CALL CLEAR8 ( QWK, 0.D0, NEL )
C
IF ( ICN.EQ.0) GO TO 90
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
QQWK = DVAL1( TEP, BE, AT, NDE, 1, NEL, I )
QWK( I) = QQWK* RA* PRN
C
100 CONTINUE
C ------------------------------------------------
C
90 CONTINUE
C
CALL MQV12 ( AT, VOR, BE, CE, FVC, UN, VN, U, V, LBC,
1 NDE, NTB, NCN, QWK, ICN )
C
CALL FVCTN ( D, VRO, FVC, NCC, KCC, NNP, NNZ, NZ2 )
C
C ----- To solve simultaneous equation ---
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( IBV( I).EQ.9999) GO TO 200
VOR( I)= VRO( I)
C
200 CONTINUE
C -----------------------------------------------
C
CALL DSLVE ( FVC, VOR, P2, NCC, KCC, IBV, NNP, NNZ,
1 NZ2, OMG(2), ERR(2), IVR )
C
RETURN
END
C **********************************************************************

SUBROUTINE LPOV ( PSI, VOR, TEP, U, V, NNP, NEL, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), U( NEL), V( NEL)
C
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
C
IF ( MOD( LOP, LOP1)) 120 ,110 ,120
C
110 LOP1 = LOP + 1
GO TO 100
C
120 LOP1 = LOP
C
100 CONTINUE
C
C ----------------------------------------------
C
CALL PRTS ( PSI, VOR, TEP, U, V, ICN )
C
C ----------------------------------------------
C
WRITE(6,2000)
2000 FORMAT(/,' *** Loop limit over ***')
WRITE(8,2000)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPST ( PSI, VOR, TEP, U, V, NNP, NEL, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), U( NEL), V( NEL)
C
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
C
IEND = 1
C
IF ( MOD( LOP, LOP1)) 120, 110, 120
C
110 LOP1 = LOP + 1
GO TO 100
C
120 LOP1 = LOP
100 CONTINUE
C
C ----------------------------------------------
C
CALL PRTS ( PSI, VOR, TEP, U, V, ICN )
C
C ----------------------------------------------
C
IF ( IEND.EQ.1) WRITE(6,2000) LOP
IF ( IEND.EQ.1) WRITE(8,2000) LOP
2000 FORMAT (/' *** Stationary condition reached *** Loop =', I5 )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MODD ( D1, P1, P2, AT, BE, CE, DX1, DY1,
1 DX2, DY2, NCW, NDE, NEL, NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NCW
C
DIMENSION D1 ( NNZ), P1 ( NNZ), P2 ( NNZ), AT ( NEL),
1 NCW( NNZ), BE ( NEL,3), CE ( NEL,3), NDE( NEL,3),
2 DX1( NEL), DY1( NEL), DX2( NEL), DY2( NEL)
C
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
C
C ----- To modify D1 ---
C
IF ( NDM.EQ.2) THEN
C
C ----------------------------------------------------------------------
C
CALL DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, D1, NCW, NEL,
1 NNZ, 2 )
C
C ----------------------------------------------------------------------
ELSE
C ----------------------------------------------------------------------
C
CALL DMAT ( AT, NDE, BE, CE, DX1, DY1, DX2, DY2, D1, NCW, NEL,
1 NNZ, 3 )
C
C ----------------------------------------------------------------------
END IF
C
W1 = 1.D0 / DT
C ------------------------------------------------
DO 100 I = 1, NNZ
C
P2( I)= ALP* D1( I) + P1( I)* W1
D1( I)= P2( I) - D1( I)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVEI ( IX, IA, IY, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION IX( N), IY( N)
C
DO 100 I = 1, N
C
IY ( I)= IX( I)* IA
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVER ( X, A, Y, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( N), Y( N)
C
DO 100 I = 1, N
C
Y ( I)= X( I)* A
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MQV12 ( AT, C, BE, CE, F, UN, VN, U, V,
1 LBC, NDE, NTB, NCN, QWK, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICN, NDE( NEL,3), NTB( NNP,10), NCN( NNP), LBC( LUV)
C
COMMON /NNL/ NNP, NEL, LUV
C
DIMENSION F ( NNP), U ( NEL), V( NEL), UN( NNP),
1 VN( NNP), AT( NEL), C( NNP), BE( NEL,3),
2 CE( NEL,3), QWK( NEL)
C
C ----------------------------------------------------------
DO 100 I = 1, NNP
C
UNA = 0
VNA = 0
C
UNB = 0
VNB = 0
C
K = NCN ( I)
C
C -----------------------------------------------
DO 150 N = 1, K
C
UNA = UNA + AT( NTB(I,N))* U( NTB(I,N))
VNA = VNA + AT( NTB(I,N))* V( NTB(I,N))
C
UNB = UNB + AT( NTB(I,N))
VNB = VNB + AT( NTB(I,N))
C
150 CONTINUE
C -----------------------------------------------
C
UN( I)= UNA / UNB
VN( I)= VNA / VNB
C
100 CONTINUE
C ----------------------------------------------------------
DO 200 J = 1, LUV
C
UN( LBC( J))= 0.D0
VN( LBC( J))= 0.D0
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( ICN.EQ.1 ) GO TO 333
C
C ------ Attention : ICN = 0 ---
C
C For Moving Wall --- 113/196 ---
C
C ----------------------------
DO 300 K = 106, 113
C
UN( K)= 1.D0
300 CONTINUE
C ----------------------------
C
333 CONTINUE
C
CALL CLEAR8 ( F, 0.D0, NNP )
C
C --------------------------------------
DO 400 NE = 1, NEL
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
C
U1 = UN( N1)
U2 = UN( N2)
U3 = UN( N3)
C
V1 = VN( N1)
V2 = VN( N2)
V3 = VN( N3)
C
AR = AT( NE)
C
B1 = BE( NE,1)
B2 = BE( NE,2)
B3 = BE( NE,3)
C
C1 = CE( NE,1)
C2 = CE( NE,2)
C3 = CE( NE,3)
C
CF1 = ( B1* C( N1) + B2* C( N2) + B3* C( N3)) / 24.D0
CF2 = ( C1* C( N1) + C2* C( N2) + C3* C( N3)) / 24.D0
C
CF3 = ( B1* U1 + C1* V1 + B2* U2
1 + C2* V2 + B3* U3 + C3* V3 )/ 24.D0
C
F( N1)= F( N1) - AR* QWK( NE)/3.D0
1 + CF1* ( 2.D0* U1 + U2 + U3) + CF2* ( 2.D0* V1 + V2 + V3)
C
F( N2)= F( N2) - AR* QWK( NE)/3.D0
1 + CF1* ( U1 + 2.D0* U2 + U3) + CF2* ( V1 + 2.D0* V2 + V3)
C
F( N3)= F( N3) - AR* QWK( NE)/3.D0
1 + CF1* ( U1 + U2 + 2.D0* U3) + CF2* ( V1 + V2 + 2.D0* V3)
C
C ----- ICNservato form ICNsidered ---
C
F( N1)= F( N1) + CF3* ( 2.D0* C( N1) + C( N2) + C( N3))
F( N2)= F( N2) + CF3* ( C( N1) + 2.D0* C( N2) + C( N3))
F( N3)= F( N3) + CF3* ( C( N1) + C( N2) + 2.D0* C( N3))
C
400 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VORW ( SH, DSH, PSI, VOR, NED, LBC, NDB, NNP,
1 LUV, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NED, LBC, NDB, ICN
C
DIMENSION SH ( LUV), DSH( LUV), PSI( NNP), VOR( NNP),
1 NED( LUV), LBC( LUV), NDB( LUV,10)
C
C ----- Vorticity at wall ---
C
IF ( LUV.LE.0) RETURN
C
C ------------------------------------------------
DO 100 I = 1, LUV
C
ND = LBC( I)
NE = NED( I)
N1 = NDB( I,1)
C
IF ( NE.EQ.0) VOR( ND)= 0.D0
IF ( NE.EQ.0) GO TO 100
C
APSI = PSI( N1)
AVOR = VOR( N1)
C
VOR( ND)= DSH( I)* ( PSI( ND) - APSI ) - AVOR* 0.5
C
IF ( ICN.EQ.1 ) GO TO 100 ! For Natural Convection Problem
C
C ----- Attention : For a Movig Wall --- 113 / 196 ---
C
IF ( ND.GE.107.AND.ND.LE.112 ) VOR( ND)= VOR( ND) - SH( I)
C -----------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTD ( PSI, VOR, TEP, UN, VN, X, Y, NDE, IRS, NRS,
1 ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, IRS, NRS, ICN
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /ITL/ IPS, ITP, IVR
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
COMMON /CST/ OMG(3), ERR(3), DLT(3), RA, RE, PRN
C
COMMON /VVT/ VARV( 2000), VART( 2000), VARP( 2000)
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), UN( NNP),
1 VN ( NNP), X ( NNP), Y ( NNP), NDE( NEL,3),
2 IRS( 616)
C
WRITE(10,2000) ICN, LOP, DT, RE, RA, PRN, NNP, NEL
2000 FORMAT(I2, I4, F6.3, 3F8.3, 2I4)
WRITE(10,2100) ( X( I), Y( I), I = 1, NNP)
2100 FORMAT(10F8.4)
WRITE(10,2200) ( ( NDE( I,J ), I = 1, NEL), J = 1, 3)
2200 FORMAT(20I4)
WRITE(10,2200) NRS, ( IRS( I), I=1, 616)
WRITE(10,2300) ( VARP( I), I = 1, LOP)
2300 FORMAT(10F12.5)
C
IF ( ICN.EQ.1) WRITE(10,2300) ( VARV( I), VART( I), I = 1, LOP)
IF ( ICN.EQ.0) WRITE(10,2300) ( VARV( I), I = 1, LOP)
C
WRITE(10,2100) ( UN( I), VN( I), I = 1, NNP)
C
IF (ICN.EQ.1)
1WRITE(10,2100) ( PSI( I), VOR( I), TEP( I), I = 1, NNP)
C
IF (ICN.EQ.0)
1WRITE(10,2100) ( PSI( I), VOR( I), I = 1, NNP)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMAT ( AT, NDE, P, NCC, NCW, KCC, NNZ, NZ2, NC2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, KCC, NC2
C
COMMON /NNL/ NNP, NEL, LUV
C
DIMENSION AT ( NEL), P ( NNZ), NDE( NEL,3), NCC( NZ2),
1 NCW( NNZ), KCC( NNP), A1 ( 3,3), N ( 6),
2 NC2( 2,NZ2)
C
W1D6 = 1.D0 / 6.D0
W1D12 = 1.D0 / 12.D0
C
A1 ( 1,1)= W1D6
A1 ( 1,2)= W1D12
A1 ( 1,3)= W1D12
C
A1 ( 2,1)= W1D12
A1 ( 2,2)= W1D6
A1 ( 2,3)= W1D12
C
A1 ( 3,1)= W1D12
A1 ( 3,2)= W1D12
A1 ( 3,3)= W1D6
C
CALL CLEAR8 ( P, 0.D0, NNZ )
C
NZK = 0
C
C ------------------------------------------------
DO 100 IEL = 1, NEL
C
N ( 1)= NDE( IEL,1)
N ( 2)= NDE( IEL,2)
N ( 3)= NDE( IEL,3)
C
ATI = AT ( IEL)
C
C --------------------------------------
DO 150 I = 1, 3
DO 150 J = 1, 3
C --------------------------------------
C
II = N( I)
JJ = N( J)
AV = ATI* A1( I,J)
C
CALL DCNL ( AV, II, JJ, P, NCC, NZK, NNZ )
C
150 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
DO 200 IJ = 1, NZ2
C
NC2( 1,IJ) = NCC( IJ) / 200
NC2( 2,IJ) = NCC( IJ) - 200* NC2( 1,IJ)
C
200 CONTINUE
C -----------------------------------------------------
C
CALL MOVEI ( NCC, 1, NCW, NNZ )
C
C -----------------------------------------------------
C
CALL DCNIJ ( KCC, NC2, NPWK, NZK, NNP, NZ2 )
C
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRTS ( PSI, VOR, TEP, U, V, ICN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICN
C
DIMENSION PSI( NNP), VOR( NNP), TEP( NNP), U( NEL), V( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /TMC/ DT, LMT, LOP1, LOP, ALP
C
ISW = MOD( LOP,LOP1)
C
IF ( ISW.NE.0) RETURN
C
WRITE(8,2000) LOP
WRITE(8,2100) ( I, PSI( I), I = 1, NNP)
2100 FORMAT(25('-'),'/ Stream Function /',33('-')//1X,
1 3(' NDE VALUE ')/,(1X,6(I5,1X,F10.5)))
WRITE(8,2000) LOP
WRITE(8,2200) ( I, VOR( I), I = 1, NNP)
2200 FORMAT(25('-'),'/ Vorticity /',33('-')//1X,
1 3(' NDE VALUE ')/,(1X,6(I5,1X,F10.5)))
C
IF ( ICN.EQ.0) GO TO 900
WRITE(8,2000) LOP
WRITE(8,2300) ( I, TEP( I), I = 1, NNP)
2300 FORMAT(25('-'),'/ Temperature /',33('-')//1X,
1 3(' NDE VALUE ')/,(1X,6(I5,1X,F10.5)))
C
900 CONTINUE
WRITE(8,2000) LOP
2000 FORMAT(/,64X,' * LOP =',I4)
WRITE(8,2400) ( I, U( I), V( I), I = 1, NEL)
2400 FORMAT(25('-'),'/ Velocity /',33('-')//1X,
1 3(' ELM U V ')/,(1X,3(I5,1X,2F10.5)))
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SORTM ( NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER WORK1, WORK2, DAT1, DAT2
C
C ---- Attention NEL*3 ---
C
DIMENSION DAT1( 588), DAT2( 588)
C
C ---- DATA 113 / 196 ---
C
KOS = NEL* 3
REWIND 1
C
READ(1,1000) ( DAT1( I), DAT2( I), I = 1, KOS )
1000 FORMAT(2I5)
C
C ------------------------------------------------
DO 100 I = 1, KOS-1
C
C --------------------------------------
DO 150 J = I+1, KOS
C
IF ( DAT1( I).EQ.DAT1( J)) GO TO 200
C
IF ( DAT1( I).GT.DAT1( J)) THEN
WORK1 = DAT1( I)
WORK2 = DAT2( I)
C
DAT1( I)= DAT1( J)
DAT2( I)= DAT2( J)
C
DAT1( J)= WORK1
DAT2( J)= WORK2
END IF
C
GO TO 210
C
200 CONTINUE
C
IF ( DAT2( I).GT.DAT2( J)) THEN
WORK1 = DAT2( I)
WORK2 = DAT1( I)
C
DAT2( I)= DAT2( J)
DAT1( I)= DAT1( J)
C
DAT2( J)= WORK1
DAT1( J)= WORK2
END IF
C
210 CONTINUE
C
150 CONTINUE
C --------------------------------------
C
100 CONTINUE
C -----------------------------------------------
C
WRITE(2,1000) ( DAT1( I), DAT2( I), I = 1, KOS )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SUBS ( AV, M, N, A, L, K, NT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( NT), L( NT)
C
K = K + 1
IF ( K.GT.NT) WRITE(6,2000)
2000 FORMAT(' * Error K.GT.NT *')
C
JE = K - N
C ----------------------------
DO 100 JJ = 1, JE
C
J1 = K - JJ + 1
A( J1)= A( J1-1)
L( J1)= L( J1-1)
C
100 CONTINUE
C ----------------------------
C
A( N)= AV
L( N)= M
C
RETURN
END
C **********************************************************************