! gethahn.f90
! program written by Tim Mitchell, Tyndall Centre
! loads the (v large) Hahn data-set and extracts the useful stuff and saves it to CRU ts files
! now does this in a two-pass process, as per Phil Jones email on 22.02.02
!   pass 1: create a clim for each month/hour/stn
!   pass 2: calc anoms, av 2+ anom --> day, av 20+ days --> month
! now calcs wind direction properly
! written 13.02.02
! last modified 25.04.02
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! f90 -o ./../obs/gethahn time.f90 filenames.f90 crutsfiles.f90 wmokey.f90 ./../obs/gethahn.f90

program GetHahn

use Time
use FileNames
use CRUtsFiles
use WMOKey

implicit none

integer, parameter :: QPass = 1		!@@@@@@@@@@@@@@@@@@ set to either 1 or 2 @@@@@@@@@@@@@@@@@

integer, parameter :: ThreshClim = 150	! min values to set a climatology (rec=?)
integer, parameter :: ThreshDay = 20	! min days to calc a month anom (rec=20)
integer, parameter :: ThreshHour = 2	! min hours to calc a day anom (rec=2)

integer, parameter :: BegYear = 1971	! 1971...1996
integer, parameter :: BegMonth = 1	! 1...12
integer, parameter :: EndYear = 1996	! 1971...1996
integer, parameter :: EndMonth = 12	! 1...12

real, pointer, dimension (:,:,:)	:: OutOct,OutSLP,OutWnd,OutDir,OutVap,OutTmp,OutWdx,OutWdy,OutAny
real, pointer, dimension (:,:,:)	:: ClimOct,ClimSLP,ClimWnd,ClimDir,ClimVap,ClimTmp,ClimWdx,ClimWdy,ClimAny
real, pointer, dimension (:)		:: Lat,Lon,Elv
real, pointer, dimension (:)		:: MyLat,MyLon,MyElv
real, dimension (8)			:: HourMean
real, dimension (8)			:: InMulti  = (/1.00,0.1,0.1,0.1,0.1,0.001,0.001,1.0/)
real, dimension (8)			:: OutMulti = (/0.01,0.1,0.1,0.1,0.1,0.001,0.001,0.1/)

integer, pointer, dimension (:,:,:)		:: InOct,InSLP,InWnd,InDir,InVap,InTmp,InWdx,InWdy,InAny
integer, pointer, dimension (:,:,:)		:: TotOct,TotSLP,TotWnd,TotDir,TotVap,TotTmp,TotWdx,TotWdy,TotAny
integer, pointer, dimension (:,:,:)		:: SaveAny,CountStnOut,CountStnTot
integer, pointer, dimension (:,:)		:: StnInfo,CountStnIn,CountAllIn,MonthLengths
integer, pointer, dimension (:)			:: Code			! WMO codes (6-digit, size=StnN)
integer, pointer, dimension (:)			:: YearAD

integer, dimension (31)				:: DayObs
integer, dimension (8)				:: HourObs
integer, dimension (8)				:: InMiss = (/-1,-1,-1,900,900,-999,-999,-1/)

character (len=3), dimension (12), parameter 	:: MonthName=(/'JAN','FEB','MAR','APR','MAY','JUN', &
							       'JUL','AUG','SEP','OCT','NOV','DEC'/)
character (len=80), pointer, dimension (:,:) 	:: RawFiles
character (len=20), pointer, dimension (:)	:: MyNameStn,NameStn
character (len=13), pointer, dimension (:)	:: MyNameCty,NameCty
character (len=09), pointer, dimension (:)	:: MyNameLoc,NameLoc
character (len=04), dimension (8)		:: OutSuffices = &
				(/'.oct','.slp','.wnd','.vap','.tmp','.wdx','.wdy','.dir'/)

real, parameter :: Kay1wet =  17.380
real, parameter :: Kay1ice =  22.440
real, parameter :: Kay2wet = 239.000
real, parameter :: Kay2ice = 272.400
real, parameter :: Kay3    =   6.107
real, parameter :: Pi      =   3.1415927

