! opday.f90
! f90 program written by Tim Mitchell on 29.01.01
! last modified on 29.04.03
! tool to operate on standard .day files
! pgf90 -o ./../obs/opday filenames.f90 annfiles.f90 synfiles.f90 dayfiles.f90 
!	stangeneral.f90 sortmod.f90 gales.f90 basicfun.f90 ./../obs/opday.f90

program OpDay

use FileNames
use ANNFiles
use SYNFiles
use DAYFiles
use StanGeneral
use SortMod
use Gales
use BasicFun

implicit none

real, pointer, dimension (:,:,:,:)	:: Data, FileFlow
real, pointer, dimension (:,:,:)	:: FileDaily
real, pointer, dimension (:,:)		:: AnnOutput
real, pointer, dimension (:)		:: ClassBounds

integer, pointer, dimension (:,:,:)	:: SeasonKey, FileLamb, FileAuto, FileJenk
integer, pointer, dimension (:,:)	:: MonthLengths, GaleThresh
integer, pointer, dimension (:)		:: YearAD, FileYearAD, BaseYearAD

character (len=9), pointer, dimension (:) 	:: AnnColTitles

real, parameter :: MissVal = -999.0

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

real :: MissAccept
real :: RealConstant, Speed, Force, ScreenVal
real :: DataMiss, DataFrac, DataMin, DataMax
real :: OpMiss, OpEn, OpFrac, OpThresh, OpDiff, OpTot, OpLen, OpBase, MaxLen, MissThresh

integer :: ReadStatus, AllocStat, MenuChoice, NextFree 
integer :: ChosenCol, ChosenColA, ChosenColB, ChosenColF, ChosenColZ, ChosenOp
integer :: StartYear, EndYear, FirstYear, YearN
integer :: BlockYearN, BlockStartAD, BlockEndAD, BlockStartX, BlockEndX
integer :: XYear, XMonth, XDay, XHeader, XFooter, XFileYear, XVari, XCol, XSevere, XClass, XBound
integer :: XGrowSeas,XJulian
integer :: HeaderN, FooterN, MonthN, DayN, VariN, StringLen, ColN, LoadColN, ClassN
integer :: FileYear0,FileYear1,Year0,Year1,Day0,Day1,Month0,Month1,Julian0,Julian1
integer :: QComp, QIntFloat, QFileComp, QUpDown
integer :: NonBlank, IntConstant, Overlap, VariCode, CurrentSeason
integer :: BegYear,BegMon,BegDay,BaseADYear0,BaseADYear1,BaseYearN,BaseYear0,BaseYear1
integer :: PrevEndYear,PrevEndMonth,PrevEndDay,PrevEndJulian,GrowSeasEnd

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
    print*, "  > Not yet written."
    print*
  else if (MenuChoice.EQ. 5)  then
    call AnomaliseA
  else if (MenuChoice.EQ. 6)  then
    call GetPerExceed
  else if (MenuChoice.EQ. 7)  then
    call ScreenExceed
  else if (MenuChoice.EQ.10)  then
    call LoadFromDAY
  else if (MenuChoice.EQ.11)  then
    call LoadFromSYN
  else if (MenuChoice.EQ.20)  then
    call SaveToDAY
  else if (MenuChoice.EQ.21)  then
    call SaveToSYN
  else if (MenuChoice.EQ.30)  then
    call CheckMaster
  else if (MenuChoice.EQ.40)  then
    call CalcFreqA
  else if (MenuChoice.EQ.41)  then
    call CalcPerKay
  else if (MenuChoice.EQ.42)  then
    call CalcGaleFreq
  else if (MenuChoice.EQ.43)  then
    call SeasonaliseA
  else if (MenuChoice.EQ.44)  then
    call SeasonaliseA
  else if (MenuChoice.EQ.45)  then
    call GrowSeasA
  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. Record periods <> k"
    print*, "  >  7. Screen out values <> k"
    print*, "  > 10. Load data from .day"
    print*, "  > 11. Load data from .syn"
    print*, "  > 20. Save data to   .day"
    print*, "  > 21. Save data to   .syn"
    print*, "  > 30. Check main arrays"
    print*, "  > 40. Classify A       --> .ann"
    print*, "  > 41. Longest per <> k --> .ann"
    print*, "  > 42. Gale freq        --> .ann"
    print*, "  > 43. Sum season of A  --> .ann"
    print*, "  > 44. Av. season of A  --> .ann"
    print*, "  > 45. Growing season A --> .ann"
    print*, "  > 99. Exit"
  end if
  
  if (MenuChoice.EQ.99) exit
