! gridtogrim.f90
! f90 main program written on 23.03.01 by Tim Mitchell
! last modification on 17.12.01
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../grim/gridtogrim time.f90 filenames.f90 grimfiles.f90 
!	./../grim/gridtogrim.f90 2> /tyn1/tim/scratch/stderr.txt
! loads data file(s) into standard grim files  

program GridToGrim

use Time
use FileNames
use GrimFiles

implicit none

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

integer, pointer, dimension (:) 		:: RowInt, YearAD, ReGrid
integer, pointer, dimension (:,:) 		:: Grid, RefGrid

character (len=80), pointer, dimension (:,:) 	:: RawFile
character (len=80), pointer, dimension (:) 	:: OutFile
character (len= 3), pointer, dimension (:) 	:: MonthsNames

real, dimension (4) 				:: Bounds

character (len=3), dimension (12),target :: MonthsFirstCap = &
	(/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'/)
character (len=3), dimension (12),target :: MonthsNoCap    = &
	(/'jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec'/)
character (len=2), dimension (12) :: MonthsNumeral  = &
	(/'01','02','03','04','05','06','07','08','09','10','11','12'/)

real, parameter 		:: MissVal = -999.0

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

real :: FileMissVal,Multiplier,Fraction

integer :: LongN,LatN,DataN, RegN, FileN,RowN,ColN,CellN, BlockN, YearN,MonthN, RawFileN, RefBoxN
integer :: FileHeadN,MonthHeadN,YearHeadN
integer :: AllocStat,ReadStatus
integer :: DateCol
integer :: XRow,XCol,XCell, XReg, XLong,XLat, XFile, XBlock, XFileHead,XMonthHead,XYearHead, XYear,XMonth
integer :: XRawFile,XBound
integer :: StringLen, PosChar, Col0Cell
integer :: Cell0,Cell1, YearAD0,YearAD1
integer :: QIntReal,QDateGreen,QZip,QDataGrim,QNorthSouth
integer :: SubLen, SubBeg, SuffixBeg,SuffixEnd,SuffixLen, GrimSubBeg, YearSubLen,YearSubBeg,YearSubEnd
integer :: MonthSubLen,MonthSubBeg,MonthSubEnd
integer :: LoadFileLen,CurrentPos,CurrentCod,TotMiss

character (len=1)  :: CharIntReal
character (len=4)  :: SaveSuffix, RefSuffix
character (len=80) :: OrigSub, ThisSub, YearSub, SpecSub
character (len=40) :: LineFormat
character (len=80) :: GridFilePath,GivenFile,FilePath,AutoLoadPath,AutoSavePath,RefFile,GrimFile,LoadFile
character (len=80) :: Info, InfoItem
character (len=80) :: Trash
character (len=20) :: TextFound

logical, parameter :: True = .TRUE.

!*******************************************************************************
! main program

call Initialise

if (QDataGrim.EQ.1) then
  call DomainFromData
else
  call DomainFromGrim
end if

call GetInfoLine
call GetDimensions
call DataSpec
call FirstExec
if (FileN.GT.1) call AutoSpec

call LoadData

call Finalise

contains

!*******************************************************************************
! initialise
! in this program we take raw data files organised as sequential spatial grids
!   (one grid per month, with 12 months (Jan,Feb,...,Dec) * YearN per file, or one month per file)
! we store the data in standard grim files

subroutine Initialise

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

print*
print*, "  > ***** GridToGrim : store gridded monthly data *****"
print*

print*, "  > Construct domain from data (=1) or from a grim file (=2) ? "
do
	read (*,*,iostat=ReadStatus), QDataGrim
	if (ReadStatus.LE.0.AND.QDataGrim.GE.1.AND.QDataGrim.LE.2) exit
end do

end subroutine Initialise

!*******************************************************************************
! specify reference domain = RefGrid (LongN,LatN), with RefBoxN valid cells

subroutine DomainFromData

print*, "  > Select the number of longitude, latitude cells: "
do
	read (*,*,iostat=ReadStatus), LongN, LatN
	if (ReadStatus.LE.0.AND.LongN.GT.0.AND.LatN.GT.0) exit
end do

print*, "  > Select the bounds of the grid (not cell centres). "
print*, "  > Enter min,max longitude, min,max latitude: "
do
	read (*,*,iostat=ReadStatus), Bounds(1),Bounds(2),Bounds(3),Bounds(4)
	if      (Bounds(1).LT.-182.OR.Bounds(2).GT.182.OR.Bounds(2).LE.Bounds(1)) then
	  print*, "  > Unacceptable longitudes. Try again."
	  ReadStatus = 1
	else if (Bounds(3).LT. -92.OR.Bounds(4).GT. 92.OR.Bounds(4).LE.Bounds(3)) then
	  print*, "  > Unacceptable latitudes. Try again."
	  ReadStatus = 1	
	end if
	
	if (ReadStatus.LE.0) exit
end do

DataN   = LongN * LatN
RefBoxN = DataN
Info    = ""

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

XCell = 0
do XLong = 1, LongN
  do XLat = 1, LatN
    XCell = XCell + 1
    RefGrid(XLong,XLat) = XCell
  end do
end do

end subroutine DomainFromData

!*******************************************************************************
! specify reference domain = RefGrid (LongN,LatN), with RefBoxN valid cells

subroutine DomainFromGrim

print*, "  > Enter the grim filepath: "
do
	read (*,*,iostat=ReadStatus), RefFile
	if (ReadStatus.LE.0.AND.RefFile.NE."") exit
end do

call LoadGrim (Data,RefGrid,YearAD,Bounds,Info,RefFile,"    ",RefSuffix)

LongN = size (RefGrid,1) ; LatN = size (RefGrid,2) ; DataN = LongN * LatN
RefBoxN = size (Data,3)
Info = ""

print "(a,3i8)",   "   > Grid dimensions and domain size: ", LongN,LatN,RefBoxN
print "(a,4f8.2)", "   > Grid bounds: ", (Bounds(XBound),XBound=1,4)

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

end subroutine DomainFromGrim

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

subroutine GetInfoLine

print*, "  > Name the grid: "
do
	read (*,*,iostat=ReadStatus), InfoItem
	if (ReadStatus.LE.0.AND.InfoItem.NE."") exit
end do
Info = trim(InfoItem)

print*, "  > Name the domain: "
do
	read (*,*,iostat=ReadStatus), InfoItem
	if (ReadStatus.LE.0.AND.InfoItem.NE."") exit
end do
Info = trim(Info) // " " // trim(InfoItem)

print*, "  > Name the scenario: "
do
	read (*,*,iostat=ReadStatus), InfoItem
	if (ReadStatus.LE.0.AND.InfoItem.NE."") exit
end do
Info = trim(Info) // " " // trim(InfoItem)

end subroutine GetInfoLine

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

subroutine GetDimensions

print*, "  > Select the first, last years to upload: "
do
	read (*,*,iostat=ReadStatus), YearAD0, YearAD1
	if (ReadStatus.LE.0.AND.YearAD1.GE.YearAD0) exit
end do

YearN  = YearAD1 - YearAD0 + 1
MonthN = 12

print*, "  > Select the number of upload executions: "
do
	read (*,*,iostat=ReadStatus), FileN
	if (ReadStatus.LE.0) exit
end do

! allocate Grid  (LongN,LatN), &
! Grid   = MissVal	  

allocate (ReGrid(DataN),                &
	  YearAD(YearN),	        stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetDimensions: Allocation failure #####"
ReGrid = MissVal

do XYear = 1, YearN
  YearAD(XYear) = YearAD0 + XYear - 1
end do

end subroutine GetDimensions

!*******************************************************************************
! specify first execution

subroutine FirstExec

if (FileN.GT.1) print*, "  > Specify the first execution. "

print "(a,i4,a4,i4,a2)", "   > Enter the number of raw files (1 or ", YearN, " or ", (YearN*12), "):"
do
	read (*,*,iostat=ReadStatus), RawFileN
	if (ReadStatus.LE.0.AND.RawFileN.NE.1.AND.RawFileN.NE.YearN.AND.RawFileN.NE.(YearN*12)) then
	  print*, "  > That number is not allowed. Try again."
	  ReadStatus = 1
	end if
	if (ReadStatus.LE.0) exit
end do

allocate (RawFile(FileN,RawFileN),&
          OutFile(FileN),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AutoSpec: Allocation failure #####"
RawFile = "missing" ; OutFile = "missing"

if (RawFileN.EQ.1) print*, "  > Enter the filepath of the raw data file: "
if (RawFileN.GT.1) print*, "  > Enter the filepath of the first raw data file: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
RawFile(1,1) = LoadPath (GivenFile,"    ")
  	  
if (RawFileN.EQ.YearN) then
  TextFound = GetTextFromInt (YearAD0)
  YearSub   = trim(adjustl(TextFound))
  
  YearSubLen = len_trim(YearSub)
  YearSubBeg = index(RawFile(1,1),YearSub(1:YearSubLen))
  YearSubEnd = YearSubBeg + YearSubLen - 1
  
  if (YearSubBeg.GT.0) then
   do XRawFile = 2, RawFileN
    GivenFile = RawFile(1,1)    
    TextFound = GetTextFromInt (YearAD0+XRawFile-1)    
    GivenFile (YearSubBeg:YearSubEnd) = trim(adjustl(TextFound))    
    RawFile(1,XRawFile) = GivenFile
   end do
  end if
else if (RawFileN.EQ.(YearN*12)) then
  TextFound = GetTextFromInt (YearAD0)
  YearSub   = trim(adjustl(TextFound))
  
  YearSubLen = len_trim(YearSub)
  YearSubBeg = index(RawFile(1,1),YearSub(1:YearSubLen))
  YearSubEnd = YearSubBeg + YearSubLen - 1
  
  MonthSubBeg = index(RawFile(1,1),'jan')
  if (MonthSubBeg.EQ.0) then
    MonthSubBeg = index(RawFile(1,1),'Jan')
    if (MonthSubBeg.EQ.0) then
      MonthSubBeg = index(RawFile(1,1),'01')
      if (MonthSubBeg.EQ.0) then
        MonthSubBeg = 0  
      else
        MonthSubLen = 2
      end if    
    else
      MonthSubLen = 3
      MonthsNames => MonthsFirstCap
    end if  
  else
    MonthSubLen = 3
    MonthsNames => MonthsNoCap
  end if
  MonthSubEnd = MonthSubBeg + MonthSubLen - 1
  
  if (YearSubBeg.GT.0.AND.MonthSubBeg.GT.0) then
   do XYear = 1, YearN
    do XMonth = 1, MonthN
     GivenFile = RawFile(1,1)    
    
     TextFound = GetTextFromInt (YearAD0+XYear-1)    
     GivenFile (YearSubBeg:YearSubEnd) = trim(adjustl(TextFound))    
    
     if      (MonthSubLen.EQ.2) then
      GivenFile (MonthSubBeg:MonthSubEnd) = MonthsNumeral (XMonth)
     else if (MonthSubLen.EQ.3) then
      GivenFile (MonthSubBeg:MonthSubEnd) = MonthsNames (XMonth)
     end if
     
     XRawFile = ((XYear-1)*12) + XMonth
     RawFile(1,XRawFile) = GivenFile
    end do
   end do
  end if
end if
  
print*, "  > Enter the filepath of the grim file to save: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
	
call ReviewCall (GivenFile,"    ",AutoSavePath,SaveSuffix,2)	! checks for file/suffix consistency
OutFile(1) = AutoSavePath

end subroutine FirstExec
    
!*******************************************************************************
! specify automatics

subroutine AutoSpec

print*, "  > Enter the substring to vary between executions: "
do
		read (*,*,iostat=ReadStatus), OrigSub
		if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
		else if (OrigSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
		end if
		if (ReadStatus.LE.0.AND.OrigSub.NE."") exit
end do
	
SubLen = len(trim(OrigSub))
  
print*, "  > Enter the substring in each execution, starting with no.2:"
do XFile = 2, FileN
	do
		read (*,*,iostat=ReadStatus), SpecSub
		if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
		else if (SpecSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
			ReadStatus = 1
		else if (len_trim(SpecSub).NE.SubLen) then
			print*, "  > The substring has a wrong length. Try again."
			ReadStatus = 1
		end if
		if (ReadStatus.LE.0) exit
	end do
        
        do XRawFile = 1, RawFileN
          GivenFile = RawFile(1,XRawFile)
          SubBeg = index(GivenFile,OrigSub(1:SubLen))			! start of orig sub
          GivenFile (SubBeg:(SubBeg+SubLen-1)) = trim(SpecSub)	        ! replace orig sub
          RawFile(XFile,XRawFile) = GivenFile
        end do
        
        GivenFile = OutFile(1)
        SubBeg = index(GivenFile,OrigSub(1:SubLen))			! start of orig sub
        GivenFile (SubBeg:(SubBeg+SubLen-1)) = trim(SpecSub)	        ! replace orig sub
        OutFile(XFile) = GivenFile        
end do
  
end subroutine AutoSpec

!*******************************************************************************
! specify data structure in original files

subroutine DataSpec

print*, "  > Enter the number of headers per file, year, month: "
do
	read (*,*,iostat=ReadStatus), FileHeadN,YearHeadN,MonthHeadN
	if (ReadStatus.LE.0.AND.FileHeadN.GE.0.AND.YearHeadN.GE.0.AND.MonthHeadN.GE.0) exit
end do
  
print*, "  > Enter the number of columns, rows in a monthly grid: "
do
	read (*,*,iostat=ReadStatus), ColN, RowN
	if (ReadStatus.LE.0.AND.RowN.GE.1.AND.ColN.GE.1) exit
end do
  
allocate (RowInt (ColN), &
	  RowReal(ColN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DataSpec: Allocation failure #####"
RowInt  = MissVal
RowReal = MissVal

CellN = RowN * ColN

if (CellN.EQ.DataN) then
  if      (RowN.EQ.LatN.AND.ColN.EQ.LongN) then
    print*, "  > Define the raw data structure."
    print*, "  > Top-bottom: N-S (=1) or S-N (=2) ? "
    do
	read (*,*,iostat=ReadStatus), QNorthSouth
	if (ReadStatus.LE.0.AND.QNorthSouth.GE.1.AND.QNorthSouth.LE.2) exit
    end do

    print*, "  > Left-right: W-E"
    print "(a,i4,a)", "   > Enter the long cell (1-", ColN, ") for the first column in a monthly grid: "
    do
	read (*,*,iostat=ReadStatus), Col0Cell
	if (ReadStatus.LE.0.AND.Col0Cell.GE.1.AND.Col0Cell.LE.ColN) exit
    end do
    
    if (QNorthSouth.EQ.1) XLat  = LatN + 1
    if (QNorthSouth.EQ.2) XLat  = 0
    XCell = 0
    
    do XRow = 1, RowN
      if (QNorthSouth.EQ.1) XLat = XLat - 1
      if (QNorthSouth.EQ.2) XLat = XLat + 1
      
      do XLong = Col0Cell, ColN
        XCell = XCell + 1
!        Grid   (XLong,XLat) = XCell
        ReGrid (XCell)      = RefGrid (XLong,XLat)
      end do
      
      if (Col0Cell.GT.1) then
       do XLong = 1, (Col0Cell-1)
        XCell = XCell + 1
!        Grid   (XLong,XLat) = XCell
        ReGrid (XCell)      = RefGrid (XLong,XLat)
       end do      
      end if
    end do
    
    BlockN = 1
  else if (ColN.EQ.LatN.AND.RowN.EQ.LongN) then
    print*, "  > ##### ERROR: rotated section X not written #####"
  else
    							! The grid is actually split by long in file
    if (Fraction.EQ.floor(Fraction)) then		!    but we imagine it as a single grid    
      print*, "  > ##### ERROR: multiple blocks per month not written #####"
    else
      print*, "  > ##### ERROR: difficult blocks not written #####"
    end if
  end if
else
  print*, "  > ##### ERROR: The cells in file do not match with grid. #####"
end if

print*, "  > Enter the data line format, as multiples of i,f,e: "
do
	read (*,*,iostat=ReadStatus), LineFormat
	
	StringLen = len (trim(LineFormat))	
        
	PosChar = max (scan(LineFormat,'I'),scan(LineFormat,'i'),scan(LineFormat,'F'),scan(LineFormat,'f'))
	PosChar = max (scan(LineFormat,'E'),scan(LineFormat,'e'),PosChar)
	
	CharIntReal = ""
	if (PosChar.GT.0) CharIntReal = LineFormat(PosChar:PosChar)
	
	if (LineFormat(1:1).NE."(" .OR. LineFormat(StringLen:StringLen).NE.")") then
	    ReadStatus = 99
	    print*, "  > Format is unacceptable. Retry."
	else 
	  if      (CharIntReal.EQ.'I'.OR.CharIntReal.EQ.'i') then
	    QIntReal = 1
	  else if (CharIntReal.EQ.'F'.OR.CharIntReal.EQ.'f') then
	    QIntReal = 2
	  else if (CharIntReal.EQ.'F'.OR.CharIntReal.EQ.'f') then
	    QIntReal = 2
	  else
	    ReadStatus = 99
	    print*, "  > Format is unacceptable. Retry."
	  end if
	end if
	
	if (ReadStatus.LE.0) exit
end do
    
print*, "  > Enter the multiplier by which to convert data: "
do
	read (*,*,iostat=ReadStatus), Multiplier
	if (ReadStatus.LE.0) exit
end do
  
print*, "  > Enter the file missing value: "
do
	read (*,*,iostat=ReadStatus), FileMissVal
	if (ReadStatus.LE.0) exit
end do
  
end subroutine DataSpec

!*******************************************************************************
! load data

subroutine LoadData

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

do XFile = 1, FileN
 XRawFile = 0
 
 if (RawFileN.EQ.1) call StartUpFile
 
 do XYear = 1, YearN
  if (RawFileN.EQ.YearN.AND.RawFileN.NE.1) call StartUpFile
 
  if (RawFileN.NE.(YearN*12).AND.YearHeadN.GT.0) then	! load year headers
       do XYearHead = 1, YearHeadN
        read (2,"(a80)"), Trash
       end do
  end if
  
  do XMonth = 1, MonthN
    if (RawFileN.EQ.(YearN*12)) call StartUpFile
 
    if (MonthHeadN.GT.0) then				! load month headers
      do XMonthHead = 1, MonthHeadN
        read (2,"(a80)"), Trash
      end do
    end if
    
    XCell = 0
    do XRow = 1, RowN					! load data by row,col into long vector
      if (QIntReal.EQ.1) then
      	read (2,fmt=LineFormat), (RowInt (XCol), XCol=1,ColN)
        
      	do XCol = 1, ColN
      	  XCell = XCell + 1
      	  if (ReGrid(XCell).NE.MissVal) Data(XYear,XMonth,ReGrid(XCell)) = real (RowInt(XCol))
      	end do
      else if (QIntReal.EQ.2) then
      	read (2,fmt=LineFormat), (RowReal(XCol), XCol=1,ColN)

      	do XCol = 1, ColN
      	  XCell = XCell + 1
      	  if (ReGrid(XCell).NE.MissVal) Data(XYear,XMonth,ReGrid(XCell)) = RowReal(XCol)
      	end do
      end if  
    end do

    if (RawFileN.EQ.(YearN*12)) call ShutDownFile
  end do  
  
  if (RawFileN.EQ.YearN.AND.RawFileN.NE.1) call ShutDownFile
 end do
 
 if (RawFileN.EQ.1) call ShutDownFile
 
 TotMiss = 0
 do XYear = 1, YearN
    do XMonth = 1, MonthN
      do XCell = 1, RefBoxN
        if (Data(XYear,XMonth,XCell).NE.FileMissVal) then
        	Data(XYear,XMonth,XCell) = Data(XYear,XMonth,XCell) * Multiplier
        else
        	Data(XYear,XMonth,XCell) = MissVal
        	TotMiss = TotMiss + 1
        end if
      end do
    end do
 end do
  
 print "(2(a,i12))", "   > Missing: ", TotMiss, " and Total: ", (YearN*MonthN*RefBoxN)
 print "(2a)", "   > Saving:  ", trim(OutFile(XFile)) 
 call SaveGrim (Data,RefGrid,YearAD,Bounds,Info,OutFile(XFile),"    ",SaveSuffix,NoZip=1)
  
 Data = MissVal
end do

end subroutine LoadData

!*******************************************************************************
! get file unzipped, open, and the headers read

subroutine StartUpFile

XRawFile = XRawFile + 1

LoadFile = LoadPath (RawFile(XFile,XRawFile),"    ")
print "(2a)", "   > Loading: ", trim(LoadFile)

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 not zipped
end if

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

if (FileHeadN.GT.0) then				! load file headers
      do XFileHead = 1, FileHeadN
        read (2,"(a80)"), Trash
      end do
end if

end subroutine StartUpFile

!*******************************************************************************
! get file shut, rezipped

subroutine ShutDownFile

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

end subroutine ShutDownFile

!*******************************************************************************
! finalise

subroutine Finalise

close (99)

print*
				! Grid
deallocate (RowInt,RowReal,Data,RefGrid,ReGrid,YearAD,RawFile,OutFile, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Finalise: Deallocation failure #####"

end subroutine Finalise

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

end program GridToGrim