气象程序设计及绘图实习报告

实习要求：
•根据所提供的实习说明和资料，按要求完成实习内容；
•将完成的实习内容按步骤写入实习报告中，包括资料说明、Fortran程序、CTL文件、gs文件、绘制的图形（要求在图形上方用draw title的方式添加标题，注明姓名拼音及学号）
•根据所绘图形及查找相关文献，说明本次实习所求物理量场的意义。

实习内容：
•Nino34海温指数与海平面气压相关图
•Nino34海温指数与中国160站降水相关图
1.绘制1月份Nino34海温指数与1月份海平面气压相关图
1.1 资料说明
➢Nino34指数：Nino34.txt是1951年1月至2013年12月Nino3.4区（热带太平洋：西经170度~西经120度、北纬5度~南纬5度）区域平均海温指数资料，共有63行数据，每行数据第一个数字为年份，后面12个数字为该年1-12月的海温指数；
➢全球海平面气压月平均数据：slp.jan.grd为1951至2013年的1月全球海平面气压场二进制数据，年数为63年，格点数144*73，水平分辨率2.5︒*2.5︒；
1.2步骤：
1用Fortran编写corr.grid.f90文件，计算1951-2013年1月Nino34海温指数与海平面气压相关系数，计算结果保存于corr.slp.grd 文件中；
corr.grid.f90文件：
program cx
implicit none
realt(63),p(73,144,63),r(73,144),relation, z(63)
character*12 c(63)
character*7 a(63)
integer i,j,k
open(5,file='D:\expriments\grads\slp.jan. grd',form='binary')
open(6,file='D:\expriments\grads\nino34. txt')
do i=1,63
read(6,'(A12)') c(i)
end do
do i=1,63
read(c(i),'(A7,F5.2)') a(i),t(i)
end do
do k=1,63
read(5) ((p(i,j,k),j=1,144),i=1,73)
end do
!print*,p
do i=1,73
do j=1,144
do k=1,63
z(k)=p(i,j,k)
end do
r(i,j)=relation(z,t,63)
end do
end do
!print*,r open(7,file='D:\expriments\grads\corr.slp. txt')
open(8,file='D:\expriments\grads\corr.slp. grd',form='binary')
write(7,'(144F5.2)')
((r(i,j),j=1,144),i=1,73)
write(8) ((r(i,j),j=1,144),i=1,73)
end
function relation(a,b,n)
implicit none
integer,intent(in)::n
real,intent(in)::a(n),b(n)
real::relation !返回的相关系数integer::i,j !循环控制变量
real::sfenzi,sfenmu1,sfenmu2,s !加法器real::amean,bmean !a,b向量的平均值
!计算平均值
s=0.
do i=1,n
s=s+a(i)
end do
amean=s/n
s=0.
do i=1,n
s=s+b(i)
end do
bmean=s/n !计算相关系数sfenzi=0.
sfenmu1=0.
sfenmu2=0.
do i=1,n
sfenzi=sfenzi+(a(i)-amean)*(b(i)-bmean)
sfenmu1=sfenmu1+(a(i)-amean)**2 sfenmu2=sfenmu2+(b(i)-bmean)**2 end do
relation=sfenzi/sqrt(sfenmu1*sfenmu2)
end function relation
保存于corr.slp.grd 文件中。

书写corr.slp.grd 的数据描述文件1.ctl：
dset e:\grads\20150706qimo\test1\corr.slp.grd
title slp
undef -9.99e33
xdef 144 linear -180 2.5
ydef 73 linear -90 2.5
zdef 1 levels 1000
tdef 63 linear jun1951 1yr
vars 1
slp 0 0 sea level pressure
endvars
编写1.gs文件：
'reinit'
'open e:\grads\20150706qimo\test1\1.ctl'
'set lev 1000'
'set lon -180 180'
'set lat -90 90'
'set t 1'
'set gxout shaded'
'set grads off'
'set black -0.2 0.2'
'd slp'
'set gxout contour'
'd slp'
'draw title Lu Yu '
'cbar 1.0 0'
'printim e:\grads\20150706qimo\test1\1.png white'
;
绘制1951-2013年1月Nino34海温指数与海平面气压相关系数图：
1951-2013年1月Nino34海温指数与海平面气压相关系数图
2 绘制1月份Nino34海温指数与7月份中国160站降水相关图：
2.1资料说明
➢Nino34指数：Nino34.txt是1951年1月至2013年12月Nino3.4区（热带太平洋：西经170度~西经120度、北纬5度~南纬5度）区域平均海温指数资料，共有63行数据，每行数据第一个数字为年份，后面12个数字为该年1-12月的海温指数；
➢160站降水资料：r1607.txt是中国160站1951~2013年共63年7月降水资料，数据是1951年7月各站降水量按照id.txt中各站顺序写入，该年160站都写入后再写下1952年7月各站降水量，...................，直至2013年7月各站降水量按照id.txt中各站顺序写入；
➢160站站点资料：id.txt是160站站号和纬度、经度资料，共有160行。

