编辑代码

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-20107月)
  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