!PROGRAM 3DWAVEMAKING

!   任意地形及港形考量岸壁反射率海底摩擦效應的多方向不規則波港池造波(時間領域)

! 本程式為公開程式，可不經授權逕行引用，但尚未實作驗証完成，有問題可電郵連絡。
!   網格元素依「應用DelaunaY三角分割建置3維海域各邊界網格元素」自動生成
!   http://www.xn--06q65y2zay3jryqs5a29hepz.tw/wave_html_VB/19ModelGen/ModelGen.html
!   將「W3GOE_TIME.IN」複製至執行資料夾(自行建立)
!   理論分析參考「時間領域3維港池造波問題」
!   http://www.xn--06q65y2zay3jryqs5a29hepz.tw/bem_html/%E6%96%B03D%E9%80%A0%E6%B3%A2/%E6%96%B03D%E9%80%A0%E6%B3%A2.pdf
!   參考數值解析用相關不規則波譜公式摘要
!   各離岸堤元素數預設為200，可自行變更

!   周宗仁  6.18.2024   

        REAL,ALLOCATABLE:: XX(:),YY(:),ZZ(:),A(:),X(:,:),Y(:,:),Z(:,:)
        REAL,ALLOCATABLE:: XN(:),YN(:),ZN(:),GB(:,:),G(:,:)
        REAL,ALLOCATABLE:: KD(:,:,:),U(:),V(:),SK(:),KKH(:),KKT(:)
        REAL,ALLOCATABLE:: CC(:,:),KR(:),FR(:),Fcn(:),Fn(:)
        REAL,ALLOCATABLE:: SPEC(:),SPEC0(:)
        REAL,ALLOCATABLE:: RNUM(:),S0(:),HKK(:),RE(:),SS(:),THETA(:)
        REAL,ALLOCATABLE:: F1(:),F1B(:),F2B(:)
        COMPLEX,ALLOCATABLE:: AFA1(:),AFA2(:)
        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(:,:),JNB(:),NEI(:,:)
        INTEGER,ALLOCATABLE:: X_SURFACE(:) ,Y_SURFACE(:) ,Z_SURFACE(:)
        INTEGER,ALLOCATABLE:: SURFACE_JNB(:),SURFACE_NEI(:,:)
        INTEGER,ALLOCATABLE:: VVSURFACE(:,:),ELEMENT(:,:)
        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,RIGHT_N,RIGHT_E,HARBOR_N,HARBOR_E 
        INTEGER::  SURFACE_N ,SURFACE_E,BED_N,BED_E
        COMPLEX::  RS,IM,TEMP
        REAL::  KH,MAIN_THETA,E25(25),F25(25),XP(3,3),P(2,2),E(2,2)

      OPEN(UNIT=21,FILE='W3.IN',FORM='FORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=25,FILE='W3GOE_TIME.IN',
     & FORM='FORMATTED',STATUS='UNKNOWN')
      OPEN(UNIT=23,FILE='CURRENT.W3',FORM='FORMATTED',
     &  STATUS='UNKNOWN')              ! 流速分佈(預設=0)，即暫不考量
