!程序簡介
!1968年1月青島驗潮站逐時水位觀測數據,進行潮汐調和分析和預報
!選入的分潮(18項，其中有兩項為隨從分潮)：  
! 1 MSf  
! 2 Q1   
! 3 O1   
!17 P1 (K1的隨從分潮)    由長期資料得振幅比 P1/K1=0.274
! 4 K1   
! 5 J1   
! 6 OO1  
! 7 N2   
! 8 M2   
! 9 S2   
!18 K2  (S2的隨從分潮)   由長期資料得振幅比 K2/S2=0.249
!10 MO3  
!11 MK3  
!12 SK3  
!13 MN4  
!14 M4      淺水分潮
!15 MS4     淺水分潮
!16 2MS6
 
 
!數據是從1968年1月1日0時開始的，即：
!          數據第一行：100  140  213  297  363  379  342  277  203  142 79  20 
!對應時間1968年1月1日： 0時 1時  2時  …………        …………   ………… 11時  依此類推
!共744個數據
!選區時十六日零時為中間時刻，需從原潮位資料中去掉后面的23個數據，剩下721個潮位數據
!回報青島1968年1月逐時潮位時用到的四大分潮：O1,K1,M2,S2，兩個淺水分潮：M4,MS4
!----------------------------------------------------------------------------------------
!總體構架
!!第一步計算分潮的角速率
!第二步計算天文要素tao，s,h,p,N,p_
!第三步計算分潮的天文初相角V0
!第四步計算基本分潮的f，u
!第五步計算天文初相角
!第六步數據讀入
!第七步計算方程
!第八步解X的賦值
!第九步計算調和常數
!第十步求744個逐時水位數據
!第十一步求每天高低潮潮高和潮時
!第十二步使用Matlab分析資料
!-------------------------------------------------------------------------------------------



!Gauss_Jordan法求解線形方程組模塊
module LinearAlgebra
  implicit none
contains
! Gauss_Jordan法
subroutine Gauss_Jordan(A,S,ANS)
  implicit none
  real    :: A(:,:)
  real    :: S(:)
  real    :: ANS(:)
  real, allocatable :: B(:,:)
  integer :: i, N
  N = size(A,1)  
  allocate(B(N,N))
  ! 保存原先的矩陣A,及數組S
  B=A 
  ANS=S
  ! 把B化成對角線矩陣(除了對角線外,都為0)
  call Upper(B,ANS,N) ! 先把B化成上三角矩陣
 
  call Lower(B,ANS,N) ! 再把B化成下三角矩陣
  ! 求解
  forall(i=1:N)
    ANS(i)=ANS(i)/B(i,i) 
  end forall
  return
end subroutine Gauss_Jordan
! 求上三角矩陣的子程序
subroutine Upper(M,S,N)
  implicit none
  integer :: N
  real    :: M(N,N)
  real    :: S(N)
  integer :: I,J
  real :: E
  do I=1,N-1
    do J=I+1,N              
      E=M(J,I)/M(I,I)
      M(J,I:N)=M(J,I:N)-M(I,I:N)*E
      S(J)=S(J)-S(I)*E
    end do
  end do
  return
end subroutine Upper
! 求下三角矩陣的子程序
subroutine Lower(M,S,N)
  implicit none
  integer :: N
  real    :: M(N,N)
  real    :: S(N)
  integer :: I,J
  real :: E
  do I=N,2,-1
    do J=I-1,1,-1           
      E=M(J,I)/M(I,I)
      M(J,1:N)=M(J,1:N)-M(I,1:N)*E
      S(J)=S(J)-S(I)*E
    end do
  end do
  return
end subroutine Lower
end module
!----------------------------------------


