! glofiles.f90
! module procedure written by Tim Mitchell on 29.03.01
! last modification on 14.11.01
! contains the relevant .glo file routines:
!    LoadGloVec and GetGloHeaders
!    LoadGloGrid is only if you want to load the grid of values as a raw thing (LongN,LatN)
!    SaveGlo
!    GetGloWeights,GetGloRef

module GloFiles

use FileNames

implicit none

contains

!*******************************************************************************
! if you don't have MapIDLReg before the call, get it
! CallGloFile may be blank

subroutine LoadGloGrid (CallGloFile, FileData, Silent)

integer, intent(in), optional			:: Silent															
character (len=80), intent (in)			:: CallGloFile			! may be blank in call

real, pointer, dimension (:,:) :: FileData

real, parameter :: MissVal = -999.0

real :: TimMin, TimMax, TimTot, GloAv
real :: TotalMiss, PerCentMiss

character (len=80) :: Title, ModelFilePath, GloFile
character (len=10) :: Format

integer :: LongN, LatN, RegN
integer :: XLong, XLat, XEight, XTen, XSix, XReg, XLetter
integer :: Long0, Long1
integer :: AllocStat, ReadStatus
integer :: RegMisMatch
integer :: GloFileLen, QZip

!***************************************
! load .glo data

GloFile = LoadPath (CallGloFile,".glo")					! get .glo file path

