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