! slicestats.f90
! module procedure written by Tim Mitchell on 14.12.99
! requires SortMod.f90

module SliceStats

implicit none

contains

!*******************************************************************************
! converts single time series of data into slice Statistic specified
! median calc after sorting data; skew also requires nedian
! other stats all require summation of powers of the data instead
! median, mean, variance are calc as sample stats in usual manner
! skewness = (3*(mean--median))/stdev  (Francis p88): <0=left skew=neg skew >0=opposite 0=normal
! kurtosis = (fourth moment/(second moment**2))--3 (Francis p88) -3<sk<0=platykurtic sk>0=leptokurtic 0=normal

subroutine SliceStatistics (JobYearN,SliceLen,Statistic,MissAccept,DataVec)

use SortMod

real, pointer, dimension (:) 	:: DataVec, SortVec
integer, pointer, dimension (:) :: Ranks

real, intent (in) 	:: MissAccept				
integer, intent (in) 	:: JobYearN, SliceLen, Statistic	! Stat:0=medi,1=mean,2=var,3=skew,4=kurt

real, allocatable, dimension (:) :: Sliced

real, parameter :: MissVal = -999.0

real :: Median
real :: Mom1, Mom2, Mom3, Mom4, MomZero1, MomZero2, MomZero3, MomZero4
real :: DataEn, DataTot, DataSqTot, DataCubTot, DataQuadTot
real :: Threshold
real :: MidReal

integer :: AllocStat
integer :: XMem, XYear, XValue
integer :: SliceStartYear, SliceMidYear, MidValue
integer :: SliceMissN

!**************************

Median=MissVal
MomZero1=MissVal ; MomZero2=MissVal ; MomZero3=MissVal ; MomZero4=MissVal
Mom1=MissVal ; Mom2=MissVal ; Mom3=MissVal ; Mom4=MissVal

