! cleanmeta.f90
! f90 program written by Tim Mitchell on 14.05.02 as accession.f90
! rewritten as CleanMeta.f90 (March 2003), 
! then rewritten again as cleanmeta.f90 (June 2003)
! program to process new files of station monthly obs
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/cleanmeta time.f90 filenames.f90 grimfiles.f90 
!	crutsfiles.f90 gridops.f90 sortmod.f90 ctyfiles.f90 countries.f90 
!	./../cruts/cleanmeta.f90 2> /tyn1/tim/scratch/stderr.txt

program CleanMeta

use Time
use FileNames
use GrimFiles
use CRUtsFiles
use GridOps
use SortMod
use CtyFiles
use Countries

implicit none

logical, pointer, dimension (:)		:: Check

real, pointer, dimension (:,:)		:: RelA,RelB,RelR
real, pointer, dimension (:)		:: LatO,LonO,ElvO, LatM,LonM,ElvM
real, pointer, dimension (:)		:: LatA,LonA,ElvA
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

integer, pointer, dimension (:,:,:)	:: DataO,DataSrcO
integer, pointer, dimension (:,:)	:: DataNmlO
integer, pointer, dimension (:,:)	:: Info,TrashInfo
integer, pointer, dimension (:,:)	:: CandidExe,CandidWye,Combo
integer, pointer, dimension (:)		:: StnOldO,StnNewM,StnOldA,StnNewA
integer, pointer, dimension (:)         :: StnSrcCodeO,StnSrcCodeA,YearADO
integer, pointer, dimension (:)		:: StnCodeBegO,StnCodeEndO,StnContinentO
integer, pointer, dimension (:)		:: StnSrcDateA,MissInts,Source
integer, pointer, dimension (:) 	:: MasterCode0,MasterCode1,MasterContinent
integer, pointer, dimension (:) 	:: Order,DecayOrder,CandidIndexM,CandidIndexA

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

real, dimension (12)	:: MonthA,MonthB,MonthR,SumR
integer, dimension (12)	:: Min,Max

real, parameter :: DataMissVal 	= -9999.0
real, parameter :: MissVal 	= -999.0
real, parameter :: MultiSource 	= -888.0
real, parameter :: SumRThresh 	= 0.25

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

logical :: Invalid

real :: OpVal,OpTot,OpEn,OpTest,OpDiff,StdErr,Lower,Upper
real :: Factor,Distance,StnODigit6,RealVal,CorrelMean,WeightMean,DecayMemory,Estimate
real :: TotNull,TotNeg,TotSmall,MeanR, CurrentKm,DecimalKM

integer :: ReadStatus, AllocStat
integer :: QOrigCode,QOrigSrc,QDigits,QIncludeOrig,QConvertedCty,QCandid,QDecay
integer :: QDodgy,QAlter,QClarify
integer :: NOYear,NDYear,NOStn,NMStn,NAStn
integer :: XOYear,XDYear,XOStn,XMStn,XAStn
integer :: NChar,NBatch,NMonth,NMasterCty,NPrior,NCandid,NDecay,NRel,NSource,NSuspect
integer :: XChar,XBatch,XMonth,XMasterCty,XPrior,XCandid,XDecay,XRel,XSource,XSuspect
integer :: SubBeg,SubLen,Ascii,Code,OrigSrcCode,AnswerNoYes,Reverse,Find,Seek
integer :: OStart,OEnd,DStart,DEnd,Butt,ButtYX,DecayDistance,StampInt,Current,RunLen
integer :: CandidAStn,CandidDStn,CandidMStn,MinusSrc,ExtraSrc,SuspectBeg,SuspectEnd
integer :: AssignAStn,AssignDStn,AssignMStn, Log1,Log2,Log3,Log4

character (len=80) :: GivenFile,Variable
character (len=80) :: OrigFileL,OrigSrcFileL,MasterFileL,AccessFileL
character (len=80) :: OrigFileS,OrigSrcFileS,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			! initial tasks and options
call SelectOriginal
call SelectVariable

call LoadingOriginal		! load sufficent for the metadata clean + check
call LoadingMaster
call LoadingAccessions

call CleanNames			! metadata clean + check
call CleanOtherMeta
QClarify=1 ; call ClarifyDegrees
QClarify=2 ; call ClarifyDegrees
call CheckLocation				! main loop
if (VariableSuffix.NE.".pre".AND.VariableSuffix.NE.".wet") call CheckOtherMiss
call CheckMinData
call PrepareSave
call SavingOriginal

call Conclude

!*******************************************************************************
! CharByChar		goes through text chracter by character
! CheckLocation		checks lat/lon of each stn
! CleanNames		makes stn,cty names caps and no-hyphens
! CleanOtherMeta	cleans meta data	
! Conclude 		conclude program
! DumpDodgyLocations	dump any stns with dodgy locations to file
! FindDatabase		identify database file
! GetCtyDist		get average distance from mid-point by country
! GetFromMaster		total up all lat/lon in master code file by country
! GetFromOriginal	total up all lat/lon in original .cts file by country
! Intro			initialises program
! LoadingMaster		load master .hdr file
! LoadingOriginal	load original .cts file
! LogMarker		put progress marker in log file
! ManualMeta		allows manual changes to metadata
! PrepareSave		prepares for saving
! SavingOriginal	saves original .cts file, with metadata tidied up
! SelectOriginal	provides user with load options and specifications
! SelectProc		processing options
! SelectVariable	select climate variable, by suffix

