潮汐调和分析

program main

real*8 h(100000),h0(100000),preh(100000),prehf(100000),ph0(100000),ph(100000) !实测数据

real*8 w(1000) !分潮角速度

real*8 x(1000),y(1000),aa(1000,1000),bb(1000,1000),f1(1000),f2(1000)

real*8 q(1000,1000),d(1000),v(1000,1000),aaa(1000,1000),bbb(1000,1000)

integer n,ntidal,khl,kh1,n2,n12 ,np,nep,np12,nbp !n为观测数据时间长度!!ntidal为分潮数目

real*8 tdif,dt,err,pai,tbi,max,min

character*15 b,char

WRITE(*,111)

WRITE(*,222)

WRITE(*,111)

111 format(20x,'********************')

222 FORMAT(20x,'**<<潮汐调和分析>>**')

pai=3.1415926

dt=1.d0

tdif=8.d0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!读取预报时段实测潮位(连云港，1979年7月)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

open(1,file='d:\Lyg1979.dat')

read(1,*)char,icck,iyyk,immk,iddk,ihhk,iccj,iyyj,immj,iddj,ihhj

call GDAY(iddk,immk,iyyk,icck,kd) !开始时间转换为18521015后的天数

nbp=(kd-1)*24+ihhk-tdif

call GDAY(iddj,immj,iyyj,iccj,kd) !开始时间转换为18521015后的天数

nep=(kd-1)*24+ihhj-tdif

np=nep-nbp+1

np12=np/12

nep=12*np12+nbp-1

np=nep-nbp+1

read(1,*)

do i=1,np12

read(1,*)(ph0(j),j=(i-1)*12+1,12*i)

enddo

close(1)

call GDAY(01,07,79,19,kd)

nmlb=(kd-1)*24+00-tdif

call GDAY(31,07,79,19,kd)

nmle=(kd-1)*24+23-tdif

nmn=nmle-nmlb+1

write(*,*)np,nmn,nmle-nmlb-1+i

do i=1,nmn

ph(i)=ph0(nmlb-nbp-1+i)

end do

write(*,*)(ph(i),i=1,nmn) end

潮汐的变化规律

潮汐的变化规律 由于太阳与月亮对地球的引力作用，我国大部分沿海地区均有一昼夜各出现海水涨落两次的潮汐现象。每月的农历初一至初五（或农历十六至二十）为大潮汐（当地人称“大活汛”）；农历初六至十二（或农历二十一至农历二十五）为小潮汐（当地人称“死汛”）；而初九或二十四为最小潮（当地人称“死汛底”）。每天的潮汐时间均后延45分钟左右，如此周而复始 有个计算公式共，仅供大家参考。 满潮时间=（农历日—1或16）乘以0.8+10:32 干潮时间=满潮时间加或减6:12 潮汐表编辑 潮汐预报表的简称。它预报沿海某些地点在未来一定时期的每天 潮汐情况。在航运方面，有些水道和港湾须在高潮前后才能航行和进出港；在军事方面，有时为了选择有利的登陆地点和时间，就必须考虑和掌握潮汐的情况；在生产方面，沿海的渔业、水产养殖业、农业、盐业、资源开发、港口工程建设、测量、环境保护和潮汐发电等，都要掌握潮汐变化的规律。潮汐表就是为这些方面服务的。 中文名 潮汐预报表 外文名

Tidal prediction table 作用 预报沿海某些地点潮汐情况 服务行业 航运，军事，生产... 最早文献 《海涛志》 包括 主港逐日预报表，附港差比数等 目录 1简介 2文献来源 3港差比数 4潮汐信息 5简便算法 6潮汐时间 1简介编辑 cháo xī biǎo 潮汐表 tide tables 潮汐表又称潮汐长期预测表，即在正常天气情况下由天文因素影响所

产生的潮汐。 2文献来源编辑 英国开尔文 中国唐代窦叔蒙在《海涛志》一文中提出了根据月相推算高潮时刻的图表法，这是保存下来的介绍潮汐预报方法的最早的文献，大约比英国的《伦敦桥潮候表》早400年。19世纪60年代末，英国开尔文和G.H.达尔文等人提出了潮汐调和分析方法，后来还设计和制造了机械的潮汐推算机，使潮汐表的编算工作得到迅速发展。自20世纪60年代以来，电子计算机已广泛应用在潮汐推算工作中。 潮汐表一般包括主港逐日预报表（通常有高潮和低潮的时间和潮高，有的港还有每小时的潮高）、附港差比数、潮信和任意时刻的潮高计算等内容。 主港逐日预报表 潮汐现象可视为由许多不同周期的分潮叠加而成，故任意时刻的潮高可表示为 图片中A为平均海平面在潮高基准面上的高度，表示分潮的圆频率，为交点因子，d为格林威治开始时的天文相角，H和为分潮的调和常数──振幅和迟角。这样，应用已求出的该港的潮汐调和常数，就能