end do

call Conclude

contains

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

subroutine Intro

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

print*
print*, "  > ##### OperateDay.f90 ##### Tool for operating on .day files #####"
print*

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
MonthN = 12 ; DayN = 31 ; VariN = 8
NextFree = 1 ; VariCode = MissVal ; MissAccept = 10.0

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

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

print*

end subroutine Intro

!*******************************************************************************
! load from .day file

subroutine LoadFromDAY

call LoadDAY (Blank, VariCode, FileYearAD, FileDaily)

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

do XFileYear = FileYear0, FileYear1
  XYear = Year0 + XFileYear - FileYear0
  
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      Data(NextFree,XYear,XMonth,XDay) = FileDaily(XFileYear,XMonth,XDay)
    end do
  end do
end do

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

call StoreMoveOn
print*

end subroutine LoadFromDAY

!*******************************************************************************
! load from .syn file

subroutine LoadFromSYN

call LoadSYN      (Blank,FileYearAD,FileLamb,FileJenk,FileAuto,FileFlow)

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

print*, "  > Variables: Lamb,Jenk,Auto,PM-1000,W,S,F,D,ZW,ZS,S"
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.11) exit
end do

do XCol = 1, LoadColN
  if (LoadColN.LT.11) then
    print "(a53,i4)", "   > Select column (1-11) from file --> array column: ", NextFree
    do
	read (*,*,iostat=ReadStatus), ChosenCol
	if (ReadStatus.LE.0 .AND. ChosenCol.GE.1 .AND. ChosenCol.LE.11) exit
    end do
  else
    ChosenCol = XCol
  end if
  
  do XFileYear = FileYear0, FileYear1
   XYear = Year0 + XFileYear - FileYear0
  
   do XMonth = 1, MonthN
    do XDay = 1, DayN
      if (ChosenCol.EQ.1) Data(NextFree,XYear,XMonth,XDay) = real(FileLamb(XFileYear,XMonth,XDay))
      if (ChosenCol.EQ.2) Data(NextFree,XYear,XMonth,XDay) = real(FileJenk(XFileYear,XMonth,XDay))
      if (ChosenCol.EQ.3) Data(NextFree,XYear,XMonth,XDay) = real(FileAuto(XFileYear,XMonth,XDay))
      if (ChosenCol.GE.4) Data(NextFree,XYear,XMonth,XDay) = FileFlow(XFileYear,XMonth,XDay,(ChosenCol-3))
    end do
   end do
  end do

  call StoreMoveOn  
end do

