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

!利用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,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，把年月份数据赋值给copen(unit=2,file=XM_July1996.dat)read(2,(a)a1print*,数据文件的第一行信息:,a1do i=1,62 read(2,(16a5)aa(i,:)end dodo i=1,62 read(aa(i,5:16),*)bb(i,:) read(aa(i,3:4),*)c(i,:)end dodo i=1,62 c(i,2)=int(real(c(i,2)/10.0)end doclose(2)!计算分潮角速率ww=(/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*wprint*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为随从分潮!计算交点订正角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 !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)=1f(10)=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(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,13 V0(i)=(14.49205212*3+180)*n(i,1)+(0.54901653*3+277.025+129.3848*96+13.1764*(220)*n(i,2)+& &(0.04106864*3+280.190-0.23872*96+0.98565*(220)*n(i,3)+(0.00464183*3+334.385+40.66249*96+& &0.11140*(220)*n(i,4)-(0.00220641*3+259.157-19.32818*96-0.05295*(220)*n(i,5)+& &(0.00000196*3+281.221+0.01718*96+0.000047*(220)*n(i,6)end doprint*print*,初始幅角V0:,V0!设caita为潮高数据do i=1,61 j=-371+(i-1)*12 caita(j:j+11)=bb(i,:)end docaita(361:371)=bb(62,1:11)!计算法方程等式右边的数据，相邻数据时间间隔为1小时f1(0)=sum(caita)!f1为A阵中除第一行外的等式右边一维数据f1(1:13)=0do i=1,13 do j=-371,371 f1(i)=f1(i)+caita(j)*cosd(j*w(i) end doend do!f2为B阵中等式右边的一维数据f2=0do i=1,13 do j=-371,371 f2(i)=f2(i)+caita(j)*sind(j*w(i) end doend do !计算A阵中的系数矩阵Aa(0,0)=743do j=1,13 a(0,j)=sind(743.0/2*w(j)/sind(0.5*w(j) a(j,0)=a(0,j)end dodo j=1,13 a(j,j)=0.5*(743+sind(743.0*w(j)/sind(w(j)end dodo i=1,13 do j=i+1,13 a(i,j)=0.5*(sind(743.0/2*(w(i)-w(j)/sind(0.5*(w(i)-w(j)+& & sind(743.0/2*(w(i)+w(j)/sind(0.5*(w(i)+w(j) a(j,i)=a(i,j) end doend doprint*print*,系数矩阵A:,a!计算B阵中的系数矩阵Bdo j=1,13b(j,j)=0.5*(743-sind(743.0*w(j)/sind(w(j)end dodo i=1,13 do j=i+1,13 b(i,j)=0.5*(sind(743.0/2*(w(i)-w(j)/sind(0.5*(w(i)-w(j)-& & sind(743.0/2*(w(i)+w(j)/sind(0.5*(w(i)+w(j) b(j,i)=b(i,j) end doend doprint*print*,系数矩阵B:,b!Guass-Seidel迭代法求解方程组h=0;g=0;i1=0doh0=hg0=g!A阵do i=0,11 s1=0 do j=0,11 s1=s1+xiaoa(j)*a(i,j) end do xiaoa(i)=-s1/a(i,i)+f1(i)/a(i,i)+xiaoa(i)-xiaoa(12)*a(i,12)/a(i,i)-xiaoa(13)*a(i,13)/a(i,i)end do!B阵do i=1,11 s1=0 do j=1,11 s1=s1+xiaob(j)*b(i,j) end do xiaob(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,11 sita(j)=atand(xiaob(j)/xiaoa(j)+180 r(j)=sqrt(xiaoa(j)*2+xiaob(j)*2) g(j)=V0(j)+u(j)+sita(j) h(j)=r(j)/f(j)end dodo i=1,11 do while(g(i)360.or.g(i)360)then do g(i)=g(i)-360 if(g(i)0.and.g(i)360)exit end do else end if if(g(i)0.and.g(i)230.or.gg1230)then do gg1=gg1-360 if(gg1-130)exit end do else end if if(gg1-130.and.gg1230.or.gg2230)then do gg2=gg2-360 if(gg2-130)exit end do else end if if(gg2-130.and.gg2230)exit end do else end ifend dog(12)=g(4)-0.075*(g(4)-g(3)g(13)=g(9)+0.081*(g(9)-g(8)sita(12)=-(u(12)+V0(12)-g(12)sita(13)=-(u(13)+V0(13)-g(13)h(12)=h(4)*0.324 h(13)=h(9)*0.282do j=12,13 r(j)=h(j)*f(j) xiaoa(j)=r(j)*cosd(sita(j) xiaob(j)=r(j)*sind(sita(j)end do g1=g-g0 h1=h-h0 i1=i1+1if(all(abs(h1)10.0).and.all(abs(g1)caita(j-1).and.caita(j)caita(j+1)then i1=i1+1 !高潮个数 end if if(caita(j)caita(j-1).and.caita(j)caita1(j-1).and.caita1(j)caita1(j+1)then i1=i1+1 do t=j-1,j+1,0.01 s=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(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)+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)then ma=sa hightt(i1)=int(t-floor(t)*60) hight(i1)=floor(t)+371 end if end do hightide(i1)=maelseend if!低潮潮位及时刻if(caita1(j)caita1(j-1).and.caita1(j)sa)then mi=sa lowtt(i2)=int(t-floor(t)*60) lowt(i2)=floor(t)+371 end if end do lowtide(i2)=mielseend if end do!将高潮位写入hightide.txt，第1列为潮位，第2、3列为时刻open(unit=2,file=hightide.txt) do i=1,i1 write(2,*)hightide(i),hight(i),hightt(i) end doclose(2)do i=1,i1 j=1 if(hight(i)23) hight(i)=hight(i)-24 j=j+1 hightrq(i)=j end 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.txt，第1列为潮位，第2、3列为时刻open(unit=2,file=lowtide.txt) do i=1,i2 write(2,*)lowtide(i),lowt(i),lowtt(i) end doclose(2)do i=1,i2j=1 if(lowt(i)23) lowt(i)=lowt(i)-24 j=j+1 lowtrq(i)=j end doend doopen(unit=2,file=lowtidexiu.xls)do i=1,i2write(2,*)lowtide(i),lowtrq(i),lowt(i),lowtt(i) !将低潮位及时