! opcruts.f90
! f90 program written by Tim Mitchell on 11.02.02
! last modified on 01.05.02
! program to manipulate CRU ts files (.cts .hdr)
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/opcruts filenames.f90 time.f90 grimfiles.f90 
!	crutsfiles.f90 saveperfiles.f90 annfiles.f90 cetgeneral.f90 basicfun.f90 
!	wmokey.f90 gridops.f90 grid.f90 ctyfiles.f90 
!	./../cruts/opcruts.f90 2> /tyn1/tim/scratch/stderr.txt

program OpCRUts

use FileNames
use Time
use GrimFiles
use CRUtsFiles
use SavePerFiles
use AnnFiles
use CETGeneral
use BasicFun
use WMOKey
use GridOps
use Grid
use CtyFiles

implicit none

real, pointer, dimension (:,:,:)		:: Stdev,Means,RealData
real, pointer, dimension (:,:)			:: StnMon,StnSea,SumsPer,NormMean,NormStdev
real, pointer, dimension (:,:)			:: RefGridLat,RefGridLon
real, pointer, dimension (:)			:: SrcEn,SrcTot,SrcTotSq
real, pointer, dimension (:)			:: StnAnn, Scatter, SumsAnn
real, pointer, dimension (:)			:: ALat,ALon,AElv,BLat,BLon,BElv,CLat,CLon,CElv

real, dimension (12)	:: ThreshLo,ThreshHi,MonthMeans,TheoryMax
real, dimension ( 4)	:: RefBounds

integer, pointer, dimension (:,:,:)		:: Data,DataA,DataB,DataC
integer, pointer, dimension (:,:,:)		:: DataSrc,DataASrc,DataBSrc,DataCSrc
integer, pointer, dimension (:,:,:)		:: LogMeans,LogStdev
integer, pointer, dimension (:,:)		:: DtbNormalsA,DtbNormalsC
integer, pointer, dimension (:,:)		:: HulmeA,Normals,RefGrid,MonthLengths
integer, pointer, dimension (:,:)		:: StnInfoA,StnInfoB,StnInfoC
integer, pointer, dimension (:,:)		:: LogPass,LogMiss,LogSkip,LogFail,LogAny,LogNorm
integer, pointer, dimension (:)			:: YearAD,StnYearAD,AYearAD,BYearAD,CYearAD,Order
integer, pointer, dimension (:)			:: CtyCode0,CtyCode1,CtyReg
integer, pointer, dimension (:)			:: AStn,BStn,CStn, AStnOld,BStnOld,CStnOld, MapCA,MapCB
integer, pointer, dimension (:)			:: StnSrcCodeA,StnSrcDateA
integer, pointer, dimension (:)			:: NmlYr0,NmlYr1,NmlSrc,NmlInc

character (len=80), pointer, dimension (:)	:: SaveBatch,DatabaseBatch
character (len=80), pointer, dimension (:)	:: Bat01,Bat02,Bat03,Bat04,Bat05,Bat06
character (len=80), pointer, dimension (:)	:: Bat07,Bat08,Bat09,Bat10,Bat11,Bat12
character (len=20), pointer, dimension (:)	:: StnNameA,StnNameB,StnNameC,AUni,Source
character (len=13), pointer, dimension (:)	:: StnCtyA,StnCtyB,StnCtyC,CtyName,CtyFinal
character (len=09), pointer, dimension (:)	:: StnLocalA,StnLocalB,StnLocalC,ColTitles
character (len=04), pointer, dimension (:)	:: StnSrcSuffixA

character (len=20), dimension (9) :: VersionKey=(/"WMO STATIONS MAY 98","WMO MAY 2002","NCDC MAY 2002", &
						"GHCN V2 TMP SERIES","HULME PRECIP","LISTER MAY 02", &
						"CRU NORMALS","CRU BIG.HDR.HYP","CRU TMP SERIES"/)
character (len=04), dimension (4) :: LogName=(/"PASS","MISS","SKIP","FAIL"/)
character (len=09), dimension(12) :: RegTitles=(/'    Globe','   Europe','  ex-USSR', &
						' Mid-East','     Asia','   Africa', &
	   					' NAmerica',' CAmerica',' SAmerica', &
	   					'   Antarc','  Oceania','   Marine'/)
character (len=02), dimension (12) :: MonthTxt=(/"01","02","03","04","05","06","07","08","09","10","11","12"/)
	   
real, dimension(12), parameter 	:: Ephemeris=(/-0.346,-0.216,-0.029, 0.173, 0.331, 0.403, &
					        0.369, 0.237, 0.044,-0.158,-0.316,-0.387/)
real, parameter 		:: MissVal = -999.0
integer, parameter 		:: DataMissVal = -9999
integer, parameter 		:: MultiSource=-888

character (len=80) :: LogMeansFile="/cru/mikeh1/f709762/scratch/log/log-opcruts-means.cts"
character (len=80) :: LogStdevFile="/cru/mikeh1/f709762/scratch/log/log-opcruts-stdev.cts"

logical :: ValidNorm

real :: Lat0,Lat1,Lon0,Lon1,Elv0,Elv1, PrecLat,PrecLon,PrecElv
real :: DataMulti,ThreshStdev,ThreshReject,ThreshPercent,ThreshOpCalc 
real :: EntryReal,Limit, PrevLat,PrevLon, Constant, SeqStatus, Sigma,Middl,Excess
real :: MissThresh,StdevThresh,DistanceThresh,Factor, ExeSpace,WyeSpace
real :: OpVal,OpTot,OpEn,OpTotSq,OpStdev,OpMean,OpDiff,OpMean0,OpMean1,OpStdev0,OpStdev1,OpStdevJ
real :: Distance,DecayDistance,Loaded,LatRad,Denom,Numer

integer :: ReadStatus,AllocStat
integer :: QMenu,QDataHeadLoad,QDataHeadSave,QTrueQuasi,QLongType,QCrossDL,QFile,QNoStdev,QNoPer
integer :: QDiff,QFake,QAlter,QManual,QReverse,QShift,QType,QSource,QOperation,QRestrict,QCombine
integer :: QManCheck,QReject,QNRMfile,QOutput,QSignCheck,QMisMatch,QStoreOld,QHomogSource,QMerge
integer :: QAnomPercent,QPosOnly
integer :: NAStn,NBStn,NCStn,NInfo,NStnYear,NAYear,NBYear,NCYear,NMonth,NVersion,NCty,NReg,NSrc
integer :: XAStn,XBStn,XCStn,XInfo,XStnYear,XAYear,XBYear,XCYear,XMonth,XVersion,XCty,XReg,XSrc
integer :: NBatch,NExe,NWye,NFile
integer :: XBatch,XExe,XWye,XFile
integer :: XAYear0,XLog,XSeason,XSeasonMonth,XSeasonYear,XMonth0,XYear0,XMonth1,XYear1,XBMonth
integer :: AddYearAD0,AddYearAD1,PrevStn,EquiStn,RegSub0,StatType,SourceCode,SuffixBeg, Digits
integer :: YearAD0,YearAD1,Month0,Month1, AStart,AEnd,BStart,BEnd,CAStart,CAEnd,CBStart,CBEnd
integer :: IntLat,IntLon,IntElv, LatMV,LonMV,ElvMV,DataMV, LatF,LonF,ElvF
integer :: FakeStn,OrigStn,NewStn, TotA,TotB,TotAB, Min,Max,MinImpose,MaxImpose
integer :: Log01,Log02,Log03,Log04,Log05,Log06,Log07,Log08,Log09,Log11,Log12,Log13
integer :: Log0,Log1,Log2,Log3,Log4,Log5,Log6,Log7
integer :: InRange, CheckYear0,CheckYear1,SaveYearAD0,SaveYearAD1, NormMin,StrLen,SubBeg
integer :: RefBoxes,SrcRemove

character (len=80) :: GivenFile,LoadFileA,LoadFileB,LoadFileC,SaveFileA,SaveFileB,SaveFileC,HeadFormat
character (len=80) :: LoadFileASrc,LoadFileBSrc,SaveFileASrc,SaveFileBSrc,SaveFileCSrc,Stem,Tip
character (len=80) :: CurrStr,PrevStr,Variable,Title,Trash
character (len=80) :: DatabaseFilter,DatabaseFilterSrc,DatabaseFileL,DatabaseSrcFileL
character (len=20) :: SourceStr,EntryString,YearTxt,LineFormat
character (len=08) :: Date
character (len=04) :: LoadSuffix,SaveSuffix,SaveFileType,FileSuffix,TimeVal

call Intro
call Menu
if (QMenu.EQ. 1) call ConvertToHead
if (QMenu.EQ. 2) call ConvertToPerSet
if (QMenu.EQ. 3) call GenerateSrcFile
if (QMenu.EQ. 4) call ConvertNewCode
if (QMenu.EQ. 5) call AddSourceHdr
if (QMenu.EQ. 6) call ConvertWMO
if (QMenu.EQ. 7) call ConvertLegacy
if (QMenu.EQ.10) call CorrectError
if (QMenu.EQ.11) call IDaddAB
if (QMenu.EQ.12) call CombineAB
if (QMenu.EQ.13) call DescribeAPer
if (QMenu.EQ.14) call DescribeAAnn
if (QMenu.EQ.18) call OperateAk
if (QMenu.EQ.19) call OperateAB
if (QMenu.EQ.21) call GetLatLongRect
if (QMenu.EQ.22) call TruncateDate
if (QMenu.EQ.23) call QCImplicitStdev
if (QMenu.EQ.25) call Anomalise
if (QMenu.EQ.26) call FindClosest
if (QMenu.EQ.27) call MultiplyCodes
if (QMenu.EQ.28) call SyncCtsSrc
if (QMenu.EQ.29) call ImpossibleElev
if (QMenu.EQ.30) call ConvertDegHun
if (QMenu.EQ.31) call MultiplyNRM
if (QMenu.EQ.34) call PruneSource
if (QMenu.EQ.40) call ConvertVariableDtb
if (QMenu.EQ.41) call ConvertVariableCts
call Finish

contains

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

subroutine Intro

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

print*
print*, "  > ***** OpCRUts: manipulates .cts .hdr .src .dtb .dts files *****"
print*

NMonth = 12 ; NVersion = 9 ; NSrc = 1000

end subroutine Intro

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

subroutine Menu

print*, "  > Select an operation to perform (0=list): "
do
	read (*,*,iostat=ReadStatus), QMenu
	
	if (QMenu.EQ.0) then
	  print*, "  >  1. Convert into a CRU ts file."
	  print*, "  >  2. Convert into a set of .per files."
	  print*, "  >  3. Generate a parallel .src file."
	  print*, "  >  4. Convert from Mark New source code file"
	  print*, "  >  5. Add source info to .hdr file"
	  print*, "  >  6. Convert WMO data file to .hdr file"
	  print*, "  >  7. Convert legacy .cts to .cts file"
	  print*, "  > 11. Identify additionality of B over A."
	  print*, "  > 12. Combine A (priority) and B."
	  print*, "  > 13. Describe A by month     (output .per file)."
	  print*, "  > 14. Describe A by continent (output .ann file)."
	  print*, "  > 18. Operate A.k"
	  print*, "  > 19. Operate A.B"
	  print*, "  > 21. Extract all data within lat/long rectangle."
	  print*, "  > 22. Truncate by date."
	  print*, "  > 23. Quality control using implicit stdev."
	  print*, "  > 25. Anomalise .cts (now anomdtb.f90)" 
	  print*, "  > 26. Find closest stns"
	  print*, "  > 27. Multiply station codes."
	  print*, "  > 28. Check .cts and .src are in sync"
	  print*, "  > 29. Remove impossible elevations."
	  print*, "  > 30. Convert deg+minutes to deg+hundredths"
	  print*, "  > 31. Multiply .nrm"
	  print*, "  > 34. Remove data from a selected source"
	  print*, "  > 40. Convert variable .dtb"
	  print*, "  > 41. Convert variable .cts"
	end if
	
	if (ReadStatus.LE.0.AND.QMenu.GE.1) exit
end do    

end subroutine Menu

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

subroutine ConvertToHead

print*, "  > Load a true (=1), quasi (=2), MikeH (=3), PhilJ (=4) file?"
do
	read (*,*,iostat=ReadStatus), QTrueQuasi
	if (ReadStatus.LE.0.AND.QTrueQuasi.GE.1.AND.QTrueQuasi.LE.4) exit
end do

if (QTrueQuasi.EQ.1) then
  print*, "  > Load a .src (=0) .cts (=1) .hdr (=2) .dtb (=3) .dts (=4) file ?"
  do
	read (*,*,iostat=ReadStatus), QDataHeadLoad
	if (ReadStatus.LE.0.AND.QDataHeadLoad.GE.0.AND.QDataHeadLoad.LE.4) exit
  end do
