! homoegeneity.f90
! written by Tim Mitchell
! module contains routines to carry out homoegenity testing
! based upon Int J Clim 17(1):25-34 (1997) Appendix 4

module Homogeneity

use Regress

implicit none

contains

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

subroutine CruTSTestMon (MonScore,QPassFail,Statistic,DataO,DataD,Order,Suffix,&
			BegN,EndN,BegO,EndO,BegD,EndD,CandO,CandD,Factor, &
			XBreakYear,SmallOrder)

real, pointer, dimension (:,:) 		:: MonScore	! NYear,NMonth
real, pointer, dimension (:,:) 		:: DataR	! NRStn,NYear
real, pointer, dimension (:) 		:: DataC	! NYear
real, pointer, dimension (:) 		:: SubOut	! NYear
real, pointer, dimension (:) 		:: AnnScore	! NYear
real, dimension (12)			:: MaxValue

integer, pointer, dimension (:,:,:)	:: DataO,DataD
integer, pointer, dimension (:)		:: Order

logical, pointer, dimension (:)		:: SourceBreak

real, intent(in), optional	:: Factor	! default=1;>1 makes it easier to pass test
real, intent(out)		:: Statistic	! the test stat at failure year

integer, intent(in), optional	:: SmallOrder	! Order size is not NDStn, but NRStn
integer, intent(in), optional	:: XBreakYear	! force exam this year only (1=first in common per)
integer, intent(in)		:: BegO,EndO	! beg/end indices of common per in O
integer, intent(in)		:: BegD,EndD	! beg/end indices of common per in D
integer, intent(in)		:: BegN,EndN    ! b/e i of nml per within the common per
integer, intent(in)		:: CandO	! candidate stn index (from DataO)
integer, intent(in), optional	:: CandD	! candidate stn index (from DataD) to merge
integer, intent(out)		:: QPassFail	! -999=no-test,0=pass,>0=fail-at-this-year
integer, dimension (12)		:: MaxIndex

character (len=4), intent(in)	:: Suffix

real, parameter, dimension (12) :: SigStat = [11.2,8.9,7.1,5.6,4.4,3.5,2.9,2.4,2.2,2.1,2.0,2.0]
real, parameter 		:: MissVal = -999.0
integer, parameter 		:: DataMissVal = -9999

real :: OpTot,OpEn,SigFactor

integer :: ReadStatus,AllocStat, QCurrent,QTestYear,QTestMonth
integer :: NYear,NDStn,NRStn,NOYear,NDYear,NMonth,NTestMonth
integer :: XYear,XDStn,XRStn,XOYear,XDYear,XMonth,XTestMonth
integer :: AnomType,Break,Comm0,Comm1

!*************************************** introduction

!write (99,"(6i6)"), BegN,EndN,BegO,EndO,BegD,EndD	! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

AnomType = 0
if (Suffix.EQ.".pre".OR.Suffix.EQ."prc".OR.Suffix.EQ.".prd" &
		.OR.Suffix.EQ.".prp".OR.Suffix.EQ.".prr") AnomType = 1

if (present(Factor)) then
  SigFactor = Factor
else
  SigFactor = 1
end if

NDStn = size(DataD,3) ; NMonth = 12

if (present(SmallOrder)) then
 NRStn = size(Order,1)
else
 NRStn = 0 ; XDStn = 0
 do
  XDStn = XDStn + 1
  if (Order(XDStn).NE.MissVal) NRStn = NRStn + 1
  if (XDStn.EQ.NDStn.OR.Order(XDStn).EQ.MissVal) exit
 end do
end if

if (present(CandD)) then
  NYear  = size(DataD,1)
  NOYear = size(DataO,1)
else
  NYear = EndO - BegO + 1
end if

if (EndN.NE.MissVal.AND.BegN.NE.MissVal) then
  Comm0=BegN ; Comm1=EndN
