! ghcndiscon.f90
! written by Tim Mitchell in April 2003
! module to correct discontinuities using GHCN method
! contains: MakeContinuous,MakeVersus,FilterOrig,etc

module GHCNdiscon

use Time
use FileNames
use PlainPerm
use Regress
use FDist_K
use FDist_BR
use FDist
use TTest

implicit none

contains

!*******************************************************************************
! Orig is the original data-set, including both stns to be checked (Check) and all
!	relevant stns already in the database
! if Adj is absent, Orig is returned untouched, but with discon in Disc
!	else if present, Adj has best estimate (missing where unchecked)
!	and Orig has the best estimate in periods (+-5yrs) unchecked
! Differ: .TRUE. means (y+1) minus y, .FALSE. (y+1) divide by y
! if set, AdjStn and AdjTS indicate the amount of adjustment made, in standard
!	deviations (from the original climate data), in the form of a RMSE
!	AdjStn gives the RMSE for each stn, AdjTS gives the RMSE for each year/mon

subroutine MakeContinuous (Orig,RefTS,Check,Disc,YearAD,Differ,Verbose, &
			   AdjStn,AdjTS,Adj)

logical, pointer, dimension (:,:)	:: Disc		! discon? (NY,NS)
logical, pointer, dimension (:)		:: Check	! check? (NS)
logical, pointer, dimension (:)		:: Suspects	! suspected discon (NY)
logical, pointer, dimension (:)		:: Breaks	! confirmed discon (NY)

integer, pointer, dimension (:,:,:)	:: Orig		! full climate data-set
integer, pointer, dimension (:,:,:)	:: RefTS	! reference time-series
integer, pointer, dimension (:,:,:), optional	:: Adj	! adjusted climate data-set
integer, pointer, dimension (:)		:: YearAD	! only for verbosity
integer, dimension (12)			:: Buffer	

real, pointer, dimension (:,:)		:: Best		! current best est (NY,NM)
real, pointer, dimension (:,:)		:: Versus	! Best minus(div-by) RefTS
real, pointer, dimension (:,:)		:: OneAdj	! adjust for stn (NY,NM)
real, pointer, dimension (:,:)		:: AllSum	! adj SE sum (NY,NM)
real, pointer, dimension (:,:)		:: AllEnn	! adj SE enn(NY,NM)
real, pointer, dimension (:,:),optional	:: AdjTS	! ts of adjustments
real, pointer, dimension (:),optional	:: AdjStn	! mean adjustment by stn

logical, intent (in) 	:: Differ	! see above
logical, intent (in) 	:: Verbose	! .TRUE. for verbosity

real, parameter		:: MonFract = 0.2	! min months for flag (0-1)
real, parameter 	:: MissVal = -999.0
integer, parameter 	:: DataMissVal = -9999

logical :: QEnd,QRatio
real	:: OneRMSE
integer :: AllocStat
integer :: NVal,NYear,NMonth,NStn,NBreak,NSuspect,NDisc
integer :: XVal,XYear,XMonth,XStn,XBreak

QRatio=.TRUE. ; if (Differ) QRatio=.FALSE.
NYear=size(Orig,1) ; NMonth=12 ; NStn=size(Orig,3)
if (size(RefTS,1).NE.NYear.OR.size(RefTS,2).NE.12.OR.size(RefTS,3).NE.NStn) &
		print*, "  > ##### MakeContinuous: RefTS* #####"
if (size(Disc,1).NE.NYear.OR.size(Disc,2).NE.NStn.OR.size(Check,1).NE.NStn) &
		print*, "  > ##### MakeContinuous: Disc/Check #####"

