戻る　

C ****************************************************************** C
C C
C Fluid Analysis Software by Finite Element Method C
C C
C ( Triangular Linear Element ) C
C C
C F 2 D - F E T C
C C
C ( Version: 3.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ****************************************************************** C
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE, INB, IEB, KC1, MTR, IBP, IBV, LBC, NED,
1 NDB, NCN, NLB, NC2, IGN, IRS, NGL, NRS
C
INTEGER*2 CON
C
C ----- DATA : 113 Nodes & 196 Elements ---
C
DIMENSION X ( 113), Y ( 113), NDE( 196,3), INB( 113),
1 IEB( 196), AT ( 196), U ( 196), V ( 196),
2 UO (196), VO ( 196), UN ( 113), VN ( 113),
3 MTR( 196), LBC( 28), P1 ( 421), P2 ( 421),
4 NC2( 2,729), KC1( 113), NCW( 421), FVC( 113)
C
DIMENSION BE ( 196,3), CE ( 196,3), DP ( 421), D1 ( 421),
1 D2 ( 421), NED( 28), DSH( 28), SH ( 28),
2 AKX( 196), AKY( 196), AX2( 196), AY2( 196),
3 PSI( 113), PSO( 113), IBP( 113),
4 VOR( 113), VRO( 113), IBV( 113)
C
DIMENSION NLB( 113,10), NDB( 28,10), NCN( 113), NC1( 729),
1 IGN( 56), IRS( 616)
C
COMMON /CNST/ OMG( 3), ERR( 3), DTA( 3), RE
C
COMMON /FLOW/ NVRC, IVRC( 50), VORC( 50), NPSC, IPSC( 50),
1 PSC( 50), COF(100)
C
COMMON /NNL / NNP, NEL, LUV
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
COMMON /ITAL/ IPS, IVR
C
C ----- Attention : < = 500 ---
C
COMMON /VRPL/ VARV( 500), VRUV( 500), VARP( 500)
C
OPEN( 1, FILE='SORTIN.DAT' )
OPEN( 2, FILE='SORTOUT.DAT' )
C
OPEN( 7, FILE='F2D_FET_7F.DAT' )
OPEN( 8, FILE='F2D_FET_OUT.DAT' )
OPEN(10, FILE='F2D_FET_FIG.DAT' )
C
CALL TTLE
C
READ(7,1000) NNP, NEL, LUV
1000 FORMAT(10X,I5,5X,I5,5X,I5)
C
CALL INPT ( X, Y, NDE, INB, IEB, AT, U, V, DSH, SH,
1 NDB, NED, MTR, LBC, NNP, NEL, LUV, NNZ, NZ2, CON )
C
CALL BILD ( NDE, IGN, IRS, NEL, NGL, NRS )
C
LOP = 0
C
CALL INIT ( IBP, IBV, LBC, X, Y, NDE, BE, CE, PSI, PSO,
1 VOR, VRO, AT, U, V, NCN, NLB, AKX, AKY, AX2,
2 AY2, MTR, CON )
C
CALL PMAT ( AT, NDE, P1, NC1, NCW, KC1, NEL, NNP, NNZ, NZ2, NC2 )
C
C ---- DP & D1 Matrix ---
C
CALL DMAT ( AT, NDE, BE, CE, AKX, AKY, AX2, AY2, DP, NCW,
1 NEL, NNZ, 1 )
C
CALL DMAT ( AT, NDE, BE, CE, AKX, AKY, AX2, AY2, D1, NCW,
1 NEL, NNZ, 2 )
C
CALL DMAT ( AT, NDE, BE, CE, AKX, AKY, AX2, AY2, D2, NCW,
1 NEL, NNZ, 3 )
C
C ----- Vorticity on a wall ---
C
CALL NBCN ( SH, DSH, PSI, VRO, NED, LBC, NDB, NNP, LUV )
C
C ***********************************************
5555 CONTINUE
C
LOP = LOP + 1
C
CALL LPSI ( PSI, PSO, IBP, VRO, DP, AT, BE, CE, AKX,
1 AKY, AX2, AY2, FVC, NC2, KC1, NDE, U, V,
2 MTR, COF, NNP, NNZ, NZ2, NEL )
C
C ----- Vorticity on a wall ---
C
CALL NBCN ( SH, DSH, PSI, VOR, NED, LBC, NDB, NNP, LUV )
C
C MODIFY D1, D2
C
CALL MODD ( D1, P1, P2, AT, BE, CE, AKX, AKY, AX2, AY2,
1 NCW, NDE, NEL, NNZ, 2 )
C
C ----- Vorticity Equation ---
C
CALL MODD ( D2, P1, P2, AT, BE, CE, AKX, AKY, AX2, AY2,
1 NCW, NDE, NEL, NNZ, 3 )
C
CALL LPVR ( VOR, VRO, IBV, D2, P2, AT, BE, CE, U, V,
1 UN, VN, FVC, NC2, KC1, NDE, NLB, NCN, LBC, LUV,
2 NNP, NNZ, NZ2, NEL, CON )
C
CALL PRSB ( PSI, VOR, U, V, NNP, NEL, CON )
C
C ---- Stationarity Check ---
C
CALL CHCK ( PSI, PSO, VOR, VRO, DTA( 3), U, V, UN, VN, UO, VO,
1 NNP, NEL, CON, *999 )
C
GO TO 5555
C *******************************************
999 CONTINUE
C
CALL FOUT ( X, Y, NDE, UN, VN, PSI, VOR, NNP, NEL, CON,
1 IRS, NRS )
C
CLOSE ( 1)
CLOSE ( 2)
C
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
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
RETURN
END
C **********************************************************************
C
SUBROUTINE BCSB5 ( DSH, SH, NDB, NED, LUV )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
C ----- Boundary condition for Cavity problem ---
C
C （ 61 Nodes / 100 Elements )
C
INTEGER*2 NDB, NED
C
DIMENSION DSH( LUV), SH( LUV), NDB( LUV,10), NED( LUV)
C
C ----- Attention NDB( LUV) --- for 61 / 100 ---
C
DO 100 I = 1, LUV
C
SH ( I) = 15.D0
DSH( I) = 75.D0
NED( I) = 4
C
100 CONTINUE
C
NED( 1) = 0
NED( 6) = 0
NED(11) = 0
NED(16) = 0
C
NDB( 1,1) = 7
NDB( 2,1) = 13
NDB( 3,1) = 24
NDB( 4,1) = 35
NDB( 5,1) = 46
NDB( 6,1) = 51
NDB( 7,1) = 46
NDB( 8,1) = 47
NDB( 9,1) = 48
NDB(10,1) = 49
C
NDB(11,1) = 55
NDB(12,1) = 49
NDB(13,1) = 38
NDB(14,1) = 27
NDB(15,1) = 16
NDB(16,1) = 11
NDB(17,1) = 16
NDB(18,1) = 15
NDB(19,1) = 14
NDB(20,1) = 13
C
RETURN
END
C *********************************************************************
C
SUBROUTINE BCSB ( DSH, SH, NDB, NED, LUV )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
C ----- New Boundary Condition ---
C
INTEGER*2 NDB, NED
C
DIMENSION NDB( LUV,10), NED( LUV), DSH( LUV), SH( LUV)
C
C ----- For 113 / 196 ---
C
DO 100 I = 1, LUV
C
SH ( I) = 21.D0
DSH( I) = 147.D0
NED( I) = 4
C
100 CONTINUE
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, IGN, IRS, NEL, NGL, NRS )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
C ----- Outer line & Mesh line ---
C
INTEGER*2 NDE, IGN, IRS, NGL
C
DIMENSION NDE( NEL,3), IGN( 56), IRS( 616), KZ5( 2,3)
C
DATA KZ5 / 1,2, 2,3, 3,1 /
C
C ----- SORT ---
C
C -------------------------------------
DO 100 I = 1, NEL
C
DO 100 J = 1, 3
C
N1 = NDE( I,KZ5( 1,J))
N2 = NDE( I,KZ5( 2,J))
C
IF( N1.LE.N2) GO TO 150
C
N3 = N1
N1 = N2
N2 = N3
C
150 CONTINUE
C
WRITE(1,2000) N1, N2
2000 FORMAT(2I5)
C
100 CONTINUE
C --------------------------------------
C
CALL SORTM ( NEL )
C
C ----- To set IGN, IRS ---
C
NGL = 0
NRS = 0
IAL = 0
C
IOY = 0
I = 0
C
C ----- IGN : OUTLINE / IRS : Component for Meshes ---
C
REWIND 2
C
190 READ(2,2000,END=240) ND1O, ND2O
I = I + 1
C
GO TO 200
C
240 IF( I.NE.0) GO TO 320
C
WRITE(6,2200)
2200 FORMAT(/,' *** NO RECORD IN SORTOUT FILE *** STOP (BILD)')
C
STOP
C
200 READ(2,2000,END=230) ND1N, ND2N
C
I = I + 1
C
GO TO 250
C
230 ND1N = 0
ND2N = 0
C
250 IF ( ND1O.EQ.ND1N .AND. ND2O.EQ.ND2N ) GO TO 300
C
IAL = IAL + 1
IOY = IOY + 1
C
IRS( IAL*2-1) = ND1O
IRS( IAL*2 ) = ND2O
IGN( IOY*2-1) = ND1O
IGN( IOY*2 ) = ND2O
C
IF( ND1N.LE.0) GO TO 320
C
ND1O = ND1N
ND2O = ND2N
C
GO TO 200
C
300 IAL = IAL + 1
C
IRS( IAL*2-1) = ND1O
IRS( IAL*2 ) = ND2O
C
GO TO 190
C
320 CONTINUE
C
NGL = IOY
NRS = IAL
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( PSI, PSO, VOR, VRO, EPS, U, V, UN, VN,
1 UO, VO, NNP, NEL, CON, * )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 CON
C
DIMENSION PSI( NNP), PSO( NNP), VOR( NNP), VRO( NNP),
1 U ( NEL), V ( NEL), UN ( NNP), VN ( NNP),
2 UO ( NEL), VO ( NEL)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
COMMON /ITAL/ IPS, IVR
C
C ******** Attention **** <= 500 ****************
C
COMMON /VRPL/ VARV( 500), VRUV( 500), VARP( 500)
C
EPS = 0.01D0 ! (%)
C
ICHK = 0
C
C ----------------------------
VRMX = -100.D0
DO 100 I = 1, NNP
VRMX = DMAX1( VRMX, DABS( PSI( I)) )
100 CONTINUE
VMX = VRMX* 0.01D0 ! 1(%) Up
C ----------------------------
C
DMAX = -100.D0
DO 200 I = 1, NNP
C
IF ( DABS( PSI( I) ).LE.VMX ) GO TO 200
C
DMAX = DMAX1( DMAX, DABS( (PSI( I)-PSO( I) ) / PSI( I) ) )
200 CONTINUE
VARP( LOP) = DMAX* 100.D0
C ----------------------------
C
VRMX = -100.D0
DO 300 I = 1, NNP
VRMX = DMAX1( VRMX, DABS( VOR( I)) )
300 CONTINUE
VMX = VRMX * 0.01D0 ! 1(%) Up
C ----------------------------
C
DMAX = -100.D0
DO 400 I = 1, NNP
IF ( DABS( VOR( I) ).LE.VMX ) GO TO 400
DMAX = DMAX1 ( DMAX, DABS( (VOR( I) - VRO( I) ) / VOR( I) ) )
400 CONTINUE
VARV( LOP) = DMAX* 100.D0
C ----------------------------
C
VMIN = 0.01D0 ! 1(%) Up
DEFR = -1000.D0
DO 600 I = 1, NEL
IF ( DABS( U( I) ).LE.VMIN ) GO TO 600
DEFR = DMAX1( DEFR, DABS( ( U( I) - UO( I) ) / U( I) ) )
600 CONTINUE
VRUV( LOP) = DEFR* 100.D0
C ----------------------------
C
WRITE(6,2000) LOP, VARP( LOP), VARV( LOP), VRUV( LOP)
2000 FORMAT(' * LOP =', I4,' Var_P=', F8.3,' (%)',' Var_Vr=', F8.3,
1 ' (%)',' Var_Vel=',F8.3,' (%)')
C
C ----- Stationarity Check ---
C
IF ( VARV( LOP).LT.EPS ) THEN
C
WRITE(6,*) ' '
WRITE(6,*) ' ** Stationary condition reached at LOP =', LOP
C
RETURN 1
C
END IF
C
C -------------------------------
C
IF ( LOP.GE.LPLM ) THEN
C
CALL LPOV ( PSI, VOR, U, V, UN, VN )
C
RETURN 1
END IF
C
C -------------------------------------
DO 700 I = 1, NNP
C
IF ( DABS( UN( I)).GE.2.D0) GO TO 900
IF ( DABS( VN( I)).GE.2.D0) GO TO 900
C
700 CONTINUE
C --------------------------------------
C
CALL MOVER ( PSI, 1.D0, PSO, NNP )
CALL MOVER ( VOR, 1.D0, VRO, NNP )
C
CALL MOVER ( U, 1.D0, UO, NEL )
CALL MOVER ( V, 1.D0, VO, NEL )
C
RETURN
C
900 WRITE(6,2100)
2100 FORMAT(/,' *** U. or V .GE.2 * STOP ***')
C
STOP
END
C *********************************************************************
C
SUBROUTINE CPRT ( C, MSG, NNP )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
CHARACTER*4 MSG
C
DIMENSION C( NNP), MSG( 4)
C
WRITE(8,2000) MSG, ( I, C( I),I = 1, NNP )
C
2000 FORMAT(26('-'),'/ ',4A4,' /',34('-')//1X,
1 5('NDE C-VALUE ',3X)/(1X,5( I3,1X,F9.5,3X)))
C
RETURN
END
C
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
C
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
C
IF ( K.EQ.0) N2 = K+1
IF ( K.EQ.0) GO TO 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
C
IF ( M-L( N)) 200, 300, 400
C
200 CONTINUE
C
N2 = N
N = N2 - ( N2-N1 )/2
C
IF ( N.NE.N2) GO TO 100
C
IF ( M-L( N1)) 110, 115, 120
C
110 N = N1
C
120 CONTINUE
C
CALL SUBS ( AV, M, N, A, L, K, NT )
C
RETURN
C
115 N = N1
C
GO TO 300
C
550 N = N2
C
300 CONTINUE
C
A( N) = A( N) + AV
C
RETURN
C
400 CONTINUE
C
N1 = N
N = N1 + ( N2-N1)/2
C
IF ( N.NE.N1) GO TO 100
C
IF ( M-L( N2)) 500, 550, 600
C
500 N = N2
GO TO 120
C
600 K = K + 1
C
IF (K.GT.NT) WRITE(6,2000)
C
2000 FORMAT(' ** ERR ** ''NT'' IN ''SOR'' IS TOO SMALL **')
C
A( K) = AV
L( K) = M
C
RETURN
C
100 CONTINUE
C --------------------------------------
C
WRITE(6,2100) K
C
2100 FORMAT(' * Error --- DO Loop Strange, K =',I10 )
C
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
C
C ... SET IA( NI)
C
IA( NI) = I - 1
NI = JA(1,I)
C
IF ( NI.GT.NIX ) GO TO 601
C
101 JA(1,I) = I
C
100 CONTINUE
C --------------------------------------
C
C ----- To set Last Row's counter ---
C
IA( NI) = NJ
C
C ----- To set ( I,J) in JA for I > J ---
C
C --------------------------------------
C
DO 200 I = 2, NI
C
I1 = I - 1
IS = IA( I1)
C
C -----------------------------
C
DO 300 J = 1, I1
C
JS = 1
C
IF ( J.GT.1 ) JS = IA(J-1) + 1
C
JL = IA( J)
C
C ----------------------
C
DO 400 K = JS, JL
C
IF ( JA(2,K).NE.I ) GO TO 400
C
K1 = K
C
GO TO 410
C
400 CONTINUE
C -----------------------
C
C ----- To set ( I,J) not exist ---
C
GO TO 300
C
410 CONTINUE
C
C ... SHIFT JA
C
NJIS = NJ - IS
C
C -----------------------
C
DO 500 L = 1, NJIS
C
JA(1,NJ-L+2) = JA(1,NJ-L+1)
C
JA(2,NJ-L+2) = JA(2,NJ-L+1)
C
500 CONTINUE
C -----------------------
C
NJ = NJ + 1
C
IF ( NJ.GT.NJX ) GO TO 602
C
C ----- To increase IA ---
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
C
STOP
C
602 WRITE(6,2000) ANJX, NJX, ANJX, NJ
2000 FORMAT(' ** ERR (DCNIJ) ',A3,'(',I5,') TOO SMALL ',A2,'=',I5)
C
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
WRITE(6,*) ' '
WRITE(6,*) ' Fluid Analysis Software by FEM '
WRITE(6,*) ' '
WRITE(6,*) ' F 2 D - F E T '
WRITE(6,*) ' '
WRITE(6,*) ' ( Modified Galerkin Method - Triangular Linear ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Version: 3.0 ) '
WRITE(6,*) ' '
WRITE(6,*) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE DMAT( AT, NDE, BE, CE, AKX, AKY, AX2, AY2,
1 D, NCCW, NEL, NNZ, NDM )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE
C
DIMENSION AT( NEL), D( NNZ), NCCW( NNZ), NDE( NEL,3),
1 BE( NEL,3), CE( NEL,3), AKX( NEL), AKY( NEL),
2 A1(3,3), N(6), AX2( NEL), AY2( NEL)
C
CALL CLEAR8 ( D, 0.D0 , NNZ )
C
NZWK = NNZ
C
DO 100 IEL = 1, NEL
C
B1 = BE ( IEL,1 )
B2 = BE ( IEL,2 )
B3 = BE ( IEL,3 )
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
C
W1 = AKX( IEL)* 0.25D0 / AT( IEL)
W2 = AKY( IEL) * 0.25D0 / AT( IEL)
C
ELSE
C
W1 = AX2( IEL) * 0.25D0 / AT( IEL)
W2 = AY2( IEL) * 0.25D0 / AT( IEL)
C
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
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
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
DO 150 I = 1, 3
C
DO 150 J = 1, 3
C
II = N ( I)
JJ = N ( J )
C
AV = A1( I,J )
C
CALL DCNL ( AV, II, JJ, D, NCCW, NZWK, NNZ )
C
150 CONTINUE
C
100 CONTINUE
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),
1 NCC(2,NZ2), NBC( NNP)
C
C---- S.O.R. METHOD ---------------------------------------------
C
CMAX = 0.D0
C
DO 100 I = 1, NNP
C
AC = DABS (C( I))
C
IF ( AC.GE.CMAX ) CMAX = AC
C
100 CONTINUE
C
EPS = CMAX * 0.01D0
** EPS = CMAX * 0.001D0
C
LIMPR = 3 * NNP
KITAL = 2 * NNP
NITAL = 0
C
900 NITAL = NITAL + 1
C
ERMAX = 0.D0
C
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 200
C
IF ( I.EQ.1) GO TO 210
C
MS = KCC( I-1) + 1
C
GO TO 220
C
210 MS = 1
C
220 ML = KCC( I)
C
IF ( MS.GT.ML) GO TO 1000
C
DK = FVC( I)
C
DO 300 M = MS, ML
C
M1 = NCC(1,M)
C
M2 = NCC(2,M)
C
IF ( M2.EQ.I) DSM = 1.D0 / TSM( M1)
C
IF ( M2.EQ.I) GO TO 300
C
DK = DK - TSM( M1)* C( M2)
C
300 CONTINUE
C
DC = - C( I) + DSM * DK
C( I) = C( I) + OMG * DC
TEMP = DABS(C( I))
C
IF (TEMP.LE.EPS) GO TO 200
C
TEMP = DABS ( DC/C( I) )
C
IF (TEMP.LE.ERMAX) GO TO 200
C
ERMAX = TEMP
C
200 CONTINUE
C
IF (ERMAX.LE.ERR) INL = NITAL

