中期水位资料对潮汐进行调和分析

!利用1996年7月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报

program work

implicit none

character*80::a1

character(len=5),dimension(62,16)::aa

integer::bb(62,12),c(62,2),caita(-371:371),i,i1,i2,j,t1

real::N0,n(13,6),a(0:13,0:13),b(1:13,1:13),s,s0,s1,s2,s3,sa,hh !n代表Doodson代码；a,b为系数矩阵

real::xiaoa(0:13),xiaob(13),gg1,gg2,pjchaocha,t,ma,mi

!计算法方程所需的参数

real,dimension(1:13)::w,u,f,V0,f1(0:13),f2

!f1和f2为法方程右边系数

real,dimension(13)::sita,h,g,r,h0(13),g0(13),h1(13),g1(13)

!调和常数参数

real,dimension(-371:371)::caita1,caita3,caita4,caita8,caita9,caita5,caita11 !主分潮、浅水分潮的潮高数值

real,dimension(:),allocatable::hightide,lowtide,chaocha

!高低潮数值

integer,dimension(:),allocatable::hightrq,lowtrq,hight,hightt,lowt,lowtt !读取数据，把潮位数据赋值给bb，把年月份数据赋值给c

open(unit=2,file='XM_July1996.dat')

read(2,'(a)')a1

print*,'数据文件的第一行信息:',a1

do i=1,62

read(2,'(16a5)')aa(i,:)

end do

do i=1,62

read(aa(i,5:16),*)bb(i,:)

read(aa(i,3:4),*)c(i,:)

end do

do i=1,62

c(i,2)=int(real(c(i,2))/10.0)

end do

close(2)

!计算分潮角速率w

w=(/0.002822,0.037219,0.038731,0.041781,0.163845,0.241534,0.078999 ,&

&0.080511,0.083333,0.122292,0.161023,0.041553,0.083561/)

w=360*w

print*

print*,'角速率w:',w

!计算N0 (middle time:1996-7-16 ; data sum:744, middle number:372 )

N0=259.157-19.32818*(1996-1900)-0.05295*(31*3+30*2+29+15+int((9 5.0)/4.0)) !初始升交点平均黄经

N0=-(0.00220641*3+N0)

print*

!转换成格林威治时间

print*,'N0:',N0

!数字序号对应选取的分潮，但将5、6(P1、K2)分别与12、13(MS4、M6)对调，其中P1、K2为随从分潮

!计算交点订正角u

u(3)=10.8*sind(N0)-1.34*sind(2*N0)+0.19*sind(3*N0)

u(4)=-8.86*sind(N0)+0.68*sind(2*N0)-0.07*sind(3*N0)

u(8)=-2.14*sind(N0)

u(13)=-17.74*sind(N0)+0.68*sind(2*N0)-0.04*sind(3*N0)

u(1)=-u(8)

u(2)=u(3)

u(7)=u(8)

u(9)=0

u(10)=u(8)+u(4)

u(11)=2*u(8)

u(5)=u(8)

u(6)=3*u(8)

u(12)=0 !

print*

print*,'交点订正角u:',u

!计算交点因子f

f(3)=1.0089+0.1871*cosd(N0)-0.147*cosd(2*N0)+0.0014*cosd(3*N0) f(4)=1.006+0.115*cosd(N0)-0.0088*cosd(2*N0)+0.0006*cosd(3*N0) f(8)=1.0004-0.0373*cosd(N0)+0.0003*cosd(2*N0)

f(13)=1.0241+0.2863*cosd(N0)+0.0083*cosd(2*N0)-0.0015*cosd(3*N0 )

f(1)=f(8)

f(2)=f(3)

f(7)=f(8)

f(9)=1

f(10)=f(8)*f(4)

f(11)=f(8)**2

f(5)=f(8)**2

f(6)=f(8)**3

f(12)=1 !

print*

print*,'交点因子f:',f

!查表得到的Doodson代码n(1,:)=(/0,2,-2,0,0,0/)

n(2,:)=(/1,-2,0,1,0,0/)

n(3,:)=(/1,-1,0,0,0,0/)

n(4,:)=(/1,1,0,0,0,0/)

n(5,:)=(/4,2,-2,0,0,0/)

n(6,:)=(/6,0,0,0,0,0/)

n(7,:)=(/2,-1,0,1,0,0/)

n(8,:)=(/2,0,0,0,0,0/)

n(9,:)=(/2,2,-2,0,0,0/)

n(10,:)=(/3,1,0,0,0,0/)

n(11,:)=(/4,0,0,0,0,0/)

n(12,:)=(/1,1,-2,0,0,0/)

n(13,:)=(/2,2,0,0,0,0/)

!计算V0

do i=1,13