! qcclimatt.f90
! program to Quality Control the CLIMAT extra Temperature data supplied directly from Hadley
! written by Dr. Tim Mitchell (Tyndall Centre) on 21.01.02
! last modified on 22.02.02
! f90 -o ./../obs/qcclimatt filenames.f90 crutsfiles.f90 wmokey.f90 ./../obs/qcclimatt.f90

program QCClimatT

use FileNames
use CRUtsFiles
use WMOkey

implicit none

real, pointer, dimension (:,:,:)		:: Raw
real, pointer, dimension (:,:)			:: PhilMeans,PhilSigma
real, pointer, dimension (:)			:: Lat,Long,Elev
real, dimension (7)				:: ClassBounds = [0.005,0.105,0.305,0.505,0.805,0.995,1.005]

integer, pointer, dimension (:,:,:)		:: DTR
integer, pointer, dimension (:,:,:)		:: CtyRaw,CtyQC,RawQC,RegQC,BitQC
integer, pointer, dimension (:,:)		:: Station,RegInfo,CheckRaw,SumQC,StnQC,StnInfo,RawTime,PreQCTime
integer, pointer, dimension (:)			:: CheckStn,CtyBound0,CtyBound1,CtyReg
integer, pointer, dimension (:)			:: PhilMeansStn,PhilMeansYr0,PhilMeansYr1
integer, pointer, dimension (:)			:: PhilSigmaStn,PhilSigmaYr0,PhilSigmaYr1
integer, pointer, dimension (:)			:: FileLines,YearAD,Code
integer, dimension (16)				:: LineInfo,LastInfo
integer, dimension (12)				:: PhilInfo

character (len=80), pointer, dimension (:)	:: ClimatTFiles,StepText,BitText
character (len=20), pointer, dimension (:)	:: NameStn
character (len=13), pointer, dimension (:)	:: NameCty,CtyName
character (len=02), pointer, dimension (:)	:: CtyAcro

real, parameter 		:: FileMissVal = -32768.0
real, parameter 		:: DataMissVal = -9999.0
real, parameter 		:: MissVal = -999.0

character (len=80), parameter	:: Blank = "  "

real :: AvMinMax,OpFraction
real :: KSigmas,KPJThresh

integer :: ReadStatus, AllocStat
integer :: RawFileN, LineN, StationN, CheckStnN, CtyN, InfoN, MonthN, YearN, RegN, StepN, BitN
integer :: XRawFile, XLine, XStation, XCheckStn, XCty, XInfo, XMonth, XYear, XReg, XStep, XBit
integer :: PhilMeansN,PhilSigmaN, ClassN
integer :: XPhilMeans,XPhilSigma, XClass
integer :: ThisCty, OpTot, OpMiss
integer :: WMOStation, WMOBlock, WMOCode
integer :: QRepeatEntry,QRepeatMonth,QRepeatReport,QOneMissVal

character (len=80) :: ClimatTFilter = "/cru/mikeh1/f709762/obs/climat/_raw/????.dat"
character (len=80) :: PhilMeansFile = "/cru/mikeh1/f709762/obs/philj/pjones-stn-tmean-norms.dat"
character (len=80) :: PhilSigmaFile = "/cru/mikeh1/f709762/obs/philj/pjones-stn-tmean-sigma.dat"
character (len=80) :: Trash

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

call Intro
call GetCtyBounds
call GetPhilJones
call GetFileSpecs
call AllocArrays

print*, "  > Loading and QC'ing data..."
do XRawFile = 1, RawFileN
  LineN = FileLines(XRawFile)
  call LoadRaw
  call QualityControl
  call SummStn
  call EndRawFile
end do

print*, "  > Summarising results..."
call DumpRawInfo
call DumpQCInfo
call DumpData
call Conclude

contains

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

subroutine Intro

open (99,file="/cru/mikeh1/f709762/scratch/log/log-qcclimatt.dat",status="replace",action="write")

print*
print*, "  > ***** QCClimatT: QC on extra CLIMAT Temperature data *****"
print*

KSigmas = 5.0 ; KPJThresh = 1.0 * 10.0			! KPJThresh is *10 to allow for degC*10
QOneMissVal = 1						! 0=allow -8000series missing values, 1=-9999 only

call GetBatch (ClimatTFilter,ClimatTFiles)		! locate raw extra CLIMAT T data files
RawFileN = size (ClimatTFiles,1)

