! signoi.f90
! f90 main program written on 03.05.00 by Tim Mitchell
! last modification on 03.05.00
! f90 -o signoi initialmod.f90 loadmod.f90 savemod.f90 transformmod.f90 signoi.f90

program SigNoi

use InitialMod
use LoadMod
use SaveMod
use TransformMod

implicit none

real, pointer, dimension (:,:,:)	:: TimLoaded
real, pointer, dimension (:,:)		:: TimProc, LinToSave, FileSeries
real, pointer, dimension (:)		:: NoiseLoaded, Vectorised, GloToSave

integer, pointer, dimension (:) 	:: WorkMapRawReg, WorkRegSizes
integer, pointer, dimension (:) 	:: WorkADYear
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw, WorkMapIDLReg

character (len=80), pointer, dimension (:) 	:: FileRegNames, WorkRegNames

real, parameter :: MissVal = -999.0

real :: YearThresh, RegThresh
real :: WorkMissAccept
real :: CalcEn, CalcTot
real :: PerCentMiss, TotalMiss

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN, WorkRegN
integer :: WorkMonth0, WorkMonth1, WorkMonthN, WorkYearN, WorkDecN
integer :: WorkTimN, FileRegN, RegValidN
integer :: AllocStat, ReadStatus
integer :: XFile, XLin, XYear, XFind, XTim, XReg
integer :: QMain, QMissAccept, QSliceify, QSave
integer :: ExitCheck
integer :: YearMissN
integer :: SelectFile, SelectADYear, ChosenADYear, ChosenFileN

character (len=10) :: WorkGridTitle
character (len=80) :: WorkGridFilePath, WorkLinTitle, WorkRegTitle, WorkGloFilePath, Blank

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

Blank = ""

call Initialise
call GrabTim

print*, "  > We load the noise (RMSE sampling error) from .glo file."
call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, NoiseLoaded)

open (99,file='/cru/u2/f709762/data/scratch/logsignoi.dat',access='sequential',status='replace',action='write')

print*, "  > Alter % acceptable missing (1=no,2=yes), currently: ", WorkMissAccept
do
	read (*,*,iostat=ReadStatus), QMissAccept
	if (ReadStatus.LE.0.AND.QMissAccept.GE.1.AND.QMissAccept.LE.2) exit
end do
if (QMissAccept.EQ.2) call SetMissAccept (WorkMissAccept)

print*, "  > Transform regional series into 30-year slices ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QSliceify
	if (ReadStatus.LE.0.AND.QSliceify.GE.1.AND.QSliceify.LE.2) exit
end do
if (QSliceify.EQ.2) call IntoSlices

call CalcSignal

print*, "  > Save signals as .tim ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QSave
	if (ReadStatus.LE.0.AND.QSave.GE.1.AND.QSave.LE.2) exit
end do
if (QSave.EQ.2) call SaveTimProc

call CalcSignalNoise

print*, "  > Save signal-to-noise ratios as .tim ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QSave
	if (ReadStatus.LE.0.AND.QSave.GE.1.AND.QSave.LE.2) exit
end do
if (QSave.EQ.2) call SaveTimProc

print*, "  > Save global av signal-to-noise ratios as .lin ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QSave
	if (ReadStatus.LE.0.AND.QSave.GE.1.AND.QSave.LE.2) exit
end do
if (QSave.EQ.2) call SaveGlobalAv

print*, "  > Save single year signal-to-noise ratios as .glo ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QSave
	if (ReadStatus.LE.0.AND.QSave.GE.1.AND.QSave.LE.2) exit
end do
if (QSave.EQ.2) call SaveOneYear

deallocate (TimProc,WorkRegNames)
print*
close (99)

contains

!*******************************************************************************
! initialise

subroutine Initialise

print*, "  > ***** Signal-to-Noise *****"
print*, "  > operates on .tim's provided noise is already in .glo"
print*

WorkMissAccept = 10.0