GloFileLen = len_trim(GloFile)
if (GloFileLen.GE.2.AND.GloFile((GloFileLen-1):GloFileLen).EQ.'.Z') then
  QZip = 2							! file is zipped
  call system ('uncompress ' // trim(GloFile))
  GloFile ((GloFileLen-1):GloFileLen) = "  "			! change filename to the unzipped file
else
  QZip = 1							! file already unzipped
end if

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

read (2, fmt="(A80)"), Title
read (2, fmt="(A10)"), Format
Format = trim(adjustl(Format))

if      (Format(2:2).EQ."8") then
  XLetter = 3
else if (Format(2:2).EQ."6") then
  XLetter = 3
else if (Format(2:3).EQ."10") then
  XLetter = 4
else
  XLetter = 3
end if
							! allows .glo without format specifier
if (Format(XLetter:XLetter).NE.'E'.AND.Format(XLetter:XLetter).NE.'F') then
    close (2)
    open (2, file=GloFile, status="old", access="sequential", form="formatted", action="read")
    read (2, fmt="(A80)"), Title
    Format = "(8F12.4)"
end if
if (.NOT.present(Silent)) print "(2a)", "   > ", trim(adjustl(Title))

read (2, fmt="(2I6)"), LongN, LatN
read (2, fmt="(A80)"), ModelFilePath

allocate (FileData  (LongN,LatN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGloGrid: Allocation failure #####"
FileData = MissVal

if      (Format(2:2).EQ."8") then
 do XLat = 1, LatN
  do XEight = 1, (LongN/8)
   Long0 = (XEight-1)*8 + 1
   Long1 = Long0 + 7
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (Format(2:3).EQ."10") then
 do XLat = 1, LatN
  do XTen = 1, (LongN/10)
   Long0 = (XTen-1)*10 + 1
   Long1 = Long0 + 9
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (Format(2:2).EQ."6") then
 do XLat = 1, LatN
  do XSix = 1, (LongN/6)
   Long0 = (XSix-1)*6 + 1
   Long1 = Long0 + 5
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
end if

close (2)

where (FileData.LT.0.1E-05.AND.FileData.GT.-0.1E-05) FileData = 0.0

if (QZip.EQ.2) call system ('compress ' // GloFile // ' &')

end subroutine LoadGloGrid

!*******************************************************************************
! if you don't have MapIDLReg before the call, get it
! CallGloFile may be blank

subroutine LoadGloVec (CallGloFile, MapIDLReg, GloVec, Silent)

real, pointer, dimension (:)			:: GloVec			! one value per region
					
integer, pointer, dimension (:,:) 		:: MapIDLReg			! maps IDL grid to regions

integer, intent(in), optional			:: Silent															
character (len=80), intent (in)			:: CallGloFile			! may be blank in call

real, allocatable, dimension (:,:) :: FileData

real, parameter :: MissVal = -999.0

real :: TimMin, TimMax, TimTot, GloAv
real :: TotalMiss, PerCentMiss
real :: TinyVal, HugeVal

character (len=80) :: Title, ModelFilePath, GloFile
character (len=10) :: Format

integer :: LongN, LatN, RegN
integer :: XLong, XLat, XEight, XSix, XReg, XTen, XLetter
integer :: Long0, Long1
integer :: AllocStat, ReadStatus
integer :: RegMisMatch
integer :: GloFileLen, QZip

!***************************************
! load .glo data

RegN = maxval (MapIDLReg)

GloFile = LoadPath (CallGloFile,".glo")					! get .glo file path

GloFileLen = len_trim(GloFile)
if (GloFileLen.GE.2.AND.GloFile((GloFileLen-1):GloFileLen).EQ.'.Z') then
  QZip = 2							! file is zipped
  call system ('uncompress ' // trim(GloFile))
  GloFile ((GloFileLen-1):GloFileLen) = "  "			! change filename to the unzipped file
else
  QZip = 1							! file already unzipped
end if

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

read (2, fmt="(A80)"), Title
read (2, fmt="(A10)"), Format
Format = trim(adjustl(Format))

if      (Format(2:2).EQ."8") then
  XLetter = 3
else if (Format(2:2).EQ."6") then
  XLetter = 3
else if (Format(2:3).EQ."10") then
  XLetter = 4
else
  XLetter = 3
end if
							! allows .glo without format specifier
if (Format(XLetter:XLetter).NE.'E'.AND.Format(XLetter:XLetter).NE.'F') then
    close (2)
    open (2, file=GloFile, status="old", access="sequential", form="formatted", action="read")
    read (2, fmt="(A80)"), Title
    Format = "(8F12.4)"
end if
if (.NOT.present(Silent)) print "(2a)", "   > ", trim(adjustl(Title))

read (2, fmt="(2I6)"), LongN, LatN
read (2, fmt="(A80)"), ModelFilePath

write (99,"(2a)"), "got headers", Format	! @@@@@@@@@@@@@@@@@@@@@@@

allocate (FileData  (LongN,LatN), &
	  GloVec    (RegN),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadGlo: Allocation failure #####"
FileData = MissVal ; GloVec = MissVal

if      (Format(2:2).EQ."8") then
 do XLat = 1, LatN
  do XEight = 1, (LongN/8)
   Long0 = (XEight-1)*8 + 1
   Long1 = Long0 + 7
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (Format(2:3).EQ."10") then
 do XLat = 1, LatN
  do XTen = 1, (LongN/10)
   Long0 = (XTen-1)*10 + 1
   Long1 = Long0 + 9
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (Format(2:2).EQ."6") then
 do XLat = 1, LatN
  do XSix = 1, (LongN/6)
   Long0 = (XSix-1)*6 + 1
   Long1 = Long0 + 5
  
   read (2, fmt=Format), (FileData (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
end if

close (2)

write (99,"(2a)"), "got data", Format	! @@@@@@@@@@@@@@@@@@@@@@@

if (QZip.EQ.2) call system ('compress ' // GloFile // ' &')

!***************************************
! process .glo data

RegMisMatch = 0

TinyVal = tiny(GloVec) ; HugeVal = huge(GloVec)

do XLong = 1, LongN
  do XLat = 1, LatN
   if (FileData(XLong,XLat).LT.TinyVal.AND.FileData(XLong,XLat).GT.(0.0-TinyVal)) &
   			FileData(XLong,XLat) = TinyVal
   
   if (FileData(XLong,XLat).GT.HugeVal) FileData(XLong,XLat) = HugeVal
   
   if (FileData(XLong,XLat).LT.(0.0-HugeVal)) FileData(XLong,XLat) = 0.0 - HugeVal
   
   				! grid box relates to a region
   if (MapIDLReg(XLong,XLat).NE.MissVal) then

     				! ...and region already has a value
     if (GloVec(MapIDLReg(XLong,XLat)).NE.MissVal) then
       				! ......but the region val does not equal the grid box val
       if (GloVec(MapIDLReg(XLong,XLat)).NE.FileData(XLong,XLat)) RegMisMatch = RegMisMatch + 1
       
     else			! ...but region does not (yet) have a value
       				! ...so we give it a value
       GloVec (MapIDLReg(XLong,XLat)) = FileData (XLong,XLat)
     end if

   else				! grid box does not relate to a region
       				! ...yet it contains a value
     if (FileData(XLong,XLat).NE.MissVal) RegMisMatch = RegMisMatch + 1

   end if
   
  end do
end do

if (RegMisMatch.GT.0) then
	print*, "  > ##### ERROR: LoadGlo: .glo Belongs to Wrong Region Set #####"
	print*, "  > Number of regions in call: ", RegN
	GloVec = MissVal
end if

deallocate (FileData)

if (.NOT.present(Silent)) then
  TimMin    =  100000.0
  TimMax    = -100000.0
  TimTot    =       0.0
  TotalMiss =       0.0
  do XReg = 1, RegN
    if      (GloVec (XReg).EQ.MissVal) then
       TotalMiss = TotalMiss + 1.0
    else
       TimTot = TimTot + GloVec (XReg)
       
       if (GloVec (XReg).GT.TimMax)  then
         TimMax    = GloVec (XReg)
       end if
       
       if (GloVec (XReg).LT.TimMin)  then
         TimMin    = GloVec (XReg)
       end if
    end if
  end do

  PerCentMiss = 100.0 * TotalMiss / (RegN)
  GloAv       = MissVal

  if (TotalMiss.LT.RegN) GloAv = TimTot / (RegN - TotalMiss)
  print "(a28,4f10.4)", "   > % miss, min, max, mean: ", PerCentMiss, TimMin, TimMax, GloAv
end if

end subroutine LoadGloVec

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

subroutine GetGloHeaders (CallGloFile,Title,Format,LongN,LatN,ModelFilePath)

character (len=80), intent (in)		:: CallGloFile			! may be blank in call
character (len=80), intent (out) 	:: Title, ModelFilePath
character (len=10), intent (out) 	:: Format

integer, intent (out) 			:: LongN, LatN

integer :: GloFileLen, QZip, XLetter

character (len=80) :: GloFile

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

GloFile = LoadPath (CallGloFile,".glo")					! get .glo file path

GloFileLen = len_trim(GloFile)
if (GloFileLen.GE.2.AND.GloFile((GloFileLen-1):GloFileLen).EQ.'.Z') then
  QZip = 2							! file is zipped
  call system ('uncompress ' // trim(GloFile))
  GloFile ((GloFileLen-1):GloFileLen) = "  "			! change filename to the unzipped file
else
  QZip = 1							! file already unzipped
end if

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

read (2, fmt="(A80)"), Title
read (2, fmt="(A10)"), Format
Format = trim(adjustl(Format))

if      (Format(2:2).EQ."8") then
  XLetter = 3
else if (Format(2:2).EQ."6") then
  XLetter = 3
else if (Format(2:3).EQ."10") then
  XLetter = 4
else
  XLetter = 3
end if
							! allows .glo without format specifier
if (Format(XLetter:XLetter).NE.'E'.AND.Format(XLetter:XLetter).NE.'F') then
    close (2)
    open (2, file=GloFile, status="old", access="sequential", form="formatted", action="read")
    read (2, fmt="(A80)"), Title
    Format = "(8F12.4)"
end if

read (2, fmt="(2I6)"), LongN, LatN
read (2, fmt="(A80)"), ModelFilePath

close (2)

if (QZip.EQ.2) call system ('compress ' // GloFile // ' &')

end subroutine GetGloHeaders

!*******************************************************************************
! altered on 12.07.00 to prompt for title; SpecTitle and SpecFilePath (len=80) can both be set to ""
! altered on 04.10.01 to allow for save as HadCM2 when calling with 97,73 grid

subroutine SaveGlo (CallLongN,CallLatN,CallRegN,CallFilePath,SpecFilePath,SpecTitle,GloAnaSlice,CallMapIDLReg)

integer, pointer, dimension (:,:) 		:: CallMapIDLReg, MapIDLReg
					! maps IDL grid to regions
real, pointer, dimension (:)			:: GloAnaSlice
					! one value per region

character (len=80), intent (in)			:: CallFilePath, SpecFilePath, SpecTitle
integer, intent (in)				:: CallLongN, CallLatN, CallRegN

real, allocatable, dimension (:,:) :: OutArray

real, parameter :: MissVal = -999.0

character (len=80) :: GivenFile, GloFilePath, GloTitle, FilePath
character (len=10) :: GloFormat

integer :: LongN, LatN
integer :: XLong, XLat, XEight, XTen, XSix
integer :: Long0, Long1, NoPolesBeg
integer :: AllocStat			! status of allocation statement
integer :: ReadStatus			! status of user input

if (SpecTitle.EQ."") then
  print*, "  > Enter the .glo file title: "
  do
	read (*,*,iostat=ReadStatus), GloTitle
	if (ReadStatus.LE.0.AND.GloTitle.NE."") exit
  end do
else
  GloTitle = SpecTitle
end if

if (SpecFilePath.EQ."") then
  print*, "  > Enter the filepath of the .glo file to save: "
  do
	do
		read (*,*,iostat=ReadStatus), GivenFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=GivenFile, name=GloFilePath)
	open (1, file=GloFilePath, status="new", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
  end do
else
  GloFilePath = SpecFilePath
end if

if (CallLongN.EQ.97.AND.CallLatN.EQ.73) then		! this goes from h2grid to HadCM2
  LatN=73 ; LongN=96
  
  NoPolesBeg = index(CallFilePath,'h2grid.ref')
  FilePath = CallFilePath(1:(NoPolesBeg-1)) // 'hadcm2.ref'
  
  allocate (MapIDLReg(96,73)) ; MapIDLReg = MissVal
  do XLong = 1, LongN
    do XLat = 2, 72
      MapIDLReg(XLong,XLat) = CallMapIDLReg(XLong,XLat)
    end do
  end do
else
  LatN=CallLatN ; LongN=CallLongN ; FilePath=CallFilePath ; MapIDLReg=>CallMapIDLReg
end if

allocate (OutArray (LongN,LatN))

do XLong = 1, LongN
  do XLat = 1, LatN
    if (MapIDLReg(XLong,XLat).NE.MissVal) then
      OutArray (XLong,XLat) = GloAnaSlice (MapIDLReg (XLong,XLat))
    else
      OutArray (XLong,XLat) = MissVal
    end if
  end do
end do 

if ((real(LongN)/8.0).EQ.nint(real(LongN)/8.0)) then
  GloFormat = "(8E12.4)"			! examine LoadGlo before altering
else if ((real(LongN)/10.0).EQ.nint(real(LongN)/10.0)) then
  GloFormat = "(10E12.4)"			! examine LoadGlo before altering
else if ((real(LongN)/6.0).EQ.nint(real(LongN)/6.0)) then
  GloFormat = "(6E12.4)"			! examine LoadGlo before altering
else
  print*, "  > ##### ERROR: SaveGlo: LongN not divisible by 10 or 8 or 6 #####"
end if

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

write (2, fmt="(A80)"), GloTitle
write (2, fmt="(A10)"), GloFormat
write (2, fmt="(2I6)"), LongN, LatN
write (2, fmt="(A80)"), FilePath

if (GloFormat.EQ."(8E12.4)") then
 do XLat = 1, LatN
  do XEight = 1, (LongN/8)
   Long0 = (XEight-1)*8 + 1
   Long1 = Long0 + 7
  
   write (2, fmt=GloFormat), (OutArray (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (GloFormat.EQ."(10E12.4)") then
 do XLat = 1, LatN
  do XTen = 1, (LongN/10)
   Long0 = (XTen-1)*10 + 1
   Long1 = Long0 + 9
  
   write (2, fmt=GloFormat), (OutArray (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else if (GloFormat.EQ."(6E12.4)") then
 do XLat = 1, LatN
  do XSix = 1, (LongN/6)
   Long0 = (XSix-1)*6 + 1
   Long1 = Long0 + 5
  
   write (2, fmt=GloFormat), (OutArray (XLong, XLat), XLong=Long0,Long1 )
  end do
 end do
else
  print*, "  > ##### ERROR: SaveGlo: No save of data #####"
end if

close (2)

deallocate (OutArray)

if (FilePath.NE.CallFilePath) then
  deallocate (MapIDLReg)
else
  nullify (MapIDLReg)
end if

end subroutine SaveGlo

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

subroutine GetGloWeights (ExeN,WyeN,GloWeights)

real, pointer, dimension (:,:) 	:: GloWeights, WeightsXtra

integer, intent(in) 		:: ExeN, WyeN

real, parameter :: MissVal = -999.0

integer :: XWye, XExe
integer :: FileSubBeg, AllocStat

character (len=80) :: GenericFile

write (99,*), "running getgloweights..."					! ########################
GenericFile  = GetGloRef (ExeN,WyeN)
FileSubBeg   = index(GenericFile,".ref")
GenericFile  = GenericFile(1:(FileSubBeg-1)) // ".weights.glo"
  
write (99,*), "calling loadglogrid..."					! ########################
call LoadGloGrid (GenericFile, GloWeights, Silent=1)

write (99,*), "checking gloweights..."					! ########################
if (size(GloWeights,1).NE.ExeN.OR.size(GloWeights,2).NE.WyeN) then
    if (ExeN.EQ.97.AND.size(GloWeights,1).EQ.96) then
      allocate (WeightsXtra(ExeN,WyeN), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: GetGloWeights: Allocation failure: WeightsXtra #####"
      
      do XWye = 1, WyeN
        do XExe = 1, 96
          WeightsXtra(XExe,XWye) = GloWeights(XExe,XWye)
        end do
        WeightsXtra(97,XWye) = GloWeights(1,XWye)
      end do
      
      deallocate (GloWeights, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: GetGloWeights: Deallocation failure: GloWeights #####"
      
      allocate (GloWeights(ExeN,WyeN), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: GetGloWeights: Allocation failure: GloWeights #####"
      
      do XWye = 1, WyeN
        do XExe = 1, ExeN
          GloWeights(XExe,XWye) = WeightsXtra(XExe,XWye)
        end do
      end do
      
      deallocate (WeightsXtra, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: GetGloWeights: Deallocation failure: WeightsXtra #####"      
    else
      GloWeights = MissVal
      print*, "  > ##### ERROR: GetGloWeights: Loaded .glo has wrong dims #####"
    end if
end if
write (99,*), "finished getgloweights"					! ########################

end subroutine GetGloWeights

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

function GetGloRef (ExeN,WyeN)

integer, intent (in) :: ExeN, WyeN

integer :: ReadStatus

character (len=80) :: GivenFile, GetGloRef

if      (ExeN.EQ.720.AND.WyeN.EQ.360) then
  GetGloRef = "/tyn1/tim/constants/goglo/half-degree.ref"
else if (ExeN.EQ.160.AND.WyeN.EQ.106) then
  GetGloRef = "/tyn1/tim/constants/goglo/half-degree-euro.ref"
else if (ExeN.EQ.480.AND.WyeN.EQ.313) then
  GetGloRef = "/tyn1/tim/constants/goglo/tenmin-euro.ref"
else if (ExeN.EQ.258.AND.WyeN.EQ.228) then
  GetGloRef = "/tyn1/tim/constants/goglo/tenmin-ateam.ref"
else if (ExeN.EQ.192.AND.WyeN.EQ.144) then
  GetGloRef = "/tyn1/tim/constants/goglo/regridx2.ref"
else if (ExeN.EQ.180.AND.WyeN.EQ.290) then
  GetGloRef = "/tyn1/tim/constants/goglo/ng5.ref"
else if (ExeN.EQ.144.AND.WyeN.EQ. 73) then
  GetGloRef = "/tyn1/tim/constants/goglo/richjones.ref"
else if (ExeN.EQ.130.AND.WyeN.EQ.242) then
  GetGloRef = "/tyn1/tim/constants/goglo/ng5cut.ref"
else if (ExeN.EQ. 96.AND.WyeN.EQ. 73) then
  GetGloRef = "/tyn1/tim/constants/goglo/hadcm2.ref"  
else if (ExeN.EQ. 97.AND.WyeN.EQ. 73) then
  GetGloRef = "/tyn1/tim/constants/goglo/h2grid.ref"  
else if (ExeN.EQ. 80.AND.WyeN.EQ. 72) then
  GetGloRef = "/tyn1/tim/constants/goglo/tenmin-brit.ref"
else if (ExeN.EQ. 64.AND.WyeN.EQ. 56) then
  GetGloRef = "/tyn1/tim/constants/goglo/csiro.ref"
else if (ExeN.EQ. 96.AND.WyeN.EQ. 48) then
  GetGloRef = "/tyn1/tim/constants/goglo/cgcm1.ref"
else if (ExeN.EQ.128.AND.WyeN.EQ. 64) then
  GetGloRef = "/tyn1/tim/constants/goglo/ncar-pcm3.ref"
else
  print*, "  > Enter the grid .ref file (/tyn1/tim/constants/goglo/): "
  do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0) exit
  end do
  GetGloRef = LoadPath (GivenFile,".ref")
end if 

end function GetGloRef

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

end module GloFiles