contains

!*******************************************************************************
! goes through text chracter by character

subroutine CharByChar

NChar = len_trim(adjustl(Dirty))
if (NChar.EQ.0) then
    Dirty = "UNKNOWN" ; NChar = 7
end if
  
do XChar = 1, NChar			
    Ascii = iachar(Dirty(XChar:XChar))
    if      (Ascii.GE.97.AND.Ascii.LE.122) then 
    	Dirty(XChar:XChar) = achar(Ascii-32)		! make uppercase
    else if (Ascii.EQ.45) then
    	Dirty(XChar:XChar) = " "			! remove hyphenations
    else if (Ascii.EQ.42) then
    	Dirty(XChar:XChar) = " "			! remove asterisks
    else if (Ascii.EQ.35) then
    	Dirty(XChar:XChar) = "%"			! remove hashes
    else if (Ascii.GE. 1.AND.Ascii.LE.31) then 
    	Dirty(XChar:XChar) = " "			! remove funny characters
    else if (Ascii.LT. 1.OR.Ascii.GT.127) then
    	Dirty(XChar:XChar) = " "			! remove funny characters
    end if

    if (XChar.EQ.NChar) then			! remove trailing  
      if      (Ascii.EQ.38) then
    	Dirty(XChar:XChar) = " "			! ampersand
      else if (Ascii.EQ.40) then
    	Dirty(XChar:XChar) = " "			! round bracket
      else if (Ascii.EQ.91) then
    	Dirty(XChar:XChar) = " "			! square bracket
      else if (Ascii.EQ.44) then
    	Dirty(XChar:XChar) = " "			! comma
      end if
    end if
end do

NChar = len_trim(adjustl(Dirty))
if (NChar.EQ.0) then
    Dirty = "UNKNOWN" ; NChar = 7
end if

end subroutine CharByChar

!*******************************************************************************
! checks lat/lon of each stn

subroutine CheckLocation

print*, "  > Checking the lat/lon of each stn in .cts file ..."
QIncludeOrig=0 ; QDodgy=0
do
  if (associated(StnCodeBegO)) then
    deallocate (StnCodeBegO,StnCodeEndO,StnContinentO, stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: CheckLocation: Dealloc #####"
  end if
  
  call ClarifyCty (StnCtyO,StnNameO,LatO,LonO,StnOldO,StnCtyM,StnNameM, &
  			LatM,LonM,StnCodeBegO,StnCodeEndO,StnContinentO)
  call GetFromMaster
  call GetFromOriginal
  call GetCtyDist
  call DumpDodgyLocations

  if (QDodgy.GT.0) then
    QIncludeOrig=1
    print "(2a)", "   > Lat/lon check: eyeball ", trim(LatLonFile)
    call ManualMeta
  end if

  deallocate (MasterCty,MasterLat,MasterLon,MasterEn, &
  	      OriginalLat,OriginalLon,OriginalEn,CtyLon,CtyLat,CtyTot,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CheckLocation: Dealloc #####"
  
  if (QDodgy.EQ.0) exit
end do

end subroutine CheckLocation

!*******************************************************************************
! checks whether a stn has at least one datum per calendar month 
! 	if not, the stn is omitted, as useless, because it can never be merged

subroutine CheckMinData

do XOStn = 1, NOStn
  Invalid=.FALSE. ; XMonth=0
  do
    XMonth=XMonth+1 ; XOYear=0
    do
      XOYear=XOYear+1
      if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) exit
      if (XOYear.EQ.NOYear) exit
    end do
    if (DataO(XOYear,XMonth,XOStn).EQ.DataMissVal) Invalid=.TRUE.
    if (XMonth.EQ.NMonth.OR.Invalid) exit
  end do
  
  if (Invalid) StnOldO(XOStn)=MissVal
end do

end subroutine CheckMinData

!*******************************************************************************
! checks whether a stn has unregistered missing values posing as valid data

subroutine CheckOtherMiss

do XOStn = 1, NOStn
 do
  Invalid=.FALSE. ; XOYear=0 ; Current=DataMissVal ; RunLen=0
  do
    XOYear=XOYear+1 ; XMonth=0
    do
      XMonth=XMonth+1
      if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) then
        if (DataO(XOYear,XMonth,XOStn).EQ.Current) then
          RunLen=RunLen+1 ; if (RunLen.EQ.4) Invalid=.TRUE.
        else
          RunLen=1
        end if
        Current=DataO(XOYear,XMonth,XOStn)
      end if
      if (XMonth.EQ.NMonth.OR.Invalid) exit
    end do
    if (XOYear.EQ.NOYear.OR.Invalid) exit
  end do
  
  if (Invalid) then
    do XOYear=1,NOYear
      do XMonth=1,NMonth
        if (DataO(XOYear,XMonth,XOStn).EQ.Current) &
        		DataO(XOYear,XMonth,XOStn)=DataMissVal
      end do
    end do
  end if
  
  if (.NOT.Invalid) exit
 end do 
