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