潮汐调和分析

2018.10.29补充

修复了一个BUG(line229)

---

潮汐调和分析是潮汐认为是不同周期的正弦分潮叠加的结果。对于每一个分潮，都是假想天体作用的结果。得知一段时间内的观测数据皆可以拟合出各个分潮。以下是利用最小二乘法将一段时间内的潮汐分潮计算得出并进行预报的FORTRAN程序。本作业只给出了主程序部分，因为是复制的end有没对齐没有检查，不保证正确，欢迎查BUG。

!传统最小二乘法调和分析程序
program main
real*8 h(100000),h0(100000),preh(100000),prehf(100000),ph0(100000),ph(100000),ph1(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  !tdif--与格林威治时差
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!读取分潮角速度
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(1,file='D:\hw3\w63.dat')
ntidal=0
read(1,*)
do i=1,1000
read(1,*) j,w(i)
IF(j==-99) exit
ntidal=ntidal+1
end do
do i=1,ntidal
w(i)=w(i)*pai/180.
enddo
write(*,22) ntidal
22 format('分潮数目=',i3)
close(1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!读取预报时段实测潮位(连云港，1979年7月)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(2,file='D:\hw3\Lyg1979.dat')
read(2,*) char,icc1,iyy1,imm1,idd1,ihh1,iccl,iyyl,imml,iddl,ihhl

call GDAY(idd1,imm1,iyy1,icc1,kd) !时间转换第几天
nbp=(kd-1)*24+ihh1-tdif  !时间转换第几小时
CALL GDAY(IDDL,IMML,IYYL,ICCL,KD)
nep=(KD-1)*24+ihhL-tdif
np=nep-nbp+1 !总共有几小时
np12=np/12

read (2,*)
do i=1,np12
read(2,*) (ph0(j),j=(i-1)*12+1,12*i)
end do
close(2)

CALL GDAY(01,09,79,19,KD) !9月数据
nmlb = (KD-1)*24+00-tdif
CALL GDAY(30,09,79,19,KD)
nmle = (KD-1)*24+23-tdif
nmn  = nmle-nmlb+1

do i=1,nmn
ph(i)=ph0(nmlb-nbp+i)
end do

CALL GDAY(01,03,79,19,KD) !3月数据
nmlb = (KD-1)*24+00-tdif
CALL GDAY(31,03,79,19,KD)
nmle = (KD-1)*24+23-tdif
nmn0  = nmle-nmlb+1
!write(*,*)nmn,np,nmlb-nbp-1+nmn

do i=1,nmn0
ph1(i)=ph0(nmlb-nbp+i)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!读取调和分析时段潮位 (连云港，1979年2-7月)
!计算起止时间n
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(5,file='Lyg1979.dat')
read(5,*) b,icc1,iyy1,imm1,idd1,ihh1,iccl,iyyl,imml,iddl,ihhl

call GDAY(01,01,IYY1,ICC1,KD) !时间转换
KH1=(KD-1)*24+ihh1-tdif
CALL GDAY(31,01,IYYL,ICCL,KD)
KHL=(KD-1)*24+23-tdif
n=khL-kh1+1 !预报使用数据长度
n00=n/12

call GDAY(01,02,IYY1,ICC1,KD) !时间转换
KH1=(KD-1)*24+ihh1-tdif
CALL GDAY(31,07,IYYL,ICCL,KD)
KHL=(KD-1)*24+23-tdif
n=khL-kh1+1 !预报使用数据长度
n12=n/12

write(*,11) n
11 format('时间长度(小时)=',i5)

read (5,*)
do i = 1, n00
read (5,*)
end do
do i=1,n12
read(5,*) (h(j),j=(i-1)*12+1,12*i)
end do

h0=h

do i=1,n
avg=avg+h(i)
end do
avg=avg/n
do i=1,n
h(i)=h(i)-avg !距平化
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!求系数矩阵a，用于计算Z0,x1,...,xJ
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
n2=(n-1)/2
n=n-1         !n为奇数
!write(*,*) n,n2
aa(1,1)=n

do i=2,ntidal+1
aa(i,1)=dsin(n*w(i-1)*dt/2.)/dsin(w(i-1)*dt/2.)   !双精度
aa(i,i)=(n+dsin(n*w(i-1)*dt)/dsin(w(i-1)*dt))/2.
end do

do i=2,ntidal+1
do j=i+1,ntidal+1
aa(i,j)=dsin(n*(w(i-1)-w(j-1))*dt/2.)/dsin((w(i-1)-w(j-1))*dt/2.)
aa(i,j)=(aa(i,j)+dsin(n*(w(i-1)+w(j-1))*dt/2.)/dsin((w(i-1)+w(j-1))*dt/2.))/2.
end do
end do

do i=2,ntidal+1       !对称
do j=1,i-1
aa(i,j)=aa(j,i)
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!求系数矩阵b，用于计算y1,y2,...,yJ
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,ntidal
bb(i,i)=(n-dsin(n*w(i)*dt)/dsin(w(i)*dt))/2.
do j=i+1,ntidal
bb(i,j)=dsin(n*(w(i)-w(j))*dt/2.)/dsin((w(i)-w(j))*dt/2.)
bb(i,j)=(bb(i,j)-dsin(n*(w(i)+w(j))*dt/2.)/dsin((w(i)+w(j))*dt/2.))/2.
end do
end do
do i=1,ntidal
do j=1,i-1
bb(i,j)=bb(j,i)
enddo
enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!求右端向量f1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
f1(1)=0.
do i=1,n  !-n2,n2等价
f1(1)=f1(1)+h(i)
end do
do i=2,ntidal+1
f1(i)=0.
do j=-n2,n2
f1(i)=f1(i)+h(j+n2+1)*dcos(j*w(i-1)*dt)
end do
end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!求右端向量f2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,ntidal
f2(i)=0.
do j=-n2,n2
f2(i)=f2(i)+h(j+n2+1)*dsin(j*w(i)*dt)
end do
end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!计算x矩阵
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call qrdcmp(aa,ntidal+1,ntidal+1,q)
call qrbksb(aa,ntidal+1,q,f1,x)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!计算y矩阵
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call qrdcmp(bb,ntidal,ntidal,q)
call qrbksb(bb,ntidal,q,f2,y)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!回报(连云港，1979年3月)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call GDAY(01,02,IYY1,ICC1,KD) !时间转换
KH10=(KD-1)*24+ihh1-tdif
call GDAY(01,03,IYY1,ICC1,KD) !时间转换
nmlb = (KD-1)*24+00-tdif
do i=1,nmn0
t=(nmlb-KH10+i-n2-1)*dt
preh(i)=x(1)+avg !x(1)是Z0
do j=1,ntidal
preh(i)=preh(i)+x(j+1)*dcos(w(j)*t)
preh(i)=preh(i)+y(j)*dsin(w(j)*t)
end do
end do

open (2001,file='lyg3.dat')
do j=1,nmn0
write(2001,'(f12.3)')ph1(j)
end do
close(2001)

open (2002,file='hindcast.dat')
do j=1,nmn0
write(2002,'(f12.3)')preh(j)
end do
close(2002)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!误差平方和
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
err=0.
do i=1,nmn0
err=err+(preh(i)-ph1(i))**2
end do

err=sqrt(err/nmn0)
write(*,2) err
write(*,3) x(1)+avg
2 format('均方差=',f12.3,'cm')
3 format('平均海平面=',f12.3,'cm')
open(2002,file='result_hindcast.dat')
write(2002,2) err
write(2002,3) x(1)+avg
close(2002)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!预报(连云港，1979年9月)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call GDAY(01,02,IYY1,ICC1,KD) !时间转换
KH10=(KD-1)*24+ihh1-tdif
call GDAY(01,09,IYY1,ICC1,KD) !时间转换
nmlb = (KD-1)*24+00-tdif
do i=1,nmn
t=(nmlb-KH10+i-n2-1)*dt
prehf(i)=x(1)+avg !x(1)是Z0
do j=1,ntidal
prehf(i)=prehf(i)+x(j+1)*dcos(w(j)*t)
prehf(i)=prehf(i)+y(j)*dsin(w(j)*t)
end do
end do

open (2003,file='lyg9.dat')
do j=1,nmn
write(2003,'(f12.3)')ph(j)
end do
close(2003)

open (2004,file='forcast.dat')
do j=1,nmn
write(2004,'(f12.3)')prehf(j)
end do
close(2004)

end

---

为了验证程序这里还给出MATLAB绘图程序

close all; clear; clc;
lyg3 = load( 'E:\物理海洋学\hw3\lyg3.dat' );
hindcast = load( 'E:\物理海洋学\hw3\hindcast.dat' );
figure
h1 = plot( 1 : 1 / 24 : 31.99, lyg3, 'r-', 'LineWidth', 1 );
hold on
h2 = plot( 1 : 1 / 24 : 31.99, hindcast, 'b--', 'LineWidth', 1 );
grid on
xlabel( '时间（天）' );
ylabel( '潮高（cm）' );
legend( '实际值', '回报值' );
axis tight

lyg9 = load( 'E:\物理海洋学\hw3\lyg9.dat' );
forcast = load( 'E:\物理海洋学\hw3\forcast.dat' );
figure
h3 = plot( 1 : 1 / 24 : 30.99, lyg9, 'r-', 'LineWidth', 1 );
hold on
h4 = plot( 1 : 1 / 24 : 30.99, forcast, 'b--', 'LineWidth', 1 );
grid on
xlabel( '时间（天）' );
ylabel( '潮高（cm）' );
legend( '实际值', '预报值' )
axis tight

做出的图像是这样的