! gloops.f90
! module in which routines to handle .glo files are held
! contains: GloToGlo

module GloOps

use FileNames
use GloFiles

implicit none

contains

!*******************************************************************************
! transforms a .glo into an equivalent .glo with a different region set 

subroutine GloToGlo (CallOrigFile,CallMintFile,OrigMapIDLReg,MintMapIDLReg,CallWeights,MagMin,MagMax)

real, pointer, dimension (:,:),optional	:: CallWeights				! (LongN,LatN)
real, pointer, dimension (:)		:: OrigVec, MintVec			! data
real, allocatable, dimension (:)	:: MintTot, MintEn			! manipulations
real, allocatable, dimension (:,:)	:: Weights				! (LongN,LatN)
					
integer, pointer, dimension (:,:) 	:: OrigMapIDLReg,MintMapIDLReg		! maps IDL grid to regions
integer, allocatable, dimension (:)	:: MintRegSizes				! manipulations

character (len=80), intent(in) 		:: CallOrigFile,CallMintFile		! blank OK

real, intent(in), optional		:: MagMin,MagMax			! limits of magnitudes to save

real, parameter :: MissVal = -999.0

integer :: AllocStat, ReadStatus
integer :: LongN, LatN, OrigRegN, MintRegN
integer :: XLong, XLat, XReg

character (len=80) :: OrigFile,MintFile,OrigTitle,MintTitle, ModelFilePath
character (len=10) :: Format

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

OrigFile = LoadPath (CallOrigFile,".glo")					! get file paths
MintFile = SavePath (CallMintFile,".glo")

call LoadGloVec (OrigFile, OrigMapIDLReg, OrigVec, Silent=1)			! load .glo vector
call GetGloHeaders (OrigFile,OrigTitle,Format,LongN,LatN,ModelFilePath)		! load .glo headers

OrigRegN = maxval (OrigMapIDLReg) ; MintRegN = maxval (MintMapIDLReg)		! get region totals

MintTitle = trim(adjustl(OrigTitle)) // " -> re-reg"				! get Mint title

allocate (MintTot      (MintRegN), &
	  MintEn       (MintRegN), &
	  MintVec      (MintRegN), &
	  MintRegSizes (MintRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GloToGlo: Allocation failure: Mint* #####"
MintTot = 0.0 ; MintEn  = 0.0 ; MintVec= MissVal ; MintRegSizes = 0

allocate (Weights (LongN,LatN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GloToGlo: Allocation failure: Weights #####"

if (present(CallWeights)) then
  do XLong = 1, LongN
    do XLat = 1, LatN
      Weights(XLong,XLat) = CallWeights(XLong,XLat)
    end do
  end do
else
  Weights = 1.0
end if

do XLong = 1, LongN
  do XLat = 1, LatN
    if (MintMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is wanted for Mint region set    
     if (MintRegSizes (MintMapIDLReg (XLong,XLat)).NE.MissVal) then	! the region is non-missing thus far

      MintRegSizes (MintMapIDLReg (XLong,XLat)) = MintRegSizes (MintMapIDLReg (XLong,XLat)) + 1
      
      if (Weights(XLong,XLat).NE.MissVal) then			! weight is available    	
       if (OrigMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is specified in old region set    	
    	if (OrigVec (OrigMapIDLReg(XLong,XLat)).NE.MissVal) then	! value is valid    	
    	  MintTot (MintMapIDLReg(XLong,XLat)) = MintTot (MintMapIDLReg(XLong,XLat)) + (Weights(XLong,XLat)* &
    	  						OrigVec (OrigMapIDLReg(XLong,XLat)))
          MintEn  (MintMapIDLReg(XLong,XLat)) = MintEn  (MintMapIDLReg(XLong,XLat)) +  Weights(XLong,XLat)
        else
         MintRegSizes (MintMapIDLReg (XLong,XLat)) = MissVal
        end if
       else
        MintRegSizes (MintMapIDLReg (XLong,XLat)) = MissVal
       end if
      else
       MintRegSizes (MintMapIDLReg (XLong,XLat)) = MissVal
      end if 
     
     end if   
    end if
  end do
end do

do XReg = 1, MintRegN
  if (MintRegSizes(XReg).NE.MissVal.AND.MintEn(XReg).GT.0) MintVec (XReg) = MintTot(XReg) / MintEn(XReg)
end do

if (present(MagMin)) then
 do XReg = 1, MintRegN
  if (MintVec(XReg).NE.MissVal) then
   if (abs(MintVec(XReg)).LT.MagMin) MintVec(XReg) = sign(MagMin,MintVec(XReg))
  end if
 end do
end if

if (present(MagMax)) then
 do XReg = 1, MintRegN
  if (MintVec(XReg).NE.MissVal) then
   if (abs(MintVec(XReg)).GT.MagMax) MintVec(XReg) = sign(MagMax,MintVec(XReg))
  end if
 end do
end if

call SaveGlo (LongN,LatN,MintRegN,ModelFilePath,MintFile,MintTitle,MintVec,MintMapIDLReg)
if (MintRegN.LE.10) then
  print "(a)", "   > Means of the regions: "
  do XReg = 1, MintRegN
	print "(a,i2,e12.4)", "   > ", XReg, MintVec(XReg)
  end do
end if

deallocate (MintTot,MintEn,MintVec,MintRegSizes,Weights, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GloToGlo: Deallocation failure #####"

end subroutine GloToGlo

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

end module GloOps