中期水位资料对潮汐进行调和分析.docx
《中期水位资料对潮汐进行调和分析.docx》由会员分享，可在线阅读，更多相关《中期水位资料对潮汐进行调和分析.docx（19页珍藏版）》请在冰豆网上搜索。
中期水位资料对潮汐进行调和分析
!
利用1996年7月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报
programwork
implicitnone
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
doi=1,62
read（2,'（16a5）'）aa（i,:
）
enddo
doi=1,62
read（aa（i,5:
16）,*）bb（i,:
）
read（aa（i,3:
4）,*）c（i,:
）
enddo
doi=1,62
c（i,2）=int（real（c（i,2））/10.0）
enddo
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（middletime:
1996-7-16;datasum:
744,middlenumber:
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
doi=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）
enddo
print*
print*,'初始幅角V0:
',V0
!
设caita为潮高数据
doi=1,61
j=-371+（i-1）*12
caita（j:
j+11）=bb（i,:
）
enddo
caita（361:
371）=bb（62,1:
11）
!
计算法方程等式右边的数据，相邻数据时间间隔为1小时
f1（0）=sum（caita）
!
f1为A阵中除第一行外的等式右边一维数据
f1（1:
13）=0
doi=1,13
doj=-371,371
f1（i）=f1（i）+caita（j）*cosd（j*w（i））
enddo
enddo
!
f2为B阵中等式右边的一维数据
f2=0
doi=1,13
doj=-371,371
f2（i）=f2（i）+caita（j）*sind（j*w（i））
enddo
enddo
!
计算A阵中的系数矩阵A
a（0,0）=743
doj=1,13
a（0,j）=sind（743.0/2*w（j））/sind（0.5*w（j））
a（j,0）=a（0,j）
enddo
doj=1,13
a（j,j）=0.5*（743+sind（743.0*w（j））/sind（w（j）））
enddo
doi=1,13
doj=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）
enddo
enddo
print*
print*,'系数矩阵A:
',a
!
计算B阵中的系数矩阵B
doj=1,13
b（j,j）=0.5*（743-sind（743.0*w（j））/sind（w（j）））
enddo
doi=1,13
doj=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）
enddo
enddo
print*
print*,'系数矩阵B:
',b
!
Guass-Seidel迭代法求解方程组
h=0;g=0;i1=0
do
h0=h
g0=g
!
A阵
doi=0,11
s1=0
doj=0,11
s1=s1+xiaoa（j）*a（i,j）
enddo
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）
enddo
!
B阵
doi=1,11
s1=0
doj=1,11
s1=s1+xiaob（j）*b（i,j）
enddo
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）
enddo
!
计算调和常数h,g
doj=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）
enddo
doi=1,11
dowhile（g（i）>360.or.g（i）<0）
if（g（i）>360）then
do
g（i）=g（i）-360
if（g（i）>0.and.g（i）<360）exit
enddo
else
endif
if（g（i）<0）then
do
g（i）=g（i）+360
if（g（i）>0.and.g（i）<360）exit
enddo
else
endif
enddo
enddo
gg1=g（4）-g（3）
gg2=g（9）-g（8）
dowhile（gg1>230.or.gg1<-130）
if（gg1>230）then
do
gg1=gg1-360
if（gg1<230.and.gg1>-130）exit
enddo
else
endif
if（gg1<-130）then
do
gg1=gg1+360
if（gg1>-130.and.gg1<230）exit
enddo
else
endif
enddo
dowhile（gg2>230.or.gg2<-130）
if（gg2>230）then
do
gg2=gg2-360
if（gg2<230.and.gg2>-130）exit
enddo
else
endif
if（gg2<-130）then
do
gg2=gg2+360
if（gg2>-130.and.gg2<230）exit
enddo
else
endif
enddo
g（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.282
doj=12,13
r（j）=h（j）*f（j）
xiaoa（j）=r（j）*cosd（sita（j））
xiaob（j）=r（j）*sind（sita（j））
enddo
g1=g-g0
h1=h-h0
i1=i1+1
if（all（abs（h1）<10.0）.and.all（abs（g1）<2.0））exit!
退出循环条件
enddo
print*
print*,'系数a:
',xiaoa
print*
print*,'系数b:
',xiaob
print*
print*,'迭代循环次数:
',i1
print*
print*,'调和常数h:
',h
print*
print*,'调和常数g:
',g
!
计算平均水位
s0=xiaoa（0）
print*
print*,'平均水位s0:
',s0
!
向文件中输入数据
!
将各分潮的调和常数写入'hg.txt'
open（unit=2,file='hg.txt'）
doi=1,13
write（2,*）h（i）,g（i）
enddo
close
（2）
!
将所有潮汐数据写入'tides.txt'
open（unit=2,file='tides.txt'）
doi=1,62
write（2,'（12i5）'）bb（i,:
）
enddo
close
（2）
!
将所有主要分潮（3,4,8,9）、浅水分潮（5,11）数据按随时间的变化情况写入向量中
doj=-371,371
caita3（j）=s0+f（3）*h（3）*cosd（w（3）*j-sita（3））
caita4（j）=s0+f（4）*h（4）*cosd（w（4）*j-sita（4））
caita8（j）=s0+f（8）*h（8）*cosd（w（8）*j-sita（8））
caita9（j）=s0+f（9）*h（9）*cosd（w（9）*j-sita（9））
caita5（j）=s0+f（5）*h（5）*cosd（w（5）*j-sita（5））
caita11（j）=s0+f（11）*h（11）*cosd（w（11）*j-sita（11））
caita1（j）=caita3（j）+caita4（j）+caita8（j）+caita9（j）+caita5（j）+caita11（j）-5*s0
enddo
!
计算高低潮个数
i1=0;i2=0
doj=-370,370
if（caita（j）>caita（j-1）.and.caita（j）>caita（j+1））then
i1=i1+1!
高潮个数
endif
if（caita（j）i2=i2+1!
低潮个数
endif
enddo
print*
print*,'高潮个数:
',i1
print*,'低潮个数:
',i2
if（allocated（hightide））deallocate（hightide）
if（allocated（lowtide））deallocate（lowtide）
if（allocated（hight））deallocate（hight）
if（allocated（lowt））deallocate（lowt）
allocate（hightide（1:
i1）,lowtide（1:
i2））
!
将观测时的高潮写入文件
open（unit=2,file='realh.txt'）
doi=1,size（hightide）
write（2,*）hightide（i）
enddo
close
（2）
!
将观测时的低潮写入文件
open（unit=2,file='reall.txt'）
doi=1,size（lowtide）
write（2,*）lowtide（i）
enddo
close
（2）
!
潮汐预报部分
!
计算高低潮对应的潮位及时刻
if（allocated（hightide））deallocate（hightide）
if（allocated（lowtide））deallocate（lowtide）
if（allocated（hight））deallocate（hight）
if（allocated（lowt））deallocate（lowt）
if（allocated（lowtrq））deallocate（lowtrq）
if（allocated（hightrq））deallocate（hightrq）
i=min（i1,i2）
allocate（hightide（1:
i1）,lowtide（1:
i2）,hightrq（i1）,lowtrq（i2）,hight（i1）,hightt（i1）,lowtt（i2）,lowt（i2）,chaocha（i））
i1=0;i2=0
doj=-370,370
!
高潮潮位及时刻
if（caita1（j）>caita1（j-1）.and.caita1（j）>caita1（j+1））then
i1=i1+1
dot=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（sma=sa
hightt（i1）=int（（t-floor（t））*60）
hight（i1）=floor（t）+371
endif
enddo
hightide（i1）=ma
else
endif
!
低潮潮位及时刻
if（caita1（j）i2=i2+1
dot=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（s>sa）then
mi=sa
lowtt（i2）=int（（t-floor（t））*60）
lowt（i2）=floor（t）+371
endif
enddo
lowtide（i2）=mi
else
endif
enddo
!
将高潮位写入'hightide.txt'，第1列为潮位，第2、3列为时刻
open（unit=2,file='hightide.txt'）
doi=1,i1
write（2,*）hightide（i）,hight（i）,hightt（i）
enddo
close
（2）
doi=1,i1
j=1
if（hight（i）<=23）then
hightrq（i）=1
endif
dowhile（hight（i）>23）
hight（i）=hight（i）-24
j=j+1
hightrq（i）=j
enddo
enddo
open（unit=2,file='hightidexiu.xls'）
doi=1,i1
hightide（i）=int（hightide（i）+0.5）*1.0
write（2,*）hightide（i）,hightrq（i）,hight（i）,hightt（i）!
将高潮位及时刻写入'hightidexiu.xls'，第1列为潮位，第2、3、4列为天数、小时、分钟
enddo
close
（2）
!
将低潮位写入'lowtide.txt'，第1列为潮位，第2、3列为时刻
open（unit=2,file='lowtide.txt'）
doi=1,i2
write（2,*）lowtide（i）,lowt（i）,lowtt（i）
enddo
close
（2）
doi=1,i2
j=1
if（lowt（i）<=23）then
lowtrq（i）=1
endif
dowhile（lowt（i）>23）
lowt（i）=lowt（i）-24
j=j+1
lowtrq（i）=j
enddo
enddo
op
升级会员 