!PROGRAM Slope ! 時間領域2維任意海底地形造波問題 ! 本程式為公開程式，可不經授權逕行引用，但尚未實作驗証完成，有問題可電郵連絡。 ! 原始論文為http://www.xn--06q65y2zay3jryqs5a29hepz.tw/chou_html/reference/1996chou_05.pdf ! 追加考量海底摩擦、岸壁消波效應，必要解複數連立方程式 ! 將斜坡水平長設立為0時，可變身為空水槽 ! 周宗仁 2024.06.14 REAL,ALLOCATABLE:: X(:),Z(:),S(:),XN(:),ZN(:),XX(:),ZZ(:) REAL,ALLOCATABLE:: H(:,:),G(:,:),WH(:),KD(:),KT(:) REAL,ALLOCATABLE:: AN(:),F2S(:),F2N(:),F2X(:),F2Z(:) REAL,ALLOCATABLE:: RNUM(:),S0(:),RE(:),HKK(:) COMPLEX,ALLOCATABLE::F1(:),F2(:),F3(:),F4(:),CC(:) COMPLEX,ALLOCATABLE::F1B(:),F2B(:),F3B(:),F4B(:) COMPLEX::RS,IM REAL::L3,E25(25),FF25(25) OPEN(UNIT=1,FILE='WH.OUT',FORM='FORMATTED',STATUS='UNKNOWN') OPEN(UNIT=2,FILE='KD.OUT',FORM='FORMATTED',STATUS='UNKNOWN') ! 以下 UNIT 預留為作業空間，勿佔用，以免造成錯誤 OPEN(UNIT=10,FILE='10.APT',FORM='UNFORMATTED',STATUS='UNKNOWN') OPEN(UNIT=11,FILE='11.APT',FORM='UNFORMATTED',STATUS='UNKNOWN') OPEN(UNIT=12,FILE='12.APT',FORM='UNFORMATTED',STATUS='UNKNOWN') OPEN(UNIT=14,FILE='14.APT',FORM='UNFORMATTED',STATUS='UNKNOWN') OPEN(UNIT=15,FILE='15.APT',FORM='UNFORMATTED',STATUS='UNKNOWN') IM=(0.,1.) PAI=3.1415927 GR=9.81 ITYPE=1 !活塞式造波板 WRITE(*,*) '入射波領域等深域水深 H0=？ 預設為1' WRITE(*,*) '輸入造波板距離 S2/H0=？' ! READ(*,*) S2 WRITE(*,*) '輸入斜坡水平距離 L3/H0=？' ! READ(*,*) L3 WRITE(*,*) '輸入斜坡摩擦係數 BETA3=？' ! READ(*,*) BETA3 WRITE(*,*) '輸入等水深海底摩擦係數 BETA4=？' ! READ(*,*) BETA4 H0=1 ! 基準水深 S2=80 L3=30 S4=S2-L3 S3=SQRT(S4**2+H0**2) BETA3=0.005 BETA4=0.005 ! 設定海底地形(可自行調整)採用1次元素逆時針積分，各節點數如下 WRITE(*,*) '輸入造波板節點數 N1=？' ! READ(*,*) N1 WRITE(*,*) '輸入水面節點數 N2=？' ! READ(*,*) N2 WRITE(*,*) '輸入斜坡點數數 N3=？' ! READ(*,*) N3 WRITE(*,*) '輸入等水深海底點數數 N4=？' ! READ(*,*) N4 N1=21 ! 自設 N2=161 N3=61 N4=101 NT=N1+N2+N3+N4 NT2=N1+N2 NT3=NT2+N3 N11=N1-1 N21=N2-1 N31=N3-1 N41=N4-1 NS=N11+N21+N31+N41 NS2=N11+N21 NS3=NS2+N31 N1X=N11-1 N2X=N21-1 N3X=N31-1 N4X=N41-1 ALLOCATE(H(NT,NT),G(NT,NT),X(NT),Z(NT),S(NT),XN(NT),ZN(NT)) ALLOCATE(KD(N2),KT(N2),WH(N2),CC(NT)) ALLOCATE(F1(N1),F2(N2),F3(N3),F4(N4)) ALLOCATE(F1B(N1),F2B(N2),F3B(N3),F4B(N4)) ALLOCATE(AN(N2),F2S(N2),F2N(N2),F2X(N2),F2Z(N2)) ! 設定元素座標及線段長(1次元素)，自行發揮 ! 以造波板與海底交接處為起始(逆時針) ! 端點座標 DO I=1,N1 X(I)=S2 Z(I)=-H0+H0/N11*(II-1) END DO DO II=1,N2 I=II+N1 X(I)=S2-S2/N21*(II-1) Z(I)=0 END DO DO II=1,N3 I=II+NT2 X(I)=L3/N31*(II-1) Z(I)=-H0/N31*(II-1) END DO DO II=1,N4 I=II+NT3 X(I)=L3+S4/N41*(II-1) Z(I)=-H0 END DO !元素長 DO I=1,N11 S(I)=H0/N11 END DO DO II=1,N21 I=II+N11 S(I)=S2/N21 END DO S3=SQRT(H0**2+L3**2) !斜坡長 DO II=1,N31 I=II+NS2 S(I)=S3/N31 END DO DO II=1,N41 I=II+NS3 S(I)=S4/N41 END DO ! 波浪資訊 WRITE(*,*) '輸入造波形式 IWAVE =1～3' WRITE(*,*) '1：規則波' WRITE(*,*) '2：孤立波 ' WRITE(*,*) '3：單方向不規則波' READ(*,*) IWAVE IF(IWAVE.EQ.1) THEN WRITE(*,*) '輸入入射波有義週期(秒) T31=？' READ(*,*) T31 WRITE(*,*) '輸入入射波有義波高(m) H31=？' READ(*,*) H31 WRITE(*,*) '輸入時間間隔數 NTT=？' READ(*,*) NTT HK=2*PAI/WAVEL(T31,H0) DT=T31/NTT ! 增分間隔 DF=2*PAI/T31 EE=-DF*H31*(SINH(HK)*COSH(HK)+HK)/2./SINH(HK)**2 !活塞式造波函數 END IF IF(IWAVE.EQ.2) THEN WRITE(*,*) '輸入孤立波波高(m) HH=？' READ(*,*) HH WRITE(*,*) '輸入時間間隔數 NTT=？' READ(*,*) NTT DT=NTT GR1=GR*H0 HH=HH/H0 CE=SQRT(GR1*(1+HH)) XO=SQRT(4.*HH/3./(1+HH)) OMEGA=SQRT(GR1)*SQRT(.75*HH*(1+HH)) TC=PAI/OMEGA !孤立波特徵值 XLEFF=9.5766/SQRT(HH) TEFF=XLEFF/CE DDT=TC/DT END IF IF(IWAVE.EQ.3) THEN WRITE(*,*) '波譜形式 ITPSPC(1～4)：？' WRITE(*,*) '1：Bretschneider波譜' WRITE(*,*) '2：Bretschneider-MitsuYasu波譜' WRITE(*,*) '3：JONSWAP譜' WRITE(*,*) '4：自訂' READ(*,*) ITPSPC WRITE(*,*) '輸入入射波有義週期(秒) T31=？' READ(*,*) T31 WRITE(*,*) '輸入入射波有義波高(m) H31=？' READ(*,*) H31 WRITE(*,*) '輸入時間間隔數 NTT=？' READ(*,*) NTT WRITE(*,*) '輸入波譜頻率分割數 NSPEC=？' READ(*,*) NSPEC ALLOCATE(RNUM(NSPEC),S0(NSPEC),RE(NSPEC),HKK(NSPEC)) HK=2*PAI/WAVEL(T31,H0) DTT=NTT IF(ITPSPC.EQ.4) THEN ! 自訂 WRITE(*,*) 'FPEAK=？' READ(*,*) FP WRITE(*,*) 'SPEAK=？' READ(*,*) SP WRITE(*,*) 'E25(I)=？' READ(*,*) (E25(I),I=2,25) END IF WRITE(*,*) '最小週期 TMIN=？' READ(*,*) TMIN WRITE(*,*) '最大週期 TMAX=？' READ(*,*) TMAX ! IXY=7 FMAX=1/TMIN FMIN=1/TMAX DF=(FMAX-FMIN)/NSPEC DT=1./(2.*NSPEC*DF) DDT=TMIN/DTT IF(ITPSPC.EQ.1) CALL BM0(S0,NSPEC,NSPEC0,H31,T31,DT) IF(ITPSPC.EQ.2) CALL BMS(S0,NSPEC,NSPEC0,H31,T31,DT) IF(ITPSPC.EQ.3) CALL JONSWAP(S0,NSPEC,NSPEC0,H31,T31,DT) IF(ITPSPC.EQ.4) CALL SPECIN(S0,NSPEC,NSPEC0,DT,E25,FF25,FP,SP) CALL WVMOD(RE,NSPEC,NSPEC0,H0,DT,TMAX,TMIN,ITYPE,HKK) DO J=1,NSPEC CALL RANDU(IXY,IY,FL) IXY=IY RNUM(J)=PAI2*FL END DO END IF ! 開始演算 WRITE(*,*) '輸入是否為第1次計算？是=0 否<>0' READ(*,*) ICASE IF(ICASE.EQ.0) THEN DO I=1,N1 F1(I)=0 F1B(I)=0 END DO DO I=1,N2 F2(I)=0 F2B(I)=0 END DO WRITE(*,*) '開始時刻START=0' START=0 WRITE(*,*) '輸入結束時刻EEND=？' READ(*,*) EEND END IF IF(ICASE.NE.0) THEN WRITE(*,*) '輸入開始時刻START=？' READ(*,*) START WRITE(*,*) '輸入結束時刻EEND=？' READ(*,*) EEND REWIND 12 READ(12) (F2(I),I=1,N2) READ(12) (X(I),Z(I),S(I),I=1,NT) END IF ! 開始反覆計算 ISTART=START*100 !因DO LOOP必要為整數 IEND=EEND*100 ISTART1=ISTART+1 IEND1=IEND IF(ICASE.NE.0) ISTART1=1 IF(ICASE.NE.0) IEND1=IEND-ISTART DO ITT=ISTART1,IEND1 IT1=1+FLOAT(NTT)/100*(ITT-1) IT2=FLOAT(NTT)/100*ITT DO IT=IT1,IT2 TTX=START+IT/FLOAT(NTT) WRITE( *,*) 'IT=',TTX IF(IWAVE.EQ.1) CALL REGULAR(IT,U,EE,DF,DT) IF(IWAVE.EQ.2) CALL SOLITON(IT,U,XO,OMEGA,TC,DDT) IF(IWAVE.EQ.3) CALL FREQ_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM) DO I=1,N1 F1B(I)=U !造波板條件 END DO ! 2維邊界積分方程式 ID=4 ! 邊界頂角數 CALL G2_ALOG_LINEAR(H,G,X,Z,S,NT,ID) !結果存於H ! REWIND 11 ! WRITE(11) ((H(I,J),J=1,NT),I=1,NT) !備用 ! 建立連立方程式 ! H：實數部，G：虛數部，依序存於15 DO I=1,NT DO J=1,NT G(I,J)=0 ! 虛數部歸零 END DO END DO DO I=1,N1 DO J=1,N1 H(I,J)=0 IF(I.EQ.J) H(I,J)=-1 END DO DO JJ=1,N3 J=JJ+NT2 H(I,J)=0 G(I,J)=BETA3*H(I,J) END DO DO JJ=1,N4 J=JJ+NT3 H(I,J)=0 G(I,J)=BETA4*H(I,J) END DO END DO DO II=1,N2 I=II+N1 DO J=1,N1 H(I,J)=0 END DO DO JJ=1,N3 J=JJ+NT2 H(I,J)=0 G(I,J)=BETA3*H(I,J) END DO DO JJ=1,N4 J=JJ+NT3 H(I,J)=0 G(I,J)=BETA4*H(I,J) END DO END DO DO II=1,N3 I=II+NT2 DO J=1,N1 H(I,J)=0 END DO DO JJ=1,N3 J=JJ+NT2 H(I,J)=0 IF(I.EQ.J) H(I,J)=-1 G(I,J)=BETA3*H(I,J) END DO DO JJ=1,N4 J=JJ+NT3 H(I,J)=0 G(I,J)=BETA4*H(I,J) END DO END DO DO II=1,N4 I=II+NT3 DO J=1,N1 H(I,J)=0 END DO DO JJ=1,N2 J=JJ+N1 H(I,J)=0 END DO DO JJ=1,N4 J=JJ+NT3 H(I,J)=0 IF(I.EQ.J) H(I,J)=-1 G(I,J)=BETA4*H(I,J) END DO END DO REWIND 15 WRITE(15) ((H(I,J),J=1,NT),I=1,NT) WRITE(15) ((G(I,J),J=1,NT),I=1,NT) CALL CMINVS(H,G,NT,NT) ! 結果：實數矩陣存於H，實數矩陣存於G ! 連立方程式右邊已知項 DO I=1,NT CC(I)=0 !歸零 END DO DO I=1,N1 RS=0 DO J=1,N1 RS=RS+H(I,J)*U END DO CC(I)=RS END DO DO II=1,N2 I=II+N1 RS=0 DO J=1,N1 RS=RS+H(I,J)*U END DO CC(I)=RS+F2(I) END DO DO II=1,N3+N4 I=II+NT2 RS=0 DO J=1,N1 RS=RS+H(I,J)*U END DO CC(I)=RS END DO ! 結果 DO I=1,N1 RS=(0.,0.) DO J=1,NT RS=RS+CMPLX(H(I,J),G(I,J))*CC(J) END DO F1(I)=RS END DO DO II=1,N2 I=II+N1 RS=(0.,0.) DO J=1,NT RS=RS+CMPLX(H(I,J),G(I,J))*CC(J) END DO F2B(I)=RS END DO DO II=1,N3 I=II+NT2 RS=(0.,0.) DO J=1,NT RS=RS+CMPLX(H(I,J),G(I,J))*CC(J) END DO F3(I)=RS END DO DO II=1,N4 I=II+NT3 RS=(0.,0.) DO J=1,NT RS=RS+CMPLX(H(I,J),G(I,J))*CC(J) END DO F4(I)=RS END DO ! 計算t+1時刻水面波形 DO I=2,N21 II=I+N1 I1=II-1 F2S(I)=(F2(I+1)-F2(I))/S(II-1) AN(I)=-ATAN2((Z(I1)-Z(II)),(X(I1)-X(II))) END DO F2S(1)=-U AN(1)=AN(2) AN(N2)=AN(N21) F2S(N21)=F2S(N2X) F2S(N2)=f2s(N21) IF(L3.EQ.0) F2S(N2)=0 DO I=1,N1 XN(I)=1 ZN(I)=0 END DO DO I=1,N2 II=I+N1 XN(II)=SIN(AN(I)) ZN(II)=COS(AN(I)) END DO DO I=1,N3 II=I+NT2 XN(II)=-H0/SQRT(H0**2+L3**2) ZN(II)=-L3/SQRT(H0**2+L3**2) END DO DO I=1,N4 II=I+NT3 XN(II)=0 ZN(II)=-1 END DO DO I=1,N2 SN=SIN(AN(I)) CS=COS(AN(I)) F2X(I)=F2B(I)*SN-F2S(I)*CS F2Z(I)=F2B(I)*CS+F2S(I)*SN END DO DO I=1,N2 II=N1+I Z(II)=Z(II)+F2Z(I)*DDT X(II)=X(II)+F2X(I)*DDT END DO IF(L3.EQ.0) X(NT2)=0 DO I=1,N2 II=I+N1 R=F2S(I) T=F2B(I) IF(I.EQ.N2) Z(II)=Z(II)/2 F2(I)=F2(I)+DDT*(0.5*(R*R+T*T)-Z(II)*GR) END DO HW1=H0+Z(N1+1) DO I=1,N1 X(I)=X(I)+U*DDT Z(I)=-H0+HW1/N11*(I-1) END DO X(N1+1)=X(N1) HW2=H0+Z(NT2) DO I=1,N3 II=I+NT2 Z(II)=Z(NT2)-HW2/N31*(I-1) IF(L3.EQ.0) X(II)=0 IF(L3.NE.0) X(II)=X(NT2)+(L3-X(NT2))/N31*(I-1) END DO DO I=1,N4 II=I+NT3 X(II)=L3+(X(1)-L3)/N41*(I-1) END DO DO I=1,N11 S(I)=Z(I+1)-Z(I) END DO DO I=1,N21 II=I+N11 I1=II+1 I2=II+2 S(II)=SQRT((X(I2)-X(I1))**2+(Z(I2)-Z(I1))**2) END DO DO I=1,N31 II=I+NS2 S(II)=SQRT(HW2**2+(L3-X(NT2))**2)/N31 END DO DO I=1,N41 II=NS3+I S(II)=(X(1)-L3)/N41 END DO S(NS + 1) = 0 !填空 S(NS + 2) = 0 S(NS + 3) = 0 S(NS + 4) = 0 DO I=1,N2 WH(I)=REAL(F2(I)) END DO WRITE(1,1) (WH(I),I=1,N2) 1 FORMAT(10F5.2) REWIND 12 WRITE(12) (F2(I),I=1,N2) WRITE(12) (X(I),Z(I),S(I),I=1,NT) END DO !IT END DO !ITT STOP END ! 副程式 SUBROUTINE G2_ALOG_LINEAR(H,G,X,Y,S,NT,ID) ! 計算[K]=[H]^[G] ! 使用 UNFORMATTED UNIT : 7,8,10 ! 結果存放在 UNIT 10 ((H(I,J),J=1,NT)I=1,NT) ! ID：邊界頂角數 REAL X(NT),Y(NT),S(NT),H(NT,NT),G(NT,NT) PAI=3.1415927 SSMAX=S(1) DO I=1,NT-ID IF(S(I).GT.SSMAX) SSMAX=S(I) END DO DO I=1,NT INDEX=0 XI=X(I) YI=Y(I) DO J=1,NT J1=J-1 J2=J+1 XJ=X(J) YJ=Y(J) IF(J.EQ.1) THEN XJ1=0 YJ1=0 ELSE XJ1=X(J1) YJ1=Y(J1) END IF IF(J.EQ.NT) THEN XJ2=0 YJ2=0 ELSE XJ2=X(J2) YJ2=Y(J2) END IF IF(J.EQ.1) R1=0 IF(J.NE.1) R1=SQRT((XJ1-XJ)**2+(YJ1-YJ)**2) IF(J.EQ.NT) R2=0 IF(J.NE.NT) R2=SQRT((XJ2-XJ)**2+(YJ2-YJ)**2) IF(XJ.EQ.XJ2.AND.YJ.EQ.YJ2) INDEX=INDEX+1 IF(R2.GT.SSMAX) INDEX=INDEX+1 IF(J.EQ.1) SJ1 =0 IF(J.NE.1) SJ1 =S(J-INDEX-1) SJ2=S(J-INDEX) IF(R2.EQ.0) SJ2=0 IF(R1.EQ.0) SJ1=0 IF(R2.GT.SSMAX) SJ2=0 IF(R1.GT.SSMAX) SJ1=0 IF(I.EQ.J) GOTO 31 CALL HG1(H1,G1,XI,YI,XJ2,YJ2,XJ,YJ,SJ2) CALL HG2(H2,G2,XI,YI,XJ,YJ,XJ1,YJ1,SJ1) GO TO 25 31 H1=0 H2=0 IF(SJ2.EQ.0) G1=0 IF(SJ2.NE.0) G1=SJ2/2/PAI*(3./4.-1./2.*ALOG(SJ2/2)) IF(SJ1.EQ.0) G2=0 IF(SJ1.NE.0) G2=SJ1/2/PAI*(5./4.-3./2.*ALOG(SJ1/2)) 25 HH=H1+H2 GG=G1+G2 IF(R1.EQ.0.OR.R2.EQ.0) HH=HH*2 IF(I.EQ.J) HH=HH+1 H(I,J)=HH G(I,J)=GG END DO ! J END DO ! I EPS=1E-7 !有關pivot精度，視必要可適確調整 CALL MINVS(H,NT,NT,EPS,ILL) IF(ILL.NE.0) WRITE(*,10) ILL 10 FORMAT(1X,'ILL MINVS.1 ='I7) CALL SEKI(H,G,NT,NT,NT,NT) RETURN END SUBROUTINE HG1(H1,G1,XI,YI,XJ2,YJ2,XJ,YJ,SJ2) IF(SJ2.EQ.0) GOTO 100 PAI=3.1415927 Z=.57735 X1=(XJ2+XJ)/2 X2=(XJ2-XJ)/2 Y1=(YJ2+YJ)/2 Y2=(YJ2-YJ)/2 XE1=X1+X2*(-Z) YE1=Y1+Y2*(-Z) XE2=X1+X2*Z YE2=Y1+Y2*Z X11=XE1-XI Y11=YE1-YI X22=XE2-XI Y22=YE2-YI RM1=SQRT(X11**2+Y11**2) RM2=SQRT(X22**2+Y22**2) A1=XI*2*Y2-YI*2*X2 A2=YJ*2*X2-XJ*2*Y2 AA=ABS(A1+A2) BB=SQRT((2*Y2)**2+(2*X2)**2) RRP=AA/BB C1=(XJ-XI)*(YJ2-YI) D1=(XJ2-XI)*(YJ-YI) IF(C1.LT.D1) RRP=-RRP H1=-SJ2/4/PAI*(1.57735*RRP/RM1**2+.42265*RRP/RM2**2) G1=-SJ2/4/PAI*(1.57735*ALOG(RM1)+.42265*ALOG(RM2)) GOTO 101 100 H1=0 G1=0 101 CONTINUE RETURN END SUBROUTINE HG2(H2,G2,XI,YI,XJ,YJ,XJ1,YJ1,SJ1) IF(SJ1.EQ.0) GOTO 100 PAI=3.1415927 Z=.57735 X1=(XJ+XJ1)/2 X2=(XJ-XJ1)/2 Y1=(YJ+YJ1)/2 Y2=(YJ-YJ1)/2 XE1=X1+X2*(-Z) YE1=Y1+Y2*(-Z) XE2=X1+X2*Z YE2=Y1+Y2*Z X11=XE1-XI Y11=YE1-YI X22=XE2-XI Y22=YE2-YI RM1=SQRT(X11**2+Y11**2) RM2=SQRT(X22**2+Y22**2) A1=XI*2*Y2-YI*2*X2 A2=YJ1*2*X2-XJ1*2*Y2 AA=ABS(A1+A2) BB=SQRT((2*Y2)**2+(2*X2)**2) RRP=AA/BB C1=(XJ1-XI)*(YJ-YI) D1=(XJ-XI)*(YJ1-YI) IF(C1.LT.D1) RRP=-RRP H2=-SJ1/4/PAI*(0.42265*RRP/RM1**2+1.57735*RRP/RM2**2) G2=-SJ1/4/PAI*(0.42265*ALOG(RM1)+1.57735*ALOG(RM2)) GOTO 101 100 H2=0 G2=0 101 CONTINUE RETURN END SUBROUTINE CMINVS(A,B,M,N) ! 解複數連立方程式(應用逆矩陣)A+iB ! 原數據依實數部虛數部序，存檔於15 ! 結果：實數矩陣存於A，實數矩陣存於B ! 14：作業檔 REAL::A(M,M),B(M,M) IM=(0.,1.) EPS=1E-7 !有關pivot精度，視必要可適確調整 REWIND 15 READ(15) ((A(I,J),J=1,N),I=1,N) CALL MINVS(A,M,N,EPS,ILL) IF(ILL.NE.0) WRITE(*,10) ILL 10 FORMAT(1X,'ILL CS.1 ='I7) READ(15) ((B(I,J),J=1,N),I=1,N) CALL SEKI(A,B,M,N,N,N) REWIND 14 WRITE(14) ((A(I,J),J=1,N),I=1,N) CALL SEKI(B,A,M,N,N,N) REWIND 15 READ(15) ((A(I,J),J=1,N),I=1,N) DO I=1,N DO J=1,N A(I,J)=A(I,J)+B(I,J) END DO END DO CALL MINVS(A,M,N,EPS,ILL) IF(ILL.NE.0) WRITE(*,11) ILL 11 FORMAT(1X,'ILL CS.2 =',I7) REWIND 14 READ(14) ((B(I,J),J=1,N),I=1,N) CALL SEKI(B,A,M,N,N,N) DO I=1,N DO J=1,N B(I,J)=-B(I,J) END DO END DO RETURN END SUBROUTINE SHG_KH(SHG,HK) ! **利用牛頓近似法，已知無因次週頻率shg從分散關係式求 kh** IF((SHG-10.).LE.0) GOTO 2 XY=SHG GO TO 6 2 IF((SHG-1.0).GE.0) GOTO 4 X=SQRT(SHG) GO TO 5 4 X=SHG 5 COTHX=1./TANH(X) XY=X-(X-SHG*COTHX)/(1.+SHG*(COTHX**2-1.)) E=1.-XY/X X=XY IF((ABS(E)-0.0005).GE.0) GOTO 5 6 HK=XY RETURN END SUBROUTINE REGULAR(IT,U,E,DF,DT) ! 規則波造波 U=E*SIN(DF*DT*(IT-1)) RETURN END SUBROUTINE SOLITON(IT,U,XO,OMEGA,TC,DDT) ! 孤立波造波 U=-XO*OMEGA/COSH(OMEGA*(DDT*IT-TC))**2 RETURN END SUBROUTINE FREQ_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM) ! 單方向不規則波造波 REAL RNUM(NSPEC),S0(NSPEC),RE(NSPEC) PAI2=6.283185308 TEMP=0 DO J=1,NSPEC SIGMAT=PAI2*DF*(J-1) AA=RE(J)*SQRT(2.*DF*S0(J)) TEMP=TEMP+AA*SIGMAT*COS(SIGMAT*(IT-1)*DDT-RNUM(J)) END DO U=TEMP RETURN END SUBROUTINE BM0(X,N,M,H3,T3,DT) ! BRETSCHNEIDER-MITSUYASU SPECTRUM ! X : BRETSCHNEIDER-MITSUYASU FREQUENCY SPECTRUM ! N : DATA NUMBER OF SPECTRUM ! H3 : SIGNIFICANT WAVE HEIGHT ! T3 : SIGNIFICANT WAVE PERIOD ! DT : TIME INTERVAL REAL X(M) DF=1.0/(2*N*DT) A=0.257*H3**2 B=T3**(-4) C=-1.03*B AB=A*B X(1)=0.0 DO I=2,N F0=DF*(I-1) CFB=F0**(-4) CF=C*CFB IF(CF .LT. -25.0) GOTO 2 X(I)=AB*CFB/F0*EXP(CF) GOTO 1 2 X(I)=0. 1 END DO RETURN END SUBROUTINE BMS(X,N,M,H3,T3,DT) ! 修正 BRETSCHNEIDER-MITSUYASU SPECTRUM ! X : 修正 BRETSCHNEIDER-MITSUYASU FREQUENCY SPECTRUM ! N : DATA NUMBER OF SPECTRUM ! H3 : SIGNIFICANT WAVE HEIGHT ! T3 : SIGNIFICANT WAVE PERIOD ! DT : TIME INTERVAL DIMENSION X(M) DF=1.0/(2*N*DT) A=0.205*H3**2 B=T3**(-4) C=-0.75*B AB=A*B X(1)=0.0 DO I=2,N F0=DF*(I-1) CFB=F0**(-4) CF=C*CFB IF(CF .LT. -25.0) GOTO 2 X(I)=AB*CFB/F0*EXP(CF) GOTO 1 2 X(I)=0. 1 END DO RETURN END SUBROUTINE JONSWAP(X,N,M,H3,T3,DT) ! X : JONSWAP FREQUENCY SPECTRUM ! N : DATA NUMBER OF SPECTRUM ! H3 : SIGNIFICANT WAVE HEIGHT ! T3 : SIGNIFICANT WAVE PERIOD ! DT : TIME INTERVAL DIMENSION X(M) DF=1.0/(2*N*DT) GANMA=3.3 ALPHA=0.0624/(0.23+0.0336*GANMA-0.185*(1.9+GANMA)**(-1)) BETA=ALPHA*(1.094-0.01915*LOG(GANMA)) TPEAK=T3/(1-0.132*(GANMA+0.2)**(-0.559)) FPEAK=1/(TPEAK) C1J=BETA*H3**2*TPEAK**(-4) C2J=-1.25*TPEAK**(-4) X(1)=0.0 DO I=2,N F=DF*(I-1) IF (F.LT.FPEAK) SIGMA=0.07 IF (F.GE.FPEAK) SIGMA=0.09 F1=C2J*F**(-4) F2=-(TPEAK*F-1.00)**2/(2.00*SIGMA**2) IF (F2.LT.-25.0) THEN F0=1.0 ELSE F0=GANMA**EXP(F2) END IF IF (F1.LT.-25.0) THEN X(I)=0.0 ELSE X(I)=C1J*F**(-5)*EXP(F1)*F0 END IF END DO RETURN END SUBROUTINE SPECIN(X,N,M,DT,E,F,FPEAK,SPEAK) REAL::X(M),E(25),F(25) ! 注意：係數隨案而異 F(1)=0;E(1)=0 DO I=2,25 F(I)=0.25*(I-1) END DO X(1)=0.0 DO I=2,N DF=1.0/(2*N*DT) DFT=DF*(I-1) IF (DFT.GE.5.0) GOTO 11 IF (DFT.GE.F(1) .AND. DFT.LT.F(2)) THEN XF= ABS((DFT-F(1))/(F(2)-F(1))) X(I)=(1.0-XF)*E(1)+XF*E(2) GOTO 20 ELSEIF (DFT.GE.F(23) .AND. DFT.LT.F(24)) THEN XF=ABS((DFT-F(23))/(F(24)-F(23))) X(I)=(1.0-XF)*E(23)+XF*E(24) GOTO 20 ELSEIF (DFT.GE.F(24) .AND. DFT.LT.F(25)) THEN XF=(DFT-F(24))/(F(25)-F(24)) X(I)=(1.0-XF)*E(24)+XF*E(25) GOTO 20 ENDIF DO J=2,22 IF (DFT.GE.F(J) .AND. DFT.LT.F(J+1)) THEN XF=(DFT-F(J))/(F(J+1)-F(J)) X(I)=(1.0-XF)*E(J)+XF*E(J+1) GOTO 20 ENDIF END DO 11 X(I)=0.0 20 END DO DO IP=2,N DF=1.0/(2*N*DT) DFT=DF*(IP-1) IF (DFT.GE.FPEAK) GOTO 221 END DO 221 IPEAK=IP DELTAF=SPEAK-X(IPEAK) NHS=IPEAK-20 NHE=IPEAK+20 IF (NHS.LE.1) NHS=2 IF (NHE.GE.N) NHE=N DO IH=NHS,NHE DF=1.0/(2*N*DT) DFT=DF*(IH-1) XDELTA=(-((DFT-FPEAK)/0.25)**2+1.0)*DELTAF X(IH)=X(IH)+XDELTA END DO RETURN END SUBROUTINE WVMOD(X,N,M,D,DT,TMAX,TMIN,ITYPE,HKK) ! 造波板特性函數活塞式或拍拉式 ! X : 特性函數 ! N : DATA NUMBER ! D : WATER DEPTH (M) ! DT : TIME INTERVAL (SEC) ! TMAX : MAXIMUM PERIOD FOR WAVE GENERATOR (SEC) ! TMIN : MINIMUM PERIOD FOR WAVE GENERATOR (SEC) ! ITYPE: 造波板型式 1--->活塞式 ; 2--->拍拉式 REAL::KH,X(M),HKK(M) PAI=3.14159265 DF=1.0/(2*N*DT) X(1)=0.0 HKK(1)=0 DO I=2,N T=1.0/(DF*(I-1)) IF(T .GT. TMAX) GOTO 5 IF(T .LT. TMIN) GOTO 5 KH=2.*PAI*D/WAVEL(T,D) HKK(I)=KH IF (KH.GT.40.0) THEN X(I)=X(I-1) GOTO 1 ENDIF AA=SINH(KH) BB=COSH(KH) IF(ITYPE.EQ.2) GOTO 50 X(I)=(AA*BB+KH)/AA**2 GOTO 1 50 X(I)=KH*(AA*BB+KH)/(AA*(1.0-BB+KH*AA)) GOTO 1 5 X(I)=0.0 1 END DO RETURN END FUNCTION WAVEL(T,D) !已知週期、水深從分散關係式求波長 ! CALCULATION OF WAVE LENGTH (M) ! T : WAVE PERIOD (SEC) ! D : WATER DEPTH (M) PAI2=6.283185 DD=PAI2*D/(9.8*T**2/PAI2) X=DD IF(DD .LT. 1.0) X=SQRT(DD) 4 COTH=1.0/TANH(X) XX=X-(X-DD*COTH)/(1.0+DD*(COTH**2-1.0)) EMIN=1.0-XX/X X=XX IF(ABS(EMIN) .GE. 0.0005) GOTO 4 WAVEL=PAI2/X RETURN END SUBROUTINE RANDU(IX,IY,FL) IY=MOD(12869*IX+6925,32768) FL=FLOAT(IY)/2.0**15 RETURN END SUBROUTINE SEKI(A,B,M,NA,NK,NB) ! 矩陣相乘 ! [A],[B]宣告(m, m) ! [A]=[A][B] ! (na,nb)=(na,nk)x(nk,nb) ! [A]資料被取代 REAL:: A(M,M),B(M,M) REAL:: C(NK) DO I=1,NA DO K=1,NK C(K)=A(I,K) END DO DO J=1,NB R=0 DO K=1,NK R=R+C(K)*B(K,J) END DO A(I,J)=R END DO END DO WRITE(*,*) 'SEKI' RETURN END SUBROUTINE MINVS(A,KO,N,EPS,ILL) REAL::A(KO,KO) INTEGER:: X(KO) LOGICAL B ILL=0 IF((KO.GE.N).AND.(N.GE.2).AND.(EPS.GT.0.0)) GOTO 1 ILL=30000 RETURN 1 DO I=1,N X(I)=I END DO DO K=1,N M=K W=0.0 DO I=K,N ABSA=ABS(A(I,K)) IF(ABSA.LE.W) GOTO 20 W=ABSA M=I 20 END DO IF(M.EQ.K) GOTO 50 L=X(M) X(M)=X(K) X(K)=L DO J=1,N W=A(K,J) A(K,J)=A(M,J) A(M,J)=W END DO 50 IF(ABS(A(K,K)).GE.EPS) GOTO 60 ILL=K RETURN 60 P=1.0/A(K,K) DO J=1,N A(K,J)=A(K,J)*P END DO DO I=1,N T=-A(I,K) B=(I.NE.K).AND.(T.NE.0.0) IF(.NOT.B) GOTO 100 DO J=1,N A(I,J)=A(I,J)+A(K,J)*T END DO 100 A(I,K)=P*T END DO A(K,K)=P END DO DO I=1,N 120 IF(X(I).EQ.I) GOTO 140 K=X(I) DO J=1,N W=A(J,I) A(J,I)=A(J,K) A(J,K)=W END DO L=X(I) X(I)=X(K) X(K)=L GOTO 120 140 END DO WRITE(*,*) 'MINVS' RETURN END ! SUBROUTINE WA(A,B,NA,N,NS) ! [A]N,NS + [B]N,NS =====> [B]N,NS REAL::A(NA,NA),B(NA,NA) DO I=1,N DO J=1,NS B(I,J)=A(I,J)+B(I,J) END DO END DO WRITE(*,*) 'WA' RETURN END ! SUBROUTINE SA(A,B,NA,N,NS) ! [A]N,NS - [B]N,NS =====> [B]N,NS REAL::A(NA,NA),B(NA,NA) DO I=1,N DO J=1,NS B(I,J)=A(I,J)-B(I,J) END DO END DO WRITE(*,*) 'SA' RETURN END