! maskpolesglo.f90
! f90 main program written on 03.03.00 by Tim Mitchell
! last modification on 03.03.00
! f90 -o maskpolesglo initialmod.f90 loadmod.f90 savemod.f90 maskpolesglo.f90

program MaskPolesGlo

use InitialMod
use LoadMod
use SaveMod

implicit none

real, pointer, dimension (:)		:: ThisGlo

integer, pointer, dimension (:) 	:: WorkMapRawReg, WorkRegSizes
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw, WorkMapIDLReg

character (len=20), pointer, dimension (:) 	:: WorkRegNames

real, parameter :: MissVal = -999.0

real :: RegSqTot, RegTot, RegEn, RegMean, RegStDev, RegRMSE

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN
integer :: WorkRegN, WorkGloN
integer :: AllocStat, ReadStatus
integer :: GloFileN, XGloFile, XLong, XReg

character (len=10) :: WorkGridTitle
character (len=40) :: WorkRegTitle
character (len=80) :: WorkGridFilePath, WorkGloTitle, Blank

!*******************************************************************************
! main loop

Blank = ""

print*, "  > ***** MaskPolesGlo *****"
print*

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

print*, "  > Select the number of .glo files with poles to be masked out: "
do
	read (*,*,iostat=ReadStatus), GloFileN
	if (ReadStatus.LE.0.AND.GloFileN.GE.1) exit
end do

do XGloFile = 1, GloFileN
  print*, "  > File number: ", XGloFile
  call LoadGlo (WorkLongN,WorkLatN,WorkRegN,WorkMapIDLReg,ThisGlo)
  
  do XLong = 1, WorkLongN
    XReg = WorkMapIDLReg (XLong,1)
    ThisGlo (XReg) = MissVal
    XReg = WorkMapIDLReg (XLong,WorkLatN)
    ThisGlo (XReg) = MissVal
  end do
  
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank, &
	      ThisGlo,WorkMapIDLReg)
  
  RegSqTot = 0.0
  RegTot   = 0.0
  RegEn    = 0.0
  RegMean  = MissVal
  RegStDev = MissVal
  RegRMSE  = MissVal
  
  do XReg = 1, WorkRegN
    if (ThisGlo(XReg).NE.MissVal) then
      RegSqTot = RegSqTot + ThisGlo(XReg) ** 2
      RegTot   = RegTot   + ThisGlo(XReg)
      RegEn    = RegEn    + 1.0
    end if
  end do
  
  if (RegEn.GE.0) then
    RegMean  = RegTot/RegEn
    RegStDev = (RegSqTot/RegEn) - (RegMean**2)
    RegStDev = sqrt (RegStDev)
    RegRMSE  = (RegSqTot/RegEn)
    RegRMSE  = sqrt (RegRMSE)
  end if
  
  print*, "  > Stats when all regions are combined into a global mean: "
  print*, "  > Regions in set, regions valid, mean, stdev, rmse: "
  print "(i6,4f12.5)", WorkRegN, RegEn, RegMean, RegStDev, RegRMSE
end do

deallocate (ThisGlo)
deallocate (WorkMapIDLReg,WorkRegSizes,WorkRegNames)

print*

end program MaskPolesGlo