program tide
use LinearAlgebra
implicit none
!定義變量
integer,parameter::P=16,Q=2,M=721,year=1968,D=15	!共16個主要分潮，2個隨從分潮,觀測有721個數據 資料的年數1968
!D為從year年1月1日算起所經過的日數
real,parameter::pi=3.14159265/180	!角度轉化為弧度
integer::yy !從1900年至year年（year年除外）間的閏年數
integer::i,j,k,o,r,flagh,flagl	!循環變量
!-------------------------------
!real,dimension(1:744)::odat	!原始數據
real,dimension(-360:360)::dat	!觀測數據,取奇數個721
real::wt,ws,wh,wp,wN,wp_	!σtao、σs、σh、σp、σN、σp'
integer,dimension(6,P+Q) ::dn	!各個分潮杜得森數
real,dimension(P+Q)::w,f,u,V0,V	!各個分潮的角速率σ,交點因子,交點訂正角,理論初相角,V為天文初相角，V=V0+u
real::fM2,fO1,fK1,fJ1,fOO1,fk2,uM2,uO1,uK1,uJ1,uOO1,uk2	!幾個重要分潮的交點因子
real::t0,s0,h0,p0,N0,p_0,N !用來計算V
!-----------------------------------------
real,dimension(P+Q)::ks,yt	!keseξ,yitaη各個分潮的ξ和η
real,dimension(0:P,0:P+Q)::A
real,dimension(P,P+Q)::B	!方程組中用到的系數A,B
real,dimension(0:P)::FF
real,dimension(P)::GG	!方程組中用到的系數F,G
real,dimension(Q)::kk,fai,sc=(/4,9/) !差比關系中分潮的振幅比和遲角差,4,9分別代表隨從分潮P1,K2所對應主分潮所在位置K1-4,S2-9
real,dimension(P+Q)::a1,b1,H,g	!待求值，各個分潮調和常數
real::S		!平均水位
real,dimension((2*(P+Q)+1),(2*(P+Q)+1))::AA
real,dimension(2*(P+Q)+1)::X,BB !通過AA*X=BB解方程組
!-----------------------------------------------------------
!real,dimension(-360:383)::hh    !每小時的潮高
!real,dimension(-21600:22980)::tm,hm	!每分鐘的潮高
!real,dimension(1:100)::hmh,hml   !hmh定義為逐分高潮水位，hml定義為逐分低潮水位
!real,dimension(1:100)::th,tl   !th定義為高潮時(min),tl定義為低潮時(min)
!real,dimension(:),allocatable::hwh,hwl !未插值求高低潮結果
!real,dimension(:),allocatable::twh,twl 
!real,dimension(150,-2:2)::ycz  
!real,dimension(:),allocatable::cz,AAA,BBB,CCC,DDD,ttt,tr,hcz,tcz   !插值法求高低潮
!real,dimension(:),allocatable::hczh,hczl !插值法求高低潮結果
!real,dimension(:),allocatable::tczh,tczl 
!integer,dimension(6)::yb=(/3,4,8,9,14,15/)	!四大分潮和兩個淺水分潮回報潮位時用到的分潮，依次為O1,K1,M2,S2,M4,MS4
!------------------------------------------------------------------------------------------------------------------


!各個分潮杜得森數賦值
dn(1:6,1)=(/0,2,-2,0,0,0/)		!Msf
dn(1:6,2)=(/1,-2,0,1,0,0/)		!Q1             
dn(1:6,3)=(/1,-1,0,0,0,0/)		!O1
dn(1:6,4)=(/1,1,0,0,0,0/)			!K1
dn(1:6,5)=(/1,2,0,-1,0,0/)		!J1
dn(1:6,6)=(/1,3,0,0,0,0/)			!OO1
dn(1:6,7)=(/2,-1,0,1,0,0/)		!N2						
dn(1:6,8)=(/2,0,0,0,0,0/)			!M2
dn(1:6,9)=(/2,2,-2,0,0,0/)		!S2
dn(1:6,10)=(/3,-1,0,0,0,0/)		!MO3
dn(1:6,11)=(/3,1,0,0,0,0/)		!MK3
dn(1:6,12)=(/3,3,-2,0,0,0/)		!SK3
dn(1:6,13)=(/4,-1,0,1,0,0/)		!MN4
dn(1:6,14)=(/4,0,0,0,0,0/)		!M4 淺水分潮
dn(1:6,15)=(/4,2,-2,0,0,0/)		!MS4 淺水分潮
dn(1:6,16)=(/6,2,-2,0,0,0/)		!2MS6
dn(1:6,17)=(/1,1,-2,0,0,0/)		!P1  (K1的隨從分潮)
dn(1:6,18)=(/2,2,0,0,0,0/)		!K2   (S2的隨從分潮)
!-----------------------------------------------------------------------------------------------------