IF (ERMAX.LE.ERR) RETURN
C
IF ( NITAL.GE.KITAL) RETURN
C
IF ( NITAL.GE.LIMPR) GO TO 900
C
GO TO 900
C
1000 WRITE (6,2000) I, MS, ML
C
2000 FORMAT ('1',' ** ERR IN " DSLVE " ** ',3I6)
C
RETURN
END
C *********************************************************************
C
FUNCTION DVAL1 ( U, BE, AT, NDE, NNP, NEL, I )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
C---- BIBUNCHI NO KEISAN ----------
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 EPRNT ( X, Y, AT, NDE, MTR, INB, IEB,
1 LBC, NLUV, NNP, NEL, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE, MTR, INB, IEB, LBC, CON
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,3), AT( NEL), MTR( NEL),
1 X0( 3), Y0( 3), A( 3), INB( NNP), IEB( NEL),
2 LBC( NLUV)
C
COMMON /CNST / OMG( 3), ERR( 3), DTA( 3), RE
COMMON /FLOW / NVRC, IVRC( 50), VORC( 50), NPSC, IPSC( 50),
1 PSC( 50), COF(100)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
C--- HEADING ----
C
WRITE(8,2001)
WRITE(8,2000)
WRITE(8,2003) NNP, NEL
WRITE(8,2000)
WRITE(8,2001)
C
2000 FORMAT(/)
2001 FORMAT(10X, 60('*'))
2003 FORMAT(10X,'*',11X,'===== F A L S - MT (V.2.5) =====',
1 11X,'*'/,10X,'*',58X,'*'/,10X,'*',11X,'Thermal Fluid Analysis ( '
2 ,I3,' /',I3,')', 13X,'*'/,10X,'*', 58X,'*'/,10X,'*',
5 10X,' Copyright 2010 : Ysuhiro MATSUDA','****')
C
C---AREA ----
C
WRITE(8,2004) NNP,NEL
2004 FORMAT(//,4X,'/----- AREA Information -----/'/' ',9X,'NNP ( NO. OF
1 NDE )---------------',I5/' ',9X,'NEL ( NO. OF ELEMENT )--
2----------',I5)
C
C-- RUN ----
C
WRITE(8,2005)
C
IF (CON .EQ. 0) THEN
C
WRITE(8,2040)
C
ELSE
C
WRITE(8,2041)
C
END IF
C
WRITE(8,2006) ALP,RE
WRITE(8,2007) DT,LPLM,NPR
C
2005 FORMAT(//,4X,'/----- RUN Information -----/'/)
2006 FORMAT(9X,' TIME SCHEME PARAMETER =',F8.2/
1 ,9X,' REYNOLDS NO. =',F8.2/)
2007 FORMAT(9X,' TIME CONSTANTS'/,9X,
1 ' DELT =',F8.2/,9X,' NO. OF LOO
2P LIMIT =',I7/,9X,' NO. OF PRINT STEP LOP =',I7)

2040 FORMAT(9X,'FLOW COMPUTATION'/)
2041 FORMAT(9X,'FLOW & HEAT COMPUTATION'/)
C
C---- CNST ----
C
WRITE(8,2008)
C
DO 100 I = 1, 3
C
GO TO ( 24, 25, 26 ), I
C
24 CONTINUE
C
WRITE(8,2009)
WRITE(8,2011) OMG( I)
WRITE(8,2012) ERR( I)
C
GO TO 100
C
25 CONTINUE
C
WRITE(8,2010)
WRITE(8,2011) OMG( I)
WRITE(8,2012) ERR( I)
WRITE(8,2013) DTA( I)
C
GO TO 100
C
26 CONTINUE
C

100 CONTINUE
C
2008 FORMAT(//,4X,'/----- CONSTANT Information -----/'/)
2009 FORMAT(9X,'NAGARE KANSU KEISAN ')
2010 FORMAT(9X,'UZUDO KEISAN ')
2011 FORMAT( 11X,'S.O.R. ACCELERATION FACTOR =',F7.2)
2012 FORMAT( 11X,'S.O.R. SYUSOKU HANTEI =',F10.5)
2013 FORMAT(/11X,'HITEIJYO SYUSOKU HANTEI =',F8.3)
C
C--- NDE ----
C
WRITE(8,2014) NNP
C
WRITE(8,2015) ( INB( I),X( I),Y( I),I=1,NNP)
C
2014 FORMAT(/,4X,'/----- NDE Information -----/' //
1 ,9X,' NO.OF NDES =',I5//
2 ,9X,3(' NDE NO.',4X,'X ( M)',4X,'Y ( M)',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 /)
C
WRITE(8,2017)
2017 FORMAT( 6X,'ELMT NO. NDE-I NDE-J NDE-K AREA',
1 6X,'MAT# VISCOSITY'/)
C
DO 200 I=1, NEL
C
DO 210 J=1, 3
C
X0( J) = X( NDE( I,J))
Y0( J) = Y( NDE( I,J))
C
210 CONTINUE
C
C--- TO CALCULATE 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))
D = (A( 1)+A( 2)+A( 3)) / 2.D0
C
K=MTR( I)
C
I1 = INB( NDE( I,1))
I2 = INB( NDE( I,2))
I3 = INB( 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) WRITE (6,2021) I
C
IF (D.LE.0) STOP
C
200 CONTINUE
C
C--- WALL NDE NUMBER ----
C
WRITE(8,2019) NLUV
2019 FORMAT(/,4X,'/--- Boundary NDE Information & ',
1 'STREAM FUNCTION VALUE ---/'//9X,'NO. OF Boundary NDES =
2',I5//3(' NO. NDE NO. PSI. ')/)
C
WRITE(8,2020) ( I, INB( LBC( I)),PSC( I), I = 1, NLUV )
C
2020 FORMAT(3(6X,I2,7X,I3,5X,F3.1,5X))
C
2021 FORMAT(' *** AREA < 0. *** ELMT NO.=',I3,'*STOP*')
C
RETURN
END
C *********************************************************************
C
SUBROUTINE FOUT ( X, Y, NDE, UN, VN, PSI, VOR, NNP, NEL,
1 CON, IRS, NRS )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE, IRS, NRS

INTEGER*2 CON
C
C ********** Attention *** <= 500 ************
C
COMMON /VRPL/ VARV( 500), VRUV( 500), VARP( 500)
C
COMMON /CNST/ OMG( 3), ERR( 3), DTA( 3), RE
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
DIMENSION X( 113), Y( 113), NDE( 196,3), IRS(616), UN( 113),
1 VN( 113), PSI( 113), VOR( 113)
C
C ------------------------------
C
WRITE(10,2040) CON, LOP, DT, RE, NNP, NEL
2040 FORMAT( I2, I4, F6.3, F8.1, 2I4)
C
WRITE(10,2050) ( X( I), Y( I), I=1, NNP )
2050 FORMAT(10F8.4)
C
WRITE(10,2060) ( ( NDE( I,J), I=1, NEL ), J=1, 3 )
2060 FORMAT(20I4)
C
WRITE(10,2060) NRS, ( IRS( I), I=1, 616 )
C
WRITE(10,2050) ( VRUV( I), VARP( I), I = 1, LOP )
WRITE(10,2050) ( VARV( I), I = 1, LOP )
C
WRITE(10,2050) ( UN( I), VN( I), I=1, NNP )
WRITE(10,2050) ( PSI( I), VOR( I), I=1, NNP )
C
WRITE(6,*) ' '
WRITE(6,*) ' **** End ****'
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
DO 100 I = 1, NNP
C
JS = 1
JL = KCC( I)
C
IF ( I.GT.1) JS= KCC( I-1) + 1
C
FW = - F( I)
C
DO 150 J=JS, JL
C
KPTR = NCC(1,J)
KCOL = NCC(2,J)
C
FW = FW + D( KPTR ) * CO( KCOL )
C
150 CONTINUE
C
F( I) = FW
C
100 CONTINUE
C
RETURN
END
C *********************************************************************
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)
1 + (2.D0 * VOR( N1) + VOR( N2) + VOR( N3) ) * S
FVC( N2) = FVC( N2)
1 + (VOR( N1) + 2.D0 * VOR( N2) + VOR( N3) ) * S
FVC( N3) = FVC( N3)
1 + (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, PSI, PSO,
1 VOR, VRO, AT, U, V, NCN, NLB, AKX, AKY,
2 AX2, AY2, MTR, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE, MTR, IBP, IBV, LBC, NCN, NLB, CON
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), AKX( NEL), AKY( NEL), PSI( NNP),
3 PSO( NNP), AX2( NEL), AY2( NEL),
4 VOR( NNP), VRO( NNP), LBC( LUV), IBP( NNP),
5 IBV( NNP), NLB( NNP,10), NCN( NNP)
C
COMMON /FLOW / NVRC, IVRC( 50), VORC( 50),
1 NPSC, IPSC( 50), PSC( 50), COF(100)
C
COMMON /NNL / NNP, NEL, LUV
C
COMMON / ITAL / IPS, IVR
C
IPS = 0
IVR = 0
C
C ... CLEAR WORK AREA
C
CALL CLEAR2 ( IBP, 1, NNP )
CALL CLEAR2 ( IBV, 1, NNP )
C
C ... CONSTRUCT THE ELEMENT 'B' & 'C' MATRICES
C
CALL BCMAT ( X, Y, NDE, BE, CE, NNP, NEL )
C
C ... SET INITIAL VALUE FOR PSI & VOR
C
CALL CLEAR8 ( PSI, 0.D0, NNP )
CALL CLEAR8 ( PSO, 0.D0, NNP )
CALL CLEAR8 ( VOR, 0.D0, NNP )
CALL CLEAR8 ( VRO, 0.D0, NNP )
C
DO 90 I = 1 ,NEL
C
U( I) = DVAL1( PSI, CE, AT, NDE, 1, NEL, I )
C
V( I) = - DVAL1( PSI, BE, AT, NDE, 1, NEL, I )
C
90 CONTINUE
C
IF ( NPSC.LE.0 ) GO TO 110
C
DO 100 I = 1, NPSC
C
IBP( IPSC( I)) = 9999
PSI ( IPSC( I)) = PSC( I)
C
100 CONTINUE
C
110 IF ( NVRC.LE.0 ) GO TO 120
C
DO 200 I = 1, NVRC
C
IBV( IVRC( I)) = 9999
VOR ( IVRC( I)) = VORC( I)
C
200 CONTINUE
C
120 CONTINUE
C
IF ( LUV.LE.0 ) GO TO 150
C
DO 400 I = 1, LUV
C
IBV( LBC( I) ) = 9999
C
400 CONTINUE
C
150 CONTINUE
C
C---- NO. OF MESH AT THE NDE
C
DO 500 I = 1, NNP
C
NCN ( I) = 0
C
500 CONTINUE
C
DO 600 I = 1, NEL
C
DO 600 J = 1, 3
C
NE = NDE( I, J )
C
NCN ( NE ) = NCN( NE ) + 1
C
IF ( NCN( NE).GT.10 ) WRITE(6,2000) NCN( NE)
C
2000 FORMAT(' *** NO. OF ELEMEN AT NDE=',I5,' GT 10 *** STOP' )
C
IF ( NCN( NE).GT.10) STOP
C
NLB( NE,NCN( NE)) = I
C
600 CONTINUE
C
DO 700 I = 1, NEL
C
AKX( I) = COF ( MTR( I) )
AKY( I) = COF ( MTR( I) )
C
700 CONTINUE
C
CALL MOVER ( PSI, 1.D0, PSO, NNP )
CALL MOVER ( VOR, 1.D0, VRO, NNP )
C
CALL PRSB ( PSI, VOR, U, V, NNP, NEL, CON )
C
CALL CDIF ( AKX, AKY, AX2, AY2, U, V, MTR, COF, NEL )
C
RETURN
END
C *********************************************************************
C
SUBROUTINE INPT ( X, Y, NDE, INB, IEB, AT, U, V, DSH, SH,
1 NDB, NED, MTR, LBC, NNP,NEL,LUV, NNZ,
2 NZ2, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
C--------------------------------------------------------------------C
C INPT- DATA READ & PRINT C
C
INTEGER*2 NDE, INB, IEB, MTR, LBC, NDB, NED, CON
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,3), AT( NEL), INB( NNP),
1 IEB( NEL), U( NEL), V( NEL), MTR( NEL), LBC( LUV),
2 DSH( LUV), SH( LUV), NDB( LUV,10), NED( LUV),
3 IWK(10), WKX(10)
C
COMMON /CNST / OMG( 3), ERR( 3), DTA( 3), RE
COMMON /TIMC / DT, LPLM, NPR, LOP, ALP

COMMON /FLOW / NVRC, IVRC( 50), VORC( 50),
1 NPSC, IPSC( 50), PSC( 50), COF(100)
C
C--- CLEAR ----------
C
NNPX = 0
NELX = 0
NVRC = 0
C
NPSC = 0
NLUV = 0
C
CALL CLEAR8 ( COF, 1.D0, 100 )
CALL CLEAR8 ( U, 0.D0, NEL )
CALL CLEAR8 ( V, 0.D0, NEL )
C
C---- Attention ---------------
C
DO 50 I = 1, NEL
50 MTR( I) = 1
C
C--- CONST. VALUE -----------
C
OMG( 1) = 1.6D0
OMG( 2) = 1.1D0
OMG( 3) = 1.1D0

ERR( 1) = 0.0001D0
ERR( 2) = 0.0001D0
ERR( 3) = 0.0001D0

DTA( 1) = 0.01D0
DTA( 2) = 0.01D0
DTA( 3) = 0.01D0

GRAV = 1000000.D0
C
C ----- Time Scheme Parameter ---
C
ALP = 1.D0
C
C ----- NORS: Number of Sennbunn ---
C
NRS = 308
C
NNZ = NRS + NNP
NZ2 = NRS* 2 + NNP
C
C
C ----- TIM ---
C
READ(7,1001) DT, LPLM, NPR
1001 FORMAT(10X,F5.0,5X,I5,5X,I5)
C
IF ( NPR.EQ.0) NPR = LPLM
C
C ----- JTN ---
C
DO 100 I = 1, NNP
C
READ(7,1003) INB( I), X( I), Y( I)
1003 FORMAT(10X,I5,5X,2F10.0)
C
100 CONTINUE
C
C ----- ELM ---
C
DO 110 I = 1, NEL
C
READ(7,1004) IEB( I), NDE( I,1), NDE( I,2), NDE( I,3)
1004 FORMAT(10X,4I5)
C
110 CONTINUE
C
C ----- PSI ---
C
DO 120 II = 1, 10
C
READ(7,1005) ( IWK( I), WKX( I), I=1,3 )
1005 FORMAT(10X,3( I5,5X,F10.0))
C
DO 125 I = 1, 3
C
IF ( IWK( I).EQ.0 ) GO TO 125
C
NPSC = NPSC+1
IPSC( NPSC) = IWK( I)
PSC ( NPSC) = WKX( I)
C
125 CONTINUE
C
120 CONTINUE
C
C ----- WAL---
C
DO 150 II = 1, 3
C
READ(7,1006) ( IWK( I), I=1, 10 )
1006 FORMAT(10X,10I5)
C
DO 160 I = 1, 10
C
IF ( IWK( I).EQ.0 ) GO TO 160
C
NLUV = NLUV+1
LBC( NLUV)= IWK( I)
C
160 CONTINUE
C
150 CONTINUE
C
CON = 0
C
WRITE(6,*) ' '
WRITE(6,2000)
2000 FORMAT(' ',' Re = ?')
C
** READ(5,*) RE
RE = 100.D0
C
COF( 1) = 1.D0 / RE
C
WRITE(6,*) ' '
WRITE(6,2100)
2100 FORMAT(' ',' DELT = ? LOP = ? ')
WRITE(6,*) ' '
C
** READ(5,*) DT, LPLM
DT = 0.1D0
LPLM = 150
C
WRITE(6,2200) RE, DT, LPLM
2200 FORMAT(' *** Re =',F7.1,' DT =',F6.3,' LPLM =', I4,/)
C
C ----- To print ---
C
CALL EPRNT ( X, Y, AT, NDE, MTR, INB, IEB, LBC, NLUV, NNP, NEL,
1 CON )
C
LUV = NLUV
C
CALL BCSB ( DSH, SH, NDB, NED, LUV )
C
RETURN
END
C *********************************************************************
C
SUBROUTINE CDIF ( AKX, AKY, AX2, AY2, U, V, MTR, COF, NEL )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 MTR
C
DIMENSION AKX( NEL), AKY( NEL), AX2( NEL), AY2( NEL),
1 U( NEL), V( NEL), COF( NEL), MTR( NEL)
C
COMMON / TIMC/ DT, LPLM, NPR, LOP, ALP
COMMON / CNST / OMG( 3), ERR( 3), DTA( 3), RE
C
DO 100 I = 1, NEL
C
C-----Attention---- MODIFIED THERMAL DIFFUSIVITY ( = 1 / ( RE ) )
C
AKX( I) = COF( MTR( I)) + 0.25D0 * DT * ( U( I) + V( I) )** 2
AKY( I) = COF( MTR( I)) + 0.25D0 * DT * ( U( I) + V( I) )** 2
C
C-----Attention---- MODIFIED DYNAMIC VISCOSITY (= 1 / RE )
C
AX2( I) = COF( MTR( I)) + 0.25D0 * DT * ( U( I) + V( I) )** 2
AY2( I) = COF( MTR( I)) + 0.25D0 * DT * ( U( I) + V( I) )** 2
C
100 CONTINUE
C
RETURN
END
C *********************************************************************
C
SUBROUTINE LPSI ( PSI, PSO, IBP, VOR, DP, AT, BE, CE,
1 AKX, AKY, AX2, AY2, FVC, NCC, KCC, NDE,
2 U, V, MTR, COF, NNP, NNZ, NZ2, NEL )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 KCC, NDE, IBP, MTR, NCC
C
DIMENSION PSI( NNP), PSO( NNP), IBP( NNP), VOR( NNP),
1 AKX( NEL), AKY( NEL), DP ( NNZ), FVC( NNP),
2 NCC(2,NZ2), KCC( NNP), U ( NEL), V ( NEL),
3 AT ( NEL), MTR( NEL), COF( NEL), NDE( NEL,3),
4 BE ( NEL,3), CE ( NEL,3), AX2( NEL), AY2( NEL)
C
COMMON /CNST/ OMG( 3), ERR( 3), DTA( 3), RE
COMMON /ITAL/ IPS, IVR
C
C ----- To construct F & solve an 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 ...TO CALCULATE 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 recompute AKX & AKY ---
C
CALL CDIF ( AKX, AKY, AX2, AY2, U, V, MTR, COF, NEL )
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 NLB, NCN, LBC, LUV, NNP, NNZ, NZ2, NEL,
3 CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 KCC, IBV, NDE, NLB( NNP,10), NCN( NNP), LBC( LUV),
1 NCC, CON
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)
C
COMMON /ITAL/ IPS, IVR
COMMON /CNST/ OMG( 3), ERR( 3), DTA( 3), RE
COMMON /FLOW/ NVRC, IVRC( 50), VORC( 50), NPSC, IPSC( 50),
1 PSC( 50), COF(100)
C
C ----- To construct F ---
C
CALL MQV12 ( AT, VOR, BE, CE, FVC, UN, VN, U, V, LBC,
1 NDE, NLB, NCN, NNP, NEL, LUV )
C
CALL FVCTN ( D, VRO, FVC, NCC, KCC, NNP, NNZ, NZ2 )
C
C ----- To solve a simultaneous equation ---
C
DO 100 I = 1, NNP
C
IF ( IBV( I).EQ.9999 ) GO TO 100
VOR( I) = VRO( I)
C
100 CONTINUE
C
CALL DSLVE ( FVC, VOR, P2, NCC, KCC, IBV, NNP, NNZ, NZ2,
1 OMG( 2), ERR( 2), IVR )
C
RETURN
END
C *********************************************************************
C
SUBROUTINE LPOV ( PSI, VOR, U, V, NNP, NEL, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
DIMENSION PSI( NNP), VOR( NNP), U( NEL), V( NEL)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
INTEGER*2 CON
C
IF ( MOD( LOP,NPR) ) 120 ,110 ,120
C
110 NPR = LOP + 1
C
GO TO 100
C
120 NPR = LOP
C
100 CONTINUE
C
CALL PRSB ( PSI, VOR, U, V, NNP, NEL, CON )
C
WRITE(6,2000)
2000 FORMAT (/,' *** LOP limit over ***')
WRITE(8,2000)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPST ( PSI, VOR, U, V, NNP, NEL, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 CON
C
DIMENSION PSI( NNP), VOR( NNP), U( NEL), V( NEL)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
IEND = 1
C
IF ( MOD( LOP,NPR) ) 120, 110, 120
C
110 NPR = LOP + 1
C
GO TO 100
C
120 NPR = LOP
C
100 CONTINUE
C
CALL PRSB ( PSI, VOR, U, V, NNP, NEL, CON )
C
IF ( IEND.EQ.1 ) WRITE(6,2000) LOP
IF ( IEND.EQ.1 ) WRITE(8,2000) LOP
C
2000 FORMAT (/' *** Stationary condition reached. *** LOP =',I6)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MODD ( D1, P1, P2, AT, BE, CE, AKX, AKY,
1 AX2, AY2, 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 AKX( NEL), AKY( NEL), AX2( NEL), AY2( NEL)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
C ----- To modify D1 ---
C
IF ( NDM .EQ. 2) THEN
C
CALL DMAT ( AT, NDE, BE, CE, AKX, AKY, AX2, AY2, D1, NCW, NEL,
1 NNZ, 2 )
C
ELSE
CALL DMAT ( AT, NDE, BE, CE, AKX, AKY, AX2, AY2, D1, NCW, NEL,
1 NNZ, 3 )
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
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
C
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
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MQV12 ( AT, C, BE, CE, F, UN, VN, U, V,
1 LBC, NDE, NLB, NCN, NNP, NEL, LUV )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE( NEL,3), NLB( NNP,10), NCN( NNP), LBC( LUV)
C
DIMENSION F ( NNP), U( NEL), V ( NEL), UN( NNP), VN ( NNP),
1 AT( NEL), C( NNP), BE( NEL,3), CE( NEL,3)
C
DO 100 I = 1, NNP
C
UNA = 0.D0
VNA = 0.D0
UNB = 0.D0
VNB = 0.D0

K = NCN ( I)
C
DO 150 N = 1, K
C
UNA = UNA + AT( NLB( I,N) )* U( NLB( I,N) )
VNA = VNA + AT( NLB( I,N) )* V( NLB( I,N) )
C
UNB = UNB + AT( NLB( I,N) )
VNB = VNB + AT( NLB( I,N) )
C
150 CONTINUE
C
UN( I)= UNA / UNB
VN( I)= VNA / VNB
C
100 CONTINUE
C
C
DO 200 J = 1, LUV
C
UN( LBC( J)) = 0.D0
VN( LBC( J)) = 0.D0
C
200 CONTINUE
C
C ----- Attention --- 113/196 ---
C
DO 300 K = 107, 112
C
UN( K) = -1.D0
C
300 CONTINUE
C
CALL CLEAR8 ( F, 0.D0, NNP )
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

CF3= ( B1*U1 + C1*V1 + B2*U2 + C2*V2 + B3*U3 + C3*V3 )/ 24.D0
C
F( N1) = F( N1) + CF1* ( 2.D0*U1+U2+U3 ) + CF2* ( 2.D0*V1+V2+V3 )
C
F( N2) = F( N2) + CF1* ( U1+2.D0*U2+U3 ) + CF2* ( V1+2.D0*V2+V3 )
C
F( N3) = F( N3) + CF1* ( U1+U2+2.D0*U3 ) + CF2* ( V1+V2+2.D0*V3 )
C
C ----- Conservation form considered ---
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
RETURN
END
C **********************************************************************
C
SUBROUTINE NBCN ( SH, DSH, PSI, VOR, NED, LBC, NDB, NNP, LUV )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NED, LBC, NDB
C
DIMENSION SH ( LUV), DSH( LUV), PSI( NNP), VOR( NNP),
1 NED( LUV), LBC( LUV), NDB( LUV,10)
C
C ----- Vorticity on a wall ---
C
IF ( LUV.LE.0 ) RETURN
C
DO 100 I = 1, LUV
C
ND = LBC( I)
NE = NED( I)
C
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.5D0
C
C ----- Attention ( 113 / 196 ) ---
C
IF ( ND.GE. 107.AND.ND.LE. 112 ) VOR( ND) = VOR( ND) - SH( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMAT ( AT, NDE, P, NCC, NCW, KCC, NEL, NNP, NNZ,
1 NZ2, NC2 )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 NDE, KCC, NC2
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
NZWK = 0
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
DO 150 I = 1, 3
C
DO 150 J = 1, 3
C
II = N ( I)
JJ = N ( J)
C
AV = ATI* A1( I,J)
C
CALL DCNL ( AV, II, JJ, P, NCC, NZWK, NNZ )
C
150 CONTINUE
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
CALL MOVEI ( NCC, 1, NCW, NNZ )
C
CALL DCNIJ ( KCC, NC2, NPWK, NZWK, NNP, NZ2 )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRSB ( PSI, VOR, U, V, NNP, NEL, CON )
C
IMPLICIT REAL*8 ( A-H, O-Z)
C
INTEGER*2 CON
C
DIMENSION PSI( NNP), VOR( NNP), U( NEL), V( NEL)
C
COMMON /TIMC/ DT, LPLM, NPR, LOP, ALP
C
ISW = MOD ( LOP,NPR )
C
IF ( ISW.NE.0) RETURN
C
WRITE(8,2000) LOP
2000 FORMAT(64X,' Loop = ',I4)
C
CALL CPRT ( PSI ,'Stream Function ', NNP )
C
CALL CPRT ( VOR ,' Vorticity ',NNP )
C
WRITE(8,2100) ( I, U( I), V( I), I = 1, NEL )
C
2100 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 KOS, WORK1, WORK2, DAT1, DAT2
C
C ----- Attention ; NEL*3 --- DATA 113 / 196 ---
C
DIMENSION DAT1( 588), DAT2( 588)
C
KOS = NEL* 3
C
REWIND 1
C
READ(1,1000) ( DAT1( I), DAT2( I), I = 1, KOS )
1000 FORMAT(2I5)
C
C --------------------------------------
C
DO 100 I = 1, KOS-1
C
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
C
IF ( K.GT.NT ) WRITE(6,2000)
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
A( N) = AV
L( N) = M
C
2000 FORMAT(' * Error K.GT.NT *')
C
RETURN
END
C ******************************* DATA ********************************
ARE NNP= 113 NEL= 196 LUV= 28
TIM DLT=0.071 LPL= 200 STP= 0
JTN 1 1.000000 0.0
JTN 2 0.857143 0.0
JTN 3 0.714286 0.0
JTN 4 0.571429 0.0
JTN 5 0.428571 0.0
JTN 6 0.285714 0.0
JTN 7 0.142857 0.0
JTN 8 0.000000 0.0
JTN 9 0.928571 0.071429
JTN 10 0.785714 0.071429
JTN 11 0.642857 0.071429
JTN 12 0.500000 0.071429
JTN 13 0.357143 0.071429
JTN 14 0.214286 0.071429
JTN 15 0.071429 0.071429
JTN 16 1.000000 0.142857
JTN 17 0.857143 0.142857
JTN 18 0.714286 0.142857
JTN 19 0.571429 0.142857
JTN 20 0.428571 0.142857
JTN 21 0.285714 0.142857
JTN 22 0.142857 0.142857
JTN 23 0.000000 0.142857
JTN 24 0.928571 0.214286
JTN 25 0.785714 0.214286
JTN 26 0.642857 0.214286
JTN 27 0.500000 0.214286
JTN 28 0.357143 0.214286
JTN 29 0.214286 0.214286
JTN 30 0.071429 0.214286
JTN 31 1.000000 0.285714
JTN 32 0.857143 0.285714
JTN 33 0.714286 0.285714
JTN 34 0.571429 0.285714
JTN 35 0.428571 0.285714
JTN 36 0.285714 0.285714
JTN 37 0.142857 0.285714
JTN 38 0.000000 0.285714
JTN 39 0.928571 0.357143
JTN 40 0.785714 0.357143
JTN 41 0.642857 0.357143
JTN 42 0.500000 0.357143
JTN 43 0.357143 0.357143
JTN 44 0.214286 0.357143
JTN 45 0.071429 0.357143
JTN 46 1.000000 0.428571
JTN 47 0.857143 0.428571
JTN 48 0.714286 0.428571
JTN 49 0.571429 0.428571
JTN 50 0.428571 0.428571
JTN 51 0.285714 0.428571
JTN 52 0.142857 0.428571
JTN 53 0.000000 0.428571
JTN 54 0.928571 0.500000
JTN 55 0.785714 0.500000
JTN 56 0.642857 0.500000
JTN 57 0.500000 0.500000
JTN 58 0.357143 0.500000
JTN 59 0.214286 0.500000
JTN 60 0.071429 0.500000
JTN 61 1.000000 0.571429
JTN 62 0.857143 0.571429
JTN 63 0.714286 0.571429
JTN 64 0.571429 0.571429
JTN 65 0.428571 0.571429
JTN 66 0.285714 0.571429
JTN 67 0.142857 0.571429
JTN 68 0.000000 0.571429
JTN 69 0.928571 0.642857
JTN 70 0.785714 0.642857
JTN 71 0.642857 0.642857
JTN 72 0.500000 0.642857
JTN 73 0.357143 0.642857
JTN 74 0.214286 0.642857
JTN 75 0.071429 0.642857
JTN 76 1.000000 0.714286
JTN 77 0.857143 0.714286
JTN 78 0.714286 0.714286
JTN 79 0.571429 0.714286
JTN 80 0.428571 0.714286
JTN 81 0.285714 0.714286
JTN 82 0.142857 0.714286
JTN 83 0.000000 0.714286
JTN 84 0.928571 0.785714
JTN 85 0.785714 0.785714
JTN 86 0.642857 0.785714
JTN 87 0.500000 0.785714
JTN 88 0.357143 0.785714
JTN 89 0.214286 0.785714
JTN 90 0.071429 0.785714
JTN 91 1.000000 0.857143
JTN 92 0.857143 0.857143
JTN 93 0.714286 0.857143
JTN 94 0.571429 0.857143
JTN 95 0.428571 0.857143
JTN 96 0.285714 0.857143
JTN 97 0.142857 0.857143
JTN 98 0.000000 0.857143
JTN 99 0.928571 0.928571
JTN 100 0.785714 0.928571
JTN 101 0.642857 0.928571
JTN 102 0.500000 0.928571
JTN 103 0.357143 0.928571
JTN 104 0.214286 0.928571
JTN 105 0.071429 0.928571
JTN 106 1.000000 1.000000
JTN 107 0.857143 1.000000
JTN 108 0.714286 1.000000
JTN 109 0.571429 1.000000
JTN 110 0.428571 1.000000
JTN 111 0.285714 1.000000
JTN 112 0.142857 1.000000
JTN 113 0.000000 1.000000
ELM 1 1 16 9 1
ELM 5 2 17 10 1
ELM 9 3 18 11 1
ELM 13 4 19 12 1
ELM 17 5 20 13 1
ELM 21 6 21 14 1
ELM 25 7 22 15 1
ELM 2 1 9 2 1
ELM 6 2 10 3 1
ELM 10 3 11 4 1
ELM 14 4 12 5 1
ELM 18 5 13 6 1
ELM 22 6 14 7 1
ELM 26 7 15 8 1
ELM 3 16 17 9 1
ELM 7 17 18 10 1
ELM 11 18 19 11 1
ELM 15 19 20 12 1
ELM 19 20 21 13 1
ELM 23 21 22 14 1
ELM 27 22 23 15 1
ELM 4 17 2 9 1
ELM 8 18 3 10 1
ELM 12 19 4 11 1
ELM 16 20 5 12 1
ELM 20 21 6 13 1
ELM 24 22 7 14 1
ELM 28 23 8 15 1
ELM 29 16 31 24 1
ELM 33 17 32 25 1
ELM 37 18 33 26 1
ELM 41 19 34 27 1
ELM 45 20 35 28 1
ELM 49 21 36 29 1
ELM 53 22 37 30 1
ELM 30 16 24 17 1
ELM 34 17 25 18 1
ELM 38 18 26 19 1
ELM 42 19 27 20 1
ELM 46 20 28 21 1
ELM 50 21 29 22 1
ELM 54 22 30 23 1
ELM 31 31 32 24 1
ELM 35 32 33 25 1
ELM 39 33 34 26 1
ELM 43 34 35 27 1
ELM 47 35 36 28 1
ELM 51 36 37 29 1
ELM 55 37 38 30 1
ELM 32 32 17 24 1
ELM 36 33 18 25 1
ELM 40 34 19 26 1
ELM 44 35 20 27 1
ELM 48 36 21 28 1
ELM 52 37 22 29 1
ELM 56 38 23 30 1
ELM 57 31 46 39 1
ELM 61 32 47 40 1
ELM 65 33 48 41 1
ELM 69 34 49 42 1
ELM 73 35 50 43 1
ELM 77 36 51 44 1
ELM 81 37 52 45 1
ELM 58 31 39 32 1
ELM 62 32 40 33 1
ELM 66 33 41 34 1
ELM 70 34 42 35 1
ELM 74 35 43 36 1
ELM 78 36 44 37 1
ELM 82 37 45 38 1
ELM 59 46 47 39 1
ELM 63 47 48 40 1
ELM 67 48 49 41 1
ELM 71 49 50 42 1
ELM 75 50 51 43 1
ELM 79 51 52 44 1
ELM 83 52 53 45 1
ELM 60 47 32 39 1
ELM 64 48 33 40 1
ELM 68 49 34 41 1
ELM 72 50 35 42 1
ELM 76 51 36 43 1
ELM 80 52 37 44 1
ELM 84 53 38 45 1
ELM 85 46 61 54 1
ELM 89 47 62 55 1
ELM 93 48 63 56 1
ELM 97 49 64 57 1
ELM 101 50 65 58 1
ELM 105 51 66 59 1
ELM 109 52 67 60 1
ELM 86 46 54 47 1
ELM 90 47 55 48 1
ELM 94 48 56 49 1
ELM 98 49 57 50 1
ELM 102 50 58 51 1
ELM 106 51 59 52 1
ELM 110 52 60 53 1
ELM 87 61 62 54 1
ELM 91 62 63 55 1
ELM 95 63 64 56 1
ELM 99 64 65 57 1
ELM 103 65 66 58 1
ELM 107 66 67 59 1
ELM 111 67 68 60 1
ELM 88 62 47 54 1
ELM 92 63 48 55 1
ELM 96 64 49 56 1
ELM 100 65 50 57 1
ELM 104 66 51 58 1
ELM 108 67 52 59 1
ELM 112 68 53 60 1
ELM 113 61 76 69 1
ELM 117 62 77 70 1
ELM 121 63 78 71 1
ELM 125 64 79 72 1
ELM 129 65 80 73 1
ELM 133 66 81 74 1
ELM 137 67 82 75 1
ELM 114 61 69 62 1
ELM 118 62 70 63 1
ELM 122 63 71 64 1
ELM 126 64 72 65 1
ELM 130 65 73 66 1
ELM 134 66 74 67 1
ELM 138 67 75 68 1
ELM 115 76 77 69 1
ELM 119 77 78 70 1
ELM 123 78 79 71 1
ELM 127 79 80 72 1
ELM 131 80 81 73 1
ELM 135 81 82 74 1
ELM 139 82 83 75 1
ELM 116 77 62 69 1
ELM 120 78 63 70 1
ELM 124 79 64 71 1
ELM 128 80 65 72 1
ELM 132 81 66 73 1
ELM 136 82 67 74 1
ELM 140 83 68 75 1
ELM 141 76 91 84 1
ELM 145 77 92 85 1
ELM 149 78 93 86 1
ELM 153 79 94 87 1
ELM 157 80 95 88 1
ELM 161 81 96 89 1
ELM 165 82 97 90 1
ELM 142 76 84 77 1
ELM 146 77 85 78 1
ELM 150 78 86 79 1
ELM 154 79 87 80 1
ELM 158 80 88 81 1
ELM 162 81 89 82 1
ELM 166 82 90 83 1
ELM 143 91 92 84 1
ELM 147 92 93 85 1
ELM 151 93 94 86 1
ELM 155 94 95 87 1
ELM 159 95 96 88 1
ELM 163 96 97 89 1
ELM 167 97 98 90 1
ELM 144 92 77 84 1
ELM 148 93 78 85 1
ELM 152 94 79 86 1
ELM 156 95 80 87 1
ELM 160 96 81 88 1
ELM 164 97 82 89 1
ELM 168 98 83 90 1
ELM 169 91 106 99 1
ELM 173 92 107 100 1
ELM 177 93 108 101 1
ELM 181 94 109 102 1
ELM 185 95 110 103 1
ELM 189 96 111 104 1
ELM 193 97 112 105 1
ELM 170 91 99 92 1
ELM 174 92 100 93 1
ELM 178 93 101 94 1
ELM 182 94 102 95 1
ELM 186 95 103 96 1
ELM 190 96 104 97 1
ELM 194 97 105 98 1
ELM 171 106 107 99 1
ELM 175 107 108 100 1
ELM 179 108 109 101 1
ELM 183 109 110 102 1
ELM 187 110 111 103 1
ELM 191 111 112 104 1
ELM 195 112 113 105 1
ELM 172 107 92 99 1
ELM 176 108 93 100 1
ELM 180 109 94 101 1
ELM 184 110 95 102 1
ELM 188 111 96 103 1
ELM 192 112 97 104 1
ELM 196 113 98 105 1
PSI 01 0. 02 0. 03 0.
PSI 04 0. 05 0. 06 0.
PSI 07 0. 08 0.
PSI 16 0. 31 0. 46 0.
PSI 61 0. 76 0. 91 0.
PSI 23 0. 38 0. 53 0.
PSI 68 0. 83 0. 98 0.
PSI 106 0. 107 0. 108 0.
PSI 109 0. 110 0. 111 0.
PSI 112 0. 113 0.
WAL 01 16 31 46 61 76 91 106 107 108
WAL 109 110 111 112 113 98 83 68 53 38
WAL 23 08 07 06 05 04 03 02
C *********************************************************