! grimfiles.f90
! module holding standard save and load GRIM (GRIdded Monthly) routines
! SaveGrim,LoadGrim,SaveGrimMonthPost,LoadGrimMonthPrep,SaveGrimMonth,LoadGrimMonth
!   SaveGrip,LoadGrip,CheckGridAB
! written by Tim Mitchell on 22.03.01
! last modified on 20.11.01
! This file has been modified by Soenke Zaehle (PIK) on 30.07.02 for use with XLF compiler
! (PIK). changes are: 
! - removed "," after read and write statements
! - replaced SCRATCH statements with REPLACE statements

module GrimFiles

use FileNames
use Time

implicit none

contains

!*******************************************************************************
! this ensures that the data is saved in the order implicit in the grid dimensions:
!   i.e. long (X) varies slowest (-180...180), and lat (Y) quickest (-90...90)

subroutine SaveGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix,NoZip,Silent)

real, pointer, dimension (:,:,:)	:: Data

integer, pointer, dimension (:,:)	:: Grid
integer, pointer, dimension (:)		:: YearAD

real, dimension (4)			:: Bounds

integer, dimension (12)			:: YearData

character (len=80), dimension (5)	:: Headers

integer, optional, intent (in)		:: NoZip,Silent

character (len=80), intent (in)		:: Info, CallFile
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: SaveSuffix

real, parameter :: MissVal = -999.0

real :: Divisor

integer :: YearN,MonthN,BoxN,ExeN,WyeN
integer :: XYear,XMonth,XBox,XHeader,XExe,XWye
integer :: AllocStat, TotExceed

character (len=80) 		:: SaveFile, Variable

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

YearN = size (Data,1) ; MonthN = size (Data,2) ; BoxN = size (Data,3) 
ExeN  = size (Grid,1) ; WyeN   = size (Grid,2) ; TotExceed = 0

if (size(YearAD).NE.YearN) print*, "  > ##### ERROR: SaveGrim: incompatible arrays #####"
if (MonthN      .NE.   12) print*, "  > ##### ERROR: SaveGrim: weird calendar #####"

!if (Bounds(1).LT.-180.OR.Bounds(2).GT.180.OR.Bounds(1).GT.Bounds(2)) &
!	print*, "  > ##### ERROR: SaveGrim: longitudes not within spherical geometry #####"
!
!if (Bounds(3).LT.-90.OR.Bounds(4).GT.90.OR.Bounds(3).GT.Bounds(4)) &
!	print*, "  > ##### ERROR: SaveGrim: latitudes  not within spherical geometry #####"

if (Bounds(2).EQ.Bounds(1)) print*, "  > ##### ERROR: SaveGrim: zero longitude range #####"
if (Bounds(4).EQ.Bounds(3)) print*, "  > ##### ERROR: SaveGrim: zero latitude range #####"

call ReviewCall (CallFile,CallSuffix,SaveFile,SaveSuffix,2)	! checks for file/suffix consistency

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

call CheckVariSuffix (SaveSuffix,Variable,Divisor)		! get Divisor

if (.not.present(Silent)) print "(2a)", "   > ", trim(Info)

open  (2,file=SaveFile,status="replace",access="sequential",form="formatted",action="write")

call WriteHeaders (SaveSuffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD(1),YearAD(YearN),"grim")

