!PROGRAM SeaBed_IRE

!     頻率領域2維任意海底地形引起不規則波變形
!     本程式為公開程式，可不經授權逕行引用，但尚未實作驗証完成，有問題可電郵連絡。
!     本論文原始發表於1976年「グリ-ン公式による2次元水面波の境界值問題の解析」。
!      http://www.xn--06q65y2zay3jryqs5a29hepz.tw/chou_html/reference/1975ijima_01.pdf
!     追加海底摩擦效應及不規則波效應
!     定義圖參考「頻率領域2維任意海底地形引起波變形」
!     http://www.xn--06q65y2zay3jryqs5a29hepz.tw/bem_html/2D_SeaBottom/Seabed.pdf

!     本程式追加利用累加法分析不規則波作用時的水面波高及週期分佈!     周宗仁  2024.06.10

      REAL,ALLOCATABLE:: X(:),Z(:),S(:),XN(:),ZN(:),XX(:),ZZ(:)
      REAL,ALLOCATABLE:: H(:,:),G(:,:)
      REAL,ALLOCATABLE:: WH(:,:),KD(:),KT(:),KP(:),KS(:),KR(:)
      REAL,ALLOCATABLE:: SPEC(:),SPEC0(:),FCN(:),FNN(:)
      COMPLEX,ALLOCATABLE::F1(:,:),F2(:,:),ZO(:)
      COMPLEX,ALLOCATABLE::FA1(:),FA2(:),FA3(:),FA4(:)   !速度勢
      COMPLEX,ALLOCATABLE::FA1B(:),FA2B(:),FA3B(:),FA4B(:)  !速度勢導函數
      COMPLEX::RS,IM,TEMP,FF0
      REAL::KH1,KH2,L1,L2,NO1,NO2

      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=14,FILE='14.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=15,FILE='15.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')

      IM=(0.,1.)
      PAI=3.1415927

      WRITE(*,*)  '入射波領域等深域水深 H1=？ 預設為10'
      H1=10        ! 基準水深 
      WRITE(*,*)  '輸入入射波領域等深域水深 H2/H1=？'
      READ(*,*)  H2
      WRITE(*,*)  '輸入入射波等深域距原點距離 L1/H1=？'
      READ(*,*)  L1
      WRITE(*,*)  '輸入通過波等深域距原點距離 L2/H1=？'
      READ(*,*)  L2
      WRITE(*,*)  '輸入通過波側等深域間距 S1/H1=？'
      READ(*,*)  S1
      WRITE(*,*)  '輸入入射波等側深域間距 S3/H1=？'
      READ(*,*)  S3
      WRITE(*,*)  '輸入海底摩擦係數 AFA=？'
      READ(*,*)  AFA

      H2=H2*H1
      L1=L1*H1
      L2=L2*H1
      SI=S1*H1
      S3=S3*H1
      SL=L1+L2
      S2=SL-S1-S3

!  設定海底地形(可自行調整)採用一定元素逆時針積分，各邊界元素數如下

      WRITE(*,*)  '輸入入射波假想邊界元素數   N1=？ '
      READ(*,*) N1  
      WRITE(*,*)  '輸入任意地形水面邊界元素數 N2=？'
      READ(*,*) N2
      WRITE(*,*)  '輸入通過波假想邊界元素數   N3=？'
      READ(*,*) N3
      WRITE(*,*)  '輸入海底左等水深邊界元素數 N41=？'
      READ(*,*) N41
      WRITE(*,*)  '輸入任意海底邊界元素數     N42=？'
      READ(*,*) N42
      WRITE(*,*)  '輸入海底右等水深邊界元素數 N43=？'
      READ(*,*) N43

      N4=N41+N42+N43
      NT=N1+N2+N3+N4
      NT1=NT+1
 
      ALLOCATE(H(NT,NT),G(NT,NT),X(NT),Z(NT),S(NT),XN(NT),ZN(NT))
      ALLOCATE(F1(N1,N1),F2(N1,N1))
      ALLOCATE(XX(NT1),ZZ(NT1),ZO(NT))
      ALLOCATE(FA1(N1),FA2(N2),FA3(N3),FA4(N4))
      ALLOCATE(FA1B(N1),FA2B(N2),FA3B(N3),FA4B(N4))
!  設定元素座標及線段長(一定元素)，自行發揮

!  以入射波領域等深域邊界與海底交接處為起始(逆時針)

