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

实习要求： ?根据所提供的实习说明和资料，按要求完成实习内容；?将完成的实习内容按步骤写入实习报告中，包括资料说明、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',f orm='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 dorelation=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 20131301070' '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(160),ti m character*12 c(63) character*7 a(63) character*8 stid(160) integer nlev,nflag open(6,file='E:\grads\20150706qimo\test2\nin o34.txt') open(5,file='E:\grads\20150706qimo\test2\r16 07.txt') open(10,file=E:\grads\20150706qimo\test1\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\cor r.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画站点图时插值时所需格点数据Gri。