每行第一个数据为8位字符型的站号，后面两个数字分别为该站点纬度和经度；
2.2 步骤：
用Fortran编写corr.station.f90文件，计算1951-2013年1月Nino34海温指数与7月160站降水相关系数，计算结果保存于corr.7.txt和corr.7.grd文件中；
program cx
implicit none
integer p(160,63),i,j,k
real
s(63),r(160),relation,t(63),lon(160),lat(16 0),tim
character*12 c(63)
character*7 a(63)
character*8 stid(160)
integer nlev,nflag
open(6,file='E:\grads\20150706qimo\test2 \nino34.txt')
open(5,file='E:\grads\20150706qimo\test2 \r1607.txt')
open(10,file=E:\grads\20150706qimo\test 1\id.txt') do i=1,160
read(10,*) stid(i),lat(i),lon(i) end do
print*,stid
do j=1,63
read (5,*) (p(i,j),i=1,160) end do
do i=1,63
read(6,'(A12)') c(i)
end do
do i=1,63
read(c(i),'(A7,F5.2)') a(i),t(i) end do
do i=1,160
do j=1,63
s(j)=p(i,j)
end do
r(i)=relation(s,t,63)
end do
!print*,r
tim=0
nlev=1
nflag=1
open(7,file='E:\grads\20150706qimo\test2 \corr.7.txt')
open(8,file='E:\grads\20150706qimo\test1 \corr.7.grd',form='binary')
do i=1,160
write(7,*)
stid(i),lon(i),lat(i),tim,nlev,nflag,r(i) end do
do i=1,160
write(8)stid(i),lat(i),lon(i),tim,nlev,nflag, r(i)
end do
nlev=0
write(8)stid(160-1),lat(160-1),lon(160-1),tim,nlev,nflag,r(160-1)
end
function relation(a,b,n)
!本程序计算两列向量的相关系数
!a,b分别是待计算的向量
!n是向量的长度，要求两列向量等长implicit none
integer,intent(in)::n
real,intent(in)::a(n),b(n) real::relation !返回的相关系数integer::i,j !循环控制变量
real::sfenzi,sfenmu1,sfenmu2,s !加法器real::amean,bmean !a,b向量的平均值!计算平均值
s=0.
do i=1,n
s=s+a(i)
end do
amean=s/n
s=0.
do i=1,n
s=s+b(i)
end do
bmean=s/n
!计算相关系数
sfenzi=0.
sfenmu1=0.
sfenmu2=0.
do i=1,n
sfenzi=sfenzi+(a(i)-amean)*(b(i)-bmean)
sfenmu1=sfenmu1+(a(i)-amean)**2
sfenmu2=sfenmu2+(b(i)-bmean)**2 end do
relation=sfenzi/sqrt(sfenmu1*sfenmu2) end function relation
书写corr.7.grd 的数据描述文件2.ctl ；
dset E:\grads\20150706qimo\test2\corr.7.dat dtype station
stnmap E:\grads\20150706qimo\test2\corr7.map UNDEF 2.5E33
TITLE corr
TDEF 1 linear jan1951 1mo
VARS 1
r 0 99 relation
endvars
准备corr.7.grd的站点映射文件corr7.map；
准备GrADS画站点图时插值时所需格点数据Grid.grd和及其描述文件Grid161.ctl（已提供）；
编写2.gs文件，绘制1951-2013年1月Nino34海温指数与7月中国160站降水相关系数图，保存于corr.7.gmf 文件中。

'reinit' ‘open E:\grads\20150706qimo\test2\grid161.ctl' 'open E:\grads\20150706qimo\test2\2.ctl'
'set lat 3 54'
'set lon 70 135'
'define a=oacres(g,r.2)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'set grads off'
'set mpdsetworld'
'set map 1 1 1'
'enable print E:\grads\20150706qimo\test1\corr7.gmf' 'set gxout shaded' 'set cmin 0'
'd aa'
'set gxout contour' 'set ccolor 2'
'd aa'
'draw title Lu Yu ' 'print'
'disable print'
1951-2013年1月Nino34海温指数与7月我国160站降水相关系数图
分析：
图1 是1951-2013年1月Nino34海温指数与海平面气压相关系数图，从图中可以看出，赤道中东太平洋地区海温指数和海平面气压成明显正相关，最大相关系数已超过0.8；而赤道锡太平洋地区海温指数和海平面气压则成负相关，最小相关系数低于-0.6。

由此可以解释当发生厄尔尼诺事件时，由于赤道中东太平洋和秘鲁沿岸等地区海温上升，使沃克环流减弱，对中东太平洋地区的下沉气流产生负的影响，使中东太平洋地区海平面气压降低；与之相反，西太平洋地区上升气流减弱，海平面气压升高。

图2为1951-2013年1月Nino34海温指数与7月我国160站降水相关系数图，可以看出在我国西北部、中部、东北部分地区以及长江流域的降水与海温指数成正相关，而在青海，西藏广州等地则与海温指数成负相关。

由此可见，当热带太平洋出现暖水位相（即厄尔尼诺事件）时, 我国东北南部、华北大部以及长江中下游地区往往会有降水较少；而当热带太平洋出现冷水位相（即拉尼娜事件）时，该地区则往往会有较充沛的降水。

当发生厄尔尼诺时，哈得莱环流增强，因此与哈得莱环流下沉支密切相关的副热带高压也会增强，使得我国夏季降水受到影响。