! annbulk.f90
! performs basic operations on .ann files in bulk
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../obs/annbulk time.f90 filenames.f90 annfiles.f90 basicfun.f90 
!	./../obs/annbulk.f90

program AnnBulk

use Time
use FileNames
use AnnFiles
use BasicFun

implicit none

!*******************************************************************************
! definitions

real, pointer, dimension (:,:) 	   		:: InnData,MidData,OutData
real, pointer, dimension (:) 	   		:: InnVec,MidVec,OutVec

integer, pointer, dimension (:) 	   	:: YearAD,InnYearAD,MidYearAD,OutYearAD
integer, pointer, dimension (:) 	   	:: Prune

character (len=80), pointer, dimension (:) 	:: InnFile,MidFile,OutFile
character (len= 9), pointer, dimension (:) 	:: Labels,InnLabels,MidLabels,OutLabels

real, parameter	:: MissVal = -999.0

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

logical :: IfPrune

real :: MagMin,MagMax,OpConstant

integer :: ReadStatus,AllocStat
integer :: QOperation,QOpPerform,QOutMag,QReplicate,QDecPlace,QKeepMiss
integer :: NCol,NInnCol,NMidCol,NOutCol,NYear,NInnYear,NMidYear,NOutYear
integer :: XCol,XInnCol,XMidCol,XOutCol,XYear,XInnYear,XMidYear,XOutYear
integer :: NFile,NMidFile,NPrune
integer :: XFile,XMidFile,XPrune
integer :: InnSubLen,InnSubBeg,InnFileLen,LastSlash,THalf
integer :: InnBeg,InnEnd,MidBeg,MidEnd,BegYearAD,EndYearAD

character (len=80) :: InnSub,OutSub,GivenFile,NameOnly
character (len= 9) :: InnText,MidText,OutText

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

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

print*
print*, "  > ***** AnnBulk: performs basic ops on .ann files in bulk *****"
print*

call SelectOp

print*

close (99)

contains

!*******************************************************************************
! main selection

subroutine SelectOp

print*, "  > Select the operation to perform (0=list): "
do
    read (*,*,iostat=ReadStatus), QOperation
    if (QOperation.EQ.0) then
      print*, "  >  7. smooth with Gaussian filter"
      print*, "  >  8. operate .ann on constant"
      print*, "  >  9. operate .ann on .ann"
      print*, "  > 10. prune columns"
      print*, "  > 11. restrict period"
    end if
    if (ReadStatus.LE.0.AND.QOperation.NE.0) exit
end do

if (QOperation.EQ. 7) call SmoothGaussian
if (QOperation.EQ. 8) call OpAnnConstant
if (QOperation.EQ. 9) call OpMultiAnn
if (QOperation.EQ.10) call PruneCols
if (QOperation.EQ.11) call RestrictPeriod

end subroutine SelectOp

!*******************************************************************************
!  7. smooth with Gaussian filter

subroutine SmoothGaussian

print*, "  > Select the periodicity at which to smooth: "
do
	read (*,*,iostat=ReadStatus), THalf
	if (ReadStatus.LE.0 .AND. THalf.GE.1) 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

call QueryOutMag
call GetInnFile
call GetOutFile