deallocate (FileYearAD,FileLamb,FileJenk,FileAuto,FileFlow, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadFromSYN: Deallocation failure #####"

print*

end subroutine LoadFromSYN

!*******************************************************************************
! save to .syn file

subroutine SaveToSYN

allocate ( FileFlow (YearN, MonthN, DayN, 8), &
	   FileLamb (YearN, MonthN, DayN),    &
	   FileJenk (YearN, MonthN, DayN),    &
	   FileAuto (YearN, MonthN, DayN),    stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveToSYN: Allocation failure #####"
FileFlow = MissVal
FileLamb = -9 ; FileJenk = -9 ; FileAuto = -9
print*, "  > Variables: Lamb,Jenk,Auto,PM-1000,W,S,F,D,ZW,ZS,S"

do XVari = 1, 11
  print "(a58,i4)", "   > Select column to use (-9=missing) for variable: ", XVari
  ChosenCol = -1
  do
	read (*,*,iostat=ReadStatus), ChosenCol
	if (ChosenCol.NE.-9) then
	  if (ChosenCol.LT.1.OR.ChosenCol.GT.ColN) ChosenCol = -1
	end if
	if (ReadStatus.LE.0 .AND. ChosenCol.NE.-1) exit
  end do

  if (ChosenCol.NE.-9) then
   do XYear = 1, YearN  
    do XMonth = 1, MonthN
     do XDay = 1, DayN
      if (XVari.EQ.1) FileLamb(XYear,XMonth,XDay) 		= nint(Data(ChosenCol,XYear,XMonth,XDay))
      if (XVari.EQ.2) FileJenk(XYear,XMonth,XDay) 		= nint(Data(ChosenCol,XYear,XMonth,XDay))
      if (XVari.EQ.3) FileAuto(XYear,XMonth,XDay) 		= nint(Data(ChosenCol,XYear,XMonth,XDay))
      if (XVari.GE.4) FileFlow(XYear,XMonth,XDay,(XVari-3)) 	= Data(ChosenCol,XYear,XMonth,XDay)
     end do
    end do
   end do    
  end if  
end do

call SaveSYN (Blank,YearAD,FileLamb,FileJenk,FileAuto,FileFlow)

deallocate (FileLamb,FileJenk,FileAuto,FileFlow, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveToSYN: Deallocation failure #####"

print*

end subroutine SaveToSYN

!*******************************************************************************
! save to .day file

subroutine SaveToDAY

call SelectColumn

allocate ( FileDaily (YearN, MonthN, DayN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveToDAY: Allocation failure #####"

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      FileDaily(XYear,XMonth,XDay) = Data(ChosenCol,XYear,XMonth,XDay)
    end do
  end do
end do

call SaveDAY (Blank,VariCode,YearAD,FileDaily)
print*

end subroutine SaveToDAY

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

subroutine OperateAk

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

call SelectColumn

if (ChosenOp.NE.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
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      Data(NextFree,XYear,XMonth,XDay) = OpAwithB (Data(ChosenCol,XYear,XMonth,XDay), &
      					RealConstant,ChosenOp)
    end do
  end do
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
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      Data(NextFree,XYear,XMonth,XDay) = OpAwithB &
      	(Data(ChosenColA,XYear,XMonth,XDay),Data(ChosenColB,XYear,XMonth,XDay),ChosenOp)
    end do
  end do
end do

call StoreMoveOn
print*

end subroutine OperateAB

!*******************************************************************************
! anomalise A against period from within itself

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

BaseYearN = BaseADYear1 - BaseADYear0 + 1
OpThresh  = real(BaseYearN) * (100.0 - MissAccept) / 100.0

allocate ( BaseYearAD (BaseYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AnomaliseA: Allocation failure #####"

do XYear = 1, BaseYearN
  BaseYearAD (XYear) = BaseADYear0 + XYear - 1
end do

call CommonVecPer (BaseYearAD,YearAD,BaseYear0,BaseYear1,Year0,Year1)		! identify common period
print*, "  > Period identified: ", BaseYearAD(BaseYear0), BaseYearAD(BaseYear1)

do XMonth = 1, MonthN
  do XDay = 1, DayN
    OpEn   = 0.0
    OpTot  = 0.0
    
    do XYear = Year0, Year1							! calc base
      if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
        OpEn  = OpEn  + 1
        OpTot = OpTot + Data(ChosenCol,XYear,XMonth,XDay)
      end if
    end do
    
    if (OpEn.GT.OpThresh) then							! calc anomalies
      OpBase = OpTot / OpEn
      
      do XYear = 1, YearN
        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) Data(NextFree,XYear,XMonth,XDay) = &
        				Data(ChosenCol,XYear,XMonth,XDay) - OpBase
      end do    
    else if (XMonth.EQ.2.AND.XDay.EQ.29.AND.OpBase.NE.MissVal) then		! 29th Feb anom v 28th Feb
      do XYear = 1, YearN
        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) Data(NextFree,XYear,XMonth,XDay) = &
        				Data(ChosenCol,XYear,XMonth,XDay) - OpBase
      end do    
    else
      OpBase = MissVal
    end if
  end do
end do

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

call StoreMoveOn
print*

end subroutine AnomaliseA

!*******************************************************************************
! locate periods where value exceeds constant

subroutine GetPerExceed

call SelectColumn
call GetMonthLengths (YearAD,MonthLengths)

print*, "  > Enter the threshold: "
do
	read (*,*,iostat=ReadStatus), OpThresh
	if (ReadStatus.LE.0) exit
end do

print*, "  > Find periods when threshold is (=1) or is not (=-1) exceeded: "
do
	read (*,*,iostat=ReadStatus), QUpDown
	if (ReadStatus.LE.0.AND.QUpDown.LE.1.AND.QUpDown.GE.-1) exit
end do

OpEn = 0 ; BegDay = 0 ; BegMon = 0 ; BegYear = 0

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, MonthLengths(XYear,XMonth)
      if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
        Data(NextFree,XYear,XMonth,XDay) = 0
        OpDiff = Data(ChosenCol,XYear,XMonth,XDay) - OpThresh
        
        if ((OpDiff*QUpDown).GE.0) then
          OpEn = OpEn + 1
          
          if (OpEn.EQ.1) then
            BegDay = XDay ; BegMon = XMonth ; BegYear = XYear
          end if
          
          Data (NextFree,BegYear,BegMon,BegDay) = OpEn
        else
          OpEn = 0 ; BegDay = 0 ; BegMon = 0 ; BegYear = 0
        end if
      else
          OpEn = 0 ; BegDay = 0 ; BegMon = 0 ; BegYear = 0
      end if
    end do
  end do
end do

print*, "  > We record the length of a period of exceedance on the initial day."
print*, "  > All other non-missing days are set to zero."
call StoreMoveOn
print*

end subroutine GetPerExceed

!*******************************************************************************
! screen out all values <> specified constant

subroutine ScreenExceed

call SelectColumn
call GetMonthLengths (YearAD,MonthLengths)

print*, "  > Screen out when threshold is (=1) or is not (=-1) exceeded: "
do
	read (*,*,iostat=ReadStatus), QUpDown
	if (ReadStatus.LE.0.AND.QUpDown.LE.1.AND.QUpDown.GE.-1) exit
end do

print*, "  > Enter the threshold: "
do
	read (*,*,iostat=ReadStatus), OpThresh
	if (ReadStatus.LE.0) exit
end do

print*, "  > Enter the value to place in screened cells: "
do
	read (*,*,iostat=ReadStatus), ScreenVal
	if (ReadStatus.LE.0) exit
end do

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, MonthLengths(XYear,XMonth)
      if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
        if      (QUpDown.EQ.-1.AND.Data(ChosenCol,XYear,XMonth,XDay).LT.OpThresh) then
        	Data(NextFree,XYear,XMonth,XDay) = ScreenVal        
        else if (QUpDown.EQ. 1.AND.Data(ChosenCol,XYear,XMonth,XDay).GT.OpThresh) then
        	Data(NextFree,XYear,XMonth,XDay) = ScreenVal
        else
                Data(NextFree,XYear,XMonth,XDay) = Data(ChosenCol,XYear,XMonth,XDay)
        end if
      end if
    end do
  end do
end do

call StoreMoveOn
print*

end subroutine ScreenExceed

!*******************************************************************************
! calculate frequency of value in specified season

subroutine CalcFreqA

call SelectColumn
call MakeSeasonKey (YearAD,SeasonKey)

print*, "  > Enter the number of classes: "
do
	read (*,*,iostat=ReadStatus), ClassN
	if (ClassN.LT.2) print*, "  > Number of classes must exceed one. Try again."
	if (ReadStatus.LE.0 .AND. ClassN.GE.2) exit
end do

allocate (ClassBounds  (ClassN-1),          &
 	  AnnOutput    (YearN, (ClassN+1)), &
	  AnnColTitles (ClassN+1),          stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcFreqA: Allocation failure #####"
AnnOutput    = 0 ; AnnColTitles = '        X' ; ClassBounds  = MissVal
AnnColTitles(ClassN+1)='  MISSING' ; AnnColTitles(1)='   LOWEST' ; AnnColTitles(ClassN)='  HIGHEST'

print*, "  > Enter the bounds between the classes, in ascending order: "
do XBound = 1, (ClassN-1)
  do
	read (*,*,iostat=ReadStatus), ClassBounds(XBound)
	if (ReadStatus.GT.0) print*, "  > Not a real. Re-enter."
	if (XBound.GT.1.AND.ClassBounds(XBound).LE.ClassBounds(XBound-1)) then
		print*, "  > Too low. Re-enter."
		ClassBounds(XBound) = MissVal
	end if
	if (ReadStatus.LE.0 .AND. ClassBounds(XBound).NE.MissVal) exit
  end do
end do

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      if (SeasonKey(XYear,XMonth,XDay).NE.MissVal) then
        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
          XBound = 0
          do
            XBound = XBound + 1
            
            if (Data(ChosenCol,XYear,XMonth,XDay).LT.ClassBounds(XBound)) &
            		AnnOutput(SeasonKey(XYear,XMonth,XDay),XBound) = &
          		AnnOutput(SeasonKey(XYear,XMonth,XDay),XBound) + 1

            if (Data(ChosenCol,XYear,XMonth,XDay).LT.ClassBounds(XBound)) exit            
            if (XBound.EQ.(ClassN-1)) exit
          end do
          
          if (Data(ChosenCol,XYear,XMonth,XDay).GE.ClassBounds(ClassN-1)) &
            		AnnOutput(SeasonKey(XYear,XMonth,XDay),ClassN) = &
          		AnnOutput(SeasonKey(XYear,XMonth,XDay),ClassN) + 1

        else
          AnnOutput(SeasonKey(XYear,XMonth,XDay),(ClassN+1)) = &
          		AnnOutput(SeasonKey(XYear,XMonth,XDay),(ClassN+1)) + 1
        end if
      end if
    end do
  end do
end do

do XYear = 1, YearN
  OpEn = 0
  
  do XClass = 1, ClassN
      OpEn = OpEn + AnnOutput (XYear,XClass)
  end do
  
  MissThresh = (MissAccept / 100.0) * OpEn
  
  if (AnnOutput(XYear,(ClassN+1)).GE.MissThresh) then
    do XClass = 1, ClassN
      AnnOutput (XYear,XClass) = MissVal
    end do  
  end if
end do

call MakeANNColTitles (AnnColTitles)
call SaveANN          (Blank, YearAD, AnnColTitles, AnnOutput)
print*

deallocate (AnnOutput,AnnColTitles,ClassBounds,SeasonKey, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcFreqA: Deallocation failure #####"

end subroutine CalcFreqA

!*******************************************************************************
! calculate length of longest period that satisfies the following requirements: 
! begin with first period where meanT>5degK for 5 consecutive days
! end with first subsequent period when meanT<5degK for 5 consecutive days

subroutine GrowSeasA

call SelectColumn
call GetMonthLengths (YearAD,MonthLengths)

allocate (AnnOutput    (YearN,5), &
	  AnnColTitles (5),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrowSeasA: Allocation failure #####"
AnnOutput = 0 
AnnColTitles(1) = '   LENGTH' ; AnnColTitles(2) = '   JULIAN' ; AnnColTitles(3) = '    BLANK'
AnnColTitles(4) = ' JULI-BEG' ; AnnColTitles(5) = ' JULI-END'

print*, "  > Enter the temperature threshold (e.g. 5.0 degC): "
do
	read (*,*,iostat=ReadStatus), OpThresh
	if (ReadStatus.LE.0) exit
end do

print*, "  > Enter the length threshold (e.g. 5 days): "
do
	read (*,*,iostat=ReadStatus), OpLen
	if (ReadStatus.LE.0) exit
end do

do XYear = 1, YearN
  if (MonthLengths(XYear,2).EQ.28) then
    AnnOutput (XYear,2) = 365
  else
    AnnOutput (XYear,2) = 366
  end if
end do

OpEn = 0.0 ; OpTot = 0.0				! current mini-per length, current grow-seas length
PrevEndYear = 1 ; PrevEndMonth = 1 ; PrevEndDay = 1	! the lastday of last growing season
PrevEndJulian = 1					! the last Julian day of last growing season
  
do XGrowSeas = 1, YearN
  GrowSeasEnd = 0							! =1 when end of grow seas is found
  XYear=PrevEndYear ; XMonth=PrevEndMonth ; XDay=PrevEndDay		! start at end of last growseas
  XJulian=PrevEndJulian
  
  do									! year loop    
    do									! month loop
      do								! day loop
        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
          if (OpTot.EQ.0) then						! no grow-seas at present
            if (Data(ChosenCol,XYear,XMonth,XDay).GT.OpThresh) then		! potential mini-per
              OpEn = OpEn + 1            
              
              if (OpEn.EQ.OpLen) then						! this = potential grow-seas
                OpTot = OpLen ; OpEn = 0
              end if
            else								! no mini-per at present
              OpEn = 0
            end if
          else								! potential grow-seas in progress
            OpTot = OpTot + 1							! len of potential grow-seas
            
            if (Data(ChosenCol,XYear,XMonth,XDay).LT.OpThresh) then		! potential mini-per
              OpEn = OpEn + 1            
              
              if (OpEn.EQ.OpLen) then						! this = end of potential g-s
                if (OpTot.GT.AnnOutput(XGrowSeas,1)) then			! this = growing season
                  AnnOutput(XGrowSeas,1) = OpTot - OpLen
                  AnnOutput(XGrowSeas,4) = XJulian + ((XYear-XGrowSeas)*AnnOutput(XGrowSeas,2)) - OpTot + 1
                  AnnOutput(XGrowSeas,5) = XJulian + ((XYear-XGrowSeas)*AnnOutput(XGrowSeas,2)) - OpLen
                  
                  PrevEndYear=XYear ; PrevEndMonth=XMonth ; PrevEndDay=XDay ; PrevEndJulian=XJulian
                end if
                                  
                if (XYear.GT.XGrowSeas) GrowSeasEnd = 1				! if end + nextyr -> exit 

                OpTot = 0 ; OpEn = 0						! reset mini-per,potential g-s
              end if
            else								! no mini-per at present
              OpEn = 0
            end if
          end if
        end if
        
        XDay = XDay + 1 ; XJulian = XJulian + 1
        if (XDay.GT.MonthLengths(XYear,XMonth)) XDay = 1		! end of month is reached
        if (XDay.EQ.1) exit
        if (GrowSeasEnd.EQ.1) exit
      end do
      
      XMonth = XMonth + 1      
      if (XMonth.GT.12) XMonth = 1					! end of year is reached    
      if (XMonth.EQ.1) exit    
      if (GrowSeasEnd.EQ.1) exit
    end do
    
    if (XYear.EQ.XGrowSeas) then					! current year = year of growing seas
      if (OpTot.EQ.0) GrowSeasEnd = 1					! not current grow seas -> exit
    end if
    
    XYear = XYear + 1
    XJulian = 1								! resets Julian day    
    
    if (XYear.GT.YearN.AND.OpTot.GT.0) then				! if end-of-file prior to end of g-s
      if (OpTot.GT.AnnOutput(XGrowSeas,1)) then				! if this was gonna be g-s
        AnnOutput(XGrowSeas,1) = MissVal
        AnnOutput(XGrowSeas,4) = MissVal
	AnnOutput(XGrowSeas,5) = MissVal
      end if
    end if
    
    if (XYear.GT.YearN) exit						! end of data file is reached					
    if (GrowSeasEnd.EQ.1) exit
  end do
end do      
      
call SaveANN (Blank, YearAD, AnnColTitles, AnnOutput)
print*

deallocate (AnnOutput,AnnColTitles,MonthLengths, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrowSeasA: Deallocation failure #####"

end subroutine GrowSeasA

!*******************************************************************************
! calculate longest period in days in specified season when value exceeds/stays-below Kay

subroutine CalcPerKay

call SelectColumn
call MakeSeasonKey (YearAD,SeasonKey)

allocate (AnnOutput    (YearN,5), &
	  AnnColTitles (5),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcPerKay: Allocation failure #####"
AnnOutput = 0 
AnnColTitles(1) = '  MAX LEN' ; AnnColTitles(2) = '   PERIOD' ; AnnColTitles(3) = '  MISSING'
AnnColTitles(4) = 'DAY0(SEA)' ; AnnColTitles(5) = 'DAY1(SEA)'

print*, "  > Enter the threshold: "
do
	read (*,*,iostat=ReadStatus), OpThresh
	if (ReadStatus.LE.0) exit
end do

print*, "  > Find longest period when threshold is (=1) or is not (=2) exceeded: "
do
	read (*,*,iostat=ReadStatus), QUpDown
	if (ReadStatus.LE.0.AND.QUpDown.GE.1.AND.QUpDown.LE.2) exit
end do

CurrentSeason = 1 ; OpEn = 0

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      if (SeasonKey(XYear,XMonth,XDay).GT.CurrentSeason.AND.SeasonKey(XYear,XMonth,XDay).NE.MissVal) then
        OpEn   = 0        
        CurrentSeason = SeasonKey(XYear,XMonth,XDay)
      end if
      
      if (SeasonKey(XYear,XMonth,XDay).EQ.CurrentSeason) then
        AnnOutput (CurrentSeason,2) = AnnOutput (CurrentSeason,2) + 1.0

        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
          if      (QUpDown.EQ.1.AND.Data(ChosenCol,XYear,XMonth,XDay).GE.OpThresh) then
          	OpEn = OpEn + 1
          else if (QUpDown.EQ.2.AND.Data(ChosenCol,XYear,XMonth,XDay).LE.OpThresh) then
          	OpEn = OpEn + 1
          else
          	OpEn = 0
          end if
        else
          AnnOutput (CurrentSeason,3) = AnnOutput (CurrentSeason,3) + 1.0
        end if        
      end if
      
      if (OpEn.GT.AnnOutput(CurrentSeason,1)) then
      	AnnOutput (CurrentSeason,1) = OpEn 
      	AnnOutput (CurrentSeason,4) = AnnOutput (CurrentSeason,2) - OpEn + 1 
      	AnnOutput (CurrentSeason,5) = AnnOutput (CurrentSeason,2) 
      end if
    end do
  end do
end do

do XYear = 1, YearN
  MissThresh = (MissAccept / 100.0) * AnnOutput(XYear,2)
  if (AnnOutput(XYear,3).GT.MissThresh) AnnOutput (CurrentSeason,1) = MissVal
end do
	
call SaveANN (Blank, YearAD, AnnColTitles, AnnOutput)
print*

deallocate (AnnOutput,AnnColTitles,SeasonKey, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcPerKay: Deallocation failure #####"

end subroutine CalcPerKay

!*******************************************************************************
! calculate gale frequencies for specified season

subroutine CalcGaleFreq

print*, "  > To measure gales we require columns with variables F, Z."
call SelectColumn
ChosenColF = ChosenCol
call SelectColumn
ChosenColZ = ChosenCol

call GetGaleThresh (YearAD,GaleThresh)
call MakeSeasonKey (YearAD,SeasonKey)

allocate (AnnOutput (YearN, 5), &
	  AnnColTitles (5),     stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcGaleFreq: Allocation failure #####"
AnnOutput = 0
AnnColTitles(1)='     GALE' ; AnnColTitles(2)='   SEVERE' ; AnnColTitles(3)=' V SEVERE' 
AnnColTitles(4)='   LENGTH' ; AnnColTitles(5)='  MISSING'

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      if (SeasonKey(XYear,XMonth,XDay).NE.MissVal) then
        AnnOutput(SeasonKey(XYear,XMonth,XDay),4) = AnnOutput(SeasonKey(XYear,XMonth,XDay),4) + 1
        
        Speed = GaleIndex (Data(ChosenColF,XYear,XMonth,XDay),Data(ChosenColZ,XYear,XMonth,XDay))
        Force = GaleScore (Speed,GaleThresh(XYear,1),GaleThresh(XYear,2),GaleThresh(XYear,3))
        
        if (Force.NE.MissVal) then
          do XSevere = 1, 3
            if (Force.GE.XSevere) AnnOutput(SeasonKey(XYear,XMonth,XDay),XSevere) = &
          		AnnOutput(SeasonKey(XYear,XMonth,XDay),XSevere) + 1
          end do
        else
          AnnOutput(SeasonKey(XYear,XMonth,XDay),5) = AnnOutput(SeasonKey(XYear,XMonth,XDay),5) + 1
        end if
      end if
    end do
  end do
end do

do XYear = 1, YearN
  MissThresh = (MissAccept / 100.0) * AnnOutput(XYear,4)
  if (AnnOutput(XYear,5).GE.MissThresh) then
  	AnnOutput (XYear,1) = MissVal
  	AnnOutput (XYear,2) = MissVal
  	AnnOutput (XYear,3) = MissVal
  end if
end do
	
call SaveANN          (Blank, YearAD, AnnColTitles, AnnOutput)
print*

deallocate (AnnOutput,AnnColTitles,GaleThresh,SeasonKey, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcGaleFreq: Deallocation failure #####"

end subroutine CalcGaleFreq

!*******************************************************************************
! calculate the sum (=43) or average (=44) of values in specified season

subroutine SeasonaliseA

if (MenuChoice.NE.43.AND.MenuChoice.NE.44) print*, "  > ##### ERROR: SeasonaliseA: Menu-subroutine mismatch"

call SelectColumn
call MakeSeasonKey (YearAD,SeasonKey)

allocate (AnnOutput    (YearN, 3), &
	  AnnColTitles (3),        stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeasonaliseA: Allocation failure #####"
AnnOutput    = 0
if (MenuChoice.EQ.43) AnnColTitles (1) = ' seas-sum' 
if (MenuChoice.EQ.44) AnnColTitles (1) = '  seas-av' 
AnnColTitles (2) = '   length' ; AnnColTitles (3) = '  missing'

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, DayN
      if (SeasonKey(XYear,XMonth,XDay).NE.MissVal) then
        AnnOutput(SeasonKey(XYear,XMonth,XDay),2)   = AnnOutput(SeasonKey(XYear,XMonth,XDay),2) + 1

        if (Data(ChosenCol,XYear,XMonth,XDay).NE.MissVal) then
          AnnOutput(SeasonKey(XYear,XMonth,XDay),1) = AnnOutput(SeasonKey(XYear,XMonth,XDay),1) + &
          		Data(ChosenCol,XYear,XMonth,XDay)
        else
          AnnOutput(SeasonKey(XYear,XMonth,XDay),3) = AnnOutput(SeasonKey(XYear,XMonth,XDay),3) + 1
        end if
      end if
    end do
  end do
end do

do XYear = 1, YearN
  MissThresh = (MissAccept / 100.0) * AnnOutput(XYear,2)
  
  if (AnnOutput(XYear,3).GE.MissThresh) then
  	AnnOutput (XYear,1) = MissVal
  else if (MenuChoice.EQ.44) then
  	AnnOutput (XYear,1) = AnnOutput (XYear,1) / (AnnOutput (XYear,2) - AnnOutput (XYear,3))
  end if
end do

call MakeANNColTitles (AnnColTitles)
call SaveANN          (Blank, YearAD, AnnColTitles, AnnOutput)
print*

deallocate (AnnOutput,AnnColTitles,SeasonKey, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeasonaliseA: Deallocation failure #####"

end subroutine SeasonaliseA

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

subroutine CheckMaster

call GetMonthLengths (YearAD,MonthLengths)

OpEn     = 0.0
do XYear = 1, YearN
  do XMonth = 1, MonthN
    OpEn = OpEn + MonthLengths(XYear,XMonth)
  end do
end do

print "(a38,f10.2)", "   > Number of days in main array:    ",   OpEn
print "(a17)",       "   > COL    MISS%     MIN     MAX"

do XCol = 1, ColN
 DataMiss =  0.0
 DataMin  =  100000.0
 DataMax  = -100000.0
 
 do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XDay = 1, MonthLengths(XYear,XMonth)
      if (Data (XCol,XYear,XMonth,XDay).EQ.MissVal) DataMiss = DataMiss + 1.0
      if (Data (XCol,XYear,XMonth,XDay).LT.DataMin) DataMin  = Data (XCol,XYear,XMonth,XDay)
      if (Data (XCol,XYear,XMonth,XDay).GT.DataMax) DataMax  = Data (XCol,XYear,XMonth,XDay)
    end do
  end do
 end do

 DataFrac = 100.0 * DataMiss /  OpEn
 
 print "(i9,3f8.2)", XCol, DataFrac, DataMin, DataMax 
end do

print*

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

end subroutine CheckMaster

!*******************************************************************************
! 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*, "  > Results stored in column: ", 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
  
  XMonth = 0
  do
    XMonth = XMonth + 1
    
    XDay = 0
    do
      XDay = XDay + 1
    
      if (Data(NextFree,XYear,XMonth,XDay).NE.MissVal) NonBlank = 1
      if (NonBlank.EQ.1) exit
      if (XDay.EQ.DayN) exit
    end do
  
    if (NonBlank.EQ.1) exit
    if (XMonth.EQ.MonthN) exit
  end do
  
  if (NonBlank.EQ.1) exit
  if (XYear.EQ.YearN) exit
end do

end subroutine CheckBlank

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

subroutine Conclude

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

close (99)

print*

end subroutine Conclude

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

end program OpDay