## Fortran 90 SUBROUTINE UNBIASED_CORRELATION

```
SUBROUTINE UNBIASED_CORRELATION( N, A, B, msg, lag, r, t, m)

! Computes the correlation of two series A and B of length
! N as a function of lag. The correlation is unbiased
! because the sums in the covariance and standard devia-
! tions are divided by m, not m-1, where m is the number
! of overlapping grid points.

! The correlation r is defined as
!
!           cov(A,B)
!     r = -------------
!          s(A) * s(B)
!
! where cov() is covariance and s() is standard deviation.

! The series A and B must be prepared such that they
! are sent to this routine with zero lag. If necessary, this
! is accomplished by padding the beginning or ends (or both)
! of each time series with missing values (msg) so that their
! elements correspond one to one. The resultant time series
! will be of equal length N. Also, missing values (msg)
! distributed thoughout each time series is perfectly
! acceptable. Time series B will be shifted by the amount
! lag relative to time series A:

! RETURNS: r, the unbiased correlation, t, the significance
! of the unbiased correlation ( t is set to msg if B = A
! at lag 0, i.e. an autocorrelation at lag 0, or if the
! correlation is 1.0 in general), and m, the number of
! overlapping grid points.

! A:      t1  t2  t3  t4        ...       tN
! B:              t1  t2  t3  t4        ...       tN

!         \_ _/
!           V
!           lag = 2 (Defined to be > 0.)

IMPLICIT NONE

INTEGER, INTENT(IN)             :: N
REAL, DIMENSION(N), INTENT(IN)  :: A
REAL, DIMENSION(N), INTENT(IN)  :: B
REAL, INTENT(IN)                :: msg
INTEGER, INTENT(IN)             :: lag
REAL, INTENT(OUT)               :: r
REAL, INTENT(OUT)               :: t
INTEGER, INTENT(OUT)            :: m

REAL, DIMENSION(:), ALLOCATABLE :: x
REAL, DIMENSION(:), ALLOCATABLE :: y
REAL, DIMENSION(:), ALLOCATABLE :: xdev
REAL, DIMENSION(:), ALLOCATABLE :: ydev
REAL, DIMENSION(:), ALLOCATABLE :: xdevydev
REAL, DIMENSION(:), ALLOCATABLE :: xdevxdev
REAL, DIMENSION(:), ALLOCATABLE :: ydevydev

REAL                            :: xmn
REAL                            :: ymn
REAL                            :: COVxy
REAL                            :: Sx
REAL                            :: Sy

ALLOCATE( x(1:3*N) )
ALLOCATE( y(1:3*N) )

x(:) = msg
y(:) = msg

x(N+1:2*N) = A(:)
y(N+1:2*N) = B(:)

y = EOSHIFT( y, SHIFT = -lag, BOUNDARY = msg )

WHERE ( x == msg ) y = msg
WHERE ( y == msg ) x = msg

m = COUNT( x /= msg )

IF ( m < 3 ) THEN

WRITE (*,'(A40, I1)') "The number of overlapping points is m = ", m
WRITE (*,'(A35)')     "The value of m should be 3 or more."
WRITE (*,'(A17)')     "Decrease the lag."
WRITE (*,'(A52)')     "Execution halted in SUBROUTINE UNBIASED_CORRELATION."
STOP

END IF

xmn = SUM( x, DIM = 1, MASK = x /= msg ) / REAL(m)
ymn = SUM( y, DIM = 1, MASK = y /= msg ) / REAL(m)

ALLOCATE( xdev(1:3*N) )
ALLOCATE( ydev(1:3*N) )

xdev(:) = x(:) - xmn
WHERE ( x == msg ) xdev(:) = msg

ydev(:) = y(:) - ymn
WHERE ( y == msg ) ydev(:) = msg

ALLOCATE( xdevydev(1:3*N) )

xdevydev(:) = xdev(:) * ydev(:)
WHERE ( x == msg ) xdevydev(:) = msg

COVxy = SUM( xdevydev, DIM = 1, MASK = xdevydev /= msg ) / REAL(m)

ALLOCATE( xdevxdev(1:3*N) )

xdevxdev(:) = xdev(:) * xdev(:)
WHERE ( x == msg ) xdevxdev(:) = msg

Sx = SQRT( SUM( xdevxdev, DIM = 1, MASK = xdevxdev /= msg ) / REAL(m) )

ALLOCATE( ydevydev(1:3*N) )

ydevydev(:) = ydev(:) * ydev(:)
WHERE ( y == msg ) ydevydev(:) = msg

Sy = SQRT( SUM( ydevydev, DIM = 1, MASK = ydevydev /= msg ) / REAL(m) )

r = COVxy / ( Sx * Sy )

IF ( r /= 1.0 ) THEN

t = r * SQRT( (m - 2) / ( 1 - r*r ) )

ELSE

t = msg

END IF

DEALLOCATE( x, y, xdev, ydev, xdevydev, xdevxdev, ydevydev )

END SUBROUTINE UNBIASED_CORRELATION

```