! transformmod.f90
! module procedure written by Tim Mitchell on 14.12.99
! last modification on 29.02.00
!   transforms temp deg C to temk deg K (only to be rarely used; see ExtractMod for vice versa)
!   transforms data into anom abs
!   transforms data into anom percent
!   detrends data
!   ### BubbleSort has moved to basicfun.f90 ###

module TransformMod

implicit none

contains

!*******************************************************************************
! detrend a region-set of time series

subroutine DeTrendTim (RegN,YearN,MissAccept,QPrint,ADYear,Data,RegAye,RegBee)

real, pointer, dimension (:,:)	:: Data
real, pointer, dimension (:)	:: RegAye, RegBee

integer, pointer, dimension (:)	:: ADYear

real, intent (in)	:: MissAccept
real, parameter		:: MissVal = -999.0
real			:: RegThresh, YearThresh
real			:: Num, Den
real			:: SumX, SumY, SumXX, SumYY, SumXY, En
real			:: ValidReg, ValidAye, ValidYear

integer, intent (in)	:: RegN, YearN, QPrint
integer			:: XReg, XYear
integer			:: AllocStat, ReadStatus

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

if (QPrint.EQ.1) print*, "  > Detrending a region-set of time series..."

allocate ( RegAye (RegN), &
	   RegBee (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DeTrendTim: Allocation failure #####"
RegAye = MissVal
RegBee = MissVal

YearThresh = YearN * (100.0-MissAccept)/100.0
RegThresh  = RegN  * (100.0-MissAccept)/100.0
ValidReg   = 0.0
ValidAye   = 0.0
ValidYear  = 0.0

do XReg = 1, RegN
  SumX  = 0.0						! initialise sums
  SumY  = 0.0
  SumXX = 0.0
  SumXY = 0.0
  SumXY = 0.0
  En    = 0.0
  
  do XYear = 1, YearN					! add together into sums
    if (Data(XReg,XYear).NE.MissVal.AND.ADYear(XYear).NE.MissVal) then
      SumX  = SumX  +  ADYear (XYear)
      SumY  = SumY  +  Data (XReg,XYear)
      SumXX = SumXX + (ADYear (XYear) ** 2)
      SumYY = SumYY + (Data (XReg,XYear) ** 2)
      SumXY = SumXY +  ADYear (XYear) * Data (XReg,XYear)
      En    = En    +  1.0
    end if
  end do
  
  if (En.GE.YearThresh) then				! calc a and b in y=ax+b
    Num = (En*SumXY)-(SumX*SumY)
    Den = (En*SumXX)-(SumX*SumX)
    if (Den.NE.0) RegAye(XReg) = Num / Den

    if (RegAye(XReg).NE.MissVal) RegBee(XReg) = (SumY-(RegAye(XReg)*SumX))/En
    
    ValidReg = ValidReg + 1.0
  end if
end do

do XReg = 1, RegN					! detrend data series
  if (RegAye(XReg).NE.MissVal) then
    ValidAye = ValidAye + 1.0
    
    do XYear = 1, YearN
      if (Data(XReg,XYear).NE.MissVal) then
      	Data(XReg,XYear) = Data(XReg,XYear) - (RegAye(XReg) * XYear)
      	ValidYear = ValidYear + 1.0
      end if
    end do
  else
    do XYear = 1, YearN
      Data(XReg,XYear) = MissVal
    end do
  end if
end do  

ValidReg  = 100.0 * (ValidReg  /  RegN)
ValidAye  = 100.0 * (ValidAye  /  RegN)
ValidYear = 100.0 * (ValidYear / (RegN*YearN))
  
write (99,"(3f12.4)"), ValidReg, ValidAye, ValidYear

if (QPrint.EQ.1) then
  print*, "  > Percent of valid regions for LSR : ", ValidReg
  print*, "  > Percent of regions with valid 'a': ", ValidAye
  print*, "  > Percent of valid data detrended  : ", ValidYear 
end if

end subroutine DeTrendTim

!*******************************************************************************
! converts single time series of data into slice means

subroutine SimpleSlice (JobYearN,SliceLen,MissAccept,DataVec)

real, pointer, dimension (:) 	:: DataVec

real, intent (in) 	:: MissAccept
integer, intent (in) 	:: JobYearN, SliceLen

real, allocatable, dimension (:) :: Sliced

real, parameter :: MissVal = -999.0

real :: DataEn, DataTot

integer :: AllocStat
integer :: XMem, XYear, SliceStartYear, SliceMidYear
integer :: TotalMiss, Threshold

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

allocate (Sliced(JobYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SimpleSlice: Allocation failure #####"
Sliced = MissVal
Threshold = nint (MissAccept*real(SliceLen)/100.0)

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  TotalMiss    = 0
  DataEn       = 0.0
  DataTot      = 0.0
    
  do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (DataVec(XYear).EQ.MissVal) then
      	TotalMiss = TotalMiss  + 1
      else
        DataTot   = DataTot + DataVec(XYear)
        DataEn    = DataEn  + 1.0
      end if
  end do
    
  if (TotalMiss.LE.Threshold.AND.DataEn.GE.1) then	! checks that no. of valid values in slice is good
    	Sliced (SliceMidYear) = DataTot / DataEn
  end if
end do

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

deallocate (Sliced, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SimpleSlice: Deallocation failure #####"

end subroutine SimpleSlice

!*******************************************************************************
! converts parallel time series of data into a time series of slice means

subroutine Sliceify (JobMemN,JobYearN,SliceLen,MissAccept,DataArray,DataVec)

real, pointer, dimension (:) 	:: DataVec
real, pointer, dimension (:,:) 	:: DataArray

real, intent (in) :: MissAccept

integer, intent (in) :: JobMemN, JobYearN, SliceLen

real, parameter :: MissVal = -999.0

real :: DataEn, DataTot

integer :: XMem, XYear, SliceStartYear, SliceMidYear
integer :: TotalMiss, Threshold

DataVec = MissVal

Threshold = nint (MissAccept*real(SliceLen*JobMemN)/100.0)

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  TotalMiss    = 0
  DataEn       = 0.0
  DataTot      = 0.0
    
  do XMem = 1, JobMemN
    do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (DataArray(XMem,XYear).EQ.MissVal) then
      	TotalMiss = TotalMiss  + 1
      else
        DataTot = DataTot + DataArray(XMem,XYear)
        DataEn  = DataEn  + 1.0
      end if
    end do
  end do
    
  if (TotalMiss.LE.Threshold.AND.DataEn.GE.1) then	! checks that no. of valid values in slice is ok
    	DataVec (SliceMidYear) = DataTot / DataEn
  end if
end do

end subroutine Sliceify

!*******************************************************************************
! converts single time series of data into slice stdev

subroutine SliceStDevCol (JobYearN,SliceLen,QSampPop,MissAccept,DataVec)

real, pointer, dimension (:) 	:: DataVec

real, intent (in) 	:: MissAccept
integer, intent (in) 	:: JobYearN, SliceLen, QSampPop 	! sample=1, pop=2

real, allocatable, dimension (:) :: Sliced

real, parameter :: MissVal = -999.0

real :: DataEn, DataTot, DataSqTot

integer :: AllocStat
integer :: XMem, XYear, SliceStartYear, SliceMidYear
integer :: TotalMiss, Threshold

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

allocate (Sliced(JobYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStDevCol: Allocation failure #####"
Sliced = MissVal
Threshold = nint (MissAccept*real(SliceLen)/100.0)

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  TotalMiss    = 0
  DataEn       = 0.0
  DataTot      = 0.0
  DataSqTot    = 0.0
    
  do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (DataVec(XYear).EQ.MissVal) then
      	TotalMiss = TotalMiss  + 1
      else
        DataSqTot = DataSqTot + (DataVec(XYear) ** 2)
        DataTot   = DataTot   +  DataVec(XYear)
        DataEn    = DataEn    +  1.0
      end if
  end do
    
  if (TotalMiss.LE.Threshold.AND.DataEn.GE.2) then	! checks that no. of valid values in slice is good    	
    	Sliced (SliceMidYear) = (DataSqTot / DataEn) - ( (DataTot / DataEn) ** 2 )
    	if (Sliced(SliceMidYear).GT.0) then
    	  if (QSampPop.EQ.2) Sliced (SliceMidYear) = Sliced (SliceMidYear) * DataEn / (DataEn - 1)
    	  Sliced (SliceMidYear) = sqrt ( Sliced (SliceMidYear) )
    	else
    	  Sliced (SliceMidYear) = MissVal
    	end if
  end if
end do

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

deallocate (Sliced, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceStDevCol: Deallocation failure #####"

end subroutine SliceStDevCol

!*******************************************************************************
! converts ensemble of data into a time series of slice stdev

subroutine SliceStDevEns (JobMemN,JobYearN,SliceLen,QSampPop,MissAccept,DataArray,DataVec)

real, pointer, dimension (:) 	:: DataVec
real, pointer, dimension (:,:) 	:: DataArray

real, intent (in) :: MissAccept

integer, intent (in) :: JobMemN, JobYearN, SliceLen, QSampPop	! 1=sample, 2=pop

real, parameter :: MissVal = -999.0

real :: DataEn, DataTot, DataSqTot

integer :: XMem, XYear, SliceStartYear, SliceMidYear
integer :: TotalMiss, Threshold

DataVec = MissVal

Threshold = nint (MissAccept*real(SliceLen*JobMemN)/100.0)

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  TotalMiss    = 0
  DataEn       = 0.0
  DataTot      = 0.0
  DataSqTot    = 0.0
    
  do XMem = 1, JobMemN
    do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (DataArray(XMem,XYear).EQ.MissVal) then
      	TotalMiss = TotalMiss  + 1
      else
        DataSqTot = DataSqTot + (DataArray(XMem,XYear) ** 2)
        DataTot   = DataTot   +  DataArray(XMem,XYear)
        DataEn    = DataEn    +  1.0
      end if
    end do
  end do
    
  if (TotalMiss.LE.Threshold.AND.DataEn.GE.2) then
    	DataVec (SliceMidYear) = (DataSqTot / DataEn) - ( (DataTot / DataEn) ** 2)
    	if (DataVec (SliceMidYear).GT.0) then
    	  if (QSampPop.EQ.2) DataVec (SliceMidYear) = DataVec (SliceMidYear) * DataEn / (DataEn - 1)
    	  DataVec (SliceMidYear) = sqrt ( DataVec (SliceMidYear) )
    	else
    	  DataVec (SliceMidYear) = MissVal
    	end if
  end if
end do

end subroutine SliceStDevEns

!*******************************************************************************
! (only to be rarely used; see ExtractMod for vice versa)
! all work to be done in temp (degC); temk only for plotting purposes

subroutine TemptoTemk (RegN,DataVec,OutVec)

real, pointer, dimension (:) :: DataVec, OutVec		! def,alloc,init in call
integer, intent (in) :: RegN

real, parameter :: MissVal = -999.0
integer :: XReg

do XReg = 1, RegN
  if (DataVec(XReg).NE.MissVal) then
    OutVec (XReg) = DataVec (XReg) + 273.15
  else
    OutVec (XReg) = MissVal
  end if
end do

end subroutine TemptoTemk

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

subroutine AnomaliseAbs (RegN,BaseVec,DataVec,OutVec)

real, pointer, dimension (:) :: BaseVec,DataVec,OutVec	! def,alloc,init in call
integer, intent (in) :: RegN

real, parameter :: MissVal = -999.0
integer :: XReg

do XReg = 1, RegN
  if (DataVec(XReg).NE.MissVal.AND.BaseVec(XReg).NE.MissVal) then
    OutVec (XReg) = DataVec (XReg) - BaseVec (XReg)
  else
    OutVec (XReg) = MissVal
  end if
end do

end subroutine AnomaliseAbs

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

subroutine AnomalisePerCent (RegN,BaseVec,DataVec,OutVec)

real, pointer, dimension (:) :: BaseVec,DataVec,OutVec	! def,alloc,init in call
integer, intent (in) :: RegN

real, parameter :: MissVal = -999.0
integer :: XReg

do XReg = 1, RegN
  if (DataVec(XReg).NE.MissVal.AND.BaseVec(XReg).NE.MissVal) then
    OutVec (XReg) = ( (DataVec(XReg)-BaseVec(XReg)) / BaseVec(XReg)) * 100.0
  else
    OutVec (XReg) = MissVal
  end if
end do

end subroutine AnomalisePerCent

!*******************************************************************************
! QRetainMiss added to call on 25.8.00

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

real, pointer, dimension (:) 	:: TSInVec,TSLowVec,TSHighVec	! def,alloc,init in call
integer, intent (in) 		:: TimeN, THalf, QRetainMiss	! QRetainMiss:1=no,2=yes

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

real, parameter :: MissVal = -999.0

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

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

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
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 
allocate (TSPad(TimeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
TSPad = 0.0

do XTime = 1, TimeN
  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
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
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
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)

end subroutine GaussSmooth

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

end module TransformMod