气象统计实习报告.docx
《气象统计实习报告.docx》由会员分享,可在线阅读,更多相关《气象统计实习报告.docx(32页珍藏版)》请在冰豆网上搜索。
![气象统计实习报告.docx](https://file1.bdocx.com/fileroot1/2022-11/22/67b1e9d1-3e26-49ef-b90b-e05194cbeb42/67b1e9d1-3e26-49ef-b90b-e05194cbeb421.gif)
气象统计实习报告
?
?
实?
习?
报?
告?
书
课程名称:
气象统计方法课程实践
姓名:
学号:
班级:
级气科班
\\*实习一求500hPa高度场气候场、距平场和均方差场
实习时间:
第9周周三1、2节
1.资料介绍
有一500hPa高度场资料,文件名,范围:
60~150E,0~40N.
时段:
~共48个月。
水平分辨率:
*,格点数:
37*17。
2.要求
编fortran程序,求500hPa高度场的
(1)气候场;
(2)距平场;
(3)均方差场。
并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图,图注要清楚,即要注明是哪个时间的图形,并做简单分析。
注:
给出了如何用fortran读取ASCII码资料.
⏹FORTRAN
programsx1
implicitnone
integernx,ny,mo,yr
parameter(nx=37,ny=17,mo=12,yr=4)
realvar(nx,ny,mo,yr)
realat(nx,ny,mo),xd(nx,ny,mo,yr),sx(nx,ny,mo)
integeri,j,m,t,it,iy,irec
open(5,file='d:
\study\form\shixione\')
doiy=1,4
dom=1,12
read(5,1000)
read(5,3000)((var(i,j,m,iy),i=1,nx),j=1,ny)
enddo
enddo
close(5)
!
计算气候场at
dot=1,12
doj=1,ny
doi=1,nx
at(i,j,t)=0
doit=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
dot=1,12
doj=1,ny
doi=1,nx
xd(i,j,t,1)=0
doit=1,4
xd(i,j,t,it)=var(i,j,t,it)-at(i,j,t)
enddo
enddo
enddo
enddo
!
求均方差场sx
dot=1,12
doj=1,ny
doi=1,nx
sx(i,j,t)=0
doit=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\',form='unformatted',access='direct',recl=nx*ny)
irec=0
dot=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\',form='unformatted',access='direct',recl=nx*ny)
irec=0
doit=1,4
dot=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\',form='unformatted',access='direct',recl=nx*ny)
irec=0
dot=1,12
irec=irec+1
write(12,rec=irec)((sx(i,j,t),i=1,nx),j=1,ny)
enddo
close(12)
1000format(2i7)
2000format
3000format
4000format
endprogramsx1
⏹运行结果:
⏹Grads文件
气候场
'reinit'
'enableprintd:
\study\form\shixione\'
'opend:
\study\form\shixione\'
'setgridoff'
'setgradsoff'
'setlat040'
'setlon60150'
'setlev500'
mon=1
while(mon<=12)
'sett'mon''
'dh'
'drawtitle1982year'mon'month'
'print'
'c'
mon=mon+1
endwhile
'disableprint'
;
距平场
'reinit'
'enableprintd:
\study\form\shixione\'
'opend:
\study\form\shixione\'
'setgridoff'
'setgradsoff'
'setlat040'
'setlon60150'
'setlev500'
year=1982
while(year<=1984)
mon=1
while(mon<=12)
'sett'mon''
'dh'
'drawtitle500hPa'year'year'mon'monthanomaly'
'print'
'c'
mon=mon+1
endwhile
year=year+1
endwhile
'disableprint';
均方差
'reinit'
'enableprintd:
\study\form\shixione\'
'opend:
\study\form\shixione\'
'setgridoff'
'setgradsoff'
'setlat040'
'setlon60150'
'setlev500'
mon=1
while(mon<=12)
'sett'mon''
'dh'
'drawtitle500hPa1982year'mon'monthMean-squareDeviation'
'print'
'c'
mon=mon+1
endwhile
'disableprint'
;
⏹运行结果:
上面图中只展示了部分而未全部添加
*实习二计算给定数据资料的简单相关系数和自相关系数
实习时间:
第10周周三1、2节
根据下表中年平均气温和冬季平均气温的等级数据进行下列计算:
1)计算两个气温之间的简单相关系数。
2)分别找出两个气温数据自相关系数绝对值最大的滞后时间长度。
(滞后长度τ最大取10)
要求:
实习报告中附出简单相关系数或自相关系数程序。
答案:
r=
年平均气温在滞后长度j=7、冬季序列在j=4最大。
⏹FORTRAN
计算相关系数r
PROGRAMEXAM
IMPLICITNONE
INTEGER,PARAMETER:
:
N=20
INTEGERi,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
REALy(N),w(N)
DATAy/,,,,,,,,,,,,,,,,,,,
DATAw/,,,,,,,,,,,,,,,,,,,
REALsyy(N),sww(N),r(N),rty(N),rtw(N),rtyw(N),rxy_ty(N),rxy_tw(N),rxy_tyw(N)
!
求两数组平均值
DOi=1,N
avr_y=avr_y+y(i)
avr_w=avr_w+w(i)
ENDDO
avr_y=avr_y/N
avr_w=avr_w/N
!
简单相关系数
DOj=1,N
syy(j)=(y(j)-avr_y)**2
sy=sy+syy(j)
sww(j)=(w(j)-avr_w)**2
sw=sw+sww(j)
ENDDO
sy=sqrt(sy/N)
sw=sqrt(sw/N)
DOj=1,N
r(j)=((y(j)-avr_y)/sy)*((w(j)-avr_w)/sw)
rxy=rxy+r(j)
ENDDO
rxy=rxy/N
PRINT"(/'1970-1989年全年平均气温与冬季平均气温的简单相关系数rxy=',",rxy
k=0
!
自相关系数
DOty=1,N/2
DOi=1,N-ty
rty(i)=((y(i)-avr_y)/sy)*((y(i+ty)-avr_y)/sy)
rxy_ty(ty)=rxy_ty(ty)+rty(i)
ENDDO
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
ENDIF
ENDDO
PRINT"('全年平均气温绝对值最大自相关系数rxy_ty=',,/,'滞后时间长度k=',I2)",rxy_ty(k),k
k=0
DOtw=1,N/2
DOi=1,N-tw
rtw(i)=((w(i)-avr_w)/sw)*((w(i+tw)-avr_w)/sw)
rxy_tw(tw)=rxy_tw(tw)+rtw(i)
ENDDO
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
ENDIF
ENDDO
PRINT"('冬季平均气温绝对值最大自相关系数rxy_tw=',,/,'滞后时间长度k=',I2)",rxy_tw(k),k
k=0
END
⏹运行成果:
*实习四求给定数据的一元线性回归方程
实习时间:
第12周周三1、2节
利用下表数据,以环流指标为预报因子,气温为预报量,计算气温和环流指标之间的一元线性回归方程,并对回归方程进行检验。
年份
气温T
环流指标
1951
32
1952
25
1953
20
1954
26
1955
27
1956
24
1957
28
1958
0
24
1959
15
1960
16
1961
24
1962
30
1963
22
1964
30
1965
24
1966
33
1967
26
1968
20
1969
32
1970
35
答案:
F=>Fα=,回归方程显着
programsx4
implicitnone
integeri
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,,,,,,,,,,,,
doi=1,20
f1=f1+x(i)*y(i)
enddo
avx=sum(x)/20!
求均值|
avy=sum(y)/20
doi=1,20
f2=f2+x(i)**2
enddo
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'
doi=1,20
sx=sx+(x(i)-avx)**2
sy=sy+(y(i)-avy)**2
sxy=sxy+(x(i)-avx)*(y(i)-avy)
enddo
sx=sqrt(sx/20)
sy=sqrt(sy/20)
sxy=sxy/20
rxy=sx*b/sy
f=rxy**2*18/(1-rxy**2)
if(f>print*,'F=',f,'>Fα=,回归方程显着
if(f<=print*,'F=',f,'End
⏹运行成果:
*实习七计算给定数据的11年滑动平均和累积距平
实习时间:
第15周周三1、2节
利用数据,编写11点滑动平均的程序,给出了阅读资料的fortran程序。
数据在文件夹中单独给出。
要求:
实习报告中附出程序,并给出原数据和滑动后数据的图形(1张图)和累积距平数据图形(1张图)
programsx7
implicitnone
integern,ih,nyear
parameter(n=85,ih=11,nyear=1922)
real:
:
x(n),x1(n-ih+1)=0,x2(n)=0,averx=0
integeri,j
open(2,file='d:
\study\form\shiyanseven\')
open(3,file='d:
\study\form\shiyanseven\')
open(4,file='d:
\study\form\shiyanseven\')
read(2,*)(x(i),i=1,n)
doj=1,n-ih+1!
movingaverage
doi=1,ih
x1(j)=x1(j)+x(i+j-1)
enddo
x1(j)=x1(j)/ih
enddo
doi=1,n-ih+1
write(3,*)x1(i)
enddo
doi=1,n!
averageofx
averx=averx+x(i)
enddo
averx=averx/n
doi=1,n!
accelarate
doj=1,i
x2(i)=x2(i)+(x(j)-averx)
enddo
enddo
doi=1,n
write(4,*)x2(i)
enddo
End
⏹图形显示:
*实习八对给定的海温数据进行EOF分析
实习时间:
第16周周三1、2节
给出海表温度距平数据资料,以及相应的数据描述文件,对其进行EOF分析,资料的时空范围可以根据获知。
数据在文件夹中单独给出,距平或者标准化距平处理后再进行EOF。
给出了如何读取资料,
为对距平或者标准化距平处理后的资料进行EOF分析。
要求:
实习报告中给出第一特征向量及其时间系数,并分析其时空特征。
programmain
parameter(mo=12,my=43,nx=18,ny=12,mt=516)
dimensionavf(mo,nx,ny),df(mo,nx,ny)
dimensionsst(nx,ny,mt),sst3(nx,ny,mt)
open(1,file='d:
\study\form\shixieight\',form='unformatted',access='direct',recl=nx*ny)
irec=1
doit=1,mt
read(1,rec=irec)((sst(i,j,it),i=1,nx),j=1,ny)
irec=irec+1
enddo
!
average
doj=1,ny
doi=1,nx
dok=1,mo
doit=k,mt,12
avf(k,i,j)=avf(k,i,j)+sst(i,j,it)
enddo
avf(k,i,j)=avf(k,i,j)/my
enddo
enddo
enddo
!
variance
doj=1,ny
doi=1,nx
dok=1,mo
doit=k,mt,12
df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j))**2
enddo
df(k,i,j)=sqrt(df(k,i,j)/my)
enddo
enddo
enddo
!
standardizing
doj=1,ny
doi=1,nx
dok=1,mo
doit=k,mt,12
if(sst(i,j,it)==then
sst3(i,j,it)=
else
sst3(i,j,it)=(sst(i,j,it)-avf(k,i,j))/df(k,i,j)
endif
enddo
enddo
enddo
enddo
!
outputfileopen(2,file='d:
\study\form\shixieight\',form='unformatted',access='direct',recl=nx*ny)
irec=1
doit=1,mt
write(2,rec=irec)((sst3(i,j,it),i=1,nx),j=1,ny)
irec=irec+1
enddo
close
(2)
close
(1)
end
⏹运行程序后:
得到文件,即图中standard
下面为老师给出程序,为运行方便将部分简略。
PROGRAMEOFPW
PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2)
PARAMETER(ff=,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)
open(11,file='d:
\study\form\shixieight\',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)
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,
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:
\study\form\shixieight\',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:
\study\form\shixieight\',form='unformatted',access='direct',recl=kvt)
irec=0
do345it=1,m
irec=irec+1
345write(21,rec=irec)(tcf(it,iik),iik=1,kvt)
close(21)
106format
open(21,file='d:
\study\form\shixieight\')
write(21,106)(er(iiii,3),iiii=1,kv)
close(21)
STOP
END
SUBROUTINETest1(n1,m,ff,f,nf)
DIMENSIONF(N1,M)
DIMENSIONnF(N1)
doi=1,n1
nf(i)=
enddo
doi=1,n1
dok=1,m
if(f(i,k).then
f(i,k)=
nf(i)=nf(i)+1
endif
enddo
enddo
return
end
SUBROUTINETRANSF(n1,m,f,nf,avf,df,ks)
DIMENSIONF(N1,M),AVF(N1)
DIMENSIONDF(N1)
DIMENSIONnF(N1)
ifgoto30
endif
doi=1,n1
avf(i)=
enddo
ifgoto5
endif
doi=1,n1
df(i)=
enddo
5continue
DO141I=1,N1
if(nf(i).goto141
do12j=1,m
12AVF(I)=AVF(I)+F(I,J)
AVF(I)=AVF(I)/M
DO14J=1,M
14F(I,J)=F(I,J)-AVF(I)
141CONTINUE
IFTHEN
RETURN
ELSE
DO241I=1,N1
if(nf(i).goto241
DO22J=1,M
22DF(I)=DF(I)+F(I,J)*F(I,J)
DF(I)=SQRT(DF(I)/M)
DO24J=1,M
24F(I,J)=F(I,J)/DF(I)
241CONTINUE
ENDIF
30CONTINUE
RETURN
END
SUBROUTINEFORMA(N,M,MNH,F,A)
DIMENSIONF(N,M),A(MNH,MNH)
IF(M-N)40,50,50
40DO44I=1,MNH
DO44J=I,MNH
A(I,J)=
DO42IS=1,N
42A(I,J)=A(I,J)+F