! updatedtb.f90
! f90 program written by Tim Mitchell on 14.05.02 as accession.f90
! rewritten as updatedtb.f90
! last modified on 25.03.03
! program to process new files of station monthly obs
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/updatedtb time.f90 filenames.f90 grimfiles.f90 crutsfiles.f90 
!	gridops.f90 sortmod.f90 plainperm.f90 regress.f90 ctyfiles.f90 countries.f90 
!	ghcnrefiter.f90 ./../cruts/updatedtb.f90 2> /tyn1/tim/scratch/stderr.txt

program UpdateDTB

use Time
use FileNames
use GrimFiles
use CRUtsFiles
use GridOps
use SortMod
use PlainPerm
use Regress
use CtyFiles
use Countries
use GHCNrefIter

implicit none

logical, pointer, dimension (:,:)	:: Trust
logical, pointer, dimension (:)		:: Normalise,GotRef,GetSingle,CandValid,PerTrust,Calc

real, pointer, dimension (:,:,:)	:: Par
real, pointer, dimension (:,:)		:: RelA,RelB,RelR, PoolWei
real, pointer, dimension (:)		:: LatO,LonO,ElvO, LatM,LonM,ElvM
real, pointer, dimension (:)		:: LatA,LonA,ElvA, LatD,LonD,ElvD
real, pointer, dimension (:)		:: MasterLat,MasterLon,MasterEn
real, pointer, dimension (:)		:: OriginalLat,OriginalLon,OriginalEn
real, pointer, dimension (:)		:: CtyLon,CtyLat,CtyTot
real, pointer, dimension (:)		:: MatchCodeA,MatchDistanceA,MatchNameA,MatchScoreA
real, pointer, dimension (:)		:: CandidScoreA,CandidScoreM
real, pointer, dimension (:)		:: Decay, VecO,VecD,VecR, Exe,Wye

integer, pointer, dimension (:,:,:)	:: DataO,DataSrcO,DataD,DataSrcD, RefTS
integer, pointer, dimension (:,:,:) 	:: FD,AbsD,ParPer
integer, pointer, dimension (:,:)	:: DataNmlO,DataNmlD,DataNmlSrcD
integer, pointer, dimension (:,:)	:: Info,TrashInfo, StnPer,Neigh,CandFD,CandAbs,Pool
integer, pointer, dimension (:,:)	:: Combo
integer, pointer, dimension (:)		:: StnOldO,StnNewM,StnOldA,StnNewA,StnOldD,StnNewD
integer, pointer, dimension (:)         :: StnSrcCodeO,StnSrcCodeA,YearADO,YearADD,TrashYearAD
integer, pointer, dimension (:)		:: StnCodeBegO,StnCodeEndO,StnContinentO
integer, pointer, dimension (:)		:: StnSrcDateA,MissInts, Per2Stn,Stn2Per,PerBeg,PerEnd
integer, pointer, dimension (:) 	:: MasterCode0,MasterCode1,MasterContinent
integer, pointer, dimension (:) 	:: CandidIndexM,CandidIndexA
integer, pointer, dimension (:) 	:: OrderA,OrderD,OrderM

character (len=80), pointer, dimension (:) :: Batch
character (len=80), pointer, dimension (:) :: MasterBatch
character (len=20), pointer, dimension (:) :: StnNameO,StnNameM,StnNameA,StnNameD,TrashName
character (len=13), pointer, dimension (:) :: StnCtyO,StnCtyM,StnCtyA,StnCtyD,TrashCty
character (len=13), pointer, dimension (:) :: MasterCty,MasterRawName,MasterFinalName
character (len=09), pointer, dimension (:) :: StnLocalO,StnLocalM,StnLocalA,StnLocalD,TrashLocal
character (len=04), pointer, dimension (:) :: StnSrcSuffixA

logical, dimension (12)	:: CheckMon
real, dimension (12)	:: MonthA,MonthB,MonthR,SumR
integer, dimension (12)	:: Min,Max,MonthPrior,MonthFresh,Adjust

real, parameter 	:: MissVal 	= -999.0
real, parameter 	:: SumRThresh 	= 0.25
integer, parameter 	:: DataMissVal 	= -9999
integer, parameter 	:: MultiSource 	= -888

character (len=80), parameter :: LatLonFile = './../../../scratch/test-latlon.txt'

logical :: ZOneDP,Differ,Success,Okay,Dense,Extrer,Overlapper,Butter,Adder,Pseudo

real :: OpVal,OpTotSq,OpTot,OpEn,OpZero,OpStDev,Aye,Bee,Are,Enn,ErrSum,ErrEnn,RMSE
real :: Factor,Distance,StnODigit6,RealVal,CorrelMean,WeightMean,DecayMemory,Estimate
real :: TotNull,TotNeg,TotSmall,TotLarge,TotPos,TotMon,MeanR,Lat1DP,Lon1DP,Dist1DP,Dist2DP
real :: PriorTot,PriorEn,FreshTot,FreshEn,NuffFresh,NuffPrior,SmallThresh,LargeThresh

integer :: ReadStatus, AllocStat
integer :: QOrigCode,QOrigSrc,QDigits,QIncludeOrig,QConvertedCty,QDecay
integer :: QDodgy,QAlter,QProceed,QRelyMeta
integer :: NOYear,NDYear,NOStn,NMStn,NAStn,NDStn,NODYear
integer :: XOYear,XDYear,XOStn,XMStn,XAStn,XDStn,XODYear
integer :: NChar,NBatch,NMonth,NMasterCty,NPrior,NCandid,NDecay,NRel,NAttempt
integer :: XChar,XBatch,XMonth,XMasterCty,XPrior,XCandid,XDecay,XRel,XAttempt
integer :: SubBeg,SubLen,Ascii,Code,OrigSrcCode,AnswerNoYes,Reverse,Find,Seek
integer :: Butt,ButtYX,DecayDistance,StampInt,Window,Increment,NeedThresh
integer :: CandidAStn,CandidDStn,CandidMStn
integer :: AssignAStn,AssignDStn,AssignMStn, Log1,Log2,Log3,Log4, Elastic
integer :: MergeBeg,MergeEnd,Beg,End,ClimBeg,ClimEnd,OStart,OEnd,DStart,DEnd

character (len=80) :: GivenFile,Variable
character (len=80) :: OrigFileL,OrigSrcFileL,DatabaseFileL,MasterFileL,AccessFileL
character (len=80) :: OrigFileS,OrigSrcFileS,DatabaseFileS,MasterFileS,AccessFileS
character (len=80) :: DatabaseSrcFileL
character (len=80) :: DatabaseSrcFileS
character (len=20) :: Dirty,Name,Stamp20
character (len=10) :: CurrentTime,FileTime,Stamp10,LogText
character (len=08) :: CurrentDate,FileDate
character (len=04) :: VariableSuffix

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

call Intro			
call SelectOriginal
call SelectVariable
call FindDatabase

call LoadingOriginal
call LoadingMaster
call LoadingAccessions
call ProcessOriginal
call IdentifyStnSrc		

call SetupAssign
call Assignment					! main loop
call ShutdownAssign

call PrepareSave
call SavingDatabase
call SavingAccession
call SavingMaster

call Conclude

!*******************************************************************************
! Assignment
! AssignStnMeta
! CheckDataDensity	checks data>2 for each month, else Dense=FALSE
! CheckOverlap
! Conclude 		conclude program
! ConcludeAssign
! FindDatabase		identify database file
! GetAdditions		identify additions to existing database from this stn 
! GetButtJoin		uses ref ts to ensure butt-join fits well
! IdentifyStnSrc	identifies a single src code for each stn
! Intro			initialises program
! LoadingAccessions	loading accessions file
! LoadingDatabase	loading database files (.dtb and .dts)
! LoadingMaster		load master .hdr file
! LoadingOriginal	load original .cts file
! LocateStn
! LogMarker		put progress marker in log file
! MakeNewDatabase	make new database filled with missing values
! MatchCandid 		identify master candidates from priors
! MatchCodeAccess	match code with accessions file
! MatchCodeMaster	match code with master file
! MatchDistanceAccess	match location with accessions file
! MatchNameAccess	match stn name with accessions file
! MatchScoreAccess	calc match score for each access file stn
! MatchScoreMaster	obtains final master score for candidate
! MergeOverlap
! MergeExisting
! NeededStretch
! PrepareAssign
! PrepareSave		prepares for saving
! ProcessOriginal
! SavingAccession
! SavingDatabase
! SavingMaster
! SelectOriginal	provides user with load options and specifications
! SelectProc		processing options
! SelectVariable	select climate variable, by suffix
! SetupAssign		prepares arrays for assignment operations
! ShutDownAssign

