! anomdtb.f90
! f90 program written by Tim Mitchell 
! originally written as part of opcruts.f90 (originated 11.02.02) as options 25
! program to convert CRU ts .dtb files to anomaly .txt files for use in gridding
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/anomdtb filenames.f90 time.f90 grimfiles.f90 
!	crutsfiles.f90 perfiles.f90 annfiles.f90 cetgeneral.f90 basicfun.f90 
!	wmokey.f90 gridops.f90 grid.f90 ctyfiles.f90 
!	./../cruts/anomdtb.f90 2> /tyn1/tim/scratch/stderr.txt

program AnomDTB

use FileNames
use Time
use GrimFiles
use CRUtsFiles
use PerFiles
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 Specify
call Anomalise
call Dumping
call Finish

contains

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

subroutine Intro

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

print*
print*, "  > ***** AnomDTB: converts .dtb to anom .txt for gridding *****"
print*

NMonth = 12 ; NVersion = 9 ; NSrc = 1000

end subroutine Intro

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

subroutine Specify

QHomogSource=0

print*, "  > Enter the suffix of the variable required:"
do
	read (*,*,iostat=ReadStatus), LoadSuffix
	if (ReadStatus.LE.0) then
		call CheckVariSuffix (LoadSuffix,Variable,Factor)
  		if (Variable.EQ."") then
  		  print*, "  > Variable unrecognised."
  		else
  		  if (LoadSuffix.EQ.".pre".OR.LoadSuffix.EQ.".wet") then
  			print*, "  > Will calculate percentage anomalies." ; QAnomPercent=1
  		  else
       			QAnomPercent=0
  		  end if
  		end if
	end if
	if (ReadStatus.LE.0.AND.Variable.NE."") exit
end do

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)

if (FileSuffix.EQ.".cts") then
 print*, "  > Use supplementary .nrm file? (1=no,2=yes) ?"
 do
	read (*,*,iostat=ReadStatus), QNRMfile
	if (ReadStatus.LE.0.AND.QNRMfile.GE.1.AND.QNRMfile.LE.2) exit
 end do
 if (QNRMfile.EQ.2) then
  print*, "  > Select the .nrm file to load:"
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  LoadFileB = LoadPath (GivenFile,".nrm")
 end if
else
 QNRMfile=1
end if

print*, "  > Specify the start,end of the normals period: "
do
	read (*,*,iostat=ReadStatus), YearAD0,YearAD1
	if (ReadStatus.LE.0.AND.YearAD1.GE.YearAD0) exit
end do

print*, "  > Specify the missing percentage permitted: "
do
	read (*,*,iostat=ReadStatus), MissThresh
	if (ReadStatus.LE.0.AND.MissThresh.GE.0.AND.MissThresh.LE.100) exit
end do
NormMin = ceiling((100.0-MissThresh)*(YearAD1-YearAD0+1)/100.0)
print*, "  > Data required for a normal: ", NormMin	! @@@@@@@@@@@@@@@@@@@@@@@

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

print*, "  > Select outputs (1=.cts,2=.ann,3=.txt,4=.stn): "
do
	read (*,*,iostat=ReadStatus), QOutput
	if (QOutput.EQ.0) then
	  print*, "  >   1: save to .cts and .src, summarise to .ann"
	  print*, "  >   2: save summary of stns to .ann only"
	  print*, "  >   3: save to pre-interpol .txt, summarise to .ann"
	  print*, "  >   4: save stns-in-range .stn, summarise to .ann"
	end if
	if (ReadStatus.LE.0.AND.QOutput.GE.1.AND.QOutput.LE.4) exit
end do

DistanceThresh=0				! DON'T ALLOW FOR THIS IN SAVE TO .CTS
						! it is only valid when divorced from station info
if (QOutput.GE.2) then	
  print*, "  > Check for duplicate stns after anomalising? (0=no,>0=km range)"
  do
	read (*,*,iostat=ReadStatus), DistanceThresh
	if (ReadStatus.LE.0.AND.DistanceThresh.GE.0) exit
  end do
end if

if (QOutput.EQ.1) then
  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")
