气象统计实习报告.docx

气象统计实习报告.docx

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

气象统计实习报告

?

?

实?

习?

报?

告?

书

课程名称：

气象统计方法课程实践

姓名：

学号：

班级：

级气科班

\\*实习一求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