! crutsfiles.f90
! module in which standard save to and load from CRU ts file routines are held
! permitted: .cts .hdr .src .nrm (anything permitted for LoadCTS only)
! contains: LoadCTS,SaveCTS,MakeStnInfo,CombineStn
! data missing values: it is very important that all missing values in the files be set to -9999
! converted to use on beo1 with Portland Group compiler 

module CRUtsFiles

use FileNames

implicit none

contains

!*******************************************************************************
! all the arrays are alloc within the routine
! StnLocal added 1.5.02 = a9 local stn identifier field that prev was in StnInfo(*,8)
! Source (legacy) added 1.5.02 for master code files: NOT with Hulme & NmlData,requires HeadOnly
! SrcCode,SrcSuffix,SrcDate added 24.10.02 for revised master code files
! if Extra is set, the NStn is expanded by the size of Extra, and arrays filled with MVs
! if PhilJ is set, 2 more metadata are used (code for homog, 1st usable year),requires HeadOnly etc
! assume LatMissVal=-9999 ; LonMissVal=-99999 ; ElvMissVal=-999 ; DataMissVal=-9999
!	to change set LatMV,LonMV,ElvMV,DataMV
! assume Lat*100,Lon*100,Elv*1 ; to change set LatF,LonF,ElvF
! Legacy, when set, allows loading of old format for headers
! NmlData may be accompanied by NmlYr0, NmlYr1, NmlSrc, NmlInc, and most other options,
!   but NOT Data, YearAD, Hulme, Legacy, HeadOnly, Source
! DtbNormals (or .dtb) implies a database file with a normals line below the metadata line

subroutine LoadCTS (StnInfo,StnLocal,StnName,StnCty,Code,Lat,Lon,Elv,OldCode,Data,YearAD,&
        NmlData,DtbNormals,CallFile,Hulme,Legacy,HeadOnly,HeadForm,LongType,Silent,Extra,PhilJ, &
        YearADMin,YearADMax,Source,SrcCode,SrcSuffix,SrcDate, &
        LatMV,LonMV,ElvMV,DataMV,LatF,LonF,ElvF,NmlYr0,NmlYr1,NmlSrc,NmlInc)

real, pointer, dimension (:), optional		:: Lat,Lon,Elv	! additional vec form

integer, pointer, dimension (:,:,:), optional	:: Data		! data integers (not if HeadOnly)
integer, pointer, dimension (:,:), optional	:: NmlData	! data integers (DataMV okay)
integer, pointer, dimension (:,:), optional	:: DtbNormals	! data integers (DataMV okay)
integer, pointer, dimension (:,:), optional	:: Hulme	! annual column in Hulme only
integer, pointer, dimension (:,:)		:: StnInfo	! station key, stored as integers
integer, pointer, dimension (:), optional	:: Code		! additional vec form
integer, pointer, dimension (:), optional	:: OldCode	! additional vec form
integer, pointer, dimension (:), optional	:: YearAD	! not if HeadOnly or NmlData
integer, pointer, dimension (:), optional	:: NmlYr0, NmlYr1 ! normals period (StnN)
integer, pointer, dimension (:), optional	:: NmlSrc, NmlInc ! normals source code, % present
integer, pointer, dimension (:), optional	:: SrcCode,SrcDate ! source info

character (len=20), pointer, dimension (:), optional	:: Source
character (len=20), pointer, dimension (:)		:: StnName
character (len=13), pointer, dimension (:)		:: StnCty
character (len=09), pointer, dimension (:)		:: StnLocal
character (len=04), pointer, dimension (:), optional	:: SrcSuffix

integer, intent(in), optional 	:: LatMV,LonMV,ElvMV,DataMV	! raw missing values
integer, intent(in), optional 	:: LatF,LonF,ElvF		! Factors
integer, intent(in), optional	:: HeadOnly	! if set, load headers-only file
integer, intent(in), optional	:: LongType	! 2=reverse-sign,3=180-360-->-180-0
integer, intent(in), optional 	:: Legacy	! if set, load old header format	 
integer, intent(in), optional 	:: Silent	 
integer, intent(in), optional 	:: Extra	 
integer, intent(in), optional 	:: PhilJ	 
integer, intent(in), optional 	:: YearADMin,YearADMax	 ! ensure YearAD covers this range

character (len=80), intent(in), optional	:: CallFile	! filepath to load (optional)
character (len=80), intent(in), optional	:: HeadForm	! format of the headers

real, parameter :: MissVal = -999.0

integer, parameter :: DataMissVal = -9999

integer :: LatFactor,LonFactor,ElvFactor, Year0,Year1,SuffixBeg
integer :: ReadStatus,AllocStat
integer :: FileStnN, StnN, MonthN, YearN, LineN, InfoN
integer :: XFileStn, XStn, XMonth, XYear, XLine, XInfo, XYearAD
integer :: HWmo,HLat,HLong,HElev,HYear,HYear0,HYear1,HGcm,HSrc,HInc,HHomog,HFirst
integer :: QDatabase
integer :: LatMissVal,LonMissVal,ElvMissVal,RawMissVal

character (len=80) :: GivenFile, LoadFile, HeaderFormat,DataFormat,Trash
character (len=80) :: ReadStnName, ReadCtyName, ReadStnLocal, HSource
character (len= 9) :: HLocal
character (len= 4) :: Suffix
character (len= 1) :: Skip

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