allocate (Best(NYear,NMonth), Versus(NYear,NMonth), &
	  OneAdj(NYear,NMonth), AllSum(NYear,NMonth), AllEnn(NYear,NMonth), &
	  Suspects(NYear), Breaks(NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeContinuous: Alloc #####"
Disc=.FALSE. ; AllSum=0 ; AllEnn=0 ; NDisc=0

do XStn = 1, NStn
  if (Check(XStn)) then
    if (Verbose) then
      write (99,*), ""
      write (99,"(a8,i4,a)"), "Station ", XStn, " --------"
    end if
    
    Best = MissVal
    do XYear = 1, NYear		! fill Best
      do XMonth = 1, NMonth
        if (	Orig(XYear,XMonth,XStn).NE.DataMissVal.AND. &
        	RefTS(XYear,XMonth,XStn).NE.DataMissVal) &
        	Best(XYear,XMonth) = real(Orig(XYear,XMonth,XStn))
      end do
    end do
    
    do
!     write (99,"(a)"), "call MakeVersus..."		! ###################
     call MakeVersus (Versus,Best,RefTS,QRatio,NYear,XStn)		! Versus based on Best
!     write (99,"(a)"), "call IdentifySuspects..."		! ###################
     call IdentifySuspects (Versus,Suspects,YearAD,Verbose,MonFract,QRatio)	! get suspects
     NSuspect = count(Suspects)
     
     if (NSuspect.GT.0) then
!       write (99,"(a)"), "call IdentifyBreaks..."		! ###################
       call IdentifyBreaks (Versus,Suspects,Breaks,YearAD,Verbose,&
       				MonFract,QRatio) ! find definite breaks
       NBreak = count(Breaks)
     
       if (NBreak.GT.0) then
!         write (99,"(a)"), "call UpdateDisc..."		! ###################
         call UpdateDisc (Disc,Breaks,XStn,QEnd,YearAD,Verbose,QRatio)	! update disc list
         if (.not.QEnd) then
!           write (99,"(a)"), "call CorrectStep..."		! ###################
           call CorrectStep(Disc,Orig,RefTS,XStn,Best,QRatio,YearAD)	! revise best estimate
           if (Verbose) write (99,*), ""
         end if
       end if
     end if
     
     if (QEnd.OR.NSuspect.EQ.0.OR.NBreak.EQ.0) exit
    end do
    
    if (present(AdjStn).OR.present(AdjTS)) then		! measure adj to ts
!      write (99,"(a)"), "call MeasureAdjOne..."		! ###################
      call MeasureAdjOne (Orig,Best,OneAdj,OneRMSE,XStn)
      
      if (present(AdjStn)) AdjStn(XStn)=OneRMSE
      if (present(AdjTS)) then
        do XYear = 1, NYear		
          do XMonth = 1, NMonth
            if (OneAdj(XYear,XMonth).NE.MissVal) then
              AllSum(XYear,XMonth)=AllSum(XYear,XMonth)+OneAdj(XYear,XMonth)
              AllEnn(XYear,XMonth)=AllEnn(XYear,XMonth)+1
            end if
          end do
        end do
      end if
    end if
    
    if (present(Adj)) then
      do XYear = 1, NYear		! fill Adj where Orig has been checked	
        do XMonth = 1, NMonth
          if      (RefTS(XYear,XMonth,XStn).NE.DataMissVal.AND. &
          	   Best(XYear,XMonth).NE.MissVal) &
            		Adj(XYear,XMonth,XStn) = nint(Best(XYear,XMonth))
        end do
      end do
    else
      	! do nowt - Disc is sufficient
    end if
  end if
end do

if (present(AdjTS)) then		! create final ts of adj made	
  do XYear=1,NYear
    do XMonth=1,NMonth
      if (AllEnn(XYear,XMonth).GT.0) then
        AdjTS(XYear,XMonth)=(AllSum(XYear,XMonth)/AllEnn(XYear,XMonth))
        if (AdjTS(XYear,XMonth).GT.0.001) then
          AdjTS(XYear,XMonth)=sqrt(AdjTS(XYear,XMonth))
        else
          AdjTS(XYear,XMonth)=0
        end if
      end if
    end do
  end do
end if

do XYear=1,NYear
  do XStn=1,NStn
    if (Disc(XYear,XStn)) NDisc=NDisc+1
  end do
end do

deallocate (Best,Versus,OneAdj,AllSum,AllEnn,Suspects,Breaks, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeContinuous: Dealloc #####"

end subroutine MakeContinuous

!*******************************************************************************
! measure the size of the adjustment made for all stns
! this is used by homogiter.f90 (directly)
! note that Raw is the array from the Raw file loaded, and =/NStn
! the adj stats are the RMSE, where the error is the adjustment made, in units
!	given by the stdev for that calendar month from the RefTS

subroutine MeasureAdjAll (Raw,Adj,RefTS,AdjMon,AdjStn,Examine)

logical, pointer, dimension (:)		:: Examine
integer, pointer, dimension (:,:,:)	:: Raw,Adj,RefTS
integer, pointer, dimension (:,:)	:: AdjMonEn
real, pointer, dimension (:,:)		:: AdjMon
real, pointer, dimension (:)		:: AdjStn
real, dimension (12)			:: StDev	! RefTS stdev

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

real    :: OpEn,OpTot,OpTotSq, Diff
integer :: NYear,NMonth,NStn,XYear,XMonth,XStn, AdjStnEn,AllocStat

NYear=size(Raw,1) ; NMonth=size(Raw,2) ; NStn=size(Raw,3)

allocate (AdjMonEn(NYear,NMonth), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### PruneRaw: alloc #####"
AdjMonEn=0 ; AdjMon=0 ; AdjStn=0

do XStn=1,NStn
 if (Examine(XStn)) then		! only include the Raw stns
  AdjStnEn=0 ; StDev=MissVal

  do XMonth=1,NMonth			! calc ref ts stdev by month
   OpEn=0 ; OpTot=0 ; OpTotSq=0		! (remember: refTS has (x,s) of Raw)
   do XYear=1,NYear
    if (RefTS(XYear,XMonth,XStn).NE.DataMissVal) then
      OpEn=OpEn+1 ; OpTot=OpTot+RefTS(XYear,XMonth,XStn)
      OpTotSq=OpTotSq+(RefTS(XYear,XMonth,XStn)**2)
    end if
   end do
   if (OpEn.GT.1) then
    StDev(XMonth)=((OpTotSq/OpEn)-((OpTot/OpEn)**2))*(OpEn/(OpEn-1))
    if (StDev(XMonth).LE.0) then
      StDev(XMonth)=MissVal
    else
      StDev(XMonth)=sqrt(StDev(XMonth))
    end if
   end if
  end do
  
  do XMonth=1,NMonth			! calc the adj made at each moment			
    if (StDev(XMonth).NE.MissVal) then		
      do XYear=1,NYear				
        if (Adj(XYear,XMonth,XStn).NE.DataMissVal) then	
          Diff = ((Adj(XYear,XMonth,XStn) - Raw(XYear,XMonth,XStn)) / &
          		StDev(XMonth)) ** 2
          
          AdjStn(XStn) = AdjStn(XStn) + Diff	! calc for stn
          AdjStnEn = AdjStnEn + 1
          					! calc for time
          AdjMon(XYear,XMonth) = AdjMon(XYear,XMonth) + Diff
          AdjMonEn(XYear,XMonth) = AdjMonEn(XYear,XMonth) + 1
        end if
      end do
    end if
  end do
  
  if (AdjStnEn.GT.0) then		! final stat for stn
    AdjStn(XStn) = sqrt (AdjStn(XStn) / real(AdjStnEn))
  else
    AdjStn(XStn) = MissVal
  end if
 end if
end do

do XYear=1,NYear			! final stat for time
  do XMonth=1,NMonth
    if (AdjMonEn(XYear,XMonth).GT.0) then
      AdjMon(XYear,XMonth) = sqrt(AdjMon(XYear,XMonth) / &
      				real(AdjMonEn(XYear,XMonth)))
    else
      AdjMon(XYear,XMonth) = MissVal
    end if
  end do
end do

end subroutine MeasureAdjAll

!*******************************************************************************
! measure the size of the adjustment made for this stn, as a time-series
! this is used by homogenise.f90 (indirectly via MakeContinuous)

subroutine MeasureAdjOne (Orig,Best,OneAdj,OneRMSE,XStn)

integer, pointer, dimension (:,:,:)	:: Orig		! full climate data-set
real, pointer, dimension (:,:)		:: OneAdj	! adjust for stn (NY,NM)
real, pointer, dimension (:,:)		:: Best		! current best est (NY,NM)
real, dimension (12)			:: StDev	! original stdev
real, intent(out)  			:: OneRMSE
integer, intent(in)  			:: XStn

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

real    :: OpEn,OpTot,OpTotSq, Sum,En
integer :: NYear,NMonth,XYear,XMonth

NYear=size(Best,1) ; NMonth=size(Best,2)
StDev=MissVal ; OneAdj=MissVal ; OneRMSE=MissVal ; Sum=0 ; En=0

do XMonth=1,NMonth			! calc original stdev by month
  OpEn=0 ; OpTot=0 ; OpTotSq=0
  do XYear=1,NYear
    if (Orig(XYear,XMonth,XStn).NE.DataMissVal) then
      OpEn=OpEn+1 ; OpTot=OpTot+Orig(XYear,XMonth,XStn)
      OpTotSq=OpTotSq+(Orig(XYear,XMonth,XStn)**2)
    end if
  end do
  if (OpEn.GT.1) then
    StDev(XMonth)=((OpTotSq/OpEn)-((OpTot/OpEn)**2))*(OpEn/(OpEn-1))
    if (StDev(XMonth).LT.0) then
      StDev(XMonth)=MissVal
    else
      StDev(XMonth)=sqrt(StDev(XMonth))
    end if
  end if
end do

do XMonth=1,NMonth				! calc OneAdj made (sq abs,in stdev)
 if (StDev(XMonth).NE.MissVal) then		! as given, calcs for homog+retry
  do XYear=1,NYear				! if filtered by Orig=/DMV & RefTS=/DMV
    if (Best(XYear,XMonth).NE.MissVal) then	! this would calc for homog only
      OneAdj(XYear,XMonth) = ((Best(XYear,XMonth)-real(Orig(XYear,XMonth,XStn))) / &
      				StDev(XMonth)) ** 2
      Sum=Sum+OneAdj(XYear,XMonth) ; En=En+1
    end if
  end do
 end if
end do

if (En.GT.0) then
  OneRMSE=(Sum/En)
  if (OneRMSE.GT.0.001) then
    OneRMSE=sqrt(OneRMSE)
  else
    OneRMSE=0
  end if
end if

end subroutine MeasureAdjOne

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

subroutine CorrectStep (Disc,Orig,RefTS,XStn,Best,QRatio,YearAD)

logical, pointer, dimension (:,:)	:: Disc		! discon? (NY,NS)
integer, pointer, dimension (:,:,:)	:: Orig		! full climate data-set
integer, pointer, dimension (:,:,:)	:: RefTS	! full ref TS
integer, pointer, dimension (:)		:: YearAD	! only for verbosity
real, pointer, dimension (:,:)		:: OrigVersus	! original Versus
real, pointer, dimension (:,:)		:: Best		! current best est (NY,NM)
real, dimension (12)			:: Correct	! correction to apply
logical, intent (in) 	:: QRatio	! .TRUE. for precip etc, .FALSE. otherwise
integer, intent(in)  	:: XStn

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

logical :: AcceptComp
real 	:: OpTot,OpEn,Before,After,Log0,Log1
integer :: NYear,NMonth,NPane,XYear,XMonth,XExam,XPane, Year0,Year1,Year2,Year3
integer :: LastDisc,ThisDisc,NextDisc,Allocstat,OldLen, YearA,YearB,YearC,YearD, Elastic

NYear=size(Disc,1) ; NMonth=12 ; LastDisc=NYear+1 
LastDisc=NYear+1 ; ThisDisc=MissVal ; NextDisc=MissVal
NPane=12 ; if (QRatio) NPane=24
Elastic=5 ; if (QRatio) Elastic=10	! min len of any calib stretch

Best = MissVal
do XYear=1,NYear
  do XMonth = 1, NMonth
    if (Orig(XYear,XMonth,XStn).NE.DataMissVal) &
          	Best(XYear,XMonth) = Orig(XYear,XMonth,XStn)
  end do
end do

allocate (OrigVersus(NYear,NMonth), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### CorrectStep: Alloc #####"

do XYear = NYear,1,-1
  if (Disc(XYear,XStn)) then	! year=disc
    if (ThisDisc.EQ.MissVal) then	! first disc
      ThisDisc=XYear
    else				! later disc
      NextDisc=XYear
    end if
  else if (XYear.EQ.1.AND.ThisDisc.NE.MissVal) then
    NextDisc=0				! got suspect, but reached end of period
  end if
  
  if (NextDisc.NE.MissVal) then	! got everything for check, so check
    Year0=NextDisc+1 ; Year1=ThisDisc ; Year2=ThisDisc+1 ; Year3=LastDisc-1

    call MakeVersus (OrigVersus,Best,RefTS,QRatio,NYear,XStn)	! original Versus
    
    Correct=MissVal		! Correct(XMonth) contains the adjustment to make
    
    do XMonth = 1, NMonth	! calc corrections to make
      Before=MissVal ; After=MissVal
      YearA=Year0 ; YearB=Year1 ; YearC=Year2 ; YearD=Year3 ; AcceptComp=.FALSE.
      
      OpTot=0 ; OpEn=0
      do XExam=YearA,YearB
        if (OrigVersus(XExam,XMonth).NE.MissVal) then
          OpTot=OpTot+OrigVersus(XExam,XMonth) ; OpEn=OpEn+1
        end if
      end do
      if (OpEn.GT.0) Before=OpTot/OpEn
      
      OpTot=0 ; OpEn=0
      do XExam=YearC,YearD
        if (OrigVersus(XExam,XMonth).NE.MissVal) then
          OpTot=OpTot+OrigVersus(XExam,XMonth) ; OpEn=OpEn+1
        end if
      end do
      if (OpEn.GT.0) After=OpTot/OpEn
       
      if (Before.NE.MissVal.AND.After.NE.MissVal) then
        if (QRatio) then
          if (Before.EQ.0) then
	    Correct(XMonth)=0.0
	  else
	    Correct(XMonth)=After/Before
	  end if
        else
          Correct(XMonth)=After-Before
        end if
      end if
    end do
    
    call SmoothCorrect (Correct,QRatio) ! smooth the corr into annual cycle
!    if (.NOT.QRatio) call SmoothCorrect (Correct,QRatio)
    
    do XMonth=1,NMonth  		! make the corrections
      if (Correct(XMonth).NE.MissVal) then
        do XExam = 1,ThisDisc
          if (Best(XExam,XMonth).NE.MissVal) then
            if (QRatio) then
              Best(XExam,XMonth) = real(Best(XExam,XMonth)) * Correct(XMonth)
            else
              Best(XExam,XMonth) = real(Best(XExam,XMonth)) + Correct(XMonth)
            end if
          end if
        end do
      end if
    end do
    
    LastDisc=ThisDisc ; ThisDisc=NextDisc ; NextDisc=MissVal
  end if
end do
  
deallocate (OrigVersus, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### CorrectStep: Dealloc #####"

end subroutine CorrectStep

!*******************************************************************************
! this subroutine smooths the corrections by calendar month
! this is advisable, because the correction to make is initially determined
!	by a rather small sample for each individual month
! by smoothing we can eliminate some of the potentially erratic corrections
! we can also add estimated corrections for months without an initial calc

subroutine SmoothCorrect (Correct,QRatio)

real, pointer, dimension (:) 	:: TSInVec,TSLowVec,TSHighVec	! vectors
real, dimension (12)		:: Correct			! corr to apply

logical, intent (in) 	:: QRatio	! .TRUE. for precip etc, .FALSE. otherwise
real, parameter 	:: MissVal = -999.0

real 	:: OpTotSq,OpTot,OpEn, OrigMean,SmooMean,OrigStDev,SmooStDev
integer :: NTriple,NMonth,NThree, AllocStat
integer :: XTriple,XMonth,XThree

NMonth=12 ; NTriple=36 ; NThree=3

allocate (TSInVec(NTriple),TSLowVec(NTriple),TSHighVec(NTriple),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SmoothCorrect: alloc #####"

do XMonth=1,NMonth		! fill TSInVec with initial corrections
  do XThree=0,2
    TSInVec((XThree*12)+XMonth) = Correct(XMonth)
  end do
end do

do XMonth=1,NMonth		! in-fill any missing values
  if (Correct(XMonth).EQ.MissVal) then
    OpTot=0 ; OpEn=0
    
    do XTriple=(10+XMonth),(14+XMonth)
      if (TSInVec(XTriple).NE.MissVal) then
        OpTot=OpTot+TSInVec(XTriple) ; OpEn=OpEn+1
      end if
    end do
    
    if (OpEn.GT.1) then
      do XThree=0,2
        TSInVec((XThree*12)+XMonth) = OpTot/OpEn
      end do
    end if
  end if
end do

OpTotSq=0 ; OpTot=0 ; OpEn=0 ; OrigMean=MissVal	! find original mean
do XMonth=1,NMonth			
  if (TSInVec(12+XMonth).NE.MissVal) then
    OpEn=OpEn+1 
    OpTot=OpTot+TSInVec(12+XMonth)
    OpTotSq=OpTotSq+(TSInVec(12+XMonth)**2)
  end if 
end do
if (OpEn.EQ.12) then
  OrigMean=OpTot/OpEn 
  OrigStDev=(OpTotSq/OpEn)-(OrigMean**2)
  if (OrigStDev.GE.0) then
    OrigStDev=sqrt(OrigStDev)
  else
    OrigStDev=1
  end if
end if

call GaussSmooth (NTriple,9,1,TSInVec,TSLowVec,TSHighVec)

if (QRatio) then
    	! not possible to ensure orig and smoo corr tally
else
 OpTotSq=0 ; OpTot=0 ; OpEn=0 ; SmooMean=MissVal	! find smoothed mean
 do XMonth=1,NMonth
  Correct(XMonth) = TSLowVec(12+XMonth)

  if (TSLowVec(12+XMonth).NE.MissVal) then
    OpEn=OpEn+1 
    OpTot=OpTot+TSLowVec(12+XMonth)
    OpTotSq=OpTotSq+(TSLowVec(12+XMonth)**2)
  end if 
 end do
 if (OpEn.EQ.12) then
   SmooMean=OpTot/OpEn
   SmooStDev=(OpTotSq/OpEn)-(SmooMean**2)
   if (SmooStDev.GE.0) then
     SmooStDev=sqrt(SmooStDev)
   else
     SmooStDev=1
   end if
 end if
					! ensure orig and smoo corr tally
 if (SmooMean.NE.MissVal.AND.OrigMean.NE.MissVal) then
    do XMonth=1,NMonth
      Correct(XMonth) = Correct(XMonth) + (OrigMean-SmooMean)
      Correct(XMonth) = OrigMean + &
      			((Correct(XMonth)-OrigMean) * (OrigStDev/SmooStDev))
    end do
 end if
end if

deallocate (TSInVec,TSLowVec,TSHighVec,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SmoothCorrect: alloc #####"

end subroutine SmoothCorrect

!*******************************************************************************
! for a given XStn, identify all additional breaks (i.e. those outside the range of an 
!	already-recognised discontinuity +- 4 years)
! additionals are entered as discontinuities
! QEnd is returned .TRUE. when no additionals are found

subroutine UpdateDisc (Disc,Breaks,XStn,QEnd,YearAD,Verbose,QRatio)

logical, pointer, dimension (:,:)	:: Disc		! discon? (NY,NS)
logical, pointer, dimension (:)		:: Breaks	! confirmed discon (NY)
integer,pointer,dimension(:)		:: YearAD	! only include for verbosity
logical, intent(in)  :: Verbose
logical, intent(in)  :: QRatio	! .TRUE. for precip etc
integer, intent(in)  :: XStn
logical, intent(out) :: QEnd

logical :: QClash
integer :: NYear,XDiscYear,XBreakYear,Year0,Year1,Elastic

Elastic=5 ; if (QRatio) Elastic=10	! min len of any calib stretch
NYear=size(Breaks) ; QEnd=.TRUE.

do XBreakYear = 1, NYear
  if      (Breaks(XBreakYear)) then
    Year0=XBreakYear-Elastic+1 ; Year1=XBreakYear+Elastic+1
    if (Year0.LT.1) Year0=1
    if (Year1.GT.NYear) Year1=NYear
    
    XDiscYear=Year0-1 ; QClash=.FALSE.
    do
      XDiscYear=XDiscYear+1
      if (Disc(XDiscYear,XStn)) QClash=.TRUE.
      if (QClash.OR.XDiscYear.EQ.Year1) exit
    end do
    
    if (.NOT.QClash) then
      QEnd = .FALSE.
      Disc(XBreakYear,XStn) = .TRUE.
      if (Verbose) write (99,"(a8,i4)"), "DISCON  ", YearAD(XBreakYear)
    else
      Breaks(XBreakYear) = .FALSE.
      if (Verbose) write (99,"(a8,i4)"), "Refuse  ", YearAD(XBreakYear)
    end if
  end if
end do

end subroutine UpdateDisc

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

subroutine IdentifyBreaks (Versus,Suspects,Breaks,YearAD,Verbose,MonFract,QRatio)

logical, pointer, dimension (:)	:: Breaks	! confirmed discon (NY)
logical, pointer, dimension (:)	:: Suspects	! suspected discon (NY)
real, pointer, dimension (:,:)	:: Perm		! NPerm,NYear
real, pointer, dimension (:,:)	:: Versus	! Best minus(div-by) RefTS
real, pointer, dimension (:)	:: Single	! monthly version of Versus
integer,pointer,dimension(:)	:: YearAD	! only include for verbosity
logical, intent(in)		:: Verbose
logical, intent (in) 		:: QRatio	! .TRUE. for precip etc
real, intent(in) 		:: MonFract	! min months for flag (0-1)

real, parameter :: MissVal = -999.0

logical :: QPassMonth
real :: TestStat,PermStat
integer :: LastSuspect,ThisSuspect,NextSuspect,Year0,Year1,AllocStat,Part
integer :: XYear,NYear,XMonth,NMonth,XLimit,NPerm,XPerm
integer :: NSample,NSmaller,NThresh,NTest,NFail,NPass

NYear=size(Versus,1) ; NMonth=12 ; NPerm=200 
Breaks = .FALSE.
LastSuspect=0 ; ThisSuspect=0 ; NextSuspect=0

Part=12 ; if (QRatio) Part=24		! max len of partition

allocate (Single(NYear),Perm(NPerm,NYear),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### IdentifyBreaks: Alloc #####"

do XYear = 1, NYear
  if (Suspects(XYear)) then	! year=suspect
    if (ThisSuspect.EQ.0) then		! no 'check' year, so fill
      ThisSuspect=XYear
    else				! got 'check' year, so make this the next
      NextSuspect=XYear
    end if
  else if (XYear.EQ.NYear.AND.ThisSuspect.NE.0) then
    NextSuspect=XYear+1		! got suspect, but reached end of period
  end if
  
  if (NextSuspect.NE.0) then	! got everything for check, so check
    Year0=(LastSuspect+1) ; Year1=(NextSuspect-1) ; NTest=0 ; NPass=0
    if (Year0.LT.(ThisSuspect-Part+1)) Year0=(ThisSuspect-Part+1) 
    if (Year1.GT.(ThisSuspect+Part))   Year1=(ThisSuspect+Part)
    
    do XMonth = 1, NMonth	! check each month individually
      Single=MissVal
      do XLimit = Year0,Year1	! fill single with relevant data
        Single(XLimit)=Versus(XLimit,XMonth)
      end do
      				! calc test statistic
      call ClusterDistance (Single,Year0,ThisSuspect,(ThisSuspect+1),Year1,QRatio,TestStat)
      
      if (TestStat.NE.MissVal) then
      				! get repartitioned sample in Perm - plainperm.f90
        call PartitionSample (Single,Single,Perm,Perm, &
      			    Year0,ThisSuspect,(ThisSuspect+1),Year1)
        
        NSample=NPerm ; NSmaller=0 ; NThresh=nint(real(NSample)*0.05) 
        QPassMonth=.FALSE. ; XPerm=0
        do			! each perm...
          XPerm=XPerm+1
          
          do XLimit = Year0,Year1	! fill single with permutation
            Single(XLimit)=Perm(XPerm,XLimit)
          end do
          				! calc mean distance within clusters
          call ClusterDistance (Single,Year0,ThisSuspect,(ThisSuspect+1),Year1,QRatio,PermStat)
          
          if      (PermStat.EQ.MissVal) then
            NSample=NSample-1
            NThresh=nint(real(NSample)*0.05)
          else if (PermStat.LT.TestStat) then
            NSmaller=NSmaller+1		! smaller than test stat
            if (NSmaller.GT.NThresh) QPassMonth=.TRUE.
          end if
          
          if (XPerm.EQ.NPerm.OR.QPassMonth) exit
        end do
        
        if (NSmaller.GT.NThresh) NPass=NPass+1
        NTest=NTest+1
      end if
    end do
    
    NFail=NTest-NPass
    if (real(NFail).GT.(real(NTest)*MonFract)) then	! confirm discon
      if (Verbose) write (99,"(a8,3(i4,a),a8,2(i2,a))"), "Confirm ", YearAD(ThisSuspect), &
        	" (", YearAD(Year0), "-", YearAD(Year1), ") ", &
        	"P-test (", NFail, "/", NTest, ")"
      Breaks(ThisSuspect) = .TRUE.		
      LastSuspect=ThisSuspect ; ThisSuspect=NextSuspect ; NextSuspect=0
    else					! reject discon
      if (Verbose) write (99,"(a8,3(i4,a),a8,2(i2,a))"), "Discard ", YearAD(ThisSuspect), &
        	" (", YearAD(Year0), "-", YearAD(Year1), ") ", &
        	"P-test (", NFail, "/", NTest, ")"
      ThisSuspect=NextSuspect ; NextSuspect=0
    end if
  end if
end do

deallocate (Single,Perm,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### IdentifyBreaks: Dealloc #####"

end subroutine IdentifyBreaks

!*******************************************************************************
! given a vector with two partitions (A,B)
! this calculates the weighted mean distance between a pair within a partition

subroutine ClusterDistance (Single,BegA,EndA,BegB,EndB,QRatio,Distance)

real, pointer, dimension (:) 	:: Single 			! (NValue)
integer, intent(in) 		:: BegA,EndA,BegB,EndB		! XValue indices
logical, intent (in) 		:: QRatio			! .TRUE. for precip etc
real, intent(out) 		:: Distance

real, parameter :: MissVal = -999.0

real :: OpTot,OpEn
integer :: XMem0,XMem1,Elastic

Elastic=5 ; if (QRatio) Elastic=10	! min len of any calib stretch
OpTot=0.0 ; OpEn=0.0 ; Distance=MissVal
        
do XMem0 = BegA,(EndA-1)	! look at each pair in partition 1
         if (Single(XMem0).NE.MissVal) then
          do XMem1 = (XMem0+1),EndA
           if (Single(XMem1).NE.MissVal) then
             OpEn=OpEn+1 ; OpTot=OpTot+abs(Single(XMem0)-Single(XMem1))
           end if
          end do
         end if
end do

if (OpEn.GE.Elastic) then        
  do XMem0 = BegB, (EndB-1)	! look at each pair in partition 2
         if (Single(XMem0).NE.MissVal) then
          do XMem1 = (XMem0+1),EndB
           if (Single(XMem1).NE.MissVal) then
             OpEn=OpEn+1 ; OpTot=OpTot+abs(Single(XMem0)-Single(XMem1))
           end if
          end do
         end if
  end do
  
  if (OpEn.GE.(Elastic*2)) Distance=OpTot/OpEn
end if

end subroutine ClusterDistance

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

subroutine IdentifySuspects (Versus,Suspects,YearAD,Verbose,MonFract,QRatio)

logical, pointer, dimension (:)	:: Suspects		! suspected discon (NY)
real, pointer, dimension (:,:)	:: Versus		! Best minus(div-by) RefTS
real, pointer, dimension (:)	:: Single		! one calendar month
integer, pointer, dimension (:)	:: QueryBeg,QueryEnd	! subsections to examine
integer,pointer,dimension(:)	:: YearAD		! only include for verbosity
real, dimension (12)		:: RSS1,RSS2		! 1=whole,2=subsections	
integer, dimension (12)		:: RSS2Enn		! n in subsections
logical, intent (in) 		:: QRatio		! .TRUE. for precip etc
logical, intent (in)		:: Verbose
real, intent (in) 		:: MonFract		! min months for flag (0-1)

real, parameter :: MissVal = -999.0

real (dp) :: FTestStat,Pea,Que
real :: Upper
integer :: Suspect,QFail,CurrentBeg,CurrentEnd,LastQuery,Err,AllocStat
integer :: XMonth,NMonth,XYear,NYear,NTest,NFail

NMonth=12 ; NYear=size(Versus,1)

allocate (Single(NYear), &
	  QueryBeg(NYear),QueryEnd(NYear),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### IdentifySuspects: Alloc #####"
QueryBeg=0 ; QueryEnd=0 ; Suspects=.FALSE.
LastQuery=1 ; QueryBeg(1)=1 ; QueryEnd(1)=size(Suspects)

do
  CurrentBeg = QueryBeg(LastQuery) ; CurrentEnd = QueryEnd(LastQuery)
  call MostSuspect (Versus,CurrentBeg,CurrentEnd,Suspect,RSS1,RSS2,RSS2Enn,MonFract,&
  			YearAD,QRatio)
  
  QFail=0
  if (Suspect.NE.MissVal) then		
    NTest=0 ; NFail=0			
    do XMonth = 1, NMonth
      if (RSS2Enn(XMonth).GT.4.AND.RSS2(XMonth).GT.0) then
        FTestStat = ((RSS1(XMonth)-RSS2(XMonth))*(real(RSS2Enn(XMonth))-4.0)) &
        	   / (RSS2(XMonth)*3.0)
        call FProb(3,(RSS2Enn(XMonth)-4),FTestStat,Pea,Que,Err)		! fdist.f90
  	if (Err.EQ.0)    NTest=NTest+1
  	if (Pea.GT.0.95) NFail=NFail+1
      end if
    end do
    
    if (real(NFail).GE.(real(NTest)*MonFract)) then	! fail under F-test
      QFail = 1				
    else					! so try t-test
      NTest=0 ; NFail=0
      do XMonth = 1, NMonth
        Single=MissVal
        do XYear = CurrentBeg,CurrentEnd
          Single(XYear)=Versus(XYear,XMonth)
        end do
        call DiffMeans (Single,Single,CurrentBeg,Suspect, &
        		(Suspect+1),CurrentEnd,Upper)
        if (Upper.NE.MissVal) NTest=NTest+1
        if (Upper.GT.0.95)    NFail=NFail+1
      end do
      
      if (real(NFail).GT.(real(NTest)*MonFract)) QFail = 2	! fail under t-test
    end if
  end if
  
  if (Verbose) then
    if      (Suspect.EQ.MissVal) then
        write (99,"(a8,4x,a,2(i4,a))"), "Sparse  ", &
        	" (", YearAD(CurrentBeg), "-", YearAD(CurrentEnd), ") " 
    else if (QFail.EQ.0) then
        write (99,"(a8,3(i4,a))"), "Ignore  ", YearAD(Suspect), &
        	" (", YearAD(CurrentBeg), "-", YearAD(CurrentEnd), ") " 
    else if (QFail.EQ.1) then
        write (99,"(a8,3(i4,a),a8,2(i2,a))"), "Suspect ", YearAD(Suspect), &
        	" (", YearAD(CurrentBeg), "-", YearAD(CurrentEnd), ") ", &
        	"F-test (", NFail, "/", NTest, ")"
    else if (QFail.EQ.2) then
        write (99,"(a8,3(i4,a),a8,2(i2,a))"), "Suspect ", YearAD(Suspect), &
        	" (", YearAD(CurrentBeg), "-", YearAD(CurrentEnd), ") ", &
        	"t-test (", NFail, "/", NTest, ")"
    end if
  end if

  if (QFail.GT.0) then
    Suspects(Suspect)=.TRUE.
    QueryBeg(LastQuery)=CurrentBeg ; QueryEnd(LastQuery)=Suspect-1
    LastQuery=LastQuery+1
    QueryBeg(LastQuery)=Suspect+1 ; QueryEnd(LastQuery)=CurrentEnd
  else
    LastQuery=LastQuery-1
  end if
    
  if (LastQuery.EQ.0) exit
end do

deallocate (Single,QueryBeg,QueryEnd,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### IdentifySuspects: Dealloc #####"

end subroutine IdentifySuspects

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

subroutine MostSuspect (Versus,Beg,End,Suspect,RSS1,RSS2,RSS2Enn,MonFract,YearAD,QRatio)

real, pointer, dimension (:,:)	:: Versus	! Best Minus(div-by) RefTS
real, pointer, dimension (:,:)	:: Full		! raw version of RSS2
real, pointer, dimension (:,:)	:: Stand	! RSS2 / RSS1
real, pointer, dimension (:)	:: Exe,Wye	! NYear vectors
real, dimension (12)		:: RSS1,RSS2	! 1=whole,2=subsections	
integer,pointer,dimension (:,:)	:: FullEnn	! full n in subsections
integer,pointer,dimension(:)	:: YearAD	! only include for verbosity
integer, dimension (12)		:: RSS2Enn	! n in subsections
logical, intent (in) 		:: QRatio	! .TRUE. for precip etc, .FALSE. otherwise
integer, intent(in) 		:: Beg,End	! period of interest to query
integer, intent(out) 		:: Suspect	! most suspicious year index
real, intent(in) 		:: MonFract	! Min months for flag (0-1)

real, parameter :: MissVal = -999.0

real :: Aye,Bee,Are,Enn,OpTotSq,OpTot,OpEn,OpMean,OpStDev,RSS,Min,EnnPart1,EnnPart2
integer :: AllocStat,MinYear,Elastic
integer :: XYear,XMonth,XQuery,XTestYear
integer :: NYear,NMonth,NQuery,NTestYear

RSS1=MissVal ; RSS2=MissVal ; RSS2Enn=MissVal ; Suspect=MissVal
NYear=size(Versus,1) ; NMonth=12 ; NQuery=End-Beg+1 
Min=0-MissVal ; MinYear=MissVal
Elastic=5 ; if (QRatio) Elastic=10	! min len of any calib stretch

if (NQuery.GE.(Elastic*2)) then
 allocate (Exe(NQuery),Wye(NQuery),FullEnn(NYear,NMonth), &
 	   Full(NYear,NMonth),Stand(NYear,NMonth), stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### MostSuspect: Alloc #####"
 FullEnn=MissVal ; Full=MissVal ; Stand=MissVal
 
 do XMonth = 1, NMonth				
  OpEn=0 ; OpTot=0 ; OpTotSq=0		! initialise RSS2 distribution
  
  Exe=MissVal ; Wye=MissVal
  do XYear = Beg, End			! fill Exe,Wye
    if (Versus(XYear,XMonth).NE.MissVal) then
      XQuery = XYear-Beg+1
      Exe(XQuery) = real(XQuery)
      Wye(XQuery) = Versus(XYear,XMonth)/100.0
    end if
  end do
  					! calc RSS1
  call RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,Elastic,RSS=RSS1(XMonth),Beg=1,End=NQuery)
  
  if (RSS1(XMonth).NE.MissVal.AND.RSS1(XMonth).NE.0) then
   do XQuery = Elastic,(NQuery-Elastic)		! test each year in turn
    XTestYear=Beg+XQuery-1			! calc RSS2 part 1
    call RegressVectors (Exe,Wye,Aye,Bee,Are,EnnPart1,Elastic,RSS=RSS,Beg=1,End=XQuery)
    
    if (RSS.NE.MissVal) then
  						! calc RSS2 part 2
      call RegressVectors (Exe,Wye,Aye,Bee,Are,EnnPart2,Elastic,RSS=Full(XTestYear,XMonth), &
      			   Beg=(XQuery+1),End=NQuery)
      
      if (Full(XTestYear,XMonth).NE.MissVal) then	! calc full RSS2
      	FullEnn(XTestYear,XMonth) = nint(EnnPart1 + EnnPart2)
      	Full(XTestYear,XMonth) = Full(XTestYear,XMonth) + RSS
      	Stand(XTestYear,XMonth) =  Full(XTestYear,XMonth) / RSS1(XMonth)
      end if
    end if
   end do
  end if
 end do
 
 do XTestYear = (Beg+Elastic-1),(End-Elastic)		! find most suspicious year
   OpEn=0 ; OpTot=0
   do XMonth = 1, NMonth
     if (Stand(XTestYear,XMonth).NE.MissVal) then
       OpEn=OpEn+1 ; OpTot=OpTot+Stand(XTestYear,XMonth)
     end if
   end do
   if (real(OpEn).GE.(MonFract*12.0)) then
     if ((OpTot/OpEn).LT.Min) then
       Min=OpTot/OpEn ; MinYear=XTestYear
     end if
   end if
 end do
 
 if (MinYear.NE.MissVal) then		! return most suspicious year
   Suspect = MinYear
   do XMonth = 1, NMonth
     RSS2(XMonth) = Full(Suspect,XMonth)
     RSS2Enn(XMonth) = FullEnn(Suspect,XMonth)
   end do
 end if

 deallocate (Exe,Wye,Full,FullEnn,Stand, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### MostSuspect: Alloc #####"
end if

end subroutine MostSuspect

!*******************************************************************************
! calcs Versus for MakeContinuous on the basis of supplied info
! I'm not convinced by the algorithm for precip, esp for dry regions, but better?

subroutine MakeVersus (Versus,Best,RefTS,QRatio,NYear,XStn)

integer, pointer, dimension (:,:,:)	:: RefTS	! reference time-series

real, pointer, dimension (:,:)		:: Best		! current best est (NY,NM)
real, pointer, dimension (:,:)		:: Versus	! Best minus(div-by) RefTS

logical, intent (in) 	:: QRatio	! .TRUE. for precip etc, .FALSE. otherwise
integer, intent (in) 	:: XStn		! stn required from RefTS
integer, intent (in) 	:: NYear	! no. years in V,B,R

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

integer :: XMonth,XYear

Versus = MissVal
if (QRatio) then
  do XYear = 1, NYear
    do XMonth = 1, 12
      if (Best(XYear,XMonth).NE.MissVal.AND.RefTS(XYear,XMonth,XStn).NE.DataMissVal) then
        if      (Best(XYear,XMonth).NE.0.AND.RefTS(XYear,XMonth,XStn).NE.0) then
          Versus(XYear,XMonth) = Best(XYear,XMonth) / real(RefTS(XYear,XMonth,XStn))
        else if (Best(XYear,XMonth).NE.0) then
          Versus(XYear,XMonth) = 10.0		! value / zero
        else if (RefTS(XYear,XMonth,XStn).NE.0) then
          Versus(XYear,XMonth) =  0.1		! zero / value
        else
          Versus(XYear,XMonth) =  1.0		! zero / zero
	end if
      end if
    end do
  end do
else
  do XYear = 1, NYear
    do XMonth = 1, 12
      if (Best(XYear,XMonth).NE.MissVal.AND.RefTS(XYear,XMonth,XStn).NE.DataMissVal) &
         Versus(XYear,XMonth) = Best(XYear,XMonth) - real(RefTS(XYear,XMonth,XStn))
    end do
  end do
end if

end subroutine MakeVersus

!*******************************************************************************
! take the three main arrays (Orig,RefTS,Adj)
! use RefTS to filter which data goes in homog and which in retry

subroutine FilterOrig (Orig,RefTS,Adj,Check,Differ)

logical, pointer, dimension (:)		:: Check,Retain

integer, pointer, dimension (:,:,:)	:: Orig		! full climate data-set
integer, pointer, dimension (:,:,:)	:: RefTS	! reference time-series
integer, pointer, dimension (:,:,:)	:: Adj		! adjusted climate data-set

logical, intent (in) 	:: Differ	! .FALSE. for precip etc

integer, parameter 	:: DataMissVal = -9999 

logical :: QRetain

integer :: NYear,NMonth,NStn,NRetain
integer :: XYear,XMonth,XStn,XRetain, Beg,End,AllocStat,Elastic

NYear=size(Adj,1) ; NMonth=size(Adj,2) ; NStn=size(Adj,3) 
Elastic=10 ; if (Differ) Elastic=5	! min len of any calib stretch

allocate (Retain(NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### PruneOrig: alloc #####"

do XStn=1,NStn
  Retain=.FALSE.
  
  do XYear=1,NYear
    XMonth=0 ; QRetain=.FALSE.
    do				! decide whether year needs retaining
          XMonth=XMonth+1
          if (	Orig (XYear,XMonth,XStn).NE.DataMissVal .AND. &
          	RefTS(XYear,XMonth,XStn).EQ.DataMissVal) QRetain=.TRUE.
          if (XMonth.EQ.NMonth.OR.QRetain) exit
    end do
    
    if (QRetain) then		! year needs retaining
      Beg=XYear-Elastic ; if (Beg.LT.1)     Beg=1
      End=XYear+Elastic ; if (End.GT.NYear) End=NYear

      do XRetain=Beg,End
        Retain(XRetain)=.TRUE.
      end do
    end if
  end do
  
  do XYear=1,NYear
    do XMonth=1,NMonth
      if (Orig (XYear,XMonth,XStn).NE.DataMissVal) then
        if (RefTS(XYear,XMonth,XStn).NE.DataMissVal) then
          Adj(XYear,XMonth,XStn) = Orig(XYear,XMonth,XStn)
        end if
        
        if (.NOT.Retain(XYear)) then
          Orig(XYear,XMonth,XStn) = DataMissVal
        end if
      end if
    end do
  end do
end do

deallocate (Retain, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### PruneOrig: dealloc #####"

end subroutine FilterOrig

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

end module GHCNdiscon