戻る

C ********************************************************************* C
C C
PROGRAM MAIN
C C
C W A V E - F D M C
C C
C ( Version 3.5 ) C
C C
C Wave Analysis by Finite Difference Method C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ********************************************************************* C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /SZE/ M, N
C
DIMENSION Z ( 600), FX ( 600), FY ( 600), U ( 600),
1 V ( 600), HX ( 600), HY ( 600), HHX( 80,10),
2 HHY( 80,10), ICD( 600), IID( 80,10)
C
OPEN ( 8, FILE='WAVE_FDM_OUT' )
OPEN ( 11, FILE='VAR_DAT')
c
CALL TTLE
C
CALL INPT ( HHX, HHY, DT, DS, IID )
C
CALL ITSB ( U, V, Z, FX, FY, HX, HY, HHX, HHY, ICD,
1 IID )
C
CALL WAVD ( FX, FY, U, V, HX, HY, DT, DS, Z, ICD )
C
CLOSE ( 8)
CLOSE ( 11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE ERCP ( T, ZZ, DS, SMZ, LOP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /SZE/ M, N
COMMON /ABC/ MM, EXL( 30), Z( 30)
C
DIMENSION ZZ( 80,10), X( 30)
C
MM = 15
IQQ = 3
C
IF ( M.EQ.80) THEN
MM = 27
IQQ = 4
END IF
C
SMZ = 0.D0
DO 100 I = 1, MM-2
C
X( I+1 ) = DFLOAT( I-1)* DS
C
100 CONTINUE
C
IF ( LOP.EQ.1) THEN
WRITE(6,*)
WRITE(6,2000) ( X( I), I = 2, MM-1 )
2000 FORMAT(' * X =',10F7.3)
WRITE(6,*)
END IF
C
DO 200 I = 2, MM-1
C
Z( I) = ZZ( I,IQQ)* 0.01D0
C
200 CONTINUE
C
C ----- Analytical Results ---
C
RRT1 = 10.D0* T
C
DO 300 I = 2, MM-1
C
EXL( I) = DSIN( T - X( I)/10.D0)
IF ( X( I).GT.RRT1) EXL( I) = 0.D0
C
300 CONTINUE
C
C ----- Relative Error ---
C
SMZ = 0.D0
DO 400 I = 2, MM-1
C
SMZ = SMZ + DABS( Z( I) - EXL( I))
C
400 CONTINUE
C
SMZ = SMZ/ DFLOAT( MM-2) ! Aver. Estimation
SMZ = SMZ/ 2.D0 ! Max.: 1. x 2 Considering "+" & "-"
SMZ = SMZ* 100.D0 ! (%) ! Rel. ERCPor
C
IF ( MOD( LOP,10).EQ.0) THEN
WRITE(6,2100) LOP, T, SMZ, Z( 2)
2100 FORMAT(' * LOP =', I4,' T=',F7.3,' Rel.Error =',F7.3,'(%)',
1 ' Z( 2)=',F7.3)
C
WRITE(8,*) LOP, T, SMZ, Z( 2)
END IF
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ITSB ( U, V, Z, FX, FY, HX, HY, HHX, HHY, ICD, IID )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION U ( 600), V ( 600), Z ( 600), FX ( 600),
1 FY ( 600), HX ( 600), HY ( 600), HHX( 80,10),
2 HHY( 80,10), ICD( 600), IID( 80,10)
C
COMMON /SZE/ M, N
C
K = 0
C --------------------------------------
DO 100 J = 1, N
DO 100 I = 1, M
C --------------------------------------
C
K = K + 1
C
HX( K) = HHX( I,J)
HY( K) = HHY( I,J)
C
Z ( K) = 0.D0
FX( K) = 0.D0
FY( K) = 0.D0
U ( K) = 0.D0
V ( K) = 0.D0
C
ICD( K) = IID( I,J)
C
100 CONTINUE
C --------------------------------------
C
RETURN
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 FDM 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 INPT ( HX, HY, DT, DS, ICD )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /SZE/ M, N
COMMON /CND/ ISR, IED, JSR( 100), JED( 100)
COMMON /FOR/ NFR( 4), IFR( 80,4), JFR( 80,4)
COMMON /FRC/ A ( 80,4), T ( 80,4), ALP( 80,4)
C
DIMENSION H ( 80,10), HX ( 80,10), HY ( 80,10),
1 ICD( 80,10), NCRD( 20), CRDN( 13 ),
2 CRD( 19 ), ICON( 40), HH ( 40)
C
CHARACTER*4 CRDN, CID, CRD
CHARACTER*1 QQ
C
DATA CRDN /'A','B','C','D','E','F','G','H','I','J','K','L','M'/
C
PAI = 3.14159 26535 89793 D0
C
WRITE(6,*) ' '
WRITE(6,*) '*** L or S ( ?) ***'
C
** READ(5,*) QQ
QQ = 'L'
** QQ = 'S'
C
IF ( QQ.EQ.'S') OPEN ( 7, FILE ='WAVE.DAT' )
IF ( QQ.EQ.'L') OPEN ( 7, FILE ='LWAVE.DAT')
C
READ(7,1500) M, N
1500 FORMAT(10X,2I5)
WRITE(6,2000) M, N
2000 FORMAT(/,' * NX =',I4,' NY =',I4)
WRITE(6,2300)
2300 FORMAT(/,' * DELT = ? ')
*** READ(5,*) DT
C
DT = 0.3D0
C ---------------
C
DS = 7.5D0* 100.D0
IF ( M.EQ.80) DS = 3.75D0* 100.D0
C
C ----------------------------
DO 120 I = 1, 20
C
NCRD( I) = 0
C
120 CONTINUE
C ----------------------------
DO 140 I = 1, 4
C
NFR( I) = 0
C ----------------------------
DO 140 J = 1, 80
C
IFR( J,I) = 0
JFR( J,I) = 0
C
140 CONTINUE
C ----------------------------
C
ISR = 0
IED = 0
C ----------------------------
DO 160 I = 1, 100
C
JSR( I) = 0
JED( I) = 0
C
160 CONTINUE
C ----------------------------
1000 CONTINUE
READ(7,1200) CID, CRD
1200 FORMAT(A1,19A4)
C
C --------------------------------------
DO 180 I = 1, 9
C
ID = I
IF ( CRDN( I).EQ.CID) GO TO 1100
C
180 CONTINUE
C --------------------------------------
GO TO 9000
C
1100 CONTINUE
BACKSPACE 7
C
NCRD( ID) = NCRD( ID) + 1
GO TO ( 1, 1, 1, 1, 1, 6, 1, 8, 9 ), ID
C
1 CONTINUE
GO TO 1000
C
6 CONTINUE
READ(7,1300) I1, I2, I3, W1, W2, W3, I4, I5, I6, W4, W5, W6
1300 FORMAT(6X,I1,2I3,3F10.3,I1,2I3,3F10.3)
C
C ------------------------------------------------
DO 200 I = 1, 4
C
IF ( I1.NE.I) GO TO 220
I7 = NFR( I) + 1
NFR( I) = I7
IFR( I7,I) = I2
JFR( I7,I) = I3
A ( I7,I) = W1
T ( I7,I) = W2
ALP( I7,I) = W3
C
220 CONTINUE
IF ( I4.NE.I) GO TO 200
C
I8 = NFR( I) + 1
NFR( I) = I8
IFR( I8,I) = I5
JFR( I8,I) = I6
A ( I8,I) = W4
T ( I8,I) = W5
ALP( I8,I) = W6
C
200 CONTINUE
C ------------------------------------------------
GO TO 1000
C
8 CONTINUE
READ(7,1400) I, N1,( ICON( J), J = 1, 35 )
1400 FORMAT(4X,2I3,35I2)
C
N2 = N1 + 34
IF ( N.LT.N2) N2 = N
C
C ----------------------------
DO 240 J = N1, N2
C
JJ = J - N1 + 1
ICD( I,J) = ICON( JJ)
C
240 CONTINUE
C ----------------------------
GO TO 1000
C
9 CONTINUE
READ(7,1600) I, N1, ( HH( J), J = 1, 5 )
1600 FORMAT(4X, 2I3, -2P14F5.0)
C
N2 = N1 + 13
IF ( N.LT.N2 ) N2 = N
C
C ----------------------------
DO 260 J = N1, N2
C
JJ = J - N1 + 1
H( I,J) = HH( JJ)
C
260 CONTINUE
C ----------------------------
GO TO 1000
C
9000 CONTINUE
C
C ------------------------------------------------
DO 270 I = 1, 4
C
NNN = NFR( I)
IF ( NNN.EQ.0) GO TO 270
C
C --------------------------------------
DO 275 J = 1, NNN
C
T ( J,I) = PAI* 2.D0/( T( J,I)* 3600.D0)
ALP( J,I) = ALP( J,I)* PAI/ 180.D0
C
275 CONTINUE
C --------------------------------------
C
270 CONTINUE
C ------------------------------------------------
C
C ----- ISR, IED, JSR, JED ---
C
C ------------------------------------------------
DO 280 I = 1, M
C
C --------------------------------------
DO 285 J = 1, N
C
IF ( ICD( I,J).LE.0) GO TO 285
JED( I) = J
IF ( JSR( I).EQ.0) JSR( I) = J
C
285 CONTINUE
C --------------------------------------
C
IF ( JSR( I).EQ.0) GO TO 280
IED = I
IF ( ISR.EQ.0) ISR = I
C
280 CONTINUE
C ------------------------------------------------
C
C ----------------------------------------------------------
DO 300 I = 1, M
DO 300 J = 1, N
C ----------------------------------------------------------
C
IC = ICD( I,J)
GO TO ( 907, 903, 905, 900, 906, 901, 900, 900, 900, 900, 903,
1 907, 905, 904, 902, 904, 902 ), IC
C
900 CONTINUE
ICD( I,J) = 0
GO TO 300
C
901 CONTINUE
ICD( I,J) = 1
GO TO 300
C
902 CONTINUE
ICD( I,J) = 2
GO TO 300
C
903 CONTINUE
ICD( I,J) = 3
GO TO 300
C
904 CONTINUE
ICD( I,J) = 4
GO TO 300
C
905 CONTINUE
ICD( I,J) = 5
GO TO 300
C
906 CONTINUE
ICD( I,J) = 6
GO TO 300
C
907 CONTINUE
ICD( I,J) = 7
C
300 CONTINUE
C ----------------------------------------------------------
C
NN = NFR( 1)
IF ( NN.EQ.0) GO TO 8000
C
C --------------------------------------
DO 310 I = 1, NN
C
II = IFR( I,1)
JJ = JFR( I,1)
H( II,JJ) = H( II,JJ+1)
C
310 CONTINUE
C --------------------------------------
C
8000 CONTINUE
NN = NFR( 2)
IF ( NN.EQ.0) GO TO 8100
C
C --------------------------------------
DO 320 I = 1, NN
C
II = IFR( I,2)
JJ = JFR( I,2)
H( II,JJ) = H( II+1,JJ)
C
320 CONTINUE
C --------------------------------------
C
8100 CONTINUE
NN = NFR( 3)
IF ( NN.EQ.0) GO TO 8200
C
C --------------------------------------
DO 340 I = 1, NN
C
II = IFR( I,3)
JJ = JFR( I,3)
H( II,JJ) = H( II,JJ-1)
H( II,JJ+1) = H( II,JJ-1)
C
340 CONTINUE
C --------------------------------------
C
8200 CONTINUE
NN = NFR( 4)
IF ( NN.EQ.0) GO TO 8300
C
C --------------------------------------
DO 360 I = 1, NN
C
II = IFR( I,4)
JJ = JFR( I,4)
H( II, JJ) = H( II-1,JJ)
H( II+1,JJ) = H( II-1,JJ)
C
360 CONTINUE
C --------------------------------------
C
8300 CONTINUE
C
C ----- Depth ---
C
C ------------------------------------------------
DO 380 I = 1, M
DO 380 J = 1, N
C
HX( I,J) = ( H( I,J) + H( I+1, J))* 0.5D0
HY( I,J) = ( H( I,J) + H( I,J +1))* 0.5D0
C
380 CONTINUE
C ------------------------------------------------
C
WRITE(6,2200) DT, ISR, IED
2200 FORMAT(' * DT =',F6.3,' ISR =',I3,' IED =',I4,/)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WAVD ( FX, FY, U, V, HX, HY, DT, DS, RZ, ICD )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /SZE/ M, N
COMMON /ABC/ MM, EXL( 30), Z( 30)
COMMON /CND/ ISR, IED, JSR( 100), JED( 100)
COMMON /FOR/ NFR( 4), IFR( 80,4), JFR( 80,4)
COMMON /FRC/ A ( 80,4), T ( 80,4), ALP( 80,4)
C
DIMENSION RZ ( 600), FX ( 600), FY( 600), U ( 600),
1 HX ( 600), HY ( 600), V ( 600), ICD( 600),
2 ZR ( 80,10)
C
WRITE(6,*)' ** Final Time ( ?)'
WRITE(6,*) ' '
C
** READ(5,*) DDX
DDX = 20.D0
C
WRITE(6,*) ' ** Correction = 0 1 2 **'
** READ(5,*) IS
IS = 0
C -----------
C
F2 = 1.D0
F1 = 0.D0
BB = 1000.D0* DT/ DS
C
IF ( IS.EQ.1) THEN
F1 = 0.5D0 - 0.5D0/ BB
F2 = 1.
END IF
C
IF ( IS.EQ.2) THEN
F1 = ( 4.D0* BB** 2 - 2.D0* BB )/4.D0 /BB** 2
F2 = 2.D0* BB/ DSQRT( 4.D0* ( BB ** 2 ))
END IF
C
IF ( IS.EQ.3) THEN
F1 = ( BB* BB - BB)/2.D0/ BB* BB
F2 = BB/ DSQRT( 2.D0* BB* BB )
END IF
C
C1 = DT / DS
C3 = 1000.D0* C1* ( -1.D0 )* F2
WKK = 1000.D0* F1* C1* C1* 1000.D0
RTM = 0.D0
CMN = 1.D0
DSS = DS* 0.01D0
C
LOP = 0
C
990 RTM = RTM + DT
LOP = LOP + 1
C
C ------ Water Level ---
C
C --------------------------------------------------------------
DO 100 I = ISR, IED
C
IF ( JSR( I).GT.99) JSR( I) = 2
JS = JSR( I)
JE = JED( I)
IJ = ( JS-2)* M + I
C
C ---------------------------------------------------------------
DO 100 J = JS, JE
C
IJ = IJ + M
IF ( MOD( ICD( IJ),2).EQ.0) GO TO 100
C
RZ( IJ) = RZ( IJ)
1 + C1 * ( FX( IJ) - FX( IJ+M) + FY( IJ+1) - FY( IJ))
2 + WKK* ( RZ( IJ+1) - 2.D0* RZ( IJ) + RZ( IJ-1))
3 + WKK* ( RZ( IJ+M) - 2.D0* RZ( IJ) + RZ( IJ-M))
C
IF ( DABS( RZ( IJ)).GT.150.D0) THEN
WRITE(6,*) ' ** Over 1.5 *** STOP * LOP =', LOP
STOP
END IF
C
100 CONTINUE
C ---------------------------------------------------------------
C
C ----- Sine Curve ---
C
C ----------------------------------------------------------
DO 200 J = 1, 4
C
IE = NFR( J)
IF ( IE.EQ.0) GO TO 200
C
C ------------------------------------------------
DO 220 I = 1, IE
C
IJ = ( JFR( I,J)-1)* DFLOAT( M) + IFR( I,J)
RZ( IJ) = A( I,J)* DSIN( RTM* T( I,J) - ALP( I,J))
C
220 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
C
C ----- FX, FY ---
C
IO = 0
C ------------------------------------------------
DO 300 I = ISR, IED
C
JS = JSR( I)
JE = JED( I)
IJ = ( JS-2)* M + I
C
C ------------------------------------------------
DO 300 J = JS, JE
C
IJ = IJ+M
ICU = ICD( IJ)
ICV = MOD( ICU - MOD( ICU,2),4)
RZIJ = RZ ( IJ)
C
C ----- FX ---
C
IF ( ICU.LT.4) GO TO 900
IJM = IJ-M
RZIJM = RZ( IJM)
DEPX = ( RZIJ + RZIJM)* 0.5D0 + HX( IJ)
C
IF ( DEPX.LE.CMN) GO TO 920
C
FXIJ = FX( IJ) + C3* DEPX* ( RZIJ - RZIJM )
1 + WKK* ( FX( IJ+1) - 2.D0* FX( IJ) + FX( IJ-1))
2 + WKK* ( FX( IJ+M) - 2.D0* FX( IJ) + FX( IJ-M))
C
IF ( DABS( FXIJ).LT.1.0E -20) FXIJ = 0.D0
C
UIJ = FXIJ/ DEPX
FX( IJ) = FXIJ
U ( IJ) = UIJ
GO TO 900
C
920 CONTINUE
FX( IJ) = 0.D0
U ( IJ) = 0.D0
C
900 CONTINUE
C
C ----- FY ---
C
IF ( ICV.EQ.0) GO TO 940
IMJ = IJ-1
RZIMJ = RZ( IMJ)
DEPY = ( RZIMJ + RZIJ)* 0.5D0 + HY( IJ)
C
IF ( DEPY.LE.CMN) GO TO 950
C
FYIJ = FY( IJ) + C3* DEPY* ( RZIMJ - RZIJ)
1 + WKK* ( FY( IJ+1) - 2.D0* FY( IJ) + FY( IJ-1))
2 + WKK* ( FY( IJ+M) - 2.D0* FY( IJ) + FY( IJ-M))
C
IF ( DABS( FYIJ).LT.1.D -20) FYIJ = 0.D0
C
VIJ = FYIJ/ DEPY
FY( IJ) = FYIJ
V ( IJ) = VIJ
GO TO 940
C
950 CONTINUE
FY( IJ) = 0.D0
V ( IJ) = 0.D0
C
940 CONTINUE
C
300 CONTINUE
C ------------------------------------------------
C
C ----------------------------------------------------------
DO 700 J = 1, 4
C
IE = NFR( J)
IF ( IE.EQ.0) GO TO 700
GO TO ( 1600, 1610, 1620, 1630), J
C
1600 CONTINUE
C ------------------------------------------------
DO 720 I = 1, IE
C
IJ = ( JFR( I,J)-1)* M + IFR( I,J)
IJP = IJ + M
FX( IJ) = FX( IJP)
FY( IJ-M) = FY( IJ)
V ( IJ-M) = V ( IJ)
C
720 CONTINUE
C ------------------------------------------------
GO TO 700
C
1610 CONTINUE
C ------------------------------------------------
DO 740 I = 1, IE
C
IJ = ( JFR( I,J)-1)* M + IFR( I,J)
IK = IJ + 1
FX( IJ-1) = FX( IJ)
FY( IJ) = FY( IK )
U ( IJ-1) = U ( IJ)
C
740 CONTINUE
C ------------------------------------------------
GO TO 700
C
1620 CONTINUE
C
C ------------------------------------------------
DO 760 I = 1, IE
C
IJ = ( JFR( I,J)-1)* M + IFR( I,J)
IJP = IJ + M
FX( IJP) = FX( IJ)
FY( IJP) = FY( IJ)
U ( IJP) = U ( IJ)
C
760 CONTINUE
C ------------------------------------------------
GO TO 700
C
1630 CONTINUE
C
C ------------------------------------------------
DO 780 I = 1, IE
C
IJ = ( JFR( I,J)-1)* M + IFR( I,J)
IK = IJ + 1
FX( IK ) = FX( IJ)
FY( IK ) = FY( IJ)
V ( IK ) = V ( IJ)
C
780 CONTINUE
C -----------------------------------------------
C
700 CONTINUE
C
C ----- Output ---
C
K = 0
C
C ---------------------------
DO 800 J = 1, N
DO 800 I = 1, M
C
K = K + 1
ZR( I,J) = RZ( K)
C
800 CONTINUE
C -------------------------------------------
C
CALL ERCP ( RTM, ZR, DSS, SMZ, LOP )
C
C -------------------------------------------
C
ET = ET + SMZ
IF ( RTM.LE.DDX) GO TO 990
ET = DSQRT( ET/ DFLOAT( LOP))
C
WRITE(6,*) ' '
WRITE(6,2100) RTM, DT, ET
2100 FORMAT(' *** Final *** T =',F7.3,' DT =',F6.3,' Et =',F8.3)
C
WRITE(11,*) MM
WRITE(11,*) ( Z ( I), I = 2, MM-1 )
WRITE(11,*) ( EXL( I), I = 2, MM-1 )
C
MM = 15
IF ( M.EQ.80) MM = 27
C
RETURN
END
C ******************************************************************