! seqtim.f90
! f90 main program written on 19.06.00 by Tim Mitchell
! last modification on 26.09.00
! f90 -o ./../goglo/seqtim initialmod.f90 loadmod.f90 savemod.f90 transformmod.f90 sortmod.f90 
!        basicfun.f90 slicestats.f90 scalemod.f90 trendsmod.f90 arraymod.f90 ./../goglo/seqtim.f90
! purpose is to permit us to avoid multiple saves to .tim and truncations to zero   

program SeqTim

use InitialMod
use LoadMod
use SaveMod
use TransformMod
use SortMod
use BasicFun
use SliceStats
use ScaleMod
use TrendsMod
use ArrayMod

implicit none

real, pointer, dimension (:,:,:)	:: TimLoaded
real, pointer, dimension (:,:)		:: TimSaved, ThisTim, LinLoaded, LinSaved
real, pointer, dimension (:,:)		:: ExeSeries, WyeSeries, OperationKey, OpTim
real, pointer, dimension (:)		:: Vectorised, GloLoaded, FunData, OperationRes, YearVec
real, pointer, dimension (:)		:: AyeKay, BeeKay, AreKay, LoPass, HiPass
real, pointer, dimension (:)		:: ExeLine, WyeLine, WyeFit, Bases, Chunks
integer, pointer, dimension (:) 	:: WorkMapRawReg, WorkRegSizes, WorkADYear, UseReg
integer, pointer, dimension (:) 	:: NewMapRawReg, NewRegSizes
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw, WorkMapIDLReg, NewMapIDLReg

character (len=20), pointer, dimension (:) 	:: WorkRegNames, NewRegNames
character (len=80), pointer, dimension (:) 	:: ThisNames

real, allocatable, dimension (:,:)	:: NewTot, NewEn

integer, allocatable, dimension (:)	:: UseTim, RegSame, YearSame

real, parameter :: MissVal = -999.0

real :: MissAccept
real :: RunTot, RunEn, RunSqTot
real :: Aye, Bee, Are, Num, Den
real :: SumX, SumY, SumXX, SumYY, SumXY, En
real :: OpTot, OpSqTot, OpEn, OpMean, OpStDev
real :: RegSqTot, RegTot, RegEn, RegMean, RegStDev, RegRMSE
real :: CutOffVal, TestVal, ChosenKay, RealConstant
real :: RegThresh
real :: ValueA, ValueB

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN
integer :: WorkMonth0,WorkMonth1,WorkMonthN,WorkYearN,WorkDecN
integer :: WorkRegN, WorkTimN, ThisRegN, NewRegN
integer :: SelectTimN, ChosenTim, BlankColN, AllColUsed, LastColUsed
integer :: AllocStat, ReadStatus, ValidTot
integer :: XReg, XGlo, XTimFree, XTimSave, XSelect, XTim, XYear, XLine, XOperation, XFile, XLong, XLat
integer :: QAnom, QMain, QSampPop, QConBee, QColEns, QStatistic, QKayTot, QConvert, QWeight, QKeepMiss
integer :: QDeTrend
integer :: CutOffSign, CutOffType, RegKept, RegCutOff
integer :: ChosenExp, Operator, IntConstant
integer :: SliceLen, OperationN
integer :: ThisLinN, XSelectLin
integer :: SelectReg, SelectRegN, SelTim
integer :: EnsRegSame, EnsYearSame
integer :: GloFileN, LinFileN
integer :: Iter, Log, PerLen, GapLen, Offset, Year0, Year1

character (len=10) :: WorkGridTitle
character (len=40) :: WorkRegTitle, NewRegTitle
character (len=80) :: WorkGridFilePath, Blank

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

open (99,file="/cru/u2/f709762/data/scratch/log-seqtim.dat",status="replace",action="write")

print*, "  > ***** SeqTim *****"
print*

Blank      = ""
MissAccept = 10
AllColUsed = 0