! 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_TIME.IN  **********************

        READ(25,*) H0                 ! 外海等水深領域水深
        READ(25,*) NDZ                ! 假想邊界面水深方向層數(預設=2)
        READ(25,*) N1                 ! 自由水面
        READ(25,*) N2                 ! 造波假想邊界面
        READ(25,*) N3                 ! 左假想邊界面
        READ(25,*) N4                 ! 右假想邊界面
        READ(25,*) N5                 ! 岸壁等
        READ(25,*) N6                 ! 海底面
        READ(25,*) NIN                ! 離岸堤數
        READ(25,*) UV_INDEX           ! 有無潮流 無：0
         
	ITYPE=1              !造波板型式
        NCH=N2-1             !造波板數
        NSPEC0=2048          !波譜分割數
        IM=(0.,1.) 
        GR=9.81
        PA2I=1.570796327 ; PAI =3.141592654 ; PAI2=6.283185308

        N=N1
        NS=N2+N3+N4+N5+N6
        NT=N+NS           !  總元素數
        NX=MAX(N,NS)      !  作業陣列大小，取大值

        ALLOCATE(GB(NX,NX),G(NX,NX))
        ALLOCATE(U(N1),V(N1),SK(N1))
        ALLOCATE(CC(NX,8),KR(N5),FR(N6),AFA1(N5),AFA2(N6))
        ALLOCATE(RNUM(NSPEC0),S0(NSPEC0))
        ALLOCATE(HKK(NSPEC0),RE(NSPEC0),SS(NSPEC0),THETA(NSPEC0))
        ALLOCATE(F1(N1),F1B(N1),F2B(NS))

        DO I=1,N5
        READ(25,*) KR(I)  !岸壁反射率，目前只適用規則波                
        END DO
        DO I=1,N6
        READ(25,*) FR(I)  !海底摩擦相關係數
        END DO
        DO I=1,N6
        FR(I)=0           !目前只適用無摩擦效應           
	AFA2(I)=IM*FR(I)                          
        END DO

        GAMMA=0           !假想邊界透過率，目前只適用透過率=0即全反射

        IF (UV_INDEX.EQ.0) THEN                    !無流場
        DO I=1,N1
        U(I)=0
        V(I)=0
        END DO
        END IF

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

        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

        TEMP=(X_START(START_N)-X_START(1))**2+
     &       (Y_START(START_N)-Y_START(1))**2
        TMP=SQRT(TEMP)       !造波機總長

        BLOD=TMP/NCH         !造波板寬

        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                                 !  有離岸堤時
        N_ISLAND=100    ! 各離岸堤最大元素數
        NMAX=10 ! 最多離岸堤數
        NJNB=5 ! 各離岸堤最大相鄰元素數
        ALLOCATE(ISLAND_N(NIN) , ISLAND_E(NIN))
        ALLOCATE(X_ISLAND(NIN,N_ISLAND),Y_ISLAND(NIN,N_ISLAND))
        ALLOCATE(Z_ISLAND(NIN,N_ISLAND),VVISLAND(NIN,N_ISLAND,4))
        ALLOCATE(ISLAND_JNB(NMAX,NJNB) ,ISLAND_NEI(NMAX,NJNB,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   !II
        END IF

!       1次元素的4頂角X,Y,Z座標

        NT_E=LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E+BED_E
        IF (NIN.NE.0) Then
        NT_E=LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E+BED_E+NMAX*NJNB
        END IF

        ALLOCATE(X(NT_E,4),Y(NT_E,4),Z(NT_E,4),A(NT_E))

        DO II=1,SURFACE_E
        DO J=1,4
        X(II,J)=X_SURFACE(VVSURFACE(II,J))
        Y(II,J)=Y_SURFACE(VVSURFACE(II,J))
        Z(II,J)=Z_SURFACE(VVSURFACE(II,J))
        END DO
        END DO
        DO II=1,START_E
        I=II+SURFACE_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,LEFT_E
        I=II+SURFACE_E+START_E
        DO J=1,4
        X(I,J)=X_LEFT(VVLEFT(II,J))
        Y(I,J)=Y_LEFT(VVLEFT(II,J))
        Z(I,J)=Z_LEFT(VVLEFT(II,J))
        END DO
        END DO
        DO II=1,RIGHT_E
        I=II+SURFACE_E+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+SURFACE_E+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,BED_E
        I=II+SURFACE_E+LEFT_E+START_E+RIGHT_E+HARBOR_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

        NT_E=I

        ALLOCATE(X(NT_E,4),Y(NT_E,4),Z(NT_E,4),A(NT_E))

!     格點座標

        NT_N=LEFT_N+START_N+RIGHT_N+HARBOR_N+SURFACE_N+BED_N
        IF (NIN.NE.0) THEN
        NT_N=LEFT_N+START_N+RIGHT_N+HARBOR_N+SURFACE_N+BED_N+NMAX*NJNB
        END IF
        ALLOCATE(XX(NT_N),YY(NT_N),ZZ(NT_N))

        DO II=1,SURFACE_N
        XX(II)=X_SURFACE(II)
        YY(II)=Y_SURFACE(II)
        ZZ(II)=Z_SURFACE(II)
        END DO
        DO II=1,START_N
        I=II+SURFACE_N
        XX(I)=X_START(II)
        YY(I)=Y_START(II)
        ZZ(I)=Z_START(II)
        END DO
        DO II=1,LEFT_N
        I=II+SURFACE_N+START_N
        XX(I)=X_LEFT(II)
        YY(I)=Y_LEFT(II)
        ZZ(I)=Z_LEFT(II)
        END DO
        DO II=1,RIGHT_N
        I=II+SURFACE_N+LEFT_N+START_N
        XX(I)=X_RIGHT(II)
        YY(I)=Y_RIGHT(II)
        ZZ(I)=Z_RIGHT(II)
        END DO
        DO II=1,HARBOR_N
        I=II+SURFACE_N+LEFT_N+START_N+RIGHT_N
        XX(I)=X_HARBOR(II)
        YY(I)=Y_HARBOR(II)
        ZZ(I)=Z_HARBOR(II)
        END DO
        DO II=1,BED_N
        I=II+SURFACN_E+LEFT_N+START_N+RIGHT_N+HARBOR_N
        XX(I)=X_BED(II)
        YY(I)=Y_BED(II)
        ZZ(I)=Z_BED(II)
        END DO

        IF (NIN.NE.0) THEN
        DO III=1,NIN
        IV=0
        DO II=1,ISLAND_N(III)
        IV=IV+ISLAND_N(III)
        I=IV+LEFT_N+START_N+RIGHT_N+HARBOR_N+SURFACE_N+BED_N
        XX(I)=X_ISLAND(III,II)
        YY(I)=Y_ISLAND(III,II)
        ZZ(I)=Z_ISLAND(III,II)
        END DO
        END DO
        END IF
        NT_N=I

        ALLOCATE(XX(NT_N),YY(NT_N),ZZ(NT_N))

!    整合各節點的相鄰元素數及編號

        ALLOCATE(JNB(NT_N),NEI(NT_N,10))

        DO II=1,SURFACE_N
        JNB(II)=SURFACE_JNB(II)
        DO J =1,JNB(II)
        NEI(II,J)=SURFACE_NEI(II,J)
        END DO
        END DO
        DO II=1,START_N
        I=II+SURFACE_N
        JNB(I)=START_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=START_NEI(II,J)
        END DO
        END DO
        DO II=1,LEFT_N
        I=II+SURFACE_N+START_N
        JNB(I)=LEFT_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=LEFT_NEI(II,J)
        END DO
        END DO
        DO II=1,RIGHT_N
        I=II+SURFACE_N+START_N+LEFT_N
        JNB(I)=RIGHT_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=RIGHT_NEI(II,J)
        END DO
        END DO
        DO II=1,HARBOR_N
        I=II+SURFACE_N+START_N+LEFT_N+RIGHT_N
        JNB(I)=HARBOR_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=HARBOR_NEI(II,J)
        END DO
        END DO
        DO II=1,BED_N
        I=II+SURFACE_N+START_N+LEFT_N+RIGHT_N+HARBOR_N
        JNB(I)=BED_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=BED_NEI(II,J)
        END DO
        END DO

        IF (NIN.NE.0) THEN
        DO III=1,NIN
        IV=0
        DO II=1,ISLAND_N(III)
        IV=IV+ISLAND_N(III)
        I=IV+LEFT_N+START_N+RIGHT_N+HARBOR_N+SURFACE_N+BED_N
        JNB(I)=BED_JNB(II)
        DO J =1,JNB(I)
        NEI(I,J)=ISLAND_NEI(III,II,J)
        END DO
        END DO
        END DO   !III
        END IF

!   水面節點與元素的關係，覓出各元素4節點的編號

        ALLOCATE(ELEMENT(SURFACE_E,4))

        DO J=1,SURFACE_E 
        DO K=1,4
        XE=X(J,K)
        YE=Y(J,K)
        DO I=1,SURFACE_N
        XNODE=XX(I)
        YNODE=YY(I)
        IF((XE.EQ.XNODE).AND.(YE.EQ.YNODE)) THEN
        ELEMENT(J,K)=I
        END IF
        END DO   !I
        END DO   !K
        END DO   !J

!******************    計算法線方向值、面積    ***************

! X,Y,Z ： 各元素格點座標
! XNN,YNN,ZNN : 元素中點法線方向值
! A ：元素面積
! NT ：格點最大值

      DO J=1,NT_E

!     起始法線方向值
      R1=X(J,3)-X(J,1)
      R2=Y(J,3)-Y(J,1)
      R3=Z(J,3)-Z(J,1)
      R4=X(J,4)-X(J,2)
      R5=Y(J,4)-Y(J,2)
      R6=Z(J,4)-Z(J,2)
      R=SQRT((R5*R3-R6*R2)**2+(R6*R1-R4*R3)**2+(R4*R2-R5*R1)**2)
      XN(J)=(R5*R3-R6*R2)/R
      YN(J)=(R6*R1-R4*R3)/R
      ZN(J)=(R4*R2-R5*R1)/R

!     元素面積
      T1=X(J,2)-X(J,1)
      T2=Y(J,2)-Y(J,1)
      T3=Z(J,2)-Z(J,1)
      T4=X(J,4)-X(J,1)
      T5=Y(J,4)-Y(J,1)
      T6=Z(J,4)-Z(J,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)
      A(J)=R+T

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

      ICHANGE=1    ! 變更ICHANGE=0即可改變法線方向(向外為正)

      IF(ICHANGE.EQ.0) THEN
      DO I=1,HARBOR_E         !  岸壁面
      II=I+SURFACE_E+LEFT_E+START_E+RIGHT_E
      XN(II)=-XN(II)
      YN(II)=-YN(II)
      ZN(II)=-ZN(II)
      END DO
      END IF

      DO J=1,BED_E
      JJ=J+LEFT_E+START_E+RIGHT_E+HARBOR_E+SURFACE_E         !  海底面，因與水面反向
      ZN(JJ)=-ZN(JJ)
      END DO
    
! **********************  波浪資訊  **********

        WRITE(*,*) '輸入造波形式 WAVE =1～5'
        WRITE(*,*) '1：規則波'
        WRITE(*,*) '2：短蜂波 '
        WRITE(*,*) '3：孤立波 '
        WRITE(*,*) '4：單方向不規則波'
        WRITE(*,*) '5：多方向不規則波 '
        READ(*,*)  WAVE
 
        IF(WAVE.EQ.1) THEN
        WRITE(*,*) '輸入入射波有義週期(秒) T31=？？'
        READ(*,*) T31
        WRITE(*,*) '輸入入射波有義波高(m) H31=？？'
        READ(*,*) H31
        WRITE(*,*)  '輸入主入射波角度(°) TETA0=？？'
        READ(*,*) TETA0
        WRITE(*,*)  '輸入時間間隔數 NTT=？？'
        READ(*,*) NTT
        HK=2*PAI/WAVEL(T31,H0)
	DT=T31/NTT			 
        RTETA0=TETA0*PAI/180.0
        DF=2*PAI/T31
        EE=-DF*H31*(SINH(HK)*COSH(HK)+HK)/2./SINH(HK)**2
        WRITE(*,*) 'T1/3=',T31,'  H1/3=',H31,'  TETA=',TETA0

        DO I=1,N5
        AFA1(I)=IM*KH*SQRT(1-KR(I)**2)            !岸壁反射率相關係數
        END DO
	END IF

        IF(WAVE.EQ.2) THEN
        WRITE(*,*) '輸入入射波有義週期(秒) T31=？？'
        READ(*,*) T31
        WRITE(*,*) '輸入入射波有義波高(m) H31=？？'
        READ(*,*) H31
        WRITE(*,*)  '輸入主入射波角度(°) TETA0=？？'
        READ(*,*) TETA01
        WRITE(*,*) '輸入入射波有義週期(秒) T32=？？'
        READ(*,*) T32
        WRITE(*,*) '輸入入射波有義波高(m) H32=？？'
        READ(*,*) H32
        WRITE(*,*)  '輸入主入射波角度(°) TETA02=？？'
        READ(*,*) TETA02
        WRITE(*,*)  '輸入時間間隔數 NTT=？？'
        READ(*,*) NTT

	DT=T31/NTT
        HK1=2*PAI/WAVEL(T31,H0)
        HK2=2*PAI/WAVEL(T32,H0)				 
        RTETA01=TETA01*PAI/180.0
        RTETA02=TETA02*PAI/180.0
        DF1=2*PAI/T31
        DF2=2*PAI/T32
        E1=-DF1*H31*(SINH(HK1)*COSH(HK1)+HK1)/2./SINH(HK1)**2
        E2=-DF2*H32*(SINH(HK2)*COSH(HK2)+HK2)/2./SINH(HK2)**2
        WRITE(*,*) 'T1=',T31,'  H1=',H31,'  TETA1=',TETA01
        WRITE(*,*) 'T2=',T32,'  H2=',H32,'  TETA2=',TETA02
        END IF
 
        IF(WAVE.EQ.3) THEN
        WRITE(*,*) '輸入孤立波波高(m) H=？？'
        READ(*,*) H
        WRITE(*,*) '輸入入射波角度(°) TETA0=？？'
        READ(*,*) TETA0
        WRITE(*,*)  '輸入時間間隔數 NTT=？？'
        READ(*,*) NTT
	DT=NTT			 
        RTETA0=TETA0*PAI/180.0
	GR1=GR*H0
	H=H/H0
	CE=SQRT(GR1*(1+H))
	XO=SQRT(4.*H/3./(1+H))
	OMEGA=SQRT(GR1)*SQRT(.75*H*(1+H))
	TC=PAI/OMEGA
	XLEFF=9.5766/SQRT(H)
	TEFF=XLEFF/CE
	DDT=TC/DT				 	 
        END IF

        IF(WAVE.EQ.4) 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(*,*)  '輸入主入射波角度(°) TETA0=？？'
        READ(*,*) TETA0
        WRITE(*,*)  '輸入時間間隔數 NTT=？？'
        READ(*,*) NTT
        WRITE(*,*)  '輸入波譜頻率分割數 NSPEC=？？'
        READ(*,*) NSPEC

        HK=2*PAI/WAVEL(T31,H0)
	DTT=NTT			 
        RTETA0=TETA0*PAI/180.0

	IF(ITPSPC.EQ.4) THEN
        WRITE(*,*)  'FPEAK=？？'
 	READ(*,*) FPEAK
        WRITE(*,*)  'SPEAK=？？'
	READ(*,*) SPEAK
        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,F25,
     &  FPEAK,SPEAK)

        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+HKK(J)*I*BLOD*(COS(RTETA0))
	END DO 
        END IF

        IF(WAVE.EQ.5) THEN
        WRITE(*,*) '波譜形式 ITPSPC(1～4)：？？'
        WRITE(*,*) '1：Bretschneider波譜'
        WRITE(*,*) '2：Bretschneider-MitsuYasu波譜'
        WRITE(*,*) '3： JONSWAP譜  4：自訂'

	READ(*,*) ITPSPC

        WRITE(*,*) '輸入入射波有義週期(秒) T31=？？'
        READ(*,*) T31
        WRITE(*,*) '輸入入射波有義波高(m) H31=？？'
        READ(*,*) H31
        WRITE(*,*)  '輸入主入射波角度(°) TETA0=？？'
        READ(*,*) TETA0
        WRITE(*,*)  '輸入方向集中度 SMAX？？'
        WRITE(*,*)  '風波          =10'
        WRITE(*,*)  '衰減距離短湧浪=25'
        WRITE(*,*)  '衰減距離長湧浪=75'
        READ(*,*) SMAX
        WRITE(*,*)  '輸入時間間隔數 NTT=？？'
        READ(*,*) NTT
        WRITE(*,*)  '輸入波譜頻率分割數 NSPEC=？？'
        READ(*,*) NSPEC

        HK=2*PAI/WAVEL(T31,H0)
	DTT=NTT			 
        RTETA0=TETA0*PAI/180.0

	IF(ITPSPC.EQ.4) THEN
        WRITE(*,*)  'FPEAK=？？'
 	READ(*,*) FPEAK
        WRITE(*,*)  'SPEAK=？'
	READ(*,*) SPEAK
        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,F25,
     &     FPEAK,SPEAK)

        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+HKK(J)*I*BLOD*(COS(RTETA0))
	END DO 