!  中點座標
      DO I=1,N1
      X(I)=L1
      Z(I)=-H1+H1/N1/2+H1/N1*(I-1)
      S(I)=H1/N1
      XN(I)=1
      ZN(I)=0
      END DO

      DO II=1,N2
      I=II+N1
      X(I)=L1-SL/N2/2-SL/N2*(II-1)
      Z(I)=0
      S(I)=SL/N2
      XN(I)=0
      ZN(I)=1
      END DO

      DO II=1,N3
      I=II+N1+N2
      X(I)=-L2
      Z(I)=-H2/N3/2-H2/N3*(II-1)
      S(I)=H3/N3
      XN(I)=-1
      ZN(I)=0
      END DO

      DO II=1,N41
      I=II+N1+N2+N3
      X(I)=-L2+S1/N41/2+S1/N41*(II-1)
      Z(I)=-H2
      S(I)=S1/N41
      XN(I)=0
      ZN(I)=-1
      END DO

      SS2=SQRT((H1-H2)**2+(S2)**2)
      DO II=1,N42
      I=II+N1+N2+N3+N41
      X(I)=-L2+S1+S2/N41/2+S2/N42*(II-1)
      Z(I)=-H2-(H1-H2)/N42/2-(H1-H2)/N42*(II-1)
      S(I)=SS2/N42
      XN(I)=-(H1-H2)/SS2
      ZN(I)=-S2/SS2
      END DO

      DO II=1,N43
      I=II+N1+N2+N3+N41+N42
      X(I)=-L2+S1+S2+S3/N43/2+S3/N43*(II-1)
      Z(I)=-H1
      S(I)=S3/N43
      XN(I)=0
      ZN(I)=-1
      END DO

!  端點座標(GUASS積分用)
      DO I=1,N1+1
      XX(I)=L1
      ZZ(I)=-H1+H1/N1*(I-1)
      END DO

      DO II=1,N2
      I=II+N1+1
      XX(I)=L1-SL/N2*II
      ZZ(I)=0
      END DO

      DO II=1,N3
      I=II+N1+1+N2
      XX(I)=-L2
      ZZ(I)=-H2/N3*II
      END DO

      DO II=1,N41
      I=II+N1+1+N2+N3
      XX(I)=-L2+S1/N41*II
      ZZ(I)=-H2
      END DO

      DO II=1,N42
      I=II+N1+1+N2+N3+N41
      XX(I)=-L2+S1+S2/N42*II
      ZZ(I)=-H2-(H1-H2)/N42*II
      END DO

      DO II=1,N43
      I=II+N1+1+N2+N3+N41+N42
      XX(I)=-L2+S1+S2+S3/N43*II
      ZZ(I)=-H1
      END DO

!  波浪資訊


!  只計算規則波時，將波譜等能量分割數 M_NUM設定為2即可

        WRITE(*,*) '輸入入射波有義週期(秒) T1/3=？？'
        READ(*,*) T13

        WRITE(*,*) '輸入入射波有義波高(m) H1/3=？？'
        READ(*,*) H13

        WRITE(*,*) '輸入波譜等能量分割數 M_NUM=？？'
        READ(*,*) M_NUM

        WRITE(*,*) '輸入集中度Smax =？'
        WRITE(*,*) '風波=10' 
        WRITE(*,*) '衰減距離短湧浪=25'
        WRITE(*,*) '衰減距離長湧浪=75' 
        READ(*,*) SMAX

        ALLOCATE(FCN(M_NUM),FNN(M_NUM),SPEC(M_NUM),SPEC0(M_NUM))
        ALLOCATE(KR(M_NUM),KP(M_NUM),WH(M_NUM,N2),KS(M_NUM))
        ALLOCATE(KD(N2),KT(N2))

        DO IT=1,M_NUM
        X2=M_NUM/IT
        FNN(IT)=1.0071*SQRT(SQRT(ALOG(X2))) !  入射波譜等能量分割
        END DO

!  開始計算

        DO IT=1,M_NUM-1   !週期回路

        IF(IT.NE.1) X1=M_NUM/(IT-1)
        IF(IT.EQ.1) X1=0

        X2=M_NUM/IT
        E1=ERF(SQRT(2*ALOG(X1)));E2=ERF(SQRT(2*ALOG(X2)))
        
        FCN(IT)=1/0.9*SQRT(2.9124*M_NUM*(E1-E2))          !無因次區間代表頻率
        SPEC0(IT)=0.257/FCN(IT)**5*EXP(-1.03/FCN(IT)**4)  !無因次

        F=FCN(IT)/T13                                     ! 因次區間代表頻率
        SPEC(IT)=0.257*H13*H13*T13/(T13*F)**5*EXP(-1.03/(T13*F)**4) !因次