do XExe = 1, ExeN
 do XWye = 1, WyeN
  if (Grid(XExe,XWye).NE.MissVal) then
   write (2,"(a9,i4,a1,i4)") "Grid-ref=", XExe, ",", XWye
  
   do XYear = 1, YearN
    do XMonth = 1, MonthN					
      if (Data(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal) then
      		YearData(XMonth) = nint(Data(XYear,XMonth,Grid(XExe,XWye)) / Divisor)
      		if (YearData(XMonth).GT.99999.OR.YearData(XMonth).LT.-9999) then
      		  YearData(XMonth) = nint(MissVal)
      		  TotExceed = TotExceed + 1
      		end if
      else
      		YearData(XMonth) = nint(MissVal)
      end if
    end do
    
    write (2,"(12i5)") (YearData(XMonth), XMonth=1,MonthN)
   end do
  end if
 end do
end do

close (2)

if (TotExceed.GT.0) print "(a,i9)", "   > Box-months outside range: ", TotExceed

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

if (present(NoZip)) then
  if (.not.present(Silent)) print*, "  > The file is not being compressed."
else
  if (.not.present(Silent)) print*, "  > Compressing in background..."
  call system ('compress ' // SaveFile // ' &')
end if

end subroutine SaveGrim

!*******************************************************************************
! when using this routine, do not assume that this Grid and your Grid are identical
! however, if the same cells are non-missing in the original saves to grim,
!   the SaveGrim procedure ensures that the Grids obtained will be the same from the two loads from grim
! Data (A) may be defined (=Y) or not (=N) in the call; MasterYearAD (B) may (=Y) or may not (=N) be specified
!   if A=Y,B=Y: loading requires sizes A=B, and common period between B and file
!   if A=Y,B=N: loading requires sizes A=file
!   if A=N,B=Y: loading requires common period between B and file, Data=(MasterYearN,12,BoxN)
!   if A=N,B=N: loaded automatically, Data=(FileYearN,12,BoxN)
! if Data is defined in the call, it must also have already been initialised if that is necessary
! FileYearAD IS the genuine FileYearAD, and may differ from Data and MasterYearAD
! use Override with caution! it fills in FileYear0 (and hence FileYear1)

subroutine LoadGrim (Data,Grid,FileYearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix,&
			MasterYearAD,Quiet,OverRide)

real, pointer, dimension (:,:,:)	:: Data

integer, pointer, dimension (:,:)		:: Grid
integer, pointer, dimension (:)			:: FileYearAD
integer, pointer, optional, dimension (:)	:: MasterYearAD	 ! may or may not be included: see top

integer, dimension (12)			:: YearData

character (len=80), dimension (5)	:: Headers

real, dimension (4)			:: Bounds

integer, intent(in), optional		:: Quiet,OverRide

character (len=80), intent (out)	:: Info
character (len=80), intent (in)		:: CallFile
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: LoadSuffix

real, parameter :: MissVal = -999.0

real :: Multiplier

integer :: FileYearN,MasterYearN,MonthN,BoxN,ExeN,WyeN
integer :: XFileYear,XMasterYear,XMonth,XBox,XHeader,XExe,XWye
integer :: AllocStat
integer :: LoadFileLen,FileYearAD0,FileYearAD1,FileYear0,FileYear1,MasterYear0,MasterYear1
integer :: QZip
integer :: DataEn,DataMiss

character (len=80) :: LoadFile, Variable, Trash

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

call ReviewCall (CallFile,CallSuffix,LoadFile,LoadSuffix,1)	! checks for file/suffix consistency

LoadFileLen = len_trim(LoadFile)
if (LoadFileLen.GT.1.AND.LoadFile((LoadFileLen-1):LoadFileLen).EQ.".Z") then
  QZip = 2							! file is zipped
  if (.not.present(Quiet)) print*, "  > Uncompressing on the fly..."
  call system ('uncompress ' // LoadFile)
  LoadFile ((LoadFileLen-1):LoadFileLen) = "  "			! change filename to the unzipped file
else
  QZip = 1							! file already unzipped
end if

call CheckVariSuffix (LoadSuffix,Variable,Multiplier)		! get multiplier

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

open  (2,file=LoadFile,status="old",access="sequential",form="formatted",action="read")

call ReadHeaders (LoadSuffix,Info,Bounds,ExeN,WyeN,BoxN,FileYearAD0,FileYearAD1,"grim")

if (Bounds(2).EQ.Bounds(1)) print*, "  > ##### ERROR: LoadGrim: zero longitude range #####"
if (Bounds(4).EQ.Bounds(3)) print*, "  > ##### ERROR: LoadGrim: zero latitude range #####"

FileYearN  = FileYearAD1 - FileYearAD0 + 1
MonthN = 12

if      (present(MasterYearAD).AND.associated(Data)) then
  if (size(Data,1).NE.size(MasterYearAD,1)) then
  	print*, "  > ##### ERROR: LoadGrim: call mismatch #####"
  	MasterYearN = MissVal
  else
        MasterYearN = size(Data,1)
  end if
else if (present(MasterYearAD)) then
        MasterYearN = size(MasterYearAD,1)
else if (associated(Data)) then
        MasterYearN = size(Data,1)
else
	MasterYearN = FileYearN
end if

if (associated(Data)) then
  if (size(Data,2).NE.MonthN)   print*, "  > ##### ERROR: LoadGrim: Months mismatch #####"
  if (size(Data,3).NE.BoxN)     print*, "  > ##### ERROR: LoadGrim: Domain mismatch #####"
else
  allocate (Data (MasterYearN,MonthN,BoxN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGRIM: Allocation failure: Data #####"
  Data = MissVal
end if

allocate (Grid       (ExeN,WyeN), &
	  FileYearAD (FileYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGRIM: Allocation failure: All #####"
Grid = MissVal

if (present(OverRide)) then
  FileYearAD0=OverRide ; FileYearAD1=OverRide+FileYearN-1
end if
do XFileYear = 1, FileYearN
  FileYearAD(XFileYear) = FileYearAD0 + XFileYear - 1
end do

if (present(MasterYearAD)) then
  call CommonVecPer (FileYearAD,MasterYearAD,FileYear0,FileYear1,MasterYear0,MasterYear1)
  
  if (FileYear0.EQ.MissVal) then
    print "(a)", "   > ##### ERROR: LoadGrim: The file and spec are mismatched in period #####"
    FileYear0 = MissVal ; FileYear1 = MissVal ; MasterYear0 = MissVal ; MasterYear1 = MissVal
  else
    if (.not.present(Quiet)) print "(a,i4,a1,i4))", "   > Period loading: ", FileYearAD(FileYear0), &
    							"-", FileYearAD(FileYear1)
  end if
else
  if (size(Data,1).EQ.FileYearN) then
    FileYear0 = 1 ; FileYear1 = FileYearN ; MasterYear0 = 1 ; MasterYear1 = FileYearN
  else
    print*, "  > ##### ERROR: LoadGrim: Period mismatch #####"
    FileYear0 = MissVal ; FileYear1 = MissVal ; MasterYear0 = MissVal ; MasterYear1 = MissVal
  end if
end if

if (FileYear0.NE.MissVal) then
 if (.not.present(Quiet)) print "(2a)", "   > ", trim(Info)
 DataEn = 0 ; DataMiss = 0
 
 do XBox = 1, BoxN
  read (2,"(a9,i4,a1,i4)") Trash, XExe, Trash, XWye
  
  Grid(XExe,XWye) = XBox
  
  do XFileYear = 1, FileYearN  
    read (2,"(12i5)") (YearData(XMonth), XMonth=1,MonthN)
    
    if (XFileYear.GE.FileYear0.AND.XFileYear.LE.FileYear1) then
     XMasterYear = MasterYear0 + XFileYear - FileYear0    
     do XMonth = 1, MonthN
      DataEn = DataEn + 1
      if (YearData(XMonth).EQ.MissVal) then
        Data(XMasterYear,XMonth,XBox) = MissVal
        DataMiss = DataMiss + 1
      else
      	Data(XMasterYear,XMonth,XBox) = real(YearData(XMonth)) * Multiplier
      end if
     end do
    end if
  end do
 end do

 if (.not.present(Quiet)) print "(a,2i12)", "   > Total and missing values: ", DataEn, DataMiss
else
 print*, "  > ##### Nothing was loaded from the file. #####"
end if

close (2)

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

if (QZip.EQ.2) then
  call system ('compress ' // LoadFile // ' &')
  if (.not.present(Quiet)) print*, "  > Recompressing in the background..."
end if

end subroutine LoadGrim

!*******************************************************************************
! saves GRIdded Seasonic files (12*month,4*season,1*annual)
! this ensures that the data is saved in the order implicit in the grid dimensions:
!   i.e. long (X) varies slowest (-180...180), and lat (Y) quickest (-90...90)

subroutine SaveGrip (Data,Grid,Bounds,Info,CallFile,CallSuffix,SaveSuffix)

real, pointer, dimension (:,:)		:: Data

integer, pointer, dimension (:,:)	:: Grid

real, dimension (4)			:: Bounds

character (len=80), dimension (5)	:: Headers

character (len=80), intent (in)		:: Info, CallFile
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: SaveSuffix

real, parameter :: MissVal = -999.0

real :: Divisor

integer :: SeasonN,BoxN,ExeN,WyeN
integer :: XSeason,XBox,XHeader,XExe,XWye
integer :: AllocStat

character (len=80) :: SaveFile, Variable, SaveBinFile

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

SeasonN = size (Data,1) ; BoxN = size (Data,2) 
ExeN  = size (Grid,1) ; WyeN   = size (Grid,2)

if (SeasonN  .NE. 17) print*, "  > ##### ERROR: SaveGrip: weird calendar #####"

!if (Bounds(1).LT.-180.OR.Bounds(2).GT.180.OR.Bounds(1).GT.Bounds(2)) &
!	print*, "  > ##### ERROR: SaveGrip: longitudes not within spherical geometry #####"
!
!if (Bounds(3).LT.-90.OR.Bounds(4).GT.90.OR.Bounds(3).GT.Bounds(4)) &
!	print*, "  > ##### ERROR: SaveGrip: latitudes  not within spherical geometry #####"

if (Bounds(2).EQ.Bounds(1)) print*, "  > ##### ERROR: SaveGrip: zero longitude range #####"
if (Bounds(4).EQ.Bounds(3)) print*, "  > ##### ERROR: SaveGrip: zero latitude range #####"

call ReviewCall (CallFile,CallSuffix,SaveFile,SaveSuffix,2)	! checks for file/suffix consistency
SaveFile = trim(SaveFile) // ".X"

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

print "(2a)", "   > ", trim(Info)

open  (2,file=SaveFile,status="replace",access="sequential",form="unformatted",action="write")

call WriteHeaders (SaveSuffix,Info,Bounds,ExeN,WyeN,BoxN,0,0,"grip")

do XExe = 1, ExeN
 do XWye = 1, WyeN
  if (Grid(XExe,XWye).NE.MissVal) then   
   write (2) XExe, XWye, (Data(XSeason,Grid(XExe,XWye)), XSeason=1,SeasonN)  
  end if
 end do
end do

close (2)

end subroutine SaveGrip

!*******************************************************************************
! when using this routine, do not assume that this Grid and your Grid are identical
! however, if the same cells are non-missing in the original saves to grim,
!   the SaveGrip procedure ensures that the Grids obtained will be the same from the two loads from grim
! Data (A) may be defined (=Y) or not (=N) in the call
! if Data is defined in the call, it must also have already been initialised if that is necessary

subroutine LoadGrip (Data,Grid,Bounds,Info,CallFile,CallSuffix,LoadSuffix)

real, pointer, dimension (:,:)		:: Data
real, dimension (17)			:: Line

integer, pointer, dimension (:,:)	:: Grid

character (len=80), dimension (5)	:: Headers

real, dimension (4)			:: Bounds

character (len=80), intent (out)	:: Info
character (len=80), intent (in)		:: CallFile
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: LoadSuffix

real, parameter :: MissVal = -999.0

integer, parameter :: SeasonN = 17

integer :: BoxN,ExeN,WyeN
integer :: XSeason,XBox,XHeader,XExe,XWye
integer :: AllocStat
integer :: LoadFileLen,QBin
integer :: Year0,Year1
integer :: DataEn,DataMiss

character (len=80) :: LoadFile, Variable, Trash

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

call ReviewCall (CallFile,CallSuffix,LoadFile,LoadSuffix,1)	! checks for file/suffix consistency

open  (2,file=LoadFile,status="old",access="sequential",form="unformatted",action="read")

call ReadHeaders (LoadSuffix,Info,Bounds,ExeN,WyeN,BoxN,Year0,Year1,"grip")

if (Bounds(2).EQ.Bounds(1)) print*, "  > ##### ERROR: LoadGrip: zero longitude range #####"
if (Bounds(4).EQ.Bounds(3)) print*, "  > ##### ERROR: LoadGrip: zero latitude range #####"

if (associated(Data)) then
  if (size(Data,1).NE.SeasonN)  print*, "  > ##### ERROR: LoadGrip: SeasonN mismatch #####"
  if (size(Data,2).NE.BoxN)     print*, "  > ##### ERROR: LoadGrip: Domain mismatch #####"
else
  allocate (Data (SeasonN,BoxN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGrip: Allocation failure: Data #####"
  Data = MissVal
end if

allocate (Grid (ExeN,WyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGrip: Allocation failure: All #####"
Grid = MissVal

if (Year1.GE.Year0) then
 print "(2a)", "   > ", trim(Info)

 do XBox = 1, BoxN
  read (2) XExe, XWye, (Line(XSeason), XSeason=1,SeasonN)
  
  Grid (XExe,XWye) = XBox
  
  do XSeason = 1, SeasonN
    Data(XSeason,XBox) = Line(XSeason)
  end do
 end do
else
 print*, "  > Nothing was loaded from the file."
end if

close (2)

end subroutine LoadGrip

!*******************************************************************************
! this saves the data from individual month/year files into one big array
! so do work, save to year/month files, and call this
! hence Save Monthly Grim files Post-analysis

subroutine SaveGrimMonthPost (DataFile,BoxN,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix)

real, pointer, dimension (:,:,:)		:: Data

real, pointer, dimension (:)			:: DataMonth

integer, pointer, dimension (:,:)		:: Grid
integer, pointer, dimension (:)			:: YearAD

character (len=80), pointer, dimension (:,:)	:: DataFile
character (len= 4), pointer, dimension (:)	:: YearText
character (len= 2), pointer, dimension (:)	:: MonthText

real, dimension (4)				:: Bounds

integer, intent (in)			:: BoxN

character (len=80), intent (in)		:: Info, CallFile
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: SaveSuffix

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: YearN,MonthN,XYear,XMonth,XBox
integer :: CurrentPos,CurrentCod

character (len= 4) :: GivenText

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

YearN  = size (DataFile,1)
MonthN = size (DataFile,2)

allocate (Data(YearN,MonthN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGrimMonthPost: Allocation failure #####"

do XYear = 1, YearN						! get Data from individual time steps
  do XMonth = 1, MonthN
    call LoadGrimMonth (DataFile(XYear,XMonth),DataMonth)	! ... and dump to file

    do XBox = 1, BoxN
      Data(XYear,XMonth,XBox) = DataMonth(XBox)  
    end do
    
    deallocate (DataMonth,stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGrimMonthPost: Deallocation failure #####"
  end do
end do

call SaveGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix)

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

end subroutine SaveGrimMonthPost

!*******************************************************************************
! this loads the data into one big array
! saves it as individual month/year files

subroutine LoadGrimMonthPrep (DataFile,BoxN,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix)

real, pointer, dimension (:,:,:)		:: Data

real, pointer, dimension (:)			:: DataMonth

integer, pointer, dimension (:,:)		:: Grid
integer, pointer, dimension (:)			:: YearAD

character (len=80), pointer, dimension (:,:)	:: DataFile
character (len= 4), pointer, dimension (:)	:: YearText
character (len= 2), pointer, dimension (:)	:: MonthText

real, dimension (4)				:: Bounds

integer, intent (out)			:: BoxN

character (len=80), intent (in)		:: CallFile
character (len=80), intent (out)	:: Info
character (len= 4), intent (in)		:: CallSuffix
character (len= 4), intent (out)	:: LoadSuffix

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: YearN,MonthN,XYear,XMonth,XBox
integer :: CurrentPos,CurrentCod

character (len= 4) :: GivenText

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

call LoadGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix)

YearN  = size (Data,1)
MonthN = size (Data,2)
BoxN   = size (Data,3)

allocate (DataFile  (YearN,MonthN), &
	  YearText  (YearN),        &
	  MonthText (MonthN),       &
	  DataMonth (BoxN),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGrimMonthPrep: Allocation failure #####"

open   (1,file="/cru/scratch2/f709762/year-scratch.txt",status="replace")  
write  (1,"(i4)") YearAD(1) 
rewind (1)
read   (1,"(a4)") GivenText
close  (1)

YearText(1) = GivenText					! get year AD in text form
do XYear = 2, YearN
  GivenText = YearText (XYear-1)
    
  CurrentPos = 4
  do
      CurrentCod = iachar(GivenText(CurrentPos:CurrentPos))
      
      if (CurrentCod.LE.56) then
        GivenText(CurrentPos:CurrentPos) = achar(CurrentCod+1)
        CurrentPos = 0
      else
        GivenText(CurrentPos:CurrentPos) = achar(48)
        CurrentPos = CurrentPos - 1
      end if
      
      if (CurrentPos.EQ.0) exit
  end do
    
  YearText (XYear) = GivenText
end do
  
MonthText(1) = "01" ; MonthText(2) = "02" ; MonthText(3) = "03" ; MonthText(4) = "04"
MonthText(5) = "05" ; MonthText(6) = "06" ; MonthText(7) = "07" ; MonthText(8) = "08"
MonthText(9) = "09" ; MonthText(10) = "10" ; MonthText(11) = "11" ; MonthText(12) = "12"

do XYear = 1, YearN						! split up Data into individual time steps
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      DataMonth(XBox) = Data(XYear,XMonth,XBox)
    end do
    
    DataFile(XYear,XMonth) = '/cru/scratch2/f709762/' // YearText(XYear) // '.' &
    				// MonthText(XMonth) // LoadSuffix
    
    call SaveGrimMonth (DataFile(XYear,XMonth),DataMonth)	! ... and dump to file
  end do
end do

deallocate (Data,DataMonth,YearText,MonthText,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGrimMonthPrep: Deallocation failure #####"

end subroutine LoadGrimMonthPrep

!*******************************************************************************
! checks the filename and suffix given in call
! returns the correct suffix, and correct filename with the correct suffix

subroutine ReviewCall (CallFile,CallSuffix,CheckFile,CheckSuffix,LoadSave)

character (len=80), intent (in)		:: CallFile
character (len= 4), intent (in)		:: CallSuffix

character (len=80), intent (out) 	:: CheckFile
character (len= 4), intent (out) 	:: CheckSuffix

integer, intent (in)			:: LoadSave		! 1=load,2=save

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

character (len=80) 			:: Variable, WorkFile
character (len= 4)		 	:: FileSuffix

real :: Factor

integer :: SuffixStart,ReadStatus,ZipStart,CallFileLen

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

CheckFile = "" ; CheckSuffix = "" ; Variable = ""

CallFileLen = len_trim(CallFile)				! remove any zip or binary expression
if (CallFileLen.GT.1) then
 if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then
  WorkFile = CallFile(1:(CallFileLen-2)) // "  "
 else
  WorkFile = CallFile
 end if
else
  WorkFile = CallFile
end if

if (CallSuffix.EQ."".OR.CallSuffix.EQ."    ") then		! suffix is not present in call
  if (WorkFile.EQ."") then					! and there is no file in call
    do
      print*, "  > Enter the grim filepath: "			! 	get filepath
      do								
		read (*,*,iostat=ReadStatus) WorkFile
		if (ReadStatus.EQ.0) exit
      end do
    
      CheckSuffix = GetFileSuffix (WorkFile)			!	get suffix from filepath
      
      if (CheckSuffix.NE."") exit
    end do
    								! 	check filepath
    if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) 
    if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix)     
  else								! and there is a file in call
    do
      CheckSuffix = GetFileSuffix (WorkFile)			!	get suffix from filepath
      
      if (CheckSuffix.EQ."") then
        print*, "  > Enter the grim filepath: "			! 	get filepath
        do								
		read (*,*,iostat=ReadStatus) WorkFile
		if (ReadStatus.EQ.0) exit
        end do
      end if
      
      if (CheckSuffix.NE."") exit
    end do
        							! 	check filepath
    if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) 
    if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix)
  end if     
else								! suffix present in call
  call CheckVariSuffix (CallSuffix,Variable,Factor)
  
  if (Variable.EQ."") then					! 	but unrecognised...
    print*, "  > Supplied suffix is invalid. Supply another. "	! 	therefore obtain another
    do
	CheckSuffix = GetFreshSuffix ()
	if (Variable.NE."") exit
    end do
  else								! 	and is recognised
    CheckSuffix = CallSuffix
  end if
    
  if (WorkFile.NE."") then					! gotta suffix, investigate file in call
    FileSuffix = GetFileSuffix (WorkFile)			!	get suffix from file
    
    if (FileSuffix.EQ.CheckSuffix) then				!	suffix good, so check filepath
      if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) 
      if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix) 
    else
      if (FileSuffix.EQ."") then				! 	suffix bad, so get filepath
      	print*, "  > File suffix is not recognised."
      else
      	print*, "  > File suffix is not the same as the supplied suffix."
      end if
      
      if (LoadSave.EQ.1) CheckFile = LoadPath (Blank,CheckSuffix) 
      if (LoadSave.EQ.2) CheckFile = SavePath (Blank,CheckSuffix) 
    end if
  else								! gotta suffix but no file in call
    if (LoadSave.EQ.1) CheckFile = LoadPath (Blank,CheckSuffix) 
    if (LoadSave.EQ.2) CheckFile = SavePath (Blank,CheckSuffix) 
  end if
end if

end subroutine ReviewCall

!*******************************************************************************
! get any legitimate suffix

function GetFreshSuffix (Verbose, NoPrompt)

integer, intent (in), optional		:: Verbose
integer, optional			:: NoPrompt

real		    	:: Factor
integer			:: ReadStatus
character (len=80)	:: Variable
character (len= 4) 	:: GetFreshSuffix, Suffix

if (.not.present(NoPrompt)) print*, "  > Enter the grim file suffix: "

Variable = ""
do
	read (*,*,iostat=ReadStatus) Suffix
	if (Suffix.EQ."") then
			print*, "  > Blank is unacceptable. Try again."
	else
			call CheckVariSuffix (Suffix,Variable,Factor)
			if (Variable.EQ."") print*, "  > Suffix unrecognised. Try again."
	end if
	if (ReadStatus.EQ.0.AND.Variable.NE."") exit
end do
if (present(Verbose)) print "(2a)", "   > Variable: ", trim(Variable)

GetFreshSuffix = Suffix

end function GetFreshSuffix

!*******************************************************************************
! extract suffix from file

function GetFileSuffix (CallFile,Verbose)

real		    			:: Factor
integer, intent (in), optional		:: Verbose
integer					:: ReadStatus, CallFileLen, SuffixStart
character (len=80)			:: Variable, CallFile, WorkFile
character (len= 4) 			:: GetFileSuffix, Suffix

CallFileLen = len_trim(CallFile)				! remove any zip or binary expression
if (CallFileLen.GT.1) then
 if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then
  WorkFile = CallFile(1:(CallFileLen-2)) // "  "
 else
  WorkFile = CallFile
 end if
else
  WorkFile = CallFile
end if

Suffix = "" ; Variable = "" ; GetFileSuffix = "    "
SuffixStart = index (WorkFile,".",.TRUE.)			! identify suffix      

if (SuffixStart.GT.0) then					! check suffix
	Suffix = WorkFile(SuffixStart:(SuffixStart+3))
	call CheckVariSuffix (Suffix,Variable,Factor)
	if (Variable.EQ."") print*, "  > Suffix unrecognised."
	if (Variable.NE."") then
		if (present(Verbose)) print "(2a)", "   > Variable: ", trim(Variable)
		GetFileSuffix = Suffix				! return suffix
	end if
else
	print*, "  > No identifiable suffix in filepath."
end if
        
end function GetFileSuffix

!*******************************************************************************
! relies on Suffix being correct !!

subroutine WriteHeaders (Suffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD0,YearAD1,FileType)

real, dimension (4)				:: Bounds

character (len=80), intent (in)			:: Info
character (len= 4), intent (in)			:: Suffix, FileType

integer, intent(in) :: BoxN,ExeN,WyeN,YearAD0,YearAD1

real :: Multi

integer :: XBound

character (len=80) :: Variable
character (len=12) :: Date,Time
character (len= 4) :: Year
character (len= 2) :: Month,Day,Hour,Minute

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

call CheckVariSuffix (Suffix,Variable,Multi)

call date_and_time (Date, Time)
Year  = Date (1:4)
Month = Date (5:6)
Day   = Date (7:8)
Hour  = Time (1:2)
Minute= Time (3:4)

if      (FileType.EQ."grim") then
  write (2,"(a15,a4,a17,a2,a1,a2,a1,a4,a4,a2,a1,a2,a20)") "Tyndall Centre ", FileType, " file created on ", &
		Day, ".", Month, ".", Year, " at ", Hour, ":", Minute, " by Dr. Tim Mitchell"
  write (2,"(a4,a3,a)") Suffix, " = ", trim(Variable)  
  write (2,"(a)") Info
  write (2,"(a6,f7.2,a1,f7.2,a8,f7.2,a1,f7.2,a12,i4,a1,i4,a1)") &
		"[Long=", Bounds(1), ",", Bounds(2), "] [Lati=", Bounds(3), ",", Bounds(4), &
		"] [Grid X,Y=", ExeN, ",", WyeN, "]"
  write (2,"(a7,i8,a9,i4,a1,i4,a9,f10.4,a16)") &
		"[Boxes=", BoxN, "] [Years=", YearAD0, "-", YearAD1, "] [Multi=", Multi, "] [Missing=-999]"
else if (FileType.EQ."grip") then
  write (2) FileType, Day, Month, Year, Hour, Minute
  write (2) Suffix, Variable
  write (2) Info
  write (2) (Bounds(XBound),XBound=1,4), ExeN, WyeN
  write (2) BoxN, 0, 0, 1.0, -999.0
end if

end subroutine WriteHeaders

!*******************************************************************************
! checks filename suffix against header suffix

subroutine ReadHeaders (FileNameSuffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD0,YearAD1,FileType)

real, dimension (4)			:: Bounds

character (len= 4), intent (in)		:: FileNameSuffix,FileType
character (len=80), intent (out)	:: Info

integer, intent(out) :: BoxN,ExeN,WyeN,YearAD0,YearAD1

real :: LoadMulti, LoadMissVal

integer :: XBound

character (len=80) :: Trash, LoadVariable
character (len= 4) :: HeaderSuffix, LoadFileType
character (len= 4) :: Year
character (len= 2) :: Month,Day,Hour,Minute

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

if      (FileType.EQ."grim") then
  read (2,"(a15,a4)") 	Trash, LoadFileType
  read (2,"(a4)")  	HeaderSuffix
  read (2,"(a80)") 	Info
  read (2,"(a6,f7.2,a1,f7.2,a8,f7.2,a1,f7.2,a12,i4,a1,i4)") &
		Trash, Bounds(1), Trash, Bounds(2), Trash, Bounds(3), Trash, Bounds(4), &
		Trash, ExeN, Trash, WyeN
  read (2,"(a7,i8,a9,i4,a1,i4)") Trash, BoxN, Trash, YearAD0, Trash, YearAD1
else if (FileType.EQ."grip") then
  read (2) LoadFileType, Day, Month, Year, Hour, Minute
!  read (2), HeaderSuffix, LoadVariable
  read (2) HeaderSuffix
  read (2) Info
  read (2) (Bounds(XBound),XBound=1,4), ExeN, WyeN
  read (2) BoxN, YearAD0,YearAD1, LoadMulti, LoadMissVal
end if

if (LoadFileType.EQ."file") LoadFileType = "grim"

if (LoadFileType.NE.FileType) then
  print*, "  > ##### ERROR: ReadHeaders: Mismatch between called and actual file type."
  YearAD1 = YearAD0 - 1
end if

if (HeaderSuffix.NE.FileNameSuffix) then
  print*, "  > ##### ERROR: ReadHeaders: Mismatch between filename and header suffixes."
  YearAD1 = YearAD0 - 1
end if

end subroutine ReadHeaders

!*******************************************************************************
! saves a single month's data in a temporary file

subroutine SaveGrimMonth (CallFile,DataMonth)

real, pointer, dimension (:) 	:: DataMonth
character (len=80), intent (in) :: CallFile
integer :: BoxN, XBox

BoxN = size (DataMonth)
open  (3,file=CallFile,status="new",access="sequential",form="unformatted",action="write")
write (3) BoxN 
do XBox = 1, BoxN
  write (3) DataMonth(XBox)
end do
close (3)

end subroutine SaveGrimMonth

!*******************************************************************************
! loads a single month's data from a temporary file

subroutine LoadGrimMonth (CallFile,DataMonth)

real, pointer, dimension (:) 	:: DataMonth
character (len=80), intent (in) :: CallFile
integer :: BoxN, XBox, AllocStat

open (3,file=CallFile,status="old",access="sequential",form="unformatted",action="read")
read (3) BoxN 
allocate (DataMonth(BoxN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGrimMonth: Allocation failure #####"
do XBox = 1, BoxN
  read (3) DataMonth(XBox)
end do
close (3)

end subroutine LoadGrimMonth

!*******************************************************************************
! BoxN : >0:grids are the same with this number of boxes, MissVal=grids are not the same

subroutine CheckGridAB (GridA,GridB,BoxN)

integer, pointer, dimension (:,:) :: GridA, GridB
integer, intent(out) :: BoxN
real, parameter :: MissVal = -999.0
integer :: ExeN,WyeN,XExe,XWye

BoxN = 0

if (size(GridA,1).NE.size(GridB,1)) then
  BoxN = MissVal
else
  ExeN = size(GridA,1)
end if

if (size(GridA,2).NE.size(GridB,2)) then
  BoxN = MissVal
else
  WyeN = size(GridA,2)
end if

if (BoxN.NE.MissVal) then
 XExe = 0
 do
  XExe = XExe + 1 ; XWye = 0
  
  do
    XWye = XWye + 1
    
    if (GridA(XExe,XWye).EQ.GridB(XExe,XWye)) then
      if (GridA(XExe,XWye).NE.MissVal) BoxN = BoxN + 1
    else
      BoxN = MissVal
    end if
    
    if (XWye.EQ.WyeN) exit
    if (BoxN.EQ.MissVal) exit
  end do
  
  if (XExe.EQ.ExeN) exit
  if (BoxN.EQ.MissVal) exit
 end do
end if

end subroutine CheckGridAB

!*******************************************************************************
! check given file suffix
! Variable is either blank (suffix unrecognised) or contains description of variable
! Factor is zero unless Variable is set

subroutine CheckVariSuffix (Suffix,Variable,Factor)

real,		    intent (out)	:: Factor
character (len=80), intent (out)	:: Variable
character (len= 4), intent (in)		:: Suffix

Variable = "undefined"

if      (Suffix.EQ.".alp") then
  	Variable="precip (mm) gamma distribution alpha"
  	Factor = 0.01
else if (Suffix.EQ.".bet") then
  	Variable="precip (mm) gamma distribution beta"
  	Factor = 0.0001
else if (Suffix.EQ.".cld") then
  	Variable="cloud cover (percentage)"
  	Factor = 0.1
else if (Suffix.EQ.".dew") then
  	Variable="dew point temperature (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".dir") then
  	Variable="wind direction (degrees)"
  	Factor = 0.1
else if (Suffix.EQ.".dpd") then
  	Variable="dew point depression (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".dsf") then
  	Variable="downward shortwave flux (W/m2)"
  	Factor = 0.1
else if (Suffix.EQ.".dtr") then
  	Variable="diurnal temperature range (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".elv") then
  	Variable="elevation (km)"
  	Factor = 0.001
else if (Suffix.EQ.".frs") then
  	Variable="ground frost frequency (days)"
  	Factor = 0.01
else if (Suffix.EQ.".oct") then
  	Variable="cloud cover (octas)"
  	Factor = 0.01
else if (Suffix.EQ.".pre") then
  	Variable="precipitation (mm)"
  	Factor = 0.1
else if (Suffix.EQ.".prc") then
  	Variable="precip change (%)"
  	Factor = 0.001
else if (Suffix.EQ.".prd") then
  	Variable="precip change (%)"
  	Factor = 0.1
else if (Suffix.EQ.".prp") then
  	Variable="precipitation rate (mm/day)"
  	Factor = 0.1
else if (Suffix.EQ.".prr") then
  	Variable="precipitation rate (mm/day)"
  	Factor = 0.001
else if (Suffix.EQ.".rad") then
  	Variable="radiation (W/(m*m))"
  	Factor = 0.1
else if (Suffix.EQ.".rhm") then
  	Variable="relative humidity (%)"
  	Factor = 0.1
else if (Suffix.EQ.".slp") then
  	Variable="sea level pressure"
  	Factor = 0.1
else if (Suffix.EQ.".spc") then
  	Variable="sunshine percentage"
  	Factor = 0.1
else if (Suffix.EQ.".spn") then
  	Variable="sunshine (percentage of normal)"
  	Factor = 0.1
else if (Suffix.EQ.".ssh") then
  	Variable="sunshine hours"
  	Factor = 0.01
else if (Suffix.EQ.".stp") then
  	Variable="station pressure"
  	Factor = 0.1
else if (Suffix.EQ.".tmp") then
  	Variable="near-surface temperature (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".tmk") then
  	Variable="near-surface temperature (degrees Kelvin)"
  	Factor = 0.1
else if (Suffix.EQ.".tmn") then
  	Variable="near-surface temperature minimum (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".tik") then
  	Variable="near-surface temperature minimum (degrees Kelvin)"
  	Factor = 0.1
else if (Suffix.EQ.".tmx") then
  	Variable="near-surface temperature maximum (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".twb") then
  	Variable="wet-bulb temperature (degrees Celsius)"
  	Factor = 0.1
else if (Suffix.EQ.".txk") then
  	Variable="near-surface temperature maximum (degrees Kelvin)"
  	Factor = 0.1
else if (Suffix.EQ.".tcl") then
  	Variable="total cloud (fraction) from long-wave radiation"
  	Factor = 0.0001
else if (Suffix.EQ.".tcs") then
  	Variable="total cloud (fraction) from short-wave radiation"
  	Factor = 0.0001
else if (Suffix.EQ.".vap") then
  	Variable="vapour pressure (hPa)"
  	Factor = 0.1
else if (Suffix.EQ.".wet") then
  	Variable="wet day frequency (days)"
  	Factor = 0.01
else if (Suffix.EQ.".wnd") then
  	Variable="wind (m/s)"
  	Factor = 0.1
else
	Variable = ""
	Factor   = 0.0
end if

end subroutine CheckVariSuffix

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

end module GRIMFiles