!========================= DIRECTION ========================
         L=(INT(0.11*SMAX+0.55))
         L2=2*L 
         KG=1
         KK=1
         LL2=L2
  100    KG=KG*LL2
         KK=KK*(LL2-1)
         LL2=LL2-2
         IF (LL2.GT.0) GOTO 100
         EXPR=PAI*FLOAT(KK)/FLOAT(KG)
         IRX=IXY

         DO J=1,NSPEC
  210    CALL RANDU(IRX,IRY,FL)
         IRX=IRY
         RCON=FL*EXPR
         THETAL=0.0
         THETAR=PAI
  200    THETAM=(THETAL+THETAR)/2.0
         THETAD=(THETAR-THETAL)
         GI=THETAM
         ALFA=THETAM-PA2I

         DO N=1,L
         GI=((SIN(ALFA)*(COS(ALFA))**(2*N-1))+(2*N-1)*GI)/(2*N)
         END DO

         DEF=GI-RCON
         IF(DEF .EQ. 0.0) GOTO 300
         IF(DEF .GT. 0.0) THETAR=THETAM
         IF(DEF .LT. 0.0) THETAL=THETAM
         IF(THETAD .GE.  0.0017) GOTO 200
  300    IF (TETA0 .LT. 90.0) THEN
         THETAM=THETAM-(PA2I-THETA0)
         IF (THETAM.LE.0.0 .OR. THETAM.GE.PA2I+THETA0) GOTO 210
         ELSE IF (TETA0 .GT. 90.0) THEN
         THETAM=THETAM-(PA2I-THETA0)
         IF (THETAM.LE.THETA0-PA2I .OR. THETAM.GE.PAI) GOTO 210
         END IF
         THETA(J)=THETAM
         END DO    ! J

	 DO J=1,NSPEC
         CALL RANDU(IXY,IY,FL)
         IXY=IY
         RNUM(J)=PAI2*FL+HKK(J)*I*BLOD*(COS(THETA(J)))
	 END DO

         END IF    ! WAVE.EQ.5