else if (QOutput.EQ.2) then
  print*, "  > Select the .ann file to save:"
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  SaveFileC = SavePath (GivenFile,".ann")
else if (QOutput.EQ.3) then
  print*, "  > Select the generic .txt file to save (yy.mm=auto):"
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  SaveFileB = SavePath (GivenFile,".txt")
else if (QOutput.EQ.4) then
  print*, "  > Select the .stn file to save:"
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
  SaveFileB = SavePath (GivenFile,".stn")

  print*, "  > Enter the correlation decay distance:"
  do
	read (*,*,iostat=ReadStatus), DecayDistance
	if (ReadStatus.LE.0) exit
  end do

  print*, "  > Submit a grim that contains the appropriate grid."
  call GrabGrid (1,RefGrid,RefBounds,RefBoxes,Quiet=1)
  NExe=size(RefGrid,1) ; NWye=size(RefGrid,2)
end if

if (QOutput.EQ.3.OR.QOutput.EQ.4) then
  print*, "  > Select the first,last years AD to save: "
  do
	read (*,*,iostat=ReadStatus), SaveYearAD0,SaveYearAD1
	if (ReadStatus.LE.0.AND.SaveYearAD1.GE.SaveYearAD0) exit
  end do
end if

end subroutine Specify

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

subroutine Anomalise

