戻る

C ********************************************************************* C
C C
C W A V E - F E M C
C C
C Wave Analysis by Finite Element Method C
C C
C ( First-order Approx. ) C
C C
C ( V.3.5 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ********************************************************************* C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE
C
DIMENSION X( 39), Y( 39), NODE( 24,4), AT( 24), LNE( 2,5,4)
C
COMMON /AREA/ NNP, NEL, SCX, SCY
C
OPEN ( 3, FILE='WFEM_INP.DAT')
OPEN ( 8, FILE='WAVE_FEM_OUT')
OPEN ( 10, FILE='WAVE_FEM_PLDAT')
OPEN ( 11, FILE='WFEM_VAR' )
C
READ(3,1000) NNP, NEL, SCX, SCY
1000 FORMAT(10X,2I5,2F10.0)
C
CALL SINP ( X, Y, NODE, AT, LNE )
C
CALL CALC ( X, Y, NODE, AT, LNE )
C
CALL EPRT ( X, Y, NODE, AT, LNE )
C
CALL WAVS ( X, Y, NODE, AT )
C
CLOSE ( 3)
CLOSE ( 8)
CLOSE (10)
CLOSE (11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C WAVE - FDM C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( V.3.5 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Wave Analysis by Finite Difference Method C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WAVS ( X, Y, NODE, AT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE, NBU, NBV, NBZ
C
COMMON /AREA/ NNP, NEL, SCX, SCY
COMMON /WVCM/ DT, RED, LP1, NLOP
C
DIMENSION X ( NNP), Y ( NNP), NODE( NEL,3), AT ( NEL),
1 AE ( 24), BE ( 24), PMT ( 39,39), DMT( 39,39),
2 RMT( 39,39), AVC( 39), U ( 39), V ( 39),
3 Z ( 39), UO ( 39), VO ( 39), ZO ( 39),
4 UA ( 39), VA ( 39), ZA ( 39), FX ( 39),
5 FY ( 39), FZ ( 39), F ( 39), G1 ( 39),
6 G2 ( 39), NBU( 39), NBV ( 39), NBZ( 39)
C
WRITE(6,2000) NNP, NEL, SCX, SCY
2000 FORMAT(/,' * NNP =',I3,' NEL =',I3,' SCX =',F8.4,
1 ' SCY =', F8.4)
WRITE(6,*) ' '
WRITE(6,2100) ( X( I), I=1, NNP)
2100 FORMAT(' * X =',10F7.3)
WRITE(6,*) ' '
WRITE(6,2200) ( Y( I), I=1, NNP)
2200 FORMAT(' * Y =',10F7.3)
WRITE(6,*) ' '
C
CALL ZVC2 ( NNP, NBU )
C
CALL ZVC2 ( NNP, NBV )
C
CALL ZVC2 ( NNP, NBZ )
C
LOP = 0
T = 0.D0
C
G = 10.D0
H = 10.D0
F2 = - 0.398D0* DT
F3 = 0.036D0* G* H* DT
C
S2 = 0.D0
S3 = - 1.D0 * DT
C
CALL MTRX1 ( X, Y, NODE, AT, AE, BE, PMT, DMT, F3, NNP, NEL )
C
C ----- Concentration ---
C
CALL CNSS ( PMT, AVC, NNP)
C
C ----- Construction of Matrix ( D.ALPHA.DT + P ) or ( P ) ---
C
AC1 = 1.D0
AC2 = - DT
C
CALL BMTRX ( PMT, DMT, RMT, NNP, AC1, AC2 )
C
DO 100 I = 1, NNP
C
UA( I) = 0.D0
VA( I) = 0.D0
ZA( I) = 0.D0
C
UO( I) = 0.D0
VO( I) = 0.D0
ZO( I) = 0.D0
C
100 CONTINUE
C
990 LOP = LOP + 1
T = T + 0.5D0* DT
N = 0
C
C1 = F2
C2 = 0.D0
GO TO 980
C
985 T = T + 0.5D0* DT
N = N + 1
C1 = S2
C2 = S3
C
980 CONTINUE
C
CALL BCVAL ( U, V, Z, NBU, NBV, NBZ, T )
C
C ----- 'U' ---
C
CALL VECU ( F, ZA, AE, NNP, NODE, NEL )
C
CALL BVEC ( FX, F, G2, NNP, C1, C2 )
C
CALL FINT ( G2, UO, RMT, G1, NNP)
C
CALL EXSL ( G1, U, AVC, NBU, NNP)
C
C ----- 'V' CALC. ---
C
CALL VECV ( F, ZA, BE, NNP, NODE, NEL )
C
CALL BVEC ( FY, F, G2, NNP, C1, C2 )
C
CALL FINT ( G2, VO, RMT, G1, NNP)
C
CALL EXSL ( G1, V, AVC, NBV, NNP)
C
CALL UVRST ( U, V )
C
C ----- 'Z' CALC. ---
C
CALL VECZ ( F, UA, VA, AE, BE, NNP, NODE, NEL )
C
CALL BVEC ( FZ, F, G2, NNP, C1, C2 )
C
CALL FINT ( G2, ZO, RMT, G1, NNP)
C
CALL EXSL ( G1, Z, AVC, NBZ, NNP)
C
CALL COPY ( U, UA, NNP)
CALL COPY ( V, VA, NNP)
CALL COPY ( Z, ZA, NNP)
C
IF ( N.LT.1) GO TO 985
C
CALL VECU ( FX, Z, AE, NNP, NODE, NEL )
C
CALL VECV ( FY, Z, BE, NNP, NODE, NEL )
C
CALL VECZ ( FZ, U, V, AE, BE, NNP, NODE, NEL )
C
CALL TOUT ( U, V, Z, UO, VO, ZO, NODE, T, IOV, X, Y,
1 LOP, *990 )
C
WRITE(6,2300)
2300 FORMAT(/,'*** END ***')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BCVAL ( U, V, Z, NBU, NBV, NBZ, T )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NBU( NNP), NBV( NNP), NBZ( NNP)
C
COMMON /AREA/ NNP, NEL, SCX, SCY
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /FRCM/ NFOR, IFTP( 63), IFOR( 63)
C
DIMENSION U( NNP), V( NNP), Z( NNP)
C
C ------------------------------------------------
DO 300 I = 1, NFOR
C
J = IFOR( I)
C
IF ( J.EQ.0) GO TO 300
IF ( IFTP( I).EQ.1) GO TO 310
IF ( IFTP( I).EQ.2) GO TO 320
C
C ----- Z KYOSEI ---
C
310 CONTINUE
C
TEMP = DSIN( T)
C --------------------
C
IF ( NBZ( J).EQ.9999) GO TO 311
Z ( J) = TEMP
NBZ( J) = 9999
GO TO 300
C
311 Z ( J) = TEMP
GOTO 300
C
C ----- U KYOSEI ---
C
320 CONTINUE
TEMP = DSIN( T )
C
IF ( NBU( J).EQ.9999) GO TO 323
U ( J) = TEMP
NBU( J) = 9999
GO TO 300
C
323 U( J) = TEMP
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BMTRX ( PMT, DMT, A, NNP, C1, C2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- B - Left ---
C
DIMENSION PMT( NNP, NNP), DMT( NNP, NNP), A( NNP, NNP)
C
DO 100 I = 1, NNP
DO 100 J = 1, NNP
C
A( I,J) = C1* PMT( I,J) + C2* DMT( I,J)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BVEC ( FOLD, FNEW, F, NNP, C1, C2 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION FOLD( NNP), FNEW( NNP), F( NNP)
C
DO 100 I = 1, NNP
C
F( I) = C1* FOLD( I) + C2* FNEW( I)
C
100 CONTINUE
C
RETURN
END
C *********************************************************************
C
SUBROUTINE CALC ( X, Y, NODE, AT, LNE )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE
C
COMMON /AREA/ NNP, NEL, SCX, SCY
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /FRCM/ NFOR, IFTP( 63), IFOR( 63)
COMMON /FLCM/ NFLW, IFLW( 10)
COMMON /LDCM/ NBDR, IBDR( 100), DBDR( 100)
C
DIMENSION X( NNP), Y( NNP), NODE( NEL,4), AT( NEL),
1 LNE( 2,5,4), WK1D( 50)
C
DATA PAI / 3.14159 26535 89793 D0 /
C
CALL ZRV ( 50, WK1D )
C
C ----------------------------------------------------------
DO 100 I = 1, 5
C
IF ( LNE( 1,I,1).EQ.0) GO TO 990
C
J1 = LNE( 1,I,3)
J2 = J1 + 1
J4 = LNE( 1,I,4)
J3 = J4 - 1
WK9 = 0.D0
IF ( J3.LE.J1) J3 = J1
C
C ------------------------------------------------
DO 120 J = J1, J3
C
WK1 = ( X( IFLW( J+1)) - X( IFLW( J)))** 2
WK2 = ( Y( IFLW( J+1)) - Y( IFLW( J)))** 2
WK1D( J+1) = DSQRT ( WK1 + WK2)
WK9 = WK9 + WK1D( J+1)
C
120 CONTINUE
C ------------------------------------------------
C
IF ( J3.LE.J2) J3 = J2
C ------------------------------------------------
DO 140 J = J2, J3
C
IFTP( J) = IFTP( J1)
IFOR ( J) = IFLW ( J )
C
140 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
990 CONTINUE
C
NFOR = NFLW
C
CALL ZRV ( 50, WK1D )
C
C ----------------------------------------------------------
DO 160 I = 1, 5
C
IF ( LNE( 2,I,1).EQ.0) GO TO 333
J1 = LNE( 2,I,3)
J4 = LNE( 2,I,4)
C
J2 = J1 + 1
J3 = J4 - 1
ISW = 0
C
IF ( J3.LE.J1) J3 = J1
C ------------------------------------------------
DO 180 J = J1, J3
C
WKX = X( IBDR( J+1)) - X( IBDR( J))
WKY = Y( IBDR( J+1)) - Y( IBDR( J))
WK1D( J) = ATAN2( WKY, WKX)
C
180 CONTINUE
C ------------------------------------------------
C
IF ( LNE( 2,I,2).LE.2) GO TO 160
IF ( IBDR( J1).EQ.IBDR( J4)) GO TO 2060
C
DBDR( J1) = WK1D( J1)
DBDR( J4) = WK1D( J3)
C
IF ( DBDR( J1).LT.0.D0) DBDR( J1) = DBDR( J1) + PAI
IF ( DBDR( J4).LT.0.D0) DBDR( J4) = DBDR( J4) + PAI
IF ( J3.LE.J2) J3 = J2
C
C ------------------------------------------------
DO 185 J = J2, J3
C
DBDR( J) = ( WK1D( J-1) + WK1D( J))/ 2.
IF ( DBDR( J).LT.0. ) DBDR( J) = DBDR( J) + PAI
C
185 CONTINUE
C ------------------------------------------------
GO TO 160
C
2060 CONTINUE
C ------------------------------------------------
DO 190 J = J1, J3
C
J5 = J-1
IF ( J.EQ.J1) J5 = J3
DBDR( J) = ( WK1D( J5 ) + WK1D( J))/ 2.
IF ( DBDR( J).LT.0.D0) DBDR( J) = DBDR( J) + PAI
C
190 CONTINUE
C ------------------------------------------------
C
DBDR( J4) = DBDR( J1)
C
160 CONTINUE
C ----------------------------------------------------------
C
333 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CNSS ( A, B, NNP)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Concetration of Matrix ---
C
DIMENSION A( NNP, NNP), B( NNP)
C
C ---------------------------------
DO 100 I = 1, NNP
C
B( I) = 0.D0
C
100 CONTINUE
C ---------------------------------
DO 200 I = 1, NNP
DO 200 J = 1, NNP
C
B( I) = B( I) + A( I,J)
C
200 CONTINUE
C ---------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COPY ( A, B, NNP)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A( NNP), B( NNP)
C
DO 100 I = 1, NNP
C
B( I) = A( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
FUNCTION DVAL1 ( U, BE, AT, NODE, NNP, NEL, I )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Differenced Value ( 1st order ) ---
C
INTEGER*2 NODE( NEL,3)
C
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
DVAL1 = ( BE( I,1)* U1 + BE( I,2 )* U2 + BE( I,3)* U3)/ 2.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DVALA ( A, Z, AE, NODE, NNP, NEL, I )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- 'F': 1st order approx.---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION A( 4), Z( NNP), AE( NEL)
C
DO 100 J = 1, 4
C
A( J) = 0.D0
C
100 CONTINUE
C
Z1 = Z( NODE( I,1))
Z2 = Z( NODE( I,2))
Z3 = Z( NODE( I,3))
Z4 = Z( NODE( I,4))
C
C13 = 0.33333 33333 33333 33333 D0
C16 = 0.16666 66666 66666 66667 D0
C
A( 1) = AE( I)* ( C13* ( Z2 - Z1) + C16* ( Z3 - Z4))
A( 2) = A ( 1)
A( 3) = AE( I)* ( C16* ( Z2 - Z1) + C13* ( Z3 - Z4))
A( 4) = A ( 3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DVALB ( A, Z, BE, NODE, NNP, NEL, I )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- 'F': 1st order approx.---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION A( 4), Z( NNP), BE( NEL)
C
DO 100 J = 1, 4
C
A( J) = 0.D0
C
100 CONTINUE
C
Z1 = Z( NODE( I,1))
Z2 = Z( NODE( I,2))
Z3 = Z( NODE( I,3))
Z4 = Z( NODE( I,4))
C
C13 = 0.33333 33333 33333 33333 D0
C16 = 0.16666 66666 66666 66667 D0
C
A( 1) = BE( I)* ( C13* ( Z4 - Z1) + C16* ( Z3 - Z2 ))
A( 2) = A ( 1)
A( 3) = BE( I)* ( C16* ( Z4 - Z1) + C13* ( Z3 - Z2 ))
A( 4) = A ( 3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EPRT ( X, Y, NODE, AT, LNE )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE
C
COMMON /RPCM/ NRP, IRP( 20)
COMMON /FLCM/ NFLW, IFLW( 10)
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /LDCM/ NBDR, IBDR( 100), DBDR( 100)
COMMON /AREA/ NNP, NEL, SCX, SCY
C
DIMENSION X( NNP), Y( NNP), NODE( NEL,4), AT( NEL),
1 LNE( 2,5,4), IW1D( 20), WK2D( 2,10)
C
C ----- HEADING ---
C
WRITE(8,2000)
WRITE(8,2100)
WRITE(8,2200)
WRITE(8,2200)
WRITE(8,2100)
C
2000 FORMAT(' ',10X, 80('*'))
2100 FORMAT(' ',10X,'*',78X,'*')
2200 FORMAT(' ',10X,'*',15X,
1 '==== I N F O R M A T I O N L I S T ====',19X,'*')
C
C ----- AREA ---
C
WRITE(8,2300) NNP, NEL, SCX, SCY
2300 FORMAT(/,'----- Area Information -----'//,
1 9X,'NNP ( SETTEN SU) ...............',I7/
2 9X,'NEL ( YOUSO SU) ...............',I7/
3 9X,'SCX ( SETTEN SCALE KEISU X) ....',G14.6/
4 9X,'SCY ( SETTEN SCALE KEISU Y) ....',G14.6 )
C
C ----- TIDE ---
C
WRITE(8,2400) DT, RED, LP1
2400 FORMAT(/,'----- TIDE INFORMATION -----'//
1 9X,'TIME STEP .....................',F9.3,'( SEC)'/
2 9X,'END TIME ....................',F9.3,'( SEC)'/
3 9X,'LP1 ( ZEN-TEN PRINT KANKAKU) ...',I9, '( - )' )
C
C ----- Reppoint ---
C
IF ( NRP.EQ.0) GO TO 110
C
WRITE(8,2500)
C
I1 = ( NRP + 4)/ 5
C
DO 100 I = 1, I1
C
J2 = I* 5
J1 = J2 - 4
IF ( J2.GT.NRP ) J2 = NRP
WRITE(8, 951) J1, J2, ( IRP( J), J = J1, J2)
C
100 CONTINUE
C
2500 FORMAT(/,'*-------- REPPOINT INFORMATION ---------*'/
1 9X,'REP-#',5X,5(' NODE NAME '))
951 FORMAT(/,9X,I2,'-',I2,5X,5( I8,2X,10X))
C
C ----- FLW ---
C
110 IF ( NFLW.EQ.0) GO TO 120
C
WRITE(8,2600)
C
DO 200 I = 1, 5
C
IF ( LNE( 1,I,1).EQ.0) GO TO 120
J1 = LNE( 1,I,3)
J2 = LNE( 1,I,4)
WRITE(8,2700) LNE( 1,I,1), ( IFLW( J),J = J1, J2 )
C
200 CONTINUE
C
2600 FORMAT(/,'*-------- FLOW INFORMATION ---------*')
2700 FORMAT(/9X,'LNE',I3,1X,'NODE =',( 10I7/))
C
C ----- BDR ---
C
120 IF ( NBDR.EQ.0) GO TO 140
C
WRITE(8,2800)
C
DO 300 I = 1, 5
C
IF ( LNE( 2,I,1).EQ.0) GO TO 140
J1 = LNE( 2,I,3)
J2 = LNE( 2,I,4)
WRITE(8,2850) LNE( 2,I,1)
C
DO 350 JJ = J1, J2, 4
C
J3 = JJ + 4
IF ( J3.GT.J2) J3 = J2
WRITE(8,2900) ( IBDR( J), DBDR( J), J = JJ, J3)
C
350 CONTINUE
C
300 CONTINUE
C
2800 FORMAT(/,'*---------- Boundary Information ---------------*'//
1 9X,' LNE-#',3X,5('NODE KAKUDO '))
2850 FORMAT(/,9X,I6)
2900 FORMAT(' ',19X,5( I4,F10.6,6X))
C
C ----- OTHER ---
C
140 CONTINUE
WRITE(8,2950)
2950 FORMAT(/,'*------ NODE Information ----------*')
WRITE(8,2970)
2970 FORMAT(/,9X,4('NODE X (M) Y (M) ')//)
C
DO 400 I = 1, NNP, 4
C
I1 = I + 3
IF ( I1.GT.NNP) I1 = NNP
J1 = 0
C
DO 450 J = I, I1
C
J1 = J1 + 1
WK2D( 1,J1) = X( J)
WK2D( 2,J1) = Y( J)
IW1D( J1) = J
C
450 CONTINUE
C
WRITE(8,2980) ( IW1D( J),' ',WK2D( 1,J),WK2D( 2, J),J = 1,J1)
2980 FORMAT(9X,4( I4,1X,A1,2F10.3,4X))
C
400 CONTINUE
C
WRITE(8,2990)
WRITE(8,2995)
2990 FORMAT(/,'*----- ELEMENT INFORMATION --------*')
2995 FORMAT(/, 4X ,'LEMT# MENSEKI NODE-1 NODE-2 NODE-3'
1 ,'NODE-4'/)
C
DO 500 I = 1, NEL
DO 550 J = 1, 4
C
IW1D( J) = NODE( I,J)
IF ( J.LE.4) IW1D( J) = IW1D( J)
C
550 CONTINUE
C
WRITE(8,2997) I, AT( I), ( IW1D( J), J = 1, 4)
2997 FORMAT(4X,I3,F13.3,1X,4I8)
C
500 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXSL ( FVEC, C, TSM, NBC, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Explicit Formulation ---
C
INTEGER*2 NBC( NNP)
C
DIMENSION FVEC( NNP), C( NNP), TSM( NNP)
C
DO 100 I = 1, NNP
C
IF ( NBC( I).EQ.9999) GO TO 100
C( I) = FVEC( I)/ TSM( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FINT ( F, C, DMT, G, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Right hand side ---
C
DIMENSION F( NNP), C( NNP), DMT( NNP,NNP), G( NNP)
C
DO 100 I = 1, NNP
C
G( I) = 0.D0
C
100 CONTINUE
C
DO 200 I = 1,NNP
DO 200 J = 1,NNP
C
G( I) = G( I) + DMT( I,J)* C( J)
C
200 CONTINUE
C
DO 300 I = 1, NNP
C
G( I) = G( I) + F( I)
C
300 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE UVRST ( U, V )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- 'U','V' restriction ---
C
COMMON /AREA / NNP, NEL, SCX, SCY
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /LDCM/ NBDR, IBDR( 100), DBDR( 100)
C
DIMENSION U( NNP), V( NNP)
C
DATA PAI / 3.14159 26535 89793 D0 /
C
DO 100 I = 1, NBDR
C
IF ( IBDR( I).LE.0) GO TO 100
J = IBDR( I)
T = DBDR( I) + PAI* 0.5D0
C
SINT = DSIN( T )
COST = DCOS( T )
C
IF ( DABS( SINT ) - DABS( COST )) 20, 20, 10
C
C ----- ( 'V' from 'U' ) ---
C
10 V( J) = - U( J)* COST/ SINT
GOTO 333
C
C ----- ( 'U' from 'V' ) ---
C
20 U( J) = - V( J)* SINT/ COST
GOTO 333
C
333 CONTINUE
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
FUNCTION LIJ ( N, I,J)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LIJ = ( N-1)* ( I-1) + ( I* ( 3-I)-2)/ 2 + J
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MTRX1 ( X, Y, NODE, AT, AE, BE, PMT, DMT,
1 F3, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE ( NEL,4)
C
DIMENSION X ( NNP), Y ( NNP), AT ( NEL), AE( NEL),
1 BE( NEL), PMT( NNP,NNP), DMT( NNP,NNP)
C
DIMENSION XL( 4), YL( 4), AM( 10), AM0( 10)
C
KKK1 = 0
C
CALL ZARY ( NNP, NNP, PMT )
C
CALL ZARY ( NNP, NNP, DMT )
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
AM0( 1) = AT( I)/ 9.D0
AM0( 2) = AT( I)/ 18.D0
AM0( 3) = AT( I)/ 36.D0
AM0( 4) = AM0( 2 )
AM0( 5) = AM0( 1)
C
AM0( 6) = AM0( 2)
AM0( 7) = AM0( 3)
AM0( 8) = AM0( 1)
AM0( 9) = AM0( 2)
AM0(10) = AM0( 1)
C
C ------------------------------------------------
DO 120 K1 = 1, 4
C
IG = NODE( I, K1)
C
C --------------------------------------
DO 140 K2 = 1, 4
C
JG = NODE( I, K2 )
IF ( K1.GE.K2 ) GO TO 15
AV1 = AM0 ( LIJ( 4, K1, K2 ))
GO TO 16
15 CONTINUE
AV1 = AM0 ( LIJ( 4, K2, K1))
16 CONTINUE
PMT( IG,JG) = PMT( IG,JG) + AV1
C
140 CONTINUE
C --------------------------------------
C
120 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
C13 = 0.33333 33333 33333 33333 D0
C16 = 0.16666 66666 66666 66667 D0
C
C ----------------------------------------------------------
DO 900 I = 1, NEL
C
C ------------------------------------------------
DO 920 K = 1, 4
C
XL( K) = X( NODE( I, K ))
YL( K) = Y( NODE( I, K ))
C
920 CONTINUE
C ------------------------------------------------
C
AE( I) = ( YL( 3) - YL( 1))/ 2.D0
BE( I) = ( XL( 2 ) - XL( 1))/ 2.D0
A1 = AE( I)/ BE( I)
B1 = BE( I)/ AE( I)
C
AM( 1) = A1* ( C13) + B1* ( C13)
AM( 2) = A1* ( - C13) + B1* ( C16)
AM( 3) = A1* ( - C16) + B1* ( - C16)
AM( 4) = A1* ( C16) + B1* ( - C13)
AM( 5) = A1* ( C13) + B1* ( C13)
C
AM( 6) = A1* ( C16) + B1* ( - C13)
AM( 7) = A1* ( - C16) + B1* ( - C16)
AM( 8) = A1* ( C13) + B1* ( C13)
AM( 9) = A1* ( - C13) + B1* ( C16)
AM(10) = A1* ( C13) + B1* ( C13)
C
C ------------------------------------------------
DO 940 K1 = 1, 4
C
IG = NODE( I, K1)
C ------------------------------------------------
DO 940 K2 = 1, 4
C
JG = NODE( I,K2)
IF ( K1.GE.K2) GO TO 131
AV2 = AM ( LIJ( 4, K1, K2 ))* F3
GO TO 132
131 CONTINUE
AV2 = AM ( LIJ( 4, K2, K1)) * F3
132 CONTINUE
DMT( IG, JG ) = DMT( IG, JG ) + AV2
C
940 CONTINUE
C ------------------------------------------------
C
900 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CNER ( U, SEC, N, X, Y, T, LOP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C U : For Rep Point SEC : Print out Time
C N : Number of Rep-points NNP : Number of Points
C
COMMON /AREA/ NNP, NEL, SCX, SCY
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /RPCM/ NRP, IRP( 20)
C
DIMENSION U( 20), EXU( 20), X( NNP), Y( NNP)
C
COMMON / STER / TSMZ
C
SMZ = 0.D0
SUMS = 0.D0
C
C ----- Anal. Value ---
C
RRT1 = 10.D0* T
C
DO 100 I = 1, N
C
EXU( I) = DSIN( T - X( IRP( I))/ 10.D0)
IF ( X( IRP( I)).GT.RRT1) EXU( I) = 0.D0
C
100 CONTINUE
C
C ----- Relative Error ---
C
DO 200 I = 1, N
C
SMZ = SMZ + DABS( U( I) - EXU( I))
C
200 CONTINUE
C
SMZ = SMZ/ DFLOAT( N) ! Aver. estimation
SMZ = SMZ/ 2.D0 ! Max.: 1. x 2 Considering "+" & "-"
SMZ = SMZ* 100.D0 ! (%) ! Rel. Error
C
WRITE(8,2000) LOP, SEC, SMZ
2000 FORMAT(' * LOP =', I4,' Time =',F9.3,' Rel. Error=',F11.3,'(%)')
IF ( MOD( LOP,10).EQ.0) THEN
WRITE(6,2100) LOP, SEC, U( 1), SMZ
2100 FORMAT(' * LOP =', I4,' T =',F7.3,' * U( 1)=', F8.3,
1 ' Rel.Error =',F8.3,' (%)')
END IF
C
WRITE(10,*) LOP, SEC, U( 1), SMZ
TSMZ = TSMZ + SMZ
C
IF ( LOP.EQ.NLOP) THEN
TSMZ = TSMZ/ DFLOAT( LOP)
WRITE(6,2200) TSMZ, N
2200 FORMAT(/,' * Aver. Error =',F8.3,' (%) N =',I3)
C
WRITE(11,*) N, NNP
WRITE(11,*) ( U ( I), I=1, N )
WRITE(11,*) ( EXU( I), I=1, N )
WRITE(11,*) ( X( I), I=1, NNP)
WRITE(11,*) ( Y( I), I=1, NNP)
C
END IF
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCAL ( AT, NODE, BE, CE, Z, PX, PY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Px & Py ----
C
INTEGER*2 NODE( NEL,3)
C
DIMENSION AT( NEL), BE( NEL,3), CE( NEL,3), Z( NNP),
1 PX( NEL), PY( NEL)
C
G = 10.D0
H = 10.D0
C
DO 100 I = 1, NEL
C
DZX = DVAL1( Z, BE, AT, NODE, NNP, NEL, I )
DZY = DVAL1( Z, CE, AT, NODE, NNP, NEL, I )
C
PX( I) = - G* DZX
PY( I) = 0.D0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE QCAL ( AT, NODE, BE, CE, U, V, Q, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE( NEL,3)
C
DIMENSION AT( NEL), BE( NEL,3), CE( NEL,3), U( NNP), V( NNP),
1 Q ( NEL)
C
H = 10.D0
C
DO 100 I = 1, NEL
C
Q( I) = ( - DVAL1( U, BE, AT, NODE, NNP, NEL, I )
1 - DVAL1( V, CE, AT, NODE, NNP, NEL, I ))* H
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SINP ( X, Y, NODE, AT, LNE )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE
CHARACTER*4 HEAD, DT80
C
COMMON /STER/ TSMZ
COMMON /RPCM/ NRP, IRP( 20)
COMMON /FLCM/ NFLW, IFLW( 10)
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /FRCM/ NFOR, IFTP( 63), IFOR( 63)
COMMON /LDCM/ NBDR, IBDR( 100), DBDR( 100)
COMMON /AREA/ NNP, NEL, SCX, SCY
C
DIMENSION X ( NNP), Y ( NNP), AT( NEL), NODE( NEL,4),
1 LNE( 2,5,4), IW1D( 20), IW2D( 2,10),
2 WK1D( 20), WK2D( 2, 10)
C
DIMENSION DT80( 20), HEAD( 10)
C
DATA HEAD( 1) /'AREA'/, HEAD( 2)/'TITL'/, HEAD( 3)/'TIDE'/
DATA HEAD( 4) /'NODE'/, HEAD( 5)/'ELEM'/, HEAD( 6)/'REPP'/
DATA HEAD( 7) /'FORC'/, HEAD( 8)/'FLW '/, HEAD( 9)/'BDR'/
DATA HEAD(10) /'END '/
C
CALL TTLE
C
REWIND 3
C
C ----- Initial Setting ---
C
DO 100 I = 1, 2
DO 100 J = 1, 5
DO 100 K = 1, 4
C
LNE( I,J,K) = 0
C
100 CONTINUE
C
TSMZ = 0.D0
C
900 CONTINUE
READ (3,1000,END=990) ( DT80( J), J = 1, 20)
WRITE(8,1000) ( DT80( J),J = 1, 20)
1000 FORMAT(20A4)
C
DO 110 I = 1, 10
C
IF ( DT80( 1).EQ.HEAD( I)) GO TO 130
C
110 CONTINUE
C
130 CONTINUE
C
BACKSPACE 3
C
C ( 1) ( 2) ( 3) ( 4) (5 ) (6) (7) (8) (9) ( 10)
C TITL TIDE NODE ELEM REPP FORC FLW BDR END
C
GO TO ( 200, 220, 240, 320, 340, 360, 520, 540, 560, 580),I
C
200 READ(3,1100) ( DT80( J), J = 1, 20)
1100 FORMAT(20A4)
GO TO 900
C
C ----- /TITL/ ---
C
220 READ(3,1200) TITL1
1200 FORMAT( 8X,A80)
GO TO 900
C
C ----- /TIDE/ ---
C
240 READ(3,1300) ( WK1D( I),I=4, 7),LP1
1300 FORMAT( 20X,4F5.0,I5 )
WRITE(6,*) ' '
WRITE(6,2000)
2000 FORMAT(' ** DT = ? ')
** READ(5,*) DT
DT = 1.15D0 ! DT = 1.15 Loop = 100
C
WRITE(6,*) '** NLOP = ? '
** READ(6,*) NLOP
NLOP = 100
RED = DT* FLOAT( NLOP )
C
IF ( LP1.LE.0) LP1 = 1
GO TO 900
C
C /NODE/
C
320 READ(3,1400) ( IW1D( I),( WK2D( I,J),J=1,2),I=1,2 )
1400 FORMAT(10X,2( I5,5X,2F10.0))
C
DO 330 I = 1, 2
C
IF ( IW1D( I).LE.0) GO TO 330
X( IW1D( I)) = WK2D( I,1)* SCX
Y( IW1D( I)) = WK2D( I,2)* SCY
C
330 CONTINUE
GO TO 900
C
C ----- /ELEM/ ---
C
340 READ(3,1500) ( ( IW2D( I,J),J = 1, 5 ),I = 1, 2 )
1500 FORMAT( 10X,2(5I5,5X))
C
CALL ZRV ( 20, WK2D )
C
DO 350 I = 1, 2
C
IF ( IW2D( I,1).LE.0) GO TO 350
C
DO 345 J = 1, 4
C
K = IW2D( I,J+1)
NODE( IW2D( I,1),J) = K
WK2D( 1,J) = X( K)
WK2D( 2,J) = Y( K)
C
345 CONTINUE
C
WK1D( 1) = WK2D( 1,2) - WK2D( 1,1)
WK1D( 2) = WK2D( 2,3) - WK2D( 2,1)
AT( IW2D( I,1)) = WK1D( 1)* WK1D( 2)
C
350 CONTINUE
GO TO 900
C
C /REPP/
C
360 READ(3,1600) ( IW1D( I),I = 1, 6 )
1600 FORMAT(10X,6(I5,5X))
C
C --------------------------------------
DO 370 I = 1, 6
C
K = IW1D( I)
IF ( K.EQ.0) GO TO 370
NRP = NRP + 1
IRP( NRP) = K
C
370 CONTINUE
C --------------------------------------
GO TO 900
C
C ----- /FORC/ ---
C
520 READ(3,1700) L, M, ( IW1D( I), I = 1, 2 )
1700 FORMAT(10X,4I5 )
C
C --------------------------------------
DO 522 I = 1, 5
C
IF ( LNE( 1,I,1).EQ.L ) GO TO 524
C
522 CONTINUE
C --------------------------------------
C
524 K1 = LNE( 1,I,3)
K2 = LNE( 1,I,4)
NFOR = NFOR + 1
C
C --------------------------------------
DO 526 I = 1, 2
C
J = K1
IF ( I.EQ.2) J = K2
IFTP( J) = M
IFOR( J) = IW1D( I)
C
526 CONTINUE
C --------------------------------------
GO TO 900
C
C ----- /FLW/ ---
C
540 READ(3,1800) L, ( IW1D( I), I = 1, 12 )
1800 FORMAT(10X,13I5 )
C
C --------------------------------------
DO 542 I = 1, 5
C
IF ( LNE( 1,I,1).EQ.0) GO TO 544
IF ( LNE( 1,I,1).EQ.L) GO TO 546
C
542 CONTINUE
C --------------------------------------
C
544 LNE( 1,I,1) = L
LNE( 1,I,3) = NFLW + 1
C
546 CONTINUE
C
C --------------------------------------
DO 550 J = 1, 12
C
K = IW1D( J)
IF ( K.EQ.0) GO TO 550
NFLW = NFLW + 1
LNE( 1,I,2) = LNE( 1,I,2) + 1
LNE( 1,I,4) = NFLW
IFLW( NFLW ) = K
C
550 CONTINUE
C --------------------------------------
GO TO 900
C
C ----- /BDR/ ---
C
560 READ(3,1900) L, ( IW1D( I), I = 1, 11)
1900 FORMAT(10X,12I5 )
C
C --------------------------------------
DO 562 I = 1, 5
C
IF ( LNE( 2,I,1).EQ.0) GO TO 564
IF ( LNE( 2,I,1).EQ.L ) GO TO 566
C
562 CONTINUE
C --------------------------------------
C
564 NBDR = NBDR + 1
IF ( NBDR.GT.100) WRITE(8,2100)
2100 FORMAT('*** NBDR.GT.100, STOP ***')
C
IF ( NBDR.GT.100) STOP
LNE( 2,I,1) = L
LNE( 2,I,3) = NBDR + 1
C
566 CONTINUE
C
C --------------------------------------
DO 570 J = 1, 11
C
K = IW1D( J)
IF ( K.EQ.0) GO TO 570
NBDR = NBDR + 1
LNE( 2,I,2) = LNE( 2,I,2 ) + 1
LNE( 2,I,4) = NBDR
IBDR( NBDR )= K
C
570 CONTINUE
C --------------------------------------
GO TO 900
C
C ----- /END/ ---
C
580 CONTINUE
C
990 CONTINUE
C
REWIND 3
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TOUT ( U, V, Z, UO, VO, ZO, NODE, T, IOV, X, Y, LOP,
1 * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NODE( NEL,3)
C
COMMON /RPCM/ NRP, IRP( 20)
COMMON /WVCM/ DT, RED, LP1, NLOP
COMMON /AREA/ NNP, NEL, SCX, SCY
C
DIMENSION U ( NNP), V( NNP), Z ( NNP), UO( NNP), VO( NNP),
1 ZO( NNP), X( NNP), UR( 20), VR( 20), ZR( 20)
C
H = 10.D0
C
C ----- Move New to Old ---
C
DO 100 I = 1, NNP
C
UO( I) = U( I)
VO( I) = V( I)
ZO( I) = Z( I)
C
100 CONTINUE
C
IOV = 0
VMX = 0.D0
C
DO 200 I = 1, NNP
C
VMX = DMAX1( VMX, DABS( Z( I)))
C
200 CONTINUE
C
IF ( VMX.GT. 1.5) IOV = 1
IF ( VMX.LT.-1.5) IOV = 1
C
C ----- REP-point value ---
C
IF ( NRP.LE.0) GO TO 35
C
DO 300 I = 1, NRP
C
J = IRP( I)
ZR( I) = Z( J)
UR( I) = U( J)
VR( I) = V( J)
C
300 CONTINUE
C
35 CONTINUE
C
C ----- Time Calc. ---
C
SEC = T
C
C ----- ALL point value ---
C
IF ( IOV.GT.0) GO TO 50
GO TO 345
C
50 WRITE(8,2100) SEC, ( I, U( I),I = 1, NNP)
WRITE(8,2200) SEC, ( I, Z( I),I = 1, NNP)
2100 FORMAT('1',5X,'VELOCITY ---- UE TIME =',
1 F5.2,' S' //(5X,5( I5,F12.6,2X)))
2200 FORMAT('1',5X,' SUII ---- ZE TIME =',
1 F5.2,' S' //(5X,5( I5,F12.6,2X)))
C
345 CONTINUE
C
C ----- REP point value ---
C
IF ( NRP.LE.0) GO TO 333
C
WRITE(8,2300) SEC, ( IRP( I),UR( I),ZR( I),I = 1, NRP )
2300 FORMAT(/,'*** REP-POINT WAVE VALUES ***',' TIME =',
1 F5.2,' S ',/,(5X,I5,2F15.6))
C
C -----------------------------------------------
C
CALL CNER ( UR, SEC, NRP, X, Y, T, LOP )
C
C -----------------------------------------------
333 CONTINUE
IF ( IOV.GT.0) WRITE(6,2400)
2400 FORMAT(/,' *** Diverged *** STOP **')
C
IF ( IOV.GT.0) STOP
IF ( LOP.LT.NLOP) RETURN 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VEC1 ( AT, Q, NODE, F, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- 1st Approximation ---
C
INTEGER*2 NODE( NEL,3)
C
DIMENSION AT( NEL), Q( NEL), F( NNP)
C
C -----------------------------------------------------
DO 100 I = 1, NNP
C
F( I) = 0.D0
C
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.D0
C
200 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VECU ( F, Z, AE, NNP, NODE, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- F1 & F2 ( 1st Approxmation ) ---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION F( NNP), A( 4), Z( NNP), AE( NEL)
C
G = 10.D0
C ------------------------------------------------
DO 100 I = 1, NNP
C
F( I) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NEL
C
CALL DVALA ( A, Z, AE, NODE, NNP, NEL, I )
C
C ------------------------------------------------
DO 200 J = 1, 4
C
F( NODE( I,J)) = F( NODE( I,J)) + G* A( J)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VECV ( F, Z, BE, NNP, NODE, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- F1 & F2 ( 1st Appr.) ---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION F( NNP), A( 4), Z( NNP), BE( NEL)
C
G = 10.D0
C ------------------------------------------------
DO 100 I = 1, NNP
C
F( I) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NEL
C
CALL DVALB ( A, Z, BE, NODE, NNP, NEL, I )
C
C ------------------------------------------------
DO 200 J = 1, 4
C
F( NODE( I,J)) = F( NODE( I,J)) + G* A( J)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VECZ ( F, U, V, AE, BE, NNP, NODE, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- F1 & F2 ( 1st Approximation ) ---
C
INTEGER*2 NODE( NEL,4)
C
DIMENSION F ( NNP), U ( NNP), V ( NNP), AE( NEL), BE( NEL),
1 AA( 4), AB( 4)
C
H = 10.D0
C ----------------------------
DO 100 I = 1, NNP
C
F( I) = 0.D0
C
100 CONTINUE
C ----------------------------
DO 200 I = 1, 4
C
AA( I) = 0.D0
AB( I) = 0.D0
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 I = 1, NEL
C
CALL DVALA ( AA, U, AE, NODE, NNP, NEL, I )
C
CALL DVALB ( AB, V, BE, NODE, NNP, NEL, I )
C
C ----------------------------------------------------------
DO 300 J = 1, 4
C
F( NODE( I,J)) = F( NODE( I,J)) + H* ( AA( J) + AB( J))
C
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZARY ( NR, NC, A )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A( NR, NC)
C
NNR = NR
NNC = NC
IF ( NR.LE.1) NNR = 1
IF ( NC.LE.1) NNC = 1
C
DO 100 J = 1, NNR
DO 100 I = 1, NNC
C
A( J, I ) = 0.D0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZRV ( NV4, AV )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION AV( NV4)
C
IF ( NV4.LE.1) NV4 = 1
C
DO 100 I = 1, NV4
C
AV( I) = 0.D0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZVC2 ( NV2, NV )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NV( NV2)
C
IF ( NV2.LE.1) NV2 = 1
C
DO 100 I = 1, NV2
C
NV( I) = 0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************