!  開始演算
!  注意：針對節點計算

        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,NS
        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 21
        READ(21) (F1(I),I=1,N1)    !OK
        READ(21) (XX(I),YY(I),ZZ(I),I=1,NT_N)
        READ(21) ((X(I,J),Y(I,J),Z(I,J),J=1,4),I=1,NT_E)
        CLOSE(21)
        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

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

!        DO I=1,N1
!        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)
!
      IF(WAVE.EQ.1) CALL REGULAR(IT,U,EE,DF,DT,HK,RTETA0,NCH,BLOD)
      IF(WAVE.EQ.2) CALL CROSS(IT, 
     &       U,E1,E2,DT,HK1,HK2,RTETA01,RTETA02,DF1,DF2,NCH,BLOD)
      IF(WAVE.EQ.3) CALL SOLITON(IT,U,XO,OMEGA,TC,RTETA0,DDT,NCH,BLOD) 
      IF(WAVE.EQ.4) CALL FREQ_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM,NCH)
      IF(WAVE.EQ.5) CALL MUTI_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM,NCH)

      DO I=1,N2
      F2B(I)=U(I)   !造波板條件
      END DO

!    2X2部分矩陣邊界積分方程式 !分成4小陣列存檔於[10]

      CALL HG_4(GB,G,NX,NT_N,NT_E,N1,NS,X,Y,Z,XX,YY,ZZ,JNB,NEI)

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