end do

end subroutine CheckOtherMiss

!*******************************************************************************
! checks whether degrees have a fraction or no. minutes after decimal point 

subroutine ClarifyDegrees

Lower=0 ; Upper=0
do XOStn = 1, NOStn
  if (LatO(XOStn).NE.MissVal.AND.LatO(XOStn).NE.0.0 .AND. &
      LonO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.0.0) then
    if (abs(mod(LatO(XOStn),1.0)).LT.0.6.AND.abs(mod(LonO(XOStn),1.0)).LT.0.6) then
      Lower=Lower+1
    else
      Upper=Upper+1
    end if
  end if
end do

if (Upper.EQ.0) then
  print*, "  > The lat/lon have minutes. Converting to decimal..."
  
  allocate (Check(NOStn),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ClarifyDegrees: alloc #####"
  Check=.TRUE.
  
  call ConvertDecimal
  
  deallocate (Check,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ClarifyDegrees: dealloc #####"
else if ((Lower+Upper).GT.0) then
  OpTest=Lower/(Lower+Upper)
  OpDiff=abs(OpTest-0.36)
  StdErr=sqrt((OpTest*(1.0-OpTest))/(Lower+Upper))
  
  if (OpDiff.GT.(3.0*StdErr)) then
    print*, "  > @@@@@ PROBLEM: Some stns have deg.min not deg.deg @@@@@"
    print "(a,3f8.2)", "   > @@@@@ test stat (stdev) =", OpDiff/StdErr,Lower,Upper

    if (QClarify.EQ.1) call SeekDegMinSrc
    if (QClarify.EQ.2) call SeekDegMinMaster
    if (QClarify.EQ.2) call SeekDegMinSection
  end if
end if

end subroutine ClarifyDegrees

!*******************************************************************************
! makes stn,cty names caps and no-hyphens

subroutine CleanNames

do XOStn = 1, NOStn
  Dirty = trim(adjustl(StnNameO(XOStn)))
  call CharByChar
  StnNameO(XOStn) = Dirty(1:20)
  
  Dirty = trim(adjustl(StnCtyO(XOStn))) ; Dirty(14:20) = "       "
  call CharByChar
  StnCtyO(XOStn) = Dirty(1:13)
end do

end subroutine CleanNames

!*******************************************************************************
! cleans meta data

subroutine CleanOtherMeta

do XOStn = 1, NOStn
  if (LatO(XOStn).LT. -90.OR.LatO(XOStn).GT.  90 .OR. &
      LonO(XOStn).LT.-180.OR.LonO(XOStn).GT. 180) then
      	LatO(XOStn) = MissVal ; LonO(XOStn) = MissVal
  end if
  if (ElvO(XOStn).LT.-400.OR.ElvO(XOStn).GT.6000) ElvO(XOStn) = MissVal
end do

end subroutine CleanOtherMeta

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

subroutine Conclude

print*
close (99)

end subroutine Conclude

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

subroutine ConvertDecimal

do XOStn = 1, NOStn
 if (Check(XOStn)) then
  if (	LatO(XOStn).NE.MissVal.AND.abs(mod(LatO(XOStn),1.0)).GE.0.6 .OR. &
  	LonO(XOStn).NE.MissVal.AND.abs(mod(LonO(XOStn),1.0)).GE.0.6) then
  		! nothing required here - already deg.deg
  else
    if (LatO(XOStn).NE.MissVal) LatO(XOStn) = &
    		LatO(XOStn)+(0.66666*mod(LatO(XOStn),1.0))
    if (LonO(XOStn).NE.MissVal) LonO(XOStn) = &
    		LonO(XOStn)+(0.66666*mod(LonO(XOStn),1.0))
  end if
 end if
end do

end subroutine ConvertDecimal

!*******************************************************************************
! dump any stns with dodgy locations to file

subroutine DumpDodgyLocations

open (1,file=LatLonFile,status="replace",action="write")
write (1,"(a7,a7,a8,a8,x,a20,x,a13,2a8)"), "7-digit","    lat","    long","     elv", &
        			"station name        ","country      ","distance","   rad*3 "
QDodgy=0
do XOStn = 1, NOStn			! check the lat/lon of each station
 if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal) then
  XMasterCty=0
  do					! find cty on master list
    XMasterCty=XMasterCty+1

    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn).AND.CtyTot(XMasterCty).GT.2.AND. &
    		CtyLat(XMasterCty).NE.MissVal.AND.CtyLon(XMasterCty).NE.MissVal) then
      
      Distance = GetDistance (CtyLat(XMasterCty),CtyLon(XMasterCty),LatO(XOStn),LonO(XOStn))
      
      if (Distance.GT.(CtyTot(XMasterCty)*3.0)) then
        StnODigit6 = real(StnOldO(XOStn))/10.0 ; QConvertedCty=1
        
        if (nint(StnODigit6).EQ.StnODigit6) then
          XMStn=0
          do
            XMStn=XMStn+1
            if (StnODigit6.EQ.StnNewM(XMStn).AND.StnNameO(XOStn).EQ.StnNameM(XMStn)) then
              StnCtyO(XOStn)=StnCtyM(XMstn) ; QConvertedCty=1
            end if
            if (StnODigit6.EQ.StnNewM(XMStn).OR.StnNewM(XMStn).EQ.MissVal) exit
          end do
        end if
        
        if (QConvertedCty.EQ.0) then
          write (1,"(i8,f7.2,f8.2,f8.1,x,a20,x,a13,2f8.1)"), &
          		StnOldO(XOStn),LatO(XOStn),LonO(XOStn),ElvO(XOStn), &
        		StnNameO(XOStn),StnCtyO(XOStn),Distance,(CtyTot(XMasterCty)*3.0)
          QDodgy=QDodgy+1
        end if
      end if
    end if
    
    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn).OR.XMasterCty.EQ.NMasterCty) exit
  end do
 end if