real :: OpTot,OpEn, HourTot,HourEn, DayTot,DayEn, MonthClim, Angle
real :: Tdry,Tdew,Twet

integer :: ReadStatus, AllocStat
integer :: NStn,NYear,NMonth,NDay,NHour,NLine,NCode,NHourValid,NDayValid,NVari
integer :: XStn,XYear,XMonth,XDay,XHour,XLine,XCode,XHourValid,XDayValid,XVari
integer :: PrevDay,PrevHour,PrevCode
integer :: QZip
integer :: LoadFileLen,HourOpp,WMOCode
integer :: IDy,IHr,IIB,ILat,ILon,IID,ICld,ISLP,IWS,IWD,IAT,IDD,IEl
integer :: OpMin,OpMax,Hour0,Hour1,HourGap,TrashInt

character (len=80) :: RawFormat = "(4x,2i2,i1,3i5,3x,i1,28x,i5,2i3,i4,i3,i4,2x)"
character (len=80) :: CodeFile,LoadFile,SaveStem,SaveTip,SaveFile
character (len=20) :: YearStr20,BegYearTxt,EndYearTxt
character (len=02) :: YearStr02

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

call Intro
call GetCode
call GetFilePaths
call AllocArrays
if (QPass.EQ.2) call LoadSynClim

do XYear = 1, NYear
  do XMonth = 1, NMonth
    if (RawFiles(XYear,XMonth).NE."blank") then
      print "(a,a3,x,i4)", "   > Operating on ", MonthName(XMonth), YearAD(XYear)
      
      call PrepRawFile
      call StoreIn
      call ShutRawFile
      if (QPass.EQ.2) call InToOut
      if (QPass.EQ.1) call InToClim
      
      call DumpYrMonStats
    end if
  end do
end do

call HarmoniseStn

if (QPass.EQ.2) then
  call RestoreAbs
  call MakeDir
  call DumpCruTS
  call DumpOpStats
else if (QPass.EQ.1) then
  call CheckClim
  call DumpSynClim
  call DumpClimStats
end if

call Conclude

contains

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

subroutine Intro

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

print*
print*, "  > ***** GetHahn.f90: loads Hahn data files and saves to CRU ts *****"
print*

end subroutine Intro

!*******************************************************************************
! these codes are 5-digit WMO codes

subroutine GetCode

CodeFile = "/tyn1/tim/data/obs/hahn/ystaty.txt"

open  (1,file=CodeFile,status="old",action="read")
read  (1,"(i6)"), NStn