!         部分陣列分割如下
!                   11(N1xN1) | 12(N1xNS)
!                  ---------------------
!                   21(NSxN1) | 22(NSxNS)

!--------------------------連立方程式右邊已知項
        DO I=1,NX
        DO J=1,8
        CC(I,J)=0                    !歸零
        END DO
        END DO

        REWIND 10
        READ(10) ((GB(I,J),J=1,N1),I=1,N1)     !空讀 
        READ(10) ((GB(I,J),J=1,N1),I=1,NS)
        DO I=1,N1
        TEMP=0
        DO J=1,N2
        TEMP=TEMP+GB(I,J)*U(J)    !k12*U
        END DO
        CC(I,1)=TEMP+F1(I)    !k12*U+F1
        END DO
        READ(10) ((GB(I,J),J=1,NS),I=1,N1)     !空讀 
        READ(10) ((GB(I,J),J=1,NS),I=1,NS)   
        DO I=1,NS
        TEMP=0
        DO J=1,N2
        TEMP=TEMP+GB(I,J)*U(J)
        END DO
        CC(I,2)=TEMP
        END DO

        REWIND 31                      !  寫入連立方程式右邊已知項
        WRITE(31) (CC(I,1),I=1,N1)     !  實數部分
        WRITE(31) (CC(I,2),I=1,NS)     !  實數部分 
        WRITE(31) (CC(I,3),I=1,N1)     !  虛數部分 = 0
        WRITE(31) (CC(I,4),I=1,NS)     !  虛數部分 = 0

!        ENDFILE(10)

        REWIND 10  ! READ
        REWIND 1   ! WRITE
        REWIND 2   ! WRITE
!   建立11部分陣列係數  實數部分
        READ(10) ((GB(I,J),J=1,N1),I=1,N1)     !  讀取[11]
        WRITE(1) ((GB(I,J),J=1,N1),I=1,N1)     !寫入[11]實數

        READ(10) ((GB(I,J),J=1,NS),I=1,N1)     !  讀取[12]
!   建立[12]部分陣列係數  實數部分
        DO I=1,N1
        DO J=1,NS
        GB(I,J)=0
        END DO
        END DO
        WRITE(1) ((GB(I,J),J=1,NS),I=1,N1)     !寫入[12]實數


!   建立[21]部分陣列係數  實數部分
        READ(10) ((GB(I,J),J=1,N1),I=1,NS)      !讀取[21]實數
        WRITE(1) ((GB(I,J),J=1,N1),I=1,NS)      !寫入[21]實數

!   建立[22]部分陣列係數  實數部分
        READ(10) ((GB(I,J),J=1,NS),I=1,NS)      !讀取[22]實數
        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)       !寫入[22]實數

!   建立[11]部分陣列係數  虛數部分
        DO I=1,N1
        DO J=1,N1
        G(I,J)=0                                !  歸零
        END DO
        END DO
        WRITE(2)  ((G(I,J),J=1,N1),I=1,N1)      !寫入[11]虛數

!   建立[12]部分陣列係數  虛數部分
        BACKSPACE 10
        BACKSPACE 10
        READ(10) ((GB(I,J),J=1,NS),I=1,N1)      !讀取[12]部分陣列
        DO I=1,N1

        DO J=1,N2
        G(I,J)=0
        END DO
        DO J=1,N3+N4
        JJ=J+N2
        G(I,JJ)=GAMMA*GB(I,JJ)
        END DO      
        DO J=1,N5
        JJ=J+N2+N3+N4
        G(I,JJ)=KH*SQRT(1-KR(J)**2) *GB(I,JJ)
        END DO
        DO J=1,N6
        JJ=J+N2+N3+N4+N5
        G(I,JJ)=FR(J)*GB(I,JJ)
        END DO

        END DO   !  I

        WRITE(2) ((G(I,J),J=1,NS),I=1,N1)   !寫入[12]虛數

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

        WRITE(2) ((G(I,J),J=1,N1),I=1,NS)    !寫入[21]虛數

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

        READ(10) ((GB(I,J),J=1,N1),I=1,NS)    !空讀
        READ(10) ((GB(I,J),J=1,NS),I=1,NS)

        DO I=1,NS
        DO J=1,N2
        G(I,J)=0
        END DO
        DO J=1,N3+N4
        JJ=J+N2
        G(I,JJ)=GAMMA*GB(I,JJ)
        END DO
        DO J=1,N5
        JJ=J+N2+N3+N4
        G(I,JJ)=KH*SQRT(1.-KR(J)**2)*GB(I,JJ)
        END DO
        DO J=1,N6
        JJ=J+N2+N3+N4+N5
        G(I,JJ)=FR(J)*GB(I,JJ)
        END DO
        END DO
        WRITE(2) ((G(I,J),J=1,NS),I=1,NS)    !寫入[22]虛數

!
        CALL CQ2(GB,G,CC,NX,N1,NS)       ! 解連立方程式

        REWIND 32
        READ(32) (CC(I,5),I=1,N1)        !實數
        READ(32) (CC(I,6),I=1,NS)        !實數
        READ(32) (CC(I,7),I=1,N1)        !虛數
        READ(32) (CC(I,8),I=1,NS)        !虛數

      DO I=1,N1
      F1B(I)=CC(I,5)
      END DO

!  T+1時刻

      DO I=1,N1

      JK=JNB(I)
      DO JKS=1,JK
      K=NEI(I,JKS)

!     法線方向值         
      R1=X(K,3)-X(K,1)
      R2=Y(K,3)-Y(K,1)
      R3=Z(K,3)-Z(K,1)
      R4=X(K,4)-X(K,2)
      R5=Y(K,4)-Y(K,2)
      R6=Z(K,4)-Z(K,2)
      R=SQRT((R5*R3-R6*R2)**2+(R6*R1-R4*R3)**2+(R4*R2-R5*R1)**2)
      XN(JKS)=(R5*R3-R6*R2)/R
      YN(JKS)=(R6*R1-R4*R3)/R
      ZN(JKS)=(R4*R2-R5*R1)/R