print*, "  > Operating..."

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,silent=1) 	! get .cts file
  SuffixBeg=index(LoadFileA,".cts")
  LoadFileC=LoadFileA(1:SuffixBeg-1) // ".src"
  SuffixBeg=index(LoadFileC,"/",.TRUE.)
  call LoadCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,Data=DataC,YearAD=CYearAD, &
  			CallFile=LoadFileC,Silent=1)          			! get .src file
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,silent=1) 	! get .dtb file
  SuffixBeg=index(LoadFileA,".dtb")
  LoadFileC=LoadFileA(1:SuffixBeg-1) // ".dts"
  SuffixBeg=index(LoadFileC,"/",.TRUE.)
  call LoadCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,Data=DataC,YearAD=CYearAD, &
  			CallFile=LoadFileC,Silent=1,DtbNormals=DtbNormalsC)     ! get .dts file
  deallocate (DtbNormalsC,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Anomalise: Deallocation failure: DtbNormalsC #####"
end if

Loaded = count(DataA.NE.DataMissVal)
NAYear = size(DataA,1) ; NAStn = size(DataA,3) 
AStart = YearAD0 - AYearAD(1) + 1 ; AEnd = AStart + YearAD1 - YearAD0		! normals period
if (AStart.LT.1.AND.AEnd.GE.1) AStart = 1
if (AEnd.GT.NAYear.AND.AStart.LE.NAYear) AEnd = NAYear
if (AStart.LT.1.OR.AStart.GT.NAYear.OR.AEnd.LT.1.OR.AEnd.GT.NAYear) then
	AStart=-999 ; AEnd=-999
  	print*, "  > ##### ERROR: Anomalise: No valid normals period #####"
end if

if (QOutput.EQ.3) then								! output period
  BStart = SaveYearAD0 - AYearAD(1) + 1 ; BEnd = BStart + SaveYearAD1 - SaveYearAD0
  if (BStart.LT.1.AND.BEnd.GE.1) BStart = 1
  if (BEnd.GT.NAYear.AND.BStart.LE.NAYear) BEnd = NAYear
  if (BStart.LT.1.OR.BStart.GT.NAYear.OR.BEnd.LT.1.OR.BEnd.GT.NAYear) then
  	BStart=-999 ; BEnd=-999
  	print*, "  > ##### ERROR: Anomalise: No valid output period #####"
!  else
!  	print "(a,2i5)", "   > The output period is:", AYearAD(BStart),AYearAD(BEnd)	! @@@@@@@@@@@@@@
  end if
end if

allocate (NormMean  (NMonth,NAStn), &
          NormStdev (NMonth,NAStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Anomalise: Allocation failure: Norm #####"
NormMean=DataMissVal ; NormStdev=DataMissVal

print "(a,6x,a)", "   > NORMALS", "      MEAN percent      STDEV percent"
									! look for normals from .dtb
! if (FileSuffix.EQ.".dtb".AND.YearAD0.EQ.1961.AND.YearAD1.EQ.1990) call NormalsDtbLine
call NormalsDtbLine		! ########################

Log1=0 ; Log2=0
do XAStn = 1, NAStn           						! look for normals from .cts
  do XMonth = 1, NMonth
   if (NormMean(XMonth,XAStn).NE.DataMissVal) then
    OpEn=0 ; OpTot=0 ; OpTotSq=0
    do XAYear = 1, NAYear
          if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
            OpEn=OpEn+1
            OpTot=OpTot+DataA(XAYear,XMonth,XAStn)
            OpTotSq=OpTotSq+(DataA(XAYear,XMonth,XAStn)**2)
          end if
    end do
				
    if (OpEn.GE.5) then
          NormStDev(XMonth,XAStn) = &
        		sqrt((OpEn/(OpEn-1))*((OpTotSq/OpEn)-((OpTot/OpEn)**2)))
          do XAYear = 1, NAYear
            if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) Log2=Log2+1
          end do				
    end if
   else
    OpEn = 0.0 ; OpTot = 0.0 ; OpTotSq = 0.0
    if (AStart.NE.-999.AND.AEnd.NE.-999) then				! get totals for stn/mon
        do XAYear = AStart, AEnd
          if (DataA(XAYear,XMonth,XAStn).NE.-9999) then
          	OpEn = OpEn + 1.0
          	OpTot = OpTot + (real(DataA(XAYear,XMonth,XAStn))/Factor)
          	OpTotSq = OpTotSq + ((real(DataA(XAYear,XMonth,XAStn))/Factor) ** 2)
          end if
        end do
    end if
    
    if (OpEn.GE.NormMin) then						! calc av & stdev if possible
      NormMean  (XMonth,XAStn) = Factor*OpTot/OpEn
      if (OpTotSq.GT.0) NormStdev (XMonth,XAStn) = Factor*sqrt((OpTotSq/OpEn)-((OpTot/OpEn)**2))
      
      OpEn = 0.0 ; OpTot = 0.0 ; OpTotSq = 0.0
      if (AStart.NE.-999.AND.AEnd.NE.-999) then				! recalc normals, no outliers
        do XAYear = AStart, AEnd
          if (DataA(XAYear,XMonth,XAStn).NE.-9999.AND. &
          		abs(real(DataA(XAYear,XMonth,XAStn))-NormMean(XMonth,XAStn)).LE. &
          		(NormStdev(XMonth,XAStn)*StdevThresh)) then
          	OpEn = OpEn + 1.0
          	OpTot = OpTot + (real(DataA(XAYear,XMonth,XAStn))/Factor)
          	OpTotSq = OpTotSq + ((real(DataA(XAYear,XMonth,XAStn))/Factor) ** 2)
          end if
        end do
      end if
      
      if (OpEn.GE.NormMin) then
        NormMean  (XMonth,XAStn) = Factor*OpTot/OpEn
        if (OpTotSq.GT.0) then
          NormStdev (XMonth,XAStn) = Factor*sqrt((OpEn/(OpEn-1))*((OpTotSq/OpEn)-((OpTot/OpEn)**2)))
        else
          NormStdev (XMonth,XAStn) = DataMissVal
        end if
	
	do XAYear=1,NAYear
	  if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
	    Log1=Log1+1
	    if (NormStdev(XMonth,XAStn).NE.DataMissVal) Log2=Log2+1
	  end if
	end do
      else
        NormMean  (XMonth,XAStn) = DataMissVal
        NormStdev (XMonth,XAStn) = DataMissVal
      end if
    end if    
   end if
  end do
end do
print "(a,8x,a,2(x,i10,x,f7.1))", "   > ", ".cts", &
		Log1,(100*(real(Log1)/Loaded)), Log2,(100*(real(Log2)/Loaded))

if (QNRMfile.EQ.2) then
  call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Code=BStn,OldCode=BStnOld, &     ! get .nrm file
                NmlData=Normals,NmlInc=NmlInc,CallFile=LoadFileB,Silent=1)
  
  Log1=0
  NBStn=size(Normals,2) ; XBStn=1 ; XAStn=1                           	! look for .nrm normals
  do
    do    
      if (AStn(XAStn).GT.BStn(XBStn)) XBStn=XBStn+1
      if (XBStn.LE.NBStn) then
       if (AStn(XAstn).EQ.BStn(XBStn)) then
        if (StnInfoB(XBStn,5).GE.YearAD0.AND.StnInfoB(XBStn,6).LE.YearAD1.AND. &
             real((StnInfoB(XBStn,6)-StnInfoB(XBStn,5)+1))*(real(NmlInc(XBStn))/100.0) &
             .GE.real(NormMin)) then
                    QMisMatch = 0 ; XMonth = 0
                    do XMonth = 1, NMonth
                       if (Normals(XMonth,XBStn).NE.DataMissVal) then
                        OpEn=0 ; OpTot=0 ; OpTotSq=0
                        do XAYear = 1, NAYear
                         if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
                             OpEn = OpEn + 1
                             OpTot = OpTot + (DataA(XAYear,XMonth,XAStn)/Factor)
                             OpTotSq = OpTotSq + (DataA(XAYear,XMonth,XAStn)/Factor) ** 2
                         end if
                        end do
                        if (OpEn.GE.5) then					! difference of means test
                             OpStDev = Factor*sqrt((OpEn/(OpEn-1))*((OpTotSq/OpEn)-((OpTot/OpEn)**2)))
                             OpMean  = Factor*(OpTot/OpEn)
                             OpDiff  = OpMean - Normals(XMonth,XBStn)		! stdev: assume .nrm = .cts
                             if ((abs(OpDiff)/(OpStdev/2.0)).GT.3.0) QMisMatch=QMisMatch+1
                        end if
                       end if                                              
                    end do
                    
                    if (QMisMatch.LT.3) then
                      do XMonth = 1, NMonth
                        if (NormMean(XMonth,XAStn).EQ.DataMissVal.AND. &
				Normals(XMonth,XBStn).NE.DataMissVal) then
                    	    NormMean(XMonth,XAStn)=Normals(XMonth,XBStn)
			    
			    do XAYear=1,NAYear
			      if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) Log1=Log1+1
			    end do
                    	end if
                      end do
                    end if
        end if
        XBStn=XBStn+1
       end if
      end if
      if (XBStn.GT.NBStn) exit
      if (AStn(XAStn).LT.BStn(XBStn)) exit
    end do
    
    XAStn=XAStn+1
    if (XAStn.GT.NAStn) exit
    if (XBStn.GT.NBStn) exit
  end do
  print "(a,8x,a,(x,i10,x,f7.1))", "   > ", ".nrm", Log1,(100*(real(Log1)/Loaded))