call GridSelect  (WorkGrid,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call RegSelect   (WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkMapIDLReg,WorkRegSizes,WorkRegNames,&
		  WorkRegTitle,WorkRegN)
call PeriodSelect(WorkYearN,WorkDecN,WorkADYear)
call GrabTim

QMain = -1

do
  if     (QMain.EQ.-1) then
  else if (QMain.EQ.0) then
    call SetMissAccept (MissAccept)
  else if (QMain.EQ.1) then
    call SelectTimOpN
    call SelectTimOp
    call CalcMean
  else if (QMain.EQ.2) then
    call SelectTimOpN
    call SelectTimOp
    call CalcStDev
  else if (QMain.EQ.3) then
    SelectTimN = 2
    call SelectTimOp
    call CalcLSRReg
  else if (QMain.EQ.4) then
    SelectTimN = 1
    call SelectTimOp
    call DeTrendA
  else if (QMain.EQ.5) then
    SelectTimN = 1
    call SelectTimOp
    call GlobalMean
  else if (QMain.EQ.6) then
    SelectTimN = 3
    call SelectTimOp
    call PosABC
  else if (QMain.EQ.7) then
    SelectTimN = 2
    call SelectTimOp
    call MaskAB
  else if (QMain.EQ.8) then
    SelectTimN = 1
    call SelectTimOp
    call OperateAConstant
  else if (QMain.EQ.9) then
    SelectTimN = 2
    call SelectTimOp
    call OperateAB
  else if (QMain.EQ.10) then
    call ClearLastCol
  else if (QMain.EQ.11) then
    call GrabGlo
  else if (QMain.EQ.12) then
    call GrabLin
  else if (QMain.EQ.13) then
    SelectTimN = 1
    call SelectTimOp
    call ExtractTrendA
  else if (QMain.EQ.14) then
    SelectTimN = 1
    call SelectTimOp
    call StDevChunks
  else if (QMain.EQ.20) then
    SelectTimN = 1
    call SelectTimOp
    call SaveToTim
  else if (QMain.EQ.21) then
    SelectTimN = 1
    call SelectTimOp
    call SaveSlices
  else if (QMain.EQ.22) then
    SelectTimN = 1
    call SelectTimOp
    call SaveRegions
  else if (QMain.EQ.23) then
    SelectTimN = 1
    call SelectTimOp
    call SaveGlobalMean
  else if (QMain.EQ.24) then
    SelectTimN = 1
    call SelectTimOp
    call ConvertA
  else if (QMain.EQ.30) then
    call PrintToScreen
  else if (QMain.EQ.40) then
    call SequenceOp
  else if (QMain.EQ.50) then
    call SelectTimOpN
    call SelectTimOp
    call SlicePerMom
  else if (QMain.EQ.51) then
    SelectTimN = 1
    call SelectTimOp
    call GaussSmoothA
  else if (QMain.EQ.52) then
    SelectTimN = 1
    call SelectTimOp
    call AnomAPerA
  else
    print*, "  > Select the comparison to make: "
    print*, "  >   0. Change % acceptable missing, currently: ", MissAccept
    print*, "  >   1. Mean  of multiple A"
    print*, "  >   2. Stdev of multiple A"
    print*, "  >   3. LSR of each region in A(x) with B(y)"
    print*, "  >   4. Detrend A"
    print*, "  >   5. Global-mean of A, filling col"
    print*, "  >   6. Position of A relative to B, C"
    print*, "  >   7. Mask A using B"
    print*, "  >   8. Operate on A with k"
    print*, "  >   9. Operate on A with B"
    print*, "  >  10. Clear last column filled"
    print*, "  >  11. Fill columns with .glo files"
    print*, "  >  12. Fill columns with .lin files"
    print*, "  >  13. Robust trend for each region in A"
    print*, "  >  14. Calc stdev of chunks from A -> .glo"
    print*, "  >  20. Save A to .tim"
    print*, "  >  21. Save year slices from A to .glo"
    print*, "  >  22. Save regions from A to .lin"
    print*, "  >  23. Save global mean from A to .lin"
    print*, "  >  24. Convert A to another region set, to .tim"
    print*, "  >  30. Print stats for all to screen"
    print*, "  >  40. Perform sequence of operations"
    print*, "  >  50. Sliceify on A,B,... for per k (any moment)"
    print*, "  >  51. Smooth A with Gaussian filter"
    print*, "  >  52. Anomalise A using period from A"
    print*, "  >  99. Exit."
  end if
  
  if (allocated(UseTim)) deallocate (UseTim, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Main: Deallocation failure #####"
  
  print*
  print*, "  > Main menu: select option: "
  do
	read (*,*,iostat=ReadStatus), QMain
	if (ReadStatus.LE.0) exit
  end do
  
  if (QMain.EQ.99) exit
end do

deallocate (TimLoaded,RegSame,YearSame)
deallocate (WorkMapIDLReg,WorkRegSizes,WorkRegNames)

print*

close (99)

contains

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

subroutine GrabTim

BlankColN = 10

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

print*,  "  > Specify the number of blank columns: "
do
   read (*,*,iostat=ReadStatus), BlankColN
   if (ReadStatus.LE.0.AND.BlankColN.GE.0) exit
end do

XTimFree = WorkTimN + 1
WorkTimN = WorkTimN + BlankColN

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

do XTim = 1, (WorkTimN-BlankColN)
  print*, "  > Load .tim number ", XTim
  
  call LoadTim (WorkYearN,WorkADYear,ThisRegN,ThisNames,ThisTim)

  if (ThisRegN.NE.WorkRegN) print*, "  > ##### ERROR: GrabTim: incorrect no. of regions in file #####"
  
  do XYear = 1, WorkYearN
   do XReg = 1, WorkRegN
    TimLoaded (XTim,XReg,XYear) = ThisTim (XReg,XYear)
   end do
  end do
  
  deallocate (ThisTim,ThisNames, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabTim: Deallocation failure #####"
end do

end subroutine GrabTim

!*******************************************************************************
! select number of .glo

subroutine SelectTimOpN

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

end subroutine SelectTimOpN

!*******************************************************************************
! select particular .tim to use

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

if (AllColUsed.EQ.0) then
  LastColUsed = XTimFree - 1
else
  LastColUsed = XTimFree
end if

print*, "  > Select particular .tim to use in op, totalling: ", SelectTimN
do XSelect = 1, SelectTimN
    print*, "  > Select .tim number: ", XSelect
    do
	read (*,*,iostat=ReadStatus), ChosenTim
	if (ReadStatus.LE.0.AND.ChosenTim.GE.1.AND.ChosenTim.LE.LastColUsed) exit
    end do
    UseTim (XSelect) = ChosenTim
end do  

end subroutine SelectTimOp

!*******************************************************************************
! calc Mean

subroutine CalcMean
   
print*, "  > Calculating..."
  
allocate (FunData(SelectTimN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcMean: Allocation failure #####"

do XReg = 1, WorkRegN    
  do XYear = 1, WorkYearN
    
    do XTim = 1, SelectTimN
      FunData (XTim) = TimLoaded(UseTim(XTim),XReg,XYear)
    end do

    TimLoaded (XTimFree,XReg,XYear) = OpCalcMean (FunData,SelectTimN,MissAccept)

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

call AlertMiss

end subroutine CalcMean

!*******************************************************************************
! calc StDev

subroutine CalcStDev
   
allocate (FunData(SelectTimN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcStDev: Allocation failure #####"

if (WorkTimN.GE.2) then
 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..."  

 do XReg = 1, WorkRegN    
  do XYear = 1, WorkYearN
    
    do XTim = 1, SelectTimN
      FunData (XTim) = TimLoaded(UseTim(XTim),XReg,XYear)
    end do

    TimLoaded (XTimFree,XReg,XYear) = OpCalcStDev (FunData,SelectTimN,MissAccept,QSampPop)

  end do
 end do
  
 call AlertMiss
else
 print*, "  > Insufficient .tim selected."
end if

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

end subroutine CalcStDev

!*******************************************************************************
! calc LSR between regions in A(x) and B(y)

subroutine CalcLSRReg
   
print*, "  > Regress for y=ax (=1) or y=ax+b (=2) ? "
do
	read (*,*,iostat=ReadStatus), QKayTot
	if (ReadStatus.LE.0.AND.QKayTot.GE.1.AND.QKayTot.LE.2) exit
end do

allocate (WyeSeries(1,WorkYearN), &
	  ExeSeries(1,WorkYearN), &
	  AyeKay   (WorkRegN),  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcLSRReg: Allocation failure #####"

if (QKayTot.Eq.2) then
  allocate (BeeKay (WorkRegN), &
  	    AreKay (WorkRegN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcLSRReg: Allocation failure #####"
end if

print*, "  > Calculating..."
  
do XReg = 1, WorkRegN    
  do XYear = 1, WorkYearN
    ExeSeries (1,XYear) = TimLoaded(UseTim(1),XReg,XYear)
    WyeSeries (1,XYear) = TimLoaded(UseTim(2),XReg,XYear)
  end do
  
  if (QKayTot.EQ.1) then
    call LinearLSRZero (1,WorkYearN,ExeSeries,WyeSeries,Aye)
    AyeKay (XReg) = Aye
  else
    call LinearLSR (1,WorkYearN,ExeSeries,WyeSeries,Aye,Bee,Are)
    AyeKay (XReg) = Aye
    BeeKay (XReg) = Bee
    AreKay (XReg) = Are
  end if
end do

print*, "  > Save 'a'."
call SaveGlo (WorkLongN, WorkLatN, WorkRegN, WorkGridFilePath, Blank, Blank, AyeKay, WorkMapIDLReg)

if (QKayTot.Eq.2) then
  print*, "  > Save 'b'."
  call SaveGlo (WorkLongN, WorkLatN, WorkRegN, WorkGridFilePath, Blank, Blank, BeeKay, WorkMapIDLReg)
  print*, "  > Save 'r'."
  call SaveGlo (WorkLongN, WorkLatN, WorkRegN, WorkGridFilePath, Blank, Blank, AreKay, WorkMapIDLReg)
end if

if (QKayTot.Eq.2) then
  deallocate (BeeKay,AreKay, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcLSRReg: Deallocation failure #####"
end if

deallocate (WyeSeries,ExeSeries,AyeKay, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcLSRReg: Deallocation failure #####"

end subroutine CalcLSRReg

!*******************************************************************************
! calc global mean and fill column with it

subroutine DeTrendA
   
allocate (ThisTim(WorkRegN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DeTrendA: Allocation failure #####"

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

call DeTrendTim (WorkRegN,WorkYearN,MissAccept,1,WorkADYear,ThisTim,AyeKay,BeeKay)

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    TimLoaded(XTimFree,XReg,XYear) = ThisTim (XReg,XYear)
  end do
end do

call AlertMiss

deallocate (ThisTim, AyeKay, BeeKay, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DeTrendA: Deallocation failure #####"

end subroutine DeTrendA

!*******************************************************************************
! extract trend (y=ax+b) from each region
! uses robust method of LSR

subroutine ExtractTrendA
   
allocate (ExeLine(WorkYearN), &
	  WyeLine(WorkYearN), &
	  AyeKay  (WorkRegN), &
	  BeeKay  (WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ExtractTrendA: Allocation failure #####"

do XYear = 1, WorkYearN					! get years
    ExeLine (XYear) = real (WorkADYear (XYear))
end do

print*, "  > Calculating..."

do XReg = 1, WorkRegN						! loop by region
    do XYear = 1, WorkYearN					! get region time Line
      WyeLine (XYear) = TimLoaded(UseTim(1),XReg,XYear)
    end do
    
    call RobustLSR (WorkYearN,ExeLine,WyeLine,WyeFit,Bee,Aye,Iter)	! reversal of Aye, Bee fits routine
    
    deallocate (WyeFit, stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: ExtractTrendA: Deallocation failure #####"

    AyeKay(XReg) = Aye
    BeeKay(XReg) = Bee
end do

print*, "  > Save 'a' from y=ax+b to .glo."
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,AyeKay,WorkMapIDLReg)

print*, "  > Save 'b' from y=ax+b to .glo."
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,BeeKay,WorkMapIDLReg)

deallocate (ExeLine,WyeLine,AyeKay,BeeKay, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ExtractTrendA: Deallocation failure #####"

end subroutine ExtractTrendA

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

subroutine StDevChunks

print*, "  > Enter the offset from the start to omit (>=0): "
do
	read (*,*,iostat=ReadStatus), Offset
	if (ReadStatus.LE.0.AND.Offset.GE.0.AND.Offset.LT.WorkYearN) exit
end do

print*, "  > Enter the length of chunk and of gap: "
do
	read (*,*,iostat=ReadStatus), PerLen, GapLen
	if (ReadStatus.LE.0.AND.PerLen.GT.0.AND.PerLen.LE.WorkYearN.AND.GapLen.GE.0) exit
end do

print*, "  > Calculate 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

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

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    Vectorised (XYear) = TimLoaded(UseTim(1),XReg,XYear)
  end do
  
  call GetChunkMeans (Vectorised,Chunks,Offset,PerLen,GapLen,MissAccept)
  
  GloLoaded(XReg) = OpCalcStDev (Chunks, nint(MissVal), MissAccept, QSampPop)
  
  deallocate (Chunks,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: StDevChunks: Deallocation failure #####"
end do

call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,GloLoaded,WorkMapIDLReg)

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

end subroutine StDevChunks

!*******************************************************************************
! calc global mean and fill column with it

subroutine GlobalMean
   
print*, "  > Calculate unweighted (=1) or weighted (=2) mean ? "
do
	read (*,*,iostat=ReadStatus), QWeight
	if (ReadStatus.LE.0.AND.QWeight.GE.1.AND.QWeight.LE.2) exit
end do

if (QWeight.EQ.2) then
  print*, "  > Load .glo file containing weights: "
  call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, GloLoaded)  
end if

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

print*, "  > Calculating..."
  
do XYear = 1, WorkYearN
  if (YearSame(UseTim(1)).NE.1.OR.XYear.EQ.1) then
   do XReg = 1, WorkRegN    
    FunData (XReg) = TimLoaded(UseTim(1),XReg,XYear)
   end do
  
   if (QWeight.EQ.1) OpMean = OpCalcMean     (FunData,WorkRegN,MissAccept)
   if (QWeight.EQ.2) OpMean = OpWeightedMean (FunData,GloLoaded,WorkRegN,MissAccept)
  end if
  
  do XReg = 1, WorkRegN    
   TimLoaded (XTimFree,XReg,XYear) = OpMean
  end do
end do

RegSame  (XTimFree) = 1
YearSame (XTimFree) = 0
 
if (QWeight.EQ.2) then
  deallocate (GloLoaded, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GlobalMean: Deallocation failure #####"
end if

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

call AlertMiss

end subroutine GlobalMean

!*******************************************************************************
! pos of A relative to B and C

subroutine PosABC

print*, "  > Calculating..."
  
do XReg = 1, WorkRegN    
 do XYear = 1, WorkYearN   
   if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal.AND.TimLoaded(UseTim(2),XReg,XYear).NE.MissVal &
   			.AND.TimLoaded(UseTim(3),XReg,XYear).NE.MissVal) then
     if (TimLoaded(UseTim(1),XReg,XYear).LT.TimLoaded(UseTim(2),XReg,XYear) &
     			.AND.TimLoaded(UseTim(1),XReg,XYear).LT.TimLoaded(UseTim(3),XReg,XYear)) then
       TimLoaded (XTimFree,XReg,XYear) = -1
     else if (TimLoaded(1,XReg,XYear).GT.TimLoaded(UseTim(2),XReg,XYear) &
     			.AND.TimLoaded(UseTim(1),XReg,XYear).GT.TimLoaded(UseTim(3),XReg,XYear)) then
       TimLoaded (XTimFree,XReg,XYear) = 1
     else
       TimLoaded (XTimFree,XReg,XYear) = 0
     end if
   else
       TimLoaded (XTimFree,XReg,XYear) = MissVal
   end if
 end do
end do

call AlertMiss

end subroutine PosABC

!*******************************************************************************
! mask A using B

subroutine MaskAB
   
 print*, "  > Is the threshold a real (=1) or a magnitude (=2): "
 do
	read (*,*,iostat=ReadStatus), CutOffType
	if (ReadStatus.LE.0.AND.CutOffType.GE.1.AND.CutOffType.LE.2) exit
 end do

 print*, "  > Enter the threshold in B beyond which to set A to missing: "
 do
	read (*,*,iostat=ReadStatus), CutOffVal
	if (ReadStatus.LE.0) exit
 end do

 print*, "  > Cut off values smaller (=1) or greater (=2) than: ", CutOffVal
 do
	read (*,*,iostat=ReadStatus), CutOffSign
	if (ReadStatus.LE.0.AND.CutOffSign.GE.1.AND.CutOffSign.LE.2) exit
 end do
 
print*, "  > Calculating..."
  
 RegKept = 0
 RegCutOff = 0
 
 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 (CutOffType.EQ.2) then
       		TestVal = abs (TimLoaded(UseTim(2),XReg,XYear))
       else
       		TestVal = TimLoaded(UseTim(2),XReg,XYear)
       end if
       
       if      (TestVal.LT.CutOffVal.AND.CutOffSign.EQ.1) then
       	  TimLoaded (XTimFree,XReg,XYear) = MissVal
       	  RegCutOff = RegCutOff + 1	
       else if (TestVal.GT.CutOffVal.AND.CutOffSign.EQ.2) then
          TimLoaded (XTimFree,XReg,XYear) = MissVal
       	  RegCutOff = RegCutOff + 1	
       else
          TimLoaded (XTimFree,XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear)
       	  RegKept = RegKept + 1	
       end if
   else
          TimLoaded (XTimFree,XReg,XYear) = MissVal
   end if
  end do
 end do
 
print*, "  > Region-years kept, cut off: ", RegKept, RegCutOff
 
call AlertMiss

end subroutine MaskAB

!*******************************************************************************
! divide by constant

subroutine OperateAConstant

call OpAKayOptions

print*, "  > Calculating..."
  
do XReg = 1, WorkRegN
 do XYear = 1, WorkYearN
    TimLoaded (XTimFree,XReg,XYear) = OpAwithB (TimLoaded(UseTim(1),XReg,XYear),RealConstant,Operator)        
 end do
end do

call AlertMiss

end subroutine OperateAConstant

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

subroutine OpAKayOptions

print*, "  > Divide =1, times =2, minus =3, add =4, sqroot =5, expon =6, abs =7"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.7) 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.5) then
  RealConstant = 0
else if (Operator.EQ.6) then
  print*, "  > Enter constant (integer):"
  do
	read (*,*,iostat=ReadStatus), RealConstant
	if (ReadStatus.LE.0) exit
  end do
else if (Operator.EQ.7) then
  RealConstant = 0
end if

end subroutine OpAKayOptions

!*******************************************************************************
! operate on A with B

subroutine OperateAB

call OpABOptions

print*, "  > Calculating..."
  
do XReg = 1, WorkRegN
 do XYear = 1, WorkYearN
    TimLoaded (XTimFree,XReg,XYear) = OpAwithB (TimLoaded(UseTim(1),XReg,XYear), &
    			TimLoaded(UseTim(2),XReg,XYear),Operator)        
 end do
end do

call AlertMiss

end subroutine OperateAB

!*******************************************************************************
! operate on A with B

subroutine OpABOptions

print*, "  > A operator B"
print*, "  > Operator is divide =1, multi =2, minus =3, add =4"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.4) exit
end do

end subroutine OpABOptions

!*******************************************************************************
! clear last column

subroutine ClearLastCol

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    TimLoaded ((XTimFree-1),XReg,XYear) = MissVal
  end do
end do

XTimFree = XTimFree - 1

print*, "  > Column cleared and made next free: ", XTimFree 

end subroutine ClearLastCol

!*******************************************************************************
! fill column with .glo

subroutine GrabGlo

print*, "  > Enter the number of .glo files to load: "
do
	read (*,*,iostat=ReadStatus), GloFileN
	if (ReadStatus.LE.0 .AND. GloFileN.GE.0) exit
end do

do XFile = 1, GloFileN
 print*, "  > Load .glo file: ", XFile
 call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, GloLoaded)

 do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    TimLoaded (XTimFree,XReg,XYear) = GloLoaded (XReg)
  end do
 end do

 YearSame (XTimFree) = 1		! denotes that entire column has identical years
 RegSame  (XTimFree) = 0

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

 call AlertMiss
end do

end subroutine GrabGlo

!*******************************************************************************
! fill column with .lin

subroutine GrabLin

print*, "  > Enter the number of .lin files to load: "
do
	read (*,*,iostat=ReadStatus), LinFileN
	if (ReadStatus.LE.0 .AND. LinFileN.GE.0) exit
end do

do XFile = 1, LinFileN
 print*, "  > Load .lin file: ", XFile
 call LoadLin (WorkYearN, WorkADYear, ThisLinN, ThisNames, LinLoaded)

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

 do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    TimLoaded (XTimFree,XReg,XYear) = LinLoaded (XSelectLin,XYear)
  end do
 end do

 RegSame  (XTimFree) = 1		! denotes that entire column has identical regions
 YearSame (XTimFree) = 0
 
 deallocate (ThisNames,LinLoaded, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabLin: Deallocation failure #####"

 call AlertMiss
end do

end subroutine GrabLin

!*******************************************************************************
! save .tim file

subroutine SaveToTim

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

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

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

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

end subroutine SaveToTim

!*******************************************************************************
! convert A to another region set and save to .tim
! WARNING: this will run, but not produce the right results, unless the old is one-region-per-box

subroutine ConvertA

call RegSelect (WorkGrid, WorkLongN, WorkLatN, WorkDataN, NewMapIDLReg, &
		NewRegSizes, NewRegNames, NewRegTitle, NewRegN)

print*, "  > Convert as totals (=1) or means (=2) ? "
do
	read (*,*,iostat=ReadStatus), QConvert
	if (ReadStatus.LE.0.AND.QConvert.GE.1.AND.QConvert.LE.2) exit
end do

print*, "  > Converting..."

allocate (NewTot   (NewRegN,WorkYearN), &
	  NewEn    (NewRegN,WorkYearN), &
	  TimSaved (NewRegN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertA: Allocation failure #####"
NewTot   = 0.0
NewEn    = 0.0
TimSaved = MissVal

do XLong = 1, WorkLongN
  do XLat = 1, WorkLatN
    if (NewMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is wanted for new region set    
      if (WorkMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is specified in old region set    	
    	
    	do XYear = 1, WorkYearN
    	  if (TimLoaded (UseTim(1),WorkMapIDLReg(XLong,XLat),XYear).NE.MissVal) then	! value is valid    	
    	    NewTot (NewMapIDLReg(XLong,XLat),XYear) = NewTot (NewMapIDLReg(XLong,XLat),XYear) + &
    	  						TimLoaded (UseTim(1),WorkMapIDLReg(XLong,XLat),XYear)
            NewEn  (NewMapIDLReg(XLong,XLat),XYear) = NewEn  (NewMapIDLReg(XLong,XLat),XYear) + 1.0
          end if
        end do
        
      end if    
    end if    
  end do
end do

do XReg = 1, NewRegN
 do XYear = 1, WorkYearN
  if (NewEn(XReg,XYear).EQ.NewRegSizes(XReg)) then
    if (QConvert.EQ.1) then
      TimSaved (XReg,XYear) = NewTot (XReg,XYear)
    else
      TimSaved (XReg,XYear) = NewTot (XReg,XYear) / NewEn (XReg,XYear)
    end if
  end if
 end do
end do

call SaveTim (NewRegN, WorkYearN, Blank, Blank, NewRegNames, WorkADYear, TimSaved)

deallocate (NewMapIDLReg, NewRegSizes, NewRegNames, NewTot, NewEn, TimSaved, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertA: Deallocation failure #####"

end subroutine ConvertA

!*******************************************************************************
! save specified slices to .glo

subroutine SaveSlices

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, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertGlo: Deallocation failure #####"

end subroutine SaveSlices

!*******************************************************************************
! save specified regions to .lin

subroutine SaveRegions

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

if (WorkRegN.GT.1) then
 SelectRegN = 0
 do
  print*, "  > Specify a region to save (0=list,-1=end): "
  do
    read (*,*,iostat=ReadStatus), SelectReg
    if (ReadStatus.NE.0.OR.SelectReg.LT.-1.OR.SelectReg.GT.WorkRegN) print*, "  > Unacceptable entry."
    if (ReadStatus.LE.0.AND.SelectReg.GE.-1.AND.SelectReg.LE.WorkRegN) exit
  end do
  
  if (SelectReg.EQ.0) then
    do XReg = 1, WorkRegN
      print "(i6,a1,a20)", XReg, " ", WorkRegNames(XReg)
    end do
  else if (SelectReg.GE.1.AND.SelectReg.LE.WorkRegN) then
    SelectRegN          = SelectRegN + 1
    UseReg (SelectRegN) = SelectReg
    print "(i6,a1,a20)", SelectReg, " ", WorkRegNames(SelectReg)
  else if (SelectReg.NE.-1) then
    print*, "  > Unacceptable entry."
  end if
  
  if (SelectReg.EQ.-1.AND.SelectRegN.GE.1) exit
 end do
else
 SelectRegN = 1
 UseReg(1)  = 1
end if

allocate (LinSaved (SelectRegN,WorkYearN), 	&
	  ThisNames(SelectRegN),		stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveRegions: Allocation failure #####"
LinSaved = MissVal

do XReg = 1, SelectRegN
  ThisNames (XReg) = WorkRegNames (UseReg(XReg))
end do

do XReg = 1, SelectRegN
  do XYear = 1, WorkYearN
    LinSaved (XReg,XYear) = TimLoaded (UseTim(1),UseReg(XReg),XYear)
  end do
end do

call SaveLin (SelectRegN, WorkYearN, ThisNames, WorkADYear, LinSaved)

deallocate (LinSaved,ThisNames,UseReg, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveRegions: Deallocation failure #####"

end subroutine SaveRegions

!*******************************************************************************
! save global mean to .lin

subroutine SaveGlobalMean

print*, "  > Calculate unweighted (=1) or weighted (=2) mean ? "
do
	read (*,*,iostat=ReadStatus), QWeight
	if (ReadStatus.LE.0.AND.QWeight.GE.1.AND.QWeight.LE.2) exit
end do

if (QWeight.EQ.2) then
  print*, "  > Load .glo file containing weights: "
  call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, GloLoaded)  
end if

allocate (LinSaved (1,WorkYearN), 	&
	  ThisNames(1),			&
	  FunData(WorkRegN),		stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGlobalMean: Allocation failure #####"
LinSaved     = MissVal
ThisNames(1) = "global mean"

do XYear = 1, WorkYearN
   do XReg = 1, WorkRegN    
    FunData (XReg) = TimLoaded(UseTim(1),XReg,XYear)
   end do
   
   if (QWeight.EQ.1) LinSaved (1,XYear) = OpCalcMean     (FunData,          WorkRegN,MissAccept)
   if (QWeight.EQ.2) LinSaved (1,XYear) = OpWeightedMean (FunData,GloLoaded,WorkRegN,MissAccept)
end do

call SaveLin (1, WorkYearN, ThisNames, WorkADYear, LinSaved)

if (QWeight.EQ.2) then
  deallocate (GloLoaded, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGlobalMean: Deallocation failure #####"
end if

deallocate (LinSaved,ThisNames,FunData, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveGlobalMean: Deallocation failure #####"

end subroutine SaveGlobalMean

!*******************************************************************************
! print details to screen

subroutine PrintToScreen

print "(a4,5a12)", " Col", "   Year-Reg.", "         Sum", "  Sum Squares", "        Mean", "       StDev" 

do XTim = 1, (XTimFree-1)
  OpTot   = 0.0
  OpSqTot = 0.0
  OpEn    = 0.0
  OpMean  = MissVal
  OpStDev = MissVal
  
  do XYear = 1, WorkYearN
   do XReg = 1, WorkRegN
    if (TimLoaded(XTim,XReg,XYear).NE.MissVal) then
      OpTot   = OpTot   +  TimLoaded(XTim,XReg,XYear)
      OpSqTot = OpSqTot + (TimLoaded(XTim,XReg,XYear) ** 2)
      OpEn    = OpEn    +  1.0
    end if
   end do
  end do
  
  if (OpEn.GE.1) then
  	OpMean  =  OpTot / OpEn
  	OpStDev = (OpSqTot / OpEn) - (OpMean ** 2)
  	OpStDev = sqrt (OpStDev)
  end if
  
  print "(i4,5e12.4)", XTim, OpEn, OpTot, OpSqTot, OpMean, OpStDev
end do

end subroutine PrintToScreen

!*******************************************************************************
! sequence of operations

subroutine SequenceOp

print*, "  > Enter the number of operations: "
do
	read (*,*,iostat=ReadStatus), OperationN
	if (ReadStatus.LE.0.AND.OperationN.GE.1) exit
end do

allocate (OperationKey (OperationN,4), 	&
	  OperationRes (OperationN),    stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SequenceOp: Allocation failure #####"

print*, "  > DEFINE OPERATIONS."
do XOperation = 1, OperationN
  print*, "  > Operation: ", XOperation
  print*, "  > Operate on A with constant (=8) or B (=9): "
  do
	read (*,*,iostat=ReadStatus), QConBee
	if (ReadStatus.LE.0.AND.QConBee.GE.8.AND.QConBee.LE.9) exit
  end do
  
  if      (QConBee.EQ.8) then
  	OperationKey (XOperation,1) = 8

  	print*, "  > Select particular .tim A from either/or: "
  	call SeqSelTim
  	OperationKey (XOperation,2) = SelTim

  	call OpAKayOptions
  	OperationKey (XOperation,3) = RealConstant
  	OperationKey (XOperation,4) = Operator
  else if (QConBee.EQ.9) then
  	OperationKey (XOperation,1) = 9

  	print*, "  > Select particular .tim A from either/or: "
  	call SeqSelTim
  	OperationKey (XOperation,2) = SelTim

  	call OpABOptions
  	OperationKey (XOperation,4) = Operator

  	print*, "  > Select particular .tim B from either/or: "
  	call SeqSelTim
  	OperationKey (XOperation,3) = SelTim
  end if  
end do

print*, "  > Calculating..."

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    OperationRes = MissVal
    
    do XOperation = 1, OperationN
      if      (OperationKey(XOperation,2).LE.WorkTimN) then				! value A
      	ValueA = TimLoaded (OperationKey(XOperation,2),XReg,XYear)
      else
      	ValueA = OperationRes ((OperationKey(XOperation,2)-WorkTimN))
      end if
      
      if (OperationKey(XOperation,1).EQ.9) then						! value B or constant
       if      (OperationKey(XOperation,3).LE.WorkTimN) then
       	ValueB = TimLoaded (OperationKey(XOperation,3),XReg,XYear)
       else
        ValueB = OperationRes ((OperationKey(XOperation,3)-WorkTimN))
       end if
      else
        ValueB = OperationKey(XOperation,3)
      end if
      
      Operator = nint(OperationKey(XOperation,4))
      OperationRes (XOperation) = OpAwithB (ValueA,ValueB,Operator)			! operation
    end do
    
    TimLoaded (XTimFree,XReg,XYear) = OperationRes (OperationN)  			! final value  
  end do
end do

print*, "  > We have stored the result in column: ", XTimFree
if (XTimFree.EQ.WorkTimN) AllColUsed = 1
if (XTimFree.LT.WorkTimN) XTimFree = XTimFree + 1

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

end subroutine SequenceOp

!*******************************************************************************
! sequence of operations

subroutine SeqSelTim

print*, "  >     .tim in main (usual no), this seq output (op no + ", WorkTimN
do
	read (*,*,iostat=ReadStatus), SelTim
	if (ReadStatus.LE.0.AND.SelTim.GE.1.AND.SelTim.LT.(WorkTimN+XOperation)) exit
end do
  
end subroutine SeqSelTim

!*******************************************************************************
! alert to all MissVal in latest column

subroutine AlertMiss

ValidTot = 0
XReg 	 = 0

do
  XReg  = XReg + 1
  XYear = 0
  
  do
    XYear = XYear + 1
    if (TimLoaded(XTimFree,XReg,XYear).NE.MissVal) ValidTot = ValidTot + 1
    if (ValidTot.GT.0) exit
    if (XYear.EQ.WorkYearN) exit
  end do
  
  if (ValidTot.GT.0) exit
  if (XReg.EQ.WorkRegN) exit
end do

if (ValidTot.GT.0) then
  print*, "  > We have filled column: ", XTimFree
  if (XTimFree.EQ.WorkTimN) AllColUsed = 1
  if (XTimFree.LT.WorkTimN) XTimFree = XTimFree + 1
else
  print*, "  > @@@@@ All column missing, so left unfilled. @@@@@"
  RegSame  (XTimFree) = 0
  YearSame (XTimFree) = 0
end if

end subroutine AlertMiss

!*******************************************************************************
! sliceify A,B,... by per k into any statistic
! median, mean, variance are calc as sample stats in usual manner
! skewness = (3*(mean--median))/stdev  (Francis p88): <0=left skew=neg skew >0=opposite 0=normal
! kurtosis = (fourth moment/(second moment**2))--3 (Francis p88) -3<sk<0=platykurtic sk>0=leptokurtic 0=normal

subroutine SlicePerMom

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

print*, "  > Enter the statistic to calc for each slice:"
print*, "  >     0=median,1=mean,2=variance,3=skewness,4=kurtosis"
do
	  read (*,*,iostat=ReadStatus), QStatistic
	  if (ReadStatus.LE.0.AND.QStatistic.GE.0.AND.QStatistic.LE.4) exit
end do

print*, "  > Detrend each slice before calculating moment ? (1=no,2=yes)"
do
	  read (*,*,iostat=ReadStatus), QDeTrend
	  if (ReadStatus.LE.0.AND.QDeTrend.GE.1.AND.QDeTrend.LE.2) exit
end do

if (SelectTimN.GT.1) then
  print*, "  > Calc for each column individually (=1), or as ensemble (=2) ?"
  do
	  read (*,*,iostat=ReadStatus), QColEns
	  if (ReadStatus.LE.0.AND.QColEns.GE.1.AND.QColEns.LE.2) exit
  end do
else
  QColEns = 1
end if

if (QColEns.EQ.1) call SlicePerCol
if (QColEns.EQ.2) call SlicePerEns

end subroutine SlicePerMom

!*******************************************************************************
! sliceify A,B,... by per k as columns

subroutine SlicePerCol

print*, "  > Sliceifying each column and region..."

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

do XTim = 1, SelectTimN
  do XReg = 1, WorkRegN
    if (XReg.EQ.1.OR.RegSame(UseTim(XTim)).EQ.0) then
     do XYear = 1, WorkYearN
      Vectorised (XYear) = TimLoaded (UseTim(XTim),XReg,XYear)
      YearVec    (XYear) = real (WorkADYear(XYear))
     end do
    
     if (QDeTrend.EQ.2) call DeTrendCol (YearVec,Vectorised)
     call SliceStatistics (WorkYearN,SliceLen,QStatistic,MissAccept,Vectorised)
     
    end if
    
    do XYear = 1, WorkYearN
      TimLoaded (XTimFree,XReg,XYear) = Vectorised (XYear)
    end do    
  end do
  
  print*, "  > Have sliced col: ", UseTim(XTim), " into: ", XTimFree
  
  if (RegSame(UseTim(XTim)).EQ.1) then
  	RegSame (XTimFree) = 1
  else
  	RegSame (XTimFree) = 0
  end if
  
  YearSame (XTimFree) = 0
 
  if (XTimFree.EQ.WorkTimN) AllColUsed = 1
  if (XTimFree.LT.WorkTimN) XTimFree = XTimFree + 1
end do

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

end subroutine SlicePerCol

!*******************************************************************************
! sliceify A,B,... by per k into stdev as ensemble

subroutine SlicePerEns

print*, "  > Sliceifying the ensemble for each region..."

allocate (OpTim (SelectTimN,WorkYearN), &
	  Vectorised (WorkYearN),	&
	  YearVec    (WorkYearN),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SlicePerEns: Allocation failure #####"

EnsRegSame = 1
do XTim = 1, SelectTimN
  if (RegSame(UseTim(XTim)).NE.1) EnsRegSame = EnsRegSame + 1
end do

do XReg = 1, WorkRegN
   if (XReg.EQ.1.OR.EnsRegSame.NE.1) then
    Vectorised = MissVal
    
    do XTim = 1, SelectTimN
     do XYear = 1, WorkYearN
      OpTim (XTim,XYear) = TimLoaded (UseTim(XTim),XReg,XYear)
      YearVec    (XYear) = real (WorkADYear(XYear))
     end do
    end do
    
    if (QDeTrend.EQ.2) call DeTrendCol (YearVec,Vectorised)
    call SliceStatisticsEns (SelectTimN,WorkYearN,SliceLen,QStatistic,MissAccept,OpTim,Vectorised)
   end if
    
   do XYear = 1, WorkYearN
      TimLoaded (XTimFree,XReg,XYear) = Vectorised (XYear)
   end do
end do
  
if (EnsRegSame.EQ.1) then
  	RegSame (XTimFree) = 1
else
  	RegSame (XTimFree) = 0
end if
  
YearSame (XTimFree) = 0
 
print*, "  > Have sliced ensemble into: ", XTimFree

if (XTimFree.EQ.WorkTimN) AllColUsed = 1
if (XTimFree.LT.WorkTimN) XTimFree = XTimFree + 1

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

end subroutine SlicePerEns

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

subroutine GaussSmoothA

print*, "  > Enter the time period (years) to smooth:"
do
	  read (*,*,iostat=ReadStatus), SliceLen
	  if (ReadStatus.LE.0.AND.SliceLen.GE.2.AND.SliceLen.LT.(WorkYearN/2)) exit
end do

print*, "  > Reinstate missing values? (1=no, 2=yes)"
do
	  read (*,*,iostat=ReadStatus), QKeepMiss
	  if (ReadStatus.LE.0.AND.QKeepMiss.GE.1.AND.QKeepMiss.LE.2) exit
end do

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

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    Vectorised (XYear) = TimLoaded (UseTim(1),XReg,XYear)
  end do
  
  call GaussSmooth (WorkYearN,SliceLen,QKeepMiss,Vectorised,LoPass,HiPass)
  
  do XYear = 1, WorkYearN
    TimLoaded (XTimFree,XReg,XYear) = LoPass (XYear)
  end do
end do

call AlertMiss

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

end subroutine GaussSmoothA

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

subroutine AnomAPerA

print*, "  > Enter the first, last years of base period:"
do
	  read (*,*,iostat=ReadStatus), Year0, Year1
	  if (ReadStatus.LE.0.AND.Year0.GE.1.AND.Year1.GE.Year0.AND.Year1.LE.WorkYearN) exit
end do
PerLen = Year1 - Year0 + 1

print*, "  > Calculate absolute (=1) or percentage (=2) anomalies ?"
do
	  read (*,*,iostat=ReadStatus), QAnom
	  if (ReadStatus.LE.0.AND.QAnom.GE.1.AND.QAnom.LE.2) exit
end do

allocate (Vectorised(PerLen), &
	  Bases(WorkRegN),    stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AnomAPerA: Allocation failure #####"

do XReg = 1, WorkRegN				! calc bases
  do XYear = 1, PerLen
    Vectorised(XYear) = TimLoaded(UseTim(1),XReg,(Year0+XYear-1))
  end do
  
  Bases(XReg) = OpCalcMean (Vectorised, PerLen, MissAccept)
end do

if (QAnom.EQ.1) then
 do XReg = 1, WorkRegN				! calc abs anomalies
  do XYear = 1, WorkYearN
    if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal.AND.Bases(XReg).NE.MissVal) &
    	TimLoaded(XTimFree,XReg,XYear) = TimLoaded(UseTim(1),XReg,XYear) - Bases(XReg)
  end do
 end do
else						! calc % anomalies
 do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (TimLoaded(UseTim(1),XReg,XYear).NE.MissVal.AND.Bases(XReg).NE.MissVal.AND.Bases(XReg).NE.0) &
        TimLoaded(XTimFree,XReg,XYear) = 100.0 * (TimLoaded(UseTim(1),XReg,XYear) - Bases(XReg)) / Bases(XReg)
  end do
 end do
end if

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

call AlertMiss

end subroutine AnomAPerA

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

end program SeqTim