! operatetim.f90
! f90 main program written on 05.05.00 by Tim Mitchell
! last modification on 16.05.00
! f90 -o operatetim initialmod.f90 loadmod.f90 savemod.f90 scalemod.f90 transformmod.f90 operatetim.f90

program OperateTim

use InitialMod
use LoadMod
use SaveMod
use ScaleMod
use TransformMod

implicit none

real, pointer, dimension (:,:,:)	:: TimLoaded
real, pointer, dimension (:,:)		:: LinToSave, FileSeries, OpTim, LinLoaded
real, pointer, dimension (:,:)		:: FileData, OpError, OpSampN, OpStDev, Scalar
real, pointer, dimension (:)		:: Vectorised, GloLoaded
real, pointer, dimension (:)		:: Pattern, OpGlo

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

integer, allocatable, dimension (:)	:: UseTim

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

real, parameter :: MissVal = -999.0

real :: YearThresh, RegThresh, MemThresh
real :: MissAccept
real :: CalcEn, CalcTot, OpTotal, OpEn, OpSqTot, OpTot
real :: PerCentMiss, TotalMiss, TimMin, TimMax
real :: Aye, RealConstant

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN, WorkRegN
integer :: WorkMonth0, WorkMonth1, WorkMonthN, WorkYearN, WorkDecN
integer :: WorkTimN, FileRegN, RegValidN, RegYearValidN, FileLineN
integer :: AllocStat, ReadStatus
integer :: XFile, XLin, XYear, XFind, XTim, XReg, XLong, XLine
integer :: QMain, QMissAccept, QSliceify, QSave, QAnother, QSampPop, QSaveErr, QCalcSig
integer :: ExitCheck, QuitProc
integer :: YearMissN
integer :: SelectFile, SelectADYear, ChosenADYear, ChosenFileN, SelectTim, UseTimN, SelectedLine
integer :: SliceLen, Operator, IntConstant

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

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

call Initialise
call GrabTim

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

QAnother = -1

do
  QuitProc = 0
  
  if      (QAnother.EQ.0) then
  	call SetMissAccept (MissAccept)
  else if (QAnother.EQ.1) then
  	SliceLen = 30
  	call IntoSlices
  else if (QAnother.EQ.2) then
  	UseTimN = 1
  	call SelectFiles
  	call ConvertGlo
  else if (QAnother.EQ.4) then
  	UseTimN = 2
  	call SelectFiles
  	if (QuitProc.EQ.0) call DivideAB
  else if (QAnother.EQ.5) then
  	UseTimN = 1
  	call SelectFiles
  	call OperateAGlo
  else if (QAnother.EQ.6) then
        call SelectUseTimN
  	call SelectFiles
  	call AverageTim
  else if (QAnother.EQ.7) then
  	UseTimN = 1
  	call SelectFiles
  	call AverageTim
  else if (QAnother.EQ.8) then
        call SelectUseTimN
  	call SelectFiles
  	call StDevTim
  else if (QAnother.EQ.9) then
        call CalcScaleEq
  else if (QAnother.EQ.10) then
        call ApplyScaling
  else if (QAnother.EQ.11) then
  	SliceLen = -1
  	call IntoSlices
  else if (QAnother.EQ.12) then
  	UseTimN = 1
  	call SelectFiles
  	call OperateALin
  else if (QAnother.EQ.13) then
  	UseTimN = 1
  	call SelectFiles
  	call OperateAConstant
  else
    print*, "  >  0. Change % acceptable missing (currently: ", MissAccept
    print*, "  >  1. Transform into 30y slices"
    print*, "  >  2. Extract year-slice(s) from A into .glo(s)"
    print*, "  >  4. Divide A by B"
    print*, "  >  5. Operate on A with specific .glo"
    print*, "  >  6. Average specified A,B,... into .tim"
    print*, "  >  7. Save specified A into .tim"
    print*, "  >  8. Stdev specified A,B,... into .tim"
    print*, "  >  9. Calc scaling eq using specified A,B,..."
    print*, "  > 10. Scale specified A,B,..."
    print*, "  > 11. Transform into slices of specified length"
    print*, "  > 12. Operate on A with specific .lin"
    print*, "  > 13. Operate on A with constant"
    print*, "  > 99. Exit"  
  end if
  
  print*
  print*, "  > Main menu: make your choice: "
  do
	read (*,*,iostat=ReadStatus), QAnother
	if (ReadStatus.LE.0 .AND. QAnother.GE.0 .AND. QAnother.LE.99) exit
  end do
  
  if (QAnother.EQ.99) exit
