!----------------------------------------------------------------------------- ! ! std_dev -- Compute mean, variance, and standard deviation for a time series. ! ! This version also includes data masking, to facilitate things like ! seasonal subsetting. ! ! 1.00 2008-dec-30 Original version. By Dave Allured. ! Adapted from corr_regress v1.00, ! also from standardize_seasonal v2.02. ! Compute mean and standard deviation for 1-D time series. ! 1.01 2009-jan-06 Improvements to match std_dev_grid v1.01. ! Change order of calling arguments. ! Change required minimum count to fractional threshold. ! Add variance output. ! Fix incorrect output for degenerate time series. ! ! ---- 2014-mar-02 Double precision version. ! Remove threshold checks, require exact missing values. ! ! Input: 1-D time series, 1-D data mask, and control parameters. See below. ! ! Output: Mean, variance, and standard deviation, scalars. ! In case of insufficient data, outputs are set to vmiss. ! ! Notes: ! ! Dimension T is simply the data indexing dimension. It is called ! the time dimension herein, but it may have any definition suitable ! to the calling program. ! ! Time subsetting for random time steps is supported, by using the ! tags(T) data mask array to select time steps to be included. ! This may be used, for example, to get seasonal statistics from ! year-round (continuous) time series. ! ! This version uses the two pass algorithm to minimize roundoff error. ! !----------------------------------------------------------------------------- module std__dev use config, only : dp implicit none ! Special public limit parameters for debugging. ! No check is done unless at least one is set positive. integer :: std_dev_min_count = -1 integer :: std_dev_max_count = -1 contains subroutine std_dev (data, tags, vmiss, thresh, mean, variance, stdev) implicit none real(dp), intent (in ) :: data(:) ! 1-D input time series (T) logical, intent (in ) :: tags(:) ! include only TRUE data pairs (T) real(dp), intent (in ) :: vmiss ! missing value in input data real(dp), intent (in ) :: thresh ! 0-1: required fraction of data present ! applies only to tag selected data real(dp), intent (out) :: mean ! computed mean real(dp), intent (out) :: variance ! computed variance real(dp), intent (out) :: stdev ! computed standard deviation ! Local variables. integer ntimes, ndata, required_count, j real(dp) val real(dp) sumx, sdd, dif, dmean, dvariance ! Derived parameters. ntimes = size (data) required_count = ceiling (thresh * count (tags)) ! # of time steps required ! for each grid point for valid output required_count = max (1, required_count) ! require at least one datum ! at each grid point ! Consistency check. if (ntimes /= size (tags)) then print * write (*, '(2(a,i0))') ' *** std_dev: Time dimensions for data' & // ' and tag arrays = ', ntimes, ', ', size (tags) print *, '*** std_dev: Abort: Data and tag time dimensions differ.' call exit (1) end if ! Compute the mean for selected time steps. sumx = 0 ndata = 0 do j = 1, ntimes if (.not. tags(j)) cycle ! include only selected times val = data(j) if (val /= vmiss) then sumx = sumx + val ndata = ndata + 1 end if end do if (ndata > 0) then ! suppress divide by zero dmean = sumx / ndata ! double precision for calcs mean = dmean ! single precision for output end if ! Compute variance and standard deviation. sdd = 0 if (ndata > 1) then ! normal points, 2 or more data: do j = 1, ntimes ! accumulate deviations squared if (.not. tags(j)) cycle val = data(j) if (val /= vmiss) then dif = val - dmean sdd = sdd + (dif * dif) end if end do dvariance = sdd / (ndata - 1) ! compute variance variance = dvariance stdev = sqrt (dvariance) ! compute standard deviation else ! for degenerate points variance = 0 ! with only one datum: stdev = 0 ! variance is zero end if ! DEBUG: Check data count. This section is completely inactive ! if debug parameters are left at their default values. if (std_dev_min_count >= 0 .or. std_dev_max_count >= 0) then if (ndata < std_dev_min_count .or. ndata > std_dev_max_count) then print * print *, '*** std_dev: tags array:' do j = 1, ntimes, 20 write (*, "(i6, ' :', 20l2)") j, tags(j:min (j+19, ntimes)) end do print * write (*, '(a, i0, a)') ' *** std_dev: Abort: Data count = ', ndata, & ', outside of debug limits.' call exit (1) end if end if ! Output missing values if data count is below threshold. ! Suppress output if data count is below threshold. if (ndata < required_count) then ! at least one datum always required mean = vmiss variance = vmiss stdev = vmiss end if end subroutine std_dev end module std__dev