tide.f90

1968年1月青岛验潮站逐时水位观测数据,进行潮汐调和分析和预报 !选入的分潮(18项) 数据是从1968年1月1日0时开始的

!程序简介
!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))