! opann.f90
! f90 program written by Tim Mitchell on 29.01.01
! last modified on 21.11.01
! tool to operate on standard .ann files
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../obs/opann filenames.f90 annfiles.f90 sortmod.f90 usesort.f90 
!	basicfun.f90 time.f90 plainperm.f90 regress.f90 pattscale.f90 
!	./../obs/opann.f90 2> /tyn1/tim/scratch/stderr.txt

program OpAnn

use FileNames
use AnnFiles
use SortMod
use UseSort
use BasicFun
use Time
use PlainPerm
use Regress
use PattScale

implicit none

real, pointer, dimension (:,:)	:: Data, FileData
real, pointer, dimension (:)	:: OpVec, ResultVecA, ResultVecB, LSRExe, LSRWye

integer, pointer, dimension (:)		  :: YearAD, FileYearAD

character (len=80), pointer, dimension (:) :: FileBatch
character (len= 9), pointer, dimension (:) :: ColTitles, FileColTitles

real, parameter :: MissVal = -999.0

character (len=80), parameter :: Blank = ""

real :: RealConstant, Speed, Force, MissAccept
real :: DataMiss, DataMean, DataTot, DataFrac, OpMiss,OpEn,OpFrac,OpTot,OpBase,OpThresh,OpMean
real :: Forc2co2,FitAlpha,FitBeta,Sens2co2Init,DTDt,Forc,EffClimSens,LastEffClimSens
real :: Trash1,Trash2, KayX,KayY

integer :: ReadStatus, AllocStat, MenuChoice, NextFree 
integer :: ChosenCol, ChosenColA, ChosenColB, ChosenColF, ChosenColZ, ChosenOp
integer :: StartYear, EndYear, FirstYear, YearN, FileN
integer :: BlockYearN, BlockStartAD, BlockEndAD, BlockStartX, BlockEndX, XFile
integer :: XYear, XFileYear, XCol, XFileCol, XLoadCol, XSaveCol, XTrendYear
integer :: StringLen, ColN, LoadColN, FileColN, SaveColN
integer :: FileYear0,FileYear1,Year0,Year1,Day0,Day1,Month0,Month1,Julian0,Julian1,TrendYear0,TrendYear1
integer :: QComp, QIntFloat, QFileComp, QKeepMiss, QAbsRel
integer :: NonBlank, IntConstant, Overlap, THalf, PerLen, GapLen,TrendTLen
integer :: BaseYear0,BaseYear1, BaseADYear0,BaseADYear1, BaseYearN

character (len=80) :: GivenFile, OrigFile, LineFormat, Waste
character (len=1)  :: IntFloat

!*******************************************************************************
! main

call Intro

do
  print*, "  > Main menu. Make your choice. (0=list)"
  do
	read (*,*,iostat=ReadStatus), MenuChoice
	if (ReadStatus.LE.0) exit
  end do
	
  if      (MenuChoice.EQ.1)  then
    print*, "  > Are you sure? (1=no,2=yes)"
    do
	read (*,*,iostat=ReadStatus), MenuChoice
	if (ReadStatus.LE.0.AND.MenuChoice.GE.1.AND.MenuChoice.LE.2) exit
    end do
    if (MenuChoice.EQ.1) print*, "  > No changes made."
    if (MenuChoice.EQ.2) call Intro
  else if (MenuChoice.EQ. 2)  then
    call OperateAk
  else if (MenuChoice.EQ. 3)  then
    call OperateAB
  else if (MenuChoice.EQ. 4)  then
    call GaussSmoothA
  else if (MenuChoice.EQ. 5)  then
    call AnomaliseA
  else if (MenuChoice.EQ. 6)  then
    call CombineAB
  else if (MenuChoice.EQ. 7)  then
    call MaskAB
  else if (MenuChoice.EQ.10)  then
    call LoadFromANN
  else if (MenuChoice.EQ.11)  then
    call LoadFromANNset
  else if (MenuChoice.EQ.20)  then
    call DumpToANN
  else if (MenuChoice.EQ.30)  then
    call CheckMaster
  else if (MenuChoice.EQ.40)  then
    call CalcEquT
  else if (MenuChoice.EQ.41)  then
    call PrintPeriodMean
  else if (MenuChoice.NE.99) then
    print*, "  >  1. Reinitialise"
    print*, "  >  2. Operate   A.k"
    print*, "  >  3. Operate   A.B"
    print*, "  >  4. Smooth    A with Gaussian filter"
    print*, "  >  5. Anomalise A relative to period"
    print*, "  >  6. Combine columns A and B"
    print*, "  >  7. Where A=x B=y"
    print*, "  > 10. Load data from .ann"
    print*, "  > 11. Load data from set of .ann files"
    print*, "  > 20. Save data to   .ann"
    print*, "  > 30. Check main arrays"
    print*, "  > 40. Calc equilib T: A(gloT),B(forc)"
    print*, "  > 41. Print period mean"
    print*, "  > 99. Exit"
  end if
  
  if (MenuChoice.EQ.99) exit
