! anomproducts.f90
! f90 main program written on 20.12.99 by Tim Mitchell
! last modification on 02.10.00
! now fireproofed against PC crashes
! f90 -o anomproducts initialmod.f90 arraymod.f90 loadmod.f90 runselectmod.f90 extractmod.f90 
! transformmod.f90 savemod.f90 anomproducts.f90

program AnomProducts

use InitialMod
use LoadMod
use RunSelectMod
use ExtractMod	
use TransformMod
use SaveMod

implicit none

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

real, pointer, dimension (:)	:: WorkGloAnaSlice, TempIn, TempOut, TempOutA, TempOutB
real, pointer, dimension (:)	:: WorkBases, WorkMultiplier
real, pointer, dimension (:,:)	:: WorkLinAnaSeries, WorkGotDecade, WorkGotFull

integer, pointer, dimension (:) 	:: WorkADYear, WorkMapRawReg, WorkRegSizes, WorkAnaSizes
integer, pointer, dimension (:) 	:: WorkDecYearN, WorkDecGetYear0, WorkDecGetYear1
integer, pointer, dimension (:) 	:: WorkDecStyleThis, WorkDecStyleNext
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw, WorkMapIDLReg
integer, pointer, dimension (:,:) 	:: WorkDecBlockKey, WorkGotMem

integer, allocatable, dimension (:,:) 	:: DecValid

character (len=20), pointer, dimension (:) 	:: WorkRegNames
character (len=80), pointer, dimension (:) 	:: WorkDecPathThis, WorkDecPathNext

real, parameter :: MissVal = -999.0

integer :: WorkMonth0,WorkMonth1,WorkMonthN,WorkYearN,WorkDecN
integer :: WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkRegN
integer :: AllocStat, ReadStatus
integer :: XDec, XReg, XYear, XMem, XLong
integer :: Year0, Year1, ADYear0, ADYear1
integer :: DegChosen, SaveChosen, AnomChosen, QDeTrend
integer :: MemN, Unitary
integer :: DataMissN, RegValidN, YearValidN
integer :: QSaveUnsmoo, QSaveSmoo, QSaveRegMeans, QMaskPoles

character (len=10) :: WorkGridTitle
character (len=40) :: WorkRegTitle
character (len=80) :: WorkGridFilePath, WorkDecTitle, SaveTitle, GivenFile
character (len=80) :: WorkDecPathA, WorkDecPathB, Blank
character (len=80) :: SmooTitle, UnsmooTitle, SmooFile, UnsmooFile

!*******************************************************************************
! preliminaries

Blank = ""

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

print*
print*, "  > ***** AnomProducts *****"
print*, "  > Calculates anomalies using .glo as a base"
print*

