tide.f90

1968年1月青岛验潮站逐时水位观测数据,进行潮汐调和分析和预报 !选入的分潮(18项) 数据是从1968年1月1日0时开始的

  end do
end do
!输出FF_GG
open(1,file='FF_GG.txt')
write(1,*),'FF'	
do j=1,P
  write(1,*),j
  write(1,*),FF(j)
enddo
write(1,*),'GG'
do j=1,P
  write(1,*),j
  write(1,*),GG(j)
enddo
close(1)

!随从分潮的和赋值
!振幅比
!P1/K1=0.274
kk(1)=0.274
!K2/S2=0.249	
kk(2)=0.249
!迟角差
fai(1)=50*(-0.075)*pi
fai(2)=50*0.081*pi



!对AA*X=BB系数矩阵AA，BB赋值
AA=0	
BB=0

!AA(a,b)

!page48第1个方程中
AA(1,1)=A(0,0)	!S0系数
do j=1,P+Q
AA(1,j+1)=ks(j)*A(0,j) !aj系数
AA(1,P+Q+j+1)=yt(j)*A(0,j)  !bj系数
end do
!AA的第一行2(P+Q)+1个元素,1个方程

!page48第2，3方程中
do i=1,P	
 AA(i+1,1)=ks(i)*A(0,i)	!第2方程S0系数	
 AA(P+i+1,1)=yt(i)*A(0,i)!第3方程S0系数
  do j=1,P+Q
    AA(i+1,j+1)=ks(i)*ks(j)*A(i,j)+yt(i)*yt(j)*B(i,j)	!第2方程aj系数	
	AA(i+1,P+Q+j+1)=ks(i)*yt(j)*A(i,j)-yt(i)*ks(j)*B(i,j)!第2方程bj系数
	!2到P+1,共P行，每行2(P+Q)+1个元素

    AA(P+i+1,j+1)=yt(i)*ks(j)*A(i,j)-ks(i)*yt(j)*B(i,j)   !第3方程aj系数
    AA(P+i+1,P+Q+j+1)=yt(i)*yt(j)*A(i,j)+ks(i)*ks(j)*B(i,j) !第3方程bj系数
	!2到P+1,共P行，每行2(P+Q)+1个元素

	!AA(1:2P+1;1:2(P+Q)+1) 已赋值
  end do
end do

!BB
!page48第1方程中
BB(1)=FF(0)

do i=1,P
!page48第2，3方程中
BB(i+1)=ks(i)*FF(i)-yt(i)*GG(i)!第2方程中
BB(P+i+1)=yt(i)*FF(i)+ks(i)*GG(i)!第3方程中
end do
!BB(1:2P+1)已赋值


!对于差比关系，可写出2Q个方程，其中的系数，BB=0
!AA的值由page49(1.4.5)方程的系数给出
!由(1.4.5)的第一个方程
k=1
do j=1,size(kk) !j=P
    AA(2*P+1+k,1+P+j)=-1 
    AA(2*P+1+k,1+sc(j))=kk(j)*cos(fai(j)) 
    AA(2*P+1+k,1+P+Q+sc(j))=-kk(j)*sin(fai(j))
	k=k+2
end do
!由(1.4.5)的第二个方程 
o=2
do j=1,size(kk) !j=P
    AA(2*P+1+o,1+P+Q+P+j)=-1
    AA(2*P+1+o,1+sc(j))=kk(j)*sin(fai(j))
    AA(2*P+1+o,1+P+Q+sc(j))=kk(j)*cos(fai(j))
	o=o+2
end do

!输出系数矩阵
!AA
open(1,file='系数矩阵AA_BB.txt')	
do i=1,2*(p+q)+1
  do j=1,2*(p+q)+1     
  write(1,'(f20.8)',advance='no'),AA(i,j)
  if (j/(2*(P+Q)+1)==1)then
  write(1,*)
  end if
  enddo
enddo
write(1,*),'BB'
do i=1,2*(p+q)+1
write(1,*),BB(i)
enddo
close(1)

!call Gauss_Jordan法求解线形方程组模块
call Gauss_Jordan(AA,BB,X)	!AA*X=BB求解
!------------------------------------------------------------------------------


!第八步
!解X的赋值
!每个分潮的a,b
S=X(1)	!平均水位	
do i=1,P+Q
a1(i)=X(i+1)
b1(i)=X(P+Q+i+1)
end do
!输出分潮的a1,b1
open(1,file='分潮的a1,b1.txt')
write(1,*),'a1'
do i=1,P+Q
   write(1,*),i,a1(i)