end if

Log0=0 ; Log1=0 ; Log2=0 ; Log3=0 ; Log4=0 ; Log5=0
do XAStn = 1, NAStn           						! anomalise
  if (ALat(XAStn).NE.MissVal.OR.ALon(XAStn).NE.MissVal) then
   do XMonth = 1, NMonth
    do XAYear = 1, NAYear
      if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
        if (NormMean(XMonth,XAStn).NE.DataMissVal.AND.(QAnomPercent.EQ.0.OR. &
        			NormMean(XMonth,XAStn).NE.0)) then
          if (NormStdev(XMonth,XAStn).NE.DataMissVal) then
            Log5=Log5+1				! stdev check performed
	    if (abs(DataA(XAYear,XMonth,XAStn)-NormMean(XMonth,XAStn)).LE. &
            			(NormStdev(XMonth,XAStn)*StdevThresh)) then
              if (QAnomPercent.EQ.0) then
                DataA(XAYear,XMonth,XAStn) = DataA(XAYear,XMonth,XAStn) - NormMean(XMonth,XAStn)
              else
                DataA(XAYear,XMonth,XAStn) = nint(1000.0*((real(DataA(XAYear,XMonth,XAStn)) / &
                				real(NormMean(XMonth,XAStn)))-1.0))
              end if
              Log1=Log1+1			! accepted for moment
            else
              DataA(XAYear,XMonth,XAStn) = DataMissVal ; DataC(XAYear,XMonth,XAStn) = DataMissVal
              Log3=Log3+1			! outside stdev range
            end if
          else
            if (QAnomPercent.EQ.0) then
                DataA(XAYear,XMonth,XAStn) = DataA(XAYear,XMonth,XAStn) - NormMean(XMonth,XAStn)
            else
                DataA(XAYear,XMonth,XAStn) = nint(1000.0*((real(DataA(XAYear,XMonth,XAStn)) / &
                				real(NormMean(XMonth,XAStn)))-1.0))
            end if
            Log1=Log1+1				! accepted for moment
          end if
        else
          DataA(XAYear,XMonth,XAStn) = DataMissVal ; DataC(XAYear,XMonth,XAStn) = DataMissVal
          Log2=Log2+1				! no normal
        end if
      end if
      
      if (DataA(XAYear,XMonth,XAStn).GT.99999) DataA(XAYear,XMonth,XAStn) = 99998
      if (DataA(XAYear,XMonth,XAStn).LT.-9999) DataA(XAYear,XMonth,XAStn) = -9998
    end do
   end do
  else
   AStn(XAStn) = -999   
   
   do XMonth = 1, NMonth
    do XAYear = 1, NAYear
     if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
       Log0=Log0+1				! no lat/lon
       DataA(XAYear,XMonth,XAStn) = DataMissVal ; DataC(XAYear,XMonth,XAStn) = DataMissVal
     end if
    end do
   end do
  end if
