! perbulk.f90
! program to handle .per files in bulk
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../obs/perbulk time.f90 filenames.f90 perfiles.f90 annfiles.f90 
!	plainperm.f90 regress.f90 
! 	cetgeneral.f90 ./../obs/perbulk.f90 2> /tyn1/tim/scratch/stderr.txt
! last modified on 13.03.02

program PerBulk

use Time
use FileNames
use PerFiles
use AnnFiles
use Regress
use CETGeneral

implicit none

real, pointer, dimension (:,:)			:: InMon, InSea, OutDec, OutAll, LSRAye, LSRBee
real, pointer, dimension (:)			:: InAnn, LSRExe, LSRWye, VecA,VecB,VecC

integer, pointer, dimension (:) 		:: YearAD

character (len=80), pointer, dimension (:) 	:: InFileA, InFileB, OutFileA, FileNameOnly
character (len= 9), pointer, dimension(:)	:: ColTitles			! can be 'X', no blanks

real, parameter :: MissVal = -999.0

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

real :: OpTot, OpEn, OpMiss, BaseTot, BaseEn
real :: LSRPea

integer :: ReadStatus, AllocStat
integer :: MenuChoice, FileVariCode
integer :: FileN, FileYearN, YearN
integer :: XFile, XFileYear, XYear, XMonth, XSeason, XDecade, XAnnual, XDecadeYear, XBegYear, XBegDecade
integer :: BegYearAD, Silent, SmooPer
integer :: InSubLen,InSubBeg,InFileALen,LastSlash,OutSubBeg,StatType,NameEnd
integer :: QVariTime,QLSRForm,QMissCalc

character (len=80) :: GivenFile, NameOnly
character (len=80) :: InSub, OutSub

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

call Intro
call Menu
call Final

contains

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

subroutine Intro

open (99,file="/tyn1/tim/scratch/log-perbulk.dat",status="replace",action="write")
print*
print*, "  > ***** PerBulk.f90 : handles .per files in bulk *****"
print*

end subroutine Intro

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

subroutine Menu

MenuChoice = 0
do
  print*, "  > MENU: make your choice (0=list):"
  do
	read (*,*,iostat=ReadStatus), MenuChoice
	if (ReadStatus.LE.0) exit
  end do
  
  if      (MenuChoice.EQ. 0) then
    print*, "  > 10. save InAnn OutDec means"
    print*, "  > 20. regression: y=at or y=at+b"
    print*, "  > 30. smooth with Gaussian filter"
    print*, "  > 99. end"    
  else if (MenuChoice.EQ.10) then
    call SaveAnnDecMeans
  else if (MenuChoice.EQ.20) then
    call RegressionYATB
  else if (MenuChoice.EQ.30) then
    call SmoothGaussFilter
  else if (MenuChoice.EQ.99) then
    print*
  else
    print*, "  > Pardon ?"
  end if
  
  if (MenuChoice.EQ.99) exit
end do

end subroutine Menu

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

subroutine SaveAnnDecMeans

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

call GetOutFileA
do XFile = 1, FileN
  GivenFile = OutFileA(XFile)
  OutSubBeg = index(GivenFile,'.per')
  if (OutSubBeg.GT.0) OutFileA(XFile) = GivenFile(1:(OutSubBeg-1)) // ".ann"
end do

print*, "  > Select the year AD to begin to calc: "
do
	read (*,*,iostat=ReadStatus), BegYearAD
	if (ReadStatus.LE.0) exit
end do

