program corr_calculation
implicit none
integer, parameter :: nlon=144, nlat=73, nyear=60
real :: nino34(nyear), precip(nlon,nlat,nyear), corr(nlon,nlat)
real :: sumx, sumy, sumxy, sumx2, sumy2, xmean, ymean
integer :: i, j, k, iyear, ios
character(80) :: filename
! 读取Nino3.4数据 (1951-2010年1月指数)
open(10, file='Nino34.txt', status='old')
do i = 1, 63
read(10, *, iostat=ios) iyear
if (iyear < 1951) cycle
if (iyear > 2010) exit
read(10, *) nino34(iyear-1950) ! 存储1月数据
do k = 1, 11 ! 跳过其他月份
read(10, *)
end do
end do
close(10)
! 读取降水数据 (1951-2010年7月)
open(20, file='pre7.grd', form='unformatted', access='direct', recl=nlon*nlat*4)
do k = 1, nyear
read(20, rec=k) ((precip(i,j,k), i=1,nlon), j=1,nlat)
end do
close(20)
! 计算相关系数
do i = 1, nlon
do j = 1, nlat
sumx = 0.0; sumy = 0.0; sumxy = 0.0; sumx2 = 0.0; sumy2 = 0.0
do k = 1, nyear
if (precip(i,j,k) < -9000) cycle ! 跳过缺测值
sumx = sumx + nino34(k)
sumy = sumy + precip(i,j,k)
sumxy = sumxy + nino34(k) * precip(i,j,k)
sumx2 = sumx2 + nino34(k)**2
sumy2 = sumy2 + precip(i,j,k)**2
end do
corr(i,j) = (nyear*sumxy - sumx*sumy) / &
sqrt( (nyear*sumx2 - sumx**2) * (nyear*sumy2 - sumy**2) )
end do
end do
! 输出结果
open(30, file='corr.txt')
do j = 1, nlat
write(30, '(144f8.3)') (corr(i,j), i=1,nlon)
end do
close(30)
open(40, file='corr.grd', form='unformatted', access='direct', recl=nlon*nlat*4)
write(40, rec=1) corr
close(40)
end program corr_calculation