end do

print "(a,6x,a,x,a)", "   > PROCESS", "  DECISION percent %of-chk"
print "(a,x,i10,2(x,f7.1))", "   > no lat/lon  ", &
	Log0,(100*(real(Log0)/Loaded)),(100*(real(Log0)/Loaded))
print "(a,x,i10,2(x,f7.1))", "   > no normal   ", &
	Log2,(100*(real(Log2)/Loaded)),(100*(real(Log2)/(Loaded-real(Log0))))
print "(a,x,i10,2(x,f7.1))", "   > out-of-range", &
	Log3,(100*(real(Log3)/Loaded)),(100*(real(Log3)/real(Log5)))

if (DistanceThresh.GT.0) then			! merge duplicate stations, judging by distance
  do XAStn=1,NAStn				! this is toggled distance/omit for save option 3
   if (AStn(XAStn).NE.MissVal) then
     do XBStn=XAStn+1,NAStn
      if (AStn(XBStn).NE.MissVal) then
        if (GetDistance (ALat(XAStn),ALon(XAStn),ALat(XBStn),ALon(XBStn)).LT.DistanceThresh) then
	  write (99,"(i7,a,i7)"), AStn(XAStn), " duplicates ", AStn(XBStn)
          AStn(XAStn) = -999

          do XAYear=1,NAYear
            do XMonth=1,NMonth
              if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) then
                if (DataA(XAYear,XMonth,XBStn).EQ.DataMissVal) then
                  DataA(XAYear,XMonth,XBStn) = DataA(XAYear,XMonth,XAStn)
              	else
              	  Log4=Log4+1
              	end if
              end if
            end do
          end do
        end if
      end if
     end do
   end if
  end do

  print "(a,x,i10,2(x,f7.1))", "   > duplicated  ", &
	Log4,(100*(real(Log4)/Loaded)),(100*(real(Log4)/real(Log1)))
  Log1=Log1-Log4
end if
print "(a,x,i10,(x,f7.1))", "   > accepted    ", &
	Log1,(100*(real(Log1)/Loaded))

end subroutine Anomalise

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

subroutine Dumping

if (QOutput.EQ.1) then				! dump to .cts and .src
  call MakeStnInfo (StnInfoB,AStn,AStnOld,ALat,ALon,AElv,1,YearAD=AYearAD,Data=DataA)
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,YearAD=AYearAD,Data=DataA, &
  		CallFile=SaveFileB,Silent=1)

  do XAYear=1,NAYear
    do XMonth=1,NMonth
      do XAStn=1,NAStn
        if (DataA(XAYear,XMonth,XAStn).EQ.-9999) DataC(XAYear,XMonth,XAStn)=-9999
      end do
    end do
  end do
  
  SuffixBeg=index(SaveFileB,".cts")
  SaveFileC=SaveFileB(1:SuffixBeg-1) // ".src"
  SuffixBeg=index(SaveFileC,"/",.TRUE.)
  print "(2a)", "   > Saving source file:   ", trim(SaveFileC(SuffixBeg+1:80))
  call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,YearAD=AYearAD,Data=DataC, &
  		CallFile=SaveFileC,Silent=1)  
else if (QOutput.EQ.3) then			! dump to .txt
  print "(a,2(i4,a))", "   > Dumping years ", AYearAD(BStart), "-", AYearAD(BEnd), " to .txt files..."
  SuffixBeg=index(SaveFileB,".txt")
  do XAYear=BStart,BEnd
    YearTxt = GetTextFromInt (AYearAD(XAYear))
    do XMonth=1,NMonth
       GivenFile = SaveFileB(1:SuffixBeg-1) // "." // trim(YearTxt) // "." // MonthTxt(XMonth) // ".txt"
       open  (9,file=GivenFile,status="replace",access="sequential",form="formatted",action="write")
       do XAStn = 1, NAStn
         if (DataA(XAYear,XMonth,XAStn).NE.-9999) write (9,"(2f8.2,f8.1,f13.5,i7)"), &
           ALat(XAStn),ALon(XAStn),AElv(XAStn),real(DataA(XAYear,XMonth,XAStn))*Factor,AStn(XAStn)
       end do
       close (9)
    end do
  end do  
