! time.f90
! module to hold standard routines for manipulating temporal things
! contains: GaussSmooth,GetMonthLengths,CommonVecPer

module Time

implicit none

contains

!*******************************************************************************
! compare year vectors to find common period

subroutine CommonVecPer (AYears,BYears,AStart,AEnd,BStart,BEnd)

integer, pointer, dimension (:) :: AYears, BYears

integer, intent (out) :: AStart,BStart,AEnd,BEnd

integer :: AYearN,BYearN,ARest,BRest
integer :: Offset,CommonLen

real, parameter :: MissVal = -999.0

AYearN = size (AYears)
BYearN = size (BYears)

Offset = BYears(1) - AYears(1) + 1
if (Offset.GE.1) then
  AStart = Offset
  BStart = 1
else
  AStart = 1
  BStart = 2 - Offset
end if

ARest = AYearN - AStart
BRest = BYearN - BStart

if (ARest.GE.0.AND.BRest.GE.0) then
  if (ARest.LT.BRest) then
    CommonLen = ARest + 1
  else
    CommonLen = BRest + 1
  end if
  
  AEnd = AStart + CommonLen - 1
  BEnd = BStart + CommonLen - 1  
else
  AStart=MissVal
  BStart=MissVal
  AEnd=MissVal
  BEnd=MissVal
end if

end subroutine CommonVecPer

!*******************************************************************************
! smooths using Gaussian filtering

subroutine GaussSmooth (TimeN,THalf,QRetainMiss,TSInVec,TSLowVec,TSHighVec,MissThresh)

real, pointer, dimension (:) 	:: TSInVec,TSLowVec,TSHighVec	! def,alloc,init in call
integer, intent (in) 		:: TimeN, THalf, QRetainMiss	! QRetainMiss:1=no,2=yes
real, intent (in), optional	:: MissThresh			! %miss=OK ; if>,output=MissVal

real, allocatable, dimension (:) 	:: Weight, TSPad, TSPadExt

real, parameter :: MissVal = -999.0

real	:: WeightFactor, WeightRoot, WeightSum
real    :: LocalTot, LocalN
real    :: MeanStart, MeanEnd
real	:: Real0, Real1
real    :: MissPercent

integer :: AllocStat
integer :: Integer0, Integer1
integer :: WeightN
integer :: XWeight, XTime, XValue
integer :: EndN
integer :: MissTot, QSetMiss

QSetMiss = 0

if (present(MissThresh)) then
 MissTot = 0
 do XTime = 1, TimeN
  if (TSInVec(XTime).EQ.MissVal) MissTot = MissTot + 1
 end do
 MissPercent = (real(MissTot) / real(TimeN)) * 100.0

 if (MissPercent.GE.MissThresh) then
  TSInVec = 999 ; QSetMiss = 1
 end if
end if

!write (99,*), "##### computing no. weights"			! ###################

WeightN = nint ((real(THalf)/2.5)+0.5)		! compute no. weights required
if (modulo(WeightN,2).EQ.0) then
	WeightN = WeightN + 2
else
	WeightN = WeightN + 1
end if
if (WeightN.LE.7) WeightN = 7