end do

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

print*
close (99)

contains

!*******************************************************************************
! calc scaling eq

subroutine CalcScaleEq

if (WorkTimN.GT.1) then
  print*, "  > Select the number of .tim files from which to calc scaling eq."
  do
	  read (*,*,iostat=ReadStatus), UseTimN
	  if (ReadStatus.LE.0.AND.UseTimN.GE.1.AND.UseTimN.LE.WorkTimN) exit
  end do
else
  UseTimN = 1
end if

call SelectFiles       ! UseTim (1...UseTimN) referring to TimLoaded (XTim,XReg,XYear)
call LoadScalar        ! Scalar (1...UseTimN,1...WorkYearN)
print*, "  > Calculating pattern..."

allocate (OpGlo    (WorkRegN),          &
	  FileData (UseTimN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcScaleEq: Allocation failure #####"
OpGlo    = MissVal
FileData = MissVal

do XReg = 1, WorkRegN
  do XTim = 1, UseTimN
    do XYear = 1, WorkYearN
      FileData (UseTim(XTim),XYear) = TimLoaded (UseTim(XTim),XReg,XYear)
    end do
  end do
  
  call LinearLSRZero (UseTimN,WorkYearN,Scalar,FileData,Aye)
  
  OpGlo (XReg) = Aye
end do

print*, "  > Enter the .glo file title to save scaling pattern: "
do
	read (*,*,iostat=ReadStatus), SaveTitle
	if (ReadStatus.LE.0 .AND. SaveTitle.NE."") exit
end do
call SaveGlo (WorkLongN, WorkLatN, WorkRegN, SaveTitle, WorkGridFilePath, &
			 OpGlo, WorkMapIDLReg)

deallocate (OpGlo,Scalar,UseTim,FileData, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcScaleEq: Deallocation failure #####"

end subroutine CalcScaleEq

!*******************************************************************************
! apply scaling

subroutine ApplyScaling

if (WorkTimN.GT.1) then
  print*, "  > Select the number of .tim files to estimate by scaling."
  do
	  read (*,*,iostat=ReadStatus), UseTimN
	  if (ReadStatus.LE.0.AND.UseTimN.GE.1.AND.UseTimN.LE.WorkTimN) exit
  end do
else
  UseTimN = 1
end if

call SelectFiles
call LoadScalar

print*, "  > Load the pattern."
call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, Pattern)

call ScalingError   

deallocate (Pattern,Scalar,UseTim, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ScalingMenu: Allocation failure #####"

print*, "  > Save scaling errors to .tim ? (1=no,2=yes)"
print*, "  > (this file cannot be used to calc sig, as n is unknown)"
do
	  read (*,*,iostat=ReadStatus), QSaveErr
	  if (ReadStatus.LE.0.AND.QSaveErr.GE.1.AND.QSaveErr.LE.2) exit
end do
if (QSaveErr.EQ.2) then
  print*, "  > Enter the .tim file title to save errors: "
  do
	read (*,*,iostat=ReadStatus), SaveTitle
	if (ReadStatus.LE.0 .AND. SaveTitle.NE."") exit
  end do
  call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpError)
end if

print*, "  > Calc sig of scaling errors ? (1=no,2=yes)"
do
	  read (*,*,iostat=ReadStatus), QCalcSig
	  if (ReadStatus.LE.0.AND.QCalcSig.GE.1.AND.QCalcSig.LE.2) exit
end do
if (QCalcSig.EQ.2) call ScalingErrorSig  

deallocate (OpError,OpSampN, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ScalingMenu: Allocation failure #####"

end subroutine ApplyScaling

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

subroutine Initialise

print*, "  > ***** Operations on .tim files *****"
print*

MissAccept = 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

if (WorkLongN.EQ.96.AND.WorkLatN.EQ.73.AND.WorkRegN.EQ.7008) then
  print*, "  > Converting regperbox into rpb..."
  
  do XLong = 1, WorkLongN
    XReg = WorkMapIDLReg (XLong,1)
    TimLoaded (1:WorkTimN,XReg,1:WorkYearN) = MissVal
    XReg = WorkMapIDLReg (XLong,WorkLatN)
    TimLoaded (1:WorkTimN,XReg,1:WorkYearN) = MissVal
  end do
end if

end subroutine GrabTim

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

subroutine IntoSlices

if (SliceLen.EQ.-1) then
  print*, "  > Enter the slice length:"
  do
	  read (*,*,iostat=ReadStatus), SliceLen
	  if (ReadStatus.LE.0.AND.SliceLen.GE.2.AND.SliceLen.LT.WorkYearN) exit
  end do
end if

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,SliceLen,MissAccept,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

!*******************************************************************************
! select one or more files from those loaded

subroutine SelectUseTimN

if (WorkTimN.EQ.1) then
  print*, "  > Insufficient files loaded to perform that operation."
  QuitProc = 1
else
  print*, "  > Enter the number of .tim to use: "
  do
	  read (*,*,iostat=ReadStatus), UseTimN
	  if (ReadStatus.LE.0.AND.UseTimN.GT.1.AND.UseTimN.LE.WorkTimN) exit
  end do
end if

end subroutine SelectUseTimN

!*******************************************************************************
! select one or more files from those loaded

subroutine SelectFiles

write (99,*), "### SelectFiles"

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

if (UseTimN.GT.WorkTimN) then
 print*, "  > Insufficient files loaded to perform that operation."
 QuitProc = 1
else
 if (WorkTimN.GT.1) then
   print*, "  > Select .tim files totalling: ", UseTimN  
   
   do XTim = 1, UseTimN
     print*, "  > Select .tim file number:", XTim 
     do
	  read (*,*,iostat=ReadStatus), SelectTim
	  if (ReadStatus.LE.0.AND.SelectTim.GE.1.AND.SelectTim.LE.WorkTimN) exit
     end do
     UseTim (XTim) = SelectTim
     write (99,"(2i4)"), XTim, UseTim (XTim)
   end do
 else
     UseTim (1) = 1
     write (99,"(2i4)"), 1, 1
 end if
end if

end subroutine SelectFiles

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

subroutine ConvertGlo

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

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear)
  end do
end do

call TimToGlo (WorkLongN, WorkLatN, WorkRegN, WorkYearN, WorkGridFilePath, &
			WorkADYear, WorkMapIDLReg, OpTim)

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

end subroutine ConvertGlo

!*******************************************************************************
! divide one .tim by another .tim

subroutine DivideAB

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

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal.AND.TimLoaded(UseTim(2),XReg,XYear).NE.MissVal) then
      if (TimLoaded(UseTim(2),XReg,XYear).NE.0) then
        OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) / TimLoaded(UseTim(2),XReg,XYear)
      end if
    end if
  end do