do XFile = 1, FileN
  GivenFile = InFileA(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print*
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadPer (InFileA(XFile), FileVariCode, YearAD, InMon, InSea, InAnn)
  deallocate (InMon, InSea)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: MonSea: Deallocation failure #####"

  FileYearN = size (YearAD,1)
  
  allocate (OutDec (FileYearN,5), &
	    ColTitles (5),         stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: OutDec: Allocation failure #####"
  OutDec = MissVal
  ColTitles = (/'  ann tot',' dec anom',' dec mean',' no years','  missing'/)

  XFileYear = 0 ; XBegYear = 0 ; XDecadeYear = 0 ; XDecade = 0
  OpTot = 0.0 ; OpEn = 0.0 ; OpMiss = 0.0 ; BaseTot = 0.0 ; BaseEn = 0.0
  do
    XFileYear = XFileYear + 1
    if (InAnn(XFileYear).NE.MissVal) OutDec(XFileYear,1) = InAnn(XFileYear)
    
    if (YearAD(XFileYear).EQ.BegYearAD) XBegYear = XFileYear
        
    if (XBegYear.GT.0) then
      XDecadeYear = XDecadeYear + 1      
      if (InAnn(XFileYear).NE.MissVal) then
        OpTot = OpTot + InAnn(XFileYear)
        OpEn  = OpEn  + 1
      else
        OpMiss = OpMiss + 1
        OpEn   = OpEn   + 1
      end if
    end if
    
    if (XDecadeYear.EQ.10.OR.XFileYear.EQ.FileYearN) then
      XDecade = XDecade + 1
      XBegDecade = XBegYear + ((XDecade-1)*10)
      
      if (OpEn.GT.0) then
        OutDec (XBegDecade,3) = OpTot / OpEn
        OutDec (XBegDecade,4) = OpEn
        OutDec (XBegDecade,5) = OpMiss
      end if
      
      XDecadeYear = 0 ; OpTot = 0.0 ; OpEn = 0.0 ; OpMiss = 0.0
    end if
    
    if (YearAD(XFileYear).GE.1961.AND.YearAD(XFileYear).LE.1990) then
      if (InAnn(XFileYear).NE.MissVal) then
        BaseTot = BaseTot + InAnn(XFileYear)
        BaseEn  = BaseEn  + 1
      end if
    end if
    
    if (XFileYear.EQ.FileYearN) exit
  end do
  
  if (BaseEn.GE.10) then
    do XFileYear = 1, FileYearN
      if (OutDec(XFileYear,3).NE.MissVal) OutDec(XFileYear,2) = OutDec(XFileYear,3) - (BaseTot/BaseEn)
    end do
  end if
  
  call SaveANN (OutFileA(XFile), YearAD, ColTitles, OutDec, DecPlaceN=1)

  deallocate (InAnn, YearAD, OutDec, ColTitles)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: AnnYAD: Deallocation failure #####"
end do

deallocate (InFileA,OutFileA)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: DecInOut: Deallocation failure #####"

end subroutine SaveAnnDecMeans

!*******************************************************************************
! regress y=at+b: inB = (outA * inA) + outB

subroutine RegressionYATB

print*, "  > Regression form: y=at (=1) or y=at+b (=2) ?" 
do
	read (*,*,iostat=ReadStatus), QLSRForm
	if (ReadStatus.LE.0.AND.QLSRForm.GE.1.AND.QLSRForm.LE.2) exit
end do
  
print*, "  > Select the 'y' files. There must be a common substring."
call GetBatch (Blank,InFileA)
FileN = size (InFileA, 1)

allocate (LSRAye (FileN,17), &
  	  LSRBee (FileN,17), &
  	  FileNameOnly (FileN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: RegressionYATB: Allocation failure: A,B #####"
LSRAye = MissVal ; LSRBee = MissVal

print*, "  > Operating in bulk..."
do XFile = 1, FileN
  GivenFile = InFileA(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  FileNameOnly(XFile)  = adjustl(trim(GivenFile((LastSlash+1):80)))
  
  call LoadPer (InFileA(XFile), FileVariCode, YearAD, InMon, InSea, InAnn, Silent=1)  
  FileYearN = size (YearAD,1)
  
  allocate (LSRExe (FileYearN), &
  	    LSRWye (FileYearN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: RegressionYATB: Allocation failure: X,Y #####"
  
  do XYear = 1, FileYearN
      LSRExe (XYear) = YearAD (XYear)
  end do

  do XMonth = 1, 12
    do XYear = 1, FileYearN
      LSRWye (XYear) = InMon (XYear,XMonth)
    end do
    
    if (QLSRForm.EQ.1) then
      call LinearLSRZeroVec (LSRWye,LSRExe,LSRAye(XFile,XMonth))
    else
      call LinearLSRVec (LSRWye,LSRExe,LSRBee(XFile,XMonth),LSRAye(XFile,XMonth),LSRPea)
    end if
  end do
  
  do XSeason = 1, 4
    do XYear = 1, FileYearN
      LSRWye (XYear) = InSea (XYear,XSeason)
    end do
    
    if (QLSRForm.EQ.1) then
      call LinearLSRZeroVec (LSRWye,LSRExe,LSRAye(XFile,(XSeason+12)))
    else
      call LinearLSRVec (LSRWye,LSRExe,LSRBee(XFile,(XSeason+12)),LSRAye(XFile,(XSeason+12)),LSRPea)
    end if
  end do
  
  do XYear = 1, FileYearN
      LSRWye (XYear) = InAnn (XYear)
  end do
    
  if (QLSRForm.EQ.1) then
      call LinearLSRZeroVec (LSRWye,LSRExe,LSRAye(XFile,17))
  else
      call LinearLSRVec (LSRWye,LSRExe,LSRBee(XFile,17),LSRAye(XFile,17),LSRPea)
  end if
  
  deallocate (LSRExe,LSRWye, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: RegressionYATB: Deallocation failure: X,Y #####"  

								!###########################
  write (99,"(a20,17e12.5)"), trim(FileNameOnly(XFile)), (LSRAye(XFile,XSeason),XSeason=1,17)
end do

write (99,*), "GRADIENTS  (12*mon,4*sea,1*ann)"
do XFile = 1, FileN
  write (99,"(a20,17e12.5)"), trim(FileNameOnly(XFile)), (LSRAye(XFile,XSeason),XSeason=1,17)
end do

if (QLSRForm.EQ.2) then
 write (99,*), "INTERCEPTS (12*mon,4*sea,1*ann)"
 do XFile = 1, FileN
  write (99,"(a20,17e12.5)"), trim(FileNameOnly(XFile)), (LSRBee(XFile,XSeason),XSeason=1,17)
 end do
end if

print*, "  > We have dumped the results to log file."

end subroutine RegressionYATB

!*******************************************************************************
! smooth with Gaussian filter

subroutine SmoothGaussFilter

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

call GetOutFileA

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

print*, "  > Keep missing seas/ann (=1), or calc from interpol monthly (=2) ?"
do
	read (*,*,iostat=ReadStatus), QMissCalc
	if (ReadStatus.LE.0.AND.QMissCalc.GE.1.AND.QMissCalc.LE.2) exit
end do

print*, "  > Enter the statistic (-1=min,0=mean,1=max,2=sum): "
do
	read (*,*,iostat=ReadStatus), StatType
	if (ReadStatus.LE.0.AND.StatType.GE.-1.AND.StatType.LE.2) exit
end do

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

do XFile = 1, FileN
 call LoadPer (InFileA(XFile), FileVariCode, YearAD, InMon, InSea, InAnn, Silent=1)  
 
 GivenFile = InFileA(XFile)
 LastSlash = index(GivenFile,"/",.TRUE.)
 NameEnd = index(GivenFile,".per",.TRUE.) ; NameEnd=NameEnd+3
 NameOnly  = trim(GivenFile((LastSlash+1):NameEnd))
  
 YearN = size (YearAD,1)
  
 if (YearN.GT.SmooPer) then
  print "(2a)", "   > Smoothing ", trim(NameOnly)
  
  allocate (VecA(YearN), &
  	    VecB(YearN), &
  	    VecC(YearN), &
  	    OutAll(YearN,17), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SmoothGaussFilter: Allocation failure: Vec* #####"
 
  do XMonth = 1, 12
    VecA=MissVal ; VecB=MissVal ; VecC=MissVal
    
    do XYear = 1, YearN
      VecA(XYear) = InMon(XYear,XMonth)
    end do
    
    call GaussSmooth (YearN,SmooPer,1,VecA,VecB,VecC,MissThresh=20.0)
    
    do XYear = 1, YearN
      OutAll(XYear,XMonth) = VecB(XYear)
      if (InMon(XYear,XMonth).EQ.MissVal) InMon(XYear,XMonth) = VecB(XYear)
    end do
  end do
  
  if (QMissCalc.EQ.2) then
    InSea=MissVal ; InAnn=MissVal
    
    if (StatType.EQ.-1) call FillSeaAnnMin  (YearAD,InMon,InSea,InAnn)
    if (StatType.EQ. 0) call FillSeaAnnMean (YearAD,InMon,InSea,InAnn)
    if (StatType.EQ. 1) call FillSeaAnnMax  (YearAD,InMon,InSea,InAnn)
    if (StatType.EQ. 2) call FillSeaAnnSum  (YearAD,InMon,InSea,InAnn)  
  end if
  
  do XSeason = 1, 4
    VecA=MissVal ; VecB=MissVal ; VecC=MissVal
    
    do XYear = 1, YearN
      VecA(XYear) = InSea(XYear,XSeason)
    end do
    
    call GaussSmooth (YearN,SmooPer,1,VecA,VecB,VecC,MissThresh=20.0)
    
    do XYear = 1, YearN
      OutAll(XYear,(XSeason+12)) = VecB(XYear)
    end do
  end do
  
  do XAnnual = 1, 1
    VecA=MissVal ; VecB=MissVal ; VecC=MissVal
    
    do XYear = 1, YearN
      VecA(XYear) = InAnn(XYear)
    end do
    
    call GaussSmooth (YearN,SmooPer,1,VecA,VecB,VecC,MissThresh=20.0)
    
    do XYear = 1, YearN
      OutAll(XYear,17) = VecB(XYear)
    end do
  end do
  
  call SavePER (OutFileA(XFile), YearAD, StatType, AllData=OutAll, NoResponse=1)
  
  deallocate (VecA,VecB,VecC,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SmoothGaussFilter: Deallocation failure: Vec* #####"
 else
  print "(3a)", "   > Rejecting ", trim(NameOnly), " due to short record."
 end if
  
 deallocate (InMon,InSea,InAnn,YearAD,stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: SmoothGaussFilter: Deallocation failure: In* #####"
end do

end subroutine SmoothGaussFilter

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

subroutine GetOutFileA

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

print "(3a)", "   > Enter the substring in the NEW files to replace '", trim(InSub), "':"
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

InSubLen = len(trim(InSub))
  
allocate (OutFileA (FileN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PerBulk: OutFileA: Allocation failure #####"
OutFileA = ""

do XFile = 1, FileN
      GivenFile = InFileA(XFile)
      InSubBeg  = index(GivenFile,trim(InSub))
      InFileALen = len(trim(GivenFile))
      
      OutFileA(XFile) = GivenFile(1:(InSubBeg-1)) // trim(OutSub) // &
      		GivenFile((InSubBeg+InSubLen):InFileALen)
end do

end subroutine GetOutFileA

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

subroutine Final

close (99)

end subroutine Final

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

end program PerBulk