!   牛頓近似法求KH
        SHG1=4*PAI*PAI*H1/T/T/9.8
        CALL SHG_KH(SHG1,KH1)

        SHG2=4*PAI*PAI*H2/T/T/9.8
        CALL SHG_KH(SHG2,KH2)

!  2維邊界積分方程式

      CALL G2_ALOG_CONST(H,G,X,Z,S,XX,ZZ,NT,NT1)   !結果存於H

!      REWIND 11
!      WRITE(11) ((H(I,J),J=1,NT),I=1,NT)  !備用

! 入射波速度勢

      DO  I= 1,NT
      ZO(I)=0       ! 歸零
      END DO
      DO  I= 1,N1
      ZO(I)=2*IM*KH1*COSH(KH1*(Z(I)+H1))/COSH(KH1)
      END DO

! 入射波邊界面速度勢與速度勢導函數間的關係式

      NO1=0.5*(1+2*KH1/SINH(2*KH1)) 
      DO I=1,N1
      DO J=1,N1
      F1(I,J)=IM*COSH(KH1*(Z(I)+H1))*COSH(KH1*(Z(J)+H1))
     & /NO1/SINH(KH1)/COSH(KH1)*S(J)
      END DO
      END DO

! 通過波邊界面速度勢與速度勢導函數間的關係式
 
      NO2=0.5*(1+2*KH2/SINH(2*KH2)) 
      DO I=1,N3
      II=I+N1+N2
      DO J=1,N3
      JJ=JJ+N1+N2
      F2(I,J)=IM*COSH(KH2*(Z(II)+H2))*COSH(KH2*(Z(JJ)+H2))
     & /NO2/SINH(KH2)/COSH(KH2)*S(JJ)
      END DO
      END DO

!  建立連立方程式

!  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)=H(I,J)-REAL(F1(II,JJ))
      G(I,J)=-AIMAG(F1(II,JJ))
      END DO
      DO JJ=1,N2
      J=JJ+N1
      H(I,J)=SHG1*H(I,J)
      END DO
      DO JJ=1,N4
      J=JJ+N1+N2+N3
      H(I,J)=0
      G(I,J)=AFA*H(I,J)
      END DO
      END DO
     
      DO II=1,N2
      I=II+N1
      DO JJ=1,N2
      J=JJ+N1
      H(I,J)=SHG1*H(I,J) 
      IF (I.EQ.J) H(I,J)= H(I,J)-1  
      END DO
      DO JJ=1,N4
      J=JJ+N1+N2+N3
      H(I,J)=0
      G(I,J)=AFA*H(I,J)
      END DO
      END DO
          
      DO II=1,N3
      I=II+N1+N2
      DO JJ=1,N2
      J=JJ+N1
      H(I,J)=SHG1*H(I,J) 
      END DO
      DO JJ=1,N3
      J=JJ+N1+N2
      H(I,J)=H(I,J)-REAL(F2(II,JJ))
      G(I,J)=-AIMAG(F2(II,JJ))    
      END DO
      DO JJ=1,N4
      J=JJ+N1+N2+N3
      H(I,J)=0
      G(I,J)=AFA*H(I,J)
      END DO
      END DO
     
      DO II=1,N4
      I=II+N1+N2+N3
      DO JJ=1,N2
      J=JJ+N1
      H(I,J)=SHG1*H(I,J)    
      END DO
      DO JJ=1,N4
      J=JJ+N1+N2+N3
      H(I,J)=0
      IF (I.EQ.J) H(I,J)=-1 
      G(I,J)=AFA*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,N1
      RS=(0.,0.) 
      DO J=1,NT
      RS=RS+CMPLX(H(I,J),G(I,J))*ZO(J)
      END DO
      FA1B(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))*ZO(J)
      END DO
      FA2(I)=RS
      END DO      

      DO II=1,N3
      I=II+N1+N2
      RS=(0.,0.) 
      DO J=1,NT
      RS=RS+CMPLX(H(I,J),G(I,J))*ZO(J)
      END DO
      FA3B(I)=RS
      END DO      

      DO II=1,N4
      I=II+N1+N2+N3
      RS=(0.,0.) 
      DO J=1,NT
      RS=RS+CMPLX(H(I,J),G(I,J))*ZO(J)
      END DO
      FA4(I)=RS
      END DO

      RS=0
      DO I=1,N1
      II=I+N1+N2
      RS=RS+FA1B(I)*COSH(Z(I)+H2)*S(I)
      END DO
      RS=IM/NO1/SINH(KH1)*RS
      KR(IT)=CABS(1+RS)                     !反射率

      RS=0
      DO I=1,N3
      II=I+N1+N2
      RS=RS+FA3B(I)*COSH(Z(II)+H2)*S(II)
      END DO
      RS=IM/NO2/SINH(KH2)*RS
      KR(IT)=CABS(RS)                       !透過率

      KS(IT)=SQRT(1-KR(IT)**2-KR(IT)**2)    !摩擦損失率