else
  if (NYear.LT.30) then
    Comm0=1 ; Comm1=NYear
  else
    Comm0=NYear-29 ; Comm1=NYear
  end if
end if

allocate (MonScore(NYear,NMonth),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CruTSTestMon: Allocation failure #####"
MonScore = MissVal

!*************************************** MAIN (don't bother if potential baseline < 10)

if ((Comm1-Comm0).GE.9) then

 allocate (DataR(NRStn,NYear), &
	   DataC(NYear), &
	   SubOut(NYear), &
	   AnnScore(NYear),stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: CruTSTestMon: Allocation failure #####"
 AnnScore = MissVal ; MaxIndex = MissVal ; MaxValue = MissVal

!*************************************** compare original and database series

 if (present(CandD)) then
!  write (99,*), "compare original and database series"	! @@@@@@@@@@@@@@@@@@@@@@@@
  do XMonth = 1, NMonth
    DataC=MissVal ; DataR=MissVal
    
    do XDYear = 1, NYear		! fill arrays
      XYear = XDYear
      XOYear=XDYear+BegO-BegD

      if      (DataD(XDYear,XMonth,CandD).NE.DataMissVal) then
        DataC(XYear) = real(DataD(XDYear,XMonth,CandD))/10
      else if (XOYear.GE.1.AND.XOYear.LE.NOYear) then
        if (DataO(XOYear,XMonth,CandO).NE.DataMissVal) &
        		DataC(XYear) = real(DataO(XOYear,XMonth,CandO))/10
      end if
      
      if (DataC(XYear).NE.MissVal) then
        do XRStn = 1, NRStn
          if (Order(XRStn).NE.CandD.AND.DataD(XDYear,XMonth,Order(XRStn)).NE.DataMissVal) &
          		DataR(XRStn,XYear) = real(DataD(XDYear,XMonth,Order(XRStn)))/10
        end do
      end if
    end do
    
    SubOut=MissVal			! carry out test (not BreakVec - see below score calc)
    if (present(XBreakYear)) then
    	call MultiSiteTest (DataC,DataR,Comm0,Comm1,AnomType,QPassFail,TestVec=SubOut,XBreakYear=XBreakYear)
    else
    	call MultiSiteTest (DataC,DataR,Comm0,Comm1,AnomType,QPassFail,TestVec=SubOut)
    end if
    MaxIndex(XMonth)=QPassFail
    if (MaxIndex(XMonth).GT.0) MaxValue(XMonth)=SubOut(MaxIndex(XMonth))
    
    do XYear = 1, NYear			! store monthly scores
      MonScore (XYear,XMonth) = SubOut (XYear)
    end do
  end do
 end if

!*************************************** examine original series only

 if (.not.present(CandD)) then
!  write (99,*), "examine original series only"	! @@@@@@@@@@@@@@@@@@@@@@@@
  do XMonth = 1, NMonth
!    write (99,*), "month=",xmonth	! @@@@@@@@@@@@@@@@@@@@@@@@
    DataC=MissVal ; DataR=MissVal
    
    do XYear = 1, NYear			! fill arrays
      XOYear=BegO+XYear-1 ; XDYear=BegD+XYear-1
      if (DataO(XOYear,XMonth,CandO).NE.DataMissVal) then
        DataC(XYear) = real(DataO(XOYear,XMonth,CandO))/10
        do XRStn = 1, NRStn
          if (DataD(XDYear,XMonth,Order(XRStn)).NE.DataMissVal) &
          		DataR(XRStn,XYear) = real(DataD(XDYear,XMonth,Order(XRStn)))/10
        end do
      end if
    end do
    
    SubOut=MissVal			! carry out test
!    write (99,*), "call MultiSiteTest"	! @@@@@@@@@@@@@@@@@@@@@@@@
    call MultiSiteTest (DataC,DataR,Comm0,Comm1,AnomType,QPassFail,TestVec=SubOut)
!    write (99,*), "done MultiSiteTest",QPassFail	! @@@@@@@@@@@@@@@@@@@@@@@@
    MaxIndex(XMonth)=QPassFail
    if (MaxIndex(XMonth).GT.0) MaxValue(XMonth)=SubOut(MaxIndex(XMonth))
    
    do XYear = 1, NYear			! store monthly scores
      MonScore (XYear,XMonth) = SubOut (XYear)
    end do
!    write (99,*), "done month=",xmonth	! @@@@@@@@@@@@@@@@@@@@@@@@
  end do
!  write (99,*), "done original series only"	! @@@@@@@@@@@@@@@@@@@@@@@@
 end if

!*************************************** calc scores

!  write (99,*), "calc scores"	! @@@@@@@@@@@@@@@@@@@@@@@@
 QPassFail = 0 ; Statistic = MissVal ; QTestMonth = 0
 do XMonth = 1, NMonth
  if (MaxIndex(XMonth).EQ.0) then	! no inhomog found for month 
   QTestMonth = QTestMonth + 1
  else if (MaxIndex(XMonth).GT.0) then	! inhomog found for month
   OpEn = 0 ; OpTot = 0 ; QTestMonth = QTestMonth + 1
   
   do XTestMonth = 1, NMonth		! find average peak within +-5 year range
    if (MaxIndex(XTestMonth).GT.0) then
     if (abs(MaxIndex(XMonth)-MaxIndex(XTestMonth)).LE.5) then
       OpEn = OpEn + 1 ; OpTot = OpTot + abs(MaxValue(XTestMonth))
     end if
    end if
   end do
   
   if (OpEn.GT.0) then			! find highest sig by calendar month
     if ((OpTot/OpEn).GT.(SigFactor*SigStat(OpEn))) then  ! @@@@ NOTE SigFactor @@@@@@@
       if (((OpTot/OpEn)/SigStat(OpEn)).GT.Statistic) then
         Statistic = (OpTot/OpEn)/SigStat(OpEn)
         QPassFail = MaxIndex(XMonth)
       end if
     end if
   end if
  end if
 end do
					! if few months tested and no inhomog found
 if (QTestMonth.LT.3.AND.QPassFail.EQ.0) QPassFail = MissVal

 deallocate (DataR,DataC,SubOut,AnnScore,stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: CruTSTestMon: Deallocation failure #####"

!*************************************** INSUFFICIENT RECORD FOR CHECK

else

 QPassFail = MissVal ; Statistic = MissVal
 
end if

end subroutine CruTSTestMon

!*******************************************************************************
! C = candidate t-s to be tested for homogeneity
! R = reference t-s from adjacent sites

subroutine MultiSiteTest (DataC,DataR,BegN,EndN,AnomType,QPassFail,Break,TestVec,BreakVec,XBreakYear)

real, pointer, dimension (:,:) 		:: DataR	! NRStn,NYear
real, pointer, dimension (:) 		:: DataC	! NYear
real, pointer, dimension (:) 		:: DifferencesR	! NYear
real, pointer, dimension (:) 		:: DifferencesC	! NYear
real, pointer, dimension (:) 		:: Anomalies	! NYear
real, pointer, dimension (:) 		:: Correlation	! NRStn
real, pointer, dimension (:) 		:: BaselineR	! NRStn
real, pointer, dimension (:), optional 	:: TestVec	! NYear: test stats: requires alloc/init

logical, pointer, dimension (:), optional :: BreakVec	! NYear: only test where .TRUE.
logical, pointer, dimension (:)           :: TestYear	! NYear: only test where .TRUE.

integer, intent(out) 		:: QPassFail		! -999=untested,0=pass,>0=fail-year
integer, intent(in), optional	:: Break		! only potential break point to test
integer, intent(in), optional	:: XBreakYear		! force exam this year only (1=first in common per)
integer, intent(in) 		:: BegN,EndN		! the ideal indices of the beg/end years
integer, intent(in) 		:: AnomType		! 0=abs,1=rel

real, parameter :: MissVal = -999.0

real :: Aye,Bee, OpTot,OpTotSq,OpEn,OpThresh, OpDenom,OpNumer, OpMean,OpStDev
real :: BaselineC, TestRatio,MaxRatio,CorrSum

integer :: AllocStat,ReadStatus,QNoTest,QPrescribed,QFirstA,QLastA,Comm0,Comm1,NormLen
integer :: NRStn,NYear
integer :: XRStn,XYear

!*************************************** introduction

NRStn = size(DataR,1) ; NYear = size(DataR,2) ; QPassFail=-999 ; QNoTest=0

allocate (Correlation (NRStn), &
	  DifferencesR(NYear), &
	  DifferencesC(NYear), &
	  BaselineR   (NRStn), &
	  Anomalies   (NYear), &
	  TestYear    (NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MultiSiteTest: Allocation failure #####"
Correlation = MissVal ; DifferencesC = MissVal ; DifferencesR = MissVal
BaselineR = MissVal ; Anomalies = MissVal ; TestYear = .FALSE.

!*************************************** calc baseline values

!write (99,*), "calc normal period"	! @@@@@@@@@@@@@@@@@@@@@@@@
BaselineC=MissVal ; QPrescribed = 0 
NormLen = EndN-BegN+1 ; OpThresh = 10+((NormLen-15)*0.8)
if (OpThresh.LT.10) OpThresh = 10
Comm0=BegN ; Comm1=EndN
!write (99,"(a,4i6)"), "BEGIN",BegN,EndN,NormLen,nint(OpThresh)	! @@@@@@@@@@@@@@@@@@@@@@

do 
 OpTot=0 ; OpEn=0
 do XYear = Comm0,Comm1			! calc stats for period
  if (DataC(XYear).NE.MissVal) then
    OpEn=OpEn+1 ; OpTot=OpTot+DataC(XYear)
  end if
 end do
 
 if (OpEn.GE.OpThresh) then		! enough valid vals, so calc baseline
   BaselineC = OpTot / OpEn
 else					! else
   if (QPrescribed.EQ.0) then	! the prescribed period fails
     Comm1 = NYear ; Comm0 = Comm1-NormLen+1
     QPrescribed = 1
   else				! the tested unprescribed period fails
     Comm1 = Comm1 - 1 ; Comm0 = Comm0 - 1
   end if
!   write (99,"(a,2i6)"), "SHIFT", Comm0,Comm1	! @@@@@@@@@@@@@@@@@@@@@@
 end if
 
 if (Comm0.LT.1.AND.BaselineC.EQ.MissVal) then	! no more periods of this length to test
   if (NormLen.GE.20) then			! but we can test shorter periods
     NormLen = NormLen-5 ; OpThresh = OpThresh-4
     Comm1=EndN ; Comm0=EndN-NormLen+1 ; QPrescribed=0
   else						! and we cannot shorten the period any further
     NormLen = MissVal
   end if
!   write (99,"(a,4i6)"), "REVIS",Comm0,Comm1,NormLen,nint(OpThresh)	! @@@@@@@@@@@@@@@@@@@@@@
 end if
 
 if (BaselineC.NE.MissVal.OR.NormLen.EQ.MissVal) exit
end do
if (BaselineC.EQ.MissVal) then
  QNoTest = 1
end if

if (QNoTest.EQ.0) then 
 do XRStn = 1, NRStn
  OpTot=0 ; OpEn=0
  do XYear = Comm0,Comm1
    if (DataR(XRStn,XYear).NE.MissVal) then
      OpEn=OpEn+1 ; OpTot=OpTot+DataR(XRStn,XYear)
    end if
  end do
  if (OpEn.GE.OpThresh) BaselineR(XRStn) = OpTot / OpEn
 end do
end if

!*************************************** calc correlation coefficients

if (QNoTest.EQ.0) then
! write (99,*), "calc correlation coefficients"	! @@@@@@@@@@@@@@@@@@@@@@@@
 do XYear = 2, NYear				! calc C difference series
    if (DataC(XYear).NE.MissVal.AND.DataC((XYear-1)).NE.MissVal) &
    		DifferencesC(XYear) = DataC(XYear)-DataC((XYear-1))
 end do

 CorrSum = 0
 do XRStn = 1, NRStn		! iterate by R stn
  if (BaselineR(XRStn).NE.MissVal) then
   DifferencesR = MissVal
   do XYear = 2, NYear				! calc R difference series
    if (DataR(XRStn,XYear).NE.MissVal.AND.DataR(XRStn,(XYear-1)).NE.MissVal) &
    		DifferencesR(XYear) = DataR(XRStn,XYear)-DataR(XRStn,(XYear-1))
   end do
   call LinearLSRVec (DifferencesC,DifferencesR,Aye,Bee,Correlation(XRStn))	! calc correlation
   if (Correlation(XRStn).GT.0) CorrSum = CorrSum + Correlation(XRStn)
  end if
 end do
 if (CorrSum.LT.2) then
   QNoTest = 1	! require decent regional series to test homogeneity
 end if
end if

!*************************************** calc weighted anomalies

if (QNoTest.EQ.0) then
! write (99,*), "calc weighted anomalies"	! @@@@@@@@@@@@@@@@@@@@@@@@
 do XYear = 1, NYear
  OpNumer=0 ; OpDenom=0 ; CorrSum=0
  
  if (DataC(XYear).NE.MissVal) then
    do XRStn = 1, NRStn
      if (DataR(XRStn,XYear).NE.MissVal.AND.BaselineR(XRStn).NE.MissVal.AND. &
      				Correlation(XRStn).GT.0) then
        if      (AnomType.EQ.0) then
          OpNumer = OpNumer + ((Correlation(XRStn) ** 2) * (DataR(XRStn,XYear) - &
          		BaselineR(XRStn)))
        else if (AnomType.EQ.1.AND.BaselineR(XRStn).NE.0) then
          OpNumer = OpNumer + (((Correlation(XRStn) ** 2) * DataR(XRStn,XYear)) / &
          		BaselineR(XRStn))
        end if
        
        OpDenom = OpDenom + (Correlation(XRStn) ** 2)
        CorrSum = CorrSum +  Correlation(XRStn)
      end if
    end do
    if (OpDenom.NE.0.AND.CorrSum.GT.2) then
      if      (AnomType.EQ.0) then
        Anomalies(XYear) = DataC(XYear) - BaselineC - (OpNumer/OpDenom)
      else if (AnomType.EQ.1.AND.BaselineC.NE.0.AND.OpNumer.NE.0) then
        Anomalies(XYear) = (DataC(XYear)/BaselineC) / (OpNumer/OpDenom)
      end if
    end if
  end if
 end do
end if

!*************************************** standardise weighted anomaly series

!if (QNoTest.EQ.0) then
!! write (99,*), "standardise"	! @@@@@@@@@@@@@@@@@@@@@@
! OpEn=0 ; OpTot=0 ; OpTotSq=0
! do XYear = 1, NYear
!  if (Anomalies(XYear).NE.MissVal) then
!    OpEn=OpEn+1 ; OpTot=OpTot+Anomalies(XYear) ; OpTotSq=OpTotSq+(Anomalies(XYear)**2)
!  end if
! end do
! if (OpEn.GT.1) then
!  OpMean=OpTot/OpEn 
!  OpStDev=sqrt((OpEn/(OpEn-1))*((OpTotSq/OpEn)-(OpMean**2)))
!  do XYear = 1, NYear
!    if (Anomalies(XYear).NE.MissVal) Anomalies(XYear) = (Anomalies(XYear) - OpMean) / OpStDev
!  end do
! else
!  QNoTest = 1
! end if
!end if

!*************************************** decide which years to test

if (QNoTest.EQ.0) then
!  write (99,*), "decide which years to test"	! @@@@@@@@@@@@@@@@@@@@@@@@
  if (present(Break)) then
    TestYear(Break) = .TRUE.
  else if (present(XBreakYear)) then
    TestYear(XBreakYear) = .TRUE.
  else if (present(BreakVec)) then
    do XYear = 1, NYear
      if (BreakVec(XYear).EQ..TRUE.) TestYear(XYear) = .TRUE.
    end do
  else
    QFirstA = MissVal ; QLastA = MissVal	! don't consider first5 / last5
    XYear = 0
    do
      XYear = XYear + 1
      if (Anomalies(XYear).NE.MissVal) QFirstA = XYear
      if (QFirstA.NE.MissVal.OR.XYear.EQ.NYear) exit
    end do
    XYear = NYear + 1
    do
      XYear = XYear - 1
      if (Anomalies(XYear).NE.MissVal) QLastA = XYear
      if (QLastA.NE.MissVal.OR.XYear.EQ.1) exit
    end do
    
    if ((QLastA-5).GE.(QFirstA+5)) then
      do XYear = (QFirstA+5), (QLastA-5)
        TestYear(XYear) = .TRUE.
      end do
    end if
  end if
end if

!*************************************** test for single shift

if (QNoTest.EQ.0) then
!  write (99,*), "test for single shift"	! @@@@@@@@@@@@@@@@@@@@@@@@
  QPassFail = MissVal ; MaxRatio = 0.0
  
  do XYear = 2, NYear-1
   if (TestYear(XYear).EQ..TRUE.) then
    call SingleShift (Anomalies,XYear,TestRatio,Simple=1)	

    if (present(TestVec)) TestVec(XYear) = TestRatio
    if (TestRatio.NE.MissVal) then
      if (QPassFail.EQ.MissVal) QPassFail = 0
      if (abs(TestRatio).GE.MaxRatio) then
      	MaxRatio = abs(TestRatio) ; QPassFail = XYear
      end if
    end if    

   end if
  end do
  
  if (MaxRatio.LT.2.AND.QPassFail.NE.MissVal) QPassFail = 0
end if

!*************************************** deallocate

deallocate (Correlation,DifferencesR,DifferencesC,BaselineR,Anomalies,TestYear,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MultiSiteTest: Deallocation failure #####"

end subroutine MultiSiteTest

!*******************************************************************************
! returns the result of the homg test for the alternative hypothesis that there is an
!	inhomogeneity in the series at the specified break point
! an inhomogeneity is detected if abs(TestRatio)>1

subroutine SingleShift (Series,Break,TestRatio,Simple)

real, pointer, dimension (:)	:: Series
integer, intent (in)		:: Break
real, intent(out)		:: TestRatio
integer, intent (in), optional	:: Simple

real, parameter :: MissVal = -999.0
real :: Sum1,Sum2,SumSq1,SumSq2,StDev1,StDev2,Miss1,Miss2,Valid1,Valid2
integer :: En,XVal

TestRatio = MissVal
En = size(Series)
Sum1=0 ; Sum2=0 ; SumSq1=0 ; SumSq2=0 ; Miss1=0 ; Miss2=0 

do XVal = 1, Break
  if (Series(XVal).NE.MissVal) then
    Sum1 = Sum1 + Series(XVal)
    SumSq1 = SumSq1 + (Series(XVal)**2)
  else
    Miss1 = Miss1 + 1
  end if
end do

do XVal = Break+1, En
  if (Series(XVal).NE.MissVal) then
    Sum2 = Sum2 + Series(XVal)
    SumSq2 = SumSq2 + (Series(XVal)**2)
  else
    Miss2 = Miss2 + 1
  end if
end do

Valid1 = Break-Miss1 ; Valid2 = En-Break-Miss2

if (Valid1.GE.5.AND.Valid2.GE.5) then
 if (present(Simple)) then
  StDev1 = (Valid1/(Valid1-1))*((SumSq1/Valid1)-((Sum1/Valid1)**2))
  StDev2 = (Valid2/(Valid2-1))*((SumSq2/Valid2)-((Sum2/Valid2)**2))
  
  if (StDev1.GT.0.AND.StDev2.GT.0) then
    TestRatio = ((Sum1/Valid1)-(Sum2/Valid2)) / sqrt((StDev1/Valid1)+(StDev2/Valid2))
  end if
 else
  StDev1 = sqrt((SumSq1-((Sum1**2)/Valid1))/Valid1)
  StDev2 = sqrt((SumSq2-((Sum2**2)/Valid2))/Valid2)
  
  if (StDev1.GT.0.AND.StDev2.GT.0) then
    TestRatio = (-2*Valid1*log(StDev1)) - (2*Valid2*log(StDev2)) - 1  
    TestRatio = TestRatio / ((0.3752*log(Valid2))+15.17)
  end if
 end if
end if

end subroutine SingleShift

!*******************************************************************************
! although this 'works' in the sense of running without execution errors, it does not
!	seem to do a very good job of detecting inhomogeneities, whether with the 'Simple'
!	option turned on or off. This appears to be because there is no indication within
!	the testing procedure itself as to whether the change is gradual or sudden
! I am now developing CruTSTestMon to see whether the same testing procedure can be made
!	more effective by utilising the multiple streams of information in the same year
!	to detect a simultaneous change across all months

subroutine CruTSTestAnn (QPassFail,DataO,DataD,Order,Suffix,&
			BegN,EndN,BegO,EndO,BegD,EndD,CandO,CandD)

real, pointer, dimension (:,:) 		:: DataR	! NRStn,NYear
real, pointer, dimension (:) 		:: DataC	! NYear

integer, pointer, dimension (:,:,:)	:: DataO,DataD
integer, pointer, dimension (:)		:: Order

integer, intent(in)		:: BegO,EndO	! beg/end indices of common per in O
integer, intent(in)		:: BegD,EndD	! beg/end indices of common per in D
integer, intent(in)		:: BegN,EndN    ! b/e i of nml per within the common per
integer, intent(in)		:: CandO	! candidate stn index (from DataO)
integer, intent(in), optional	:: CandD	! candidate stn index (from DataD) to merge
integer, intent(out)		:: QPassFail	! -999=no-test,0=pass,>0=fail-at-this-year

character (len=4), intent(in)	:: Suffix

real, parameter 	:: MissVal = -999.0
integer, parameter 	:: DataMissVal = -9999

integer :: ReadStatus,AllocStat, QFirstO,QFirstD,QCurrent
integer :: NYear,NDStn,NRStn,NOYear,NDYear,NMonth
integer :: XYear,XDStn,XRStn,XOYear,XDYear,XMonth
integer :: OmitD,AnomType,Break

!*************************************** introduction

if (present(CandD)) then
  OmitD = CandD
else
  OmitD = MissVal
end if

NDStn = size(DataD,3) ; NMonth = 12

NRStn = 0 ; XDStn = 0
do
  XDStn = XDStn + 1
  if (Order(XDStn).NE.MissVal) NRStn = NRStn + 1
  if (XDStn.EQ.NDStn.OR.Order(XDStn).EQ.MissVal) exit
end do

NYear = EndO - BegO + 1

allocate (DataR(NRStn,NYear), &
	  DataC(NYear),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CruTSTest: Allocation failure #####"
DataR=0 ; DataC=0

!*************************************** create series (annual averages)

XYear = 0				! get series for original candidate
do XOYear = BegO, EndO
  XYear = XYear + 1 ; XMonth = 0
  do
    XMonth = XMonth + 1
    if (DataO(XOYear,XMonth,CandO).NE.DataMissVal) then
      DataC(XYear)=DataC(XYear)+(real(DataO(XOYear,XMonth,CandO))/10.0)
      if (XMonth.EQ.12) DataC(XYear)=DataC(XYear)/12.0
    else
      DataC(XYear)=MissVal
    end if
    if (XMonth.EQ.12.OR.DataC(XYear).EQ.MissVal) exit
  end do
end do

if (present(CandD)) then		! create merged series
 XYear=0 ; Break=0 ; QFirstO=MissVal ; QFirstD=MissVal ; QCurrent=MissVal
 do XDYear = BegD, EndD
  XYear = XYear + 1  ; XMonth = 0			
  if (DataC(XYear).EQ.MissVal) then	! don't overwrite original   
   do
    XMonth = XMonth + 1
    if (DataD(XDYear,XMonth,CandD).NE.DataMissVal) then
      DataC(XYear)=DataC(XYear)+(real(DataD(XDYear,XMonth,CandD))/10.0)
      if (XMonth.EQ.12) then
      	DataC(XYear)=DataC(XYear)/12.0 ; QCurrent=1
      end if
    else
      DataC(XYear)=MissVal
    end if
    if (XMonth.EQ.12.OR.DataC(XYear).EQ.MissVal) exit
   end do
  else
   QCurrent = 0
  end if
  
  if (DataC(XYear).NE.MissVal.AND.Break.NE.MissVal) then	! monitor break point
    if (QCurrent.EQ.0) then
      if (QFirstO.EQ.MissVal) then
        QFirstO = XYear
      else
        if (QFirstD.GT.QFirstO) Break = MissVal
      end if
    else
      if (QFirstD.EQ.MissVal) then
        QFirstD = XYear
      else
        if (QFirstD.LT.QFirstO) Break = MissVal
      end if
    end if
  end if
  
 end do
 
 if (Break.NE.MissVal) then	! find break point (MissVal indicates multiple breaks)
   if (QFirstD.EQ.MissVal.OR.QFirstO.EQ.MissVal) then
     Break = 0		! no merged series - no point testing
   else if (QFirstD.GT.QFirstO) then
     Break = QFirstD-1	! D element follows O element
   else   
     Break = QFirstO-1	! opposite is true
   end if
 end if
else
 Break = MissVal
end if

do XRStn = 1, NRStn
 if (Order(XRStn).NE.OmitD) then	! don't use merged stn as ref stn
  XYear = 0
  do XDYear = BegD, EndD
    XYear = XYear + 1 ; XMonth = 0
    do
      XMonth = XMonth + 1
      if (DataD(XDYear,XMonth,Order(XRStn)).NE.DataMissVal) then
    	DataR(XRStn,XYear)=DataR(XRStn,XYear)+(real(DataD(XDYear,XMonth,Order(XRStn)))/10.0)
    	if (XMonth.EQ.12) DataR(XRStn,XYear)=DataR(XRStn,XYear)/12.0
      else
        DataR(XRStn,XYear)=MissVal
      end if
      if (XMonth.EQ.12.OR.DataR(XRStn,XYear).EQ.MissVal) exit
    end do
  end do
 else
  do XYear = 1, NYear
    DataR(XRStn,XYear) = MissVal
  end do
 end if
end do

!*************************************** run test

AnomType = 0
if (Suffix.EQ.".pre".OR.Suffix.EQ."prc".OR.Suffix.EQ.".prd" &
		.OR.Suffix.EQ.".prp".OR.Suffix.EQ.".prr") AnomType = 1

if (BegN.NE.MissVal.AND.EndN.NE.MissVal.AND.(EndN-BegN).GE.21) then
 if (Break.EQ.MissVal) then
  call MultiSiteTest (DataC,DataR,BegN,EndN,AnomType,QPassFail)
 else if (Break.NE.0) then
  call MultiSiteTest (DataC,DataR,BegN,EndN,AnomType,QPassFail,Break=Break)
 else
  QPassFail = MissVal
 end if
else
 QPassFail = MissVal
end if

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

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

end subroutine CruTSTestAnn

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

end module Homogeneity