! homogiter.f90
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/homogiter time.f90 filenames.f90 crutsfiles.f90 
!	perfiles.f90 cetgeneral.f90 gridops.f90 sortmod.f90 plainperm.f90 regress.f90 
!	fdist_k.f90 fdist_br.f90 fdist.f90 ttest.f90 ghcnrefiter.f90 ghcndiscon.f90 
!	./../cruts/homogiter.f90 2> /tyn1/tim/scratch/stderr.txt

program Homogiter

use Time
use FileNames
use CRUtsFiles
use PerFiles
use CETGeneral
use GridOps
use SortMod
use PlainPerm
use Regress
use FDist_K
use FDist_BR
use FDist
use TTest
use GHCNrefIter
use GHCNdiscon

implicit none

logical, pointer, dimension (:,:)	:: Disc,DiscCrude,Trust,Believe,Got,Need
logical, pointer, dimension (:)		:: Examine,Sought,GotRef

real, pointer, dimension (:,:)		:: AdjMon,AdjSea
real, pointer, dimension (:)		:: OLat,OLon,OElv, AdjStn,AdjAnn
real, pointer, dimension (:)		:: Lat,Lon,Elv, DLat,DLon,DElv

integer, pointer, dimension (:,:,:)	:: Data,DData,OData,PData,RefTS,DumpTarget,Adj
integer, pointer, dimension (:,:,:)	:: DataSrc,ODataSrc,AdjSrc
integer, pointer, dimension (:,:)	:: Info,DInfo,OInfo,DNormals,TrashInfo
integer, pointer, dimension (:)		:: YearAD,DYearAD,PYearAD,TrashYearAD
integer, pointer, dimension (:)		:: Code,OCode,DCode,PCode,Order,OldCode,OOldCode
integer, pointer, dimension (:)		:: EarlBeg,EarlEnd,LateBeg,LateEnd,ThisBeg,ThisEnd
integer, allocatable, dimension (:)	:: StnLen
integer, dimension (12)			:: Climatology,DumpData

character (len=80), pointer, dimension (:) :: Batch
character (len=20), pointer, dimension (:) :: Name,DName,OName,TrashName
character (len=13), pointer, dimension (:) :: Cty,DCty,OCty,TrashCty
character (len=09), pointer, dimension (:) :: Local,DLocal,OLocal,TrashLocal

real, parameter 	:: Missval = -999.0
integer, parameter 	:: DataMissval = -9999

logical :: Verbose,PlotAllManip,Differ,Plotted,QBlank,Match,Clash,ForceLate

real :: Aye,Bee,Enn,Are,Pea,Decay,TTestResult,OpTot,OpEn,OpMean,OpWin,OpAdd,OpValid
real :: OpOData,OpRefTS,OpAdj,OpClim
real :: PlotThresh,Req,LatTot,LonTot,DistTot,DistTotSq,Dist,DistMax
real :: CentroidLat,CentroidLon,CentroidRange, BoxLatDeg,BoxLonDeg
real :: RealTarget,RealClim,RealOData,RealRefTS,RealAdj

integer :: AllocStat,ReadStatus, Beg,End,DBeg,DEnd,KBeg,KEnd,PBeg,PEnd
integer :: Dropped,IterInt,Omit,ExtraTerm
integer :: QDump,QAbs,QAnom,QStep,QDumpView,QOneCont,QUseDStn
integer :: QCountOnly,QType,QRetry
integer :: NYear,NDYear,NMonth,NStn,NDStn,NOStn,NIterate,NDisc,NTrust,NGotRef,NSought
integer :: XYear,XDYear,XMonth,XStn,XDStn,XOstn,XIterate,XDisc,XTrust,XGotRef,XSought
integer :: NBatch,NInfo,NWindow,NMiss,NExam,NOYear,NLat,NLon,NGot,NNeed,NPYear,NPStn
integer :: XBatch,XInfo,XWindow,XMiss,XExam,XOYear,XLat,XLon,XGot,XNeed,XPYear,XPStn
integer :: BegO,EndO,BegR,EndR, SubBeg,SubLen, Log0,Log1,LogAdj,LogNon, OBeg,OEnd
integer :: OuterLen,Elastic

character (len=99) :: Command
character (len=80) :: DataFile,DatabaseFile,Title,GivenFile,DumpFile,PerFile,HdrFile
character (len=80) :: DataSrcFile,SaveSrcFile,SaveFileOne,SaveFileTwo,ParallelFile
character (len=20) :: String20,IterStr,CallForm
character (len=10) :: FileTime
character (len= 8) :: StnText,DumpText,FileDate,FileText
character (len= 7) :: ScreenText
character (len= 4) :: IterText,VariableSuffix

call Main

contains

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

subroutine Main

call Initialise
call SpecifyCommon

print*, "  > Loading..."
if (DatabaseFile.NE."none") then
  call LoadDatabase
else
  call LoadDatabaseDummy
end if  
call LoadOriginal
call GetCoverage

QCountOnly=1 ; NStn=NOStn ; call Screen
call SpecComb
QcountOnly=0 ; XStn=NOStn ; call Screen
call MakeArrays
if (QRetry.GT.0) call LoadParallel

call Process
call SaveTwo

if      (QDumpView.EQ.1) then
  print*, "  > Plotting stns with discon..."
  call PlotDiscStn
else if (QDumpView.EQ.2) then
  print*, "  > Dumping stns with discon..."
  call DumpDiscStn
end if

if (ExtraTerm.EQ.1) call system ('printf "end" >> ./../../../scratch/f90toidl.pipe')
print*
close (99)

end subroutine Main

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

subroutine Process