contains

!*******************************************************************************
! carries out the assignment (iterating by stn, main loop)

subroutine Assignment

do XOStn = 1, NOStn		! assign signed 6-digit stn code
 call CheckDataDensity			! checks data>2 for each month, else Dense=FALSE
 
 if (Dense) then
  call PrepareAssign
  
  call MatchCodeAccess				! match accession entries
  call MatchDistanceAccess
  call MatchNameAccess
  call MatchScoreAccess
  call MatchCandid				! identify master candidates
  
  if (NCandid.GT.0) then
    NAttempt=NCandid ; OrderM=MissVal
    
    do XCandid=1,NCandid			! generates the metadata score for each cand
      CandidAStn = CandidIndexA(XCandid)	!	stored in CandidScoreM(XCandid)
      CandidMStn = CandidIndexM(XCandid)
      
      call MatchCodeMaster
      call MatchScoreMaster
      						! for qualifying matches, seek dtb stn
      if (CandidScoreM(XCandid).LT.0.OR.CandidScoreA(XCandid).LT.0) then
        if (LocateStn(StnNewM(CandidMStn),StnNewD).NE.MissVal) &
		CandidScoreM(XCandid)=CandidScoreM(XCandid)-100
      else
        NAttempt=NAttempt-1 ; CandidScoreM(XCandid)=MissVal
      end if
    end do
    
    if (NAttempt.GT.0) then
     XAttempt=0
     do						! attempt a match, by candidate
      XAttempt=XAttempt+1 ; OpTot=100000 ; OpEn=MissVal
      do XCandid=1,NCandid
        if (CandidScoreM(XCandid).LT.OpTot .AND. &
	    CandidScoreM(XCandid).NE.MissVal) OpEn=XCandid
      end do
      if (OpEn.NE.MissVal) then
        XCandid=OpEn
        CandidAStn = CandidIndexA(XCandid)	
        CandidMStn = CandidIndexM(XCandid)
        CandidDStn = LocateStn (StnNewM(CandidMStn),StnNewD)
	CandidScoreM(XCandid) = MissVal
      end if
      
      if (CandidDStn.NE.MissVal) then	! dtb stn exists
	  call GetAdditions				! look for extra data

	  call CheckNinety				! carry out 90% check
	  if (.NOT.Overlapper) call CheckOverlap	! look for overlap
	  
	  if ((.NOT.Extrer).AND.Butter) then	! if no-extra but good-match
	    Butter=.FALSE. ; Overlapper=.TRUE.		! make overlapping
	  end if
	  
	  if      (Overlapper) then		! accept overlap
  	    call AssignStnMeta				! assign metadata
  	    call MergeOverlap				! merge data
          else if (Butter) then			! attempt butt join
	    AssignDStn=CandidDStn			! assume a good match
	    call MergeExisting				! merge data (not sources)
            call GetButtJoin				! seek butt-join (correct if nec)
          end if
      else				! dtb stn does not exist
      	  Adder=.TRUE.
	  call AssignStnMeta				! assign metadata
          call MergeExisting				! add data
      end if
      
      if (XAttempt.EQ.NAttempt.OR.Overlapper.OR.Butter.OR.Adder) exit
     end do
    end if
  end if
    
  if ((.NOT.Overlapper).AND.(.NOT.Butter).AND.(.NOT.Adder)) then
    Pseudo=.TRUE.
    call AssignStnMeta
    call MergeExisting
  end if
    
  call ConcludeAssign
 end if
end do

end subroutine Assignment

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

subroutine AssignStnMeta
					! ensure M stn index
if (Pseudo) then			! new pseudo
  AssignMStn = LocateStn (nint(MissVal),StnNewM)	! find next free master stn
  call FindNextPseudo					! place pseudo in StnNewM
else					! any other assignment
  AssignMStn = CandidMStn
end if
					! ensure D stn index
if (Pseudo.OR.Adder) then		! new pseudo or dtb entry
  AssignDStn = LocateStn (nint(MissVal),StnNewD)	! find next free dtb stn
else					! any other assignment
  AssignDStn = CandidDStn
end if
					! ensure A stn index
if (Pseudo.OR.Adder) then		! new pseudo or dtb entry
  AssignAStn = LocateStn (nint(MissVal),StnNewA)	! find next free dtb stn
else					! any other assignment
  if (StnOldA(CandidAStn).EQ.StnOldO(XOStn)   .AND. &
      LatA(CandidAStn).EQ.LatO(XOStn)         .AND. &
      LonA(CandidAStn).EQ.LonO(XOStn)         .AND. &
      ElvA(CandidAStn).EQ.ElvO(XOStn)         .AND. &
      StnNameA(CandidAStn).EQ.StnNameO(XOStn) .AND. &
      StnCtyA(CandidAStn).EQ.StnCtyO(XOStn)   .AND. &
      StnSrcCodeA(CandidAStn).EQ.StnSrcCodeO(XOStn)) then
    AssignAStn = CandidAStn
  else
    AssignAStn = LocateStn (nint(MissVal),StnNewA)
  end if
end if

if (LatM(AssignMStn).EQ.MissVal) LatM(AssignMStn) = LatO(XOStn)
if (LonM(AssignMStn).EQ.MissVal) LonM(AssignMStn) = LonO(XOStn)
if (ElvM(AssignMStn).EQ.MissVal) ElvM(AssignMStn) = ElvO(XOStn)
if (trim(StnNameM(AssignMStn)).EQ."UNKNOWN") StnNameM(AssignMStn) = StnNameO(XOStn)
if (trim(StnCtyM(AssignMStn)).EQ."UNKNOWN") StnCtyM(AssignMStn) = StnCtyO(XOStn)
if (StnLocalM(AssignMStn).EQ."   nocode") StnLocalM(AssignMStn) = StnLocalO(XOStn)

if (StnNewD(AssignDStn).EQ.MissVal) StnNewD(AssignDStn) = StnNewM(AssignMStn)
if (LatD(AssignDStn).EQ.MissVal) LatD(AssignDStn) = LatM(AssignMStn)
if (LonD(AssignDStn).EQ.MissVal) LonD(AssignDStn) = LonM(AssignMStn)
if (ElvD(AssignDStn).EQ.MissVal) ElvD(AssignDStn) = ElvM(AssignMStn)
if (trim(StnNameD(AssignDStn)).EQ."UNKNOWN") StnNameD(AssignDStn) = StnNameM(AssignMStn)
if (trim(StnCtyD(AssignDStn)).EQ."UNKNOWN") StnCtyD(AssignDStn) = StnCtyM(AssignMStn)
if (StnLocalD(AssignDStn).EQ."   nocode") StnLocalD(AssignDStn) = StnLocalM(AssignMStn)
if (.NOT.Normalise(AssignDStn)) Normalise(AssignDStn)=.TRUE.	! calc normal

if (StnOldA(AssignAStn).EQ.MissVal) then	! new accession entry required
  StnOldA(AssignAStn) = StnOldO(XOStn) 
  StnNewA(AssignAStn) = StnNewM(AssignMStn)
  LatA(AssignAStn) = LatO(XOStn)
  LonA(AssignAStn) = LonO(XOStn)
  ElvA(AssignAStn) = ElvO(XOStn)
  StnNameA(AssignAStn)  = StnNameO(XOStn)
  StnCtyA(AssignAStn)   = StnCtyO(XOStn)
  StnLocalA(AssignAStn) = StnLocalO(XOStn)
  StnSrcCodeA(AssignAStn)   = StnSrcCodeO(XOStn)
  StnSrcSuffixA(AssignAStn) = VariableSuffix
  StnSrcDateA(AssignAStn)   = StampInt
end if

end subroutine AssignStnMeta