else if (QTrueQuasi.EQ.2) then
  print*, "  > Load a data (=1) or header (=2) file ?"
  do
	read (*,*,iostat=ReadStatus), QDataHeadLoad
	if (ReadStatus.LE.0.AND.QDataHeadLoad.GE.1.AND.QDataHeadLoad.LE.2) exit
  end do
else
  QDataHeadLoad = 1
end if

if      (QTrueQuasi.EQ.1) then
  	if (QDataHeadLoad.EQ.0) LoadSuffix = ".src"
  	if (QDataHeadLoad.EQ.1) LoadSuffix = ".cts"
  	if (QDataHeadLoad.EQ.2) LoadSuffix = ".hdr"
  	if (QDataHeadLoad.EQ.3) LoadSuffix = ".dtb"
  	if (QDataHeadLoad.EQ.4) LoadSuffix = ".dts"
else if (QTrueQuasi.EQ.4) then
  	LoadSuffix = "    "
  	HeadFormat = "(i6,i4,i5,i5,x,a20,x,a13,x,4i4)"
  	QLongType  = 2
  	LatF=10 ; LonF=10 ; ElvF=1 ; DataMV=-999
  	LatMV=-999 ; LonMV=-1999 ; ElvMV=-999
else
  	LoadSuffix = "    "
  	
  	if (QTrueQuasi.EQ.2) then
	  print*, "  > Enter the header line format (i?,i?,i?,i?,a??,a??,i?,i?,i?,a?): "
	  do
		read (*,*,iostat=ReadStatus), HeadFormat
		if (ReadStatus.LE.0.AND.HeadFormat.NE."") exit
	  end do

	  print*, "  > Are the longs 0:180:-180:0 (=1), 0:-180:180:0 (=2), 0:360 (=3): "
	  do
		read (*,*,iostat=ReadStatus), QLongType
		if (ReadStatus.LE.0.AND.QLongType.GE.1.AND.QLongType.LE.3) exit
	  end do
	  
	  print*, "  > Enter the lat,lon,elv factors (integers):"
	  do
		read (*,*,iostat=ReadStatus), LatF,LonF,ElvF
		if (ReadStatus.LE.0) exit
	  end do

	  print*, "  > Enter the lat,lon,elv missing values (integers):"
	  do
		read (*,*,iostat=ReadStatus), LatMV,LonMV,ElvMV
		if (ReadStatus.LE.0) exit
	  end do

  	  if (QDataHeadLoad.EQ.1) then
  	    print*, "  > Enter the data missing value in the file: "
	    do
		read (*,*,iostat=ReadStatus), DataMV
		if (ReadStatus.LE.0) exit
	    end do
  	  end if
  	end if
end if
print*, "  > Select the file to load, with suffix: ", LoadSuffix
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,LoadSuffix)

print*, "  > Save as .src (=0) .cts (=1) .hdr (=2) .dtb (=3) .dts (=4) file ?"
do
	read (*,*,iostat=ReadStatus), QDataHeadSave
	if (ReadStatus.LE.0.AND.QDataHeadSave.GE.0.AND.QDataHeadSave.LE.4) exit
end do

if      (QDataHeadSave.EQ.0) then
  SaveSuffix = ".src" ; QDataHeadSave = 1
else if (QDataHeadSave.EQ.1) then
  SaveSuffix = ".cts"
else if (QDataHeadSave.EQ.2) then
  SaveSuffix = ".hdr"
else if (QDataHeadSave.EQ.3) then
  SaveSuffix = ".dtb"
else if (QDataHeadSave.EQ.4) then
  SaveSuffix = ".dts"
end if

print*, "  > Select the file to save, with suffix: ", SaveSuffix
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,SaveSuffix)

print*, "  > Loading..."
if      (QTrueQuasi.EQ.1) then
  if (QDataHeadLoad.EQ.0.OR.QDataHeadLoad.EQ.1) then
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA)
  else if (QDataHeadLoad.EQ.2) then
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			CallFile=LoadFileA,HeadOnly=1)
  else if (QDataHeadLoad.EQ.3.OR.QDataHeadLoad.EQ.4) then
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,DtbNormals=DtbNormalsA,YearAD=YearAD,CallFile=LoadFileA)
  end if
else if (QTrueQuasi.EQ.2) then
  if (QDataHeadLoad.EQ.1) then
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA,HeadForm=HeadFormat, &
			LatF=LatF,LonF=LonF,ElvF=ElvF, &
			LongType=QLongType,LatMV=LatMV,LonMV=LonMV,ElvMV=ElvMV,DataMV=DataMV)
  else
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			CallFile=LoadFileA,HeadOnly=1,HeadForm=HeadFormat, &
			LatF=LatF,LonF=LonF,ElvF=ElvF, &
			LongType=QLongType,LatMV=LatMV,LonMV=LonMV,ElvMV=ElvMV)
  end if
else if (QTrueQuasi.EQ.3) then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA,Hulme=HulmeA)
else if (QTrueQuasi.EQ.4) then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA,HeadForm=HeadFormat, &
			LatF=LatF,LonF=LonF,ElvF=ElvF,PhilJ=1, &
			LongType=QLongType,LatMV=LatMV,LonMV=LonMV,ElvMV=ElvMV,DataMV=DataMV)
end if
NAStn = size(AStn,1)				

print*, "  > Saving..."

if (QDataHeadLoad.LE.1) then		! YearAD and DataA already exist
  QNoPer = 1					! only save stations with genuine data
else if (QDataHeadLoad.EQ.2) then	! YearAD and DataA do not already exist
  allocate (DataA    (1,NMonth,NAStn), &
	    YearAD   (1), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToHead: Allocation failure: DataA #####"
  DataA=-9999 ; YearAD=-999 
  QNoPer = 0					! save all stns regardless of missing data

else if (QDataHeadLoad.GT.2) then	! YearAD and DataA already exist
  QNoPer = 1
end if

if (QDataHeadLoad.LE.2.AND.QDataHeadSave.GT.2) then
  allocate (DtbNormalsA(NMonth,NAStn), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToHead: Allocation failure: DtbNA #####"
  DtbNormalsA=-9999
end if

if (QDataHeadSave.LE.2) then
  call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,QNoPer,YearAD=YearAD,Data=DataA)
else
  call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,QNoPer,YearAD=YearAD,Data=DataA, &
  			DtbNormals=DtbNormalsA)
end if

allocate (Order(NAStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToHead: Alloc: Order #####"
do XAStn=1,NAStn
    Order(XAStn)=XAStn
end do

if (QDataHeadSave.LE.1) then
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD, &
  		Order=Order,CallFile=SaveFileA)
else if (QDataHeadSave.EQ.2) then
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,HeadOnly=1, &
  		Order=Order,CallFile=SaveFileA)
else
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD, &
  		CallFile=SaveFileA,Order=Order, DtbNormals=DtbNormalsA)
end if

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

end subroutine ConvertToHead

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

subroutine ConvertToPerSet

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Enter the multiplier to apply to the data in the file:"
do
	read (*,*,iostat=ReadStatus), DataMulti
	if (ReadStatus.LE.0) exit
end do

print*, "  > Enter the statistic (-1=min,0=mean,1=max,2=sum): "
do
	read (*,*,iostat=ReadStatus), StatType
	if (ReadStatus.LE.0.AND.StatType.GE.-1.AND.StatType.LE.2) exit
end do

print*, "  > Select the generic .per file to save (include 'regcode'):"
do
	read (*,*,iostat=ReadStatus), GivenFile
	
	if (ReadStatus.LE.0.AND.GivenFile.NE."") then
	  RegSub0 = index(GivenFile,"regcode")
	else
	  RegSub0 = 0
	end if
	
	if (RegSub0.GT.0) exit
end do

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

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=Data,YearAD=YearAD,CallFile=LoadFileA)

NAStn = size(StnInfoA,1)

