From f20a1e73135ccad0a3d3f774e053888631e1de4c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 7 Jan 2020 13:37:22 +0100 Subject: [PATCH] Fortran-Lib: remove NR references --- pre-proc/process_soilgrid/mo_kind.f90 | 128 +- pre-proc/process_soilgrid/mo_moment.f90 | 1223 +++++++++--------- src/lib/mo_boxcox.f90 | 857 ++----------- src/lib/mo_corr.f90 | 1564 +---------------------- src/lib/mo_kind.f90 | 32 +- src/lib/mo_moment.f90 | 94 +- 6 files changed, 892 insertions(+), 3006 deletions(-) mode change 100755 => 100644 pre-proc/process_soilgrid/mo_moment.f90 mode change 100644 => 100755 src/lib/mo_kind.f90 diff --git a/pre-proc/process_soilgrid/mo_kind.f90 b/pre-proc/process_soilgrid/mo_kind.f90 index 7a6cd616..f6e37651 100755 --- a/pre-proc/process_soilgrid/mo_kind.f90 +++ b/pre-proc/process_soilgrid/mo_kind.f90 @@ -1,51 +1,93 @@ -MODULE mo_kind +!> \file mo_kind.f90 + +!> \brief Define number representations + +!> \details This module declares the desired ranges and precisions of the number representations, +!> such as single precision or double precision, 32-bit or 46-bit integer, etc. +!> It confirms mostly with the nrtype module of Numerical Recipes in f90. + +!> \authors Juliane Mai, Matthias Cuntz, Nov 2011 +!> \date 2011-2014 +!> \copyright GNU Lesser General Public License http://www.gnu.org/licenses/ + +! Number model from which the SELECTED_REAL_KIND are requested: +! 4 byte REAL 8 byte REAL +! IEEE: precision = 6 precision = 15 +! exponent = 37 exponent = 307 +! CRAY: - precision = 13 +! exponent = 2465 + +! Written Juliane Mai, Matthias Cuntz, Nov 2011 +! Modified Matthias Cuntz, Nov 2011 - private/public +! - documentation +! - removed tab characters +! Matthias Cuntz, May 2014 - iso_fortran_env and iso_c_binding +! S. Mueller, Dec 2019 - remove NR specific types - ! This module declares the desired ranges and precisions of the number representations, - ! such as single precision or double precision, 32-bit or 46-bit integer, etc. - ! It confirms mostly with the nrtype module of Numerical Recipes in f90. +! License +! ------- +! This file is part of the UFZ Fortran library. - ! Number model from which the SELECTED_REAL_KIND are requested: - ! 4 byte REAL 8 byte REAL - ! IEEE: precision = 6 precision = 15 - ! exponent = 37 exponent = 307 - ! CRAY: - precision = 13 - ! exponent = 2465 +! The UFZ Fortran library is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. + +! The UFZ Fortran library is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. + +! You should have received a copy of the GNU Lesser General Public License +! along with the UFZ Fortran library (cf. gpl.txt and lgpl.txt). +! If not, see . + +! Copyright 2011-2014 Matthias Cuntz, Juliane Mai + +MODULE mo_kind - ! Written Juliane Mai, Nov 2011 - ! Modified Matthias Cuntz, Nov 2011 - private/public - ! - documentation - ! - removed tab characters + ! Does not work with compilers intel v11 and sun v12.2 + ! use, intrinsic :: iso_fortran_env, only: & + ! int8!, int16, int32, int64, real32, real64 + use, intrinsic :: iso_c_binding, only: & + c_short, c_int, c_long_long, c_float, c_double, c_float_complex, c_double_complex, c_bool IMPLICIT NONE - PRIVATE - - PUBLIC :: i1, i2, i4, i8, sp, dp, spc, dpc, lgt, sprs2_sp, sprs2_dp - - INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND(18) ! 8 Byte Integer - INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND(9) ! 4 Byte Integer - INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND(4) ! 2 Byte Integer - INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND(2) ! 1 Byte Integer - INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6,37) ! Single Precision Real - INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,307) ! Double Precision Real - ! Note: Intrinsic KIND is an extension to the Fortran 90 standard. - ! Matthias C thinks that spc=sp and dpc=dp would be enough. - INTEGER, PARAMETER :: spc = KIND((1.0_sp,1.0_sp)) ! Single Precision Complex - INTEGER, PARAMETER :: dpc = KIND((1.0_dp,1.0_dp)) ! Double Precision Complex - INTEGER, PARAMETER :: lgt = KIND(.true.) ! Logical - - ! Numerical recipes types for sparse arrays - TYPE sprs2_sp - INTEGER(I4) :: n, len - REAL(SP), DIMENSION(:), POINTER :: val - INTEGER(I4), DIMENSION(:), POINTER :: irow - INTEGER(I4), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_sp - TYPE sprs2_dp - INTEGER(I4) :: n, len - REAL(DP), DIMENSION(:), POINTER :: val - INTEGER(I4), DIMENSION(:), POINTER :: irow - INTEGER(I4), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_dp + !> 1 Byte Integer Kind + INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND(2) + ! INTEGER, PARAMETER :: i1 = int8 ! c_word does not exist; should be c_bool, probably; but see above + !> 2 Byte Integer Kind + ! INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND(4) + ! INTEGER, PARAMETER :: i2 = int16 + INTEGER, PARAMETER :: i2 = c_short + !> 4 Byte Integer Kind + ! INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND(9) + ! INTEGER, PARAMETER :: i4 = int32 + INTEGER, PARAMETER :: i4 = c_int + !> 8 Byte Integer Kind + ! INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND(18) + ! INTEGER, PARAMETER :: i8 = int64 + INTEGER, PARAMETER :: i8 = c_long_long + !> Single Precision Real Kind + ! INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6,37) + ! INTEGER, PARAMETER :: sp = real32 + INTEGER, PARAMETER :: sp = c_float + !> Double Precision Real Kind + ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,307) + ! INTEGER, PARAMETER :: dp = real64 + INTEGER, PARAMETER :: dp = c_double + !> Single Precision Complex Kind + ! INTEGER, PARAMETER :: spc = KIND((1.0_sp,1.0_sp)) + ! INTEGER, PARAMETER :: spc = sp + INTEGER, PARAMETER :: spc = c_float_complex + !> Double Precision Complex Kind + ! INTEGER, PARAMETER :: dpc = KIND((1.0_dp,1.0_dp)) + ! INTEGER, PARAMETER :: dpc = dp + INTEGER, PARAMETER :: dpc = c_double_complex + !> Logical Kind + INTEGER, PARAMETER :: lgt = KIND(.true.) + !INTEGER, PARAMETER :: lgt = c_bool + ! Types have to be in a public section for doxygen END MODULE mo_kind diff --git a/pre-proc/process_soilgrid/mo_moment.f90 b/pre-proc/process_soilgrid/mo_moment.f90 old mode 100755 new mode 100644 index 6e6fc98d..50bdcac3 --- a/pre-proc/process_soilgrid/mo_moment.f90 +++ b/pre-proc/process_soilgrid/mo_moment.f90 @@ -7,16 +7,14 @@ MODULE mo_moment ! i.e. they are normally divided by n and not (n-1) ! Literature - ! Moments (incl. mean, stddev, etc.) - Chapter 14 of - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 ! Central moments and error variances ! LH Benedict & RD Gould, Towards better uncertainty estimates for turbulence statistics, ! Experiments in Fluids 22, 129-136, 1996 ! Written Nov 2011, Matthias Cuntz ! Modified, MC, Dec 2011 - mod. correlation, covariance + ! Modified by M. Schroen, Sep 2015, average/mean for single value + ! Modified by S. Mueller, Dec 2019, no citations for common sence ! License ! ------- @@ -38,7 +36,7 @@ MODULE mo_moment ! Copyright 2011-2012 Matthias Cuntz - USE mo_kind, ONLY: i4, sp, dp + USE mo_kind, ONLY : i4, sp, dp IMPLICIT NONE @@ -99,15 +97,10 @@ MODULE mo_moment ! m = absdev(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE absdev - MODULE PROCEDURE absdev_sp, absdev_dp + MODULE PROCEDURE absdev_sp, absdev_dp END INTERFACE absdev ! ------------------------------------------------------------------ @@ -152,15 +145,10 @@ MODULE mo_moment ! m = average(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE average - MODULE PROCEDURE average_sp, average_dp + MODULE PROCEDURE average_sp, average_dp END INTERFACE average ! ------------------------------------------------------------------ @@ -214,7 +202,7 @@ MODULE mo_moment ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE central_moment - MODULE PROCEDURE central_moment_sp, central_moment_dp + MODULE PROCEDURE central_moment_sp, central_moment_dp END INTERFACE central_moment ! ------------------------------------------------------------------ @@ -268,7 +256,7 @@ MODULE mo_moment ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE central_moment_var - MODULE PROCEDURE central_moment_var_sp, central_moment_var_dp + MODULE PROCEDURE central_moment_var_sp, central_moment_var_dp END INTERFACE central_moment_var ! ------------------------------------------------------------------ @@ -315,16 +303,11 @@ MODULE mo_moment ! m = correlation(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.))) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 ! Modified, MC, Dec 2011 - covariance as <(x-)(y-)> instead of - INTERFACE correlation - MODULE PROCEDURE correlation_sp, correlation_dp + MODULE PROCEDURE correlation_sp, correlation_dp END INTERFACE correlation ! ------------------------------------------------------------------ @@ -371,16 +354,11 @@ MODULE mo_moment ! m = covariance(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.))) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 ! Modified, MC, Dec 2011 - covariance as <(x-)(y-)> instead of - INTERFACE covariance - MODULE PROCEDURE covariance_sp, covariance_dp + MODULE PROCEDURE covariance_sp, covariance_dp END INTERFACE covariance ! ------------------------------------------------------------------ @@ -426,15 +404,10 @@ MODULE mo_moment ! m = kurtosis(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE kurtosis - MODULE PROCEDURE kurtosis_sp, kurtosis_dp + MODULE PROCEDURE kurtosis_sp, kurtosis_dp END INTERFACE kurtosis ! ------------------------------------------------------------------ @@ -479,15 +452,10 @@ MODULE mo_moment ! m = mean(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE mean - MODULE PROCEDURE mean_sp, mean_dp + MODULE PROCEDURE mean_sp, mean_dp END INTERFACE mean ! ------------------------------------------------------------------ @@ -544,7 +512,7 @@ MODULE mo_moment ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE mixed_central_moment - MODULE PROCEDURE mixed_central_moment_sp, mixed_central_moment_dp + MODULE PROCEDURE mixed_central_moment_sp, mixed_central_moment_dp END INTERFACE mixed_central_moment ! ------------------------------------------------------------------ @@ -602,7 +570,7 @@ MODULE mo_moment ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE mixed_central_moment_var - MODULE PROCEDURE mixed_central_moment_var_sp, mixed_central_moment_var_dp + MODULE PROCEDURE mixed_central_moment_var_sp, mixed_central_moment_var_dp END INTERFACE mixed_central_moment_var ! ------------------------------------------------------------------ @@ -611,10 +579,10 @@ MODULE mo_moment ! moment ! PURPOSE - ! Calculates the first four sample moments of a vector, i.e. mean, variance, skewness, and kurtosis, + ! Calculates the first four sample/population moments of a vector, i.e. mean, variance, skewness, and kurtosis, ! as well as standard deviation and mean absolute devation. ! mean = sum(x)/n - ! variance = sum((x-mean(x))**2)/(n-1) + ! variance = sum((x-mean(x))**2)/(n-1) or sum((x-mean(x))**2)/n if sample is False ! skewness = sum(((x-mean(x))/stddev(x))**3)/n ! kurtosis = sum(((x-mean(x))/stddev(x))**4)/n - 3 ! stddev = sqrt(variance) @@ -626,7 +594,7 @@ MODULE mo_moment ! CALLING SEQUENCE ! call moment(dat, average=average, variance=variance, skewness=skewness, kurtosis=kurtosis, & - ! mean=mean, stddev=stddev, absdev=absdev, mask=mask) + ! mean=mean, stddev=stddev, absdev=absdev, mask=mask, sample=sample) ! INTENT(IN) ! real(sp/dp) :: dat(:) 1D-array with input numbers @@ -640,17 +608,19 @@ MODULE mo_moment ! INTENT(IN), OPTIONAL ! logical :: mask(:) 1D-array of logical values with size(dat). ! If present, only those locations in vec corresponding to the true values in mask are used. + ! logical :: sample logical value + ! If present and False, the population variance and std-dev will be calculated (divide by n). ! INTENT(INOUT), OPTIONAL ! None ! INTENT(OUT), OPTIONAL ! real(sp/dp) :: average average of input vector - ! real(sp/dp) :: variance sample variance of input vector + ! real(sp/dp) :: variance sample variance of input vector (either a sample or pupulation moment) ! real(sp/dp) :: skewness skewness of input vector ! real(sp/dp) :: kurtosis excess kurtosis of input vector ! real(sp/dp) :: mean same as average - ! real(sp/dp) :: stddev sqaure root of variance + ! real(sp/dp) :: stddev sqaure root of variance (either a sample or pupulation moment) ! real(sp/dp) :: absdev mean absolute deviations from average ! RESTRICTIONS @@ -661,15 +631,11 @@ MODULE mo_moment ! call moment(vec, mask=(vec >= 0.), mean=m, stddev=s) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 + ! Modified, S. Mueller, Dec 2019 added optional sample input-para to switch sample to population variance and std-dev. INTERFACE moment - MODULE PROCEDURE moment_sp, moment_dp + MODULE PROCEDURE moment_sp, moment_dp END INTERFACE moment ! ------------------------------------------------------------------ @@ -714,15 +680,10 @@ MODULE mo_moment ! m = skewness(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE skewness - MODULE PROCEDURE skewness_sp, skewness_dp + MODULE PROCEDURE skewness_sp, skewness_dp END INTERFACE skewness ! ------------------------------------------------------------------ @@ -769,15 +730,10 @@ MODULE mo_moment ! m = stddev(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE stddev - MODULE PROCEDURE stddev_sp, stddev_dp + MODULE PROCEDURE stddev_sp, stddev_dp END INTERFACE stddev ! ------------------------------------------------------------------ @@ -824,15 +780,10 @@ MODULE mo_moment ! m = variance(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE variance - MODULE PROCEDURE variance_sp, variance_dp + MODULE PROCEDURE variance_sp, variance_dp END INTERFACE variance ! ------------------------------------------------------------------ @@ -849,9 +800,9 @@ FUNCTION absdev_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: absdev_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: absdev_dp REAL(dp) :: n @@ -859,19 +810,19 @@ FUNCTION absdev_dp(dat, mask) LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error absdev_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error absdev_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'absdev_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'absdev_dp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n + ave = sum(dat(:), mask = maske) / n ! Sum of absolute deviation - absdev_dp = sum(abs(dat(:)-ave), mask=maske)/n + absdev_dp = sum(abs(dat(:) - ave), mask = maske) / n END FUNCTION absdev_dp @@ -880,9 +831,9 @@ FUNCTION absdev_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: absdev_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: absdev_sp REAL(sp) :: n @@ -890,19 +841,19 @@ FUNCTION absdev_sp(dat, mask) LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error absdev_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error absdev_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'absdev_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'absdev_sp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n + ave = sum(dat(:), mask = maske) / n ! Sum of absolute deviation - absdev_sp = sum(abs(dat(:)-ave), mask=maske)/n + absdev_sp = sum(abs(dat(:) - ave), mask = maske) / n END FUNCTION absdev_sp @@ -912,25 +863,24 @@ FUNCTION average_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: average_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: average_dp REAL(dp) :: n LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error average_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error average_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'average_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if ! Average - average_dp = sum(dat(:), mask=maske)/n + average_dp = sum(dat(:), mask = maske) / n END FUNCTION average_dp @@ -939,25 +889,24 @@ FUNCTION average_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: average_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: average_sp REAL(sp) :: n LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error average_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error average_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'average_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if ! Average - average_sp = sum(dat(:), mask=maske)/n + average_sp = sum(dat(:), mask = maske) / n END FUNCTION average_sp @@ -967,37 +916,37 @@ FUNCTION central_moment_dp(x, r, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: r - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: central_moment_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + INTEGER(i4), INTENT(IN) :: r + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: central_moment_dp - REAL(dp) :: n, mx + REAL(dp) :: n, mx LOGICAL, DIMENSION(size(x)) :: maske if (r .lt. 0) then - central_moment_dp = 0.0_dp - return - endif + central_moment_dp = 0.0_dp + return + end if if (r .eq. 0) then - central_moment_dp = 1.0_dp - return - endif + central_moment_dp = 1.0_dp + return + end if if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error central_moment_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error central_moment_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (2.0_dp+tiny(2.0_dp))) stop 'central_moment_dp: n must be at least 3' + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'central_moment_dp: n must be at least 3' ! average - mx = sum(x, mask=maske) / n + mx = sum(x, mask = maske) / n ! central moment - central_moment_dp = sum((x-mx)**r, mask=maske) / n + central_moment_dp = sum((x - mx)**r, mask = maske) / n END FUNCTION central_moment_dp @@ -1006,37 +955,37 @@ FUNCTION central_moment_sp(x, r, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: r - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: central_moment_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + INTEGER(i4), INTENT(IN) :: r + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: central_moment_sp REAL(sp) :: n, mx LOGICAL, DIMENSION(size(x)) :: maske if (r .lt. 0) then - central_moment_sp = 0.0_sp - return - endif + central_moment_sp = 0.0_sp + return + end if if (r .eq. 0) then - central_moment_sp = 1.0_sp - return - endif + central_moment_sp = 1.0_sp + return + end if if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error central_moment_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error central_moment_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (2.0_sp+tiny(2.0_sp))) stop 'central_moment_sp: n must be at least 3' + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'central_moment_sp: n must be at least 3' ! average - mx = sum(x, mask=maske) / n + mx = sum(x, mask = maske) / n ! central moment - central_moment_sp = sum((x-mx)**r, mask=maske) / n + central_moment_sp = sum((x - mx)**r, mask = maske) / n END FUNCTION central_moment_sp @@ -1046,36 +995,36 @@ FUNCTION central_moment_var_dp(x, r, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: r - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: central_moment_var_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + INTEGER(i4), INTENT(IN) :: r + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: central_moment_var_dp REAL(dp) :: n, rr, u2r, ur, urm1, urp1, u2 LOGICAL, DIMENSION(size(x)) :: maske if (r.le.1) then - central_moment_var_dp = 0.0_dp - return - endif + central_moment_var_dp = 0.0_dp + return + end if if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error central_moment_var_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error central_moment_var_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (2.0_dp+tiny(2.0_dp))) stop 'central_moment_var_dp: n must be at least 3' - - u2r = central_moment(x, 2*r, mask=maske) - ur = central_moment(x, r, mask=maske) - urm1 = central_moment(x, r-1, mask=maske) - u2 = central_moment(x, 2, mask=maske) - urp1 = central_moment(x, r+1, mask=maske) - rr = real(r,dp) - central_moment_var_dp = (u2r - ur*ur + rr*rr*urm1*urm1*u2 - 2.0_dp*rr*urp1*urm1) / n + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'central_moment_var_dp: n must be at least 3' + + u2r = central_moment(x, 2 * r, mask = maske) + ur = central_moment(x, r, mask = maske) + urm1 = central_moment(x, r - 1, mask = maske) + u2 = central_moment(x, 2, mask = maske) + urp1 = central_moment(x, r + 1, mask = maske) + rr = real(r, dp) + central_moment_var_dp = (u2r - ur * ur + rr * rr * urm1 * urm1 * u2 - 2.0_dp * rr * urp1 * urm1) / n END FUNCTION central_moment_var_dp @@ -1084,36 +1033,36 @@ FUNCTION central_moment_var_sp(x, r, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: r - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: central_moment_var_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + INTEGER(i4), INTENT(IN) :: r + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: central_moment_var_sp REAL(sp) :: n, rr, u2r, ur, urm1, urp1, u2 LOGICAL, DIMENSION(size(x)) :: maske if (r.le.1) then - central_moment_var_sp = 0.0_sp - return - endif + central_moment_var_sp = 0.0_sp + return + end if if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error central_moment_var_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error central_moment_var_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (2.0_sp+tiny(2.0_sp))) stop 'central_moment_var_sp: n must be at least 3' - - u2r = central_moment(x, 2*r, mask=maske) - ur = central_moment(x, r, mask=maske) - urm1 = central_moment(x, r-1, mask=maske) - u2 = central_moment(x, 2, mask=maske) - urp1 = central_moment(x, r+1, mask=maske) - rr = real(r,sp) - central_moment_var_sp = (u2r - ur*ur + rr*rr*urm1*urm1*u2 - 2.0_sp*rr*urp1*urm1) / n + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'central_moment_var_sp: n must be at least 3' + + u2r = central_moment(x, 2 * r, mask = maske) + ur = central_moment(x, r, mask = maske) + urm1 = central_moment(x, r - 1, mask = maske) + u2 = central_moment(x, 2, mask = maske) + urp1 = central_moment(x, r + 1, mask = maske) + rr = real(r, sp) + central_moment_var_sp = (u2r - ur * ur + rr * rr * urm1 * urm1 * u2 - 2.0_sp * rr * urp1 * urm1) / n END FUNCTION central_moment_var_sp @@ -1123,33 +1072,33 @@ FUNCTION correlation_dp(x, y, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: correlation_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + REAL(dp), DIMENSION(:), INTENT(IN) :: y + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: correlation_dp - REAL(dp) :: n - REAL(dp) :: mx, my - REAL(dp) :: sx, sy, covar + REAL(dp) :: n + REAL(dp) :: mx, my + REAL(dp) :: sx, sy, covar LOGICAL, DIMENSION(size(x)) :: maske if (size(x) .ne. size(y)) stop 'Error correlation_dp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error correlation_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error correlation_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'correlation_dp: n must be at least 2' + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'correlation_dp: n must be at least 2' ! Mean and Stddev of x and y - call moment(x, mx, stddev=sx, mask=maske) - call moment(y, my, stddev=sy, mask=maske) - covar = sum((x-mx)*(y-my), mask=maske) / n + call moment(x, mx, stddev = sx, mask = maske) + call moment(y, my, stddev = sy, mask = maske) + covar = sum((x - mx) * (y - my), mask = maske) / n ! correlation - correlation_dp = covar / (sx*sy) + correlation_dp = covar / (sx * sy) END FUNCTION correlation_dp @@ -1158,33 +1107,33 @@ FUNCTION correlation_sp(x, y, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: correlation_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + REAL(sp), DIMENSION(:), INTENT(IN) :: y + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: correlation_sp - REAL(sp) :: n - REAL(sp) :: mx, my - REAL(sp) :: sx, sy, covar + REAL(sp) :: n + REAL(sp) :: mx, my + REAL(sp) :: sx, sy, covar LOGICAL, DIMENSION(size(x)) :: maske if (size(x) .ne. size(y)) stop 'Error correlation_sp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error correlation_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error correlation_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'correlation_sp: n must be at least 2' + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'correlation_sp: n must be at least 2' ! Mean and Stddev of x and y - call moment(x, mx, stddev=sx, mask=maske) - call moment(y, my, stddev=sy, mask=maske) - covar = sum((x-mx)*(y-my), mask=maske) / n + call moment(x, mx, stddev = sx, mask = maske) + call moment(y, my, stddev = sy, mask = maske) + covar = sum((x - mx) * (y - my), mask = maske) / n ! correlation - correlation_sp = covar / (sx*sy) + correlation_sp = covar / (sx * sy) END FUNCTION correlation_sp @@ -1194,30 +1143,30 @@ FUNCTION covariance_dp(x, y, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: covariance_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + REAL(dp), DIMENSION(:), INTENT(IN) :: y + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: covariance_dp - REAL(dp) :: n - REAL(dp) :: mx, my + REAL(dp) :: n + REAL(dp) :: mx, my LOGICAL, DIMENSION(size(x)) :: maske if (size(x) .ne. size(y)) stop 'Error covariance_dp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error covariance_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error covariance_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'covariance_dp: n must be at least 2' + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'covariance_dp: n must be at least 2' ! Mean of x and y - mx = mean(x, mask=maske) - my = mean(y, mask=maske) - covariance_dp = sum((x-mx)*(y-my), mask=maske) / n + mx = mean(x, mask = maske) + my = mean(y, mask = maske) + covariance_dp = sum((x - mx) * (y - my), mask = maske) / n END FUNCTION covariance_dp @@ -1226,30 +1175,30 @@ FUNCTION covariance_sp(x, y, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: covariance_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + REAL(sp), DIMENSION(:), INTENT(IN) :: y + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: covariance_sp - REAL(sp) :: n - REAL(sp) :: mx, my + REAL(sp) :: n + REAL(sp) :: mx, my LOGICAL, DIMENSION(size(x)) :: maske if (size(x) .ne. size(y)) stop 'Error covariance_sp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error covariance_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error covariance_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'covariance_sp: n must be at least 2' + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'covariance_sp: n must be at least 2' ! Mean of x and y - mx = mean(x, mask=maske) - my = mean(y, mask=maske) - covariance_sp = sum((x-mx)*(y-my), mask=maske) / n + mx = mean(x, mask = maske) + my = mean(y, mask = maske) + covariance_sp = sum((x - mx) * (y - my), mask = maske) / n END FUNCTION covariance_sp @@ -1259,39 +1208,39 @@ FUNCTION kurtosis_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: kurtosis_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: kurtosis_dp REAL(dp) :: n REAL(dp) :: ep, ave, var REAL(dp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error kurtosis_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error kurtosis_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'kurtosis_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'kurtosis_dp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_dp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_dp) if (abs(var) .lt. tiny(0.0_dp)) stop 'kurtosis_dp: no kurtosis when zero variance' ! Kurtosis - p(:) = p(:)*s(:)*s(:) - kurtosis_dp = sum(p(:), mask=maske) - kurtosis_dp = kurtosis_dp/(n*var*var) - 3.0_dp + p(:) = p(:) * s(:) * s(:) + kurtosis_dp = sum(p(:), mask = maske) + kurtosis_dp = kurtosis_dp / (n * var * var) - 3.0_dp END FUNCTION kurtosis_dp @@ -1300,39 +1249,39 @@ FUNCTION kurtosis_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: kurtosis_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: kurtosis_sp REAL(sp) :: n REAL(sp) :: ep, ave, var REAL(sp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error kurtosis_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error kurtosis_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'kurtosis_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'kurtosis_sp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_sp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_sp) if (abs(var) .lt. tiny(0.0_sp)) stop 'kurtosis_sp: no kurtosis when zero variance' ! Kurtosis - p(:) = p(:)*s(:)*s(:) - kurtosis_sp = sum(p(:), mask=maske) - kurtosis_sp = kurtosis_sp/(n*var*var) - 3.0_sp + p(:) = p(:) * s(:) * s(:) + kurtosis_sp = sum(p(:), mask = maske) + kurtosis_sp = kurtosis_sp / (n * var * var) - 3.0_sp END FUNCTION kurtosis_sp @@ -1342,26 +1291,25 @@ FUNCTION mean_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: mean_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: mean_dp REAL(dp) :: n LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error mean_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error mean_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'mean_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if ! Mean - mean_dp = sum(dat(:), mask=maske)/n + mean_dp = sum(dat(:), mask = maske) / n END FUNCTION mean_dp @@ -1370,26 +1318,25 @@ FUNCTION mean_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: mean_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: mean_sp REAL(sp) :: n LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error mean_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error mean_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'mean_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if ! Mean - mean_sp = sum(dat(:), mask=maske)/n + mean_sp = sum(dat(:), mask = maske) / n END FUNCTION mean_sp @@ -1399,48 +1346,48 @@ FUNCTION mixed_central_moment_dp(x, y, r, s, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: r - INTEGER(i4), INTENT(IN) :: s - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: mixed_central_moment_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + REAL(dp), DIMENSION(:), INTENT(IN) :: y + INTEGER(i4), INTENT(IN) :: r + INTEGER(i4), INTENT(IN) :: s + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: mixed_central_moment_dp - REAL(dp) :: n, mx, my + REAL(dp) :: n, mx, my REAL(dp), DIMENSION(size(x)) :: xx, yy - LOGICAL, DIMENSION(size(x)) :: maske + LOGICAL, DIMENSION(size(x)) :: maske if (r.lt.0 .or. s.lt.0) then - mixed_central_moment_dp = 0.0_dp - return - endif + mixed_central_moment_dp = 0.0_dp + return + end if if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_dp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (2.0_dp+tiny(2.0_dp))) stop 'mixed_central_moment_dp: n must be at least 3' + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'mixed_central_moment_dp: n must be at least 3' ! Averages of x and y - mx = sum(x, mask=maske) / n - my = sum(y, mask=maske) / n + mx = sum(x, mask = maske) / n + my = sum(y, mask = maske) / n ! Mixed central moment if (r>0) then - xx = (x-mx)**r + xx = (x - mx)**r else - xx = 1._dp - endif + xx = 1._dp + end if if (s>0) then - yy = (y-my)**s + yy = (y - my)**s else - yy = 1._dp - endif - mixed_central_moment_dp = sum(xx*yy, mask=maske) / n + yy = 1._dp + end if + mixed_central_moment_dp = sum(xx * yy, mask = maske) / n END FUNCTION mixed_central_moment_dp @@ -1449,48 +1396,48 @@ FUNCTION mixed_central_moment_sp(x, y, r, s, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: r - INTEGER(i4), INTENT(IN) :: s - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: mixed_central_moment_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + REAL(sp), DIMENSION(:), INTENT(IN) :: y + INTEGER(i4), INTENT(IN) :: r + INTEGER(i4), INTENT(IN) :: s + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: mixed_central_moment_sp - REAL(sp) :: n, mx, my + REAL(sp) :: n, mx, my REAL(sp), DIMENSION(size(x)) :: xx, yy - LOGICAL, DIMENSION(size(x)) :: maske + LOGICAL, DIMENSION(size(x)) :: maske if (r.lt.0 .or. s.lt.0) then - mixed_central_moment_sp = 0.0_sp - return - endif + mixed_central_moment_sp = 0.0_sp + return + end if if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_sp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (2.0_sp+tiny(2.0_sp))) stop 'mixed_central_moment_sp: n must be at least 3' + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'mixed_central_moment_sp: n must be at least 3' ! Averages of x and y - mx = sum(x, mask=maske) / n - my = sum(y, mask=maske) / n + mx = sum(x, mask = maske) / n + my = sum(y, mask = maske) / n ! Mixed central moment if (r>0) then - xx = (x-mx)**r + xx = (x - mx)**r else - xx = 1._sp - endif + xx = 1._sp + end if if (s>0) then - yy = (y-my)**s + yy = (y - my)**s else - yy = 1._sp - endif - mixed_central_moment_sp = sum(xx*yy, mask=maske) / n + yy = 1._sp + end if + mixed_central_moment_sp = sum(xx * yy, mask = maske) / n END FUNCTION mixed_central_moment_sp @@ -1500,12 +1447,12 @@ FUNCTION mixed_central_moment_var_dp(x, y, r, s, mask) ! Error variance of mixed central moment (Benedict & Gould 1996) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: r - INTEGER(i4), INTENT(IN) :: s - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: mixed_central_moment_var_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: x + REAL(dp), DIMENSION(:), INTENT(IN) :: y + INTEGER(i4), INTENT(IN) :: r + INTEGER(i4), INTENT(IN) :: s + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: mixed_central_moment_var_dp REAL(dp) :: u2r2s, urs, urm1s, u20, urp1s, ursm1, u02, ursp1, u11 REAL(dp) :: n, rr, ss @@ -1513,31 +1460,31 @@ FUNCTION mixed_central_moment_var_dp(x, y, r, s, mask) if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_var_dp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_dp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_dp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(x),dp) - endif - if (n .le. (2.0_dp+tiny(2.0_dp))) stop 'mixed_central_moment_var_dp: n must be at least 3' - - u2r2s = mixed_central_moment(x, y, 2*r, 2*s, mask=maske) - urs = mixed_central_moment(x, y, r, s, mask=maske) - urm1s = mixed_central_moment(x, y, r-1, s, mask=maske) - u20 = mixed_central_moment(x, y, 2, 0, mask=maske) - urp1s = mixed_central_moment(x, y, r+1, s, mask=maske) - ursm1 = mixed_central_moment(x, y, r, s-1, mask=maske) - u02 = mixed_central_moment(x, y, 0, 2, mask=maske) - ursp1 = mixed_central_moment(x, y, r, s+1, mask=maske) - u11 = mixed_central_moment(x, y, 1, 1, mask=maske) - rr = real(r,dp) - ss = real(s,dp) - - mixed_central_moment_var_dp = (u2r2s - urs*urs & - + rr*rr**u20*urm1s*urm1s + ss*ss*u02*ursm1*ursm1 & - + 2.0_dp*rr*ss*u11*urm1s*ursm1 & - - 2.0_dp*rr*urp1s*urm1s - 2.0_dp*ss*ursp1*ursm1) / n + maske(:) = .true. + n = real(size(x), dp) + end if + if (n .le. (2.0_dp + tiny(2.0_dp))) stop 'mixed_central_moment_var_dp: n must be at least 3' + + u2r2s = mixed_central_moment(x, y, 2 * r, 2 * s, mask = maske) + urs = mixed_central_moment(x, y, r, s, mask = maske) + urm1s = mixed_central_moment(x, y, r - 1, s, mask = maske) + u20 = mixed_central_moment(x, y, 2, 0, mask = maske) + urp1s = mixed_central_moment(x, y, r + 1, s, mask = maske) + ursm1 = mixed_central_moment(x, y, r, s - 1, mask = maske) + u02 = mixed_central_moment(x, y, 0, 2, mask = maske) + ursp1 = mixed_central_moment(x, y, r, s + 1, mask = maske) + u11 = mixed_central_moment(x, y, 1, 1, mask = maske) + rr = real(r, dp) + ss = real(s, dp) + + mixed_central_moment_var_dp = (u2r2s - urs * urs & + + rr * rr**u20 * urm1s * urm1s + ss * ss * u02 * ursm1 * ursm1 & + + 2.0_dp * rr * ss * u11 * urm1s * ursm1 & + - 2.0_dp * rr * urp1s * urm1s - 2.0_dp * ss * ursp1 * ursm1) / n END FUNCTION mixed_central_moment_var_dp @@ -1546,12 +1493,12 @@ FUNCTION mixed_central_moment_var_sp(x, y, r, s, mask) ! Error variance of mixed central moment (Benedict & Gould 1996) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: r - INTEGER(i4), INTENT(IN) :: s - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: mixed_central_moment_var_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: x + REAL(sp), DIMENSION(:), INTENT(IN) :: y + INTEGER(i4), INTENT(IN) :: r + INTEGER(i4), INTENT(IN) :: s + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: mixed_central_moment_var_sp REAL(sp) :: u2r2s, urs, urm1s, u20, urp1s, ursm1, u02, ursp1, u11 REAL(sp) :: n, rr, ss @@ -1559,171 +1506,185 @@ FUNCTION mixed_central_moment_var_sp(x, y, r, s, mask) if (size(x) .ne. size(y)) stop 'Error mixed_central_moment_var_sp: size(x) .ne. size(y)' if (present(mask)) then - if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_sp: size(mask) .ne. size(x)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(x)) stop 'Error mixed_central_moment_var_sp: size(mask) .ne. size(x)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(x),sp) - endif - if (n .le. (2.0_sp+tiny(2.0_sp))) stop 'mixed_central_moment_var_sp: n must be at least 3' - - u2r2s = mixed_central_moment(x, y, 2*r, 2*s, mask=maske) - urs = mixed_central_moment(x, y, r, s, mask=maske) - urm1s = mixed_central_moment(x, y, r-1, s, mask=maske) - u20 = mixed_central_moment(x, y, 2, 0, mask=maske) - urp1s = mixed_central_moment(x, y, r+1, s, mask=maske) - ursm1 = mixed_central_moment(x, y, r, s-1, mask=maske) - u02 = mixed_central_moment(x, y, 0, 2, mask=maske) - ursp1 = mixed_central_moment(x, y, r, s+1, mask=maske) - u11 = mixed_central_moment(x, y, 1, 1, mask=maske) - rr = real(r,sp) - ss = real(s,sp) - - mixed_central_moment_var_sp = (u2r2s - urs*urs & - + rr*rr**u20*urm1s*urm1s + ss*ss*u02*ursm1*ursm1 & - + 2.0_sp*rr*ss*u11*urm1s*ursm1 & - - 2.0_sp*rr*urp1s*urm1s - 2.0_sp*ss*ursp1*ursm1) / n + maske(:) = .true. + n = real(size(x), sp) + end if + if (n .le. (2.0_sp + tiny(2.0_sp))) stop 'mixed_central_moment_var_sp: n must be at least 3' + + u2r2s = mixed_central_moment(x, y, 2 * r, 2 * s, mask = maske) + urs = mixed_central_moment(x, y, r, s, mask = maske) + urm1s = mixed_central_moment(x, y, r - 1, s, mask = maske) + u20 = mixed_central_moment(x, y, 2, 0, mask = maske) + urp1s = mixed_central_moment(x, y, r + 1, s, mask = maske) + ursm1 = mixed_central_moment(x, y, r, s - 1, mask = maske) + u02 = mixed_central_moment(x, y, 0, 2, mask = maske) + ursp1 = mixed_central_moment(x, y, r, s + 1, mask = maske) + u11 = mixed_central_moment(x, y, 1, 1, mask = maske) + rr = real(r, sp) + ss = real(s, sp) + + mixed_central_moment_var_sp = (u2r2s - urs * urs & + + rr * rr**u20 * urm1s * urm1s + ss * ss * u02 * ursm1 * ursm1 & + + 2.0_sp * rr * ss * u11 * urm1s * ursm1 & + - 2.0_sp * rr * urp1s * urm1s - 2.0_sp * ss * ursp1 * ursm1) / n END FUNCTION mixed_central_moment_var_sp ! ------------------------------------------------------------------ - SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask) + SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - REAL(dp), OPTIONAL, INTENT(OUT) :: average - REAL(dp), OPTIONAL, INTENT(OUT) :: variance - REAL(dp), OPTIONAL, INTENT(OUT) :: skewness - REAL(dp), OPTIONAL, INTENT(OUT) :: kurtosis - REAL(dp), OPTIONAL, INTENT(OUT) :: mean - REAL(dp), OPTIONAL, INTENT(OUT) :: stddev - REAL(dp), OPTIONAL, INTENT(OUT) :: absdev - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + REAL(dp), OPTIONAL, INTENT(OUT) :: average + REAL(dp), OPTIONAL, INTENT(OUT) :: variance + REAL(dp), OPTIONAL, INTENT(OUT) :: skewness + REAL(dp), OPTIONAL, INTENT(OUT) :: kurtosis + REAL(dp), OPTIONAL, INTENT(OUT) :: mean + REAL(dp), OPTIONAL, INTENT(OUT) :: stddev + REAL(dp), OPTIONAL, INTENT(OUT) :: absdev + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + LOGICAL, OPTIONAL, INTENT(IN) :: sample - REAL(dp) :: n + REAL(dp) :: n, div_n ! divisor depending if sample or population moments REAL(dp) :: ep, ave, var REAL(dp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error moment_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error moment_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'moment_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_dp: n must be at least 2' + + ! set divisor for population (n) or sample (n-1) variance + div_n = n - 1.0_dp + if (present(sample)) then + if (.not. sample) div_n = n + end if ! Any optional argument if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return + present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return ! Average - ave = sum(dat(:), mask=maske)/n + ave = sum(dat(:), mask = maske) / n if (present(average)) average = ave - if (present(mean)) mean = ave + if (present(mean)) mean = ave if (.not. (present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(stddev) .or. present(absdev))) return + present(kurtosis) .or. present(stddev) .or. present(absdev))) return ! Absolute deviation - s(:) = dat(:)-ave - if (present(absdev)) absdev = sum(abs(s(:)), mask=maske)/n + s(:) = dat(:) - ave + if (present(absdev)) absdev = sum(abs(s(:)), mask = maske) / n ! Variance / Standard deviation if (.not. (present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(stddev))) return - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_dp) + present(kurtosis) .or. present(stddev))) return + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / div_n if (present(variance)) variance = var ! Standard deviation - if (present(stddev)) stddev = sqrt(var) + if (present(stddev)) stddev = sqrt(var) if (.not. (present(skewness) .or. present(kurtosis))) return ! Skewness if (abs(var) .lt. tiny(0.0_dp)) stop 'moment_dp: no skewness or kurtosis when zero variance' - p(:) = p(:)*s(:) + p(:) = p(:) * s(:) if (present(skewness)) then - skewness = sum(p(:), mask=maske) - skewness = skewness/(n*stddev*stddev*stddev) - endif + skewness = sum(p(:), mask = maske) + skewness = skewness / (n * stddev * stddev * stddev) + end if ! Kurtosis if (present(kurtosis)) then - p(:) = p(:)*s(:) - kurtosis = sum(p(:), mask=maske) - kurtosis = kurtosis/(n*variance*variance) - 3.0_dp + p(:) = p(:) * s(:) + kurtosis = sum(p(:), mask = maske) + kurtosis = kurtosis / (n * variance * variance) - 3.0_dp end if END SUBROUTINE moment_dp - SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask) + SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - REAL(sp), OPTIONAL, INTENT(OUT) :: average - REAL(sp), OPTIONAL, INTENT(OUT) :: variance - REAL(sp), OPTIONAL, INTENT(OUT) :: skewness - REAL(sp), OPTIONAL, INTENT(OUT) :: kurtosis - REAL(sp), OPTIONAL, INTENT(OUT) :: mean - REAL(sp), OPTIONAL, INTENT(OUT) :: stddev - REAL(sp), OPTIONAL, INTENT(OUT) :: absdev - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + REAL(sp), OPTIONAL, INTENT(OUT) :: average + REAL(sp), OPTIONAL, INTENT(OUT) :: variance + REAL(sp), OPTIONAL, INTENT(OUT) :: skewness + REAL(sp), OPTIONAL, INTENT(OUT) :: kurtosis + REAL(sp), OPTIONAL, INTENT(OUT) :: mean + REAL(sp), OPTIONAL, INTENT(OUT) :: stddev + REAL(sp), OPTIONAL, INTENT(OUT) :: absdev + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + LOGICAL, OPTIONAL, INTENT(IN) :: sample - REAL(sp) :: n + REAL(sp) :: n, div_n ! divisor depending if sample or population moments REAL(sp) :: ep, ave, var REAL(sp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error moment_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error moment_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'moment_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_sp: n must be at least 2' + + ! set divisor for population (n) or sample (n-1) variance + div_n = n - 1.0_sp + if (present(sample)) then + if (.not. sample) div_n = n + end if ! Any optional argument if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return + present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return ! Average - ave = sum(dat(:), mask=maske)/n + ave = sum(dat(:), mask = maske) / n if (present(average)) average = ave - if (present(mean)) mean = ave + if (present(mean)) mean = ave if (.not. (present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(stddev) .or. present(absdev))) return + present(kurtosis) .or. present(stddev) .or. present(absdev))) return ! Absolute deviation - s(:) = dat(:)-ave - if (present(absdev)) absdev = sum(abs(s(:)), mask=maske)/n + s(:) = dat(:) - ave + if (present(absdev)) absdev = sum(abs(s(:)), mask = maske) / n ! Variance / Standard deviation if (.not. (present(variance) .or. present(skewness) .or. & - present(kurtosis) .or. present(stddev))) return - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_sp) + present(kurtosis) .or. present(stddev))) return + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / div_n if (present(variance)) variance = var ! Standard deviation - if (present(stddev)) stddev = sqrt(var) + if (present(stddev)) stddev = sqrt(var) if (.not. (present(skewness) .or. present(kurtosis))) return ! Skewness if (abs(var) .lt. tiny(0.0_sp)) stop 'moment_sp: no skewness or kurtosis when zero variance' - p(:) = p(:)*s(:) + p(:) = p(:) * s(:) if (present(skewness)) then - skewness = sum(p(:), mask=maske) - skewness = skewness/(n*stddev*stddev*stddev) - endif + skewness = sum(p(:), mask = maske) + skewness = skewness / (n * stddev * stddev * stddev) + end if ! Kurtosis if (present(kurtosis)) then - p(:) = p(:)*s(:) - kurtosis = sum(p(:), mask=maske) - kurtosis = kurtosis/(n*variance*variance) - 3.0_sp + p(:) = p(:) * s(:) + kurtosis = sum(p(:), mask = maske) + kurtosis = kurtosis / (n * variance * variance) - 3.0_sp end if END SUBROUTINE moment_sp @@ -1734,34 +1695,34 @@ FUNCTION stddev_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: stddev_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: stddev_dp REAL(dp) :: n REAL(dp) :: ep, ave, var REAL(dp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error stddev_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error stddev_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'stddev_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'stddev_dp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_dp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_dp) stddev_dp = sqrt(var) END FUNCTION stddev_dp @@ -1771,34 +1732,34 @@ FUNCTION stddev_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: stddev_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: stddev_sp REAL(sp) :: n REAL(sp) :: ep, ave, var REAL(sp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error stddev_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error stddev_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'stddev_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'stddev_sp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_sp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_sp) stddev_sp = sqrt(var) END FUNCTION stddev_sp @@ -1809,40 +1770,40 @@ FUNCTION skewness_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: skewness_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: skewness_dp REAL(dp) :: n REAL(dp) :: ep, ave, var, stddev REAL(dp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error skewness_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error skewness_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'skewness_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'skewness_dp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_dp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_dp) stddev = sqrt(var) ! Skewness if (abs(var) .lt. tiny(0.0_dp)) stop 'skewness_dp: no skewness when zero variance' - p(:) = p(:)*s(:) - skewness_dp = sum(p(:), mask=maske) - skewness_dp = skewness_dp/(n*stddev*stddev*stddev) + p(:) = p(:) * s(:) + skewness_dp = sum(p(:), mask = maske) + skewness_dp = skewness_dp / (n * stddev * stddev * stddev) END FUNCTION skewness_dp @@ -1851,40 +1812,40 @@ FUNCTION skewness_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: skewness_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: skewness_sp REAL(sp) :: n REAL(sp) :: ep, ave, var, stddev REAL(sp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error skewness_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error skewness_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'skewness_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'skewness_sp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - var = (var-ep*ep/n)/(n-1.0_sp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + var = (var - ep * ep / n) / (n - 1.0_sp) stddev = sqrt(var) ! Skewness if (abs(var) .lt. tiny(0.0_sp)) stop 'skewness_sp: no skewness when zero variance' - p(:) = p(:)*s(:) - skewness_sp = sum(p(:), mask=maske) - skewness_sp = skewness_sp/(n*stddev*stddev*stddev) + p(:) = p(:) * s(:) + skewness_sp = sum(p(:), mask = maske) + skewness_sp = skewness_sp / (n * stddev * stddev * stddev) END FUNCTION skewness_sp @@ -1894,34 +1855,34 @@ FUNCTION variance_dp(dat, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: variance_dp + REAL(dp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(dp) :: variance_dp REAL(dp) :: n REAL(dp) :: ep, ave, var REAL(dp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error variance_dp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),dp) + if (size(mask) .ne. size(dat)) stop 'Error variance_dp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), dp) else - maske(:) = .true. - n = real(size(dat),dp) - endif - if (n .le. (1.0_dp+tiny(1.0_dp))) stop 'variance_dp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), dp) + end if + if (n .le. (1.0_dp + tiny(1.0_dp))) stop 'variance_dp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - variance_dp = (var-ep*ep/n)/(n-1.0_dp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + variance_dp = (var - ep * ep / n) / (n - 1.0_dp) END FUNCTION variance_dp @@ -1930,34 +1891,34 @@ FUNCTION variance_sp(dat, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(IN) :: dat - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: variance_sp + REAL(sp), DIMENSION(:), INTENT(IN) :: dat + LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + REAL(sp) :: variance_sp REAL(sp) :: n REAL(sp) :: ep, ave, var REAL(sp), DIMENSION(size(dat)) :: p, s - LOGICAL, DIMENSION(size(dat)) :: maske + LOGICAL, DIMENSION(size(dat)) :: maske if (present(mask)) then - if (size(mask) .ne. size(dat)) stop 'Error variance_sp: size(mask) .ne. size(dat)' - maske = mask - n = real(count(maske),sp) + if (size(mask) .ne. size(dat)) stop 'Error variance_sp: size(mask) .ne. size(dat)' + maske = mask + n = real(count(maske), sp) else - maske(:) = .true. - n = real(size(dat),sp) - endif - if (n .le. (1.0_sp+tiny(1.0_sp))) stop 'variance_sp: n must be at least 2' + maske(:) = .true. + n = real(size(dat), sp) + end if + if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'variance_sp: n must be at least 2' ! Average - ave = sum(dat(:), mask=maske)/n - s(:) = dat(:)-ave + ave = sum(dat(:), mask = maske) / n + s(:) = dat(:) - ave ! Variance / Standard deviation - ep = sum(s(:), mask=maske) - p(:) = s(:)*s(:) - var = sum(p(:), mask=maske) - variance_sp = (var-ep*ep/n)/(n-1.0_sp) + ep = sum(s(:), mask = maske) + p(:) = s(:) * s(:) + var = sum(p(:), mask = maske) + variance_sp = (var - ep * ep / n) / (n - 1.0_sp) END FUNCTION variance_sp diff --git a/src/lib/mo_boxcox.f90 b/src/lib/mo_boxcox.f90 index 6ea272a9..a3adfadd 100644 --- a/src/lib/mo_boxcox.f90 +++ b/src/lib/mo_boxcox.f90 @@ -4,14 +4,12 @@ MODULE mo_boxcox ! as well as estimating the best exponent for the Box-Cox transformation ! Usage: - ! USE mo_boxcox, ONLY: boxcox, get_boxcox - ! lmbda = get_boxcox(data) + ! USE mo_boxcox, ONLY: boxcox, invboxcox ! new_data = boxcox(data, lmbda) ! data = invboxcox(new_data, lmbda) ! Written March 2011, Matthias Cuntz - ! - modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox, get_boxcox - ! - modified numerical recipes: brent, mnbrak, swap, shft + ! - modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox ! License ! ------- @@ -33,44 +31,15 @@ MODULE mo_boxcox ! Copyright 2011-2012 Matthias Cuntz, Juliane Mai + ! History: + ! Dec 2019: Robert Schweppe: - removed NR code (get_boxcox) - ! Note on Numerical Recipes License - ! --------------------------------- - ! Be aware that some code is under the Numerical Recipes License 3rd - ! edition - - ! The Numerical Recipes Personal Single-User License lets you personally - ! use Numerical Recipes code ("the code") on any number of computers, - ! but only one computer at a time. You are not permitted to allow anyone - ! else to access or use the code. You may, under this license, transfer - ! precompiled, executable applications incorporating the code to other, - ! unlicensed, persons, providing that (i) the application is - ! noncommercial (i.e., does not involve the selling or licensing of the - ! application for a fee), and (ii) the application was first developed, - ! compiled, and successfully run by you, and (iii) the code is bound - ! into the application in such a manner that it cannot be accessed as - ! individual routines and cannot practicably be unbound and used in - ! other programs. That is, under this license, your application user - ! must not be able to use Numerical Recipes code as part of a program - ! library or "mix and match" workbench. - - ! Businesses and organizations that purchase the disk or code download, - ! and that thus acquire one or more Numerical Recipes Personal - ! Single-User Licenses, may permanently assign those licenses, in the - ! number acquired, to individual employees. Such an assignment must be - ! made before the code is first used and, once made, it is irrevocable - ! and can not be transferred. - - ! If you do not hold a Numerical Recipes License, this code is only for - ! informational and educational purposes but cannot be used. - - USE mo_kind, ONLY: i4, sp, dp - USE mo_utils, only: eq, le, ge + USE mo_kind, ONLY : i4, sp, dp + USE mo_utils, only : eq, le, ge IMPLICIT NONE PUBLIC :: boxcox ! Calculate Box-Cox power transformed values given the original values and the exponent lambda - PUBLIC :: get_boxcox ! Find the exponent that maximises the log-likelihood function PUBLIC :: invboxcox ! Calculate the inverse Box-Cox given the transformed values and the exponent lambda ! ------------------------------------------------------------------ @@ -124,66 +93,11 @@ MODULE mo_boxcox ! - Modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox, get_boxcox ! - Modified numerical recipes: brent, mnbrak, swap, shft INTERFACE boxcox - MODULE PROCEDURE boxcox_sp, boxcox_dp + MODULE PROCEDURE boxcox_sp, boxcox_dp END INTERFACE boxcox ! ------------------------------------------------------------------ - ! NAME - ! get_boxcox - - ! PURPOSE - ! Get lambda for Box-Cox transformation. - ! - ! Transform a positive dataset with a Box-Cox power transformation. - ! boxcox(x) = (x**lambda - 1)/lambda if lambda <> 0 - ! ln(x) if lambda = 0 - ! - ! If an optinal mask is given, then the Box-Cox transformation is only performed on - ! those locations that correspond to true values in the mask. - ! x can be single or double precision. The result will have the same numerical precision. - - ! CALLING SEQUENCE - ! out = get_boxcox(x, mask=mask) - - ! INTENT(IN) - ! real(sp/dp) :: x(:) 1D-array with input numbers (>0.) - - ! INTENT(INOUT) - ! None - - ! INTENT(OUT) - ! real(sp/dp) :: lambda Lambda for power transformed values (at mask=.True.) - - ! INTENT(IN), OPTIONAL - ! logical :: mask(:) 1D-array of logical values with size(x). - ! If present, only those locations in vec corresponding to the true values in mask are used. - - ! INTENT(INOUT), OPTIONAL - ! None - - ! INTENT(OUT), OPTIONAL - ! None - - ! RESTRICTIONS - ! Input values must be positive, i.e. > 0. - - ! EXAMPLE - ! lmbda = get_boxcox(data, mask=(data>0.)) - ! new_data = boxcox(data, lmbda, mask=(data>0.)) - ! idata = invboxcox(new_data, lmbda, mask=(data>0.)) - ! -> see also example in test directory - - ! HISTORY - ! Written, Matthias Cuntz, March 2011 - ! - Modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox, get_boxcox - ! - Modified numerical recipes: brent, mnbrak, swap, shft - INTERFACE get_boxcox - MODULE PROCEDURE get_boxcox_sp, get_boxcox_dp - END INTERFACE get_boxcox - - ! ------------------------------------------------------------------ - ! NAME ! invboxcox @@ -256,29 +170,29 @@ FUNCTION boxcox_sp(x, lmbda, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(in) :: x - REAL(sp), INTENT(in) :: lmbda - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(sp), DIMENSION(size(x)) :: boxcox_sp + REAL(sp), DIMENSION(:), INTENT(in) :: x + REAL(sp), INTENT(in) :: lmbda + LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask + REAL(sp), DIMENSION(size(x)) :: boxcox_sp LOGICAL, DIMENSION(size(x)) :: maske - REAL(sp) :: lmbda1 + REAL(sp) :: lmbda1 maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error boxcox_sp: size(mask) /= size(x)' - maske = mask + if (size(mask) /= size(x)) stop 'Error boxcox_sp: size(mask) /= size(x)' + maske = mask endif - if (any((le(x,0.0_sp)) .and. maske)) stop 'Error boxcox_sp: x <= 0' - if (abs(lmbda) .lt. tiny(0.0_sp)) then - where (maske) - boxcox_sp = log(x) - elsewhere - boxcox_sp = x - end where + if (any((le(x, 0.0_sp)) .and. maske)) stop 'Error boxcox_sp: x <= 0' + if (abs(lmbda) < tiny(0.0_sp)) then + where (maske) + boxcox_sp = log(x) + elsewhere + boxcox_sp = x + end where else - lmbda1 = 1.0_sp / lmbda - boxcox_sp = merge((exp(lmbda*log(x)) - 1.0_sp) * lmbda1, x, maske) + lmbda1 = 1.0_sp / lmbda + boxcox_sp = merge((exp(lmbda * log(x)) - 1.0_sp) * lmbda1, x, maske) endif END FUNCTION boxcox_sp @@ -288,347 +202,70 @@ FUNCTION boxcox_dp(x, lmbda, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(in) :: x - REAL(dp), INTENT(in) :: lmbda - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(dp), DIMENSION(size(x)) :: boxcox_dp + REAL(dp), DIMENSION(:), INTENT(in) :: x + REAL(dp), INTENT(in) :: lmbda + LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask + REAL(dp), DIMENSION(size(x)) :: boxcox_dp LOGICAL, DIMENSION(size(x)) :: maske - REAL(dp) :: lmbda1 + REAL(dp) :: lmbda1 maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error boxcox_dp: size(mask) /= size(x)' - maske = mask + if (size(mask) /= size(x)) stop 'Error boxcox_dp: size(mask) /= size(x)' + maske = mask endif - if (any((le(x,0.0_dp)) .and. maske)) then - stop 'Error boxcox_dp: x <= 0' + if (any((le(x, 0.0_dp)) .and. maske)) then + stop 'Error boxcox_dp: x <= 0' end if - if (abs(lmbda) .lt. tiny(0.0_dp)) then - where (maske) - boxcox_dp = log(x) - elsewhere - boxcox_dp = x - end where + if (abs(lmbda) < tiny(0.0_dp)) then + where (maske) + boxcox_dp = log(x) + elsewhere + boxcox_dp = x + end where else - lmbda1 = 1.0_dp / lmbda - boxcox_dp = merge((exp(lmbda*log(x)) - 1.0_dp) * lmbda1, x, maske) + lmbda1 = 1.0_dp / lmbda + boxcox_dp = merge((exp(lmbda * log(x)) - 1.0_dp) * lmbda1, x, maske) endif END FUNCTION boxcox_dp ! ------------------------------------------------------------------ - ! Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is - ! between ax and cx, and f(bx) is less than both f(ax) and f(cx)), this routine isolates - ! the minimum to a fractional precision of about tol using Brent''s method. The abscissa of - ! the minimum is returned. - ! Parameters: Maximum allowed number of iterations; goldenratio; and a small number that - ! protects against trying to achieve fractional accuracy for a minimum that happens to be - ! exactly zero. - Function brent_sp(ax, bx, cx, func_sp, tol, dat) - - IMPLICIT NONE - - REAL(sp), INTENT(in) :: ax, bx, cx - INTERFACE - FUNCTION func_sp(x, da) - USE mo_kind, ONLY: sp - IMPLICIT NONE - REAL(sp), INTENT(in) :: x - REAL(sp), DIMENSION(:), INTENT(in) :: da - REAL(sp) :: func_sp - END FUNCTION func_sp - END INTERFACE - REAL(sp), INTENT(in) :: tol - REAL(sp), DIMENSION(:), INTENT(in) :: dat - REAL(sp) :: brent_sp - - INTEGER(i4), PARAMETER :: ITMAX = 500 - REAL(sp), PARAMETER :: CGOLD = 0.3819660_sp - REAL(sp), PARAMETER :: ZEPS = 1.0e-3_sp*epsilon(ax) - INTEGER(i4) :: iter - REAL(sp) :: a, b, d, e, etemp, fu, fv, fw, fx - REAL(sp) :: p, q, r, tol1, tol2, u, v, w, x, xm - - d = 0.0_sp - - a = min(ax,cx) - b = max(ax,cx) - v = bx - w = v - x = v - e = 0.0_sp - fx = func_sp(x, dat) - fv = fx - fw = fx - do iter=1, ITMAX - xm = 0.5_sp*(a+b) - tol1 = tol*abs(x) + ZEPS - tol2 = 2.0_sp*tol1 - if (le(abs(x-xm),(tol2-0.5_sp*(b-a)))) then - brent_sp = x - return - end if - if (abs(e) > tol1) then - r = (x-w)*(fx-fv) - q = (x-v)*(fx-fw) - p = (x-v)*q-(x-w)*r - q = 2.0_sp*(q-r) - if (q > 0.0_sp) p = -p - q = abs(q) - etemp = e - e = d - if (ge(abs(p),abs(0.5_sp*q*etemp)) .or. & - le(p,q*(a-x)) .or. ge(p,q*(b-x))) then - e = merge(a-x, b-x, ge(x,xm)) - d = CGOLD*e - else - d = p/q - u = x+d - if (u-a < tol2 .or. b-u < tol2) d = sign(tol1,xm-x) - end if - else - e = merge(a-x, b-x, ge(x,xm)) - d = CGOLD*e - end if - u = merge(x+d, x+sign(tol1,d), ge(abs(d),tol1)) - fu = func_sp(u, dat) - if (le(fu,fx)) then - if (ge(u,x)) then - a = x - else - b = x - end if - call shft_sp(v,w,x,u) - call shft_sp(fv,fw,fx,fu) - else - if (u < x) then - a = u - else - b = u - end if - if (le(fu,fw) .or. eq(w,x)) then - v = w - fv = fw - w = u - fw = fu - else if (le(fu,fv) .or. eq(v,x) .or. eq(v,w)) then - v = u - fv = fu - end if - end if - end do - ! If it comes to here, there were too much iterations - STOP 'brent_sp: exceed maximum iterations' - - END FUNCTION brent_sp - - - FUNCTION brent_dp(ax, bx, cx, func_dp, tol, dat) - - IMPLICIT NONE - - REAL(dp), INTENT(in) :: ax, bx, cx - INTERFACE - FUNCTION func_dp(x, da) - USE mo_kind, ONLY: dp - IMPLICIT NONE - REAL(dp), INTENT(in) :: x - REAL(dp), DIMENSION(:), INTENT(in) :: da - REAL(dp) :: func_dp - END FUNCTION func_dp - END INTERFACE - REAL(dp), INTENT(in) :: tol - REAL(dp), DIMENSION(:), INTENT(in) :: dat - REAL(dp) :: brent_dp - - INTEGER(i4), PARAMETER :: ITMAX=500 - REAL(dp), PARAMETER :: CGOLD=0.3819660_dp - REAL(dp), PARAMETER :: ZEPS=1.0e-3_dp*epsilon(ax) - INTEGER(i4) :: iter - REAL(dp) :: a, b, d, e, etemp, fu, fv, fw, fx - REAL(dp) :: p, q, r, tol1, tol2, u, v, w, x, xm - - d = 0.0_dp - - a = min(ax,cx) - b = max(ax,cx) - v = bx - w = v - x = v - e = 0.0_dp - fx = func_dp(x, dat) - fv = fx - fw = fx - do iter=1, ITMAX - xm = 0.5_dp*(a+b) - tol1 = tol*abs(x) + ZEPS - tol2 = 2.0_dp*tol1 - if (le(abs(x-xm),(tol2-0.5_dp*(b-a)))) then - brent_dp = x - return - end if - if (abs(e) > tol1) then - r = (x-w)*(fx-fv) - q = (x-v)*(fx-fw) - p = (x-v)*q-(x-w)*r - q = 2.0_dp*(q-r) - if (q > 0.0_dp) p = -p - q = abs(q) - etemp = e - e = d - if (ge(abs(p),abs(0.5_dp*q*etemp)) .or. & - le(p,q*(a-x)) .or. ge(p,q*(b-x))) then - e = merge(a-x, b-x, ge(x,xm)) - d = CGOLD*e - else - d = p/q - u = x+d - if (u-a < tol2 .or. b-u < tol2) d = sign(tol1,xm-x) - end if - else - e = merge(a-x, b-x, ge(x,xm)) - d = CGOLD*e - end if - u = merge(x+d, x+sign(tol1,d), ge(abs(d),tol1)) - fu = func_dp(u, dat) - if (le(fu,fx)) then - if (ge(u,x)) then - a = x - else - b = x - end if - call shft_dp(v,w,x,u) - call shft_dp(fv,fw,fx,fu) - else - if (u < x) then - a = u - else - b = u - end if - if (le(fu,fw) .or. eq(w,x)) then - v = w - fv = fw - w = u - fw = fu - else if (le(fu,fv) .or. eq(v,x) .or. eq(v,w)) then - v = u - fv = fu - end if - end if - end do - ! If it comes to here, there were too much iterations - STOP 'brent_dp: exceed maximum iterations' - - END FUNCTION brent_dp - - ! ------------------------------------------------------------------ - - ! Find the lambda that maximizes the log-likelihood function - FUNCTION get_boxcox_sp(x, mask) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(in) :: x - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(sp) :: get_boxcox_sp - - REAL(sp), PARAMETER :: tol = 1.48e-08 - REAL(sp) :: ax, bx, cx - REAL(sp) :: fa, fb, fc - LOGICAL, DIMENSION(size(x)) :: maske - REAL(sp), DIMENSION(:), ALLOCATABLE :: xx - INTEGER(i4) :: nn - - - maske(:) = .true. - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error get_boxcox_sp: size(mask) /= size(x)' - maske = mask - endif - nn = count(maske) - if (allocated(xx)) deallocate(xx) - allocate(xx(nn)) - xx = pack(x, mask=maske) - ax = -2.0_sp - bx = 2.0_sp - call mnbrak_sp(ax, bx, cx, fa, fb, fc, llf_boxcox_sp, xx) - get_boxcox_sp = brent_sp(ax, bx, cx, llf_boxcox_sp, tol, xx) - - END FUNCTION get_boxcox_sp - - - FUNCTION get_boxcox_dp(x, mask) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(in) :: x - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(dp) :: get_boxcox_dp - - REAL(dp), PARAMETER :: tol = 1.48e-08 - REAL(dp) :: ax, bx, cx - REAL(dp) :: fa, fb, fc - LOGICAL, DIMENSION(size(x)) :: maske - REAL(dp), DIMENSION(:), ALLOCATABLE :: xx - INTEGER(i4) :: nn - - - maske(:) = .true. - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error get_boxcox_dp: size(mask) /= size(x)' - maske = mask - endif - nn = count(maske) - if (allocated(xx)) deallocate(xx) - allocate(xx(nn)) - xx = pack(x, mask=maske) - if (any(le(xx,0.0_dp))) then - stop 'Error get_boxcox_dp: x <= 0' - end if - ax = -5.0_dp - bx = 5.0_dp - call mnbrak_dp(ax, bx, cx, fa, fb, fc, llf_boxcox_dp, xx) - get_boxcox_dp = brent_dp(ax, bx, cx, llf_boxcox_dp, tol, xx) - - END FUNCTION get_boxcox_dp - - ! ------------------------------------------------------------------ - FUNCTION invboxcox_1d_sp(x, lmbda, mask) IMPLICIT NONE - REAL(sp), DIMENSION(:), INTENT(in) :: x - REAL(sp), INTENT(in) :: lmbda - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(sp), DIMENSION(size(x)) :: invboxcox_1d_sp + REAL(sp), DIMENSION(:), INTENT(in) :: x + REAL(sp), INTENT(in) :: lmbda + LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask + REAL(sp), DIMENSION(size(x)) :: invboxcox_1d_sp - LOGICAL, DIMENSION(size(x)) :: maske - REAL(sp) :: lmbda1 + LOGICAL, DIMENSION(size(x)) :: maske + REAL(sp) :: lmbda1 REAL(sp), DIMENSION(size(x)) :: temp maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error invboxcox_1d_sp: size(mask) /= size(x)' - maske = mask + if (size(mask) /= size(x)) stop 'Error invboxcox_1d_sp: size(mask) /= size(x)' + maske = mask endif - if (abs(lmbda) .lt. tiny(0.0_sp)) then - where (maske) - invboxcox_1d_sp = exp(x) - elsewhere - invboxcox_1d_sp = x - end where + if (abs(lmbda) < tiny(0.0_sp)) then + where (maske) + invboxcox_1d_sp = exp(x) + elsewhere + invboxcox_1d_sp = x + end where else - lmbda1 = 1.0_sp / lmbda - temp = lmbda*x+1.0_sp - where (temp .gt. 0.0_sp) - temp = exp(lmbda1*log(temp)) - elsewhere - temp = 0.0_sp - end where - invboxcox_1d_sp = merge(temp, x, maske) + lmbda1 = 1.0_sp / lmbda + temp = lmbda * x + 1.0_sp + where (temp > 0.0_sp) + temp = exp(lmbda1 * log(temp)) + elsewhere + temp = 0.0_sp + end where + invboxcox_1d_sp = merge(temp, x, maske) endif END FUNCTION invboxcox_1d_sp @@ -637,35 +274,35 @@ FUNCTION invboxcox_1d_dp(x, lmbda, mask) IMPLICIT NONE - REAL(dp), DIMENSION(:), INTENT(in) :: x - REAL(dp), INTENT(in) :: lmbda - LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask - REAL(dp), DIMENSION(size(x)) :: invboxcox_1d_dp + REAL(dp), DIMENSION(:), INTENT(in) :: x + REAL(dp), INTENT(in) :: lmbda + LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask + REAL(dp), DIMENSION(size(x)) :: invboxcox_1d_dp - LOGICAL, DIMENSION(size(x)) :: maske - REAL(dp) :: lmbda1 + LOGICAL, DIMENSION(size(x)) :: maske + REAL(dp) :: lmbda1 REAL(dp), DIMENSION(size(x)) :: temp maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error invboxcox_1d_dp: size(mask) /= size(x)' - maske = mask + if (size(mask) /= size(x)) stop 'Error invboxcox_1d_dp: size(mask) /= size(x)' + maske = mask endif - if (abs(lmbda) .lt. tiny(0.0_dp)) then - where (maske) - invboxcox_1d_dp = exp(x) - elsewhere - invboxcox_1d_dp = x - end where + if (abs(lmbda) < tiny(0.0_dp)) then + where (maske) + invboxcox_1d_dp = exp(x) + elsewhere + invboxcox_1d_dp = x + end where else - lmbda1 = 1.0_dp / lmbda - temp = lmbda*x+1.0_dp - where (temp .gt. 0.0_dp) - temp = exp(lmbda1*log(temp)) - elsewhere - temp = 0.0_dp - end where - invboxcox_1d_dp = merge(temp, x, maske) + lmbda1 = 1.0_dp / lmbda + temp = lmbda * x + 1.0_dp + where (temp > 0.0_dp) + temp = exp(lmbda1 * log(temp)) + elsewhere + temp = 0.0_dp + end where + invboxcox_1d_dp = merge(temp, x, maske) endif END FUNCTION invboxcox_1d_dp @@ -674,24 +311,24 @@ FUNCTION invboxcox_0d_sp(x, lmbda) IMPLICIT NONE - REAL(sp), INTENT(in) :: x - REAL(sp), INTENT(in) :: lmbda - REAL(sp) :: invboxcox_0d_sp + REAL(sp), INTENT(in) :: x + REAL(sp), INTENT(in) :: lmbda + REAL(sp) :: invboxcox_0d_sp - REAL(sp) :: lmbda1 - REAL(sp) :: temp + REAL(sp) :: lmbda1 + REAL(sp) :: temp - if (abs(lmbda) .lt. tiny(0.0_sp)) then - invboxcox_0d_sp = exp(x) + if (abs(lmbda) < tiny(0.0_sp)) then + invboxcox_0d_sp = exp(x) else - lmbda1 = 1.0_sp / lmbda - temp = lmbda*x+1.0_sp - if (temp .gt. 0.0_sp) then - temp = exp(lmbda1*log(temp)) - else - temp = 0.0_sp - end if - invboxcox_0d_sp = temp + lmbda1 = 1.0_sp / lmbda + temp = lmbda * x + 1.0_sp + if (temp > 0.0_sp) then + temp = exp(lmbda1 * log(temp)) + else + temp = 0.0_sp + end if + invboxcox_0d_sp = temp endif END FUNCTION invboxcox_0d_sp @@ -700,292 +337,26 @@ FUNCTION invboxcox_0d_dp(x, lmbda) IMPLICIT NONE - REAL(dp), INTENT(in) :: x - REAL(dp), INTENT(in) :: lmbda - REAL(dp) :: invboxcox_0d_dp + REAL(dp), INTENT(in) :: x + REAL(dp), INTENT(in) :: lmbda + REAL(dp) :: invboxcox_0d_dp - REAL(dp) :: lmbda1 - REAL(dp) :: temp + REAL(dp) :: lmbda1 + REAL(dp) :: temp - if (abs(lmbda) .lt. tiny(0.0_dp)) then - invboxcox_0d_dp = exp(x) + if (abs(lmbda) < tiny(0.0_dp)) then + invboxcox_0d_dp = exp(x) else - lmbda1 = 1.0_dp / lmbda - temp = lmbda*x+1.0_dp - if (temp .gt. 0.0_dp) then - temp = exp(lmbda1*log(temp)) - else - temp = 0.0_dp - end if - invboxcox_0d_dp = temp + lmbda1 = 1.0_dp / lmbda + temp = lmbda * x + 1.0_dp + if (temp > 0.0_dp) then + temp = exp(lmbda1 * log(temp)) + else + temp = 0.0_dp + end if + invboxcox_0d_dp = temp endif END FUNCTION invboxcox_0d_dp - ! ------------------------------------------------------------------ - - ! The boxcox (negative) log-likelihood function. - FUNCTION llf_boxcox_sp(lmbda, x) - - IMPLICIT NONE - - REAL(sp), INTENT(in) :: lmbda - REAL(sp), DIMENSION(:), INTENT(in) :: x - REAL(sp) :: llf_boxcox_sp - - REAL(sp), DIMENSION(size(x)) :: y - REAL(sp) :: N, my, f, s - - N = real(size(x),sp) - y = boxcox_sp(x, lmbda) - my = sum(y) / N - f = (lmbda-1.0_sp) * sum(log(x)) - - ! f = f - 0.5_dp*N * log(sum((y-my)*(y-my))/N) - s = Min( huge(1.0_sp)/N , Max( N*tiny(1.0_sp) , sum((y-my)*(y-my)) ) ) - f = f - 0.5_sp*N * log(s/N) - - llf_boxcox_sp = -f - - END FUNCTION llf_boxcox_sp - - - FUNCTION llf_boxcox_dp(lmbda, x) - - IMPLICIT NONE - - REAL(dp), INTENT(in) :: lmbda - REAL(dp), DIMENSION(:), INTENT(in) :: x - REAL(dp) :: llf_boxcox_dp - - REAL(dp), DIMENSION(size(x)) :: y - REAL(dp) :: N, my, f, s - - N = real(size(x),dp) - y = boxcox_dp(x, lmbda) - my = sum(y) / N - f = (lmbda-1.0_dp) * sum(log(x)) - - ! f = f - 0.5_dp*N * log(sum((y-my)*(y-my))/N) - s = Min( huge(1.0_dp)/N , Max( N*tiny(1.0_dp) , sum((y-my)*(y-my)) ) ) - f = f - 0.5_dp*N * log(s/N) - - llf_boxcox_dp = -f - - END FUNCTION llf_boxcox_dp - - ! ------------------------------------------------------------------ - - ! Given a function func, and given distinct initial points ax and bx, this routine searches - ! in the downhill direction (defined by the functionas evaluated at the initial points) and - ! returns new points ax, bx, cx that bracket a minimum of the function. - ! Parameters: GOLD is the default ratio by which successive intervals are magnified; GLIMIT - ! is the maximum magnification allowed for a parabolic-fit step. - SUBROUTINE mnbrak_sp(ax, bx, cx, fa, fb, fc, func_sp, dat) - - IMPLICIT NONE - - REAL(sp), INTENT(inout) :: ax, bx - REAL(sp), INTENT(out) :: cx - REAL(sp), INTENT(out) :: fa, fb, fc - INTERFACE - FUNCTION func_sp(x, da) - USE mo_kind, ONLY: sp - IMPLICIT NONE - REAL(sp), INTENT(in) :: x - REAL(sp), DIMENSION(:), INTENT(in) :: da - REAL(sp) :: func_sp - END FUNCTION func_sp - END INTERFACE - REAL(sp), DIMENSION(:), INTENT(in) :: dat - - REAL(sp), PARAMETER :: GOLD = 1.618034_sp - REAL(sp), PARAMETER :: GLIMIT = 100.0_sp - REAL(sp), PARAMETER :: TINY = 1.0e-20_sp - REAL(sp) :: fu, q, r, u, ulim - - fa = func_sp(ax, dat) - fb = func_sp(bx, dat) - if (fb > fa) then - call swap_sp(ax,bx) - call swap_sp(fa,fb) - end if - cx = bx + GOLD*(bx-ax) - fc = func_sp(cx, dat) - do - if (fb < fc) return - r = (bx-ax)*(fb-fc) - q = (bx-cx)*(fb-fa) - u = bx-((bx-cx)*q-(bx-ax)*r)/(2.0_sp*sign(max(abs(q-r),TINY),q-r)) - ulim = bx + GLIMIT*(cx-bx) - if ((bx-u)*(u-cx) > 0.0_sp) then - fu = func_sp(u, dat) - if (fu < fc) then - ax = bx - fa = fb - bx = u - fb = fu - return - else if (fu > fb) then - cx = u - fc = fu - return - end if - u = cx + GOLD*(cx-bx) - fu = func_sp(u, dat) - else if ((cx-u)*(u-ulim) > 0.0_sp) then - fu = func_sp(u, dat) - if (fu < fc) then - bx = cx - cx = u - u = cx + GOLD*(cx-bx) - call shft_sp(fb,fc,fu,func_sp(u, dat)) - end if - else if (ge((u-ulim)*(ulim-cx),0.0_sp)) then - u = ulim - fu = func_sp(u, dat) - else - u = cx + GOLD*(cx-bx) - fu = func_sp(u, dat) - end if - call shft_sp(ax,bx,cx,u) - call shft_sp(fa,fb,fc,fu) - end do - - END SUBROUTINE mnbrak_sp - - - SUBROUTINE mnbrak_dp(ax, bx, cx, fa, fb, fc, func_dp, dat) - - IMPLICIT NONE - - REAL(dp), INTENT(inout) :: ax, bx - REAL(dp), INTENT(out) :: cx - REAL(dp), INTENT(out) :: fa, fb, fc - INTERFACE - FUNCTION func_dp(x, da) - USE mo_kind, ONLY: dp - IMPLICIT NONE - REAL(dp), INTENT(in) :: x - REAL(dp), DIMENSION(:), INTENT(in) :: da - REAL(dp) :: func_dp - END FUNCTION func_dp - END INTERFACE - REAL(dp), DIMENSION(:), INTENT(in) :: dat - - REAL(dp), PARAMETER :: GOLD = 1.618034_dp - REAL(dp), PARAMETER :: GLIMIT = 100.0_dp - REAL(dp), PARAMETER :: TINY = 1.0e-20_dp - REAL(dp) :: fu, q, r, u, ulim - - fa = func_dp(ax, dat) - fb = func_dp(bx, dat) - if (fb > fa) then - call swap_dp(ax,bx) - call swap_dp(fa,fb) - end if - cx = bx + GOLD*(bx-ax) - fc = func_dp(cx, dat) - do - if (fb < fc) return - r = (bx-ax)*(fb-fc) - q = (bx-cx)*(fb-fa) - u = bx-((bx-cx)*q-(bx-ax)*r)/(2.0_dp*sign(max(abs(q-r),TINY),q-r)) - ulim = bx + GLIMIT*(cx-bx) - if ((bx-u)*(u-cx) > 0.0_dp) then - fu = func_dp(u, dat) - if (fu < fc) then - ax = bx - fa = fb - bx = u - fb = fu - return - else if (fu > fb) then - cx = u - fc = fu - return - end if - u = cx + GOLD*(cx-bx) - fu = func_dp(u, dat) - else if ((cx-u)*(u-ulim) > 0.0_dp) then - fu = func_dp(u, dat) - if (fu < fc) then - bx = cx - cx = u - u = cx + GOLD*(cx-bx) - call shft_dp(fb,fc,fu,func_dp(u, dat)) - end if - else if (ge((u-ulim)*(ulim-cx),0.0_dp)) then - u = ulim - fu = func_dp(u, dat) - else - u = cx + GOLD*(cx-bx) - fu = func_dp(u, dat) - end if - call shft_dp(ax,bx,cx,u) - call shft_dp(fa,fb,fc,fu) - end do - - END SUBROUTINE mnbrak_dp - - ! ------------------------------------------------------------------ - - ! Helper for brent and mnbrak - SUBROUTINE shft_sp(a,b,c,d) - - IMPLICIT NONE - - REAL(sp), INTENT(out) :: a - REAL(sp), INTENT(inout) :: b,c - REAL(sp), INTENT(in) :: d - - a = b - b = c - c = d - - END SUBROUTINE shft_sp - - - SUBROUTINE shft_dp(a,b,c,d) - - IMPLICIT NONE - - REAL(dp), INTENT(out) :: a - REAL(dp), INTENT(inout) :: b,c - REAL(dp), INTENT(in) :: d - - a = b - b = c - c = d - - END SUBROUTINE shft_dp - - ! ------------------------------------------------------------------ - - ! Helper for mnbrak - SUBROUTINE swap_sp(a,b) - - IMPLICIT NONE - - REAL(sp), INTENT(inout) :: a, b - REAL(sp) :: dum - dum = a - a = b - b = dum - - END SUBROUTINE swap_sp - - - SUBROUTINE swap_dp(a,b) - - IMPLICIT NONE - - REAL(dp), INTENT(inout) :: a, b - REAL(dp) :: dum - dum = a - a = b - b = dum - - END SUBROUTINE swap_dp - END MODULE mo_boxcox diff --git a/src/lib/mo_corr.f90 b/src/lib/mo_corr.f90 index 826c6247..35850a61 100644 --- a/src/lib/mo_corr.f90 +++ b/src/lib/mo_corr.f90 @@ -1,14 +1,8 @@ MODULE mo_corr - ! This module provides auto- and crosscorrelation function calculations + ! This module provides autocorrelation function calculations - ! Literature - ! Corr, FFT - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - - ! Written March 2011, Matthias Cuntz + ! Rewritten December 2019, Sebastian Mueller ! License ! ------- @@ -28,106 +22,14 @@ MODULE mo_corr ! along with the UFZ Fortran library (cf. gpl.txt and lgpl.txt). ! If not, see . - ! Copyright 2011 Matthias Cuntz - - - ! Note on Numerical Recipes License - ! --------------------------------- - ! Be aware that some code is under the Numerical Recipes License 3rd - ! edition - - ! The Numerical Recipes Personal Single-User License lets you personally - ! use Numerical Recipes code ("the code") on any number of computers, - ! but only one computer at a time. You are not permitted to allow anyone - ! else to access or use the code. You may, under this license, transfer - ! precompiled, executable applications incorporating the code to other, - ! unlicensed, persons, providing that (i) the application is - ! noncommercial (i.e., does not involve the selling or licensing of the - ! application for a fee), and (ii) the application was first developed, - ! compiled, and successfully run by you, and (iii) the code is bound - ! into the application in such a manner that it cannot be accessed as - ! individual routines and cannot practicably be unbound and used in - ! other programs. That is, under this license, your application user - ! must not be able to use Numerical Recipes code as part of a program - ! library or "mix and match" workbench. + ! Copyright 2019 Sebastian Mueller - ! Businesses and organizations that purchase the disk or code download, - ! and that thus acquire one or more Numerical Recipes Personal - ! Single-User Licenses, may permanently assign those licenses, in the - ! number acquired, to individual employees. Such an assignment must be - ! made before the code is first used and, once made, it is irrevocable - ! and can not be transferred. - ! If you do not hold a Numerical Recipes License, this code is only for - ! informational and educational purposes but cannot be used. - - USE mo_kind, ONLY : i4, sp, dp, spc, dpc - USE mo_constants, ONLY : TWOPI_sp, TWOPI_dp, PI_dp + USE mo_kind, ONLY : i4, sp, dp Implicit NONE - PUBLIC :: autocoeffk ! coeff_k so that autocorr = coeff_k/coeff_0 PUBLIC :: autocorr ! Autocorrelation coefficient at lag k = autocoeffk(k)/autocoeffk(0) - PUBLIC :: corr ! Correlation (function) with optional highpass filtering: covariance=correlation(1)/n - PUBLIC :: crosscoeffk ! coeff_k so that crosscorr = coeff_k/coeff_0, crosscoeffk(0) = covariance - PUBLIC :: crosscorr ! Crosscorrelation coefficient at lag k = crosscoeffk(k)/crosscoeffk(0) - - ! ------------------------------------------------------------------ - - ! NAME - ! autocoeffk - - ! PURPOSE - ! Coefficient at lag k so that autocorrelation coefficient at lag k - ! is autocoeffk(x,k)/autocoeffk(x,0). - ! - ! If an optinal mask is given, the calculations are only over those locations that correspond - ! to true values in the mask. - ! x can be single or double precision. The result will have the same numerical precision. - - ! CALLING SEQUENCE - ! ak = autocoeffk(x, k, mask=mask) - - ! INTENT(IN) - ! real(sp/dp) :: x(:) Time series - ! integer(i4) :: k[(:)] Lag for autocorrelation - - ! INTENT(INOUT) - ! None - - ! INTENT(OUT) - ! real(sp/dp) :: ak[(:)] coefficient so that ak/autocoeffk(x,0) is the autocorrelation coefficient at lag k - - ! INTENT(IN), OPTIONAL - ! logical :: mask(:) 1D-array of logical values with size(vec). - ! If present, only those locations in vec corresponding to the true values in mask are used. - - ! INTENT(INOUT), OPTIONAL - ! None - - ! INTENT(OUT), OPTIONAL - ! None - - ! RESTRICTIONS - ! None - - ! EXAMPLE - ! ! Variance = autocoeffk(0) - ! var = autocoeffk(x,0) - ! -> see also example in test directory - - ! LITERATURE - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - - ! HISTORY - ! Written, Matthias Cuntz, Nov 2011 - ! Modified, Stephan Thober, Nov 2012 - added 1d version - INTERFACE autocoeffk - MODULE PROCEDURE autocoeffk_sp, autocoeffk_dp, & - autocoeffk_1d_dp, autocoeffk_1d_sp - END INTERFACE autocoeffk ! ------------------------------------------------------------------ @@ -173,388 +75,95 @@ MODULE mo_corr ! acorr = autocorr(x,size(x)/2) ! -> see also example in test directory - ! LITERATURE - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 ! Modified, Stephan Thober, Nov 2012 - added 1d version + ! Modified, Sebastian Mueller, Dec 2019 - rewritten INTERFACE autocorr - MODULE PROCEDURE autocorr_sp, autocorr_dp, & - autocorr_1d_sp, autocorr_1d_dp + MODULE PROCEDURE autocorr_sp, autocorr_dp, autocorr_1d_sp, autocorr_1d_dp END INTERFACE autocorr ! ------------------------------------------------------------------ - ! NAME - ! corr - - ! PURPOSE - ! Computes the correlation of two real data sets data1 and data2 of length N (includ- - ! ing any user-supplied zero padding) with Fast Fourier Transform (FFT). N must be an integer - ! power of 2 for the FFT routine. Corr takes only the elements up to the last power of 2 and - ! returns the number of elements used in nadjust. - ! The answer is returned as the function corr, an array of length N (nadjust). - ! The answer is stored in wrap-around order, i.e., correlations at increasingly negative lags - ! are in corr(N) on down to corr(N/2+1), while correlations at increasingly positive lags are - ! in correl(1) (zero lag) on up to correl(N/2). Sign convention of this routine: if data1 lags - ! data2, i.e., is shifted to the right of it, then correl will show a peak at positive lags. - ! - ! Optional high-pass filtering of the time series in Fourier space is implemented. - ! - ! Note covariance(x,y) = corr(x,y)/n - - ! CALLING SEQUENCE - ! cfunc = corr(data1,data2,nadjust=nadjust,nhigh=nhigh,nwin=nwin) - - ! INTENT(IN) - ! real(sp/dp) :: data1(:) 1st time series - ! real(sp/dp) :: data2(:) 2nd time series - - ! INTENT(INOUT) - ! None - - ! INTENT(OUT) - ! real(sp/dp) :: corr(size(data1)) Correlation function between data1 and data2 in wrap-around order - - ! INTENT(IN), OPTIONAL - ! integer(i4) :: nhigh If >0 then nhigh upper frequencies are filtered. nwin then defines - ! the used window function for filtering. (default: 0) - ! integer(i4) :: nwin Window function for highpass filtering (default: 1) - ! 0: no filtering - ! 1: ideal highpass, i.e. cut out the nhigh upper frequencies - ! 2: linear interpolation of 0 to 1 from highest to highest-nhigh - ! frequency; similar Bartlett window - - ! INTENT(INOUT), OPTIONAL - ! None - - ! INTENT(OUT), OPTIONAL - ! integer(i4) :: nadjust Actual used number of elements. - - ! RESTRICTIONS - ! None - - ! EXAMPLE - ! ! Covariance function: covariance(x,y) = corr(x,y)/n - ! wT = corr(w, T, nadjust=nwT) - ! wT(1:nwT) = wT(1:nwT) / real(nwT,dp) - ! -> see also example in test directory - - ! LITERATURE - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - - ! HISTORY - ! Written, Matthias Cuntz, Nov 2011 - INTERFACE corr - MODULE PROCEDURE corr_sp, corr_dp - END INTERFACE corr - - ! ------------------------------------------------------------------ - - ! NAME - ! crosscoeffk - - ! PURPOSE - ! Coefficient at lag k so that crosscorrelation coefficient at lag k - ! is crosscoeffk(x,y,k)/crosscoeffk(x,y,0). - ! -> crosscoeffk(x,y,0) = covariance(x,y) - ! - ! If an optinal mask is given, the calculations are only over those locations that correspond - ! to true values in the mask. - ! x can be single or double precision. The result will have the same numerical precision. - - ! CALLING SEQUENCE - ! ck = crosscoeffk(x, y, k, mask=mask) - - ! INTENT(IN) - ! real(sp/dp) :: x(:) 1st time series - ! real(sp/dp) :: y(:) 2nd time series - ! integer(i4) :: k Lag for crosscorrelation - - ! INTENT(INOUT) - ! None - - ! INTENT(OUT) - ! real(sp/dp) :: ck coefficient so that ck/crosscoeffk(x,0) is the crosscorrelation coefficient at lag k - - ! INTENT(IN), OPTIONAL - ! logical :: mask(:) 1D-array of logical values with size(vec). - ! If present, only those locations in vec corresponding to the true values in mask are used. - - ! INTENT(INOUT), OPTIONAL - ! None - - ! INTENT(OUT), OPTIONAL - ! None - - ! RESTRICTIONS - ! None - - ! EXAMPLE - ! ! covariance = crosscoeffk(0) - ! cov = crosscoeffk(x,y,0) - ! -> see also example in test directory - - ! LITERATURE - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - - ! HISTORY - ! Written, Matthias Cuntz, Nov 2011 - INTERFACE crosscoeffk - MODULE PROCEDURE crosscoeffk_sp, crosscoeffk_dp - END INTERFACE crosscoeffk - - ! ------------------------------------------------------------------ - - ! NAME - ! crosscorr - - ! PURPOSE - ! Element at lag k of crosscorrelation function - ! crosscorr(x,y,k) = crosscoeffk(x,y,k)/crosscoeffk(x,y,0). - ! - ! If an optinal mask is given, the calculations are only over those locations that correspond - ! to true values in the mask. - ! x can be single or double precision. The result will have the same numerical precision. - - ! CALLING SEQUENCE - ! ck = crosscorr(x, y, k, mask=mask) - - ! INTENT(IN) - ! real(sp/dp) :: x(:) 1st time series - ! real(sp/dp) :: y(:) 2nd time series - ! integer(i4) :: k Lag for crosscorrelation - - ! INTENT(INOUT) - ! None - - ! INTENT(OUT) - ! real(sp/dp) :: ck Coefficient of crosscorrelation function at lag k - - ! INTENT(IN), OPTIONAL - ! logical :: mask(:) 1D-array of logical values with size(vec). - ! If present, only those locations in vec corresponding to the true values in mask are used. - - ! INTENT(INOUT), OPTIONAL - ! None - - ! INTENT(OUT), OPTIONAL - ! None - - ! RESTRICTIONS - ! None - - ! EXAMPLE - ! ! Last crosscorrelation element - ! ccorr = crosscorr(x,y,size(x)/2) - ! -> see also example in test directory - - ! LITERATURE - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 - - ! HISTORY - ! Written, Matthias Cuntz, Nov 2011 - INTERFACE crosscorr - MODULE PROCEDURE crosscorr_sp, crosscorr_dp - END INTERFACE crosscorr - - ! ------------------------------------------------------------------ - - PRIVATE - - ! ------------------------------------------------------------------ - - ! Private routines, mostly from numerical recipes - INTERFACE arth - MODULE PROCEDURE arth_sp, arth_dp, arth_i4 - END INTERFACE arth - INTERFACE four1 - MODULE PROCEDURE four1_sp, four1_dp - END INTERFACE four1 - INTERFACE fourrow - MODULE PROCEDURE fourrow_sp, fourrow_dp - END INTERFACE fourrow - INTERFACE realft - MODULE PROCEDURE realft_sp, realft_dp - END INTERFACE realft - INTERFACE swap - MODULE PROCEDURE & ! swap_i4, & - !swap_sp, swap_1d_sp, & - !swap_dp, & !swap_1d_dp, & - ! swap_spc, & - swap_1d_spc, & - ! swap_dpc, & - swap_1d_dpc !, & - ! masked_swap_sp, masked_swap_1d_sp, masked_swap_2d_sp, & - ! masked_swap_dp, masked_swap_1d_dp, masked_swap_2d_dp, & - ! masked_swap_spc, masked_swap_1d_spc, masked_swap_2d_spc, & - ! masked_swap_dpc, masked_swap_1d_dpc, masked_swap_2d_dpc - END INTERFACE swap - !INTERFACE zroots_unity - ! MODULE PROCEDURE zroots_unity_sp, zroots_unity_dp - !END INTERFACE zroots_unity - - INTEGER(i4), PARAMETER :: NPAR_ARTH = 16, NPAR2_ARTH = 8 - - ! ------------------------------------------------------------------ - CONTAINS ! ------------------------------------------------------------------ - ! From numerical recipes documentation - ! Returns an array of length n containing an arithmetic progression whose - ! first value is first and whose increment is increment. If first and - ! increment have rank greater than zero, returns an array of one larger rank, - ! with the last subscript having size n and indexing the progressions. - FUNCTION arth_sp(first, increment, n) - - IMPLICIT NONE - - REAL(sp), INTENT(IN) :: first, increment - INTEGER(i4), INTENT(IN) :: n - REAL(sp), DIMENSION(n) :: arth_sp - - INTEGER(i4) :: k, k2 - REAL(sp) :: temp - - if (n > 0) arth_sp(1) = first - if (n <= NPAR_ARTH) then - do k = 2, n - arth_sp(k) = arth_sp(k - 1) + increment - end do - else - do k = 2, NPAR2_ARTH - arth_sp(k) = arth_sp(k - 1) + increment - end do - temp = increment * NPAR2_ARTH - k = NPAR2_ARTH - do - if (k >= n) exit - k2 = k + k - arth_sp(k + 1 : min(k2, n)) = temp + arth_sp(1 : min(k, n - k)) - temp = temp + temp - k = k2 - end do - end if - - END FUNCTION arth_sp - - - FUNCTION arth_dp(first, increment, n) - - IMPLICIT NONE - - REAL(dp), INTENT(IN) :: first, increment - INTEGER(i4), INTENT(IN) :: n - REAL(dp), DIMENSION(n) :: arth_dp - - INTEGER(i4) :: k, k2 - REAL(dp) :: temp - - if (n > 0) arth_dp(1) = first - if (n <= NPAR_ARTH) then - do k = 2, n - arth_dp(k) = arth_dp(k - 1) + increment - end do - else - do k = 2, NPAR2_ARTH - arth_dp(k) = arth_dp(k - 1) + increment - end do - temp = increment * NPAR2_ARTH - k = NPAR2_ARTH - do - if (k >= n) exit - k2 = k + k - arth_dp(k + 1 : min(k2, n)) = temp + arth_dp(1 : min(k, n - k)) - temp = temp + temp - k = k2 - end do - end if - - END FUNCTION arth_dp - + FUNCTION autocorr_dp(x, k, mask) result(acf) - FUNCTION arth_i4(first, increment, n) - - IMPLICIT NONE - - INTEGER(i4), INTENT(IN) :: first, increment, n - INTEGER(i4), DIMENSION(n) :: arth_i4 - - INTEGER(i4) :: k, k2, temp - - if (n > 0) arth_i4(1) = first - if (n <= NPAR_ARTH) then - do k = 2, n - arth_i4(k) = arth_i4(k - 1) + increment - end do - else - do k = 2, NPAR2_ARTH - arth_i4(k) = arth_i4(k - 1) + increment - end do - temp = increment * NPAR2_ARTH - k = NPAR2_ARTH - do - if (k >= n) exit - k2 = k + k - arth_i4(k + 1 : min(k2, n)) = temp + arth_i4(1 : min(k, n - k)) - temp = temp + temp - k = k2 - end do - end if - - END FUNCTION arth_i4 - - ! ------------------------------------------------------------------ - - FUNCTION autocoeffk_dp(x, k, mask) + USE mo_moment, ONLY: moment IMPLICIT NONE REAL(dp), DIMENSION(:), INTENT(IN) :: x INTEGER(i4), INTENT(IN) :: k LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: autocoeffk_dp + REAL(dp) :: acf + + INTEGER(i4) :: kk ! absolute value of lag k + INTEGER(i4) :: nn ! number of true values in mask + REAL(dp) :: n ! real of nn + REAL(dp) :: mean, var ! moments + REAL(dp), DIMENSION(size(x)) :: x_shift + LOGICAL, DIMENSION(size(x)) :: maske + maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocoeffk_dp: size(mask) /= size(x)' - autocoeffk_dp = crosscoeffk(x, x, k, mask) - else - autocoeffk_dp = crosscoeffk(x, x, k) + if (size(mask) /= size(x)) stop 'Error autocorr_dp: size(mask) /= size(x)' + maske = mask end if - END FUNCTION autocoeffk_dp + ! calculate mean and population variance of x + call moment(x, mean=mean, variance=var, mask=maske, sample=.false.) + ! shift x to 0 mean + x_shift = x - mean + ! use absolute value of k + kk = abs(k) + nn = size(x) + n = real(count(maske(1 : nn - kk).and.maske(1 + kk : nn)), dp) + ! calculate the auto-correlation function + acf = sum(x_shift(1 : nn - kk) * x_shift(1 + kk : nn), mask = (maske(1 : nn - kk) .and. maske(1 + kk : nn))) / n / var + + END FUNCTION autocorr_dp + FUNCTION autocorr_sp(x, k, mask) result(acf) - FUNCTION autocoeffk_sp(x, k, mask) + USE mo_moment, ONLY: moment IMPLICIT NONE REAL(sp), DIMENSION(:), INTENT(IN) :: x INTEGER(i4), INTENT(IN) :: k LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: autocoeffk_sp + REAL(sp) :: acf + + INTEGER(i4) :: kk ! absolute value of lag k + INTEGER(i4) :: nn ! number of true values in mask + REAL(sp) :: n ! real of nn + REAL(sp) :: mean, var ! moments + REAL(sp), DIMENSION(size(x)) :: x_shift + LOGICAL, DIMENSION(size(x)) :: maske + maske(:) = .true. if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocoeffk_sp: size(mask) /= size(x)' - autocoeffk_sp = crosscoeffk(x, x, k, mask) - else - autocoeffk_sp = crosscoeffk(x, x, k) + if (size(mask) /= size(x)) stop 'Error autocorr_sp: size(mask) /= size(x)' + maske = mask end if - END FUNCTION autocoeffk_sp + ! calculate mean and population variance of x + call moment(x, mean=mean, variance=var, mask=maske, sample=.false.) + ! shift x to 0 mean + x_shift = x - mean + ! use absolute value of k + kk = abs(k) + nn = size(x) + n = real(count(maske(1 : nn - kk).and.maske(1 + kk : nn)), sp) + ! calculate the auto-correlation function + acf = sum(x_shift(1 : nn - kk) * x_shift(1 + kk : nn), mask = (maske(1 : nn - kk) .and. maske(1 + kk : nn))) / n / var + + END FUNCTION autocorr_sp - FUNCTION autocoeffk_1d_dp(x, k, mask) result(acf) + FUNCTION autocorr_1d_dp(x, k, mask) result(acf) IMPLICIT NONE @@ -565,20 +174,19 @@ FUNCTION autocoeffk_1d_dp(x, k, mask) result(acf) REAL(dp), DIMENSION(size(k)) :: acf if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocoeffk_1d_dp: size(mask) /= size(x)' do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i), mask) + acf(i) = autocorr(x, k(i), mask) end do else do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i)) + acf(i) = autocorr(x, k(i)) end do end if - END FUNCTION autocoeffk_1d_dp + END FUNCTION autocorr_1d_dp - FUNCTION autocoeffk_1d_sp(x, k, mask) result(acf) + FUNCTION autocorr_1d_sp(x, k, mask) result(acf) IMPLICIT NONE @@ -589,1061 +197,15 @@ FUNCTION autocoeffk_1d_sp(x, k, mask) result(acf) REAL(sp), DIMENSION(size(k)) :: acf if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocoeffk_1d_sp: size(mask) /= size(x)' do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i), mask) + acf(i) = autocorr(x, k(i), mask) end do else do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i)) - end do - end if - - END FUNCTION autocoeffk_1d_sp - - ! ------------------------------------------------------------------ - - FUNCTION autocorr_dp(x, k, mask) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: autocorr_dp - - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocorr_1d_dp: size(mask) /= size(x)' - autocorr_dp = crosscoeffk(x, x, k, mask) / crosscoeffk(x, x, 0, mask) - else - autocorr_dp = crosscoeffk(x, x, k) / crosscoeffk(x, x, 0) - end if - - END FUNCTION autocorr_dp - - - FUNCTION autocorr_sp(x, k, mask) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: autocorr_sp - - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocorr_1d_sp: size(mask) /= size(x)' - autocorr_sp = crosscoeffk(x, x, k, mask) / crosscoeffk(x, x, 0, mask) - else - autocorr_sp = crosscoeffk(x, x, k) / crosscoeffk(x, x, 0) - end if - - END FUNCTION autocorr_sp - - FUNCTION autocorr_1d_dp(x, k, mask) result(acf) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), DIMENSION(:), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - INTEGER(i4) :: i - REAL(dp), DIMENSION(size(k)) :: acf - REAL(dp) :: c0 - - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocorr_dp: size(mask) /= size(x)' - c0 = crosscoeffk(x, x, 0, mask) - do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i), mask) / c0 - end do - else - c0 = crosscoeffk(x, x, 0) - do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i)) / c0 - end do - end if - - END FUNCTION autocorr_1d_dp - - FUNCTION autocorr_1d_sp(x, k, mask) result(acf) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(IN) :: x - INTEGER(i4), DIMENSION(:), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - INTEGER(i4) :: i - REAL(sp), DIMENSION(size(k)) :: acf - REAL(sp) :: c0 - - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error autocorr_sp: size(mask) /= size(x)' - c0 = crosscoeffk(x, x, 0, mask) - do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i), mask) / c0 - end do - else - c0 = crosscoeffk(x, x, 0) - do i = 1, size(k) - acf(i) = crosscoeffk(x, x, k(i)) / c0 + acf(i) = autocorr(x, k(i)) end do end if END FUNCTION autocorr_1d_sp - ! ------------------------------------------------------------------ - - FUNCTION corr_dp(data1, data2, nadjust, nhigh, nwin) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(IN) :: data1, data2 - INTEGER(i4), OPTIONAL, INTENT(OUT) :: nadjust - INTEGER(i4), OPTIONAL, INTENT(IN) :: nhigh - INTEGER(i4), OPTIONAL, INTENT(IN) :: nwin - REAL(dp), DIMENSION(size(data1)) :: corr_dp - - REAL(dp), DIMENSION(:), ALLOCATABLE :: dat1, dat2, corrout - COMPLEX(dpc), DIMENSION(:), ALLOCATABLE :: cdat1, cdat2 - REAL(dp), DIMENSION(:), ALLOCATABLE :: win1 - !COMPLEX(dpc), DIMENSION(:), ALLOCATABLE :: cwin1 - INTEGER(i4) :: i, n, no2, iwin, ihigh - REAL(dp) :: no2_1, ihigh1 - - n = size(data1) - if (size(data2) /= n) stop 'Error corr_dp: size(data1) /= size(data2)' - - if (present(nwin)) then - iwin = nwin - else - iwin = 1 - end if - - if (present(nhigh)) then - ihigh = nhigh - else - ihigh = 0 - end if - - if (iand(n, n - 1) /= 0) then - if (present(nadjust)) then - n = 2**floor(log(real(n, dp)) / log(2.0_dp)) - nadjust = n - else - stop 'Error corr_dp: size(data1) must be a power of 2' - end if - else - if (present(nadjust)) then - nadjust = n - end if - end if - - allocate(dat1(n)) - allocate(dat2(n)) - dat1 = data1(1 : n) - dat2 = data2(1 : n) - - no2 = n / 2 - no2_1 = 1.0_dp / real(no2, dp) - allocate(cdat1(no2)) - allocate(cdat2(no2)) - allocate(corrout(n)) - - ! FFT - call realft(dat1, 1, cdat1) - call realft(dat2, 1, cdat2) - - ! Highpass - if (ihigh > 0) then - ! FxH - allocate(win1(no2)) - !allocate(cwin1(no2)) - select case(iwin) - case(0) ! no window - win1(1 : no2) = 1.0_dp - case(1) ! ideal high pass filter - win1(1 : ihigh) = 0.0_dp - win1(ihigh + 1 : no2) = 1.0_dp - case(2) ! similar Bartlett window - ihigh1 = 1.0_dp / real(ihigh, dp) - forall(i = 1 : ihigh) win1(i) = real(i - 1, dp) * ihigh1 - win1(ihigh + 1 : no2) = 1.0_dp - case default - stop 'Unimplemented window option in corr_dp' - end select - !cwin1 = cmplx(win1, win1, kind=dpc) - ! low pass - ! cdat1(1:no2) = cdat1(1:no2)*cwin1(1:no2) - ! cdat2(1:no2) = cdat2(1:no2)*cwin1(1:no2) - cdat1(1 : no2) = cdat1(1 : no2) * win1(1 : no2) - cdat2(1 : no2) = cdat2(1 : no2) * win1(1 : no2) - end if - - ! FxF* - cdat1(1) = cmplx(real(cdat1(1)) * real(cdat2(1)) * no2_1, & - aimag(cdat1(1)) * aimag(cdat2(1)) * no2_1, kind = dpc) - cdat1(2 :) = cdat1(2 :) * conjg(cdat2(2 :)) * no2_1 - - ! IFFT - call realft(corrout, -1, cdat1) - corr_dp(1 : n) = corrout(1 : n) - if (size(corr_dp) > n) corr_dp(n + 1 :) = 0.0_dp - - deallocate(dat1) - deallocate(dat2) - deallocate(cdat1) - deallocate(cdat2) - deallocate(corrout) - if (ihigh > 0) then - deallocate(win1) - !deallocate(cwin1) - end if - - END FUNCTION corr_dp - - - FUNCTION corr_sp(data1, data2, nadjust, nhigh, nwin) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(IN) :: data1, data2 - INTEGER(i4), OPTIONAL, INTENT(OUT) :: nadjust - INTEGER(i4), OPTIONAL, INTENT(IN) :: nhigh - INTEGER(i4), OPTIONAL, INTENT(IN) :: nwin - REAL(sp), DIMENSION(size(data1)) :: corr_sp - - REAL(sp), DIMENSION(:), ALLOCATABLE :: dat1, dat2, corrout - COMPLEX(spc), DIMENSION(:), ALLOCATABLE :: cdat1, cdat2 - REAL(sp), DIMENSION(:), ALLOCATABLE :: win1 - !COMPLEX(spc), DIMENSION(:), ALLOCATABLE :: cwin1 - INTEGER(i4) :: i, n, no2, iwin, ihigh - REAL(sp) :: no2_1, ihigh1 - - n = size(data1) - if (size(data2) /= n) stop 'Error corr_sp: size(data1) /= size(data2)' - - if (present(nwin)) then - iwin = nwin - else - iwin = 1 - end if - - if (present(nhigh)) then - ihigh = nhigh - else - ihigh = 0 - end if - - if (iand(n, n - 1) /= 0) then - if (present(nadjust)) then - n = 2**floor(log(real(n, sp)) / log(2.0_sp)) - nadjust = n - else - stop 'Error corr_sp: size(data1) must be a power of 2' - end if - else - if (present(nadjust)) then - nadjust = n - end if - end if - - allocate(dat1(n)) - allocate(dat2(n)) - dat1 = data1(1 : n) - dat2 = data2(1 : n) - - no2 = n / 2 - no2_1 = 1.0_sp / real(no2, sp) - allocate(cdat1(no2)) - allocate(cdat2(no2)) - allocate(corrout(n)) - - ! FFT - call realft(dat1, 1, cdat1) - call realft(dat2, 1, cdat2) - - ! Highpass - if (ihigh > 0) then - ! FxH - allocate(win1(no2)) - !allocate(cwin1(no2)) - select case(iwin) - case(0) ! no window - win1(1 : no2) = 1.0_sp - case(1) ! ideal high pass filter - win1(1 : ihigh) = 0.0_sp - win1(ihigh + 1 : no2) = 1.0_sp - case(2) ! similar Bartlett window - ihigh1 = 1.0_sp / real(ihigh, sp) - forall(i = 1 : ihigh) win1(i) = real(i - 1, sp) * ihigh1 - win1(ihigh + 1 : no2) = 1.0_sp - case default - stop 'Unimplemented window option in corr_sp' - end select - !cwin1 = cmplx(win1, win1, kind=spc) - ! low pass - ! cdat1(1:no2) = cdat1(1:no2)*cwin1(1:no2) - ! cdat2(1:no2) = cdat2(1:no2)*cwin1(1:no2) - cdat1(1 : no2) = cdat1(1 : no2) * win1(1 : no2) - cdat2(1 : no2) = cdat2(1 : no2) * win1(1 : no2) - end if - - ! FxF* - cdat1(1) = cmplx(real(cdat1(1)) * real(cdat2(1)) * no2_1, & - aimag(cdat1(1)) * aimag(cdat2(1)) * no2_1, kind = spc) - cdat1(2 :) = cdat1(2 :) * conjg(cdat2(2 :)) * no2_1 - - ! IFFT - call realft(corrout, -1, cdat1) - corr_sp(1 : n) = corrout(1 : n) - if (size(corr_sp) > n) corr_sp(n + 1 :) = 0.0_sp - - deallocate(dat1) - deallocate(dat2) - deallocate(cdat1) - deallocate(cdat2) - deallocate(corrout) - if (ihigh > 0) then - deallocate(win1) - !deallocate(cwin1) - end if - - END FUNCTION corr_sp - - ! ------------------------------------------------------------------ - - FUNCTION crosscoeffk_dp(x, y, k, mask) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: crosscoeffk_dp - - INTEGER(i4) :: nn ! number of true values in mask - INTEGER(i4) :: nnn ! number of true values in mask .and. shifted mask by lag k - REAL(dp) :: n ! real of nn or nnn - INTEGER(i4) :: kk ! absolute value of lag k - REAL(dp) :: ave - - REAL(dp), DIMENSION(size(x)) :: xdash - REAL(dp), DIMENSION(size(x)) :: ydash - LOGICAL, DIMENSION(size(x)) :: maske - - maske(:) = .true. - if (present(mask)) then - if (size(x) /= size(y)) stop 'Error crosscoeffk_dp: size(x) /= size(y)' - if (size(mask) /= size(x)) stop 'Error crosscoeffk_dp: size(mask) /= size(x)' - maske = mask - end if - - nn = count(maske) - n = real(nn, dp) - ! crosscoeffk(x, y, k) = crosscoeffk(y, x, -k) - if (k >= 0) then - ave = sum(x(:), mask = maske) / n - xdash = x - ave - ave = sum(y(:), mask = maske) / n - ydash = y - ave - else - ave = sum(y(:), mask = maske) / n - xdash = y - ave - ave = sum(x(:), mask = maske) / n - ydash = x - ave - end if - kk = abs(k) - nnn = size(x, 1) - n = real(count(maske(1 : nnn - kk).and.maske(1 + kk : nnn)), dp) - crosscoeffk_dp = sum(xdash(1 : nnn - kk) * ydash(1 + kk : nnn), mask = (maske(1 : nnn - kk).and.maske(1 + kk : nnn))) / n - - END FUNCTION crosscoeffk_dp - - FUNCTION crosscoeffk_sp(x, y, k, mask) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: crosscoeffk_sp - - INTEGER(i4) :: nn ! number of true values in mask - INTEGER(i4) :: nnn ! number of true values in mask .and. shifted mask by lag k - REAL(sp) :: n ! real of nn or nnn - INTEGER(i4) :: kk ! absolute value of lag k - REAL(sp) :: ave - - REAL(sp), DIMENSION(size(x)) :: xdash - REAL(sp), DIMENSION(size(x)) :: ydash - LOGICAL, DIMENSION(size(x)) :: maske - - maske(:) = .true. - if (present(mask)) then - if (size(x) /= size(y)) stop 'Error crosscoeffk_sp: size(x) /= size(y)' - if (size(mask) /= size(x)) stop 'Error crosscoeffk_sp: size(mask) /= size(x)' - maske = mask - end if - - nn = count(maske) - n = real(nn, sp) - ! crosscoeffk(x, y, k) = crosscoeffk(y, x, -k) - if (k >= 0) then - ave = sum(x(:), mask = maske) / n - xdash = x - ave - ave = sum(y(:), mask = maske) / n - ydash = y - ave - else - ave = sum(y(:), mask = maske) / n - xdash = y - ave - ave = sum(x(:), mask = maske) / n - ydash = x - ave - end if - kk = abs(k) - nnn = size(x, 1) - n = real(count(maske(1 : nnn - kk).and.maske(1 + kk : nnn)), sp) - crosscoeffk_sp = sum(xdash(1 : nnn - kk) * ydash(1 + kk : nnn), mask = (maske(1 : nnn - kk).and.maske(1 + kk : nnn))) / n - - END FUNCTION crosscoeffk_sp - - ! ------------------------------------------------------------------ - - FUNCTION crosscorr_dp(x, y, k, mask) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(IN) :: x - REAL(dp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(dp) :: crosscorr_dp - - if (size(x) /= size(y)) stop 'Error crosscorr_dp: size(x) /= size(y)' - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error crosscorr_dp: size(mask) /= size(x)' - crosscorr_dp = crosscoeffk(x, y, k, mask) / crosscoeffk(x, y, 0, mask) - else - crosscorr_dp = crosscoeffk(x, y, k) / crosscoeffk(x, y, 0) - end if - - END FUNCTION crosscorr_dp - - FUNCTION crosscorr_sp(x, y, k, mask) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(IN) :: x - REAL(sp), DIMENSION(:), INTENT(IN) :: y - INTEGER(i4), INTENT(IN) :: k - LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask - REAL(sp) :: crosscorr_sp - - if (size(x) /= size(y)) stop 'Error crosscorr_sp: size(x) /= size(y)' - if (present(mask)) then - if (size(mask) /= size(x)) stop 'Error crosscorr_sp: size(mask) /= size(x)' - crosscorr_sp = crosscoeffk(x, y, k, mask) / crosscoeffk(x, y, 0, mask) - else - crosscorr_sp = crosscoeffk(x, y, k) / crosscoeffk(x, y, 0) - end if - - END FUNCTION crosscorr_sp - - ! ------------------------------------------------------------------ - - ! From numerical recipes documentation - ! Replaces a complex array data by its discrete Fourier transform, if isign is input as 1; - ! or replaces data by its inverse discrete Fourier transform times the size of data, if isign - ! is input as -1. The size of data must be an integer power of 2. Parallelismis achieved - ! by internally reshaping the input array to two dimensions. (Use this version if fourrow is - ! faster than fourcol on your machine.) - - SUBROUTINE four1_sp(data, isign) - - IMPLICIT NONE - - COMPLEX(spc), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - - COMPLEX(spc), DIMENSION(:, :), ALLOCATABLE :: dat, temp - COMPLEX(dpc), DIMENSION(:), ALLOCATABLE :: w, wp - REAL(dp), DIMENSION(:), ALLOCATABLE :: theta - INTEGER(i4) :: n, m1, m2, j - - n = size(data) - if (iand(n, n - 1) /= 0) stop 'Error four1_sp: size(data1) must be a power of 2' - - m1 = 2**ceiling(0.5_sp * log(real(n, sp)) / 0.693147_sp) - m2 = n / m1 - allocate(dat(m1, m2), theta(m1), w(m1), wp(m1), temp(m2, m1)) - dat = reshape(data, shape(dat)) - call fourrow(dat, isign) - theta = arth(0, isign, m1) * TWOPI_dp / real(n, dp) - wp = cmplx(-2.0_dp * sin(0.5_dp * theta)**2, sin(theta), kind = dpc) - w = cmplx(1.0_dp, 0.0_dp, kind = dpc) - do j = 2, m2 - w = w * wp + w - dat(:, j) = dat(:, j) * cmplx(w, kind = spc) - end do - temp = transpose(dat) - call fourrow(temp, isign) - data = reshape(temp, shape(data)) - deallocate(dat, w, wp, theta, temp) - - END SUBROUTINE four1_sp - - - SUBROUTINE four1_dp(data, isign) - - IMPLICIT NONE - - COMPLEX(dpc), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - - COMPLEX(dpc), DIMENSION(:, :), ALLOCATABLE :: dat, temp - COMPLEX(dpc), DIMENSION(:), ALLOCATABLE :: w, wp - REAL(dp), DIMENSION(:), ALLOCATABLE :: theta - INTEGER(i4) :: n, m1, m2, j - - n = size(data) - if (iand(n, n - 1) /= 0) stop 'Error four1_dp: size(data1) must be a power of 2' - - m1 = 2**ceiling(0.5_sp * log(real(n, sp)) / 0.693147_sp) - m2 = n / m1 - allocate(dat(m1, m2), theta(m1), w(m1), wp(m1), temp(m2, m1)) - dat = reshape(data, shape(dat)) - call fourrow(dat, isign) - theta = arth(0, isign, m1) * TWOPI_dp / real(n, dp) - wp = cmplx(-2.0_dp * sin(0.5_dp * theta)**2, sin(theta), kind = dpc) - w = cmplx(1.0_dp, 0.0_dp, kind = dpc) - do j = 2, m2 - w = w * wp + w - dat(:, j) = dat(:, j) * w - end do - temp = transpose(dat) - call fourrow(temp, isign) - data = reshape(temp, shape(data)) - deallocate(dat, w, wp, theta, temp) - - END SUBROUTINE four1_dp - - ! ------------------------------------------------------------------ - - ! From numerical recipes documentation - ! Replaces each row (constant first index) of data(1:M,1:N) by its discrete Fourier trans- - ! form (transform on second index), if isign is input as 1; or replaces each row of data - ! by N times its inverse discrete Fourier transform, if isign is input as -1. N must be an - ! integer power of 2. Parallelism is M-fold on the first index of data. - - SUBROUTINE fourrow_sp(data, isign) - - IMPLICIT NONE - - COMPLEX(spc), DIMENSION(:, :), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - - INTEGER(i4) :: n, i, istep, j, m, mmax, n2 - REAL(dp) :: theta - COMPLEX(spc), DIMENSION(size(data, 1)) :: temp - COMPLEX(dpc) :: w, wp - COMPLEX(spc) :: ws - - n = size(data, 2) - if (iand(n, n - 1) /= 0) stop 'Error fourrow_sp: size(data,2) must be a power of 2' - n2 = n / 2 - j = n2 - do i = 1, n - 2 - if (j > i) call swap(data(:, j + 1), data(:, i + 1)) - m = n2 - do - if (m < 2 .or. j < m) exit - j = j - m - m = m / 2 - end do - j = j + m - end do - mmax = 1 - do - if (n <= mmax) exit - istep = 2 * mmax - theta = PI_dp / real(isign * mmax, dp) - wp = cmplx(-2.0_dp * sin(0.5_dp * theta)**2, sin(theta), kind = dpc) - w = cmplx(1.0_dp, 0.0_dp, kind = dpc) - do m = 1, mmax - ws = cmplx(w, kind = spc) - do i = m, n, istep - j = i + mmax - temp = ws * data(:, j) - data(:, j) = data(:, i) - temp - data(:, i) = data(:, i) + temp - end do - w = w * wp + w - end do - mmax = istep - end do - - END SUBROUTINE fourrow_sp - - - SUBROUTINE fourrow_dp(data, isign) - - IMPLICIT NONE - - COMPLEX(dpc), DIMENSION(:, :), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - - INTEGER(i4) :: n, i, istep, j, m, mmax, n2 - REAL(dp) :: theta - COMPLEX(dpc), DIMENSION(size(data, 1)) :: temp - COMPLEX(dpc) :: w, wp - COMPLEX(dpc) :: ws - - n = size(data, 2) - if (iand(n, n - 1) /= 0) stop 'Error fourrow_dp: size(data,2) must be a power of 2' - n2 = n / 2 - j = n2 - do i = 1, n - 2 - if (j > i) call swap(data(:, j + 1), data(:, i + 1)) - m = n2 - do - if (m < 2 .or. j < m) exit - j = j - m - m = m / 2 - end do - j = j + m - end do - mmax = 1 - do - if (n <= mmax) exit - istep = 2 * mmax - theta = PI_dp / real(isign * mmax, dp) - wp = cmplx(-2.0_dp * sin(0.5_dp * theta)**2, sin(theta), kind = dpc) - w = cmplx(1.0_dp, 0.0_dp, kind = dpc) - do m = 1, mmax - ws = w - do i = m, n, istep - j = i + mmax - temp = ws * data(:, j) - data(:, j) = data(:, i) - temp - data(:, i) = data(:, i) + temp - end do - w = w * wp + w - end do - mmax = istep - end do - - END SUBROUTINE fourrow_dp - - ! ------------------------------------------------------------------ - - ! From numerical recipes documentation - ! When isign=1, calculates the Fourier transform of a set of N real-valued data points, - ! input in the array data. If the optional argument zdata is not present, the data are replaced - ! by the positive frequency half of its complex Fourier transform. There al-valued first and - ! last components of the complex transform are returned as elements data(1) and data(2), - ! respectively. If the complex array zdata of lengthN/2 ispresent, data is unchanged and - ! the transform is returned in zdata. N must be a power of 2. If isign=-1, this routine - ! calculates the inverse transform of a complex data array if it is the transform of real data. - ! (Result in this case must be multiplied by 2/N.) The data can be supplied either in data, - ! with zdata absent, or inzdata. - SUBROUTINE realft_dp(data, isign, zdata) - - IMPLICIT NONE - - REAL(dp), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - COMPLEX(dpc), DIMENSION(:), OPTIONAL, TARGET :: zdata - - INTEGER(i4) :: n, nh, nq - COMPLEX(dpc), DIMENSION(size(data) / 4) :: w - COMPLEX(dpc), DIMENSION(size(data) / 4 - 1) :: h1, h2 - COMPLEX(dpc), DIMENSION(:), POINTER :: cdata - COMPLEX(dpc) :: z - REAL(dp) :: c1 = 0.5_dp, c2 - - n = size(data) - if (iand(n, n - 1) /= 0) stop 'Error realft_dp: size(data) must be a power of 2' - nh = n / 2 - nq = n / 4 - if (present(zdata)) then - if (n / 2 /= size(zdata)) stop 'Error realft_dp: size(zdata) /= size(data)/2' - cdata => zdata - if (isign == 1) cdata = cmplx(data(1 : n - 1 : 2), data(2 : n : 2), kind = dpc) - else - allocate(cdata(n / 2)) - cdata = cmplx(data(1 : n - 1 : 2), data(2 : n : 2), kind = dpc) - end if - if (isign == 1) then - c2 = -0.5_dp - call four1(cdata, +1) - else - c2 = 0.5_dp - end if - w = zroots_unity_dp(sign(n, isign), n / 4) - w = cmplx(-aimag(w), real(w), kind = dpc) - h1 = c1 * (cdata(2 : nq) + conjg(cdata(nh : nq + 2 : -1))) - h2 = c2 * (cdata(2 : nq) - conjg(cdata(nh : nq + 2 : -1))) - cdata(2 : nq) = h1 + w(2 : nq) * h2 - cdata(nh : nq + 2 : -1) = conjg(h1 - w(2 : nq) * h2) - z = cdata(1) - if (isign == 1) then - cdata(1) = cmplx(real(z) + aimag(z), real(z) - aimag(z), kind = dpc) - else - cdata(1) = cmplx(c1 * (real(z) + aimag(z)), c1 * (real(z) - aimag(z)), kind = dpc) - call four1(cdata, -1) - end if - if (present(zdata)) then - if (isign /= 1) then - data(1 : n - 1 : 2) = real(cdata) - data(2 : n : 2) = aimag(cdata) - end if - else - data(1 : n - 1 : 2) = real(cdata) - data(2 : n : 2) = aimag(cdata) - deallocate(cdata) - end if - - END SUBROUTINE realft_dp - - - SUBROUTINE realft_sp(data, isign, zdata) - - IMPLICIT NONE - - REAL(sp), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(i4), INTENT(IN) :: isign - COMPLEX(spc), DIMENSION(:), OPTIONAL, TARGET :: zdata - - INTEGER(i4) :: n, nh, nq - COMPLEX(spc), DIMENSION(size(data) / 4) :: w - COMPLEX(spc), DIMENSION(size(data) / 4 - 1) :: h1, h2 - COMPLEX(spc), DIMENSION(:), POINTER :: cdata - COMPLEX(spc) :: z - REAL(sp) :: c1 = 0.5_sp, c2 - - n = size(data) - if (iand(n, n - 1) /= 0) stop 'Error realft_sp: size(data) must be a power of 2' - nh = n / 2 - nq = n / 4 - if (present(zdata)) then - if (n / 2 /= size(zdata)) stop 'Error realft_sp: size(zdata) /= size(data)/2' - cdata => zdata - if (isign == 1) cdata = cmplx(data(1 : n - 1 : 2), data(2 : n : 2), kind = spc) - else - allocate(cdata(n / 2)) - cdata = cmplx(data(1 : n - 1 : 2), data(2 : n : 2), kind = spc) - end if - if (isign == 1) then - c2 = -0.5_sp - call four1(cdata, +1) - else - c2 = 0.5_sp - end if - w = zroots_unity_sp(sign(n, isign), n / 4) - w = cmplx(-aimag(w), real(w), kind = spc) - h1 = c1 * (cdata(2 : nq) + conjg(cdata(nh : nq + 2 : -1))) - h2 = c2 * (cdata(2 : nq) - conjg(cdata(nh : nq + 2 : -1))) - cdata(2 : nq) = h1 + w(2 : nq) * h2 - cdata(nh : nq + 2 : -1) = conjg(h1 - w(2 : nq) * h2) - z = cdata(1) - if (isign == 1) then - cdata(1) = cmplx(real(z) + aimag(z), real(z) - aimag(z), kind = spc) - else - cdata(1) = cmplx(c1 * (real(z) + aimag(z)), c1 * (real(z) - aimag(z)), kind = spc) - call four1(cdata, -1) - end if - if (present(zdata)) then - if (isign /= 1) then - data(1 : n - 1 : 2) = real(cdata) - data(2 : n : 2) = aimag(cdata) - end if - else - data(1 : n - 1 : 2) = real(cdata) - data(2 : n : 2) = aimag(cdata) - deallocate(cdata) - end if - - END SUBROUTINE realft_sp - - ! ------------------------------------------------------------------ - - ! From numerical recipes documentation - ! Swaps the corresponding elements of a and b. If mask is present, performs - ! the swap only where mask is true. (Following code is the unmasked case. - ! For speed at runtime, the masked case is implemented by overloading, not - ! by testing for the optional argument.) - ! SUBROUTINE swap_i4(a,b) - ! INTEGER(i4), INTENT(INOUT) :: a,b - ! INTEGER(i4) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_i4 - - - ! SUBROUTINE swap_sp(a,b) - ! REAL(sp), INTENT(INOUT) :: a,b - ! REAL(sp) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_sp - - - ! SUBROUTINE masked_swap_sp(a,b,mask) - ! REAL(sp), INTENT(INOUT) :: a,b - ! LOGICAL, INTENT(IN) :: mask - ! REAL(sp) :: swp - ! if (mask) then - ! swp=a - ! a=b - ! b=swp - ! end if - ! END SUBROUTINE masked_swap_sp - - - ! SUBROUTINE masked_swap_1d_sp(a,b,mask) - ! REAL(sp), DIMENSION(:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:), INTENT(IN) :: mask - ! REAL(sp), DIMENSION(size(a)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_1d_sp - - - ! SUBROUTINE masked_swap_2d_sp(a,b,mask) - ! REAL(sp), DIMENSION(:,:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:,:), INTENT(IN) :: mask - ! REAL(sp), DIMENSION(size(a,1),size(a,2)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_2d_sp - - - ! SUBROUTINE swap_1d_sp(a,b) - ! REAL(sp), DIMENSION(:), INTENT(INOUT) :: a,b - ! REAL(sp), DIMENSION(SIZE(a)) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_1d_sp - - - ! SUBROUTINE swap_dp(a,b) - ! REAL(dp), INTENT(INOUT) :: a,b - ! REAL(dp) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_dp - - - ! SUBROUTINE swap_1d_dp(a,b) - ! REAL(dp), DIMENSION(:), INTENT(INOUT) :: a,b - ! REAL(dp), DIMENSION(SIZE(a)) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_1d_dp - - - ! SUBROUTINE masked_swap_dp(a,b,mask) - ! REAL(dp), INTENT(INOUT) :: a,b - ! LOGICAL, INTENT(IN) :: mask - ! REAL(dp) :: swp - ! if (mask) then - ! swp=a - ! a=b - ! b=swp - ! end if - ! END SUBROUTINE masked_swap_dp - - - ! SUBROUTINE masked_swap_1d_dp(a,b,mask) - ! REAL(dp), DIMENSION(:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:), INTENT(IN) :: mask - ! REAL(dp), DIMENSION(size(a)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_1d_dp - - - ! SUBROUTINE masked_swap_2d_dp(a,b,mask) - ! REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:,:), INTENT(IN) :: mask - ! REAL(dp), DIMENSION(size(a,1),size(a,2)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_2d_dp - - - ! SUBROUTINE swap_spc(a,b) - ! COMPLEX(spc), INTENT(INOUT) :: a,b - ! COMPLEX(spc) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_spc - - - ! SUBROUTINE masked_swap_spc(a,b,mask) - ! COMPLEX(spc), INTENT(INOUT) :: a,b - ! LOGICAL, INTENT(IN) :: mask - ! COMPLEX(spc) :: swp - ! if (mask) then - ! swp=a - ! a=b - ! b=swp - ! end if - ! END SUBROUTINE masked_swap_spc - - - ! SUBROUTINE masked_swap_1d_spc(a,b,mask) - ! COMPLEX(spc), DIMENSION(:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:), INTENT(IN) :: mask - ! COMPLEX(spc), DIMENSION(size(a)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_1d_spc - - - ! SUBROUTINE masked_swap_2d_spc(a,b,mask) - ! COMPLEX(spc), DIMENSION(:,:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:,:), INTENT(IN) :: mask - ! COMPLEX(spc), DIMENSION(size(a,1),size(a,2)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_2d_spc - - - SUBROUTINE swap_1d_spc(a, b) - COMPLEX(spc), DIMENSION(:), INTENT(INOUT) :: a, b - COMPLEX(spc), DIMENSION(SIZE(a)) :: dum - dum = a - a = b - b = dum - END SUBROUTINE swap_1d_spc - - - ! SUBROUTINE swap_dpc(a,b) - ! COMPLEX(dpc), INTENT(INOUT) :: a,b - ! COMPLEX(dpc) :: dum - ! dum=a - ! a=b - ! b=dum - ! END SUBROUTINE swap_dpc - - - ! SUBROUTINE masked_swap_dpc(a,b,mask) - ! COMPLEX(dpc), INTENT(INOUT) :: a,b - ! LOGICAL, INTENT(IN) :: mask - ! COMPLEX(dpc) :: swp - ! if (mask) then - ! swp=a - ! a=b - ! b=swp - ! end if - ! END SUBROUTINE masked_swap_dpc - - - ! SUBROUTINE masked_swap_1d_dpc(a,b,mask) - ! COMPLEX(dpc), DIMENSION(:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:), INTENT(IN) :: mask - ! COMPLEX(dpc), DIMENSION(size(a)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_1d_dpc - - - ! SUBROUTINE masked_swap_2d_dpc(a,b,mask) - ! COMPLEX(dpc), DIMENSION(:,:), INTENT(INOUT) :: a,b - ! LOGICAL, DIMENSION(:,:), INTENT(IN) :: mask - ! COMPLEX(dpc), DIMENSION(size(a,1),size(a,2)) :: swp - ! where (mask) - ! swp=a - ! a=b - ! b=swp - ! end where - ! END SUBROUTINE masked_swap_2d_dpc - - - SUBROUTINE swap_1d_dpc(a, b) - COMPLEX(dpc), DIMENSION(:), INTENT(INOUT) :: a, b - COMPLEX(dpc), DIMENSION(SIZE(a)) :: dum - dum = a - a = b - b = dum - END SUBROUTINE swap_1d_dpc - - ! ------------------------------------------------------------------ - - ! From numerical recipes documentation - ! Returns a complex array containing nn consecutive powers of the nth - ! complex root of unity. - FUNCTION zroots_unity_dp(n, nn) - - INTEGER(i4), INTENT(IN) :: n, nn - COMPLEX(dpc), DIMENSION(nn) :: zroots_unity_dp - INTEGER(i4) :: k - REAL(dp) :: theta - - zroots_unity_dp(1) = 1.0_dp - theta = TWOPI_dp / n - k = 1 - do - if (k >= nn) exit - zroots_unity_dp(k + 1) = cmplx(cos(k * theta), sin(k * theta), kind = dpc) - zroots_unity_dp(k + 2 : min(2 * k, nn)) = zroots_unity_dp(k + 1) * & - zroots_unity_dp(2 : min(k, nn - k)) - k = 2 * k - end do - - END FUNCTION zroots_unity_dp - - - FUNCTION zroots_unity_sp(n, nn) - - INTEGER(i4), INTENT(IN) :: n, nn - COMPLEX(spc), DIMENSION(nn) :: zroots_unity_sp - INTEGER(i4) :: k - REAL(sp) :: theta - - zroots_unity_sp(1) = 1.0_sp - theta = TWOPI_sp / n - k = 1 - do - if (k >= nn) exit - zroots_unity_sp(k + 1) = cmplx(cos(k * theta), sin(k * theta), kind = spc) - zroots_unity_sp(k + 2 : min(2 * k, nn)) = zroots_unity_sp(k + 1) * & - zroots_unity_sp(2 : min(k, nn - k)) - k = 2 * k - end do - - END FUNCTION zroots_unity_sp - - ! ------------------------------------------------------------------ - END MODULE mo_corr diff --git a/src/lib/mo_kind.f90 b/src/lib/mo_kind.f90 old mode 100644 new mode 100755 index 992cd69c..f6e37651 --- a/src/lib/mo_kind.f90 +++ b/src/lib/mo_kind.f90 @@ -22,6 +22,7 @@ ! - documentation ! - removed tab characters ! Matthias Cuntz, May 2014 - iso_fortran_env and iso_c_binding +! S. Mueller, Dec 2019 - remove NR specific types ! License ! ------- @@ -48,34 +49,34 @@ MODULE mo_kind ! Does not work with compilers intel v11 and sun v12.2 ! use, intrinsic :: iso_fortran_env, only: & ! int8!, int16, int32, int64, real32, real64 - use, intrinsic :: iso_c_binding, only : & - c_short, c_int, c_long_long, c_float, c_double, c_float_complex, c_double_complex, c_bool + use, intrinsic :: iso_c_binding, only: & + c_short, c_int, c_long_long, c_float, c_double, c_float_complex, c_double_complex, c_bool IMPLICIT NONE !> 1 Byte Integer Kind - INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND(2) + INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND(2) ! INTEGER, PARAMETER :: i1 = int8 ! c_word does not exist; should be c_bool, probably; but see above !> 2 Byte Integer Kind ! INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND(4) ! INTEGER, PARAMETER :: i2 = int16 - INTEGER, PARAMETER :: i2 = c_short + INTEGER, PARAMETER :: i2 = c_short !> 4 Byte Integer Kind ! INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND(9) ! INTEGER, PARAMETER :: i4 = int32 - INTEGER, PARAMETER :: i4 = c_int + INTEGER, PARAMETER :: i4 = c_int !> 8 Byte Integer Kind ! INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND(18) ! INTEGER, PARAMETER :: i8 = int64 - INTEGER, PARAMETER :: i8 = c_long_long + INTEGER, PARAMETER :: i8 = c_long_long !> Single Precision Real Kind ! INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6,37) ! INTEGER, PARAMETER :: sp = real32 - INTEGER, PARAMETER :: sp = c_float + INTEGER, PARAMETER :: sp = c_float !> Double Precision Real Kind ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,307) ! INTEGER, PARAMETER :: dp = real64 - INTEGER, PARAMETER :: dp = c_double + INTEGER, PARAMETER :: dp = c_double !> Single Precision Complex Kind ! INTEGER, PARAMETER :: spc = KIND((1.0_sp,1.0_sp)) ! INTEGER, PARAMETER :: spc = sp @@ -89,19 +90,4 @@ MODULE mo_kind !INTEGER, PARAMETER :: lgt = c_bool ! Types have to be in a public section for doxygen - !> Single Precision Numerical Recipes types for sparse arrays - TYPE sprs2_sp - INTEGER(I4) :: n, len - REAL(SP), DIMENSION(:), POINTER :: val - INTEGER(I4), DIMENSION(:), POINTER :: irow - INTEGER(I4), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_sp - !> Double Precision Numerical Recipes types for sparse arrays - TYPE sprs2_dp - INTEGER(I4) :: n, len - REAL(DP), DIMENSION(:), POINTER :: val - INTEGER(I4), DIMENSION(:), POINTER :: irow - INTEGER(I4), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_dp - END MODULE mo_kind diff --git a/src/lib/mo_moment.f90 b/src/lib/mo_moment.f90 index 729dd590..50bdcac3 100644 --- a/src/lib/mo_moment.f90 +++ b/src/lib/mo_moment.f90 @@ -7,10 +7,6 @@ MODULE mo_moment ! i.e. they are normally divided by n and not (n-1) ! Literature - ! Moments (incl. mean, stddev, etc.) - Chapter 14 of - ! WH Press, SA Teukolsky, WT Vetterling, BP Flannery, - ! Numerical Recipes in Fortran 90 - The Art of Parallel Scientific Computing, 2nd Edition - ! Volume 2 of Fortran Numerical Recipes, Cambridge University Press, UK, 1996 ! Central moments and error variances ! LH Benedict & RD Gould, Towards better uncertainty estimates for turbulence statistics, ! Experiments in Fluids 22, 129-136, 1996 @@ -18,6 +14,7 @@ MODULE mo_moment ! Written Nov 2011, Matthias Cuntz ! Modified, MC, Dec 2011 - mod. correlation, covariance ! Modified by M. Schroen, Sep 2015, average/mean for single value + ! Modified by S. Mueller, Dec 2019, no citations for common sence ! License ! ------- @@ -100,11 +97,6 @@ MODULE mo_moment ! m = absdev(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE absdev @@ -153,11 +145,6 @@ MODULE mo_moment ! m = average(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE average @@ -316,11 +303,6 @@ MODULE mo_moment ! m = correlation(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.))) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 ! Modified, MC, Dec 2011 - covariance as <(x-)(y-)> instead of - @@ -372,11 +354,6 @@ MODULE mo_moment ! m = covariance(vec1, vec2, mask=((vec1 >= 0.) .and. (vec2 >= 0.))) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 ! Modified, MC, Dec 2011 - covariance as <(x-)(y-)> instead of - @@ -427,11 +404,6 @@ MODULE mo_moment ! m = kurtosis(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE kurtosis @@ -480,11 +452,6 @@ MODULE mo_moment ! m = mean(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE mean @@ -612,10 +579,10 @@ MODULE mo_moment ! moment ! PURPOSE - ! Calculates the first four sample moments of a vector, i.e. mean, variance, skewness, and kurtosis, + ! Calculates the first four sample/population moments of a vector, i.e. mean, variance, skewness, and kurtosis, ! as well as standard deviation and mean absolute devation. ! mean = sum(x)/n - ! variance = sum((x-mean(x))**2)/(n-1) + ! variance = sum((x-mean(x))**2)/(n-1) or sum((x-mean(x))**2)/n if sample is False ! skewness = sum(((x-mean(x))/stddev(x))**3)/n ! kurtosis = sum(((x-mean(x))/stddev(x))**4)/n - 3 ! stddev = sqrt(variance) @@ -627,7 +594,7 @@ MODULE mo_moment ! CALLING SEQUENCE ! call moment(dat, average=average, variance=variance, skewness=skewness, kurtosis=kurtosis, & - ! mean=mean, stddev=stddev, absdev=absdev, mask=mask) + ! mean=mean, stddev=stddev, absdev=absdev, mask=mask, sample=sample) ! INTENT(IN) ! real(sp/dp) :: dat(:) 1D-array with input numbers @@ -641,17 +608,19 @@ MODULE mo_moment ! INTENT(IN), OPTIONAL ! logical :: mask(:) 1D-array of logical values with size(dat). ! If present, only those locations in vec corresponding to the true values in mask are used. + ! logical :: sample logical value + ! If present and False, the population variance and std-dev will be calculated (divide by n). ! INTENT(INOUT), OPTIONAL ! None ! INTENT(OUT), OPTIONAL ! real(sp/dp) :: average average of input vector - ! real(sp/dp) :: variance sample variance of input vector + ! real(sp/dp) :: variance sample variance of input vector (either a sample or pupulation moment) ! real(sp/dp) :: skewness skewness of input vector ! real(sp/dp) :: kurtosis excess kurtosis of input vector ! real(sp/dp) :: mean same as average - ! real(sp/dp) :: stddev sqaure root of variance + ! real(sp/dp) :: stddev sqaure root of variance (either a sample or pupulation moment) ! real(sp/dp) :: absdev mean absolute deviations from average ! RESTRICTIONS @@ -662,13 +631,9 @@ MODULE mo_moment ! call moment(vec, mask=(vec >= 0.), mean=m, stddev=s) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 + ! Modified, S. Mueller, Dec 2019 added optional sample input-para to switch sample to population variance and std-dev. INTERFACE moment MODULE PROCEDURE moment_sp, moment_dp END INTERFACE moment @@ -715,11 +680,6 @@ MODULE mo_moment ! m = skewness(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE skewness @@ -770,11 +730,6 @@ MODULE mo_moment ! m = stddev(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE stddev @@ -825,11 +780,6 @@ MODULE mo_moment ! m = variance(vec, mask=(vec >= 0.)) ! -> see also example in test directory - ! LITERATURE - ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 - - ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes, - ! Cambridge University Press, UK, 1996 - ! HISTORY ! Written, Matthias Cuntz, Nov 2011 INTERFACE variance @@ -1586,7 +1536,7 @@ END FUNCTION mixed_central_moment_var_sp ! ------------------------------------------------------------------ - SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask) + SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample) IMPLICIT NONE @@ -1599,8 +1549,9 @@ SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, a REAL(dp), OPTIONAL, INTENT(OUT) :: stddev REAL(dp), OPTIONAL, INTENT(OUT) :: absdev LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + LOGICAL, OPTIONAL, INTENT(IN) :: sample - REAL(dp) :: n + REAL(dp) :: n, div_n ! divisor depending if sample or population moments REAL(dp) :: ep, ave, var REAL(dp), DIMENSION(size(dat)) :: p, s @@ -1616,6 +1567,12 @@ SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, a end if if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_dp: n must be at least 2' + ! set divisor for population (n) or sample (n-1) variance + div_n = n - 1.0_dp + if (present(sample)) then + if (.not. sample) div_n = n + end if + ! Any optional argument if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. & present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return @@ -1634,7 +1591,7 @@ SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, a ep = sum(s(:), mask = maske) p(:) = s(:) * s(:) var = sum(p(:), mask = maske) - var = (var - ep * ep / n) / (n - 1.0_dp) + var = (var - ep * ep / n) / div_n if (present(variance)) variance = var ! Standard deviation if (present(stddev)) stddev = sqrt(var) @@ -1656,7 +1613,7 @@ SUBROUTINE moment_dp(dat, average, variance, skewness, kurtosis, mean, stddev, a END SUBROUTINE moment_dp - SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask) + SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, absdev, mask, sample) IMPLICIT NONE @@ -1669,8 +1626,9 @@ SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, a REAL(sp), OPTIONAL, INTENT(OUT) :: stddev REAL(sp), OPTIONAL, INTENT(OUT) :: absdev LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask + LOGICAL, OPTIONAL, INTENT(IN) :: sample - REAL(sp) :: n + REAL(sp) :: n, div_n ! divisor depending if sample or population moments REAL(sp) :: ep, ave, var REAL(sp), DIMENSION(size(dat)) :: p, s @@ -1686,6 +1644,12 @@ SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, a end if if (n .le. (1.0_sp + tiny(1.0_sp))) stop 'moment_sp: n must be at least 2' + ! set divisor for population (n) or sample (n-1) variance + div_n = n - 1.0_sp + if (present(sample)) then + if (.not. sample) div_n = n + end if + ! Any optional argument if (.not. (present(average) .or. present(variance) .or. present(skewness) .or. & present(kurtosis) .or. present(mean) .or. present(stddev) .or. present(absdev))) return @@ -1704,7 +1668,7 @@ SUBROUTINE moment_sp(dat, average, variance, skewness, kurtosis, mean, stddev, a ep = sum(s(:), mask = maske) p(:) = s(:) * s(:) var = sum(p(:), mask = maske) - var = (var - ep * ep / n) / (n - 1.0_sp) + var = (var - ep * ep / n) / div_n if (present(variance)) variance = var ! Standard deviation if (present(stddev)) stddev = sqrt(var)