! interpol.f90
! module to hold routines for spatial interpolation
! contains: CrudeInterpol,FlatSurf,CrudeSmoo,Dupli90s
! requires higher up in call: SortMod, FileNames, GrimFiles, Grid, Regress 

module Interpol

implicit none

contains

!*******************************************************************************
! manages crude interpolation from one grid to another, minimising array sizes while loading all at once
! approach:
! 1. initially loads all data and defines grids
! 2. iterates through each time step and transforms spatial field into high-res spatial field
! 3. saves high-res data set as new grim
! steps in iteration:
! a. places all raw data from RAW vector into larger vector (same LOw res) padded with missing boxes
! b. by growing the valid domain we give values to all the LOw res missing boxes (using NearNeigh)
! c. we smooth out the LOw res data into HIgh res data (using CrudeSmoo)
! d. we SHRINK the bounds of the HIgh res data to match the target (InterPOLated) data
! e. we reduce the domain of the SHRINK data into that of the target (InterPOLated) data

subroutine CrudeInterpol (LoadFile,SaveFile,ExeMulti,WyeMulti,Extent,IpolGrid,IpolBounds,IpolName)

use FileNames
use GrimFiles
use Grid

real, pointer, dimension (:,:,:)	:: BriefRawData, RawData,IpolData
real, pointer, dimension (:)		:: LoVec,IterVec,HiVec,ShrinkVec

integer, pointer, dimension (:,:)	:: RawGrid,IpolGrid
integer, pointer, dimension (:,:)	:: LoGrid,HiGrid,LoToHi,ShrinkGrid,HiToShrink
integer, pointer, dimension (:)		:: BriefYearAD,YearAD

real, dimension (4) 			:: RawBounds,IpolBounds

real, intent(in)			:: ExeMulti, WyeMulti
integer, intent(in)			:: Extent

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: YearN,MonthN,XMonth,XYear
integer :: RawExeN, RawWyeN, RawBoxN, RawCellN, XRawExe, XRawWye
integer :: LoExeN, LoWyeN, LoBoxN, XLoExe, XLoWye
integer :: IpolExeN,IpolWyeN,IpolBoxN,XIpolExe,XIpolWye
integer :: ShrinkBoxN,HiExeN,HiWyeN,HiCellN,XIter,CountMiss
integer :: ValidShrink

character (len=80) :: RawInfo,IpolInfo,IpolName
character (len=80) :: LoadFile,SaveFile
character (len= 4) :: RawSuffix,IpolSuffix

!***************************************
! get initial data and all grids
								! load Raw data+grid
call LoadGrim (RawData,RawGrid,YearAD,RawBounds,RawInfo,LoadFile,"    ",RawSuffix)

!call LoadGrim (BriefRawData,RawGrid,BriefYearAD,RawBounds,RawInfo,LoadFile,"    ",RawSuffix)
!print*, "  > Duplicating into 1996-2000..."
!call Dupli90s (BriefRawData,RawData,BriefYearAD,YearAD)
!deallocate (BriefRawData,BriefYearAD,stat=AllocStat)
!if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Allocation failure: Brief* #####"

print*, "  > Initialising..."

IpolInfo = trim(RawInfo) // ' ->' // trim(IpolName) 

YearN   = size(RawData,1) ; MonthN  = size(RawData,2) ; RawBoxN = size(RawData,3)
RawExeN = size(RawGrid,1) ; RawWyeN = size(RawGrid,2) ; RawCellN = RawExeN * RawWyeN
LoExeN = RawExeN ; LoWyeN = RawWyeN ; LoBoxN = RawCellN

IpolExeN = size (IpolGrid,1) ; IpolWyeN = size (IpolGrid,2) ; IpolBoxN = maxval(IpolGrid)