allocate (YearAD(RawFileN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure #####"
do XRawFile = 1, RawFileN
  YearAD(XRawFile) = 1995 + XRawFile - 1
end do

MonthN = 12 ; RegN = 12 ; StepN = 11 ; BitN = 16 ; ClassN = 7

end subroutine Intro

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

subroutine GetCtyBounds

print*, "  > Getting country information..."

call LoadWMOCty (CtyBound0,CtyBound1,CtyAcro,CtyReg,CtyName,ExtraCty=1)
CtyN = size(CtyBound0) - 1
CtyName(CtyN+1) = "unidentified" ; CtyReg(CtyN+1) = 12

allocate (CtyRaw (RawFileN,MonthN,CtyN+1), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetCtyBounds: Allocation failure: CtyRaw #####"
CtyRaw=0

end subroutine GetCtyBounds

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

subroutine GetPhilJones

print*, "  > Getting Phil Jones' information..."

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

allocate (PhilMeansYr0 (PhilMeansN),   &
	  PhilMeansYr1 (PhilMeansN),   &
	  PhilMeansStn (PhilMeansN),   &
	  PhilMeans    (PhilMeansN,12),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetPhilJones: Allocation failure: PhilMeans* #####"

open  (3,file=PhilMeansFile,status="old",access="sequential",form="formatted",action="read")
do XPhilMeans = 1, PhilMeansN
  read (3,"(i6,a10,2i5,12i5)"),PhilMeansStn(XPhilMeans),Trash,PhilMeansYr0(XPhilMeans), &
  			PhilMeansYr1(XPhilMeans), (PhilInfo(XInfo),XInfo=1,MonthN)
  do XMonth = 1, MonthN
    PhilMeans(XPhilMeans,XMonth) = real(PhilInfo(XMonth))
  end do
end do
close (3)

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

allocate (PhilSigmaYr0 (PhilSigmaN),   &
	  PhilSigmaYr1 (PhilSigmaN),   &
	  PhilSigmaStn (PhilSigmaN),   &
	  PhilSigma    (PhilSigmaN,12),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetPhilJones: Allocation failure: PhilSigma* #####"

open  (3,file=PhilSigmaFile,status="old",access="sequential",form="formatted",action="read")
do XPhilSigma = 1, PhilSigmaN
  read (3,"(i6,a10,2i5,12i5)"),PhilSigmaStn(XPhilSigma),Trash,PhilSigmaYr0(XPhilSigma), &
  			PhilSigmaYr1(XPhilSigma), (PhilInfo(XInfo),XInfo=1,MonthN)
  do XMonth = 1, MonthN
    PhilSigma(XPhilSigma,XMonth) = real(PhilInfo(XMonth))
  end do
end do
close (3)

end subroutine GetPhilJones

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

subroutine GetFileSpecs

print*, "  > Identifying file specifications..."

CheckStnN = 1000000
allocate (CheckStn  (CheckStnN),  &
	  FileLines (RawFileN),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFileSpecs: Allocation failure: CheckStn #####"
CheckStn = 0 ; FileLines = 0 ; StationN = 0

do XRawFile = 1, RawFileN
  call system ('wc -l ' // trim(ClimatTFiles(XRawFile)) // ' > trashme-gfs.txt')
  open  (3,file='trashme-gfs.txt',status="old",access="sequential",form="formatted",action="read")
  read  (3,"(i10)"), FileLines(XRawFile)			! get number of lines
  close (3)
  call system ('rm trashme-gfs.txt')
  
  open  (3,file=ClimatTFiles(XRawFile),status="old",access="sequential",form="formatted",action="read")
  do XLine = 1, FileLines(XRawFile)
    read (3,"(a12,2i6)"), Trash, WMOBlock, WMOStation
    WMOCode = (WMOBlock * 10000) + WMOStation
    CheckStn(WMOCode) = CheckStn(WMOCode) + 1			! store existence of station/month
    if (CheckStn(WMOCode).EQ.1) StationN = StationN + 1
  end do
  close (3)
end do
print*, "  > WMO stations in files total: ", StationN

allocate (Station (StationN,4), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFileSpecs: Allocation failure: Station #####"
Station = MissVal		! 1=WMO-stn,2=cty-index,3=PhilJ-mean-index,4=PhilJ-sigma-index

XStation = 0
do XCheckStn = 1, CheckStnN
  if (CheckStn(XCheckStn).GT.0) then
    XStation = XStation + 1

    Station(XStation,1) = XCheckStn

    call LocateStn
    Station(XStation,2) = ThisCty
    
    XPhilMeans = 0
    do
      XPhilMeans = XPhilMeans + 1
      if (PhilMeansStn(XPhilMeans).EQ.XCheckStn) Station(XStation,3) = XPhilMeans      
      if (PhilMeansStn(XPhilMeans).GT.XCheckStn) exit
      if (XPhilMeans.EQ.PhilMeansN) exit
    end do
    
    XPhilSigma = 0
    do
      XPhilSigma = XPhilSigma + 1
      if (PhilSigmaStn(XPhilSigma).EQ.XCheckStn) Station(XStation,4) = XPhilSigma      
      if (PhilSigmaStn(XPhilSigma).GT.XCheckStn) exit
      if (XPhilSigma.EQ.PhilSigmaN) exit
    end do
  end if
end do

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

end subroutine GetFileSpecs

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

subroutine AllocArrays

print*, "  > Allocating arrays..."

allocate (RegQC     (RawFileN,RegN,StepN),    &
	  BitQC     (RawFileN,StepN,0:BitN-1),&
	  StepText  (StepN),                  &
	  BitText   (StepN),                  &
	  CheckRaw  (RawFileN,5),             &
	  StnQC     (RawFileN+1,ClassN+1),    &
	  RawTime   (RawFileN,MonthN),	      &
	  PreQCTime (RawFileN,MonthN+1),      &
	  DTR       (RawFileN,MonthN,StationN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure #####"
RegQC=0 ; BitQC=0 ; CheckRaw=0 ; StepText=Blank ; BitText=Blank ; StnQC = 0 ; DTR = DataMissVal
RawTime=0 ; PreQCTime=0

end subroutine AllocArrays

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

subroutine LoadRaw

allocate (Raw (MonthN,StationN,16), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadRaw: Allocation failure: Raw #####"
Raw = MissVal

open  (3,file=ClimatTFiles(XRawFile),status="old",access="sequential",form="formatted",action="read")
XStation = 1
do XLine = 1, LineN
  read (3,"(20i6)"), XYear,XMonth,WMOBlock,WMOStation,(LineInfo(XInfo),XInfo=1,16)
  CheckRaw (XRawFile,1) = CheckRaw (XRawFile,1) + 1
  
  if (XMonth.GE.1.AND.XMonth.LE.12) then
   PreQCTime (XRawFile,XMonth) = PreQCTime (XRawFile,XMonth) + 1
   WMOCode = (WMOBlock * 10000) + WMOStation		! identify WMO code
   do							! identify XStation
    if (WMOCode.EQ.Station(XStation,1)) exit
    if (WMOCode.NE.Station(XStation,1)) XStation = XStation + 1
    if (XStation.GT.StationN) XStation = 1
   end do
  
   if (Raw(XMonth,XStation,1).NE.MissVal) then
   	CheckRaw (XRawFile,3) = CheckRaw (XRawFile,3) + 1
   	QRepeatEntry = 0
   	QRepeatReport = 1
   else
   	QRepeatEntry = 1 ; QRepeatReport = 0 ;  XInfo = 0
   	do
     	  XInfo = XInfo + 1
     	  if (LineInfo(XInfo).NE.LastInfo(XInfo)) QRepeatEntry = 0
     	  if (XInfo.EQ.16.OR.QRepeatEntry.EQ.0) exit
   	end do
   end if
    
   if (XMonth.GT.1) then
     	  QRepeatMonth = 1 ; XInfo = 0
     	  do
      	    XInfo = XInfo + 1
      	    if (real(LineInfo(XInfo)).NE.Raw((XMonth-1),XStation,XInfo)) QRepeatMonth = 0
      	    if (XInfo.EQ.16.OR.QRepeatMonth.EQ.0) exit
     	  end do
   else
     	  QRepeatMonth = 0
   end if
    
   if (QRepeatEntry.EQ.0.AND.QRepeatMonth.EQ.0) then
      do XInfo = 1, 16					! store data for station-month
        Raw (XMonth,XStation,XInfo) = real(LineInfo(XInfo))
      end do
  							! register station-month
      if (QRepeatReport.EQ.0) then
        RawTime(XRawFile,XMonth) = RawTime(XRawFile,XMonth) + 1
      	CtyRaw (XRawFile,XMonth,Station(XStation,2)) = CtyRaw (XRawFile,XMonth,Station(XStation,2)) + 1
      	CheckRaw (XRawFile,5) = CheckRaw (XRawFile,5) + 1
      end if
   else
      if (QRepeatReport.EQ.0) CheckRaw (XRawFile,4) = CheckRaw (XRawFile,4) + 1
   end if
  else
   PreQCTime (XRawFile,13) = PreQCTime (XRawFile,13) + 1
   CheckRaw (XRawFile,2) = CheckRaw (XRawFile,2) + 1
  end if
  
  LastInfo = LineInfo
end do
close (3)

allocate (RawQC (MonthN,StationN,StepN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadRaw: Allocation failure: RawQC #####"
RawQC = MissVal

end subroutine LoadRaw

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

subroutine SummStn

do XStation = 1, StationN
  OpTot = 0 ; OpMiss = 0
  
  do XMonth = 1, MonthN
    if (RawQC(XMonth,XStation,1).NE.MissVal) then
     OpTot = OpTot + 1

!     if (RawQC(XMonth,XStation,1).GT.0.OR.RawQC(XMonth,XStation,3).GT.0.OR.RawQC(XMonth,XStation,4).GT.0.OR.&
!      	 RawQC(XMonth,XStation,7).GT.0.OR.RawQC(XMonth,XStation,9).GT.0.OR.RawQC(XMonth,XStation,11).GT.0) then
!        	OpMiss = OpMiss + 1
!     else
!       		DTR (XRawFile,XMonth,XStation) = nint((Raw(XMonth,XStation,9)-Raw(XMonth,XStation,13)))
!     end if

     if      (RawQC(XMonth,XStation, 1).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8100 - RawQC(XMonth,XStation, 1)
     else if (RawQC(XMonth,XStation, 3).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8200 - RawQC(XMonth,XStation, 3)
     else if (RawQC(XMonth,XStation, 4).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8300 - RawQC(XMonth,XStation, 4)
     else if (RawQC(XMonth,XStation, 7).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8400 - RawQC(XMonth,XStation, 7)
     else if (RawQC(XMonth,XStation, 9).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8500 - RawQC(XMonth,XStation, 9)
     else if (RawQC(XMonth,XStation,11).GT.0) then
       		DTR (XRawFile,XMonth,XStation) = -8600 - RawQC(XMonth,XStation,11)
     else
       		DTR (XRawFile,XMonth,XStation) = nint((Raw(XMonth,XStation,9)-Raw(XMonth,XStation,13)))
     end if
     
     if (QOneMissVal.EQ.1.AND.DTR(XRawFile,XMonth,XStation).LT.-8000) DTR(XRawFile,XMonth,XStation)=DataMissVal
    end if
  end do
  
  if (OpTot.GT.0) then
    if (OpMiss.GT.0) then
    	OpFraction = real(OpMiss) / real(OpTot)
    else
    	OpFraction = 0.0
    end if
    
    XClass = 0
    do
      XClass = XClass + 1
      if (OpFraction.LT.ClassBounds(XClass)) StnQC(XRawFile,XClass) = StnQC(XRawFile,XClass) + 1
      if (OpFraction.LT.ClassBounds(XClass)) exit
      if (XClass.EQ.ClassN) exit
    end do
  end if
end do

end subroutine SummStn

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

subroutine EndRawFile

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

end subroutine EndRawFile

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

subroutine QualityControl

do XMonth = 1, MonthN
 do XStation = 1, StationN
  if (Raw(XMonth,XStation,1).NE.MissVal) then
   call Step1
   StepText(1) = "rejected if any of Tmean,Tmin,Tmax are missing"
   BitQC(XRawFile,1,RawQC(XMonth,XStation,1)) = BitQC(XRawFile,1,RawQC(XMonth,XStation,1)) + 1
   BitText(1) = "missing obs: 1=mean,2=min,4=max [fatal]"
   
   if (RawQC(XMonth,XStation,1).GT.0) then
    RegQC(XRawFile,CtyReg(Station(XStation,2)),1) = RegQC(XRawFile,CtyReg(Station(XStation,2)),1) + 1
   else
    call Step2
    StepText(2) = "no rejections"
    BitQC(XRawFile,2,RawQC(XMonth,XStation,2)) = BitQC(XRawFile,2,RawQC(XMonth,XStation,2)) + 1
    BitText(2) = "missing days: 1=mean,2=min,4=max [not fatal]"
    
    if (RawQC(XMonth,XStation,2).LT.7) then
     call Step3
     StepText(3) = "rejected if any 'days count' for Tmean,Tmin,Tmax is present and <15"
     BitQC(XRawFile,3,RawQC(XMonth,XStation,3)) = BitQC(XRawFile,3,RawQC(XMonth,XStation,3)) + 1
     BitText(3) = "days<15: 1=mean,2=min,4=max [>0 is fatal]"
     
     if (RawQC(XMonth,XStation,3).GT.0) then
      RegQC(XRawFile,CtyReg(Station(XStation,2)),3) = RegQC(XRawFile,CtyReg(Station(XStation,2)),3) + 1
     end if
    end if
    
    if (RawQC(XMonth,XStation,3).LE.0) then
     call Step4
     StepText(4) = "rejected unless min<mean<max and abs[((Tmin+Tmax)/2)-Tmean]=<2.0degC"
     BitQC(XRawFile,4,RawQC(XMonth,XStation,4)) = BitQC(XRawFile,4,RawQC(XMonth,XStation,4)) + 1
     BitText(4) = "failed tests: 1:min<mean<max 2:abs[((Tmin+Tmax)/2)-Tmean]=<2.0degC [fatal]"
      
     if (RawQC(XMonth,XStation,4).GT.0) then
      RegQC(XRawFile,CtyReg(Station(XStation,2)),4) = RegQC(XRawFile,CtyReg(Station(XStation,2)),4) + 1
     else
      call Step5       
      StepText(5) = "no rejections"
      BitQC(XRawFile,5,RawQC(XMonth,XStation,5)) = BitQC(XRawFile,5,RawQC(XMonth,XStation,5)) + 1
      BitText(5) = "missing normals: 1=mean,2=min,4=max,8=sigma [not fatal]"

      if (RawQC(XMonth,XStation,5).LT.7) then
       call Step6        
       StepText(6) = "no rejections"
       BitQC(XRawFile,6,RawQC(XMonth,XStation,6)) = BitQC(XRawFile,6,RawQC(XMonth,XStation,6)) + 1
       BitText(6) = "normals<15yr: 1=mean,2=min,4=max,8=period [not fatal]"

       if (RawQC(XMonth,XStation,6).LT.7) then
        call Step7         
        StepText(7) = "rejected unless (normal-?*sigma)<Tstat<(normal+?*sigma) (CLIMAT+,?=KSigmas)"
        BitQC(XRawFile,7,RawQC(XMonth,XStation,7)) = BitQC(XRawFile,7,RawQC(XMonth,XStation,7)) + 1
        BitText(7) = "failed test: (X-?s)<Tstat<(X+?s):1=mean,2=min,4=max (CLIMAT+,?=KSigmas))[fatal]"

        if (RawQC(XMonth,XStation,7).GT.0) then
          RegQC(XRawFile,CtyReg(Station(XStation,2)),7) = RegQC(XRawFile,CtyReg(Station(XStation,2)),7) + 1
        end if
       end if
      end if
       
      if (RawQC(XMonth,XStation,7).LE.0) then
       call Step8
       StepText(8) = "no rejections"
       BitQC(XRawFile,8,RawQC(XMonth,XStation,8)) = BitQC(XRawFile,8,RawQC(XMonth,XStation,8)) + 1
       BitText(8) = "PhilJ v CLIMAT+ normals: no comparison for 1=mean,2=sigma [not fatal]"
       
       if (RawQC(XMonth,XStation,8).LT.3) then
        call Step9
        StepText(9) = "rejected if PhilJ minus CLIMAT+ normals > KPJThresh: 1=mean,2=sigma [fatal]"
        BitQC(XRawFile,9,RawQC(XMonth,XStation,9)) = BitQC(XRawFile,9,RawQC(XMonth,XStation,9)) + 1
        BitText(9) = "failed normals test: abs[PhilJ--CLIMAT+]>KPJThresh: 1=mean,2=sigma [fatal]"
        
        if (RawQC(XMonth,XStation,9).GT.0) then
          RegQC(XRawFile,CtyReg(Station(XStation,2)),9) = RegQC(XRawFile,CtyReg(Station(XStation,2)),9) + 1
        end if
       end if
       
       if (RawQC(XMonth,XStation,9).LE.0) then
        call Step10
        StepText(10) = "no rejections"
        BitQC(XRawFile,10,RawQC(XMonth,XStation,10)) = BitQC(XRawFile,10,RawQC(XMonth,XStation,10)) + 1
        BitText(10) = "no test with PhilJ sigma: (X-?s)<Tstat<(X+?s): 1=mean,2=min,4=max [not fatal]"
       
        if (RawQC(XMonth,XStation,10).LT.7) then
         call Step11
         StepText(11) = "rejected unless (X-?s)<Tstat<(X+?s) (with PhilJ sigma: ?=KSigmas)"
         BitQC(XRawFile,11,RawQC(XMonth,XStation,11)) = BitQC(XRawFile,11,RawQC(XMonth,XStation,11)) + 1
         BitText(11) = "failed PhilJ sigma test: (X-?s)<Tstat<(X+?s): 1=mean,2=min,4=max [fatal]"
        
         if (RawQC(XMonth,XStation,11).GT.0) then
          RegQC(XRawFile,CtyReg(Station(XStation,2)),11) = RegQC(XRawFile,CtyReg(Station(XStation,2)),11) + 1
         end if         
        end if
       end if
      end if
     end if
    end if
   end if
  end if
 end do
end do

end subroutine QualityControl

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

subroutine Step1

RawQC(XMonth,XStation,1) = 0
    
if (Raw(XMonth,XStation, 4).EQ.FileMissVal) RawQC(XMonth,XStation,1) = RawQC(XMonth,XStation,1) +  1
if (Raw(XMonth,XStation,13).EQ.FileMissVal) RawQC(XMonth,XStation,1) = RawQC(XMonth,XStation,1) +  2
if (Raw(XMonth,XStation, 9).EQ.FileMissVal) RawQC(XMonth,XStation,1) = RawQC(XMonth,XStation,1) +  4

end subroutine Step1

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

subroutine Step2

RawQC(XMonth,XStation,2) = 0

if (Raw(XMonth,XStation, 5).EQ.FileMissVal) RawQC(XMonth,XStation,2) = RawQC(XMonth,XStation,2) + 1
if (Raw(XMonth,XStation,10).EQ.FileMissVal) RawQC(XMonth,XStation,2) = RawQC(XMonth,XStation,2) + 2
if (Raw(XMonth,XStation,10).EQ.FileMissVal) RawQC(XMonth,XStation,2) = RawQC(XMonth,XStation,2) + 4

end subroutine Step2

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

subroutine Step3

RawQC(XMonth,XStation,3) = 0

if (Raw(XMonth,XStation, 5).NE.FileMissVal.AND.Raw(XMonth,XStation, 5).LT.15) &
    					RawQC(XMonth,XStation,3) = RawQC(XMonth,XStation,3) + 1
if (Raw(XMonth,XStation,10).NE.FileMissVal.AND.Raw(XMonth,XStation, 5).LT.15) &
    					RawQC(XMonth,XStation,3) = RawQC(XMonth,XStation,3) + 2
if (Raw(XMonth,XStation,10).NE.FileMissVal.AND.Raw(XMonth,XStation, 5).LT.15) &
    					RawQC(XMonth,XStation,3) = RawQC(XMonth,XStation,3) + 4

end subroutine Step3

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

subroutine Step4

RawQC(XMonth,XStation,4) = 0

if (Raw(XMonth,XStation,13).GE.Raw(XMonth,XStation,4).OR. &
    			Raw(XMonth,XStation, 9).LE.Raw(XMonth,XStation,4)) &
    		RawQC(XMonth,XStation,4) = RawQC(XMonth,XStation,4) + 1

AvMinMax = (Raw(XMonth,XStation,13)+Raw(XMonth,XStation,9)) / 2.0

if (abs(AvMinMax-Raw(XMonth,XStation,4)).GT.20.0) &
    		RawQC(XMonth,XStation,4) = RawQC(XMonth,XStation,4) + 2

end subroutine Step4

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

subroutine Step5

RawQC(XMonth,XStation,5) = 0

if (Raw(XMonth,XStation, 7).EQ.FileMissVal) RawQC(XMonth,XStation,5) = RawQC(XMonth,XStation,5) +  1
if (Raw(XMonth,XStation,14).EQ.FileMissVal) RawQC(XMonth,XStation,5) = RawQC(XMonth,XStation,5) +  2
if (Raw(XMonth,XStation,11).EQ.FileMissVal) RawQC(XMonth,XStation,5) = RawQC(XMonth,XStation,5) +  4
if (Raw(XMonth,XStation, 6).EQ.FileMissVal) RawQC(XMonth,XStation,5) = RawQC(XMonth,XStation,5) +  8

end subroutine Step5

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

subroutine Step6

RawQC(XMonth,XStation,6) = 0

if (Raw(XMonth,XStation, 8).NE.FileMissVal.AND.Raw(XMonth,XStation, 8).LT.15) &
		RawQC(XMonth,XStation,6) = RawQC(XMonth,XStation,6) + 1
if (Raw(XMonth,XStation,12).NE.FileMissVal.AND.Raw(XMonth,XStation,12).LT.15) &
		RawQC(XMonth,XStation,6) = RawQC(XMonth,XStation,6) + 2
if (Raw(XMonth,XStation,12).NE.FileMissVal.AND.Raw(XMonth,XStation,12).LT.15) &
		RawQC(XMonth,XStation,6) = RawQC(XMonth,XStation,6) + 4

if (Raw(XMonth,XStation,1).NE.FileMissVal.AND.Raw(XMonth,XStation,2).NE.FileMissVal &
	.AND.(Raw(XMonth,XStation,2)-Raw(XMonth,XStation,1)).LT.14) & 
		RawQC(XMonth,XStation,6) = RawQC(XMonth,XStation,6) + 8

end subroutine Step6

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

subroutine Step7

RawQC(XMonth,XStation,7) = 0

if (RawQC(XMonth,XStation,6).NE.1.AND.RawQC(XMonth,XStation,6).NE.3.AND.RawQC(XMonth,XStation,6).NE.5) then
 if (Raw(XMonth,XStation, 7).NE.FileMissVal) then
  if (Raw(XMonth,XStation, 4).LT.(Raw(XMonth,XStation, 7)-(Raw(XMonth,XStation,6)*KSigmas)) .OR. &
      Raw(XMonth,XStation, 4).GT.(Raw(XMonth,XStation, 7)+(Raw(XMonth,XStation,6)*KSigmas))) &
      		RawQC(XMonth,XStation,7) = RawQC(XMonth,XStation,7) +  1
 end if
end if

if (RawQC(XMonth,XStation,6).NE.2.AND.RawQC(XMonth,XStation,6).NE.3.AND.RawQC(XMonth,XStation,6).NE.6) then
 if (Raw(XMonth,XStation,14).NE.FileMissVal) then
  if (Raw(XMonth,XStation,13).LT.(Raw(XMonth,XStation,14)-(Raw(XMonth,XStation,6)*KSigmas)) .OR. &
      Raw(XMonth,XStation,13).GT.(Raw(XMonth,XStation,14)+(Raw(XMonth,XStation,6)*KSigmas))) &
      		RawQC(XMonth,XStation,7) = RawQC(XMonth,XStation,7) +  2
 end if
end if

if (RawQC(XMonth,XStation,6).NE.4.AND.RawQC(XMonth,XStation,6).NE.5.AND.RawQC(XMonth,XStation,6).NE.6) then
 if (Raw(XMonth,XStation,11).NE.FileMissVal) then
  if (Raw(XMonth,XStation, 9).LT.(Raw(XMonth,XStation,11)-(Raw(XMonth,XStation,6)*KSigmas)) .OR. &
      Raw(XMonth,XStation, 9).GT.(Raw(XMonth,XStation,11)+(Raw(XMonth,XStation,6)*KSigmas))) &
      		RawQC(XMonth,XStation,7) = RawQC(XMonth,XStation,7) +  4
 end if
end if

end subroutine Step7

!*******************************************************************************
! can we make a comparison between PhilJ and CLIMAT+ normal or sigma ?

subroutine Step8

!print*, "  > step8..."				! #################################

RawQC(XMonth,XStation,8) = 0

if (Raw(XMonth,XStation,1).NE.FileMissVal.AND.Raw(XMonth,XStation,2).NE.FileMissVal.AND. &
	(Raw(XMonth,XStation,2)-Raw(XMonth,XStation,1)).GE.14) then
 
 if (Station(XStation,3).NE.MissVal) then
  if (PhilMeans(Station(XStation,3),XMonth).EQ.MissVal.OR. &
	PhilMeansYr0(Station(XStation,3)).EQ.MissVal.OR.PhilMeansYr1(Station(XStation,3)).EQ.MissVal) then
   RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
  else
   if ((PhilMeansYr1(Station(XStation,3))-PhilMeansYr0(Station(XStation,3))).LT.14) then
    RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
   else
    if (Raw(XMonth,XStation,7).NE.FileMissVal) then
     if (Raw(XMonth,XStation,8).EQ.FileMissVal.OR.Raw(XMonth,XStation,8).GE.15) then
      if ((Raw(XMonth,XStation,1)+1900.0).NE.PhilMeansYr0(Station(XStation,3)).OR. &
      		(Raw(XMonth,XStation,2)+1900.0).NE.PhilMeansYr1(Station(XStation,3))) then
        RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
      end if
     else
      RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
     end if
    else
     RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
    end if
   end if
  end if
 else
  RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 1
 end if
 
 if (Station(XStation,4).NE.MissVal) then
  if (PhilSigma(Station(XStation,4),XMonth).EQ.MissVal.OR. &
	PhilSigmaYr0(Station(XStation,4)).EQ.MissVal.OR.PhilSigmaYr1(Station(XStation,4)).EQ.MissVal) then
   RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
  else
   if ((PhilSigmaYr1(Station(XStation,4))-PhilSigmaYr0(Station(XStation,4))).LT.14) then
    RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
   else
    if (Raw(XMonth,XStation,6).NE.FileMissVal) then
     if (Raw(XMonth,XStation,8).EQ.FileMissVal.OR.Raw(XMonth,XStation,8).GE.15) then
      if ((Raw(XMonth,XStation,1)+1900.0).NE.PhilSigmaYr0(Station(XStation,4)).OR. &
      		(Raw(XMonth,XStation,2)+1900.0).NE.PhilSigmaYr1(Station(XStation,4))) then
        RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
      end if
     else
      RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
     end if
    else
     RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
    end if
   end if
  end if
 else
  RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 2
 end if
else
 RawQC(XMonth,XStation,8) = RawQC(XMonth,XStation,8) + 3
end if

!print*, "  > end"				! #################################

end subroutine Step8

!*******************************************************************************
! make the comparison between between PhilJ and CLIMAT+ normal or sigma

subroutine Step9

!print*, "  > step9..."				! #################################

RawQC(XMonth,XStation,9) = 0

if (RawQC(XMonth,XStation,8).NE.1.AND. &
		abs(Raw(XMonth,XStation,7)-PhilMeans(Station(XStation,3),XMonth)).GT.KPJThresh) &
		RawQC(XMonth,XStation,9) = RawQC(XMonth,XStation,9) + 1

if (RawQC(XMonth,XStation,8).NE.2.AND. &
		abs(Raw(XMonth,XStation,6)-PhilSigma(Station(XStation,4),XMonth)).GT.KPJThresh) &
		RawQC(XMonth,XStation,9) = RawQC(XMonth,XStation,9) + 2

!print*, "  > end"				! #################################

end subroutine Step9

!*******************************************************************************
! can we use PhilJ sigma to compare CLIMAT+ values with normals ?

subroutine Step10

!print*, "  > step10..."				! #################################
!print*, XMonth,XStation
!print*, Station(XStation,4),PhilSigma(Station(XStation,4),XMonth)
!print*, PhilSigmaYr1(Station(XStation,4)),PhilSigmaYr0(Station(XStation,4))

RawQC(XMonth,XStation,10) = 0

if (Station(XStation,4).EQ.MissVal) then
 RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 7
else
 if (PhilSigma(Station(XStation,4),XMonth).EQ.MissVal.OR. &
			(PhilSigmaYr1(Station(XStation,4))-PhilSigmaYr0(Station(XStation,4))).LT.14) then
  RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 7
 else
  if (Station(XStation,3).EQ.MissVal) then
   RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 1
  else
   if (PhilMeans(Station(XStation,3),XMonth).EQ.MissVal.OR. &
			(PhilMeansYr1(Station(XStation,3))-PhilMeansYr0(Station(XStation,3))).LT.14) then
    RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 1
   end if
  end if
  
  if (Raw(XMonth,XStation,1).EQ.FileMissVal.OR.Raw(XMonth,XStation,2).EQ.FileMissVal.OR. &
	(Raw(XMonth,XStation,2)-Raw(XMonth,XStation,1)).GE.14) then
    if (Raw(XMonth,XStation,14).NE.FileMissVal) then
      if (Raw(XMonth,XStation,12).NE.FileMissVal.AND.Raw(XMonth,XStation,12).LT.15) &
   		RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 2
    else
      RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 2
    end if
    
    if (Raw(XMonth,XStation,11).NE.FileMissVal) then
      if (Raw(XMonth,XStation,12).NE.FileMissVal.AND.Raw(XMonth,XStation,12).LT.15) &
   		RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 4
    else
      RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 4
    end if    
  else
    RawQC(XMonth,XStation,10) = RawQC(XMonth,XStation,10) + 6
  end if
 end if
end if

!print*, "  > end"				! #################################

end subroutine Step10

!*******************************************************************************
! we use PhilJ sigma to compare CLIMAT+ values with normals

subroutine Step11

!print*, "  > step11..."				! #################################

RawQC(XMonth,XStation,11) = 0

if (RawQC(XMonth,XStation,10).NE.1.AND.RawQC(XMonth,XStation,10).NE.3.AND.RawQC(XMonth,XStation,10).NE.5) then
 if (Raw(XMonth,XStation, 4).LT.(PhilMeans(Station(XStation,3),XMonth)-(PhilSigma(Station(XStation,4),XMonth)*KSigmas)) .OR. &
     Raw(XMonth,XStation, 4).GT.(PhilMeans(Station(XStation,3),XMonth)+(PhilSigma(Station(XStation,4),XMonth)*KSigmas))) &
      		RawQC(XMonth,XStation,11) = RawQC(XMonth,XStation,11) +  1
end if

if (RawQC(XMonth,XStation,10).NE.2.AND.RawQC(XMonth,XStation,10).NE.3.AND.RawQC(XMonth,XStation,10).NE.6) then
 if (Raw(XMonth,XStation,13).LT.(Raw(XMonth,XStation,14)-(PhilSigma(Station(XStation,4),XMonth)*KSigmas)) .OR. &
     Raw(XMonth,XStation,13).GT.(Raw(XMonth,XStation,14)+(PhilSigma(Station(XStation,4),XMonth)*KSigmas))) &
      		RawQC(XMonth,XStation,11) = RawQC(XMonth,XStation,11) +  2
end if

if (RawQC(XMonth,XStation,10).NE.4.AND.RawQC(XMonth,XStation,10).NE.5.AND.RawQC(XMonth,XStation,10).NE.6) then
 if (Raw(XMonth,XStation, 9).LT.(Raw(XMonth,XStation,11)-(PhilSigma(Station(XStation,4),XMonth)*KSigmas)) .OR. &
     Raw(XMonth,XStation, 9).GT.(Raw(XMonth,XStation,11)+(PhilSigma(Station(XStation,4),XMonth)*KSigmas))) &
      		RawQC(XMonth,XStation,11) = RawQC(XMonth,XStation,11) +  4
end if

!print*, "  > end"				! #################################

end subroutine Step11

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

subroutine DumpRawInfo

write (99,*), "Summary table of raw data in file, by year/month: "
write (99,"(a)"), " Year  Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec  ???"
do XRawFile = 1, RawFileN
  write (99,"(i5,13i5)"), XRawFile, (PreQCTime(XRawFile,XMonth),XMonth=1,13)
end do
write (99,*)

write (99,*), "Summary table of raw data in file, by pre-QC category: "
write (99,"(a)"), " Year Lines NoMon n*Mon n*Inf Valid"
do XRawFile = 1, RawFileN
  write (99,"(i5,5i6)"), XRawFile, (CheckRaw(XRawFile,XInfo),XInfo=1,5)
end do
write (99,*)

write (99,*), "Summary table of valid data in file, by year/month: "
write (99,"(a)"), " Year  Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec"
do XRawFile = 1, RawFileN
  write (99,"(i5,12i5)"), XRawFile, (RawTime(XRawFile,XMonth),XMonth=1,12)
end do
write (99,*)

allocate (RegInfo (RawFileN,RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpRawInfo: Allocation failure: RegInfo #####"
RegInfo = 0

write (99,*), "Summary table of raw loaded data: "
write (99,*), "Year Euro USSR MidE Asia Afri NAme CAme SAme Anta Ocea Buoy Misc   TOT"
do XRawFile = 1, RawFileN
  do XCty = 1, (CtyN+1)
    do XMonth = 1, MonthN
      RegInfo(XRawFile,CtyReg(XCty)) = RegInfo(XRawFile,CtyReg(XCty)) + CtyRaw(XRawFile,XMonth,XCty)
    end do
  end do
  
  OpTot = 0
  do XReg = 1, RegN
    OpTot = OpTot + RegInfo(XRawFile,XReg)
  end do
  
  write (99,"(i5,12i5,i6)"), XRawFile, (RegInfo(XRawFile,XReg),XReg=1,RegN), OpTot
end do
write (99,*)

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

end subroutine DumpRawInfo

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

subroutine DumpQCInfo

allocate (SumQC (RegN,StepN+1), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpQCInfo: Allocation failure: SumQC 1 #####"
SumQC = 0

write (99,*), "Summary table of all data rejected."
write (99,"(a)"), " Step  Euro  USSR  MidE  Asia  Afri  NAme  CAme  SAme  Anta  Ocea  Buoy  Misc   TOT"
do XStep = 1, StepN
  OpTot = 0
  do XReg = 1, RegN
    do XRawFile = 1, RawFileN
      SumQC(XReg,XStep) = SumQC(XReg,XStep) + RegQC(XRawFile,XReg,XStep)
    end do
    OpTot = OpTot + SumQC(XReg,XStep)
  end do
  write (99,"(i5,13i6)"), XStep, (SumQC(XReg,XStep),XReg=1,RegN), OpTot
end do

OpTot = 0
do XReg = 1, RegN
    do XStep = 1, StepN
      SumQC(XReg,StepN+1) = SumQC(XReg,StepN+1) + SumQC(XReg,XStep)
    end do
    OpTot = OpTot + SumQC(XReg,StepN+1)
end do
write (99,"(a5,13i6)"), "TOTAL", (SumQC(XReg,StepN+1),XReg=1,RegN), OpTot
write (99,*)

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

allocate (SumQC (StepN,BitN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpQCInfo: Allocation failure: SumQC 2 #####"
SumQC = 0

write (99,*), "Summary table of all data processing."
write (99,"(a)"), "       0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15   SUM"
do XStep = 1, StepN
  OpTot = 0
  do XBit = 0, BitN-1
    do XRawFile = 1, RawFileN
      SumQC(XStep,XBit) = SumQC(XStep,XBit) + BitQC(XRawFile,XStep,XBit)
    end do
    OpTot = OpTot + SumQC(XStep,XBit)
  end do
  write (99,"(i2,17i6)"), XStep, (SumQC(XStep,XBit),XBit=0,BitN-1), OpTot
end do
write (99,*)

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

write (99,*), "Summary table of all stations, classed by % of months rejected."
write (99,"(a)"), " Year     0  1-10 11-30 31-50 51-80 81-99   100   TOT"
do XRawFile = 1, RawFileN
  do XClass = 1, ClassN
    StnQC(XRawFile,ClassN+1) = StnQC(XRawFile,ClassN+1) + StnQC(XRawFile,XClass)
  end do
  
  write (99,"(i5,8i6)"), XRawFile, (StnQC(XRawFile,XClass),XClass=1,ClassN+1)
end do

do XClass = 1, ClassN+1
  do XRawFile = 1, RawFileN
    StnQC(RawFileN+1,XClass) = StnQC(RawFileN+1,XClass) + StnQC(XRawFile,XClass)
  end do
end do
write (99,"(a5,8i6)"), "  TOT", (StnQC(RawFileN+1,XClass),XClass=1,ClassN+1)
write (99,*)

do XStep = 1, StepN
 if (trim(StepText(XStep)).NE."no rejections") then
  write (99,*), "Summary table of data rejected at step: ", XStep
  write (99,*), trim(StepText(XStep))
  write (99,*), "Year Euro USSR MidE Asia Afri NAme CAme SAme Anta Ocea Buoy Misc   TOT"
 
  do XRawFile = 1, RawFileN
   OpTot = 0
   do XReg = 1, RegN
    OpTot = OpTot + RegQC(XRawFile,XReg,XStep)
   end do
  
   write (99,"(i5,12i5,i6)"), XRawFile, (RegQC(XRawFile,XReg,XStep),XReg=1,RegN), OpTot
  end do
  write (99,*)
 end if
end do

write (99,*), "Key to bits in data processing tables: "
do XStep = 1, StepN
  write (99,"(i2,a1,a)"), XStep, " ", trim(BitText(XStep))
end do
write (99,*)
 
do XRawFile = 1, RawFileN
  write (99,*), "Summary table of data processing for year: ", XRawFile
  write (99,"(a)"), "      0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15  SUM"
    
  do XStep = 1, StepN
    OpTot = 0
    do XBit = 0, BitN-1
      OpTot = OpTot + BitQC(XRawFile,XStep,XBit)
    end do
    
    write (99,"(i2,17i5)"), XStep, (BitQC(XRawFile,XStep,XBit),XBit=0,BitN-1), OpTot
  end do
  
  write (99,*)
end do

end subroutine DumpQCInfo

!*******************************************************************************
! this is the subroutine that will finally dump the QC'd data into CRU ts files
! when I have got it working, temporarily reverse the process above in which DTR is filled
!	so that it fills with rejected data instead, so that I can eyeball the rejected data

subroutine DumpData

print*, "  > Saving the data as a CRU ts file..."

allocate (Code(StationN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpData: Allocation failure #####"
do XStation = 1, StationN
  Code(XStation)=Station(XStation,1)
end do

call GetCRUtsInfo (Code,Lat,Long,Elev,NameStn,NameCty,6)

call MakeStnInfo (StnInfo,Code,Lat,Long,Elev,1,YearAD=YearAD,Data=DTR,DataMissVal=DataMissVal)

call SaveCTS (StnInfo,NameStn,NameCty,Data=DTR,YearAD=YearAD,Silent=1)

deallocate (Code,Lat,Long,Elev,NameStn,NameCty,StnInfo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpData: Deallocation failure #####"

end subroutine DumpData

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

subroutine LocateStn

XCty = 0 ; ThisCty = 0
do
  XCty = XCty + 1
  
  if (XCheckStn.GE.CtyBound0(XCty).AND.XCheckStn.LE.CtyBound1(XCty)) ThisCty = XCty
  if (XCty.EQ.CtyN.AND.ThisCty.EQ.0) ThisCty = CtyN + 1
  
  if (XCty.EQ.CtyN) exit
end do

end subroutine LocateStn

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

subroutine Conclude

print*

close (99)

end subroutine Conclude

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

end program QCClimatT