在氣象研究中，有時需要計算兩個實體量時間序列的相關系數，例如塔西提站海平面氣壓與全球海平面氣壓的相關系數。這裡，以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)如提供的程式在使用中有任何問題歡迎私信交流。