气象统计方法实习BD
实习一：气候场、距平场、均方差场 编程如下:
parameter(ii=37,jj=17,mon=12,year=4)
real var(ii,jj,mon,year),ave(ii,jj,mon),jp(ii,jj,mon,year) real s(ii,jj,mon) integer i,j,iy,m
open(5,file='d:\ex1\h500.dat')
open(6,file='d:\ex1\ave.grd',form='binary') open(7,file='d:\ex1\jp.grd',form='binary') open(8,file='d:\ex1\s.grd',form='binary')
open(12,file='d:\ex1\outall.grd',form='binary' open(9,file='d:\ex1\ave.txt') open(10,file='d:\ex1\jp.txt') open(11,file='d:\ex1\s.txt') ！读数据
DO iy=1,4 do m=1,12 !ccc read h500
read(5,1000)
read(5,2000) ((var(i,j,m,iy),i=1,ii),j=1,jj) enddo enddo ！计算气候场 do j=1,jj do i=1,ii do m=1,12
ave(i,j,m)=var(i,j,m,1)+var(i,j,m,2)+var(i,j,m,3)+var(i,j,m,4) ave(i,j,m)=ave(i,j,m)/4.0 enddo enddo enddo ！计算距平场 do iy=1,4 do m=1,12 do j=1,jj do i=1,ii
jp(i,j,m,iy)=var(i,j,m,iy)-ave(i,j,m) enddo enddo enddo enddo ！计算均方差场
do j=1,jj do i=1,ii do m=1,12
s(i,j,m)=jp(i,j,m,1)*jp(i,j,m,1)+jp(i,j,m,2)*jp(i,j,m,2)+jp(i,j /,m,3)*jp(i,j,m,3)+jp(i,j,m,4)*jp(i,j,m,4) s(i,j,m)=s(i,j,m)/4.0 s(i,j,m)=sqrt(s(i,j,m)) enddo enddo enddo
do iy=1,4 do m=1,12
write(6)((ave(i,j,m),i=1,ii),j=1,jj) write(7)((jp(i,j,m,iy),i=1,ii),j=1,jj) write(8)((s(i,j,m),i=1,ii),j=1,jj)
write(9,2000)((ave(i,j,m),i=1,ii),j=1,jj) write(10,2000)((jp(i,j,m,iy),i=1,ii),j=1,jj) write(11,2000)((s(i,j,m),i=1,ii),j=1,jj) write(12)((ave(i,j,m),i=1,ii),j=1,jj) write(12)((jp(i,j,m,iy),i=1,ii),j=1,jj) write(12)((s(i,j,m),i=1,ii),j=1,jj)
enddo enddo
1000 format(2i7) 2000 format(37f8.1) close(5) close(6) close(7) close(8) close(9) close(10) close(11) close(12) end
给ave配的ctl文件： dset ^d:\ex1\ave.grd undef -9.99E+33
title NCEP/NCAR REANALYSIS PROJECT xdef 37 linear 60.000 2.500 ydef 17 linear 0.000 2.500 zdef 1 levels 500
tdef 12 linear JAN1982 12mo vars 1
ave 1 99 H500 endvars
给ave配的gs文件： 'reinit'
'open d:\ex1\ave.ctl'
'enable print d:\ex1\ave.gmf' mon=1
while(mon
'draw title qihouchang of 'mon' ' 'print' 'c'
mon=mon+1 endwhile
'disable print' ;
气候场图：
一月份高度的气候场呈现南高北低的状态，陆地上的高度场比较稀疏，而在西太平洋上高度场比较密集。
八月份高度的气候场呈现东高西低的状态，在我国东北部以北以及印度东北部出现低压中心，而在赤道西太平洋地区出现高压中心。35°N以北高度分布很密集，而35°N以南比较稀疏。
给jp配的ctl文件： dset ^d:\ex1\jp.grd undef -9.99E+33
title NCEP/NCAR REANALYSIS PROJECT xdef 37 linear 60.000 2.500 ydef 17 linear 0.000 2.500 zdef 1 levels 500
tdef 48 linear JAN1982 1mo vars 1
jp 1 99 H500 endvars
给jp配的gs文件： 'reinit'
'open d:\ex1\jp.ctl'
'enable print d:\ex1\jp.gmf' year=1982
while(year
while(mon
'draw title jupingchang of 'year'.'mon'' 'print' 'c'
mon=mon+1
endwhile year=year+1 endwhile
'disable print' ;
距平场图：
1983年6月距平场在日本地区出现低压中心，在我国南部出现高压中心,在亚洲西北部也有高压中心。赤道至25°N间以及25°N-40°N，60°E-100°E间基本都是正距平，而在25°N-40°N，100°E-150°E间基本都是负距平。
1984年7月距平场在亚洲大陆西部、日本地区、赤道西太平洋地区形成低压中心，太平洋西北部形成高压中心。
给s配的ctl文件： dset ^d:\ex1\s.grd undef -9.99E+33
title NCEP/NCAR REANALYSIS PROJECT xdef 37 linear 60.000 2.500 ydef 17 linear 0.000 2.500 zdef 1 levels 500
tdef 12 linear JAN1982 12mo vars 1
s 1 99 H500 endvars
给s配的gs文件： 'reinit'
'open d:\ex1\s.ctl'
'enable print d:\ex1\s.gmf' mon=1
while(mon
'draw title junfangchachang of 'mon' ' 'print' 'c'
mon=mon+1 endwhile
'disable print' ;
均方差场图：
一月份高度的均方差场整体呈现南小北大的状态。说明低纬地区高度的波动幅度比较小，而中高纬地区高度的波动比较大。
八月份高度的均方差场在亚洲大陆西部有极大值，在30°N处包括赤道-30°N 、60°E-85°E这些区域高度的波动幅度比较小，30°N以南以北地区高度的波动幅度较大。
实习二：相关系数 Fortran程序如下：
program ex2
integer,parameter::n=20,p=10 integer i,j,t1,t2,t3 real
a(n),b(n),jpa(n),jpb(n),zxfc1(p),zxgxs1(p),zxfc2(p),zxgxs2(p),lhxfc(p),lhxgxs(p)
real::s1=0.0,s2=0.0,sum1=0.0,sum2=0.0,sum3=0.0,ave1,ave2,r,fc1,fc2 data
a/3.40,3.30,3.20,2.90,3.40,2.80,3.60,3.00,2.80,3.00,3.10,3.00,2.90,2.70,3.50,3.20,3.10,2.80,2.90,2.90/ data
b/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/
data jpa/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ data jpb/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ data zxfc1/0,0,0,0,0,0,0,0,0,0/ data zxgxs1/0,0,0,0,0,0,0,0,0,0/ data zxfc2/0,0,0,0,0,0,0,0,0,0/ data zxgxs2/0,0,0,0,0,0,0,0,0,0/ data lhxfc/0,0,0,0,0,0,0,0,0,0/ data lhxgxs/0,0,0,0,0,0,0,0,0,0/
!求均值 do i=1,n s1=s1+a(i) enddo ave1=s1/n do i=1,n s2=s2+b(i) enddo ave2=s2/n !求距平 do i=1,n
jpa(i)=a(i)-ave1 enddo do i=1,n
jpb(i)=b(i)-ave2 enddo
!求相关系数 do i=1,n
sum1=sum1+jpa(i)*jpa(i) sum2=sum2+jpb(i)*jpb(i) sum3=sum3+jpa(i)*jpb(i) enddo
r=sum3/(sqrt(sum1*sum2))
print*,'中国1970-1989年年平均和冬季平均气温的相关系数为r=',r !求方差 fc1=sum1/n fc2=sum2/n
!求自协方差 do j=1,p do i=1,n-j
zxfc1(j)=zxfc1(j)+jpa(i)*jpa(i+j) enddo
zxfc1(j)=zxfc1(j)/(n-j) enddo do j=1,p do i=1,n-j
zxfc2(j)=zxfc2(j)+jpb(i)*jpb(i+j) enddo
zxfc2(j)=zxfc2(j)/(n-j) enddo
!自相关系数 do i=1,p
zxgxs1(i)=zxfc1(i)/fc1 zxgxs2(i)=zxfc2(i)/fc2 enddo
!落后交叉协方差 do j=1,p do i=1,n-j
lhxfc(j)=lhxfc(j)+jpa(i)*jpb(i+j) enddo
lhxfc(j)=lhxfc(j)/(n-j) enddo
!落后相关系数 do i=1,p
lhxgxs(i)=lhxfc(i)/(sqrt(fc1*fc2)) enddo
print*,'年平均气温不同滞后时刻所对应的自相关系数为：' print *,((zxgxs1(i),','),i=1,P)
print *
print*,'冬季平均气温不同滞后时刻所对应的自相关系数为：' print *,((zxgxs2(i),','),i=1,P) print *
print*,'年平均气温和冬季平均气温之间不同滞后时刻所对应的落后交叉相关系数为：' print *,((lhxgxs(i),','),i=1,P) print * !最大绝对值 t1=1 t2=1 t3=1 do i=2,p
if(abs(zxgxs1(i))>abs(zxgxs1(t1)))t1=i if(abs(zxgxs2(i))>abs(zxgxs2(t2)))t2=i if(abs(lhxgxs(i))>abs(lhxgxs(t3)))t3=i enddo
print *,'中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度为：',t1,zxgxs1(t1)
print *,'中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度为：',t2,zxgxs2(t2)
print *,'中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度为：',t3,lhxgxs(t3) end
程序运行如下：
中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度是7，自相关系数是-0.3724，说明在这些年的数据中，滞后7年的序列与原序列的相关系数绝对值最大，呈反相关。
中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度是4，自相关系数是-0.3678，说明在这些年的数据中，滞后4年的序列与原序列的相关系数绝对值最大，呈反相关。
中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度是3，自相关系数是-0.4066，说明在这些年的数据中，滞后3年的冬季平均气温序列与原平均气温序列的相关系数绝对值最大，呈反相关。
实习三：一元线性回归 Fortran程序如下： program ex3
integer,parameter::n=17 integer i,j
real::suma1=0,suma2=0,avea1,avea2,s1=0,s2=0,b0,b,x=14.5,y integer a1(n) real a2(n),a3(n)
data a1/0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/
data a2/30.0,29.1,28.4,28.1,28.0,27.7,27.5,27.2,27.0,26.8,26.5,26.3,26.1,25.7,25.3,24.8,24.0/ data a3/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ open(5,file='d:\3\ys.txt')
open(15,file='d:\3\ys.grd',form='binary') open(6,file='d:\3\hs.txt')
open(16,file='d:\3\hs.grd',form='binary')
!求平均值 do i=1,n
suma1=suma1+a1(i) suma2=suma2+a2(i) enddo
avea1=suma1/n avea2=suma2/n !求b,b0值 do i=1,n
s1=s1+a1(i)*a2(i) s2=s2+a1(i)*a1(i) enddo
b=(s1-suma1*suma2/n)/(s2-suma1*suma1/n) b0=avea2-avea1*b y=b0+b*x
print*,'y=',b0,'+',b,'x' print* print*,y
!求回归数据 do i=1,n
a3(i)=b0+b*a1(i) enddo
write(5,1000)(a1(i),i=1,n) write(5,2000)(a2(i),i=1,n) write(15)(a1(i),i=1,n) write(15)(a2(i),i=1,n)
write(6,1000)(a1(i),i=1,n) write(6,2000)(a3(i),i=1,n) write(16)(a1(i),i=1,n) write(16)(a3(i),i=1,n)
1000 format(i7) 2000 format(f8.1) end
程序运行：
所给数据做出的线性回归曲线的斜率是-0.3012，截距是29.3804，说明所给数据y随着x递减。
实习四：滑动平均 Fortran程序如下： PROGRAM MA integer i integer,parameter::n=85,ih=11,nyear=1922 integer::yr(n)=0 real X(n),X1(n) real::s(75)=0
C ********************************************** C * N: SAMPLE SIZE OF THE TIME SERIES * C * IH: MOVING LENGTH * C * NYEAR: FIRST YEAR OF THE SERIES * C * X(N): OROGINAL TIME SERIES * C * X1(N-IH+1): MOVED SERIES * C ********************************************** OPEN(2,FILE='d:\4\MA.DAT') open(12,file='d:\4\ma.grd',form='binary') OPEN(3,FILE='d:\4\hd.DAT') open(13,file='d:\4\hd.grd',form='binary') !年份 do i=1,n yr(i)=1922-1+i enddo !读入数据 READ(2,*)(X(I),I=1,N) !计算滑动平均 do i=1,11 s(1)=s(1)+x(i) enddo x1(6)=s(1)/ih do i=7,80 s(i-5)=s(i-6)+x(i+ih-6)-x(i-6) x1(i)=s(i-5)/ih enddo !写入数据
write(3,'(1x,i5,1X,f5.1,1X,f5.1)')((yr(i),x(i),x1(i)),i=1,n) write(12)((x(i)),i=1,n) write(13)((x1(i)),i=1,n) close(2) close(3) close(12) close(13) end
原数据和滑动后数据的图形：
由图可知，所给数据在1922-1955年之间呈波动下降趋势，在1955-1968年呈波动上升趋势，上升幅度较大，而1968-2006年之间大致在同一水平上波动，没有升降趋势。
实习五：eof Fortran程序： C$large
PROGRAM EOFPW
C THIS PROGRAM USES EOF FOR ANALYSING TIME SERIES C OF METEOROLOGICAL FIELD C M:LENTH OF TIME SERIES C N:NUMBER OF GRID-POINTS
C KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDERDLIZED DEPATURE C KV:NUMBER OF EIGENVALUES WILL BE OUTPUT
C KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT C MNH=MIN(M,N)
C Evf=EIGENVACTORS, tcF=TIME COEFFICIENTS FOR EGVT. C ER(KV,1)=LAMDA,LAMDA EIGENVALUE C ER(KV,2)=ACCUMULATE LAMDA
C ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTO c EIGENVACTOR.
C ER(KV,4)=ACCUMULATE ER(KV,3) C
PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2) PARAMETER(ff=-999.0,nx=18,ny=12) C
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)
CCCCCCCCCCCCCCCCINPUT DATA
open(11,file='d:\5\sstpx.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)
CCCCCCCCCCCCCCCCINPUT DATA CCCCCCCCCCCCCCCCCCC ccccccccccccccccccccccccccccccccccccc 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:\5\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:\5\tcf.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)
c CALL OUTER(KV,ER,mnh) open(21,file='d:\5\dats.dat')
write(21,106)(er(iiii,3),iiii=1,kv) close(21)
c CALL OUTVT(KVT,N,M,MNH,S,F,Evf,tcF) 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)
C THIS SUBROUTINE PROVIDES INITIAL F BY KS and kv. 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
if(ks.eq.0)then goto 5 endif do i=1,n1 df(i)=0.0 enddo
cccccccccccccccccc 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
IF(KS.EQ.0) THEN 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) C THIS SUBROUTINE FORMS A BY F 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)
C THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A 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(U.EQ.0.0) G=1.
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 IF(S3.GT.S2) GOTO 50 RETURN END
SUBROUTINE ARRANG(KV,MNH,A,ER,S)
C THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES C FROM MAX TO MIN
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)
C THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS (M.GE.N,SAVED IN S;
C M.LT.N,SAVED IN F) AND ITS TIME COEFFICENTS SERIES (M.GE.N, C SAVED IN F; M.LT.N,SAVED IN S)
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
cccccccccccccccccccccccccccccccccccccccccc c t=0.0
c do 365 i=1,mnh c365 t=t+s(i,1)*s(i,2) c write(*,*)t
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF(N.LE.M) THEN 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
ccccccccccccccccccccccccccccccccccccccccccccccccccccc 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) c this subroutine sent undefine value ff to evf 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)
C THIS SUBROUTINE PRINTS ARRAY ER
C ER(KV,1) FOR SEQUENCE OF EIGENVALUE FROM BIG TO SMALL C ER(KV,2) FOR EIGENVALUE FROM BIG TO SMALL
C ER(KV,3) FOR SMALL LO=(LAMDA/TOTAL VARIANCE) C ER(KV,4) FOR BIG LO=SUM OF SMALL LO) 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) C THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS C AND ITS TIME-COEFFICENT SERIES
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
IF(M.GE.N) THEN
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 IF(M.GE.N) THEN
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
程序运行：
1949、1950、1951、1953、1955、1956、1960-1963、1968、1971、1972、1974、1976、1979、1981、1985、1989这些年份，赤道太平洋中部至东部的海温都是较正常水平偏低的，而在1952、1954、1958、1964、1966、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这些年份，赤道太平洋地区的海温是偏高的。