潮汐调和分析
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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