else if (QOutput.EQ.4) then			! dump to .stn
  print*, "  > Calculating station coverages..."
  NBYear=SaveYearAD1-SaveYearAD0+1
  
  allocate (RefGridLat(NExe,NWye),RefGridLon(NExe,NWye), &
	    BYearAD(NBYear),RealData(NBYear,NMonth,RefBoxes), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### WithinRange: Alloc: DataB #####"
  RefGridLat=MissVal ; RefGridLon=MissVal ; RealData=0
  ExeSpace=(RefBounds(2)-RefBounds(1))/real(NExe)
  WyeSpace=(RefBounds(4)-RefBounds(3))/real(NWye)

  do XBYear = 1, NBYear
    BYearAD(XBYear) = SaveYearAD0+XBYear-1
  end do
  call CommonVecPer (AYearAD,BYearAD,AStart,AEnd,BStart,BEnd)

  do XExe=1,NExe
   do XWye=1,NWye
    if (RefGrid(XExe,XWye).GT.0) then
      RefGridLon(XExe,XWye) = RefBounds(1)+((real(XExe)-0.5)*ExeSpace)
      RefGridLat(XExe,XWye) = RefBounds(3)+((real(XWye)-0.5)*WyeSpace)
    end if
   end do
  end do
  
  do XAStn=1,NAStn
    XMonth=0 ; ValidNorm=.FALSE.
    do
      XMonth=XMonth+1
      if (NormMean(XMonth,XAStn).NE.DataMissVal) ValidNorm=.TRUE.
      if (ValidNorm.OR.XMonth.EQ.12) exit
    end do
    
    if (ValidNorm) then
     do XExe = 1, NExe
      do XWye = 1, NWye
       if (RefGrid(XExe,XWye).GT.0) then
        Distance = GetDistance(ALat(XAStn),ALon(XAStn), &
        		RefGridLat(XExe,XWye),RefGridLon(XExe,XWye))	! gridops.f90
        if (Distance.LT.DecayDistance) then
          do XAYear = AStart,AEnd
            XBYear=(XAYear-AStart)+BStart
            do XMonth = 1, NMonth
              if (DataA(XAYear,XMonth,XAStn).NE.-9999) &
              	  RealData(XBYear,XMonth,RefGrid(XExe,XWye)) = &
              	  RealData(XBYear,XMonth,RefGrid(XExe,XWye)) + 1
            end do
          end do
        end if
       end if
      end do
     end do
    end if
  end do
  
  print*, "  > Saving..."
  Title="station count"
  call SaveGrim (RealData,RefGrid,BYearAD,RefBounds,Title,SaveFileB,".stn", &
  		 SaveSuffix,NoZip=1,Silent=1)
end if

if (QOutput.EQ.1) then
    SuffixBeg=index(SaveFileB,".cts")
    SaveFileC=SaveFileB(1:SuffixBeg-1) // ".ann"
else if (QOutput.EQ.3) then
    SuffixBeg=index(SaveFileB,".txt")
    SaveFileC=SaveFileB(1:SuffixBeg-1) // ".ann"
else if (QOutput.EQ.4) then
    SuffixBeg=index(SaveFileB,".stn")
    SaveFileC=SaveFileB(1:SuffixBeg-1) // ".ann"
end if
  
SuffixBeg=index(SaveFileC,"/",.TRUE.)
! print "(2a)", "   > Saving summary file:  ", trim(SaveFileC(SuffixBeg+1:80))

call ClassifyContinent
call SaveANN (SaveFileC,AYearAD,ColTitles,SumsPer,DecPlaceN=1)

end subroutine Dumping

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

subroutine Finish

print*

close (99)

end subroutine Finish

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

subroutine NormalsDtbLine

Log1=0
do XAStn = 1, NAStn
  do XMonth = 1, NMonth
    if (NormMean(XMonth,XAStn).EQ.DataMissVal) then
      if (DtbNormalsA(XMonth,XAStn).NE.DataMissVal) then
        NormMean(XMonth,XAStn) = DtbNormalsA(XMonth,XAStn)
        
        do XAYear = 1, NAYear
          if (DataA(XAYear,XMonth,XAStn).NE.DataMissVal) Log1=Log1+1
        end do				
      end if
    end if
  end do
end do

print "(a,8x,a,(x,i10,x,f7.1))", "   > ", ".dtb", &
		Log1,(100*(real(Log1)/Loaded))

end subroutine NormalsDtbLine

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

subroutine ClassifyContinent

NCty=MasterSize()						! ctyfiles.f90
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)	! ctyfiles.f90

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 

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

end program AnomDTB