! comparetim.f90
! f90 main program written on 08.02.00 by Tim Mitchell
! last modification on 28.04.00
! f90 -o comparetim initialmod.f90 loadmod.f90 savemod.f90 scalemod.f90 transformmod.f90 comparetim.f90

program CompareTim

use InitialMod
use LoadMod
use SaveMod
use ScaleMod
use TransformMod

implicit none

real, pointer, dimension (:)		:: VecSaved, LinVec, LoPass, HiPass 
real, pointer, dimension (:,:)		:: LinLoaded, FileLin, LinSaved, TimSaved, OpYear, OpData

integer, pointer, dimension (:) 	:: WorkADYear
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw

character (len=80), pointer, dimension (:)	:: FileLinNames, WorkLinNames, SaveLinNames

integer, allocatable, dimension (:)	:: LineSelect
integer, allocatable, dimension (:,:)	:: PairSelect

real, parameter :: MissVal = -999.0

real :: RunTot, RunEn, RunSqTot
real :: AllSumX, AllSumY, AllSumXX, AllSumYY, AllSumXY, AllEn
real :: Numer, Denom, AllAye, AllBee, AllAre
real :: KaySelect, RealConstant
real :: MissAccept
real :: OpAye, OpBee, OpPea

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN
integer :: WorkYearN, WorkDecN
integer :: WorkFileTimN, WorkFileLinN, WorkLinN, FileLinN, SaveTimN
integer :: AllocStat, ReadStatus
integer :: XFile, XLin, XYear, XYear0, XFind, XSelect
integer :: QAnom, QMain, QSuffix
integer :: SelectLin, SelectLinA, SelectLinB, SelectLinN, QIntercept
integer :: FirstFree, CurrentLin, SliceLen
integer :: BandSide, BandWidth
integer :: Operator, IntConstant

character (len=10) :: WorkGridTitle
character (len=80) :: WorkGridFilePath, WorkTimTitle, Blank

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

Blank = ""

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

call Initialise

QMain = -1
do
  QSuffix = 0
  
  if      (QMain.EQ.1) then
    QSuffix = 1
    call SelectTwo
    call TwoFromOne
  else if (QMain.EQ.2) then
    call SelectMany
    call SelectSuffix
    call CalcMean
  else if (QMain.EQ.3) then
    QSuffix = 1
    call SelectMany
    call CalcStDev
  else if (QMain.EQ.4) then
    QSuffix = 0
    call SelectMany
    call CalcBestFitTime
  else if (QMain.EQ.5) then
    QSuffix = 0
    call CalcBestFitLines
  else if (QMain.EQ.6) then
    QSuffix = 0
    call SelectMany
    call CalcSqRoot
  else if (QMain.EQ.7) then
    QSuffix = 2
    call SelectMany
    call OpWithKay
  else if (QMain.EQ.8) then
    QSuffix = 2
    call SelectMany
    call CalcMean
    call SliceMean
  else if (QMain.EQ.9) then
    QSuffix = 2
    call SelectMany
    call SliceMany
  else if (QMain.EQ.10) then
    QSuffix = 2
    call SelectMany
    call GradientMany
  else if (QMain.EQ.11) then
    QSuffix = 3
    call SelectMany
    call GaussMany
  else
    print*, "  >   1. First line minus second line"
    print*, "  >   2. Mean of a number of lines"
    print*, "  >   3. St-dev of a number of lines"
    print*, "  >   4. LSR best-fit of lines v time"
    print*, "  >   5. LSR best-fit of lines v lines"
    print*, "  >   6. Square root a number of lines"
    print*, "  >   7. Operate on a number of lines with k"
    print*, "  >   8. Average lines and slice into 30y means" 
    print*, "  >   9. Slice lines into 30y means"
    print*, "  >  10. Obtain gradient of individual lines"
    print*, "  >  11. Gauss-smooth"
    print*, "  >  99. Exit"
  end if
  
  if (QSuffix.EQ.1) call DumpToLin
  if (QSuffix.EQ.2) call DumpToTim
  if (QSuffix.EQ.3) call DumpMultiToLin
  
  print*
  print*, "  > Main menu:"
  do
	read (*,*,iostat=ReadStatus), QMain
	if (ReadStatus.LE.0.AND.QMain.GE.1) exit
  end do
  
  if (QMain.EQ.99) exit