if (present(Hulme).AND.present(HeadOnly)) print*, &
		"  > ##### ERROR: LoadCRUts: Call Clash: Hulme, HeadOnly #####"
if (present(NmlData)) then
  if (present(Data).OR.present(YearAD).OR.present(Hulme).OR.present(Legacy).OR. &
  		present(HeadOnly).OR.present(Source)) &
  		print*, "  > ##### ERROR: LoadCRUts: Call Clash: NmlData present #####"
else
  if (present(NmlYr0).OR.present(NmlYr1).OR.present(NmlSrc).OR.present(NmlInc)) &
  		print*, "  > ##### ERROR: LoadCRUts: Call Clash: NmlData missing #####"
end if

if (present(PhilJ).AND.(present(NmlData).OR.present(Source).OR.present(SrcCode).OR.present(Hulme))) &
  		print*, "  > ##### ERROR: LoadCRUts: Call Clash: PhilJ present #####"

QDatabase = 0
if (present(DtbNormals)) QDatabase = 1

MonthN = 12

LatMissVal=-9999 ; LonMissVal=-99999 ; ElvMissVal=-999		! usual/hulme hdr miss vals
if (present(Legacy)) then
  LatMissVal=-999 ; LonMissVal=-9999 ; ElvMissVal=-999		! legacy hdr miss vals
end if
if (present(LatMV)) LatMissVal = LatMV				! custom hdr miss vals
if (present(LonMV)) LonMissVal = LonMV
if (present(ElvMV)) ElvMissVal = ElvMV

RawMissVal = -9999
if (present(Hulme))  RawMissVal = -10
if (present(DataMV)) RawMissVal = DataMV

LatFactor=100 ; LonFactor=100 ; ElvFactor=1			! usual/hulme hdr factors
if (present(Legacy)) then
  LatFactor=10 ; LonFactor=10 ; ElvFactor=1			! legacy hdr factors
end if
if (present(LatF)) LatFactor = LatF				! custom hdr factors
if (present(LonF)) LonFactor = LonF
if (present(ElvF)) ElvFactor = ElvF

DataFormat="(i4,12i5)"						! standard data format
if (present(Hulme)) DataFormat="(i4,12i5,i6)"			! hulme data format

if (present(HeadForm)) then
  HeaderFormat = trim(HeadForm)
  Suffix = "    "
else if (present(NmlData)) then
  HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,i2,x,i3)"
  Suffix = ".nrm"
else if (present(Hulme)) then
  HeaderFormat = "(i7,i5,i6,i5,a20,a13,i4,i4,i7,a9)"
  Suffix = "    "
else if (present(Legacy)) then
  HeaderFormat = "(i7,i5,i6,i5,x,a20,x,a13,x,i4,x,i4,i7,a9)"
  Suffix = ".cts"
else
  if (present(Source)) then
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,a20)"
  else if (present(SrcCode)) then
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,i4,x,a4,x,i10)"
  else
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9)"
  end if
  
  if (present(HeadOnly)) then
   Suffix = ".hdr"
  else
   if (present(CallFile)) then
    SuffixBeg = index(CallFile,".",.TRUE.)
    if      (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".src") then
      Suffix = ".src"
    else if (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".dtb") then
      Suffix = ".dtb" ; QDatabase = 1
    else if (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".dts") then
      Suffix = ".dts" ; QDatabase = 1
    else
      Suffix = ".cts"
    end if
   else
    Suffix = ".cts"
   end if
  end if
end if

GivenFile = ""								! file to load
if (present(CallFile)) GivenFile = CallFile
if (GivenFile.EQ."") then
  print "(2a)", "  > Enter the file to load with suffix: ", Suffix  
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
end if
LoadFile = LoadPath (GivenFile,Suffix)

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

if (present(HeadOnly)) then
 StnN = LineN
else if (present(NmlData)) then
 StnN = LineN / 2
else    
 open  (2,file=LoadFile,status="old",access="sequential",form="formatted",action="read")    
 XLine = 1 ; Year0 = 1000000 ; Year1 = -1000000	; StnN = 0		! get StationN,YearN
 do
  if      (present(Source)) then
    read (2,HeaderFormat),HWmo,HLat,HLong,HElev,ReadStnName,ReadCtyName,HYear0,HYear1,HGcm,HLocal,HSource
  else if (present(PhilJ)) then
    read (2,HeaderFormat),HWmo,HLat,HLong,HElev,ReadStnName,ReadCtyName,HYear0,HYear1,HHomog,HFirst
  else
    read (2,HeaderFormat),HWmo,HLat,HLong,HElev,ReadStnName,ReadCtyName,HYear0,HYear1,HGcm,HLocal
  end if
  
  if (QDatabase.EQ.1) then
  	read (2,"(a1)"), Skip			! skip normals line
  	XLine = XLine + 1
  end if
  if (HYear0.LT.Year0) Year0 = HYear0
  if (HYear1.GT.Year1) Year1 = HYear1
  do XYear = HYear0,HYear1
    read (2,"(a1)"), Skip			! skip time-series lines
  end do
  XLine = XLine + 1 + (HYear1 - HYear0 + 1)
  StnN = StnN + 1
  if (XLine.GT.LineN) exit
 end do

 if (present(YearADMin)) then
   if (YearADMin.LT.Year0) Year0 = YearADMin
 end if
 if (present(YearADMax)) then
   if (YearADMax.GT.Year1) Year1 = YearADMax
 end if
 
 YearN = Year1 - Year0 + 1
 close (2)