end do

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

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

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

end subroutine DivideAB

!*******************************************************************************
! divide or subtract one .tim by a .glo

subroutine OperateAGlo

print*, "  > Divide by (=1), multiply by (=2), subtract (=3), or add (=4) .glo ?"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.4) exit
end do

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

call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, GloLoaded)

do XReg = 1, WorkRegN
  if (GloLoaded(XReg).NE.MissVal) then
    do XYear = 1, WorkYearN
      if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal) then
        if (Operator.EQ.1) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) / GloLoaded (XReg)
        if (Operator.EQ.2) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) * GloLoaded (XReg)
        if (Operator.EQ.3) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) - GloLoaded (XReg)
        if (Operator.EQ.4) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) + GloLoaded (XReg)
      end if
    end do
  end if
end do

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

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

deallocate (GloLoaded, OpTim, UseTim, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OperateAGlo: Deallocation failure #####"

end subroutine OperateAGlo

!*******************************************************************************
! divide or subtract one .tim by a .lin

subroutine OperateALin

print*, "  > Divide by (=1), multiply by (=2), subtract (=3), or add (=4) .lin ?"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.4) exit
end do

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

call LoadLin (WorkYearN, WorkADYear, FileLineN, FileLineNames, LinLoaded)

if (FileLineN.EQ.1) then
  SelectedLine = 1
else
  print*, "  > Enter the index of the line to adopt as the scalar (0=list):"
  do
    do
	read (*,*,iostat=ReadStatus), SelectedLine
	if (ReadStatus.LE.0 .AND. SelectedLine.GE.0 .AND. SelectedLine .LE. FileLineN) exit
    end do
  
    if (SelectedLine.EQ.0) then
      do XLine = 1, FileLineN
        print*, XLine, FileLineNames(XLine)
      end do  
    end if
  
    if (SelectedLine.GT.0) exit
  end do    
end if