end do

close (1)

end subroutine DumpDodgyLocations

!*******************************************************************************
! get average distance from mid-point by country

subroutine GetCtyDist

allocate (CtyLat   (NMasterCty), &
	  CtyLon   (NMasterCty), &
	  CtyTot   (NMasterCty), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetCountryDistComb: Allocation failure #####"
CtyLat = 0.0 ; CtyLon = 0.0 ; CtyTot = 0.0

do XMasterCty=1,NMasterCty	! get lat/lon mid-point by country
  if      (QIncludeOrig.EQ.0.AND.MasterEn(XMasterCty).GT.0) then
    CtyLat(XMasterCty) = MasterLat(XMasterCty) / MasterEn(XMasterCty)
    CtyLon(XMasterCty) = MasterLon(XMasterCty) / MasterEn(XMasterCty)
  else if (QIncludeOrig.EQ.1.AND.(MasterEn(XMasterCty)+OriginalEn(XMasterCty)).GT.0) then
    CtyLat(XMasterCty) = (MasterLat(XMasterCty)+OriginalLat(XMasterCty)) / &
    				(MasterEn(XMasterCty)+OriginalEn(XMasterCty))
    CtyLon(XMasterCty) = (MasterLon(XMasterCty)+OriginalLon(XMasterCty)) / &
    				(MasterEn(XMasterCty)+OriginalEn(XMasterCty))
  end if
end do

do XMStn = 1, NMStn		! total up all distances from mid-point by country (master)
  XMasterCty=0
  do
    XMasterCty=XMasterCty+1
    if (MasterCty(XMasterCty).EQ.StnCtyM(XMStn)) then
      if (LatM(XMStn).NE.MissVal.AND.LonM(XMStn).NE.MissVal) CtyTot(XMasterCty) = &
      		CtyTot(XMasterCty) + &
      		GetDistance (CtyLat(XMasterCty),CtyLon(XMasterCty),LatM(XMStn),LonM(XMStn))
    end if 
    if (MasterCty(XMasterCty).EQ.StnCtyM(XMStn).OR.XMasterCty.EQ.NMasterCty) exit
  end do
end do

if (QIncludeOrig.EQ.1) then
 do XOStn = 1, NOStn		! total up all distances from mid-point by country (original)
  XMasterCty=0
  do
    XMasterCty=XMasterCty+1
    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn)) then
      if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal) CtyTot(XMasterCty) = &
      		CtyTot(XMasterCty) + &
      		GetDistance (CtyLat(XMasterCty),CtyLon(XMasterCty),LatO(XOStn),LonO(XOStn))
    end if 
    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn).OR.XMasterCty.EQ.NMasterCty) exit
  end do
 end do
end if

do XMasterCty=1,NMasterCty	! get average distance from mid-point by country
  if      (QIncludeOrig.EQ.0.AND.MasterEn(XMasterCty).GT.0) then
    CtyTot(XMasterCty) = CtyTot(XMasterCty) / MasterEn(XMasterCty)
  else if (QIncludeOrig.EQ.1.AND.(MasterEn(XMasterCty)+OriginalEn(XMasterCty)).GT.0) then
    CtyTot(XMasterCty) = CtyTot(XMasterCty) / (MasterEn(XMasterCty)+OriginalEn(XMasterCty))
  end if
end do

end subroutine GetCtyDist

!*******************************************************************************
! total up all lat/lon in master code file by country

subroutine GetFromMaster

