气候变化及其诊断.docx
- 文档编号:29077744
- 上传时间:2023-07-20
- 格式:DOCX
- 页数:32
- 大小:478.77KB
气候变化及其诊断.docx
《气候变化及其诊断.docx》由会员分享,可在线阅读,更多相关《气候变化及其诊断.docx(32页珍藏版)》请在冰豆网上搜索。
气候变化及其诊断
气候变化及其诊断
实习报告
学院数学与统计学院
专业09大气科学
学号20091338017
姓名谷天宇
任课教师张春莹
实习一气候场、距平场、标准化场,标准差场
提供500hPa位势高度场资料(1982年1月-1985年12月),选取1月和7月
1气候场(.gs)
'reinit'
'openf:
\gu\eof\tcf.ctl'
'enableprintf:
\gu\eof\tcf1.gmf'
'sett1516'
'dt1'
'print'
'disableprint'
pulldummy
'enableprintf:
\gu\eof\tcf2.gmf'
'c'
'dt2'
'print'
'disableprint'
;
图.一月七月气候场
2结果分析:
一月份显示,500hpa等位势高度线基本与纬度线平行,且随着纬度的增加而增加,但总体气压梯度较小,这是因为,地面接收的太阳辐射随着纬度增加而减少,由南往北温度逐渐降低。
如上图所示,7月份,西太平洋副高位于我国东南部25N附近,形成闭合区间,气压梯度强度较大。
2距平场
'reinit'
'openf:
\gu\h500\h500.ctl'
'enableprintf:
\gu\h500\jp.gmf'
'setgradsoff'
'setgridoff'
'setlat040'
'setlon60150'
'setlev500'
'sett112'
'definehave=ave(h,t+0,t=48,12)'
'modifyhaveseasonal'
'sett148'
'definehano=h-have'
'sett148'
it=1
while(it<=48)
'sett'it''
'dhano'
'drawtitlet='it''
'print'
'c'
it=it+1
endwhile
'disableprint'
;
图.一月七月距平场
结果分析如下:
如图所示,一月份所显示,北半球由于海陆相间,冬季陆地作为冷源,容易形成高压中心,在35N附近形成的正的气压距平中心(西伯利亚高压);海洋相对而言就是热源,易形成低压中心,所以在西太平洋区存在负的气压距平中心。
七月份图,在夏季北半球大陆作为热源,易形成低压中心,所以存在负的气压距平中心,海洋作为冷源,容易形成高压中心,存在正的气压距平中心。
3均方差场
'reinit'
'openf:
\gu\h500\h500.ctl'
'setgridoff'
'setgradsoff'
'enableprintf:
\gu\h500\jfc.gmf'
'sett112'
'definehave=ave(h,t+0,t=48,12)'
'modifyhaveseasonal'
'sett148'
'definehano=h-have'
'sett112'
'defineSD=ave(hano*hano,t+0,t=48,12)'
'sett148'
'definestd=sqrt(SD)'
'sett112'
it=1
while(it<=12)
'sett'it''
'dstd'
'drawtitlejunfangcha(t='it')'
'print'
'c'
it=it+1
endwhile
'disableprint'
;
图.一月七月均方差图
结果分析:
图一月份的显示,北半球在陆上存在两个最低气压均方差中心(105E,18N),(115E,23N),说明这两个地方的气压值较为稳定,活动幅度较小。
因为副热带高压的移动范围为10°N-40°N之间,所以这个两个地区一直处在副高控制下,气象一直较高,变化幅度不大。
七月份图,夏季北半球西太平洋上存在较大的气压均方差中心,说明此处气压值变化程度较大。
4标准场
'reinit'
'openf:
\gu\h500\h500.ctl'
'setgridoff'
'setgradsoff'
'enableprintf:
\gu\h500\bzc.gmf'
'sett112'
'definehave=ave(h,t+0,t=48,12)'
'modifyhaveseasonal'
'sett148'
'definehano=h-have'
'sett112'
'defineSD=ave(hano*hano,t+0,t=48,12)'
'definestd=sqrt(SD)'
'modifystdseasonal'
'sett148'
'definestdhano=hano/std'
'sett148'
it=1
while(it<=48)
'sett'it''
'dstdhano'
'drawtitlebianzhun(t='it')'
'print'
'c'
it=it+1
endwhile
'disableprint'
图.一月七月标准差场
结果分析:
如图所示,1月份时,亚欧大陆被西伯利亚高压控制,形成一个气压高值中心。
而海洋相对是热源,有利于低压的形成,如北太平洋的阿留申低压,所以形成负的气压中心。
7月份(夏季)相反,北半球大陆是热源,形成低压。
如(30°N,115°E)形成东西较长的负的气压距平中心。
而我国的青藏高原地区存在高压中心。
实习二相关系数
提供中国160站7月降水资料(1951年-2010年),选取自己所在地和邻近点
1提取石家庄天津(1951年-2010年)七月份降水
parameter(n=160,m=60)
integerit(n,m),in,im
open(11,file='f:
\gu\160\R16007.txt')
open(12,file='f:
\gu\160\r160sjz.txt')
open(13,file='f:
\gu\160\r160tj.txt')
read(11,*)((it(in,im),in=1,n),im=1,m)
in=37!
天津
doim=1,m
write(13,*)it(in,im)
enddo
in=35!
石家庄
doim=1,m
write(12,*)it(in,im)
enddo
end
2求天津与石家庄的降水的相关系数
programxianguanxishu
parameter(m=60)
integer:
:
i,j
dimension:
:
t2(m),t3(m)
real:
:
sum2=0.0,sum3=0.0,ave2,ave3,s2,s3,temp2,temp3,sum12,s12,r12
open(2,file='f:
\gu\160\r160sjz.txt')
open(3,file='f:
\gu\160\r160tj.txt')
doi=1,m
read(2,*)t2(i)
enddo
doi=1,m
read(3,*)t3(i)
enddo
doi=1,m
sum2=sum2+t2(i)
ave2=sum2/m
sum3=sum3+t3(i)
ave3=sum3/m
enddo
sum2=0.0
sum3=0.0
doi=1,m
temp2=(t2(i)-ave2)**2
sum2=sum2+temp2
temp3=(t3(i)-ave3)**2
sum3=sum3+temp3
enddo
s2=sqrt(sum2/m)!
石家庄方差
s3=sqrt(sum3/m)!
天津方差
doi=1,m
sum12=sum12+(t2(i)-ave2)*(t3(i)-ave3)
enddo
s12=sum12/m!
石家庄与天津的协方差
r12=s12/(s2*s3)
print*,'石家庄和天津降水的相关系数:
',r12
end
图.石家庄和天津1951-2010年7月份降水分布图
利用相关系数的t检验方法,r=0.3396631.
得统计量t=
=2.750311,给定信度范围为α=0.01对应可查的tα,比较t与tα,则t>tα,总体相关,即石家庄降水与天津地区降水存在显著的正相关关系。
实习三一元线形回归
提供数据资料JJA(1948-2008年)
1求取一元线性回归方程
programyiyuanxianxing3
parameter(N=61)
realx(N),y(N),a,b,avex,avey,sxy,sxx,syy,sx,sy,F,xi,yi
open(1,file='f:
\gu\yyxx\jja.txt')
open(2,file='f:
\gu\yyxx\xianxing.txt')
doi=1,N
read(1,*)x(i),y(i)
write(2,*)x(i),y(i)
enddo
avex=0.0
avey=0.0
doi=1,N
avex=avex+x(i)
enddo
avex=avex/N
doi=1,N
avey=avey+y(i)
enddo
avey=avey/N
sxy=0.0
doi=1,N!
协方差
sxy=sxy+(x(i)-avex)*(y(i)-avey)
enddo
sxy=sxy/N
sxx=0.0
syy=0.0
doi=1,N!
方差
sxx=sxx+(x(i)-avex)**2
syy=syy+(y(i)-avey)**2
enddo
sxx=sxx/N
syy=syy/N
b=sxy/sxx
a=avey-b*avex
print*,a,b
print*,'一元回归方程为:
''y=',a,b,'x'!
计算一元线性回归方程
end
程序运行可得:
图.源数据分布图与一元线性回归图
根据
解得一元线性回归方程为:
y=42.29835-2.1384394x,可知,jja的数值随年份的增加逐渐减小(如图一元线性回归图)。
实习四滑动t检验
提供中国160站7月降水资料(1951年-2010年),选取自己所在地
1fortran程序编写
parameter(n=60,n1=10,n2=10,ta=2.85)
parameter(undef=-9999.)
realx(n),t(n)
open(1,file='f:
\gu\hdt\sjz.txt')
open(2,file='f:
\gu\hdt\MovingttestT.txt')
doi=1,n
read(1,*)x(i)
enddo
doi=1,n
t(i)=undef
enddo
do10i=n1,n-n2
x1=0.0
doj1=i-n1+1,i
x1=x1+x(j1)
enddo
x1=x1/real(n1)
s1=0.0
doj1=i-n1+1,i
s1=s1+(x(j1)-x1)**2
enddo
s1=s1/real(n1)
x2=0.0
doj2=i,i+n2-1
x2=x2+x(j2)
enddo
x2=x2/real(n2)
s2=0.0
doj2=i,i+n2-1
s2=s2+(x(j2)-x2)**2
enddo
s2=s2/real(n2)
s=sqrt((n1*s1+n2*s2)/real(n1+n2-2))
t(i)=(x1-x2)/s/sqrt(real(n1+n2)/real(n1*n2))
10continue
doi=n1,n-n2
write(2,*)t(i)
enddo
close
(1)
close
(2)
end
我们利用滑动T检验原理
编写Fortran程序,进行计算。
原始数据有60个,滑动长度取10,利用左边两个公式得到41个滑动统计量,给定显著性水平0.01,查t分布表得到临界值
=2.85,若
,认为基准点前后的两子序列两子序列均值无显著差异,否则认为在基准点时刻出现了突变。
发现41个滑动统计量都在0.01显著性水平临界值范围以内,石家庄在1951-2010年间的7月降水量均值分布正常,没有突变现象发生。
实习五EOF
提供数据资料sstpx.grd
1Fortran程序编写
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='f:
\gu\eof\sstpx.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='f:
\gu\eof\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='f:
\gu\eof\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='f:
\gu\eof\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.-1)then
goto30
endif
doi=1,n1
avf(i)=0.0
enddo
if(ks.eq.0)then
goto5
endif
doi=1,n1
df(i)=0.0
enddo
5continue
DO141I=1,N1
if(nf(i).ne.0)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
IF(KS.EQ.0)THEN
RETURN
ELSE
DO241I=1,N1
if(nf(i).ne.0)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)=0.0
DO42IS=1,N
42A(I,J)=A(I,J)+F(IS,I)*F(IS,J)
A(J,I)=A(I,J)
44CONTINUE
RETURN
50DO54I=1,MNH
DO54J=I,MNH
A(I,J)=0.0
DO52JS=1,M
52A(I,J)=A(I,J)+F(I,JS)*F(J,JS)
A(J,I)=A(I,J)
54CONTINUE
RETURN
END
SUBROUTINEJCB(N,A,S,EPS)
DIMENSIONA(N,N),S(N,N)
DO30I=1,N
DO30J=1,I
IF(I-J)20,10,20
10S(I,J)=1.
GOTO30
20S(I,J)=0.
S(J,I)=0.
30CONTINUE
G=0.
DO40I=2,N
I1=I-1
DO40J=1,I1
40G=G+2.*A(I,J)*A(I,J)
S1=SQRT(G)
S2=EPS/FLOAT(N)*S1
S3=S1
L=0
50S3=S3/FLOAT(N)
60DO130IQ=2,N
IQ1=IQ-1
DO130IP=1,IQ1
IF(ABS(A(IP,IQ)).LT.S3)GOTO130
L=1
V1=A(IP,IP)
V2=A(IP,IQ)
V3=A(IQ,IQ)
U=0.5*(V1-V3)
IF(U.EQ.0.0)G=1.
IF(ABS(U).GE.1E-10)G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)
ST=G/SQRT(2.*(1.+SQRT(1.-G*G)))
CT=SQRT(1.-ST*ST)
DO110I=1,N
G=A(I,IP)*CT-A(I,IQ)*ST
A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
A(I,IP)=G
G=S(I,IP)*CT-S(I,IQ)*ST
S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT
110S(I,IP)=G
DO120I=1,N
A(IP,I)=A(I,IP)
120A(IQ,I)=A(I,IQ)
G=2.*V2*ST*CT
A(IP,IP)=V1*CT*CT+V3*ST*ST-G
A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
A(IQ,IP)=A(IP,IQ)
130CONTINUE
IF(L-1)150,140,150
140L=0
GOTO60
150IF(S3.GT.S2)GOTO50
RETURN
END
SUBROUTINEARRANG(KV,MNH,A,ER,S)
DIMENSIONA(MNH,MNH),ER(mnh,4),S(MNH,MNH)
TR=0.0
DO200I=1,MNH
TR=TR+A(I,I)
200ER(I,1)=A(I,I)
MNH1=MNH-1
DO210K1=MNH1,1,-1
DO210K2=K1,MNH1
IF(ER(K2,1).LT.ER(K2+1,1))THEN
C=ER(K2+1,1)
ER(K2+1,1)=ER(K2,1)
ER(K2,1)=C
DO205I=1,MNH
C=S(I,K2+1)
S(I,K2+1)=S(I,K2)
S(I,K2)=C
205CONTINUE
ENDIF
210CONTINUE
ER(1,2)=ER(1,1)
DO220I=2,KV
ER(I,2)=ER(I-1,2)+ER(I,1)
220CONTINUE
DO230I=1,KV
ER(I,3)=ER(I,1)/TR
ER(I,4)=ER(I,2)/TR
230CONTINUE
WRITE(6,250)TR
250FORMAT(/5X,'TOTALSQUAREERROR=',F20.5)
RETURN
END
SUBROUTINETCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)
DIMENSIONS(MNH,MNH),F(N,M),V(MNH),ER(mnh,4),evf(n,kvt),tcf(m,kvt)
DO360J=1,KVT
C=0.
DO350I=1,MNH
350C=C+S(I,J)*S(I,J)
C=SQRT(C)
DO160I=1,MNH
s(I,J)=S(I,J)/C
160evf(I,J)=S(I,J)/C
360CONTIN
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 气候变化 及其 诊断