print*, "  > Select parameters that will govern operations."
call GridSelect   (WorkGrid,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call PeriodSelect (WorkYearN,WorkDecN,WorkADYear)
call RegSelect    (WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkMapIDLReg,WorkRegSizes,WorkRegNames,&
		   WorkRegTitle,WorkRegN)

end subroutine Initialise

!*******************************************************************************
! load .tim

subroutine GrabTim

print*, "  > Select the number of .tim files to load."
do
	read (*,*,iostat=ReadStatus), WorkTimN
	if (ReadStatus.LE.0.AND.WorkTimN.GE.1) exit
end do

allocate (TimLoaded    (WorkTimN,WorkRegN,WorkYearN), &
	  WorkRegNames (WorkRegN),	              stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabTim: Allocation failure #####"
TimLoaded    = MissVal
WorkRegNames = "blank"

do XTim = 1, WorkTimN
  print*, "  > File number : ", XTim
  
  do
    call LoadTim (WorkYearN,WorkADYear,FileRegN,FileRegNames,FileSeries)
    
    if (FileRegN.NE.WorkRegN) then
      print*, "  > ##### Incorrect no. of regions in file. Try again #####"
      ExitCheck = 1
      
      deallocate (FileRegNames,FileSeries,stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
    else
      ExitCheck = 0
    end if
    
    if (ExitCheck.EQ.0) exit
  end do
  
  do XReg = 1, WorkRegN
    WorkRegNames (XReg) = FileRegNames (XReg)
    
    do XYear = 1, WorkYearN
      TimLoaded (XTim,XReg,XYear) = FileSeries (XReg,XYear)
    end do
  end do
  
  deallocate (FileRegNames,FileSeries,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabTim: Deallocation failure #####"
 end do

TotalMiss = 0.0
do XTim = 1, WorkTimN
  do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      if (TimLoaded(XTim,XReg,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1.0
    end do
  end do
end do
PerCentMiss = 100.0 * TotalMiss / (WorkTimN*WorkRegN*WorkYearN)
print "(a38,f8.2)", "   > After loading: percentage missing: ", PerCentMiss

end subroutine GrabTim

!*******************************************************************************
! sliceifying into 30y slices

subroutine IntoSlices

print*, "  > Sliceifying each file and region individually..."

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

do XTim = 1, WorkTimN
  do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      Vectorised (XYear) = TimLoaded (XTim,XReg,XYear)
    end do
    
    call SimpleSlice (WorkYearN,30,WorkMissAccept,Vectorised)
    
    do XYear = 1, WorkYearN
      TimLoaded (XTim,XReg,XYear) = Vectorised (XYear)
    end do
  end do
  
  if (WorkTimN.GT.1) print*, "  > Have slicified file: ", XTim
end do

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

TotalMiss = 0.0
do XTim = 1, WorkTimN
  do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      if (TimLoaded(XTim,XReg,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1.0
    end do
  end do
end do
PerCentMiss = 100.0 * TotalMiss / (WorkTimN*WorkRegN*WorkYearN)
print "(a42,f8.2)", "   > After sliceifying: percentage missing: ", PerCentMiss

end subroutine IntoSlices

!*******************************************************************************
! calcs signal RMSE

subroutine CalcSignal

print*, "  > Calculating signals..."

allocate (TimProc (WorkRegN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcSignal: Allocation failure #####"
TimProc = MissVal

RegThresh  = real(WorkTimN)*(100.0-WorkMissAccept)/100.0

do XYear = 1, WorkYearN
  do XReg = 1, WorkRegN
    CalcEn  = 0.0
    CalcTot = 0.0
    
    do XTim = 1, WorkTimN
      if (TimLoaded(XTim,XReg,XYear).NE.MissVal) then
        CalcEn  = CalcEn  + 1.0
        CalcTot = CalcTot + (TimLoaded (XTim,XReg,XYear) ** 2)
      endif
    end do
    
    if (CalcEn.GE.RegThresh) then
    	TimProc(XReg,XYear) = CalcTot / CalcEn
    	TimProc(XReg,XYear) = sqrt (TimProc(XReg,XYear))
    end if
  end do
end do

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

TotalMiss = 0.0
do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      if (TimProc(XReg,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1.0
    end do
end do
PerCentMiss = 100.0 * TotalMiss / (WorkRegN*WorkYearN)
print "(a41,f8.2)", "   > After signalling: percentage missing: ", PerCentMiss

end subroutine CalcSignal

!*******************************************************************************
! calcs signal-to-noise ratios

subroutine CalcSignalNoise

print*, "  > Dividing signals by noise..."

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (TimProc(XReg,XYear).NE.MissVal) then
      if (NoiseLoaded(XReg).NE.MissVal) then
        TimProc(XReg,XYear) = TimProc(XReg,XYear) / NoiseLoaded(XReg)
      else
        TimProc(XReg,XYear) = MissVal
      end if
    end if
  end do
end do

TotalMiss = 0.0
do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      if (TimProc(XReg,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1.0
    end do
end do
PerCentMiss = 100.0 * TotalMiss / (WorkRegN*WorkYearN)
print "(a38,f8.2)", "   > After noising: percentage missing: ", PerCentMiss

end subroutine CalcSignalNoise

!*******************************************************************************
! save TimProc to .tim

subroutine SaveTimProc

print*, "  > Enter the .tim file title: "
do
	read (*,*,iostat=ReadStatus), WorkLinTitle
	if (ReadStatus.LE.0.AND.WorkLinTitle.NE."") exit
end do

call SaveTim (WorkRegN, WorkYearN, Blank, Blank, WorkRegNames, WorkADYear, TimProc)

end subroutine SaveTimProc

!*******************************************************************************
! save global average to .lin

subroutine SaveGlobalAv

print*, "  > Globally averaging..."

allocate (LinToSave (1,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGlobalAv: Allocation failure #####"
LinToSave = MissVal

YearThresh = real(WorkYearN)*(100.0-WorkMissAccept)/100.0

do XYear = 1, WorkYearN
 CalcEn  = 0.0
 CalcTot = 0.0
  
 do XReg = 1, WorkRegN
  if (TimProc(XReg,XYear).NE.MissVal) then
    CalcEn  = CalcEn  + 1.0
    CalcTot = CalcTot + TimProc(XReg,XYear)
  end if
  
  if (CalcEn.GE.YearThresh) LinToSave (1,XYear) = CalcTot / CalcEn
 end do
end do

TotalMiss = 0.0
do XYear = 1, WorkYearN
      if (LinToSave(1,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1.0
end do
PerCentMiss = 100.0 * TotalMiss / WorkYearN
print "(a41,f8.2)", "   > After globalling: percentage missing: ", PerCentMiss

call SaveLin (1,WorkYearN,WorkRegNames,WorkADYear,LinToSave)

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

end subroutine SaveGlobalAv

!*******************************************************************************
! save global slice for single year

subroutine SaveOneYear

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

ChosenADYear =  0

do
  do
    print*, "  > Enter the year AD (-99=end): "
    do
        read (*,*,iostat=ReadStatus), SelectADYear
	if (ReadStatus.LE.0) exit
    end do
    
    if (SelectADYear.NE.-99) then
      do XYear = 1, WorkYearN
        if (WorkADYear(XYear).EQ.SelectADYear) ChosenADYear = XYear
      end do
    else
      ChosenADYear = 1
    end if
    
    if (ChosenADYear.EQ.0) print*, "  > Year out of range. Try again."
    if (ChosenADYear.GE.1) exit
  end do
  
  if (SelectADYear.NE.-99) then
    do XReg = 1, WorkRegN
      GloToSave(XReg) = TimProc(XReg,ChosenADYear)
    end do
    
    print*, "  > Enter the .glo file title: "
    do
	read (*,*,iostat=ReadStatus), WorkLinTitle
	if (ReadStatus.LE.0.AND.WorkLinTitle.NE."") exit
    end do

    call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,GloToSave,WorkMapIDLReg)
  end if
  
  if (SelectADYear.EQ.-99) exit
end do

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

end subroutine SaveOneYear

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

end program SigNoi