在气象研究中，有时需要计算两个物理量时间序列的相关系数，例如塔西提站海平面气压与全球海平面气压的相关系数。这里，以FORTRAN语言和绘图软件GRADS为例进行说明。
一、计算相关系数的子程序
(一)相关系数的计算公式
这里给出的是皮尔逊相关系数的计算公式：
(二)程子序
!求两个一维时间序列的相关系数子程序
!n为时间长度，x,y为两个时间序列，r为相关系数
!auther: Yue Shen 2020/12/20
subroutine cor(n,x,y,r)
real x(n),y(n)
real avex,avey,sumx,sumy,varx,vary,covxy
sumx=0.0
sumy=0.0
varx=0.0
vary=0.0
covxy=0.0
!分别计算x和y两个序列的总和与平均值
do i=1,n
sumx=sumx+x(i)
sumy=sumy+y(i)
end do
avex=sumx/n
avey=sumy/n
!分别计算x，y序列的方差
do i=1,n
varx=varx+(x(i)-avex)**2
vary=vary+(y(i)-avey)**2
end do
varx=varx/n
vary=vary/n
!计算x和y的协方差
do i=1,n
covxy=covxy+(x(i)-avex)*(y(i)-avey)
end do
covxy=covxy/n
!计算相关系数
r=covxy/(sqrt(varx)*sqrt(vary))
end subroutine cor
(3)其他的相关系数
在一些文章的研究中，我们还会碰到其他的相关系数，例如偏相关系数、滑动相关系数等。这里就不在作详细介绍。感兴趣的读者可搜索相关文章进行了解，如：胡景高，周兵，陶丽.南亚高压特征参数与我国夏季降水的关系分析[J].气象，2010，36(4)：51-56. 一文的研究中就用到了偏相关系数和滑动相关系数。
二、塔西提站海平面气压与全球海平面气压的相关系数图的绘制
(一)资料
NCAR/NCEP月平均海平面气压再分析资料1979-2019
塔西提站海平面气压1951-2019
所用均可在NCAR/NCEP网站下载得到。
(二)数据处理程序
!auther: Yue Shen 2020/12/20
program so
real mslp(144,73,1979:2019,12),tahiti(1951:2019,12),r(144,73,12)
!read data mslp.grd and tahiti_slp.txt
open(1,file='E:\myresearch\so\mslp.grd',form='binary')
do iy=1979,2019
do im=1,12
do j=1,73
do i=1,144
read(1) mslp(i,j,iy,im)
end do
end do
end do
end do
close(1)
open(1,file='E:\myresearch\so\tahiti_slp.txt')
read(1,*)
read(1,*)
read(1,*)
do iy=1951,2019
read(1,'(i4,12(f6.1))') iyear,(tahiti(iy,im),im=1,12)
end do
close(1)
!caculate the correlation between mslp and tahiti_mslp 1979-2019
do im=1,12
do j=1,73
do i=1,144
call cor(41,mslp(i,j,1979:2019,im),tahiti(1979:2019,im),r(i,j,im))
end do
end do
end do
!save r to r,grd
open(1,file='E:\myresearch\so\r.grd',form='binary')
do im=1,12
do j=1,73
do i=1,144
write(1) r(i,j,im)
end do
end do
end do
close(1)
end program
(三)数据描述文件r.ctl以及绘图命令集dr.gs
1.r.ctl
dset E:\myresearch\so\mslp.grd
title Monthly NCEP/DOE Reanalysis 2
undef -9.96921e+36
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 linear 0 1
tdef 502 linear 00Z01JAN1979 1mo
vars 1
mslp 0 99 t,y,x Monthly Mean Sea Level Pressure
endvars
2.dr.gs
*auther: Yue Shen 2020/12/20
'reinit'
'open E:\myresearch\so\r.ctl'
i=1
while(i<=12)
'set t 'i''
'set grads off'
'set gxout contour'
'set cint 0.1'
'set black -0.3 0.3'
'd r'
'draw title correlation coefficient between slp of tahiti and other place 'i' '
'gxprint E:\myresearch\so\'i'.png png white'
'c'
i=i+1
endwhile
;
(四)部分结果
1月份
2月份
可以看到，在冬季塔西提站海平面气压与印度洋东部地区海平面气压存在明显的负相关关系。但是在某些月份(5月、6月、7月和10月)这种相关关系则不明显。
【说明】
(1)建立《气象程序集》的目的是为了减少大家在学习中编程时遇到的麻烦，将更多的时间用在思考问题上，毕竟个人认为编程绘图仅是探究问题或结果展示的一种途径(模式开发另说)。
(2)《气象程序集》中提供的所有代码有需要的均可复制使用，图像我想也没人会用。
(3)转载本文须注明出处。
(4)如提供的程序在使用中有任何问题欢迎私信交流。