气象统计方法实习BD文档格式.docx
《气象统计方法实习BD文档格式.docx》由会员分享，可在线阅读，更多相关《气象统计方法实习BD文档格式.docx（37页珍藏版）》请在冰点文库上搜索。
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
计算距平场
doiy=1,4
jp（i,j,m,iy）=var（i,j,m,iy）-ave（i,j,m）
计算均方差场
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
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）
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'
'
print'
c'
mon=mon+1
endwhile
disableprint'
;
气候场图：
一月份高度的气候场呈现南高北低的状态，陆地上的高度场比较稀疏，而在西太平洋上高度场比较密集。
八月份高度的气候场呈现东高西低的状态，在我国东北部以北以及印度东北部出现低压中心，而在赤道西太平洋地区出现高压中心。
35°
N以北高度分布很密集，而35°
N以南比较稀疏。
给jp配的ctl文件：
\ex1\jp.grd
tdef48linearJAN19821mo
jp199H500
给jp配的gs文件：
\ex1\jp.ctl'
\ex1\jp.gmf'
year=1982
while（year<
=1985）
djp'
drawtitlejupingchangof'
year'
.'
year=year+1
距平场图：
1983年6月距平场在日本地区出现低压中心，在我国南部出现高压中心,在亚洲西北部也有高压中心。
赤道至25°
N间以及25°
N-40°
N，60°
E-100°
E间基本都是正距平，而在25°
N，100°
E-150°
E间基本都是负距平。
1984年7月距平场在亚洲大陆西部、日本地区、赤道西太平洋地区形成低压中心，太平洋西北部形成高压中心。
给s配的ctl文件：
\ex1\s.grd
s199H500
给s配的gs文件：
\ex1\s.ctl'
\ex1\s.gmf'
ds'
drawtitlejunfangchachangof'
均方差场图：
一月份高度的均方差场整体呈现南小北大的状态。
说明低纬地区高度的波动幅度比较小，而中高纬地区高度的波动比较大。
八月份高度的均方差场在亚洲大陆西部有极大值，在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
s2=s2+b（i）
ave2=s2/n
求距平
jpa（i）=a（i）-ave1
jpb（i）=b（i）-ave2
求相关系数
sum1=sum1+jpa（i）*jpa（i）
sum2=sum2+jpb（i）*jpb（i）
sum3=sum3+jpa（i）*jpb（i）
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）
zxfc1（j）=zxfc1（j）/（n-j）
zxfc2（j）=zxfc2（j）+jpb（i）*jpb（i+j）
zxfc2（j）=zxfc2（j）/（n-j）
自相关系数
doi=1,p
zxgxs1（i）=zxfc1（i）/fc1
zxgxs2（i）=zxfc2（i）/fc2
落后交叉协方差
lhxfc（j）=lhxfc（j）+jpa（i）*jpb（i+j）
lhxfc（j）=lhxfc（j）/（n-j）
落后相关系数
lhxgxs（i）=lhxfc（i）/（sqrt（fc1*fc2））
年平均气温不同滞后时刻所对应的自相关系数为：
print*,（（zxgxs1（i）,'
'
）,i=1,P）
print*
冬季平均气温不同滞后时刻所对应的自相关系数为：
print*,（（zxgxs2（i）,'
年平均气温和冬季平均气温之间不同滞后时刻所对应的落后交叉相关系数为：
print*,（（lhxgxs（i）,'
最大绝对值
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
print*,'
中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度为：
t1,zxgxs1（t1）
中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度为：
t2,zxgxs2（t2）
中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度为：
t3,lhxgxs（t3）
end
程序运行如下：
中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度是7，自相关系数是-0.3724，说明在这些年的数据中，滞后7年的序列与原序列的相关系数绝对值最大，呈反相关。
中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度是4，自相关系数是-0.3678，说明在这些年的数据中，滞后4年的序列与原序列的相关系数绝对值最大，呈反相关。
中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度是3，自相关系数是-0.4066，说明在这些年的数据中，滞后3年的冬季平均气温序列与原平均气温序列的相关系数绝对值最大，呈反相关。
实习三：
一元线性回归
programex3
n=17
integeri,j
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='
\3\ys.txt'
open（15,file='
\3\ys.grd'
open（6,file='
\3\hs.txt'
open（16,file='
\3\hs.grd'
求平均值
suma1=suma1+a1（i）
suma2=suma2+a2（i）
avea1=suma1/n
avea2=suma2/n
求b,b0值
s1=s1+a1（i）*a2（i）
s2=s2+a1（i）*a1（i）
b=（s1-suma1*suma2/n）/（s2-suma1*suma1/n）
b0=avea2-avea1*b
y=b0+b*x
y='
b0,'
+'
b,'
x'
print*
print*,y
求回归数据
a3（i）=b0+b*a1（i）
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）
程序运行：
所给数据做出的线性回归曲线的斜率是-0.3012，截距是29.3804，说明所给数据y随着x递减。
实习四：
滑动平均
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*
OPEN（2,FILE='
\4\MA.DAT'
\4\ma.grd'
OPEN（3,FILE='
\4\hd.DAT'
open（13,file='
\4\hd.grd'
!
年份
doi=1,n
yr（i）=1922-1+i
读入数据
READ（2,*）（X（I）,I=1,N）
计算滑动平均
doi=1,11
s
（1）=s
（1）+x（i）
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
写入数据
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（13）
原数据和滑动后数据的图形：
由图可知，所给数据在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）
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
\5\sstpx.grd'
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
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）
ok4'
CALLJCB（MNH,A,S,0.00001）
ok5'
CALLARRANG（KV,MNH,A,ER,S）
ok6'
CALLTCOEFF（KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER）
ok7'
calltest3（N,ff,nf,evf,kvt）
ok8'
open（21,file='
\5\evf.grd'
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）
\5\tcf.grd'
recl
升级会员 