allocate (LoGrid    (LoExeN,LoWyeN),     &
	  ShrinkGrid(IpolExeN,IpolWyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Allocation failure: LoGrid #####"

do XLoExe = 1, LoExeN						! get Lo grid
  do XLoWye = 1, LoWyeN
    LoGrid(XLoExe,XLoWye) = ((XLoExe-1)*LoWyeN) + XLoWye
  end do
end do

call MakeHiResGrid (ExeMulti,WyeMulti,LoGrid,HiGrid,LoToHi)	! get Hi and LoToHi grids
HiExeN = size (HiGrid,1) ; HiWyeN = size (HiGrid,2) ; HiCellN = HiExeN * HiWyeN

do XIpolExe = 1, IpolExeN					! get Shrink grid
  do XIpolWye = 1, IpolWyeN
    ShrinkGrid(XIpolExe,XIpolWye) = ((XIpolExe-1)*IpolWyeN) + XIpolWye
  end do 
end do
ShrinkBoxN = IpolExeN * IpolWyeN
    								! get HiToShrink grid
call DefineShrinking (RawBounds,IpolBounds,HiGrid,ShrinkGrid,HiToShrink,ValidShrink)
    
allocate (IpolData (YearN,MonthN,IpolBoxN), &
	  HiVec    (HiCellN),               &
	  ShrinkVec(ShrinkBoxN),            stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Allocation failure: ID,HV,SV #####"
IpolData = MissVal ; HiVec = MissVal ; ShrinkVec = MissVal
    
!***************************************
! iterate a month at a time

if (ValidShrink.EQ.0) then
 print*, "  > Iterating by time step..."

 do XYear = 1, YearN
  do XMonth = 1, MonthN
    allocate (LoVec(RawCellN),  &
    	      IterVec(RawCellN),stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Allocation failure: LV #####"
    LoVec = MissVal ; IterVec = MissVal
    
    do XRawExe = 1, RawExeN					! get all Lo vector info available from Raw
      do XRawWye = 1, RawWyeN
        if (RawGrid(XRawExe,XRawWye).NE.MissVal) &
        		LoVec(LoGrid(XRawExe,XRawWye)) = RawData(XYear,XMonth,RawGrid(XRawExe,XRawWye))
      end do
    end do
    
    do								! get rest of Lo vector
      if (maxval(LoVec).GT.MissVal) then			! if there are non-missing values present...
        if (minval(LoVec).EQ.MissVal) then			! and there are also missing values present
          call NearNeigh (1,LoGrid,LoVec,IterVec)      		! method is to iterate through
          LoVec = IterVec					! on each iteration we grow the edges of valid
        end if
      else
        print "(a,2i4)", "   > *** WARNING: month is missing: year,month = ", YearAD(XYear), XMonth
      end if
      
      if (maxval(LoVec).EQ.MissVal) exit			! unless the entire vector is missing	
      if (minval(LoVec).GT.MissVal) exit			! or until the entire vector is valid	
    end do
    
    deallocate (IterVec,stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Deallocation failure: IV #####"

    HiVec = MissVal
    call CrudeSmoo (Extent,HiGrid,LoToHi,LoVec,HiVec)		! get Hi vector by local smoothing
    
    nullify (LoVec)
    
    ShrinkVec = MissVal
    do XIpolExe = 1, IpolExeN					! shrink Hi vector into Shrink vector
      do XIpolWye = 1, IpolWyeN					
        if (ShrinkGrid(XIpolExe,XIpolWye).NE.MissVal) then
          if (HiToShrink(XIpolExe,XIpolWye).NE.MissVal) then
              ShrinkVec(ShrinkGrid(XIpolExe,XIpolWye)) = HiVec(HiToShrink(XIpolExe,XIpolWye))
          end if
        end if
      end do
    end do
    
    do XIpolExe = 1, IpolExeN					! redomain into IpolData
      do XIpolWye = 1, IpolWyeN
        if (IpolGrid(XIpolExe,XIpolWye).NE.MissVal) then
          if (ShrinkGrid(XIpolExe,XIpolWye).NE.MissVal) then
            IpolData(XYear,XMonth,IpolGrid(XIpolExe,XIpolWye)) = ShrinkVec(ShrinkGrid(XIpolExe,XIpolWye))
          end if
        end if
      end do
    end do
  end do
 end do
 
 deallocate (RawData,stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Deallocation failure: RawData #####"    

 print*, "  > Saving..."

 call SaveGrim (IpolData,IpolGrid,YearAD,IpolBounds,IpolInfo,SaveFile,RawSuffix,IpolSuffix)
 
 deallocate (LoGrid,HiGrid,LoToHi,ShrinkGrid,HiToShrink,IpolData,YearAD,stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Deallocation failure: final #####"
end if

deallocate (HiVec,    &
	    ShrinkVec,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CrudeInterpol: Deallocation failure: SV #####"    

end subroutine CrudeInterpol

!*******************************************************************************
! (ideally) takes 2D array with missing values and returns flat surface without missing values
! appears to produce the desired output in some circumstances, but not all
! the potential problem may lie in v. large data sets, sparse data sets, or high contrast data sets
! Data is supplied in the call, Surf is returned from the call

subroutine FlatSurf (Data,Surf)

use Regress

real, pointer, dimension (:,:) 	:: Data,Surf
real, pointer, dimension (:) 	:: DataVec,DataRef,DataFit,VecCept,VecGrad

real, parameter :: MissVal = -999.0

real :: VecTot,VecEn
real :: ExeCept,ExeGrad,WyeCept,WyeGrad
real :: Correl

integer :: AllocStat
integer :: Iter
integer :: ExeN, WyeN
integer :: XExe, XWye

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

ExeN = size(Data,1) ; WyeN = size(Data,2)
ExeCept = MissVal ; ExeGrad = MissVal ; WyeCept = MissVal ; WyeGrad = MissVal 

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

allocate (DataVec(WyeN), &			! calc equation for each set of Wye
	  DataRef(WyeN), &
	  VecCept(ExeN), &
	  VecGrad(ExeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FlatSurf: Allocation failure: Exe #####"

do XWye = 1, WyeN
  DataRef(XWye) = XWye
end do

do XExe = 1, ExeN  
  do XWye = 1, WyeN
    DataVec(XWye) = Data(XExe,XWye)    
  end do
  
  call LinearLSRVec (DataRef,DataVec,VecCept(XExe),VecGrad(XExe),Correl)
  
!  call RobustLSR (WyeN,DataRef,DataVec,DataFit,VecCept(XExe),VecGrad(XExe),Iter)
  
!  deallocate(DataFit,stat=AllocStat)
!  if (AllocStat.NE.0) print "(a,i8)", "   > ##### ERROR: FlatSurf: Deallocation failure: ", XExe
end do

VecTot = 0.0 ; VecEn = 0.0			! calc average gradient from all sets of Wye
do XExe = 1, ExeN
  if (VecGrad(XExe).NE.MissVal) then
    VecTot = VecTot + VecGrad(XExe)
    VecEn  = VecEn  + 1.0
  end if
end do
if (VecEn.GE.1) ExeGrad = VecTot / VecEn

VecTot = 0.0 ; VecEn = 0.0			! calc average intercept from all sets of Wye
do XExe = 1, ExeN
  if (VecCept(XExe).NE.MissVal) then
    VecTot = VecTot + VecCept(XExe)
    VecEn  = VecEn  + 1.0
  end if
end do
if (VecEn.GE.1) ExeCept = VecTot / VecEn

deallocate (DataVec,DataRef,VecCept,VecGrad, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FlatSurf: Deallocation failure: Exe #####"

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

allocate (DataVec(ExeN), &			! calc equation for each set of Exe
	  DataRef(ExeN), &
	  VecCept(WyeN), &
	  VecGrad(WyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FlatSurf: Allocation failure: Wye #####"

do XExe = 1, ExeN
  DataRef(XExe) = XExe
end do

do XWye = 1, WyeN  
  do XExe = 1, ExeN
    DataVec(XExe) = Data(XExe,XWye)    
  end do
  
  call LinearLSRVec (DataRef,DataVec,VecCept(XExe),VecGrad(XExe),Correl)
  
!  call RobustLSR (ExeN,DataRef,DataVec,DataFit,VecCept(XWye),VecGrad(XWye),Iter)
  
!  deallocate(DataFit,stat=AllocStat)
!  if (AllocStat.NE.0) print "(a,i8)", "   > ##### ERROR: FlatSurf: Deallocation failure: ", XWye
end do

VecTot = 0.0 ; VecEn = 0.0			! calc average gradient from all sets of Exe
do XWye = 1, WyeN
  if (VecGrad(XWye).NE.MissVal) then
    VecTot = VecTot + VecGrad(XWye)
    VecEn  = VecEn  + 1.0
  end if
end do
if (VecEn.GE.1) WyeGrad = VecTot / VecEn

VecTot = 0.0 ; VecEn = 0.0			! calc average intercept from all sets of Exe
do XWye = 1, WyeN
  if (VecCept(XWye).NE.MissVal) then
    VecTot = VecTot + VecCept(XWye)
    VecEn  = VecEn  + 1.0
  end if
end do
if (VecEn.GE.1) WyeCept = VecTot / VecEn

deallocate (DataVec,DataRef,VecCept,VecGrad, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FlatSurf: Deallocation failure: Wye #####"

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

allocate (Surf(ExeN,WyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FlatSurf: Allocation failure: Surf #####"

if (ExeGrad.NE.MissVal.AND.ExeCept.NE.MissVal.AND.WyeGrad.NE.MissVal.AND.WyeCept.NE.MissVal) then
  do XExe = 1, ExeN
    do XWye = 1, WyeN
      Surf (XExe,XWye) = (ExeCept+(ExeGrad*real(XExe))+WyeCept+(WyeGrad*real(XWye))) / 2.0
    end do
  end do
else
  Surf = MissVal
end if

end subroutine FlatSurf

!*******************************************************************************
! calculates the average of surrounding grid boxes
! all are already defined and filled
! nothing in HiVec is made MissVal, all is either filled or else returned as it is

subroutine CrudeSmoo (Extent,HiGrid,LoToHi,LoVec,HiVec)

real, pointer, dimension (:)		:: LoVec,HiVec

integer, pointer, dimension (:,:) 	:: HiGrid,LoToHi

integer, intent (in) :: Extent			! distance (boxes) away from box-to-fill to incorporate

real, parameter :: MissVal = -999.0

real :: OpTot, OpEn

integer :: AllocStat
integer :: ThisHiExe,ThisHiWye,SquareHiExe,SquareHiWye
integer :: HiExeN,HiWyeN,HiCellN

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

HiExeN = size (HiGrid,1) ; HiWyeN = size (HiGrid,2) ; HiCellN = HiExeN * HiWyeN

do ThisHiExe = 1, HiExeN
  do ThisHiWye = 1, HiWyeN
    OpTot = 0.0 ; OpEn  = 0.0
    
    do SquareHiExe = max(1,(ThisHiExe-Extent)), min(HiExeN,(ThisHiExe+Extent))
      do SquareHiWye = max(1,(ThisHiWye-Extent)), min(HiWyeN,(ThisHiWye+Extent))
        if (LoToHi(SquareHiExe,SquareHiWye).NE.MissVal) then
          if (LoVec(LoToHi(SquareHiExe,SquareHiWye)).NE.MissVal) then
            OpTot = OpTot + LoVec(LoToHi(SquareHiExe,SquareHiWye))
            OpEn  = OpEn  + 1.0
          end if
        end if
      end do
    end do
    
    if (OpEn.GT.0) HiVec(HiGrid(ThisHiExe,ThisHiWye)) = OpTot / OpEn
  end do
end do

end subroutine CrudeSmoo

!*******************************************************************************
! calculates the average of surrounding grid boxes
! all are already defined and filled
! if a cell in InVec is...
! ...filled, OutVec = InVec
! ...unfilled, we try to fill it based on surrounding cells, otherwise set to missing

subroutine NearNeigh (Extent,Grid,InVec,OutVec)

real, pointer, dimension (:)		:: InVec,OutVec

integer, pointer, dimension (:,:) 	:: Grid

integer, intent (in) :: Extent			! distance (boxes) away from box-to-fill to incorporate

real, parameter :: MissVal = -999.0

real :: OpTot, OpEn

integer :: AllocStat
integer :: ThisExe,ThisWye,SquareExe,SquareWye
integer :: ExeN,WyeN
integer :: ExeSqMin,ExeSqMax,WyeSqMin,WyeSqMax

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

ExeN = size (Grid,1) ; WyeN = size (Grid,2)

do ThisExe = 1, ExeN
  do ThisWye = 1, WyeN
   if (Grid(ThisExe,ThisWye).NE.MissVal) then
    if (InVec(Grid(ThisExe,ThisWye)).EQ.MissVal) then
     OpTot = 0.0 ; OpEn  = 0.0
     ExeSqMin = max(1,(ThisExe-Extent)) ; ExeSqMax = min(ExeN,(ThisExe+Extent))
     WyeSqMin = max(1,(ThisWye-Extent)) ; WyeSqMax = min(WyeN,(ThisWye+Extent))
     
     do SquareExe = ExeSqMin, ExeSqMax
      do SquareWye = WyeSqMin, WyeSqMax
        if (Grid(SquareExe,SquareWye).NE.MissVal) then
          if (InVec(Grid(SquareExe,SquareWye)).NE.MissVal) then
            OpTot = OpTot + InVec(Grid(SquareExe,SquareWye))
            OpEn  = OpEn  + 1.0
          end if
        end if
      end do
     end do
     
     if (OpEn.GE.1) then
     	OutVec(Grid(ThisExe,ThisWye)) = OpTot / OpEn
     else
        OutVec(Grid(ThisExe,ThisWye)) = MissVal
     end if 
    else
        OutVec(Grid(ThisExe,ThisWye)) = InVec(Grid(ThisExe,ThisWye))
    end if
   end if
   
  end do
end do

end subroutine NearNeigh

!*******************************************************************************
! extends a grim data set (1901-1995) to cover 20th century (1901-2000) by duplication

subroutine Dupli90s (DataIn,DataOut,YearADIn,YearADOut)

real, pointer, dimension (:,:,:) 	:: DataIn,DataOut

integer, pointer, dimension (:) 	:: YearADIn,YearADOut

integer :: AllocStat
integer :: YearInN, YearOutN, MonthN, BoxN
integer :: XYearIn,XYearOut,XMonth,XBox

YearInN = size (DataIn,1) ; MonthN = size (DataIn,2) ; BoxN = size (DataIn,3)

if (YearInN.EQ.95.AND.size(YearADIn).EQ.95.AND.YearADIn(1).EQ.1901) then
  YearOutN = 100
  
  allocate (DataOut  (YearOutN,MonthN,BoxN), &
  	    YearADOut(YearOutN),             stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Dupli90s: Allocation failure #####"

  do XYearOut = 1, YearOutN
    YearADOut(XYearOut) = 1900 + XYearOut

    if (XYearOut.LE.95) XYearIn = XYearOut
    if (XYearOut.GT.95) XYearIn = XYearOut - 5
    
    do XMonth = 1, MonthN
      do XBox = 1, BoxN
         DataOut(XYearOut,XMonth,XBox) = DataIn(XYearIn,XMonth,XBox)
      end do
    end do
  end do
  
  print*, "  > Duplicated 1991-1995 into 1996-2000." 
else
  print*, "  > Dupli90s: wrong period length. No duplication."
end if

end subroutine Dupli90s

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

end module Interpol