!*******************************************************************************
! checks data>2 for each month, else Dense=FALSE

subroutine CheckDataDensity

Dense=.TRUE. ; XMonth=0
do
  XMonth=XMonth+1 ; XOYear=0 ; OpEn=0
  do
    XOYear=XOYear+1
    if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) OpEn=OpEn+1
    if (XOYear.EQ.NOYear.OR.OpEn.EQ.2) exit
  end do
  if (OpEn.LT.2) Dense=.FALSE.
  if (XMonth.EQ.12.OR.(.NOT.Dense)) exit
end do

end subroutine CheckDataDensity

!*******************************************************************************
! carry out 90% test = if 90% data in overlap are identical,
!   with at least 3 months checked, assume overlap is identical

subroutine CheckNinety

TotNeg=0 ; TotPos=0 ; TotMon=0 ; CheckMon=.FALSE.
do XMonth=1,NMonth
  do XODYear = 1, NODYear
    XDYear = DStart + XODYear - 1 ; XOYear = OStart + XODYear - 1
    if (DataD(XDYear,XMonth,CandidDStn).NE.DataMissVal.AND.&
    	DataO(XOYear,XMonth,XOStn).NE.DataMissVal) then
      if (.NOT.CheckMon(XMonth)) then
        CheckMon(XMonth)=.TRUE. ; TotMon=TotMon+1
      end if
      
      if (DataD(XDYear,XMonth,CandidDStn).EQ.DataO(XOYear,XMonth,XOStn)) then
        TotPos=TotPos+1
      else
        TotNeg=TotNeg+1
      end if
    end if
  end do
end do

if (TotMon.GE.3) then
  if ((TotPos/(TotPos+TotNeg)).GE.0.9) then
    Overlapper=.TRUE. ; Butter=.FALSE.
  end if
end if

end subroutine CheckNinety

!*******************************************************************************
! checks overlap between original and candidate database entry

subroutine CheckOverlap

TotNull=0 ; TotNeg=0 ; TotSmall=0 ; TotLarge=0 ; MeanR=0
do XMonth=1,NMonth
  Exe=MissVal ; Wye=MissVal
  do XODYear = 1, NODYear
    XDYear = DStart + XODYear - 1 ; XOYear = OStart + XODYear - 1
    if (DataD(XDYear,XMonth,CandidDStn).NE.DataMissVal) &
    		Exe(XODYear) = DataD(XDYear,XMonth,CandidDStn)
    if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) &
    		Wye(XODYear) = DataO(XOYear,XMonth,XOStn)
  end do

  if (Differ) then
    call MergeForDiff  (Exe,Wye,1,NODYear)
  else
    call MergeForRatio (Exe,Wye,1,NODYear)
  end if
  
  call RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,4)	! 4 is deliberate
!  write (99,"(a,i3,4f10.2)"), "overlap month",XMonth,Aye,Bee,Are,Enn	! ####################
  
  SmallThresh=0.7-(Enn*0.05) ; if (SmallThresh.LT.0.1) SmallThresh=0.1
  LargeThresh=SmallThresh+0.2
  
  if      (Are.EQ.MissVal.OR.(Are.EQ.0.0.AND.Aye.EQ.0.0)) then
      TotNull=TotNull+1
  else if (Are.LT.SmallThresh) then
      TotNeg=TotNeg+1     ; MeanR=MeanR+Are
  else if (Are.LT.LargeThresh) then
      TotSmall=TotSmall+1 ; MeanR=MeanR+Are
  else
      TotLarge=TotLarge+1 ; MeanR=MeanR+Are
  end if  
end do

if (TotNull.LT.12) MeanR=MeanR/(12-TotNull)

! write (99,"(a,4f5.1,f10.2)"), "overlap total",TotNull,TotNeg,TotSmall,TotLarge,MeanR ! ####################

if      (TotNull.GE.6) then
  Overlapper=.FALSE. ; Butter=.TRUE.	! nulls - butt join
else if ((TotSmall+(2*TotNeg)).GE.TotLarge) then
  Overlapper=.FALSE. ; Butter=.FALSE.	! bad overlap - no match!
  write (99,"(a10,2i8,f7.2,f8.2,i5,x,a20,x,a13)"), "  uncorrel", StnOldO(XOStn), &
	StnNewD(CandidDStn),LatD(CandidDStn),LonD(CandidDStn),nint(ElvD(CandidDStn)), &
	StnNameD(CandidDStn),StnCtyD(CandidDStn)
else if (MeanR.LT.0.2) then
  Overlapper=.FALSE. ; Butter=.FALSE.	! bad overlap - no match!
  write (99,"(a10,2i8,f7.2,f8.2,i5,x,a20,x,a13)"), "  uncorrel", StnOldO(XOStn), &
	StnNewD(CandidDStn),LatD(CandidDStn),LonD(CandidDStn),nint(ElvD(CandidDStn)), &
	StnNameD(CandidDStn),StnCtyD(CandidDStn)
else if (MeanR.LT.0.4) then
  Overlapper=.FALSE. ; Butter=.TRUE.	! doubtful overlap - butt join
else
  Overlapper=.TRUE.  ; Butter=.FALSE.	! excellent overlap - accept!
end if

end subroutine CheckOverlap

!*******************************************************************************
! conclude program

subroutine Conclude

print*
close (99)

end subroutine Conclude

!*******************************************************************************
! conclude stn assignment

subroutine ConcludeAssign

if      (Pseudo) then
  LogText = "pseudo"	; Log1=Log1+1
else if (Adder) then
  LogText = "new-entry"	; Log2=Log2+1
else if (Overlapper) then
  LogText = "overlap"	; Log3=Log3+1
else if (Butter) then
  LogText = "butt-join"	; Log4=Log4+1
end if

write (99,"(a10,2i8,f7.2,f8.2,i5,x,a20,x,a13)"), LogText, StnOldO(XOStn), &
	StnNewD(AssignDStn),LatO(XOStn),LonO(XOStn),nint(ElvO(XOStn)), &
	StnNameO(XOStn),StnCtyO(XOStn)
	
end subroutine ConcludeAssign

!*******************************************************************************
! identify database file

subroutine FindDatabase

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

end subroutine FindDatabase

!*******************************************************************************
! finds next available pseudo code, and place in StnNewM(AssignMStn)

subroutine FindNextPseudo

Seek = 0-StnCodeBegO(XOStn)
do
  Find = LocateStn(Seek,StnNewM)
  if (Find.NE.MissVal) Seek = Seek - 1
  if (Find.EQ.MissVal) exit
end do
StnNewM(AssignMStn) = Seek

end subroutine FindNextPseudo

!*******************************************************************************
! examine new entry; if there is nothing to be learned, make this overlap

subroutine GetAdditions

Extrer=.FALSE. ; XODYear=0
do
  XODYear=XODYear+1 ; XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
  do XMonth = 1, NMonth
    if (	DataD(XDYear,XMonth,CandidDStn).EQ.DataMissVal .AND. &
    		DataO(XOYear,XMonth,XOStn).NE.DataMissVal) Extrer=.TRUE.
  end do
  if (XODYear.EQ.NODYear.OR.Extrer) exit
end do

end subroutine GetAdditions

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

subroutine GetButtJoin