allocate (AUni(NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToPerSet: Allocation failure: A* #####"
do XAStn = 1, NAStn
  AUni(XAStn) = GetTextFromInt(StnInfoA(XAStn,1))
end do

Stem = GivenFile(01:(RegSub0-1))
Tip  = GivenFile((RegSub0+7):80)
call MakeBatch (Stem,Tip,AUni,SaveBatch)

do XAStn = 1, NAStn
  NStnYear = StnInfoA(XAStn,6) - StnInfoA(XAStn,5) + 1		! no. of years in stn record
  
  allocate (StnMon   (NStnYear,12), &				! arrays for saving to .per
  	    StnSea   (NStnYear, 4), &
  	    StnAnn   (NStnYear),    &
  	    StnYearAD(NStnYear),    stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToPerSet: Allocation failure: Stn* #####"
  StnMon=MissVal ; StnSea=MissVal ; StnAnn=MissVal		! initialise arrays
  
  XAYear0 = 0							! find index of first stn year in Data
  do
    XAYear0 = XAYear0 + 1
    if (YearAD(XAYear0).EQ.StnInfoA(XAStn,5)) exit
  end do
  
  XAYear = XAYear0 - 1
  do XStnYear = 1, NStnYear					! iterate by stn year
    XAYear = XAYear + 1						! identify corresponding year in Data
    
    StnYearAD(XStnYear) = StnInfoA(XAStn,5) + XStnYear - 1	! fill stn YearAD
    
    do XMonth = 1, 12						! fill stn Monthly data file
      if (Data(XAYear,XMonth,XAStn).NE.DataMissVal) &
      		StnMon(XStnYear,XMonth) = Data(XAYear,XMonth,XAStn) * DataMulti
    end do
  end do
  
  if (StatType.EQ.-1) call FillSeaAnnMin  (StnYearAD,StnMon,StnSea,StnAnn)
  if (StatType.EQ. 0) call FillSeaAnnMean (StnYearAD,StnMon,StnSea,StnAnn)
  if (StatType.EQ. 1) call FillSeaAnnMax  (StnYearAD,StnMon,StnSea,StnAnn)
  if (StatType.EQ. 2) call FillSeaAnnSum  (StnYearAD,StnMon,StnSea,StnAnn)
  
  call SavePER (SaveBatch(XAStn),StnYearAD,StatType,Monthly=StnMon,Seasonal=StnSea,Annual=StnAnn,NoResponse=1)
    
  deallocate (StnMon,StnSea,StnAnn,StnYearAD,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToPerSet: Deallocation failure: Stn* #####"
end do

end subroutine ConvertToPerSet

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

subroutine GenerateSrcFile

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")
SuffixBeg = index(LoadFileA,".cts")
SaveFileB = LoadFileA
SaveFileB(SuffixBeg:(SuffixBeg+3)) = ".src"

print*, "  > Enter the source code:"
do
	read (*,*,iostat=ReadStatus), SourceCode
	if (ReadStatus.LE.0) exit
end do

print*, "  > Loading from .cts and saving to .src ..."

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD,CallFile=LoadFileA,Silent=1)

NAYear=size(DataA,1) ; NAStn=size(DataA,3)

allocate (DataB(NAYear,NMonth,NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertToPerSet: Allocation failure: A* #####"
DataB = -9999

do XAYear = 1, NAYear
  do XMonth = 1, NMonth
    do XAStn = 1, NAStn
      if (DataA(XAYear,XMonth,XAStn).NE.-9999) DataB(XAYear,XMonth,XAStn) = SourceCode
    end do
  end do
end do

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataB,YearAD=YearAD,CallFile=SaveFileB,Silent=1)

end subroutine GenerateSrcFile

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

subroutine ConvertNewCode

print*, "  > Select the Mark New source code master file:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,"    ")

print*, "  > Select the .hdr file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".hdr")

call system ('wc -l ' // LoadFileA // ' > trashme-loadcts.txt')		! get number of lines
open  (3,file='trashme-loadcts.txt',status="old",access="sequential",form="formatted",action="read")
read  (3,"(i10)"), NAStn
close (3)
call system ('rm trashme-loadcts.txt')

print*, "  > Stations in file total: ", NAStn

allocate (AStn     (NAStn), &
          AStnOld  (NAStn), &
	  StnNameA (NAStn), &
	  StnLocalA(NAStn), &
	  ALat     (NAStn), &
	  ALon     (NAStn), &
	  AElv     (NAStn), &
	  StnCtyA  (NAStn), &
	  Source   (NAStn), &
	  YearAD   (1),     &
	  DataA    (1,NMonth,NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertNewCode: Allocation failure #####"
LatMV=-9999 ; LonMV=-99999 ; ElvMV=-999 ; AStnOld=-999
YearAD=-999 ; DataA=-9999 ; ALat=MissVal ; ALon=MissVal ; AElv=MissVal

print*, "  > Loading... "

open  (2,file=LoadFileA,status="old",access="sequential",form="formatted",action="read")    
do XAStn = 1, NAStn
  read (2,"(i7,x,a20,4x,a9,3i6,x,a13,a19)"), AStn(XAStn),StnNameA(XAStn), &
  		StnLocalA(XAStn),IntLat,IntLon,IntElv,StnCtyA(XAStn),Source(XAStn)

  if (IntLon.EQ.-9999.AND.IntLat.EQ.-9999) IntLon=LonMV
  if (IntLon.EQ.19990) IntLon=LonMV
  if (IntLat.LE.-9990) IntLat=LatMV
  if (IntElv.LE. -999) IntElv=ElvMV
  
!  if (AStn(XAStn).EQ.9739017) print "(3i9)", IntLat,IntLon,IntElv	! @@@@@@@@@@@@@@@

  if (IntLat.NE.LatMV) ALat(XAStn)=real(IntLat)/100.0
  if (IntLon.NE.LonMV) ALon(XAStn)=real(IntLon)/100.0
  if (IntElv.NE.ElvMV) AElv(XAStn)=real(IntElv)
  
!  if (AStn(XAStn).EQ.9739017) print "(3f9.2)", ALat(XAStn),ALon(XAStn),AElv(XAStn)	! @@@@@@@@@@@@@@@
end do
close (2)

print*, "  > Saving... "

call MakeStnInfo (StnInfoA,AStn,AStnOld,ALat,ALon,AElv,0,Data=DataA,YearAD=YearAD)

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileB,Silent=1,HeadOnly=1,Source=Source)

end subroutine ConvertNewCode

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

subroutine AddSourceHdr

print*, "  > Select the .hdr file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".hdr")

print*, "  > Enter the source string (up to 20 char): "
do
	read (*,*,iostat=ReadStatus), SourceStr
	if (ReadStatus.LE.0.AND.SourceStr.NE."") exit
end do

print*, "  > Select the .hdr file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".hdr")

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=LoadFileA,HeadOnly=1)

NAStn = size(StnInfoA,1)
allocate (Source(NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AddSourceHdr: Allocation failure #####"
Source = trim(SourceStr)

print*, "  > Saving..."
call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileB,HeadOnly=1,Source=Source)

end subroutine AddSourceHdr

!*******************************************************************************
! take a (processed?) WMO code file and convert to .hdr with 7-digit codes

subroutine ConvertWMO

print*, "  > Select the WMO code file:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,"    ")

print*, "  > Select the .hdr file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".hdr")

print*, "  > Operating..."

call system ('wc -l ' // LoadFileA // ' > trashme-loadcts.txt')		! get number of lines
open  (3,file='trashme-loadcts.txt',status="old",access="sequential",form="formatted",action="read")
read  (3,"(i10)"), NAStn
close (3)
call system ('rm trashme-loadcts.txt')

allocate (AStn     (NAStn), &
	  AStnOld  (NAStn), &
	  StnLocalA(NAStn), &
	  DataA    (1,NMonth,NAStn), &
	  YearAD   (1), &
	  Source   (NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertNewCode: Allocation failure #####"
StnLocalA="   nocode" ; DataA=-9999 ; YearAD=-999 ; Source="NCDC MAY 2002" ; AStnOld=-999

open  (2,file=LoadFileA,status="old",access="sequential",form="formatted",action="read")    
do XAStn = 1, NAStn
  read (2,"(i6)"), AStn(XAStn)
end do
close (2)

call GetCRUtsInfo (AStn,ALat,ALon,AElv,StnNameA,StnCtyA,6)

AStn = AStn * 10			! convert from 6-digit code to 7-digit code

call MakeStnInfo (StnInfoA,AStn,AStnOld,ALat,ALon,AElv,0,Data=DataA,YearAD=YearAD)

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileB,Silent=1,HeadOnly=1,Source=Source)

end subroutine ConvertWMO

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

subroutine ConvertLegacy

print*, "  > Select the legacy .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Enter the no. digits in the WMO code:"
do
	read (*,*,iostat=ReadStatus), Digits
	if (ReadStatus.LE.0.AND.Digits.GE.5.AND.Digits.LE.7) exit
end do

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".cts")

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA,Legacy=1)

call GetCRUtsInfo (AStn,BLat,BLon,BElv,StnNameB,StnCtyB,Digits)

NAStn = size(AStn,1)

do XAStn = 1, NAStn
  if      (ALat(XAStn).NE.MissVal.AND.BLat(XAStn).NE.MissVal .AND. &
  	   ALon(XAStn).NE.MissVal.AND.BLon(XAStn).NE.MissVal) then
      if (nint(ALat(XAStn)*10.0).EQ.nint(BLat(XAStn)*10.0) .AND. &
      	  nint(ALon(XAStn)*10.0).EQ.nint(BLon(XAStn)*10.0)) then
        ALat(XAStn)=BLat(XAStn) ; ALon(XAStn)=BLon(XAStn)
      end if
  end if
end do

call MakeStnInfo (StnInfoA,AStn,AStnOld,ALat,ALon,AElv,1,Data=DataA,YearAD=YearAD)

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileB,Data=DataA,YearAD=YearAD,Silent=1)

end subroutine ConvertLegacy

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

subroutine CorrectError

print*, "  > Select the .hdr file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".hdr")

print*, "  > Select the .hdr file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".hdr")

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			HeadOnly=1,CallFile=LoadFileA,Silent=1,Source=Source)
NAStn = size(AStn,1)
do XAStn = 1, NAStn
  if (StnInfoA(XAStn,7).EQ.-999) AStn(XAStn)=-999.0
  if (StnInfoA(XAStn,7).GT.0.AND.StnInfoA(XAStn,7).LT.10000) AStn(XAStn)=-999.0
end do

allocate (YearAD (1), &
	  DataA  (1,NMonth,NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CorrectError: Allocation failure #####"
DataA = -9999 ; YearAD = -999

print*, "  > Saving..."

call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,0,YearAD=YearAD,Data=DataA)

call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileB, &
		HeadOnly=1,Source=Source)

end subroutine CorrectError

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

subroutine IDaddAB

print*, "  > Select the A (master) header file (.hdr) to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".hdr")

print*, "  > Select the B (extra)  header file (.hdr) to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileB = LoadPath (GivenFile,".hdr")

print*, "  > Select the header file (.hdr) to save."
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".hdr")

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=LoadFileA,HeadOnly=1)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,CallFile=LoadFileB,HeadOnly=1)

NAStn = size (StnNameA,1) ; NBStn = size (StnNameB,1) ; NCStn = NBStn

allocate (StnInfoC (NCStn,8), &
	  StnLocalC(NCStn),   &
	  StnNameC (NCStn),   &
	  StnCtyC  (NCStn),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: IDaddAB: Allocation failure: Data #####"
StnInfoC = MissVal ; StnNameC = "" ; StnCtyC = "" ; StnLocalC="   nocode"

PrevStn = 1
do XBStn = 1, NBStn
  if (StnInfoB(XBStn,1).GT.0) then
    XAStn = PrevStn-1 ; EquiStn = MissVal ; AddYearAD0 = MissVal ; AddYearAD1 = MissVal
    do      
      XAStn = XAStn + 1
      if (StnInfoA(XAStn,1).LT.StnInfoB(XBStn,1)) PrevStn = XAStn
      if (StnInfoA(XAStn,1).EQ.StnInfoB(XBStn,1)) EquiStn = XAStn
      if (StnInfoA(XAStn,1).GT.StnInfoB(XBStn,1)) exit
      if (XAStn.EQ.NAStn) exit
    end do
    
    if (EquiStn.EQ.MissVal) then
      AddYearAD0 = StnInfoB(XBStn,5) ; AddYearAD1 = StnInfoB(XBStn,6)
    else
      if (StnInfoA(XAStn,5).GT.StnInfoB(XBStn,5)) then
        AddYearAD0 = StnInfoB(XBStn,5)
        if (StnInfoA(XAStn,6).LT.StnInfoB(XBStn,6)) then
          AddYearAD1 = StnInfoB(XBStn,6)
        else
          if (StnInfoA(XAStn,5).GT.StnInfoB(XBStn,6)) then
            AddYearAD1 = StnInfoB(XBStn,6)
          else
            AddYearAD1 = StnInfoA(XAStn,5) - 1
          end if
        end if
      else if (StnInfoA(XAStn,6).LT.StnInfoB(XBStn,6)) then
        AddYearAD1 = StnInfoB(XBStn,6)
        if (StnInfoA(XAStn,6).LT.StnInfoB(XBStn,5)) then
            AddYearAD0 = StnInfoB(XBStn,5)
        else
            AddYearAD0 = StnInfoA(XAStn,6) + 1
        end if
      end if
    end if
    
    if (AddYearAD0.NE.MissVal) then
      StnLocalC(XBStn) = StnLocalB(XBStn)
      StnNameC (XBStn) = StnNameB (XBStn)
      StnCtyC  (XBStn) = StnCtyB  (XBStn)
      do XInfo = 1, 8
        StnInfoC(XBStn,XInfo) = StnInfoB(XBStn,XInfo)
      end do
      StnInfoC(XBStn,5)=AddYearAD0 ; StnInfoC(XBStn,6)=AddYearAD1
    end if
  end if
end do

print*, "  > Saving..."
call SaveCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,CallFile=SaveFileC,HeadOnly=1)

end subroutine IDaddAB

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

subroutine CombineAB

print*, "  > Combine: C=A+B (=0), C=A where B has data (=1) ?"
do
	read (*,*,iostat=ReadStatus), QCombine
	if (ReadStatus.LE.0.AND.QCombine.GE.0.AND.QCombine.LE.1) exit
end do

print*, "  > Select .cts file A (priority) to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select .cts file B to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileB = LoadPath (GivenFile,".cts")

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".cts")
SuffixBeg = index(LoadFileB,".cts")
SaveFileB = LoadFileB ; SaveFileB(SuffixBeg:SuffixBeg+3) = ".ann"

print*, "  > Process parallel source files (0=no,1=yes) ?"
do
	read (*,*,iostat=ReadStatus), QSource
	if (ReadStatus.LE.0.AND.QSource.GE.0.AND.QSource.LE.1) exit
end do
if (QSource.EQ.1) then
  SuffixBeg=index(LoadFileA,".cts") ; LoadFileASrc=LoadFileA ; LoadFileASrc(SuffixBeg:SuffixBeg+3)=".src"
  SuffixBeg=index(LoadFileB,".cts") ; LoadFileBSrc=LoadFileB ; LoadFileBSrc(SuffixBeg:SuffixBeg+3)=".src"
  SuffixBeg=index(SaveFileC,".cts") ; SaveFileCSrc=SaveFileC ; SaveFileCSrc(SuffixBeg:SuffixBeg+3)=".src"
end if

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=AYearAD,CallFile=LoadFileA,Silent=1)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Code=BStn,OldCode=BStnOld,Lat=BLat,Lon=BLon,Elv=BElv, &
			Data=DataB,YearAD=BYearAD,CallFile=LoadFileB,Silent=1)

call AllocateC

print*, "  > Combining stations..."
XAStn=1 ; XBStn = 1 ; XCStn = 0
do						! combine at the station level
  XCStn = XCStn + 1
  
  if (XAStn.LE.NAStn.AND.XBStn.LE.NBStn) then
    if (AStn(XAStn).EQ.BStn(XBStn)) QFile = 3
    if (AStn(XAStn).LT.BStn(XBStn)) QFile = 1
    if (AStn(XAStn).GT.BStn(XBStn)) QFile = 2
  else if (XAStn.LE.NAStn) then
    QFile = 1
  else if (XBStn.LE.NBStn) then
    QFile = 2
  else
    QFile = 0
  end if
  
  if      (QFile.EQ.3) then
    CStn(XCStn)=AStn(XAStn) ; CStnOld(XCStn)=AStnOld(XAStn) ; CLat(XCStn)=ALat(XAStn) ; CLon(XCStn)=ALon(XAStn) ; CElv(XCStn)=AElv(XAStn)    
    StnNameC(XCStn)=StnNameA(XAStn) ; StnCtyC(XCStn)=StnCtyA(XAStn)
    if (len_trim(StnLocalA(XAStn)).EQ.0) StnLocalC(XCStn)=StnLocalB(XBStn)
    if (len_trim(StnLocalC(XCStn)).EQ.0) StnLocalC(XCStn)="   nocode"
    if (len_trim(StnNameA(XAStn)) .EQ.0) StnNameC(XCStn)=StnNameB(XBStn)
    if (len_trim(StnCtyA (XAStn)) .EQ.0) StnCtyC (XCStn)=StnCtyB (XBStn)
    if (AElv(XAStn).EQ.-999.AND.BElv(XBStn).NE.-999) CElv(XCStn) = BElv(XBStn)    
    MapCA(XCStn) = XAStn ; MapCB(XCStn) = XBStn    
    XAStn = XAStn + 1 ; XBStn = XBStn + 1
  else if (QFile.EQ.1) then
    CStn(XCStn)=AStn(XAStn) ; CStnOld(XCStn)=AStnOld(XAStn) ; CLat(XCStn)=ALat(XAStn) ; CLon(XCStn)=ALon(XAStn) ; CElv(XCStn)=AElv(XAStn)
    StnNameC(XCStn)=StnNameA(XAStn) ; StnCtyC(XCStn)=StnCtyA(XAStn)
    if (len_trim(StnLocalA(XAStn)).NE.0) StnLocalC(XCStn)=StnLocalA(XAStn)
    MapCA(XCStn) = XAStn    
    XAStn = XAStn + 1
  else if (QFile.EQ.2) then
    CStn(XCStn)=BStn(XBStn) ; CStnOld(XCStn)=BStnOld(XBStn) ; CLat(XCStn)=BLat(XBStn) ; CLon(XCStn)=BLon(XBStn) ; CElv(XCStn)=BElv(XBStn)
    StnNameC(XCStn)=StnNameB(XBStn) ; StnCtyC(XCStn)=StnCtyB(XBStn)
    if (len_trim(StnLocalB(XBStn)).NE.0) StnLocalC(XCStn)=StnLocalB(XBStn)
    MapCB(XCStn) = XBStn    
    XBStn = XBStn + 1
  end if
  
  if (QFile.EQ.0) print "(a,i6)", "   > Total stations available: ", (XCStn-1)
  if (QFile.EQ.0) exit
