气象统计实习报告.docx

气象统计实习报告.docx

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

气象统计实习报告

实 习 报 告 书

课程名称：

气象统计方法课程实践

姓名:

学号:

班级：

级气科班

\\＊实习一求500hPa高度场气候场、距平场和均方差场

实习时间：

第9周周三1、2节

1.资料介绍

有一500hPa高度场资料,文件名h500。

dat,范围：

60～150E，0~40N。

时段：

1982。

1~1985。

12共48个月。

水平分辨率：

2.5＊2.5,格点数:

37*17。

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

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\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

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

32

1952

1。

2

25

1953

2。

2

20

1954

2.4

26

1955

—0.5

27

1956

2.5

24

1957

—1.1

28

1958

0

24

1959

6.2

15

1960

2.7

16

1961

3.2

24

1962

-1。

1

30

1963

2.5

22

1964

1。

2

30

1965

1。

8

24

1966

0.6

33

1967

2.4

26

1968

2。

5

20

1969

1。

2

32

1970

-0.8

35

答案:

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，’

41，回归方程不显著

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

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

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

enddo

enddo

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

close

（2）

close

（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

enddo

return

end

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

DIMENSIONF（N1,M），AVF（N1）

DIMENSIONDF（N1）

DIMENSIONnF（N1）

if（ks.eq。

-