end do

call Conclude

contains

!*******************************************************************************
! intro

subroutine Intro

open (99,file="./../../../scratch/log-opann.dat",status="replace",action="write")

print*
print*, "  > ##### OpAnn.f90 ##### Tool for operating on .ann files #####"
print*

MissAccept = 10.0

print*, "  > Enter the number of columns in the master array:"
do
	read (*,*,iostat=ReadStatus), ColN
	if (ReadStatus.LE.0 .AND. ColN.GE.1) exit
end do

print*, "  > Enter the start and end years AD of the master array:"
do
	read (*,*,iostat=ReadStatus), StartYear, EndYear
	if (ReadStatus.LE.0 .AND. StartYear.LE.EndYear) exit
end do

YearN  = EndYear - StartYear + 1
NextFree = 1

allocate ( Data      (YearN, ColN), &
	   ColTitles (ColN), &
	   YearAD    (YearN),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure #####"
Data      = MissVal
ColTitles = 'X'

do XYear = 1, YearN
  YearAD (XYear) = StartYear + XYear - 1
end do

print*

end subroutine Intro

!*******************************************************************************
! load from .ann file

subroutine LoadFromANN

call LoadANN      (Blank, FileYearAD, FileColTitles, FileData)
call CommonVecPer (FileYearAD,YearAD,FileYear0,FileYear1,Year0,Year1)

if (FileYear0.EQ.MissVal) then
  print*, "  > The loaded file has no period in common with the master array."
else
  print*, "  > The common period is: ", FileYearAD(FileYear0), FileYearAD(FileYear1)
end if

FileColN = size (FileColTitles)

print*, "  > COL    TITLE"
do XCol = 1, FileColN
  print "(i8,a9)", XCol, FileColTitles(XCol)
end do

if (FileColN.GT.1) then
  print*, "  > Enter the number of columns to include in array: "
  do
	read (*,*,iostat=ReadStatus), LoadColN
	if (ReadStatus.LE.0 .AND. LoadColN.GE.1 .AND. LoadColN.LE.FileColN) exit
  end do
else
  LoadColN = 1
end if

do XLoadCol = 1, LoadColN
  if (LoadColN.LT.FileColN) then
    print "(a47,i4)", "   > Select column from file --> array column: ", NextFree
    do
	read (*,*,iostat=ReadStatus), ChosenCol
	if (ReadStatus.LE.0 .AND. ChosenCol.GE.1 .AND. ChosenCol.LE.FileColN) exit
    end do
  else
    ChosenCol = XLoadCol
  end if

  ColTitles(NextFree) = FileColTitles(ChosenCol)
  
  do XFileYear = FileYear0, FileYear1
   XYear = Year0 + XFileYear - FileYear0
  
   Data(XYear,NextFree) = FileData(XFileYear,ChosenCol)
  end do

  call StoreMoveOn  
end do

deallocate (FileYearAD,FileColTitles,FileData, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadFromANN: Deallocation failure #####"

print*

end subroutine LoadFromANN

!*******************************************************************************
! load from set of .ann files

subroutine LoadFromANNset

GivenFile = ""
call GetBatch (GivenFile,FileBatch)
FileN = size(FileBatch)

do XFile = 1, FileN
  call LoadANN      (FileBatch(XFile), FileYearAD, FileColTitles, FileData)
  call CommonVecPer (FileYearAD,YearAD,FileYear0,FileYear1,Year0,Year1)
  LoadColN = size (FileColTitles)
  
  do XLoadCol = 1, LoadColN
    do XFileYear = FileYear0, FileYear1
      XYear = Year0 + XFileYear - FileYear0
      Data(XYear,NextFree) = FileData(XFileYear,XLoadCol)
    end do
    ColTitles(NextFree) = FileColTitles(XLoadCol)
    call StoreMoveOn  
  end do
  
  deallocate (FileYearAD,FileColTitles,FileData, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadFromANN: Deallocation failure #####"
end do

end subroutine LoadFromANNset

!*******************************************************************************
! smooth a column with a Gaussian filter

subroutine GaussSmoothA

call SelectColumn

print*, "  > Select the periodicity at which to smooth: "
do
	read (*,*,iostat=ReadStatus), THalf
	if (ReadStatus.LE.0 .AND. THalf.GE.1 .AND. THalf.LE.YearN) exit
end do

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

allocate (OpVec      (YearN), &
	  ResultVecA (YearN), &
	  ResultVecB (YearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GaussSmoothA: Allocation failure #####"
OpVec = MissVal

do XYear = 1, YearN
    OpVec (XYear) = Data (XYear,ChosenCol)
end do
  
call GaussSmooth (YearN,THalf,QKeepMiss,OpVec,ResultVecA,ResultVecB)

do XYear = 1, YearN
    Data (XYear,NextFree) = ResultVecA (XYear)
end do

deallocate (OpVec,ResultVecA,ResultVecB, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GaussSmoothA: Deallocation failure #####"

call StoreMoveOn
print*

end subroutine GaussSmoothA

!*******************************************************************************
! anomalise a period against a column within it

subroutine AnomaliseA

call SelectColumn

print*, "  > Select first, last years AD as base period: "
do
	read (*,*,iostat=ReadStatus), BaseADYear0,BaseADYear1
	if (ReadStatus.LE.0.AND.BaseADYear0.GE.YearAD(1).AND.BaseADYear1.LE.YearAD(YearN) &
				.AND.BaseADYear1.GE.BaseADYear0) exit
end do

print*, "  > Use absolute (=0, minus) or relative (=1, divide) anomalies? "
do
	read (*,*,iostat=ReadStatus), QAbsRel
	if (ReadStatus.LE.0.AND.QAbsRel.GE.0.AND.QAbsRel.LE.1) exit
end do

BaseYearN = BaseADYear1 - BaseADYear0 + 1
BaseYear0 = BaseADYear0 - YearAD(1)   + 1
BaseYear1 = BaseYear0   + BaseYearN   - 1

OpThresh  = real(BaseYearN) * (MissAccept / 100.0)

OpTot  = 0.0
OpMiss = 0.0
do XYear = BaseYear0, BaseYear1
    if (Data (XYear,ChosenCol) .NE. MissVal) then
    	OpTot = OpTot + Data (XYear,ChosenCol)
    else
        OpMiss = OpMiss + 1
    end if
end do

if (OpMiss.GT.OpThresh) then
  print*, "  > Too many missing values in base period to calc anomalies."
else
  OpBase = OpTot / (real(BaseYearN) - OpMiss)
end if
  
do XYear = 1, YearN
  if (Data (XYear,ChosenCol).NE.MissVal) then
    if (QAbsRel.EQ.0) then
      Data (XYear,NextFree) = Data (XYear,ChosenCol) - OpBase
    else
      if (OpBase.NE.0) Data (XYear,NextFree) = Data (XYear,ChosenCol) / OpBase
    end if
  else
    Data (XYear,NextFree) = MissVal
  end if
end do

call StoreMoveOn
print*

end subroutine AnomaliseA

!*******************************************************************************
! mask A . B

subroutine MaskAB

print*, "  > Select columns A and B in turn: "
call SelectColumn
ChosenColA = ChosenCol
call SelectColumn
ChosenColB = ChosenCol

print*, "  > Select x and y:"
do
	read (*,*,iostat=ReadStatus), KayX,KayY
	if (ReadStatus.LE.0) exit
end do

do XYear = 1, YearN
    Data(XYear,NextFree) = Data(XYear,ChosenColB)
    if (Data(XYear,ChosenColA).EQ.KayX) Data(XYear,NextFree) = KayY
end do

call StoreMoveOn
print*

end subroutine MaskAB

!*******************************************************************************
! combine A . B

subroutine CombineAB

print*, "  > Select columns A and B in turn: "
call SelectColumn
ChosenColA = ChosenCol
call SelectColumn
ChosenColB = ChosenCol

do XYear = 1, YearN
    if      (Data(XYear,ChosenColA).NE.MissVal.AND.Data(XYear,ChosenColB).EQ.MissVal) then
      Data(XYear,NextFree) = Data(XYear,ChosenColA)
    else if (Data(XYear,ChosenColA).EQ.MissVal.AND.Data(XYear,ChosenColB).NE.MissVal) then
      Data(XYear,NextFree) = Data(XYear,ChosenColB)
    else
      Data(XYear,NextFree) = MissVal
    end if
end do

call StoreMoveOn
print*

end subroutine CombineAB

!*******************************************************************************
! operate A . constant

subroutine OperateAk

print*, "  > Select the operation: A/k=1, A*k=2, A-k=3, A+k=4, sqrt(A)=5, A**k=6,"
print*, "  > 			   abs(A)=7, log(A)=8, e(A)=9, ln(A)=10:"
do
	read (*,*,iostat=ReadStatus), ChosenOp
	if (ReadStatus.LE.0 .AND. ChosenOp.GE.1 .AND. ChosenOp.LE.10) exit
end do

call SelectColumn

if (ChosenOp.LT.7.AND.ChosenOp.NE.5) then
  print*, "  > Select the real constant: "
  do
	read (*,*,iostat=ReadStatus), RealConstant
	if (ReadStatus.LE.0) exit
  end do
end if

do XYear = 1, YearN
      Data(XYear,NextFree) = OpAwithB (Data(XYear,ChosenCol),RealConstant,ChosenOp)
end do

call StoreMoveOn
print*

end subroutine OperateAk

!*******************************************************************************
! operate A . B

subroutine OperateAB

print*, "  > Select the operation (A/B=1, A*B=2, A-B=3, A+B=4, A**B=6: "
do
	read (*,*,iostat=ReadStatus), ChosenOp
	if (ReadStatus.LE.0 .AND. ChosenOp.GE.1 .AND. ChosenOp.LE.6 .AND. ChosenOp.NE.5) exit
end do

print*, "  > Select columns A and B in turn: "
call SelectColumn
ChosenColA = ChosenCol
call SelectColumn
ChosenColB = ChosenCol

do XYear = 1, YearN
      Data(XYear,NextFree) = OpAwithB (Data(XYear,ChosenColA),Data(XYear,ChosenColB),ChosenOp)
end do

call StoreMoveOn
print*

end subroutine OperateAB

!*******************************************************************************
! check master array

subroutine CheckMaster

print "(a25)", "   > COL    MISS%    MEAN"

do XCol = 1, ColN
 DataMiss = 0.0 ; DataTot = 0.0 ; DataMean = MissVal

 do XYear = 1, YearN
      if (Data (XYear,XCol) .EQ. MissVal) then
      	DataMiss = DataMiss + 1.0
      else
        DataTot = DataTot + Data (XYear,XCol)
      end if
 end do

 DataFrac = 100.0 * DataMiss / YearN
 if ((YearN-DataMiss).GT.0) DataMean = DataTot / (YearN-DataMiss)
 
 print "(i9,f8.2,f8.2,a1,a9)", XCol, DataFrac, DataMean, " ", ColTitles(XCol) 
end do

print*

end subroutine CheckMaster

!*******************************************************************************
! calc equilibrium temperature

subroutine CalcEquT

print*, "  > Select cols A (global T anom.) and B (rad.forc.anom.) in turn: "
call SelectColumn
ChosenColA = ChosenCol
call SelectColumn
ChosenColB = ChosenCol

call GetKaySet (0,Forc2co2,FitAlpha,FitBeta,Sens2co2Init,TrendTLen)

print*, "  > Select the period length: "
do
	read (*,*,iostat=ReadStatus), PerLen
	if (ReadStatus.LE.0 .AND. PerLen.GE.1) exit
end do

allocate (LSRExe (TrendTLen), &
	  LSRWye (TrendTLen), &
	  OpVec  (YearN),     &
	  ResultVecA(YearN),  &
	  ResultVecB(YearN),  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcEquT: Allocation failure #####"

do XYear = 1, YearN								! calculates Forc
   OpVec(XYear) = Data(XYear,ChosenColB)
end do 
call GaussSmooth (YearN,PerLen,1,OpVec,ResultVecA,ResultVecB)
 
LastEffClimSens = Sens2co2Init
do XYear = 1, YearN
 TrendYear0 = nint (real(XYear)-(real(TrendTLen)/2.0))
 TrendYear1 = TrendYear0 + TrendTLen - 1

 LSRExe=MissVal ; LSRWye=MissVal
 do XTrendYear = TrendYear0, TrendYear1						! calculates DTDt
  if (XTrendYear.GE.1.AND.XTrendYear.LE.YearN) then          
          LSRExe (XTrendYear-TrendYear0+1) = XTrendYear-TrendYear0+1
          LSRWye (XTrendYear-TrendYear0+1) = Data(XTrendYear,ChosenColA)
  end if
 end do
 call LinearLSRVec (LSRExe,LSRWye,Trash1,DTDt,Trash2)
 
 Forc = ResultVecA(XYear)
 
 if (LastEffClimSens.NE.MissVal.AND.DTDt.NE.MissVal) then
   EffClimSens = LastEffClimSens + (DTDt * FitAlpha * exp (FitBeta * DTDt))	! calculates eff.clim.sens.
   LastEffClimSens = EffClimSens
 end if
 
 if (EffClimSens.NE.MissVal.AND.Forc.NE.MissVal) &
   Data(XYear,NextFree) = EffClimSens * (Forc / Forc2co2)			! calcs equilibrium T       
end do

deallocate (LSRExe,LSRWye,OpVec,ResultVecA,ResultVecB,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcEquT: Deallocation failure #####"

call StoreMoveOn
print*

end subroutine CalcEquT

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

subroutine PrintPeriodMean

print*, "  > Select first, last years AD in period: "
do
	read (*,*,iostat=ReadStatus), BaseADYear0,BaseADYear1
	if (ReadStatus.LE.0.AND.BaseADYear0.GE.YearAD(1).AND.BaseADYear1.LE.YearAD(YearN) &
				.AND.BaseADYear1.GE.BaseADYear0) exit
end do

do XCol = 1, ColN
  OpTot = 0.0 ; OpEn = 0.0 ; OpMean = MissVal
  
  do XYear = 1, YearN
    if (YearAD(XYear).GE.BaseADYear0.AND.YearAD(XYear).LE.BaseADYear1) then
      if (Data(XYear,XCol).NE.MissVal) then
        OpTot = OpTot + Data(XYear,XCol) ; OpEn = OpEn + 1
      end if
    end if
  end do
  
  if (OpEn.GE.1) OpMean = OpTot / OpEn
  
  print "(a,i3,2e12.4)", "     ", XCol, OpMean, OpEn
end do

end subroutine PrintPeriodMean

!*******************************************************************************
! select column

subroutine SelectColumn

print*, "  > Select the column: "
do
	read (*,*,iostat=ReadStatus), ChosenCol
	if (ReadStatus.LE.0 .AND. ChosenCol.GE.1 .AND. ChosenCol.LE.ColN) exit
end do

end subroutine SelectColumn

!*******************************************************************************
! store data in column and move on

subroutine StoreMoveOn

call CheckBlank

if (NonBlank.EQ.0) then
  print*, "  > Operation has no valid results. No results stored."
else
  print "(a13,i4,a9)", "   > Column: ", NextFree, adjustr(ColTitles(NextFree))
  
  NextFree = NextFree + 1
  if (NextFree.GT.ColN) then
    NextFree = NextFree - 1
    print*, "  > All columns are being used. Next free set to: ", NextFree
  end if
end if

end subroutine StoreMoveOn

!*******************************************************************************
! check whether column is (=0) or is not (=1) blank

subroutine CheckBlank

NonBlank = 0

XYear    = 0
do
  XYear = XYear + 1
  
  if (Data(XYear,NextFree).NE.MissVal) NonBlank = 1
  if (NonBlank.EQ.1) exit
  if (XYear.EQ.YearN) exit
end do

end subroutine CheckBlank

!*******************************************************************************
! save selected columns to .ann file

subroutine DumpToAnn

call CheckMaster

if (ColN.GT.1) then
  print*, "  > Enter the number of columns to save to .ann: "
  do
	read (*,*,iostat=ReadStatus), SaveColN
	if (ReadStatus.LE.0 .AND. SaveColN.GE.1 .AND. SaveColN.LE.ColN) exit
  end do
else
  SaveColN = 1
end if

allocate ( FileData      (YearN, SaveColN), &
	   FileColTitles (SaveColN),        stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpToAnn: Allocation failure #####"

if (SaveColN.LT.ColN) print "(a,i4)", "   > Select each column from array in turn: "
do XSaveCol = 1, SaveColN
  if (SaveColN.LT.ColN) then
    do
	read (*,*,iostat=ReadStatus), ChosenCol
	if (ReadStatus.GT.0 .OR. ChosenCol.LT.1 .OR. ChosenCol.GT.ColN) print*, "  > Bad entry. Try again."
	if (ReadStatus.LE.0 .AND. ChosenCol.GE.1 .AND. ChosenCol.LE.ColN) exit
    end do
  else
    ChosenCol = XSaveCol
  end if

  FileColTitles(XSaveCol) = ColTitles(ChosenCol)
  
  do XYear = 1, YearN
    FileData(XYear,XSaveCol) = Data(XYear,ChosenCol) 
  end do
end do

call MakeANNColTitles (FileColTitles)
call SaveANN (Blank, YearAD, FileColTitles, FileData)
print*

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

end subroutine DumpToAnn

!*******************************************************************************
! conclude

subroutine Conclude

deallocate (YearAD,Data,ColTitles,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Conclude: Deallocation failure #####"

close (99)

print*

end subroutine Conclude

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

end program OpAnn