end do

allocate (CYearAD  (NCYear),&
          DataC    (NCYear,NMonth,NCStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineAB: Allocation failure: CYearAD,DataC #####"
DataC = -9999
do XCYear = 1, NCYear
  CYearAD(XCYear) = YearAD0 + XCYear - 1
end do

call CommonVecPer (AYearAD,CYearAD,AStart,AEnd,CAStart,CAEnd)
call CommonVecPer (BYearAD,CYearAD,BStart,BEnd,CBStart,CBEnd)
print "(a,i4,a1,i4)", "   > Period from file A: ", AYearAD(AStart),"-",AYearAD(AEnd)
print "(a,i4,a1,i4)", "   > Period from file B: ", BYearAD(BStart),"-",BYearAD(BEnd)
print "(a,i4,a1,i4)", "   > Combined period:    ", CYearAD(1),"-",CYearAD(NCYear)

if (QCombine.EQ.0) call CombineTheTwoAdd
if (QCombine.EQ.1) call CombineTheTwoIfData
  
if (QCombine.EQ.0) print "(a,3i10)", "   > Data unique to A, B, common: ", TotA,TotB,TotAB
if (QCombine.EQ.1) print "(a,3i10)", "   > Okay, no-B-stn, no-B-data:   ", TotAB,TotA,TotB

print*, "  > Saving combined .cts file..."
call MakeStnInfo (StnInfoC,CStn,CStnOld,CLat,CLon,CElv,1,YearAD=CYearAD,Data=DataC)	! save .cts file
call SaveCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,Data=DataC,YearAD=CYearAD,CallFile=SaveFileC,Silent=1)

deallocate (DataA,DataB,AYearAD,BYearAD,StnInfoA,StnLocalA,StnNameA,StnCtyA,StnInfoB,StnLocalB,StnNameB,StnCtyB, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineAB: Deallocation failure: AB part 1 #####"

if (QSource.EQ.1) then
  DataC = -9999
  
  print*, "  > Loading parallel source files..."
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,YearAD=AYearAD,Data=DataA,CallFile=LoadFileASrc,Silent=1)
  call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,YearAD=BYearAD,Data=DataB,CallFile=LoadFileBSrc,Silent=1)
  
  deallocate (StnInfoA,StnLocalA,StnNameA,StnCtyA,StnInfoB,StnLocalB,StnNameB,StnCtyB, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineAB: Deallocation failure: AB part 2 #####"

  if (QCombine.EQ.0) call CombineTheTwoAdd
  if (QCombine.EQ.1) call CombineTheTwoIfData
  
  print*, "  > Saving parallel source file..."
  call SaveCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,Data=DataC,YearAD=CYearAD,CallFile=SaveFileCSrc,Silent=1)
end if

end subroutine CombineAB

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

subroutine AllocateC

NAYear = size(DataA,1) ; NAStn = size(DataA,3)
NBYear = size(DataB,1) ; NBStn = size(DataB,3)
print "(a,i6)", "   > Stations from file A:     ", NAStn
print "(a,i6)", "   > Stations from file B:     ", NBStn

YearAD0 = AYearAD(1) ; if (YearAD0.GT.BYearAD(1)) YearAD0 = BYearAD(1)
YearAD1 = AYearAD(NAYear) ; if (YearAD1.LT.BYearAD(NBYear)) YearAD1 = BYearAD(NBYear)
NCYear = YearAD1 - YearAD0 + 1 ; NCStn = NAStn + NBStn

allocate (CStn     (NCStn), &			! make station-specific allocations
	  CStnOld  (NCStn), &
	  CLat     (NCStn), &
	  CLon     (NCStn), &
	  CElv     (NCStn), &
	  StnNameC (NCStn), &
	  StnLocalC(NCStn), &
	  StnCtyC  (NCStn), &
	  MapCA    (NCStn), &
	  MapCB    (NCStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocateC: Allocation failure: CStn #####"
CStn=-999 ; CLat=-999 ; CLon=-9999 ; CElv=-999 ; MapCA=-999 ; MapCB=-999 
StnNameC="" ; StnCtyC="" ; StnLocalC="   nocode"

end subroutine AllocateC

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

subroutine CombineTheTwoAdd

TotA=0 ; TotB=0 ; TotAB=0

do XCStn = 1, NCStn
  if (MapCA(XCStn).NE.-999) then
    do XAYear = AStart, AEnd
      XCYear = CAStart + XAYear - AStart
      do XMonth = 1, NMonth
        if (DataA(XAYear,XMonth,MapCA(XCStn)).NE.-9999) then
        	TotA = TotA + 1
        	DataC   (XCYear,XMonth,XCStn) = DataA   (XAYear,XMonth,MapCA(XCStn))
        end if
      end do
    end do
  end if

  if (MapCB(XCStn).NE.-999) then
    do XBYear = BStart, BEnd
      XCYear = CBStart + XBYear - BStart
      do XMonth = 1, NMonth
        if (DataB(XBYear,XMonth,MapCB(XCStn)).NE.-9999) then
        	if (DataC(XCYear,XMonth,XCStn).EQ.-9999) then
        	  DataC   (XCYear,XMonth,XCStn) = DataB   (XBYear,XMonth,MapCB(XCStn))
        	  TotB = TotB + 1
        	else
        	  TotA = TotA - 1 ; TotAB = TotAB + 1
        	end if
        end if
      end do
    end do
  end if
end do

end subroutine CombineTheTwoAdd

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

subroutine CombineTheTwoIfData

TotA=0 ; TotB=0 ; TotAB=0

do XCStn = 1, NCStn
  if (MapCA(XCStn).NE.-999) then
   if (MapCB(XCStn).NE.-999) then
    do XAYear = AStart, AEnd
      XCYear = CAStart + XAYear -  AStart
      XBYear =  BStart + XCYear - CBStart
      
      do XMonth = 1, NMonth
        if (DataA(XAYear,XMonth,MapCA(XCStn)).NE.-9999) then
          if (DataB(XBYear,XMonth,MapCB(XCStn)).NE.-9999) then
        	TotAB = TotAB + 1
        	DataC (XCYear,XMonth,XCStn) = DataA (XAYear,XMonth,MapCA(XCStn))
          else
          	TotB = TotB + 1
          end if
        end if
      end do
    end do
   else
    TotA = TotA + 1
   end if
  end if
end do

end subroutine CombineTheTwoIfData

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

subroutine DescribeAPer

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select the .per file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".per")

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=Data,YearAD=YearAD,CallFile=LoadFileA)

NAStn = size(StnInfoA,1) ; NAYear = size(YearAD,1)

allocate (SumsPer(NAYear,17), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DescribeAPer: Allocation failure #####"
SumsPer = 0

do XAYear = 1, NAYear
  do XMonth = 1, NMonth
    do XAStn = 1, NAStn
      if (Data(XAYear,XMonth,XAStn).NE.-9999) SumsPer(XAYear,XMonth)=SumsPer(XAYear,XMonth)+1
    end do
  end do
end do

do XSeasonYear = 1, NAYear
 do XSeason = 1, 4
  do XSeasonMonth = (XSeason*3),(XSeason*3)+2
    XMonth=XSeasonMonth ; XAYear=XSeasonYear
    if (XMonth.GT.12) then
    	XMonth=XMonth-12 ; XAYear=XAYear+1
    end if
    if (XAYear.LE.NAYear) then
    	SumsPer(XSeasonYear,12+XSeason)=SumsPer(XSeasonYear,12+XSeason)+SumsPer(XAYear,XMonth)
    else
    	SumsPer(XSeasonYear,12+XSeason)=-999.0
    end if
  end do
  if (SumsPer(XSeasonYear,12+XSeason).NE.-999.0) &
  				SumsPer(XSeasonYear,12+XSeason)=SumsPer(XSeasonYear,12+XSeason)/3.0
 end do
end do

do XAYear = 1, NAYear
  do XMonth = 1, 12
    SumsPer(XAYear,17)=SumsPer(XAYear,17)+SumsPer(XAYear,XMonth)
  end do
  if (SumsPer(XAYear,17).NE.-9999) SumsPer(XAYear,17)=SumsPer(XAYear,17)/12.0
end do

LineFormat='(i5,12f9.1,5f9.1)'
call SavePER (SaveFileB,YearAD,2,AllData=SumsPer,NoResponse=1,CallLineFormat=LineFormat)

end subroutine DescribeAPer

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

subroutine DescribeAAnn

print*, "  > Select the .cts or .dtb file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".cts") ; if (SuffixBeg.GT.0) FileSuffix = ".cts"
        SuffixBeg = index(GivenFile,".dtb") ; if (SuffixBeg.GT.0) FileSuffix = ".dtb"
	if (ReadStatus.LE.0.AND.(FileSuffix.EQ.".cts".OR.FileSuffix.EQ.".dtb")) exit
end do
LoadFileA  = LoadPath (GivenFile,FileSuffix)

print*, "  > Select the .ann file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".ann")

print*, "  > Operating ..."
if (FileSuffix.EQ.'.cts') then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA)
else
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld, &
  			Lat=ALat,Lon=ALon,Elv=AElv,DtbNormals=DtbNormalsA, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA,silent=1)
end if

NAStn = size(StnInfoA,1) ; NAYear = size(YearAD,1)

call ClassifyContinent

print*, "  > Total valid data: ", TotA
call SaveANN (SaveFileB,YearAD,ColTitles,SumsPer,DecPlaceN=1)

end subroutine DescribeAAnn

!*******************************************************************************
! CAREFUL! this is called both from DescribeAAnn and Anomalise

subroutine ClassifyContinent