allocate (Weight(WeightN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
Weight = 0.0
						! compute weights
!write (99,*), "##### computing weights"				! ###################

WeightFactor = -18.0 / (real(THalf*THalf))
WeightRoot   = 1.0 / sqrt (2.0 * 3.1415927)
Weight (1)   = WeightRoot
WeightSum    = Weight (1)

do XWeight = 2, WeightN
  Weight (XWeight) = WeightRoot * exp (WeightFactor * real(XWeight) * real(XWeight))
  WeightSum = WeightSum + 2.0 * Weight (XWeight)
end do

do XWeight = 1, WeightN
  Weight (XWeight) = Weight (XWeight) / WeightSum
end do
						! pad the ts with local mean where missing 
!write (99,*), "##### padding where missing"			! ###################

allocate (TSPad(TimeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
TSPad = 0.0

do XTime = 1, TimeN
!  write (99,"(i4,f8.2)"), XTime,TSInVec(XTime)			! ###################
  if (TSInVec(XTime).EQ.MissVal) then
    Real0  = max ((XTime-WeightN+1),1)
    Real1  = min ((XTime+WeightN-1),TimeN)
    
    Integer0 = nint (Real0)
    Integer1 = nint (Real1)
    
    LocalTot = 0.0
    LocalN   = 0.0

    do XValue = Integer0, Integer1
      if (TSInVec(XValue).NE.MissVal) then
        LocalTot = LocalTot + TSInVec(XValue)
        LocalN   = LocalN   + 1.0
      end if
    end do
    
    if (LocalN.GT.0.AND.LocalTot.NE.0) TSPad (XTime) = LocalTot / LocalN
  else
    TSPad (XTime) = TSInVec (XTime) 
  end if
end do
						! extend ends of ts by mean from each end
!write (99,*), "##### extending ends"			! ###################

EndN = WeightN - 1

LocalTot = 0.0
LocalN   = 0.0
do XTime = 1, EndN
        LocalTot = LocalTot + TSPad(XTime)
        LocalN   = LocalN   + 1.0
end do
MeanStart = LocalTot / LocalN

LocalTot = 0.0
LocalN   = 0.0
do XTime = (TimeN-EndN+1), TimeN
        LocalTot = LocalTot + TSPad(XTime)
        LocalN   = LocalN   + 1.0
end do
MeanEnd = LocalTot / LocalN

allocate (TSPadExt(TimeN+2*EndN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
TSPadExt = 0.0

do XValue = 1, EndN
  TSPadExt (XValue) = MeanStart
end do
do XValue = (EndN+1), (EndN+TimeN)
  TSPadExt (XValue) = TSPad (XValue-EndN)
end do
do XValue = (EndN+TimeN+1), (TimeN+2*EndN)
  TSPadExt (XValue) = MeanEnd
end do
						! apply the filter
!write (99,*), "##### applying filter"			! ###################
! if it fails below, it is probably because you haven't alloc all the arrays in call

do XTime = 1, TimeN
  WeightSum = Weight (1) * TSPadExt (EndN+XTime)
  
  do XWeight = 2, WeightN
    WeightSum = WeightSum + Weight(XWeight) * &
    			(TSPadExt(XTime+EndN-XWeight+1)+TSPadExt(XTime+EndN+XWeight-1))
  end do
  
  TSLowVec (XTime) = WeightSum
end do
						! compute the high-residual
!write (99,*), "##### computing high residual"		! ###################

do XTime = 1, TimeN
  TSHighVec (XTime) = TSInVec (XTime) - TSLowVec (XTime)
end do
						! reinsert MissVal unless called otherwise
if (QRetainMiss.EQ.2) then
 do XTime = 1, TimeN
  if (TSInVec(XTime).EQ.MissVal) then
   TSLowVec (XTime) = MissVal
   TSHighVec(XTime) = MissVal
  end if
 end do
end if

deallocate (Weight, TSPad, TSPadExt)

if (QSetMiss.EQ.1) then
  TSLowVec = MissVal ; TSHighVec = MissVal
end if

end subroutine GaussSmooth

!*******************************************************************************
! get month lengths appropriate to years AD

subroutine GetMonthLengths (YearAD,MonthLengths)

integer, pointer, dimension (:,:) 	:: MonthLengths
integer, pointer, dimension (:) 	:: YearAD

real    :: Rem004,Rem100,Rem400
integer :: XYear, XMonth, YearN

YearN = size (YearAD)

if (.NOT.associated(MonthLengths)) print*, &
	"  > ##### GetMonthLengths: must be alloc in call #####"

MonthLengths = 31

do XYear = 1, YearN
  MonthLengths(XYear,4 ) = 30
  MonthLengths(XYear,6 ) = 30
  MonthLengths(XYear,9 ) = 30
  MonthLengths(XYear,11) = 30
  
  Rem004 = mod (YearAD(XYear),  4)
  Rem100 = mod (YearAD(XYear),100)
  Rem400 = mod (YearAD(XYear),400)
  
  if (Rem004.EQ.0) then
    if (Rem100.EQ.0) then
      if (Rem400.EQ.0) then
        MonthLengths(XYear,2) = 29
      else
        MonthLengths(XYear,2) = 28
      end if
    else
        MonthLengths(XYear,2) = 29
    end if
  else
        MonthLengths(XYear,2) = 28
  end if
end do

end subroutine GetMonthLengths

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

end module Time