!print*,"如果原數據文件和tideanalysis.exe在同一個文件夾中，請輸入ok"
!read(*,*)

!第一步
!計算分潮的角速率

!給σt、σs、σh、σp、σN、σp'賦值,角度轉化為弧度
!單位:弧度/平太陽時
wt=14.49205212*pi	
ws=0.54901653*pi	
wh=0.04106864*pi
wp=0.00464183*pi
wN=0.00220641*pi
wp_=0.00000196*pi

!各個分潮的角速率σ
do i=1,P+Q	
w(i)=dn(1,i)*wT+dn(2,i)*ws+dn(3,i)*wh+dn(4,i)*wp+dn(5,i)*wN+dn(6,i)*wp_
end do
!----------------------------------------------------------------------------


!第二步
!計算天文要素tao，s,h,p,N,p_
!yy為1900年至year年（year年除外）間的閏年數，若year年為閏年，則該年也算在內。、
yy=int((year-1901)/4.)

!選16日零時為中間時刻，t=0
!由于程序計算三角函數化弧度*pi
s0=(277.025+129.38481*(year-1900)+13.17640*(D+yy))*pi	!書上y＝1968，D＝15
h0=(280.190-0.23872*(year-1900)+0.98565*(D+yy))*pi	!課本公式給出的這幾個參量時間取的格林威治時
p0=(334.385+40.66249*(year-1900)+0.11140*(D+yy))*pi	
N0=(259.157-19.32818*(year-1900)-0.05295*(D+yy))*pi
p_0=(281.221+0.01718*(year-1900)+0.000047*(D+yy))*pi
N=N0
t0=h0-s0
!------------------------------------------------------------------------------------------------



!第三步
!計算分潮的天文初相角V0
do i=1,P+Q	
V0(i)=dn(1,i)*t0+dn(2,i)*s0+dn(3,i)*h0+dn(4,i)*p0+dn(5,i)*N0+dn(6,i)*p_0
end do
!------------------------------------------------------------------------------------------------



!第四步計算基本分潮的f，u 選擇區時時間系統

!遲角差為區時專用遲角

!計算用到的幾個重要分潮交點因子
fM2=1.0004-0.0373*cos(N)+0.0003*cos(2*N)	
fO1=1.0089+0.1871*cos(N)-0.0147*cos(2*N)+0.0014*cos(3*N)
fK1=1.0060+0.1150*cos(N)-0.0088*cos(2*N)+0.0006*cos(3*N)
fJ1=1.0129+0.1676*cos(N)-0.0170*cos(2*N)+0.0016*cos(3*N)
fOO1=1.1027+0.6504*cos(N)+0.0317*cos(2*N)-0.0014*cos(3*N)
fk2=1.0241-0.2863*cos(N)+0.0083*cos(2*N)-0.0015*cos(3*N)

!計算用到的幾個重要分潮交點訂正角
uM2=(-2.14*sin(N))*pi				
uO1=(10.80*sin(N)-1.34*sin(2*N)+0.19*sin(3*N))*pi
uK1=(-8.86*sin(N)+0.68*sin(2*N)-0.07*sin(3*N))*pi
uJ1=(-12.94*sin(N)+1.34*sin(2*N)-0.19*sin(3*N))*pi
uOO1=(-36.68*sin(N)+4.02*sin(2*N)-0.57*sin(3*N))*pi
uk2=(-17.74*sin(N)+0.68*sin(2*N)-0.04*sin(3*N))*pi