allocate (Code(NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetCode: Allocation failure #####"

do XStn = 1, NStn
  read  (1,"(i5)"), Code(XStn)
end do
close (1)

end subroutine GetCode

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

subroutine GetFilePaths

NYear = EndYear - BegYear + 1
NMonth = 12 ; NDay = 31 ; NHour = 8 ; NVari = 8

allocate (RawFiles(NYear,NMonth), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetFilePaths: Allocation failure: RawFiles* #####"
RawFiles = "blank"

XYear = 1 ; XMonth = BegMonth
do
  YearStr20 = GetTextFromInt (BegYear+XYear-1)
  YearStr02 = trim(adjustl(YearStr20(3:4)))
    
  do
    RawFiles(XYear,XMonth) = "/tyn1/tim/data/obs/hahn/_raw/" // MonthName(XMonth) // YearStr02 // "L"
    
    XMonth = XMonth + 1
    if (XYear.EQ.NYear.AND.XMonth.GT.EndMonth) exit
    if (XMonth.GT.NMonth) exit
  end do
  
  XMonth = 1
  XYear  = XYear + 1
  if (XYear.GT.NYear) exit
end do

end subroutine GetFilePaths

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

subroutine AllocArrays

allocate (InOct(NDay,NHour,NStn), &		! for initi see PrepRawFile
	  InSLP(NDay,NHour,NStn), &
	  InWdx(NDay,NHour,NStn), &
	  InWdy(NDay,NHour,NStn), &
	  InWnd(NDay,NHour,NStn), &
	  InVap(NDay,NHour,NStn), &
	  InTmp(NDay,NHour,NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: In* #####"

allocate (ClimOct(NMonth,NHour,NStn), &
	  ClimSLP(NMonth,NHour,NStn), &
	  ClimWdx(NMonth,NHour,NStn), &
	  ClimWdy(NMonth,NHour,NStn), &
	  ClimWnd(NMonth,NHour,NStn), &
	  ClimVap(NMonth,NHour,NStn), &
	  ClimTmp(NMonth,NHour,NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: Clim* #####"
ClimOct=-9999 ; ClimSLP=-9999 ; ClimWnd=-9999 ; ClimVap=-9999 ; ClimTmp=-9999
ClimWdx=-9999 ; ClimWdy=-9999

if (QPass.EQ.2) then
  allocate (OutOct(NYear,NMonth,NStn), &
	    OutSLP(NYear,NMonth,NStn), &
	    OutWdx(NYear,NMonth,NStn), &
	    OutWdy(NYear,NMonth,NStn), &
	    OutWnd(NYear,NMonth,NStn), &
	    OutDir(NYear,NMonth,NStn), &
	    OutVap(NYear,NMonth,NStn), &
	    OutTmp(NYear,NMonth,NStn), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: Out* #####"
  OutOct=-9999 ; OutSLP=-9999 ; OutWnd=-9999 ; OutDir=-9999 ; OutVap=-9999 ; OutTmp=-9999
  OutWdx=-9999 ; OutWdy=-9999
else if (QPass.EQ.1) then
  allocate (TotOct(NMonth,NHour,NStn), &
	    TotSLP(NMonth,NHour,NStn), &
	    TotWdx(NMonth,NHour,NStn), &
	    TotWdy(NMonth,NHour,NStn), &
	    TotWnd(NMonth,NHour,NStn), &
	    TotVap(NMonth,NHour,NStn), &
	    TotTmp(NMonth,NHour,NStn), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: Tot* #####"
  TotOct=0 ; TotSLP=0 ; TotWdx=0 ; TotWdy=0 ; TotWnd=0 ; TotVap=0 ; TotTmp=0
end if

allocate (Lat(NStn), &
	  Lon(NStn), &
	  Elv(NStn), &
	  NameStn(NStn), &
	  NameCty(NStn), &
	  NameLoc(NStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: Stn* #####"
Lat=-999.0 ; Lon=-999.0 ; Elv=-999.0 ; NameStn="" ; NameCty="unidentified" ; NameLoc="   nocode"

allocate (CountAllIn (NYear,NMonth),       &
	  CountStnIn (NDay,NHour),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: Count* #####"
CountAllIn = 0 ; CountStnIn = 0

if (QPass.EQ.2) then
  allocate (CountStnOut(NYear,NMonth,NVari), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: CountStnOut #####"
  CountStnOut = 0
else if (QPass.EQ.1) then
  allocate (CountStnTot(NMonth,NHour,NVari), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: CountStnTot #####"
  CountStnTot = 0
end if

allocate (YearAD(NYear),MonthLengths(NYear,NMonth), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AllocArrays: Allocation failure: YearAD #####"
do XYear = 1, NYear
  YearAD(XYear) = BegYear + XYear - 1
end do

call GetMonthLengths (YearAD,MonthLengths)

end subroutine AllocArrays

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

subroutine PrepRawFile

InOct=-999 ; InSLP=-999 ; InWdx=-999 ; InWdy=-999 ; InWnd=-999 ; InVap=-999 ; InTmp=-999

LoadFile = LoadPath(RawFiles(XYear,XMonth),"    ")
LoadFileLen = len_trim(LoadFile)
if (LoadFileLen.GT.1.AND.LoadFile((LoadFileLen-1):LoadFileLen).EQ.".Z") then
  QZip = 2							! file is zipped
  call system ('uncompress ' // LoadFile)
  LoadFile ((LoadFileLen-1):LoadFileLen) = "  "			! change filename to the unzipped file
else
  QZip = 1							! file already unzipped
end if

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

end subroutine PrepRawFile

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

subroutine StoreIn

CountStnIn = 0 ; PrevCode = 0 ; PrevDay = 0 ; PrevHour = 0

open  (1,file=LoadFile,status="old",action="read")
do XLine = 1, NLine
  read (1,RawFormat), IDy,IHr,IIB,ILat,ILon,IID,ICld,ISLP,IWS,IWD,IAT,IDD,IEl
  
  XDay = IDy ; XHour = (IHr/3)+1 ; WMOCode = IID ; XStn = -999  ;  XCode = PrevCode
  
  if (PrevDay.NE.XDay.OR.PrevHour.NE.XHour.OR.WMOCode.LT.Code(PrevCode)) XCode = 0
  do								! find station
    XCode = XCode + 1
    
    if (WMOCode.EQ.Code(XCode)) then
    	PrevCode = XCode
    	XStn = XCode
    end if
    
    if (WMOCode.LT.Code(XCode)) exit    
    if (XCode.EQ.NStn) 		exit
    if (XStn.NE.-999) 		exit
  end do
  PrevDay = XDay ; PrevHour = XHour
  
  if (XStn.NE.-999) then
    CountStnIn(XDay,XHour) = CountStnIn(XDay,XHour) + 1
    
    if (IDD.NE.700.AND.IDD.NE.900.AND.IAT.NE.900) then
      Tdry = real(IAT) * 0.1
      Tdew = Tdry - (real(IDD) * 0.1)
      Twet = (0.6 * Tdew) + (0.4 * Tdry)

      if (Twet.GE.0) then
        InVap(XDay,XHour,XStn) = nint(10 * Kay3 * exp((Kay1wet * Tdew) / (Kay2wet + Tdew)))
      else
        InVap(XDay,XHour,XStn) = nint(10 * Kay3 * exp((Kay1ice * Tdew) / (Kay2ice + Tdew)))
      end if
    end if
    
    InSLP(XDay,XHour,XStn)=ISLP
    InTmp(XDay,XHour,XStn)=IAT 
    
    if (IIB.EQ.  1) InOct(XDay,XHour,XStn)=ICld
    if (IWS.NE.999) InWnd(XDay,XHour,XStn)=IWS 
      
    if (IWD.GE.0.AND.IWD.LE.360) then
      InWdx(XDay,XHour,XStn) = nint(1000 * sin(real(IWD)*Pi/180.0))
      InWdy(XDay,XHour,XStn) = nint(1000 * cos(real(IWD)*Pi/180.0))
    end if
    
    if (Elv(XStn).EQ.-999) then
      Lat(XStn)=real(ILat)/100.0 ; Lon(XStn)=real(ILon)/100.0 ; Elv(XStn)=real(IEl)
      if (Lon(XStn).GT.180.0) Lon(XStn) = Lon(XStn) - 360.0
    end if
  end if
end do
close (1)

end subroutine StoreIn

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

subroutine ShutRawFile

if (QZip.EQ.2) call system ('compress ' // LoadFile // ' &')
  
end subroutine ShutRawFile

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

subroutine InToClim

do XVari = 1, 7					! iterate by variable
 call SelectArraysIn ; call SelectArraysClim ; call SelectArraysTot

 do XStn = 1, NStn				! iterate by station  
   do XHour = 1, NHour				! iterate by hour
     HourTot = 0.0 ; HourEn = 0.0
     
     do XDay = 1, NDay								! check for datum
       if (InAny(XDay,XHour,XStn).NE.InMiss(XVari).AND.InAny(XDay,XHour,XStn).NE.-999) then
         HourEn = HourEn + 1.0
         HourTot = HourTot + (InAny(XDay,XHour,XStn) * InMulti(XVari))
       end if
     end do
     
     if (HourEn.GT.0) then							! if valid data
      if (TotAny(XMonth,XHour,XStn).EQ.0) then					! if no previous data
       ClimAny(XMonth,XHour,XStn) = HourTot/HourEn
       TotAny(XMonth,XHour,XStn) = nint(HourEn)
      else									! else update average
       ClimAny(XMonth,XHour,XStn) = ((ClimAny(XMonth,XHour,XStn)*real(TotAny(XMonth,XHour,XStn))) &
       					+HourTot)/(real(TotAny(XMonth,XHour,XStn))+HourEn)
       TotAny(XMonth,XHour,XStn) = TotAny(XMonth,XHour,XStn) + nint(HourEn)
      end if
     end if
   end do
 end do

 nullify(InAny,TotAny,ClimAny)
end do

end subroutine InToClim

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

subroutine InToOut

do XVari = 1, 7					! iterate by variable
 call SelectArraysIn ; call SelectArraysClim ; call SelectArraysOut
 
 do XStn = 1, NStn				! iterate by station  
    DayTot = 0.0 ; DayEn = 0.0
    
    do XDay = 1, NDay				! iterate by day
      HourTot = 0.0 ; HourEn = 0.0
      
      do XHour = 1, NHour			! iterate by hour
        if (ClimAny(XMonth,XHour,XStn).NE.-9999) then				! check for clim
          if (InAny(XDay,XHour,XStn).NE.InMiss(XVari).AND.InAny(XDay,XHour,XStn).NE.-999) then
            HourEn = HourEn + 1.0
            HourTot = HourTot + ((InAny(XDay,XHour,XStn)*InMulti(XVari)) - ClimAny(XMonth,XHour,XStn))
          end if
        end if
      end do
      
      if (HourEn.GE.ThreshHour) then
        DayEn = DayEn + 1.0
        DayTot = DayTot + (HourTot / HourEn)
      end if
    end do
    
    if (DayEn.GE.ThreshDay) OutAny(XYear,XMonth,XStn) = DayTot / DayEn
 end do
 
 nullify(InAny,OutAny,ClimAny)
end do

end subroutine InToOut

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

subroutine HarmoniseStn

call GetCRUtsInfo (Code,MyLat,MyLon,MyElv,MyNameStn,MyNameCty,Digits=5) 

do XStn = 1, NStn
  if (MyLat(XStn).NE.-999.0.AND.Lat(XStn).NE.-999.0.AND.abs(MyLat(XStn)-Lat(XStn)).LE.0.1) then
    if (MyLon(XStn).NE.-999.0.AND.Lon(XStn).NE.-999.0.AND.abs(MyLon(XStn)-Lon(XStn)).LE.0.1) then
      NameStn(XStn) = MyNameStn(XStn)		! NameLoc is already set to "   nocode"
    end if
  end if
  
  NameCty(XStn) = MyNameCty(XStn)
end do

deallocate (MyLat,MyLon,MyElv,MyNameStn,MyNameCty,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: HarmoniseStn: Deallocation failure #####"

Code = Code * 10

end subroutine HarmoniseStn

!*******************************************************************************
! if pass 1, check values going into clim total at least 30

subroutine CheckClim

do XVari = 1, 7
 call SelectArraysClim ; call SelectArraysTot
 
 do XMonth = 1, NMonth
  do XStn = 1, NStn				! iterate by station  
   do XHour = 1, NHour				! iterate by hour
     if (TotAny(XMonth,XHour,XStn).LT.ThreshClim) then
     	ClimAny(XMonth,XHour,XStn) = -9999					! sets missing
     else
        CountStnTot(XMonth,XHour,XVari)=CountStnTot(XMonth,XHour,XVari)+1	! counts valid
     end if
   end do
  end do
 end do
 
 nullify (TotAny,ClimAny)
end do

end subroutine CheckClim

!*******************************************************************************
! if pass 2, restore anomalies back to absolute values

subroutine RestoreAbs

do XVari = 1, 7
 call SelectArraysClim ; call SelectArraysOut
 
 do XStn = 1, NStn					! iterate by station  
   do XMonth = 1, NMonth				! iterate by month
     MonthClim=-9999 ; HourObs=0 ; Hour0=-9999 ; HourGap=-9999 ; Hour1=-9999 ; HourTot=0
     do XHour = 1, NHour				! find which hours have clims
       if (ClimAny(XMonth,XHour,XStn).NE.-9999) then
       	 HourObs(XHour) = 1
       	 HourTot = HourTot + 1
       end if
     end do
     
     if (HourTot.EQ.8) then
        				Hour0=1 ; HourGap=1 ; Hour1=8
     else if (HourObs(1).EQ.1.AND.HourObs(3).EQ.1.AND.HourObs(5).EQ.1.AND.HourObs(7).EQ.1) then
        				Hour0=1 ; HourGap=2 ; Hour1=7
     else if (HourObs(2).EQ.1.AND.HourObs(4).EQ.1.AND.HourObs(6).EQ.1.AND.HourObs(8).EQ.1) then
        				Hour0=2 ; HourGap=2 ; Hour1=8
     else if (HourObs(1).EQ.1.AND.HourObs(5).EQ.1) then
        				Hour0=1 ; HourGap=4 ; Hour1=5
     else if (HourObs(2).EQ.1.AND.HourObs(6).EQ.1) then
        				Hour0=2 ; HourGap=4 ; Hour1=6
     else if (HourObs(3).EQ.1.AND.HourObs(7).EQ.1) then
        				Hour0=3 ; HourGap=4 ; Hour1=7
     else if (HourObs(4).EQ.1.AND.HourObs(8).EQ.1) then
        				Hour0=4 ; HourGap=4 ; Hour1=8
     end if
     
     if (Hour0.NE.-9999) then
       HourTot=0.0 ; HourEn=0.0				! calc monthly clim
       do XHour = Hour0,Hour1,HourGap
       	 HourTot = HourTot + ClimAny(XMonth,XHour,XStn)
         HourEn = HourEn + 1.0
       end do
       
       MonthClim = HourTot / HourEn
       
       do XYear = 1, NYear				! add clim back to anomalies
         if (OutAny(XYear,XMonth,XStn).NE.-9999) then
         	OutAny(XYear,XMonth,XStn) = OutAny(XYear,XMonth,XStn) + MonthClim
         	CountStnOut(XYear,XMonth,XVari)=CountStnOut(XYear,XMonth,XVari)+1 
         end if
       end do	
     else
       do XYear = 1, NYear				! reject all data for this var/stn/cal-month
         OutAny(XYear,XMonth,XStn) = -9999
       end do
     end if     
   end do
 end do
 
 nullify (OutAny,ClimAny)
end do

end subroutine RestoreAbs

!*******************************************************************************
! if pass 2, make wind direction from x,y components

subroutine MakeDir

do XStn = 1, NStn					! iterate by station  
  do XYear = 1, NYear					! iterate by year
    do XMonth = 1, NMonth				! iterate by month
      if (OutWdx(XYear,XMonth,XStn).NE.-9999.AND.OutWdy(XYear,XMonth,XStn).NE.-9999) then
        Angle = atan(abs(OutWdx(XYear,XMonth,XStn)/OutWdy(XYear,XMonth,XStn)))
        Angle = Angle * 180.0 / Pi
        
        if      (OutWdx(XYear,XMonth,XStn).GE.0.AND.OutWdy(XYear,XMonth,XStn).GE.0) then
          OutDir(XYear,XMonth,XStn) = Angle
        else if (OutWdx(XYear,XMonth,XStn).GE.0) then
          OutDir(XYear,XMonth,XStn) = 180.0 - Angle
        else if (OutWdy(XYear,XMonth,XStn).LT.0) then
          OutDir(XYear,XMonth,XStn) = 180.0 + Angle
        else
          OutDir(XYear,XMonth,XStn) = 360.0 - Angle
        end if
        
        CountStnOut(XYear,XMonth,8)=CountStnOut(XYear,XMonth,8)+1 
      end if
    end do
  end do
end do

end subroutine MakeDir

!*******************************************************************************
! I know all the variales are named Save, but it reduces variable declarations!

subroutine LoadSynClim

SaveStem = "/tyn1/tim/data/obs/hahn/" 
do XVari = 1, 7
  SaveTip  = "clim" // OutSuffices(XVari) // ".txt"
  SaveFile = trim(SaveStem) // trim(SaveTip)  
  print "(2a)", "   > Loading from: ", trim(SaveTip)
  
  call SelectArraysClim
  
  open  (2,file=SaveFile,status="old",access="sequential",form="formatted",action="read")
  do XStn = 1, NStn
    read (2,"(i8)"), TrashInt
    do XMonth = 1, NMonth
      read (2,"(8e13.5)"), (ClimAny(XMonth,XHour,XStn),XHour=1,NHour)
    end do
  end do
  close (2)
  
  nullify (ClimAny)
end do

end subroutine LoadSynClim

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

subroutine DumpSynClim

SaveStem = "/tyn1/tim/data/hahn/" 
do XVari = 1, 7
  SaveTip  = "clim" // OutSuffices(XVari) // ".txt"
  SaveFile = trim(SaveStem) // trim(SaveTip)  
  print "(2a)", "   > Saving to: ", trim(SaveTip)
  
  call SelectArraysClim
  
  open  (2,file=SaveFile,status="replace",access="sequential",form="formatted",action="write")
  do XStn = 1, NStn
    write (2,"(i8)"), Code(XStn)
    do XMonth = 1, NMonth
      write (2,"(8e13.5)"), (ClimAny(XMonth,XHour,XStn),XHour=1,NHour)
    end do
  end do
  close (2)
  
  nullify (ClimAny)
end do

end subroutine DumpSynClim

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

subroutine DumpCRUts

allocate (SaveAny(NYear,NMonth,NStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpCRUts: Allocation failure: SaveAny #####"

BegYearTxt = GetTextFromInt(BegYear)
EndYearTxt = GetTextFromInt(EndYear)

SaveStem = "/tyn1/tim/data/hahn/" 
do XVari = 1, 8
 if (XVari.NE.6.AND.XVari.NE.7) then
  SaveTip  = "hahn." // trim(BegYearTxt) // "-" // trim(EndYearTxt) // OutSuffices(XVari) // ".cts"
  SaveFile = trim(SaveStem) // trim(SaveTip)
  
  print "(2a)", "   > Saving to: ", trim(SaveTip)
  
  call SelectArraysOut
  
  SaveAny = -9999
  do XYear = 1, NYear
    do XMonth = 1, NMonth
      do XStn = 1, NStn
       if (OutAny(XYear,XMonth,XStn).NE.-9999) &
       			SaveAny(XYear,XMonth,XStn)=nint(OutAny(XYear,XMonth,XStn)/OutMulti(XVari))
      end do
    end do
  end do
  
  call MakeStnInfo (StnInfo,Code,Code,Lat,Lon,Elv,1,YearAD=YearAD,Data=SaveAny,Silent=1)
  
  call SaveCTS (StnInfo,NameLoc,NameStn,NameCty,SaveAny,YearAD,CallFile=SaveFile,Silent=1)
  
  deallocate (StnInfo, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpCRUts: Deallocation failure: StnInfo #####"

  nullify(OutAny)
 end if
end do

end subroutine DumpCRUts

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

subroutine DumpYrMonStats

write (99,"(i4,a4,a)"), YearAD(XYear),MonthName(XMonth),": station totals"
write (99,"(a)"), " DAY    00    03    06    09    12    15    18    21"
do XDay = 1, MonthLengths(XYear,XMonth)
    write (99,"(i4,8i6)"), XDay,(CountStnIn(XDay,XHour),XHour=1,NHour)
end do
write (99,*)

CountAllIn(XYear,XMonth) = sum(CountStnIn)

end subroutine DumpYrMonStats

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

subroutine DumpOpStats

write (99,"(a)"), "The inputs (hours*days*stations), and outputs (stations)."
write (99,"(2a4,a8,8a6)"), "YEAR"," MON","  INPUTS",(OutSuffices(XVari),XVari=1,NVari)
do XYear = 1, NYear
  do XMonth = 1, NMonth
    write (99,"(i4,a4,i8,8i6)"), YearAD(XYear),MonthName(XMonth),CountAllIn(XYear,XMonth), &
    				(CountStnOut(XYear,XMonth,XVari),XVari=1,NVari)
  end do
end do
write (99,*)

end subroutine DumpOpStats

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

subroutine DumpClimStats

write (99,"(a)"), "The climatologies (stations)."
write (99,"(2a4,8a6)"), " MON"," SYN",(OutSuffices(XVari),XVari=1,NVari)
do XMonth = 1, NMonth
  do XHour = 1, NHour
    write (99,"(a4,i4,8i6)"), MonthName(XMonth),XHour,(CountStnTot(XMonth,XHour,XVari),XVari=1,NVari)
  end do
end do
write (99,*)

end subroutine DumpClimStats

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

subroutine Conclude

print*

close (99)

end subroutine Conclude

!*******************************************************************************
! ['.oct','.slp','.wnd','.vap','.tmp','.wdx','.wdy','.dir']

subroutine SelectArraysIn

if      (XVari.EQ.1) then
 	InAny => InOct
else if (XVari.EQ.2) then
 	InAny => InSLP
else if (XVari.EQ.3) then 
 	InAny => InWnd
else if (XVari.EQ.4) then 
 	InAny => InVap
else if (XVari.EQ.5) then 
 	InAny => InTmp
else if (XVari.EQ.6) then 
 	InAny => InWdx
else if (XVari.EQ.7) then 
 	InAny => InWdy
end if
  
end subroutine SelectArraysIn

!*******************************************************************************
! ['.oct','.slp','.wnd','.vap','.tmp','.wdx','.wdy','.dir']

subroutine SelectArraysClim

if      (XVari.EQ.1) then
 	ClimAny => ClimOct
else if (XVari.EQ.2) then
 	ClimAny => ClimSLP
else if (XVari.EQ.3) then 
 	ClimAny => ClimWnd
else if (XVari.EQ.4) then 
 	ClimAny => ClimVap
else if (XVari.EQ.5) then 
 	ClimAny => ClimTmp
else if (XVari.EQ.6) then 
 	ClimAny => ClimWdx
else if (XVari.EQ.7) then 
 	ClimAny => ClimWdy
end if
  
end subroutine SelectArraysClim

!*******************************************************************************
! ['.oct','.slp','.wnd','.vap','.tmp','.wdx','.wdy','.dir']

subroutine SelectArraysTot

if      (XVari.EQ.1) then
 	TotAny => TotOct
else if (XVari.EQ.2) then
 	TotAny => TotSLP
else if (XVari.EQ.3) then 
 	TotAny => TotWnd
else if (XVari.EQ.4) then 
 	TotAny => TotVap
else if (XVari.EQ.5) then 
 	TotAny => TotTmp
else if (XVari.EQ.6) then 
 	TotAny => TotWdx
else if (XVari.EQ.7) then 
 	TotAny => TotWdy
end if
  
end subroutine SelectArraysTot

!*******************************************************************************
! ['.oct','.slp','.wnd','.vap','.tmp','.wdx','.wdy','.dir']

subroutine SelectArraysOut

if      (XVari.EQ.1) then
 	OutAny => OutOct
else if (XVari.EQ.2) then
 	OutAny => OutSLP
else if (XVari.EQ.3) then 
 	OutAny => OutWnd
else if (XVari.EQ.4) then 
 	OutAny => OutVap
else if (XVari.EQ.5) then 
 	OutAny => OutTmp
else if (XVari.EQ.6) then 
 	OutAny => OutWdx
else if (XVari.EQ.7) then 
 	OutAny => OutWdy
else if (XVari.EQ.8) then 
 	OutAny => OutDir
end if
  
end subroutine SelectArraysOut

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

end program GetHahn