Calculate the Pearson correlation coefficient of two arrays

pearsonr.f90

real(kind=8) function pearsonr(x, y) result(r)
    ! given two arrays x and y, this function computes their Pearson correlation coefficient r
    implicit none
    real(kind=8), dimension(:) :: x, y
    real(kind=8), dimension(size(x)) :: xt, yt
    real(kind=8) :: ax, ay, df, sxx, sxy, syy
    integer :: n
    
    if (size(x) /= size(y)) STOP 'Dimension mismatch in pearsonr'
    
    n = size(x)
    ! find the means and subtract them
    ax = sum(x)/n
    ay = sum(y)/n
    xt = x - ax
    yt = y - ay

    sxx = dot_product(xt,xt)
    syy = dot_product(yt,yt)
    sxy = dot_product(xt,yt)
    r = sxy / sqrt(sxx*syy)
    ! abs(r) cannot be > 1, except for artifacts of floating point arithmetic
    r = max(min(r, 1.d0), -1.d0)

end function pearsonr