end do  

deallocate (LinLoaded,WorkLinNames)

print*

close (99)

contains

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

subroutine Initialise

print*, "  > ***** CompareTim *****"
print*, "  > Operates on .tim and .lin files; saves .tim file"
print*

MissAccept = 10.0

print*, "  > Select parameters that will govern operations."
call GridSelect   (WorkGrid,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call PeriodSelect (WorkYearN,WorkDecN,WorkADYear)

print*, "  > Select the number of .tim files to load."
do
	read (*,*,iostat=ReadStatus), WorkFileTimN
	if (ReadStatus.LE.0.AND.WorkFileTimN.GE.0) exit
end do
print*, "  > Select the number of .lin files to load."
do
	read (*,*,iostat=ReadStatus), WorkFileLinN
	if (ReadStatus.LE.0.AND.WorkFileLinN.GE.0) exit
end do
print*, "  > Enter the total number of lines to load."
do
	read (*,*,iostat=ReadStatus), WorkLinN
	if (ReadStatus.LE.0.AND.WorkLinN.GE.1) exit
end do

allocate (LinLoaded(WorkLinN,WorkYearN), &
	  WorkLinNames (WorkLinN),	 stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

LinLoaded = MissVal
WorkLinNames = "blank"

FirstFree = 0

if (WorkFileTimN.GE.1) call GrabTim
if (WorkFileLinN.GE.1) call GrabLin

end subroutine Initialise

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

subroutine GrabTim

 print*, "  > We load the .tim files."
 do XFile = 1, WorkFileTimN
  print*, "  > File number : ", XFile
  
  call LoadTim (WorkYearN,WorkADYear,FileLinN,FileLinNames,FileLin)
  
  do XLin = 1, FileLinN
    FirstFree = FirstFree + 1
    
    if (FirstFree.GT.WorkLinN) then
      print*, "  > Not enough lines allocated to load line: ", FileLinNames (XLin)
    else
      WorkLinNames (FirstFree) = FileLinNames (XLin)
      do XYear = 1, WorkYearN
        LinLoaded (FirstFree,XYear) = FileLin (XLin,XYear)
      end do
    end if 
  end do
  
  deallocate (FileLinNames,FileLin,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
 end do

end subroutine GrabTim

!*******************************************************************************
! load .lin

subroutine GrabLin

 print*, "  > We load the .lin files."
 do XFile = 1, WorkFileLinN
  print*, "  > File number : ", XFile
  
  call LoadLin (WorkYearN,WorkADYear,FileLinN,FileLinNames,FileLin)
  
  do XLin = 1, FileLinN
    FirstFree = FirstFree + 1
    
    if (FirstFree.GT.WorkLinN) then
      print*, "  > Not enough lines allocated to load line: ", FileLinNames (XLin)
    else
      WorkLinNames (FirstFree) = FileLinNames (XLin)
      do XYear = 1, WorkYearN
        LinLoaded (FirstFree,XYear) = FileLin (XLin,XYear)
      end do
    end if 
  end do
  
  deallocate (FileLinNames,FileLin,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
 end do

end subroutine GrabLin

!*******************************************************************************
! select two lines

subroutine SelectSuffix

print*, "  > Save as a .lin (=1) or .tim (=2) ? "
do
  read (*,*,iostat=ReadStatus), QSuffix
  if (ReadStatus.LE.0.AND.QSuffix.GE.1.AND.QSuffix.LE.2) exit
end do

end subroutine SelectSuffix

!*******************************************************************************
! select two lines

subroutine SelectTwo

print*, "  > Operation A minus B: select line A (0=list): "
do
  read (*,*,iostat=ReadStatus), SelectLinA
  if (ReadStatus.LE.0.AND.SelectLinA.EQ.0) then
    do XLin = 1, WorkLinN
      print "(I6,A3,A40)", XLin, " : ", WorkLinNames (XLin)
    end do
  end if
  if (ReadStatus.LE.0.AND.SelectLinA.GE.1.AND.SelectLinA.LE.WorkLinN) exit
end do

print*, "  > Operation A minus B: select line B (0=list): "
do
  read (*,*,iostat=ReadStatus), SelectLinB
  if (ReadStatus.LE.0.AND.SelectLinB.EQ.0) then
    do XLin = 1, WorkLinN
      print "(I6,A3,A40)", XLin, " : ", WorkLinNames (XLin)
    end do
  end if
  if (ReadStatus.LE.0.AND.SelectLinB.GE.1.AND.SelectLinB.LE.WorkLinN) exit
end do

end subroutine SelectTwo

!*******************************************************************************
! select any number of lines

subroutine SelectMany

if (allocated(LineSelect)) deallocate (LineSelect)
allocate (LineSelect(WorkLinN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
LineSelect = 0

print*, "  > Select lines (0=list ; -99=end): "
do
  read (*,*,iostat=ReadStatus), SelectLin
  if (ReadStatus.LE.0.AND.SelectLin.EQ.0) then
    do XLin = 1, WorkLinN
      print "(I6,A3,A40)", XLin, " : ", WorkLinNames (XLin)
    end do
  end if
  if (ReadStatus.LE.0.AND.SelectLin.GE.1.AND.SelectLin.LE.WorkLinN) LineSelect(SelectLin) = 1
  if (ReadStatus.LE.0.AND.SelectLin.EQ.-99) exit
end do

SelectLinN = 0
do XLin = 1, WorkLinN
  if (LineSelect(XLin).EQ.1) SelectLinN = SelectLinN + 1
end do

end subroutine SelectMany

!*******************************************************************************
! subtract Two from One

subroutine TwoFromOne

allocate (LinSaved(1,WorkYearN), &
	  SaveLinNames(1)      , stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: TwoFromOne: Allocation failure #####"
LinSaved = MissVal
SaveLinNames (1) = "one line subtracted from the other"

  print*, "  > Calculate abs diff (=1), % diff (=2), fraction (=3) ? "
  do
	read (*,*,iostat=ReadStatus), QAnom
	if (ReadStatus.LE.0.AND.QAnom.GE.1.AND.QAnom.LE.3) exit
  end do

  if      (QAnom.EQ.1) then
   do XYear = 1, WorkYearN
    if (LinLoaded(SelectLinA,XYear).NE.MissVal.AND.LinLoaded(SelectLinB,XYear).NE.MissVal) then
      LinSaved (1,XYear) = LinLoaded(SelectLinA,XYear) - LinLoaded(SelectLinB,XYear)
    end if
   end do
  else if (QAnom.EQ.2) then
   do XYear = 1, WorkYearN
    if (LinLoaded(SelectLinA,XYear).NE.MissVal.AND.LinLoaded(SelectLinB,XYear).NE.MissVal) then
     if (LinLoaded(SelectLinB,XYear).NE.0) then
      LinSaved (1,XYear) = 100.0 * (LinLoaded(SelectLinA,XYear) - LinLoaded(SelectLinB,XYear)) &
      			/ abs (LinLoaded(SelectLinB,XYear))
     end if
    end if
   end do
  else if (QAnom.EQ.3) then
   do XYear = 1, WorkYearN
    if (LinLoaded(SelectLinA,XYear).NE.MissVal.AND.LinLoaded(SelectLinB,XYear).NE.MissVal) then
     if (LinLoaded(SelectLinB,XYear).NE.0) then
      LinSaved (1,XYear) = LinLoaded(SelectLinA,XYear) / LinLoaded(SelectLinB,XYear)
     end if
    end if
   end do
  end if

end subroutine TwoFromOne

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

subroutine CalcMean

if (QSuffix.EQ.1) then 
  allocate (LinSaved(1,WorkYearN),  &
	    SaveLinNames(1)      ,  stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcMean: Allocation failure #####"
  LinSaved = MissVal
  SaveLinNames (1) = "mean of a number of lines"
else if (QSuffix.EQ.2) then
  allocate (TimSaved(1,WorkYearN),  &
	    SaveLinNames(1)      ,  stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcMean: Allocation failure #####"
  TimSaved = MissVal
  SaveLinNames (1) = "mean of a number of lines"
  SelectLinN = 1
end if

do XYear = 1, WorkYearN
   RunTot = 0.0
   RunEn  = 0.0
   
   do XLin = 1, WorkLinN
     if (LineSelect(XLin).EQ.1.AND.LinLoaded(XLin,XYear).NE.MissVal) then
      RunTot = RunTot + LinLoaded(XLin,XYear)
      RunEn  = RunEn  + 1.0
     end if
   end do
   
   if (RunEn.GE.1) then
     if (QSuffix.EQ.1) LinSaved (1,XYear) = RunTot / RunEn
     if (QSuffix.EQ.2) TimSaved (1,XYear) = RunTot / RunEn
   end if
end do

write (99,*), "have calc mean"
 
end subroutine CalcMean

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

subroutine CalcStDev

allocate (LinSaved(1,WorkYearN), &
	  SaveLinNames(1)      , stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcStDev: Allocation failure #####"
LinSaved = MissVal
SaveLinNames (1) = "stdev of a number of lines"

 do XYear = 1, WorkYearN
   RunSqTot = 0.0
   RunTot   = 0.0
   RunEn    = 0.0
   
   do XLin = 1, WorkLinN
     if (LineSelect(XLin).EQ.1.AND.LinLoaded(XLin,XYear).NE.MissVal) then
      RunSqTot = RunSqTot + (LinLoaded(XLin,XYear) * LinLoaded(XLin,XYear))
      RunTot   = RunTot   + LinLoaded(XLin,XYear)
      RunEn    = RunEn    + 1.0
     end if
   end do
   
   if (RunEn.GE.1) then
     LinSaved (1,XYear) = (RunSqTot/RunEn) - ((RunTot/RunEn)*(RunTot/RunEn))
     LinSaved (1,XYear) = sqrt (LinSaved (1,XYear))
   end if
 end do
 
end subroutine CalcStDev

!*******************************************************************************
! calc SqRoot

subroutine CalcSqRoot

allocate (TimSaved(SelectLinN,WorkYearN), &
	  SaveLinNames(SelectLinN)      , stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcSqRoot: Allocation failure #####"
TimSaved = MissVal

SelectLinA = 0

do XLin = 1, WorkLinN
 if (LineSelect(XLin).EQ.1) then
  SelectLinA = SelectLinA + 1
  
  SaveLinNames (SelectLinA) = WorkLinNames (XLin)
  
  do XYear = 1, WorkYearN
   if (LinLoaded(XLin,XYear).NE.MissVal.AND.LinLoaded(XLin,XYear).GE.0) &
     TimSaved (SelectLinA,XYear) = sqrt ( LinLoaded(XLin,XYear) )
  end do
 end if
end do
   
end subroutine CalcSqRoot

!*******************************************************************************
! operate with constant

subroutine OpWithKay

allocate (TimSaved(SelectLinN,WorkYearN), &
	  SaveLinNames(SelectLinN)      , stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcSqRoot: Allocation failure #####"
TimSaved = MissVal

SelectLinA = 0

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

do XLin = 1, WorkLinN
 if (LineSelect(XLin).EQ.1) then
  SelectLinA = SelectLinA + 1
  
  SaveLinNames (SelectLinA) = WorkLinNames (XLin)
  
  do XYear = 1, WorkYearN
   if (LinLoaded(XLin,XYear).NE.MissVal) then
     if (Operator.EQ.1) TimSaved (SelectLinA,XYear) = LinLoaded(XLin,XYear) / RealConstant
     if (Operator.EQ.2) TimSaved (SelectLinA,XYear) = LinLoaded(XLin,XYear) * RealConstant
     if (Operator.EQ.3) TimSaved (SelectLinA,XYear) = LinLoaded(XLin,XYear) - RealConstant
     if (Operator.EQ.4) TimSaved (SelectLinA,XYear) = LinLoaded(XLin,XYear) + RealConstant
     if (Operator.EQ.5) TimSaved (SelectLinA,XYear) = sqrt (LinLoaded(XLin,XYear))
     if (Operator.EQ.6) TimSaved (SelectLinA,XYear) = LinLoaded(XLin,XYear) / IntConstant
   end if
  end do
 end if
end do
      
end subroutine OpWithKay

!*******************************************************************************
! sliceify into 30y means

subroutine SliceMean

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

do XYear = 1, WorkYearN
  VecSaved (XYear) = LinSaved (1,XYear)
end do

call SimpleSlice (WorkYearN,30,MissAccept,VecSaved)

do XYear = 1, WorkYearN
  LinSaved (1,XYear) = VecSaved (XYear)
end do

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

end subroutine SliceMean

!*******************************************************************************
! sliceify into 30y means

subroutine SliceMany

allocate (TimSaved(SelectLinN,WorkYearN),  &
	  SaveLinNames(SelectLinN)      ,  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SliceMany: Allocation failure #####"
TimSaved = MissVal

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

XSelect = 0							! selected line 1...SelectLinN

do XLin = 1, WorkLinN						! loop through every line
  if (LineSelect(XLin).EQ.1) then				! if it is a chosen line...
    
    XSelect = XSelect + 1					! ...then increment chosen line
    
    SaveLinNames (XSelect) = WorkLinNames (XLin)		! ...save name of line
    
    do XYear = 1, WorkYearN					! ...load up vector with data
      VecSaved (XYear) = LinLoaded (XLin,XYear)
    end do
    
    call SimpleSlice (WorkYearN,30,MissAccept,VecSaved)		! ...sliceify it

    do XYear = 1, WorkYearN					! ...put it in saved lines
      TimSaved (XSelect,XYear) = VecSaved (XYear)
    end do
    
  end if
end do

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

end subroutine SliceMany

!*******************************************************************************
! calculate gradients

subroutine GradientMany

BandSide  = 15				! in calc the gradient we use BandSide years either side of Year
BandWidth = (BandSide * 2) + 1

allocate (TimSaved(SelectLinN,WorkYearN),  &
	  SaveLinNames(SelectLinN)      ,  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GradientMany: Allocation failure #####"
TimSaved = MissVal

allocate (OpYear(1,BandWidth), &
	  OpData(1,BandWidth), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GradientMany: Allocation failure #####"
OpYear = MissVal
OpData = MissVal

XSelect = 0							! selected line 1...SelectLinN

do XLin = 1, WorkLinN						! loop through every line
  if (LineSelect(XLin).EQ.1) then				! if it is a chosen line...
    
    XSelect = XSelect + 1					! ...then increment chosen line
    
    SaveLinNames (XSelect) = WorkLinNames (XLin)		! ...save name of line
    
    do XYear0 = 1, (WorkYearN-BandWidth+1)			! ...iterate by start year of band
      do XYear = 1, BandWidth					! ......iterate by band year
        OpYear (1,XYear) = XYear				! .........fill op vector with years
        OpData (1,XYear) = LinLoaded (XLin,(XYear+XYear0-1))	! .........fill op vector with data
      end do

      call LinearLSR (1,BandWidth,OpYear,OpData,OpAye,OpBee,OpPea)	! ......calc gradient

      TimSaved (XSelect,(XYear0+BandSide)) = OpAye			! ......put it in saved lines
    end do
        
  end if
end do

deallocate (OpYear, &
	    OpData, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GradientMany: Deallocation failure #####"

end subroutine GradientMany

!*******************************************************************************
! calc best fit by LSR for lines v times

subroutine CalcBestFitTime

print*, "  > Set intercept to zero (1=no,2=yes) ?"
do
	read (*,*,iostat=ReadStatus), QIntercept
	if (ReadStatus.LE.0.AND.QIntercept.GE.1.AND.QIntercept.LE.2) exit
end do

AllSumX  = 0.0
AllSumY  = 0.0
AllSumXX = 0.0
AllSumYY = 0.0
AllSumXY = 0.0
AllEn    = 0.0
   
do XYear = 1, WorkYearN
   do XLin = 1, WorkLinN
     if (LineSelect(XLin).EQ.1.AND.LinLoaded(XLin,XYear).NE.MissVal) then
      AllSumX  = AllSumX  + WorkADYear(XYear)
      AllSumY  = AllSumY  + LinLoaded(XLin,XYear)
      AllSumXX = AllSumXX + (WorkADYear(XYear)*WorkADYear(XYear))
      AllSumYY = AllSumYY + (LinLoaded(XLin,XYear)*LinLoaded(XLin,XYear))
      AllSumXY = AllSumXY + (WorkADYear(XYear)*LinLoaded(XLin,XYear))
      AllEn    = AllEn    + 1.0
     end if
   end do
end do

Numer  = 0.0
Denom  = 0.0
AllAye = MissVal
AllBee = MissVal
AllAre = MissVal
 
if (QIntercept.EQ.1) then
  Numer = (AllEn * AllSumXY) - (AllSumX * AllSumY)	! y=ax+b 'a'
  Denom = (AllEn * AllSumXX) - (AllSumX * AllSumX)
  if (Denom.NE.0) AllAye = Numer / Denom
  
  Numer = AllSumY - (AllAye * AllSumX)			! y=ax+b 'b'
  if (AllEn.NE.0) AllBee = Numer / AllEn
  							
  if (AllEn.NE.0) then					! y=ax+b 'r'
  	Numer  = (AllSumXY - ((AllSumX / AllEn) * (AllSumY / AllEn))) / AllEn
  	Denom  = sqrt ((AllSumXX - ((AllSumX / AllEn) ** 2)) / AllEn) * &
  		 sqrt ((AllSumYY - ((AllSumY / AllEn) ** 2)) / AllEn)
  	if (Denom.NE.0) AllAre = Numer / Denom
  end if  
  
  print*, "  > y=ax+b : values of a, b, r"
  print "(3F12.6)", AllAye, AllBee, AllAre
else
  if (AllSumXX.NE.0) AllAye = AllSumXY / AllSumXX	! y=ax   'a'
  
  if (AllEn.NE.0) then					! y=ax+b 'r'
  	Numer  = (AllSumXY - ((AllSumX / AllEn) * (AllSumY / AllEn))) / AllEn
  	Denom  = sqrt ((AllSumXX - ((AllSumX / AllEn) ** 2)) / AllEn) * &
  		 sqrt ((AllSumYY - ((AllSumY / AllEn) ** 2)) / AllEn)
  	if (Denom.NE.0) AllAre = Numer / Denom
  end if  
  print*, "  > y=ax : value of a ; y=ax+b : value of r"
  print "(2F12.6)", AllAye, AllAre
end if

end subroutine CalcBestFitTime

!*******************************************************************************
! calc best fit by LSR for lines v lines

subroutine CalcBestFitLines

if (allocated(PairSelect)) deallocate (PairSelect)
allocate (PairSelect(WorkLinN,2), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
PairSelect = 0
FirstFree  = 0

print*, "  > Enter the set of x-var lines (0=list ; -99=end): "
do
  read (*,*,iostat=ReadStatus), SelectLin
  if (ReadStatus.LE.0.AND.SelectLin.EQ.0) then
    do XLin = 1, WorkLinN
      print "(I6,A3,A40)", XLin, " : ", WorkLinNames (XLin)
    end do
  end if
  if (ReadStatus.LE.0.AND.SelectLin.GE.1.AND.SelectLin.LE.WorkLinN) then
    FirstFree = FirstFree + 1
    PairSelect(SelectLin,1) = FirstFree
  end if
  if (ReadStatus.LE.0.AND.SelectLin.EQ.-99) exit
end do
SelectLinN = FirstFree

print*, "  > Enter the set of y-var lines. "
do XLin = 1, SelectLinN
  print*, "  > Enter the y-var line corresponding to x-var line: ", XLin
  do
    read (*,*,iostat=ReadStatus), SelectLin
    if (ReadStatus.LE.0.AND.SelectLin.GE.1.AND.SelectLin.LE.WorkLinN) exit
  end do
  
  PairSelect(SelectLin,2) = XLin
end do

print*, "  > Set intercept to zero (1=no,2=yes) ?"
do
	read (*,*,iostat=ReadStatus), QIntercept
	if (ReadStatus.LE.0.AND.QIntercept.GE.1.AND.QIntercept.LE.2) exit
end do

AllSumX  = 0.0
AllSumY  = 0.0
AllSumXX = 0.0
AllSumYY = 0.0
AllSumXY = 0.0
AllEn    = 0.0
   
do XLin = 1, SelectLinN
  SelectLinA = MissVal
  SelectLinB = MissVal
  
  do XFind = 1, WorkLinN
    if (PairSelect(XFind,1).EQ.XLin) SelectLinA = XFind
    if (PairSelect(XFind,2).EQ.XLin) SelectLinB = XFind
  end do
  
  if (SelectLinA.NE.MissVal.AND.SelectLinB.NE.MissVal) then
   do XYear = 1, WorkYearN
     if (LinLoaded(SelectLinA,XYear).NE.MissVal.AND.LinLoaded(SelectLinB,XYear).NE.MissVal) then
      AllSumX  = AllSumX  + LinLoaded(SelectLinA,XYear)
      AllSumY  = AllSumY  + LinLoaded(SelectLinB,XYear)
      AllSumXX = AllSumXX + (LinLoaded(SelectLinA,XYear)*LinLoaded(SelectLinA,XYear))
      AllSumYY = AllSumYY + (LinLoaded(SelectLinB,XYear)*LinLoaded(SelectLinB,XYear))
      AllSumXY = AllSumXY + (LinLoaded(SelectLinA,XYear)*LinLoaded(SelectLinB,XYear))
      AllEn    = AllEn    + 1.0
     end if
   end do
  end if
end do

Numer  = 0.0
Denom  = 0.0
AllAye = MissVal
AllBee = MissVal
AllAre = MissVal
 
if (QIntercept.EQ.1) then
  Numer = (AllEn * AllSumXY) - (AllSumX * AllSumY)	! y=ax+b 'a'
  Denom = (AllEn * AllSumXX) - (AllSumX * AllSumX)
  if (Denom.NE.0) AllAye = Numer / Denom
  
  Numer = AllSumY - (AllAye * AllSumX)			! y=ax+b 'b'
  if (AllEn.NE.0) AllBee = Numer / AllEn
  							
  if (AllEn.NE.0) then					! y=ax+b 'r'
  	Numer  = (AllSumXY - ((AllSumX / AllEn) * (AllSumY / AllEn))) / AllEn
  	Denom  = sqrt ((AllSumXX - ((AllSumX / AllEn) ** 2)) / AllEn) * &
  		 sqrt ((AllSumYY - ((AllSumY / AllEn) ** 2)) / AllEn)
  	if (Denom.NE.0) AllAre = Numer / Denom
  end if  
  
  print*, "  > y=ax+b : values of a, b, r"
  print "(3F12.6)", AllAye, AllBee, AllAre
else
  if (AllSumXX.NE.0) AllAye = AllSumXY / AllSumXX	! y=ax   'a'
  
  if (AllEn.NE.0) then					! y=ax+b 'r'
  	Numer  = (AllSumXY - ((AllSumX / AllEn) * (AllSumY / AllEn))) / AllEn
  	Denom  = sqrt ((AllSumXX - ((AllSumX / AllEn) ** 2)) / AllEn) * &
  		 sqrt ((AllSumYY - ((AllSumY / AllEn) ** 2)) / AllEn)
  	if (Denom.NE.0) AllAre = Numer / Denom
  end if  
  print*, "  > y=ax : value of a ; y=ax+b : value of r"
  print "(2F12.6)", AllAye, AllAre
end if

deallocate (PairSelect)

end subroutine CalcBestFitLines

!*******************************************************************************
! Gauss smooth

subroutine GaussMany

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

allocate (LinSaved(SelectLinN,WorkYearN), &
	  SaveLinNames(SelectLinN)      , &
	  LinVec(WorkYearN)             , &
	  LoPass(WorkYearN)             , &
	  HiPass(WorkYearN)             , stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcStDev: Allocation failure #####"
LinSaved     = MissVal
SaveLinNames = ""

CurrentLin = 0
do XLin = 1, WorkLinN
  if (LineSelect(XLin).EQ.1) then
    CurrentLin = CurrentLin + 1
    LoPass = MissVal
    HiPass = MissVal
    
    do XYear = 1, WorkYearN
      LinVec (XYear) = LinLoaded(XLin,XYear)
    end do
    
    call GaussSmooth (WorkYearN,SliceLen,1,LinVec,LoPass,HiPass)
    
    do XYear = 1, WorkYearN
      LinSaved (CurrentLin,XYear) = LoPass(XYear)
    end do
  end if
end do

end subroutine GaussMany

!*******************************************************************************
! save .lin file (QSuffix=1)

subroutine DumpToLin

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

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

end subroutine DumpToLin

!*******************************************************************************
! save .tim file (QSuffix=2)

subroutine DumpToTim

call SaveTim (SelectLinN,WorkYearN,Blank,Blank,SaveLinNames,WorkADYear,TimSaved)

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

end subroutine DumpToTim

!*******************************************************************************
! save .lin file (QSuffix=3)

subroutine DumpMultiToLin

call SaveLin (SelectLinN,WorkYearN,SaveLinNames,WorkADYear,LinSaved)

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

end subroutine DumpMultiToLin

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

end program CompareTim