2014年气象统计实习报告
实 习 报 告 书
课程名称：气象统计方法课程实践
姓名： 学号：
班级： 级气科 班
\\\\*实习一 求500hPa高度场气候场、距平场和均方差场
实习时间：第9周周三1、2节
1. 资料介绍
有一500hPa高度场资料，文件名h500.dat，范围：60~150E，0~40N.
2．要求
编fortran程序，求500hPa高度场的 （1）气候场； （2）距平场； （3）均方差场。
并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。
注：h500.For给出了如何用fortran读取ASCII码资料h500.dat. FORTRAN
program sx1 implicit none integer nx,ny,mo,yr
parameter(nx=37,ny=17,mo=12,yr=4) real var(nx,ny,mo,yr) real
at(nx,ny,mo),xd(nx,ny,mo,yr),sx(nx,ny
,mo)
integer i,j,m,t,it,iy,irec
open(5,file='d:\\study\\form\\shixione\\h500.dat') do iy=1,4 do m=1,12 read(5,1000)
read(5,3000)((var(i,j,m,iy),i=1,nx),j=1,ny) enddo enddo close(5) !计算气候场at do t=1,12 do j=1,ny do i=1,nx at(i,j,t)=0 do it=1,4
at(i,j,t)=at(i,j,t)+var(i,j,t,it) enddo
at(i,j,t)=at(i,j,t)/4 enddo enddo enddo !求距平场xd do t=1,12 do j=1,ny do i=1,nx xd(i,j,t,1)=0 do it=1,4
xd(i,j,t,it)=var(i,j,t,it)-at(i,j,t) enddo enddo enddo enddo
!求均方差场sx do t=1,12 do j=1,ny do i=1,nx sx(i,j,t)=0 do it=1,4
sx(i,j,t)=sx(i,j,t)+(var(i,j,t,it)-at(i,j,t))**2 enddo
sx(i,j,t)=sqrt(sx(i,j,t)/4) enddo enddo enddo !写入气候场
open(10,file='d:\\study\\form\\shixione\\at.grd',form='unformatted',access='direct',recl=nx*ny) irec=0 do t=1,12 irec=irec+1
write(10,rec=irec)((at(i,j,t),i=1,nx),j=1,ny) enddo close(10) !写入距平场
open(11,file='d:\\study\\form\\shixione\\xd.grd',form='unformatted',access='di
rect',recl=nx*ny) irec=0 do it=1,4 do t=1,12 irec=irec+1
write(11,rec=irec)((xd(i,j,t,it),i=1,nx),j=1,ny) enddo enddo close(11) ！写入均方差场
open(12,file='d:\\study\\form\\shixione\\sx.grd',form='unformatted',access='direct',recl=nx*ny) irec=0 do t=1,12 irec=irec+1
write(12,rec=irec)((sx(i,j,t),i=1,nx),j=1,ny) enddo close(12) 1000 format(2i7)
2000 format(37F6.2) 3000 format(37f8.1) 4000 format(37f7.2) end program sx1
运行结果：
Grads文件 气候场
'reinit' 'enable print
d:\\study\\form\\shixione\\at.gmf' 'open d:\\study\\form\\shixione\\at.ctl' 'set grid off' 'set grads off' 'set lat 0 40'
'set lon 60 150' 'set lev 500' mon=1
while(mon<=12) 'set t 'mon'' 'd h'
'draw title 1982year'mon'month' 'print' 'c' mon=mon+1 endwhile 'disable print' ;
距平场
'reinit' 'enable print
d:\\study\\form\\shixione\\sx.gmf' 'open d:\\study\\form\\shixione\\sx.ctl' 'set grid off' 'set grads off' 'set lat 0 40' 'set lon 60 150' 'set lev 500' year=1982
while(year<=1984) mon=1
while(mon<=12) 'set t 'mon'' 'd h'
'draw title 500hPa 'year'year'mon'month anomaly' 'print' 'c' mon=mon+1 endwhile year=year+1 endwhile
'disable print';
均方差
'reinit' 'enable print
d:\\study\\form\\shixione\\sx.gmf' 'open d:\\study\\form\\shixione\\sx.ctl' 'set grid off' 'set grads off' 'set lat 0 40' 'set lon 60 150' 'set lev 500' mon=1
while(mon<=12) 'set t 'mon'' 'd h'
'draw title 500hPa 1982year'mon'month Mean-square Deviation' 'print' 'c' mon=mon+1 endwhile 'disable print' ;
运行结果：
上面图中只展示了部分而未全部添加
*实习二 计算给定数据资料的简单相关系数和自相关系数 实习时间：第10周周三1、2节
根据下表中年平均气温和冬季平均气温的等级数据进行下列计算： 1）计算两个气温之间的简单相关系数。
2）分别找出两个气温数据自相关系数绝对值最大的滞后时间长度。（滞后长度τ最大取10）
要求：实习报告中附出简单相关系数或自相关系数程序。 答案：r=0.47
年平均气温在滞后长度j=7、冬季序列在j=4最大。 FORTRAN
计算相关系数r PROGRAM EXAM IMPLICIT NONE
INTEGER,PARAMETER::N=20 INTEGER i,j,k,ty,tw,tyw
REAL::avr_y=0,avr_w=0,sy=0,sw=0,rxy=0,max_y=0,max_w=0,max_yw=0 REAL y(N),w(N) DATA
y/3.4,3.3,3.2,2.9,3.4,2.8,3.6,3.0,2.8,3.0,3.1,3.0,2.9,
2.7,3.5,3.2,3.1,2.8,2.9,2.9/ DATA
w/3.24,3.14,3.26,2.38,3.32,2.71,2.84,3.94,2.75,1.83,2.80,2.81,2.63,3.20,3.60,3.40,3.07,1.87,2.63,2.47/ REAL
syy(N),sww(N),r(N),rty(N),rtw(N),rtyw(N),rxy_ty(N),rxy_tw(N),rxy_tyw(N)
!求两数组平均值 DO i=1,N avr_y=avr_y+y(i) avr_w=avr_w+w(i) END DO avr_y=avr_y/N avr_w=avr_w/N !简单相关系数 DO j=1,N
syy(j)=(y(j)-avr_y)**2 sy=sy+syy(j)
sww(j)=(w(j)-avr_w)**2 sw=sw+sww(j) END DO
sy=sqrt(sy/N) sw=sqrt(sw/N) DO j=1,N
r(j)=((y(j)-avr_y)/sy)*((w(j)-avr_w)/sw) rxy=rxy+r(j) END DO rxy=rxy/N
PRINT \"(/'1970-1989年全年平均气温与冬季平均气温的简单相关系数rxy=',f5.2)\k=0 !自相关系数 DO ty=1,N/2 DO i=1,N-ty
rty(i)=((y(i)-avr_y)/sy)*((y(i+ty)-avr_y)/sy) rxy_ty(ty)=rxy_ty(ty)+rty(i) END DO
rxy_ty(ty)=rxy_ty(ty)/(N-ty) rxy_ty(ty)=ABS(rxy_ty(ty)) IF(rxy_ty(ty)>max_y) THEN max_y=rxy_ty(ty) k=ty END IF
END DO
PRINT \"('全年平均气温绝对值最大自相关系数rxy_ty=',f7.4,/,'k=0 DO tw=1,N/2 DO i=1,N-tw
rtw(i)=((w(i)-avr_w)/sw)*((w(i+tw)-avr_w)/sw) rxy_tw(tw)=rxy_tw(tw)+rtw(i) END DO
rxy_tw(tw)=rxy_tw(tw)/(N-tw) rxy_tw(tw)=ABS(rxy_tw(tw)) IF(rxy_tw(tw)>max_w) THEN max_w=rxy_tw(tw) k=tw END IF END DO
PRINT \"('冬季平均气温绝对值最大自相关系数rxy_tw=',f7.4,/,'k=0
滞
后
时
间
长
度
k=',I2)\
滞
后
时
间
长
度
k=',I2)\
END
运行成果：
*实习四 求给定数据的一元线性回归方程 实习时间：第12周周三1、2节
利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。
年份 1951 1952 1953 1954 气温T 0.9 1.2 2.2 2.4 环流指标 32 25 20 26
1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970
-0.5 2.5 -1.1 0 6.2 2.7 3.2 -1.1 2.5 1.2 1.8 0.6 2.4 2.5 1.2 -0.8 27 24 28 24 15 16 24 30 22 30 24 33 26 20 32 35 答案：
ˆ 7.5-0.23x yF=20.18>Fα=4.41，回归方程显著
program sx4 implicit none integer i
real:: x(20),y(20),
f1=0,f2=0,b,avx,avy,b0,sx=0,sy=0,sxy=0,f,rxy
data(x(i),i=1,20)/32,25,20,26,27,24,28,24,15,16,24,30,22,30,24,33,26,20,32,35/ data(y(i),i=1,20)/0.9,1.2,2.2,2.4,-0.5,2.
5,-1.1,0,6.2,2.7,3.2,-1.1,2.5,1.2,1.8,0.6,2.4,2.5,1.2,-0.8/ do i=1,20 f1=f1+x(i)*y(i) end do
avx=sum(x)/20 !求均值¦ avy=sum(y)/20 do i=1,20 f2=f2+x(i)**2 end do
b=(20*f1-sum(x)*sum(y))/(20*f2-sum(x)**2) !求¨®b b0=avy-b*avx print *,'b=',b print *,'b0=',b0 print*,'y=',b0,b,'x' sx=sx+(x(i)-avx)**2 sy=sy+(y(i)-avy)**2
sxy=sxy+(x(i)-avx)*(y(i)-avy) end do sx=sqrt(sx/20)
do i=1,20
sy=sqrt(sy/20) sxy=sxy/20 rxy=sx*b/sy
f=rxy**2*18/(1-rxy**2)
if(f>4.41) print *,'F=',f,'>Fα=4.41，回归方程显著
if(f<=4.41) print *,'F=',f,' *实习七 计算给定数据的11年滑动平均和累积距平 实习时间：第15周周三1、2节 利用数据ma.dat，编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。 要求：实习报告中附出程序，并给出原数据和滑动后数据的图形（1 张图）和累积距平数据图形（1张图） program sx7 implicit none integer n,ih,nyear parameter(n=85,ih=11,nyear=1922) real::x(n),x1(n-ih+1)=0,x2(n)=0,averx=0 integer i,j open(2,file='d:\\study\\form\\shiyanseven\\ma.dat') open(3,file='d:\\study\\form\\shiyanseven\\maresult.dat') open(4,file='d:\\study\\form\\shiyanseven\\accresult.dat') read(2,*)(x(i),i=1,n) do j=1,n-ih+1 !moving average do i=1,ih x1(j)=x1(j)+x(i+j-1) end do x1(j)=x1(j)/ih end do do i=1,n-ih+1 write(3,*) x1(i) end do do i=1,n !average of x averx=averx+x(i) end do averx=averx/n do i=1,n !accelarate do j=1,i x2(i)=x2(i)+(x(j)-averx) end do end do do i=1,n write(4,*) x2(i) end do End 图形显示： *实习八 对给定的海温数据进行EOF分析 实习时间：第16周周三1、2节 给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF分析，资料的时空范围可以根据sstpx.ctl获知。 数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。 Zhunsst.for给出了如何读取资料， Ssteof.for为对距平或者标准化距平处理后的资料进行EOF分析。 要求：实习报告中给出第一特征向量及其时间系数，并分析其时空特征。 program main parameter(mo=12,my=43,nx=18,ny=12,mt=516) dimension avf(mo,nx,ny),df(mo,nx,ny) dimension sst(nx,ny,mt),sst3(nx,ny,mt) open(1,file='d:\\study\\form\\shixieight\\sstpx.grd',form='unformatted',access='direct',recl=nx*ny) irec=1 do it=1,mt read(1,rec=irec)((sst(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do ! average do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 avf(k,i,j)=avf(k,i,j)+sst(i,j,it) end do avf(k,i,j)=avf(k,i,j)/my end do end do end do ! variance do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j))**2 end do df(k,i,j)=sqrt(df(k,i,j)/my) end do end do end do ! standardizing do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 if(sst(i,j,it)==-999.0)then sst3(i,j,it)=-999.0 else sst3(i,j,it)=(sst(i,j,it)-avf(k,i,j))/df(k,i,j) end if end do end do end do end do ! output file open(2,file='d:\\study\\form\\shixieight\\standard.grd',form='unformatted',access='direct',recl=nx*ny) irec=1 do it=1,mt write(2,rec=irec)((sst3(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do close(2) close(1) end 运行程序后： 得到zhunsst.grd文件，即图中standard 下面为老师给出程序，为运行方便将部分简略。 PROGRAM EOFPW PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2) PARAMETER(ff=-999.0,nx=18,ny=12) DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(mnh,4),DF(N),V(MNH),AVF(N),evf(N,KVT),tCF(M,KVT),data(Nx,ny),nf(N) open(11,file='d:\\study\\form\\shixieight\\zhunsst.grd',form='unformatted',access='direct',recl=nx*ny) do 132 it=1,m read(11,rec=it)((data(i,j),i=1,nx),j=1,ny ) do 132 jj=1,ny do 132 ii=1,nx kkkk=nx*(jj-1)+ii f(kkkk,it)=data(ii,jj) 132 continue close(11) CALL Test1(n,m,ff,f,nf) write(*,*)'ok2' CALL TRANSF(N,M,F,nf,AVF,DF,KS) write(*,*)'ok3' CALL FORMA(N,M,MNH,F,A) write(*,*)'ok4' CALL JCB(MNH,A,S,0.00001) write(*,*)'ok5' CALL ARRANG(KV,MNH,A,ER,S) write(*,*)'ok6' CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER) write(*,*)'ok7' call test3(N,ff,nf,evf,kvt) write(*,*)'ok8' open(21,file='d:\\study\\form\\shixieight\\evf.grd',form='unformatted',access='direct',recl=nx*ny) irec=0 do 668 kk=1,kvt irec=irec+1 668 write(21,rec=irec)((evf(nx*(j-1)+i,kk),i=1,nx),j=1,ny) close(21) open(21,file='d:\\study\\form\\shixieight\cf.grd',form='unformatted',access='direct',recl=kvt) irec=0 do 345 it=1,m irec=irec+1 345 write(21,rec=irec) (tcf(it,iik),iik=1,kvt) close(21) 106 format(10f8.4) open(21,file='d:\\study\\form\\shixieight\\dats.dat') write(21,106)(er(iiii,3),iiii=1,kv) close(21) STOP END SUBROUTINE Test1(n1,m,ff,f,nf) DIMENSION F(N1,M) DIMENSION nF(N1) do i=1,n1 nf(i)=0.0 enddo do i=1,n1 do k=1,m if(f(i,k).eq.ff)then f(i,k)=0.0 nf(i)=nf(i)+1 endif enddo enddo return end SUBROUTINE TRANSF(n1,m,f,nf,avf,df,ks) DIMENSION F(N1,M),AVF(N1) DIMENSION DF(N1) DIMENSION nF(N1) if(ks.eq.-1)then goto 30 endif do i=1,n1 avf(i)=0.0 enddo ifthen goto 5 endif do i=1,n1 df(i)=0.0 enddo 5 continue DO 141 I=1,N1 if(nf(i).ne.0) goto 141 do 12 j=1,m 12 AVF(I)=AVF(I)+F(I,J) AVF(I)=AVF(I)/M DO 14 J=1,M 14 F(I,J)=F(I,J)-AVF(I) 141 CONTINUE IFTHEN RETURN ELSE DO 241 I=1,N1 if(nf(i).ne.0) goto 241 DO 22 J=1,M 22 DF(I)=DF(I)+F(I,J)*F(I,J) DF(I)=SQRT(DF(I)/M) DO 24 J=1,M 24 F(I,J)=F(I,J)/DF(I) 241 CONTINUE ENDIF 30 CONTINUE RETURN END SUBROUTINE FORMA(N,M,MNH,F,A) DIMENSION F(N,M),A(MNH,MNH) IF(M-N) 40,50,50 40 DO 44 I=1,MNH DO 44 J=I,MNH A(I,J)=0.0 DO 42 IS=1,N 42 A(I,J)=A(I,J)+F(IS,I)*F(IS,J) A(J,I)=A(I,J) 44 CONTINUE RETURN 50 DO 54 I=1,MNH DO 54 J=I,MNH A(I,J)=0.0 DO 52 JS=1,M 52 A(I,J)=A(I,J)+F(I,JS)*F(J,JS) A(J,I)=A(I,J) 54 CONTINUE RETURN END SUBROUTINE JCB(N,A,S,EPS) DIMENSION A(N,N),S(N,N) DO 30 I=1,N DO 30 J=1,I IF(I-J) 20,10,20 10 S(I,J)=1. GO TO 30 20 S(I,J)=0. S(J,I)=0. 30 CONTINUE G=0. DO 40 I=2,N I1=I-1 DO 40 J=1,I1 40 G=G+2.*A(I,J)*A(I,J) S1=SQRT(G) S2=EPS/FLOAT(N)*S1 S3=S1 L=0 50 S3=S3/FLOAT(N) 60 DO 130 IQ=2,N IQ1=IQ-1 DO 130 IP=1,IQ1 IF(ABS(A(IP,IQ)).LT.S3) GOTO 130 L=1 V1=A(IP,IP) V2=A(IP,IQ) V3=A(IQ,IQ) U=0.5*(V1-V3) IF IF(ABS(U).GE.1E-10) G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U) ST=G/SQRT(2.*(1.+SQRT(1.-G*G))) CT=SQRT(1.-ST*ST) DO 110 I=1,N G=A(I,IP)*CT-A(I,IQ)*ST A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT A(I,IP)=G G=S(I,IP)*CT-S(I,IQ)*ST S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT 110 S(I,IP)=G DO 120 I=1,N A(IP,I)=A(I,IP) 120 A(IQ,I)=A(I,IQ) G=2.*V2*ST*CT A(IP,IP)=V1*CT*CT+V3*ST*ST-G A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST) A(IQ,IP)=A(IP,IQ) 130 CONTINUE IF(L-1) 150,140,150 140 L=0 GO TO 60 150 IFGOTO 50 RETURN END SUBROUTINE ARRANG(KV,MNH,A,ER,S) DIMENSION A(MNH,MNH),ER(mnh,4),S(MNH,MNH) TR=0.0 DO 200 I=1,MNH TR=TR+A(I,I) 200 ER(I,1)=A(I,I) MNH1=MNH-1 DO 210 K1=MNH1,1,-1 DO 210 K2=K1,MNH1 IF(ER(K2,1).LT.ER(K2+1,1)) THEN C=ER(K2+1,1) ER(K2+1,1)=ER(K2,1) ER(K2,1)=C DO 205 I=1,MNH C=S(I,K2+1) S(I,K2+1)=S(I,K2) S(I,K2)=C 205 CONTINUE ENDIF 210 CONTINUE ER(1,2)=ER(1,1) DO 220 I=2,KV ER(I,2)=ER(I-1,2)+ER(I,1) 220 CONTINUE DO 230 I=1,KV ER(I,3)=ER(I,1)/TR ER(I,4)=ER(I,2)/TR 230 CONTINUE WRITE(6,250) TR 250 FORMAT(/5X,'TOTAL SQUARE ERROR=',F20.5) RETURN END SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER) DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(mnh,4),evf(n,kvt),tcf(m,kvt) DO 360 J=1,KVT C=0. DO 350 I=1,MNH 350 C=C+S(I,J)*S(I,J) C=SQRT(C) DO 160 I=1,MNH s(I,J)=S(I,J)/C 160 evf(I,J)=S(I,J)/C 360 CONTINUE IFTHEN DO 390 J=1,M DO 370 I=1,N V(I)=F(I,J) F(I,J)=0. 370 CONTINUE do 371 is=1,kvt tcf(j,is)=0. 371 continue DO 380 IS=1,KVT DO 380 I=1,N f(is,j)=F(IS,J)+V(I)*S(I,IS) 380 tcf(j,is)=tcf(J,is)+V(I)*S(I,IS) 390 CONTINUE ELSE DO 410 I=1,N DO 400 J=1,M V(J)=F(I,J) F(I,J)=0. 400 CONTINUE DO 410 JS=1,KVT DO 410 J=1,M f(I,JS)=F(I,JS)+V(J)*S(J,JS) 410 CONTINUE DO 430 JS=1,KVT DO 420 J=1,M tcf(J,JS)=S(J,JS)*SQRT(ER(JS,1)) 420 CONTINUE DO 430 I=1,N evf(I,JS)=F(I,JS)/SQRT(ER(JS,1)) 430 CONTINUE t=0.0 do 3650 i=1,m 3650 t=t+tcf(i,1)*tcf(i,2) write(*,*)t t=0.0 do 3651 i=1,n 3651 t=t+evf(i,1)*evf(i,2) write(*,*)t ENDIF RETURN END SUBROUTINE test3(N1,ff,nf,evf,kvt) dimension nf(n1),evf(n1,kvt) do i=1,n1 if(nf(i).ne.0)then do k=1,kvt evf(i,k)=ff enddo endif enddo return end SUBROUTINE OUTER(KV,ER,mnh) DIMENSION ER(mnh,4) WRITE(16,510) 510 FORMAT(/30X,'EIGENVALUE AND ANALYSIS ERROR',/) WRITE(16,520) 520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH) WRITE(16,530) (IS,(ER(IS,J),J=1,4),IS=1,KV) 530 FORMAT(1X,I10,4F15.5) WRITE(16,540) 540 FORMAT(//) RETURN END SUBROUTINE OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF) DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT) WRITE(16,560) 560 FORMAT(30X,'STANDARD EIGENVECTORS',/) WRITE(16,570) (IS,IS=1,KVT) 570 FORMAT(3X,80I7) DO 550 I=1,N IFTHEN WRITE(16,580) I,(S(I,JS),JS=1,KVT) 580 FORMAT(1X,I3,80F7.3,/) DO 11 JS=1,KVT EGVT(I,JS)=S(I,JS) 11 CONTINUE ELSE WRITE(16,590) I,(F(I,JS),JS=1,KVT) 590 FORMAT(1X,I3,80F7.3) DO 12 JS=1,KVT EGVT(I,JS)=F(I,JS) 12 CONTINUE ENDIF 550 CONTINUE WRITE(16,720) 720 FORMAT(//) WRITE(16,610) 610 FORMAT(30X,'TIME-COEFFICENT SERIES OF S. E.'/) WRITE(16,620) (IS,IS=1,KVT) 620 FORMAT(3X,80I7) DO 600 J=1,M IFTHEN WRITE(16,630) J,(F(IS,J),IS=1,KVT) 630 FORMAT(1X,I3,80F7.1) DO 13 IS=1,KVT ECOF(J,IS)=F(IS,J) 13 CONTINUE ELSE WRITE(16,640) J,(S(J,IS),IS=1,KVT) 640 FORMAT(1X,I3,80F7.1) DO 14 IS=1,KVT ECOF(J,IS)=S(J,IS) 14 CONTINUE ENDIF 600 CONTINUE RETURN END 运行结果 编写CTL文件 Evf.ctl dset d:\\study\\form\\shixieight\\evf.grd title Coads SSTA E undef -999.0 xdef 18 linear 120 10 ydef 12 linear -27.5 5 zdef 1 linear 1000 1 tdef 2 linear 1jan1948 1month vars 1 S 0 99 Coads SST anomaly interperated using Endvars Tcf.ctl dset d:\\study\\form\\shixieight\cf.grd title Coads SSTA T undef -999.0 xdef 1 linear 120 10 ydef 1 linear -27.5 5 zdef 1 linear 1000 1 tdef 516 linear 1jan1948 1month vars 2 a 0 99 time coefficient 1 b 0 99 time coefficient 2 endvars 根据TCL文件编写相应的GS文件 Evf.gs 'reinit' 'enable print d:\\study\\form\\shixieight\\evf.gmf' 'open d:\\study\\form\\shixieight\\evf.ctl' day=1 while(day<=2) 'set t 'day'' 'd s' 'draw title 'day'evf' 'print' 'c' day=day+1 endwhile 'disable print' ; Tcfa.gs 'reinit' 'enable print d:\\study\\form\\shixieight\cf.gmf' 'open d:\\study\\form\\shixieight\cf.ctl' 'set t 1 516' 'd a' 'draw title tcfa' 'print' 'disable print' ; Tcfb.gs 'reinit' 'enable print d:\\study\\form\\shixieight\cfb.gmf' 'open d:\\study\\form\\shixieight\cf.ctl' 'set t 1 516' 'd b' 'draw title tcfb' 'print' 'disable print' ; 效果图 空间场 时间系数 tcfa 由图中可以看出： 1949-1951、1953、1955-1956、1960-1963、1968、1971-1972、1974、1976、1979、1981、1985、1989这些年份，赤道太平洋中部至东部的海温都是较正常水平偏低，而在1952、1954、1958、1964、1969、1973、1977-1978、1980、1982-1983、1987、1991这些年份，赤道太平洋中部至东部的海温都是偏高的，尤其在1982-1983年，据记载，这段时间出现了厄尔尼诺现象。 空间场 时间系数 从图中可以看出： 1948-1954、1955-1958、1965、1969、1972、1974-1977、1979、1983-1985、1987、1989 这些年份，赤道太平洋地区的海温偏低，其余年份尤其是1954、1959、1964、1967、1968、1973、1978、1980、1982、1986、1988、1991这些年份，赤道太平洋地区的海温是偏高。 因篇幅问题不能全部显示，请点此查看更多更全内容