print*, "  > Processing..."
write (99,"(i4,a,i4)"), YearAD(1), "-", YearAD(NYear)
do XStn=1,NStn
  write (99,"(i4,x,a20,x,a13,2f8.2)"), XStn,Name(XStn),Cty(XStn),Lat(XStn),Lon(XStn)
end do

call MakeRefTS (RefTS,Data,Lat,Lon,Sought,Trust,Decay,Differ, &
  			.FALSE.,nint(MissVal),GotRef)
call Trustworthy (Data,RefTS,Sought,Trust)	! ghcnrefiter.f90
call MakeContinuous (Data,RefTS,Sought,DiscCrude,YearAD,Differ,Verbose)

NGotRef=count(GotRef) ; NSought=NOStn ; NDisc=count(DiscCrude)
OuterLen=maxval(StnLen,Sought)
print "(a,4a7)", "   > ITERATIONS (iter gaps)","    GOT","  TO-DO","  DISCS"," MAXLEN"
print "(a,4i7)", "   > readied               ",NGotRef,NSought,NDisc,OuterLen

Omit=0 ; XIterate=0 ; RefTS=DataMissVal ; ForceLate=.TRUE.
do			! MAIN ITERATION LOOP
  XIterate=XIterate+1			! calc quality ref ts, but with restrictions
  
  if (ForceLate) then
    call MakeRefTS (RefTS,Data,Lat,Lon,Sought,Trust,Decay,Differ, &
  		.TRUE.,Omit,GotRef,Disc=DiscCrude, &
		ForceBeg=LateBeg,ForceEnd=LateEnd,MakeMore=.TRUE.)
  else
    call MakeRefTS (RefTS,Data,Lat,Lon,Sought,Trust,Decay,Differ, &
  		.TRUE.,Omit,GotRef,Disc=DiscCrude, &
		ForceBeg=EarlBeg,ForceEnd=EarlEnd,MakeMore=.TRUE.)
  end if
  
  NGotRef=count(GotRef) ; NSought=NSought-NGotRef 
  if (count(DiscCrude).NE.NDisc) print*, "  > @@@@@ ERROR: Disc error @@@@@" 
  write (99,"(a,2i4,3i7)"), "CHECK",XIterate,Omit,NGotRef,NSought,NDisc
  
  if (NGotRef.GT.0) then	! have got more stns with ref ts
    					! for these stns, correct original Data
    call MakeContinuous (Data,RefTS,GotRef,Disc,YearAD,Differ,Verbose,Adj=Data)
    					! ... and make checked part trusted
    call Trustworthy    (Data,RefTS,GotRef,Believe)
	
    do XStn=1,NStn
      if (GotRef(XStn)) then		! stn has ref ts
        Sought(XStn)=.FALSE.		! so no longer needs checking
        
	do XYear=1,NYear
          if (Believe(XYear,XStn)) then	! for the years with a valid ref ts
	    Trust(NYear,NStn)=.TRUE.	! the series may enter other ref ts
	    if (DiscCrude(XYear,XStn)) then
              DiscCrude(XYear,XStn)=.FALSE.	! any discon have been healed	
              NDisc=NDisc-1
            end if
	  end if
        end do
      end if
    end do
  
    if (ForceLate) then
      ForceLate=.FALSE.
    else
      ForceLate=.TRUE.
    end if
    
    ScreenText="checked"
    print "(2a,2(a,i4),a,4i7)", "   > ", ScreenText, "    (", &	
    		XIterate,",",Omit,")",NGotRef,NSought,NDisc,OuterLen
  else
    Omit=Omit+5				! relax restrictions one notch
    ForceLate=.TRUE.
  
    do XStn=1,NStn
      if (Sought(XStn)) then
        if ((StnLen(XStn)-Omit).LE.(2*Elastic)) then
          Sought(XStn)=.FALSE. ; NSought=NSought-1
        end if
      end if
    end do
    
    ScreenText="relaxed" ; NGotRef=MissVal
  end if
  
  if (NSought.GT.0) then
      OuterLen=maxval(StnLen,Sought)
  else
      OuterLen=0
  end if
    
  if (NSought.EQ.0.OR.Omit.GE.OuterLen) exit
end do