enddo
write(1,*),'b1'
do i=1,P+Q
   write(1,*),i,b1(i)
enddo
close(1)
!---------------------------------------------------------


!第九步

!计算调和常数
do i=1,P+Q
H(i)=sqrt((a1(i))**2+(b1(i))**2)
g(i)=atan(b1(i)/a1(i))
if(a1(i)*cos(g(i))<0)then	!atan值域为(-π/2,π/2)，if语句用来判断g是为arctan(b/a)还是arctan(b/a)+π
g(i)=g(i)+3.14159265
end if
if(g(i)<0)then
g(i)=g(i)+2*3.14159265
end if				
end do
!输出分潮的调和常数H,g
open(1,file='调和常数.txt')	
write(1,*)'H(cm)           g(角度)            '
do i=1,P+Q
write(1,'(2f12.6)')H(i),g(i)/pi
enddo
close(1)



!用四大主要分潮和两个浅水分潮回报1968年1月的逐时水位

!第十步
!求744个逐时水位数据 -360,383
!do i=-360,383	
!hh(i)=S  !S为平均水位s0
  !do j=1,6  !六个分潮 
   ! hh(i)=hh(i)+f(yb(j))*H(yb(j))*cos(w(yb(j))*i+V(yb(j))-g(yb(j)))
  !end do
!end do
!输出逐时水位数据，间隔1hour,潮高四舍五入到cm
!open(1,file='逐时水位数据.txt')	
!write(1,*)'每小时一次的潮高（cm）'
!do i=1,31
   !do j=0,23
     ! write(1,'(a3,i2,a2,i2,a2,a4,f4.0)'),'1月',i,'日',j,'时','潮高',hh((i-1)*24+j-360)
   !enddo
!enddo
!close(1)
!open(1,file='分析用逐时水位.txt')
!do i=-360,383
!write(1,fmt="(F4.0)"),hh(i)
!end do
!close(1)

!open(1,file='分析比较原始水位逐时水位.txt')
!do i=-360,383
!write(1,fmt="(F4.0,5X,F4.0)"),odat(i+361),hh(i)
!end do
!close(1)




!第十一步
!求每天高低潮潮高和潮时

!求44580个逐分水位数据 -21600,22980
!do i=-21600,22980	
!hm(i)=S  !S为平均水位s0
 ! do j=1,6  !六个分潮 
  !  hm(i)=hm(i)+f(yb(j))*H(yb(j))*cos(w(yb(j))*(real(i)/60.)+V(yb(j))-g(yb(j)))
 ! end do
!end do


!j=1		
!k=1
!do i=-21599,22979	
  ! if(hm(i)>hm(i-1).and.hm(i)>hm(i+1))then
   !  hmh(j)=hm(i)	!高潮附近水位
!	 th(j)=i+21600	!高潮附近时(min)
!	 j=j+1
   !else if(hm(i)<hm(i-1).and.hm(i)<hm(i+1))then
    ! hml(k)=hm(i)	!低潮附近水位
	 !tl(k)=i+21600	!低潮附近时(min) 
	 !k=k+1
   !end if
!end do


!高潮个数
!flagh=0 
!do i=1,100
! if (th(i)/=0) then
 ! flagh=flagh+1
! endif
!enddo


!低潮个数
!flagl=0 
!do i=1,100
! if (tl(i)/=0) then
 ! flagl=flagl+1
 !endif
!enddo



!输出未插值高潮时，高潮水位
!allocate(twh(flagh),hwh(flagh),twl(flagl),hwl(flagl))
!do i=1,flagh
  ! twh(i)=th(i)
   !hwh(i)=hmh(i)
!enddo
!do i=1,flagl
   !twl(i)=tl(i)
   !hwl(i)=hml(i)
!enddo


!open(1,file='未插值高潮时 高潮水位.txt')	
!write(1,*)'高潮时             ','潮高(cm)'
!时间转化为日期
!k=24*60
!do i=1,31
   !do j=1,flagh
     ! if(twh(j)>=0.and.twh(j)<k) then
     ! write(1,'(a3,i2,a2,i2,a2,i2,a2,a5,f4.0)'),'1月',i,'日',floor(twh(j)/60),'时',int(twh(j)-floor(twh(j)/60)*60),'分','   ',hwh(j)
      !endif
   !enddo
   !twh=twh-24*60
!enddo
!close(1)

