! 二维潮波数值模拟与预报! 计算海域包括渤海和黄海北部（117.50 -126.67 E,34.00 - 41.00 N）! 空间分辨率是10'x10'! 开边界处的h和g由TMD TOOLBOX求得!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!program hw4implicit noneinterfacesubroutine h_get(h,dx,dy,dt,u,v,d,cd)integer::i,jreal,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cdreal,dimension(43,56),intent(in)::d,u,vreal,dimension(43,56),intent(inout)::hreal,dimension(43),intent(in)::dxend subroutine h_getsubroutine v_get(v,u,cu,cv,dx,dy,dt,d,h,f)integer::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,ureal,dimension(43,56),intent(inout)::vreal,dimension(43),intent(in)::dx,freal,dimension(43,56)::uu,ssend subroutine v_getsubroutine u_get(u,v,cu,cv,dx,dy,dt,d,h,f)integer::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,vreal,dimension(43,56),intent(inout)::ureal,dimension(43),intent(in)::dx,freal,dimension(43,56)::vv,rrend subroutine u_getend interface! < PART 1 >! dx, dy : 网格宽度，单位（m）! dt : 时间步长，取为M2分潮周期的1/240，单位（s）! d : 各个网格点的水深! h : 海水相对静止海面的高度! u, v : 海水水平流速! uu, vv : u、v的平均值! f : 科氏参量! cd : 水深控制场! cu, cv : 速度控制场! h0, g0 : 开边界处的振幅和迟角! h_check : 用于检验迭代是否达到要求的精度! eps : 迭代精度!< PART 2 >! h_3d : 迭代完成后得到的水位! a,b,g,hh : 调和常数! l1,l2,l3,l4,y1,y2 :最小二乘法正规方程组的系数!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!real,parameter::pi=3.14159625,dy=6371004*2*pi/360/6real,parameter::dt=24.8412/2*3600/240,w=360/(240*dt),eps=0.01 !M2 !real,parameter::dt=24.8412*3600/240,w=360/(240*dt),eps=0.01 !K1 integer::i,j,k,m,n,tinteger,dimension(43,56)::cd,cu,cvreal,dimension(43,56)::d,h,u,v,uu,vv,h_checkreal,dimension(43)::dx,freal,dimension(19:54)::h0,g0real,dimension(43,56,240)::h_3dreal::l1,l2,l3,l4,cos_g,sin_greal,dimension(43,56)::a,b,g,hh,y1,y2real,dimension(240)::temp1,temp2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP1:读取水深数据、控制场和开边界点的振幅与迟角!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!open(1,file='BHhai.dat')read(1,'(56f5.0)') (d(i,:),i=1,43)close(1)open(2,file='Bctrlh.dat')read(2,'(56i1)') (cd(i,:),i=1,43)close(2)open(3,file='data.out_M2')!open(3,file='data.out_K1')do i=19,54read(3,'(42x,f6.4,f11.2)') h0(i),g0(i)end doclose(3)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP2: 计算各个纬度的dx和f,计算各个网格点的流速控制场cu和cv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!do i=1,43dx(i)=6371004*cosd(41.0-(i-1)/6)*2*pi/360/6f(i)=2*7.292115e-5*sind(41.0-(i-1)/6)end dodo i=2,42do j=2,55cu(i,j)=cd(i,j)*cd(i,j+1)end doend dodo i=2,42do j=2,55cv(i,j)=cd(i,j)*cd(i-1,j)end doend do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP3:迭代循环求解!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!h_check=0do j=19,54h(42,j)=h0(j)*cosd(-g0(j))end doprint*,'迭代控制误差',epsprint*,'迭代次数'do while(maxval(abs(h-h_check))>eps)t=1h_check=hdo while(t<=240)do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)t=t+1do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)t=t+1end dok=k+1print*,kend do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !STEP 4:迭代收敛后输出h，u，v!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!t=1do while(t<=240)do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)h_3d(:,:,t)=hcall v_get(v,u,cu,cv,dx,dy,dt,d,h,f)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)t=t+1do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)h_3d(:,:,t)=hcall u_get(u,v,cu,cv,dx,dy,dt,d,h,f)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)t=t+1end do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP 5 : 用最小二乘法求解各个网格点的调和常数hh和g !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 最小二乘计算原理：《计算方法引论》P52! 调和分析原理：《海洋水文环境要素的分析方法和预报》P97! 求正规方程组的系数l1=0do t=1,240l1=l1+cosd(w*t*dt)**2end dol2=0do t=1,240l2=l2+sind(w*t*dt)*cosd(w*t*dt)end dol3=l2l4=0do t=1,240l4=l4+sind(w*t*dt)**2end dodo i=1,43do j=1,56do t=1,240temp1(t)=cosd(w*t*dt)*h_3d(i,j,t)temp2(t)=sind(w*t*dt)*h_3d(i,j,t)end doy1(i,j)=sum(temp1)y2(i,j)=sum(temp2)end doend do!求解调和常数do i=1,43do j=1,56if(cd(i,j)==0)cyclea(i,j)=(l4*y1(i,j)-l2*y2(i,j))/(l1*l4-l2*l3)b(i,j)=(l3*y1(i,j)-l1*y2(i,j))/(l2*l3-l1*l4)hh(i,j)=sqrt(a(i,j)**2+b(i,j)**2)if(hh(i,j)==0)cyclecos_g=a(i,j)/hh(i,j)sin_g=b(i,j)/hh(i,j)if(cos_g==1.and.sin_g==0)theng(i,j)=0else if(cos_g>0.and.sin_g>0)theng(i,j)=atand(b(i,j)/a(i,j))else if(cos_g==0.and.sin_g==1)theng(i,j)=90else if(cos_g<0.and.sin_g>0)theng(i,j)=atand(b(i,j)/a(i,j))+180else if(cos_g==-1.and.sin_g==0)theng(i,j)=180else if(cos_g<0.and.sin_g<0)theng(i,j)=atand(b(i,j)/a(i,j))+180else if(cos_g==0.and.sin_g==-1)theng(i,j)=270elseg(i,j)=atand(b(i,j)/a(i,j))+360end ifend doend doopen(4,file='H_M2.dat')!open(4,file='H_K1.dat')write(4,'(56f8.3)') ((hh(i,j),j=1,56),i=1,43) close(4)open(5,file='g_M2.dat')!open(5,file='g_K1.dat')write(5,'(56f8.3)') ((g(i,j),j=1,56),i=1,43) close(5)end program hw4!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !子程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!计算海水相对于静止海面的高度hsubroutine h_get(h,dx,dy,dt,u,v,d,cd)implicit noneinteger::i,jreal,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cdreal,dimension(43,56),intent(in)::d,u,vreal,dimension(43),intent(in)::dxreal,dimension(43,56),intent(inout)::hdo i=2,41do j=2,55if(cd(i,j)/=0)thenh(i,j)=h(i,j)&&-dt/dx(i)*((d(i,j)+h(i,j)+d(i,j+1)+h(i,j+1))/2*u(i,j)-&&(d(i,j)+h(i,j)+d(i,j-1)+h(i,j-1))/2*u(i,j-1))&&-dt/dy*((d(i,j)+h(i,j)+d(i-1,j)+h(i-1,j))/2*v(i,j)-&&(d(i,j)+h(i,j)+d(i+1,j)+h(i+1,j))/2*v(i+1,j))elsecycleend ifend doend doend subroutine h_get! 计算vsubroutine v_get(v,u,cu,cv,dx,dy,dt,d,h,f)implicit noneinteger::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,ureal,dimension(43,56),intent(inout)::vreal,dimension(43),intent(in)::dx,freal,dimension(43,56)::uu,ssdo i=2,42do j=2,55if(cu(i,j)/=0)thenif(i==42)thenuu(i,j)=(uu(i,j-1)+uu(i,j+1))/2elseuu(i,j)=(u(i,j)+u(i,j-1)+u(i-1,j)+u(i-1,j-1))/4end ifelsecycleend ifend doend dodo i=2,42do j=2,55ss(i,j)=sqrt(uu(i,j)**2+v(i,j)**2)if(cv(i,j)/=0)thenv(i,j)=&&1.0/(1.0/dt+0.5*kk*ss(i,j)/((d(i,j)+d(i-1,j))/2+(h(i,j)+h(i-1,j))/2))& &*(v(i,j)/dt&&-v(i,j)*0.5*kk*ss(i,j)/((d(i,j)+d(i-1,j))/2+(h(i,j)+h(i-1,j))/2)&&-0.5*uu(i,j)*(v(i,j+1)-v(i,j-1))/dx(i)&&-0.5*v(i,j)*(v(i-1,j)-v(i+1,j))/dy&&-f(i)*uu(i,j)&&-gg*(h(i-1,j)-h(i,j))/dy)elsecycleend ifend doend doend subroutine v_get!计算usubroutine u_get(u,v,cu,cv,dx,dy,dt,d,h,f)implicit noneinteger::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,vreal,dimension(43,56),intent(inout)::ureal,dimension(43),intent(in)::dx,freal,dimension(43,56)::vv,rrdo i=2,42do j=2,55if(cv(i,j)/=0)thenif(i==42)thenvv(i,j)=v(i,j)elsevv(i,j)=(v(i,j)+v(i,j+1)+v(i+1,j)+v(i+1,j+1))/4end ifelsecycleend ifend doend dodo i=2,42do j=2,55rr(i,j)=sqrt(u(i,j)**2+vv(i,j)**2)if(cu(i,j)/=0)thenu(i,j)=&&1.0/(1.0/dt+0.5*kk*rr(i,j)/((d(i,j)+d(i,j+1))/2+(h(i,j)+h(i,j+1))/2))& &*(u(i,j)/dt&&-u(i,j)*0.5*kk*rr(i,j)/((d(i,j)+d(i,j+1))/2+(h(i,j)+h(i,j+1))/2)&&-0.5*u(i,j)*(u(i,j+1)-u(i,j-1))/dx(i)&&-0.5*vv(i,j)*(u(i-1,j)-u(i+1,j))/dy&&+f(i)*vv(i,j)&&-gg*(h(i,j+1)-h(i,j))/dx(i))elsecycleend ifend doend doend subroutine u_get。