潮汐要素复习整理

潮汐原理复习思考题整理 （第四章~第五章） 第四章 1.什么是中期观测资料分析和短期观测资料分析，以及调和常数求解的实际步骤 中期观测资料分析：属于不同群的分潮的会合周期最长为1个月，因此把长度长于一个月但不足一年的观测记录称为中期观测资料 短期观测资料分析：观测的时间长度只有一天或几天 调和常数求解的实际步骤： ?中期观测资料分析(TB P103-107) 1）区分主分潮和随从分潮2）取L 段观测记录，式（4.4）可以写为（4.6） 3）将式（4.6）的余弦函数展开得到（4.7） 4）式（4.7）是包含2（P+Q）+1个未知数的由 () 1 L l l M = ∑ 个方程组成的矛盾方程组 5）通过最小二乘法得到矛盾方程组的法方程(4.10) 6）当L=1时，法方程(4.10)变为TB P106 7）引入Q个随从分潮与相应的主分潮的差比关系后，将给出另外2Q个方程(4.11) 8）进一步求得(4.12) ?短期观测资料分析(TB P116-119) 1）潮汐调和常数的初算2）潮流调和常数的计算 3）噪声方差的估计4）不合理数据的舍弃 5）调和常数和余流的计算6）潮流椭圆要素的计算 2.短期资料观测引入的参数D 和d 代表什么含义，具有什么作用？ 振幅系数D 和迟角订正d 用准调和分潮表达式比用调和分潮表达式要简单的多，不但可以简化许多分析过程，而 且对分析实际潮汐特征也能使得问题变得更容易。 3.什么是准调和分潮，它和调和分潮有什么区别 ?实际准调和分潮的振幅和相角与A 小时前的引潮力准调和分潮相应量有关，与其余时刻，特别是与当时引潮力则没有关系，故A 叫做准调和分潮的潮龄 ?区别 4.了解潮汐和潮流的自报TB P119 第五章 1.潮汐特征值的含义TB P120-121 2.对于不同潮汐类型港口潮汐特征值的计算方法

海洋要素计算(潮汐)

海洋要素计算作业之二——潮汐（威海2013年五月份） 一．本次潮汐调和分析共选取了十三个分潮： MSf，Q1，O1，K1，P1，K2，N2，M2，S2，MK3，M4，MS4，M6 为使您查看方便，将本次大作业的放在本文件夹各文件内，具体参考如下： 1.原数据为：qd.dat; 2.Fortran编程见该文件夹内：tide.f90文件； 3.求各分潮调和常数H、g的值及其中间过程得到的各值见：qd_tide.dat文件；二．对比回报值和实测值： 1. 回报1968年一月份的水位值见：huibao.dat； 2. 用matlab绘制的潮汐过程曲线见：潮汐过程曲线.bmp 3. 用给定的六个分潮求得的高潮和低潮发生的时刻及潮位值见—：gaodichao.dat; 运行tide.f90后求得威海地区2013年5月份的平均潮差。 由图可知：由于只计算了一个月的潮汐数据，所以回报值和实测值相符的不是很好，如果计算一年的数据，应该会取得比较良好的结果。

