1、!利用 1996 年 7 月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报program workimplicit nonecharacter*80:a1character(len=5),dimension(62,16):aainteger:bb(62,12),c(62,2),caita(-371:371),i,i1,i2,j,t1real: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,g
2、g2,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
3、!高低潮数值integer,dimension(:),allocatable:hightrq,lowtrq,hight,hightt,lowt,lowtt!读取数据，把潮位数据赋值给 bb，把年月份数据赋值给 copen(unit=2,file=XM_July1996.dat)read(2,(a)a1print*,数据文件的第一行信息:,a1do i=1,62read(2,(16a5)aa(i,:)end dodo i=1,62read(aa(i,5:16),*)bb(i,:)read(aa(i,3:4),*)c(i,:)end dodo i=1,62c(i,2)=int(real(c(i,2
4、)/10.0)end doclose(2)!计算分潮角速率 ww=(/0.002822,0.037219,0.038731,0.041781,0.163845,0.241534,0.078999, 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(
5、MS4 、M6) 对调，其中 P1、K2 为随从分潮!计算交点订正角 uu(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)=0u(10)=u(8)+u(4)u(11)=2*u(8)u(5)=u(8)u(6)=3*u(8)u(12)=0 !p
6、rint*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)=1f(1
7、0)=f(8)*f(4)f(11)=f(8)*2f(5)=f(8)*2f(6)=f(8)*3f(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(
8、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/)!计算 V0do i=1,13V0(i)=(14.49205212*3+180)*n(i,1)+(0.54901653*3+277.025+129.3848*96+13.1764*(220)*n(i,2)+g=0;i1=0doh0=hg0=g!A 阵do i=0,11s1=0do j=0,11s1=s1+xiaoa(j)*a(i,j)end doxiaoa(i)=-s1/a(i,i)+f1(i)/a(i,i)+xi
9、aoa(i)-xiaoa(12)*a(i,12)/a(i,i)-xiaoa(13)*a(i,13)/a(i,i)end do!B 阵do i=1,11s1=0do j=1,11s1=s1+xiaob(j)*b(i,j)end doxiaob(i)=-s1/b(i,i)+f2(i)/b(i,i)+xiaob(i)-xiaob(12)*b(i,12)/b(i,i)-xiaob(13)*b(i,13)/b(i,i)end do!计算调和常数 h,gdo j=1,11sita(j)=atand(xiaob(j)/xiaoa(j)+180r(j)=sqrt(xiaoa(j)*2+xiaob(j)*2)g
10、(j)=V0(j)+u(j)+sita(j)h(j)=r(j)/f(j)end dodo i=1,11do while(g(i)360.or.g(i)360)thendo g(i)=g(i)-360if(g(i)0.and.g(i)0.and.g(i)230.or.gg1230)thendogg1=gg1-360if(gg1-130)exitend doelseend ifif(gg1-130.and.gg1230.or.gg2230)thendogg2=gg2-360if(gg2-130)exitend doelseend ifif(gg2-130.and.gg2caita(j-1).and
11、.caita(j)caita(j+1)theni1=i1+1 !高潮个数end ifif(caita(j)caita1(j-1).and.caita1(j)caita1(j+1)theni1=i1+1do t=j-1,j+1,0.01s=s0+f(3)*h(3)*cosd(w(3)*(t)+V0(3)+u(3)-g(3)+&f(4)*h(4)*cosd(w(4)*(t)+V0(4)+u(4)-g(4)+&f(5)*h(5)*cosd(w(5)*(t)+V0(5)+u(5)-g(5)+&f(8)*h(8)*cosd(w(8)*(t)+V0(8)+u(8)-g(8)+&f(9)*h(9)*cosd
12、(w(9)*(t)+V0(9)+u(9)-g(9)+&f(11)*h(11)*cosd(w(11)*(t)+V0(11)+u(11)-g(11)sa=s0+f(3)*h(3)*cosd(w(3)*(t+1/60.0)+V0(3)+u(3)-g(3)+&f(4)*h(4)*cosd(w(4)*(t+1/60.0)+V0(4)+u(4)-g(4)+&f(5)*h(5)*cosd(w(5)*(t+1/60.0)+V0(5)+u(5)-g(5)+&f(8)*h(8)*cosd(w(8)*(t+1/60.0)+V0(8)+u(8)-g(8)+&f(9)*h(9)*cosd(w(9)*(t+1/60.0)
13、+V0(9)+u(9)-g(9)+&f(11)*h(11)*cosd(w(11)*(t+1/60.0)+V0(11)+u(11)-g(11)if(ssa)thenmi=salowtt(i2)=int(t-floor(t)*60)lowt(i2)=floor(t)+371end ifend dolowtide(i2)=mielseend if end do!将高潮位写入hightide.txt，第 1 列为潮位，第 2、3 列为时刻open(unit=2,file=hightide.txt) do i=1,i1write(2,*)hightide(i),hight(i),hightt(i)end
14、 doclose(2)do i=1,i1j=1if(hight(i)23)hight(i)=hight(i)-24j=j+1hightrq(i)=jend doend doopen(unit=2,file=hightidexiu.xls)do i=1,i1hightide(i)=int(hightide(i)+0.5)*1.0write(2,*)hightide(i),hightrq(i),hight(i),hightt(i) !将高潮位及时刻写入hightidexiu.xls，第 1 列为潮位，第 2、3、4 列为天数、小时、分钟end doclose(2)!将低潮位写入lowtide.tx
15、t，第 1 列为潮位，第 2、3 列为时刻open(unit=2,file=lowtide.txt)do i=1,i2write(2,*)lowtide(i),lowt(i),lowtt(i)end doclose(2)do i=1,i2j=1if(lowt(i)23)lowt(i)=lowt(i)-24j=j+1lowtrq(i)=jend doend doopen(unit=2,file=lowtidexiu.xls)do i=1,i2write(2,*)lowtide(i),lowtrq(i),lowt(i),lowtt(i) !将低潮位及时刻写入lowtidexiu.xls，第 1 列
16、为潮位，第 2、3、4 列为天数、小时、分钟end doclose(2)!将四大主分潮潮位数据写入zhuyaotides.txt中，第 1 列为 O1，第 2列为 K1，第 3 列为 M2，第 4 列为 S2open(unit=2,file=zhuyaotides.txt)do i=-371,371write(2,(4f7.1)caita3(i),caita4(i),caita8(i),caita9(i)end doclose(2)!将浅水分潮潮位数据写入qianshuitides.txt中 ,第 1 列为 MS4，第 2列为 M4open(unit=2,file=qianshuitides.
17、txt)do i=-371,371write(2,(2f7.1)caita5(i),caita11(i)end doclose(2)!计算平均潮差s=0do i=1,size(chaocha)s=s+(hightide(i)-lowtide(i)end dopjchaocha=s/size(chaocha)print*print*,平均潮差 :,pjchaochaprint*deallocate(hightide,lowtide,hight,lowt,chaocha,hightrq,lowtrq) !释放掉潮位数据!判断 XM 站潮汐类型hh=(h(3)+h(4)/h(8)if(hh=0.5.and.hh=2.0.and.hh=4.0)then print*,潮型数： ,hh,所以 XM 站潮汐类型为全日潮elseprint*,errorend ifend program work