end if

FileStnN = StnN
if (present(Extra)) StnN = StnN + Extra

if (present(NmlData)) then
  allocate (NmlData (MonthN,StnN),  stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: NmlData #####"
  NmlData = DataMissVal

  if (present(NmlSrc)) then
    allocate (NmlSrc  (StnN),  stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: NmlSrc #####"
    NmlSrc = DataMissVal
  end if
  if (present(NmlInc)) then
    allocate (NmlInc  (StnN),  stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: NmlInc #####"
    NmlInc = DataMissVal
  end if
else if (.not.present(HeadOnly)) then
  allocate (Data   (YearN,MonthN,StnN), &
	    YearAD (YearN),             stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Data #####"
  Data = DataMissVal
  
  do XYear = 1, YearN
    YearAD(XYear) = Year0 + XYear - 1
  end do
  
  if (present(DtbNormals)) then
    allocate (DtbNormals(MonthN,StnN), stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: DtbNormals #####"
    DtbNormals = DataMissVal
  end if
  
  if (present(Hulme)) then
    allocate (Hulme (YearN,StnN), stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Hulme #####"
    Hulme = DataMissVal
  end if
end if

allocate (StnInfo (StnN,8),            &
	  StnLocal(StnN),              &
	  StnName (StnN),              &
	  StnCty  (StnN),              stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Stn #####"
StnInfo = MissVal ; StnName = "UNKNOWN" ; StnCty = "UNKNOWN" ; StnLocal = "   nocode"

if (present(Source)) then
  allocate (Source(StnN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Source #####"
  Source = "SOURCE INFO"
else if (present(SrcCode)) then
  allocate (SrcCode(StnN), &
    	    SrcSuffix(StnN), &
    	    SrcDate(StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Src* #####"
  SrcCode = MissVal ; SrcSuffix = "NONE" ; SrcDate = MissVal
end if

open  (2,file=LoadFile,status="old",access="sequential",form="formatted",action="read")    
do XStn = 1, FileStnN
  if (present(NmlData)) then
    read (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
  			  (StnInfo(XStn,XInfo),XInfo=5,7),ReadStnLocal,HSrc,HInc
    if (present(NmlSrc)) NmlSrc(XStn)=HSrc
    if (present(NmlInc)) NmlInc(XStn)=HInc
  else if (present(Source)) then
    read (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
  			  (StnInfo(XStn,XInfo),XInfo=5,7),ReadStnLocal,Source(XStn)
  else if (present(SrcCode)) then
    read (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
  			  (StnInfo(XStn,XInfo),XInfo=5,7),ReadStnLocal, &
  			  SrcCode(XStn),SrcSuffix(XStn),SrcDate(XStn)
  else if (present(PhilJ)) then
    read (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
  			  (StnInfo(XStn,XInfo),XInfo=5,6),HHomog,HFirst
  else
    read (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
  			  (StnInfo(XStn,XInfo),XInfo=5,7),ReadStnLocal
!    write (99,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),ReadStnName,ReadCtyName, &
!  			  (StnInfo(XStn,XInfo),XInfo=5,7),ReadStnLocal	! @@@@@@@@@@@@@@@@@@@@@@@@
  end if
  
  StnLocal(XStn) = trim(ReadStnLocal)
  StnName(XStn)  = trim(ReadStnName)
  StnCty (XStn)  = trim(ReadCtyName)
  
  if (.not.present(HeadOnly)) then
    if (present(NmlData)) then
        read (2,DataFormat), HYear,(NmlData(XMonth,XStn),XMonth=1,MonthN)
    else if (QDatabase.EQ.1) then
      if (present(DtbNormals)) then
        read (2,"(i4,12i5)"), HYear,(DtbNormals(XMonth,XStn),XMonth=1,MonthN)
      else
  	read (2,"(a80)"), Trash			! skip normals line
      end if
      do XYear = (StnInfo(XStn,5)-Year0+1), (StnInfo(XStn,6)-Year0+1)
        read (2,DataFormat), HYear,(Data(XYear,XMonth,XStn),XMonth=1,MonthN)
      end do
    else if (present(Hulme)) then
      do XYear = (StnInfo(XStn,5)-Year0+1), (StnInfo(XStn,6)-Year0+1)
        read (2,DataFormat), HYear,(Data(XYear,XMonth,XStn),XMonth=1,MonthN),Hulme(XYear,XStn)
      end do
    else if (present(PhilJ)) then
      XYear = (StnInfo(XStn,5)-Year0+1)-1 ; XYearAD = StnInfo(XStn,5)-1
      do				
        XYear=XYear+1 ; XYearAD=XYearAD+1
        if (XYearAD.LT.HFirst) then
          read (2,"(a1)"), Skip			! ignore unsafe early period
        else					! read safe later period
          read (2,DataFormat), HYear,(Data(XYear,XMonth,XStn),XMonth=1,MonthN)
        end if
        if (XYear.EQ.(StnInfo(XStn,6)-Year0+1)) exit
      end do
      StnInfo(XStn,5) = HFirst
    else
      do XYear = (StnInfo(XStn,5)-Year0+1), (StnInfo(XStn,6)-Year0+1)
        read (2,DataFormat), HYear,(Data(XYear,XMonth,XStn),XMonth=1,MonthN)
      end do
    end if
  end if
end do
close (2)

if (present(Data))    where (Data.EQ.RawMissVal)       Data = DataMissVal
if (present(Hulme))   where (Hulme.EQ.RawMissVal)      Hulme = DataMissVal
if (present(NmlData)) where (NmlData.EQ.RawMissVal)    NmlData = DataMissVal
if (present(DtbNormals))   where (DtbNormals.EQ.RawMissVal) DtbNormals = DataMissVal

if (present(LongType)) then		! convert non-standard longitudes
  if      (LongType.EQ. 2) then
    do XStn = 1, FileStnN
      if (StnInfo(XStn,3).NE.LonMissVal) StnInfo(XStn,3) = StnInfo(XStn,3) * (0-1)
    end do
  else if (LongType.EQ. 3) then
    do XStn = 1, FileStnN
      if (StnInfo(XStn,3).GT.(180*LonFactor).AND.StnInfo(XStn,3).NE.LonMissVal) &
      				StnInfo(XStn,3) = -1 * ((360*LonFactor) - StnInfo(XStn,3))
    end do
  end if
end if
  
if (present(Code)) then
  allocate (Code (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Code #####"
  Code = MissVal
  do XStn = 1, FileStnN
    Code (XStn) = StnInfo (XStn,1)
  end do
end if

if (present(OldCode)) then
  allocate (OldCode (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: OldCode #####"
  OldCode = MissVal
  do XStn = 1, FileStnN
    OldCode (XStn) = StnInfo (XStn,7)
  end do
end if

if (present(Lat)) then
  allocate (Lat (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Lat #####"
  Lat = MissVal
  do XStn = 1, FileStnN
    if (StnInfo(XStn,2).NE.LatMissVal) Lat (XStn) = real(StnInfo(XStn,2)) / real(LatFactor)
  end do
end if

if (present(Lon)) then
  allocate (Lon (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Lon #####"  
  Lon = MissVal
  do XStn = 1, FileStnN
    if (StnInfo(XStn,3).NE.LonMissVal) Lon (XStn) = real(StnInfo(XStn,3)) / real(LonFactor)
  end do
end if

if (present(Elv)) then
  allocate (Elv (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: Elv #####"
  Elv = MissVal
  do XStn = 1, FileStnN
    if (StnInfo(XStn,4).NE.ElvMissVal) Elv (XStn) = real(StnInfo(XStn,4)) / real(ElvFactor)
  end do
end if

if (present(NmlYr0)) then
  allocate (NmlYr0 (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: NmlYr0 #####"
  NmlYr0 = MissVal
  do XStn = 1, FileStnN
    NmlYr0 (XStn) = StnInfo(XStn,5)
  end do
end if

if (present(NmlYr1)) then
  allocate (NmlYr1 (StnN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadCRUts: Allocation failure: NmlYr1 #####"
  NmlYr1 = MissVal
  do XStn = 1, FileStnN
    NmlYr1 (XStn) = StnInfo(XStn,6)
  end do
end if

if (.not.present(Silent)) print*, "  > Stations in file total: ", FileStnN

end subroutine LoadCTS

!*******************************************************************************
! to save a stn, it must have: StnInfo(XStn,XInfo), where XInfo=1
! SaveCTS will only save the station if it has a valid column 1 value 
! 	if valid, this routine expects cols 5,6 to also be filled
! input missing values MUST be = -9999 ; to save miss val with different code, use DataMV
! StnLocal added 1.5.02 = a9 local stn identifier field that prev was in StnInfo(*,8)
! Source (optional) added 1.5.02 = a20 for master code files: requires HeadOnly
! NmlData must be accompanied by NmlSrc, NmlInc, and may take most other options,
!   but NOT Data, YearAD, Hulme, HeadOnly, Source
! to save in code order, fill Order with the stn indices to use
! DtbNormals (NOT .dtb) implies a database file with a normals line below the metadata line

subroutine SaveCTS (StnInfo,StnLocal,StnName,StnCty,Data,YearAD,NmlData,DtbNormals,CallFile, &
		    Hulme,HeadOnly,HeadForm,Silent,Source,Order, &
		    SrcCode,SrcSuffix,SrcDate,DataMV,NmlSrc,NmlInc)

integer, pointer, dimension (:,:,:), optional	:: Data		! data, stored as integers (not if HeadOnly)
integer, pointer, dimension (:,:), optional	:: NmlData	! data integers (DataMV okay)
integer, pointer, dimension (:,:), optional	:: DtbNormals	! data integers (DataMV okay)
integer, pointer, dimension (:,:), optional	:: Hulme	! annual column in Hulme
integer, pointer, dimension (:,:)		:: StnInfo	! station key, stored as integers
integer, pointer, dimension (:), optional	:: YearAD	! not if HeadOnly
integer, pointer, dimension (:), optional	:: NmlSrc, NmlInc ! normals source code, % present
integer, pointer, dimension (:), optional	:: SrcCode,SrcDate ! source info
integer, pointer, dimension (:), optional	:: Order 	! num order of stn codes

character (len=20), pointer, dimension (:), optional	:: Source
character (len=20), pointer, dimension (:)		:: StnName
character (len=13), pointer, dimension (:)		:: StnCty
character (len=09), pointer, dimension (:)		:: StnLocal
character (len=04), pointer, dimension (:), optional	:: SrcSuffix

integer, intent(in), optional 			:: Silent	 
integer, intent(in), optional 			:: DataMV	 
integer, intent(in), optional			:: HeadOnly	! if set, save headers-only file

character (len=80), intent(in), optional	:: HeadForm	! format of the headers
character (len=80), intent(in), optional	:: CallFile	! filepath to save (optional)

real, parameter :: MissVal = -999.0

integer, parameter :: DataMissVal = -9999

integer :: ReadStatus,AllocStat
integer :: StnN, MonthN, YearN, LineN, InfoN, SaveStnN, AllStnN
integer :: XStn, XMonth, XYear, XLine, XInfo, XSaveStn, XAllStn
integer :: HWmo,HLat,HLong,HElev,HYear,HYear0,HYear1,HGcm,HLocal
integer :: Year0,Year1,SuffixBeg

character (len=80) :: GivenFile, SaveFile, HeaderFormat
character (len=20) :: HName, HCty, DataFormat
character (len= 4) :: Suffix
character (len= 1) :: Skip

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

if (present(Hulme).AND.present(HeadOnly)) print*, &
		"  > ##### ERROR: SaveCRUts: Call Clash: Hulme, HeadOnly #####"

if (present(NmlData)) then
  if (present(Data).OR.present(YearAD).OR.present(Hulme).OR. &
  		present(HeadOnly).OR.present(Source)) &
  		print*, "  > ##### ERROR: SaveCRUts: Call Clash: NmlData present #####"
  if (present(NmlSrc).AND.present(NmlInc)) then
  else
  		print*, "  > ##### ERROR: SaveCRUts: Call Clash: NmlSrc or NmlInc missing #####"
  end if
else
  if (present(NmlSrc).OR.present(NmlInc)) &
  		print*, "  > ##### ERROR: SaveCRUts: Call Clash: NmlData missing #####"
end if

MonthN = 12 ; StnN = size (StnInfo,1)

DataFormat="(i4,12i5)"
if (present(Hulme)) DataFormat="(i4,12i5,i6)"

if (present(HeadForm)) then
  HeaderFormat = trim(HeadForm)
else if (present(NmlData)) then
  HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,i2,x,i3)"
else if (present(Hulme)) then
  HeaderFormat = "(i7,i5,i6,i5,a20,a13,i4,i4,i7,a9)"
else
  if (present(Source)) then
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,a20)"
  else if (present(SrcCode)) then				! note '0' + date = 10 char !!
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,i4,x,a4,x,a1,i9)"
  else
    HeaderFormat = "(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9)"
  end if
end if

if (present(HeadOnly)) then
  Suffix = ".hdr"
else if (present(NmlData)) then
  Suffix = ".nrm"
else
  if (present(CallFile)) then
    SuffixBeg = index(CallFile,".",.TRUE.)
    if      (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".src") then
      Suffix = ".src"
    else if (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".dtb") then
      Suffix = ".dtb" ; if (.not.present(DtbNormals)) print*, &
      			"  > ##### ERROR: SaveCRUts: DtbNormals not in save to .dtb #####"
    else if (CallFile(SuffixBeg:(SuffixBeg+3)).EQ.".dts") then
      Suffix = ".dts" ; if (.not.present(DtbNormals)) print*, &
      			"  > ##### ERROR: SaveCRUts: DtbNormals not in save to .dtb #####"
    else
      Suffix = ".cts"
    end if
  else
    Suffix = ".cts"
  end if
  
  YearN = size(Data,1)
  
  if (present(DataMV)) then
    if (present(Data)) where (Data.EQ.DataMissVal) Data = DataMV
    if (present(Hulme)) where (Hulme.EQ.DataMissVal) Hulme = DataMV
    if (present(NmlData)) where (NmlData.EQ.DataMissVal) NmlData = DataMV
    if (present(DtbNormals)) where (DtbNormals.EQ.DataMissVal) DtbNormals = DataMV
  end if
end if

GivenFile = ""								! file to save
if (present(CallFile)) GivenFile = CallFile
if (trim(GivenFile).EQ."") then
  print "(2a)", "  > Enter the file to save with suffix: ", Suffix  
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
  end do
end if
SaveFile = SavePath (GivenFile,Suffix)

SaveStnN = StnN

open  (2,file=SaveFile,status="replace",access="sequential",form="formatted",action="write")    
do XAllStn = 1, StnN
 if (present(Order)) then
   XStn = Order(XAllStn)
 else
   XStn = XAllStn
 end if
 
 if (XStn.GE.1.AND.XStn.LE.StnN) then
  if (StnInfo(XStn,1).NE.MissVal) then
   if (present(NmlData)) then
    write (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),StnName(XStn),StnCty(XStn), &
  			 (StnInfo(XStn,XInfo),XInfo=5,7),StnLocal(XStn),NmlSrc(XStn),NmlInc(XStn)
   else if (present(Source)) then
    write (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),StnName(XStn),StnCty(XStn), &
  			 (StnInfo(XStn,XInfo),XInfo=5,7),StnLocal(XStn),Source(XStn)
   else if (present(SrcCode)) then
    write (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),StnName(XStn),StnCty(XStn), &
  			 (StnInfo(XStn,XInfo),XInfo=5,7),StnLocal(XStn), &
  			 SrcCode(XStn),SrcSuffix(XStn),'0',SrcDate(XStn)
   else
    write (2,HeaderFormat),(StnInfo(XStn,XInfo),XInfo=1,4),StnName(XStn),StnCty(XStn), &
  			 (StnInfo(XStn,XInfo),XInfo=5,7),StnLocal(XStn)
   end if
   
   if (present(DtbNormals)) write (2,"(a4,12i5)"), "6190",(DtbNormals(XMonth,XStn),XMonth=1,MonthN)

   if (.not.present(HeadOnly)) then
    if (present(NmlData)) then
        HYear = (100*mod(StnInfo(XStn,5),100)) + mod(StnInfo(XStn,6),100) 
        write (2,DataFormat), HYear,(NmlData(XMonth,XStn),XMonth=1,MonthN)
    else if (present(Hulme)) then
      do XYear = (StnInfo(XStn,5)-YearAD(1)+1), (StnInfo(XStn,6)-YearAD(1)+1)
        write (2,DataFormat), YearAD(XYear), (Data(XYear,XMonth,XStn),XMonth=1,MonthN), &
        			Hulme(XYear,XStn)
      end do
    else
      do XYear = (StnInfo(XStn,5)-YearAD(1)+1), (StnInfo(XStn,6)-YearAD(1)+1)
        write (2,DataFormat), YearAD(XYear), (Data(XYear,XMonth,XStn),XMonth=1,MonthN)
      end do
    end if
   end if
  else
   SaveStnN = SaveStnN - 1
  end if
 else
   SaveStnN = SaveStnN - 1
 end if
end do
close (2)

if (.not.present(Silent)) print*, "  > Stations in file total: ", SaveStnN

if (present(DataMV)) then
    if (present(Data)) where (Data.EQ.DataMV) Data = DataMissVal
    if (present(Hulme)) where (Hulme.EQ.DataMV) Hulme = DataMissVal
    if (present(NmlData)) where (NmlData.EQ.DataMV) NmlData = DataMissVal
    if (present(DtbNormals)) where (DtbNormals.EQ.DataMV) DtbNormals = DataMissVal
end if

end subroutine SaveCTS

!*******************************************************************************
! MAY BE NECESSARY TO ALLOC STNINFO PRIOR TO CALL
! make StnInfo for CRU ts file to be saved, using call info
! needs: WMOcode,Lat,Long,Elev (all vectors, StnN)
! needs: at least one of [Data,YearAD],[YearAD0,YearAD1]: [YearAD0,YearAD1] takes precedence
! needs: QNoPer (if no period obtainable: 0=save-all, 1=no-save)
! returns: StnInfo 
! SaveCTS will only save the station if it has a valid column 1 value 
! 	if valid, that routine expects cols 5,6 to also be filled
! data missing value MUST be -9999
! assume convert to Lat*100,Lon*100,Elv*1 ; to change set LatF,LonF,ElvF
! for .nrm files, supply NmlYr0, NmlYr1 for YearAD0,YearAD1
! for .dtb files, supply DtbNormals if supplying Data
! QBlank=1 if there are no stns to be saved

subroutine MakeStnInfo (StnInfo,WMOcode,OldCode,Lat,Long,Elev,QNoPer,QLongType, &
			YearAD0,YearAD1,YearAD,Data,DtbNormals,&
			Silent,QBlank,LatMV,LonMV,ElvMV,LatF,LonF,ElvF)

real, pointer, dimension (:)			:: Lat,Long,Elev

integer, pointer, dimension (:,:,:), optional	:: Data
integer, pointer, dimension (:,:), optional	:: DtbNormals
integer, pointer, dimension (:,:)		:: StnInfo
integer, pointer, dimension (:)			:: WMOcode,OldCode
integer, pointer, dimension (:), optional	:: YearAD0,YearAD1,YearAD

logical, allocatable, dimension (:)		:: Mask

character (len=20), pointer, dimension (:)	:: StnName
character (len=13), pointer, dimension (:)	:: CtyName

logical, intent(out), optional  :: QBlank

integer, intent(in), optional 	:: LatMV,LonMV,ElvMV	! missing values
integer, intent(in) 		:: QNoPer
integer, intent(in), optional	:: QLongType
integer, intent(in), optional 	:: Silent	 
integer, intent(in), optional 	:: LatF,LonF,ElvF	! Factors

real, parameter :: MissVal = -999.0

integer, parameter :: DataMissVal = -9999

integer :: LatFactor,LonFactor,ElvFactor
integer :: LatMissVal,LonMissVal,ElvMissVal
integer :: ReadStatus, AllocStat
integer :: StnN,XStn, YearN,XYear, MonthN,XMonth, SaveStnN,XSaveStn,BlankStnN,NoCodeStnN
integer :: YearADMin,YearADMax,Year

StnN = size (WMOCode,1) ; MonthN = 12 ; BlankStnN = 0 ; NoCodeStnN = 0 ; SaveStnN = 0

LatMissVal=-9999 ; LonMissVal=-99999 ; ElvMissVal=-999		! usual/hulme hdr miss vals
if (present(LatMV)) LatMissVal = LatMV				! custom hdr miss vals
if (present(LonMV)) LonMissVal = LonMV
if (present(ElvMV)) ElvMissVal = ElvMV

LatFactor=100 ; LonFactor=100 ; ElvFactor=1			! usual/hulme hdr factors
if (present(LatF)) LatFactor = LatF				! custom hdr factors
if (present(LonF)) LonFactor = LonF
if (present(ElvF)) ElvFactor = ElvF

if (.NOT.associated(StnINfo)) then
  allocate (StnInfo(StnN,8), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeStnInfo: alloc StnInfo #####"
end if
allocate (Mask   (StnN),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeStnInfo: alloc Mask #####"

StnInfo = -999
do XStn = 1, StnN
  StnInfo(XStn,2)=LatMissVal ; StnInfo(XStn,3)=LonMissVal ; StnInfo(XStn,4)=ElvMissVal
end do

do XStn = 1, StnN
 if (WMOCode(XStn).NE.MissVal) StnInfo(XStn,1)=WMOCode(XStn)
 if (OldCode(XStn).NE.MissVal) StnInfo(XStn,7)=OldCode(XStn)
 if (Lat (XStn).NE.MissVal)    StnInfo(XStn,2)=nint(Lat (XStn)*real(LatFactor))
 if (Long(XStn).NE.MissVal)    StnInfo(XStn,3)=nint(Long(XStn)*real(LonFactor))
 if (Elev(XStn).NE.MissVal)    StnInfo(XStn,4)=nint(Elev(XStn)*real(ElvFactor))
end do

if (present(QLongType)) then		! convert non-standard longitudes
  if      (QLongType.EQ. 2) then
    do XStn = 1, StnN
      if (StnInfo(XStn,3).NE.LonMissVal) StnInfo(XStn,3) = StnInfo(XStn,3) * (0-1)
    end do
  else if (QLongType.EQ. 3) then
    do XStn = 1, StnN
      if (StnInfo(XStn,3).LT.0.AND.StnInfo(XStn,3).NE.LonMissVal) &
      				StnInfo(XStn,3) = 360 - StnInfo(XStn,3)
    end do
  end if
end if
  
if (present(YearAD0).AND.present(YearAD1)) then
 Mask=.FALSE.
 where (YearAD0.NE.MissVal) Mask=.TRUE.
 YearADMin = minval(YearAD0,Mask)
 YearADMin = YearAD0(YearADMin)
 
 Mask=.FALSE.
 where (YearAD1.NE.MissVal) Mask=.TRUE.
 YearADMax = maxval(YearAD1,Mask)
 YearADMax = YearAD1(YearADMax)
 
 do XStn = 1, StnN
  if (StnInfo(XStn,1).NE.MissVal) then
   if (YearAD0(XStn).NE.MissVal.AND.YearAD1(XStn).NE.MissVal) then
  	StnInfo(XStn,5)=YearAD0(XStn)
  	StnInfo(XStn,6)=YearAD1(XStn)
  	SaveStnN = SaveStnN + 1
   else
  	if (QNoPer.EQ.0) then
  		StnInfo(XStn,5)=YearAD0(XStn)
  		StnInfo(XStn,6)=YearAD1(XStn)
  		BlankStnN = BlankStnN + 1
  		SaveStnN = SaveStnN + 1
  	else
  		StnInfo(XStn,1)=MissVal
  		StnInfo(XStn,5)=MissVal
  		StnInfo(XStn,6)=MissVal
  		BlankStnN = BlankStnN + 1
  	end if
   end if
  else
   NoCodeStnN = NoCodeStnN + 1
  end if
 end do
 
 if (.not.present(Silent)) print "(a,3i6)", "   > Totals: no-code,no-data,save: ", &
 				NoCodeStnN,BlankStnN,SaveStnN 
else if (present(Data).AND.present(YearAD)) then
 YearN = size(YearAD,1) 
 YearADMin = YearAD(1) ; YearADMax = YearAD(YearN)
 
 do XStn = 1, StnN
  if (StnInfo(XStn,1).NE.MissVal) then
   XYear = 0					! find station's first year with data
   do
     XYear = XYear + 1

     XMonth = 0
     do
       XMonth = XMonth + 1
       if (Data(XYear,XMonth,XStn).NE.DataMissVal) StnInfo(XStn,5)=YearAD(XYear)
       if (StnInfo(XStn,5).NE.MissVal) exit
       if (XMonth.EQ.MonthN) exit
     end do
     
     if (StnInfo(XStn,5).NE.MissVal) exit
     if (XYear.EQ.YearN) exit
   end do
   
   if (StnInfo(XStn,5).NE.MissVal) then		! got first yr, find last yr
    SaveStnN = SaveStnN + 1
    XYear = YearN + 1				
    do
     XYear = XYear - 1

     XMonth = 0
     do
       XMonth = XMonth + 1
       if (Data(XYear,XMonth,XStn).NE.DataMissVal) StnInfo(XStn,6)=YearAD(XYear)
       if (StnInfo(XStn,6).NE.MissVal) exit
       if (XMonth.EQ.MonthN) exit
     end do
     
     if (StnInfo(XStn,6).NE.MissVal) exit
     if (XYear.EQ.1) exit
    end do   
   else if (present(DtbNormals)) then		! no ts data, seek normals
     XMonth = 0
     do
       XMonth = XMonth + 1
       if (DtbNormals(XMonth,XStn).NE.DataMissVal) then
       		StnInfo(XStn,5)=YearAD(1) ; StnInfo(XStn,6)=YearAD(1)
       end if
       if (StnInfo(XStn,5).NE.MissVal) exit
       if (XMonth.EQ.MonthN) exit
     end do
   end if
   
   if (StnInfo(XStn,5).EQ.MissVal) then		! no data whatsoever
    if (QNoPer.EQ.0) then				! set to full range
  		StnInfo(XStn,5)=YearAD(1)
  		StnInfo(XStn,6)=YearAD(YearN)
  		BlankStnN = BlankStnN + 1
  		SaveStnN = SaveStnN + 1
    else						! set to missing
  		BlankStnN = BlankStnN + 1
  		StnInfo(XStn,1)=MissVal
  		StnInfo(XStn,5)=MissVal
  		StnInfo(XStn,6)=MissVal
    end if    
   end if
  else
   NoCodeStnN = NoCodeStnN + 1
  end if
 end do
 
 if (.not.present(Silent)) print "(a,3i6)", "   > Totals: no-code,no-data,save: ", &
 				NoCodeStnN,BlankStnN,SaveStnN 
 if (present(QBlank)) then
   QBlank=.TRUE. ; XStn=0
   do
     XStn=XStn+1
     if (StnInfo(XStn,1).NE.MissVal) QBlank=.FALSE.
     if (XStn.EQ.StnN) exit
     if (.NOT.QBlank) exit
   end do
 end if
else
 print*, "  > ***** PROBLEM: MakeStnInfo: Insufficient call: no save *****"
 StnInfo=MissVal
end if

end subroutine MakeStnInfo

!*******************************************************************************
! takes two independent station code vectors, 
!	1. combines them (.OR.) into CStn (size NCStn)
! 	2. provides vectors CfromA,CfromB (size NCStn) of indices from AStn,BStn

subroutine CombineStn (AStn,BStn,CStn,CfromA,CfromB)

integer, pointer, dimension (:) 	:: AStn,BStn,CStn,CfromA,CfromB
integer, allocatable, dimension (:) 	:: CalcCStn,CalcCfromA,CalcCfromB

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: XAStn,XBStn,XCStn
integer :: NAStn,NBStn,NCStn

NAStn=size(AStn) ; NBStn=size(BStn)
allocate (CalcCStn   (NAStn+NBStn), &
	  CalcCfromA (NAStn+NBStn), &
	  CalcCfromB (NAStn+NBStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineStn: Allocation failure: Calc #####"
CalcCStn=MissVal ; CalcCfromA=MissVal ; CalcCfromB=MissVal

XAStn=1 ; XBStn=1 ; XCStn=0
do  
  if (XAStn.LE.NAStn) then
    if (XBStn.LE.NBStn) then
      if      (AStn(XAStn).LT.BStn(XBStn)) then
        XCStn = XCStn + 1
	CalcCStn(XCStn)=AStn(XAStn) ; CalcCfromA(XCStn)=XAStn 
	XAStn = XAStn + 1
      else if (AStn(XAStn).EQ.BStn(XBStn)) then
        XCStn = XCStn + 1
	CalcCStn(XCStn)=AStn(XAStn) ; CalcCfromA(XCStn)=XAStn ; CalcCfromB(XCStn)=XBStn 
        XBStn = XBStn + 1
        XAStn = XAStn + 1
      else if (AStn(XAStn).GT.BStn(XBStn)) then
        XCStn = XCStn + 1
	CalcCStn(XCStn)=BStn(XBStn) ; CalcCfromB(XCStn)=XBStn 
        XBStn = XBStn + 1
      end if
    else
      XCStn = XCStn + 1
      CalcCStn(XCStn)=AStn(XAStn) ; CalcCfromA(XCStn)=XAStn 
      XAStn = XAStn + 1
    end if
  else
    if (XBStn.LE.NBStn) then
      XCStn = XCStn + 1
      CalcCStn(XCStn)=BStn(XBStn) ; CalcCfromB(XCStn)=XBStn 
      XBStn = XBStn + 1
    end if
  end if
  
  if (XAStn.GT.NAStn.AND.XBStn.GT.NBStn) exit
end do

NCStn=XCStn
allocate (CStn   (NCStn), &
	  CfromA (NCStn), &
	  CfromB (NCStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineStn: Allocation failure: C* #####"
do XCStn = 1, NCStn
  CStn   (XCStn) = CalcCStn   (XCStn)
  CfromA (XCStn) = CalcCfromA (XCStn)
  CfromB (XCStn) = CalcCfromB (XCStn)
end do

deallocate (CalcCStn,CalcCfromA,CalcCfromB, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CombineStn: Deallocation failure: Calc* #####"

end subroutine CombineStn

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

end module CRUtsFiles