Fortran 读取程序和求解潮汐各分潮的调和常数 范例

!程序中数组的各元素依次为分潮MSf,Q1,O1,K1,P1,K2,N2,M2,S2,MK3,M4,MS6,M6program tiaohechangshuimplicit noneinteger,dimension(78):n=(/0,2,-2,0,0,0,1,-2,0,1,0,0,1,-1,0,0,0,0,&1,1,0,0,0,0,1,1,-2,0,0,0,2,2,0,0,0,0,2,-1,0,1,0,0,2,0,0,0,0,0,2,2,-2,0,0,0,3,&1,0,0,0,0,4,0,0,0,0,0,4,2,-2,0,0,0,6,0,0,0,0,0/) !变量n中从左开始每隔六个数为一个分潮的杜德森数（n1,n2,n3,n4,n5,n6real,dimension(6):o=(/14.49205212,0.54901653,0.04106864,&0.00464183,0.00220641,0.00000196/)integer:y=1996,Y1=int(1996-1901)/4),D=198,t=3,i,j,ot,N2,g,k,g1 !经计算观测数据的中间时刻为1996年7月16日11时,转化成格林时间为t-8real,dimension(13):OO,V0,f,u,b,b1,R,C,H1,gg !其中变量OO、V0、f、u、H1、gg分别为各分潮的角速率、幅角、交点因子、交点订正角、调和函数（H1和gg）real:T0,s0,h0,p0,N0,PP0,TT,s,h,p,NN,PP,N1,fO1,fK1,fM2,fK2,uO1,uK1,uM2,uK2real,dimension(0:13,0:13):AA,AA1real,dimension(0:13,0:14):AAAreal,dimension(13,13):BB,BB1real,dimension(13,14):BBBreal,dimension(0:13):a1,areal,dimension(0:13,1):F1,F11real,dimension(13,1):F2,F22integer,dimension(62):year,month,dayinteger,dimension(1:24):time1=(/00,01,02,03,04,05,06,07,08,09,&10,11,12,13,14,15,16,17,18,19,20,21,22,23/)integer,dimension(62,12):hr1integer,dimension(-371:371):hr !经计算N=(743-1)/2=371!从(XM_July1996.dat)中读数据并创建.txt格式的文件open(8,file=F:Fortran90XM_July1996.dat)read(8,*)do i=1,62 read(8,(11x,i4,i2,i2,x,12i5),year(i),month(i),day(i),hr1(i,:)end doclose(8)open(4,file=F:Fortran90XM_July1996lie.txt)do i=1,62do j=1,12if (mod(i,2)=0) then write(4,(i4,i1,i2,i2,i5),year(i),month(i),day(i),time1(12+j),hr1(i,j)else write(4,(i4,i1,i2,i2,i5),year(i),month(i),day(i),time1(j),hr1(i,j)end ifend doend doclose(4)open(5,file=F:Fortran90XM_July1996lie.txt)do j=-371,371read(5,(10x,i5)hr(j)end doclose(5)!各天文参量每日Greenwish零时的值T0=180.0s0=277.025+129.38481*(y-1900)+13.17640*(D+Y1)h0=280.190-0.23872*(y-1900)+0.98565*(D+Y1)p0=334.385+40.66249*(y-1900)+0.11140*(D+Y1)N0=259.157-19.32818*(y-1900)-0.05295*(D+Y1)PP0=281.221+0.01718*(y-1900)+0.000047*(D+Y1)!计算某时刻的角速率TT=o(1)*t+T0s=o(2)*t+s0h=o(3)*t+h0p=o(4)*t+p0NN=o(5)*t+N0PP=o(6)*t+PP0!求中间时刻各分潮的f,uN1=-NNfO1=1.0089+0.1871*cosd(N1)-0.0147*cosd(2*N1)+0.0014*cosd(3*N1)fK1=1.0060+0.1150*cosd(N1)-0.0088*cosd(2*N1)+0.0006*cosd(3*N1)fM2=1.0004-0.0373*cosd(N1)+0.0003*cosd(2*N1)fK2=1.0241+0.2863*cosd(N1)+0.0083*cosd(2*N1)-0.0015*cosd(3*N1)uO1=10.80*sind(N1)-1.34*sind(2*N1)+0.19*sind(3*N1)uK1=-8.86*sind(N1)+0.68*sind(2*N1)-0.07*sind(3*N1)uM2=-2.14*sind(N1)uK2=-17.74*sind(N1)+0.68*sind(2*N1)-0.04*sind(3*N1)f=(/fM2,fO1,fO1,fK1,1.0,fK2,fM2,fM2,1.0,fM2*fK1,fM2*2,fM2,fM2*3/)u=(/-uM2,uO1,uO1,uK1,0.0,uK2,uM2,uM2,0.0,uM2+uK1,2*uM2,uM2,3*fM2/)!求各分潮的角速率和幅角do i=1,13OO(i)=n(6*i-5)*o(1)+n(6*i-4)*o(2)+n(6*i-3)*o(3)+n(6*i-2)*o(4)+n(6*i-1)*o(5)+n(6*i)*o(6) V0(i)=n(6*i-5)*T0+n(6*i-4)*s0+n(6*i-3)*h0+n(6*i-2)*p0+n(6*i-1)*N0+n(6*i)*PP0 end do!求解方程组的系数矩阵和增广矩阵AA=0;AA(0,0)=743;ot=1;N2=743;F1=0;F2=0do i=-371,371F1(0,1)=F1(0,1)+hr(i)end dodo i=1,13do j=-371,371F1(i,1)=F1(i,1)+hr(j)*cosd(j*OO(i)*ot)F2(i,1)=F2(i,1)+hr(j)*sind(j*OO(i)*ot)end doend dodo j=1,13AA(j,0)=sind(N2*OO(j)*ot/2.0)/sind(OO(j)*ot/2.0)AA(0,j)=sind(N2*OO(j)*ot/2.0)/sind(OO(j)*ot/2.0)AA(j,j)=(N2+sind(N2*OO(j)*ot)/sind(OO(j)*ot)/2.0BB(j,j)=(N2-sind(N2*OO(j)*ot)/sind(OO(j)*ot)/2.0do i=1,13if (i/=j) thenAA(i,j)=(sind(N2*(OO(i)-OO(j)*ot/2.0)/sind(OO(i)-OO(j)*ot/2.0)+sind(N2*(OO(i)+OO(j)*ot/2.0)/sind(OO(i)+OO(j)*ot/2.0)/2.0BB(i,j)=(sind(N2*(OO(i)-OO(j)*ot/2.0)/sind(OO(i)-OO(j)*ot/2.0)-sind(N2*(OO(i)+OO(j)*ot/2.0)/sind(OO(i)+OO(j)*ot/2.0)/2.0end ifend doend doAAA(:,0:13)=AA(:,0:13)AAA(:,14)=F1(:,1)BBB(:,1:13)=BB(:,:)BBB(:,14)=F2(:,1)!用高斯消元法求解g=size(AAA,dim=1)-1;g1=size(BBB,dim=1)!求解各分潮a的值do j=0,g-1 do i=j+1,g AAA(i,:)=AAA(i,:)-AAA(j,:)*(AAA(i,j)/AAA(j,j) end doend doAA1(:,:)=AAA(:,0:g)F11(:,1)=AAA(:,g+1)a(g)=F11(g,1)/AA1(g,g)do i=g-1,0,-1a(i)=(F11(i,1)-sum(AA1(i,(i+1):g)*a(i+1):g)/AA1(i,i)end do!求解各分潮b的值do j=1,g1-1 do i=j+1,g1 BBB(i,:)=BBB(i,:)-BBB(j,:)*(BBB(i,j)/BBB(j,j) end doend doBB1(:,:)=BBB(:,1:g1)F22(:,1)=BBB(:,g1+1)b(g1)=F22(g1,1)/BB1(g1,g1)do i=g1-1,1,-1b(i)=(F