!查表,各個分潮交點因子
f(1)=fM2	
f(2)=fO1
f(3)=fO1
f(4)=fK1
f(5)=fJ1
f(6)=fOO1
f(7)=fM2
f(8)=fM2
f(9)=1
f(10)=fM2*fO1
f(11)=fM2*fK1
f(12)=fK1
f(13)=fM2**2
f(14)=fM2**2
f(15)=fM2
f(16)=fM2**2

!查表，各個分潮交點訂正角
u(1)=-uM2	
u(2)=uO1
u(3)=uO1
u(4)=uK1
u(5)=uJ1
u(6)=uOO1
u(7)=uM2
u(8)=uM2
u(9)=0
u(10)=uM2+uO1
u(11)=uM2+uK1
u(12)=uK1
u(13)=uM2+uM2
u(14)=uM2+uM2
u(15)=uM2
u(16)=2*uM2
f(18)=fk2                                     
u(18)=uk2

!p1,k2分潮，用升交點黃經表示f，u
f(17)=sqrt(((0.0112296*sin(N)-0.00798*sin(2*N))**2)+((1-0.0112296*cos(N)+0.000798*cos(2*N))**2))                                                 !不對，需要自己計算   !P1   P1                                                                      
u(17)=atan((0.0112296*sin(N)-0.00798*sin(2*N))/(1-0.0112296*cos(N)+0.000798*cos(2*N)))
!-----------------------------------------------------------------------------------------------





!第五步
!計算天文初相角	 
do j=1,18
	 V(j)=V0(j)+u(j)
end do
!--------------------------------------



!第六步
!讀入原始數據
!open(1,file='odat.txt')
!read(1,fmt="(F3.0,11F4.0)"),((odat(i+j),j=0,11),i=1,744,12)
!close(1)

!open(1,file='分析用原始數據.txt')
!write(1,fmt="(F4.0)"),((odat(i)),i=1,744)
!close(1)


! data.txt是原始數據去掉23個后的
!輸出原始數據去掉23個后的data.txt
!open(1,file='data.txt')
!do i=1,M
!write(1,fmt="(F5.1)")odat(i)
!end do
!close(1)

!潮位資料數據讀入數組dat.txt中用來計算方程中的系數

open(unit=2,file='data.txt',status='old')	

read(2,*) dat

close(2)

!第七步
!計算方程

!各個分潮的ξ和η賦值
do i=1,P+Q	
ks(i)=f(i)*cos(V(i))
yt(i)=f(i)*sin(V(i))
end do



!A,B賦值
do i=0,P	
   do j=0,P+Q
   !A00
      if(i==0.and.j==0) then
	  A(i,j)=M
	  !A0j	   
      else if(i==0.and.j/=0)	then
	  A(i,j)=sin(0.5*M*w(j))/sin(0.5*w(j))
	  !Aj0		
      else if(i/=0.and.j==0) then
	  A(i,j)=sin(0.5*M*w(i))/sin(0.5*w(i))
	  !Ajj,Bjj		
      else if(i==j) then
         A(i,j)=0.5*(M+sin(M*w(j))/sin(w(j))) 
         B(i,j)=0.5*(M-sin(M*w(j))/sin(w(j)))					
      else
         A(i,j)=0.5*sin(0.5*M*(w(i)-w(j)))/sin(0.5*(w(i)-w(j)))+0.5*sin(0.5*M*(w(i)+w(j)))/sin(0.5*(w(i)+w(j))) 									!Aij,Bij
         B(i,j)=0.5*sin(0.5*M*(w(i)-w(j)))/sin(0.5*(w(i)-w(j)))-0.5*sin(0.5*M*(w(i)+w(j)))/sin(0.5*(w(i)+w(j))) 
      endif
   end do
end do


!F0,Fi,Gi賦值
!F0
FF(0)=0
do i=-360,360	
FF(0)=FF(0)+dat(i)
end do
!Fi,Gi
do i=1,P	
FF(i)=0
GG(i)=0
  do j=-360,360
   FF(i)=FF(i)+dat(j)*cos(j*w(i))
   GG(i)=GG(i)+dat(j)*sin(j*w(i))