allocate (Sliced(JobYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatistics: Allocation failure #####"
Sliced = MissVal

if (Statistic.EQ.3.OR.Statistic.EQ.0) then
      allocate (SortVec(SliceLen),     &
      		Ranks(SliceLen), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatistics: Allocation failure #####"
      Ranks = MissVal
end if
      
Threshold = SliceLen*(100-MissAccept)/100

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  DataEn       = 0.0

  if (Statistic.EQ.3.OR.Statistic.EQ.0) then				! for median or skewness calc
      do XValue = 1, SliceLen						! set up sorting array
        XYear = SliceStartYear + XValue - 1
        SortVec (XValue) = DataVec (XYear)
      end do
      
      Median = MedianValue (SliceLen,SortVec)
  end if
  
  if (Statistic.GE.1) then  						! for moments calc
    DataTot      = 0.0
    DataSqTot    = 0.0
    DataCubTot   = 0.0
    DataQuadTot  = 0.0
    
    do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (DataVec(XYear).NE.MissVal) then
        if (Statistic.GE.4) DataQuadTot = DataQuadTot + (DataVec(XYear) ** 4)
        if (Statistic.GE.3) DataCubTot  = DataCubTot  + (DataVec(XYear) ** 3)
        if (Statistic.GE.2) DataSqTot   = DataSqTot   + (DataVec(XYear) ** 2)
        if (Statistic.GE.1) DataTot     = DataTot     +  DataVec(XYear)

        DataEn      = DataEn      +  1.0
      end if
    end do
  end if
  
  if (DataEn.GE.Threshold.AND.DataEn.GE.1) then
    if (Statistic.GE.1) MomZero1 = DataTot     / DataEn			! moments about zero
    if (Statistic.GE.2) MomZero2 = DataSqTot   / DataEn
    if (Statistic.GE.3) MomZero3 = DataCubTot  / DataEn
    if (Statistic.GE.4) MomZero4 = DataQuadTot / DataEn
  
    if (Statistic.GE.1) Mom1 = 0					! moments about the mean
    if (Statistic.GE.2) Mom2 = MomZero2 - (MomZero1 ** 2)
    if (Statistic.GE.3) Mom3 = MomZero3 - (3 * MomZero1 * MomZero2) + (2 * (MomZero1 ** 3))
    if (Statistic.GE.4) Mom4 = MomZero4 - (4 * MomZero1 * MomZero3) + (6 * (MomZero1 ** 2) * MomZero2) &
  				- (3 * (MomZero1 ** 4))

    if (Statistic.EQ.1) Sliced (SliceMidYear)   = MomZero1					! mean
    if (Statistic.EQ.2) Sliced (SliceMidYear)   = Mom2						! variance
    
    if (Mom2.GT.0) then
      if (Statistic.EQ.3) Sliced (SliceMidYear) = (3 * (MomZero1 - Median)) / sqrt (Mom2)	! skewness
      if (Statistic.EQ.4) Sliced (SliceMidYear) = (Mom4 / (Mom2 ** 2)) - 3			! kurtosis
    end if
  end if

  if (Statistic.EQ.0.AND.Median.NE.MissVal) &			! median
  		     Sliced (SliceMidYear) = Median
end do

do XYear = 1, JobYearN
  DataVec(XYear) = Sliced(XYear)
end do

if (Statistic.EQ.3.OR.Statistic.EQ.0) then
      deallocate (SortVec,Ranks, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatistics: Deallocation failure #####"
end if
      
deallocate (Sliced, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatistics: Deallocation failure #####"

end subroutine SliceStatistics

!*******************************************************************************
! converts ensemble time series of data into slice statistic specified
! median calc after sorting data; skew also requires median
! other stats all require summation of powers of the data instead
! median, mean, variance are calc as sample stats in usual manner
! skewness = (3*(mean--median))/stdev  (Francis p88): <0=left skew=neg skew >0=opposite 0=normal
! kurtosis = (fourth moment/(second moment**2))--3 (Francis p88) -3<sk<0=platykurtic sk>0=leptokurtic 0=normal

subroutine SliceStatisticsEns (JobMemN,JobYearN,SliceLen,Statistic,MissAccept,InArray,OutVec)

use SortMod

real, pointer, dimension (:,:) 	:: InArray				! JobMemN,JobYearN
real, pointer, dimension (:) 	:: OutVec, SortVec
integer, pointer, dimension (:) :: Ranks

real, intent (in) 	:: MissAccept					
integer, intent (in) 	:: JobMemN, JobYearN, SliceLen
integer, intent (in) 	:: Statistic				! 0=medi,1=mean,2=var,3=skew,4=kurt

real, parameter :: MissVal = -999.0

real :: Median
real :: MidReal
real :: Mom1, Mom2, Mom3, Mom4, MomZero1, MomZero2, MomZero3, MomZero4
real :: DataEn, DataTot, DataSqTot, DataCubTot, DataQuadTot
real :: Threshold

integer :: AllocStat
integer :: XMem, XYear, XValue
integer :: SliceStartYear, SliceMidYear, MidValue
integer :: SliceMissN, EnsDataN, EnsValue

!**************************

Median=MissVal
MomZero1=MissVal ; MomZero2=MissVal ; MomZero3=MissVal ; MomZero4=MissVal
Mom1=MissVal ; Mom2=MissVal ; Mom3=MissVal ; Mom4=MissVal

OutVec = MissVal

EnsDataN = SliceLen*JobMemN
      
if (Statistic.EQ.3.OR.Statistic.EQ.0) then
      allocate (SortVec(EnsDataN),     &
      		Ranks(EnsDataN), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatisticsEns: Allocation failure #####"
      Ranks = MissVal
end if
      
Threshold = JobMemN*SliceLen*(100-MissAccept)/100

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  DataEn       = 0.0

  if (Statistic.EQ.3.OR.Statistic.EQ.0) then				! for median or skewness calc
      do XValue = 1, SliceLen						! set up sorting array
        XYear = SliceStartYear + XValue - 1
        
        do XMem = 1, JobMemN
          EnsValue  = ((XValue-1)*JobMemN) + XMem
          SortVec (EnsValue) = InArray (XMem,XYear)
        end do
      end do
      
      Median = MedianValue (EnsDataN,SortVec)
  end if
  
  if (Statistic.GE.1) then  						! for moment calc
    DataTot      = 0.0
    DataSqTot    = 0.0
    DataCubTot   = 0.0
    DataQuadTot  = 0.0
    
    do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
     do XMem = 1, JobMemN
      if (InArray(XMem,XYear).NE.MissVal) then
        if (Statistic.GE.4) DataQuadTot = DataQuadTot + (InArray(XMem,XYear) ** 4)
        if (Statistic.GE.3) DataCubTot  = DataCubTot  + (InArray(XMem,XYear) ** 3)
        if (Statistic.GE.2) DataSqTot   = DataSqTot   + (InArray(XMem,XYear) ** 2)
        if (Statistic.GE.1) DataTot     = DataTot     +  InArray(XMem,XYear)

        DataEn      = DataEn      +  1.0
      end if
     end do
    end do
  end if
  
  if (DataEn.GE.Threshold.AND.DataEn.GE.1) then
    if (Statistic.GE.1) MomZero1 = DataTot     / DataEn		! moments about zero
    if (Statistic.GE.2) MomZero2 = DataSqTot   / DataEn
    if (Statistic.GE.3) MomZero3 = DataCubTot  / DataEn
    if (Statistic.GE.4) MomZero4 = DataQuadTot / DataEn
  
    if (Statistic.GE.1) Mom1 = 0					! moments about the mean
    if (Statistic.GE.2) Mom2 = MomZero2 - (MomZero1 ** 2)
    if (Statistic.GE.3) Mom3 = MomZero3 - (3 * MomZero1 * MomZero2) + (2 * (MomZero1 ** 3))
    if (Statistic.GE.4) Mom4 = MomZero4 - (4 * MomZero1 * MomZero3) + (6 * (MomZero1 ** 2) * MomZero2) &
  				- (3 * (MomZero1 ** 4))

    if (Statistic.EQ.1) OutVec (SliceMidYear)   = MomZero1					! mean
    if (Statistic.EQ.2) OutVec (SliceMidYear)   = Mom2						! variance
    
    if (Mom2.GT.0) then
      if (Statistic.EQ.3) OutVec (SliceMidYear) = (3 * (MomZero1 - Median)) / sqrt (Mom2)	! skewness
      if (Statistic.EQ.4) OutVec (SliceMidYear) = (Mom4 / (Mom2 ** 2)) - 3			! kurtosis
    end if
  end if

  if (Statistic.EQ.0.AND.Median.NE.MissVal) &							! median
  		     OutVec (SliceMidYear) = Median
    
end do

if (Statistic.EQ.3.OR.Statistic.EQ.0) then
      deallocate (SortVec,Ranks, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStatisticsEns: Deallocation failure #####"
end if
      
end subroutine SliceStatisticsEns

!*******************************************************************************

end module SliceStats