NMasterCty = MasterSize()
allocate (MasterLat   (NMasterCty), &
	  MasterLon   (NMasterCty), &
	  MasterEn    (NMasterCty), &
	  MasterRawName(NMasterCty),MasterCty(NMasterCty),MasterCode0(NMasterCty), &
	  MasterCode1(NMasterCty),MasterContinent(NMasterCty),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFromMaster: Alloc #####"
MasterLat = 0.0 ; MasterLon = 0.0 ; MasterEn = 0.0

call LoadMasterCty (MasterRawName,MasterCty,MasterCode0,MasterCode1,MasterContinent)

deallocate (MasterRawName,MasterCode0,MasterCode1,MasterContinent, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFromMaster: Dealloc #####"

do XMStn = 1, NMStn
  XMasterCty=0
  do
    XMasterCty=XMasterCty+1
    if (MasterCty(XMasterCty).EQ.StnCtyM(XMStn)) then
      if (LatM(XMStn).NE.MissVal.AND.LonM(XMStn).NE.MissVal) then
        MasterLat(XMasterCty) = MasterLat(XMasterCty) + LatM(XMStn)
        MasterLon(XMasterCty) = MasterLon(XMasterCty) + LonM(XMStn)
        MasterEn (XMasterCty) = MasterEn (XMasterCty) + 1.0
      end if
    end if 
    if (MasterCty(XMasterCty).EQ.StnCtyM(XMStn).OR.XMasterCty.EQ.NMasterCty) exit
  end do
end do

end subroutine GetFromMaster

!*******************************************************************************
! total up all lat/lon in Original code file by country

subroutine GetFromOriginal

allocate (OriginalLat   (NMasterCty), &
	  OriginalLon   (NMasterCty), &
	  OriginalEn    (NMasterCty), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFromOriginal: Alloc #####"
OriginalLat = 0.0 ; OriginalLon = 0.0 ; OriginalEn = 0.0

do XOStn = 1, NOStn	
  XMasterCty=0
  do
    XMasterCty=XMasterCty+1
    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn)) then
      if (LatO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.MissVal) then
        OriginalLat(XMasterCty) = OriginalLat(XMasterCty) + LatO(XOStn)
        OriginalLon(XMasterCty) = OriginalLon(XMasterCty) + LonO(XOStn)
        OriginalEn (XMasterCty) = OriginalEn (XMasterCty) + 1.0
      end if
    end if 
    if (MasterCty(XMasterCty).EQ.StnCtyO(XOStn).OR.XMasterCty.EQ.NMasterCty) exit
  end do
end do

end subroutine GetFromOriginal

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

subroutine Intro

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

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

NMonth = 12

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)
NAStn = size(StnOldA,1)

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

end subroutine LoadingAccessions

!*******************************************************************************
! 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

if (QOrigCode.EQ.0) then		! code in main code column
  call LoadCTS (TrashInfo,StnLocalO,StnNameO,StnCtyO,Code=StnOldO, &
    			Lat=LatO,Lon=LonO,Elv=ElvO,Data=DataO, &
			YearAD=YearADO,CallFile=OrigFileL,Silent=1)
  if (QDigits.EQ.6) where (StnOldO.GT.0) StnOldO = StnOldO * 10
else					! code in 'grid-box' column
  call LoadCTS (TrashInfo,StnLocalO,StnNameO,StnCtyO,OldCode=StnOldO, &
    			Lat=LatO,Lon=LonO,Elv=ElvO,Data=DataO, &
			YearAD=YearADO,CallFile=OrigFileL,Silent=1)
end if

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

NOStn = size(StnOldO,1) ; NOYear = size(DataO,1)
print "(a,i7,2(a,i4))", "   > Original file has ", NOStn, " stns for ", &
			YearADO(1), "-", YearADO(NOYear)

