! module ArrayMod
! deals with basic array operations
! contains: RandomiseSeed, UnBlankReal, UnBlankInteger, GetChunkMeans

module ArrayMod

implicit none

contains

!*******************************************************************************
! get random value and make it act as seed, using time
! once this has been called, simply use random_number

subroutine RandomiseSeed

integer, dimension(1) :: Seed
real, parameter :: MissVal = -999.0
integer :: Minute,Second
character (len=80) :: ListDump, Trash, CommandLine

ListDump = '/cru/u2/f709762/data/scratch/listdump.txt'
ListDump = trim (ListDump)

CommandLine = 'ps -o etime > ' // ListDump
CommandLine = trim (CommandLine)
call system (CommandLine)

open  (2,file=ListDump,status="old",action="read")
read  (2,*), Trash
read  (2,"(a6,i2,a1,i2,a1,i2)"), Trash, Minute, Trash, Second
close (2)

CommandLine = 'rm ' // ListDump
CommandLine = trim (CommandLine)
call system (CommandLine)

Seed(1) = MissVal

if (Minute.NE.0) then
  if (Second.NE.0) then
	Seed(1) = (1000*Minute) / Second
  else
	Seed(1) = Minute
  end if
else
  if (Second.NE.0) &
	Seed(1) = Second
end if

! call random_seed (put=Seed)	! compaq fortran version
call random_seed

end subroutine RandomiseSeed

!*******************************************************************************
! weed out missing values and reform array so that it has the right size
!   to contain all, and only, valid values

subroutine UnBlankReal (InVector,OutVector)

real, pointer, dimension (:) 	:: InVector,OutVector

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: InN, ValidN
integer :: XIn, XValid

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

InN = size (InVector,1)

ValidN = 0
do XIn = 1, InN
  if (InVector(XIn).NE.MissVal) ValidN = ValidN + 1
end do

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

XValid = 0
do XIn = 1, InN
  if (InVector(XIn).NE.MissVal) then
  	XValid = XValid + 1
  	OutVector (XValid) = InVector (XIn)
  end if
end do

end subroutine UnBlankReal

!*******************************************************************************
! weed out missing values and reform array so that it has the right size
!   to contain all, and only, valid values

subroutine UnBlankInteger (InVector,OutVector)

integer, pointer, dimension (:) :: InVector,OutVector

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: InN, ValidN
integer :: XIn, XValid

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

InN = size (InVector,1)

ValidN = 0
do XIn = 1, InN
  if (InVector(XIn).NE.MissVal) ValidN = ValidN + 1
end do

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

XValid = 0
do XIn = 1, InN
  if (InVector(XIn).NE.MissVal) then
  	XValid = XValid + 1
  	OutVector (XValid) = InVector (XIn)
  end if
end do

end subroutine UnBlankInteger

!*******************************************************************************
! cut up Data into chunks, calc means, and return means in Chunks
! offset >= 0 and is the number of values to omit at the beginning

subroutine GetChunkMeans (Data,Chunks,Offset,PerLen,GapLen,MissAccept)

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

integer, intent(in) :: Offset, PerLen, GapLen

real, intent(in) :: MissAccept

real, parameter :: MissVal = -999.0

real :: SectN, Thresh, OpTot, OpEn

integer :: InN, ChunkN, LastBit
integer :: XChunk, XYear, XIn
integer :: Abandon, AllocStat, ValidTot

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

Abandon = 0					! checks arguments are acceptable
if (Offset.LT.0.OR.PerLen.LT.1.OR.GapLen.LT.0.OR.MissAccept.GT.100.OR.MissAccept.LT.0) Abandon = 1

if (Abandon.EQ.0) then				! calcs number of chunks
  InN     = size (Data,1)
  SectN   = (real(InN)-real(Offset)) / (real(PerLen)+real(GapLen))
  LastBit = InN - Offset - (floor(SectN) * (PerLen+GapLen))
  ChunkN  = floor(SectN)
  if (LastBit.GE.PerLen) ChunkN = ChunkN + 1
  if (ChunkN.LT.1) Abandon = 1
  Thresh = (100.0 - MissAccept) * real(PerLen) / 100.0  
end if

if (Abandon.EQ.0) then				! chunkify into means
  allocate (Chunks(ChunkN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GetChunkMeans: Allocation failure #####"
  Chunks = MissVal
  
  ValidTot = 0
  XIn      = Offset
  
  do XChunk = 1, ChunkN
    OpTot = 0.0
    OpEn  = 0.0
    
    do XYear = 1, PerLen			! go through period to obtain total
      XIn = XIn + 1
      if (Data(XIn).NE.MissVal) then
        OpTot = OpTot + Data(XIn)
        OpEn  = OpEn  + 1.0
      end if
    end do
    						! obtain chunk mean
    if (OpEn.GE.Thresh.AND.OpEn.GT.0) then
    	Chunks(XChunk) = OpTot / OpEn
    	ValidTot = ValidTot + 1
    end if
    
    if (XChunk.LT.ChunkN) XIn = XIn + GapLen	! miss out gap
  end do
else						! return without chunkifying
  allocate (Chunks(1), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GetChunkMeans: Allocation failure #####"
  Chunks = MissVal
end if

end subroutine GetChunkMeans

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

end module ArrayMod