allocate (Adj(NYear,NMonth,NStn),AdjMon(NYear,NMonth),AdjStn(NStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### Process: alloc Adj* #####"
Adj=DataMissVal ; AdjMon=MissVal ; AdjStn=MissVal

call FilterOrig (Data,RefTS,Adj,Examine,Differ)
call MeasureAdjAll (OData,Adj,RefTS,AdjMon,AdjStn,Examine)

end subroutine Process

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

subroutine Initialise

open (99,file="./../../../scratch/log-homogiter.dat", &
		status="replace",action="write")

print*
print*, "  > ***** Homogiter: homogenises new stn files *****"
print*, "  > cleanmeta.f90 must be run first, and outputs used here"
print*

NMonth=12

end subroutine Initialise

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

subroutine SpecifyCommon

print*, "  > Extra terminals (log,idl) ? (0=no,1=yes,2=dump) "
do
  read (*,*,iostat=ReadStatus), ExtraTerm
  if (ReadStatus.LE.0.AND.ExtraTerm.GE.0.AND.ExtraTerm.LE.2) exit
end do

if (ExtraTerm.EQ.1) then
  Command='xterm -iconic -geometry 90x20+500+0 -T log -e tail -f' // &
	' ./../../../scratch/log-homogenise.dat &'
  call system (Command)
  call system ('xterm -iconic -geometry 80x40+0+20 -T idl -e idl startup &')
end if

print*, "  > Select the .cts file to load: "
do
  read (*,*,iostat=ReadStatus), GivenFile
  if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
DataFile = LoadPath(GivenFile,".cts")
SubBeg = index(DataFile,".cts") - 4
VariableSuffix = DataFile(SubBeg:SubBeg+3)
QRetry = index(DataFile,"retry")

if (QRetry.GT.0) then
  SubBeg = index(DataFile,"/",.TRUE.) ; SubLen = len_trim(DataFile)
  DataFile((SubLen-8):(SubLen-8)) = "?"
  call GetBatch (DataFile,Batch,Silent=1) 
  NBatch=size(Batch) ; DataFile = Batch(NBatch)
  deallocate (Batch,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### SpecifyCommon: dealloc: Batch #####"
  IterStr = DataFile((SubLen-8):(SubLen-8))
  IterInt = GetIntFromText (IterStr)
  IterInt = IterInt + 1
  IterStr = GetTextFromInt (IterInt)
  GivenFile = "./../../../data/cruts/homog" // DataFile(SubBeg:(SubLen-9)) &
  		// trim(adjustl(IterStr)) // DataFile((SubLen-7):SubLen)
else
  SubBeg = index(DataFile,"/",.TRUE.) ; SubLen = len_trim(DataFile)
  GivenFile = "./../../../data/cruts/homog" // DataFile(SubBeg:(SubLen-8)) &
  		// ".1" // DataFile((SubLen-7):SubLen)
end if

print "(2a)", "   > Save file: ", trim(GivenFile)
SaveFileTwo = SavePath (GivenFile,".cts")

if      (VariableSuffix.EQ.'.tmp') then
  Decay=1200.0 ; Differ=.TRUE. ; Elastic=5
else if (VariableSuffix.EQ.'.dtr') then
  Decay= 750.0 ; Differ=.TRUE. ; Elastic=5
else if (VariableSuffix.EQ.'.pre') then
  Decay= 450.0 ; Differ=.FALSE. ; Elastic=10
else if (VariableSuffix.EQ.'.vap') then
  Decay=1000.0 ; Differ=.TRUE. ; Elastic=5
else if (VariableSuffix.EQ.'.cld') then
  Decay= 600.0 ; Differ=.TRUE. ; Elastic=5
else if (VariableSuffix.EQ.'.spc') then
  Decay= 600.0 ; Differ=.TRUE. ; Elastic=5
else
  print*, "  > @@@@@ ERROR: variable suffix not included: ", VariableSuffix
end if

DatabaseFile = "./../../../data/cruts/database/" // VariableSuffix(2:4) // &
		".??????????.dtb"
call GetBatch (DatabaseFile,Batch,Silent=1,ReturnUnalloc=1)
 
if (.not.associated(Batch)) then
  print*, "  > There are no prior accessions for this variable."
  DatabaseFile = "none"
else
  NBatch = size(Batch,1)
  DatabaseFile = Batch(NBatch)
  deallocate (Batch,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Specify: Dealloc #####"
end if

if (ExtraTerm.EQ.1) call system ('printf "' // trim(VariableSuffix) // &
		'" >> ./../../../scratch//f90toidl.pipe')

if (ExtraTerm.EQ.1) then
 do
  print*, "  > Plot corrected stns ? (0=no,1=yes,2=dump, -1=help) "
  read (*,*,iostat=ReadStatus), QDumpView
  if (QDumpView.EQ.-1) then
    print*, "  > 1. allows interactive selection of plots, stn by stn"
    print*, "  > 2. dumps a set of plots to ./../../../scratch/"
  end if
  if (ReadStatus.LE.0.AND.QDumpView.GE.0.AND.QDumpView.LE.2) exit
 end do
else if (ExtraTerm.EQ.2) then
 QDumpView=2
else
 QDumpView=0
end if

if (QDumpView.GT.0) then
 do
  print*, "  > Enter the adjustment threshold, above which to dump stn. (-1=help)"
  read (*,*,iostat=ReadStatus), PlotThresh
  if (PlotThresh.EQ.-1) then
    print*, "  > The adj.thr. is RMSE of changes, in st.dev. of original."
    print*, "  > To dump all, set to zero; to plot most-adjusted, set >1."
  end if
  if (ReadStatus.LE.0.AND.PlotThresh.NE.-1) exit
 end do
end if

end subroutine SpecifyCommon

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

subroutine LoadDatabaseDummy

NDYear=1 ; NDStn=1
allocate (DInfo(NDStn,8), DData(NDYear,NMonth,NDStn), DYearAD(NDYear), &
	  DLocal(NDStn), DName(NDStn), DCty(NDStn), DCode(NDStn), &
	  DLat(NDStn), DLon(NDStn), DElv(NDStn), DNormals(NMonth,NDStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### LoadDatabaseDummy: alloc #####"
DData=DataMissVal ; DNormals=DataMissVal ; DYearAD=1990
DInfo=MissVal ; DCode=MissVal ; DLat=MissVal ; DLon=MissVal ; DElv=MissVal
DLocal="   nocode" ; DName="UNKNOWN" ; DCty="UNKNOWN"

end subroutine LoadDatabaseDummy

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

subroutine LoadDatabase

SubBeg=index(DatabaseFile,"/",.TRUE.)
FileDate=DatabaseFile((SubBeg+05):(SubBeg+10)) 
FileTime=DatabaseFile((SubBeg+11):(SubBeg+14))

call LoadCTS (DInfo,DLocal,DName,DCty,Data=DData,YearAD=DYearAD, &
		Lat=DLat,Lon=DLon,Elv=DElv,Code=DCode, &
		DtbNormals=DNormals,Silent=1,CallFile=DatabaseFile)
NDYear = size(DYearAD,1) ; NDStn = size(DData,3)
print "(a9,2(a2,a1),a2,a3,a2,a1,a2,a,i7,2(a,i4))", "   > Dtb ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," - ", &
		FileTime(1:2),":",FileTime(3:4),": ",NDStn," stns ", &
		DYearAD(1), "-", DYearAD(NDYear)

end subroutine LoadDatabase

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

subroutine LoadOriginal

call LoadCTS (OInfo,OLocal,OName,OCty,Data=OData,Code=OCode,OldCode=OOldCode,Silent=1, &
	      Lat=OLat,Lon=OLon,Elv=OElv,YearAD=YearAD,CallFile=DataFile, &
	      YearADMin=DYearAD(1),YearADMax=DYearAD(NDYear))
NYear=size(OData,1) ; NOStn=size(OData,3) ; OBeg=MissVal ; OEnd=MissVal

do XYear=1,NYear
  XMonth=0
  do
    XMonth=XMonth+1 ; XOStn=0
    do
      XOStn=XOStn+1
      if (OData(XYear,XMonth,XOStn).NE.DataMissVal) then
        if (OBeg.EQ.MissVal) OBeg=XYear
        OEnd=XYear
      end if
      if (XOStn.EQ.NOStn.OR.OEnd.EQ.XYear) exit
    end do
    if (XMonth.EQ.NMonth.OR.OEnd.EQ.XYear) exit
  end do
end do
NOYear=OEnd-OBeg+1

print "(a,11x,a2,i7,2(a,i4))", "   > New file ", ": ",NOStn, &
		" stns ",YearAD(OBeg), "-", YearAD(OEnd)

SubBeg=index(DataFile,".cts") 
DataSrcFile=DataFile ; DataSrcFile(SubBeg:(SubBeg+3))=".src"
call LoadCTS (TrashInfo,TrashLocal,TrashName,TrashCty, &
		Data=ODataSrc,YearAD=TrashYearAD,Silent=1,CallFile=DataSrcFile, &
		YearADMin=DYearAD(1),YearADMax=DYearAD(NDYear))
deallocate (TrashInfo,TrashLocal,TrashName,TrashCty,TrashYearAD, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### LoadBoth: dealloc #####"

call CommonVecPer (YearAD,DYearAD,Beg,End,DBeg,DEnd)

end subroutine LoadOriginal

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

subroutine LoadParallel

ParallelFile=DataFile ; ParallelFile(QRetry:(QRetry+4))="homog"

call LoadCTS (TrashInfo,TrashLocal,TrashName,TrashCty,Code=PCode, &
		Data=PData,YearAD=PYearAD,Silent=1,CallFile=ParallelFile)
NPYear=size(PData,1) ; NPStn=size(PData,3)
print "(a,9x,a2,i7,2(a,i4))", "   > Homog file ", ": ",NPStn, &
		" stns ",PYearAD(1), "-", PYearAD(NPYear)

call CommonVecPer (YearAD,PYearAD,KBeg,KEnd,PBeg,PEnd)

if (PBeg.NE.MissVal) then	! common period between retry (O) and homog (P)
  XOStn=1 ; XPStn=1
  do				! search by stn
    if      (Code(XOStn) .EQ.MissVal) then
      XOStn=XOStn+1
    else if (PCode(XPStn).EQ.MissVal) then
      XPStn=XPStn+1
    else if (Code(XOStn).EQ.PCode(XPStn)) then		! same stn in O and P
      XYear=KBeg-1 ; ThisBeg => EarlBeg ; ThisEnd => EarlEnd
      do					! search by year
        XYear= XYear+1 ; XPYear= XYear-KBeg+PBeg 
	XMonth=0 ; Match=.FALSE. ; Clash=.FALSE.
	
	do					! search by month
	  XMonth=XMonth+1
	  					! year is present in O and P
	  if 	   ( Data( XYear,XMonth,XOStn).NE.DataMissVal.AND. &
	      	    PData(XPYear,XMonth,XPStn).NE.DataMissVal) then
	    Match=.TRUE.
	  else if (( Data( XYear,XMonth,XOStn).NE.DataMissVal.AND. &
	      	    PData(XPYear,XMonth,XPStn).EQ.DataMissVal).OR. &
		   ( Data( XYear,XMonth,XOStn).EQ.DataMissVal.AND. &
	      	    PData(XPYear,XMonth,XPStn).NE.DataMissVal)) then
	    Clash=.TRUE.
	  end if
	  
	  if (XMonth.EQ.NMonth.OR.Match.OR.Clash) exit
	end do
	
	if      (Match) then			! found match per
	  if      (ThisBeg(XOStn).EQ.MissVal) then	! unrecorded, so start
	    ThisBeg(XOStn)= XYear ; ThisEnd(XOStn)= XYear
	  else						! started
	    ThisEnd(XOStn)= XYear
	  end if
	else if (Clash) then			! found clash per
	  if      (ThisBeg(XOStn).EQ.MissVal) then	! mere continuation
	  	! nothing required
	  else						! indicates end of per
	    ThisBeg => LateBeg ; ThisEnd => LateEnd
	  end if
	end if
	
	if (XYear.EQ.KEnd) exit
      end do
      
      if (EarlBeg(XOStn).NE.MissVal.AND.LateBeg(XOStn).EQ.MissVal) then
        LateBeg(XOStn)=EarlBeg(XOStn) ; LateEnd(XOStn)=EarlEnd(XOStn)
      end if
      
      XOStn=XOStn+1 ; XPStn=XPStn+1
    else if (Code(XOStn).GT.PCode(XPStn)) then		! this P stn is not in O
      XPStn=XPStn+1
    else						! this O stn is not in P
      XOStn=XOStn+1
    end if
    
    if (XOStn.GT.NOStn.OR.XPStn.GT.NPStn) exit
  end do
end if

deallocate (TrashInfo,TrashLocal,TrashName,TrashCty,PCode,PData,PYearAD, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### LoadParallel: dealloc #####"

end subroutine LoadParallel

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

subroutine GetCoverage

NLat=180 ; NLon=360 ; BoxLatDeg=1.0 ; BoxLonDeg=1.0

allocate (Got(NLat,NLon),Need(NLat,NLon),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### GetCoverage: alloc #####"
Got=.FALSE. ; Need=.FALSE.

do XOStn=1,NOStn
  if (OLat(XOStn).NE.MissVal.AND.OLon(XOStn).NE.MissVal) then
    Enn=Enn+1
    XLat=floor(( 90.0+OLat(XOStn))/BoxLatDeg)+1 ; if (XLat.GT.NLat) XLat=NLat
    XLon=floor((180.0+OLon(XOStn))/BoxLonDeg)+1 ; if (XLon.GT.NLon) XLon=NLon
    Got(XLat,XLon)=.TRUE.
  end if
end do

call IdentifyRelevant (Got,Need,Decay)	! gridops.f90

NGot=count(Got) ; NNeed=count(Need)
print "(a,2i7)", "   > 1o cells: +stn, need: ",NGot,NNeed

end subroutine GetCoverage

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

subroutine Screen

do XDStn=1,NDStn
  if (DLat(XDStn).EQ.MissVal.OR.DLon(XDStn).EQ.MissVal) then
    QUseDStn = 0
  else
    XLat=floor(( 90.0+DLat(XDStn))/BoxLatDeg)+1 ; if (XLat.GT.NLat) XLat=NLat
    XLon=floor((180.0+DLon(XDStn))/BoxLonDeg)+1 ; if (XLon.GT.NLon) XLon=NLon
    
    if (Need(XLat,XLon)) then
      QUseDStn = 1
    else
      QUseDStn = 0
    end if
  end if
  
  if (QUseDStn.EQ.1) then
    if (QCountOnly.EQ.1) then
      NStn=NStn+1
    else
      XStn=XStn+1
      Local(XStn)=DLocal(XDStn) ; Name(XStn)=DName(XDStn) ; Cty(XStn)=DCty(XDStn)
      Lat(XStn)=DLat(XDStn) ; Lon(XStn)=DLon(XDStn) ; Elv(XStn)=DElv(XDStn)
  
      do XDYear=DBeg,DEnd
        XYear=Beg+(XDYear-DBeg)			
        do XMonth=1,NMonth
          Data(XYear,XMonth,XStn) = DData(XDYear,XMonth,XDStn)
        end do
      end do
    end if
  end if
end do

if (QCountOnly.EQ.0) then
  print "(a,2x,a2,i7,a)", "   > Restrict database ", ": ",(NStn-NOStn)," stns "
  deallocate (DInfo,DData,DLocal,DName,DCty,DCode,DLat,DLon,DElv, &
  	      Got,Need,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### Screen: dealloc #####"
end if

end subroutine Screen

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

subroutine SpecComb

allocate (Data(NYear,NMonth,NStn), DataSrc(NYear,NMonth,NStn), &
	  Local(NStn), Name(NStn), Cty(NStn), Code(NStn), &
	  Lat(NStn), Lon(NStn), Elv(NStn), OldCode(NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SpecComb: alloc #####"

Data=DataMissVal ; DataSrc=DataMissVal
Local="   nocode" ; Name="UNKNOWN" ; Cty="UNKNOWN"
Code=MissVal ; OldCode=MissVal ; Lat=MissVal ; Lon=MissVal ; Elv=MissVal

do XOStn=1,NOStn
    Code(XOStn)=OCode(XOStn) ; OldCode(XOStn)=OOldCode(XOStn) ; Lat(XOStn)=OLat(XOStn)
    Lon(XOStn)=OLon(XOStn) ; Elv(XOStn)=OElv(XOStn)
    Local(XOStn)=OLocal(XOStn) ; Name(XOstn)=OName(XOStn) ; Cty(XOStn)=OCty(XOStn)
    do XYear=1,NYear
      do XMonth=1,NMonth
        Data   (XYear,XMonth,XOStn)=OData   (XYear,XMonth,XOStn)
        DataSrc(XYear,XMonth,XOStn)=ODataSrc(XYear,XMonth,XOStn)
      end do
    end do
end do

deallocate (OInfo,OLocal,OName,OCty,OCode,OOldCode,OLat,OLon,OElv,ODataSrc, &
		stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SpecComb: dealloc 1 #####"

end subroutine SpecComb

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

subroutine MakeArrays

allocate (Disc(NYear,NStn),DiscCrude(NYear,NStn),Trust(NYear,NStn), &
	  Believe(NYear,NStn),RefTS(NYear,NMonth,NStn),Examine(NStn), &
	  Sought(NStn),GotRef(NStn),StnLen(NStn), &
	  EarlBeg(NStn),EarlEnd(NStn),LateBeg(NStn),LateEnd(NStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### MakeArrays: alloc #####"
Disc=.FALSE. ; DiscCrude=.FALSE. ; Trust=.TRUE. ; Believe=.FALSE.
Verbose=.TRUE. ; PlotAllManip=.FALSE. ; RefTS=DataMissVal
EarlBeg=MissVal ; EarlEnd=MissVal ; LateBeg=MissVal ; LateEnd=MissVal

Examine=.FALSE.
do XStn=1,NOStn
    Examine(XStn)=.TRUE.
    
    BegO=0 ; EndO=0
    do XYear=1,NYear
      do XMonth=1,NMonth
        if (Data(XYear,XMonth,XStn).NE.DataMissVal) then
	  if (BegO.EQ.0) BegO=XYear
	  EndO=XYear
	end if
      end do
    end do
    StnLen(XStn)=EndO-BegO+1
end do
Sought=Examine

end subroutine MakeArrays

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

subroutine SaveTwo

allocate (AdjSrc(NYear,NMonth,NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SaveTwo: alloc: AdjSrc #####"
AdjSrc=DataMissVal ; LogAdj=0 ; LogNon=0

do XYear=1,NYear
  do XMonth=1,NMonth
    do XStn=1,NStn
      if (Examine(XStn)) then
        if (Adj(Xyear,XMonth,XStn).NE.DataMissVal) then
      	  AdjSrc(Xyear,XMonth,XStn) = DataSrc(Xyear,XMonth,XStn)
      	  LogAdj=LogAdj+1
        else if (Data(Xyear,XMonth,XStn).NE.DataMissVal) then
          LogNon=LogNon+1
        end if
      end if
    end do
  end do
end do 
print "(a,2i7)", "   > Data to homog, retry: ",LogAdj,LogNon
write (99,"(a,2i7)"), "   > Data to homog, retry: ",LogAdj,LogNon

call MakeStnInfo (Info,Code,OldCode,Lat,Lon,Elv,1, &
		  YearAD=YearAD,Data=Adj,Silent=1,QBlank=QBlank)
if (QBlank) then
  print*, "  > It has not been possible to check any data."
else
  print "(2a)", "   > Save file: ", trim(SaveFileTwo)
  SubBeg=index(SaveFileTwo,".cts")  ; SaveSrcFile=SaveFileTwo(1:(SubBeg-1)) // ".src"
  call SaveCTS (Info,Local,Name,Cty,Data=Adj,YearAD=YearAD, &
		CallFile=SaveFileTwo,Silent=1)
  call SaveCTS (Info,Local,Name,Cty,Data=AdjSrc,YearAD=YearAD, &
		CallFile=SaveSrcFile,Silent=1)
end if

call DumpAdjHdr
call DumpAdjPer

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

do XStn=1,NStn
  if (.NOT.Examine(XStn)) Code(XStn)=MissVal	! remove dtb stns from retry
  
  do XYear=1,NYear
    do XMonth=1,NMonth
      if (Data(Xyear,XMonth,XStn).EQ.DataMissVal) &
      	  DataSrc(Xyear,XMonth,XStn) = DataMissVal
    end do
  end do
end do 

call MakeStnInfo (Info,Code,OldCode,Lat,Lon,Elv,1, &
		  YearAD=YearAD,Data=Data,Silent=1,QBlank=QBlank)
if (QBlank) then
	! nothing here
else
  SubBeg=index(SaveFileTwo,"homog") ; SaveFileTwo(SubBeg:(SubBeg+4))="retry"
  SubBeg=index(SaveFileTwo,".cts")  ; SaveSrcFile=SaveFileTwo(1:(SubBeg-1)) // ".src"
  call SaveCTS (Info,Local,Name,Cty,Data=Data,YearAD=YearAD, &
		CallFile=SaveFileTwo,Silent=1)
  call SaveCTS (Info,Local,Name,Cty,Data=DataSrc,YearAD=YearAD, &
		CallFile=SaveSrcFile,Silent=1)
end if

deallocate (AdjSrc,DataSrc,Info, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### SaveTwo: dealloc: *Src #####"

end subroutine SaveTwo

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

subroutine DumpAdjHdr

SubBeg=index(SaveFileTwo,".cts") ; HdrFile=SaveFileTwo(1:(SubBeg-1)) // ".adj.hdr"

call QuickSort (Reals=AdjStn,Order=Order,NMiss=NMiss)

CallForm="(f7.2,13x)"
do XStn=1,NStn
    String20 = GetTextFromReal (AdjStn(XStn),CallForm=CallForm)
    Local(XStn) = adjustr(String20(1:9))
end do

call SaveCTS (Info,Local,Name,Cty,CallFile=HdrFile,HeadOnly=1,Silent=1,Order=Order)
		    
deallocate (Order, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### DumpAdjHdr: dealloc #####"

end subroutine DumpAdjHdr

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

subroutine DumpAdjPer

SubBeg=index(SaveFileTwo,".cts") ; PerFile=SaveFileTwo(1:(SubBeg-1)) // ".adj.per"

allocate (AdjSea(NYear,4),AdjAnn(NYear), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### DumpAdjPer: alloc #####"

call FillSeaAnnMean (YearAD,AdjMon,AdjSea,AdjAnn)
call SavePER (PerFile,YearAD,0,Monthly=AdjMon,Seasonal=AdjSea,Annual=AdjAnn,NoResponse=1)

deallocate (AdjSea,AdjAnn, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### DumpAdjPer: dealloc #####"

end subroutine DumpAdjPer

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

subroutine PlotDiscStn

XIterate=1

do XStn = 1, NStn
  if (Examine(XStn).AND.AdjStn(XStn).GE.PlotThresh) then
    XYear=0 ; QAbs=0 ; QAnom=0 ; QStep=0 ; QDump=99
    do
     XYear=XYear+1
     
     if (Disc(XYear,XStn).OR.PlotAllManip) then
      call CalcClimatology
      QDump=7 ; call DumpForIDL ; QStep=1
      QDump=8 ; call DumpForIDL ; QDump=9 ; call DumpForIDL
      
      do
        do
          print "(5a)", "   > Plot ", trim(Name(XStn)), " (", trim(Cty(XStn)), &
          		") (0=exit,99=help)"
          read (*,*,iostat=ReadStatus), QDump
          if (ReadStatus.LE.0) exit
        end do
        
        if (QDump.GE.1.AND.QDump.LE.9) then
          call DumpForIDL

          if      (QDump.LT.4) then
            QAbs=1
          else if (QDump.LT.7) then
            QAnom=1
          else
            QStep=1
          end if
        else if (QDump.EQ.99) then
          print*, "  > ABS:  1=orig,2=refTS,3=okay"
          print*, "  > ANOM: 4=orig,5=refTS,6=okay"
          print*, "  > DIFF: 7=1-2, 8=1-3,  9=3-2"
        end if
        
        if (QDump.EQ.0) exit
      end do
     end if
     
     if (XYear.EQ.NYear.OR.QDump.EQ.0) exit
    end do
  end if 
end do

end subroutine PlotDiscStn

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

subroutine DumpDiscStn

XIterate=1 ; write (99,*)
write (99,"(a)"), "   PLOT      STAT"

do XStn = 1, NStn
  if (Examine(XStn).AND.AdjStn(XStn).GT.PlotThresh) then
    QAbs=0 ; QAnom=0 ; QStep=0 ; QDump=99
    write (99,"(i7,f10.2)"), Code(XStn), AdjStn(XStn)
    call CalcClimatology
      
    do QDump = 1, 9
          call DumpForIDL

          if      (QDump.LT.4) then
            QAbs=1
          else if (QDump.LT.7) then
            QAnom=1
          else
            QStep=1
          end if
          
	  if (ExtraTerm.EQ.1) then
            SubBeg = index(DumpFile,".txt")
            DumpFile = DumpFile(1:SubBeg) // "eps"
            do
              inquire (file=DumpFile, exist=Plotted)
              if (Plotted) exit
            end do
	  end if
    end do
  end if 
end do

end subroutine DumpDiscStn

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

subroutine CalcClimatology

Climatology=0
do XMonth = 1, NMonth
  OpTot=0 ; OpEn=0
  do XYear = 1, NYear
    if (RefTS(XYear,XMonth,XStn).NE.DataMissVal) then
      OpEn=OpEn+1 ; OpTot=OpTot+real(RefTS(XYear,XMonth,XStn))
    end if
  end do
  if (OpEn.Gt.0) then
    Climatology(XMonth)=nint(OpTot/OpEn)
    if ((.NOT.Differ).AND.OpTot.EQ.0.0) Climatology(XMonth)=0.1
  end if
end do

end subroutine CalcClimatology

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

subroutine DumpForIDL

DumpText(5:8) = "Orig"

if      (QDump.EQ.1) then
	FileText = "a.oriabs" ; DumpText(1:4) = "Abso" ; DumpTarget => OData
	Title = "Original" ; if (QAbs.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.2) then
	FileText = "b.refabs" ; DumpText(1:4) = "Abso" ; DumpTarget => RefTS
	Title = "Reference" ; if (QAbs.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.3) then
	FileText = "c.yesabs" ; DumpText(1:4) = "Abso" ; DumpTarget => Adj
	Title = "Revised" ; if (QAbs.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.4) then
	FileText = "d.oriano" ; DumpText(1:4) = "Anom" ; DumpTarget => OData
	if (Differ) then
	  Title = "Original (anomalies)" 
	else
	  Title = "Original : normal (ratio)" 
	end if
	if (QAnom.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.5) then
	FileText = "e.refano" ; DumpText(1:4) = "Anom" ; DumpTarget => RefTS
	if (Differ) then
	  Title = "Reference (anomalies)" 
	else
	  Title = "Reference : normal (ratio)"
	end if
	if (QAnom.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.6) then
	FileText = "f.yesano" ; DumpText(1:4) = "Anom" ; DumpTarget => Adj
	if (Differ) then
	  Title = "Revised (anomalies)" 
	else
	  Title = "Revised : normal (ratio)"
	end if
	if (QAnom.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.7) then
	FileText = "g.shifts" ; DumpText(1:4) = "Step" ; DumpTarget => OData
	if (Differ) then
	  Title = "Original minus Reference"
	else
	  Title = "Original : Reference (ratio)"
	end if
	if (QStep.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.8) then
	FileText = "h.change" ; DumpText(1:4) = "Step" ; DumpTarget => OData
	if (Differ) then
	  Title = "Original minus Revised" 
	else
	  Title = "Original : Revised (ratio)" 
	end if
	if (QStep.EQ.1) DumpText(5:8) = "Revi"
else if (QDump.EQ.9) then
	FileText = "i.cleand" ; DumpText(1:4) = "Step" ; DumpTarget => RefTS
	if (Differ) then
	  Title = "Revised minus Reference" 
	else
	  Title = "Revised : Reference (ratio)" 
	end if
	if (QStep.EQ.1) DumpText(5:8) = "Revi"
end if

StnText  = GetTextFromInt(Code(XStn))
IterText = GetTextFromInt(XIterate)
DumpFile = './../../../scratch/' // trim(adjustl(StnText)) // '.' // FileText // '.txt'

open (80,file=DumpFile,status="replace",action="write")
write (80,"(i8,x,a20,x,a13,x,i4,x,a8,x,i2)"), Code(XStn),Name(XStn), &
			Cty(XStn),NYear,DumpText,QDumpView
write (80,"(a80)"), Title
do XYear = 1,NYear
    DumpData=DataMissVal
    
    if (Differ) then
      call GrabMeans
    else
      call GrabTotals
    end if
    
    do XMonth=1,NMonth
      if (DumpData(XMonth).GT.99999.OR.DumpData(XMonth).LT.-9999) &
      		 DumpData(XMonth) = DataMissVal
    end do
    
    write (80,"(i4,12i5,f10.2)"), YearAD(XYear), &
    	(DumpData(XMonth), XMonth=1,NMonth), OpMean
end do
close (80)
  
if (ExtraTerm.EQ.1) call system &
	('printf "' // trim(DumpFile) // '" >> ./../../../scratch/f90toidl.pipe')

end subroutine DumpForIDL

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

subroutine GrabMeans

OpTot=0 ; OpEn=0 ; OpMean=MissVal ; OpValid=0			

do XMonth=1,NMonth
      if (DumpTarget(XYear,XMonth,XStn).NE.DataMissVal) then
      	if      (QDump.LT.4) then
	  DumpData(XMonth) = DumpTarget(XYear,XMonth,XStn)
        else if (QDump.LT.7) then
	  DumpData(XMonth) = DumpTarget(XYear,XMonth,XStn)-Climatology(XMonth)
	else if (QDump.EQ.7) then
          if (RefTS(XYear,XMonth,XStn).NE.DataMissVal) &
              DumpData(XMonth) = OData(XYear,XMonth,XStn)-RefTS(XYear,XMonth,XStn)
        else if (QDump.EQ.8) then
          if (Adj(XYear,XMonth,XStn).NE.DataMissVal) &
              DumpData(XMonth) = OData(XYear,XMonth,XStn)-Adj(XYear,XMonth,XStn)
	else if (QDump.EQ.9) then
          if (Adj(XYear,XMonth,XStn).NE.DataMissVal) &
              DumpData(XMonth) = Adj(XYear,XMonth,XStn)-RefTS(XYear,XMonth,XStn)
	end if
      end if
      
      if (DumpData(XMonth).NE.DataMissVal) then		! month = valid
        OpTot=OpTot+DumpData(XMonth) ; OpEn=OpEn+1 ; OpValid=OpValid+1
      else if (QDump.LT.7) then				! est month from yr-5, yr-1
        OpAdd=0 ; OpWin=0
        do XWindow=XYear-5,XYear-1
          if (XWindow.GT.0.AND.XWindow.LE.NYear) then
            if (DumpTarget(XWindow,XMonth,XStn).NE.DataMissVal) then
              OpAdd=OpAdd+DumpTarget(XWindow,XMonth,XStn) ; OpWin=OpWin+1
            end if
          end if
        end do
        if (OpWin.GT.1) then
          if (QDump.LT.4) then
            OpTot=OpTot+(OpAdd/OpWin) ; OpEn=OpEn+1
          else
            OpTot=OpTot+((OpAdd/OpWin)-Climatology(XMonth)) 
	    OpEn=OpEn+1
          end if
        end if
      else						! est month from mon-2,mon+2
        OpAdd=0 ; OpWin=0 ; XWindow=XMonth-3
        do
          XWindow=XWindow+1
          if (Xwindow.GT.0.AND.Xwindow.LE.NMonth) then
            if (DumpData(XWindow).NE.DataMissVal) then
              OpAdd=OpAdd+DumpData(XWindow) ; OpWin=OpWin+1
            end if
          end if
          if (XWindow.EQ.(XMonth+2).OR.XWindow.EQ.12) exit
        end do
        if (OpWin.GT.0) then
            OpTot=OpTot+(OpAdd/OpWin) ; OpEn=OpEn+1
        end if
      end if
end do
    
if (OpValid.GT.6) OpMean=OpTot/opEn
    
end subroutine GrabMeans

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

subroutine GrabTotals

OpTot=0 ; OpEn=0 ; OpMean=MissVal ; OpValid=0 ; OpClim=0

do XMonth=1,NMonth
        RealTarget = real(DumpTarget(XYear,XMonth,XStn))
        RealClim   = real(Climatology(XMonth))
        RealOData  = real(OData(XYear,XMonth,XStn))
        RealRefTS  = real(RefTS(XYear,XMonth,XStn))
        RealAdj    = real(Adj(XYear,XMonth,XStn))
        
      	if      (QDump.LT.4) then			! fill monthly data array
      	  if (RealTarget.NE.DataMissVal) then
	    DumpData(XMonth) = DumpTarget(XYear,XMonth,XStn)
	    OpTot=OpTot+RealTarget 
	    OpValid=OpValid+1
	  else
	    OpTot=OpTot+RealClim
	  end if
        else if (QDump.LT.7) then
	  if (RealTarget.NE.DataMissVal) then
	    DumpData(XMonth) = nint(100.0*RealTarget/RealClim)
	    OpTot=OpTot+RealTarget 
	    OpValid=OpValid+1
	  else
	    OpTot=OpTot+RealClim
	  end if
	else if (QDump.EQ.7) then
          if (RealOData.NE.DataMissVal.AND.RealRefTS.NE.DataMissVal.AND. &
	  		RealRefTS.NE.0) then
	    DumpData(XMonth) = nint(100.0*RealOData/RealRefTS)
	    OpTot=OpTot+((RealOData/RealRefTS)*RealClim)
	    OpValid=OpValid+1
	  else
	    OpTot=OpTot+RealClim
	  end if
        else if (QDump.EQ.8) then
          if (RealOData.NE.DataMissVal.AND.RealAdj.NE.DataMissVal.AND. &
	  		RealAdj.NE.0) then
	    DumpData(XMonth) = nint(100.0*RealOData/RealAdj)
	    OpTot=OpTot+((RealOData/RealAdj)*RealClim)
	    OpValid=OpValid+1
	  else
	    OpTot=OpTot+RealClim
	  end if
	else if (QDump.EQ.9) then
          if (RealRefTS.NE.DataMissVal.AND.RealAdj.NE.DataMissVal.AND. &
	  		RealRefTS.NE.0) then
	    DumpData(XMonth) = nint(100.0*RealAdj/RealRefTS)
	    OpTot=OpTot+((RealAdj/RealRefTS)*RealClim)
	    OpValid=OpValid+1
	  else
	    OpTot=OpTot+RealClim
	  end if
	end if
	
	OpClim=OpClim+RealClim
end do
    
if (OpValid.GT.6) then
      if      (QDump.LT.4) then
        OpMean=OpTot
      else
        OpMean=nint(100.0*OpTot/OpClim)
      end if
end if
    
end subroutine GrabTotals

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

end program Homogiter