NCty=MasterSize()
allocate (CtyName  (NCty), &
	  CtyFinal (NCty), &
	  CtyCode0 (NCty), &
	  CtyCode1 (NCty), &
	  CtyReg   (NCty), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyContinent: Allocation failure #####"

call LoadMasterCty (CtyName,CtyFinal,CtyCode0,CtyCode1,CtyReg)

allocate (SumsPer(NAYear,12), &
	  ColTitles(12),      stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyContinent: Allocation failure #####"
SumsPer = 0 ; ColTitles = RegTitles

do XAStn = 1, NAStn
 if (AStn(XAStn).NE.MissVal) then
  XCty=0 ; XReg=0
  do
    XCty=XCty+1
    if (trim(StnCtyA(XAStn)).EQ.trim(CtyName(XCty))) XReg=CtyReg(XCty)
    if (XReg.NE.0) exit
    if (XCty.EQ.NCty) exit
  end do
  
  do XAYear = 1, NAYear
    do XMonth = 1, NMonth
      if (DataA(XAYear,XMonth,XAStn).NE.-9999) then
      	SumsPer(XAYear,1)=SumsPer(XAYear,1)+1
      	if (XReg.GE.1.AND.XReg.LE.11) SumsPer(XAYear,(XReg+1))=SumsPer(XAYear,(XReg+1))+1
      end if
    end do
  end do
 end if
end do

TotA = 0
do XAYear = 1, NAYear
  TotA = TotA + SumsPer(XAYear,1)
  do XReg = 1, 12
    SumsPer(XAYear,XReg)=SumsPer(XAYear,XReg)/12.0
  end do
end do

deallocate (CtyName,CtyFinal,CtyCode0,CtyCode1,CtyReg, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyContinent: Deallocation failure #####"

end subroutine ClassifyContinent 

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

subroutine OperateAk

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select the operation (A/k=1, A*k=2, A-k=3, A+k=4, sqrt=5, **k=6, abs=7): "
do
	read (*,*,iostat=ReadStatus), QOperation
	if (ReadStatus.LE.0 .AND. QOperation.GE.1 .AND. QOperation.LE.7) exit
end do

Constant = -999.0
if (QOperation.NE.5.AND.QOperation.NE.7) then
  print*, "  > Select the constant: "
  do
	read (*,*,iostat=ReadStatus), Constant
	if (ReadStatus.LE.0) exit
  end do
end if

call EnquireRestrict

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".cts")

print*, "  > Operating ..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=Data,YearAD=YearAD,CallFile=LoadFileA)

NAStn = size(StnInfoA,1) ; NAYear = size(YearAD,1)

do XAStn = 1, NAStn
  do XAYear = 1, NAYear
    do XMonth = 1, 12
     if (Data(XAYear,XMonth,XAStn).NE.-9999) then
      Data(XAYear,XMonth,XAStn) = &
      		OpAwithB (real(Data(XAYear,XMonth,XAStn)),Constant,QOperation)
      if (QRestrict.Eq.1.OR.QRestrict.Eq.3) then
        if (Data(XAYear,XMonth,XAStn).LT.Min) Data(XAYear,XMonth,XAStn) = MinImpose
      end if
      if (QRestrict.Eq.2.OR.QRestrict.Eq.3) then
        if (Data(XAYear,XMonth,XAStn).GT.Max) Data(XAYear,XMonth,XAStn) = MaxImpose
      end if
     end if
    end do
  end do
end do

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=Data,YearAD=YearAD,CallFile=SaveFileC)

end subroutine OperateAk

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

subroutine OperateAB

print*, "  > Select .cts file A to load (priority for headers):"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select the operation (A/B=1, A*B=2, A-B=3, A+B=4): "
do
	read (*,*,iostat=ReadStatus), QOperation
	if (ReadStatus.LE.0 .AND. QOperation.GE.1 .AND. QOperation.LE.4) exit
end do

print*, "  > Select .cts file B to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileB = LoadPath (GivenFile,".cts")

call EnquireRestrict

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".cts")

print*, "  > Manipulate parallel .src files ? (0=no,1=yes)"
do
	read (*,*,iostat=ReadStatus), QSource
	if (ReadStatus.LE.0.AND.QSource.GE.0.AND.QSource.LE.1) exit
end do
if (QSource.EQ.1) then
  LoadFileASrc = LoadFileA
  SuffixBeg = index(LoadFileASrc,".cts") ; LoadFileASrc(SuffixBeg:SuffixBeg+3) = ".src"
  SaveFileCSrc = SaveFileC
  SuffixBeg = index(SaveFileCSrc,".cts") ; SaveFileCSrc(SuffixBeg:SuffixBeg+3) = ".src"
end if

print*, "  > Loading ..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Code=BStn,OldCode=BStnOld,Lat=BLat,Lon=BLon,Elv=BElv, &
		Data=DataB,YearAD=BYearAD,CallFile=LoadFileB)
if (QSource.EQ.1) call LoadCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC, &
		Data=DataC,YearAD=CYearAD,CallFile=LoadFileASrc)

NAStn = size(StnInfoA,1) ; NAYear = size(AYearAD,1)
NBStn = size(StnInfoB,1) ; NBYear = size(BYearAD,1)
call CommonVecPer (AYearAD,BYearAD,AStart,AEnd,BStart,BEnd)
print "(a,2i5)", "   > Common period:", AYearAD(AStart),AYearAD(AEnd)

XBStn = 1 ; XAStn = 1
do
!  write (99,"(a,2i8)"), "ITERATION:", AStn(XAStn),BStn(XBStn)	! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
  if (AStn(XAStn).EQ.BStn(XBStn)) then
    do XAYear = 1, NAYear
      if (XAYear.LT.AStart.OR.XAYear.GT.AEnd) then
!        write (99,"(a,i5)"), "OUTRANGE:", AYearAD(XAYear)	! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
        do XMonth = 1, NMonth
          DataA(XAYear,XMonth,XAStn) = -9999 ; if (QSource.EQ.1) DataC(XAYear,XMonth,XAStn) = -9999
        end do 
      else
        XBYear = XAYear - AStart + BStart
!        write (99,"(a,2i5)"), "IN-RANGE:", AYearAD(XAYear),BYearAD(XBYear) ! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
        do XMonth = 1, NMonth
          if (DataA(XAYear,XMonth,XAStn).NE.-9999.AND.DataB(XBYear,XMonth,XBStn).NE.-9999) then
            OpVal = OpAwithB (real(DataA(XAYear,XMonth,XAStn)), &
            					   real(DataB(XBYear,XMonth,XBStn)),QOperation)
            DataA(XAYear,XMonth,XAStn) = nint(OpVal)
!            write (99,"(a,2i5)"), XMonth,DataA(XAYear,XMonth,XAStn),DataB(XBYear,XMonth,XBStn), &
!            				 ! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
      	    if (QRestrict.Eq.1.OR.QRestrict.Eq.3) then
              if (DataA(XAYear,XMonth,XAStn).LT.Min) DataA(XAYear,XMonth,XAStn) = MinImpose
            end if
            if (QRestrict.Eq.2.OR.QRestrict.Eq.3) then
              if (DataA(XAYear,XMonth,XAStn).GT.Max) DataA(XAYear,XMonth,XAStn) = MaxImpose
            end if
          else
            DataA(XAYear,XMonth,XAStn) = -9999 ; if (QSource.EQ.1) DataC(XAYear,XMonth,XAStn) = -9999
          end if
        end do 
      end if
    end do
    XAStn = XAStn + 1 ; XBStn = XBStn + 1
  else if (AStn(XAStn).LT.BStn(XBStn)) then  
!    write (99,"(a,i8)"), "REMOVING:", AStn(XAStn)	! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
    do XAYear = 1, NAYear
      do XMonth = 1, NMonth
        DataA(XAYear,XMonth,XAStn) = -9999 ; if (QSource.EQ.1) DataC(XAYear,XMonth,XAStn) = -9999
      end do
    end do
    XAStn = XAStn + 1
  else if (AStn(XAStn).GT.BStn(XBStn)) then
!    write (99,"(a,i8)"), "IGNORING:", BStn(XBStn)	! @@@@@@@@@@@@@@@@@@@@@@@@@@@@
    XBStn = XBStn + 1
  end if
  
  if (XAStn.GT.NAStn.OR.XBStn.GT.NBStn) exit
end do  
  
print*, "  > Saving ..."
call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=AYearAD,Data=DataA)	! save .cts file
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileC)
if (QSource.EQ.1) call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataC,YearAD=AYearAD,CallFile=SaveFileCSrc)

end subroutine OperateAB

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

subroutine EnquireRestrict

print*, "  > Impose min (=1), max (=2), both (=3), neither (=0):"
do
	read (*,*,iostat=ReadStatus), QRestrict
	if (ReadStatus.LE.0.AND.QRestrict.GE.0.AND.QRestrict.LE.3) exit
end do
if (QRestrict.EQ.1.OR.QRestrict.EQ.3) then
  print*, "  > Enter the minimum value to accept, and value to impose:"
  do
	read (*,*,iostat=ReadStatus), Min,MinImpose
	if (ReadStatus.LE.0) exit
  end do
end if
if (QRestrict.EQ.2.OR.QRestrict.EQ.3) then
  print*, "  > Enter the maximum value to accept, and value to impose:"
  do
	read (*,*,iostat=ReadStatus), Max,MaxImpose
	if (ReadStatus.LE.0) exit
  end do
end if

end subroutine EnquireRestrict

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

subroutine GetLatLongRect

print*, "  > Select the .cts or .dtb file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".cts") ; if (SuffixBeg.GT.0) FileSuffix = ".cts"
        SuffixBeg = index(GivenFile,".dtb") ; if (SuffixBeg.GT.0) FileSuffix = ".dtb"
	if (ReadStatus.LE.0.AND.(FileSuffix.EQ.".cts".OR.FileSuffix.EQ.".dtb")) exit
end do
LoadFileA  = LoadPath (GivenFile,FileSuffix)

print*, "  > Enter the min,max lat and min,max long:"
do
	read (*,*,iostat=ReadStatus), Lat0,Lat1,Lon0,Lon1

	if (ReadStatus.LE.0.AND.Lat0.GE.-90.AND.Lat1.LE.90.AND.Lat1.GE.Lat0) then
	  if (Lon0.GE.-180.AND.Lon0.LE.180.AND.Lon1.GE.-180.AND.Lon1.LE.180) then
	    if (Lon0.GT.Lon1) then
	      print*, "  > We note that the selected rectangle crosses the Date Line."
	      QCrossDL = 1
	    else
	      QCrossDL = 0
	    end if
	  else
	    print*, "  > Lat/long boundaries are unacceptable. Re-try."
	    ReadStatus = 1
	  end if
	else
	  print*, "  > Lat/long boundaries are unacceptable. Re-try."
	  ReadStatus = 1
	end if
	
	if (ReadStatus.LE.0) exit
end do

print*, "  > Enter the min,max elevation (-999,-999: no restriction): "
do
	read (*,*,iostat=ReadStatus), Elv0,Elv1

	if (ReadStatus.LE.0) then
	  if (Elv0.EQ.-999.AND.Elv1.EQ.-999) then
	    Elv0 = -100000 ; Elv1 = 100000
	  else
	    if (Elv0.GT.Elv1) then	    
	      print*, "  > Elevation boundaries are unacceptable. Re-try."
	      ReadStatus = 1
	    end if
	  end if
	end if
	
	if (ReadStatus.LE.0) exit
end do

if (FileSuffix.EQ.".cts") print*, "  > Select the .cts file to save:"
if (FileSuffix.EQ.".dtb") print*, "  > Select the .dtb file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,FileSuffix)

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

print*, "  > Loading..."
if (FileSuffix.EQ.".cts") then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld, &
  		Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=Data,YearAD=AYearAD,CallFile=LoadFileA)
else
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld, &
  		Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=Data,DtbNormals=DtbNormalsA,YearAD=AYearAD,CallFile=LoadFileA)
end if

NAStn = size(StnInfoA,1)
allocate (StnInfoB(NAStn,8), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetLatLongRect: alloc StnInfoB #####"

where (ALat.LT.Lat0.OR.ALat.GT.Lat1) AStn = -999
if (Elv0.NE.-999.OR.Elv1.EQ.-999) then
  where (AElv.LT.Elv0.OR.AElv.GT.Elv1) AStn = -999
end if

if (QCrossDL.EQ.0) then
  where (ALon.LT.(Lon0).OR.ALon.GT.(Lon1)) AStn = -999
else
  where (ALon.GT.(Lon1).AND.ALon.LT.(Lon0)) AStn = -999
end if

print*, "  > Saving..."

if (FileSuffix.EQ.".cts") then
  call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,1, &
  			YearAD=AYearAD,Data=Data)
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=Data, &
  		YearAD=AYearAD,CallFile=SaveFileA)
else
  call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,1, &
  			YearAD=AYearAD,Data=Data,DtbNormals=DtbNormalsA)
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=Data, &
  		YearAD=AYearAD,CallFile=SaveFileA,DtbNormals=DtbNormalsA)
end if

end subroutine GetLatLongRect

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

subroutine TruncateDate

print*, "  > Select the .cts or .dtb file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".cts") ; if (SuffixBeg.GT.0) FileSuffix = ".cts"
        SuffixBeg = index(GivenFile,".dtb") ; if (SuffixBeg.GT.0) FileSuffix = ".dtb"
	if (ReadStatus.LE.0.AND.(FileSuffix.EQ.".cts".OR.FileSuffix.EQ.".dtb")) exit
end do
LoadFileA  = LoadPath (GivenFile,FileSuffix)

print*, "  > Loading..."

if (FileSuffix.EQ.".cts") then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
else
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,DtbNormals=DtbNormalsA,YearAD=AYearAD,CallFile=LoadFileA)
end if

NAYear = size(AYearAD,1)
NAStn  = size(DataA,  3)
print "(a,i4,a1,i4)", "   > The loaded file extends: ", AYearAD(1), "-", AYearAD(NAYear)

print*, "  > Select the initial year, month in the truncated file:"
do
	read (*,*,iostat=ReadStatus), YearAD0,Month0
	if (ReadStatus.LE.0.AND.Month0.GE.1.AND.Month0.LE.12) exit
end do