do XFile = 1, NFile
  call GetNameOnly
  call LoadANN (InnFile(XFile),InnYearAD,InnLabels,InnData)
  
  NCol=size(InnLabels,1) ; NYear=size(InnYearAD,1)
  allocate (OutData(NYear,NCol), &
  	    InnVec(NYear),MidVec(NYear),OutVec(NYear),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SmoothGaussian: Alloc #####"
  OutData=MissVal
  
  do XCol=1,NCol
    do XYear=1,NYear
      InnVec(XYear) = InnData(XYear,XCol)
    end do
    
    call GaussSmooth (NYear,THalf,QKeepMiss,InnVec,OutVec,MidVec)
    
    do XYear=1,NYear
      OutData(XYear,XCol) = OutVec(XYear)
    end do
  end do
  
  if (QOutMag.EQ.1) call CheckOutMag
    
  call SaveANN (OutFile(XFile),InnYearAD,InnLabels,OutData,DecPlaceN=QDecPlace)
    
  print "(2a,3(a,i4))", "   > ", trim(NameOnly), &
    		" for", NCol, " in ", InnYearAD(1), "-", InnYearAD(NYear)

  deallocate (OutData,InnVec,MidVec,OutVec,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SmoothGaussian: Dealloc #####"
end do

end subroutine SmoothGaussian

!*******************************************************************************
!  8. operate .ann on constant

subroutine OpAnnConstant

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

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

call QueryOutMag
call GetInnFile
call GetOutFile

do XFile = 1, NFile
  call GetNameOnly
  call LoadANN (InnFile(XFile),YearAD,Labels,InnData)
  
  NCol=size(Labels,1) ; NYear=size(YearAD,1)
  allocate (OutData(NYear,NCol),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpAnnConstant: Alloc #####"
  OutData=MissVal
  
  do XCol=1,NCol
      do XYear=1,NYear
	OutData(XYear,XCol) = OpAwithB (InnData(XYear,XCol), &
					OpConstant,QOpPerform)
      end do
  end do
  
  if (QOutMag.EQ.1) call CheckOutMag
    
  call SaveANN (OutFile(XFile),YearAD,Labels,OutData,DecPlaceN=QDecPlace)
    
  print "(2a,3(a,i4))", "   > ", trim(NameOnly), &
    		" for", NCol, " in ", YearAD(1), "-", YearAD(NYear)

  deallocate (OutData,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpAnnConstant: Dealloc #####"
end do

end subroutine OpAnnConstant

!*******************************************************************************
!  9. operate .ann on .ann

subroutine OpMultiAnn

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

call QueryOutMag
call GetParallelFiles
call GetOutFile

do XFile = 1, NFile
  call GetMidFileIndex
  call GetNameOnly
  call LoadANN (InnFile(XFile),   InnYearAD,InnLabels,InnData)
  call LoadANN (MidFile(XMidFile),MidYearAD,MidLabels,MidData)
  
  NInnCol=size(InnLabels,1) ; NMidCol=size(MidLabels,1) 
  NOutCol=NInnCol ; if (NMidCol.LT.NOutCol) NOutCol=NMidCol
  
  call CommonVecPer (InnYearAD,MidYearAD,InnBeg,InnEnd,MidBeg,MidEnd)
  
  if (InnBeg.NE.MissVal) then
    NOutYear=InnEnd-InnBeg+1
    
    allocate (OutYearAD(NOutYear),OutLabels(NOutCol),OutData(NOutYear,NOutCol), &
    	      stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiAnn: Alloc #####"
    OutYearAD=MissVal ; OutLabels="" ; OutData=MissVal
    
    do XOutYear=1,NOutYear
      XInnYear=InnBeg+XOutYear-1
      OutYearAD(XOutYear)=InnYearAD(XInnYear)
    end do
    
    do XOutCol=1,NOutCol
      InnText=trim(adjustl(InnLabels(XOutCol))) 
      MidText=trim(adjustl(MidLabels(XOutCol)))
      OutLabels(XOutCol) = " " // InnText(1:4) // MidText(1:4)
      OutLabels(XOutCol) = adjustr(OutLabels(XOutCol))
    end do
    
    do XOutCol=1,NOutCol
      do XOutYear=1,NOutYear
        XInnYear=InnBeg+XOutYear-1 ; XMidYear=MidBeg+XOutYear-1
	OutData(XOutYear,XOutCol) = OpAwithB (InnData(XInnYear,XOutCol), &
					MidData(XMidYear,XOutCol),QOpPerform)
      end do
    end do
    
    NYear=NOutYear ; NCol=NOutCol
    if (QOutMag.EQ.1) call CheckOutMag
    
    call SaveANN (OutFile(XFile),OutYearAD,OutLabels,OutData,DecPlaceN=QDecPlace)
    
    print "(2a,5(a,i4))", "   > ", trim(NameOnly), &
    		" for", NOutCol, " from", NInnCol, ",", NMidCol, &
		" in ", OutYearAD(1), "-", OutYearAD(NOutYear)

    deallocate (OutYearAD,OutLabels,OutData,stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiAnn: Dealloc #####"
  else
    print "(3a)", "   > ", trim(NameOnly), " has no common period"
  end if
end do

end subroutine OpMultiAnn

!*******************************************************************************
! 10. prune columns

subroutine PruneCols

print*, "  > Select the number of columns to prune: "
do
      read (*,*,iostat=ReadStatus), NPrune
      if (ReadStatus.LE.0.AND.NPrune.GE.1) exit
end do

allocate (Prune(NPrune),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PruneCols: Alloc #####"

do XPrune=1,NPrune
  print "(a,i2)", "   > Select a column to prune: ", XPrune
  do
      read (*,*,iostat=ReadStatus), Prune(XPrune)
      if (ReadStatus.LE.0.AND.Prune(XPrune).GE.1) exit
  end do
end do

call GetInnFile
call GetOutFile

do XFile = 1, NFile
  call GetNameOnly
  call LoadANN (InnFile(XFile),InnYearAD,InnLabels,InnData)
  
  NInnCol=size(InnLabels,1) ; NYear=size(InnYearAD,1) ; NOutCol=NInnCol
  do XInnCol=1,NInnCol
    XPrune=0
    do
      XPrune=XPrune+1
      if (Prune(XPrune).EQ.XInnCol) NOutCol=NOutCol-1
      if (XPrune.EQ.NPrune.OR.Prune(XPrune).EQ.XInnCol) exit
    end do
  end do
  
  allocate (OutData(NYear,NOutCol),OutLabels(NOutCol),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PruneCols: Alloc 2 #####"
  OutData=MissVal ; OutLabels=""
  
  XOutCol=0
  do XInnCol=1,NInnCol
    XPrune=0 ; IfPrune=.FALSE.
    do
      XPrune=XPrune+1
      if (Prune(XPrune).EQ.XInnCol) IfPrune=.TRUE.
      if (XPrune.EQ.NPrune.OR.Prune(XPrune).EQ.XInnCol) exit
    end do
    
    if (.NOT.IfPrune) then
      XOutCol=XOutCol+1
      OutLabels(XOutCol)=InnLabels(XInnCol)
      do XYear=1,NYear
        OutData(XYear,XOutCol)=InnData(XYear,XInnCol)
      end do
    end if
  end do
  
  call SaveANN (OutFile(XFile),InnYearAD,OutLabels,OutData,DecPlaceN=QDecPlace)
    
  print "(2a,4(a,i4))", "   > ", trim(NameOnly), &
    		" for", NOutCol, " of", NInnCol, " in ", InnYearAD(1), "-", InnYearAD(NYear)

  deallocate (OutData,OutLabels,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PruneCols: Dealloc 2 #####"
end do

deallocate (Prune,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PruneCols: Dealloc #####"

end subroutine PruneCols

!*******************************************************************************
! 11. restrict period

subroutine RestrictPeriod

print*, "  > Select the year AD beg,end (-999 for no limit): "
do
      read (*,*,iostat=ReadStatus), BegYearAD,EndYearAD
      if (ReadStatus.LE.0) exit
end do

call GetInnFile
call GetOutFile

do XFile = 1, NFile
  call GetNameOnly
  call LoadANN (InnFile(XFile),InnYearAD,Labels,InnData)
  
  NInnYear=size(InnYearAD,1) ; NCol=size(Labels,1)
  
  if      (BegYearAD.EQ.MissVal.AND.EndYearAD.NE.MissVal) then
    if (InnYearAD(1).LE.EndYearAD) 		BegYearAD=InnYearAD(1)
  else if (BegYearAD.NE.MissVal.AND.EndYearAD.EQ.MissVal) then
    if (InnYearAD(NInnYear).GE.BegYearAD) 	EndYearAD=InnYearAD(NInnYear)
  else if (BegYearAD.EQ.MissVal.AND.EndYearAD.EQ.MissVal) then
    BegYearAD=InnYearAD(1) ; EndYearAD=InnYearAD(NInnYear)
  end if
  
  if (BegYearAD.NE.MissVal.AND.EndYearAD.NE.MissVal) then
    NMidYear=EndYearAD-BegYearAD+1
    allocate (MidYearAD(NMidYear),stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: RestrictPeriod: Alloc 1 #####"
    do XMidYear=1,NMidYear
      MidYearAD(XMidYear) = BegYearAD + XMidYear - 1
    end do
    call CommonVecPer (InnYearAD,MidYearAD,InnBeg,InnEnd,MidBeg,MidEnd)
    deallocate (MidYearAD,stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: RestrictPeriod: Dealloc 1 #####"
  else
    InnBeg=MissVal
  end if
  
  if (InnBeg.NE.MissVal) then
    NOutYear=MidEnd-MidBeg+1
    allocate (OutData(NOutYear,NCol),OutYearAD(NOutYear), stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: RestrictPeriod: Alloc 2 #####"
    
    do XOutYear=1,NOutYear
      OutYearAD(XOutYear)=BegYearAD+MidBeg+XOutYear-2
    end do
      
    do XOutYear=1,NOutYear
      XInnYear=InnBeg+XOutYear-1
      do XCol=1,NCol
        OutData(XOutYear,XCol)=InnData(XInnYear,XCol)
      end do
    end do
    
    call SaveANN (OutFile(XFile),OutYearAD,Labels,OutData,DecPlaceN=QDecPlace)
    
    print "(2a,4(a,i4))", "   > ", trim(NameOnly), &
    		" has ",  OutYearAD(1), "-", OutYearAD(NOutYear), &
		" from ", InnYearAD(1), "-", InnYearAD(NInnYear)

    deallocate (OutData,OutYearAD, stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: RestrictPeriod: Dealloc 2 #####"
  else
    print "(3a)", "   > ", trim(NameOnly), " has no common period"
  end if
end do

end subroutine RestrictPeriod

!*******************************************************************************
! query whether to check output magnitudes

subroutine QueryOutMag

print*, "  > Restrict output magnitudes (0=no,1=yes) ?"
do
      read (*,*,iostat=ReadStatus), QOutMag
      if (ReadStatus.LE.0.AND.QOutMag.GE.0.AND.QOutMag.LE.1) exit
end do
if (QOutMag.EQ.1) then
  print*, "  > Magnitudes smaller than X are set to X. Enter X."
  do
      read (*,*,iostat=ReadStatus), MagMin
      if (ReadStatus.LE.0.AND.MagMin.GE.0) exit
  end do
  
  print*, "  > Magnitudes larger  than Y are set to Y. Enter Y."
  do
      read (*,*,iostat=ReadStatus), MagMax
      if (ReadStatus.LE.0.AND.MagMax.GE.0) exit
  end do  
end if

end subroutine QueryOutMag

!*******************************************************************************
! check output magnitudes

subroutine CheckOutMag

do XCol=1,NOutCol
  do XYear=1,NOutYear
    if (OutData(XYear,XCol).NE.MissVal) then
      if (abs(OutData(XYear,XCol)).LT.MagMin) &
      		OutData(XYear,XCol) = sign(MagMin,OutData(XYear,XCol))
      if (abs(OutData(XYear,XCol)).GT.MagMax) &
      		OutData(XYear,XCol) = sign(MagMax,OutData(XYear,XCol))
    end if
  end do
end do

end subroutine CheckOutMag

!*******************************************************************************
! get parallel files

subroutine GetParallelFiles

print*, "  > Select OLD files 'A'. There must be a common substring."
call GetBatch (Blank,InnFile)
NFile  = size (InnFile, 1)

print*, "  > Select OLD files 'B'. There must be a common substring."
do
  call GetBatch (Blank,MidFile)
  NMidFile = size(MidFile, 1)
  if (NFile.NE.NMidFile) then
    if (mod(real(NFile),real(NMidFile)).EQ.0) then
      print*, "  > Replicate the files to match A ? (0=no,1,2=yes)"
      print*, "  >   1=(a:a,n+a,...) 2=(a:n*a+1,n*a+2,...)"
      do
        read (*,*,iostat=ReadStatus), QReplicate
        if (ReadStatus.LE.0.AND.QReplicate.GE.0.AND.QReplicate.LE.2) exit
      end do
    end if

    if (QReplicate.EQ.0) then
      print*, "  > Incorrect number of files. Try again."
      deallocate (MidFile, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: GetParallelFiles: Dealloc #####"
    end if
  end if
  if (NFile.EQ.NMidFile.OR.QReplicate.GT.0) exit
end do

end subroutine GetParallelFiles

!*******************************************************************************
! get input file vector

subroutine GetInnFile

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InnFile)
NFile  = size (InnFile, 1)

end subroutine GetInnFile

!*******************************************************************************
! get output file vector

subroutine GetOutFile

print*, "  > Enter the substring in the OLD A files to alter in the NEW files: "
do
	read (*,*,iostat=ReadStatus), InnSub
	if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
	else if (InnSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
	end if
	if (ReadStatus.LE.0.AND.InnSub.NE."") exit
end do
InnSub=trim(adjustl(InnSub))

print "(3a)", "   > Enter the substring in the NEW files to replace '", trim(InnSub), "':"
do
	read (*,*,iostat=ReadStatus), OutSub
	if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
	else if (OutSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
	end if
	if (ReadStatus.LE.0.AND.OutSub.NE."") exit
end do

InnSubLen = len(trim(InnSub))
  
allocate (OutFile (NFile), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GloBulk: OutFile: Allocation failure #####"
OutFile = ""

do XFile = 1, NFile
      GivenFile = InnFile(XFile)
      InnSubBeg  = index(GivenFile,InnSub(1:InnSubLen))
      InnFileLen = len(adjustl(trim(GivenFile)))
      
      OutFile(XFile) = GivenFile(1:(InnSubBeg-1)) // trim(adjustl(OutSub)) // &
      		GivenFile((InnSubBeg+InnSubLen):InnFileLen)
end do

print*, "  > Enter the decimal places to save (1-3):"
do
      read (*,*,iostat=ReadStatus), QDecPlace
      if (ReadStatus.LE.0.AND.QDecPlace.GE.1.AND.QDecPlace.LE.3) exit
end do

end subroutine GetOutFile

!*******************************************************************************
! get 'mid' file index, depending on replication

subroutine GetMidFileIndex

if      (QReplicate.EQ.0) then		! work out 'mid' file index
    XMidFile = XFile
else if (QReplicate.EQ.1) then
    XMidFile = nint(mod(real(XFile),real(NMidFile)))
    if (XMidFile.EQ.0) XMidFile = NMidFile
else if (QReplicate.EQ.2) then
    if (XFile.EQ.1) XMidFile=0
    if (mod(real(XFile),(real(NFile)/real(NMidFile))).EQ.1) XMidFile=XMidFile+1
end if

GivenFile = OutFile(XFile)
LastSlash = index(GivenFile,"/",.TRUE.)
NameOnly  = trim(GivenFile((LastSlash+1):80))
  
end subroutine GetMidFileIndex

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

subroutine GetNameOnly

GivenFile = OutFile(XFile)
LastSlash = index(GivenFile,"/",.TRUE.)
NameOnly  = trim(GivenFile((LastSlash+1):80))
  
end subroutine GetNameOnly

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

end program AnnBulk