calc_covariance Function

    function calc_covariance(this, that)
        ! Covariance between two arrays
        ! General routine, no global data

        real(dp), intent(in) :: this(:), that(:)
        integer :: length1, length2
        integer :: i
        real(dp) :: sx, sy, sxy ! sums
        real(dp) :: meanx, meany
        real(dp) :: calc_covariance
        character(len=*), parameter :: t_r = 'calc_covariance'

        length1 = size(this)
        length2 = size(that)
        if (length1 /= length2) call stop_all(t_r, "ERROR: Something has gone very wrong indeed...")
        sx = 0.0_dp
        sy = 0.0_dp
        do i = 1, length1
            sx = sx + this(i)
            sy = sy + that(i)
            if (Errordebug > 0) write(stdout, *) this(i), that(i)
        end do
        meanx = sx / length1
        meany = sy / length1
        sxy = 0.0_dp
        do i = 1, length1
            sxy = sxy + (this(i) - meanx) * (that(i) - meany)
        end do
        calc_covariance = sxy / (length1)

    end function calc_covariance