print*, "  > Select the final year, month in the truncated file:"
do
	read (*,*,iostat=ReadStatus), YearAD1,Month1
	if (ReadStatus.LE.0.AND.Month1.GE.1.AND.Month1.LE.12.AND.YearAD1.GE.YearAD0) exit
end do

if (FileSuffix.EQ.".cts") print*, "  > Select the .cts file to save:"
if (FileSuffix.EQ.".dtb") print*, "  > Select the .dtb file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,FileSuffix)
print*, "  > Truncating..."

NBYear = YearAD1 - YearAD0 + 1

allocate (BYearAD(NBYear), &
	  DataB  (NBYear,NMonth,NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: TruncateDate: Allocation failure: BYearAD,DataB #####"
DataB = -9999

do XBYear = 1, NBYear
  BYearAD(XBYear) = YearAD0 + XBYear - 1
end do

call CommonVecPer (AYearAD,BYearAD,AStart,AEnd,BStart,BEnd)

do XBYear = BStart, BEnd
  XAYear = AStart + XBYear - BStart
  
  do XMonth = 1, NMonth
   if      (XBYear.EQ.BStart.AND.XMonth.LT.Month0) then
   			! leave DataB blank
   else if (XBYear.EQ.BEnd  .AND.XMonth.GT.Month1) then
   			! leave DataB blank
   else
     do XAStn = 1, NAStn
      DataB(XBYear,XMonth,XAStn) = DataA(XAYear,XMonth,XAStn)
     end do
   end if
  end do
end do

print*, "  > Saving..."

call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=BYearAD,Data=DataB)

if (FileSuffix.EQ.".cts") then
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataB,YearAD=BYearAD,CallFile=SaveFileB)
else
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataB,DtbNormals=DtbNormalsA, &
  		YearAD=BYearAD,CallFile=SaveFileB)
end if

end subroutine TruncateDate

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

subroutine QCImplicitStdev

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Process parallel source file? (1=no,2=yes) ?"
do
	read (*,*,iostat=ReadStatus), QSource
	if (ReadStatus.LE.0.AND.QSource.GE.1.AND.QSource.LE.2) exit
end do

print*, "  > Enter the min. no. of values required to calc a stdev:"
do
	read (*,*,iostat=ReadStatus), ThreshStdev
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do

print*, "  > If no stdev, keep (=1) or reject (=2) data ?"
do
	read (*,*,iostat=ReadStatus), QNoStdev
	if (ReadStatus.LE.0.AND.QNoStdev.GE.1.AND.QNoStdev.LE.2) exit
end do

print*, "  > Enter the no. of stdevs at which to flag data:"
do
	read (*,*,iostat=ReadStatus), ThreshReject
	if (ReadStatus.LE.0.AND.ThreshReject.GT.0) exit
end do

print*, "  > Reject (=1) or manually check (=2) flagged data ?"
do
	read (*,*,iostat=ReadStatus), QManCheck
	if (ReadStatus.LE.0.AND.QManCheck.GE.1.AND.QManCheck.LE.2) exit
end do

if (QManCheck.EQ.2) then
  print*, "  > Set earliest, latest years to check manually (-999=no limit):"
  do
	read (*,*,iostat=ReadStatus), CheckYear0, CheckYear1
!	if (CheckYear0.NE.-999.AND.(CheckYear0.LT.YearAD(1).OR. &
!			CheckYear0.GT.YearAD(NAYear))) ReadStatus=1
!	if (CheckYear1.NE.-999.AND.(CheckYear1.LT.YearAD(1).OR. &
!			CheckYear1.GT.YearAD(NAYear))) ReadStatus=1
	if (CheckYear0.NE.-999.AND.CheckYear1.NE.-999.AND.CheckYear0.GT.CheckYear1) &
			ReadStatus=1
	if (ReadStatus.LE.0) exit
  end do
end if

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileB = SavePath (GivenFile,".cts")

if (QSource.EQ.2) then
  LoadFileASrc=LoadfileA
  SuffixBeg=index(LoadFileASrc,".cts")
  LoadFileASrc(SuffixBeg:SuffixBeg+3)=".src"
  SaveFileBSrc=SaveFileB
  SuffixBeg=index(SaveFileBSrc,".cts")
  SaveFileBSrc(SuffixBeg:SuffixBeg+3)=".src"
end if

print*, "  > Loading, processing ..."

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Oldcode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=YearAD,CallFile=LoadFileA)
if (QSource.EQ.2) then
  call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB, &
			Data=DataASrc,YearAD=BYearAD,CallFile=LoadFileASrc)
  deallocate (StnInfoB,StnLocalB,StnNameB,StnCtyB,BYearAD, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: QCImplicitStdev: Allocation failure: St* #####"
end if

NAYear = size(DataA,1) ; NAStn = size(DataA,3)

allocate (StnAnn  (NAYear), &
	  LogStdev(1,NMonth,NAStn), &
	  LogMeans(1,NMonth,NAStn), &
	  Stdev   (1,NMonth,NAStn), &
	  Means   (1,NMonth,NAStn), &
	  CYearAD (1), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QCImplicitStdev: Allocation failure: St* #####"

ThreshOpCalc = 100.0 * (real(NAYear) - ThreshStdev) / real(NAYear)
    
do XMonth = 1, NMonth
  do XAStn = 1, NAStn
    do XAYear = 1, NAYear
      StnAnn(XAYear) = DataA(XAYear,XMonth,XAStn)
    end do
    
    Means(1,XMonth,XAStn) = OpCalcMean  (StnAnn,NAYear,ThreshOpCalc,CallMissVal=-9999.0)
    Stdev(1,XMonth,XAStn) = OpCalcStDev (StnAnn,NAYear,ThreshOpCalc,2,CallMissVal=-9999.0)
  end do
end do

allocate (LogPass(NAYear+1,NMonth+1), &
	  LogMiss(NAYear+1,NMonth+1), &
	  LogSkip(NAYear+1,NMonth+1), &
	  LogFail(NAYear+1,NMonth+1), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QCImplicitStdev: Allocation failure: DataB* #####"
LogPass=0 ; LogFail=0 ; LogSkip=0 ; LogMiss=0

do XAStn = 1, NAStn
  ThreshLo=-9999 ; ThreshHi=-9999 ; SeqStatus=0.0
  
  do XMonth = 1, NMonth
    if (Stdev(1,XMonth,XAStn).NE.-9999.AND.Means(1,XMonth,XAStn).NE.-9999) then
      ThreshLo(XMonth)=Means(1,XMonth,XAStn)-(Stdev(1,XMonth,XAStn)*ThreshReject)
      ThreshHi(XMonth)=Means(1,XMonth,XAStn)+(Stdev(1,XMonth,XAStn)*ThreshReject)
    end if
  end do
  
  do XAYear = 1, NAYear
    InRange=1
    if (CheckYear0.NE.-999.AND.YearAD(XAYear).LT.CheckYear0) InRange=0
    if (CheckYear1.NE.-999.AND.YearAD(XAYear).GT.CheckYear1) InRange=0
    
    do XMonth = 1, NMonth
      if (DataA(XAYear,XMonth,XAStn).NE.-9999) then
        if (ThreshLo(XMonth).NE.-9999.AND.Stdev(1,XMonth,XAStn).GT.0.0) then
	  Excess = (DataA(XAYear,XMonth,XAStn)-Means(1,XMonth,XAStn))/Stdev(1,XMonth,XAStn)
          SeqStatus=SeqStatus+abs(Excess)-ThreshReject+1
          
       	  if (DataA(XAYear,XMonth,XAStn).GE.ThreshLo(XMonth) &
       				.AND.DataA(XAYear,XMonth,XAStn).LE.ThreshHi(XMonth)) then
            LogPass(XAYear,XMonth) = LogPass(XAYear,XMonth) + 1
          else
            write(99,"(i7,x,a20,i3,i5,i6,2f6.2,2f8.2)"),AStn(XAStn),StnNameA(XAStn),XMonth,&
            		YearAD(XAYear),DataA(XAYear,XMonth,XAStn),Excess,SeqStatus, &
            		ThreshLo(XMonth),ThreshHi(XMonth)

            if (QManCheck.EQ.1.OR.InRange.EQ.0.OR.abs(Excess).GT.real(ThreshReject+1)) then
              DataA (XAYear,XMonth,XAStn) = -9999
              if (QSource.EQ.2) DataASrc(XAYear,XMonth,XAStn) = -9999
              LogFail(XAYear,XMonth) = LogFail(XAYear,XMonth) + 1
            else if (SeqStatus.GT.4.0) then
              call ManualCheck
            else
              LogPass(XAYear,XMonth) = LogPass(XAYear,XMonth) + 1              
            end if
          end if
          
          if (SeqStatus.GT.0) SeqStatus=SeqStatus-0.5
          if (SeqStatus.LT.0) SeqStatus=0.0
        else						! no stdev calculable
          if (QNoStdev.EQ.1) then
            LogSkip(XAYear,XMonth) = LogSkip(XAYear,XMonth) + 1
          else
            DataA   (XAYear,XMonth,XAStn) = -9999
            if (QSource.EQ.2) DataASrc(XAYear,XMonth,XAStn) = -9999
            LogSkip(XAYear,XMonth) = LogSkip(XAYear,XMonth) + 1
          end if
        end if
      else						! data is missing
        LogMiss(XAYear,XMonth) = LogMiss(XAYear,XMonth) + 1
      end if
    end do
  end do
end do

print*, "  > Saving QC'd data ..."
call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=YearAD,Data=DataA)
call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD, &
		CallFile=SaveFileB,Silent=1)
if (QSource.EQ.2) call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,Data=DataASrc, &
		YearAD=YearAD,CallFile=SaveFileBSrc,Silent=1)

print*, "  > Saving means and stdev to extra log files ..."
CYearAD = YearAD(1)

Means=Means*10.0 ; LogMeans=Means
Stdev=Stdev*10.0 ; LogStdev=Stdev

call system ('rm /cru/mikeh1/f709762/scratch/log/log-opcruts-means.cts')
call system ('rm /cru/mikeh1/f709762/scratch/log/log-opcruts-stdev.cts')

call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=CYearAD,Data=LogMeans)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=LogMeans,YearAD=CYearAD,CallFile=LogMeansFile,Silent=1)

deallocate (StnInfoC, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QCImplicitStdev: Deallocation failure: StnInfoC* #####"

call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=CYearAD,Data=LogStdev)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=LogStdev,YearAD=CYearAD,CallFile=LogStdevFile,Silent=1)

print*, "  > Saving QC statistics to main log file ..."

do XLog = 1, 4
 if (XLog.EQ.1) LogAny => LogPass
 if (XLog.EQ.2) LogAny => LogMiss
 if (XLog.EQ.3) LogAny => LogSkip
 if (XLog.EQ.4) LogAny => LogFail
 
 do XAYear = 1, NAYear
   do XMonth = 1, NMonth
      LogAny(XAYear,NMonth+1) = LogAny(XAYear,NMonth+1) + LogAny(XAYear,XMonth)
   end do
 end do
  
 do XMonth = 1, NMonth
   do XAYear = 1, NAYear
      LogAny(NAYear+1,XMonth) = LogAny(NAYear+1,XMonth) + LogAny(XAYear,XMonth)
   end do
 end do  
 
 do XMonth = 1, NMonth
      LogAny(NAYear+1,NMonth+1) = LogAny(NAYear+1,NMonth+1) + LogAny(NAYear+1,XMonth)
 end do
 
 write (99,"(a4,a)"),LogName(XLog), &
 	"   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec    TOT"
 do XAYear = 1, NAYear
  write (99,"(i4,12i6,i7)"), YearAD(XAYear),(LogAny(XAYear,XMonth),XMonth=1,NMonth+1)
 end do
 write (99,"(a4,12i6,i7)"), "TOT ",(LogAny(NAYear+1,XMonth),XMonth=1,NMonth+1)
 write (99,"(a)"), " "
 
 nullify(LogAny)
end do

end subroutine QCImplicitStdev

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

subroutine ManualCheck

print "(a,i8,i6,x,a20,x,a13,i3,i5,f8.2)", "   > ", AStn(XAStn),XAStn,StnNameA(XAStn), &
			StnCtyA(XAStn),XMonth,YearAD(XAYear),Excess
!print "(a4,12f7.2)", "LO  ", (ThreshLo(XBMonth),XBMonth=1,12)
!print "(a4,12f7.2)", "HI  ", (ThreshHi(XBMonth),XBMonth=1,12)
print*, "  > Reject single (=2), sequence (=1), or none (=0) ?"
do
	    read (*,*,iostat=ReadStatus), QReject
	    if (ReadStatus.LE.0.AND.QReject.GE.0.AND.QReject.LE.2) exit
end do
if      (QReject.EQ.0) then
              LogPass(XAYear,XMonth) = LogPass(XAYear,XMonth) + 1