do XYear = 1, WorkYearN
  if (LinLoaded(SelectedLine,XYear).NE.MissVal) then
    do XReg = 1, WorkRegN
      if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal) then
        if (Operator.EQ.1) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) / LinLoaded(SelectedLine,XYear)
        if (Operator.EQ.2) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) * LinLoaded(SelectedLine,XYear)
        if (Operator.EQ.3) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) - LinLoaded(SelectedLine,XYear)
        if (Operator.EQ.4) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) + LinLoaded(SelectedLine,XYear)
      end if
    end do
  end if
end do

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

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

deallocate (LinLoaded, FileLineNames, OpTim, UseTim, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OperateALin: Deallocation failure #####"

end subroutine OperateALin

!*******************************************************************************
! divide or subtract one .tim by a .glo

subroutine OperateAConstant

print*, "  > Divide =1, multiply =2, minus =3, add =4, sqroot =5, exponentiate =6"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.6) exit
end do

if (Operator.LT.5) then
  print*, "  > Enter constant (real):"
  do
	read (*,*,iostat=ReadStatus), RealConstant
	if (ReadStatus.LE.0) exit
  end do
else if (Operator.EQ.6) then
  print*, "  > Enter constant (positive integer):"
  do
	read (*,*,iostat=ReadStatus), IntConstant
	if (ReadStatus.LE.0.AND.IntConstant.GE.0) exit
  end do
end if

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

do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal) then
        if (Operator.EQ.1) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) /  RealConstant
        if (Operator.EQ.2) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) *  RealConstant
        if (Operator.EQ.3) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) -  RealConstant
        if (Operator.EQ.4) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) +  RealConstant
        if (Operator.EQ.5) OpTim (XReg,XYear) = sqrt (TimLoaded(UseTim(1),XReg,XYear))
        if (Operator.EQ.6) OpTim (XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) ** IntConstant
      end if
    end do
end do

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

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

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

end subroutine OperateAConstant

!*******************************************************************************
! average together .tim and save them

subroutine AverageTim

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

MemThresh     = (100.0 - MissAccept) * UseTimN / 100.0
RegYearValidN = 0.0

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    OpTotal = 0.0
    OpEn    = 0.0
    
    do XTim = 1, UseTimN
      if (TimLoaded(UseTim(XTim),XReg,XYear).NE.MissVal) then
        OpTotal = OpTotal + TimLoaded(UseTim(1),XReg,XYear)
        OpEn    = OpEn    + 1
      end if
    end do
    
    if (OpEn.GE.MemThresh) then
    	OpTim (XReg,XYear) = OpTotal / OpEn
    	RegYearValidN = RegYearValidN + 1
    end if
  end do
end do

print*, "  > Total region-years valid: ", RegYearValidN
print*, "  > Enter the .tim file title to save: "
do
	read (*,*,iostat=ReadStatus), SaveTitle
	if (ReadStatus.LE.0 .AND. SaveTitle.NE."") exit
end do

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

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

end subroutine AverageTim

!*******************************************************************************
! stdev of .tim and save result

subroutine StDevTim

print*, "  > Calc sample (=1) or population (=2) stdev ? "
do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0 .AND. QSampPop.GE.1 .AND. QSampPop.LE.2) exit
end do
print*, "  > Calculating standard deviations..."

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

MemThresh     = (100.0 - MissAccept) * UseTimN / 100.0
RegYearValidN = 0.0

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    OpTotal = 0.0
    OpSqTot = 0.0
    OpEn    = 0.0
    
    do XTim = 1, UseTimN
      if (TimLoaded(UseTim(XTim),XReg,XYear).NE.MissVal) then
        OpTotal = OpTotal +  TimLoaded(UseTim(XTim),XReg,XYear)
        OpSqTot = OpSqTot + (TimLoaded(UseTim(XTim),XReg,XYear) ** 2)
        OpEn    = OpEn    +  1
      end if
    end do
    
    if (OpEn.GE.MemThresh) then
      if (QSampPop.EQ.1) then
    	  OpTim (XReg,XYear) = (OpSqTot / OpEn) - ( (OpTotal / OpEn) ** 2)
    	  OpTim (XReg,XYear) =  sqrt (OpTim (XReg,XYear))
    	  RegYearValidN      =  RegYearValidN + 1
      else
    	if (OpEn.GE.2) then
    	  OpTim (XReg,XYear) = (OpSqTot / OpEn) - ( (OpTotal / OpEn) ** 2)
    	  OpTim (XReg,XYear) =  OpTim (XReg,XYear) * OpEn / (OpEn - 1)
    	  OpTim (XReg,XYear) =  sqrt (OpTim (XReg,XYear))
    	  RegYearValidN      =  RegYearValidN + 1
    	end if
      end if  
    end if
    
  end do
