! ghcnrefiter.f90
! written by Tim Mitchell in September 2003
! module to construct reference time-series using an improved GHCN method
!	with iteration added to have as much data as possible with ref ts
! the unimproved method (ghcnref.f90) gives discons, and should not be used
! the improved method (ghcnrefplus.f90) risks leaving most without ref ts
! MAIN MODULES = MakeRefTS (top), Trustworthy (bottom)

module GHCNrefIter

use FileNames
use GridOps
use SortMod
use PlainPerm
use Regress

implicit none

contains

!*******************************************************************************
! a wrapper that calls most of the other subroutines in turn
! returns the reference ts
! for the stns for which Calc is true, the ref ts is reinit and calc
! Decay contains the correlation decay distance
! Differ: .TRUE. means (y+1) minus y, .FALSE. (y+1) divide by y
! if Disc is included, any discon lead to different per for same stn
! only those portions OK in Trust are allowed to enter the ref ts
! if Careful=.TRUE. quality is emphasised relative to coverage in calc ref ts
! if Omit is not -999, this number of years may be missing from the ref ts
!	for a particular stn; if -999, there is no constraint
! where GotRef=.TRUE. a ref ts (whether complete or not) has been calc
! NeedBeg/End (indices, not AD) denote a stretch that must be incl in any ref ts
! ForceBeg/End does the same, but varying the period by stn (-999=no compulsion)
! if Omit AND Need/Force are set, the default is to allow Omit missing values in
!	the needed/forced period; if this is undesirable set MakeAll
! if MakeMore AND Need/Force are set, no ref ts will be calc unless > than N/F

subroutine MakeRefTS (RefTS,Orig,Lat,Lon,Calc,Trust,Decay,Differ, &
		      Careful,Omit,GotRef, Disc,KPar,KNeigh,KLevel, &
		      NeedBeg,NeedEnd,ForceBeg,ForceEnd,MakeAll,MakeMore,Verbose)

logical,pointer,dimension(:,:),optional	:: Disc
logical,pointer,dimension(:,:)		:: Trust,StnPerTrust
logical,pointer,dimension(:)		:: GotRef,Calc,CandValid	! stn
logical,pointer,dimension(:)		:: PerTrust
real, pointer, dimension (:,:,:)	:: Par
real, pointer, dimension (:,:)		:: PoolWei
real, pointer, dimension (:)		:: Lat,Lon
integer, pointer, dimension (:,:,:) 	:: RefTS,Orig,FD,Abs,ParPer
integer, pointer, dimension (:,:) 	:: StnPer,Neigh,CandFD,CandAbs,Pool
integer, pointer, dimension (:) 	:: Per2Stn,Stn2Per,PerBeg,PerEnd
integer, pointer, dimension (:) 	:: InclBeg,InclEnd
integer, pointer ,dimension (:),optional:: ForceBeg,ForceEnd
integer, dimension (12) 		:: ParMonTot

logical, intent(in)		:: Differ,Careful
logical, intent(in), optional	:: MakeAll,MakeMore
real, intent(in)		:: Decay
integer, intent(in)		:: Omit
integer, intent(in), optional	:: KPar,KNeigh,KLevel,NeedBeg,NeedEnd,Verbose
integer, parameter 		:: DataMissVal = -9999
real, parameter 		:: MissVal = -999.0

logical :: InclAll,InclMore
real    :: OpValid,OpEn, AreThresh,Progress
integer :: NYear,NMonth,NStn,NPar,NNeigh,NPer,NLevel,NTen
integer :: XYear,XMonth,XStn,XPAr,XNeigh,XPer,XLevel,XTen
integer :: AllocStat, CandBeg,CandEnd, DumpTot,DumpEn, RefBeg,RefEnd
integer :: DataCheck,DataBlind,Stn0Per
character (len=10) :: CurrentTime
character (len=08) :: CurrentDate

NYear=size(Orig,1) ; NMonth=12 ; NStn=size(Orig,3) ; Progress=0.0
NPar=5 ; if (present(KPar)) NPar=KPar
NNeigh=100 ; if (present(KNeigh)) NNeigh=KNeigh
NLevel=100 ; if (present(KLevel)) NLevel=KLevel
InclAll =.FALSE. ; if (present(MakeAll))  InclAll =.TRUE.
InclMore=.FALSE. ; if (present(MakeMore)) InclMore=.TRUE.

if (present(ForceBeg).AND.present(ForceEnd)) then
  InclBeg => ForceBeg ; InclEnd => ForceEnd