allocate (RefTS(NDYear,NMonth,NDStn),Trust(NDYear,NDStn),Calc(NDStn),GotRef(NDStn), &
	  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### GetButtJoin: alloc #####"
RefTS=DataMissVal ; Trust=.TRUE. ; Calc=.FALSE. ; GotRef=.FALSE.
Calc(CandidDStn)=.TRUE. ; NeedThresh=2

Increment= 1 ; call NeededStretch
if (ClimBeg.NE.MissVal.AND.ClimEnd.NE.MissVal) call MakeRefTS &
		(RefTS,DataD,LatD,LonD,Calc,Trust,real(DecayDistance),Differ, &
		.TRUE.,nint(MissVal),GotRef,NeedBeg=ClimBeg,NeedEnd=ClimEnd)

if (.NOT.GotRef(CandidDStn)) then
  Increment=-1 ; call NeededStretch
  if (ClimBeg.NE.MissVal.AND.ClimEnd.NE.MissVal) call MakeRefTS &
		(RefTS,DataD,LatD,LonD,Calc,Trust,real(DecayDistance),Differ, &
		.TRUE.,nint(MissVal),GotRef,NeedBeg=ClimBeg,NeedEnd=ClimEnd)
end if

if (GotRef(CandidDStn)) then
  do XMonth=1,NMonth
    VecO=MissVal ; VecD=MissVal ; VecR=MissVal
  
    do XODYear=1,NODYear				! set up vectors
      XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
      if (RefTS(XDYear,XMonth,CandidDStn).NE.DataMissVal) &
    		VecR(XODYear)=real(RefTS(XDYear,XMonth,CandidDStn))
      if (DataD(XDYear,XMonth,CandidDStn).NE.DataMissVal) &
    		VecD(XODYear)=real(DataD(XDYear,XMonth,CandidDStn))
      if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) &
    		VecO(XODYear)=real(DataO(XOYear,XMonth,XOStn))
    end do
  
    if (Differ) then	! adjust RefTS to match Database, then Orig to match RefTS
      call MergeForDiff  (VecD,VecR,1,NODYear)		! in ghcnrefiter.f90
      call MergeForDiff  (VecR,VecO,1,NODYear)
    else
      call MergeForRatio (VecD,VecR,1,NODYear)
      call MergeForRatio (VecR,VecO,1,NODYear)
    end if
    
    do XODYear=1,NODYear				! in-fill data
      XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
      if (	DataD   (XDYear,XMonth,CandidDStn).NE.DataMissVal .AND. &
      		DataSrcD(XDYear,XMonth,CandidDStn).EQ.DataMissVal) then
        if (VecO(XODYear).NE.MissVal) then
          if      (VecO(XODYear).LT.Min(XMonth)) then
            VecO(XODYear)=Min(XMonth)
          else if (VecO(XODYear).GT.Max(XMonth)) then
            VecO(XODYear)=Max(XMonth)
	  end if
        end if
      
        DataD    (XDYear,XMonth,CandidDStn) = nint(VecO(XODYear))
        DataSrcD (XDYear,XMonth,CandidDStn) = DataSrcO(XOYear,XMonth,XOStn)
      end if
    end do
  end do
  
  call AssignStnMeta
else if (QRelyMeta.EQ.1) then
  do XMonth=1,NMonth
    do XODYear=1,NODYear				! in-fill data
      XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
      if (	DataD   (XDYear,XMonth,CandidDStn).NE.DataMissVal .AND. &
      		DataSrcD(XDYear,XMonth,CandidDStn).EQ.DataMissVal) then
        DataSrcD (XDYear,XMonth,CandidDStn) = DataSrcO(XOYear,XMonth,XOStn)
      end if
    end do
  end do
  
  call AssignStnMeta
else
  write (99,"(a10,2i8,f7.2,f8.2,i5,x,a20,x,a13)"), "  no-refts", StnOldO(XOStn), &
	StnNewD(CandidDStn),LatD(CandidDStn),LonD(CandidDStn),nint(ElvD(CandidDStn)), &
	StnNameD(CandidDStn),StnCtyD(CandidDStn)

  do XMonth=1,NMonth
    do XDYear=1,NDYear
      if (	DataD   (XDYear,XMonth,CandidDStn).NE.DataMissVal .AND. &
      		DataSrcD(XDYear,XMonth,CandidDStn).EQ.DataMissVal) then
        DataD(XDYear,XMonth,CandidDStn) = DataMissVal
      end if
    end do
  end do
  
  Butter=.FALSE.
end if