! 計算水面波高

      DO I=1,N2
      WH(IT,I)=CABS(FA2(I))
      END DO

      WRITE(1,2) KR(IT),KP(IT),KS(IT),T13 ,IT   
      WRITE(1,1) (WH(IT,I),I=1,N2)

      END DO  !IT

!  不規則波利用波譜疊加法

        A0=0
        B0=0
        C0=0
        T0=0
        DO IT=1,M_NUM-1
        DF=FNN(IT+1)-FNN(IT)
        A0=A0+KR(IT)*KR(IT)*SPEC0(IT)*DF
        B0=B0+KP(IT)*KP(IT)*SPEC0(IT)*DF
        C0=C0+KS(IT)*KS(IT)*SPEC0(IT)*DF
        END DO
        AA0=4*SQRT(A0)  !反射率
        BB0=4*SQRT(B0)  !透過率
        CC0=4*SQRT(C0)  !消耗率

        WRITE(1,3) AA0,BB0,CC0


        DO I=1,N2
        HK=0
        TK=0
        DO IT=1,M_NUM-1
        DF=FNN(IT+1)-FNN(IT)
        DE=WH(IT,I)*WH(IT,I)*SPEC0(IT)*DF
        HK=HK+DE
        TK=TK+FCN(IT)*FCN(IT)*DE
        END DO
        KD(I)=4*SQRT(HK)                           !水面波高比
        KT(I)=1./SQRT(SQRT(1.03*PAI))*SQRT(HK/TK)  !水面波高週期比
        END DO

1     FORMAT(10F5.2)
2     FORMAT(4F5.2,I3)
3     FORMAT('反射率= ',F5.2,'透過率= ',F5.2,'消耗率= ',F5.2)
1003  FORMAT(4F10.4,I5)      
      STOP  
      END

!             副程式

      SUBROUTINE G2_ALOG_CONST(H,G,X,Z,S,XX,ZZ,NT,NT1)

!     計算[K]=[H]^[G]

      REAL::X(NT),Z(NT),S(NT),H(NT,NT),G(NT,NT),XX(NT1),ZZ(NT1)

      PAI=3.1415927

      DO I=1,NT
      DO J=1,NT

      IF(I.EQ.J) THEN
      H(I,J)=0
      R=S(J)/2
      G(I,J)=R/PAI*(ALOG(1./R)+1)
      END IF

      IF(I.NE.J) THEN
        ZM=.57735
        X1=(XX(J+1)+XX(J))/2;Y1=(ZZ(J+1)+ZZ(J))/2
        X2=(XX(J+1)-XX(J))/2;Y2=(ZZ(J+1)-ZZ(J))/2
        XE1=X1+X2*(-ZM);YE1=Y1+Y2*(-ZM)
        XE2=X1+X2*ZM;   YE2=Y1+Y2*ZM
        X11=XE1-X(I);Y11=YE1-Z(I)
        X22=XE2-X(I);Y22=YE2-Z(I)
        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=(XX(J)-X(I))*(ZZ(J+1)-Z(I))
        D1=(XX(J+1)-X(I))*(ZZ(J)-Z(I))
        IF(C1.LT.D1) RRP=-RRP
        H(I,J)=-S(J)/4/PAI*(1.57735*RRP/RM1**2+.42265*RRP/RM2**2)   !GUASS積分，直接計算亦可
        G(I,J)=-S(J)/4/PAI*(1.57735*ALOG(RM1)+.42265*ALOG(RM2))
      END IF
      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 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 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