call GridSelect   (WorkGrid,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call RegSelect    (WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkMapIDLReg,WorkRegSizes,WorkRegNames,&
		   WorkRegTitle,WorkRegN)

Unitary = 0				! finds out whether (=0) or not (=1) all regions are a single box
XReg    = 0
do
  XReg = XReg + 1
  if (WorkRegSizes(XReg).GT.1) Unitary = 1
  if (Unitary.GT.0) exit
  if (XReg.EQ.WorkRegN) exit
end do

QMaskPoles = 1				! decides whether or not to mask poles
if (WorkGrid.EQ.3.AND.Unitary.EQ.0) then
  print*, "  > Mask poles ? (1=no,2=yes)"
  do
	read (*,*,iostat=ReadStatus), QMaskPoles
	if (ReadStatus.LE.0 .AND. QMaskPoles.GE.1 .AND. QMaskPoles.LE.2) exit
  end do
end if

call PeriodSelect (WorkYearN,WorkDecN,WorkADYear)
call SeasonSelect (WorkMonth0,WorkMonth1,WorkMonthN)
call RawSelect    (WorkGrid,WorkLongN,WorkLatN,WorkMapIDLReg,WorkMapIDLRaw,WorkMapRawReg)

print*, "  > Detrend time series ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QDeTrend
	if (ReadStatus.LE.0 .AND. QDeTrend.GE.1 .AND. QDeTrend.LE.2) exit
end do

if (QDeTrend.EQ.2) then
	print*, "  > Select the .glo file with the year-multiplier to detrend."
	call LoadGlo (WorkLongN,WorkLatN,WorkRegN,WorkMapIDLReg,WorkMultiplier)
end if

print*, "  > Construct abs (=1), percent (=2) anomalies, or do not anomalise (=3) ? "
do
	read (*,*,iostat=ReadStatus), AnomChosen
	if (ReadStatus.LE.0 .AND. AnomChosen.GE.1 .AND. AnomChosen.LE.3) exit
end do

if (AnomChosen.NE.3) then
  print*, "  > Select the .glo file with the region bases."
  call LoadGlo      (WorkLongN,WorkLatN,WorkRegN,WorkMapIDLReg,WorkBases)
end if

allocate ( WorkAnaSizes    (WorkRegN), &
	   WorkGloAnaSlice (WorkRegN), &
	   WorkLinAnaSeries(WorkRegN, WorkYearN), &
	   WorkGotDecade   (WorkRegN, 10), &
	   WorkGotFull     (WorkRegN, WorkYearN), &
	   WorkGotMem      (WorkRegN, WorkYearN), &
	   DecValid        (WorkRegN, 10), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

print*, "  > Enter the number of ensemble members over which to average: "
do
	read (*,*,iostat=ReadStatus), MemN
	if (ReadStatus.LE.0 .AND. MemN.GE.1) exit
end do

WorkGotFull = 0.0
WorkGotMem  = 0

print "(a31,3i8)", "   > Members, regions, decades: ", MemN, WorkRegN, WorkDecN

!*******************************************************************************
! save options

print*, "  > Save unsmoothed region t-s to .tim ? (1=no,2=yes,fire-proof)"
do
	read (*,*,iostat=ReadStatus), QSaveUnsmoo
	if (ReadStatus.LE.0 .AND. QSaveUnsmoo.GE.1 .AND. QSaveUnsmoo.LE.2) exit
end do

if (QSaveUnsmoo.EQ.2) then
  print*, "  > Enter the .tim file title: "
  do
	read (*,*,iostat=ReadStatus), UnsmooTitle
	if (ReadStatus.LE.0.AND.UnsmooTitle.NE."") exit
  end do
  
  print*, "  > Enter the filepath of the .tim file to save: "
  do
	do
		read (*,*,iostat=ReadStatus), GivenFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=GivenFile, name=UnsmooFile)
	open (1, file=UnsmooFile, status="new", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
  end do
end if  

print*, "  > Save smoothed (Gauss,30y) region t-s to .tim ? (1=no,2=yes,fire-proof)"
do
	read (*,*,iostat=ReadStatus), QSaveSmoo
	if (ReadStatus.LE.0 .AND. QSaveSmoo.GE.1 .AND. QSaveSmoo.LE.2) exit
end do

if (QSaveSmoo.EQ.2) then
  print*, "  > Enter the .tim file title: "
  do
	read (*,*,iostat=ReadStatus), SmooTitle
	if (ReadStatus.LE.0.AND.SmooTitle.NE."") exit
  end do
  
  print*, "  > Enter the filepath of the .tim file to save: "
  do
	do
		read (*,*,iostat=ReadStatus), GivenFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=GivenFile, name=SmooFile)
	open (1, file=SmooFile, status="new", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
  end do
end if  

print*, "  > Save region means to .glo ? (1=no,2=yes,not fire-proof)"
do
	read (*,*,iostat=ReadStatus), QSaveRegMeans
	if (ReadStatus.LE.0 .AND. QSaveRegMeans.GE.1 .AND. QSaveRegMeans.LE.2) exit
end do

!*******************************************************************************
! extract data into WorkGotFull

do XMem = 1, MemN 

  if (MemN.GT.1) print*, "  > Select ensemble member: ", XMem
  
  call RunSelect    (WorkGrid,WorkMonth0,WorkMonth1,WorkYearN,WorkDecN,WorkADYear,&
		     WorkDecTitle,WorkDecStyleThis,WorkDecStyleNext,&
		     WorkDecPathThis,WorkDecPathNext,WorkDecYearN,WorkDecGetYear0,WorkDecGetYear1)

  print*, "  > Non-missing data loaded thus far, by decade: "
  print*, "  > Year AllData Regions   Years"

  do XDec = 1, WorkDecN

    DataMissN  = 0
    YearValidN = 0
    RegValidN  = 0
    DecValid   = 0
    
    Year0 = (XDec-1)*10 + 1
    Year1 = Year0 + 9
  
    WorkGotDecade = 0.0

    call BlockKey    (WorkDecYearN(XDec),WorkMonthN,WorkDecGetYear0(XDec),WorkDecGetYear1(XDec),&
  		      WorkMonth0,WorkMonth1,WorkDecPathThis(XDec),WorkDecPathNext(XDec),&
  		      WorkDecStyleThis(XDec),WorkDecStyleNext(XDec),WorkDecBlockKey)
    
    call ExtractFile (WorkLongN,WorkLatN,WorkDataN,WorkDecPathThis(XDec),WorkDecPathNext(XDec),&
    		      WorkDecStyleThis(XDec),WorkDecStyleNext(XDec),WorkRegN,WorkMonthN,&
  		      WorkDecYearN(XDec),WorkDecGetYear0(XDec),&
  		      WorkMapRawReg,WorkRegSizes,WorkDecBlockKey,WorkGotDecade)
  
    do XReg = 1, WorkRegN
     do XYear = 1, 10
      if (WorkGotDecade (XReg,XYear).NE.MissVal) then
       	WorkGotFull (XReg,(Year0+XYear-1)) = WorkGotFull (XReg,(Year0+XYear-1)) + &
       						WorkGotDecade (XReg,XYear)
       	WorkGotMem  (XReg,(Year0+XYear-1)) = WorkGotMem  (XReg,(Year0+XYear-1)) + 1
        DecValid    (XReg,XYear)           = 1
      else
        DataMissN = DataMissN + 1
      end if
     end do
    end do
    
    do XReg = 1, WorkRegN			! finds number of valid regions
      XYear = 0
      do
        XYear = XYear + 1
        if (DecValid(XReg,XYear).EQ.1) RegValidN = RegValidN + 1
        if (DecValid(XReg,XYear).EQ.1) exit
        if (XYear.EQ.10) exit
      end do      
    end do
    
    do XYear = 1, 10				! finds number of valid years
      XReg = 0
      do
        XReg = XReg + 1
        if (DecValid(XReg,XYear).EQ.1) YearValidN = YearValidN + 1
        if (DecValid(XReg,XYear).EQ.1) exit
        if (XReg.EQ.WorkRegN) exit
      end do      
    end do
    
    XYear = Year0 - 1				! finds first valid year in decade to report to screen
    do
      XYear = XYear + 1
      if (WorkADYear(XYear).NE.MissVal) exit
      if (XYear.EQ.Year1) exit
    end do    
    
    print "(a1,4i8)", " ", WorkADYear(XYear), ((10*WorkRegN)-DataMissN), RegValidN, YearValidN
  end do

end do

deallocate (DecValid, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: AnomProducts: deallocation failure #####"

!*******************************************************************************
! average out members into mean

do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (WorkGotMem(XReg,XYear).NE.0) then
       WorkGotFull(XReg,XYear) = WorkGotFull(XReg,XYear) / real(WorkGotMem(XReg,XYear))
    else
       WorkGotFull(XReg,XYear) = MissVal
    end if
  end do
end do

!*******************************************************************************
! detrend (by removing multiplier * yearAD / 100)

if (QDeTrend.EQ.2) then 

 do XYear = 1, WorkYearN
  if (WorkADYear(XYear).NE.MissVal) then
   do XReg = 1, WorkRegN
    if (WorkGotFull(XReg,XYear).NE.MissVal.AND.WorkMultiplier(XReg).NE.MissVal) then
       WorkGotFull(XReg,XYear) = WorkGotFull(XReg,XYear) - (WorkMultiplier(XReg)*real(WorkADYear(XYear))/100.0)
    else
       WorkGotFull(XReg,XYear) = MissVal
    end if
   end do
  else
   do XReg = 1, WorkRegN
       WorkGotFull(XReg,XYear) = MissVal
   end do  
  end if
 end do

end if

!*******************************************************************************
! anomalise

if (AnomChosen.EQ.1) then 

 do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (WorkGotFull(XReg,XYear).NE.MissVal.AND.WorkBases(XReg).NE.MissVal) then
       WorkGotFull(XReg,XYear) = WorkGotFull(XReg,XYear) - WorkBases(XReg)
    else
       WorkGotFull(XReg,XYear) = MissVal
    end if
  end do
 end do

else if (AnomChosen.EQ.2) then 

 do XReg = 1, WorkRegN
  do XYear = 1, WorkYearN
    if (WorkGotFull(XReg,XYear).NE.MissVal.AND.WorkBases(XReg).NE.MissVal) then
       WorkGotFull(XReg,XYear) = (WorkGotFull(XReg,XYear) - WorkBases(XReg)) &
       				 * 100.0 / WorkBases(XReg)
    else
       WorkGotFull(XReg,XYear) = MissVal
    end if
  end do
 end do

end if

!*******************************************************************************
! mask poles

if (QMaskPoles.EQ.2) then
  do XLong = 1, WorkLongN
    do XYear = 1, WorkYearN
    	WorkGotFull (WorkMapIDLReg(XLong,1),       XYear) = MissVal
    	WorkGotFull (WorkMapIDLReg(XLong,WorkLatN),XYear) = MissVal
    end do
  end do
end if

!*******************************************************************************
! save unsmoothed time series to .tim

if (QSaveUnsmoo.EQ.2) then
  call SaveTim (WorkRegN,WorkYearN,UnsmooTitle,UnsmooFile,WorkRegNames,WorkADYear,WorkGotFull)
end if

!*******************************************************************************
! generate and save smoothed time series to .tim

if (QSaveSmoo.EQ.2) then

  WorkLinAnaSeries = 0.0

  allocate (TempIn   (WorkYearN), &
  	  TempOutA (WorkYearN), &
  	  TempOutB (WorkYearN), stat=AllocStat)

  do XReg = 1, WorkRegN
    if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  
    TempIn (1:WorkYearN) = WorkGotFull (XReg,1:WorkYearN)
    TempOutA = MissVal
    TempOutB = MissVal
  
    call GaussSmooth (WorkYearN,30,1,TempIn,TempOutA,TempOutB)
  
    WorkLinAnaSeries (XReg,1:WorkYearN) = TempOutA (1:WorkYearN)
  end do

  deallocate (TempIn, TempOutA, TempOutB)

  call SaveTim (WorkRegN,WorkYearN,SmooTitle,SmooFile,WorkRegNames,WorkADYear,WorkLinAnaSeries)

end if

!*******************************************************************************
! generate and save means to .glo

if (QSaveRegMeans.EQ.2) then
 
 ADYear0 = -1
 ADYear1 = -1
 
 do
  
  if (ADYear1.GT.ADYear0) then
   
   Year0 = -1			! finds index year matching first AD year of slice
   XYear = 0
   do
     XYear = XYear + 1
     if (WorkADYear(XYear).EQ.ADYear0) Year0 = XYear
     if (Year0.NE.-1) exit
     if (XYear.EQ.WorkYearN) exit
   end do
   
   Year1 = -1			! finds index year matching last AD year of slice
   XYear = 0
   do
     XYear = XYear + 1
     if (WorkADYear(XYear).EQ.ADYear1) Year1 = XYear
     if (Year1.NE.-1) exit
     if (XYear.EQ.WorkYearN) exit
   end do
   
   if (Year0.NE.-1.AND.Year1.NE.-1) then	! if first and last index years are valid
   
    WorkGloAnaSlice  = 0.0
    WorkAnaSizes     = 0

    do XReg = 1, WorkRegN
     do XYear = Year0, Year1			! average over period selected
      if (WorkGotFull(XReg,XYear).NE.MissVal) then
        WorkGloAnaSlice (XReg) = WorkGloAnaSlice (XReg) + WorkGotFull(XReg,XYear)
        WorkAnaSizes (XReg)    = WorkAnaSizes (XReg) + 1
      end if
     end do
  
     if (WorkAnaSizes(XReg).GT.0) then
      WorkGloAnaSlice (XReg) = WorkGloAnaSlice (XReg) / real (WorkAnaSizes (XReg))
     else
      WorkGloAnaSlice (XReg) = MissVal
     end if
     
     if (WorkRegN.LE.5) print "(i4,a1,a20,f12.4)", XReg, " ", WorkRegNames(XReg), WorkGloAnaSlice (XReg)
    end do

    call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,&
  		WorkGloAnaSlice,WorkMapIDLReg)
  
   end if
  
  end if
  
  do
        print*, "  > Enter first, last years AD of slice to save (99,99=end): "
	read (*,*,iostat=ReadStatus), ADYear0, ADYear1
	if (ReadStatus.LE.0) exit
  end do
  
  if (ADYear0.EQ.99.AND.ADYear1.EQ.99) exit
 
 end do
 
end if

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

print*

deallocate (WorkAnaSizes,WorkGloAnaSlice,WorkLinAnaSeries,WorkGotDecade,WorkGotFull,WorkGotMem)
deallocate (WorkRegNames)

close (99)

end program AnomProducts