气象统计实习报告.docx

资源描述

气象统计实习报告.docx

《气象统计实习报告.docx》由会员分享，可在线阅读，更多相关《气象统计实习报告.docx（38页珍藏版）》请在冰点文库上搜索。

气象统计实习报告.docx

气象统计实习报告

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

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\h500.dat'）

doiy=1,4

dom=1,12

read（5,1000）

read（5,3000）（（var（i,j,m,iy）,i=1,nx）,j=1,ny）

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

求距平场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

求均方差场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

写入气候场

open（10,file='d:

\study\form\shixione\at.grd',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\xd.grd',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

close（11）

！

写入均方差场

open（12,file='d:

\study\form\shixione\sx.grd',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（37F6.2）

3000format（37f8.1）

4000format（37f7.2）

endprogramsx1

⏹运行结果：

⏹Grads文件

气候场

'reinit'

'enableprintd:

\study\form\shixione\at.gmf'

'opend:

\study\form\shixione\at.ctl'

'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\sx.gmf'

'opend:

\study\form\shixione\sx.ctl'

'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\sx.gmf'

'opend:

\study\form\shixione\sx.ctl'

'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=0.47

年平均气温在滞后长度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/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/

DATAw/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/

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=',f5.2）",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=',f7.4,/,'滞后时间长度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=',f7.4,/,'滞后时间长度k=',I2）",rxy_tw（k）,k

k=0

END

⏹运行成果：

*实习四求给定数据的一元线性回归方程

实习时间：

第12周周三1、2节

利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。

年份

气温T

环流指标

1951

0.9

1952

1.2

1953

2.2

1954

2.4

1955

-0.5

1956

2.5

1957

-1.1

1958

1959

6.2

1960

2.7

1961

3.2

1962

-1.1

1963

2.5

1964

1.2

1965

1.8

1966

0.6

1967

2.4

1968

2.5

1969

1.2

1970

-0.8

答案：

F=20.18>Fα=4.41，回归方程显著

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.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/

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>4.41）print*,'F=',f,'>Fα=4.41，回归方程显著

if（f<=4.41）print*,'F=',f,'

End

⏹运行成果：

*实习七计算给定数据的11年滑动平均和累积距平

实习时间：

第15周周三1、2节

利用数据ma.dat，编写11点滑动平均的程序，ma.for给出了阅读资料的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\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）

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

doi=1,n

write（4,*）x2（i）

enddo

End

⏹图形显示：

*实习八对给定的海温数据进行EOF分析

实习时间：

第16周周三1、2节

给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF分析，资料的时空范围可以根据sstpx.ctl获知。

数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。

Zhunsst.for给出了如何读取资料，

Ssteof.for为对距平或者标准化距平处理后的资料进行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\sstpx.grd',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

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

standardizing

doj=1,ny

doi=1,nx

dok=1,mo

doit=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）

endif

enddo

outputfileopen（2,file='d:

\study\form\shixieight\standard.grd',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

（2）

（1）

end

⏹运行程序后：

得到zhunsst.grd文件，即图中standard

下面为老师给出程序，为运行方便将部分简略。

PROGRAMEOFPW

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）

open（11,file='d:

\study\form\shixieight\zhunsst.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）

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:

\study\form\shixieight\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:

\study\form\shixieight\tcf.grd',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（10f8.4）

open（21,file='d:

\study\form\shixieight\dats.dat'）

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）=0.0

enddo

doi=1,n1

dok=1,m

if（f（i,k）.eq.ff）then

f（i,k）=0.0

nf（i）=nf（i）+1

endif

enddo

return

end

SUBROUTINETRANSF（n1,m,f,nf,avf,df,ks）

DIMENSIONF（N1,M）,AVF（N1）

DIMENSIONDF（N1）

DIMENSIONnF（N1）

if（ks.eq.-1）then

goto30

endif

doi

展开阅读全文