气象统计方法实习BD.docx
《气象统计方法实习BD.docx》由会员分享，可在线阅读，更多相关《气象统计方法实习BD.docx（37页珍藏版）》请在冰豆网上搜索。
气象统计方法实习BD
实习一：
气候场、距平场、均方差场
编程如下:
parameter（ii=37,jj=17,mon=12,year=4）
realvar（ii,jj,mon,year）,ave（ii,jj,mon）,jp（ii,jj,mon,year）
reals（ii,jj,mon）
integeri,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'）
！
读数据
DOiy=1,4
dom=1,12
!
cccreadh500
read（5,1000）
read（5,2000）（（var（i,j,m,iy）,i=1,ii）,j=1,jj）
enddo
enddo
！
计算气候场
doj=1,jj
doi=1,ii
dom=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
！
计算距平场
doiy=1,4
dom=1,12
doj=1,jj
doi=1,ii
jp（i,j,m,iy）=var（i,j,m,iy）-ave（i,j,m）
enddo
enddo
enddo
enddo
！
计算均方差场
doj=1,jj
doi=1,ii
dom=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
doiy=1,4
dom=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
1000format（2i7）
2000format（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
titleNCEP/NCARREANALYSISPROJECT
xdef37linear60.0002.500
ydef17linear0.0002.500
zdef1levels500
tdef12linearJAN198212mo
vars1
ave199H500
endvars
给ave配的gs文件：
'reinit'
'opend:
\ex1\ave.ctl'
'enableprintd:
\ex1\ave.gmf'
mon=1
while（mon<=12）
'sett'mon''
'dave'
'drawtitleqihouchangof'mon''
'print'
'c'
mon=mon+1
endwhile
'disableprint'
;
气候场图：
一月份高度的气候场呈现南高北低的状态，陆地上的高度场比较稀疏，而在西太平洋上高度场比较密集。
八月份高度的气候场呈现东高西低的状态，在我国东北部以北以及印度东北部出现低压中心，而在赤道西太平洋地区出现高压中心。
35°N以北高度分布很密集，而35°N以南比较稀疏。
给jp配的ctl文件：
dset^d:
\ex1\jp.grd
undef-9.99E+33
titleNCEP/NCARREANALYSISPROJECT
xdef37linear60.0002.500
ydef17linear0.0002.500
zdef1levels500
tdef48linearJAN19821mo
vars1
jp199H500
endvars
给jp配的gs文件：
'reinit'
'opend:
\ex1\jp.ctl'
'enableprintd:
\ex1\jp.gmf'
year=1982
while（year<=1985）
mon=1
while（mon<=12）
'sett'mon''
'djp'
'drawtitlejupingchangof'year'.'mon''
'print'
'c'
mon=mon+1
endwhile
year=year+1
endwhile
'disableprint'
;
距平场图：
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
titleNCEP/NCARREANALYSISPROJECT
xdef37linear60.0002.500
ydef17linear0.0002.500
zdef1levels500
tdef12linearJAN198212mo
vars1
s199H500
endvars
给s配的gs文件：
'reinit'
'opend:
\ex1\s.ctl'
'enableprintd:
\ex1\s.gmf'
mon=1
while（mon<=12）
'sett'mon''
'ds'
'drawtitlejunfangchachangof'mon''
'print'
'c'
mon=mon+1
endwhile
'disableprint'
;
均方差场图：
一月份高度的均方差场整体呈现南小北大的状态。
说明低纬地区高度的波动幅度比较小，而中高纬地区高度的波动比较大。
八月份高度的均方差场在亚洲大陆西部有极大值，在30°N处包括赤道-30°N、60°E-85°E这些区域高度的波动幅度比较小，30°N以南以北地区高度的波动幅度较大。
实习二：
相关系数
Fortran程序如下：
programex2
integer,parameter:
:
n=20,p=10
integeri,j,t1,t2,t3
reala（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
dataa/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/
datab/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/
datajpa/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
datajpb/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
datazxfc1/0,0,0,0,0,0,0,0,0,0/
datazxgxs1/0,0,0,0,0,0,0,0,0,0/
datazxfc2/0,0,0,0,0,0,0,0,0,0/
datazxgxs2/0,0,0,0,0,0,0,0,0,0/
datalhxfc/0,0,0,0,0,0,0,0,0,0/
datalhxgxs/0,0,0,0,0,0,0,0,0,0/
!
求均值
doi=1,n
s1=s1+a（i）
enddo
ave1=s1/n
doi=1,n
s2=s2+b（i）
enddo
ave2=s2/n
!
求距平
doi=1,n
jpa（i）=a（i）-ave1
enddo
doi=1,n
jpb（i）=b（i）-ave2
enddo
!
求相关系数
doi=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
!
求自协方差
doj=1,p
doi=1,n-j
zxfc1（j）=zxfc1（j）+jpa（i）*jpa（i+j）
enddo
zxfc1（j）=zxfc1（j）/（n-j）
enddo
doj=1,p
doi=1,n-j
zxfc2（j）=zxfc2（j）+jpb（i）*jpb（i+j）
enddo
zxfc2（j）=zxfc2（j）/（n-j）
enddo
!
自相关系数
doi=1,p
zxgxs1（i）=zxfc1（i）/fc1
zxgxs2（i）=zxfc2（i）/fc2
enddo
!
落后交叉协方差
doj=1,p
doi=1,n-j
lhxfc（j）=lhxfc（j）+jpa（i）*jpb（i+j）
enddo
lhxfc（j）=lhxfc（j）/（n-j）
enddo
!
落后相关系数
doi=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
doi=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程序如下：
programex3
integer,parameter:
:
n=17
integeri,j
real:
:
suma1=0,suma2=0,avea1,avea2,s1=0,s2=0,b0,b,x=14.5,y
integera1（n）
reala2（n）,a3（n）
dataa1/0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/
dataa2/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/
dataa3/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'）
!
求平均值
doi=1,n
suma1=suma1+a1（i）
suma2=suma2+a2（i）
enddo
avea1=suma1/n
avea2=suma2/n
!
求b,b0值
doi=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
!
求回归数据
doi=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）
1000format（i7）
2000format（f8.1）
end
程序运行：
所给数据做出的线性回归曲线的斜率是-0.3012，截距是29.3804，说明所给数据y随着x递减。
实习四：
滑动平均
Fortran程序如下：
PROGRAMMA
integeri
integer,parameter:
:
n=85,ih=11,nyear=1922
integer:
:
yr（n）=0
realX（n）,X1（n）
real:
:
s（75）=0
C**********************************************
C*N:
SAMPLESIZEOFTHETIMESERIES*
C*IH:
MOVINGLENGTH*
C*NYEAR:
FIRSTYEAROFTHESERIES*
C*X（N）:
OROGINALTIMESERIES*
C*X1（N-IH+1）:
MOVEDSERIES*
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'）
!
年份
doi=1,n
yr（i）=1922-1+i
enddo
!
读入数据
READ（2,*）（X（I）,I=1,N）
!
计算滑动平均
doi=1,11
s
（1）=s
（1）+x（i）
enddo
x1（6）=s
（1）/ih
doi=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
PROGRAMEOFPW
CTHISPROGRAMUSESEOFFORANALYSINGTIMESERIES
COFMETEOROLOGICALFIELD
CM:
LENTHOFTIMESERIES
CN:
NUMBEROFGRID-POINTS
CKS=-1:
SELF;KS=0:
DEPATURE;KS=1:
STANDERDLIZEDDEPATURE
CKV:
NUMBEROFEIGENVALUESWILLBEOUTPUT
CKVT:
NUMBEROFEIGENVECTORSANDTIMESERIESWILLBEOUTPUT
CMNH=MIN（M,N）
CEvf=EIGENVACTORS,tcF=TIMECOEFFICIENTSFOREGVT.
CER（KV,1）=LAMDA,LAMDAEIGENVALUE
CER（KV,2）=ACCUMULATELAMDA
CER（KV,3）=THESUMOFCOMPONENTSVECTORSPROJECTEDONTO
cEIGENVACTOR.
CER（KV,4）=ACCUMULATEER（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
DIMENSIONF（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）
CCCCCCCCCCCCCCCCINPUTDATA
open（11,file='d:
\5\sstpx.grd',form='unformatted',
!
access='direct',recl=nx*ny）
do132it=1,m
read（11,rec=it）（（data（i,j）,i=1,nx）,j=1,ny）
do132jj=1,ny
do132ii=1,nx
kkkk=nx*（jj-1）+ii
f（kkkk,it）=data（ii,jj）
132continue
close（11）
CCCCCCCCCCCCCCCCINPUTDATACCCCCCCCCCCCCCCCCCC
ccccccccccccccccccccccccccccccccccccc
CALLTest1（n,m,ff,f,nf）
write（*,*）'ok2'
CALLTRANSF（N,M,F,nf,AVF,DF,KS）
write（*,*）'ok3'
CALLFORMA（N,M,MNH,F,A）
write（*,*）'ok4'
CALLJCB（MNH,A,S,0.00001）
write（*,*）'ok5'
CALLARRANG（KV,MNH,A,ER,S）
write（*,*）'ok6'
CALLTCOEFF（KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER）
write（*,*）'ok7'
calltest3（N,ff,nf,evf,kvt）
write（*,*）'ok8'
open（21,file='d:
\5\evf.grd',form='unformatted',
!
access='direct',recl=nx*ny）
irec=0
do668kk=1,kvt
irec=irec+1
668write（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
升级会员 