!     元素面積
      T1=X(K,2)-X(K,1)
      T2=Y(K,2)-Y(K,1)
      T3=Z(K,2)-Z(K,1)
      T4=X(K,4)-X(K,1)
      T5=Y(K,4)-Y(K,1)
      T6=Z(K,4)-Z(K,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)
      A(JKS)=R+T
      END DO

      XT1=X(I,1);XT2=X(I,2);XT3=X(I,3);XT4=X(I,4)
      YT1=Y(I,1);YT2=Y(I,2);YT3=Y(I,3);YT4=Y(I,4)
      ZT1=Z(I,1);ZT2=Z(I,2);ZT3=Z(I,3);ZT4=Z(I,4)

      P1=-1;P2=-1

      XP(1,1)=.25*((-1+P2)*XT1+( 1-P2)*XT2+(1+P2)*XT3+(-1-P2)*XT4)
      XP(2,1)=.25*((-1+P1)*XT1+(-1-P1)*XT2+(1+P1)*XT3+( 1-P1)*XT4)
      XP(3,1)=XNN

      XP(1,2)=.25*((-1+P2)*YT1+( 1-P2)*YT2+(1+P2)*YT3+(-1-P2)*YT4)
      XP(2,2)=.25*((-1+P1)*YT1+(-1-P1)*YT2+(1+P1)*YT3+( 1-P1)*YT4)
      XP(3,2)=YNN

      XP(1,3)=.25*((-1+P2)*ZT1+( 1-P2)*ZT2+(1+P2)*ZT3+(-1-P2)*ZT4)
      XP(2,3)=.25*((-1+P1)*ZT1+(-1-P1)*ZT2+(1+P1)*ZT3+( 1-P1)*ZT4)
      XP(3,3)=ZNN

      CALL MINVS(XP,3,3,1.E-5,ILL)  ! Jacobian matrix inverse
      IF(ILL.NE.0) WRITE(*,11) ILL
11    FORMAT(1X,'ILL XP=',I7)

!   全體及無因次座標轉換
      IE1=ELEMENT(I,1)
      IE2=ELEMENT(I,2)
      IE3=ELEMENT(I,3)
      IE4=ELEMENT(I,4)
      FP1=(F1(IE1)-F1(IE4))/2.
      FP2=(F1(IE1)-F1(IE2))/2.
      FP3=F1B(I)   !近似

      FX=XP(1,1)*FP1+XP(1,2)*FP2+XP(1,3)*FP3
      FY=XP(2,1)*FP1+XP(2,2)*FP2+XP(2,3)*FP3
      FZ=XP(3,1)*FP1+XP(3,2)*FP2+XP(3,3)*FP3

      XX(I)=XX(I)+FX*DT
      YY(I)=YY(I)+FY*DT
      ZZ(I)=ZZ(I)+FZ*DT

      F1(I)=F1(I)+DT*(.5*(FX**2+FY**2+FZ**2)-GR*ZZ(I))  !T+1時刻F1

      END DO    !N1

! T+1時刻XX,YY,ZZ代入X,Y,Z
      DO I=1,N1
      DO K=1,4
      IE=ELEMENT(I,K)
      X(I,K)=XX(IE)
      Y(I,K)=YY(IE)
      Z(I,K)=ZZ(IE)
      END DO          
      END DO

      END DO  !IT
      END DO  !ITT

      REWIND 21
      WRITE(21) (F1(I),I=1,N1)
      WRITE(21) (XX(I),YY(I),ZZ(I),I=1,NT_N)
      WRITE(21) ((X(I,J),Y(I,J),Z(I,J),J=1,4),I=1,NT_E)
      CLOSE(21)

      STOP
      END


!****                   副程式                ****

      SUBROUTINE HG(GB,G,NX,X,Y,Z,NE,XX,YY,ZZ,NN,JNB,NEI,I,J,II,JJ)

      REAL::GB(NX,NX),G(NX,NX),X(NE,4),Y(NE,4),Z(NE,4)
      REAL::XX(NN),YY(NN),ZZ(NN)  
      REAL::XNN(10),YNN(10),ZNN(10),AA(10)
      REAL::GS(10),XG(10,4),YG(10,4),ZG(10,4)
      REAL::P(2,2),E(2,2)
      INTEGER::JNB(NN),NEI(NN,10)
  
!     確定節點各相鄰元素的編號，覓出該節點的K值，供HG等使用
!     I：源點  J：被積分點  JK：相鄰元素數  JKS：GUASS積分點
!     計算各元素基本解G及其法線方向導函數係數GB
!     NN=NT_N，NE=NT_E

      JK=JNB(J)

      DO JKS=1,JK
      K=NEI(J,JKS)

!     法線方向值         

      R1=X(K,3)-X(K,1)
      R2=Y(K,3)-Y(K,1)
      R3=Z(K,3)-Z(K,1)
      R4=X(K,4)-X(K,2)
      R5=Y(K,4)-Y(K,2)
      R6=Z(K,4)-Z(K,2)
      R=SQRT((R5*R3-R6*R2)**2+(R6*R1-R4*R3)**2+(R4*R2-R5*R1)**2)
      XNN(JKS)=(R5*R3-R6*R2)/R
      YNN(JKS)=(R6*R1-R4*R3)/R
      ZNN(JKS)=(R4*R2-R5*R1)/R

!     元素面積
      T1=X(K,2)-X(K,1)
      T2=Y(K,2)-Y(K,1)
      T3=Z(K,2)-Z(K,1)
      T4=X(K,4)-X(K,1)
      T5=Y(K,4)-Y(K,1)
      T6=Z(K,4)-Z(K,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)
      AA(JKS)=R+T

