MODULE stats
        USE nrtype
        IMPLICIT NONE

        INTERFACE avg_std_
                MODULE PROCEDURE avg_std_mask,avg_std
        END INTERFACE

        CONTAINS

        SUBROUTINE avg_std_mask(dat,mask,mean,var,sdev)
        USE nrtype
        IMPLICIT NONE

        REAL(DP), DIMENSION(:), INTENT(IN) :: dat
        REAL(DP), INTENT(out) :: mean, var, sdev
        INTEGER(I2B) :: n
        REAL(DP), DIMENSION(size(dat)) :: p,s
        REAL(DP) :: ep
        LOGICAL, DIMENSION(size(dat)) :: mask

        n=INT(COUNT(mask),I2B)
        mean=sum(dat,mask)/n
        s=dat-mean
        ep=sum(s,mask)
        p=s*s
        var=sum(p,mask)
        var=(var-ep**2/n)/(n-1)
        sdev=sqrt(var)

        END SUBROUTINE avg_std_mask

        SUBROUTINE avg_std(dat,mean,var,sdev)
        USE nrtype
        IMPLICIT NONE

        REAL(DP), DIMENSION(:), INTENT(IN) :: dat
        REAL(DP), INTENT(out) :: mean, var, sdev
        INTEGER(I2B) :: n
        REAL(DP), DIMENSION(size(dat)) :: p,s
        REAL(DP) :: ep

        n=INT(SIZE(dat,1),I2B)
        mean=sum(dat(:))/n
        s(:)=dat(:)-mean
        ep=sum(s(:))
        p(:)=s(:)*s(:)
        var=sum(p(:))
        var=(var-ep**2/n)/(n-1)
        sdev=sqrt(var)

        END SUBROUTINE avg_std


END MODULE stats