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

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

!利用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),g 1(13) !调和常数参数

real,dimension(-

371:371)::caita1,caita3,caita4,caita8,caita9,caita5,caita 11 !主分潮、浅水分潮的潮高数值

real,dimension(:),allocatable::hightide,lowtide,chaoc ha !高低潮数值

integer,dimension(:),allocatable::hightrq,lowtrq,high t,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.24 1534,0.078999,&

&0.080511,0.083333,0.122292,0.161023,0.041553,0.08356 1/) 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((95.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