! 全體座標系與局部無因次座標系轉換|G|

      DO KP=1,4

      IF(KP.EQ.1) PC=-1
      IF(KP.EQ.2) PC=1
      IF(KP.EQ.3) PC=1
      IF(KP.EQ.4) PC=-1

      RX=.25*((-1+PC)*X(K,1)-(1+PC)*X(K,2)+(1+PC)*X(K,3)+(1-PC)*X(K,4))
      RY=.25*((-1+PC)*Y(K,1)-(1+PC)*Y(K,2)+(1+PC)*Y(K,3)+(1-PC)*Y(K,4))
      RZ=.25*((-1+PC)*Z(K,1)-(1+PC)*Z(K,2)+(1+PC)*Z(K,3)+(1-PC)*Z(K,4))

      IF(KP.EQ.1.OR.KP.EQ.2) THEN
      RSX=.5*(X(K,3)-X(K,4))
      RSY=.5*(Y(K,3)-Y(K,4))
      RSZ=.5*(Z(K,3)-Z(K,4))
      END IF

      IF(KP.EQ.3.OR.KP.EQ.4) THEN
      RSX=.5*(X(K,2)-X(K,1))
      RSY=.5*(Y(K,2)-Y(K,1))
      RSZ=.5*(Z(K,2)-Z(K,1))
      END IF

      GS1=RSY*RZ-RSZ*RY
      GS2=RSZ*RX-RSX*RZ
      GS3=RSX*RY-RSY*RX

      END DO

      GS(JKS)=SQRT(GS1**2+GS2**2+GS3**2)

! GUASS 數值積分計算點座標XG,YG,ZG

      N=2
      WL=1;WM=1
      QR=1./SQRT(3.)

      P(1,1)=-QR ; E(1,1)=-QR
      P(1,2)= QR ; E(1,2)=-QR
      P(2,1)= QR ; E(2,1)= QR
      P(2,2)=-QR ; E(2,2)= QR

      DO L=1,N
      DO M=1,N
      F1=(1-P(L,M)-E(L,M)+P(L,M)*E(L,M))/4.
      F2=(1+P(L,M)-E(L,M)-P(L,M)*E(L,M))/4
      F3=(1+P(L,M)+E(L,M)+P(L,M)*E(L,M))/4.
      F4=(1-P(L,M)+E(L,M)-P(L,M)*E(L,M))/4.
      LM=M+N*(L-1)
      XG(JKS,LM)=F1*X(K,1)+F2*X(K,2)+F3*X(K,3)+F4*X(K,4)
      YG(JKS,LM)=F1*Y(K,1)+F2*Y(K,2)+F3*Y(K,3)+F4*Y(K,4)
      ZG(JKS,LM)=F1*Z(K,1)+F2*Z(K,2)+F3*Z(K,3)+F4*Z(K,4)
      END DO
      END DO

      END DO    !JKS

      IF(I.EQ.J)  THEN

      GB(II,JJ)=0