deallocate (RefTS,Trust,Calc,GotRef,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### GetButtJoin: dealloc #####"

end subroutine GetButtJoin

!*******************************************************************************
! identifies a single src code for each stn

subroutine IdentifyStnSrc

allocate (StnSrcCodeO(NOStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: IdentifyStnSrc: Alloc #####"

if (QOrigSrc.EQ.0) then
  StnSrcCodeO = OrigSrcCode
else
  StnSrcCodeO = DataMissVal
  do XOStn = 1, NOStn		! if one source in a stn, make StnSrcCodeO = src
   do XOYear = 1, NOYear	! 	else make StnSrcCodeO = MultiSource
    do XMonth = 1, NMonth
      if (DataSrcO(XOYear,XMonth,XOStn).NE.DataMissVal) then
        if      (StnSrcCodeO(XOStn).EQ.DataMissVal) then
          StnSrcCodeO(XOStn) = DataSrcO(XOYear,XMonth,XOStn)
        else if (StnSrcCodeO(XOStn).NE.DataSrcO(XOYear,XMonth,XOStn)) then
          StnSrcCodeO(XOStn) = MultiSource
        end if
      end if
    end do
   end do
  end do
end if
  
end subroutine IdentifyStnSrc

!*******************************************************************************
! initialises program

subroutine Intro

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

print*
print*, "  > ***** UpdateDTB: processes new stn files *****"
print*

NMonth = 12 ; Log1=0 ; Log2=0 ; Log3=0 ; Log4=0

call date_and_time (CurrentDate, CurrentTime)
Stamp10  = CurrentDate(3:8) // CurrentTime(1:4) 
Stamp20  = adjustr(Stamp10)
StampInt = GetIntFromText(Stamp20)

end subroutine Intro

!*******************************************************************************
! loading accessions file

subroutine LoadingAccessions

AccessFileL="./../../../data/cruts/accession/??????????.hdr"
call GetBatch (AccessFileL,Batch,Silent=1)
NBatch = size(Batch,1)
AccessFileL = Batch(NBatch) 
SubBeg=index(AccessFileL,"/",.TRUE.)
FileDate=AccessFileL((SubBeg+1):(SubBeg+6)) 
FileTime=AccessFileL((SubBeg+7):(SubBeg+10))

call LoadCTS (TrashInfo,StnLocalA,StnNameA,StnCtyA,Code=StnOldA,OldCode=StnNewA, &
		Silent=1,Lat=LatA,Lon=LonA,Elv=ElvA,HeadOnly=1, &
		CallFile=AccessFileL,Extra=NOStn,SrcCode=StnSrcCodeA, &
		SrcSuffix=StnSrcSuffixA,SrcDate=StnSrcDateA)
NAStn = size(StnOldA,1)
where (StnSrcCodeA.EQ.MissVal) StnSrcCodeA=DataMissVal

print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a)", "   > Accessions file on ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," at ", &
		FileTime(1:2),":",FileTime(3:4)," has ",(NAStn-NOStn)," stns"
deallocate (TrashInfo,Batch,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingAccessions: Dealloc #####"

end subroutine LoadingAccessions

!*******************************************************************************
! loading database files (.dtb and .dts)

subroutine LoadingDatabase

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

call LoadCTS (TrashInfo,TrashLocal,TrashName,TrashCty,Data=DataD,YearAD=YearADD, &
		DtbNormals=DataNmlD,YearADMin=YearADO(1),YearADMax=YearADO(NOYear), &
		Silent=1,CallFile=DatabaseFileL,Extra=NOStn)
NDStn = size(DataD,3) ; NDYear = size(YearADD,1)
deallocate (TrashInfo,TrashLocal,TrashName,TrashCty,YearADD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingDatabase: Dealloc 1 #####"

print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a,i10,a)", "   > Database   file on ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," at ", &
		FileTime(1:2),":",FileTime(3:4)," has ",(NDStn-NOStn)," stns;", &
		count(DataD.NE.DataMissVal), " data"

SubBeg = index(DatabaseFileL,".dtb")
DatabaseSrcFileL = DatabaseFileL(1:SubBeg) // "dts"
call LoadCTS (TrashInfo,StnLocalD,StnNameD,StnCtyD,Code=StnNewD, &
		Lat=LatD,Lon=LonD,Elv=ElvD,DtbNormals=DataNmlSrcD,Data=DataSrcD, &
		YearAD=YearADD,CallFile=DatabaseSrcFileL,Extra=NOStn,Silent=1, &
		YearADMin=YearADO(1),YearADMax=YearADO(NOYear))
deallocate (TrashInfo,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingDatabase: Dealloc 2 #####"

end subroutine LoadingDatabase

!*******************************************************************************
! load master file

subroutine LoadingMaster

MasterFileL="./../../../data/cruts/master/??????????.hdr"
call GetBatch (MasterFileL,Batch,Silent=1)
NBatch = size(Batch,1)
MasterFileL = Batch(NBatch)
SubBeg=index(MasterFileL,"/",.TRUE.)
FileDate=MasterFileL((SubBeg+1):(SubBeg+6)) 
FileTime=MasterFileL((SubBeg+7):(SubBeg+10))

call LoadCTS (TrashInfo,StnLocalM,StnNameM,StnCtyM,Code=StnNewM,Lat=LatM,Lon=LonM,Elv=ElvM, &
			HeadOnly=1,CallFile=MasterFileL,Extra=NOStn,Silent=1)
NMStn = size(StnNewM,1)

print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a)", "   > Master     file at ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," at ", &
		FileTime(1:2),":",FileTime(3:4)," has ",(NMStn-NOStn)," stns"

deallocate (TrashInfo,Batch,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingMaster: Dealloc #####"

end subroutine LoadingMaster

!*******************************************************************************
! load original data file

subroutine LoadingOriginal

call LoadCTS (TrashInfo,StnLocalO,StnNameO,StnCtyO,OldCode=StnOldO, &
    			Lat=LatO,Lon=LonO,Elv=ElvO,Data=DataO, &
			YearAD=YearADO,CallFile=OrigFileL,Silent=1)

deallocate (TrashInfo,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingOriginal: Dealloc 1 #####"

NOStn = size(StnOldO,1) ; NOYear = size(DataO,1)
print "(a,22x,a,i7,a,i10,a)", "   > Original   file","has ", NOStn, " stns;", &
		count(DataO.NE.DataMissVal), " data"

call LoadCTS (TrashInfo,TrashLocal,TrashName,TrashCty,Data=DataSrcO, &
			YearAD=TrashYearAD,CallFile=OrigSrcFileL,Silent=1)
deallocate (TrashInfo,TrashLocal,TrashName,TrashCty,TrashYearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingOriginal: Dealloc 2 #####"

end subroutine LoadingOriginal

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

function LocateStn (WantedCode,Stn)

integer :: LocateStn,WantedCode,XStn,NStn
integer, pointer, dimension (:) :: Stn
real, parameter :: MissVal=-999.0

NStn=size(Stn,1) ; XStn=0 ; LocateStn=0
do
  XStn = XStn + 1
  if (Stn(XStn).EQ.MissVal) LocateStn=MissVal
  if (Stn(XStn).EQ.WantedCode) LocateStn=XStn
  if (XStn.EQ.NStn.OR.LocateStn.NE.0) exit
end do

end function LocateStn

!*******************************************************************************
! make new database filled with missing values

subroutine MakeNewDatabase

NDStn = NOStn ; NDYear = NOYear
allocate (StnLocalD(NDStn), StnNameD(NDStn), StnCtyD(NDStn), StnNewD(NDStn), &
	  LatD(NDStn), LonD(NDStn), ElvD(NDStn), YearADD(NDYear), &
	  DataD(NDYear,NMonth,NDStn),DataSrcD(NDYear,NMonth,NDStn), &
	  DataNmlD(NMonth,NDStn),DataNmlSrcD(NMonth,NDStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeNewDatabase: Alloc #####"

StnLocalD="   nocode" ; StnNameD="UNKNOWN" ; StnCtyD="UNKNOWN" ; StnNewD=MissVal
LatD=MissVal ; LonD=MissVal ; ElvD=MissVal
DataD=DataMissVal ; DataSrcD=DataMissVal ; DataNmlD=DataMissVal ; DataNmlSrcD=DataMissVal

do XDYear = 1, NDYear
  YearADD(XDYear) = YearADO(XDYear)
end do

end subroutine MakeNewDatabase

!*******************************************************************************
! identify master candidates from priors

subroutine MatchCandid

CandidIndexA=MissVal ; CandidIndexM=MissVal ; CandidScoreA=MissVal ; NCandid=0

if (NPrior.GE.1) then					! relevant prior accessions
  call QuickSort (Reals=MatchScoreA,OrderValid=OrderA)	! sort the access scores
  
  do XPrior = 1, NPrior					! iterate by access score
    XMStn=0
    do							! iterate by master stn
     XMStn=XMStn+1
      
     if (StnNewM(XMStn).EQ.StnNewA(OrderA(XPrior))) then	! access and master stns match
        XCandid=0
        do
          XCandid=XCandid+1
          if (CandidIndexM(XCandid).EQ.XMStn)   exit	! candidate already selected
          if (CandidIndexM(XCandid).EQ.MissVal) exit	! new candidate
        end do
        
        if (CandidIndexM(XCandid).EQ.MissVal) then	! new candidate
          CandidScoreA(XCandid)=MatchScoreA(OrderA(XPrior))	! candidate score
          CandidIndexA(XCandid)=OrderA(XPrior)			! stn-index in accession
          CandidIndexM(XCandid)=XMStn 				! stn-index in master
          NCandid=NCandid+1
        end if
     end if
      
     if (XMStn.EQ.NMStn.OR.StnNewM(XMStn).EQ.StnNewA(OrderA(XPrior))) exit
    end do
  end do
end if

end subroutine MatchCandid

!*******************************************************************************
! match code with accessions file

subroutine MatchCodeAccess

MatchCodeA = MissVal					! unmatched
if (StnOldO(XOStn).NE.MissVal) then			
  do XAStn = 1, NAStn
    if (StnOldA(XAStn).EQ.StnOldO(XOStn)) then		
    	if (abs(StnSrcCodeA(XAStn)).EQ.StnSrcCodeO(XOStn) &
    			.AND.StnSrcCodeO(XOStn).NE.MultiSource) then
    	  MatchCodeA(XAStn) = -25			! acc-code + src-code match
    	else
    	  MatchCodeA(XAStn) = -10 			! accession code match
        end if
    else if (StnOldA(XAStn).NE.MissVal) then
      MatchCodeA(XAStn) = 30				! bad match
    end if
  end do
end if

end subroutine MatchCodeAccess

!*******************************************************************************
! match code with master file

subroutine MatchCodeMaster

if (StnOldO(XOStn).NE.MissVal) then				! orig stn code OK
  if (StnOldA(CandidAStn).EQ.StnOldO(XOStn)) then			! matches acc code
    if (StnSrcCodeO(XOStn).NE.MultiSource.AND. &		! matches acc source
	abs(StnSrcCodeA(CandidAStn)).EQ.StnSrcCodeO(XOStn)) then
      CandidScoreM(XCandid) = -25
    else
      CandidScoreM(XCandid) = -10
    end if
  else if ((StnNewM(CandidMStn)*10).EQ.StnOldO(XOStn)) then		
      CandidScoreM(XCandid) = -10				! matches master
  else if (StnOldO(XOStn).LT.0.OR.mod(StnOldO(XOStn),10).NE.0) then
      CandidScoreM(XCandid) =   0				! match impossible
  else if (StnNewM(CandidMStn).LT.0) then
      CandidScoreM(XCandid) =   0				! match impossible
  else
      CandidScoreM(XCandid) =  30				! match poss, but none
  end if
else
      CandidScoreM(XCandid) =   0				! match impossible
end if

end subroutine MatchCodeMaster

!*******************************************************************************
! match location with accessions file

subroutine MatchDistanceAccess

MatchDistanceA = MissVal
if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal) then	! distance OK
 if (ZOneDP .AND.mod((100*LatO(XOStn)),10.0).EQ.0 &
 			.AND.mod((100*LonO(XOStn)),10.0).EQ.0) then
   do XAStn = 1, NAStn
     if (LatA(XAStn).NE.MissVal.AND.LonA(XAStn).NE.MissVal) then
         Dist2DP = GetDistance (LatA(XAStn),LonA(XAStn),LatO(XOStn),LonO(XOStn))
         Lat1DP  = real(nint(LatA(XAStn)*10))/10.0
         Lon1DP  = real(nint(LonA(XAStn)*10))/10.0
         Dist1DP = GetDistance (Lat1DP,Lon1DP,LatO(XOStn),LonO(XOStn))
       
         if (Dist1DP.LT.Dist2DP) then
           MatchDistanceA(XAStn) = Dist1DP
         else
           MatchDistanceA(XAStn) = Dist2DP
         end if
     end if
   end do
 else
   do XAStn = 1, NAStn
     if (LatA(XAStn).NE.MissVal.AND.LonA(XAStn).NE.MissVal) then
       if (mod((100*LatA(XAStn)),10.0).EQ.0 &
 			.AND.mod((100*LonA(XAStn)),10.0).EQ.0) then
         Dist2DP = GetDistance (LatA(XAStn),LonA(XAStn),LatO(XOStn),LonO(XOStn))
         Lat1DP  = real(nint(LatO(XOStn)*10))/10.0
         Lon1DP  = real(nint(LonO(XOStn)*10))/10.0
         Dist1DP = GetDistance (Lat1DP,Lon1DP,LatA(XAStn),LonA(XAStn))
       
         if (Dist1DP.LT.Dist2DP) then
           MatchDistanceA(XAStn) = Dist1DP
         else
           MatchDistanceA(XAStn) = Dist2DP
         end if
       else
         MatchDistanceA(XAStn) = GetDistance (LatA(XAStn),LonA(XAStn), &
         					LatO(XOStn),LonO(XOStn))
       end if
     end if
   end do
 end if
end if

end subroutine MatchDistanceAccess

!*******************************************************************************
! match stn name with accessions file

subroutine MatchNameAccess

MatchNameA = MissVal
if (trim(StnNameO(XOStn)).NE.'UNKNOWN') then			! name match OK
  do XAStn = 1, NAStn
    if (trim(StnNameA(XAStn)).NE.'UNKNOWN') then
     if (trim(StnNameA(XAStn)).EQ.trim(StnNameO(XOStn))) then
      MatchNameA(XAStn) = 0.0
     else
      Name = StnNameO(XOStn)
      NChar = len_trim(adjustl(Name)) ; XChar = 0 ; SubLen = 0
      do
        XChar = XChar + 1
        SubBeg = index (StnNameA(XAStn),Name(1:XChar))
        if (SubBeg.NE.0) SubLen = SubLen + 1
        if (XChar.EQ.NChar.OR.SubBeg.EQ.0) exit
      end do
      MatchNameA(XAStn) = 10 - SubLen
     end if
    end if
  end do
end if

end subroutine MatchNameAccess

!*******************************************************************************
! calc match score for each access file stn

subroutine MatchScoreAccess

MatchScoreA = MissVal ; NPrior = 0
do XAStn = 1, NAStn
 if (MatchCodeA(XAStn).NE.MissVal) then
  MatchScoreA(XAStn) = MatchCodeA(XAStn)
  
  if (MatchDistanceA(XAStn).NE.MissVal) then
        MatchScoreA(XAStn) = MatchScoreA(XAStn) + (MatchDistanceA(XAStn)-5)
  else
  	MatchScoreA(XAStn) = MatchScoreA(XAStn) + 5
  end if
  if (MatchNameA(XAStn).NE.MissVal) then
  	MatchScoreA(XAStn) = MatchScoreA(XAStn) + (MatchNameA(XAStn)-5)
  else
  	MatchScoreA(XAStn) = MatchScoreA(XAStn) + 5
  end if
  
  if (ElvA(XAStn).NE.MissVal.AND.ElvO(XOStn).NE.MissVal) then
	OpVal = sqrt(abs(ElvA(XAStn)-ElvO(XOStn))+1)-3 
	if (OpVal.GT.5) OpVal=5
	MatchScoreA(XAStn) = MatchScoreA(XAStn) + OpVal
  end if
  if (StnCtyA(XAStn).EQ.StnCtyO(XOStn)) &
    	MatchScoreA(XAStn) = MatchScoreA(XAStn) - 1
  if (StnLocalA(XAStn).EQ.StnLocalO(XOStn)) &
    	MatchScoreA(XAStn) = MatchScoreA(XAStn) - 1

  if (MatchScoreA(XAStn).GE.30) then
      MatchScoreA(XAStn) = MissVal
  else
      NPrior = NPrior + 1
  end if
 end if
end do

end subroutine MatchScoreAccess

!*******************************************************************************
! obtains final master score for candidate

subroutine MatchScoreMaster

if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal.AND. &
		LatM(CandidMStn).NE.MissVal.AND.LonM(CandidMStn).NE.MissVal) then
 if (ZOneDP .AND.mod((100*LatO(XOStn)),10.0).EQ.0 &
 			.AND.mod((100*LonO(XOStn)),10.0).EQ.0) then
         Dist2DP = GetDistance (LatM(CandidMStn),LonM(CandidMStn), &
         			LatO(XOStn),LonO(XOStn))
         Lat1DP  = real(nint(LatM(CandidMStn)*10))/10.0
         Lon1DP  = real(nint(LonM(CandidMStn)*10))/10.0
         Dist1DP = GetDistance (Lat1DP,Lon1DP,LatO(XOStn),LonO(XOStn))
       
         if (Dist1DP.LT.Dist2DP) then
           Distance = Dist1DP
         else
           Distance = Dist2DP
         end if
 else
       if (mod((100*LatM(CandidMStn)),10.0).EQ.0 &
 			.AND.mod((100*LonM(CandidMStn)),10.0).EQ.0) then
         Dist2DP = GetDistance (LatM(CandidMStn),LonM(CandidMStn), &
         			LatO(XOStn),LonO(XOStn))
         Lat1DP  = real(nint(LatO(XOStn)*10))/10.0
         Lon1DP  = real(nint(LonO(XOStn)*10))/10.0
         Dist1DP = GetDistance (Lat1DP,Lon1DP,LatM(CandidMStn),LonM(CandidMStn))
       
         if (Dist1DP.LT.Dist2DP) then
           Distance = Dist1DP
         else
           Distance = Dist2DP
         end if
       else
         Distance = GetDistance (LatM(CandidMStn),LonM(CandidMStn), &
         			 LatO(XOStn),LonO(XOStn))
       end if
 end if

 CandidScoreM(XCandid) = CandidScoreM(XCandid) + Distance - 5
else
 if (CandidScoreM(XCandid).GT.-10) CandidScoreM(XCandid) = 3
end if

if (trim(StnNameO(XOStn)).NE.'UNKNOWN'.AND.trim(StnNameM(CandidMStn)).NE.'UNKNOWN') then
  if (trim(StnNameO(XOStn)).EQ.trim(StnNameM(CandidMStn))) then
    SubLen = 0.0
  else
    Name = StnNameO(XOStn)
    NChar = len_trim(adjustl(Name)) ; XChar = 0 ; SubLen = 0
    do
        XChar = XChar + 1
        SubBeg = index (StnNameM(CandidMStn),Name(1:XChar))
        if (SubBeg.NE.0) SubLen = SubLen + 1
        if (XChar.EQ.NChar.OR.SubBeg.EQ.0) exit
    end do
    
    SubLen = 10 - SubLen
    if (SubLen.LT.0) SubLen = 1.0
  end if
      
  CandidScoreM(XCandid) = CandidScoreM(XCandid) + SubLen - 5
end if

if (ElvM(CandidMStn).NE.MissVal.AND.ElvO(XOStn).NE.MissVal) then
  OpVal = sqrt (abs(ElvM(CandidMStn)-ElvO(XOStn)) + 1) - 3
  if (OpVal.GT.5) OpVal=5
  CandidScoreM(XCandid) = CandidScoreM(XCandid) + OpVal
end if

if (StnCtyM(CandidMStn).EQ.StnCtyO(XOStn)) &
    		CandidScoreM(XCandid) = CandidScoreM(XCandid) - 1
if (StnLocalM(CandidMStn).EQ.StnLocalO(XOStn)) &
    		CandidScoreM(XCandid) = CandidScoreM(XCandid) - 1

end subroutine MatchScoreMaster

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

subroutine MergeOverlap

do XMonth = 1, NMonth		! iterate by month
  VecO=MissVal ; VecD=MissVal
  
  do XODYear=1,NODYear			! set up vectors for this month
    XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
    if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) &
    		VecO(XODYear)=real(DataO(XOYear,XMonth,XOStn))
    if (DataD(XDYear,XMonth,AssignDStn).NE.DataMissVal) &
    		VecD(XODYear)=real(DataD(XDYear,XMonth,AssignDStn))
  end do
  
  if (Differ) then
    call MergeForDiff  (VecD,VecO,1,NODYear)		! in ghcnrefiter.f90
  else
    call MergeForRatio (VecD,VecO,1,NODYear)
  end if
  
  do XODYear=1,NODYear
    XDYear=DStart+XODYear-1 ; XOYear=OStart+XODYear-1
    if (DataD(XDYear,XMonth,AssignDStn).EQ.DataMissVal.AND. &
        DataO(XOYear,XMonth,XOStn)     .NE.DataMissVal) then
      if (VecO(XODYear).NE.MissVal) then
        if      (VecO(XODYear).LT.Min(XMonth)) then
          VecO(XODYear)=Min(XMonth)
        else if (VecO(XODYear).GT.Max(XMonth)) then
          VecO(XODYear)=Max(XMonth)
	end if
      end if
      
      DataD    (XDYear,XMonth,AssignDStn) = nint(VecO(XODYear))
      DataSrcD (XDYear,XMonth,AssignDStn) = DataSrcO(XOYear,XMonth,XOStn)
    end if
  end do
end do

end subroutine MergeOverlap

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

subroutine MergeExisting

do XODYear = 1, NODYear
  XDYear = DStart + XODYear - 1 ; XOYear = OStart + XODYear - 1
  do XMonth = 1, NMonth
    if (DataD(XDYear,XMonth,AssignDStn).EQ.DataMissVal) then
      if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) then
        if      (Butter) then				! butt-join
	  DataD(XDYear,XMonth,AssignDStn) = DataO(XOYear,XMonth,XOStn)
	else if (Adder) then				! new stn for database
	  DataD(XDYear,XMonth,AssignDStn) = DataO(XOYear,XMonth,XOStn)
	else if (Pseudo) then				! new pseudo code
	  DataD(XDYear,XMonth,AssignDStn) = DataO(XOYear,XMonth,XOStn)
	end if
	
	if (DataD(XDYear,XMonth,AssignDStn).LT.Min(XMonth)) &
			DataD(XDYear,XMonth,AssignDStn) = Min(XMonth)
	if (DataD(XDYear,XMonth,AssignDStn).GT.Max(XMonth)) &
			DataD(XDYear,XMonth,AssignDStn) = Max(XMonth)
	
	if (.NOT.Butter) DataSrcD(XDYear,XMonth,AssignDStn) = DataSrcO(XOYear,XMonth,XOStn)
      end if
    end if
  end do
end do

end subroutine MergeExisting

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

subroutine NeededStretch

MonthPrior=0 ; MonthFresh=0 ; NuffFresh=0 ; NuffPrior=0
if (Increment.GT.0) XDYear=0
if (Increment.LT.0) XDYear=NDYear+1

do						
      XDYear=XDYear+Increment ; XMonth=0
      do
        XMonth=XMonth+1
        if (	DataD(XDYear,XMonth,CandidDStn).NE.DataMissVal.AND. &
        	DataSrcD(XDYear,XMonth,CandidDStn).EQ.DataMissVal) then
          MonthFresh(XMonth)=MonthFresh(XMonth)+1
          if (MonthFresh(XMonth).EQ.NeedThresh) NuffFresh=NuffFresh+1
        end if
        if (XMonth.EQ.NMonth) exit
      end do
      if (Increment.GT.0.AND.XDYear.EQ.NDYear) exit
      if (Increment.LT.0.AND.XDYear.EQ.1) exit
      if (NuffFresh.EQ.12) exit
end do
    
if (NuffFresh.EQ.12) then			
     if (Increment.GT.0) ClimEnd=XDYear
     if (Increment.LT.0) ClimBeg=XDYear
     
     do
      XDYear=XDYear-Increment ; XMonth=0
      do
        XMonth=XMonth+1
        if (	DataD(XDYear,XMonth,CandidDStn).NE.DataMissVal.AND. &
        	DataSrcD(XDYear,XMonth,CandidDStn).NE.DataMissVal) then
          MonthPrior(XMonth)=MonthPrior(XMonth)+1
          if (MonthPrior(XMonth).EQ.NeedThresh) NuffPrior=NuffPrior+1
        end if
        if (XMonth.EQ.NMonth) exit
      end do
      if (Increment.LT.0.AND.XDYear.EQ.NDYear) exit
      if (Increment.GT.0.AND.XDYear.EQ.1) exit
      if (NuffPrior.EQ.12) exit
     end do
end if
    
if (NuffPrior.EQ.12) then
     if (Increment.GT.0) ClimBeg=XDYear
     if (Increment.LT.0) ClimEnd=XDYear
else
     ClimBeg=MissVal ; ClimEnd=MissVal
end if

end subroutine NeededStretch

!*******************************************************************************
! prepares stn for assignment

subroutine PrepareAssign

if (mod(XOStn,10).EQ.0.0) then		! progress marker
  write (99,*), ""
  write (99,"(2(a,i5),a,i3,a)"), "----- STATION ", XOStn, " OF ", NOStn, " (", &
   				nint(100.0*(real(XOStn)/real(NOStn))), "%) -----"
end if

QDecay=0				! no decay dist calc yet

Overlapper=.FALSE. ; Butter=.FALSE. ; Adder=.FALSE. ; Pseudo=.FALSE.
      						
end subroutine PrepareAssign

!*******************************************************************************
! prepares for saving

subroutine PrepareSave

print*, "  > Saving consolidated files ..."
MasterFileS      = "./../../../data/cruts/master/"    // Stamp10 // ".hdr"
AccessFileS      = "./../../../data/cruts/accession/" // Stamp10 // ".hdr"
DatabaseFileS    = "./../../../data/cruts/database/"  // VariableSuffix(2:4) // &
			"." // Stamp10 // ".dtb"
DatabaseSrcFileS = "./../../../data/cruts/database/"  // VariableSuffix(2:4) // &
			"." // Stamp10 // ".dts"

end subroutine PrepareSave

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

subroutine ProcessOriginal

call ClarifyCty (StnCtyO,StnNameO,LatO,LonO,StnOldO,StnCtyM,StnNameM, &
  			LatM,LonM,StnCodeBegO,StnCodeEndO,StnContinentO)

OpZero=0 ; OpTot=0
do XOStn=1,NOStn
  if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal) then
    OpTot=OpTot+1
    if (mod((100*LatO(XOStn)),10.0).EQ.0.AND.mod((100*LonO(XOStn)),10.0).EQ.0) &
    		OpZero=OpZero+1
  end if
end do
if (OpZero.GT.(OpTot*0.10)) then
  ZOneDP=.TRUE.
  print "(2(a,i4),a)", "   > Many stns have lat/lon to only 1 d.p. (", &
  	nint(OpZero), "/", nint(OpTot), ")"
else
  ZOneDP=.FALSE.
  print "(2(a,i4),a)", "   > Few  stns have lat/lon to only 1 d.p. (", &
  	nint(OpZero), "/", nint(OpTot), ")"
end if

end subroutine ProcessOriginal

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

subroutine SavingAccession

allocate (MissInts(NAStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingAccession: Alloc #####"
MissInts = MissVal
where (StnSrcCodeA.EQ.DataMissVal) StnSrcCodeA=MissVal

call MakeStnInfo (Info,StnOldA,StnNewA,LatA,LonA,ElvA,0,YearAD0=MissInts,YearAD1=MissInts,Silent=1)
call SaveCTS (Info,StnLocalA,StnNameA,StnCtyA,CallFile=AccessFileS, &
			HeadOnly=1,Silent=1,SrcCode=StnSrcCodeA, &
			SrcSuffix=StnSrcSuffixA,SrcDate=StnSrcDateA)

call system ('compress ' // trim(AccessFileL) // ' &')	! compress old accessions file
call system ('chmod ug=rw ' // trim(AccessFileS))	! set permissions for new acc file
! call system ('sort -1n -o ' // trim(AccessFileS) // &	! compaq f90 version
!		' ' // trim(AccessFileS) // ' &') ! sort
call system ('sort -g -o ' // trim(AccessFileS) // &	! portland group f90 version
		' ' // trim(AccessFileS) // ' &') ! sort

deallocate (MissInts,Info,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingAccession: Dealloc #####"

end subroutine SavingAccession

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

subroutine SavingDatabase

allocate (MissInts(NDStn),OrderD(NDStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingDatabase: Alloc #####"
MissInts = MissVal

call QuickSort (Ints=StnNewD,OrderValid=OrderD)

SubBeg=index(DatabaseFileS,"/",.TRUE.)
FileDate=DatabaseFileS((SubBeg+5):(SubBeg+10)) 
FileTime=DatabaseFileS((SubBeg+11):(SubBeg+14))
print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a,i10,a)", "   > Database   file on ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," at ", &
		FileTime(1:2),":",FileTime(3:4)," has ",count(StnNewD.NE.MissVal), &
		" stns;", count(DataD.NE.DataMissVal), " data"

call MakeStnInfo (Info,StnNewD,MissInts,LatD,LonD,ElvD,1, &
		Data=DataD,DtbNormals=DataNmlD,YearAD=YearADD,Silent=1)
call SaveCTS (Info,StnLocalD,StnNameD,StnCtyD,CallFile=DatabaseFileS, &
		Data=DataD,DtbNormals=DataNmlD,YearAD=YearADD,Silent=1,Order=OrderD)
call SaveCTS (Info,StnLocalD,StnNameD,StnCtyD,CallFile=DatabaseSrcFileS, &
		Data=DataSrcD,DtbNormals=DataNmlSrcD,YearAD=YearADD,Silent=1,Order=OrderD)

if (DatabaseFileL.NE."") then
	call system ('compress ' // trim(DatabaseFileL) // ' &')
	call system ('compress ' // trim(DatabaseSrcFileL) // ' &')
end if
call system ('chmod ug=rw ' // trim(DatabaseFileS))
call system ('chmod ug=rw ' // trim(DatabaseSrcFileS))

deallocate (MissInts,Info,OrderD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingDatabase: Dealloc #####"

end subroutine SavingDatabase

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

subroutine SavingMaster

allocate (MissInts(NMStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingMaster: Alloc #####"
MissInts = MissVal

call MakeStnInfo (Info,StnNewM,MissInts,LatM,LonM,ElvM,0,YearAD0=MissInts,YearAD1=MissInts,Silent=1)
call SaveCTS (Info,StnLocalM,StnNameM,StnCtyM,CallFile=MasterFileS,HeadOnly=1,Silent=1)

call system ('compress ' // trim(MasterFileL) // ' &')		! compress old master file
call system ('chmod ug=rw ' // trim(MasterFileS))		! set perm for new master file
! call system ('sort -1n -o ' // trim(MasterFileS) // &		! compaq f90 version
!		' ' // trim(MasterFileS) // ' &')	
call system ('sort -g -o ' // trim(MasterFileS) // &		! portland group f90 version
		' ' // trim(MasterFileS) // ' &')	

deallocate (MissInts,Info,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SavingMaster: Delloc #####"

end subroutine SavingMaster

!*******************************************************************************
! provides user with load options and specifications

subroutine SelectOriginal

do
 do
   print*, "  > Select the .cts file to load: "
   read (*,*,iostat=ReadStatus), GivenFile
   if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
 end do
 OrigFileL = LoadPath (GivenFile,".cts")
 
 SubBeg = index(OrigFileL,"homog")
 if (SubBeg.LT.1) then
   print*, "  > File appears unhomogenised. Proceed ? (0=no,1=yes)"
   do
     read (*,*,iostat=ReadStatus), QProceed
     if (QProceed.EQ.1) SubBeg=1000
     if (ReadStatus.LE.0) exit
   end do
 end if
 
 if (SubBeg.GT.0) exit
end do

QOrigSrc  = 1			! i.e. parallel src file exists
SubBeg = index(OrigFileL,".cts")
OrigSrcFileL = OrigFileL(1:(SubBeg-1)) // ".src"
QOrigCode = 1			! i.e. original 7-digit code in the 'grid box' col

			! this is necessary for late sources like CLIMAT and MCDW
			!    for variables where there are not likely to be overlaps
			!    with stns in the dtb, because most stns end prior to 1990
			! so in this case, set to 1, where the presumption is that if
			!    a metadata match can be made, the data will match
			! this is valid in the context of correcting stns to match latest
			!    values because earlier-ending series will have been corrected
			!    to match the latest values in those series
do
  print*, "  > Allow non-overlap matches with dtb stns using metadata ? (0=no,1=yes)"
  read (*,*,iostat=ReadStatus), QRelyMeta
  if (ReadStatus.LE.0.AND.QRelyMeta.GE.0.AND.QRelyMeta.LE.1) exit
end do

end subroutine SelectOriginal

!*******************************************************************************
! select the climate variable, by suffix

subroutine SelectVariable

SubBeg = index(OrigFileL,".cts")
VariableSuffix = OrigFileL((SubBeg-4):(SubBeg-1))

do
   if (VariableSuffix.EQ."    ") then
    do
      print*, "  > Identify the variable suffix (form:'.???'): "
      read (*,*,iostat=ReadStatus), VariableSuffix
      if (ReadStatus.LE.0) exit
    end do
   end if
  
   call CheckVariSuffix (VariableSuffix,Variable,Factor)
   if (Variable.EQ."") VariableSuffix = "    "
   if (Variable.NE."") print "(2a)", "   > The variable is: ", trim(Variable)
   if (Variable.NE."") exit
end do

if      (VariableSuffix.EQ.".pre") then
  DecayDistance = 450  ; Min=0 ; Max=99900 ; Differ=.FALSE.
else if (VariableSuffix.EQ.".tmp") then
  DecayDistance = 1200 ; Min=-2730 ; Max=9000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".dtr") then 
  DecayDistance = 750  ; Min=0 ; Max=10000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".cld") then 
  DecayDistance = 750  ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".spc") then	  	! assumed, after cld
  DecayDistance = 750  ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".vap") then 
  DecayDistance = 1200 ; Min=0 ; Max=10000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".rhm") then 	! assumed, after vap, tmp
  DecayDistance = 1200 ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".wet") then 
  DecayDistance = 450  ; Min=0  ; Differ=.FALSE.
  Max=(/3100,2900,3100,3000,3100,3000,3100,3100,3000,3100,3000,3100/)
else
  print*, "  > @@@@@ ERROR: SelectVariable: no var @@@@@"
end if

Elastic=10 ; if (Differ) Elastic=5

end subroutine SelectVariable

!*******************************************************************************
! prepares arrays for assignment operations

subroutine SetupAssign

if (DatabaseFileL.EQ."") call MakeNewDatabase
if (DatabaseFileL.NE."") call LoadingDatabase
call CommonVecPer (YearADO,YearADD,OStart,OEnd,DStart,DEnd)
print*, "  > Assigning each stn in .cts file ..."

NODYear = OEnd-OStart+1

allocate (MatchCodeA(NAStn),MatchDistanceA(NAStn),MatchNameA(NAStn), &
	  MatchScoreA(NAStn),OrderA(NAStn),OrderM(NMStn), &
	  CandidIndexA(NAStn),CandidIndexM(NMStn), &
	  CandidScoreA(NAStn),CandidScoreM(NMStn),Exe(NODYear),Wye(NODYear), &
	  Combo(NODYear,NMonth),Decay(NDStn),VecO(NODYear),VecD(NODYear), &
	  VecR(NODYear),Normalise(NDStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SetupAssign: Alloc #####"
Normalise=.FALSE.

end subroutine SetupAssign

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

subroutine ShutDownAssign

deallocate (MatchCodeA,MatchDistanceA,MatchNameA,MatchScoreA,CandidIndexA, &
	    CandidIndexM,CandidScoreA,CandidScoreM, &
	    Combo,Decay,VecO,VecD,VecR,Normalise,OrderA,OrderM,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ShutDownAssign: Dealloc #####"

print "(a,4i5)", "   > Stns: pseudo,new-entry,overlap,butt-join:", Log1,Log2,Log3,Log4

end subroutine ShutDownAssign

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

end program UpdateDTB