!PROGRAM 3DHARBOR_MULTPILE

!  本程式尚未驗証

!   任意地形及港形考量岸壁反射率海底摩擦效應潮流效應的多方向不規則波港池波高分布(頻率領域)
!   網格元素依「應用Delaunay三角分割建置3維海域各邊界網格元素」自動生成
!   http://www.xn--06q65y2zay3jryqs5a29hepz.tw/wave_html_VB/19ModelGen/ModelGen.html
!   將「W3GOE.IN」複製至執行資料夾(自行建立)
!   理論分析參考「頻率領域空間3維問題」
!   http://www.xn--06q65y2zay3jryqs5a29hepz.tw/bem_html/bem/bem_8.0.pdf
!   參考數值解析用相關不規則波譜公式摘要
!   各離岸堤元素數預設為200，可自行變更

!       周宗仁  5. 6.2024   改版GFORTRAN  增加海底摩擦效應，可增潮流效應
!               5.26.2024  更正
        REAL,ALLOCATABLE:: XX(:),YY(:),ZZ(:),S(:),X(:,:),Y(:,:),Z(:,:)
        REAL,ALLOCATABLE:: XN(:),YN(:),ZN(:),GB(:,:),G(:,:)
        REAL,ALLOCATABLE:: XM(:),YM(:),SM(:),XNM(:),YNM(:)
        REAL,ALLOCATABLE:: XMM(:),YMM(:),ZMM(:),GG(:)
        REAL,ALLOCATABLE:: KD(:,:,:),U(:),V(:),SK(:),KKH(:),KKT(:)
        REAL,ALLOCATABLE:: CC(:,:),KR(:),FR(:),Fcn(:),Fn(:)
        REAL,ALLOCATABLE:: SPEC(:),SPEC0(:)
        COMPLEX,ALLOCATABLE:: AFA1(:),AFA2(:)
        COMPLEX,ALLOCATABLE:: H(:,:),HB(:,:),FO(:),FOB(:)

        INTEGER,ALLOCATABLE:: X_LEFT(:),Y_LEFT(:),Z_LEFT(:)
        INTEGER,ALLOCATABLE:: VVLEFT(:,:),LEFT_JNB(:),LEFT_NEI(:,:)
        INTEGER,ALLOCATABLE:: X_START(:) , Y_START(:) , Z_START(:)
        INTEGER,ALLOCATABLE:: VVSTART(:,:),START_JNB(:),START_NEI(:,:)
        INTEGER,ALLOCATABLE:: X_RIGHT(:) , Y_RIGHT(:) , Z_RIGHT(:)
        INTEGER,ALLOCATABLE:: VVRIGHT(:, :),RIGHT_JNB(:),RIGHT_NEI(:,:)
        INTEGER,ALLOCATABLE:: X_HARBOR(:) , Y_HARBOR(:) , Z_HARBOR(:)
        INTEGER,ALLOCATABLE:: HARBOR_JNB(:),HARBOR_NEI(:,:)
        INTEGER,ALLOCATABLE:: VVHARBOR(:,:)
        INTEGER,ALLOCATABLE:: X_SURFACE(:) ,Y_SURFACE(:) ,Z_SURFACE(:)
        INTEGER,ALLOCATABLE:: SURFACE_JNB(:),SURFACE_NEI(:,:)
        INTEGER,ALLOCATABLE:: VVSURFACE(:,:)
        INTEGER,ALLOCATABLE:: X_BED(:) ,Y_BED(:) ,Z_BED(:)
        INTEGER,ALLOCATABLE:: VVBED(:,:),BED_JNB(:),BED_NEI(:,:)
        INTEGER,ALLOCATABLE:: X_ISLAND(:,:),Y_ISLAND(:,:),Z_ISLAND(:,:)
        INTEGER,ALLOCATABLE:: ISLAND_N(:) ,ISLAND_E(:),VVISLAND(:,:,:)
        INTEGER,ALLOCATABLE:: ISLAND_JNB(:,:) ,ISLAND_NEI(:,:,:)

        INTEGER::  LEFT_N ,LEFT_E,O_NUM,M_NUM
        INTEGER::  START_N ,START_E
        INTEGER::  RIGHT_N ,RIGHT_E
        INTEGER::  HARBOR_N ,HARBOR_E
        INTEGER::  SURFACE_N ,SURFACE_E
        INTEGER::  BED_N ,BED_E
        COMPLEX::  RS,IM,TEMP
        REAL::  KH,MAIN_THETA

      OPEN(UNIT=21,FILE='W3.IN',FORM='FORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=25,FILE='W3GOE.IN',FORM='FORMATTED',STATUS='UNKNOWN')
! 流速分佈(預設=0)，即暫不考量
      OPEN(UNIT=23,FILE='CURRENT.W3',FORM='FORMATTED',STATUS='UNKNOWN')
! OUTPUT用
      OPEN(UNIT=22,FILE='W3.OUT',FORM='FORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=28,FILE='COO.W3',FORM='FORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=29,FILE='KD.OUT',FORM='FORMATTED',STATUS='UNKNOWN')
! 以下 UNIT 預留為作業空間，勿佔用，以免造成錯誤
      OPEN(UNIT=1,FILE='1.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=2,FILE='2.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=3,FILE='3.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=4,FILE='4.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=8,FILE='8.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=9,FILE='9.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      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')
      OPEN(UNIT=30,FILE='30.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=31,FILE='31.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=32,FILE='32.APT',FORM='UNFORMATTED',STATUS='UNKNOWN')

!  ************  地理資訊  W3GOE.IN  **********************

        READ(25,*) Ho                              ! 外海等水深領域水深
        READ(25,*) NDZ                             ! 假想邊界面水深方向層數(預設=2)
        READ(25,*) M1                              ! 左假想邊界線
        READ(25,*) M2                              ! 造波假想邊界線
        READ(25,*) M3                              ! 右假想邊界線
        READ(25,*) N1                              ! 假想邊界面=(M1+M2+M3)*NDZ
        READ(25,*) N2                              ! 自由水面
        READ(25,*) N3                              ! 岸壁等
        READ(25,*) N4                              ! 海底面
        READ(25,*) NIN                             ! 離岸堤數
        READ(25,*) UV_INDEX                        ! 有無潮流 無： 0

        WRITE(*,*) 'Ho=',Ho,'NDZ=',NDZ
        WRITE(*,*) 'M1=',M1,'M2=',M2,'M3=',M3
        WRITE(*,*) 'N1=',N1,'N2=',N2,'N3=',N3,'N4=',N4

        N12=N1/NDZ
        M=M1+M2+M3
        N=N1+N2
        NS=N3+N4
        NT=N1+N2+N3+N4    !  總元素數
        NX=MAX(N,NS)      !  作業陣列大小，取大值
        WRITE(*,*) 'N=',N,'   NS=',NS

        ALLOCATE(XX(NT),YY(NT),ZZ(NT),S(NT),X(NT,4),Y(NT,4),Z(NT,4))
        ALLOCATE(XN(NT),YN(NT),ZN(NT),GB(NX,NX),G(NX,NX))
        ALLOCATE(XM(N12),YM(N12),SM(N12),XNM(N12),YNM(N12))
        ALLOCATE(XMM(N1),YMM(N1),ZMM(N1))
        ALLOCATE(U(N2),V(N2),SK(N2),GG(N2))
        ALLOCATE(CC(NX,8),KR(N3),FR(N4))
        ALLOCATE(AFA1(N3),AFA2(N4))
        ALLOCATE(H(N1,N1),HB(N1,N1),FO(N12),FOB(N12))

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

        DO I=1,N3
        READ(25,*) KR(I)
        END DO
        DO I=1,N4
        READ(25,*) FR(I)
        END DO
        DO  I=1,N3
        AFA1(I)=IM*KH*SQRT(1-KR(I)**2)            !岸壁反射率相關係數
	END DO
	DO I=1,N4
	AFA2(I)=IM*FR(I)                          !海底摩擦相關係數
	END DO

        IF (UV_INDEX.EQ.0) THEN                    !無流場
        DO I=1,N2
        U(I)=0
        V(I)=0
        END DO
        END IF
        IF (UV_INDEX.NE.0) THEN                    !有流場CURRENT.W3載入
        DO I=1,N2
        READ(23,*) U(I),V(I)
        END DO
        END IF

        DO I=1,N1
        READ(25,*) XMM(I),YMM(I),ZMM(I)              ! 假想邊界線格點座標
        END DO

        READ(25,*) LEFT_N,LEFT_E                           ! LEFT Node,Element
        ALLOCATE(X_LEFT(LEFT_N),Y_LEFT(LEFT_N),Z_LEFT(LEFT_N))
        ALLOCATE(VVLEFT(LEFT_E,4))
        ALLOCATE(LEFT_JNB(LEFT_N),LEFT_NEI(LEFT_N,10))
        DO I = 1,LEFT_N
        READ(25,*)  X_LEFT(I),Y_LEFT(I),Z_LEFT(I)
        END DO
        DO I = 1,LEFT_N
        READ(25,*)  LEFT_JNB(I)
        END DO
        DO I = 1,LEFT_N
        DO J=1,LEFT_JNB(I)
        READ(25,*)  LEFT_NEI(I,J)
        END DO
        END DO
        DO I = 1,LEFT_E
        DO J= 1,4
        READ(25,*)  VVLEFT(I, J)
        END DO
        END DO

        READ(25,*) START_N,START_E                     !   START Node,Element
        ALLOCATE(X_START(START_N) , Y_START(START_N) )
        ALLOCATE(VVSTART(START_E, 4), Z_START(START_N))
        ALLOCATE(START_JNB(START_N),START_NEI(START_N,10))
        DO I = 1,START_N
        READ(25,*) X_START(I),Y_START(I),Z_START(I)
        END DO
        DO I = 1,START_N
        READ(25,*)  START_JNB(I)
        END DO
        DO I = 1,START_N
        DO J=1,START_JNB(I)
        READ(25,*)  START_NEI(I,J)
        END DO
        END DO
        DO I = 1,START_E
        DO J= 1,4
        READ(25,*) VVSTART(I, J)
        END DO
        END DO

        READ(25,*)  RIGHT_N , RIGHT_E                   !  RIGHT Node,Element
        ALLOCATE(X_RIGHT(RIGHT_N),Y_RIGHT(RIGHT_N))
        ALLOCATE(VVRIGHT(RIGHT_E, 4),Z_RIGHT(RIGHT_N))
        ALLOCATE(RIGHT_JNB(RIGHT_N),RIGHT_NEI(RIGHT_N,10))
        DO I = 1,RIGHT_N
        READ(25,*)  X_RIGHT(I) , Y_RIGHT(I) , Z_RIGHT(I)
        END DO
        DO I = 1,RIGHT_N
        READ(25,*)  RIGHT_JNB(I)
        END DO
        DO I = 1,RIGHT_N
        DO J=1,RIGHT_JNB(I)
        READ(25,*)  RIGHT_NEI(I,J)
        END DO
        END DO
        DO I = 1,RIGHT_E
        DO J= 1,4
        READ(25,*) VVRIGHT(I,J)
        END DO
        END DO

        READ(25,*)  HARBOR_N , HARBOR_E                         !HARBOR Node,Element
        ALLOCATE(X_HARBOR(HARBOR_N) , Y_HARBOR(HARBOR_N))
        ALLOCATE(Z_HARBOR(HARBOR_N),VVHARBOR(HARBOR_E,4))
        ALLOCATE(HARBOR_JNB(HARBOR_N),HARBOR_NEI(HARBOR_N,10))
        DO I = 1,HARBOR_N
        READ(25,*)  X_HARBOR(I),Y_HARBOR(I),Z_HARBOR(I)
        END DO
        DO I = 1,HARBOR_N
        READ(25,*)  HARBOR_JNB(I)
        END DO
        DO I = 1,HARBOR_N
        DO J=1,HARBOR_JNB(I)
        READ(25,*)  HARBOR_NEI(I,J)
        END DO
        END DO
        DO I = 1,HARBOR_E
        DO J= 1,4
        READ(25,*) VVHARBOR(I,J)
        END DO
        END DO

        READ(25,*)  SURFACE_N,SURFACE_E                  !SURFACE Node,Element
        ALLOCATE(X_SURFACE(SURFACE_N) ,Y_SURFACE(SURFACE_N))
        ALLOCATE(Z_SURFACE(SURFACE_N),VVSURFACE(SURFACE_E,4))
        ALLOCATE(SURFACE_JNB(SURFACE_N),SURFACE_NEI(SURFACE_N,10))
        DO I = 1,SURFACE_N
        READ(25,*)  X_SURFACE(I),Y_SURFACE(I), Z_SURFACE(I)
        END DO 
        DO I = 1,SURFACE_N
        READ(25,*)  SURFACE_JNB(I)
        END DO
        DO I = 1,SURFACE_N
        DO J=1,SURFACE_JNB(I)
        READ(25,*)  SURFACE_NEI(I,J)
        END DO
        END DO               
        DO I = 1,SURFACE_E
        DO J= 1,4
        READ(25,*)  VVSURFACE(I,J)
        END DO
        END DO

        READ(25,*)  BED_N , BED_E                      !BED Node,Element
        ALLOCATE(X_BED(BED_N) ,Y_BED(BED_N) ,Z_BED(BED_N))
        ALLOCATE(VVBED(BED_E,4))
        ALLOCATE(BED_JNB(BED_N),BED_NEI(BED_N,10))
        DO I = 1,BED_N
        READ(25,*)  X_BED(I) , Y_BED(I) , Z_BED(I)
        END DO
        DO I = 1,BED_N
        READ(25,*)  BED_JNB(I)
        END DO
        DO I = 1,BED_N
        DO J=1,BED_JNB(I)
        READ(25,*)  BED_NEI(I,J)
        END DO
        END DO        
        DO I = 1,BED_E
        DO J= 1,4
        READ(25,*) VVBED(I,J)
        END DO
        END DO

        If (NIN.NE.0) Then                                 !  有離岸堤時
        ALLOCATE(ISLAND_N(NIN) , ISLAND_E(NIN))
        ALLOCATE(X_ISLAND(NIN,200),Y_ISLAND(NIN,200))
        ALLOCATE(Z_ISLAND(NIN,200),VVISLAND(NIN,200,4))
        ALLOCATE(ISLAND_JNB(10,200) ,ISLAND_NEI(10,200,4))
        DO II = 1,NIN
        READ(25,*)  ISLAND_N(II),ISLAND_E(II)   !  ISLAND Node,Element
        DO I = 1,ISLAND_N(II)
        READ(25,*)  X_ISLAND(II,I),Y_ISLAND(II,I),Z_ISLAND(II,I)
        END DO
        DO I = 1,ISLAND_N(II)
        READ(25,*)  ISLAND_JNB(II,I)
        END DO
        DO I = 1,ISLAND_N(II)
        DO J=1,ISLAND_JNB(II,I)
        READ(25,*)  ISLAND_NEI(II,I,J)
        END DO
        END DO        
        DO I = 1,ISLAND_E(II)
        DO J= 1,4
        READ(25,*) VVISLAND(II,I,J)
        END DO
        END DO
        END DO
        End If

       !                   為數值計算，取得各1次元素的4頂角X,Y,Z座標

        DO II=1,LEFT_E
        DO J=1,4
        X(II,J)=X_LEFT(VVLEFT(II,J))
        Y(II,J)=Y_LEFT(VVLEFT(II,J))
        Z(II,J)=Z_LEFT(VVLEFT(II,J))
        END DO
        END DO
        DO II=1,START_E
        I=II+LEFT_E
        DO J=1,4
        X(I,J)=X_START(VVSTART(II,J))
        Y(I,J)=Y_START(VVSTART(II,J))
        Z(I,J)=Z_START(VVSTART(II,J))
        END DO
        END DO
        DO II=1,RIGHT_E
        I=II+LEFT_E+START_E
        DO J=1,4
        X(I,J)=X_RIGHT(VVRIGHT(II,J))
        Y(I,J)=Y_RIGHT(VVRIGHT(II,J))
        Z(I,J)=Z_RIGHT(VVRIGHT(II,J))
        END DO
        END DO
        DO II=1,HARBOR_E
        I=II+LEFT_E+START_E+RIGHT_E
        DO J=1,4
        X(I,J)=X_HARBOR(VVHARBOR(II,J))
        Y(I,J)=Y_HARBOR(VVHARBOR(II,J))
        Z(I,J)=Z_HARBOR(VVHARBOR(II,J))
        END DO
        END DO
        DO II=1,SURFACE_E
        I=II+LEFT_E+START_E+RIGHT_E+HARBOR_E
        DO J=1,4
        X(I,J)=X_SURFACE(VVSURFACE(II,J))
        Y(I,J)=Y_SURFACE(VVSURFACE(II,J))
        Z(I,J)=Z_SURFACE(VVSURFACE(II,J))
        END DO
        END DO
        DO II=1,BED_E
        I=II+LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E
        DO J=1,4
        X(I,J)=X_BED(VVBED(II,J))
        Y(I,J)=Y_BED(VVBED(II,J))
        Z(I,J)=Z_BED(VVBED(II,J))
        END DO
        END DO

        IF (NIN.NE.0) THEN
        DO III=1,NIN
        IV=0
        DO II=1,ISLAND_E(III)
        IV=IV+ISLAND_E(III)
        I=IV+LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E+BED_E
        DO J=1,4
        X(I,J)=X_ISLAND(III,VVISLAND(III,II,J))
        Y(I,J)=Y_ISLAND(III,VVISLAND(III,II,J))
        Z(I,J)=Z_ISLAND(III,VVISLAND(III,II,J))
        END DO
        END DO
        END DO
        END IF
        IF (NIN.EQ.0) THEN
        NT=LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E+BED_E
        END IF
        IF (NIN.NE.0) THEN
        NT=I
        END IF
!******************    計算網格中點、法線方向值、面積    ***************
        WRITE(*,*) ' BEGIN XYZ  計算邊界面網格點法線'

        CALL XYZ(X,Y,Z,XN,YN,ZN,XX,YY,ZZ,S,NT,NT)
        WRITE(*,*) 'END OF XYZ'
        CALL XY(XMM,YMM,XM,YM,XNM,YNM,SM,NT,N12,NDZ)
        WRITE(*,*) 'END OF XY  計算假想邊界線網格點法線'

!******* 必要時變更法線方向，確認外掛軟體格式(預設為不變更)

        ICHANGE=1    ! 變更ICHANGE=0即可改變法線方向(向外為正)
        IF(ICHANGE.EQ.0) THEN
        DO I=1,N3         !  岸壁面
        II=I+N1+N2
        XN(II)=-XN(II)
        YN(II)=-YN(II)
        ZN(II)=-ZN(II)
        END DO
        END IF

	DO I=1,N4         !  海底面，因與水面反向
	II=N1+N2+N3+I
	ZN(II)=-ZN(II)
	END DO
! **********************  波浪資訊  **********

        WRITE(*,*) '輸入入射波有義週期(秒) T1/3=？？'
        READ(*,*) T13
        WRITE(*,*) '輸入入射波有義波高(m) H1/3=？？'
        READ(*,*) H13
        WRITE(*,*) '輸入波譜等能量分割數 M_NUM=？？'
        READ(*,*) M_NUM
        WRITE(*,*) 'T1/3=',T13,'  H1/3=',H13,'  M_NUM=',M_NUM
        WRITE(*,*)  '輸入主入射波角度(度)'
        READ(*,*) MAIN_THETA
        WRITE(*,*)   '輸入方向函數分割數'
        READ(*,*) O_NUM

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

        ALLOCATE(Fcn(M_NUM),Fn(M_NUM),SPEC(M_NUM),SPEC0(M_NUM))
        ALLOCATE(KKH(N2),KKT(N2),KD(M_NUM,O_NUM,N2))

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

    !  方向函數等角分割
        L=INT(0.11*SMAX)     ! L大於2
        L2=2*L
        KG=1
        KK=1
        LL2=L2
  100   KG=KG*LL2
        KK=KK*(LL2-1)
        LL2=LL2-2
        IF (LL2.GT.0) GO TO 100
    !   方向函數
        DOI=180/(O_NUM-1)
        DO I=1,O_NUM-1
        THETA=MAIN_THETA+DOI/2+DOI*(I-1)-90
        IF(THETA.LE.0.OR.THETA.GE.180) THEN 
        GG(I)=0
        ELSE
        theta = theta * pai / 180
        GG(I)=FLOAT(KG)/PAI/KK*COS(THETA)**L2
        END IF    
        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) !因次

        SHG=4*PAI*PAI*Ho*F*F/9.8    !數值計算用有因次
        CALL SHG_KH(SHG,KH)
       
        WRITE(*,*) 'T1/3=',T13, '  H1/3',H13, '  主入射波角度=',MAIN_THETA

        DO IMO=1,O_NUM-1     !波向回路 

        DOI=180/(O_NUM-1)        
        OO=DOI/2+DOI*(IMO-1)
        OI=OO/180.*PAI

!  入射波勢函數及法線方向導函數

        DO I=1,N12
        TEMP=CEXP(-IM*KH*(XM(I)*COS(OI)+YM(I)*SIN(OI)))
        FO(I)=-IM*TEMP
        FOB(I)=-TEMP*(COS(OI)*YNM(I)-SIN(OI)*XNM(I))*KH
        END DO
        WRITE(*,*) ' END OF FO'

!   考量流存在時波頻率變化

        DO I=1,N2
        SK(I)=SHG                                                     !不考量流(預設)
        IF(UV_INDEX.NE.0) THEN
        SK(I)=(SQRT(SHG*9.8)-KH*COS(OI)*U(I)-KH*SIN(OI)*V(I))**2/9.8  !  考量流
        END IF
        END DO

        WRITE(28,6546) (XX(I),YY(I),ZZ(I),XN(I),YN(I),ZN(I),I=1,NT)
6546    FORMAT(6F10.4)

!    邊界積分方程式

        CALL G3_4(GB,G,NX,NT,N,NS,XX,YY,ZZ,XN,YN,ZN,S,10,11)  !2x2部分矩陣3維邊界積分方程式
        WRITE(*,*) 'END OF G3_4'      !分成4小陣列存檔

        CALL G2(GB,G,NX,N12,XM,YM,SM,XNM,YNM,KH,2,H,HB,N1,FO,FOB,N12)  !假想邊界線2維邊界積分方程式
        WRITE(*,*) 'END OF G2'               !存於HB及H   HB及H為複數陣列

        CALL RQ(GB,G,NX,N1,N12,Z,ZZ,NT,KH)   ! 外海領域假想邊界面(3D)水深函數
        WRITE(*,*) 'END OF RQ'                ! GB→R    G→Q

        C=KH/SINH(KH)/0.5/(1+2*KH/SINH(2.*KH))

!   建立連立方程式(本部分程式依解析問題有不同)       實數部存於GB，虛數部存於G

! ----------------- c[R][K*][Q]

        DO I=1,N1
        DO J=1,N12
        RS=0
        DO K=1,N12
        RS=RS+GB(I,K)*H(K,J)
        END DO
        HB(I,J)=RS                 !   [R][K*]
        END DO
        END DO

        DO I=1,N1
        DO J=1,N1
        RS=0
        DO K=1,N12
        RS=RS+HB(I,K)*G(K,J)
        END DO
        H(I,J)=RS*C                 !   c[R][K*][Q]
        END DO
        END DO
!--------------------------連立方程式右邊已知項
        REWIND 31

        DO I=1,NX
        CC(I,1)=0                    !歸零
        CC(I,2)=0
        CC(I,3)=0
        CC(I,4)=0
        END DO

        DO I=1,N1
        RS=0
        DO J=1,N12
        RS=RS+GB(I,J)*FO(J)
        END DO
        CC(I,1)=REAL(RS)              !    [R][Fo]
        CC(I,3)=AIMAG(RS)
        END DO

        DO I=1,N1
        RS=0
        DO J=1,N12
        RS=RS+HB(I,J)*FOB(J)           !   [R][K*][FBo]
        END DO
        CC(I,1)=CC(I,1)-REAL(RS)       !   [R]{[Fo]-[K*][FBo]}
        CC(I,3)=CC(I,3)-AIMAG(RS)
        END DO

        WRITE(31) (CC(I,1),I=1,N)      !  實數部分
        WRITE(31) (CC(I,2),I=1,NS)     !  實數部分 = 0
        WRITE(31) (CC(I,3),I=1,N)      !  虛數部分
        WRITE(31) (CC(I,4),I=1,NS)     !  虛數部分 = 0

!         部分陣列分割如下
!                   11(NxN)  | 12(NxNS)
!                  ---------------------
!                   21(NSxN) | 22(NSxNS)

!                                 建立11部分陣列係數  實數部分
        REWIND 10
        WRITE(*,*) '1'
        READ(10) ((GB(I,J),J=1,N),I=1,N)             !  讀取[K11]

        DO J=1,N1
        DO I=1,N1
        GB(I,J)=GB(I,J)-REAL(H(I,J))                 !  [K11]-c[R][K*][Q]
        END DO
        END DO

        DO J=1,N2
        JJ=J+N1
        DO I=1,N
        GB(I,JJ)=SK(J)*GB(I,JJ)                      !  SHG[K12]
        END DO
        DO I=1,N2
        II=I+N1
        IF(II.EQ.JJ) GB(II,JJ)=GB(II,JJ)-1           !  SHG[K22]-I
        END DO
        END DO

        REWIND 1
        REWIND 2
        WRITE(1) ((GB(I,J),J=1,N),I=1,N)

!   建立[K12]部分陣列係數  實數部分
        DO I=1,N
        DO J=1,NS
        G(I,J)=0
        END DO
        END DO

        WRITE(1) ((G(I,J),J=1,NS),I=1,N)

!   建立[K21]部分陣列係數  實數部分
        WRITE(*,*) '2'
        READ(10) ((GB(I,J),J=1,NS),I=1,N)      !讀取[K12]，跳過
        WRITE(*,*) '3'
        READ(10) ((GB(I,J),J=1,N),I=1,NS)      !讀取[K21]

        DO I=1,NS
        DO J=1,N1
        G(I,J)=GB(I,J)
        END DO
        END DO
        DO I=1,NS
        DO J=1,N2
        JJ=J+N1
        G(I,JJ)=GB(I,JJ)*SK(J)
        END DO
        END DO

        WRITE(1) ((G(I,J),J=1,N),I=1,NS)

!   建立[K22]部分陣列係數  實數部分
        DO I=1,NS
        DO J=1,NS
        G(I,J)=0
        IF(I.EQ.J) G(I,J)=-1
        END DO
        END DO

        WRITE(1) ((G(I,J),J=1,NS),I=1,NS)
!   建立[K11]部分陣列係數  虛數部分
        DO I=1,N
        DO J=1,N
        G(I,J)=0                !  歸零
        END DO
        END DO

        DO I=1,N1
        DO J=1,N1
        G(I,J)=-AIMAG(H(I,J))
        END DO
        END DO

        WRITE(2)  ((G(I,J),J=1,N),I=1,N)

!   建立[K12]部分陣列係數  虛數部分
        BACKSPACE 10
        BACKSPACE 10
        WRITE(*,*) '4'
        READ(10) ((GB(I,J),J=1,NS),I=1,N)  !讀取[K12]部分陣列

        DO I=1,N
        DO J=1,N3
        G(I,J)=KH*SQRT(1.-KR(J)**2)*GB(I,J)
        END DO
        END DO
        DO I=1,N
        DO J=1,N4
        JJ=J+N3
        G(I,JJ)=KH*FR(J)*GB(I,JJ)
        END DO
        END DO
        WRITE(2) ((G(I,J),J=1,NS),I=1,N)

!   建立[K21]部分陣列係數  虛數部分
        DO I=1,NS
        DO J=1,N
        G(I,J)=0                !  歸零
        END DO
        END DO

        WRITE(2) ((G(I,J),J=1,N),I=1,NS)

!    建立[K22]部分陣列係數  虛數部分

        WRITE(*,*) '5'
        READ(10) ((GB(I,J),J=1,N),I=1,NS)     !空讀
        WRITE(*,*) '6'
        READ(10) ((GB(I,J),J=1,NS),I=1,NS)

        DO I=1,NS
        DO J=1,N3
        G(I,J)=KH*SQRT(1.-KR(J)**2)*GB(I,J)
        END DO
        END DO
        DO I=1,NS
        DO J=1,N4
        JJ=J+N3
        G(I,JJ)=FR(J)*GB(I,JJ)
        END DO
        END DO
        WRITE(2) ((G(I,J),J=1,NS),I=1,NS)
!
        CALL CQ2(GB,G,CC,NX,N,NS)       ! 解連立方程式

        REWIND 32
        READ(32) (CC(I,5),I=1,N)
        READ(32) (CC(I,6),I=1,NS)
        READ(32) (CC(I,7),I=1,N)
        READ(32) (CC(I,8),I=1,NS)

        DO  I=1,N2
        II=I+N1
        KD(IT,IMO,I)=CABS(CMPLX(CC(II,5),CC(II,7)))*SQRT(SK(I)/SHG)
        END DO

        WRITE(22,1110) T,OO
        DO I=1,N2
        WRITE(22,1002) XX(I+N1),YY(I+N1),KD(IT,IMO,I),IT,IMO,I
        END DO

        END DO  !IMO
        END DO  !IT

1002    FORMAT(3F10.4,3I5)
1110    FORMAT('T=',F6.3,5X,'OO=',F6.3)

!  不規則波波高比及週期比

        DO I=1,N2

        HK=0
        TK=0
        DO IT=1,M_NUM-1
        DO IMO=1,O_NUM-1
        DF=Fn(IT+1)-Fn(IT)
        DOO=DOI/180*PAI
        DE=KD(IT,IMO,I)*KD(IT,IMO,I)*GG(IMO)*SPEC0(IT)*DF*DOO
        HK=HK+DE
        TK=TK+Fcn(IT)*Fcn(IT)*DE
        END DO
        END DO

        KKH(I)=4*SQRT(HK)
        KKT(I)=1./SQRT(SQRT(1.03*PAI))*SQRT(HK/TK)
        END DO
        DO I=1,N2
        WRITE(29,1003) XX(I+N1),YY(I+N1),KKH(I),KKT(I),I,T,MAIN_THETA
        END DO


1003    FORMAT(4F10.4,I5,2F10.4)
        STOP
        END !PROGRAM 3DHARBOR_MULTPILE


!****                   副程式                ****
!****
!
      SUBROUTINE XYZ(X,Y,Z,XN,YN,ZN,XX,YY,ZZ,S,M,N)

! X,Y,Z ： 各元素格點座標
! XN,YN,ZN : 一定元素中點法線方向值
! XX,YY,ZZ : 一定元素中點座標
! S ：一定元素面積
! M ：矩陣預設最大值
! N ：矩陣實際值

      REAL::X(M,4),Y(M,4),Z(M,4),XN(M),YN(M),ZN(M)
      REAL::S(M),YY(M),XX(M),ZZ(M)

!     一定元素中點法線方向值
      DO I=1,N
      R1=X(I,3)-X(I,1)
      R2=Y(I,3)-Y(I,1)
      R3=Z(I,3)-Z(I,1)
      R4=X(I,4)-X(I,2)
      R5=Y(I,4)-Y(I,2)
      R6=Z(I,4)-Z(I,2)
      R=SQRT((R5*R3-R6*R2)**2+(R6*R1-R4*R3)**2+(R4*R2-R5*R1)**2)
      XN(I)=(R5*R3-R6*R2)/R
      YN(I)=(R6*R1-R4*R3)/R
      ZN(I)=(R4*R2-R5*R1)/R

!     一定元素中點座標
      XX(I)=.25*(X(I,1)+X(I,2)+X(I,3)+X(I,4))
      YY(I)=.25*(Y(I,1)+Y(I,2)+Y(I,3)+Y(I,4))
      ZZ(I)=.25*(Z(I,1)+Z(I,2)+Z(I,3)+Z(I,4))

!     一定元素面積
      T1=X(I,2)-X(I,1)
      T2=Y(I,2)-Y(I,1)
      T3=Z(I,2)-Z(I,1)
      T4=X(I,4)-X(I,1)
      T5=Y(I,4)-Y(I,1)
      T6=Z(I,4)-Z(I,1)
      R=.5*SQRT((T2*R3-T3*R2)**2+(T3*R1-T1*R3)**2+(T1*R2-T2*R1)**2)
      T=.5*SQRT((R2*T6-R3*T5)**2+(R3*T4-R1*T6)**2+(R1*T5-R2*T4)**2)
      S(I)=R+T

      END DO
      RETURN
      END

!****
      SUBROUTINE XY(X,Y,XM,YM,XNM,YNM,SM,M,N,NDZ)

! XM ,YM : 假想邊界線線分中點座標
! XNM , YNM : 假想邊界線線分中點法線方向
! SM : 假想邊界線線分
! M ：矩陣預設最大值
! N ：矩陣實際值
! NDZ：假想邊界線水深方向層數(預設=2)
!
      REAL::X(M,4),Y(M,4),XM(N),YM(N),XNM(N),YNM(N),SM(N)
      DO I=1,N
      II=I*NDZ
      XM(I)=.5*(X(II,1)+X(II,3))
      YM(I)=.5*(Y(II,1)+Y(II,3))
      SM(I)=SQRT((X(II,1)-X(II,3))**2+(Y(II,1)-Y(II,3))**2)
      XNM(I)=(X(II,3)-X(II,1))/SM(I)
      YNM(I)=(Y(II,3)-Y(II,1))/SM(I)
      END DO
      RETURN
      END
!===============
      SUBROUTINE G3_4(GB,G,M,NT,N,NS,XX,YY,ZZ,XN,YN,ZN,S,L,L1)

!  2x2部份矩陣邊界積分方程式(一定元素)
!  將邊界積分方程式分割成2x2部份矩陣，減少矩陣容量並提高計算精度

! G ：基本解作業矩陣
! GB：基本解對法線方向導函作業矩陣
! XN,YN,ZN : 一定元素中點法線方向值
! XX,YY,ZZ : 一定元素中點座標
! M ：矩陣預設最大值
! NT：全部空間網格數=N1+N2+N3+N4
! N1：假想邊界面網格數
! N2：靜水面網格數
! N3：岸壁(含海岸及防波堤)網格數
! N4：海底面網格數
! N=N1+N2
! NS=N3+N4

! 分割部分矩陣

!     NxN  |  NxNS
!    ------|-------
!     NSxN |  NSxNS

      REAL::GB(M,M),G(M,M)
      REAL::XX(NT),YY(NT),ZZ(NT),XN(NT),YN(NT),ZN(NT),S(NT)

      REWIND L
      REWIND L1

      DO II=1,N
      I=II
      DO JJ=1,N
      J=JJ
      CALL G3G(GB,G,M,NT,XX,YY,ZZ,XN,YN,ZN,S,I,J,II,JJ)
      END DO
      END DO

      WRITE(*,*) '1'
      WRITE(L) ((GB(II,JJ),JJ=1,N),II=1,N)
      WRITE(L1) ((G(II,JJ),JJ=1,N),II=1,N)

      DO II=1,N
      I=II
      DO JJ=1,NS
      J=JJ+N
      CALL G3G(GB,G,M,NT,XX,YY,ZZ,XN,YN,ZN,S,I,J,II,JJ)
      END DO
      END DO

      WRITE(*,*) '2'
      WRITE(L) ((GB(II,JJ),JJ=1,NS),II=1,N)
      WRITE(L1) ((G(II,JJ),JJ=1,NS),II=1,N)

      DO II=1,NS
      I=II+N
      DO JJ=1,N
      J=JJ
      CALL G3G(GB,G,M,NT,XX,YY,ZZ,XN,YN,ZN,S,I,J,II,JJ)
      END DO
      END DO

      WRITE(*,*) '3'
      WRITE(L) ((GB(II,JJ),JJ=1,N),II=1,NS)
      WRITE(L1) ((G(II,JJ),JJ=1,N),II=1,NS)

      DO II=1,NS
      I=II+N
      DO JJ=1,NS
      J=JJ+N
      CALL G3G(GB,G,M,NT,XX,YY,ZZ,XN,YN,ZN,S,I,J,II,JJ)
      END DO
      END DO

      WRITE(*,*) '4'
      WRITE(L) ((GB(II,JJ),JJ=1,NS),II=1,NS)
      WRITE(L1) ((G(II,JJ),JJ=1,NS),II=1,NS)

      CALL MIR2(GB,G,M,N,NS,L,31)       ! 2x2部份矩陣逆矩陣
      CALL MTR2(GB,G,M,N,NS,31,L1,L)    ! 2x2部份矩陣相乘

      RETURN
      END
!===============
      SUBROUTINE G3G(GB,G,M,NT,XX,YY,ZZ,XN,YN,ZN,S,I,J,II,JJ)

!              建立各元素基本解及其法線方向導函數係數值(2x2矩陣用)

! G ：基本解作業矩陣
! GB：基本解對法線方向導函數作業矩陣
! XN,YN,ZN : 一定元素中點法線方向值
! XX,YY,ZZ : 一定元素中點座標
! S：一定元素面積
! M ：矩陣預設最大值
! NT ：全部空閭間網格數=N1+N2+N3+N4
! N1：假想邊界面網格數
! N2：靜水面網格數
! N3：岸壁(含海岸及防波堤)網格數
! N4：海底面網格數

      REAL::GB(M,M),G(M,M)
      REAL::XX(NT),YY(NT),ZZ(NT),XN(NT),YN(NT),ZN(NT),S(NT)
      PAI=3.1415927

      IF(I.EQ.J)  THEN
      GB(II,JJ)=1.0
      G(II,JJ)=SQRT(0.9945/PAI*S(J))
      END IF

      IF(I.NE.J)  THEN
      R=SQRT((XX(I)-XX(J))**2+(YY(I)-YY(J))**2+(ZZ(I)-ZZ(J))**2)
      RN=(XX(J)-XX(I))*XN(J)+(YY(J)-YY(I))*YN(J)+(ZZ(J)-ZZ(I))*ZN(J)
      RN=-RN/R
      GB(II,JJ)=0.5/PAI/R**2*RN*S(J)
      G(II,JJ)=0.5/PAI/R*S(J)
      END IF

      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
!
        SUBROUTINE MTR2 (A,B,NX,N,NS,L,M,K)
        REAL::A(NX,NX),B(NX,NX)
        WRITE(*,*) 'BEGIN MTR2'
        REWIND L
        REWIND M
        READ(L) ((A(I,J),J=1,N),I=1,N)
        READ(M) ((B(I,J),J=1,N),I=1,N)
        CALL SEKI(A,B,NX,N,N,N)
        REWIND 30
        WRITE(30) ((A(I,J),J=1,N),I=1,N)
        DO I=1,N
        DO J=1,N
        B(I,J)=A(I,J)
        END DO
        END DO
        BACKSPACE L
        READ(L) ((A(I,J),J=1,N),I=1,N)
        READ(M) ((B(I,J),J=1,NS),I=1,N)
        CALL SEKI (A,B,NX,N,N,NS)
        WRITE(30) ((A(I,J),J=1,NS),I=1,N)
        READ(L) ((A(I,J),J=1,NS),I=1,N)
        READ(M) ((B(I,J),J=1,N ),I=1,NS)
        CALL SEKI (A,B,NX,N,NS,N)
        REWIND 30
        READ(30) ((B(I,J),J=1,N),I=1,N)
        CALL WA (A,B,NX,N,N)
        REWIND K
        WRITE(K) ((B(I,J),J=1,N),I=1,N)
        READ(L) ((A(I,J),J=1,N ),I=1,NS)
        READ(L) ((A(I,J),J=1,NS),I=1,NS)
        REWIND L
        READ(L) ((A(I,J),J=1,N ),I=1,N)
        READ(L) ((A(I,J),J=1,NS),I=1,N)
        READ(M) ((B(I,J),J=1,NS),I=1,NS)
        CALL SEKI (A,B,NX,N,NS,NS)
        READ(30) ((B(I,J),J=1,NS),I=1,N)
        CALL WA (A,B,NX,N,NS)
        WRITE(K) ((B(I,J),J=1,NS),I=1,N)
        REWIND M
        READ(M) ((B(I,J),J=1,N),I=1,N)
        READ(L) ((A(I,J),J=1,N),I=1,NS)
        CALL SEKI (A,B,NX,NS,N,N)
        REWIND 30
        WRITE(30) ((A(I,J),J=1,N),I=1,NS)
        READ(L) ((B(I,J),J=1,NS),I=1,NS)
        BACKSPACE L
        BACKSPACE L
        READ(L) ((A(I,J),J=1,N),I=1,NS)
        READ(M) ((B(I,J),J=1,NS),I=1,N)
        CALL SEKI (A,B,NX,NS,N,NS)
        WRITE(30) ((A(I,J),J=1,NS),I=1,NS)
        READ(L) ((A(I,J),J=1,NS),I=1,NS)
        READ(M) ((B(I,J),J=1,N ),I=1,NS)
        CALL SEKI (A,B,NX,NS,NS,N)
        REWIND 30
        READ(30) ((B(I,J),J=1,N ),I=1,NS)
        CALL WA (A,B,NX,NS,N)
        WRITE(K) ((B(I,J),J=1,N ),I=1,NS)
        BACKSPACE L
        READ(L) ((A(I,J),J=1,NS),I=1,NS)
        READ(M) ((B(I,J),J=1,NS),I=1,NS)
        CALL SEKI (A,B,NX,NS,NS,NS)
        READ(30) ((B(I,J),J=1,NS),I=1,NS)
        CALL WA (B,A,NX,NS,NS)
        WRITE(K) ((A(I,J),J=1,NS),I=1,NS)
        WRITE(*,*) 'END OF MTR2'
        RETURN
        END
!
        SUBROUTINE MIR2(A,B,NX,N,NS,L,M)
!**************************************************************************
!     CALCULATE LARGE MATRIX A(N+NS,N+NS)'S INVERSE   [A1]N,N  [A2]N,NS   *
!  READ A FROM UNIT L FILE, AFTER CALCULATING,        [A3]NS,N [A4]NS,NS  *
!  WRITE,UNIT M FILE WITH :                                            *
!               [M1] = [ [A1] - [A2] [A4]^ [A3] ]^ ,  ^ MEANS INVERSE     *
!               [M3] = - [A4]^ [A3] [M1]                                  *
!               [M4] = [ [A4] - [A3] [A1]^ [A2] ]^                        *
!               [M2] = - [A1]^ [A2] [M4]                                  *
!**************************************************************************
        REAL::A(NX,NX),B(NX,NX)
        WRITE(*,*) 'BEGIN MIR2'
        REWIND L
        READ(L) ((A(I,J),J=1,N ),I=1,N)
        WRITE(*,*) '1'
        READ(L) ((A(I,J),J=1,NS),I=1,N)
        WRITE(*,*) '2'
        READ(L) ((A(I,J),J=1,N ),I=1,NS)
        WRITE(*,*) '3'
        READ(L) ((B(I,J),J=1,NS),I=1,NS)
        WRITE(*,*) '4'
        CALL MINVS(B,NX,NS,1E-7,ILL)
        ILLD=1
        IF(ILL.NE.0) WRITE (6,10) ILLD,ILL
        CALL SEKI(B,A,NX,NS,NS,N)
        REWIND 30
        WRITE(30) ((B(I,J),J=1,N),I=1,NS)
        BACKSPACE L
        BACKSPACE L
        BACKSPACE L
        READ(L) ((A(I,J),J=1,NS),I=1,N)
        CALL SEKI (A,B,NX,N,NS,N)
        READ(L) ((B(I,J),J=1,N ),I=1,NS)
        READ(L) ((B(I,J),J=1,NS),I=1,NS)
        REWIND L
        READ(L) ((B(I,J),J=1,N ),I=1,N )
        CALL SA (B,A,NX,N,N)
        CALL MINVS (A,NX,N,1E-7,ILL)
        ILLD=2
        IF(ILL.NE.0) WRITE (6,10) ILLD,ILL
        READ(L) ((B(I,J),J=1,NS),I=1,N)
        READ(L) ((B(I,J),J=1,N ),I=1,NS)
        READ(L) ((B(I,J),J=1,NS),I=1,NS)
        REWIND M
        WRITE(M) ((A(I,J),J=1,N ),I=1,N)
        REWIND 30
        READ(30) ((B(I,J),J=1,N ),I=1,NS)
        CALL SEKI (B,A,NX,NS,N,N)
        DO I=1,NS
        DO J=1,N
        B(I,J)=-B(I,J)
        END DO
        END DO
        REWIND 9
        WRITE(9) ((B(I,J),J=1,N ),I=1,NS)
        REWIND L
        READ(L) ((A(I,J),J=1,N ),I=1,N)
        CALL MINVS (A,NX,N,1E-7,ILL)
        ILLD=3
        IF(ILL.NE.0) WRITE(*,10) ILLD,ILL
        READ(L) ((B(I,J),J=1,NS),I=1,N )
        CALL SEKI (A,B,NX,N,N,NS)
        REWIND 30
        WRITE(30) ((A(I,J),J=1,NS),I=1,N)
        READ(L) ((B(I,J),J=1,N ),I=1,NS)
        CALL SEKI (B,A,NX,NS,N,NS)
        READ(L) ((A(I,J),J=1,NS),I=1,NS)
        CALL SA (A,B,NX,NS,NS)
        CALL MINVS (B,NX,NS,1E-7,ILL)
        ILLD=4
        IF(ILL.NE.0) WRITE(*,10) ILLD,ILL
        WRITE(9) ((B(I,J),J=1,NS),I=1,NS)
        REWIND 30
        READ(30) ((A(I,J),J=1,NS),I=1,N)
        CALL SEKI (A,B,NX,N,NS,NS)
        DO I=1,N
        DO J=1,NS
        A(I,J)=-A(I,J)
        END DO
        END DO
        WRITE(M) ((A(I,J),J=1,NS),I=1,N)
        REWIND 9
        READ(9) ((A(I,J),J=1,N ),I=1,NS)
        WRITE(M) ((A(I,J),J=1,N),I=1,NS)
        READ(9) ((A(I,J),J=1,NS),I=1,NS)
        WRITE(M) ((A(I,J),J=1,NS),I=1,NS)
10      FORMAT(1X,'ILL(',I1,')=',I7)
        WRITE(*,*) 'END OF MIR2'
        RETURN
        END

      SUBROUTINE CS(A,EB,M,N,L,F,FS,KI)

! **  DATA INPUT : UNIT 15
! **  TEMP. FILE : UNIT 14

      REAL::A(M,M),EB(M,M)
      COMPLEX::F(L),FS(L)
      COMPLEX IM,RS
      IM=(0.,1.)
      REWIND 15
      READ(15) ((A(I,J),J=1,N),I=1,N)
      CALL MINVS(A,M,N,1E-7,ILL)
      IF(ILL.NE.0) WRITE(*,10) ILL
10    FORMAT(1X,'ILL CS.1 ='I7)
      READ(15) ((EB(I,J),J=1,N),I=1,N)
      CALL SEKI(A,EB,M,N,N,N)
      REWIND 14
      WRITE(14)  ((A(I,J),J=1,N),I=1,N)
      CALL SEKI(EB,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)+EB(I,J)
      END DO
      END DO

      CALL MINVS(A,M,N,1E-7,ILL)
      IF(ILL.NE.0) WRITE(*,11) ILL
11    FORMAT(1X,'ILL CS.2 =',I7)
      REWIND 14
      READ(14) ((EB(I,J),J=1,N),I=1,N)
      CALL SEKI(EB,A,M,N,N,N)

      IF(KI.NE.1) THEN
      DO I=1,N
      RS=(0.,0.)
      DO J=1,N
      RS=RS+(A(I,J)-IM*EB(I,J))*F(J)
      END DO
      FS(I)=RS
      END DO
      END IF

      RETURN
      END

      SUBROUTINE G2(GB,G,M,N,XM,YM,SM,XN,YN,KH,KI,H,HB,NX,F1,F2,L)

!               2維邊界積分方程式

      REAL::GB(M,M),G(M,M),XM(N),YM(N),SM(N),XN(N),YN(N)
      COMPLEX::H(NX,NX),HB(NX,NX),F1(L),F2(L)
      REAL::KH
      COMPLEX::IM,RS,RT

      IM=(0.,1.)
      PAI=3.1415927
      DO I=1,N
      DO J=1,N
      IF(I.NE.J) THEN
      R=SQRT((XM(I)-XM(J))**2+(YM(I)-YM(J))**2)
      RK=KH*R
      CALL BESJNS(RK,0,BJO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
1     FORMAT(1X,'ILL G2 =',I7)
      CALL BESJNS(RK,1,BJ1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,0,BYO,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      CALL BESYNS(RK,1,BY1,ILL)
      IF(ILL.NE.0) WRITE(*,1) ILL
      RN=((XM(J)-XM(I))*YN(J)-(YM(J)-YM(I))*XN(J))/R
      RS=0.5*IM*KH*(BJ1+IM*BY1)*SM(J)*RN
      RT=-0.5*IM*(BJO+IM*BYO)*SM(J)/HO
      END IF
      IF(I.EQ.J) THEN
      IF(KI.EQ.1) RS=(-1.,0.)
      IF(KI.EQ.2) RS=(1.,0.)
      RT=(ALOG(KH*SM(J)/4.)-0.42278-IM*PAI/2.)/PAI*SM(J)
      END IF
      GB(I,J)=REAL(RS)
      G(I,J)=AIMAG(RS)
      HB(I,J)=RT
      END DO
      END DO

      REWIND 15
      WRITE(15) ((GB(I,J),J=1,N),I=1,N)
      WRITE(15) ((G(I,J),J=1,N),I=1,N)

      CALL CS(GB,G,M,N,L,F1,F2,1)

      DO I=1,N
      DO J=1,N
      RS=(0.,0.)
      DO K=1,N
      RS=RS+(GB(I,K)-IM*G(I,K))*HB(K,J)
      END DO
      H(I,J)=RS
      END DO
      END DO

      RETURN
      END

!************************************************************************
      SUBROUTINE BESJNS(X,N,BJN,ILL)

      BJN=0.0
      IF(N.GE.0) GOTO 20
      ILL=1
      RETURN
20    IF(X.LT.0) GOTO 30
      IF(X.EQ.0) GOTO 21
      IF(X.LT.0) GOTO 31
21    IF(N.GT.1) GOTO 30
      ILL=0
      IF(N.EQ.0) BJN=1.0
      RETURN
30    ILL=2
      RETURN
31    IF((X-15.0).GT.0) GOTO 34
      NTEST=20.0+10.0*X-X*X/3.0
      GOTO 36
34    NTEST=110.0+X/2.0
36    IF((N-NTEST).LT.0) GOTO 40
38    ILL=3
      RETURN
40    ILL=0
!
! -- *** COMPUTE STARTING VALUE OF M *** --
!
      N1=0
      IF(X.LT.1.0)   GOTO 54
      IF(X.LT.10.0)  GOTO 53
      IF(X.LT.50.0)  GOTO 52
      IF(X.LT.100.0) GOTO 51
      M=X+35.0
      GOTO 55
51    M=1.1*X+25.0
      GOTO 55
52    M=1.3*X+15.0
      N1=N+10
      GOTO 55
53    M=2.0*X+8.0
      N1=N+7
      GOTO 55
54    M=4.0*X+6.0
55    M=MAX0(M,N1)
!
! --- ** SET F(M),F(M-1) ** ---
!
      FM1=1.0E-30
      FM=0.0
      ALPHA=0.0
      IF((M-(M/2)*2).NE.0) GOTO 120
      JT=-1
      GOTO 130
120   JT=1
130   M2=M-2
      DO K=1,M2
      MK=M-K
      FMK=FLOAT(MK)
      BMK=2.0*FMK*FM1/X-FM
      FM=FM1
      FM1=BMK
      IF((MK-N-1).NE.0) GOTO 150
      BJN=BMK
150   JT=-JT
      S=1+JT
      ALPHA=ALPHA+BMK*S
      END DO
      BMK=2.0*FM1/X-FM
      IF(N.NE.0) GOTO 180
       BJN=BMK
180   ALPHA=ALPHA+BMK
      BJN=BJN/ALPHA
      RETURN
      END
!***************************************
      SUBROUTINE BESYNS(X,N,BYN,ILL)
! **************************************
      DOUBLE PRECISION TERM
      BYN=-1.0E30
      IF(N.LT.0) GOTO 180
      ILL=0
      IF(X.LT.0) GOTO 190
      IF(X.EQ.0) GOTO 11
      IF(X.GT.0) GOTO 20
11    IF(N.GT.1) GOTO 190
      RETURN
20    PI=3.141592653
      IF((X-4.0).LE.0) GOTO 40
      T=4.0/X
      T=T*T
      P0=((((-0.0000037043*T+0.0000173565)*T-0.0000487613)*T
     &  +0.0001734300)*T-0.0017530620)*T+0.3989422793

      Q0=((((0.0000032312*T-0.0000142078)*T+0.0000342468)*T
     &  -0.0000869791)*T+0.0004564324)*T-0.0124669441

      P1=((((0.0000042414*T-0.0000200920)*T+0.0000580759)*T
     &  -0.0002232030)*T+0.0029218256)*T+0.3989422819

      Q1=((((-0.0000036594*T+0.0000162200)*T-0.0000398708)*T
     &  +0.0001064741)*T-0.0006390400)*T+0.0374008364

      A=SQRT(2.0*PI)
      B=4.0*A
      P0=A*P0
      Q0=B*Q0/X
      P1=A*P1
      Q1=B*Q1/X
      A=X-PI/4.0
      B=SQRT(2.0/(PI*X))
      Y0=B*(P0*SIN(A)+Q0*COS(A))
      Y1=B*(-P1*COS(A)+Q1*SIN(A))
      GOTO 90
!
! --- ** COMPUTE Y0,Y1  DO X <= 4 ** ---
!
40    XX=X/2.0
      X2=XX*XX
      T=ALOG(XX)+0.5772156649
      SUM=0.0
      TERM=T
      Y0=T

      DO L=1,15
      IF((L-1).EQ.0) GOTO 60
      SUM=SUM+1.0/FLOAT(L-1)
60    FLL=L
      TS=T-SUM
      TERM=(TERM*(-X2)/(FLL*FLL))*(1.0-1.0/(FLL*TS))
      Y0=Y0+TERM
      END DO
      TERM=XX*(T-0.5)

      SUM=0.0
      Y1=TERM
      DO L=2,16
      SUM=SUM+1.0/FLOAT(L-1)
      FLL=L
      FL1=FLL-1.0
      TS=T-SUM
      TERM=(TERM*(-X2)/(FL1*FLL))*((TS-0.5/FLL)/(TS+0.5/FL1))
      Y1=Y1+TERM
      END DO


      PI2=2.0/PI
      Y0=PI2*Y0
      Y1=-PI2/X+PI2*Y1
!
! --- ** CHECK IF ONLY Y0 OR Y1 IS DESIRED ** ---
! ** ----------------------------------------- **
90    IF((N-1).GT.0) GOTO 130
! ** ----------------------------------------- **
! --- ** RETURN EITHER Y0 OR Y1 AS REQUIRED ** ---
! ** ------------------------------------------ **
      IF(N.EQ.0) GOTO 120
      BYN=Y1
      RETURN
120   BYN=Y0
      RETURN
! ----------------------------------------------  S.45.03.03
! --- ** PERFORM RECURRENCE OPERATIONS,FIND YN(X) ** ----
! ----------------------------------------------------------
130   YA=Y0
      YB=Y1
      K=1
140   T=FLOAT(2*K)/X
      YC=T*YB-YA
      K=K+1
      IF((K-N).EQ.0) GOTO 160
      YA=YB
      YB=YC
      GOTO 140
160   BYN=YC
      RETURN
180   ILL=1
      RETURN
190   ILL=2
      RETURN
      END
!**********************************************************************
        SUBROUTINE RQ(RR,QQ,NX,MR,M,Z,ZZ,MT,KH)
        REAL::RR(NX,NX),QQ(NX,NX),KH,Z(MT,4),ZZ(MT)

        NDZ=MR/M
        DO I=1,MR
        DO J=1,M
        RR(I,J)=0
        END DO
        END DO
        DO I=1,M
        DO J=1,MR
        QQ(I,J)=0
        END DO
        END DO
        DO I=1,M
        I2=(I-1)*NDZ
        DO I1=1,NDZ
        I3=I2+I1
        DO J=1,M
        IF (J.EQ.I) THEN
        RR(I3,J)=COSH(KH*(ZZ(I3)+1))/COSH(KH)
        END IF
        END DO
        END DO
        END DO
        DO J=1,M
        J2=(J-1)*NDZ
        DO J1=1,NDZ
        J3=J2+J1
        DO I=1,M
        IF (I.EQ.J) THEN
        DS=ABS(Z(J3,1)-Z(J3,3))
        QQ(I,J3)=COSH(KH*(ZZ(J3)+1))*DS
        END IF
        END DO
        END DO
        END DO

        RETURN
        END
!*******************************************************************
        SUBROUTINE MIC2(A,B,NX,N,NS)

!  2X2矩陣的複數逆矩陣

        REAL A(NX,NX),B(NX,NX)
        CALL MIR2(A,B,NX,N,NS,1,3)
        WRITE(*,*) 'MIR2'
        CALL MTR2(A,B,NX,N,NS,3,2,8)
        WRITE(*,*) 'MTR2'
        CALL MTR2(A,B,NX,N,NS,2,8,3)
        WRITE(*,*) 'MTR2'
        CALL WAR2(A,B,NX,N,NS,1,3,4)
        WRITE(*,*) 'WAR2'
        CALL MIR2(A,B,NX,N,NS,4,3)
        WRITE(*,*) 'MIR2'
        CALL MTR2(A,B,NX,N,NS,8,3,4)
        WRITE(*,*) 'MTR2'
        RETURN
        END
!*****************************************************
        SUBROUTINE WAR2(A,B,NX,N,NS,L,M,K)
!  2X2矩陣加法

        REAL::A(NX,NX),B(NX,NX)
        REWIND L
        REWIND M
        REWIND K
        READ(L) ((A(I,J),J=1,N),I=1,N)
        READ(M) ((B(I,J),J=1,N),I=1,N)
        CALL WA(B,A,NX,N,N)
        WRITE(K) ((A(I,J),J=1,N),I=1,N)
        READ(L) ((A(I,J),J=1,NS),I=1,N)
        READ(M) ((B(I,J),J=1,NS),I=1,N)
        CALL WA(B,A,NX,N,NS)
        WRITE(K) ((A(I,J),J=1,NS),I=1,N)
        READ(L) ((A(I,J),J=1,N),I=1,NS)
        READ(M) ((B(I,J),J=1,N),I=1,NS)
        CALL WA(B,A,NX,NS,N)
        WRITE(K) ((A(I,J),J=1,N),I=1,NS)
        READ(L) ((A(I,J),J=1,NS),I=1,NS)
        READ(M) ((B(I,J),J=1,NS),I=1,NS)
        CALL WA(B,A,NX,NS,NS)
        WRITE(K) ((A(I,J),J=1,NS),I=1,NS)
        RETURN
        END
!
!8**************************************************
        SUBROUTINE CQ2(A,B,C,NX,N,NS)

!  2X2矩陣的複數連立方程式

        REAL::A(NX,NX),B(NX,NX),C(NX,8)
        CALL MIC2(A,B,NX,N,NS)
        REWIND 31
        REWIND 32

        REWIND 3
        REWIND 4
        READ(3) ((A(I,J),J=1,N),I=1,N)
        READ(4) ((B(I,J),J=1,N),I=1,N)
        READ(31) (C(I,1),I=1,N)         ! 實數部分
        READ(31) (C(I,2),I=1,NS)        ! 實數部分
        READ(31) (C(I,3),I=1,N)         ! 虛數部分
        READ(31) (C(I,4),I=1,NS)        ! 虛數部分

        DO I=1,N
        R=0.
        S=0.
        T=0.
        U=0.
        DO J=1,N
        R=R+A(I,J)*C(J,1)
        S=S+B(I,J)*C(J,3)
        T=T+A(I,J)*C(J,3)
        U=U+B(I,J)*C(J,1)
        END DO
        C(I,5)=R+S
        C(I,7)=T-U
        END DO

        READ(3) ((A(I,J),J=1,NS),I=1,N)
        READ(4) ((B(I,J),J=1,NS),I=1,N)

        DO I=1,N
        U=0.
        T=0.
        S=0.
        R=0.
        DO J=1,NS
        R=R+A(I,J)*C(J,2)
        S=S+B(I,J)*C(J,4)
        T=T+A(I,J)*C(J,4)
        U=U+B(I,J)*C(J,2)
        END DO
        C(I,5)=C(I,5)+R+S
        C(I,7)=C(I,7)+T-U
        END DO

        READ(3) ((A(I,J),J=1,N),I=1,NS)
        READ(4) ((B(I,J),J=1,N),I=1,NS)

        DO I=1,NS
        T=0.
        U=0.
        R=0.
        S=0.
        DO J=1,N
        R=R+A(I,J)*C(J,1)
        S=S+B(I,J)*C(J,3)
        T=T+A(I,J)*C(J,3)
        U=U+B(I,J)*C(J,1)
        END DO
        C(I,6)=R+S
        C(I,8)=T-U
        END DO
        READ(3) ((A(I,J),J=1,NS),I=1,NS)
        READ(4) ((B(I,J),J=1,NS),I=1,NS)
        DO I=1,NS
        R=0.
        S=0.
        T=0.
        U=0.
        DO J=1,NS
        R=R+A(I,J)*C(J,2)
        S=S+B(I,J)*C(J,4)
        T=T+A(I,J)*C(J,4)
        U=U+B(I,J)*C(J,2)
        END DO
        C(I,6)=C(I,6)+R+S
        C(I,8)=C(I,8)+T-U
        END DO
        WRITE(32) (C(I,5),I=1,N)          ! 實數部分
        WRITE(32) (C(I,6),I=1,NS)         ! 實數部分
        WRITE(32) (C(I,7),I=1,N)          ! 虛數部分
        WRITE(32) (C(I,8),I=1,NS)         ! 虛數部分

        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