else if (QReject.EQ.2) then
              DataA (XAYear,XMonth,XAStn) = -9999
              if (QSource.EQ.2) DataASrc(XAYear,XMonth,XAStn) = -9999
              LogFail(XAYear,XMonth) = LogFail(XAYear,XMonth) + 1
else if (QReject.EQ.1) then
              print*, "  > Enter the first mon+yr, last mon+yr to reject: "
	      do
	        read (*,*,iostat=ReadStatus), XMonth0,XYear0,XMonth1,XYear1
	        if (ReadStatus.LE.0.AND.XYear1.GE.XYear0.AND.XMonth0.GE.1.AND.XMonth0.LE.12 &
	        		.AND.XMonth1.GE.1.AND.XMonth1.LE.12 &
	        		.AND.XYear0.GE.YearAD(1).AND.XYear1.LE.YearAD(NAYear)) exit
	      end do
	      do XBYear=(XYear0-YearAD(1)+1),(XYear1-YearAD(1)+1)
	        do XBMonth=1,NMonth
	          if (XBYear.EQ.(XYear0-YearAD(1)+1).AND.XBMonth.LT.XMonth0) then
	          else if (XBYear.EQ.(XYear1-YearAD(1)+1).AND.XBMonth.GT.XMonth1) then
	          else
	            if (DataA(XBYear,XBMonth,XAStn).NE.-9999) then
	              DataA(XBYear,XBMonth,XAStn) = -9999
	              if (QSource.EQ.2) DataASrc(XBYear,XBMonth,XAStn) = -9999
            	      LogFail(XAYear,XMonth) = LogFail(XAYear,XMonth) + 1
	              
	              if (XBYear.LT.XAYear.OR.(XBYear.EQ.XAYear.AND.XBMonth.LT.XMonth)) then
            	        LogPass(XAYear,XMonth) = LogPass(XAYear,XMonth) - 1
	              else if (XBYear.EQ.XAYear.AND.XBMonth.EQ.XMonth) then
	              else
            	        LogMiss(XAYear,XMonth) = LogMiss(XAYear,XMonth) - 1
	              end if            
	            end if
	          end if
	        end do
	      end do
end if
! LastCheck = XAStn

end subroutine ManualCheck

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

subroutine Anomalise

print*, "  > This software is now a separate program: anomdtb.f90"

end subroutine Anomalise

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

subroutine FindClosest

print*, "  > Select the .cts or .dtb file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".cts") ; if (SuffixBeg.GT.0) FileSuffix = ".cts"
        SuffixBeg = index(GivenFile,".dtb") ; if (SuffixBeg.GT.0) FileSuffix = ".dtb"
	if (ReadStatus.LE.0.AND.(FileSuffix.EQ.".cts".OR.FileSuffix.EQ.".dtb")) exit
end do
LoadFileA  = LoadPath (GivenFile,FileSuffix)

print*, "  > Enter the lat and long: "
do
	read (*,*,iostat=ReadStatus), Lat0,Lon0
	if (ReadStatus.LE.0.AND.Lat0.GE.-90.AND.Lat0.LE.90.AND.Lon0.GE.-180.AND.Lon0.LE.180) exit
end do

print*, "  > Find the closest N stations. Specify N."
do
	read (*,*,iostat=ReadStatus), NBStn
	if (ReadStatus.LE.0.AND.NBStn.GE.1) exit
end do

print*, "  > Restrict to stns with data in period: beg,end (-999=don't restrict):"
do
	read (*,*,iostat=ReadStatus), YearAD0,YearAD1
	if (ReadStatus.LE.0) exit
end do

print*, "  > Loading, processing ..."

if (FileSuffix.EQ.".cts") then
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,oldcode=AStnOld, &
  		Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
else
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,oldcode=AStnOld, &
  		Lat=ALat,Lon=ALon,Elv=AElv,DtbNormals=DtbNormalsA, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
end if
NAYear = size(DataA,1) ; NAStn = size(DataA,3)

allocate (Scatter(NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FindClosest: Allocation failure: Scatter #####"
Scatter = MissVal

do XAStn = 1, NAStn
  if ( (YearAD0.EQ.MissVal.AND.YearAD1.EQ.MissVal).OR. &
  		(StnInfoA(XAStn,5).LE.YearAD0.AND.StnInfoA(XAStn,6).GE.YearAD0) .OR. &
  		(StnInfoA(XAStn,5).LE.YearAD1.AND.StnInfoA(XAStn,6).GE.YearAD1) .OR. & 
  		(StnInfoA(XAStn,5).GE.YearAD0.AND.StnInfoA(XAStn,6).LE.YearAD1) ) then
    if (ALat(XAStn).NE.MissVal.AND.ALon(XAStn).NE.MissVal) &
    			Scatter(XAStn) = GetDistance (ALat(XAStn),ALon(XAStn),Lat0,Lon0)
  end if
end do

do XBStn = 1, NBStn
  Limit = 100000.0
  do XAStn = 1, NAStn
    if (Scatter(XAStn).NE.MissVal.AND.Scatter(XAStn).LT.Limit) then
      Limit = Scatter(XAStn) ; XCStn = XAStn
    end if
  end do
  print "(i7,2f8.2,i6,x,a20,x,a13,2(x,i4),f8.2)", AStn(XCStn),ALat(XCStn),ALon(XCStn),nint(AElv(XCStn)), &
  		StnNameA(XCStn),StnCtyA(XCStn),StnInfoA(XCStn,5),StnInfoA(XCStn,6),Scatter(XCStn)
  Scatter(XCStn) = MissVal
end do

end subroutine FindClosest

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

subroutine MultiplyCodes

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".cts")

print*, "  > Enter the factor by which to multiply codes: "
do
	read (*,*,iostat=ReadStatus), Factor
	if (ReadStatus.LE.0) exit
end do

print*, "  > Multiply positive (>0) codes only ? (0=no,1=yes) "
do
	read (*,*,iostat=ReadStatus), QPosOnly
	if (ReadStatus.LE.0.AND.QPosOnly.GE.0.AND.QPosOnly.LE.1) exit
end do

print*, "  > Store the original code as 'old code' ? (0=no,1=yes)"
do
	read (*,*,iostat=ReadStatus), QStoreOld
	if (ReadStatus.LE.0.AND.QStoreOld.GE.0.AND.QStoreOld.LE.1) exit
end do

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,oldcode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
			Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
NAStn = size(DataA,3)

if (QStoreOld.EQ.1) then
 do XAStn = 1, NAStn
  AStnOld(XAStn) = AStn(XAStn)
 end do
end if

do XAStn = 1, NAStn
  if (QPosOnly.EQ.0.OR.AStn(XAStn).GT.0) AStn(XAStn) = nint(real(AStn(XAStn)) * Factor)
end do

call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,0,YearAD=AYearAD,Data=DataA)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileA)

end subroutine MultiplyCodes

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

subroutine SyncCtsSrc

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")
SubBeg=index(LoadFileA,".cts")
LoadFileB=LoadFileA ; LoadFileB(SubBeg:(SubBeg+3))=".src"

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA, &
			Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB, &
			Data=DataB,YearAD=BYearAD,CallFile=LoadFileB)
NAYear = size(DataA,1) ; NAStn = size(DataA,3) ; Log1=0 ; Log2=0 ; XInfo=MultiSource

do XAStn = 1, NAStn
  do XAYear = 1, NAYear
    do XMonth = 1, NMonth
      if      (DataA(XAYear,XMonth,XAStn).EQ.DataMissVal.AND. &
      	       DataB(XAYear,XMonth,XAStn).NE.DataMissVal) then
        Log1=Log1+1
	write (99,"(6i8)"), 1,AYearAD(XAYear),XMonth,StnInfoA(XAStn,1), &
				DataA(XAYear,XMonth,XAStn),DataB(XAYear,XMonth,XAStn)
	DataB(XAYear,XMonth,XAStn) = DataMissVal
      else if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal.AND. &
      	       DataB(XAYear,XMonth,XAStn).EQ.DataMissVal) then
        Log2=Log2+1
	write (99,"(6i8)"), 2,AYearAD(XAYear),XMonth,StnInfoA(XAStn,1), &
				DataA(XAYear,XMonth,XAStn),DataB(XAYear,XMonth,XAStn)
	DataB(XAYear,XMonth,XAStn) = XInfo
      else if (DataB(XAYear,XMonth,XAStn).NE.DataMissVal) then
        XInfo=DataB(XAYear,XMonth,XAStn)
      end if
    end do
  end do
end do

if (Log1.EQ.0.AND.Log2.EQ.0) then
  print*, "  > The .cts and .src files are already in sync."
else
  if (Log1.GT.0) print "(a,i10)", "   > X. Miss data but valid src:", Log1
  if (Log2.GT.0) print "(a,i10)", "   > Y. Valid data but miss src:", Log2
end if

if (Log1.GT.0) then
  print*, "  > Select the corrected .cts file to save:"
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  SaveFileA = SavePath (GivenFile,".cts")
  SubBeg=index(SaveFileA,".cts")
  SaveFileB=SaveFileA ; SaveFileB(SubBeg:(SubBeg+3))=".src"

  call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileA)
  call SaveCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Data=DataB,YearAD=BYearAD,CallFile=SaveFileB)
end if

end subroutine SyncCtsSrc

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

subroutine ImpossibleElev

print*, "  > Select the file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
StrLen=len_trim(GivenFile)
FileSuffix=GivenFile((StrLen-3):StrLen)
LoadFileA = LoadPath (GivenFile,FileSuffix)
StrLen=len_trim(LoadFileA)

print*, "  > Select the file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
FileSuffix=LoadFileA((StrLen-3):StrLen)
SaveFileA = SavePath (GivenFile,FileSuffix)

if      (LoadFileA((StrLen-3):StrLen).EQ.'.hdr') then
  if (index(LoadFileA,"master").NE.0) then
    print*, "  > Loading master .hdr file..."
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=LoadFileA,HeadOnly=1,Elv=AElv)
  else
    print*, "  > Loading accessions .hdr file..."
    call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=LoadFileA,HeadOnly=1,Elv=AElv, &
			SrcCode=StnSrcCodeA,SrcSuffix=StnSrcSuffixA,SrcDate=StnSrcDateA)
  end if
else
  print*, "  > Loading database file..."
  call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD,Elv=AElv, &
  			CallFile=LoadFileA,DtbNormals=DtbNormalsA)
end if 

NAStn=size(StnInfoA,1)
do XAStn = 1, NAStn
  if (AElv(XAStn).NE.-999.AND.(AElv(XAStn).LT.-500.OR.AElv(XAStn).GT.6000)) then
    AElv(XAStn)=MissVal ; StnInfoA(XAStn,4)=MissVal ; OpEn=OpEn+1
  end if
end do
print "(a,i5)", "   > Stns corrected total: ", nint(OpEn)

if      (LoadFileA((StrLen-3):StrLen).EQ.'.hdr') then
  if (index(LoadFileA,"master").NE.0) then
    print*, "  > Saving master .hdr file..."
    call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA,HeadOnly=1)
  else
    print*, "  > Saving accessions .hdr file..."
    call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA,HeadOnly=1, &
			SrcCode=StnSrcCodeA,SrcSuffix=StnSrcSuffixA,SrcDate=StnSrcDateA)
  end if
else
  print*, "  > Saving database file..."
  call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=YearAD, &
  			CallFile=SaveFileA,DtbNormals=DtbNormalsA)
end if 

end subroutine ImpossibleElev

!*******************************************************************************
! convert degrees and minutes in .cts header lines to degrees and hundedths

subroutine ConvertDegHun

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".cts")

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA, &
		Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
NAStn = size(DataA,3)

do XAStn = 1, NAStn
  if (ALat(XAStn).NE.-999.0) then
     OpDiff = mod (ALat(XAStn),1.0)
     OpVal  = ALat(XAStn) - OpDiff
     OpDiff = OpDiff / 0.60
     ALat(XAStn) = OpVal + OpDiff
  end if
  if (ALon(XAStn).NE.-999.0) then
     OpDiff = mod (ALon(XAStn),1.0)
     OpVal  = ALon(XAStn) - OpDiff
     OpDiff = OpDiff / 0.60
     ALon(XAStn) = OpVal + OpDiff
  end if
end do

call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,0,YearAD=AYearAD,Data=DataA)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileA)

end subroutine ConvertDegHun

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

subroutine MultiplyNRM

print*, "  > We also looks out for dodgy data missing values (-999)."
print*, "  > Select the .nrm file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the .nrm file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".nrm")

print*, "  > Select the constant: "
do
	read (*,*,iostat=ReadStatus), Factor
	if (ReadStatus.LE.0) exit
