戻る

C ******************************************************************* C
C C
C F 2 D - F E R C
CC C
C 2D - Fluid Analysis Software by F.E.M. C
C C
C ( Rectangular Linear Element ) C
C C
C ( Version 3.5 ) C
C C
C ( 121 / 100 ) 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, KC1, MTR, IBP, IBV, LBC, NBD, NTB, NCN,
1 NEB
C
C ----- 121 NDEs / 100 elements ---
C
DIMENSION X ( 121), Y ( 121), NDE( 100,4), U ( 100),
1 V ( 100), UO ( 100), VO ( 100), UN ( 121),
2 VN ( 121), ARS( 100), MTR( 100), LBC( 40),
3 P1 ( 541), P2 ( 541), NC1( 961), NC2( 2,961),
4 NCW( 541), FVC( 121), DP ( 541), D ( 541),
5 RKX( 100), RKY( 100)
C
DIMENSION PSI( 121), PSO( 121), IBP( 121), VOR( 121),
1 VRO( 121), IBV( 121), NEB( 121,10), NBD( 40)
C
DIMENSION DSH( 40), SH ( 40), ITB( 121), IEB( 100),
1 NTB( 40), NCN( 121), KC1( 121), HXX( 12),
2 HYY( 12), FF ( 121), ARW( 121)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF( 100)
C
OPEN ( 7, FILE='F2D_FER_INP10.DAT')
OPEN (10, FILE='F2D_FER_PLT.DAT' )
C
READ(7,1000) NNP, NEL, LUV
1000 FORMAT(10X,I5,5X,I5,5X,I5)
C
CALL TTLE
C
CALL INPT ( X, Y, U, V, DSH, SH, NDE, ITB, IEB,
1 NTB, NBD, MTR, LBC, NNZ, NZ2, NR2, HXX, HYY,
2 HKE, HN1, HN2 )
C
T = 0.D0
LOP = 0
C
CALL INSB ( X, Y, PSI, PSO, VOR, VRO, U, V,
1 UN, VN, UO, VO, MTR, RKX, RKY, NC1,
2 NDE, IBP, IBV, LBC, NCN, NEB, NNZ, NR2,
3 HXX, HYY, FF, ARW )
C
C ----- P1 - Matrix ---
C
CALL PMAT ( P1, NC1, NCW, X, Y, NDE, KC1, NC2, NNZ, NZ2,
1 ARS )
C
C ----- DP & D - Matrix ---
C
CALL DMAT ( X, Y, RKX, RKY, DP, NDE, NCW, NNZ, 1 )
C
CALL DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, 2 )
C
C ----- VOR on a wall ---
C
CALL NBCN ( SH, DSH, PSI, VRO, NBD, LBC, NTB, NNP, LUV )
C
9999 LOP = LOP + 1
T = T + DT
C
C ----- Stream function ---
C
CALL LPPS ( U, V, UN, VN, PSI, PSO, IBP, VRO,
1 DP, RKX, RKY, NCN, FVC, NC2, KC1, NDE,
2 MTR, X, Y, IBV, NNZ, NZ2, NR2, HXX,
3 HYY, FF, ARW )
C
C ----- VOR on walls ---
C
CALL NBCN ( SH, DSH, PSI, VOR, NBD, LBC, NTB, NNP, LUV )
C
C ----- To modify D & D2 ---
C
CALL MODD ( D, P1, P2, X, Y, RKX, RKY, NCW, NDE, NNZ, 2)
C
CALL LPVR ( UN, VN, VOR, VRO, IBV, D, P2, FVC, NC2,
1 KC1, NDE, X, Y, NNZ, NZ2, DMX, DN )
C
C ----- Stationarity check ---
C
CALL CHCK ( PSI, PSO, VOR, VRO, U, V, UN, VN, ICK,
1 UO, VO, ARS, DFR, *999 )
C
IF ( LOP.LT.NLM) GO TO 9999
C
999 CONTINUE
C
CALL PLTD ( PSI, VOR, UN, VN, X, Y, NDE )
C
CLOSE ( 7)
CLOSE (10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE PLTD ( PSI, VOR, UN, VN, X, Y, NDE )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /ITL/ IPSI, ITEP, IVOR
COMMON /CS2/ DT, NLM, LOP, ALP
C
COMMON /VAR/ VARV( 500), VRUV( 500), VARP( 500)
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF(100)
C
DIMENSION PSI( NNP), VOR( NNP), UN( NNP), VN( NNP),
1 X ( NNP), Y ( NNP), NDE( NEL,4)
C
WRITE(10,2000) LOP, DT, RE, NNP, NEL
2000 FORMAT(I4, F6.3, F8.1, 2I4)
WRITE(10,2100) ( X( I), Y( I), I = 1, NNP)
2100 FORMAT(10F8.3)
WRITE(10,2200) ( ( NDE( I,J), I = 1, NEL), J = 1, 4 )
2200 FORMAT(20I4)
C
WRITE(10,2100) ( VRUV( I), VARP( I), I = 1, LOP)
WRITE(10,2100) ( VARV( I), I = 1, LOP)
WRITE(10,2100) ( UN ( I), VN( I), I = 1, NNP)
WRITE(10,2100) ( PSI ( I), VOR( I), I = 1, NNP)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ARST ( NDE, NNZ, NZ2 )
C
C ----- To set NNZ , NZ2 etc. ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION NDE( NEL,4), INDE2( 600), M( 6)
C
COMMON /NNL/ NNP, NEL, LUV
C
C ----- Attention --- INDE2 : NEL*6 ---
C
DIMENSION KZ3( 2,6)
C
DATA KZ3 / 1,2, 1,3, 1,4, 2,3, 2,4, 3,4 /
C
IF ( NNP.EQ. 2601) L= 10100
IF ( NNP.EQ. 2601) GO TO 110
C
IF ( NNP.EQ.10201) L= 40200
IF ( NNP.EQ.10201) GO TO 110
C
IF ( NNP.EQ. 6561) L = 25760
IF ( NNP.EQ. 6561) GO TO 110
C
L = 0
C ------------------------------------------------
DO 100 I = 1, NEL
C
C --------------------------------------
DO 150 J = 1, 6
C
N1 = NDE( I,KZ3( 1,J))
N2 = NDE( I,KZ3( 2,J))
M1 = MIN0( N1,N2)
M2 = MAX0( N1,N2)
M( J) = M1* 10300 + M2
C
150 CONTINUE
C --------------------------------------
C
C --------------------------------------
DO 300 J = 1, 6
C
IL = 1
20 IF ( IL.GT.L) GO TO 10
IF ( M( J).EQ.INDE2( IL)) GO TO 300
IL = IL + 1
GO TO 20
C
10 L = L + 1
INDE2( L) = M( J)
C
300 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
110 CONTINUE
NNZ = NNP + L
NZ2 = NNP + L* 2
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BC10 ( DSH, SH, NTB, NBD, Y )
C
C ----- Attention : Boundary condition for Cavity problem - 121/100 -
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NTB, NBD
C
DIMENSION DSH( LUV), SH( LUV), NTB( LUV), NBD( LUV), Y (NNP)
C
C ----- Attention NTB ( LUV) ---
C
COMMON /NNL/ NNP, NEL, LUV
C
C ----- Attention For 100 / 121 case Only ---
C
DMS = Y( 116) - Y( 105)
SHV = 3.D0/ DMS
DSW = 3.D0/( DMS** 2 )
C
*** WRITE(6,*) ' * DMS, SH, DSH =', DMS, SHV, DSW
C
C ------------------------------------------------
DO 100 I = 1, LUV
C
* SH ( I) = 30.D0 DSH( I) = 300.D0
C
SH ( I) = SHV
DSH( I) = DSW
C ===================
C
C ----- NBD = 0 == > VORTICITY = 0 ---
C
NBD( I)= 1
C
100 CONTINUE
C ------------------------------------------------
C
NBD( 1) = 0
NBD(11) = 0
NBD(21) = 0
NBD(31) = 0
C
NTB( 1) = 13
NTB( 2) = 13
NTB( 3) = 24
NTB( 4) = 35
NTB( 5) = 46
NTB( 6) = 57
NTB( 7) = 68
NTB( 8) = 79
NTB( 9) = 90
NTB(10) =101
C
NTB(11) =101
NTB(12) =101
NTB(13) =102
NTB(14) =103
NTB(15) =104
NTB(16) =105
NTB(17) =106
NTB(18) =107
NTB(19) =108
NTB(20) =109
C
NTB(21) =109
NTB(22) =109
NTB(23) = 98
NTB(24) = 87
NTB(25) = 76
NTB(26) = 65
NTB(27) = 54
NTB(28) = 43
NTB(29) = 32
NTB(30) = 21
C
NTB(31) = 21
NTB(32) = 21
NTB(33) = 20
NTB(34) = 19
NTB(35) = 18
NTB(36) = 17
NTB(37) = 16
NTB(38) = 15
NTB(39) = 14
NTB(40) = 13
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( PSI, PSO, VOR, VRO, U, V, UN, VN, ICK,
1 UO, VO, ARS, DFR, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PSI( NNP), PSO( NNP), VOR( NNP), VRO( NNP),
1 U ( NEL), V ( NEL), UN ( NNP), VN ( NNP),
2 UO ( NEL), VO ( NEL), ARS( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /VAR/ VARV( 500), VRUV( 500), VARP( 500)
C
EPS = 0.01D0 ! (%)
C -------------------------
ICK = 0
VRMX = -100.D0
C
DO 100 I = 1, NNP
C
VRMX = DMAX1( VRMX, DABS( PSI( I)))
100 CONTINUE
VMX = VRMX* 0.01D0 ! 1(%) Up
C ------------------------------------------------
C
DMX = -100.D0
DO 200 I = 1, NNP
C
IF ( DABS( PSI( I)).LE.VMX) GO TO 200
DMX = DMAX1( DMX, DABS( (PSI( I) - PSO( I))/ PSI( I)))
200 CONTINUE
VARP( LOP) = DMX* 100.D0
C ------------------------------------------------
C
VRMX = -100.D0
DO 300 I = 1, NNP
C
VRMX = DMAX1( VRMX, DABS( VOR( I)))
300 CONTINUE
VMX = VRMX* 0.01D0 ! 1(%) Up
C ------------------------------------------------
C
DMX = -100.D0
DO 400 I = 1, NNP
C
IF ( DABS( VOR( I)).LE.VMX) GO TO 400
DMX = DMAX1( DMX, DABS( (VOR( I) - VRO( I))/ VOR( I)))
C
400 CONTINUE
VARV( LOP) = DMX* 100.D0
C ------------------------------------------------
C
VMIN = 0.01D0 ! 1(%) Up
DFR = -1000.D0
DO 500 I = 1, NEL
C
IF ( DABS( U( I)).LE.VMIN ) GO TO 500
DFR = DMAX1( DFR, DABS( ( U( I) - UO( I))/ U( I)))
C
500 CONTINUE
VRUV( LOP) = DFR* 100.D0
C ------------------------------------------------
C
NPR = 10
C =============
C
IF ( MOD( LOP,NPR).EQ.0)
1WRITE(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
WRITE(6,*) ' '
WRITE(6,2100) LOP
2100 FORMAT(' ** Stationary condition reached at Loop =', I4)
WRITE(6,2000) LOP, VARP( LOP), VARV( LOP), VRUV( LOP)
RETURN 1
END IF
C
C -------------------------------
IF ( LOP.GE.NLM) THEN
WRITE(6,2200)
2200 FORMAT (/,' Loop Limit Over')
RETURN 1
END IF
C
C ------------------------------------------------
DO 600 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
600 CONTINUE
C ------------------------------------------------
C
CALL MOVR ( PSI, 1.D0, PSO, NNP)
C
CALL MOVR ( VOR, 1.D0, VRO, NNP)
C
CALL MOVR ( U, 1.D0, UO, NEL)
C
CALL MOVR ( V, 1.D0, VO, NEL)
C
RETURN
C
900 WRITE(6,2300)
2300 FORMAT(/,' *** U. or V .GE.2 * STOP ***')
STOP
END
C **********************************************************************
C
SUBROUTINE VXG ( 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 VX2 ( 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
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 NDE.LE.32768 ---
C
IF ( I.GT.J) RETURN
C
M = I* 10300 + J
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 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
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
C
500 N = N2
GO TO 120
C
600 K = K + 1
IF ( K.GT.NT) WRITE(6,2000)
2000 FORMAT(' ** Error ** ''NT'' IN ''SOR'' IS TOO SMALL **')
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)
C
STOP
END
C **********************************************************************
C
SUBROUTINE DCIJ ( IA, JA, NI, NJ, NIX, NJX)
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IA( NIX)
CHARACTER*4 ANIX, ANJX
C
DIMENSION JA( 2,NJX)
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 -------------------------------------------
DO 100 I = 1, NJ
C
IF ( NI.EQ.JA( 1,I)) GO TO 101
C
C ----- To set IA( NI) ---
C
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 ----- To set 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 550 L = I, NI
C
IA( L) = IA( L) + 1
550 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 ------------------------------------------------
RETURN
C
601 WRITE(6,2000) ANIX, NIX, ANIX, NI
STOP
C
602 WRITE(6,2000) ANJX, NJX, ANJX, NJ
2000 FORMAT(' ** Error ( DCIJ) ',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 R '
WRITE(6,*) ' '
WRITE(6,*) ' ( Version : 3.5 ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Modified Galerkin Method - Rectangular Linear ) '
WRITE(6,*) ' '
WRITE(6,*) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION X( NNP), Y ( NNP), RKX( NNP), RKY( NNP),
1 D( NNZ), NDE( NEL,4), NCW( NNZ), A1 ( 4,4),
2 N( 4)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF( 100)
C
CALL VXG ( D, 0.D0, NNZ )
C
NZK = NNZ
C
C --------------------------------------
DO 100 IEL = 1, NEL
C
N( 1) = NDE( IEL,1)
N( 2) = NDE( IEL,2 )
N( 3) = NDE( IEL,3 )
N( 4) = NDE( IEL,4)
A = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))/ 2.D0
B = DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))/ 2.D0
C
IF ( NDM - 2) 2000, 2100, 2200
C
2000 CONTINUE
W1 = A / B
W2 = B / A
GO TO 2500
C
C ----- Attention : Dynamic Viscosity = Thermal Conductivity ---
C
2100 CONTINUE
W1 = RKX( IEL)* A/ B
W2 = RKY( IEL)* B/ A
GOTO 2500
C
2200 CONTINUE
C
TKAK = RKX( IEL) - COF( 1) + COF( 1)
W1 = TKAK* A / B
W2 = TKAK* B / A
C
C ----- W1 = A / B, W2 = B / A ---
C
2500 CONTINUE
A1( 1,1) = W1/ 3.D0 + W2 / 3.D0
A1( 1,2) = - W1/ 3.D0 + W2/ 6.D0
A1( 1,3) = - W1/ 6.D0 - W2/ 6.D0
A1( 1,4) = W1/ 6.D0 - W2/ 3.D0
C
A1( 2,1) = - W1/ 3.D0 + W2/ 6.D0
A1( 2,2) = W1/ 3.D0 + W2/ 3.D0
A1( 2,3) = W1/ 6.D0 - W2/ 3.D0
A1( 2,4) = - W1/ 6.D0 - W2/ 6.D0
C
A1( 3,1) = - W1/ 6.D0 - W2/ 6.D0
A1( 3,2) = W1/ 6.D0 - W2/ 3.D0
A1( 3,3) = W1/ 3.D0 + W2/ 3.D0
A1( 3,4) = - W1/ 3.D0 + W2/ 6.D0
C
A1( 4,1) = W1/ 6.D0 - W2/ 3.D0
A1( 4,2) = - W1/ 6.D0 - W2/ 6.D0
A1( 4,3) = - W1/ 3.D0 + W2/ 6.D0
A1( 4,4) = W1/ 3.D0 + W2/ 3.D0
C
C ----------------------------
DO 150 I = 1, 4
DO 150 J = 1, 4
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 DSLV ( FVC, C, TSM, NCC, KCC, NBC, NNZ, NZ2,
1 OMG, ERR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, NBC
C
DIMENSION FVC( NNP), C( NNP), TSM( NNZ), NCC( 2,NZ2), KCC( NNP),
1 NBC( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
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
LMPR = 3* NNP
KTAL = 2* NNP
NTAL = 0
C
900 NTAL = NTAL + 1
ERMX = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 200
IF ( I.EQ.1) GO TO 210
C
MS = KCC( I-1) + 1
C
GOTO 220
C
210 MS = 1
220 ML = KCC( I)
C
IF ( MS.GT.ML) GO TO 1000
C
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
DK = DK - TSM( M1)* C( M2 )
C
300 CONTINUE
C ----------------------------
C
DC = - C( I) + DSM* DK
C( I) = C( I) + OMG* DC
C
TEMP = DABS( C( I))
IF ( TEMP.LE.EPS) GO TO 200
TEMP = DABS( DC/ C( I))
IF ( TEMP.LE.ERMX) GO TO 200
ERMX = TEMP
C
200 CONTINUE
C ------------------------------------------------
C
IF ( ERMX.LE.ERR ) RETURN
IF ( NTAL.GE.KTAL) RETURN
IF ( NTAL.GE.LMPR) GO TO 900
C
GOTO 900
C
1000 WRITE(6,2000) I, MS, ML
2000 FORMAT ('1',' ** Error in "DSLV" ** ',3I6)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVCTN ( D, CO, F, NC2, KCC, NNZ, NZ2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC
C
DIMENSION D( NNZ), CO( NNP), F( NNP), NC2( 2,NZ2), KCC( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
C --------------------------------------
DO 100 I = 1, NNP
C
JS = 1
JL = KCC( I)
IF ( I.GT.1) JS = KCC( I-1) + 1
FW = - F( I)
C ----------------------------
DO 150 J = JS, JL
C
KPR = NC2( 1,J)
KCL = NC2( 2,J)
FW = FW + D( KPR)* CO( KCL)
C
150 CONTINUE
C ----------------------------
F( I) = FW
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVOR ( VOR, NDE, FVC, X, Y )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION VOR( NNP), FVC( NNP), X( NNP), Y( NNP), NDE( NEL,4)
C
COMMON /NNL/ NNP, NEL, LUV
C
CALL VXG ( FVC, 0.D0, NNP)
C
DO 100 NE = 1, NEL
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
S = A* B / 9.D0
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
C
FVC( N1) = FVC( N1)
C
1 + ( 4.D0*VOR( N1) + 2.D0*VOR( N2) + VOR( N3) + 2.D0*VOR( N4))*S
C
FVC( N2) = FVC( N2)
C
1 + ( 2.D0*VOR( N1) + 4.D0*VOR( N2) + 2.D0*VOR( N3) + VOR( N4))*S
C
FVC( N3) = FVC( N3)
C
1 + ( VOR( N1) + 2.D0*VOR( N2) + 4.D0*VOR( N3) + 2.D0*VOR( N4))*S
C
FVC( N4) = FVC( N4)
C
1 + ( 2.D0*VOR( N1) + VOR( N2) + 2.D0*VOR( N3) + 4.D0*VOR( N4))*S
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INSB ( X, Y, PSI, PSO, VOR, VRO, U, V,
1 UN, VN, UO, VO, MTR, RKX, RKY, NC1,
2 NDE, IBP, IBV, LBC, NCN, NEB, NNZ, NR2,
3 HXX, HYY, FF, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR, IBP, IBV, LBC, NCN, NEB
C
DIMENSION X ( NNP), Y ( NNP), PSI( NNP), PSO( NNP),
1 VOR( NNP), VRO( NNP), U ( NEL), V ( NEL),
2 UO ( NEL), VO ( NEL), UN ( NNP), VN ( NNP),
3 RKX( NNP), RKY( NNP), NC1( NNZ), NEB( NNP,10),
4 HXX( NR2), HYY( NR2 ), FF ( NNP)
C
DIMENSION NDE( NEL,4), IBP( NNP), IBV( NNP), LBC( LUV),
1 NCN( NNP), MTR( NEL), ARW( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF(100)
C
CALL VXG ( ARW, 0.D0, NNP)
CALL VXG ( UO, 0.D0, NEL )
CALL VXG ( VO, 0.D0, NEL )
C
CALL VX2 ( IBP, 1, NNP)
CALL VX2 ( IBV, 1, NNP)
C
CALL VXG ( PSI, 0.D0, NNP)
CALL VXG ( PSO, 0.D0, NNP)
CALL VXG ( VOR, 0.D0, NNP)
CALL VXG ( VRO, 0.D0, NNP)
C
IF ( NPC.LE.0) GO TO 110
C
C --------------------------------------
DO 100 I = 1, NPC
C
IBP( IPC( I)) = 9999
PSI( IPC( I)) = PSC( I)
C
100 CONTINUE
C --------------------------------------
C
110 IF ( NVC.LE.0) GO TO 121
C
C --------------------------------------
DO 200 I = 1, NVC
C
IBV( IVR( I)) = 9999
VOR( IVR( I)) = VRC( I)
C
200 CONTINUE
C --------------------------------------
121 GO TO 130
C
130 IF ( LUV.LE.0) GO TO 150
C
C --------------------------------------
DO 400 I = 1, LUV
C
IBV( LBC( I)) = 9999
400 CONTINUE
C --------------------------------------
150 CONTINUE
C
C ----- No. of meshes at each NDE ---
C
C --------------------------------------
DO 500 I = 1, NNP
C
NCN ( I) = 0
500 CONTINUE
C --------------------------------------
DO 600 I = 1, NEL
DO 600 J = 1, 4
C
NE = NDE( I,J)
NCN ( NE ) = NCN( NE ) + 1
IF ( NCN ( NE).GT.10) WRITE(6,2000) NCN ( NE)
2000 FORMAT(' *** NO. OF ELEMENTS AT NDE=',I5,' GT 10 *** STOP' )
IF ( NCN ( NE).GT.10) STOP
NEB( NE,NCN(NE)) = I
C
600 CONTINUE
C --------------------------------------
C
CALL VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV, NCN, ARW )
C
IF ( NNP.EQ.63 ) CALL VXG ( U, 0.5D0, NEL)
C
C --------------------------------------
DO 700 I = 1, NEL
C
RKX( I) = COF( MTR( I))
RKY( I) = COF( MTR( I))
C
700 CONTINUE
C -------------------------------------
C
CALL MOVR ( PSI, 1.D0, PSO, NNP)
C
CALL MOVR ( VOR, 1.D0, VRO, NNP)
C
CALL CKAK ( RKX, RKY, U, V, MTR, NR2, HXX, HYY, UN, VN, NDE, FF )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( X, Y, U, V, DSH, SH, NDE, ITB,
1 IEB, NTB, NBD, MTR, LBC, NNZ, NZ2, NR2,
2 HXX, HYY, HKE, HN1, HN2 )
C
C ----- Input DATA ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, ITB, IEB, MTR, LBC, NTB, NBD
C
CHARACTER*3 HEAD
C
DIMENSION X ( NNP), Y ( NNP), U ( NEL), V ( NEL),
1 DSH( LUV), SH ( LUV), NDE( NEL,4), ITB( NNP),
2 IEB( NEL), NTB( LUV), NBD( LUV), MTR( NEL),
3 LBC( LUV), HXX( NR2), HYY( NR2)
C
DIMENSION IWK( 10), WKX( 10), IWKX( 10,10), WK( 10,10)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
C
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF(100)
C
NSR = INT( DSQRT( DFLOAT( NEL)))
NR2 = NSR + 2
C
NPX = 0
NEX = 0
NVC = 0
NPC = 0
NLUV = 0
C
CALL VXG ( COF, 1.D0, 100 )
CALL VXG ( U, 0.D0, NEL )
CALL VXG ( V, 0.D0, NEL )
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.001D0
DLT( 2) = 0.01D0
DLT( 3) = 0.01D0
C
C ----- Time scheme parameter ---
C
ALP = 1.D0
C
C ----- TIM ---
C
READ(7,1000) DT, NLM, LOP1, ALPA, DLT1
1000 FORMAT( 10X,F12.0,5X,I5,5X,I5,5X,F10.0,5X,F10.0)
C
IF ( DABS( DLT1).GE.0.0001D0) DLT( 1) = DLT1
C
C ----- MTR ---
C
READ(7,1100) J, COF( 1), RA, ICN
1100 FORMAT( 10X,I5,5X,F10.0,5X,F10.0,4X,I1)
C
C ----- Attention Heat & Fluid Only ---
C
WRITE(6,*) ' '
WRITE(6,2000) NNP, NEL, LUV
2000 FORMAT(' ( NNP =',I5,' NEL =',I5,' LUV =',I4,' )')
WRITE(6,*) ' '
C
ICN = 0
C
C ----- Attention ----
C
WRITE(6,*) ' ** Re = ?'
** READ(5,*) RE
C
RE = 100.D0
C =================
C
COF( 1) = 1.D0/ RE
C =======================
C
WRITE(6,*) ' '
WRITE(6,*) ' ** DT = ? * LOP = ?'
** READ( 5,*) DT, NLM
DT = 0.1D0
NLM = 150
C
WRITE(6,2100) RE, DT, NLM, DLT1
2100 FORMAT(/,' * Re =',F7.1,' DT =',F6.3,' LOP =',I4,' EPS =',
1 F8.5,/ )
C
C ----- JTN ---
C
DO 100 II = 1, NNP
C
READ(7,1200) IWK( 1), WK( 1,1), WK( 2,1)
1200 FORMAT( 10X,I5,5X,2F10.0)
C
NPX = NPX + 1
ITB( NPX) = IWK( 1)
X ( NPX) = WK( 1,1)
Y ( NPX) = WK( 2,1)
C
100 CONTINUE
C
KKK = 1
C =============
C
IF ( KKK.NE.1) GO TO 3344
C
C ----- Mesh with equal length ---
C ========================
C
C --------------------------------------
DO 200 I = 1, 11
C
X( I) = 1.D0 - 0.1D0*(I-1)
Y( I) = 0.D0
C
200 CONTINUE
C --------------------------------------
DO 300 I = 12, 22
C
X( I) = 1.D0 - 0.1D0*(I-12)
Y( I) = 0.1D0
C
300 CONTINUE
C --------------------------------------
DO 400 I = 23, 33
C
X( I) = 1.D0 - 0.1D0* (I-23)
Y( I) = 0.2D0
C
400 CONTINUE
C --------------------------------------
DO 500 I = 34, 44
C
X( I) = 1.D0 - 0.1D0* (I-34)
Y( I) = 0.3D0
C
500 CONTINUE
C --------------------------------------
DO 600 I = 45, 55
C
X( I) = 1.D0 - 0.1D0* (I-45)
Y( I) = 0.4D0
C
600 CONTINUE
C --------------------------------------
DO 700 I = 56, 66
C
X( I) = 1.D0 - 0.1D0* (I-56)
Y( I) = 0.5D0
C
700 CONTINUE
C --------------------------------------
DO 800 I = 67, 77
C
X( I) = 1.D0 - 0.1D0* (I-67)
Y( I) = 0.6D0
C
800 CONTINUE
C --------------------------------------
DO 900 I = 78, 88
C
X( I) = 1.D0 - 0.1D0* (I-78)
Y( I) = 0.7D0
C
900 CONTINUE
C --------------------------------------
DO 920 I = 89, 99
C
X( I) = 1.D0 - 0.1D0* (I-89)
Y( I) = 0.8D0
C
920 CONTINUE
C --------------------------------------
DO 940 I = 100, 110
C
X( I) = 1.D0 - 0.1D0* (I-100)
Y( I) = 0.9D0
C
940 CONTINUE
C --------------------------------------
DO 960 I = 111, 121
C
X( I) = 1.D0 - 0.1D0* (I-111)
Y( I) = 1.D0
C
960 CONTINUE
C --------------------------------------
3344 CONTINUE
C
C ----- ELM ---
C
HKE = ( X( 3)-X( 2))/ ( X( 2)-X( 1))
C
HN1 = ( HKE + 1.D0)* ( HKE + 1.D0)* (HKE + 1.D0)
HN2 = HKE* ( HKE + 1.D0)* ( HKE + 1.D0)
C
C --------------------------------------
DO 980 II = 1, NEL
C
READ(7,1300) ( IWKX( 1,J), J = 1, 6 )
1300 FORMAT(10X,6I5)
IWKX( 1,6) = MAX0 ( IWKX( 1,6),1)
NEX = NEX + 1
IEB( NEX) = IWKX( 1,1)
C
NDE( NEX,1)= IWKX( 1,2)
NDE( NEX,2)= IWKX( 1,3)
NDE( NEX,3)= IWKX( 1,4)
NDE( NEX,4)= IWKX( 1,5)
MTR( NEX) = IWKX( 1,6)
C
980 CONTINUE
C --------------------------------------
C
C ----- PSI ---
C
121 READ(7,1400) HEAD
1400 FORMAT(A3)
C
IF ( HEAD.NE.'PSI' ) GO TO 140
BACKSPACE 7
C
READ(7,1500) ( IWK( I), WKX( I), I = 1, 3)
1500 FORMAT( 10X,3( I5,5X,F10.0))
C
C --------------------------------------
DO 985 I = 1, 3
C
IF ( IWK( I).EQ.0 .AND. DABS(WKX( I)).LE.0.0001D0) GO TO 985
NPC = NPC + 1
IPC( NPC) = IWK( I)
PSC( NPC) = WKX( I)
C
985 CONTINUE
C --------------------------------------
GOTO 121
C
C ----- TEP ---
C
140 BACKSPACE 7
141 READ(7,1400) HEAD
C
IF ( HEAD.NE.'TEP' ) GO TO 150
BACKSPACE 7
C
READ(7,1500) ( IWK( I), WKX( I), I = 1, 3 )
C
C --------------------------------------
DO 990 I = 1, 3
C
IF ( IWK( I).EQ.0 .AND. DABS(WKX( I)).LE.0.0001D0) GO TO 990
990 CONTINUE
C --------------------------------------
GO TO 141
C
C ----- WALL ---
C
150 CONTINUE
BACKSPACE 7
C
151 READ(7,1400) HEAD
C
IF ( HEAD.NE.'WAL' ) GO TO 160
BACKSPACE 7
C
READ(7,1600) ( IWK( I),I=1,10)
1600 FORMAT( 10X,10I5)
C --------------------------------------
DO 995 I = 1, 10
C
IF ( IWK( I).EQ.0) GO TO 995
NLUV = NLUV + 1
LBC( NLUV)= IWK( I)
C
995 CONTINUE
C --------------------------------------
GOTO 151
C
160 CONTINUE
C
C ----- To set NNZ, NZ2 --- Attention ---
C
NNZ = 541
NZ2 = 961
C ----------------------------------------
C
LUV = NLUV
C
IF ( NNP.EQ.121) CALL BC10 ( DSH, SH, NTB, NBD, Y )
C
C --------------------------------------
DO 997 I = 2, NR2-1
C
HXX( I) = X( I-1) - X( I)
HYY( I) = Y( ( I-1)* ( NR2-1)+1) - Y( ( I-2)* ( NR2-1)+1)
C
997 CONTINUE
C --------------------------------------
C
HXX( 1) = HXX( 2)
HYY( 1) = HYY( 2)
HXX( NR2) = HXX( NR2-1)
HYY( NR2) = HYY( NR2-1)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CKAK ( RKX, RKY, U, V, MTR, NR2, HXX, HYY,
1 UN, VN, NDE, FF )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, MTR
C
COMMON / NNL/ NNP, NEL, LUV
C
DIMENSION RKX( NNP), RKY( NNP), U ( NEL), V ( NEL),
1 MTR( NEL), HXX( NR2), HYY( NR2), UN( NNP),
2 VN ( NNP), NDE( NEL,4), FF ( NNP)
C
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF(100)
C
L = 0
C ------------------------------------------------
DO 100 J = 1, NR2-1
C
DO 100 I = 1, NR2-1
C ------------------------------------------------
L = L + 1
C
AA = ( HXX( I+1)* HXX( I+1)
1 / ( HXX( I+1)* HXX( I+1) + HYY( J)* HYY( J)))* HYY( J)* HYY( J)
C
BB = DT* ( UN( L)/ HXX( I+1) + VN( L)/ HYY( J))** 2/ 2.D0
C
CC = ( ( HXX( I)/ HXX( I+1) - 1.D0 )* UN( L)/ HXX( I+1) +
1 ( HYY( J+1)/ HYY( J) - 1.D0 )* VN( L)/ HYY( J))/ 6.D0
C
FF( L) = AA* ( BB + CC )
C
100 CONTINUE
C ------------------------------------------------
C
DO 200 I = 1, NEL
C
C ----- Attention --- Modified Diffusivity ---
C
RKX( I) = COF( MTR( I))
C
1 + 0.25D0* ( FF( NDE( I,1)) + FF( NDE( I,2))
2 + FF( NDE( I,3)) + FF( NDE( I,4)))
C
RKY( I) = RKX( I)
C
200 CONTINUE
C -------------------------------------------------
C
DO 300 I = 1, NNP
C
FF( I) = FF( I)/ COF( 1) + 1.D0
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPPS ( U, V, UN, VN, PSI, PSO, IBP,
1 VOR, DP, RKX, RKY, NCN, FVC, NC2,
2 KCC, NDE, MTR, X, Y, IBV, NNZ,
3 NZ2, NR2, HXX, HYY, FF, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 IBP, IBV, KCC, NDE, MTR, NCN
C
DIMENSION U ( NEL), V ( NEL), UN ( NNP), VN ( NNP),
1 PSI( NNP), PSO( NNP), IBP( NNP), VOR( NNP),
2 DP ( NNZ), RKX( NNP), RKY( NNP), NCN( NNP),
3 FVC( NNP), NC2( 2,NZ2), KCC( NNP), NDE( NEL,4),
4 MTR( NEL), X ( NNP), Y ( NNP), IBV( NNP),
5 FF ( NNP), ARW( NNP)
C
DIMENSION HXX( NR2), HYY( NR2)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF(100)
C
C ----- To construct F & solve equation -----
C
CALL FVOR ( VOR, NDE, FVC, X, Y )
C
C -------------------------------------------
C
CALL MOVR ( PSO, 1.D0, PSI, NNP )
C
C ----------------------------------------------------------------
C
CALL DSLV ( FVC, PSI, DP, NC2, KCC, IBP, NNZ, NZ2, OMG( 1),
1 ERR( 1) )
C
C ----------------------------------------------------------------
C
CALL VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV, NCN, ARW )
C
C ----- To compute RKX & RKY ---
C
CALL CKAK ( RKX, RKY, U, V, MTR, NR2, HXX, HYY, UN, VN, NDE, FF )
C
C ----------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LPVR ( UN, VN, VOR, VRO, IBV, D, P2, FVC, NC2,
1 KCC, NDE, X, Y, NNZ, NZ2, DMX, DN )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 KCC, IBV, NDE
C
DIMENSION UN ( NNP), VN ( NNP), X ( NNP), Y ( NNP),
1 VOR( NNP), VRO( NNP), IBV( NNP), D ( NNZ),
2 P2 ( NNZ), FVC( NNP), NC2( 2,NZ2), KCC( NNP),
3 NDE( NEL,4), DS ( 121)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
COMMON /FLW/ NVC, IVR( 50), VRC( 50), NPC, IPC( 50), PSC( 50),
1 COF( 100)
C
C ---- To construct F ---
C
CALL MQV12 ( VOR, FVC, UN, VN, NDE, X, Y, DMX, DS)
C
C ----------------------------------------------------------
DN = 0.D0
DO 100 I = 1, NNP
C
DN = DN + DABS( DS( I))
100 CONTINUE
DN = DN/ NNP
C ----------------------------------------------------------
C
CALL FVCTN ( D, VRO, FVC, NC2, KCC, NNZ, NZ2 )
C
C ----- To solve simultaneous equation ---------------------
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 DSLV ( FVC, VOR, P2, NC2, KCC, IBV, NNZ, NZ2, OMG( 2),
1 ERR( 2))
C
C -----------------------------------------------------------------
RETURN
END
C **********************************************************************
C
SUBROUTINE MODD ( D, P1, P2, X, Y, RKX, RKY, NCW, NDE,
1 NNZ, NDM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION D ( NNZ), P1 ( NNZ), P2 ( NNZ), X( NNP), Y( NNP),
1 RKX( NNP), RKY( NNP), NCW( NNZ), NDE( NEL,4)
C
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /NNL/ NNP, NEL, LUV
C
C ----- To modify D -----------------------------------------
C
CALL DMAT ( X, Y, RKX, RKY, D, NDE, NCW, NNZ, NDM )
C
C -------------------------------------------
W1 = 1.D0 / DT
DO 100 I = 1, NNZ
C
P2( I) = ALP* D( I) + P1( I)* W1
D ( I) = P2( I) - D( I)
C
100 CONTINUE
C -------------------------------------------
RETURN
END
C **********************************************************************
C
SUBROUTINE MOVI ( 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 MOVR ( 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 ( C, F, UN, VN, NDE, X, Y, DMX, DS)
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE( NEL,4)
C
DIMENSION C( NNP), F( NNP), UN( NNP), VN( NNP), X( NNP),
1 Y( NNP), DS( NNP)
C
C ----- Attention : '100' elements ---
C
DIMENSION WNU( 100,16), WNV( 100,16)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
COMMON /CST/ OMG( 3), ERR( 3), DLT( 3), ACF, RE
C
C --------------------------------------
C
CALL VXG ( F, 0.D0, NNP)
C
CALL VXG ( DS, 0.D0, NNP)
C
C --------------------------------------
C
C ----- PSI & OMG : Linear ---
C
C --------------------------------------------------------------------
DO 100 NE = 1, NEL
C
C1 = C( NDE( NE,1))
C2 = C( NDE( NE,2))
C3 = C( NDE( NE,3))
C4 = C( NDE( NE,4))
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
C
WNU( NE, 1) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 2) = - A/12.D0*C1 + A/12.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 3) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 4) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
C
WNU( NE, 5) = WNU( NE,2)
WNU( NE, 6) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 7) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 8) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
C
WNU( NE, 9) = WNU( NE,3)
WNU( NE,10) = WNU( NE,7)
WNU( NE,11) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
WNU( NE,12) = - A/36.D0*C1 + A/36.D0*C2 + A/12.D0*C3 - A/12.D0*C4
C
WNU( NE,13) = WNU( NE,4)
WNU( NE,14) = WNU( NE,8)
WNU( NE,15) = WNU( NE,12)
WNU( NE,16) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
C
100 CONTINUE
C --------------------------------------------------------------------
DO 200 NE = 1, NEL
C
C1 = C( NDE( NE,1))
C2 = C( NDE( NE,2))
C3 = C( NDE( NE,3))
C4 = C( NDE( NE,4))
C
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
WNV( NE, 1) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
WNV( NE, 2) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
WNV( NE, 3) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
WNV( NE, 4) = - B/12.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/12.D0*C4
C
WNV( NE, 5) = WNV( NE, 2)
WNV( NE, 6) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3+B/18.D0*C4
WNV( NE, 7) = - B/36.D0*C1 - B/12.D0*C2 + B/12.D0*C3 + B/36.D0*C4
WNV( NE, 8) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
C
WNV( NE, 9) = WNV( NE, 3)
WNV( NE,10) = WNV( NE, 7)
WNV( NE,11) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE,12) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
C
WNV( NE,13) = WNV( NE, 4)
WNV( NE,14) = WNV( NE, 8)
WNV( NE,15) = WNV( NE,12)
C
WNV( NE,16) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
C
200 CONTINUE
C --------------------------------------------------------------------
DO 300 IEL = 1, NEL
C
C ----------------------------------------------------------
DO 320 I = 1, 4
C
II = NDE( IEL,I)
FU = 0.D0
FV = 0.D0
C
C ------------------------------------------------
DO 350 J = 1, 4
C
JJ = NDE( IEL,J)
FU = FU + WNU( IEL,4*( I-1)+J)* UN( JJ)
FV = FV + WNV( IEL,4*( I-1)+J)* VN( JJ)
C
350 CONTINUE
C ------------------------------------------------
F( II) = F( II) + FU + FV
C
320 CONTINUE
C ----------------------------------------------------------
C
300 CONTINUE
C --------------------------------------------------------------------
DO 400 NE = 1, NEL
C
C1 = UN( NDE( NE,1))
C2 = UN( NDE( NE,2))
C3 = UN( NDE( NE,3))
C4 = UN( NDE( NE,4))
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))/ 2.D0
C
WNU( NE, 1) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 2) = - A/12.D0*C1 + A/12.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 3) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
WNU( NE, 4) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
C
WNU( NE, 5) = WNU( NE,2)
WNU( NE, 6) = - A/ 6.D0*C1 + A/ 6.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 7) = - A/18.D0*C1 + A/18.D0*C2 + A/18.D0*C3 - A/18.D0*C4
WNU( NE, 8) = - A/36.D0*C1 + A/36.D0*C2 + A/36.D0*C3 - A/36.D0*C4
C
WNU( NE, 9) = WNU( NE,3)
WNU( NE,10) = WNU( NE,7)
WNU( NE,11) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
WNU( NE,12) = - A/36.D0*C1 + A/36.D0*C2 + A/12.D0*C3 - A/12.D0*C4
C
WNU( NE,13) = WNU( NE,4)
WNU( NE,14) = WNU( NE,8)
WNU( NE,15) = WNU( NE,12)
WNU( NE,16) = - A/18.D0*C1 + A/18.D0*C2 + A/ 6.D0*C3 - A/ 6.D0*C4
C
400 CONTINUE
C --------------------------------------------------------------------
C
DO 500 NE = 1, NEL
C
C1 = VN( NDE( NE,1))
C2 = VN( NDE( NE,2))
C3 = VN( NDE( NE,3))
C4 = VN( NDE( NE,4))
C
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))/ 2.D0
C
WNV( NE, 1) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
WNV( NE, 2) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
WNV( NE, 3) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
WNV( NE, 4) = - B/12.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/12.D0*C4
C
WNV( NE, 5) = WNV( NE,2)
WNV( NE, 6) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE, 7) = - B/36.D0*C1 - B/12.D0*C2 + B/12.D0*C3 + B/36.D0*C4
WNV( NE, 8) = - B/36.D0*C1 - B/36.D0*C2 + B/36.D0*C3 + B/36.D0*C4
C
WNV( NE, 9) = WNV( NE,3)
WNV( NE,10) = WNV( NE,7)
WNV( NE,11) = - B/18.D0*C1 - B/ 6.D0*C2 + B/ 6.D0*C3 + B/18.D0*C4
WNV( NE,12) = - B/18.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/18.D0*C4
C
WNV( NE,13) = WNV( NE,4)
WNV( NE,14) = WNV( NE,8)
WNV( NE,15) = WNV( NE,12)
WNV( NE,16) = - B/ 6.D0*C1 - B/18.D0*C2 + B/18.D0*C3 + B/ 6.D0*C4
C
500 CONTINUE
C --------------------------------------------------------------------
DO 600 IEL = 1, NEL
C
C ----------------------------------------------------------
DO 620 I = 1, 4
C
II = NDE( IEL,I)
FU = 0.D0
FV = 0.D0
DUV = 0.D0
C
C ------------------------------------------------
DO 650 J = 1, 4
C
JJ = NDE( IEL,J)
FU = FU + WNU( IEL,4*( I-1)+J)* C( JJ)
FV = FV + WNV( IEL,4*( I-1)+J)* C( JJ)
DUV = DUV + WNU( IEL,4*( I-1)+J) + WNV( IEL,4*( I-1)+J)
C
650 CONTINUE
C -----------------------------------------------
F( II) = F ( II) + FU + FV
DS( II) = DS( II) + DUV
C
620 CONTINUE
C ----------------------------------------------------------
C
600 CONTINUE
C --------------------------------------------------------------------
DO 700 I = 1, NEL
C
A = DABS( Y( NDE( I,1)) - Y( NDE( I,3)))
B = DABS( X( NDE( I,1)) - X( NDE( I,3)))
N1 = NDE( I,1)
N2 = NDE( I,2)
N3 = NDE( I,3)
N4 = NDE( I,4)
C
700 CONTINUE
C --------------------------------------------------------------------
C
DMX = 0.D0
DO 800 I = 1, NNP
C
DMX = DMAX1( DMX, DABS( DS( I)))
800 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE NBCN ( SH, DSH, PSI, VOR, NBD, LBC, NTB, NNP, LUV )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NBD, LBC, NTB
C
DIMENSION SH ( LUV), DSH( LUV), PSI( NNP), VOR( NNP),
1 NBD( LUV), LBC( LUV), NTB( 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 = NBD( I)
N1 = NTB( I,1)
C
IF ( NE.EQ.0) VOR( ND) = 0.D0
IF ( NE.EQ.0) GO TO 100
APSI = PSI( N1)
AVOR = VOR( N1)
C
VOR( ND) = DSH( I)* ( PSI( ND) - APSI) - AVOR* 0.5D0
C
C ----- Attention : Temp. Routine --- 121/ 100 ---
C
C --- Velocity : Right direction ==> ---
C
IF ( ND.GE.112.AND.ND.LE.120 ) VOR( ND) = VOR( ND) - SH( I)
C -----------------------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMAT ( P, NCC, NCW, X, Y, NDE, KCC, NC2, NNZ,
1 NZ2, ARS)
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, KCC
C
DIMENSION P ( NNZ), NCC( NZ2), NCW( NNZ), X( NNP), Y( NNP),
1 NDE( NEL,4), KCC( NNP), NC2( 2,NZ2)
C
DIMENSION A1( 4,4), N( 4), ARS( NEL)
C
COMMON /NNL/ NNP, NEL, LUV
C
W1D9 = 1.D0/ 9.D0
W1D18 = 1.D0/ 18.D0
W1D36 = 1.D0/ 36.D0
C
A1( 1,1) = W1D9
A1( 1,2) = W1D18
A1( 1,3) = W1D36
A1( 1,4) = W1D18
C
A1( 2,1) = W1D18
A1( 2,2) = W1D9
A1( 2,3) = W1D18
A1( 2,4) = W1D36
C
A1( 3,1) = W1D36
A1( 3,2) = W1D18
A1( 3,3) = W1D9
A1( 3,4) = W1D18
C
A1( 4,1) = W1D18
A1( 4,2) = W1D36
A1( 4,3) = W1D18
A1( 4,4) = W1D9
C
C ---------------------------------
C
CALL VXG ( P, 0.D0, NNZ )
C
C ---------------------------------
C
NZK = 0
C ----------------------------------------------------------
DO 100 IEL = 1, NEL
C
N ( 1) = NDE( IEL,1)
N ( 2) = NDE( IEL,2)
N ( 3) = NDE( IEL,3)
N ( 4) = NDE( IEL,4)
C
ARS( IEL) = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))
1 * DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))
C
C ------------------------------------------------
DO 150 I = 1, 4
DO 150 J = 1, 4
C ------------------------------------------------
C
II = N ( I)
JJ = N ( J)
A = DABS( Y( NDE( IEL,1)) - Y( NDE( IEL,3)))
B = DABS( X( NDE( IEL,1)) - X( NDE( IEL,3)))
AV = A* B* A1( I,J)
C
C ------------------------------------------------
C
CALL DCNL ( AV, II, JJ, P, NCC, NZK, NNZ )
C
C ------------------------------------------------
C
150 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C --------------------------------------------------
DO 200 II = 1, NZ2
C
NC2( 1,II) = NCC( II)/ 10300
NC2( 2,II) = NCC( II) - 10300* NC2( 1,II)
C
200 CONTINUE
C --------------------------------------------------
C
CALL MOVI ( NCC, 1, NCW, NNZ)
C
C --------------------------------------------------
C
CALL DCIJ ( KCC, NC2, NPWK, NZK, NNP, NZ2 )
C
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
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
RETURN
END
C **********************************************************************
C
SUBROUTINE TOUT ( PSI, VOR, UN, VN, X1, Y1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PSI( NNP), VOR( NNP), UN( NNP), VN( NNP), X1( NNP),
1 Y1 ( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
COMMON /CS2/ DT, NLM, LOP, ALP
C
PMX = 0.D0
VMA = 0.D0
NMP = 0
NVP = 0
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( PMX.LE. PSI( I)) THEN
PMX = PSI ( I)
NMP = I
END IF
C
IF ( VMA .LE. DABS( VOR( I))) THEN
VMA = VOR ( I)
NVP = I
END IF
C
100 CONTINUE
C ------------------------------------------------
C
M = NNP - INT( DSQRT( DFLOAT( NEL)))/ 2
N = INT( DSQRT( DFLOAT( NEL))) + 1
C ------------------------------------------------
DO 200 I = M, 1, -N
C
IF ( MOD( NEL,2)) 200, 145, 140
140 UX = ( UN( I)+UN( I-1))/2.D0
WRITE(8,2000) - UX, Y1( I), I
GOTO 200
145 WRITE(8,2000) - UN( I), Y1( I), I
2000 FORMAT( E15.7,1X,E15.7,2X,I5 )
C
200 CONTINUE
C ------------------------------------------------
C
K = ( INT( DSQRT( DFLOAT( NEL)))/ 2.D0 )* N + 1
L = K + INT( DSQRT( DFLOAT( NEL)))
C
C ------------------------------------------------
DO 300 I = K, L
C
XX1 = 1.D0 - X1( I)
IF ( MOD( NEL,2)) 300, 245, 240
240 VX = ( VN( I) + VN( I+N))/ 2.D0
GO TO 300
245 CONTINUE
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DQTM ( A, B, N, S )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( NNP), B( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
S = 0.D0
DO 100 I = 1, N-1
C
H = DABS( A( I) - A( I+1))
S = S + H* ( B( I) + B( I+1))* 0.5D0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VELC ( PSI, U, V, NDE, UN, VN, X, Y, IBV,
1 NCN, ARW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NCN, IBV
C
DIMENSION PSI( NNP), U ( NEL), V ( NEL), NDE( NEL,4),
1 UN ( NNP), VN ( NNP), X ( NNP), Y ( NNP),
2 NCN( NNP), IBV( NNP), ARW( NNP)
C
COMMON /NNL/ NNP, NEL, LUV
C
CALL VXG ( UN, 0.D0, NNP)
C
CALL VXG ( VN, 0.D0, NNP)
C
CALL VXG ( ARW, 0.D0, NNP)
C
C ----- To compute UN & VN ---
C
C ----------------------------------------------------------
DO 100 NE = 1, NEL
C
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))
AB = A* B
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
C
ARW( N1) = ARW( N1) + AB
ARW( N2) = ARW( N2) + AB
ARW( N3) = ARW( N3) + AB
ARW( N4) = ARW( N4) + AB
C
SN1 = 0.D0
SN2 = X( N2) - X( N1)
SN3 = X( N2) - X( N1)
SN4 = 0.D0
C
TN1 = 0.D0
TN2 = 0.D0
TN3 = Y( N3) - Y( N2)
TN4 = Y( N3) - Y( N2)
C
UN( N1) = UN( N1) + ( ( SN1-B )* PSI( N1) - SN1* PSI( N2)
1 + SN1* PSI( N3) + ( B-SN1)* PSI( N4))
C
UN( N2) = UN( N2) + ( ( SN2-B )* PSI( N1) - SN2* PSI( N2)
1 + SN2* PSI( N3) + ( B-SN2 )* PSI( N4))
C
UN( N3) = UN( N3) + ( ( SN3-B )* PSI( N1) - SN3* PSI( N2)
1 + SN3* PSI( N3) + ( B-SN3 )* PSI( N4))
C
UN( N4) = UN( N4) + ( ( SN4-B )* PSI( N1) - SN4* PSI( N2)
1 + SN4* PSI( N3) + ( B-SN4)* PSI( N4))
C
C --------------------------------------------------------------------
C
VN( N1) = VN( N1) - ( ( TN1-A )* PSI( N1) + ( A-TN1)* PSI( N2)
1 + TN1* PSI( N3) - TN1* PSI( N4))
C
VN( N2) = VN( N2) - ( ( TN2-A )* PSI( N1) + ( A-TN2 )* PSI( N2)
1 + TN2* PSI( N3) - TN2* PSI( N4))
C
VN( N3) = VN( N3) - ( ( TN3-A )* PSI( N1) + ( A-TN3 )* PSI( N2)
1 + TN3* PSI( N3) - TN3* PSI( N4))
C
VN( N4) = VN( N4) - ( ( TN4-A )* PSI( N1) + ( A-TN4)* PSI( N2)
1 + TN4* PSI( N3) - TN4* PSI( N4))
C
100 CONTINUE
C ----------------------------------------------------------
DO 200 I = 1, NNP
C
UN( I) = UN( I)/ ARW( I)
VN( I) = VN( I)/ ARW( I)
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 NE = 1, NEL
C
N1 = NDE( NE,1)
N2 = NDE( NE,2)
N3 = NDE( NE,3)
N4 = NDE( NE,4)
A = DABS( Y( NDE( NE,1)) - Y( NDE( NE,3)))
B = DABS( X( NDE( NE,1)) - X( NDE( NE,3)))
S = B* 0.5D0
T = A* 0.5D0
AB = A* B
C
U( NE) = ( ( S - B)* PSI( N1) - S* PSI( N2)
1 + S* PSI( N3) + ( B-S)* PSI( N4))/ AB
C
V( NE) = - ( ( T - A)* PSI( N1) + ( A - T)* PSI( N2)
1 + T* PSI( N3) - T* PSI( N4))/ AB
C
300 CONTINUE
C ----------------------------------------------------------
DO 400 NE = 1, NEL
DO 400 II = 1, 4
C
JJ = NDE( NE,II)
IF ( IBV( JJ).EQ.9999) UN( JJ) = 0.D0
IF ( IBV( JJ).EQ.9999) VN( JJ) = 0.D0
C
400 CONTINUE
C ----------------------------------------------------------
C
C ----- Attention : Moving boundary : Temp. Routine --- 121/ 100 ---
C
C ----------------------------------------------------------
DO 500 ND = 1, NNP
C
IF ( ND.GE.112.AND.ND.LE.120) UN( ND) = 1.D0
500 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ********************************************************************