!输出未插值低潮时，低潮水位
!open(1,file='未插值低潮时 低潮水位.txt')	
!write(1,*)'低潮时             ','潮高(cm)'
!k=24*60
!do i=1,31
  ! do o=1,flagl
    !  if(twl(o)>0.and.twl(o)<k) then
   !   write(1,'(a3,i2,a2,i2,a2,i2,a2,a5,f4.0)'),'1月',i,'日',floor(twl(o)/60),'时',int(twl(o)-floor(twl(o)/60)*60),'分','   ',hwl(o)
    !  endif
   !enddo
   !twl=twl-24*60
!enddo
!close(1)

!deallocate(twh,hwh,twl,hwl)
!拉格朗日四次差值多项式求高低潮潮高和潮时p105-106

!分配数组
!allocate(tczh(flagh),tczl(flagl),hczh(flagh),hczl(flagl))
!allocate(cz(flagl+flagh))

!高低潮合并到一起
!do i=1,flagh
  ! cz(i)=th(i)
!enddo
!do i=j,flagl+flagh
  ! cz(i)=tl(i-j+1)
!enddo
!cz=cz-21600

!牛顿法求解
!allocate(AAA(size(cz)),BBB(size(cz)),CCC(size(cz)),DDD(size(cz)),ttt(size(cz)),tr(size(cz)),hcz(size(cz)),tcz(size(cz))) 
!do i=1,size(cz)
  !do j=-2,2
   ! ycz(i,j)=hm(cz(i)-j)
  !enddo

 ! AAA(i)=(ycz(i,1)+ycz(i,-1)-2*ycz(i,0))/6.
  !BBB(i)=(ycz(i,2)+ycz(i,-2)-2*ycz(i,0))/24.
  !CCC(i)=(ycz(i,1)-ycz(i,-1))/6.
  !DDD(i)=(ycz(i,2)-ycz(i,-2))/12.
     !  ttt(i)=ycz(i,0)
	 !  tr(i)=ttt(i)+1
  !do while (abs(ttt(i)-tr(i))>0.1e-20)
   !  tr(i)=ttt(i)
   !  ttt(i)=tr(i)-(4*AAA(i)*tr(i)**3+3*BBB(i)*tr(i)**2+2*CCC(i)*tr(i)+DDD(i))/(12*AAA(i)*tr(i)**2+6*BBB(i)*tr(i)+2*CCC(i))
  !enddo
 ! tcz(i)=cz(i)+ttt(i)

 ! hcz(i)=AAA(i)*ttt(i)**4+BBB(i)*ttt(i)**3+CCC(i)*ttt(i)**2+DDD(i)*ttt(i)+hm(cz(i))
!enddo


!将高低潮分开
!tcz=(tcz+21600)
!do i=1,flagh
  ! tczh(i)=tcz(i)
   !hczh(i)=hcz(i)
!enddo
!do i=flagh+1,flagl+flagh
  ! tczl(i-flagh)=tcz(i)
  ! hczl(i-flagh)=hcz(i)
!enddo


!潮位数据，潮时精确到min,潮高四舍五入到cm



!输出高潮时，高潮水位

!open(1,file='高潮时 高潮水位.txt')	
!write(1,*)'高潮时             ','潮高(cm)'
!时间转化为日期
!k=24*60
!do i=1,31
   !do j=1,flagh
     ! if(tczh(j)>=0.and.tczh(j)<k) then
      !write(1,'(a3,i2,a2,i2,a2,i2,a2,a5,f4.0)'),'1月',i,'日',floor(tczh(j)/60),'时',int(tczh(j)-floor(tczh(j)/60)*60),'分','   ',hczh(j)
      !endif
   !enddo
   !tczh=tczh-24*60
!enddo
!close(1)

!输出低潮时，低潮水位
!open(1,file='低潮时 低潮水位.txt')	
!write(1,*)'低潮时             ','潮高(cm)'
!k=24*60
!do i=1,31
  ! do j=1,flagl
    !  if(tczl(j)>0.and.tczl(j)<k) then
    !  write(1,'(a3,i2,a2,i2,a2,i2,a2,a5,f4.0)'),'1月',i,'日',floor(tczl(j)/60),'时',int(tczl(j)-floor(tczl(j)/60)*60),'分','   ',hczl(j)
    !  endif
  ! enddo
   !tczl=tczl-24*60
!enddo
!close(1)




!第十二步
!使用Matlab分析资料

print*,"请到源程序文件夹查看结果"
print*,"作者 04 海洋科学 单士亮 [e-mail]"
print*,"程序我已检查了多次，鉴于此程序比较长，其中难免有错误，请指正，谢谢!"
read(*,*)
end program tide
!AUTHOR 04 海洋科学 单士亮
!ALL RIGHT RESERVED