end do

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=LoadFileA,Silent=1, &     
		NmlSrc=NmlSrc,NmlInc=NmlInc,NmlData=Normals)
NAStn = size(StnNameA,1)
do XAStn = 1, NAStn
  do XMonth = 1, NMonth
    if (Normals(XMonth,XAStn).EQ.-999) then
        Normals(XMonth,XAStn) = DataMissVal
    else if (Normals(XMonth,XAStn).NE.DataMissVal) then
    	Normals(XMonth,XAStn) = nint(real(Normals(XMonth,XAStn))*Factor)
    end if
  end do
end do
call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA,Silent=1, &     
		NmlSrc=NmlSrc,NmlInc=NmlInc,NmlData=Normals)

end subroutine MultiplyNRM

!*******************************************************************************
! convert different variable databases (.dtb to .cts)
! option 1 = .rhm to .vap

subroutine ConvertVariableDtb

print*, "  > Select the conversion (1=rhm->vap):"
do
	read (*,*,iostat=ReadStatus), QMerge
	if (ReadStatus.LE.0.AND.QMerge.GE.1.AND.QMerge.LE.1) exit
end do

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".cts")

if (QMerge.EQ.1) DatabaseFilter="/tyn1/tim/data/cruts/database/rhm.??????????.dtb"
call GetBatch (DatabaseFilter,DatabaseBatch,Silent=1)
NBatch = size(DatabaseBatch,1)
DatabaseFileL = DatabaseBatch(NBatch)			! i.e. most recently dated file
Date=DatabaseFileL(36:43) ; TimeVal=DatabaseFileL(44:47)
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD, &
		CallFile=DatabaseFileL,Silent=1)
NAStn = size(StnNameA,1) ; NAYear = size(AYearAD,1)
print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a)", "   > Database file at   ",Date(7:8),".",Date(5:6),&
		".",Date(3:4)," at ",TimeVal(1:2),":",TimeVal(3:4)," has ",(NAstn)," stns."
deallocate (StnInfoA,StnLocalA,StnNameA,StnCtyA,AYearAD,DatabaseBatch,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableDtb: Deallocation failure: 1 #####"

DatabaseSrcFileL = DatabaseFileL
SubBeg = index(DatabaseSrcFileL,".dtb")
DatabaseSrcFileL(SubBeg:SubBeg+3) = ".dts"
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataC,YearAD=AYearAD,CallFile=DatabaseSrcFileL,Silent=1)
deallocate (StnInfoA,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableDtb: Deallocation failure: 2 #####"

if (QMerge.EQ.1) DatabaseFilter="/tyn1/tim/data/cruts/database/tmp.??????????.dtb"
call GetBatch (DatabaseFilter,DatabaseBatch,Silent=1)
NBatch = size(DatabaseBatch,1)
DatabaseFileL = DatabaseBatch(NBatch)			! i.e. most recently dated file
Date=DatabaseFileL(36:43) ; TimeVal=DatabaseFileL(44:47)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Code=BStn,Lat=BLat,Lon=BLon,Elv=BElv, &
		Data=DataB,YearAD=BYearAD,CallFile=DatabaseFileL,Silent=1)
NBStn = size(StnNameB,1) ; NBYear = size(BYearAD,1)
print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a)", "   > Database file at   ",Date(7:8),".",Date(5:6),&
		".",Date(3:4)," at ",TimeVal(1:2),":",TimeVal(3:4)," has ",(NBStn)," stns."
deallocate (StnInfoB,DatabaseBatch,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableDtb: Deallocation failure: 3 #####"

call CommonVecPer (AYearAD,BYearAD,AStart,AEnd,BStart,BEnd)
print "(a,2i5)", "   > Common period:", AYearAD(AStart),AYearAD(AEnd)

XAStn=1 ; XBStn=1
do
  if      (AStn(XAStn).EQ.BStn(XBStn)) then
    do XAYear = 1, NAYear
      XBYear = XAYear-AStart+BStart
      if (XBYear.LT.0.OR.XBYear.GT.NBYear) then
        do XMonth = 1, NMonth
          DataA(XAYear,XMonth,XAStn) = DataMissVal
          DataC(XAYear,XMonth,XAStn) = DataMissVal
        end do
      else
        do XMonth = 1, NMonth
          if (DataA(XAYear,XMonth,XAStn).GE.0.AND.DataA(XAYear,XMonth,XAStn).LE.1000.AND. &
          				DataB(XBYear,XMonth,XBStn).NE.DataMissVal) then
            if (QMerge.EQ.1) then			! rhm + tmp --> vap
              if (DataB(XBYear,XMonth,XBStn).GE.2) then		! assume Tw>0
                OpVal = 17.38 * (real(DataB(XBYear,XMonth,XBStn))/10.0)
                OpVal = OpVal / (239.0 + (real(DataB(XBYear,XMonth,XBStn))/10.0))
              else						! assume Tw<0
                OpVal = 22.44 * (real(DataB(XBYear,XMonth,XBStn))/10.0)
                OpVal = OpVal / (272.4 + (real(DataB(XBYear,XMonth,XBStn))/10.0))
              end if
              OpVal = 6.107 * exp(OpVal)
              OpVal = OpVal * real(DataA(XAYear,XMonth,XAStn)) / 1000.0
              DataA(XAYear,XMonth,XAStn) = nint(OpVal*10.0)	! we also retain DataC
            end if
          else
            DataA(XAYear,XMonth,XAStn) = DataMissVal
            DataC(XAYear,XMonth,XAStn) = DataMissVal
          end if
        end do
      end if
    end do
    XAStn=XAStn+1 ; XBStn=XBStn+1
  else if (AStn(XAStn).LT.BStn(XBStn)) then
    AStn(XAStn)=MissVal
    XAStn=XAStn+1
  else if (AStn(XAStn).GT.BStn(XBStn)) then
    XBStn=XBStn+1
  end if

  if (XAStn.GT.NAStn.OR.XBStn.GT.NBStn) exit
end do

print*, "  > Saving to .cts and .src ..."
call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=AYearAD,Data=DataA)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileC)

SuffixBeg=index(SaveFileC,".cts")
SaveFileC=SaveFileC(1:SuffixBeg-1) // ".src"
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataC,YearAD=AYearAD,CallFile=SaveFileC)

end subroutine ConvertVariableDtb

!*******************************************************************************
! convert different variable .cts file (.cts to .cts)
! option 1 = .snh to .spc

subroutine ConvertVariableCts

print*, "  > Select the conversion (1=snh->spc):"
do
	read (*,*,iostat=ReadStatus), QMerge
	if (ReadStatus.LE.0.AND.QMerge.GE.1.AND.QMerge.LE.1) exit
end do

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")

print*, "  > Process parallel source file? (1=no,2=yes) ?"
do
	read (*,*,iostat=ReadStatus), QSource
	if (ReadStatus.LE.0.AND.QSource.GE.1.AND.QSource.LE.2) exit
end do

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".cts")

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA, &
		Code=AStn,OldCode=AStnOld,Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA)
NAYear=size(DataA,1) ; NAStn=size(DataA,3)

if (QSource.EQ.2) then
  LoadFileASrc=LoadfileA
  SuffixBeg=index(LoadFileASrc,".cts")
  LoadFileASrc(SuffixBeg:SuffixBeg+3)=".src"
  SaveFileASrc=SaveFileA
  SuffixBeg=index(SaveFileA,".cts")
  SaveFileASrc(SuffixBeg:SuffixBeg+3)=".src"

  call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB, &
			Data=DataASrc,YearAD=BYearAD,CallFile=LoadFileASrc)
			
  deallocate (StnInfoB,StnLocalB,StnNameB,StnCtyB,BYearAD, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableCts: dealloc: B* #####"
end if

if (QMerge.EQ.1) then
  allocate (MonthLengths(NAYear,NMonth), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableCts: alloc: ML #####"

  do XAStn=1,NAStn
    if (ALat(XAStn).NE.MissVal) then
      LatRad=ALat(XAStn)*3.1415927/180.0
      do XMonth=1,NMonth		! calc max snh per day
        Numer=sin(Ephemeris(XMonth))*sin(LatRad)*(0.0-1.0)
	Denom=cos(Ephemeris(XMonth))*cos(LatRad)
	
	if ((Numer/Denom).GE.-1.AND.(Numer/Denom).LE.1) then
	  TheoryMax(XMonth) =  8.0 * acos (Numer/Denom)
	else if (Numer.GT.0) then
	  TheoryMax(XMonth) =  0.0
	else
	  TheoryMax(XMonth) = 24.0
	end if
      end do
      
      call GetMonthLengths (AYearAD,MonthLengths)
      
      do XAYear=1,NAYear
        do XMonth=1,NMonth
	  if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
	    if (TheoryMax(XMonth).GT.0) then
	      Numer=real(DataA(XAYear,XMonth,XAStn))/10.0
	      Denom=TheoryMax(XMonth)*real(MonthLengths(XAYear,XMonth))
	    
	      DataA(XAYear,XMonth,XAStn) = nint(1000.0 * Numer / Denom)
	    
	      if      (DataA(XAYear,XMonth,XAStn).LT.0.OR. &
	    		DataA(XAYear,XMonth,XAStn).GT.1010) then
	        DataA   (XAYear,XMonth,XAStn) = DataMissVal
	        if (QSource.EQ.2) DataASrc(XAYear,XMonth,XAStn) = DataMissVal
	      else if (DataA(XAYear,XMonth,XAStn).GT.1000) then
	        DataA(XAYear,XMonth,XAStn) = 1000
	      end if
	    else
	      DataA   (XAYear,XMonth,XAStn) = DataMissVal
	      if (QSource.EQ.2) DataASrc(XAYear,XMonth,XAStn) = DataMissVal
	    end if
	  end if
	end do
      end do
      
!      write (99,"(i7,x,a20,f8.2,12f6.1)"), AStn(XAStn),StnNameA(XAStn),ALat(XAStn), &
!				(TheoryMax(XMonth),XMonth=1,NMonth)
    else				! if no lat cannot calc theoretical max snh
      AStn(XAStn)=MissVal
    end if
  end do
  
  deallocate (MonthLengths, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertVariableCts: dealloc: ML #####"
end if

call MakeStnInfo (StnInfoC,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=AYearAD,Data=DataA)
call SaveCTS (StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD,CallFile=SaveFileA)
if (QSource.EQ.2) call SaveCTS &
	(StnInfoC,StnLocalA,StnNameA,StnCtyA,Data=DataASrc,YearAD=AYearAD,CallFile=SaveFileASrc)

end subroutine ConvertVariableCts

!*******************************************************************************
! option 34

subroutine PruneSource

print*, "  > Select the .cts file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".cts")
SubBeg = index(LoadFileA,".cts",.TRUE.)
LoadFileB = LoadFileA ; LoadFileB(SubBeg:(SubBeg+3)) = ".src"

print*, "  > Select the .cts file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".cts")
SubBeg = index(SaveFileA,".cts",.TRUE.)
SaveFileB = SaveFileA ; SaveFileB(SubBeg:(SubBeg+3)) = ".src"

print*, "  > Select the source code to remove:"
do
	read (*,*,iostat=ReadStatus), SrcRemove
	if (ReadStatus.LE.0) exit
end do

print*, "  > Loading..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataASrc,YearAD=AYearAD, &
		CallFile=LoadFileB,Silent=1)
deallocate (StnInfoA,StnLocalA,StnNameA,StnCtyA,AYearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### PruneSource: dealloc load src #####"

call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,OldCode=AStnOld, &
		Lat=ALat,Lon=ALon,Elv=AElv, &
		Data=DataA,YearAD=AYearAD,CallFile=LoadFileA,Silent=1)
deallocate (StnInfoA,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### PruneSource: dealloc load cts #####"

NAYear = size(DataA,1) ; NMonth=12 ; NAStn = size(DataA,3) ; Log01=0

print*, "  > Processing..."
do XAYear=1,NAYear
  do XMonth=1,NMonth
    do XAStn=1,NAStn
      if (DataASrc(XAYear,XMonth,XAStn).EQ.SrcRemove) then
        DataASrc(XAYear,XMonth,XAStn) = DataMissVal
        DataA   (XAYear,XMonth,XAStn) = DataMissVal
	Log01=Log01+1
      end if
    end do
  end do
end do
print "(a,i10)", "   > Made changes to data: ", Log01

print*, "  > Saving..."
call MakeStnInfo (StnInfoA,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=AYearAD, &
			Data=DataA,Silent=1)
call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataA,YearAD=AYearAD, &
			CallFile=SaveFileA,Silent=1)
call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Data=DataASrc,YearAD=AYearAD, &
			CallFile=SaveFileB,Silent=1)

end subroutine PruneSource

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

subroutine Finish

print*

close (99)

end subroutine Finish

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

end program OpCRUts