! extractmod.f90
! module procedure written by Tim Mitchell in Dec 1999
! last modification on 16.02.00
! contains all the routines necessary to extract a given run from raw files
! 	BlockKey, ExtractMonth, ExtractFile

module ExtractMod

use ArrayMod

implicit none

contains

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

subroutine BlockKey (FileYearN, JobMonthN, FileYear0, FileYear1, JobMonth0, JobMonth1, &
			FilePathThis, FilePathNext, FileStyleThis, FileStyleNext, FileBlockKey)

integer, pointer, dimension (:,:) :: FileBlockKey

integer, intent (in)    :: FileYearN, FileYear0, FileYear1, JobMonthN, JobMonth0, JobMonth1
integer, intent (inout) :: FileStyleThis, FileStyleNext

character (len=80), intent (inout)  :: FilePathThis, FilePathNext

real, parameter :: MissVal = -999.0

integer :: AllocStat				! status of allocation
integer :: OverLap				! 0=none 1=overlap
integer :: XYear, XMonth
integer :: LastSeasonA				! last season contained within fileA

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

if (FileYearN.NE.10) print*, "  > ##### Cannot use files other than 10 years in length #####"

allocate (FileBlockKey (FileYearN, 12), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
FileBlockKey = 0

if ( (JobMonth1-JobMonth0).GE.0 ) then
  OverLap = 0
else
  OverLap = 1
end if

if (OverLap.EQ.0) then
  
  FilePathNext  = "blank"
  FileStyleNext = MissVal
  
  do XYear = FileYear0, FileYear1
    FileBlockKey (XYear, JobMonth0) = 1
    
    if (JobMonthN.GE.3) then
      do XMonth = (JobMonth0+1),(JobMonth1-1)
        FileBlockKey (XYear, XMonth) = 2
      end do
    end if

    FileBlockKey (XYear, JobMonth1) = 3
  end do
  
else

  if (FileYear1.EQ.FileYearN .AND. FilePathNext.NE."blank") then
    LastSeasonA = FileYear1 - 1
  else
    LastSeasonA = FileYear1
  end if
  
  do XYear = FileYear0, LastSeasonA
    FileBlockKey (XYear, JobMonth0) = 1
    
    if ((12-JobMonth0).GE.1) then
      do XMonth = (JobMonth0+1),12
        FileBlockKey (XYear, XMonth) = 2
      end do
    end if

    if (JobMonth1.GT.1) then
      do XMonth = 1,(JobMonth1-1)
        FileBlockKey ((XYear+1), XMonth) = 2
      end do
    end if

    FileBlockKey ((XYear+1), JobMonth1) = 3
  end do
  
  if (FileYear1.EQ.FileYearN .AND. FilePathNext.NE."blank") then
    
    FileBlockKey (FileYear1, JobMonth0) = 1
    
    if ((12-JobMonth0).GE.1) then
      do XMonth = (JobMonth0+1),12
        FileBlockKey (10, XMonth) = 2
      end do
    end if

    if (JobMonth1.GT.1) then
      do XMonth = 1,(JobMonth1-1)
        FileBlockKey (1, XMonth) = -2
      end do
    end if

    FileBlockKey (1, JobMonth1) = -3
  
  else

    FilePathNext  = "blank"
    FileStyleNext = MissVal
  
  end if
  
end if 

end subroutine BlockKey

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

subroutine DefineStyle (StyleIndex, FileFormat, FileMissVal, FileRowN, FileColN, FileHeaderN)

integer, intent (in)  :: StyleIndex

integer, intent (out) 			:: FileRowN, FileColN, FileHeaderN
character (len=10), intent (out) 	:: FileFormat
real, intent (out)			:: FileMissVal

real, parameter :: MissVal = -999.0

FileRowN    = MissVal		! the standard series is 101...
FileColN    = MissVal		! the series 201 and 301... are = 101... (and must be maintained thus)
FileHeaderN = MissVal		! the rule is that 101... is used unless the data is in temk (deg K) form
FileFormat  = "blank"		!	in which case, the equivalent from the 201...series must be used
FileMissVal = MissVal		! ...or unless the data is in degC * 10 form,
				!	in which case, the equivalent from the 301...series must be used

if 	(StyleIndex.EQ.101.OR.StyleIndex.EQ.201.OR.StyleIndex.EQ.301) then
  FileFormat  = "e12.5"
  FileMissVal = 10000
  FileHeaderN = 6
  FileRowN    = 1168
  FileColN    = 6
else if (StyleIndex.EQ.102.OR.StyleIndex.EQ.202.OR.StyleIndex.EQ.302) then
  FileFormat  = "f8.2"
  FileMissVal = 9999.99
  FileHeaderN = 6
  FileRowN    = 701
  FileColN    = 10
else if (StyleIndex.EQ.103.OR.StyleIndex.EQ.203.OR.StyleIndex.EQ.302) then
  FileFormat  = "e13.6"
  FileMissVal = -999.0
  FileHeaderN = 1
  FileRowN    = 876
  FileColN    = 8
else if (StyleIndex.EQ.104.OR.StyleIndex.EQ.204.OR.StyleIndex.EQ.304) then
  FileFormat  = "e13.6"
  FileMissVal = -999.0
  FileHeaderN = 1
  FileRowN    = 1024
  FileColN    = 8
else if (StyleIndex.EQ.105.OR.StyleIndex.EQ.205.OR.StyleIndex.EQ.305) then
  FileFormat  = "e16.7"
  FileMissVal = -999.9
  FileHeaderN = 0
  FileRowN    = 73
  FileColN    = 96
else if (StyleIndex.EQ.106.OR.StyleIndex.EQ.206.OR.StyleIndex.EQ.306) then
  FileFormat  = "e16.7"
  FileMissVal = -9999.0
  FileHeaderN = 0
  FileRowN    = 32400
  FileColN    = 8
else
  print*, "  > ##### DefineStyle: style not recognised #####"
end if			

end subroutine DefineStyle

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

subroutine TemktoTemp (RegN,DataVec)

real, pointer, dimension (:) :: DataVec		! def,alloc,init in call
integer, intent (in) :: RegN

real, parameter :: MissVal = -999.0
integer :: XReg

do XReg = 1, RegN
  if (DataVec(XReg).NE.MissVal) then
    DataVec (XReg) = DataVec (XReg) - 273.15
  end if
end do

end subroutine TemktoTemp

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

subroutine TempTens (RegN,DataVec)

real, pointer, dimension (:) :: DataVec		! def,alloc,init in call
integer, intent (in) :: RegN

real, parameter :: MissVal = -999.0
integer :: XReg

do XReg = 1, RegN
  if (DataVec(XReg).NE.MissVal) then
    DataVec (XReg) = DataVec (XReg) / 10.0
  end if
end do

end subroutine TempTens

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

subroutine ExtractMonth (FileUnit, FileFormat, FileMissVal, FileRowN, FileColN, FileHeaderN, &
			 ModelDataN, RegN, MapRawReg, RegSizes, GotMonth)

integer, pointer, dimension (:) :: MapRawReg, RegSizes

real, pointer, dimension (:) :: GotMonth		! def,alloc,init=0 in call

integer, intent (in) :: FileUnit, FileRowN, FileColN, FileHeaderN, &
			ModelDataN, RegN

real, intent (in) :: FileMissVal

character (len=10), intent (in) :: FileFormat

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

real, allocatable, dimension (:) :: LoadVec

integer, allocatable, dimension (:) :: MissReg

character (len=80), parameter :: ScratchFile = "/cru/u2/f709762/data/scratch/deleteme.txt"

character (len=10) :: RecFormat, RecFormatLast
character (len=4)  :: FormatPrefix, FormatPrefixLast
character (len=1)  :: HeaderTrash

integer :: FileColNLast				! no of columns in last row
integer :: XRec, XDatum, XReg
integer :: AllocStat				! status of allocation attempt
integer :: RecDatum0, RecDatum1			! record first, last data
integer :: RecTotal
integer :: ReadStat
integer :: RegMissTot

real, parameter :: MissVal = -999.0

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

allocate (LoadVec  (ModelDataN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
LoadVec = MissVal

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

FileColNLast = ModelDataN - (FileRowN-1) * FileColN

open (unit=(FileUnit+2), file=ScratchFile)
  	write (unit=(FileUnit+2),fmt="(2(I4))"), FileColN, FileColNLast
close (FileUnit+2)
open (unit=(FileUnit+2), file=ScratchFile)
  	read  (unit=(FileUnit+2),fmt="(2(A4))"), FormatPrefix, FormatPrefixLast
close (FileUnit+2)

RecFormat     = '(' // trim(adjustl(FormatPrefix))     // trim(adjustl(FileFormat)) // ')'
RecFormatLast = '(' // trim(adjustl(FormatPrefixLast)) // trim(adjustl(FileFormat)) // ')'

RecFormat     = trim(adjustl(RecFormat))
RecFormatLast = trim(adjustl(RecFormatLast))

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

RecDatum0 = 1 - FileColN
RecDatum1 = 0

do XRec = 1, FileHeaderN
  read (FileUnit, fmt="(A1)"), HeaderTrash
end do

do XRec = 1, (FileRowN-1)
  RecDatum0 = RecDatum0 + FileColN
  RecDatum1 = RecDatum1 + FileColN
  
  read (FileUnit, fmt=RecFormat), (LoadVec (XDatum), XDatum=RecDatum0,RecDatum1)
end do

RecDatum0 = RecDatum0 + FileColN
RecDatum1 = ModelDataN

read (FileUnit, fmt=RecFormatLast), (LoadVec (XDatum), XDatum=RecDatum0,RecDatum1)

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

do XDatum = 1, ModelDataN
  if (MapRawReg(XDatum).NE.nint(MissVal)) then  
    if (MissReg(MapRawReg(XDatum)).LE.0) then

      if (LoadVec(XDatum).NE.FileMissVal) then
        GotMonth(MapRawReg(XDatum)) = GotMonth(MapRawReg(XDatum)) + LoadVec (XDatum)
        MissReg(MapRawReg(XDatum))  = MissReg(MapRawReg(XDatum)) - 1
      else
        MissReg(MapRawReg(XDatum))  = 1
      end if

    end if
  end if
end do

do XReg = 1, RegN
  if (MissReg(XReg).LT.0) then
    if ((MissReg(XReg)+RegSizes(XReg)).EQ.0) then
  	GotMonth(XReg) = GotMonth(XReg) / real(RegSizes(XReg))
    else
        print*, "  > ##### ERROR: ExtractMonth: region size mismatch #####"
    end if
  else
  	GotMonth(XReg) = MissVal
  end if
end do

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

deallocate (LoadVec)
deallocate (MissReg)

end subroutine ExtractMonth

!*******************************************************************************
! requires use arraymod in main program

subroutine ExtractFile (ModelLongN, ModelLatN, ModelDataN, &
			FilePathThis, FilePathNext, FileStyleThis, FileStyleNext, &
			RegN, JobMonthN, FileYearN, FileYear0, &
			MapRawReg, RegSizes, FileBlockKey, GotDecade)

integer, pointer, dimension (:,:) :: FileBlockKey
integer, pointer, dimension (:)   :: MapRawReg, RegSizes

real, pointer, dimension (:,:)  :: GotDecade	! allocated and dealloc in call,init=here
real, pointer, dimension (:)    :: GotMonth	! allocated and deallocated here

character (len=80), intent (in) :: FilePathThis, FilePathNext

integer, intent (in) :: ModelLongN, ModelLatN, ModelDataN, RegN, JobMonthN, FileYearN, FileYear0
integer, intent (in) :: FileStyleThis, FileStyleNext
			
!***************************************

integer, allocatable, dimension (:,:) :: MissDecade

integer :: XYear, XMonth, XReg, XRec, XChar
integer :: RandomNumber
integer :: AllocStat
integer :: FileUnit
integer :: ReadStatus
integer :: ThisSeason			! season/year 1...10
integer :: FileHeaderN, FileRowN, FileColN

character (len=80) :: UnZippedDump
character (len=80) :: GetFile
character (len=10) :: FileFormat
character (len=6)  :: RandomCha
character (len=1)  :: UnwantedMonth

real :: FileMissVal, RandomReal			

real, parameter :: MissVal = -999.0

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

call DefineStyle (FileStyleThis,FileFormat,FileMissVal,FileRowN,FileColN,FileHeaderN)

allocate (GotMonth  (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
GotMonth = 0.0

allocate (MissDecade(RegN,FileYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
MissDecade = 0

GotDecade = 0.0		! deliberately initialised to zero, because this is accumulated into
UnZippedDump = "blank"

open (1, file=FilePathThis, status="old", iostat=ReadStatus)
if (ReadStatus .EQ. 0) then
	close (1)
	GetFile = FilePathThis
else
        call RandomiseSeed
        call random_number (RandomReal)
        
        RandomNumber = RandomReal * 999999
        
!        RandomNumber = 1038
!        do XChar = 1, len(trim(adjustr(FilePathThis)))
!          RandomNumber = RandomNumber + iachar (FilePathThis(XChar:XChar))
!        end do

        open  (1,file='/cru/u2/f709762/data/scratch/temp.dat',status='replace')
        write (1,"(I6)"), RandomNumber
        close (1)
        open  (1,file='/cru/u2/f709762/data/scratch/temp.dat',status='old')
        read  (1,"(A6)"), RandomCha
        close (1)
        
        UnZippedDump = '/cru/u2/f709762/data/scratch/eraseme' // trim(adjustl(RandomCha))
        
	call system ('uncompress -c ' // FilePathThis // ' > ' // UnZippedDump)
	GetFile = UnZippedDump
end if

ThisSeason = FileYear0 - 1
FileUnit = 7

open (unit=FileUnit,file=GetFile,status="old",access="sequential",action="read",form="formatted") 

do XYear = 1, FileYearN
  do XMonth = 1, 12
    
    if (FileBlockKey(XYear,XMonth) .LE. 0) then      
      do XRec = 1, FileHeaderN
        read (FileUnit, "(A1)"), UnwantedMonth
      end do
      
      do XRec = 1, FileRowN
        read (FileUnit, "(A1)"), UnwantedMonth
      end do
    end if
    
    if (FileBlockKey(XYear,XMonth) .EQ. 1) ThisSeason = ThisSeason + 1
    
    if (FileBlockKey(XYear,XMonth) .GE. 1) then
      GotMonth = 0.0
      call ExtractMonth (FileUnit,FileFormat,FileMissVal,FileRowN,FileColN,FileHeaderN, &
      			 ModelDataN,RegN,MapRawReg,RegSizes,GotMonth)
      if (FileStyleThis.GE.200.AND.FileStyleThis.LE.299) call TemkToTemp (RegN,GotMonth)
      if (FileStyleThis.GE.300.AND.FileStyleThis.LE.399) call TempTens   (RegN,GotMonth)
      
      do XReg = 1, RegN
        if (MissDecade(XReg,ThisSeason).EQ.0) then 
           if (GotMonth(XReg).NE.MissVal) then
        	GotDecade(XReg,ThisSeason) = GotDecade(XReg,ThisSeason) + GotMonth(XReg)
           else
                MissDecade(XReg,ThisSeason) = 1
           end if
        end if
      end do      
    end if
    
    if (FileBlockKey(XYear,XMonth) .EQ. 3) then
      do XReg = 1, RegN
       if (MissDecade(XReg,ThisSeason).EQ.0) then 
         GotDecade(XReg,ThisSeason) = GotDecade(XReg,ThisSeason) / real(JobMonthN)
       else
         GotDecade(XReg,ThisSeason) = MissVal
       end if
      end do
    end if
  end do
end do

close (FileUnit)

if (UnZippedDump.NE."blank") call system ('rm ' // UnZippedDump)

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

if (FilePathNext.NE."blank") then

  call DefineStyle (FileStyleNext,FileFormat,FileMissVal,FileRowN,FileColN,FileHeaderN)
  
  UnZippedDump = "blank"
  
  open (1, file=FilePathNext, status="old", iostat=ReadStatus)
  
  if (ReadStatus .EQ. 0) then
	close (1)
	GetFile = FilePathNext
  else
        RandomNumber = 1123
        do XChar = 1, len(trim(adjustr(FilePathNext)))
          RandomNumber = RandomNumber + iachar (FilePathNext(XChar:Xchar))
        end do

        open  (2,file='/cru/u2/f709762/data/scratch/temp.dat',status='replace')
        write (2,"(I6)"), RandomNumber
        close (2)
        open  (2,file='/cru/u2/f709762/data/scratch/temp.dat',status='old')
        read  (2,"(A6)"), RandomCha
        close (2)
        
        UnZippedDump = '/cru/u2/f709762/data/scratch/eraseme' // trim(adjustl(RandomCha))
        
	call system ('uncompress -c ' // FilePathNext // ' > ' // UnZippedDump)
	GetFile = UnZippedDump    
  end if

  FileUnit = 7
  
  open (unit=FileUnit,file=GetFile,status="old",access="sequential",action="read",form="formatted") 

  do XMonth = 1, 12
    
    if (FileBlockKey(1,XMonth) .LE. -2) then
      GotMonth = 0.0
      call ExtractMonth (FileUnit,FileFormat,FileMissVal,FileRowN,FileColN,FileHeaderN, &
      			 ModelDataN,RegN,MapRawReg,RegSizes,GotMonth)
      
      if (FileStyleThis.GE.200.AND.FileStyleThis.LE.299) call TemkToTemp (RegN,GotMonth)
      if (FileStyleThis.GE.300.AND.FileStyleThis.LE.399) call TempTens   (RegN,GotMonth)

      if (MissDecade(XReg,ThisSeason).EQ.0) then 
           if (GotMonth(XReg).NE.MissVal) then
        	GotDecade(XReg,ThisSeason) = GotDecade(XReg,ThisSeason) + GotMonth(XReg)
           else
                MissDecade(XReg,ThisSeason) = 1
           end if
      end if
    end if
    
    if (FileBlockKey(1,XMonth) .EQ. -3) then
      if (MissDecade(XReg,ThisSeason).EQ.0) then 
         GotDecade(XReg,ThisSeason) = GotDecade(XReg,ThisSeason) / real(JobMonthN)
      else
         GotDecade(XReg,ThisSeason) = MissVal
      end if
    end if
    
  end do

  close (FileUnit)
  
  if (UnZippedDump.NE."blank") call system ('rm ' // UnZippedDump)

end if

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

deallocate (GotMonth)
deallocate (MissDecade)

end subroutine ExtractFile

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

end module ExtractMod