三．程序 %% 潮汐过程曲线图 clear,clc %% huibao=load('G:\chaoxi\huibao.dat'); % huibao=fread(fhuibao); shice=load('G:\chaoxi\qd.dat'); % shice=fread(fshice); %huibao_y=zeros(1,12*62); %shice_y=zeros(1,12*62); huibao=double(huibao'); huibao_y=double(huibao(:)); %shice_y=reshape(shice',1,[]) %for i=1:12; % for j=1:62 % huibao_y(i)=huibao(i,j) % shice_y(i)=shice(i,j) %end %end shice=double(shice'); shice_y=double(shice(:)); x=linspace(1,31,length(huibao_y)); plot(x,huibao_y,'r-') hold on plot(x,shice_y,'b-') title('威海(37°31′N ,122°08′E)2013年五月潮汐调和分析图') legend('回报值','实测值') xlabel('时间（2013年五月份）') ylabel('水位（m）')

潮汐调和分析-实验报告

中期水位观测资料的最小二乘分析报告 摘要： 本次实验中采用了KM站（28.05N,121.17E）1997年8月的逐时潮位数据，运用中期水位观测资料的最小二乘分析方法，通过奇异值剔除、调和分析、逐时潮位回报、高低潮及余水位计算等工作，对此验潮站的数据进行了初步分析，并了解了中期水位资料分析的基本流程。 报告主要内容： （1）11个主要分潮（MSf Q1 O1 K1 N2 M2 S2 MK3 M4 MS4）及两个随从分潮（P1 K2）的调和常数H和g （2）图像和数据文件的基本信息 （3）平均潮差和潮汐类型 （4）余水位特征分析 （5）误差分析 （6）程序的相关说明 （1）潮汐的调和常数： 利用最小二乘原理，通过引入差比关系的方法，我们可以成功得到11个主要分潮和2个随从分潮的调和常数如下： 分潮名称调和常数H 调和常数g MSf 121.2222 -32.38747 Q1 62.95736 233.5120 O1 225.5294 238.7111 K1 266.1612 113.2537 N2 420.5689 167.2492 M2 1922.772 174.8581 S2 679.3940 197.3759 MK3 33.19594 252.1002 M4 32.43390 121.7806 MS4 33.60584 198.8826 M6 3.762754 94.29744 P1 73.46050 109.5160 K2 192.9479 201.4156 程序运行结果如图：

其中H关系到分潮的振幅，g关系到分潮的相位。从表中可以看出，M2分潮的振幅最强，对当地潮位的贡献最大，这与实际情况相符，但K1分潮的调和常数H仅有266.1412，结果偏小。 （2）图像和数据文件的基本信息： 本次报告中包含以下数据文件： 1. KM9708new.dat 数据原始文件。 2. KM9708new_02.dat 经过奇异值订正的数据文件，为方便画图时读取，没有输入数据质量信息。 3. 调和常数.txt 保存了调和常数的相关数据 4.回报值.txt 保存了利用六个主要分潮进行数据回报得到的结果，同样为方便读取，没有输入质量控制信

潮汐知识点

潮汐知识点汇总： 1.中版潮汐表一共有六册，其中第1-3册是我国海区的，4-6册是国外海区的。潮汐表中的潮时是当地标准时。 2.中版潮汐表的主要内容：主港潮汐预报表、潮流预报表、差比数和潮信表、格林尼治月中天时刻表、梯形图卡。梯形图卡作用是求取任意潮时的潮高和任意潮高的潮时。我国一些重要港口每小时潮高可以在主港潮汐预报表里面查找。 3.差比数是指主附港之间的潮时差、潮差比和改正值；潮差比是附港平均潮差与主港平均潮差之比（对于日潮港来说，潮差比是附港回归潮大的潮差与主港回归潮大的潮差之比）；潮时差是主港和附港高（低）潮时之差，正号表示附港比主港发生的晚，负号表示早。 4.下面因素的急剧变化会导致潮汐变化的反常：台风（增水）、寒潮（减水）、春季气旋入海、降水、气压、结冰。（注意：海啸不会导致潮汐变化的反常） 5.中版潮汐表预报潮高的误差是20-30CM，预报潮时的误差是20-30分钟。 6.利用中版潮汐表查主港潮汐，可从目录查主港所在页码。 7.附港潮时的计算公式：附港潮时=主港潮时+潮时差 附港潮高计算，对于中国海区的1-3册，可以利用下面公式: 当主附港季节改正数较大时（SC>10CM）： 附港潮高=（主港潮高-（主港平均海面+主港平均海面季节改正））*潮差比+（附港平均海面+附港平均海面季节改正） 当主附港季节改正不大时(SC<10CM)：附港潮高=主港潮高*潮差比+改正值 对于外国海区的4-6册，可以利用下面公式: 附港潮高=主港潮高*潮差比+改正数+潮高季节改正数 8.利用潮信资料算潮汐：潮信资料包括平均大（小）潮升，平均高（低）潮间隙，平均海面。当地高（低）潮潮时=格林威治月上（下）中天时+当地高（低）潮间隙 对于半日潮港，可用下列方法近似求月中天时间： 上半月：月上中天时=（农历日期-1）*0.8+1200月下中天时=月上中天时±1225 下半月：月上中天时=（农历日期-16）*0.8月下中天时=月上中天时±1225 平均大潮高潮高(MHWS)=大潮升(SR) 平均大潮低潮高MLWS=2*MSL-SR 平均小潮高潮高（MHWN）=小潮升(SN) 平均小潮低潮高MLWN=2*MSL-SN 大潮日指的是初三、十八；小潮日是十一、二十五。 其他日潮高估算： 高潮潮高=大潮升-（大潮升-小潮升）/7.5*（所求日与大潮日间隔天数） 低潮潮高=2*MSL-高潮高 9.求任意潮时的潮高和任意潮高的潮时 T=涨落潮时间间隔t是任意时与高低潮的时间间隔 10.实际水深=海图水深+潮高+（CD-TD）=测水仪测的水深+船舶吃水 过浅滩最小潮高：海图水深+潮高+（CD-TD）>=船舶吃水+富余水深 过高架最大潮高：船舶水面是最大高度+潮高+安全余量>=大潮升+净空高度 其中船舶水面最大高度=主大桅高+型深—吃水

一个简化的潮汐预报准调和分析方法

一个简化的潮汐预报准调和分析方法 王如云1，2,李慧娟1，2，蒋风芝2 （1 河海大学水文水资源与水利工程科学国家重点实验室，南京210098； 2 河海大学海洋学院，南京210098） 摘要：在用现有的浅水港日潮汐准调和分析预报方法进行潮汐分析预报时，发现最小二乘法的法方 程组的系数矩阵条件数很大，数量级在108，因此矩阵是坏条件的(或为病态的)，算法不稳定。根据 潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建 立了一个简化的浅水准调和分析模型。利用连云港的多年实测数据检验，简化的准调和分析模型相 对于原准调和分析模型来讲，最小二乘法的法方程组系数矩阵条件数小很多，因此简化后的模型计 算更为稳定。在实测数据时间较长的情况下，简化前后的模型预报精度相当。但当实测数据较短时， 简化前的原模型却没有传统的调和分析模型的预报结果精度高，而简化后的模型却能保持比传统的 调和分析模型的预报结果有一定的改善。特别是简化后比简化前的模型计算时间减少了68%。 关键词：浅水潮汐；准调和分析；潮汐预报 1引言 在潮汐预报方面，一般采用调和分析方法，在深水区域此方法可以获得很好的预报效果，但在浅水区域尤其是河口区域，由于浅水潮汐的复杂性，采用此方法往往不能获得满意的效果。例如杜德森提出的60个分潮[1]，其结果不能令人满意。为此，杜德森后来又提出了一个直接对高低潮进行浅水改正的方法[2]，该方法虽然使高低潮的预报精度有了提高，但把它应用到逐时潮位预报上则有许多困难和不便之处。在浅水区域由于非线性效应的加大，潮波往往产生畸变。此时，高频振动的作用必须予以充分考虑。为了提高浅水区域潮汐预报的精度，从调和分析方法来讲就必须增加高频的浅水分潮。在水深不太浅的区域，浅水分潮的振幅会随着阶数的增高而迅速减小，所以在一般港口采用较少数目的主要浅水分潮即可满足潮汐预报的要求。但在浅水区，常常需要考虑到六阶甚至更高阶的相互作用，才能满足潮汐预报的要求。上个世纪六十年代，一些潮汐学者试图通过扩充高频分潮的数目以使预报结果获得改进，如Zelter and Cumimngs[3]以及Rossiter and Lennon[4]曾将分潮的数目扩充到110多个，但效果并不理想。方国洪等人认为，不理想的原因在于随着频率的增加，高频分潮的数目极速的增加，不可能从中挑选出少数分潮近似代替所有分潮，难以用有限数目的浅水分潮来体现总的浅水效应。可以认为通过增加浅水分潮以改进潮汐预报，其效果可能是比较有限的。 基于以上分析，方国洪等提出了一个浅水潮汐预报的准调和分析方法[5]，可以用来推算任意时刻的潮高，也可以用来推算高、低潮，效果比传统的调和分析法有了显著的改进。但我们使用此方法对连云港的多年潮位实测数据进行分析预报时，发现最小二乘法的法方程组系数矩阵条件数很大，算法不稳定。为此，我们对浅水准调和分析模型进行了简化，简化后的模型计算更为稳定，计算时间大为减少。 2 准调和分析方法介绍 方国洪等人[5]提出的浅水预报准调和方法思路是把潮高分做两部分，一部分为低频部分，由 基金项目：水文水资源与水利工程科学国家重点实验室开放研究基金(2005407411)；中国教育部科学技术研究重点项目（104104）；中国江苏省普通高等学校高新技术产业发展项目(JH03-010) 作者简介：王如云(1963-)，教授，男，安徽芜湖人，从事计算物理学研究，E-mail:wangry@https://www.wendangku.net/doc/806611098.html,

潮汐调和分析

program main real*8 h(100000),h0(100000),preh(100000),prehf(100000),ph0(100000),ph(100000) !实测数据 real*8 w(1000) !分潮角速度 real*8 x(1000),y(1000),aa(1000,1000),bb(1000,1000),f1(1000),f2(1000) real*8 q(1000,1000),d(1000),v(1000,1000),aaa(1000,1000),bbb(1000,1000) integer n,ntidal,khl,kh1,n2,n12 ,np,nep,np12,nbp !n为观测数据时间长度!!ntidal为分潮数目 real*8 tdif,dt,err,pai,tbi,max,min character*15 b,char WRITE(*,111) WRITE(*,222) WRITE(*,111) 111 format(20x,'********************') 222 FORMAT(20x,'**<<潮汐调和分析>>**') pai=3.1415926 dt=1.d0 tdif=8.d0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !读取预报时段实测潮位(连云港，1979年7月) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! open(1,file='d:\Lyg1979.dat') read(1,*)char,icck,iyyk,immk,iddk,ihhk,iccj,iyyj,immj,iddj,ihhj call GDAY(iddk,immk,iyyk,icck,kd) !开始时间转换为18521015后的天数 nbp=(kd-1)*24+ihhk-tdif call GDAY(iddj,immj,iyyj,iccj,kd) !开始时间转换为18521015后的天数 nep=(kd-1)*24+ihhj-tdif np=nep-nbp+1 np12=np/12 nep=12*np12+nbp-1 np=nep-nbp+1 read(1,*) do i=1,np12 read(1,*)(ph0(j),j=(i-1)*12+1,12*i) enddo close(1) call GDAY(01,07,79,19,kd) nmlb=(kd-1)*24+00-tdif call GDAY(31,07,79,19,kd) nmle=(kd-1)*24+23-tdif nmn=nmle-nmlb+1 write(*,*)np,nmn,nmle-nmlb-1+i do i=1,nmn

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

中期水位资料对潮汐进行调和分析 !利用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

潮汐调和分析-实验报告

中期水位观测资料的最小二乘分析报告摘要： 本次实验中采用了KM站（28.05N,121.17E）1997年8月的逐时潮位数据，运用中期水位观测资料的最小二乘分析方法，通过奇异值剔除、调和分析、逐时潮位回报、高低潮及余水位计算等工作，对此验潮站的数据进行了初步分析，并了解了中期水位资料分析的基本流程。 报告主要内容： （1）11 个主要分潮（MSf Q1 O1 K1 N2 M2 S2 MK3 M4 MS4）及两个随从分潮（P1 K2）的调和常数H和g （2）图像和数据文件的基本信息 （3）平均潮差和潮汐类型 （4）余水位特征分析 （5）误差分析 （6）程序的相关说明 （1）潮汐的调和常数： 利用最小二乘原理，通过引入差比关系的方法，我们可以成功得到11个主要分潮和2个随从分潮的调和常数如下： 分潮名称调和常数H调和常数g MSf 121.2222 -32.38747 Q1 62.95736 233.5120 O1 225.5294 238.7111 K1 266.1612 113.2537 N2 420.5689 167.2492 M2 1922.772 174.8581 S2 679.3940 197.3759 MK3 33.19594 252.1002 M4 32.43390 121.7806 MS4 33.60584 198.8826 M6 3.762754 94.29744 P1 73.46050 109.5160 K2 192.9479 201.4156 程序运行结果如图：

其中H关系到分潮的振幅，g关系到分潮的相位。从表中可以看出，M2 分潮的振幅最强，对当地潮位的贡献最大，这与实际情况相符，但K1分潮的调和常数H 仅有266.1412，结果偏小。 （2）图像和数据文件的基本信息： 本次报告中包含以下数据文件： 1. KM9708new.dat 数据原始文件。 2. KM9708new_02.dat 经过奇异值订正的数据文件，为方便画图时读取，没有输入数据质量信息。 3. 调和常数.txt 保存了调和常数的相关数据 4.回报值.txt 保存了利用六个主要分潮进行数据回报得到的结果，同样为方便读取，没有输入质量控制

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

!利用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