! grid.f90
! module to hold standard routines for manipulating grids
! contains: GrabGrid,MakeHiResGrid,DefineShrinking
! FileNames,GrimFiles, must be further up the call line

module Grid

implicit none

contains

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

subroutine GrabGrid (CallMan,Grid,Bounds,BoxN,Quiet)

use FileNames
use GrimFiles

real, pointer, dimension (:,:,:) 	:: Data

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

real, dimension(4) 			:: Bounds

integer, intent (in), optional 	:: Quiet		! minimises verbosity
integer, intent (in) 		:: CallMan		! 1,2,MissVal
integer, intent (out) 		:: BoxN			! will be MissVal if QMan=2

real, parameter :: MissVal = -999.0

integer :: ReadStatus, AllocStat
integer :: QMan
integer :: ExeN,WyeN,XBound

character (len=80) :: Path,Info
character (len= 4) :: Suffix

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

if (CallMan.EQ.MissVal) then
  print*, "  > Define from grim (=1), or manually (=2) ?"
  do
	read (*,*,iostat=ReadStatus), QMan
	if (ReadStatus.LE.0.AND.QMan.GE.1.AND.QMan.LE.2) exit
  end do
else
  QMan = CallMan
end if

if (QMan.EQ.1) then
  print*, "  > Enter the grim filepath: "
  do
	read (*,*,iostat=ReadStatus), Path
	if (ReadStatus.LE.0.AND.Path.NE."") exit
  end do

  if (present(Quiet)) then
  	call LoadGrim (Data,Grid,YearAD,Bounds,Info,Path,"    ",Suffix,Quiet=1)
  else
  	call LoadGrim (Data,Grid,YearAD,Bounds,Info,Path,"    ",Suffix)
  end if
  
  ExeN = size(Grid,1) ; WyeN = size (Grid,2) ; BoxN = size (Data,3)

  print "(a,3i8)", "   > Grid dimensions and domain size: ", ExeN, WyeN, BoxN

  deallocate (Data,YearAD,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGrid: Deallocation failure #####"
else
  print*, "  > Enter the long and lat dimensions (boxes): "
  do
	read (*,*,iostat=ReadStatus), ExeN, WyeN
	if (ReadStatus.LE.0.AND.ExeN.GE.1.AND.WyeN.GE.1) exit
  end do
  
  allocate (Grid(ExeN,WyeN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGrid: Allocation failure #####"
  Grid   = MissVal
  Bounds = MissVal
  BoxN   = MissVal

  print*, "  > Enter the min,max long and min,max lat: "
  do
	read (*,*,iostat=ReadStatus), (Bounds(XBound),XBound=1,4)
	
	if (Bounds(1).LT.-180.OR.Bounds(2).GT.180.OR.Bounds(1).GE.Bounds(2)) then
	  print*, "  > There is a problem with the longitudes. Try again."
	  ReadStatus = 1
	end if
	
	if (Bounds(3).LT. -90.OR.Bounds(4).GT. 90.OR.Bounds(3).GE.Bounds(4)) then
	  print*, "  > There is a problem with the latitudes.  Try again."
	  ReadStatus = 1
	end if
	
	if (ReadStatus.LE.0) exit
  end do
  
end if

end subroutine GrabGrid

!*******************************************************************************
! the Call integers refer to the magnification and may be set to MissVal
! GridLo must be allocated in the call
! GridHi and LoToHi are returned from the subroutine
! GridHi (ExeHiN,WyeHiN) has no missing values, containing indices 1...(ExeHiN*WyeHiN)
! LoToHi (ExeHiN,WyeHiN) contains the most-contributing box from GridLo
!   and will be set to MissVal if there is no single most-contributing box

subroutine MakeHiResGrid (CallExeMulti,CallWyeMulti,GridLo,GridHi,LoToHi)

integer, pointer, dimension (:,:) 	:: GridLo,GridHi,LoToHi

real, intent (in) :: CallExeMulti,CallWyeMulti

real, parameter :: MissVal = -999.0

real :: ExeMulti,WyeMulti

integer :: AllocStat,ReadStatus
integer :: ExeLoN,WyeLoN,ExeHiN,WyeHiN
integer :: XExeLo,XWyeLo,XExeHi,XWyeHi,XExeAdd,XWyeAdd

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

ExeLoN = size (GridLo,1) ; WyeLoN = size (GridLo,2)

if (CallExeMulti.EQ.MissVal.OR.CallWyeMulti.EQ.MissVal) then
  print*, "  > Enter the magnification ratio for the long,lat axes: "
  do
	read (*,*,iostat=ReadStatus), ExeMulti, WyeMulti
	if (ReadStatus.LE.0.AND.ExeMulti.GE.1.AND.WyeMulti.GE.1) exit
  end do
else
  ExeMulti = CallExeMulti
  WyeMulti = CallWyeMulti
end if

ExeHiN = nint (ExeMulti * ExeLoN) ; WyeHiN = nint (WyeMulti * WyeLoN)

if (real(ExeHiN).EQ.(ExeMulti*real(ExeLoN)).AND.real(WyeHiN).EQ.(WyeMulti*real(WyeLoN))) then
  allocate (GridHi(ExeHiN,WyeHiN), &
  	    LoToHi(ExeHiN,WyeHiN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeHiResGrid: Allocation failure #####"
  
  do XExeHi = 1, ExeHiN
    do XWyeHi = 1, WyeHiN
      GridHi(XExeHi,XWyeHi) = ((XExeHi-1) * WyeHiN) + XWyeHi 
    end do
  end do

  do XExeLo = 1, ExeLoN
    do XExeAdd = 1, floor(ExeMulti)
      XExeHi = nint(real(XExeAdd) + real(XExeLo-1)*real(ExeMulti))
      
      do XWyeLo = 1, WyeLoN
        do XWyeAdd = 1, floor(WyeMulti)
          XWyeHi = nint(real(XWyeAdd) + real(XWyeLo-1)*real(WyeMulti))
          
          LoToHi(XExeHi,XWyeHi) = GridLo(XExeLo,XWyeLo)
        end do
      end do
    end do
  end do
else
  print*, "  > MakeHiResGrid: no grid created: incompatible with bounds."
end if

end subroutine MakeHiResGrid

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

subroutine DefineShrinking (OrigBounds,ShrinkBounds,OrigGrid,ShrinkGrid,ReGrid,ValidShrink)

integer, pointer, dimension (:,:)	:: OrigGrid,ShrinkGrid,ReGrid

real, dimension (4) 			:: OrigBounds,ShrinkBounds

integer, intent (out) 	:: ValidShrink

real, parameter :: MissVal = -999.0

real :: OrigLonPerBox,OrigLatPerBox
real :: ShrinkLonPerBox,ShrinkLatPerBox

integer :: AllocStat,ReadStatus 
integer :: TestBoxN,OrigExeN,OrigWyeN,OrigGlobalN
integer :: ShrinkExeN,ShrinkWyeN,ShrinkGlobalN,ShrinkExeCept,ShrinkWyeCept,ShrinkBoxN
integer :: XShrinkExe,XShrinkWye,XShrinkBox
integer :: QChangeLat,QChangeLon

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

OrigExeN = size (OrigGrid,1) ; OrigWyeN = size (OrigGrid,2)
ShrinkExeN = size (ShrinkGrid,1) ; ShrinkWyeN = size (ShrinkGrid,2)
ShrinkBoxN = maxval (ShrinkGrid)

OrigLonPerBox = (OrigBounds(2) - OrigBounds(1)) / real(OrigExeN)
OrigLatPerBox = (OrigBounds(4) - OrigBounds(3)) / real(OrigWyeN)

ShrinkLonPerBox = (ShrinkBounds(2) - ShrinkBounds(1)) / real(ShrinkExeN)
ShrinkLatPerBox = (ShrinkBounds(4) - ShrinkBounds(3)) / real(ShrinkWyeN)

OrigGlobalN   = nint (360.0 / OrigLonPerBox)   * nint (180.0 / OrigLatPerBox)  
ShrinkGlobalN = nint (360.0 / ShrinkLonPerBox) * nint (180.0 / ShrinkLatPerBox)

ValidShrink = 0

if (OrigGlobalN.NE.ShrinkGlobalN) then
  QChangeLat = 2 ; QChangeLon = 2
  print*, "  > The original and shrunken grids do not precisely match."
  
  if (OrigLonPerBox.NE.ShrinkLonPerBox) then
    print "(a,f12.4)", "   > The shrunken longitude box dimension is: ", ShrinkLonPerBox
    print "(a,f12.4)", "   > Decide (1=no,2=yes) whether to make it:  ", OrigLonPerBox
    do
	read (*,*,iostat=ReadStatus), QChangeLon
	if (ReadStatus.LE.0.AND.QChangeLon.GE.1.AND.QChangeLon.LE.2) exit
    end do
    if (QChangeLon.EQ.2) ShrinkLonPerBox = OrigLonPerBox
  end if
  
  if (OrigLatPerBox.NE.ShrinkLatPerBox) then
    print "(a,f12.4)", "   > The shrunken latitude box dimension is:  ", ShrinkLatPerBox
    print "(a,f12.4)", "   > Decide (1=no,2=yes) whether to make it:  ", OrigLatPerBox
    do
	read (*,*,iostat=ReadStatus), QChangeLat
	if (ReadStatus.LE.0.AND.QChangeLat.GE.1.AND.QChangeLat.LE.2) exit
    end do
    if (QChangeLat.EQ.2) ShrinkLatPerBox = OrigLatPerBox
  end if
  
  if (QChangeLat.EQ.1.OR.QChangeLon.EQ.1) then 
  	ValidShrink = 1
  	print*, "  > The original and shrunken grids do not match. The End."
  end if
end if

if (ValidShrink.NE.1) then
  if ((ShrinkBounds(1)-OrigBounds(1)).NE.0) then
  	ShrinkExeCept = nint ((ShrinkBounds(1)-OrigBounds(1)) / OrigLonPerBox)
  else
  	ShrinkExeCept = 0
  end if

  if ((ShrinkBounds(3)-OrigBounds(3)).NE.0) then
  	ShrinkWyeCept = nint ((ShrinkBounds(3)-OrigBounds(3)) / OrigLatPerBox)
  else
  	ShrinkWyeCept = 0
  end if
  
  print "(a,i8)", "   > Globalised original and shrunk grid sizes are both: ", OrigGlobalN

  allocate (ReGrid (ShrinkExeN,ShrinkWyeN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: DefineReGrid: Allocation failure #####"

  TestBoxN = 0
  do XShrinkExe = 1, ShrinkExeN
   do XShrinkWye = 1, ShrinkWyeN
    ReGrid(XShrinkExe,XShrinkWye) = OrigGrid((XShrinkExe+ShrinkExeCept),(XShrinkWye+ShrinkWyeCept))
    if (ReGrid(XShrinkExe,XShrinkWye).NE.MissVal) TestBoxN = TestBoxN + 1
   end do
  end do

  if (ShrinkBoxN.EQ.MissVal) then					! shrink grid not yet defined
   ShrinkBoxN = TestBoxN
  
   XShrinkBox = 0							! ...so we define it
   do XShrinkExe = 1, ShrinkExeN
    do XShrinkWye = 1, ShrinkWyeN
      if (ReGrid(XShrinkExe,XShrinkWye).NE.MissVal) then
        XShrinkBox = XShrinkBox + 1
        ShrinkGrid(XShrinkExe,XShrinkWye) = XShrinkBox
      end if
    end do
   end do
  
   print "(a,i8)", "   > Shrink grids defined. Shrunken domain size: ", ShrinkBoxN
  else									! shrink grid already defined
   if (ShrinkBoxN.NE.TestBoxN) then					! ...and domain size does not match 
    print*, ShrinkBoxN,TestBoxN
    
    print*, "  > Mismatch in size of shrunken domain. File not shrunk."    
    ValidShrink = 1
   else									! ...and domain size does match
    XShrinkBox = 0							! ...so we check the domain
    do XShrinkExe = 1, ShrinkExeN
      do XShrinkWye = 1, ShrinkWyeN
        if (ReGrid(XShrinkExe,XShrinkWye).NE.MissVal) then
          XShrinkBox = XShrinkBox + 1
          if (ShrinkGrid(XShrinkExe,XShrinkWye).NE.XShrinkBox) ValidShrink = 1
        else
          if (ShrinkGrid(XShrinkExe,XShrinkWye).NE.MissVal)    ValidShrink = 1
        end if
      end do
    end do
    
    if (ValidShrink.EQ.1) then
    	print*, "  > Mismatch in cells of shrunken domain. File not shrunk."    	
    else
        print "(a,i8)", "   > Shrink grids checked v. original grid. Shrunken domain size: ", ShrinkBoxN
    end if
   end if
  end if
end if

end subroutine DefineShrinking

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

end module Grid