if (QOrigSrc.EQ.1) then
  call LoadCTS (TrashInfo,TrashLocal,TrashName,TrashCty,Data=DataSrcO, &
			CallFile=OrigSrcFileL,Silent=1)
  deallocate (TrashInfo,TrashLocal,TrashName,TrashCty,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingOriginal: Dealloc 2 #####"
  
  MinusSrc = count(DataO.NE.DataMissVal.AND.DataSrcO.EQ.DataMissVal)
  ExtraSrc = count(DataO.EQ.DataMissVal.AND.DataSrcO.NE.DataMissVal)
  
  if (MinusSrc.GT.0) print "(a,i10)", "   > @@@@@ PROBLEM: data without src code:   ", MinusSrc
  if (ExtraSrc.GT.0) print "(a,i10)", "   > @@@@@ PROBLEM: src codes without datum: ", ExtraSrc
  
  if (MinusSrc.GT.0.OR.ExtraSrc.GT.0) then
!    print*, "  > Seeking problem..."		! @@@@@@@@@@@@@@@@@@@
    
    do XOStn=1,NOStn
      do XOYear=1,NOYear
        do XMonth=1,NMonth
	  if      (DataO   (XOYear,XMonth,XOStn).NE.DataMissVal.AND. &
	           DataSrcO(XOYear,XMonth,XOStn).EQ.DataMissVal) then
!	    print "(a,3i7)", "SRC ",StnOldO(XOStn),YearADO(XOYear),XMonth	! @@@@@@@@@@@@
	    write (99,"(a,3i7)"), "SRC ",StnOldO(XOStn),YearADO(XOYear),XMonth
	  else if (DataO   (XOYear,XMonth,XOStn).EQ.DataMissVal.AND. &
	           DataSrcO(XOYear,XMonth,XOStn).NE.DataMissVal) then
!	    print "(a,3i7)", "DAT ",StnOldO(XOStn),YearADO(XOYear),XMonth	! @@@@@@@@@@@@
	    write (99,"(a,3i7)"), "DAT ",StnOldO(XOStn),YearADO(XOYear),XMonth
	  end if
	end do
      end do
    end do
  end if
else
  allocate (DataSrcO(NOYear,NMonth,NOStn),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadingOriginal: Alloc 3 #####"
  DataSrcO=DataMissVal
  do XOYear = 1, NOYear
    do XMonth = 1, NMonth
      do XOStn = 1, NOStn
        if (DataO(XOYear,XMonth,XOStn).NE.DataMissVal) &
        		DataSrcO(XOYear,XMonth,XOStn) = OrigSrcCode
      end do
    end do
  end do
end if

end subroutine LoadingOriginal

!*******************************************************************************
! allows manual changes to metadata

subroutine ManualMeta

print*, "  > Manually change station meta information."
do
  print*, "  > Enter a station code (0=exit): "
  do
    read (*,*,iostat=ReadStatus), Code
    if (ReadStatus.LE.0) exit
  end do
  if (Code.NE.0) then
    XOStn = 0
    do
        XOStn = XOStn + 1
        if (XOStn.EQ.NOStn.AND.StnOldO(XOStn).NE.Code) Code = -1
        if (StnOldO(XOStn).EQ.Code) exit
    end do
      
    if (Code.NE.-1) then
      print "(a,i8,2f8.2,f8.1,x,a20,x,a13)", "   > ", StnOldO(XOStn),LatO(XOStn), &
      		LonO(XOStn),ElvO(XOStn),StnNameO(XOStn),StnCtyO(XOStn)
      print*, "  > Alter: code(=0), lat(=1), lon(=2), elv(=3) stn-name(=4), cty-name(=5): "
      do
        read (*,*,iostat=ReadStatus), QAlter
        if (ReadStatus.LE.0.AND.QAlter.GE.0.AND.QAlter.LE.5) exit
      end do
      if      (QAlter.EQ.0) then
        print*, "  > Enter the new station code: "
        do
          read (*,*,iostat=ReadStatus), StnOldO(XOStn)
          if (ReadStatus.LE.0) exit
        end do
      else if (QAlter.GE.1.AND.QAlter.LE.3) then
        print*, "  > Enter the new real value: "
        do
          read (*,*,iostat=ReadStatus), RealVal
          if (ReadStatus.LE.0) exit
        end do
        if (QAlter.EQ.1.AND.((RealVal.GE. -90.AND.RealVal.LE. 90).OR.RealVal.EQ.MissVal)) &
        	 	 LatO(XOStn)=RealVal
        if (QAlter.EQ.2.AND.((RealVal.GE.-180.AND.RealVal.LE.180).OR.RealVal.EQ.MissVal)) &
        		 LonO(XOStn)=RealVal
        if (QAlter.EQ.3) ElvO(XOStn)=RealVal
      else if (QAlter.EQ.4) then
        print*, "  > Enter the new station name: "
        do
          read (*,*,iostat=ReadStatus), StnNameO(XOStn)
          if (ReadStatus.LE.0) exit
        end do
      else if (QAlter.EQ.5) then
        if (associated(MasterCty)) then
          print*, "  > Select the new country (>0: see cty/master.txt): "
          do
            read (*,*,iostat=ReadStatus), XMasterCty
            if (ReadStatus.LE.0.AND.XMasterCty.GE.1.AND.XMasterCty.LE.NMasterCty) exit
          end do
          StnCtyO(XOStn) = MasterCty(XMasterCty)
        else
          print*, "  > Enter the new country name: "
          do
            read (*,*,iostat=ReadStatus), StnCtyO(XOStn)
            if (ReadStatus.LE.0) exit
          end do
        end if
      end if
      print "(a,i8,2f8.2,f8.1,x,a20,x,a13)", "   > ", StnOldO(XOStn),LatO(XOStn),&
      			LonO(XOStn),ElvO(XOStn),StnNameO(XOStn),StnCtyO(XOStn)
    end if
  end if
  if (Code.EQ.0) exit
end do

end subroutine ManualMeta

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

subroutine PrepareSave

print*, "  > The cleaned file goes in ./../../../data/cruts/clean/ ..."
SubBeg = index(OrigFileL,"/",.TRUE.)
SubLen = len_trim(OrigFileL)
GivenFile = "./../../../data/cruts/clean" // OrigFileL(SubBeg:SubLen)
OrigFileS = SavePath (GivenFile,".cts")

SubBeg = index(OrigFileS,".cts")
OrigSrcFileS = OrigFileS(1:SubBeg) // "src"

end subroutine PrepareSave

!*******************************************************************************
! saves original .cts file, with metadata tidied up

subroutine SavingOriginal

call QuickSort (Ints=StnOldO,Order=Order)

call MakeStnInfo (Info,StnOldO,StnOldO,LatO,LonO,ElvO,1,YearAD=YearADO,Data=DataO,Silent=1)
call SaveCTS (Info,StnLocalO,StnNameO,StnCtyO,Data=DataO,YearAD=YearADO, &
		CallFile=OrigFileS,Silent=1,Order=Order)		
call SaveCTS (Info,StnLocalO,StnNameO,StnCtyO,Data=DataSrcO,YearAD=YearADO, &
		CallFile=OrigSrcFileS,Silent=1,Order=Order)		

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

end subroutine SavingOriginal

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

subroutine SeekDegMinSection

Lower=0 ; SuspectBeg=0 ; SuspectEnd=0
do XOStn = 1, NOStn
  if (LatO(XOStn).NE.MissVal.AND.LatO(XOStn).NE.0.0 .AND. &
      LonO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.0.0) then
    if (abs(mod(LatO(XOStn),1.0)).LT.0.6.AND.abs(mod(LonO(XOStn),1.0)).LT.0.6) then
      if (Lower.LT. 5) Lower=Lower+1
    else 
      if (Lower.GT. 0) Lower=Lower-1
    end if
    
    if      (SuspectBeg.EQ.0) then
      if (Lower.GE.5) SuspectBeg=XOStn-4
    else if (Lower.LE.0) then
      SuspectEnd=XOStn-4
    else if (XOstn.EQ.NOStn) then
      SuspectEnd=NOStn
    end if
    
    if (SuspectEnd.GT.0) then
      print "(2(a,i7))", "   > @@@@@ CHECK stns between ", StnOldO(SuspectBeg), &
      			 " and ", StnOldO(SuspectEnd)
      SuspectBeg=0 ; SuspectEnd=0 ; Lower=0
    end if
  end if
end do

end subroutine SeekDegMinSection

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

subroutine SeekDegMinMaster

allocate (Check(NOStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeekDegMinMaster: alloc #####"
Check=.FALSE. ; Log1=0 ; Log2=0 ; Log3=0 ; Log4=0

do XOStn = 1, NOStn
  if (LatO(XOStn).NE.MissVal.AND.LatO(XOStn).NE.0.0 .AND. &
      LonO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.0.0) then
   if (abs(mod(LatO(XOStn),1.0)).LT.0.6.AND.abs(mod(LonO(XOStn),1.0)).LT.0.6) then

    XAStn=0 ; Log1=Log1+1
    do
      XAStn=XAStn+1
      if (StnOldA(XAStn).EQ.StnOldO(XOStn)) then
        Distance = GetDistance (LatA(XAStn),LonA(XAStn),LatO(XOStn),LonO(XOStn))
	Log2=Log2+1
	
	if (Distance.NE.MissVal.AND.Distance.LT.40) then

	  XMStn=0 ; Log3=Log3+1
	  do
	    XMStn=XMstn+1
	    if (StnNewM(XMStn).EQ.StnNewA(XAStn).AND. &
	    	LatM(XMStn).NE.MissVal.AND.LonM(XMStn).NE.MissVal) then
	      Log4=Log4+1
	      CurrentKM = GetDistance (LatM(XMStn),LonM(XMStn),LatO(XOStn),LonO(XOStn))
	      DecimalKM = GetDistance (LatM(XMStn),LonM(XMStn), &
	      			(LatO(XOStn)+(0.66666*mod(LatO(XOStn),1.0))), &
				(LonO(XOStn)+(0.66666*mod(LonO(XOStn),1.0))) )
						! @@@@@@@@@@@@@@@@@@@@@@@
	      write (99,"(2i8,6f8.2)"), StnOldO(XOStn),StnNewM(XMStn),CurrentKM,DecimalKM, &
	      		LatO(XOStn),LonO(XOStn),LatM(XMStn),LonM(XMStn)
	      
	      if (CurrentKM.GT.DecimalKM) Check(XOStn)=.TRUE.
	    end if
	    if (XMstn.EQ.NMStn) exit
	  end do

	end if
      end if
      if (StnOldA(XAStn).GT.StnOldO(XOStn)) exit
      if (XAStn.EQ.NAStn) exit
    end do
    
   end if
  end if
end do

NSuspect = count(Check)
! print "(4i10)", Log1,Log2,Log3,Log4		! @@@@@@@@@@@@@@@@@@@@@@@@
print "(a,i10)", "   > @@@@@ Decimalise from master file: ", NSuspect
call ConvertDecimal

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

end subroutine SeekDegMinMaster

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

subroutine SeekDegMinSrc

allocate (Source(NOStn),Check(NOStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeekDegMinSrc: alloc #####"
Source=MissVal ; Check=.TRUE.

do XOStn=1,NOStn
  do XOYear=1,NOYear
    if (Source(XOStn).NE.MultiSource) then
      do XMonth=1,NMonth
        if (DataSrcO(XOYear,XMonth,XOStn).NE.DataMissVal) then
	  if (Source(XOStn).EQ.MissVal) then
	    Source(XOStn) = DataSrcO(XOYear,XMonth,XOStn)
	  else if (DataSrcO(XOYear,XMonth,XOStn).NE.Source(XOStn)) then
	    Source(XOStn) = MultiSource
	  end if
	end if
      end do
    end if
  end do
end do

XSource=0
do
  XSource=XSource+1
  
  if (Source(XSource).NE.MissVal.AND.Check(XSource)) then
    Lower=0 ; Upper=0
    
    do XOStn = XSource, NOStn
      if (Source(XOStn).EQ.Source(XSource)) then
        Check(XOStn)=.FALSE.
	if (LatO(XOStn).NE.MissVal.AND.LatO(XOStn).NE.0.0 .AND. &
            LonO(XOStn).NE.MissVal.AND.LonO(XOStn).NE.0.0) then
          if (abs(mod(LatO(XOStn),1.0)).LT.0.6.AND.abs(mod(LonO(XOStn),1.0)).LT.0.6) then
	    Lower=Lower+1
	  else
	    Upper=Upper+1
	  end if
	end if
      end if
    end do
    
    OpTest=Lower/(Lower+Upper)
    OpDiff=abs(OpTest-0.36)
    StdErr=sqrt((OpTest*(1.0-OpTest))/(Lower+Upper))
  
    print "(a,i8,a,3f8.2)", "   > @@@@@ src code ", Source(XSource), " =", &
    			OpDiff/StdErr,Lower,Upper
  end if
  
  if (XSource.EQ.NOStn) exit
end do

do
  print*, "  > Select a src code to convert to decimal (-999=end): "
  read (*,*,iostat=ReadStatus), OpTest
  
  if (ReadStatus.LE.0) then
   if (OpTest.NE.MissVal) then
    Check=.FALSE.
    do XOStn=1,NOStn
      if (Source(XOStn).EQ.OpTest) Check(XOstn)=.TRUE.
    end do
    call ConvertDecimal
   end if
  end if
  
  if (OpTest.EQ.MissVal) exit
end do


deallocate (Source,Check,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeekDegMinSrc: dealloc #####"

end subroutine SeekDegMinSrc

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

subroutine SelectOriginal

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")

do
  print*, "  > Does the file have a parallel .src file (0=no,1=yes) ?"
  read (*,*,iostat=ReadStatus), QOrigSrc
  if (ReadStatus.LE.0.AND.QOrigSrc.GE.0.AND.QOrigSrc.LE.1) exit
end do

if (QOrigSrc.EQ.1) then
  do
   print*, "  > Identify the parallel .src file: "
   read (*,*,iostat=ReadStatus), GivenFile
   if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  OrigSrcFileL = LoadPath (GivenFile,".src")
else
  do
   print*, "  > Identify the source code for the file: "
   read (*,*,iostat=ReadStatus), OrigSrcCode
   if (ReadStatus.LE.0) exit
  end do
end if

do
  print*, "  > Is the original 7-digit code in the 'grid box' col? (0=no,1=yes)"
  read (*,*,iostat=ReadStatus), QOrigCode
  if (ReadStatus.LE.0.AND.QOrigCode.GE.0.AND.QOrigCode.LE.1) exit
end do

if (QOrigCode.EQ.0) then
 do
  print*, "  > Does the file have 7 (=7) or 6 (=6) digit station codes (i7) ?"
  read (*,*,iostat=ReadStatus), QDigits
  if (ReadStatus.LE.0.AND.QDigits.GE.6.AND.QDigits.LE.7) exit
 end do
end if

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
else if (VariableSuffix.EQ.".tmp") then
  DecayDistance = 1200 ; Min=-2730 ; Max=9000
else if (VariableSuffix.EQ.".dtr") then
  DecayDistance = 750  ; Min=0 ; Max=10000
else if (VariableSuffix.EQ.".cld") then
  DecayDistance = 750  ; Min=0 ; Max=1000
else if (VariableSuffix.EQ.".spc") then		! assumed, on the basis of cld
  DecayDistance = 750  ; Min=0 ; Max=1000
else if (VariableSuffix.EQ.".vap") then
  DecayDistance = 1200 ; Min=0 ; Max=10000
else if (VariableSuffix.EQ.".rhm") then		! assumed, on the basis of vap and tmp
  DecayDistance = 1200 ; Min=0 ; Max=1000
else if (VariableSuffix.EQ.".wet") then
  DecayDistance = 450  ; Min=0 
  Max=(/3100,2900,3100,3000,3100,3000,3100,3100,3000,3100,3000,3100/)
else
  print*, "  > @@@@@ ERROR: SelectVariable: no var @@@@@"
end if

end subroutine SelectVariable

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

end program CleanMeta