!  I=J 時，應用GUASS積分對各相鄰元素的g作數值積分  (3.8)式

      N=2
      GX=0
      DO KK=1,JNB(J) 
      DO L=1,N
      DO M=1,N
      LM=M+N*(L-1)

      RILM=(XX(I)-XG(KK,LM))**2+(YY(I)-YG(KK,LM))**2+
     &(ZZ(I)-ZG(KK,LM))**2
      RILM2=SQRT(RILM)

      IF(LM.EQ.1) THEN
      GX=GX+WL*WM*(1-P(L,M)-E(L,M)+P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF
      IF(LM.EQ.2) THEN
      GX=GX+WL*WM*(1+P(L,M)-E(L,M)-P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF
      IF(LM.EQ.3) THEN
      GX=GX+WL*WM*(1+P(L,M)+E(L,M)+P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF
      IF(LM.EQ.4) THEN
      GX=GX+WL*WM*(1-P(L,M)+E(L,M)-P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF
      
      END DO
      END DO

      END DO  !KK

      G(II,JJ)=1./8/PAI*GX

      END IF

      IF(I.NE.J)  THEN

!  i<>j時，應用GUASS積分對各相鄰元素的h及g作數值積分 (3.4)式
!  X,Y,Z:節點座標
!  Xg,Yg,Zg:被積分相鄰元素GUASS積分點座標
!  XNN,YNN,ZNN被積分相鄰元素法線向量
!  AA：無因次座標轉換全體座標
      N=2
      
      HX=0
      GX=0
      DO KK=1,JNB(J)

      DO L=1,N
      DO M=1,N
      LM=M+N*(L-1)

      RILM=(XX(I)-XG(KK,LM))**2+(YY(I)-YG(KK,LM))**2
     &+(ZZ(I)-ZG(KK,LM))**2
      RILM2=SQRT(RILM)
      RN=(XG(KK,LM)-XX(I))*XNN(KK)+(YG(KK,LM)-YY(I))*YNN(KK)+
     &  (ZG(KK,LM)-ZZ(I))*ZNN(KK)
      RN=RN/RILM2

      IF(LM.EQ.1) THEN
      HX=HX+WL*WM*(1-P(L,M)-E(L,M)
     &+P(L,M)*E(L,M))/RILM*RN*GS(KK)
      GX=GX+WL*WM*(1-P(L,M)-E(L,M)+P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF

      IF(LM.EQ.2) THEN
      HX=HX+WL*WM*(1+P(L,M)-E(L,M)
     &-P(L,M)*E(L,M))/RILM*RN*GS(KK)
      GX=GX+WL*WM*(1+P(L,M)-E(L,M)-P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF

      IF(LM.EQ.3) THEN
      HX=HX+WL*WM*(1+P(L,M)+E(L,M)
     &+P(L,M)*E(L,M))/RILM*RN*GS(KK)
      GX=GX+WL*WM*(1+P(L,M)+E(L,M)+P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF

      IF(LM.EQ.4) THEN
      HX=HX+WL*WM*(1-P(L,M)+E(L,M)
     &-P(L,M)*E(L,M))/RILM*RN*GS(KK)
      GX=GX+WL*WM*(1-P(L,M)+E(L,M)-P(L,M)*E(L,M))/RILM2*GS(KK)
      END IF

      END DO   !M
      END DO   !L

      END DO   !KK

      GB(II,JJ)=-1./8/PAI*HX
      G(II,JJ) = 1./8/PAI*GX

      END IF

      RETURN
      END

      SUBROUTINE HG_4(GB,G,NX,NN,NE,N1,NS,X,Y,Z,XX,YY,ZZ,JNB,NEI)

      REAL::GB(NX,NX),G(NX,NX),X(NE,4),Y(NE,4),Z(NE,4)
      REAL::XX(NN),YY(NN),ZZ(NN)
      INTEGER::JNB(NN),NEI(NN,10)
  
!  NN=NT_N，NE=NT_E
!  2x2部份矩陣邊界積分方程式(1次元素)
!  將邊界積分方程式分割成2x2部份矩陣，減少矩陣容量並提高計算精度
!  G ：基本解作業矩陣存於L1
!  GB：基本解對法線方向導函作業矩陣存於L
!  結果存於L=10

! 分割部分矩陣
!     N1xN1  |  N1xNS
!    -------------
!     NSxN1  |  NSxNS

      L=10
      L1=11

      REWIND L
      REWIND L1

      DO II=1,N1
      I=II
      DO JJ=1,N1
      J=JJ
      CALL HG(GB,G,NX,X,Y,Z,NE,XX,YY,ZZ,NN,JNB,NEI,I,J,II,JJ)
      END DO
      END DO

      WRITE(L) ((GB(II,JJ),JJ=1,N1),II=1,N1)
      WRITE(L1) ((G(II,JJ),JJ=1,N1),II=1,N1)

      DO II=1,N1
      I=II
      DO JJ=1,NS
      J=JJ+N1
      CALL HG(GB,G,NX,X,Y,Z,NE,XX,YY,ZZ,NN,JNB,NEI,I,J,II,JJ)
      END DO
      END DO

      WRITE(L) ((GB(II,JJ),JJ=1,NS),II=1,N1)
      WRITE(L1) ((G(II,JJ),JJ=1,NS),II=1,N1)

      DO II=1,NS
      I=II+N1
      DO JJ=1,N1
      J=JJ
      CALL HG(GB,G,NX,X,Y,Z,NE,XX,YY,ZZ,NN,JNB,NEI,I,J,II,JJ)
      END DO
      END DO

      WRITE(L) ((GB(II,JJ),JJ=1,N1),II=1,NS)
      WRITE(L1) ((G(II,JJ),JJ=1,N1),II=1,NS)

      DO II=1,NS
      I=II+N1
      DO JJ=1,NS
      J=JJ+N1
      CALL HG(GB,G,NX,X,Y,Z,NE,XX,YY,ZZ,NN,JNB,NEI,I,J,II,JJ)
      END DO
      END DO

      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,NX,N1,NS,L,31)       ! 2x2部份矩陣逆矩陣
      CALL MTR2(GB,G,NX,N1,NS,31,L1,L)    ! 2x2部份矩陣相乘，存於L=10

      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 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 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 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矩陣的複數連立方程式，結果存於32

        REAL::A(NX,NX),B(NX,NX),C(NX,8)

        REWIND 31    ! INPUT
        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)        ! 虛數部分

        CALL MIC2(A,B,NX,N,NS)   !2X2矩陣的逆矩陣，結果實數存於3，虛數存於4

        REWIND 3
        REWIND 4

        READ(3) ((A(I,J),J=1,N),I=1,N)  ![11]
        READ(4) ((B(I,J),J=1,N),I=1,N)

        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)  ![12]
        READ(4) ((B(I,J),J=1,NS),I=1,N)

        DO I=1,N
        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,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)  ![21]
        READ(4) ((B(I,J),J=1,N),I=1,NS)

        DO I=1,NS
        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,6)=R+S
        C(I,8)=T-U
        END DO

        READ(3) ((A(I,J),J=1,NS),I=1,NS)  ![22]
        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


        REWIND 32    ! OUTPUT

        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
	GOTO 6
 2	IF((SHG-1.0).GE.0) GOTO 4
	X=SQRT(SHG)
	GOTO 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 RANDU(IX,IY,FL)

      IY=MOD(12869*IX+6925,32768)
      FL=FLOAT(IY)/2.0**15

      RETURN
      END


      SUBROUTINE REGULAR(IT,U,E,DF,DT,HK,RTETA0,NCH,BLOD)

!       規則波造波

	REAL U(NCH)
        DO  I=1,NCH
        U(I)=E*SIN(DF*DT*(IT-1)-HK*(I-1)*BLOD*COS(RTETA0))
        END DO
        RETURN
        END

        SUBROUTINE CROSS
     &  (IT,U,E1,E2,DT,HK1,HK2,RTETA01,RTETA02,DF1,DF2,NCH,BLOD)

!        短峰波造波

	REAL U(NCH)
        DO I=1,NCH
        U(I)=E1*SIN(DF1*DT*(IT-1)-HK1*(I-1)*BLOD*COS(RTETA01))
     &      +E2*SIN(DF2*DT*(IT-1)+HK2*(I-1)*BLOD*COS(RTETA02))
        END DO
        RETURN
        END

      SUBROUTINE SOLITON(IT,U,XO,OMEGA,TC,RTETA0,DDT,NCH,BLOD)

!       孤立波造波

	REAL U(NCH)
        DO I=1,NCH
        U(I)=-XO*OMEGA/COSH(OMEGA*(DDT*IT-TC)-(I-1)*BLOD*COS(RTETA0))**2
        END DO
        RETURN
        END

      SUBROUTINE FREQ_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM,NCH)

!       單方向不規則波造波

      REAL RNUM(NSPEC),U(NCH),S0(NSPEC),RE(NSPEC)
      PAI2=6.283185308
      DO I=1,NCH
      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))!差異在RNUM
      END DO
      U(I)=TEMP
      END DO
      RETURN
      END


      SUBROUTINE MUTI_WAVE(IT,U,S0,DF,NSPEC,RE,DDT,RNUM,NCH)

!         多方向不規則波造波

      REAL RNUM(NSPEC),U(NCH),S0(NSPEC),RE(NSPEC)
      PAI2=6.283185308
      DO I=1,NCH
      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))!差異在RNUM
      END DO
      U(I)=TEMP
      END DO
      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