–ß‚é
C
PROGRAM MAIN
C ******************************************************************** C
C C
C Tidal Flow & Diffusion Analysis C
C C
C by Finite Element Method C
C C
C TIFDAS ( V.5.4 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
INTEGER*2 NODE, NX1, NW1, MX1, MT1, KTR, KCD, NBU, NBV,
1 NBZ, LCL, NEU, NMU
C
COMMON /TEST/ KEY
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST,
1 NL3, NP12, NL12
C
DIMENSION U ( 81), V ( 81), Z ( 81), UO( 81),
1 VO( 81), ZO( 81), UE( 116), VE( 116),
2 ZE( 116), UA( 81), VA( 81), ZA( 81),
3 F ( 81)
C
DIMENSION X ( 81), Y ( 81), NODE( 116,4), AT ( 116),
1 DPH( 116), CFX( 116), CFY ( 116), NX1( 81),
2 MX1( 116), NW1( 98), MT1 ( 140), KTR( 116,3)
C
C ----- Tidal Comp. ---
C
DIMENSION WXX( 116), WYY( 116), BE ( 116,3), CE ( 116,3),
1 PMT( 278), DMT( 278), PX ( 116), PY ( 116),
2 Q ( 116), AVC( 81), BVC( 81), UD ( 81),
3 VD ( 81), ZD ( 81), DEP( 81), FX ( 81),
4 FY ( 81), FZ ( 81), KCD( 81), NBU( 81),
5 NBV( 81), NBZ( 81), NCD( 2,475)
C
DIMENSION FS2( 81), F2 ( 116), F3 ( 116), NC1( 81),
1 NC2( 81,8), B ( 116), PKM( 278), AC1( 278),
2 DL ( 116)
C
C ----- Diffusion Comp. ---
C
DIMENSION UM ( 81), VM( 81), SXR( 116), SYR( 116),
1 FVC( 81), CT( 81), LCL( 81), NEU( 81),
2 NMU( 116), FU( 81), DPM( 81,81)
C
C -------------------------
C
CALL TTLE
C
OPEN ( 9, FILE='REPV.DAT' )
OPEN ( 10, FILE='UV.DAT' )
OPEN ( 11, FILE='PLOT.DAT' )
OPEN ( 12, FILE='SMPL7.DAT' )
C
CALL INIS ( IZ )
C
C ----------------------------------------------------------------------
C
CALL INPT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1,
1 MX1, KTR, NW1, MT1 )
C
C -----------------------------------------------------------------------
C
IF ( IZ.EQ.1) GO TO 333
IF ( IZ.EQ.2) GO TO 555
C
333 CONTINUE
C ---------------------------------------------------------------------
C
CALL TIDL ( X, Y, NODE, AT, DPH, WXX, WYY,
1 BE, CE, PMT, DMT, PX, PY, Q,
2 AVC, BVC, UD, VD, ZD, DEP, FX,
3 FY, FZ, NCD, KCD, NBU, NBV, NBZ,
C
4 NX1, MX1, FS2, F2, F3, NC1, NC2,
5 B, PKM, AC1, DL, F, U, V,
6 Z, UO, VO, ZO, UE, VE, ZE,
7 UA, VA, ZA )
C ---------------------------------------------------------------------
GO TO 900
C
555 CONTINUE
C ---------------------------------------------------------------------
C
CALL SDIF ( X, Y, NODE, AT, DPH, CFX, CFY,
1 UE, VE, U, V, UM, VM, SXR,
2 SYR, NX1, MX1, Q, FVC, CT, NBV,
3 LCL, NEU, NMU, FU, DPM )
C
C ---------------------------------------------------------------------
C
900 CONTINUE
C-----------------------
C
CLOSE ( 9)
CLOSE ( 10)
CLOSE ( 11)
CLOSE ( 12)
C
STOP
END
C ***********************************************************************
C
SUBROUTINE BCVL ( DEP, UA, VA, ZA, U, V, Z, NBU, NBV,
1 NBZ, NNP, T )
C
C ----- Fixed Value at Boundary ---
C
INTEGER*2 NBU( NNP), NBV( NNP), NBZ( NNP)
C
DIMENSION DEP( NNP), U ( NNP), V( NNP), Z( NNP), UA( NNP),
1 VA ( NNP), ZA( NNP)
C
COMMON /FRCM/ NFR, IFTP( 150), IFOR( 150), FA( 150),
1 FT( 150), FR ( 150), FB ( 150)
C
DATA PAI /3.141592/
C
CALL ZVC2 ( NNP, NBU )
CALL ZVC2 ( NNP, NBV )
CALL ZVC2 ( NNP, NBZ )
C
C ----------------------------
DO 100 I = 1, NNP
C
U( I) = UA( I)
V( I) = VA( I)
Z( I) = ZA( I)
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NFR
C
J = IFOR( I)
IF ( J.EQ.0) GO TO 200
IF ( IFTP( I).EQ.1) GO TO 210
IF ( IFTP( I).EQ.2) GO TO 220
IF ( IFTP( I).EQ.3) GO TO 230
GOTO 1000
C
C ----- Forced "Z" ---
C
210 CONTINUE
TEMP = FA(I)* COS( 2.* PAI* T/ FT(I) - PAI* FR(I)/ 180.) + FB(I)
C
IF( NBZ( J).EQ.9999) GOTO 211
Z ( J) = TEMP
NBZ( J) = 9999
GOTO 200
C
211 CONTINUE
Z( J) = Z( J) + TEMP
GOTO 200
C
C ----- Forced "U" ---
C
220 CONTINUE
TEMP = ( FA( I)* COS( 2.* PAI* T/ FT( I)
1 - PAI* FR( I)/180.) + FB( I))* ( DEP( J) + Z( J))
IF ( NBU( J).EQ.9999 ) GOTO 223
U ( J) = TEMP
NBU( J) = 9999
GOTO 200
C
223 CONTINUE
U( J) = U( J) + TEMP
GO TO 200
C
C ----- Forced "V" ---
C
230 CONTINUE
TEMP = ( FA( I)* COS( 2.* PAI* T/ FT( I)
1 - PAI* FR( I)/180.) + FB( I))* ( DEP( J) + Z( J))
C
IF ( NBV( J).EQ.9999 ) GOTO 233
V ( J) = TEMP
NBV( J) = 9999
GOTO 200
C
233 V( J) = V( J) + TEMP
C
200 CONTINUE
C ------------------------------------------------
RETURN
C
1000 WRITE (6,*) ' *** Error *** BCVL ***'
STOP
END
C ***********************************************************************
C
SUBROUTINE BMTX ( PMT, DMT, A, NNZ, C1, C2 )
C
C ----- Left hand side for Expression "B" ---
C
DIMENSION PMT( NNZ), DMT( NNZ), A( NNZ), C1( NNZ)
C
DO 100 I = 1, NNZ
C
A( I) = C1( I)* PMT( I) + C2* DMT( I)
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CALC ( X, Y, NODE, CFX, CFY, KTR, LNE )
C
INTEGER*2 NODE, KTR
C
CHARACTER*4 RPNM
C
COMMON /RPC2/ RPNM( 50)
COMMON /TITE/ TITL1( 10)
COMMON /SECM/ NSA, ISEA( 100)
COMMON /FXCM/ NFXC, IFXC(200), FXC(200)
COMMON /LDCM/ NLND, ILND(400), DLND(400)
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
C
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
COMMON /DPLT/ XMIN, XMAX, YMIN, YMAX, HIT, DX, DY, DC, WID
COMMON /RPC1/ NRP, IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
C
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT,ZSTT, NPLT, NP3
COMMON /MRCM/ NMTR,COFV( 50),COFX( 50),COFY( 50), WX( 50),WY( 50),
1 RMDA( 50)
C
COMMON /FRCM/ NFR, IFTP( 150), IFOR( 150), FA( 150), FT( 150),
1 FR( 150), FB( 150)
COMMON /CMPF/ MFOR, JFTP( 150), JFOR( 150), PA( 150),
1 PT( 150), PR( 150), PB( 150)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,4), CFX ( NEL),
1 CFY ( NEL), KTR( NEL,3), IND2( 900), LNE( 2,40,4),
2 WK1D( 400)
C
PAI = 3.14159
C
K = 4
C ----------------------------
DO 100 I = 1, NEL
C
L = NODE( I,K)
CFX( I) = COFX( L)
CFY( I) = COFY( L)
C
100 CONTINUE
C ----------------------------
C
L = 0
C ------------------------------------------------
DO 200 I = 1, NEL
C
C --------------------------------------
DO 250 J = 1, 3
C
M1 = NODE( I,J)
J1 = J+1
IF ( J1.GT.3) J1 = 1
M2 = NODE( I,J1)
IF ( M1.LT.M2) GO TO 255
M3 = M2
M2 = M1
M1 = M3
255 M3 = M1* 10000 + M2
IF ( L.EQ.0) GO TO 257
C
C ----------------------------
DO 270 K = 1,L
C
K1 = L - K + 1
IF ( IND2( K1).EQ.M3) GO TO 275
C
270 CONTINUE
C ----------------------------
C
257 CONTINUE
L = L + 1
IND2( L) = M3
C
275 CONTINUE
C
250 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
NNZ = NNP + L
C --------------------------------------
C
CALL TRAN ( NEL, NODE, KTR )
C
C --------------------------------------
MB = 0
K2 = 3
J2 = K2 - 1
C ------------------------------------------------
DO 300 I = 1, NEL
C
IF ( J2.LE.1) J2 = 1
C --------------------------------------
DO 300 J = 1, J2
C
K1 = J+1
IF ( K2.LE.K1) K2 = K1
C --------------------------------------
DO 300 K = K1,K2
C
M1 = NODE( I,J) - NODE( I,K)
IF ( M1.LT.0 ) M1 = - M1
IF ( M1.GT.MB) MB = M1
C
300 CONTINUE
C --------------------------------------
C
MB = MB + 1
NZ2 = NNZ* 2 - NNP
C --------------------------------------
C
CALL ZVC4 ( 400, WK1D )
C
C --------------------------------------
C
C --------------------------------------
DO 400 I = 1, 40
C
IF ( LNE( 1,I,1).EQ.0) GO TO 550
J1 = LNE( 1,I,3)
J2 = J1 + 1
J4 = LNE( 1,I,4)
J3 = J4-1
WK9 = 0.
IF ( J3.LE.J1) J3 = J1
C
C ----------------------------
DO 450 J = J1, J3
C
WK1 = ( X( ISEA( J+1)) - X( ISEA( J)))** 2
WK2 = ( Y( ISEA( J+1)) - Y( ISEA( J)))** 2
WK1D( J+1) = SQRT( WK1 + WK2 )
WK9 = WK9 + WK1D( J+1)
C
450 CONTINUE
C ----------------------------
WK1 = FA( J4) - FA( J1)
WK2 = FT( J4) - FT( J1)
WK3 = FR( J4) - FR( J1)
WK4 = FB( J4) - FB( J1)
C
IF ( J3.LE.J2) J3 = J2
C
C ----------------------------
DO 470 J = J2, J3
C
WK8 = WK1D( J)/ WK9
IFTP( J) = IFTP( J1)
IFOR( J) = ISEA( J )
C
FA( J) = FA( J-1) + WK1* WK8
FT( J) = FT( J-1) + WK2* WK8
FR( J) = FR( J-1) + WK3* WK8
FB( J) = FB( J-1) + WK4* WK8
C
470 CONTINUE
C ----------------------------
C
400 CONTINUE
C --------------------------------------
550 CONTINUE
C
NFR = NSA
C
C ----- ADD /PFORCE/ to /FORCE/ ---
C
IF ( MFOR.LE.0) GO TO 570
C --------------------------------------
DO 600 I = 1, MFOR
C
IF ( NFR.GE.150) THEN
WRITE(6,*) 'ERR,***CALC***'
STOP
END IF
C
NFR = NFR + 1
IFTP( NFR) = JFTP( I)
IFOR( NFR) = JFOR ( I)
C
FA( NFR) = PA( I)
FT( NFR) = PT( I)
FR( NFR) = PR( I)
FB( NFR) = PB( I)
C
600 CONTINUE
C --------------------------------------
C
570 CONTINUE
C
C ----- LAND ---
C
CALL ZVC4 ( 400, WK1D )
C
C --------------------------------------
DO 700 I = 1, 40
C
IF ( LNE( 2,I,1).EQ.0) GO TO 730
J1 = LNE( 2,I,3)
J4 = LNE( 2,I,4)
J2 = J1 + 1
J3 = J4 - 1
C
IF ( J3.LE.J1) J3 = J1
C
C -----------------------------
DO 750 J = J1, J3
C
WKX = X( ILND( J+1)) - X( ILND( J))
WKY = Y( ILND( J+1)) - Y( ILND( J))
WK1D( J) = ATAN2( WKY, WKX )
C
750 CONTINUE
C ----------------------------
C
IF ( LNE( 2, I, 2 ).LE.2) GO TO 700
IF ( ILND( J1).EQ.ILND( J4)) GO TO 850
C
DLND( J1) = WK1D( J1)
DLND( J4) = WK1D( J3)
C
IF ( DLND( J1).LT.0.) DLND( J1) = DLND( J1) + PAI
IF ( DLND( J4).LT.0.) DLND( J4) = DLND( J4) + PAI
IF ( J3.LE.J2 ) J3 = J2
C
C ----------------------------
DO 800 J = J2, J3
C
DLND( J) = ( WK1D( J-1) + WK1D( J))/ 2.
IF ( DLND( J).LT.0.) DLND( J) = DLND( J) + PAI
C
800 CONTINUE
C ----------------------------
GO TO 700
C
C ----------------------------
850 DO 900 J = J1, J3
C
J5 = J - 1
IF ( J.EQ.J1) J5 = J3
DLND( J) = ( WK1D( J5 ) + WK1D( J))/2.
IF ( DLND( J).LT.0.) DLND( J) = DLND( J) + PAI
C
900 CONTINUE
C ----------------------------
C
DLND( J4) = DLND( J1)
C
700 CONTINUE
C --------------------------------------
C
730 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CFAS ( U, V, UAV, VAV, UM, VM, SXR, SYR,
1 NODE, MX1, NNP, NEL )
C
C ----- MEAN Velocity & VARIANCE SET ---
C
INTEGER*2 NODE( NEL,4), MX1( NEL)
REAL*4 NTIM, MTM
C
DIMENSION U ( NNP), V ( NNP), UAV( NEL), VAV( NEL), UM( NNP),
1 VM( NNP), SXR( NEL), SYR( NEL)
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
C
C ----- CLEAR ---
C
CALL ZVC4 ( NEL, UAV )
CALL ZVC4 ( NEL, VAV )
CALL ZVC4 ( NNP, UM )
CALL ZVC4 ( NNP, VM )
C
CALL ZVC4 ( NNP, U )
CALL ZVC4 ( NNP, V )
CALL ZVC4 ( NEL, SXR )
CALL ZVC4 ( NEL, SYR )
C
MTM = 0
IF ( KYP.EQ.9) GO TO 990
C
MTM = MTM - 43200
IF ( KYP.EQ.2 ) MTM = MTM - 43200
C
C ----- Mean Velocity ---
C
NTIM = 0
C
10 CONTINUE
C
IF( NTIM.LE.MTM ) GO TO 10
C
NTIM = NTIM + 1
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
N1 = NODE( I,1)
N2 = NODE( I,2)
N3 = NODE( I,3)
UAV( I) = UAV( I) + ( U( N1) + U( N2 ) + U( N3))/ 3.
VAV( I) = VAV( I) + ( V( N1) + V( N2 ) + V( N3))/ 3.
C
100 CONTINUE
C ----------------------------------------------------------
DO 110 I = 1, NNP
C
UM( I) = UM( I) + U( I)
VM( I) = VM( I) + V( I)
C
110 CONTINUE
C ----------------------------------------------------------
GO TO 10
C
C ----- Variance CALC. ---
C
300 CONTINUE
C
IF( NTIM.LE.MTM ) GO TO 300
C
C --------------------------------------
DO 310 I = 1, NEL
C
N1 = NODE( I, 1)
N2 = NODE( I, 2 )
N3 = NODE( I, 3)
C
UE = ( U( N1) + U( N2 ) + U( N3))/ 3.
VE = ( V( N1) + V( N2 ) + V( N3))/ 3.
SXR( I) = SXR( I) + ( UE - UAV( I))* ( UE - UAV( I))
SYR( I) = SYR( I) + ( VE - VAV( I))* ( VE - VAV( I))
C
310 CONTINUE
C --------------------------------------
GO TO 300
C
990 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CLCF ( NEL, AT, NODE, SXR, SYR, MX1, CFX, CFY,
1 CFL, NMC )
C
C ----- Calc. of Coefficients ---
C
INTEGER*2 NODE( NEL,4), MX1( NEL)
C
DIMENSION AT( NEL), SXR( NEL), SYR( NEL), CFX( NEL), CFY( NEL)
C
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
C
COMMON /TEST/ KEY
C
CALL ZVC4 ( NEL,CFX)
C
CALL ZVC4 ( NEL,CFY)
C
IF ( ABS( CFL).LE.0.001) GO TO 150
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
CARL = 2.* SQRT( AT( I)/1.7320508 )
CFX( I) = CFL* CARL* SXR( I)
CFY( I) = CFL* CARL* SYR( I)
C
100 CONTINUE
C ------------------------------------------------
C
150 CONTINUE
C
C ------------------------------------------------
DO 200 I = 1, NEL
C
CFX( I) = CFX( I) + COFX( NODE( I,NMC))
CFY( I) = CFY( I) + COFY( NODE( I,NMC))
C
200 CONTINUE
C ------------------------------------------------
C
IF ( KEY.NE.1) GO TO 250
C
C ----- Attention ---
C
AL = 0.86
C -----------------------------------
DO 300 I = 1, NEL
C
CFX( I) = CFX( I) + AL* 3.75
CFY( I) = CFY( I) + AL* 3.75
C
300 CONTINUE
C -----------------------------------
C
250 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CNDE ( FS2, F2, NC2, NC1 )
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
DIMENSION FS2( NNP), F2( NEL), NC1( NNP), NC2( NNP,8)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
NQ = NC1( I)
FNCN = 1./ FLOAT( NQ )
CTT = 0.
C
C --------------------------------------
DO 200 J = 1, NQ
C
JJ = NC2( I,J)
CTT = CTT + F2( JJ)
C
200 CONTINUE
C --------------------------------------
FS2( I) = CTT* FNCN
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE COEF ( NNP, NEL, NODE, G, H, DT, DL, FS2, F3,
1 F2, NC1, NC2, B )
C
INTEGER*2 NODE
DIMENSION NODE( NEL,4), H ( NEL), DL ( NEL), FS2( NNP),
1 F3 ( NEL), F2( NEL), NC1( NNP), NC2( NNP, 8),
2 B ( NEL)
C
C ------------------------------------------------
DO 100 I = 1, NEL
IF ( H( I).GT.RHMX) RHMX = H( I)
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NEL
C
B( I) = ( SQRT( G* RHMX ))* DT/ DL( I)
F = 1./ 6. - 1./( 3.* B( I)* B( I))
C
F3( I) = 0.5* F
F2( I) = 1./ 3. + 1./( 6.* B( I)* B( I))
C
200 CONTINUE
C ------------------------------------------------
C
CALL STNN ( NNP, NEL, NODE, NC1, NC2 )
C
CALL CNDE ( FS2, F2, NC2, NC1 )
C
C ------------------------------------------------
C
BMX = -100.
BMN = 100.
F3MX = -100.
F3MN = 100.
C
FS2X = -100.
FS2I = 100.
F2MX = -100.
F2MN = 100.
C
C ------------------------------------------------
DO 300 I = 1, NEL
C
IF ( B( I).GT.BMX) BMX = B( I)
IF ( B( I).LT.BMN) BMN = B( I)
C
IF ( F3( I).GT.F3MX) F3MX = F3( I)
IF ( F3( I).LT.F3MN) F3MN = F3( I)
C
IF ( F2( I).GT.F2MX) F2MX = F2( I)
IF ( F2( I).LT.F2MN) F2MN = F2( I)
C
300 CONTINUE
C ------------------------------------------------
DO 400 I = 1, NNP
C
IF ( FS2( I).GT.FS2X ) FS2X = FS2( I)
IF ( FS2( I).LT.FS2I.AND.FS2( I).GT.0.) FS2I = FS2( I)
C
400 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DCNIJ ( IA, JA, NJ, NIX, NJX )
C
INTEGER*2 IA( NIX), JA( 2,NJX)
C
CHARACTER*4 ANIX, ANJX
C
C ----- NI, NJ : Current length of IA, JA
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( 2,I)) GO TO 101
IA( NI ) = I-1
NI = JA( 2,I)
IF ( NI.GT.NIX ) GO TO 601
101 JA( 2, 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( 1,K).NE.I ) GO TO 400
K1 = K
GOTO 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 --------------------------------------
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 --------------------------------------
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 ----- Insert ( I,J) element ---
C
IS = IS + 1
JA( 2,IS) = JA( 2,K1)
JA( 1,IS) = J
C
300 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
RETURN
C
601 WRITE(6,2000) ANIX, NIX, ANIX, NI
STOP
602 WRITE(6,2000) ANJX, NJX, ANJX, NJ
STOP
2000 FORMAT(' ** Error ( DCNIJ) ',A3,'(',I5,') too small',A2,'=',I5)
C
END
C ***********************************************************************
C
SUBROUTINE DCNL ( AV, I, J, A, L, K, NT )
C
INTEGER*2 M2( 2)
DIMENSION A( NT), L( NT)
EQUIVALENCE ( M, M2 )
C
C ----- Attention --- NO. OF NODE.LE.32768 ---
C
IF ( I.GT.J) RETURN
C
M2( 1) = J
M2( 2) = I
IE = K
N1 = 1
C
IF ( K.EQ.0) N2 = K + 1
IF ( K.EQ.0) GO TO 10
N2 = K
10 N = N2
IF ( K.EQ.0) IE = 1
C
C ------------------------------------------------
DO 100 II = 1, IE
C
IF ( K.EQ.0) GO TO 520
IF ( M - L( N)) 200, 300, 400
200 CONTINUE
N2 = N
N = N2 - ( N2 - N1)/ 2
IF ( N.NE.N2 ) GO TO 100
IF ( M - L( N1)) 110, 298, 120
C
110 N = N1
120 CONTINUE
K = K + 1
IF ( K.GT.NT ) WRITE(6,805)
805 FORMAT (' * Error K.GT.NT *')
JE = K - N
C --------------------------------------
DO 201 JJ = 1, JE
C
J1 = K - JJ+1
A( J1) = A( J1 - 1)
L( J1) = L( J1 - 1)
C
201 CONTINUE
C --------------------------------------
A( N) = AV
L( N) = M
RETURN
C
298 N = N1
GO TO 300
C
299 N = N2
300 CONTINUE
C
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)) 510, 299, 520
510 N = N2
GOTO 120
C
520 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
RETURN
C
100 CONTINUE
C ------------------------------------------------
C
WRITE(6,2100) K
2100 FORMAT(' * Error --- DO Loop Strange, K =',I10 )
STOP
END
C ***********************************************************************
C
SUBROUTINE DEPW ( H, NODE, DEP, WXX, WYY, NEL, NNP)
C
C ----- Nodal DPH & Wind vector ---
C
INTEGER*2 NODE( NEL,4)
DIMENSION H( NEL), DEP( NNP), WXX( NEL), WYY( NEL)
C
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
N = 0
DEP( I) = 0.0
C
C --------------------------------------
DO 200 J = 1, NEL
DO 200 K = 1, 3
C --------------------------------------
C
IF( NODE( J,K).NE.I ) GO TO 200
DEP( I) = DEP( I) + H( J)
N = N + 1
C
200 CONTINUE
C --------------------------------------
DEP( I) = DEP( I)/ N
C
100 CONTINUE
C ------------------------------------------------
DO 300 I = 1, NEL
C
J = NODE( I,4)
C
WXX( I) = ROA/ ROW* GM2* WX( J)* SQRT( WX( J)** 2 + WY( J)** 2 )
1 * ( 1.+ FRIC )
C
WYY( I) = ROA/ ROW* GM2* WY( J)* SQRT( WX( J)** 2 + WY( J)** 2 )
1 * ( 1.+ FRIC )
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
FUNCTION DVL1 ( U, BE, AT, NODE, NNP, NEL, I )
C
C ----- Differentiated Value ---
C
INTEGER*2 NODE( NEL,4)
DIMENSION U( NNP), BE( NEL,3), AT( NEL)
C
U1 = U( NODE( I,1))
U2 = U( NODE( I,2))
U3 = U( NODE( I,3))
C
DVL1 = ( BE(I,1)* U1 + BE(I,2)* U2 + BE(I, 3)* U3)/( 2.* AT( I))
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE EPRT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1,
1 NW1, MX1, MT1, LNE )
C
INTEGER*2 NODE, NX1, NW1, MX1, MT1
C
CHARACTER*4 RPNM
CHARACTER*1 WD, WKD( 200), BLN, ASTR
C
COMMON /RPC2/ RPNM( 50)
COMMON /TITE/ TITL1( 10)
COMMON /SECM/ NSA, ISEA( 100)
COMMON /FXCM/ NFXC,IFXC(200),FXC(200)
COMMON /LDCM/ NLND, ILND(400), DLND(400)
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
COMMON /RPC1/ NRP, IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
C
COMMON /DPLT/ XMIN, XMAX, YMIN, YMAX, HIT, DX, DY, DC, WID
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
COMMON /FRCM/ NFR, IFTP( 150), IFOR( 150), FA( 150),
1 FT( 150), FR( 150), FB( 150)
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,4), AT ( NEL),
1 DPH ( NEL), CFX ( NEL), CFY ( NEL), LNE( 2,40,4),
2 NX1 ( NNP), NW1 ( NP12), MX1 ( NEL), MT1 ( NL12 ),
3 IWKD( 20), WK2D( 2, 10)
C
DATA BLN /' '/, ASTR /'*'/
C
CALL EPRT1 ( NX1 )
C
IF ( NSR.EQ.0) GO TO 90
C
WRITE(6,2000)
2000 FORMAT (' ---- Source Inf. ---'/
1 4X,'SOR-#',1X,2(' Node Drain-Quant. Soiled-Quant.'))
C
I1 = ( NSR + 3)/ 4
C
C ----------------------------------------------------------
DO 100 I = 1, I1
C
J2 = I* 4
J1 = J2 - 3
IF ( J2.GT.NSR) J2 = NSR
WRITE (6,2100) J1, J2,( NX1( ISOR( J)), SORW( J),
1 SORQ( J),J = J1, J2 )
C
100 CONTINUE
C ----------------------------------------------------------
C
2100 FORMAT (/,' ',I3,'-',I3,2X,2(2X,I4,2F11.4,2X))
C
90 IF ( NFXC.EQ.0) GO TO 333
WRITE(6,2200)
2200 FORMAT(/' ----- FXC Inf. ---'/,4X,'FXC-#',4X,
1 3(' Node Fixed Conc.'))
C
I1 = ( NFXC + 4)/ 5
C
C ------------------------------------------------
DO 200 I = 1, I1
C
J2 = I * 5
J1 = J2 - 4
IF ( J2.GT.NFXC) J2 = NFXC
K = 0
C
C --------------------------------------
DO 250 J = J1, J2
C
K = K + 1
WK2D( 2,K) = FXC( J)
C
IF ( IFXC( J).GT.NNP) GO TO 93
IWKD( K) = NX1( IFXC( J))
WK2D(1,K) = FLOAT( ICHAR( BLN))
GO TO 250
C
93 IWKD( K) = IFXC( J)
WK2D(1,K) = FLOAT( ICHAR ( ASTR))
C
250 CONTINUE
C --------------------------------------
C
WRITE(6,2300) J1, J2, ( IWKD( J),WK2D(1,J),WK2D(2,J),J = 1,K)
2300 FORMAT(/,3X,I3,' -',I3,/,3(5X,I4,1X,F1.0,F7.4,5X))
C
200 CONTINUE
C ------------------------------------------------
C
333 IF ( NFR.EQ.0) GO TO 110
C
C --------------------------------------
DO 106 I = 1, NFR
C
IF ( IFOR( I).GT.NNP) GO TO 102
J = NX1( IFOR( I))
WD = BLN
GO TO 103
102 J = IFOR( I)
WD = ASTR
103 CONTINUE
C
106 CONTINUE
C --------------------------------------
C
110 IF ( NSA.EQ.0) GO TO 120
C
C --------------------------------------
DO 112 I = 1,40
C
IF ( LNE( 1,I,1).EQ.0) GO TO 120
J1 = LNE( 1,I,3)
J2 = LNE( 1,I,4)
C
112 CONTINUE
C --------------------------------------
C
120 IF ( NLND.EQ.0) GO TO 140
WRITE(6,2700)
C
C --------------------------------------
DO 122 I = 1, 40
C
IF ( LNE( 2,I,1).EQ.0) GO TO 140
J1 = LNE( 2,I,3)
J2 = LNE( 2,I,4)
C
122 CONTINUE
C --------------------------------------
C
2700 FORMAT(/' ----- LAND Inf. ---'//9X,' LNE-#',4X,
12('NODE KAKUDO '))
C
140 CONTINUE
C
IF ( NLST.EQ.1.OR.NLST.GE. 3) GO TO 900
C
C ------------------------------------------------
DO 142 I = 1, NNP, 4
C
I1 = I + 3
IF ( I1.GT.NNP) I1 = NNP
J1 = 0
C
C --------------------------------------
DO 145 J = I, I1
C
J1 = J1 + 1
WK2D( 1,J1) = X( J)
WK2D( 2,J1) = Y( J)
IF ( J.GT.NNP) GO TO 146
IWKD( J1) = NX1( J)
WKD ( J1) = BLN
GO TO 145
146 IWKD( J1) = J
WKD ( J1) = ASTR
C
145 CONTINUE
C --------------------------------------
C
142 CONTINUE
C ------------------------------------------------
C
J1 = 3
J2 = 4
C
C ----------------------------------------------------------
DO 210 I = 1, NEL
C
C ------------------------------------------------
DO 202 J = 1, J1
C
IWKD( J) = NODE( I,J)
IF ( J.LE.3) IWKD( J) = NX1( IWKD( J))
C
202 CONTINUE
C ------------------------------------------------
C
IF ( AT( I).LE.0.) THEN
WRITE(6,*) ' Error *** EPRT *** '
STOP
END IF
C
210 CONTINUE
C ----------------------------------------------------------
C
I1 = NL12
IF ( NP12.GT.NL12 ) I1 = NP12
C
C --------------------------------------
DO 222 I = 1,I1
C
C ----------------------------
DO 220 J = 1,4
C
IWKD( J) = 0
C
220 CONTINUE
C ----------------------------
C
IF ( I.LE.NNP ) IWKD( 1) = NX1( I)
IF ( I.LE.NP12) IWKD( 2) = NW1( I)
IF ( I.LE.NEL ) IWKD( 3) = MX1( I)
IF ( I.LE.NL12) IWKD( 4) = MT1( I)
C
222 CONTINUE
C --------------------------------------
C
900 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE EPRT1 ( NX1 )
C
C ----- To edit List Print ---
C
INTEGER*2 NX1
C
CHARACTER*4 RPNM
C
COMMON /TITE/ TITL1( 10)
COMMON /RPC2/ RPNM( 50)
COMMON /SECM/ NSA, ISEA( 100)
COMMON /FXCM/ NFXC,IFXC(200),FXC(200)
COMMON /LDCM/ NLND, ILND(400), DLND(400)
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
C
COMMON /SRCM/ NSR,ISOR( 100), SORW( 100), SORQ( 100)
COMMON /RPC1/ NRP,IRP( 50),UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
COMMON /FRCM/ NFR, IFTP( 150), IFOR( 150), FA( 150), FT( 150),
1 FR( 150), FB( 150)
C
DIMENSION NX1( NNP)
C
WRITE(6,2000) NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST
2000 FORMAT(/' ----- Area Inf. ---',/,
1 2X,'NNP =',I7,' NEL =',I7,' MB =',I7,' NNZ =',I7,
2 ' NZ2 =',I7,/,' SCX =',G14.6,' SCY =',G14.6,
3 ' NLST =',I5)
C
C ----- Tidal Comp.---
C
IF ( KTDE.NE.1) GO TO 40
C
WRITE(6,2100) RDT, NSTR, NED, NP1, NP2, NT1, USTT, ZSTT, NPLT
2100 FORMAT(/' ----- Tidal Computation Inf. ---'/,
1 2X,'Time Step =',F6.1,' (sec)',' Start_T =',I3,' (sec)',
2 ' End_T =',I7,' (sec)'/,2X,'NP1 =',I7,' (sec)',
3 ' NP2 =',I7,' (sec)',' NT1 =',I7,
4 ' (sec)',/,' UV-Start =',F7.3,'(m/s)',1X,' Z-Start =',
5 F7.3,'(m)',' Print Tim.( Yes= 0: No= 1)',I2 )
C
C ---- Diffusion Comp. -------
C
40 IF ( MDF.NE.1) GO TO 50
WRITE(6,2200) KYP, CB, ( I, CFL( I), I = 1, 5 )
2200 FORMAT (/' ----- Diffusion Inf. ---'//
1 9X,'MODEL .......................',I7/
2 9X,'Back ground concentration =',F7.2,'(PPM)'/
3 5(9X,' Rep-Length Coeff.(',I1,') ........',F7.2/))
C
C ----- REP_POINT ---
C
50 IF ( NRP.EQ.0) GO TO 60
WRITE(6,2300)
I1 = ( NRP + 4)/ 5
C
C ----------------------------------------------------------
DO 100 I = 1, I1
C
J2 = I* 5
J1 = J2 - 4
IF ( J2.GT.NRP) J2 = NRP
WRITE(6,2400) J1,J2,( NX1( IRP( J)), RPNM( J),J = J1, J2 )
C
100 CONTINUE
C ----------------------------------------------------------
C
2300 FORMAT(/' ----- Rep. Point Inf. ---'/,
1 9X,'REP-No.',5X,3(' Node Name '))
2400 FORMAT(7X,2I5,5X,3(I5,2X,A4,6X))
C
C ----- MATERIAL ---
C
60 IF ( NMTR.EQ.0) GO TO 70
C
WRITE(6,2500)
C
C ----------------------------------------------------------
DO 200 I = 1, 50
C
IF ( ABS( COFV( I)).LE.0.0001) GO TO 200
WRITE(6,2600) I, COFV( I), COFX( I), COFY( I), WX( I), WY( I),
1 RMDA( I)
C
200 CONTINUE
C ----------------------------------------------------------
C
2500 FORMAT(/' ----- Material Inf. ---'/,
1 4X,'MTR_No COFV COFX COFY',
2 ' WX WY RMDA')
2600 FORMAT(/,5X,I3,3X,6F10.3)
C
70 WRITE(6,2700) FC,GM2,FRIC,VC,ROW,ROA,ALF
2700 FORMAT(/' ----- Constant Inf. ---'/
1 1X,' Coriolis Coeff.=',G14.6,' Bottom Friction Coeff.=',
2 G14.6,/,' Fr. Hirei Coeff.=',G14.6,' Div. Check Vel.=',
3 F5.1,' Sea Density =',F8.2,/,' Air Dens. =',G8.2,
4 ' S.O.R. Par.=',G14.6,' Time Sch. Coeff.=',G10.2,
5 /,' Error const.in S.O.R.=',G14.6,/)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE EXPL ( FVC, C, TSM, NBC, NNP)
C
INTEGER*2 NBC( NNP)
C
DIMENSION FVC( NNP), C( NNP), TSM( NNP)
C
DO 100 I = 1, NNP
C
IF ( NBC( I).EQ.9999 ) GOTO 100
C( I) = FVC( I)/ TSM( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FINT ( F, C, DMT, NCD, KCD, G, NNP, NNZ, NZ2, C1,
1 C2 )
C
INTEGER*2 KCD( NNP), NCD( 2, NZ2 )
C
DIMENSION F( NNP), C( NNP), DMT( NNZ), G( NNP)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( I.EQ.1) GO TO 101
MS = KCD( I-1) + 1
GO TO 102
101 MS = 1
102 ML = KCD( I)
C
IF ( MS.GT.ML) GOTO 900
G( I) = C1 * F( I)
C
C --------------------------------------
DO 150 M = MS, ML
C
M1 = NCD( 2, M )
M2 = NCD(1, M )
G( I) = G( I) + C2* DMT( M1)* C( M2)
C
150 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
RETURN
C
900 WRITE(6,2000)
2000 FORMAT(' *** Error ( FINT )*** MS.GT.ML')
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE FLPN
C
OPEN ( 9, FILE='REPV.DAT' )
OPEN ( 10, FILE='UV.DAT' )
OPEN ( 11, FILE='PLOT.DAT' )
OPEN ( 12, FILE='SMPL7.DAT' )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE FMT1 ( AT, H, Q, NODE, FVC, NNP, NEL)
C
DIMENSION AT( NEL), H( NEL), Q( NEL), FVC( NNP)
C
INTEGER*2 NODE( NEL,4)
C
CALL ZVC4 ( NNP, FVC)
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
ATD3 = AT( I)/ 3.
C
C ------------------------------------------------
DO 150 J = 1, 3
C
FVC( NODE( I,J)) = FVC( NODE( I,J)) + H( I)* Q( I)* ATD3
150 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE GSOL ( A, B, N, MB, MW )
C
DIMENSION A( N, 33), B( N)
C
N1 = N - 1
C
C ----------------------------------------------------------
DO 100 K = 1, N1
C
IF ( K.LE.( N - MB )) THEN
KMAX = MB - 1
ELSE
KMAX = MB - ( K - ( N - MB ))
END IF
C
K1 = K + 1
IF ( K.LE.MB ) THEN
AKK = 1. / A( K, K)
ELSE
AKK = 1. / A( K, MB )
END IF
C
BK = B( K)
KMAX1 = K + MB - 1
IF ( KMAX1.GT.N ) KMAX1 = N
C
C ------------------------------------------------
DO 200 I = K1, KMAX1
C
IF( I.GT.MB ) THEN
KK = K - ( I - MB )
ELSE
KK = K
END IF
C
AIK =A( I, KK)* AKK
B( I) = B( I) - AIK* BK
C
IF( I.LE.MB ) THEN
C --------------------------------------
DO 300 J = K, KMAX + K
C
A( I,J) = A( I,J) - AIK* A( K,J)
300 CONTINUE
C --------------------------------------
END IF
C
IF( I.GT.MB ) THEN
C
DO 400 J = KK, KMAX + KK
C
IF ( K.LT.MB ) THEN
A( I,J) = A( I,J) - AIK* A( K, K + J - KK)
ELSE
A( I,J) = A( I,J) - AIK* A( K,( MB + J - KK))
END IF
C
400 CONTINUE
C
END IF
C
200 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
C ----- BACK Substitution ---
C
K = N
B( K) = B( K)/ A( K, MB )
C
450 K = K - 1
C
IF ( K.LE.0) RETURN
C
K1 = K + 1
SUM = 0.
C
IF ( K.GE.MB ) THEN
MM = MIN0 ( ( N-K+MB), MW )
C
DO 500 J = MB + 1, MM
C
K2 = K1 + J - MB - 1
SUM = SUM + A( K,J)* B( K2)
C
500 CONTINUE
B( K) = ( B( K) - SUM)/ A( K,MB )
C
ELSE
C
DO 600 J = K1, MB+K-1
C
SUM = SUM + A( K,J)* B( J)
600 CONTINUE
C
B( K) = ( B( K) - SUM )/ A( K,K)
END IF
C
GO TO 450
C
END
C ***********************************************************************
C
SUBROUTINE INIS ( IZ )
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
CHARACTER*4 CID
C
WRITE(6,*) '*** 1 : Tidal Comp. 2 : Stationary Diffusion Comp.
1***'
** READ( 5,*) IZ
IZ = 1
REWIND 12
READ(12,1000) CID, NCLC, NNP, NEL, SCX, SCY, NLST
1000 FORMAT(A4,5X,I1,2I5,2F10.0,24X,I1)
REWIND 12
C
IF ( ABS( SCX).LE.0.001) SCX = 1.
IF ( ABS( SCY).LE.0.001) SCY = 1.
C
NL3 = NEL* 3
NP12 = NNP* 1.2
NL12 = NEL* 1.2
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE INIT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1,
1 NW1, MX1, MT1, KTR, LNE )
C
INTEGER*2 NODE, NX1, NW1, MX1, MT1, KTR
C
CHARACTER*4 BLN, RPNM
C
COMMON /TITE/ TITL1( 10)
COMMON /RPC2/ RPNM( 50)
COMMON /SECM/ NSA,ISEA( 100)
COMMON /LDCM/ NLND,ILND(400),DLND(400)
COMMON /FXCM/ NFXC,IFXC(200),FXC(200)
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /SRCM/ NSR,ISOR( 100), SORW( 100), SORQ( 100)
COMMON /RPC1/ NRP,IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /MRCM/ NMTR,COFV( 50),COFX( 50),COFY( 50),
1 WX( 50),WY( 50),RMDA( 50)
C
COMMON /FRCM/ NFR,IFTP( 150),IFOR( 150),FA( 150),
1 FT( 150),FR( 150),FB( 150)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,4), AT ( NEL),
1 DPH( NEL), CFX( NEL), CFY ( NEL), NX1( NNP),
2 NW1( NP12 ), MX1( NEL), MT1( NL12 ), KTR ( NEL,3),
3 LNE( 2,40,4)
C
DATA BLN /' '/
DATA PAI /3.14159/
C
MB = 0
NNZ = 0
NZ2 = 0
C
USTT = 0.
ZSTT = 0.
C
MDF = 0
KYP = 0
CB = 0.
NFAL = 0
C
C ----------------------------
DO 100 I = 1, 5
C
CFL( I) = 0.
100 CONTINUE
C ----------------------------
DO 200 I = 1, 50
C
IRP ( I) = 0
RPNM( I) = BLN
C
200 CONTINUE
C ----------------------------
C
NMTR = 0
C
CALL ZVC4 ( 50, COFX )
CALL ZVC4 ( 50, COFY )
C
CALL ZVC4 ( 50, WX )
CALL ZVC4 ( 50, WY )
CALL ZVC4 ( 50, RMDA )
C
FC = 2.* 0.00007292* SIN( 35./180.* PAI )
C
GM2 = 0.0026
FRIC = 0.25
VC = 10.
ROW = 1000.
C
ROA = 1.
ALF = 0.5
NSR = 0
C
C ----------------------------
DO 300 I = 1, 100
C
ISOR( I) = 0
SORW( I) = 0.
SORQ( I) = 0.
C
300 CONTINUE
C ----------------------------
DO 400 I = 1, 200
C
IFXC( I) = 0
FXC ( I) = 0.
C
400 CONTINUE
C ----------------------------
DO 500 I = 1, 150
C
IFTP( I) = 0
IFOR( I) = 0
C
FA( I) = 0.
FT( I) = 0.
FR( I) = 0.
FB( I) = 0.
C
500 CONTINUE
C ----------------------------
DO 600 I = 1, 100
C
ISEA( I) = 0
C
600 CONTINUE
C ----------------------------
DO 700 I = 1, 400
C
ILND( I) = 0
DLND( I) = 0.
C
700 CONTINUE
C ----------------------------
C
CALL ZVC4 ( NNP, X )
C
CALL ZVC4 ( NNP, Y )
C
I = NEL* 4
C
CALL ZVC2 ( I, NODE )
C
CALL ZVC4 ( NEL, AT )
CALL ZVC4 ( NEL, DPH )
CALL ZVC4 ( NEL, CFX )
CALL ZVC4 ( NEL, CFY )
C
CALL ZVC2 ( NNP, NX1 )
CALL ZVC2 ( NEL, MX1 )
CALL ZVC2 ( NP12, NW1 )
CALL ZVC2 ( NL12, MT1 )
C
I = NEL* 3
C
CALL ZVC2 ( I, KTR )
C
C ----------------------------
DO 800 I = 1, 2
DO 800 J = 1, 40
DO 800 K = 1, 4
C
LNE( I,J,K) = 0
C
800 CONTINUE
C ----------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE INIB ( U, V, Z, UD, VD, ZD, F, FX,
1 FY, FZ, AT, BE, CE, H, DEP, Q,
2 DMT, AVC, WXX, WYY, PX, PY, NODE,
3 NCD, KCD, NBU, NBV, NBZ, N4)
C
INTEGER*2 NODE,KCD, NBU, NBV, NBZ
C
DIMENSION U ( 81), V ( 81), Z ( 81), UD( NNP),
1 VD ( NNP), ZD ( NNP), F ( 81), FX( 81),
2 FY ( 81), FZ ( 81), AT ( NEL), BE( NEL,3),
C
3 CE ( NEL,3), H ( NEL), DEP( NNP), Q ( NEL),
4 DMT( NNZ), AVC( NNP), WXX( NEL), WYY( NEL),
5 PX ( NEL), PY ( NEL), NODE( NEL,4),
C
6 NCD ( 2,NZ2), KCD( NNP), NBU( NNP), NBV( NNP),
7 NBZ ( NNP)
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2,
1 NT1, USTT, ZSTT, NPLT, NP3
C
C ----- To set Initial Value ---
C
C ----------------------------
DO 100 I = 1, NNP
C
U( I) = USTT
V( I) = USTT
Z( I) = ZSTT
C
100 CONTINUE
C ----------------------------
C
CALL BCVL ( DEP, U, V, Z, U, V, Z, NBU, NBV, NBZ, NNP, 0.)
C
CALL PCAL ( H, AT, NODE, BE, CE, U, V, Z, WXX, WYY, PX,
1 PY, NNP, NEL )
C
CALL QCAL ( AT, NODE, BE, CE, U, V, Q, NNP, NEL )
C
C ----- UD ( DUDT ) ---
C
CALL VEC1 ( AT, PX, NODE, F, NNP, NEL, 1. )
C
C ----------------------------
DO 200 I = 1, NNP
C
FX( I) = F( I)
200 CONTINUE
C ----------------------------
C
CALL FINT ( F, U, DMT, NCD, KCD, F, NNP, NNZ, NZ2, -1.,-1. )
C
CALL EXPL ( F, UD, AVC, NBU, NNP)
C
C ----- VD ( = DVDT ) ---
C
CALL VEC1( AT, PY, NODE, F, NNP, NEL, 1. )
C
C ----------------------------
DO 300 I = 1, NNP
C
FY( I) = F( I)
300 CONTINUE
C ----------------------------
C
CALL FINT ( F, V, DMT, NCD, KCD, F, NNP, NNZ, NZ2, -1., -1. )
C
CALL EXPL ( F, VD, AVC, NBV, NNP )
C
C ----- DZ ( = DZDT ) ---
C
CALL VEC1( AT, Q, NODE, F, NNP, NEL, 1. )
C
C ----------------------------
DO 400 I = 1, NNP
C
FZ( I) = F( I)
400 CONTINUE
C ----------------------------
C
CALL FINT ( F, Z, DMT, NCD, KCD, F, NNP, NNZ, NZ2, -1., 0. )
C
CALL EXPL ( F, ZD, AVC, NBZ, NNP )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE INPT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1,
1 MX1, KTR, NW1, MT1)
C
INTEGER*2 NODE, NX1, NW1, MX1, MT1,KTR
C
CHARACTER*5 HEAD, EOF
CHARACTER*4 DTC, CID, RPNM, BLN
C
COMMON /TEST/ KEY
COMMON /TITE/ TITL1( 10)
COMMON /RPC2/ RPNM( 50)
COMMON /SECM/ NSA,ISEA( 100)
COMMON /FXCM/ NFXC,IFXC(200),FXC(200)
COMMON /LDCM/ NLND,ILND( 400),DLND( 400)
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /SRCM/ NSR,ISOR( 100), SORW( 100), SORQ( 100)
COMMON /RPC1/ NRP, IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
C
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
COMMON /FRCM/ NFR,IFTP( 150),IFOR( 150),FA( 150), FT( 150),
1 FR( 150),FB( 150)
COMMON /CMPF/ MFOR, JFTP( 150), JFOR( 150), PA( 150), PT( 150),
1 PR( 150), PB( 150)
C
DIMENSION X ( NNP), Y ( NNP), NODE ( NEL,4),
1 AT ( NEL), DPH( NEL), CFX ( NEL),
2 CFY ( NEL), NX1( NNP), NW1 ( NP12),
3 MX1 ( NEL), MT1( NL12), KTR ( NEL,3),
4 LNE ( 2,40,4), ICRD(10), DTC( 19 ),
5 HEAD( 25 )
C
C ( 1) ( 2) ( 3) ( 4) ( 5)
DATA HEAD /'AREA' , 'TITL' , 'TIDE' , 'DIFF' , 'NODE' ,
1 'ELEM' , 'REPP' , 'MDEP' , 'DEPT' , 'MATE' ,
2 'CONS' , 'SOUR' ,
3 'FXC ' , 'FORC' , 'SEA ' , 'LAND' , 'END ' ,
4 'CNST' , 'SOR ' , 'MTR ' ,
5 'JTN ' , 'ELM ' , 'CONT' , ' ' , ' ' /
C
DATA BLN /' ' /
DATA EOF /'EOF' /
DATA PAI /3.14159/
C
ICNTP = 0
ICNTE = 0
ICNTD = 0
C
MFOR = 0
ADDH = 0.
C
DO 111 I = 1, 10
C
ICRD( I) = 0
C
111 CONTINUE
C
CALL INIT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1, NW1, MX1,
1 MT1, KTR, LNE )
C
100 CONTINUE
C
READ(12,1000) CID
1000 FORMAT(A4)
BACKSPACE 12
READ(12,2000, END = 990) DTC
2000 FORMAT( 19A4)
C
ICNTD = ICNTD + 1
C
C ------------------------------------------------
DO 110 I = 1, 30
C
IF ( CID.EQ.HEAD( I)) GO TO 130
110 CONTINUE
C ------------------------------------------------
STOP
C
130 CONTINUE
C
BACKSPACE 12
C
C AREA TITL TIDE DIFF NODE ELEM REPP MDEP DEPT MATE CONS
C SOUR FXC FORC SEA LAND END CNST SOR MTR JTN ELM
C CONT
C
GO TO ( 200, 220, 240, 260, 320, 340, 360, 380, 400, 420, 440,
1 460, 500, 520, 540, 560, 990, 440, 460, 420, 320, 340,
2 620, 200, 200), I
C
200 READ(12,1100) DTC
1100 FORMAT(19A4)
GO TO 100
C
220 READ(12,1200) TITL1
1200 FORMAT( 10X,10A4)
GO TO 100
C
240 CONTINUE
C
C ----------------------------
C
CALL RTDE
C
C ----------------------------
C
GO TO 100
C
260 READ(12,1400) KYP, CB, CFL
1400 FORMAT( 10X,I5,6F5.0)
C
NFAL = 0
C
C ------------------------------------------------
DO 270 I = 1,5
C
IF ( CFL( I).GE.0.) NFAL = NFAL + 1
270 CONTINUE
C ------------------------------------------------
C
MDF = 0
GO TO 100
C
320 CONTINUE
C
C ------------------------------------------------
C
CALL RNDE ( NW1, NX1, X, Y, ICRD )
C
C ------------------------------------------------
GO TO 100
C
340 CONTINUE
C
C ----------------------------------------------------------
C
CALL RELM ( MT1, MX1, NW1, NODE, X, Y, AT, ICRD )
C
C ----------------------------------------------------------
GO TO 100
C
360 CONTINUE
C
C --------------------------------------
C
CALL RRPP ( NW1)
C
C --------------------------------------
GO TO 100
C
380 CONTINUE
C
C --------------------------------------
C
CALL RMDEP ( DPH, ADDH )
C
C --------------------------------------
GO TO 100
C
400 CONTINUE
C
C --------------------------------------
C
CALL RDPT ( MT1, DPH, ADDH )
C
C --------------------------------------
GO TO 100
C
420 CONTINUE
C
C ----------------------------
C
CALL RMTE
C
C ----------------------------
GO TO 100
C
440 CONTINUE
C
CALL RCONS ( PAI )
GO TO 100
C
460 CONTINUE
C
C ----------------------------
C
CALL RSOR ( NW1)
C
C ----------------------------
GO TO 100
C
500 CONTINUE
C
C ----------------------------
C
CALL RFXC ( NW1)
C
C ----------------------------
GO TO 100
C
520 CONTINUE
C
C -----------------------------------
C
CALL RFORC ( NX1, NW1,LNE )
C
C -----------------------------------
GO TO 100
C
540 CONTINUE
C
C -----------------------------------
C
CALL RSEA ( NW1, LNE )
C
C -----------------------------------
GO TO 100
C
560 CONTINUE
C
C -----------------------------------
C
CALL RLND ( NW1, LNE )
C
C ------------------------------------
GO TO 100
C
620 CONTINUE
C
C ----------------------------
C
CALL RCONT ( BLN)
C
C ----------------------------
GO TO 100
C
990 CONTINUE
C
IF ( ICRD(1).LT.NNP) THEN
WRITE(6,*) 'ERR,***INPT 2***'
STOP
END IF
C
IF ( ICRD( 2).LT.NEL) THEN
WRITE(6,*) 'ERR,***INPT 3***'
STOP
END IF
C
C ----- 'NODE' DATA relocation ---
C
L = ICRD ( 2)
M = 0
C
C --------------------------------------
DO 700 J = 1, 4
DO 700 I = 1, L
C
M = M + 1
NODE( M,1) = NODE( I,J)
C
700 CONTINUE
C --------------------------------------
C
NNP = ICRD( 1)
NEL = ICRD( 2)
DTC( 1) = EOF
C
NLND = NLND + 1
C
C --------------------------------------------------------------------
C
CALL CALC ( X, Y, NODE, CFX, CFY, KTR, LNE )
C
C --------------------------------------------------------------------
C
CALL EPRT ( X, Y, NODE, AT, DPH, CFX, CFY, NX1, NW1,
1 MX1, MT1, LNE )
C
C --------------------------------------------------------------------
C
C ----- Plotting Data ---
C
WRITE(6,2100) NSA, NLND
2100 FORMAT(/,' ** NSA, NLND =',2I5)
WRITE(11,*) NNP, NEL
C
C ------------------------------------------------
C
DO 800 I = 1, NNP
C
WRITE(11,*) X( I), Y( I)
800 CONTINUE
C ------------------------------------------------
DO 900 I = 1, NEL
C
WRITE(11,*) ( NODE( I,J), J = 1, 3)
900 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE INT1 ( NNP, CT, NBV )
C
C Boundary condution
C
INTEGER*2 NBV( NNP)
C
DIMENSION CT( NNP)
C
COMMON /FXCM/ NFXC, IFXC( 200), FXC( 200)
C
CALL ZVC2 ( NNP, NBV )
C
CALL ZVC4 ( NNP, CT )
C
DO 100 I = 1, NFXC
C
J = IFXC( I)
IF ( J.LE.0) GO TO 100
NBV( J) = 9999
CT ( J) = FXC( I)
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE KEMT ( NNP, NEL, NODE, NEU, NMU, NBV, LCL, UAV,
1 VAV, CFX, CFY, H, AT, X, Y, FU,
2 DPM, FVC, CT )
C
C ----- Global Matrix 'DPM' Construction ---
C
DIMENSION X ( NNP), Y( NNP), CFX( NEL), CFY( NEL),
1 UAV( NEL), VAV( NEL), FVC( NNP), CT ( NNP),
2 H ( NEL), AT( NEL), DPM( 81, 81), FU ( NNP),
3 LL ( 3)
C
REAL*8 XL( 3), YL( 3), B( 3), C( 3), COFX, COFY, UE, VE
INTEGER*2 NEU( NNP), NMU( NEL), NODE( NEL,4), LCL( NNP),
1 NBV( NNP)
C
CALL ZARY ( NNP, NNP, DPM )
C
C ----------------------------
DO 100 I = 1, NNP
C
NEU( I) = I
LCL( I) = I
C
100 CONTINUE
C ----------------------------
DO 150 J = 1, NEL
C
NMU( J) = J
150 CONTINUE
C ----------------------------
DO 200 K= 1, NEL
C
NE = NMU ( K)
C
C -----------------------
DO 210 I = 1, 3
C
XL( I) = X( NODE( NE, I ))
YL( I) = Y( NODE( NE, I ))
C
210 CONTINUE
C -----------------------
C
B( 1) = YL( 2) - YL( 3)
B( 2) = YL( 3) - YL( 1)
B( 3) = YL( 1) - YL( 2)
C
C( 1) = XL( 3) - XL( 2)
C( 2) = XL( 1) - XL( 3)
C( 3) = XL( 2) - XL( 1)
C
LL( 1) = LCL( NODE( NE,1))
LL( 2) = LCL( NODE( NE,2 ))
LL( 3) = LCL( NODE( NE,3))
C
COFX = CFX( NE)/ (4.* AT( NE))* H( NE )
COFY = CFY( NE)/ (4.* AT( NE))* H( NE )
C
UE = H( NE)* UAV( NE)/6.
VE = H( NE)* VAV( NE)/6.
C
C --------------------------------------------------------------------
DO 220 I = 1, 3
DO 220 J = 1, 3
C
II = LL( I)
JJ = LL( J)
C
DPM ( II,JJ) = DPM ( II,JJ)
1 + COFX* B( I)* B( J) + COFY* C( I)* C( J)
2 + UE* B( J) + VE* C( J)
C
220 CONTINUE
C --------------------------------------------------------------------
C
200 CONTINUE
C
C ----- Vector & Boundary condition ---
C
DO 300 I = 1, NNP
C
FU( I) = FVC ( NEU ( I))
300 CONTINUE
C -------------------------------------
C
DO 400 I = 1, NNP
C
IF ( NBV( NEU ( I)).NE.9999) GO TO 400
DPM( I, I ) = DPM( I,I )* 1.E10
FU ( I) = CT( NEU( I))* DPM( I,I )
C
400 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE LDVL ( U, V, NNP)
C
C ---- LAND BOUNDARY - 'U','V'- reCNSTruction
C
COMMON /LDCM/ NLND, ILND( 400), DLND( 400)
C
DIMENSION U( 81), V( 81)
C
DATA PIE /3.141592/
C
C ----------------------------------------------------------
DO 100 I = 1, NLND
C
IF ( ILND( I).LE.0) GO TO 100
J = ILND( I)
T = DLND( I) + PIE * 0.5
SINT = SIN( T)
COST = COS( T)
C
IF ( ABS( SINT) - ABS( COST)) 20, 20, 10
C
C ----- 'V' from 'U' ---
C
10 V( J) = - U( J)* COST / SINT
GOTO 200
C
C ----- 'U' from 'V' ---
C
20 U( J) = - V( J)* SINT/ COST
C
200 CONTINUE
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE MTX1 ( X, Y, NODE, AT, BE, CE, PMT, DMT,
1 WYY, NCD2, KCD, NNP, NEL, NNZ, NZ2 )
C
REAL*8 XL, YL, B, C, ELMT, AT4
INTEGER*2 KCD( NNP), NODE( NEL,4), NCD2( 2, NZ2 )
DIMENSION X ( NNP), Y ( NNP), AT ( NEL), BE ( NEL,3),
1 CE( NEL,3), PMT( NNZ), DMT( NNZ), XL ( 3),
2 YL( 3), B ( 3), C ( 3), WYY( NEL)
C
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
C
CALL ZVC4 ( NNZ, PMT )
C
CALL ZVC4 ( NNZ, DMT )
C
CALL ZVC2 ( 2* NZ2, NCD2 )
C
CALL ZVC2 ( NNP, KCD )
C
KKK1 = 0
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
DO 100 K1 = 1, 3
C ------------------------------------------------
C
IG = NODE ( I, K1)
C
C --------------------------------------
DO 170 K2 = 1, 3
C
JG = NODE( I, K2 )
ELMT = 1./12.
IF ( K1.EQ.K2 ) ELMT = 1./ 6.
C AV1 = AT( I)* ELMT
AV1 = ELMT
C
CALL DCNL ( AV1, IG, JG, PMT, NCD2, KKK1, NNZ )
C
170 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NEL
C
C --------------------------------------
DO 250 K = 1, 3
C
XL( K) = X ( NODE( I, K))
YL( K) = Y ( NODE( I, K))
C
250 CONTINUE
C --------------------------------------
C
B( 1) = YL( 2) - YL( 3)
B( 2) = YL( 3) - YL( 1)
B( 3) = YL( 1) - YL( 2)
C( 1) = XL( 3) - XL( 2)
C( 2) = XL( 1) - XL( 3)
C( 3) = XL( 2) - XL( 1)
C
C --------------------------------------
DO 270 L = 1, 3
C
BE( I, L) = B( L)
CE( I, L) = C( L)
C
270 CONTINUE
C --------------------------------------
C
AT4 = 0.25/ AT( I)
C
C --------------------------------------
DO 290 K1 = 1, 3
C
IG = NODE( I, K1)
C --------------------------------------
DO 290 K2 = 1, 3
C
JG = NODE( I,K2)
AV2 = AT4* ( B( K1)* B( K2) + C( K1)* C( K2))* WYY( I)
AV2 = AV2/ AT( I)
C
IF ( ABS( AV2).LE.0.0001) GO TO 333
C
CALL DCNL ( AV2, IG, JG, DMT, NCD2, KKK1, NNZ)
C
333 CONTINUE
C
290 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
CALL DCNIJ ( KCD, NCD2, KKK1, NNP, NZ2 )
C
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE OUT3 ( CT, CFX, CFY, NX1, CFL, CB, NNP, NEL )
C
C ----- NODAL VALUE OUTPUT ---
C
INTEGER*2 NX1( NNP)
DIMENSION CT( NNP), CFX( NEL), CFY( NEL)
C
DO 100 I = 1, NNP
C
CT( I) = CT( I) + CB
100 CONTINUE
C
WRITE(6,2000) CFL, ( NX1( I), CT( I),I = 1, NNP )
2000 FORMAT(' *** NODE CONCENTRATION *** CFL =',F7.3/
1 ( 5X,5( I5,F10.4,5X )) )
WRITE(6,*) ' ** CFX *', ( CFX( I),I = 1, NEL)
WRITE(6,*) ' ** CFY *', ( CFY( I),I = 1, NEL)
WRITE(30,*) ( CT ( I), I = 1, NNP)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PCAL ( H, AT, NODE, BE, CE, U, V, Z, WXX,
1 WYY, PX, PY, NNP, NEL)
C
C ----- PX, PY ---
C
INTEGER*2 NODE( NEL,4)
DIMENSION H ( NEL), AT( NEL), BE ( NEL,3), CE ( NEL,3), U ( 81),
1 V ( 81), Z ( 81), WXX( NEL), WYY( NEL), PX( NEL),
2 PY( NEL)
C
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
C
DATA G / 9.81274/
C
DO 100 I = 1, NEL
C
N1 = NODE( I,1)
N2 = NODE( I,2 )
N3 = NODE( I,3)
C
ZE = ( Z( N1) + Z( N2 ) + Z( N3))/ 3.
C
UEL = ( U( N1) + U( N2 ) + U( N3))/ 3.
VEL = ( V( N1) + V( N2 ) + V( N3))/ 3.
C
UES = ( U( N1) + U( N2 ) + U( N3))/ ( 3.* ( H( I) + ZE ))
VES = ( V( N1) + V( N2 ) + V( N3))/ ( 3.* ( H( I) + ZE ))
C
FX = GM2* UES* SQRT( UES** 2 + VES** 2 )
FY = GM2* VES* SQRT( UES** 2 + VES** 2 )
C
DUX = DVL1( U, BE, AT, NODE, NNP, NEL, I )
DUY = DVL1( U, CE, AT, NODE, NNP, NEL, I )
DVX = DVL1( V,BE,AT, NODE, NNP, NEL,I )
DVY = DVL1( V,CE,AT, NODE, NNP, NEL,I )
DZX = DVL1( Z,BE,AT, NODE, NNP, NEL,I )
DZY = DVL1( Z,CE,AT, NODE, NNP, NEL,I )
C
PX( I) = - G* ( H( I) + ZE )* DZX - UES* DUX
1 - VES* DUY + FC* VEL - FX + WXX( I)
C
PY( I) = - G* ( H( I) + ZE )* DZY - UES* DVX
1 - VES* DVY - FC* UEL - FY + WYY( I)
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PRSS ( MB, MW, NNP, A, AA, AAA, FU, FFU )
C
DIMENSION A ( 81,81), AA( 81,81), AAA( 81,33), FU( NNP),
1 FFU( NNP)
C
MW = 2* MB - 1
C
CALL ZARY ( NNP, NNP, AA )
C
DO 50 I = 1, NNP
DO 50 J = 1, MW
C
AAA( I,J) = 0.
50 CONTINUE
C
DO 100 I = 1, MB
DO 100 J = 1, NNP
C
AA( I,J) = A( I,J)
100 CONTINUE
C
JJ = 1
C
DO 200 I = MB + 1, NNP
DO 300 J = 1, NNP
C
IF ( J - JJ.LE.0) GO TO 250
AA( I,J - JJ) = A( I,J)
C
250 CONTINUE
IF( J.EQ.NNP) JJ = JJ+1
C
300 CONTINUE
200 CONTINUE
C
DO 400 I = 1, NNP
DO 400 J = 1, MW
C
AAA( I,J) = AA( I,J)
FFU( I) = FU( I)
C
400 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE QCAL ( AT, NODE, BE, CE, U, V, Q, NNP, NEL )
C
C ----- Q - Comp.---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION AT( NEL), BE( NEL,3), CE( NEL,3), U( 81), V( 81),
1 Q ( NEL)
C
DO 100 I = 1, NEL
C
Q( I) = - DVL1( U, BE, AT, NODE, NNP, NEL, I )
1 - DVL1( V, CE, AT, NODE, NNP, NEL, I )
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE QCLC ( NEL, NODE, H, AT, Q )
C
C ----- Source quantity ---
C
DIMENSION H( NEL), AT( NEL), Q( NEL)
C
INTEGER*2 NODE( NEL,4), NDFE( 20)
C
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
C
CALL ZVC4 ( NEL,Q )
C
IF ( NSR.LE.0) RETURN
C
T = 1./( 24.*3600. )
C ------------------------------------------------
DO 100 I = 1, NSR
C
NO = ISOR ( I)
N = 0
C
C --------------------------------------
DO 120 J = 1, NEL
DO 120 K = 1, 3
C
NG = NODE( J, K)
IF ( NG.NE.NO ) GO TO 120
C
N = N + 1
IF ( N.GT.20) THEN
WRITE(6,*) 'ERR,***QCLC***'
STOP
END IF
NDFE( N) = J
C
120 CONTINUE
C --------------------------------------
C
IF ( N.LE.0) GO TO 100
C
C --------------------------------------
DO 140 J = 1, N
C
NE = NDFE( J)
SW = AT( NE)* H( NE) + SORW( I)* T/ N
QW = SORQ( I)* T / N
C
C ----- UNIT : PPM ---
C
C Q( NE) = Q( NE) + QW / SW* 1.0E6
C
140 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RCONS ( PAI )
C
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
C
DIMENSION WK1D ( 20)
C
READ(12,1000) ( WK1D( I), I = 1, 9 )
1000 FORMAT(10X,9F5.0)
C
IF ( ABS( WK1D( 1)).GE.0.0001)
1 FC = 2.* 7.292E-5* SIN( WK1D( 1)/180.* PAI)
C
IF ( ABS( WK1D( 2)).GE.0.0001) GM2 = WK1D( 2)
IF ( ABS( WK1D( 3)).GE.0.0001) FRIC = WK1D( 3)* 0.001
IF ( ABS( WK1D( 4)).GE.0.0001) VC = WK1D( 4)
IF ( ABS( WK1D( 5)).GE.0.0001) ROW = WK1D( 5)
C
IF ( ABS( WK1D( 6)).GE.0.0001) ROA = WK1D( 6)
IF ( ABS( WK1D( 8)).GE.0.0001) ALF = WK1D( 8)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE RCONT ( BLN )
C
DIMENSION WK1D(20), CONT1(24)
C
CHARACTER*4 BLN, WK1D, CONT1
C
NCNT= 0.
C
READ(12,1000) ( WK1D( I),I = 1,12),( WK1D( I),I=21, 32)
1000 FORMAT( 10X,12A4,10X,( 12( 1X,A4)))
C
DO 100 I = 1, 12
C
C ----- Attention ? ---
C
IF ( WK1D( I+20).EQ.BLN) GO TO 100
C
NCNT = NCNT + 1
CONT1( NCNT) = WK1D( I)
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RDPT ( MT1, DPH, ADDH )
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
INTEGER*2 MT1
C
DIMENSION WK1D( 20), MT1( NL12), DPH( NEL)
C
READ(12,1000) J, ( WK1D( I), I = 1,10)
1000 FORMAT(10X,I5,5X,10F5.0)
C
C ------------------------------------------------
DO 100 I = 1, 10
C
IF ( ABS( WK1D( I)).LE.0.0001) THEN
GO TO 150
ELSE IF( J.GT.NL12) THEN
WRITE(6,*) '* Error *** RDPT 1 ***'
STOP
END IF
C
K = MT1( J)
IF ( K.EQ.0) THEN
WRITE(6,*) ' ERR,***RDPT 2***'
ELSE IF ( ABS( DPH( K) - ADDH).GE.0.0001) THEN
WRITE(6,*) ' ERR,***RDPT 3***'
STOP
END IF
C
DPH( K) = WK1D( I) + ADDH
C
150 J = J + 1
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RELM ( MT1, MX1, NW1, NODE, X, Y, AT, ICRD )
C
INTEGER*2 NODE, NW1, MX1, MT1
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
C
DIMENSION IWK2D( 2,10), WK2D( 2,10), ICRD(10), MT1( NL12),
1 MX1 ( NEL), NW1 ( NP12), NODE( NEL,4), X ( NNP),
2 Y ( NNP), WK1D( 20), AT ( NEL)
C
READ(12,1000) (( IWK2D( I,J),J = 1,5),I = 1,2 )
1000 FORMAT(10X,2(5I5,5X))
C
CALL ZVC4 ( 20, WK2D )
C
C ----------------------------------------------------------
DO 100 I = 1, 2
C
IF ( IWK2D( I,1).LE.0) THEN
C
GO TO 100
C
ELSE IF ( IWK2D( I,1).GT.NL12) THEN
WRITE(6,*) 'ERR,***RELM 1***'
ELSE IF ( ICRD( 2).GE.NEL) THEN
WRITE(6,*) 'ERR,***RELM 2***'
ELSE IF ( MT1( IWK2D( I,1)).NE.0) THEN
WRITE(6,*) 'ERR,***RELM 3***'
STOP
C
END IF
C
ICRD( 2) = ICRD( 2) + 1
MX1( ICRD ( 2)) = IWK2D( I,1)
MT1( IWK2D( I,1)) = ICRD( 2)
C
C ------------------------------------------------
DO 150 J = 1, 3
C
K = IWK2D( I,J+1)
IF ( K.GT.NP12 ) THEN
WRITE(6,*) 'ERR,***RELM 4***'
STOP
END IF
C
K = NW1( K)
IF ( K.EQ.0) THEN
WRITE(6,*) 'ERR,***RELM 5***'
STOP
END IF
C
NODE( ICRD( 2),J) = K
WK2D( 1,J) = X( K)
WK2D( 2,J) = Y( K)
C
150 CONTINUE
C ------------------------------------------------
C
IF ( IWK2D( I,5).EQ.0) IWK2D( I,5) = 1
C
NODE( ICRD( 2),4) = IWK2D( I,5)
WK1D( 1) = WK2D(1,1)* ( WK2D( 2,2) - WK2D( 2,3))
WK1D( 2) = WK2D(1,2)* ( WK2D( 2,3) - WK2D( 2,1))
WK1D( 3) = WK2D(1,3)* ( WK2D( 2,1) - WK2D( 2,2))
C
AT( ICRD( 2)) = ( W K1D(1) + WK1D( 2) + WK1D( 3))/ 2.
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RFORC ( NX1, NW1,LNE )
C
COMMON /SECM/ NSA,ISEA( 100)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
COMMON /FRCM/ NFR, IFTP( 150), IFOR( 150), FA( 150), FT( 150),
1 FR( 150), FB( 150)
C
INTEGER*2 NX1, NW1
C
DIMENSION IWKD( 20), WK2D( 2,10), LNE( 2,40,4), NX1( NNP),
1 NW1 ( NP12)
C
READ(12,1000) L, M,( IWKD( I),I = 1, 2 ),
1 (( WK2D( I,J),I = 1,2),J = 1,4 )
1000 FORMAT(10X,4I5,8F5.0)
C
C ------------------------------------------------
DO 100 I = 1,40
C
IF ( LNE( 1, I, 1).EQ.L) GO TO 110
100 CONTINUE
C ------------------------------------------------
C
WRITE(6,*) '* Error *** RFORC 1 ***'
STOP
C
110 K1 = LNE( 1,I,3)
K2 = LNE( 1,I,4)
IF ( NX1( ISEA( K1)).NE.IWKD(1) .OR.
1NX1( ISEA( K2)).NE.IWKD( 2)) THEN
WRITE(6,*) 'ERR,*** RFORC 2 ***'
STOP
END IF
C
NFR = NFR + 1
C ------------------------------------------------
DO 200 I = 1, 2
C
J = K1
IF ( I.EQ.2) J = K2
IFTP( J) = M
IFOR( J) = NW1( IWKD( I))
FA ( J) = WK2D( I,1)
FT ( J) = WK2D( I,2)
FR ( J) = WK2D( I,3)
FB ( J) = WK2D( I,4)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RFXC( NW1 )
C
COMMON /FXCM/ NFXC, IFXC( 200), FXC( 200)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
INTEGER*2 NW1
C
DIMENSION IWKD( 20), WK1D( 20), NW1( NP12)
C
READ(12,1000) ( IWKD( I), WK1D( I), I = 1, 4 )
1000 FORMAT(10X,4(I5,F10.0))
C
C ------------------------------------------------
DO 100 I = 1, 4
C
IF ( IWKD( I).EQ.0) THEN
GO TO 100
END IF
C
K = IWKD( I)
IF ( K.GT.NP12) THEN
WRITE(6,*) 'ERR,***RFXC 1***'
STOP
END IF
C
K = NW1( K)
IF ( K.EQ.0) THEN
WRITE(6,*) 'ERR,***RFXC 2***'
ELSE IF( NFXC.GE.200) THEN
WRITE(6,*) 'ERR,***RFXC 3***'
STOP
END IF
C
NFXC = NFXC + 1
IFXC( NFXC) = K
FXC ( NFXC) = WK1D( I)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RLND ( NW1, LNE )
C
COMMON /LDCM/ NLND, ILND( 400), DLND( 400)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
INTEGER*2 NW1
C
DIMENSION IWKD( 20), LNE( 2,40,4), NW1( NP12)
C
READ(12,1000) L, ( IWKD( I),I = 1, 11 )
1000 FORMAT(10X,12I5)
C
C ------------------------------------------------
DO 100 I = 1, 40
C
IF ( LNE( 2,I,1).EQ.0) GO TO 120
IF ( LNE( 2,I,1).EQ.L) GO TO 140
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,*) '* ERR,***RLND 1***'
STOP
C
120 NLND = NLND + 1
C
IF ( NLND.GT.400)
1WRITE(6,2000)
2000 FORMAT(' *** NLND.GT.400 , STOP ***')
IF ( NLND.GT.400) STOP
LNE( 2,I,1) = L
LNE( 2,I,3) = NLND + 1
C
140 CONTINUE
C
C ------------------------------------------------
DO 200 J = 1, 11
C
K = IWKD( J)
IF ( K.EQ.0) GO TO 200
IF ( K.GT.NP12) THEN
WRITE(6,*) ' * Error *** RLND 2 ***'
STOP
END IF
C
K = NW1( K)
IF ( K.EQ.0) THEN
WRITE(6,*) ' * Error *** RLND 3 ***'
C
ELSE IF( NLND.GE.400) THEN
WRITE(6,*) ' * Error *** RLND 4 ***'
STOP
END IF
C
NLND = NLND + 1
LNE( 2,I,2) = LNE( 2,I,2) + 1
LNE( 2,I,4) = NLND
ILND( NLND) = K
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RMDEP ( DPH, ADDH )
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
DIMENSION DPH( NEL)
C
READ(12,1000) ADDH
1000 FORMAT(10X,F10.0)
C
C --------------------------------------
DO 100 I = 1, NEL
C
DPH( I) = DPH( I) + ADDH
100 CONTINUE
C --------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RMTE
C
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
C
DIMENSION WK1D( 20)
C
READ(12,1000) IWK1, ( WK1D( I),I = 1, 6)
1000 FORMAT( 10X,I5,5F10.0,F5.0)
C
IF ( IWK1.LT.1 .OR. IWK1.GT.50) THEN
WRITE(6,*) ' * Error *** RMTE 1 ***'
ELSE IF ( ABS( COFV( IWK1)).GE.0.0001) THEN
WRITE(6,*) ' * Error *** RMTE 2 ***'
ELSE IF ( NMTR.GE.50) THEN
WRITE(6,*) ' * Error *** RMTE 3 ***'
STOP
END IF
C
NMTR = NMTR + 1
C
COFV( IWK1) = WK1D( 1)
COFX( IWK1) = WK1D( 4)
COFY( IWK1) = WK1D( 5)
C
WX ( IWK1) = WK1D( 2)
WY ( IWK1) = WK1D( 3)
RMDA( IWK1) = WK1D( 6)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RNDE ( NW1, NX1, X, Y, ICRD )
C
COMMON /TEST/ KEY
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
INTEGER*2 NX1, NW1
C
DIMENSION IWKD( 20), WK2D( 2,10), ICRD( 10), X( NNP),
1 Y ( NNP), NW1 ( NP12), NX1( NNP)
C
IF ( KEY.NE.2) GO TO 322
C
READ(12,1000) ( IWKD( I), ( WK2D( I,J),J = 1, 2), I = 1, 2 )
1000 FORMAT(10X,2( I5,5X,2F10.0))
GO TO 324
C
322 CONTINUE
READ(12,1100) ( IWKD( I),( WK2D( I,J),J = 1,2),I = 1,2)
1100 FORMAT( 10X,2( I5,5X,2F10.0))
324 CONTINUE
C
C ------------------------------------------------
DO 100 I = 1, 2
C
IF( IWKD( I).LE.0) THEN
GO TO 100
C
ELSE IF ( IWKD( I).GT.NP12) THEN
WRITE(6,*) ' Error *** RNDE 1 ***'
C
ELSE IF ( ICRD(1).GE.NNP) THEN
WRITE(6,*) ' Error *** RNDE 2 ***'
C
ELSE IF ( NW1( IWKD( I)).NE.0) THEN
WRITE(6,*) ' Error *** RNDE 3 ***'
STOP
END IF
C
ICRD( 1) = ICRD( 1) + 1
NX1( ICRD( 1)) = IWKD( I)
NW1( IWKD( I)) = ICRD( 1)
C
X( ICRD( 1)) = WK2D( I,1)* SCX
Y( ICRD( 1)) = WK2D( I,2)* SCY
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RRPP ( NW1 )
C
INTEGER*2 NW1
C
CHARACTER*4 WK1D,RPNM
C
COMMON /RPC2/ RPNM( 50)
COMMON /RPC1/ NRP, IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
DIMENSION IWKD( 20), WK1D( 20), NW1( NP12)
C
READ(12,1000) ( IWKD( I), WK1D( I), I = 1, 6)
1000 FORMAT(10X,6( I5,A4,1X ))
C
C ------------------------------------------------
DO 100 I = 1, 6
C
K = IWKD( I)
IF ( K.EQ.0) GO TO 100
IF ( K.GT.NP12) THEN
WRITE(6,*) 'ERR,***RRPP 1***'
STOP
END IF
C
K = NW1( K)
IF ( K.EQ.0) THEN
WRITE(6,*) 'ERR,***RRPP 2***'
ELSE IF ( NRP.GE.50) THEN
WRITE(6,*) 'NRP=', NRP
WRITE(6,*) 'ERR,***RRPP 3***'
STOP
END IF
C
NRP = NRP + 1
IRP ( NRP) = K
RPNM( NRP) = WK1D( I)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RSEA ( NW1, LNE )
C
COMMON /SECM/ NSA, ISEA( 100)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
C
INTEGER*2 NW1
C
DIMENSION IWKD( 20), LNE( 2, 40, 4), NW1( NP12 )
C
READ(12,1000) L, ( IWKD( I),I = 1,12)
1000 FORMAT( 10X,13I5)
C
C ------------------------------------------------
DO 100 I = 1, 40
C
IF( LNE( 1,I,1).EQ.0) GO TO 120
IF( LNE( 1,I,1).EQ.L) GO TO 140
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,*) ' Error *** RSEA 1 ***'
STOP
C
120 CONTINUE
C
LNE( 1,I,1) = L
LNE( 1,I,3) = NSA + 1
C
140 CONTINUE
C
C ------------------------------------------------
DO 200 J = 1, 12
C
K = IWKD( J)
IF ( K.EQ.0) GO TO 200
IF ( K.GT.NP12) THEN
WRITE(6,*) 'ERR,***RSEA 2***'
STOP
END IF
C
K = NW1( K)
IF ( K.EQ.0) THEN
WRITE(6,*) ' Error *** RSEA 3***'
ELSE IF( NSA.GE.150) THEN
WRITE(6,*) 'ERR,*** RSEA 4***'
STOP
END IF
C
NSA = NSA + 1
LNE( 1,I,2) = LNE( 1,I,2) + 1
LNE( 1,I,4) = NSA
ISEA( NSA) = K
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RSOR ( NW1 )
C
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
INTEGER*2 NW1
C
DIMENSION NW1( NP12), IWKD( 20), WK2D( 2,10)
C
READ(12,1000) ( IWKD( I), ( WK2D( I,J),J = 1, 2 ),I = 1, 2 )
1000 FORMAT(10X,2( I5, 5X, 2F10.0))
C
C ------------------------------------------------
DO 100 I = 1, 2
C
IF ( IWKD( I).EQ.0) RETURN
IF ( IWKD( I).GT.NP12 ) THEN
WRITE(6,*) 'ERR,***RSOR 1***'
STOP
END IF
C
IWKD( I) = NW1( IWKD( I))
IF( IWKD( I).EQ.0) THEN
WRITE(6,*) ' Error *** RSOR 2***'
ELSE IF( NSR.GE.200) THEN
WRITE(6,*) ' Error *** RSOR 3***'
STOP
END IF
C
NSR = NSR + 1
ISOR( NSR) = IWKD( I)
SORW( NSR) = WK2D( I,1)
SORQ( NSR) = WK2D( I,2 )
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RTDE
C
COMMON /TEST/ KEY
COMMON /OUTM/ LP1, LP2, LP3, LT1
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
C
DIMENSION WK1D( 200)
C
READ(12,1000) ( WK1D( I),I= 3, 11), USTT, ZSTT, NPLT, KEY
1000 FORMAT(10X,9F5.0,2F5.0,2I1)
C
RDT = WK1D( 3)
C
IF ( KEY.EQ.1) GO TO 2000
NSTR = ( ( WK1D( 4) - 1)* 12 + WK1D( 5))* 3600
NED = ( ( WK1D( 6) - 1)* 12 + WK1D( 7))* 3600
GO TO 3000
C
C ----- TEST ---
C
2000 NSTR = ( WK1D( 4)-1)* 12 + WK1D( 5)
NED = ( WK1D( 6)-1)* 12 + WK1D( 7)
C
3000 CONTINUE
NP1 = WK1D( 8)
NP2 = WK1D( 9)
C
IF ( ABS( WK1D(11)).LE.0.0001) WK1D(11) = 900.
NP3 = WK1D(10)
NT1 = WK1D(11)
IF ( NSTR.LE.0) NSTR = 0
IF ( MTHD.LE.0) MTHD = 3
IF ( JTP.LE.0) JTP = 5
C
IF ( WK1D( 8).GT.0.) LP1 = WK1D( 8)/ RDT + 0.5
IF ( WK1D( 9).GT.0.) LP2 = WK1D( 9)/ RDT + 0.5
IF ( WK1D(10).GT.0.) LP3 = WK1D(10)/ RDT + 0.5
IF ( WK1D(11).GT.0.) LT1 = WK1D(11)/ RDT + 0.5
C
KTDE = 1
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RVCT ( B, P, D, CO, NCC, KCC, NNP, NNZ, NZ2, NTIM )
C
C ----- 'RVEC'- Vector ---
C
INTEGER*2 NCC, KCC
C
DIMENSION B( NNP), P( NNZ), D( NNZ), CO( NNP), KCC( NNP),
1 NCC( 2,475)
C
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
JS = 1
IF( I.GT.1) JS = KCC( I-1)+1
JL = KCC( I)
FW = 0.
C
C --------------------------------------
DO 150 J = JS, JL
C
KPR = NCC( 2,J)
KCL = NCC( 1,J)
BB = P( KPR) - NTIM* ( 1.- ALF )* D( KPR)
FW = FW + BB* CO( KCL)
C
150 CONTINUE
C --------------------------------------
C
B( I) = FW
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE SDIF ( X, Y, NODE, AT, DPH, CFX, CFY,
1 UE, VE, U, V, UM, VM, SXR,
2 SYR, NX1, MX1, Q, FVC, CT, NBV,
3 LCL, NEU, NMU, FU, DPM )
C
C ----- Stationary Diffusion Analysis ---
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,4),
1 AT ( NEL), DPH( NEL), CFX ( NEL),
2 CFY( NEL), UE ( NEL), VE ( NEL),
3 U ( NNP), V ( NNP), UM ( NNP),
4 VM ( NNP), SXR( NEL), SYR ( NEL),
5 NX1( NNP), Q ( NEL), FVC ( NNP),
6 CT ( NNP), NBV( NNP), MX1 ( NEL),
7 LCL( NNP), NEU( NNP), NMU ( NEL),
8 FU ( NNP), DPT( 81,81), DP2M( 81,33),
9 FFU( 81), DPM( 81,81)
C
INTEGER*2 NODE, NX1, MX1,LCL, NEU, NMU, NBV
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
KSR = 1
C
IF ( KYP.EQ.3) KSR = NSR
C
OPEN ( 20,FILE='USER')
C
CALL INT1 ( NNP, CT, NBV )
C
CALL CFAS ( U, V, UE, VE, UM, VM, SXR, SYR, NODE, MX1,
1 NNP, NEL )
C
C ------------------------------------------------
DO 100 J = 1,KSR
C
CALL SRMD ( KSR, J, NX1 )
C
CALL QCLC ( NEL, NODE, DPH, AT, Q )
C
CALL FMT1 ( AT, DPH, Q, NODE, FVC, NNP, NEL )
C
IF ( KYP.NE.9) GO TO 150
NFAL = 1
CFL( 1) = 0.
C
150 CONTINUE
C
C --------------------------------------
DO 200 I = 1, NFAL
C
CFLX = CFL( I)
C
CALL CLCF ( NEL, AT, NODE, SXR, SYR, MX1, CFX, CFY, CFLX, 4 )
C
CALL KEMT ( NNP, NEL, NODE, NEU, NMU, NBV, LCL, UE,
1 VE, CFX, CFY, DPH, AT, X, Y, FU,
2 DPM, FVC, CT )
C
MW = 2* MB - 1
C
CALL PRSS ( MB, MW, NNP, DPM, DPT, DP2M, FU, FFU )
C
CALL GSOL ( DP2M, FFU, NNP, MB, MW )
C
CALL OUT3 ( FFU, CFX, CFY, NX1, CFLX, CB, NNP, NEL )
C
200 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,2000)
2000 FORMAT(20X,'*** END ***',/)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE SRMD ( KSR, J, NX1 )
C
C ----- KIYORITSU DATA SETUP for each 'SOR' ---
C
INTEGER*2 NX1( 1), KKSR
C
COMMON /DFCM/ MDF, KYP, CB, NFAL, CFL( 5)
COMMON /SRCM/ NSR, ISOR( 100), SORW( 100), SORQ( 100)
C
IF ( KYP.NE.3) RETURN
C
NSR = 1
C
IF( J.EQ.1) GO TO 90
ISRX = ISOR( 1)
SRWX = SORW( 1)
SRQX = SORQ( 1)
C
KKSR = KSR
IF ( KSR.LE.2 ) KKSR = 2
C
C --------------------------------------
DO 100 I = 2, KKSR
C
ISOR( I-1) = ISOR( I)
SORW( I-1) = SORW( I)
SORQ( I-1) = SORQ( I)
C
100 CONTINUE
C --------------------------------------
C
ISOR( KSR) = ISRX
SORW( KSR) = SRWX
SORQ( KSR) = SRQX
C
90 CONTINUE
WRITE(6,2000) NX1( ISOR( 1))
2000 FORMAT(/' *** Point Source Analysis ***',5X,'SOR - Point =',I5 )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE STNN ( NNP, NEL, NODE, NC1, NC2 )
C
INTEGER*2 NODE
C
DIMENSION NODE( NEL,4), NC1( NNP), NC2( NNP,8)
C
C --------------------------------------
DO 100 I = 1, NNP
C
NC1( I) = 0
100 CONTINUE
C --------------------------------------
DO 200 J = 1, 8
DO 200 I = 1, NNP
C
NC2( I,J) = 0
200 CONTINUE
C --------------------------------------
DO 300 J = 1, 3
DO 300 I = 1, NEL
C
N = NODE( I,J)
NC1( N) = NC1( N) + 1
NC2( N,NC1( N)) = I
C
300 CONTINUE
C --------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TIDL ( X, Y, NODE, AT, H, WXX, WYY,
1 BE, CE, PMT, DMT, PX, PY, Q,
2 AVC, BVC, UD, VD, ZD, DEP, FX,
3 FY, FZ, NCD, KCD, NBU, NBV, NBZ,
4 NX1, MX1, FS2, F2, F3, NC1, NC2,
5 B, PKM, AC1, DL, F, U, V,
6 Z, UO, VO, ZO, UE, VE, ZE,
7 UA, VA, ZA )
C
INTEGER*2 NODE, KCD, NBU, NBV, NBZ, NX1, MX1
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,4), AT ( NEL),
1 H ( NEL), WXX( NEL), WYY ( NEL), BE ( NEL,3),
2 CE ( NEL,3), PMT( NNZ), DMT ( NNZ), PX ( NEL),
3 PY ( NEL), Q ( NEL), AVC ( NNP), BVC( NNP),
4 UD ( NNP), VD ( NNP), ZD ( NNP), DEP( NNP),
5 FX ( NNP), FY ( NNP), FZ ( NNP), NCD( 2,NZ2),
8 KCD( NNP), NBU( NNP), NBV ( NNP), NBZ( NNP),
9 NX1( NNP), MX1( NEL)
C
DIMENSION PKM( NNZ), AC1( NNZ), DL( NEL), FS2 ( NNP),
1 F2 ( NEL), F3 ( NEL), B ( NEL), NC1( NNP),
2 NC2( NNP,8)
C
DIMENSION F ( NNP), U ( NNP), V ( NNP), Z ( NNP),
1 UO( NNP), VO( NNP), ZO( NNP), UE( NEL),
2 VE( NEL), ZE( NEL), UA( NNP), VA( NNP),
3 ZA( NNP)
C
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3,
1 NP12, NL12
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2,
1 NT1, USTT, ZSTT, NPLT, NP3
COMMON /MRCM/ NMTR, COFV( 50), COFX( 50), COFY( 50), WX( 50),
1 WY( 50), RMDA( 50)
C
ALF = 0.
C
CALL ZVC4 ( NNP, U )
CALL ZVC4 ( NNP, V )
CALL ZVC4 ( NNP, Z )
C
CALL ZVC4 ( NNP, UO )
CALL ZVC4 ( NNP, VO )
CALL ZVC4 ( NNP, ZO )
C
CALL ZVC4 ( NEL, UE )
CALL ZVC4 ( NEL, VE )
CALL ZVC4 ( NEL, ZE )
C
CALL ZVC4 ( NNP, F )
C
CALL ZVC4 ( NNP, UA )
CALL ZVC4 ( NNP, VA )
CALL ZVC4 ( NNP, ZA )
C
CALL ZVC4 ( NNP, UD )
CALL ZVC4 ( NNP, VD )
CALL ZVC4 ( NNP, ZD )
C
CALL ZVC2 ( NNP, NBU )
CALL ZVC2 ( NNP, NBV )
CALL ZVC2 ( NNP, NBZ )
C
DT = RDT
C
IF ( NSTR.LE.0) GO TO 5
GO TO 999
C
C ----- Nodal Depth & Wind Vector ---
C
5 CONTINUE
C
CALL DEPW ( H, NODE, DEP, WXX, WYY, NEL, NNP )
C
C ----- Attention : COFV( NEL) --- WYY( NEL) ---
C
G = 9.8
C
C ---------------------------------
DO 100 I = 1, NEL
C
DL( I) = SQRT( AT( I))
100 CONTINUE
C ---------------------------------
C
CALL COEF ( NNP, NEL, NODE, G, H, DT, DL, FS2, F3, F2,
1 NC1, NC2, B )
C
C ------------------------------------------------
DO 200 NENO = 1, NEL
C
WYY( NENO) = F3( NENO )* G* H( NENO)* DT
200 CONTINUE
C ------------------------------------------------
C
CALL MTX1 ( X, Y, NODE, AT, BE, CE, PMT, DMT, WYY,
1 NCD, KCD, NNP, NEL, NNZ, NZ2 )
C
CALL ZVC4 ( NEL, WYY )
C
C ----- First order --- Consistency of 'P'-Matrix ---
C
CALL ZVC4 ( NNP, AVC )
C
C --------------------------------------
DO 400 NL = 1, NEL
DO 400 NM = 1, 3
C
NK = NODE( NL,NM)
AVC( NK) = AVC( NK) + 1./ 3.
C
400 CONTINUE
C --------------------------------------
C
CALL INIB ( U, V, Z, UD, VD, ZD, F, FX, FY,
1 FZ, AT, BE, CE, H, DEP, Q, DMT, AVC,
2 WXX, WYY, PX, PY, NODE, NCD, KCD, NBU, NBV,
3 NBZ, 4 )
C
C ----- Matrix ( D.ALF.DT + P) or ( P) ---
C
DO 500 I = 1, NNZ
C
AC1( I) = 1.
500 CONTINUE
C
C2 = - ( 1.- ALF)* DT
C
CALL BMTX ( PMT, DMT, PKM, NNZ, AC1, C2 )
C
C ----- AVC --- For ( U), ( V) eq. : BVC --- for ( Z) eq.
C
C -----------------------------
DO 600 I = 1, NNP
C
UA( I) = U( I)
VA( I) = V( I)
ZA( I) = Z( I)
C
UO( I) = U( I)
VO( I) = V( I)
ZO( I) = Z( I)
C
600 CONTINUE
C ----------------------------
C
OPEN ( 20,FILE ='USER' )
C
C ----- Loop start ---
C
999 CONTINUE
C
CALL TLOP ( NODE, AT, H, WXX, WYY, BE, CE, PX,
1 PY, Q, AVC, FX, FY, FZ, PKM, DEP,
2 NCD, KCD, NBU, NBV, NBZ, NX1, MX1, FS2,
3 DT, IVR, F, U, V, Z, UO, VO,
4 ZO, UE, VE, ZE, UA, VA, ZA )
C
IF ( IVR.GT.0) THEN
WRITE(6,*) 'ERR,***TIDL***'
STOP
END IF
WRITE(6,2000)
2000 FORMAT(20X,'*** END ***')
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TLOP ( NODE, AT, H, WXX, WYY, BE, CE,
1 PX, PY, Q, AVC, FX, FY, FZ,
2 PKM, DEP, NCD, KCD, NBU, NBV, NBZ,
3 NX1, MX1, FS2, DT, IVR, F, U,
4 V, Z, UO, VO, ZO, UE, VE,
5 ZE, UA, VA, ZA )
C
INTEGER*2 NODE, KCD, NBU, NBV, NBZ, NX1, MX1
C
DIMENSION NODE( NEL,4), AT ( NEL), H ( NEL), WXX( NEL),
1 WYY ( NEL), BE ( NEL,3), CE ( NEL,3), PX ( NEL),
2 PY ( NEL), Q ( NEL), AVC( NNP), DEP( NNP),
3 NCD ( 2,NZ2), FX ( NNP), FY ( NNP), FZ ( NNP),
4 KCD ( NNP), NBU( NNP), NBV( NNP), NBZ( NNP),
5 NX1 ( NNP), MX1( NEL), PKM( NNZ), C1 ( 1200),
6 FS2 ( NNP)
C
COMMON /AREA/ NNP, NEL, MB, NNZ, NZ2, SCX, SCY, NLST, NL3, NP12,
1 NL12
C
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
DIMENSION F ( NNP), U ( NNP), V ( NNP), Z ( NNP), UO( NNP),
1 VO( NNP), ZO( NNP), UE( NEL), VE( NEL), ZE( NEL),
2 UA( NNP), VA( NNP), ZA( NNP)
C
LOP = 0
NTIM = NSTR
T = NTIM
C
900 CONTINUE
LOP = LOP + 1
T = T + 0.5* DT
C
NTIM = T
N = 0
C
C --------------------------------------
DO 100 I = 1, NNP
C
C1( I) = - FS2( I)* DT
100 CONTINUE
C --------------------------------------
C
C9 = 1.
C2 = 0.
GO TO 980
C
985 CONTINUE
T = T + 0.5* DT
N = N + 1
C
C --------------------------------------
DO 200 I = 1, NNP
C
C1( I) = - 0.2* DT
200 CONTINUE
C --------------------------------------
C
C2 = - 0.8* DT
C9 = 1.
980 CONTINUE
C
CALL BCVL ( DEP, U, V, Z, U, V, Z, NBU, NBV, NBZ, NNP, T )
C
CALL PCAL ( H, AT, NODE, BE, CE, UA, VA, ZA, WXX, WYY,
1 PX, PY, NNP, NEL )
C
C ----- 'U' CALC.---
C
CALL VEC1( AT, PX, NODE, F, NNP, NEL, 1. )
C
CALL BMTX ( FX, F, F, NNP, C1, C2 )
C
CALL FINT ( F, UO, PKM, NCD, KCD, F, NNP, NNZ, NZ2, 1., C9 )
C
CALL EXPL ( F, U, AVC, NBU, NNP)
C
C ----- 'V' CALC.---
C
CALL VEC1 ( AT, PY, NODE, F, NNP, NEL, 1. )
C
CALL BMTX ( FY, F, F, NNP, C1, C2 )
C
CALL FINT ( F, VO, PKM, NCD, KCD, F, NNP, NNZ, NZ2, 1., C9 )
C
CALL EXPL ( F, V, AVC, NBV, NNP )
C
CALL LDVL ( U, V, NNP)
C
C ----- 'Z' CALC.---
C
CALL QCAL ( AT, NODE, BE, CE, UA, VA, Q, NNP, NEL )
C
CALL VEC1 ( AT, Q, NODE, F, NNP, NEL, 1. )
C
CALL BMTX ( FZ, F, F, NNP, C1, C2 )
C
CALL FINT ( F, ZO, PKM, NCD, KCD, F, NNP, NNZ, NZ2, 1., C9 )
C
CALL EXPL ( F, Z, AVC, NBZ, NNP)
C
C --------------------------------------
DO 300 I = 1, NNP
C
UA( I) = U( I)
VA( I) = V( I)
ZA( I) = Z( I)
C
300 CONTINUE
C --------------------------------------
IF ( N.LT.1) GO TO 985
C
CALL PCAL ( H, AT, NODE, BE, CE, U, V, Z, WXX, WYY,
1 PX, PY, NNP, NEL)
C
CALL QCAL ( AT, NODE, BE, CE, U, V, Q, NNP, NEL)
C
CALL VEC1 ( AT, PX, NODE, FX, NNP, NEL, 1. )
C
CALL VEC1 ( AT, PY, NODE, FY, NNP, NEL, 1. )
C
CALL VEC1 ( AT, Q, NODE, FZ, NNP, NEL, 1. )
C
C ----- Output ---
C
CALL TOUT ( U, V, Z, UO, VO, ZO, UE, VE, ZE, NODE,
1 H, DEP, MX1, NX1, NNP, NEL, T, IVR, LOP )
C
IF ( IVR.GT.0) GO TO 999
IF ( NTIM - NED) 900, 999, 999
C
999 CONTINUE
C
C ----- Output for Plotting ---
C
WRITE( 10,*) LOP
C
C --------------------------------------
DO 400 I = 1, NEL
C
WRITE( 10,* ) UE( I), VE( I)
400 CONTINUE
C --------------------------------------
REWIND 10
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TOUT ( U, V, Z, UO, VO, ZO, UE, VE,
1 ZE, NODE, H, DEP, MX1, NX1, NNP, NEL,
2 T, IVR, LOP )
C
INTEGER*2 NODE( NEL,4), NX1( NNP), MX1( NEL)
C
CHARACTER*4 RPNM
C
DIMENSION U ( NNP), V( NNP), Z ( NNP), UO( NNP), VO( NNP),
1 ZO( NNP), H( NEL), DEP( NNP), UE( NEL), VE( NEL),
2 ZE( NEL)
C
COMMON /TTME/ NTM( 3)
COMMON /RPC2/ RPNM( 50)
COMMON /OUTM/ LP1,LP2,LP3, LT1
COMMON /RPC1/ NRP, IRP( 50), UR( 50), VR( 50), ZR( 50)
COMMON /CNST/ FC, GM2, FRIC, VC, ROW, ROA, ALF
COMMON /TDCM/ KTDE, MTHD, JTP, RDT, NSTR, NED, NP1, NP2, NT1,
1 USTT, ZSTT, NPLT, NP3
C
C --------------------------------------
DO 100 I = 1, NNP
C
UO( I) = U( I)
VO( I) = V( I)
ZO( I) = Z( I)
C
100 CONTINUE
C --------------------------------------
C
IVR = 0
VMAX = 0.
C
C ------------------------------------------------
DO 200 I = 1, NEL
C
ZW = 0.
UW = 0.
VW = 0.
C
C --------------------------------------
DO 250 J = 1, 3
C
ZW = ZW + Z( NODE( I,J))
UW = UW + U( NODE( I,J))
VW = VW + V( NODE( I,J))
C
250 CONTINUE
C --------------------------------------
C
ZE( I) = ZW/ 3.
UE( I) = UW/( 3.* ( H( I) + ZE( I)))
VE( I) = VW/( 3.* ( H( I) + ZE( I)))
C
VMAX = AMAX1( VMAX, ABS( UE( I)), ABS( VE( I)))
C
200 CONTINUE
C ------------------------------------------------
C
IF ( VMAX.GT.VC) IVR = 1
C
C ----- REP-POINT VALUE ---
C
IF ( NRP.LE.0) GO TO 350
C
C ------------------------------------------------
DO 300 I = 1, NRP
C
J = IRP( I)
ZR( I) = Z( J)
UR( I) = U( J)/ ( DEP( J) + ZR( I))
VR( I) = V( J)/ ( DEP( J) + ZR( I))
C
300 CONTINUE
C ------------------------------------------------
C
WRITE(9,*) LOP, IRP( 1), UR( 1), VR( 1), ZR( 1)
C
350 CONTINUE
NTIM = T + 0.5
NTM( 1) = T/ 43200.+ 1
NTM( 2) = AMOD( T,43200.)/ 3600.
NTM( 3) = AMOD( T,3600. )/ 60.
SEC = AMOD( T,60.)
C
C ----- To print all ppoints value ---
C
IF ( IVR.GT.0) GO TO 50
IF ( NP1.EQ.0) GO TO 90
IF ( NTIM - NED + 43200.LT.0) GO TO 90
IF ( MOD( LOP,LP1).NE.0) GO TO 90
C
50 CONTINUE
C
WRITE(6,*) ' '
WRITE(6,*) ' * Loop=', LOP
WRITE(6,2200) NTM, SEC, ( MX1( I), ZE( I), I = 1, NEL )
2200 FORMAT(' * Water level --- ZE Time =',
1 I3,' C ',I3,' H ',I3,' M ',F5.2,' S' //(3X,5( I4,F8.3,2X )))
C
90 CONTINUE
C
C ----- To print Rep. point value ---
C
IF ( NT1.EQ.0) GO TO 500
IF ( MOD( LOP, LT1).NE.0) GO TO 500
C
C ------------------------------------------------
DO 400 I = 1, NNP
C
U( I) = U( I)/ ( DEP( I) + Z( I))
V( I) = V( I)/ ( DEP( I) + Z( I))
C
400 CONTINUE
C ------------------------------------------------
C
500 CONTINUE
IF ( NP3.EQ.0) GO TO 900
IF ( MOD( LOP,LP3).NE.0) GO TO 900
C
900 CONTINUE
IF ( IVR.GT.0) WRITE(6,*) ' ** Error ** TOUT **'
IF ( IVR.GT.0) STOP
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TRAN ( NEL, NODE, KTR )
C
C ----- Re-Construction for Element Data ---
C
INTEGER*2 NODE, KTR
C
DIMENSION NODE( NEL,4), KTR( NEL,3), KZ5( 2,3)
C
DATA KZ5 / 1,2, 2,3, 3,1 /
C
C --------------------------------------
DO 100 I = 1, NEL
DO 100 J = 1, 3
C
KTR( I,J) = 9999
100 CONTINUE
C -----------------------------------------------
DO 200 I = 1, NEL
C
N = 0
C --------------------------------------
DO 300 J = 1, NEL
C
IF ( N.EQ.3) GO TO 200
IF ( I.EQ.J) GO TO 300
C
C ----------------------------
DO 400 K = 1, 3
DO 400 L = 1, 3
C
IF ( NODE( I,KZ5( 1,K)).NE.NODE( J, KZ5( 2,L))) GO TO 400
IF ( NODE( I,KZ5( 2,K)).NE.NODE( J, KZ5( 1,L))) GO TO 400
N = N + 1
KTR( I,K) = J
GO TO 300
C
400 CONTINUE
C ----------------------------
C
300 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,101)'C ************************************************ C'
WRITE(6,101)'C C'
WRITE(6,101)'C T I F D A S ( V.5.4 ) C'
WRITE(6,101)'C C'
WRITE(6,101)'C Tidal Flow & Diffusion Analysis C'
WRITE(6,101)'C C'
WRITE(6,101)'C by Finite Element Method C'
WRITE(6,101)'C C'
WRITE(6,101)'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,101)'C C'
WRITE(6,101)'C ************************************************ C'
WRITE(6,*) ' '
101 FORMAT(1X,A60)
C
RETURN
END
C ************************************************************************
C
SUBROUTINE VEC1 ( AT, Q, NODE, F, NNP, NEL, ALF )
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION AT( NEL), Q( NEL), F( NNP)
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
F( I) = 0.
100 CONTINUE
C ------------------------------------------------------
DO 200 I = 1, NEL
DO 200 J = 1, 3
C
F( NODE( I,J)) = F( NODE( I,J)) - Q( I)/3.* ALF
200 CONTINUE
C ------------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE ZARY ( NR, NC, A )
C
DIMENSION A( NR, NC )
C
NNR = NR
NNC = NC
C
IF( NR.LE.1) NNR = 1
IF( NC.LE.1) NNC = 1
C
C ---------------------------------
DO 100 J = 1, NNR
DO 100 I = 1, NNC
C
A( J,I ) = 0.
100 CONTINUE
C --------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE ZVC2 ( NV2, NV )
C
INTEGER*2 NV ( NV2 )
C
IF ( NV2.LE.1) NV2 = 1
C
DO 100 I = 1, NV2
C
NV( I) = 0
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE ZVC4 ( NV4, AV )
C
DIMENSION AV ( NV4)
C
IF ( NV4.LE.1) NV4 = 1
C
DO 100 I = 1, NV4
C
AV( I) = 0.
100 CONTINUE
C
RETURN
END
*************************************************************************