else
  allocate (InclBeg(NStn),InclEnd(NStn), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Alloc: 0 #####"
  
  if (present(NeedBeg).AND.present(NeedEnd)) then
    InclBeg=NeedBeg ; InclEnd=NeedEnd
  else
    InclBeg=MissVal ; InclEnd=MissVal
  end if
end if

if (Careful) then
  AreThresh=0.4
else
  AreThresh=0.2		! could be 0.2 is one is daring
end if

allocate (StnPer(NStn,NYear),StnPerTrust(NStn,NYear), &
	  Neigh(NStn,NNeigh), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Alloc: 1 #####"
if (present(Verbose)) write (99,"(a,i10)"), "MakeRefTS is operating on stns: ", count(Calc)

		! break up original stns into consistent per, with trusted data
		! one stn may contain 0-(NOYear/2) pers, spec in StnPer
		! total of pers from all stns = NPer
if (present(Disc)) then
  call FindPer (Orig,StnPer,StnPerTrust,NPer,Lat,Lon,Trust,Differ,Decay=Decay,Calc=Calc,Disc=Disc)
else
  call FindPer (Orig,StnPer,StnPerTrust,NPer,Lat,Lon,Trust,Differ,Decay=Decay,Calc=Calc)
end if
if (present(Verbose)) write (99,"(a,i10)"), "FindPer has identified pers: ", NPer

allocate (FD(NYear,NMonth,NPer), Per2Stn(NPer), Stn2Per(NStn), PerTrust(NPer), &
	  Abs(NYear,NMonth,NPer), PerBeg(NPer), PerEnd(NPer), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Alloc: 2 #####", NPer

		! generates set of absolute series for all pers (Abs)
		! also records beg,end of each per (PerBeg,PerEnd)
		! also records the stn pertaining to each period (Per2Stn)
		! also records the FIRST per pertaining to each stn (Stn2Per)
		! records the pers that are trusted
call InitialAbs (Orig,StnPer,StnPerTrust,Abs,Per2Stn,Stn2Per, &
		 PerTrust,PerBeg,PerEnd)
deallocate (StnPer,StnPerTrust, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: delloc: 2 #####"
if (present(Verbose)) write (99,"(a,i10)"), "InitialAbs has trusted pers: ", count(PerTrust)

		! fills Neigh with the indices of the nearest stns, in order
		! only stns within corr decay dist are included
		! the upper limit to the no. stns is NNeigh (use KNeigh in call)
		! only stns with trusted periods are allowed as neighbours
call OrderNeigh (Neigh,Stn2Per,Per2Stn,PerTrust,PerBeg,PerEnd,Lat,Lon,Decay,Differ)
if (present(Verbose)) write (99,"(a)"), "OrderNeigh has ID neighbours"

		! uses neighbours to infill gaps in Abs
		! any remaining gaps are set to the local mean
call FillGaps (Abs,Per2Stn,Stn2Per,PerBeg,PerEnd,Neigh,Differ)
if (present(Verbose)) write (99,"(a)"), "FillGaps has filled any gaps"

		! convert Abs into genuine first-difference series
		! PerBeg is changed accordingly (+1)
		! Differ=.TRUE. means (y+1) minus y, .FALSE. (y+1) divide by y
call ConvertFD (Abs,FD,PerBeg,PerEnd,Differ)
if (present(Verbose)) write (99,"(a)"), "ConvertFD has created FD series"

allocate (CandFD(NYear,NMonth), CandAbs(NYear,NMonth), CandValid(NYear), &
	  Pool(NMonth,NLevel),PoolWei(NMonth,NLevel), &
	  ParPer(NMonth,NPar,NLevel),Par(NYear,NMonth,NPar), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Alloc: 3 #####"

do XStn=1,NStn
 if (present(Verbose)) then
   if (mod(XStn,10).EQ.0.0) then
     write (99,"(2(a,i5),a,i3,a)"), "----- STATION ", XStn, " OF ", NStn, " (", &
   		nint(100.0*(real(XStn)/real(NStn))), "%) -----"
   else if ((real(XStn)/real(NStn)).GT.Progress) then
     call date_and_time (CurrentDate, CurrentTime)
     if (Progress.EQ.0.0) print "(a)", "   >  % hh:mm"
     print "(a,i2,x,a2,a1,a2)", "   > ", nint(Progress*100), &
     			CurrentTime(1:2), ":", CurrentTime(3:4)
     Progress=Progress+0.1
   end if
 end if
 
 if (Calc(XStn).AND.Stn2Per(XStn).NE.MissVal) then

  		! get the candidate series, possibly including missing values
  if (InclBeg(XStn).NE.MissVal.AND.InclEnd(XStn).NE.MissVal) then
    call GetCand  (Abs,FD,CandAbs,CandFD,Stn2Per,Per2Stn,PerBeg,PerEnd, &
  	XStn,CandBeg,CandEnd,CandValid,NeedBeg=InclBeg(XStn),NeedEnd=InclEnd(XStn))
  else
    call GetCand  (Abs,FD,CandAbs,CandFD,Stn2Per,Per2Stn,PerBeg,PerEnd, &
  	XStn,CandBeg,CandEnd,CandValid)
  end if
  
  		! build the pool of neighbouring pers on which to draw, with weights
  		! the max size of the pool is governed by NLevel (use KLevel in call)
  call BuildPool  (FD,CandFD,Stn2Per,Per2Stn,PerTrust,PerBeg,PerEnd, &
  			XStn,CandBeg,CandEnd,Neigh,Pool,PoolWei,AreThresh,Differ)
		
		! select the pers to form into each parallel FD series (ParPer)
		! the number of parallels is governed by NPar (use KPar in call)
		! each parallel must be filled completely; omissions shorten RefTS
		! each per can only contribute to one parallel
		! RefBeg,End define the ref ts: applies to all calendar months
  if (InclBeg(XStn).NE.MissVal.AND.InclEnd(XStn).NE.MissVal) then
    call FormPar (Pool,PoolWei,PerBeg,PerEnd,RefBeg,RefEnd,ParPer,CandBeg,CandEnd,&
    			CandValid,Careful,Omit,Differ,InclAll,InclMore, &
			NeedBeg=InclBeg(XStn),NeedEnd=InclEnd(XStn))
  else
    call FormPar (Pool,PoolWei,PerBeg,PerEnd,RefBeg,RefEnd,ParPer,CandBeg,CandEnd,&
    			CandValid,Careful,Omit,Differ,InclAll,InclMore)
  end if
  
  		! make each parallel ts (Par), as specified in ParPer
  call MakePar    (Abs,CandAbs,PerBeg,PerEnd,RefBeg,RefEnd,ParPer,Par,Per2Stn,Differ)
  
  		! combine the parallels into a ref ts
  call CombinePar (CandAbs,Par,ParPer,RefBeg,RefEnd,RefTS,XStn,Differ)  
 end if
end do

deallocate (FD,Abs,Per2Stn,Stn2Per,PerBeg,PerEnd,Neigh,PerTrust, &
	    CandFD,CandAbs,Pool,PoolWei,ParPer,Par,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Dealloc: 4 #####"

if (present(ForceBeg).AND.present(ForceEnd)) then
  	! nothing required
else
  deallocate (InclBeg,InclEnd,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### MakeRefTS: Dealloc: 5 #####"
end if

GotRef=.FALSE.
do XStn=1,NStn			! fill GotRef
  if (Calc(XStn)) then
    XMonth=0
    do				! quicker to iterate month,year
      XMonth=XMonth+1 ; XYear=0
      do
        XYear=XYear+1
        if (	Orig(XYear,XMonth,XStn).NE.DataMissVal.AND. &
        	RefTS(XYear,XMonth,XStn).NE.DataMissVal) &
        	GotRef(XStn)=.TRUE.
        if (XYear.EQ.NYear.OR.GotRef(XStn)) exit
      end do
      if (XMonth.EQ.NMonth.OR.GotRef(XStn)) exit
    end do
  end if
end do

end subroutine MakeRefTS

!*******************************************************************************
! uses the ts info in Par to construct a ref ts for that stn (which must be 
!	identified through CandStn)
! ParPer is only used to identify the parallels in Par that are actually filled

subroutine CombinePar (CandAbs,Par,ParPer,RefBeg,RefEnd,RefTS,CandStn,Differ)

real, pointer, dimension (:,:,:)	:: Par
real, pointer, dimension (:)		:: Exe,Wye,Weight,Comb

integer, pointer, dimension (:,:,:)	:: ParPer
integer, pointer, dimension (:,:,:)	:: RefTS
integer, pointer, dimension (:,:)	:: CandAbs

logical, intent(in)	:: Differ
integer, intent(in)	:: CandStn,RefBeg,RefEnd

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

logical :: Okay
real 	:: Aye,Bee,Enn,ParTot,ParTotSq,CandTot,CandTotSq,CandMean,ParMean,RegEnn
real 	:: AdjMean,AdjStDev,Numer,Denom,TotWeight,Actual
integer :: NYear,NMonth,NPar, AllocStat,Elastic
integer :: XYear,XMonth,XPar

Elastic=10 ; if (Differ) Elastic=5
NYear=size(Par,1) ; NMonth=12 ; NPar=size(Par,3) 

allocate (Exe(NYear),Wye(NYear),Comb(NYear),Weight(NPar), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombinePar: alloc #####"

do XMonth = 1, NMonth
 do XYear=1,NYear
   RefTS(XYear,XMonth,CandStn)=DataMissVal
 end do
 
 if (RefBeg.NE.MissVal.AND.RefEnd.NE.MissVal) then
   Exe=MissVal ; TotWeight=0 ; Weight=0 ; Comb=MissVal
   
   do XYear=RefBeg,RefEnd
     if (CandAbs(XYear,XMonth).NE.DataMissVal) Exe(XYear)=real(CandAbs(XYear,XMonth))
   end do
  
   do XPar=1,NPar				! prepare each parallel for merge
    if (ParPer(XMonth,XPar,1).NE.MissVal) then
     Wye=MissVal
     do XYear=RefBeg,RefEnd
       Wye(XYear)=Par(XYear,XMonth,XPar)
     end do
     
     if (Differ) then
       call MergeForDiff  (Exe,Wye,RefBeg,RefEnd)	! adjust par to match cand
     else
       call MergeForRatio (Exe,Wye,RefBeg,RefEnd)	
     end if
     
     do XYear=RefBeg,RefEnd
         Par(XYear,XMonth,XPar)=Wye(XYear)
     end do
     call RegressVectors (Exe,Wye,Aye,Bee,Weight(XPar),RegEnn,Elastic,&
     				Beg=RefBeg,End=RefEnd)	! calc weight of par
     if (Weight(XPar).GT.0) then
         Weight(XPar)=Weight(XPar)**2
     else
         Weight(XPar)=0.0001
     end if
     
     TotWeight=TotWeight+Weight(XPar)
    end if
   end do
   
   do XYear=RefBeg,RefEnd			! calc weighted mean -> Comb
    ParTot=0
    do XPar=1,NPar
     if (ParPer(XMonth,XPar,1).NE.MissVal) &
      		ParTot=ParTot+(Par(XYear,XMonth,XPar)*Weight(XPar))
    end do
    Comb(XYear)=ParTot/TotWeight
   end do
  
   if (Differ) then				! adjust comb to match cand
     call MergeForDiff  (Exe,Comb,RefBeg,RefEnd)	
   else
     call MergeForRatio (Exe,Comb,RefBeg,RefEnd)
   end if
   
   do XYear=RefBeg,RefEnd	
     RefTS(XYear,XMonth,CandStn)=nint(Comb(XYear))
   end do
 end if
end do

deallocate (Exe,Wye,Weight,Comb, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombinePar: dealloc #####"

end subroutine CombinePar

!*******************************************************************************
! uses ParPer to construct the parallel Abs series, returned in Par

subroutine MakePar (Abs,CandAbs,PerBeg,PerEnd,RefBeg,RefEnd,ParPer,Par,Per2Stn,Differ)

real, pointer, dimension (:)		:: Stand,Addit
real, pointer, dimension (:,:,:)	:: Par

integer, pointer, dimension (:,:,:)	:: ParPer,Abs
integer, pointer, dimension (:,:)	:: CandAbs
integer, pointer, dimension (:)		:: PerBeg,PerEnd
integer, pointer, dimension (:) 	:: Per2Stn		! for info only

logical, intent(in)			:: Differ
integer, intent(in)			:: RefBeg,RefEnd

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

logical :: Okay
real 	:: ParTot,ParTotSq,AddTot,AddTotSq,AddVal,MergeEn
real 	:: AddMean,ParMean,AdjMean,AdjStDev,Numer,Denom
integer :: NYear,NMonth,NPar,NLevel,NPer
integer :: XYear,XMonth,XPar,XLevel,XPer
integer :: Joint0,Joint1,ParBeg,ParEnd, AllocStat,Elastic

Elastic=10 ; if (Differ) Elastic=5
NYear=size(Abs,1) ; NMonth=12 ; NPer=size(Abs,3) 
NPar=size(ParPer,2) ; NLevel=size(ParPer,3) ; Par=MissVal

allocate (Stand(NYear),Addit(NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakePar: alloc #####"

do XMonth = 1, NMonth
  do XPar = 1, NPar
    XPer=ParPer(XMonth,XPar,1) ; Stand=MissVal
    if (XPer.NE.MissVal) then				! initialise parallel
      ParBeg=PerBeg(XPer) ; ParEnd=PerEnd(XPer)
      do XYear = PerBeg(XPer),PerEnd(XPer)
        Stand(XYear)=Abs(XYear,XMonth,XPer)
      end do
    end if
    
    XLevel=1
    do							! find mergers
      XLevel=XLevel+1
      
      if (ParPer(XMonth,XPar,XLevel).NE.MissVal) then
        XPer=ParPer(XMonth,XPar,XLevel)
        Joint0=PerBeg(XPer) ; if (Joint0.LT.ParBeg) Joint0=ParBeg
        Joint1=PerEnd(XPer) ; if (Joint1.GT.ParEnd) Joint1=ParEnd
        MergeEn=Joint1-Joint0+1
        if (Differ.AND.MergeEn.GT.10) then	! for diffs (not ratios), limit to 10
          Joint0=Joint1-9 ; MergeEn=10		!   to restrict incorp undetected inhomog
        end if
        
        if (MergeEn.GE.Elastic) then
          Addit=MissVal
	  do XYear=PerBeg(XPer),PerEnd(XPer)		! initialise per
	    Addit(XYear)=real(Abs(XYear,XMonth,XPer))
	  end do
	  
	  if (Differ) then
	    call MergeForDiff  (Stand,Addit,Joint0,Joint1)	! adjust per
	  else
	    call MergeForRatio (Stand,Addit,Joint0,Joint1)	! adjust per
	  end if
	  
          do XYear = (Joint1+1),PerEnd(XPer)		! append adj data
              Stand(XYear) = Addit(XYear)
          end do
          ParEnd=PerEnd(XPer)				! update par end
        end if
      end if

      if (XLevel.GT.NLevel) exit
      if (ParPer(XMonth,XPar,XLevel).EQ.MissVal) exit
    end do
    
    do XYear=1,NYear
      Par(XYear,XMonth,XPar)=Stand(XYear)		! fill with parallel
      if (XYear.LT.RefBeg.OR.XYear.GT.RefEnd) Par(XYear,XMonth,XPar)=MissVal
    end do
  end do
end do

deallocate (Stand,Addit, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakePar: dealloc #####"

end subroutine MakePar

!*******************************************************************************
! adjusts the Addit(ional) vector to match the characteristics of the corresponding 
!	Stand(ard), for ratio-based (not difference-based) data, on the
!	assumption that both are gamma-distributed
! the old (x) becomes the new (y) through y=ax**b
! b is calc iteratively, so that the shape parameters match
! then a is deduced, such that the scale parameters match

subroutine MergeForRatio (Stand,Addit,Overlap0,Overlap1)

real, pointer, dimension (:)	:: Stand,Addit
integer, intent(in)		:: Overlap0,Overlap1
real, parameter :: MissVal = -999.0

real 	:: ParTot,ParTotSq,AddTot,AddTotSq,AddMean,En,Add,Power,Multi,Iter,Cheat,New
real 	:: StandShape,StandScale,StandMean,StandVar
real	:: LateShape,LateScale,LateMean,LateVar
real	:: LastShape,LastScale,LastMean,LastVar
integer :: XYear,NYear

NYear=size(Addit) ; Power=MissVal ; Multi=MissVal ; Cheat=0.0
AddTot=0 ; AddTotSq=0 ; ParTot=0 ; ParTotSq=0 ; En=0
do XYear = Overlap0,Overlap1			
  if (Stand(XYear).NE.MissVal.AND.Addit(XYear).NE.MissVal) then
            AddTot	= AddTot	+ Addit(XYear)
            AddTotSq	= AddTotSq	+(Addit(XYear)**2)
            ParTot	= ParTot	+ Stand(XYear)
            ParTotSq	= ParTotSq	+(Stand(XYear)**2)
	    En=En+1
  end if
end do

if      (En.EQ.0) then
    Power=1.0 ; Multi=1.0
else if (En.EQ.1) then
    if (AddTot.GT.0) then
      Power=1.0 ; Multi=ParTot/AddTot
    else
      Power=1.0 ; Multi=1.0 ; Cheat=ParTot-AddTot
    end if
else if (ParTot.EQ.0) then
    Power=1.0 ; Multi=0.0
else if (AddTot.EQ.0) then
    Power=1.0 ; Multi=1.0 ; Cheat=(ParTot/En)
else
    StandMean=ParTot/En 	; StandVar=(ParTotSq/En)-(StandMean**2)
    LateMean=AddTot/En 		; LateVar=(AddTotSq/En)-(LateMean**2)
    
    if (StandVar.EQ.0.OR.LateVar.EQ.0) then
      Power=1.0 ; Multi=StandMean/LateMean
    end if
end if

if (Power.Eq.MissVal) then
  LastShape=StandShape+((LateShape-StandShape)*2)
  Power=1.0 ; Iter=0
  Add=0.4 ; if (LateShape.LT.StandShape) Add=0.0-0.2

  do				! seek power giving match to shape parameter 
    if (mod(nint(100.0*(LateShape-StandShape)),1).LT.1) exit
    if (Iter.GE.10) exit
    
    if      (LateShape.EQ.LastShape) then
      Power=Power-Add			! go back to last power
      Add=Add*0.5			! no shape, so try again with half add
    else if (((LateShape-StandShape)*(LastShape-StandShape)).LE.0) then	
      Add=Add*(0.0-0.5) 		! overshoot
    else if (abs(LateShape-StandShape).GE.abs(LastShape-StandShape)) then
      Add=Add*(0.0-2.0)			! wrong direction
    end if				! else undershoot
    
    Iter=Iter+1
    LastShape=LateShape
    Power=Power+Add
    
    AddTot=0 ; AddTotSq=0
    do XYear = Overlap0,Overlap1			
      if (Stand(XYear).NE.MissVal.AND.Addit(XYear).NE.MissVal) then
            AddTot   = AddTot   + (Addit(XYear)** Power)
            AddTotSq = AddTotSq + (Addit(XYear)**(Power*2.0))
      end if
    end do
    LateVar=(AddTotSq/En)-((AddTot/En)**2)
    
    if (LateVar.GT.0) LateShape=((AddTot/En)**2)/LateVar
  end do
  
  Multi=StandMean/(AddTot/En)
end if
  
do XYear = 1,NYear			! adjust addit data
  if (Addit(XYear).NE.MissVal) then
    New=Cheat+(Multi*(Addit(XYear)**Power))
!    write (99,"(i4,3f10.2)"), XYear,Stand(XYear),Addit(XYear),New ! #################
    Addit(XYear)=New
  end if
end do

end subroutine MergeForRatio

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

subroutine MergeForDiff (Stand,Addit,Overlap0,Overlap1)

real, pointer, dimension (:)	:: Stand,Addit
integer, intent(in)		:: Overlap0,Overlap1
real, parameter :: MissVal = -999.0
real 	:: ParTot,ParTotSq,AddTot,AddTotSq,AddMean,ParMean,Numer,Denom,AdjStDev,En
integer :: XYear,NYear

NYear=size(Addit)
!if (NYear.NE.size(Stand)) print "(a,2i6)", &
!	"   > @@@@@ ERROR: MergeForDiff: Stand v Addit", size(Stand),NYear
!if (Overlap0.LT.1.OR.Overlap1.GT.NYear) print "(a,3i6)", &
!	"   > @@@@@ ERROR: MergeForDiff: beg,end,NYear", Overlap0,Overlap1,NYear

ParTot=0 ; ParTotSq=0 ; AddTot=0 ; AddTotSq=0 ; En=0
do XYear = Overlap0,Overlap1			
  if (Stand(XYear).NE.MissVal.AND.Addit(XYear).NE.MissVal) then
            ParTot	= ParTot	+ Stand(XYear)
            ParTotSq	= ParTotSq	+(Stand(XYear)**2)
            AddTot	= AddTot	+ Addit(XYear)
            AddTotSq	= AddTotSq	+(Addit(XYear)**2)
	    En=En+1
  end if
end do

if (En.GT.0) then          
  AddMean = AddTot/En			! normal for addit data
  ParMean = ParTot/En			! normal for existing data
  Numer = ParTotSq-(ParTot*ParMean)
  Denom = AddTotSq-(AddTot*AddMean)
  AdjStDev = 1 ; if (Numer.GT.0.AND.Denom.GT.0) AdjStDev = sqrt(Numer/Denom)
	  
  do XYear = 1,NYear			! adjust addit data
    if (Addit(XYear).NE.MissVal) Addit(XYear)=ParMean+((Addit(XYear)-AddMean)*AdjStDev)
  end do
end if

end subroutine MergeForDiff

!*******************************************************************************
! this routine decides which per to use to build each parallel
! returns ParPer (Nmonth,NPar,NLevel): for each month and parallel, the index of 
!	each per to use in order; unused slots are -999
! also returns RefBeg,End - year indices for the ref ts: apply to all months
! NeedBeg,NeedEnd may be used to force a refts to include 1961-1990
! if Need is set, AllNeeded clarifies whether all the needed period
!	is actually required, or whether Omit missing values may be permitted
! AllNeeded may be overridden by MoreNeeded, which ensures that the needed period
!	is not only included, but exceeded; if so, missing values are not allowed

subroutine FormPar (Pool,PoolWei,PerBeg,PerEnd,RefBeg,RefEnd,ParPer,CandBeg, &
		    CandEnd,CandValid,Careful,Omit,Differ,AllNeeded,MoreNeeded, &
		    NeedBeg,NeedEnd)

logical, pointer, dimension (:)		:: CandValid
real, pointer, dimension (:,:)		:: PoolWei
integer, pointer, dimension (:,:,:)	:: ParPer	
integer, pointer, dimension (:,:)	:: Pool
integer, pointer, dimension (:)		:: PerBeg,PerEnd
logical, intent(in)			:: Careful,Differ,AllNeeded,MoreNeeded
integer, intent(out)			:: RefBeg,RefEnd
integer, intent(in)			:: CandBeg,CandEnd,Omit	
integer, intent(in), optional		:: NeedBeg,NeedEnd 	! stretch must be incl

real, parameter :: MissVal = -999.0

logical :: Abandon,Complete
real    :: TestQuality,BestQuality,Quality,IdealQuality,Weight
integer :: NYear,NMonth,NLevel,NPar
integer :: XYear,XMonth,XLevel,XPar
integer :: TestBeg,TestEnd,BestBeg,BestEnd,Latest,ValidTot,Elastic

Elastic=10 ; if (Differ) Elastic=5
NYear=size(CandValid,1)
NMonth=size(Pool,1) ; NLevel=size(Pool,2) ; NPar=size(ParPer,2)
if (size(ParPer,1).NE.NMonth.OR.size(ParPer,3).NE.NLevel) &
		print*, "  > ##### FormPar: bad-len: ParPer #####"

Weight=1.0 ; if (Careful) Weight=sqrt(real(NPar))
Complete=.FALSE. ; RefBeg=MissVal ; RefEnd=MissVal
TestQuality=0 ; TestBeg=CandBeg ; TestEnd=CandEnd 
BestQuality=0 ; BestBeg=MissVal ; BestEnd=MissVal
IdealQuality = real(Weight*(TestEnd-TestBeg+1))	

do
  Abandon=.FALSE.
  
  do
    XMonth=0 ; TestQuality=0 ; ParPer=MissVal
    
    do
      XMonth=XMonth+1
      call EnginePar (Pool,PoolWei,PerBeg,PerEnd,ParPer, &
		      XMonth,TestBeg,TestEnd,Quality,Latest,Abandon,Careful,Differ)
      if (Quality.NE.MissVal) then
        if (.NOT.Careful) Quality=TestEnd-TestBeg+1	! no weight from correl
        TestQuality=TestQuality+Quality
      end if
      if (Abandon.OR.XMonth.EQ.12) exit
    end do
    
    if (Abandon.AND.Latest.NE.MissVal) then	! could try with earlier end
      TestEnd=Latest		! if permitted below, we reiterate
      
      if (Omit.EQ.MissVal.OR.((TestBeg-CandBeg)+(CandEnd-TestEnd)).LE.Omit) then
        XYear=TestBeg-1 ; ValidTot=0
        do			! ensure at least 5/10 yrs of original data
          XYear=XYear+1
          if (CandValid(XYear)) ValidTot=ValidTot+1
          if (XYear.EQ.TestEnd.OR.ValidTot.GE.Elastic) exit
        end do
	
        if (ValidTot.LT.Elastic) then	! insufficient valid data
          Latest=MissVal
        else if (present(NeedEnd)) then	! a certain period required
	  if (Omit.EQ.MissVal) then		! no omissions constraint
	    if (TestEnd.LT.NeedEnd) then		! ends too early
	      Latest=MissVal
	    else if (MoreNeeded) then			! must have more
	      if (TestBeg.GE.NeedBeg.OR.TestEnd.LE.NeedEnd) Latest=MissVal
	    end if
	  else					! omissions constraint
	    if (MoreNeeded) then			! must have more
	      if (TestBeg.GE.NeedBeg.OR.TestEnd.LE.NeedEnd) Latest=MissVal
	    else if (AllNeeded) then			! must have at least
	      if (TestEnd.LT.NeedEnd) Latest=MissVal
	    else					! must have, minus Omit
	      if ((TestEnd+Omit).LT.NeedEnd) Latest=MissVal
	    end if
	  end if
        end if
      else
        Latest=MissVal		! not permitted by Omit
      end if
    end if
    
    if (Abandon.AND.Latest.EQ.MissVal) exit
    if (.NOT.Abandon) exit
  end do
  
  if (.NOT.Abandon) then				! find best score
    if ((TestQuality/12.0).GT.BestQuality) then
      BestQuality=(TestQuality/12.0) ; BestBeg=TestBeg ; BestEnd=TestEnd
    end if
  end if
  
  TestBeg=TestBeg+5 ; TestEnd=CandEnd
  IdealQuality = real(Weight*(TestEnd-TestBeg+1))	! still possible
  
  if (IdealQuality.LT.BestQuality.OR.TestBeg.GT.TestEnd) then			
    RefBeg=BestBeg ; RefEnd=BestEnd ; Complete=.TRUE.	! ideal < achieved
  else if (Omit.EQ.MissVal.OR.((TestBeg-CandBeg)+(CandEnd-TestEnd)).LE.Omit) then
    XYear=TestBeg-1 ; ValidTot=0			! ideal achievable
    do							! 5 yrs of orig?
          XYear=XYear+1
          if (CandValid(XYear)) ValidTot=ValidTot+1
          if (XYear.EQ.TestEnd.OR.ValidTot.GE.Elastic) exit
    end do
    if (ValidTot.LT.Elastic) then			! section useless
      RefBeg=BestBeg ; RefEnd=BestEnd ; Complete=.TRUE.
    else if (present(NeedBeg)) then
      if (Omit.EQ.MissVal) then			! no omissions constraint
	    if (TestBeg.GT.NeedBeg) then		! starts too late
	      Complete=.TRUE.
	    else if (MoreNeeded) then			! must have more
	      if (TestBeg.GE.NeedBeg.OR.TestEnd.LE.NeedEnd) Complete=.TRUE.
	    end if
      else					! omissions constraint
	    if (MoreNeeded) then			! must have more
	      if (TestBeg.GE.NeedBeg.OR.TestEnd.LE.NeedEnd) Complete=.TRUE.
	    else if (AllNeeded) then			! must have at least
	      if (TestBeg.GT.NeedBeg) Complete=.TRUE.
	    else					! must have, minus Omit
	      if ((TestBeg-Omit).GT.NeedBeg) Complete=.TRUE.
	    end if
      end if
      
      if (Complete) then
        RefBeg=BestBeg ; RefEnd=BestEnd
      end if
    end if
  else
    RefBeg=BestBeg ; RefEnd=BestEnd ; Complete=.TRUE.	! must end (Omit)
  end if
  
  if (Complete) exit
end do

ParPer=MissVal

if (RefBeg.NE.MissVal) then				! restore the best score
  do XMonth=1,NMonth
    call EnginePar (Pool,PoolWei,PerBeg,PerEnd,ParPer, &
		      XMonth,RefBeg,RefEnd,Quality,Latest,Abandon,Careful,Differ)
    if (Abandon) print*, "  > ##### FormPar: unexpectedly failed parallels #####"
  end do
end if

end subroutine FormPar

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

subroutine EnginePar (Pool,PoolWei,PerBeg,PerEnd,ParPer, &
		      QMonth,Beg,End,Quality,Latest,Abandon,Careful,Differ)

real, pointer, dimension (:,:)		:: PoolWei
real, pointer, dimension (:)		:: Score,PoolUse,ParWei
integer, pointer, dimension (:,:,:)	:: ParPer
integer, pointer, dimension (:,:)	:: Pool
integer, pointer, dimension (:)		:: PerBeg,PerEnd,Order,LooseEnd,ParBeg,ParEnd

logical, intent(in)	:: Differ,Careful	! TRUE:weight by sqrt(NPar)
logical, intent(out)	:: Abandon
real, intent(out)	:: Quality
integer, intent(out)	:: Latest
integer, intent(in) 	:: QMonth,Beg,End

real, parameter :: MissVal = -999.0

real    :: ScoreTot,BestTot,TargetTot,AttemptTot,Weight
integer :: BestBeg,BestEnd,CandLen,LastEnd,AllPar,Elastic
integer :: XMonth,XLevel,XScore,XPar,XParInit,XPer,XLooseEnd, XYear
integer :: NMonth,NLevel,NScore,NPar,NParInit,NPer,NLooseEnd, AllocStat

Elastic=10 ; if (Differ) Elastic=5
NMonth=size(Pool,1) ; NLevel=size(Pool,2) ; NPar=size(ParPer,2)
if (size(ParPer,1).NE.NMonth.OR.size(ParPer,3).NE.NLevel) &
		print*, "  > ##### EnginePar: bad-len: ParPer #####"
allocate (Score(NLevel),LooseEnd(NPar),PoolUse(NLevel), &
	  ParBeg(NPar),ParEnd(NPar),ParWei(NPar), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: alloc #####"

Abandon=.FALSE. ; Quality=MissVal ; Latest=MissVal ; XMonth=QMonth ; AllPar=NPar

do
  PoolUse=MissVal ; ParWei=MissVal ; ParBeg=MissVal ; ParEnd=MissVal ; LooseEnd=MissVal
  NLooseEnd=0 ; Score=MissVal ; NScore=0
  
  do XPar=1,AllPar			! this is vital when NPar=NPar-1, reiterate
    do XLevel=1,NLevel
      ParPer(XMonth,XPar,XLevel) = MissVal
    end do
  end do
  
  do XLevel=1,NLevel			! find coeff for periods starting at right time
    XPer=Pool(XMonth,XLevel)
    if (XPer.NE.MissVal) then
      if (PerBeg(XPer).LE.Beg.AND.PerEnd(XPer).GE.(Beg+Elastic-1)) then
        Score(XLevel) = 0-PoolWei(XMonth,XLevel)	! just makes it easier
        NScore=NScore+1					! to treat sorted results
      end if 
    end if						
  end do
  
  if      (NScore.LT.2) then		! insufficient pers for start year
    Abandon=.TRUE.
  else if (NScore.LT.NPar) then		! adjust NPar to match available NPer
    NPar=NScore
  end if
  
  if (.NOT.Abandon) then
    allocate (Order(NLevel),stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: alloc Order 1 #####"

    call QuickSort (Reals=Score,OrderValid=Order)
    XScore=0
    do					! initialise the parallels
      XScore=XScore+1 ; XPar=XScore ; XLevel=Order(XScore) ; XPer=Pool(XMonth,XLevel)
      PoolUse(XLevel) = real(XPar)		! shows avail of per
      ParPer(XMonth,XPar,1) = XPer		! sequence per into par
      ParBeg(XPar) = PerBeg(XPer)		! parallel beg/end
      ParEnd(XPar) = PerEnd(XPer)
      ParWei(XPar) = PoolWei(XMonth,XLevel) * (PerEnd(XPer)-Beg+1)
      if (PerEnd(XPer).LT.End) then
        LooseEnd(XPar)=PerEnd(XPer) ; NLooseEnd=NLooseEnd+1
      end if
      if (XScore.EQ.NScore.OR.XPar.EQ.NPar) exit
    end do
    
    deallocate (Order,stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: dealloc Order 1 #####"
    
    do					! improve the parallels
      if (NLooseEnd.EQ.0) exit
      
      allocate (Order(NPar),stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: alloc Order 2 #####"
      call QuickSort (Ints=LooseEnd,OrderValid=Order)	! find next parallel
      XPar=Order(1)
      deallocate (Order,stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: dealloc Order 2 #####"
    
      Score=MissVal ; NScore=0
      do XLevel=1,NLevel				! calc score for all good per
       if (PoolUse(XLevel).EQ.MissVal.AND. &		! per is unused and here
       		Pool(XMonth,XLevel).NE.MissVal) then 
        XPer=Pool(XMonth,XLevel)
        if (PerBeg(XPer).LE.(ParEnd(XPar)-Elastic+1).AND. &
            PerEnd(XPer).GT.ParEnd(XPar)) then		! per is relevant (r assumed)
          NScore=NScore+1
          if (PerEnd(XPer).LT.End) then
            Score(XLevel) = 0 - (PoolWei(XMonth,XLevel) * (PerEnd(XPer)-ParEnd(XPar)))
          else
            Score(XLevel) = 0 - (PoolWei(XMonth,XLevel) * (End-ParEnd(XPar)))
          end if
        end if
       end if
      end do
    
      if (NScore.GT.0) then			! there are relevant pers
       allocate (Order(NLevel),stat=AllocStat)
       if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: alloc Order 3 #####"
       call QuickSort (Reals=Score,OrderValid=Order)	! find best per to add
       XLevel=Order(1) ; XPer=Pool(XMonth,XLevel)
       deallocate (Order,stat=AllocStat)
       if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: dealloc Order 3 #####"
      
       ParWei(XPar) = ParWei(XPar) + (PoolWei(XMonth,XLevel)*(PerEnd(XPer)-ParEnd(XPar)))
       ParEnd(XPar) = PerEnd(XPer)		! no change to ParBeg

       XScore=0
       do
        XScore=XScore+1
        if (ParPer(XMonth,XPar,XScore).EQ.MissVal) exit
       end do
       ParPer(XMonth,XPar,XScore) = XPer	! sequence per into par
      
       PoolUse(XLevel) = real(XPar)+(real(XScore-1)*0.01)

       if (PerEnd(XPer).GE.End) then 		! parallel complete
        LooseEnd(XPar)=MissVal ; NLooseEnd=NLooseEnd-1
       else if (XScore.EQ.NLevel) then		! no more room
        NLooseEnd=0
       else					! incomplete and more room			
        LooseEnd(XPar)=PerEnd(XPer)
       end if
      end if
      
      if (NScore.EQ.0) exit
    end do
    
    if (NLooseEnd.GT.0) then	! cannot complete all parallels
      allocate (Order(NPar),stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: alloc Order 4 #####"
      call QuickSort (Ints=LooseEnd,OrderValid=Order)	! find first-finishing parallel
      XPar=Order(1)
      if (Latest.LT.LooseEnd(XPar)) Latest=LooseEnd(XPar)
      deallocate (Order,stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: EnginePar: dealloc Order 4 #####"

      if (NPar.GT.2) then
        NPar=NPar-1
      else
        Abandon=.TRUE.
      end if
    else
      Weight=NPar ; if (Careful) Weight=sqrt(real(NPar))
      Quality=0
      do XPar = 1, NPar
        Quality=Quality+ParWei(XPar)
      end do
      Quality=Quality/Weight
    end if
  end if
  
  if (Abandon.OR.Quality.NE.MissVal) exit
end do

end subroutine EnginePar

!*******************************************************************************
! returns the Pool of neighbours from which parallels may be built
! also the weight to attach to each neighbour (PoolWei)

subroutine BuildPool (FD,CandFD,Stn2Per,Per2Stn,PerTrust,PerBeg,PerEnd, &
			CandStn,CandBeg,CandEnd,Neigh,Pool,PoolWei,AreThresh,Differ)

logical, pointer, dimension (:)		:: PerTrust
real, pointer, dimension (:,:)		:: PoolWei
real, pointer, dimension (:)		:: Exe,Wye
integer, pointer, dimension (:,:,:)	:: FD
integer, pointer, dimension (:,:)	:: CandFD,Neigh,Pool
integer, pointer, dimension (:)		:: Stn2Per,Per2Stn,PerBeg,PerEnd
integer, dimension (12)			:: Score,Level

logical, intent(in)	:: Differ
integer, intent(in)	:: CandStn,CandBeg,CandEnd
real, intent(in) 	:: AreThresh
real, parameter 	:: MissVal = -999.0
integer, parameter 	:: DataMissVal = -9999
integer, parameter 	:: ScoreInit=10

real    :: Aye,Bee,Enn,Are,Pea
integer :: NStn,NPer,NYear,NMonth,NHalt,NNeigh,NLevel, AllocStat
integer :: XStn,XPer,XYear,XMonth,XHalt,XNeigh,XLevel, Joint0,Joint1,Elastic

Elastic=10 ; if (Differ) Elastic=5
NYear=size(FD,1) ; NMonth=12 ; NPer=size(FD,3) 
NNeigh=size(Neigh,2) ; NLevel=size(Pool,2)
if (size(Pool,1).NE.NMonth) &
		print*, "  > ##### BuildPool: bad-len: Pool #####"
if (size(PoolWei,1).NE.NMonth.OR.size(PoolWei,2).NE.NLevel) &
		print*, "  > ##### BuildPool: bad-len: PoolWei #####"
Pool=MissVal ; PoolWei=MissVal

allocate (Exe(NYear),Wye(NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: BuildPool: alloc #####"

XNeigh=1 ; NHalt=0 ; Score=ScoreInit ; Level=0

if (Neigh(CandStn,1).NE.MissVal) then
 do
  XStn=Neigh(CandStn,XNeigh) ; XPer=Stn2Per(XStn)
  
  do
    if (PerTrust(XPer)) then			! period is trustworthy
     Joint0=CandBeg ; if (Joint0.LT.PerBeg(XPer)) Joint0=PerBeg(XPer)
     Joint1=CandEnd ; if (Joint1.GT.PerEnd(XPer)) Joint1=PerEnd(XPer)
    
     if ((Joint1-Joint0).GE.(Elastic-1)) then	! there is a worthwhile overlap
      do XMonth = 1, NMonth			! iterate by month
       if (Score(XMonth).GT.0) then		! ignore month if halted
        Exe=MissVal ; Wye=MissVal		! fill exe,wye
        do XYear = Joint0, Joint1
          if (FD(XYear,XMonth,XPer).NE.DataMissVal.AND. &
	      CandFD(XYear,XMonth).NE.DataMissVal) then
            Exe(XYear)=FD(XYear,XMonth,XPer)
            Wye(XYear)=CandFD(XYear,XMonth)
          end if
        end do
        
        call RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,(Elastic-1), &
				Beg=Joint0,End=Joint1)
	
	if      (Are.GT.AreThresh) then		! correlated, above threshold
	  Level(XMonth)=Level(XMonth)+1
	  Pool    (XMonth,Level(XMonth)) = XPer
	  PoolWei (XMonth,Level(XMonth)) = Are**2
	  
	  if (Score(XMonth).LT.ScoreInit) Score(XMonth)=Score(XMonth)+1
	  if (Level(XMonth).EQ.NLevel) then	! no more levels available
	    Score(XMonth)=0 ; NHalt=NHalt+1
	  end if
	else if (Enn.GT.(Elastic*2)) then	! despite good support
	  Score(XMonth)=Score(XMonth)-1
	  if (Score(XMonth).EQ.0) NHalt=NHalt+1	! bad prospect of sig coming
	end if
	
       end if
      end do
     end if
    end if
    
    XPer=XPer+1
    if (XPer.GT.NPer) exit			! no more periods
    if (Per2Stn(XPer).NE.XStn) exit		! no more periods for this neighbour
    if (NHalt.EQ.12) exit			! all months halted
  end do
  
  XNeigh=XNeigh+1
  if (XNeigh.GT.NNeigh) exit			! no more neighbours
  if (Neigh(CandStn,XNeigh).EQ.MissVal) exit	! no more valid neighbours
  if (NHalt.EQ.12) exit				! all months halted
 end do
end if

deallocate (Exe,Wye, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: BuildPool: dealloc #####"

end subroutine BuildPool

!*******************************************************************************
! returns the Abs,FD time-series for the candidate (NYear,NMonth), 
!   plus beg/end indices in CandBeg,CandEnd
! there may be missing stretches between the start and end

subroutine GetCand (Abs,FD,CandAbs,CandFD,Stn2Per,Per2Stn,PerBeg,PerEnd, &
		      CandStn,CandBeg,CandEnd,CandValid,NeedBeg,NeedEnd)

logical, pointer, dimension (:)		:: CandValid
integer, pointer, dimension (:,:,:)	:: FD,Abs
integer, pointer, dimension (:,:)	:: CandFD,CandAbs
integer, pointer, dimension (:)		:: Stn2Per,Per2Stn,PerBeg,PerEnd

integer, intent(in), optional		:: NeedBeg,NeedEnd 	! stretch must be incl
integer, intent(in)			:: CandStn
integer, intent(out)			:: CandBeg,CandEnd
integer, parameter 			:: DataMissVal = -9999

integer :: NPer,NYear,NMonth,XPer,XYear,XMonth

NYear=size(FD,1) ; NMonth=12 ; NPer=size(FD,3)
if (size(CandFD,1).NE.NYear.OR.size(CandFD,2).NE.NMonth) &
		print*, "  > ##### GetCand: bad-len: CandFD #####"

CandBeg=0 ; CandEnd=0 ; CandFD=DataMissVal ; CandAbs=DataMissVal
CandValid=.FALSE.
XPer=Stn2Per(CandStn)			! first period for candidate stn

do
  if (CandBeg.EQ.0) CandBeg=PerBeg(XPer)
  if (CandEnd.LT.PerEnd(XPer)) CandEnd=PerEnd(XPer)
  
  do XYear=PerBeg(XPer),PerEnd(XPer)
    CandValid(XYear)=.TRUE.
    do XMonth = 1, NMonth
      CandAbs(XYear,XMonth) = Abs(XYear,XMonth,XPer)
      CandFD (XYear,XMonth) = FD (XYear,XMonth,XPer)	! PerBeg=DMV, so ok
    end do
  end do
  
  XPer=XPer+1
  if (XPer.GT.NPer) exit		! no more periods
  if (Per2Stn(XPer).NE.CandStn) exit	! next period is for different stn
end do
	! this forces (e.g.) 1961-90 to be included in the cand series
	! if the same option is set in FormPar, only a RefBeg/End is allowed
	!   that includes 1961-90
	! setting CandBeg/End here helps FormPar, but more importantly
	!   it ensures that throughout 1961-90 is treated as required
if (present(NeedBeg).AND.present(NeedEnd)) then
  if (CandBeg.GT.NeedBeg) CandBeg=NeedBeg
  if (CandEnd.LT.NeedEnd) CandEnd=NeedEnd
end if

end subroutine GetCand

!*******************************************************************************
! convert absolute array into first-difference array
! the beg/end vectors are adjusted accordingly
! there really should not be any missing values

subroutine ConvertFD (Abs,FD,PerBeg,PerEnd,Differ)

integer, pointer, dimension (:,:,:)	:: Abs,FD
integer, pointer, dimension (:)		:: PerBeg,PerEnd

logical, intent(in)	:: Differ		! true for (B-A), false for (100*B/A)
real, parameter 	:: MissVal = -999.0
integer, parameter 	:: DataMissVal = -9999

integer :: NYear,NMonth,NPer
integer :: XYear,XMonth,XPer

NYear=size(Abs,1) ; NMonth=12 ; NPer=size(Abs,3) ; FD=Abs
if (size(PerBeg,1).NE.NPer.OR.size(PerEnd,1).NE.NPer) &
		print*, "  > ##### ConvertFD: bad-len: beg/end #####"

if (Differ) then			! (B-A)
  do XPer=1,NPer
    do XMonth=1,NMonth				! plain difference
      do XYear=PerEnd(XPer),(PerBeg(XPer)+1),-1	! loop back is essential
        FD(XYear,XMonth,XPer) = FD(XYear,XMonth,XPer)-FD((XYear-1),XMonth,XPer)
      end do
      FD(PerBeg(XPer),XMonth,XPer)=DataMissVal	! else a spurious abs value
    end do
  end do
else					! (100*B/A)
  do XPer=1,NPer				! else divisions by zero
    do XMonth=1,NMonth
      do XYear=PerEnd(XPer),PerBeg(XPer),-1	! not essential, but matches
        if (FD(XYear,XMonth,XPer).EQ.0) FD(XYear,XMonth,XPer)=1
      end do
    end do
  end do 
  
  do XPer=1,NPer				! ratio*100
    do XMonth=1,NMonth
      do XYear=PerEnd(XPer),(PerBeg(XPer)+1),-1	! loop back is essential
        FD(XYear,XMonth,XPer) = nint(100.0*real(FD(XYear,XMonth,XPer)) &
        			/real(FD((XYear-1),XMonth,XPer)))
        if (FD(XYear,XMonth,XPer).LT.    1) FD(XYear,XMonth,XPer)=    1
        if (FD(XYear,XMonth,XPer).GT.10000) FD(XYear,XMonth,XPer)=10000
      end do
      FD(PerBeg(XPer),XMonth,XPer)=DataMissVal	! else a spurious abs value
    end do
  end do
end if
 
end subroutine ConvertFD

!*******************************************************************************
! this presumes all arrays are already filled
! it uses neighbours to infill gaps in Abs
! it iterates by neighbour, looking for one with a value, and sig. correl
! it uses a window of 25 values to calc the correlation
! if it finds that there are no close enough relationships, it infills with
!	the local (window) mean from the same station
! note that if the period includes missing values in the first,last years, these
!	will be in-filled; if this is not desirable, fix it in the beg/end vectors,
!	which probably means adjusting FindPer
! any subsequent and in-range estimates for the same calendar month will have 
!	any in-filled estimate included in the regression

subroutine FillGaps (Abs,Per2Stn,Stn2Per,PerBeg,PerEnd,Neigh,Differ)

integer, pointer, dimension (:,:,:)	:: Abs
integer, pointer, dimension (:,:)	:: Neigh
integer, pointer, dimension (:)		:: Per2Stn,Stn2Per,PerBeg,PerEnd
real, pointer, dimension (:)		:: Exe,Wye

logical, intent(in)	:: Differ

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

logical :: Success,Okay
real :: Aye,Bee,Enn,Are,Pea, OpTot,OpEn, AyeNo,BeeNo
integer :: AllocStat,Lower,Upper,Elastic
integer :: NYear,NMonth,NStn,NPer,NNeigh,NPane,NWindow,NFail
integer :: XYear,XMonth,XStn,XPer,XNeigh,XPane,XWindow,XFail,XNeighPer

Elastic=10 ; if (Differ) Elastic=5
NPane=24 ; if (Differ) NPane=12
NYear=size(Abs,1) ; NMonth=12 ; NPer=size(Abs,3) 
NStn=size(Neigh,1) ; NNeigh=size(Neigh,2) ; NWindow=(NPane*2)+1
Success = .TRUE.
! write (99,"(6i8)"), NYear,NMonth,NStn,NPer,NNeigh,NPane	

if (size(Per2Stn,1).NE.NPer.OR.size(Stn2Per).NE.NStn) &
		print*, "  > ##### FillGaps: bad-len: cross #####"
if (size(PerBeg,1).NE.NPer.OR.size(PerEnd,1).NE.NPer) &
		print*, "  > ##### FillGaps: bad-len: beg/end #####"

allocate (Exe(NWindow),Wye(NWindow), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FillGaps: alloc #####"

do XPer=1,NPer
  do XYear=PerBeg(XPer),PerEnd(XPer)
    do XMonth=1,NMonth
      if (Abs(XYear,XMonth,XPer).EQ.DataMissVal) then
        Wye=MissVal ; Lower=XYear-NPane ; Upper=XYear+NPane			
        if (Upper.GT.NYear) Upper=NYear
        if (Lower.LT.    1) Lower=1
        
        do XPane=Lower,Upper		! fill Wye
          XWindow=XPane-(XYear-NPane)+1
          if (Abs(XPane,XMonth,XPer).NE.DataMissVal) &
                    		Wye(XWindow) = real(Abs(XPane,XMonth,XPer))
        end do
              
        XNeigh=1 ; NFail=0
        do
          XNeighPer=MissVal ; XStn=MissVal
          if (Per2Stn(XPer).NE.MissVal) XStn=Neigh(Per2Stn(XPer),XNeigh)
          if (XStn.NE.MissVal) XNeighPer=Stn2Per(XStn)
          
          if (XNeighPer.NE.MissVal) then
           do
            if (Abs(XYear,XMonth,XNeighPer).NE.DataMissVal) then
              Exe=MissVal
              do XPane=Lower,Upper	! fill Exe
                XWindow=XPane+1-XYear+NPane
                if (Abs(XPane,XMonth,XNeighPer).NE.DataMissVal) &
                      Exe(XWindow) = real(Abs(XPane,XMonth,XNeighPer))
              end do
              
	      if (Differ) then		! y=a+bx for temperature
		call MergeForDiff  (Wye,Exe,1,NWindow)
	      else			! adjust neighbour to y=ax**b for precip
		call MergeForRatio (Wye,Exe,1,NWindow)
	      end if

              call RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,Elastic)
	      
              if      (Are.GE.0.2) then		! sig rel
		Abs(XYear,XMonth,XPer) = nint(Exe(NPane+1))
              else if (Enn.GT.(Elastic*2)) then	! insig, despite good opport
                NFail=NFail+1
              end if
            end if
            
            XNeighPer=XNeighPer+1
            if (XNeighPer.GT.NPer) exit			! last per for stn
            if (Per2Stn(XNeighPer).NE.XStn) exit 		! last per for stn
           end do
          end if
          
          XNeigh=XNeigh+1 
          if (NFail.GE.5) exit					! hope is gone
          if (XNeigh.GT.NNeigh) exit				! no more neigh
          if (XStn.EQ.MissVal) exit				! no more valid neigh
          if (Abs(XYear,XMonth,XPer).NE.DataMissVal) exit	! filled
        end do
        
        if (Abs(XYear,XMonth,XPer).EQ.DataMissVal) then	! not possible to use neigh
          OpTot=0 ; OpEn=0
          do XWindow=1,NWindow				! try to take local mean
            if (Wye(XWindow).NE.MissVal) then
              OpTot=OpTot+Wye(XWindow) ; OpEn=OpEn+1
            end if
          end do
          
          if (OpEn.GT.0) then				! take local mean
            Abs(XYear,XMonth,XPer)=nint(OpTot/OpEn)	
          else						! try to take full per mean
            OpTot=0 ; OpEn=0
            do XWindow=PerBeg(XPer),PerEnd(XPer)
              if (Abs(XWindow,XMonth,XPer).NE.DataMissVal) then
                OpTot=OpTot+Abs(XWindow,XMonth,XPer) ; OpEn=OpEn+1
              end if
            end do
            
            if (OpEn.GT.0) then				! take full per mean
              Abs(XYear,XMonth,XPer)=nint(OpTot/OpEn)	
            else
              write (99,"(a,6i6)"), "alg-fail ", &
            		Per2Stn(XPer),XYear,XMonth,XPer,PerBeg(XPer),PerEnd(XPer)
              print*, "  > ##### ERROR: FillGaps: alg failed #####"
              Success = .FALSE.
            end if
          end if
        end if
      end if
    end do
  end do
end do

deallocate (Exe,Wye, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FillGaps: dealloc #####"

end subroutine FillGaps

!*******************************************************************************
! fills Neigh with the indices of the nearest stns, in order, within corr decay dist
! the maximum number of neighbours is given by the 2nd dimension in Neigh
! only stns with valid pers (given by Stn2Per) are filled, or allowed to be neighbours

subroutine OrderNeigh (Neigh,Stn2Per,Per2Stn,PerTrust,PerBeg,PerEnd,Lat,Lon,Decay,Differ)

logical, pointer, dimension (:)		:: PerTrust,StnTrust
real, pointer, dimension (:)		:: Lat,Lon,Horiz
integer, pointer, dimension (:,:)	:: Neigh
integer, pointer, dimension (:)		:: Stn2Per,Per2Stn,PerBeg,PerEnd
integer, pointer, dimension (:)		:: Order,StnBeg,StnEnd
logical, intent(in)			:: Differ
real, intent(in)			:: Decay
real, parameter 			:: MissVal = -999.0

logical :: Overlap
integer :: NStn,NPer,NNeigh,NHoriz
integer :: XStn,XPer,XNeigh,XHoriz, AllocStat,Joint0,Joint1,Elastic

Elastic=10 ; if (Differ) Elastic=5
NStn=size(Neigh,1) ; NNeigh=size(Neigh,2) ; NPer=size(PerTrust,1)
if (size(Lat,1).NE.NStn.OR.size(Lon,1).NE.NStn) &
		print*, "  > ##### OrderNeigh: bad-len: Lat/Lon #####"
if (size(Stn2Per,1).NE.NStn) print*, "  > ##### OrderNeigh: bad-len: StnPer #####"

allocate (Horiz(NStn),Order(NStn),StnTrust(NStn), &
	  StnBeg(NStn),StnEnd(NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OrderNeigh: Alloc: Horiz #####"
Neigh=MissVal ; Order=MissVal ; StnTrust=.FALSE. ; stnBeg=MissVal ; StnEnd=MissVal

do XStn=1,NStn						! iterate by station
  XPer=Stn2Per(XStn)
  do
    if (XPer.EQ.MissVal.OR.XPer.GT.NPer) exit
    if (Per2Stn(XPer).NE.XStn) exit
    if (StnBeg(XStn).EQ.MissVal) StnBeg(XStn)=PerBeg(XPer)	! stn beg
    StnEnd(XStn)=PerEnd(XPer)					! stn end
    if (PerTrust(XPer)) StnTrust(XStn)=.TRUE.		! stn has >0 trusted pers
    XPer=XPer+1
  end do
end do

do XStn = 1, NStn
 Horiz=MissVal ; NHoriz=0
 
 if (Stn2Per(XStn).NE.MissVal.AND. &
 		Lat(XStn).NE.MissVal.AND.Lon(XStn).NE.MissVal) then
  do XHoriz=1,NStn					! iterate by potential neighbour
    if (StnTrust(XHoriz).AND.XHoriz.NE.XStn.AND. &	! must be >0 trusted periods
    		Lat(XHoriz).NE.MissVal.AND.Lon(XHoriz).NE.MissVal) then
      XPer=Stn2Per(XHoriz) ; Overlap=.FALSE.
      do						! look for overlap + trusted per
        if (XPer.EQ.MissVal.OR.XPer.GT.NPer.OR.Overlap) exit
        if (Per2Stn(XPer).NE.XHoriz) exit
        if (PerTrust(XPer)) then
          Joint0=PerBeg(XPer) ; if (Joint0.LT.StnBeg(XStn)) Joint0=StnBeg(XStn)
          Joint1=PerEnd(XPer) ; if (Joint1.GT.StnEnd(XStn)) Joint1=StnEnd(XStn)
          if ((Joint1-Joint0).GE.(Elastic-1)) Overlap=.TRUE.
        end if
        XPer=XPer+1
      end do
      
      if (Overlap) then					! examine distance
       Horiz(XHoriz) = GetDistance (Lat(XStn),Lon(XStn),Lat(XHoriz),Lon(XHoriz))
      
       if (Horiz(XHoriz).GT.Decay) then
        Horiz(XHoriz)=MissVal
       else						! if <decay, potential neigh
        NHoriz=NHoriz+1
       end if
      end if
    end if
  end do
  
  if (NHoriz.GT.0) then
    call QuickSort (Reals=Horiz,OrderValid=Order)
  
    do XHoriz=1,NHoriz
      if (XHoriz.LE.NNeigh) Neigh(XStn,XHoriz)=Order(XHoriz)
    end do
  end if
 end if
end do

deallocate (Horiz,Order,StnTrust,StnBeg,StnEnd, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OrderNeigh: Dealloc: Horiz #####"

end subroutine OrderNeigh

!*******************************************************************************
! initialise Abs and associated vectors (already allocated)
! fills with absolute values
! can have gaps filled and turned into genuine first diffs later

subroutine InitialAbs (Orig,StnPer,StnPerTrust,Abs,Per2Stn,Stn2Per, &
		       PerTrust,PerBeg,PerEnd)

logical, pointer, dimension (:,:)	:: StnPerTrust
logical, pointer, dimension (:)		:: PerTrust

integer, pointer, dimension (:,:,:)	:: Orig,Abs
integer, pointer, dimension (:,:)	:: StnPer
integer, pointer, dimension (:)		:: Per2Stn,Stn2Per,PerBeg,PerEnd

integer, parameter 	:: DataMissVal = -9999
real, parameter 	:: MissVal = -999.0
integer :: NYear,NMonth,NStn,NPer,NBreak,XYear,XMonth,XStn,XPer,XBreak

NYear=size(Orig,1) ; NMonth=12 ; NStn=size(Orig,3) ; NPer=size(Abs,3)
if (size(StnPer,1).NE.NStn.OR.size(StnPer,2).NE.NYear) &
		print*, "  > ##### InitialAbs: bad-len: StnPer #####"
if (size(Abs,1).NE.NYear.OR.size(Abs,2).NE.12) &
		print*, "  > ##### InitialAbs: bad-len: Abs #####"
if (size(Per2Stn,1).NE.NPer.OR.size(Stn2Per,1).NE.NStn) &
		print*, "  > ##### InitialAbs: bad-len: crossover #####"
if (size(PerBeg,1).NE.NPer.OR.size(PerEnd,1).NE.NPer) &
		print*, "  > ##### InitialAbs: bad-len: Beg/End #####"
Abs=DataMissVal ; Stn2Per=MissVal ; Per2Stn=MissVal 
PerBeg=MissVal ; PerEnd=MissVal ; PerTrust=.FALSE.

XPer=0
do XStn = 1, NStn
  XBreak=1
  do
    if (StnPer(XStn,XBreak).EQ.MissVal) exit

    if (StnPer(XStn,XBreak).NE.MissVal) then
      XPer=XPer+1				! set new period
      Per2Stn(XPer)=XStn			! identify stn for this period
      if (XBreak.EQ.1) Stn2Per(XStn)=XPer	! identify first period for this stn
      PerBeg(XPer)=StnPer(XStn,XBreak)		! identify period beg
      PerEnd(XPer)=StnPer(XStn,(XBreak+1))	! identify period end
      PerTrust(XPer)=StnPerTrust(XStn,XBreak)	! identify period trustworthiness
!      write (99,"(7i5)"), XStn,XPer,Per2Stn(XPer),Stn2Per(XStn), &
!      			PerBeg(XPer),PerEnd(XPer),PerTrust(XPer)	! @@@@@@@@@@@@@@@
      
      do XYear=PerBeg(XPer),PerEnd(XPer)	! fill Abs with absolute values
        do XMonth=1,NMonth
          Abs(XYear,XMonth,XPer) = Orig(XYear,XMonth,XStn)
        end do
      end do
    end if

    XBreak=XBreak+2
  end do
end do

end subroutine InitialAbs

!*******************************************************************************
! Orig is the absolute climate data-set of interest
! StnPer contains the beg/end points between independent periods for each station
! Lat and Lon are checked because without them any per is useless
! if Disc is included, any discon years are prohibited from being in mid-period 
! if CandStn+Decay are set, only stns within Decay range of CandStn have pers
! 	or if Calc+Decay are set, ditto for stns within range of .TRUE. in Calc

subroutine FindPer (Orig,StnPer,StnPerTrust,NPer,Lat,Lon,Trust,Differ, &
		    CandStn,Decay,Disc,Calc)

logical,pointer,dimension(:,:)		:: Trust	! Nyear,NStn
logical,pointer,dimension(:,:)		:: StnPerTrust
logical,pointer,dimension(:,:),optional	:: Disc		! Nyear,NStn
logical,pointer,dimension(:),  optional	:: Calc		! NStn
logical,allocatable,dimension(:)	:: Need		! NStn
real, pointer, dimension (:)		:: Lat,Lon
integer, pointer, dimension (:,:,:)	:: Orig
integer, pointer, dimension (:,:)	:: StnPer
integer, dimension (12)			:: Current,Total,Pre,Post
logical, intent(in)			:: Differ
real, intent(in), optional 		:: Decay
integer, intent(in), optional 		:: CandStn
integer, intent(out) 	:: NPer
integer, parameter 	:: DataMissVal = -9999
real, parameter 	:: MissVal = -999.0
real    :: Distance
logical :: ForceEnd,Reset,Revisit,Gambit
integer :: NYear,NMonth,NStn,XYear,XMonth,XStn,XPer,XCalc
integer :: Sum,Beg,End,NuffMonths,WaryPerCount,TrustPerCount,AllocStat
integer :: TrustBeg,TrustEnd,CheckBeg,CheckEnd, Elastic
integer :: PerBeg,PerEnd,CheckPer,Break,PreMin,PostMin

Elastic=10 ; if (Differ) Elastic=5
NYear=size(Orig,1) ; NMonth=12 ; NStn=size(Orig,3) 
if (size(StnPer,1).NE.NStn.OR.size(StnPer,2).NE.NYear) &
		print*, "  > ##### FindGaps: bad-len #####"
StnPer=MissVal ; ForceEnd=.FALSE. ; StnPerTrust=.FALSE.
WaryPerCount=0 ; TrustPerCount=0 ; NPer=0

allocate (Need(NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FindPer: alloc: Need #####"
Need=.TRUE.

if (present(Decay)) then
  if      (present(CandStn)) then
    Need=.FALSE.
    do XStn=1,NStn
      if (Lat(XStn).NE.MissVal.AND.Lon(XStn).NE.MissVal) then
        Distance = GetDistance(Lat(XStn),Lon(XStn),Lat(CandStn),Lon(CandStn))
	if (Distance.NE.MissVal.AND.Distance.LE.Decay) Need(XStn)=.TRUE.
      end if
    end do
  else if (present(Calc)) then
    Need=.FALSE.
    do XCalc=1,NStn
      if (Calc(XCalc)) then
        Need(XCalc)=.TRUE.
	do XStn=1,NStn
	  if (.NOT.Need(XStn)) then
	    Distance = GetDistance(Lat(XStn),Lon(XStn),Lat(XCalc),Lon(XCalc))
	    if (Distance.NE.MissVal.AND.Distance.LE.Decay) Need(XStn)=.TRUE.
	  end if
	end do
      end if
    end do
  end if
end if

do XStn = 1, NStn
 if (Need(XStn)) then
  TrustBeg=MissVal ; TrustEnd=MissVal
  do XYear=1,NYear
    if (Trust(XYear,XStn)) then
      if (TrustBeg.EQ.MissVal) TrustBeg=XYear
      TrustEnd=XYear
    end if
  end do
  
  XPer=1 ; Beg=0 ; End=0 ; Current=0; Total=0 ; NuffMonths=0
  Reset=.FALSE. ; XYear=0
  
  do XYear=1,NYear
    if (present(Disc)) ForceEnd=Disc(XYear,XStn)
    
    Sum=0
    do XMonth = 1, NMonth
      if (Orig(XYear,XMonth,XStn).NE.DataMissVal) then
        if (Current(XMonth).LT.Elastic) Current(XMonth) = Current(XMonth) + 1
        Sum = Sum + Current(XMonth)
        Total(XMonth)=Total(XMonth)+1
        if (Total(XMonth).EQ.Elastic) NuffMonths=NuffMonths+1
      else if (Current(XMonth).GT.0) then
        Current(XMonth) = Current(XMonth) - 1
        Sum = Sum + Current(XMonth)
      end if
    end do
    				! (un)trusted period ended
    if (ForceEnd) then	! discontinuity
     if (Beg.EQ.0) then
       Reset=.TRUE.				! no current per, so just restart
     else if (XYear.LT.(Beg+Elastic-1).OR.NuffMonths.LT.12) then
       Reset=.TRUE.				! hardly begun, so just restart
     else
       End=XYear				! decent per, so make this a per
     end if
    else			! not a discontinuity
     if (Beg.EQ.0) then
      if (Sum.GT.0) then
        Beg=XYear
      end if
     else if (Sum.EQ.0) then
      if (XYear.LT.(Beg+Elastic-1).OR.NuffMonths.LT.12) then	
          Reset=.TRUE.
      else
          End=XYear-Elastic					
      end if
     else if (XYear.EQ.NYear) then
      if (XYear.GT.(Beg+Elastic-2).AND.NuffMonths.EQ.12) then				
          End=XYear
      end if
     end if
    end if
    
    if (End.GT.0) then		! store this as a period
      StnPer(XStn,XPer)=Beg ; StnPer(XStn,(XPer+1))=End
      if (Beg.GE.TrustBeg.AND.End.LE.TrustEnd) then	! all of per is trusty
        StnPerTrust(XStn,XPer)=.TRUE. ; StnPerTrust(XStn,XPer+1)=.TRUE.
        TrustPerCount=TrustPerCount+1
      else
        WaryPerCount=WaryPerCount+1
      end if
      XPer=XPer+2 ; NPer=NPer+1 ; Reset=.TRUE.
    end if
    
    if (Reset) then		! reset for new period search
      Beg=0 ; End=0 ; Current=0 ; Total=0 ; NuffMonths=0 ; Reset=.FALSE.
    end if
  end do
  				
  if (TrustBeg.NE.MissVal) then	! deal with any mid-per change-of-trust
    CheckPer=1 ; Revisit=.FALSE.
    do				! loop by per
      if (StnPer(XStn,CheckPer).EQ.MissVal) exit	! no more pers for stn
      PerBeg=StnPer(XStn,CheckPer) ; PerEnd=StnPer(XStn,CheckPer+1)
      
      				! examine if change-of-trust within per
      if (	(TrustBeg.GT.PerBeg.AND.TrustBeg.LT.PerEnd) .OR. &
          	(TrustEnd.GT.PerBeg.AND.TrustEnd.LT.PerEnd) ) then
        CheckBeg=PerBeg ; CheckEnd=PerEnd	! range to check for valid per
        if (TrustBeg.GT.PerBeg.AND.TrustBeg.LT.PerEnd &
        	.AND.(.NOT.Revisit)) then
          Break=TrustBeg ; Gambit=.TRUE.	! Gambit: break is beg not end
          if (TrustEnd.GT.PerBeg.AND.TrustEnd.LT.PerEnd) then
            CheckEnd=TrustEnd ; Revisit=.TRUE.	! Revisit: return for TrustEnd
          end if
        else
          Break=TrustEnd+1 ; Gambit=.FALSE. ; Revisit=.FALSE.
        end if
        
        Pre=0 ; Post=0	
        do XYear=CheckBeg,(Break-1)
            do XMonth=1,NMonth		! find info pre-break
              if (Orig(XYear,XMonth,XStn).NE.DataMissVal) &
            		Pre(XMonth)=Pre(XMonth)+1
            end do
        end do
        do XYear=Break,CheckEnd
            do XMonth=1,NMonth		! find info post-break
              if (Orig(XYear,XMonth,XStn).NE.DataMissVal) &
            		Post(XMonth)=Post(XMonth)+1
            end do
        end do
        PreMin=minval(Pre) ; PostMin=minval(Post)
          		
        if      (PreMin.GE.Elastic.AND.PostMin.GE.Elastic) then
            StnPer(XStn,CheckPer)=PerBeg ; StnPer(XStn,CheckPer+1)=Break-1
            StnPer(XStn,XPer)=Break  ; StnPer(XStn,XPer+1)=PerEnd
            
            if (Gambit) then
              StnPerTrust(XStn,XPer)=.TRUE. 
              StnPerTrust(XStn,XPer+1)=.TRUE.
            else
              StnPerTrust(XStn,CheckPer)=.TRUE. 
              StnPerTrust(XStn,CheckPer+1)=.TRUE.
            end if

            NPer=NPer+1 ; XPer=XPer+2  		! separate pers
            TrustPerCount=TrustPerCount+1	! one more trusted per
        else if (PreMin.LT.Elastic.AND.PostMin.LT.Elastic) then
	    					! single, already untrusted per
        else if (PreMin.GE.Elastic.AND.PostMin.LT.Elastic) then
            if (Gambit) then
              					! single, already untrusted per
            else
              StnPerTrust(XStn,CheckPer)=.TRUE. ! just trust the extra data
              StnPerTrust(XStn,CheckPer+1)=.TRUE.
              TrustPerCount=TrustPerCount+1
              WaryPerCount=WaryPerCount-1
            end if
        else if (PreMin.LT.Elastic.AND.PostMin.GE.Elastic) then
            if (Gambit) then
              StnPerTrust(XStn,CheckPer)=.TRUE. ! just trust the extra data
              StnPerTrust(XStn,CheckPer+1)=.TRUE.
              TrustPerCount=TrustPerCount+1
              WaryPerCount=WaryPerCount-1
            else
              					! single, already untrusted per
            end if
        end if
      else
        Revisit=.FALSE.
      end if

      if (Break.EQ.MissVal) exit

      if (.NOT.Revisit) CheckPer=CheckPer+2 
    end do
  end if
 end if
end do

deallocate (Need, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FindPer: dealloc: Need #####"

! write (99,"(a,2i8)"), "periods trusted, untrusted: ",TrustPerCount,WaryPerCount
! print "(a,2i7)", "   > Periods trusted, not: ",TrustPerCount,WaryPerCount

end subroutine FindPer

!*******************************************************************************
! use the refTS to decide which portions of which stns can be trusted sufficiently
!   to be allowed into the refTS in the future.
! only stns for which a RefTS has been calculated are checked (i.e. Calc=true)
! stn-yrs without a refTS become untrusted (i.e. Trust=false)
! if Verbose: print the number of stn-yrs that are available to be trusted
!   and (a) have a ref ts and are therefore trustworthy, and (b) do not and are not

subroutine Trustworthy (Orig,RefTS,Calc,Trust)

logical,pointer,dimension(:,:)		:: Trust	! year,stn
logical,pointer,dimension(:)		:: Calc		! stn
integer, pointer, dimension (:,:,:) 	:: Orig,RefTS

integer, parameter :: DataMissVal = -9999
integer :: NYear,NMonth,NStn, Beg,End
integer :: XYear,XMonth,XStn

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

! write (99,"(a)"), "TRUSTWORTHY"	! @@@@@@@@@@@@@@@@@@@@@@@@
do XStn=1,NStn
 if (Calc(XStn)) then
  XYear=0 ; Beg=0 ; End=0
  do
    XYear=XYear+1
    if (RefTS(XYear,1,XStn).NE.DataMissVal) then	! checked
      if (Beg.EQ.0) then
        Beg=XYear			! start of per
      else if (XYear.EQ.NYear) then
        End=XYear			! end of per
      end if
    else if (Beg.GT.0) then				! unchecked
      End=XYear-1			! end of per
    end if
    
    if (XYear.EQ.NYear.OR.End.GT.0) exit
  end do
   
  do XYear=1,NYear
    if (XYear.LT.Beg.OR.XYear.GT.End) then
      Trust(XYear,XStn)=.FALSE.
    else
      Trust(XYear,XStn)=.TRUE.
    end if
  end do
  
!  write (99,"(3i6)"), XStn,Beg,End	! @@@@@@@@@@@@@@@@@@@@@@@@
 end if
end do

end subroutine Trustworthy

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

end module GHCNrefIter