end do

print*, "  > Total region-years valid: ", RegYearValidN
print*, "  > Enter the .tim file title to save: "
do
	read (*,*,iostat=ReadStatus), SaveTitle
	if (ReadStatus.LE.0 .AND. SaveTitle.NE."") exit
end do

call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpTim)

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

end subroutine StDevTim

!*******************************************************************************
! load scalar from .lin

subroutine LoadScalar

allocate (Scalar (UseTimN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadPattScal: Allocation failure #####"
Scalar = MissVal

print*, "  > Load the scalar for each .tim selected."
do XTim = 1, UseTimN
 print*, "  > Load the scalar for .tim number: ", XTim
 call LoadLin (WorkYearN, WorkADYear, FileLineN, FileLineNames, FileData)

 if (FileLineN.EQ.1) then
  SelectedLine = 1
 else
  print*, "  > Enter the index of the line to adopt as the scalar (0=list):"
  do
    do
	read (*,*,iostat=ReadStatus), SelectedLine
	if (ReadStatus.LE.0 .AND. SelectedLine.GE.0 .AND. SelectedLine .LE. FileLineN) exit
    end do
  
    if (SelectedLine.EQ.0) then
      do XLine = 1, FileLineN
        print*, XLine, FileLineNames(XLine)
      end do  
    end if
  
    if (SelectedLine.GT.0) exit
  end do    
 end if

 do XYear = 1, WorkYearN
    Scalar (XTim,XYear) = FileData (SelectedLine,XYear)
 end do

 deallocate (FileLineNames, FileData, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadPattScal: Deallocation failure #####"
end do

end subroutine LoadScalar

!*******************************************************************************
! scale and calc error

subroutine ScalingError

allocate (OpError (WorkRegN,WorkYearN), &
	  OpSampN (WorkRegN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ScalingError: Allocation failure #####"
OpError = MissVal
OpSampN = MissVal

MemThresh     = (100.0 - MissAccept) * UseTimN / 100.0

do XReg = 1, WorkRegN
 if (Pattern(XReg).NE.MissVal) then
  do XYear = 1, WorkYearN
    OpTot = 0.0
    OpEn  = 0.0
    
    do XTim = 1, UseTimN
      if (TimLoaded(UseTim(XTim),XReg,XYear).NE.MissVal.AND.Scalar(XTim,XYear).NE.MissVal) then
        OpTot = OpTot + (Scalar(XTim,XYear) * Pattern(XReg)) - TimLoaded(UseTim(XTim),XReg,XYear) 
        OpEn  = OpEn  + 1 
      end if
    end do
    
    if (OpEn.GE.MemThresh) then
        OpError (XReg,XYear) = OpTot / OpEn
        OpSampN (XReg,XYear) = OpEn
    end if
  end do
 end if
end do

end subroutine ScalingError

!*******************************************************************************
! calc sig of scaling error

subroutine ScalingErrorSig

print*, "  > Load stdev from .tim"
call LoadTim (WorkYearN, WorkADYear, WorkRegN, FileLineNames, OpStDev)

print*, "  > We divide the stdev by sqrt(n), since error = ensemble mean"
do XReg = 1, WorkRegN		
  do XYear = 1, WorkYearN
    if (OpStDev(XReg,XYear).NE.MissVal.AND.OpSampN(XReg,XYear).NE.MissVal) then
      OpStDev (XReg,XYear) = OpStDev (XReg,XYear)  / sqrt (OpSampN (XReg,XYear) )
    else
      OpStDev (XReg,XYear) = MissVal
    end if
  end do
end do

print*, "  > We divide the error by the adjusted stdev"
do XReg = 1, WorkRegN		
  do XYear = 1, WorkYearN
    if (OpStDev(XReg,XYear).NE.MissVal.AND.OpError(XReg,XYear).NE.MissVal) then
      OpError (XReg,XYear) = OpError (XReg,XYear) / OpStDev (XReg,XYear)
    else
      OpError (XReg,XYear) = MissVal
    end if
  end do
end do

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

print*, "  > Enter the .tim file title to save significances: "
do
	read (*,*,iostat=ReadStatus), SaveTitle
	if (ReadStatus.LE.0 .AND. SaveTitle.NE."") exit
end do
call SaveTim (WorkRegN, WorkYearN, SaveTitle, WorkRegNames, WorkADYear, OpError)

end subroutine ScalingErrorSig

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

end program OperateTim