! program written by Tim Mitchell on 21.8.00
! last modified on 4.9.00
! compares various region sets in .glo form
! f90 -o comparereg initialmod.f90 loadmod.f90 savemod.f90 comparereg.f90

program CompareReg

use InitialMod
use LoadMod
use SaveMod

implicit none

real, pointer,  dimension (:,:,:)		:: GloData
real, pointer,  dimension (:)			:: ThisGlo, BoxValid, BoxMiss

integer, pointer, dimension (:,:,:)		:: GloMap
integer, pointer, dimension (:,:)		:: MapIDLRaw, MapIDLReg, ThisMap
integer, pointer,  dimension (:)		:: MapRawReg, RegSizes, ADYear

character (len=20), dimension (:), pointer 	:: RegNames,LineNames

real, allocatable, dimension (:)		:: ScoreExp, ScoreAct, RegTot, RegSqTot

integer, allocatable, dimension (:,:,:,:)	:: Ring
integer, allocatable, dimension (:)		:: GloRegN

real, parameter :: MissVal = -999.0

real :: MissAccept
real :: RegThresh
real :: OpTot, OpMiss, OpEn, OpSqTot, OpMagTot, OpMean, OpStDev, OpMagAv, OpMin, OpDiff, OpAdjMin

integer :: LatN, LongN, GridChosen, GridDataN, RegN
integer :: Month0, Month1, MonthN, YearN, DecN
integer :: AllocStat, ReadStatus, MenuChoice
integer :: XReg, XYear, XGlo, XLat, XLong, XAdjacent
integer :: YLat, YLong
integer :: FirstFree, SelectA, SelectB, GotThisReg
integer :: AllGloN, LoadGloN, ThisRegN
integer :: NearLong, NearLat

character (len=80) :: GridFilePath, RegTitle, Blank
character (len=10) :: GridTitle

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

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

print*
print*, "  > ***** CompareReg: compares various region sets in .glo form *****"
print*

call Initialise
call MainMenu
call Finalise

close (99)

contains

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

subroutine MainMenu

MenuChoice = -1

do
  if      (MenuChoice.EQ.0)  then
    call SetMissAccept (MissAccept)
  else if (MenuChoice.EQ.1)  then
    call IdentifyTwo
    call CalcRMSE
  else if (MenuChoice.EQ.2)  then
    call IdentifyTwo
    call CompareAB
  else if (MenuChoice.EQ.3)  then
    call IdentifyTwo
    call ClosestNeighbour
  else if (MenuChoice.EQ.4)  then
    call IdentifyTwo
    call HypABEqual
  else if (MenuChoice.EQ.5)  then
    call IdentifyTwo
    call CalcCoeffVar
  else if (MenuChoice.EQ.6)  then
    call IdentifyOne
  else if (MenuChoice.EQ.30) then
    call PrintBoxStats
  else if (MenuChoice.EQ.31) then
    call PrintRegStats
  else
    print*
    print*, "  >  0. Change MissAccept, currently: ", MissAccept
    print*, "  >  1. RMSE in estimate of A boxes from B, saving using A regions"
    print*, "  >  2. A-B, box by box, printing stats"
    print*, "  >  3. Closest neighbour analysis: boxes=A, regions=B"
    print*, "  >  4. Examine A=B hypothesis by box"
    print*, "  >  5. Coeff of var in est of A boxes from B, saving with A reg"
    print*, "  > 30. Print box stats"
    print*, "  > 31. Print reg stats"
    print*, "  > 99. Exit"
  end if

  print*
  print*, "  > Main menu. Select your choice."
  do
	read (*,*,iostat=ReadStatus), MenuChoice
	if (ReadStatus.LE.0.AND.MenuChoice.GE.1.AND.MenuChoice.LE.99) exit
  end do
  
  if (MenuChoice.EQ.99) exit
end do

end subroutine MainMenu

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

subroutine Initialise

call GridSelect  (GridChosen,GridTitle,LongN,LatN,GridDataN,GridFilePath)

print*, "  > Select the number of .glo to load."
do
	read (*,*,iostat=ReadStatus), LoadGloN
	if (ReadStatus.LE.0.AND.LoadGloN.GE.1) exit
end do  

MissAccept = 10.0
Blank      = ""
FirstFree  = 1
AllGloN    = LoadGloN + 30

allocate (GloData    (AllGloN,LongN,LatN), &
	  GloMap     (AllGloN,LongN,LatN), &
	  GloRegN    (AllGloN),	           stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure #####"
GloData    = MissVal
GloMap     = MissVal
GloRegN    = MissVal

do XGlo = 1, LoadGloN
  call GrabGlo
end do

end subroutine Initialise

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

subroutine Finalise

deallocate (GloData,GloMap,GloRegN,MapIDLRaw,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Finalise: Deallocation failure #####"

print*

end subroutine Finalise

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

subroutine GrabGlo

print*, "  > Load .glo file into column: ", FirstFree

call RegSelect (GridChosen, LongN, LatN, GridDataN, MapIDLReg, RegSizes, RegNames, RegTitle, RegN)
call LoadGlo   (LongN, LatN, RegN, MapIDLReg, ThisGlo) 

GloRegN (FirstFree) = RegN

do XLong = 1, LongN
  do XLat = 1, LatN
    GloMap (FirstFree,XLong,XLat) = MapIDLReg (XLong,XLat)
    if (MapIDLReg (XLong,XLat).NE.MissVal) &
	    GloData (FirstFree,XLong,XLat) = ThisGlo (MapIDLReg (XLong,XLat))
  end do
end do

FirstFree = FirstFree + 1

deallocate (MapRawReg,MapIDLReg,RegSizes,RegNames,ThisGlo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGlo: Deallocation failure #####"

end subroutine GrabGlo

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

subroutine IdentifyOne

print*, "  > Select a column:"
do
	read (*,*,iostat=ReadStatus), SelectA
	if (ReadStatus.LE.0.AND.SelectA.GE.1.AND.SelectA.LT.FirstFree) exit
end do  

end subroutine IdentifyOne

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

subroutine IdentifyTwo

print*, "  > Select column A:"
do
	read (*,*,iostat=ReadStatus), SelectA
	if (ReadStatus.LE.0.AND.SelectA.GE.1.AND.SelectA.LT.FirstFree) exit
end do  

print*, "  > Select column B:"
do
	read (*,*,iostat=ReadStatus), SelectB
	if (ReadStatus.LE.0.AND.SelectB.GE.1.AND.SelectB.LT.FirstFree) exit
end do  

end subroutine IdentifyTwo

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

subroutine CalcRMSE

allocate (ThisMap  (LongN,LatN),       &
	  ThisGlo  (GloRegN(SelectA)), &
	  BoxValid (GloRegN(SelectA)), &
	  BoxMiss  (GloRegN(SelectA)), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcRMSE: Allocation failure #####"
ThisGlo  = 0.0
BoxValid = 0.0
BoxMiss  = 0.0

do XLong = 1, LongN
  do XLat = 1, LatN
    ThisMap (XLong,XLat) = GloMap (SelectA,XLong,XLat)
    
    if (GloMap(SelectA,XLong,XLat).NE.MissVal) then
      if (GloData(SelectA,XLong,XLat).NE.MissVal.AND.GloData(SelectB,XLong,XLat).NE.MissVal) then
        ThisGlo  (GloMap(SelectA,XLong,XLat)) = ThisGlo  (GloMap(SelectA,XLong,XLat)) + &
        					((GloData(SelectB,XLong,XLat)-GloData(SelectA,XLong,XLat))**2)
        BoxValid (GloMap(SelectA,XLong,XLat)) = BoxValid (GloMap(SelectA,XLong,XLat)) + 1.0
      else
        BoxMiss  (GloMap(SelectA,XLong,XLat)) = BoxMiss  (GloMap(SelectA,XLong,XLat)) + 1.0
      end if
    end if
  end do
end do

do XReg = 1, GloRegN(SelectA)
  RegThresh = (100.0 - MissAccept) * (BoxValid(XReg) + BoxMiss(XReg)) / 100.0
  
  if (BoxValid(XReg).GE.RegThresh) then
    ThisGlo (XReg) = ThisGlo (XReg) / BoxValid(XReg)
    ThisGlo (XReg) = sqrt (ThisGlo (XReg))
  else
    ThisGlo (XReg) = MissVal
  end if
end do

call SaveGlo (LongN,LatN,GloRegN(SelectA),GridFilePath,Blank,Blank,ThisGlo,ThisMap)

deallocate (ThisGlo,ThisMap,BoxValid,BoxMiss, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcRMSE: Deallocation failure #####"

end subroutine CalcRMSE

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

subroutine CalcCoeffVar

allocate (ThisMap  (LongN,LatN),       &
	  ThisGlo  (GloRegN(SelectA)), &
	  BoxValid (GloRegN(SelectA)), &
	  BoxMiss  (GloRegN(SelectA)), &
	  RegTot   (GloRegN(SelectA)), &
	  RegSqTot (GloRegN(SelectA)), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcCoeffVar: Allocation failure #####"
ThisGlo  = 0.0
BoxValid = 0.0
BoxMiss  = 0.0
RegTot   = 0.0
RegSqTot = 0.0

do XLong = 1, LongN
  do XLat = 1, LatN
    ThisMap (XLong,XLat) = GloMap (SelectA,XLong,XLat)
    
    if (GloMap(SelectA,XLong,XLat).NE.MissVal) then
      if (GloData(SelectB,XLong,XLat).NE.MissVal) then
        RegSqTot (GloMap(SelectA,XLong,XLat)) = RegSqTot (GloMap(SelectA,XLong,XLat)) + &
        						(GloData(SelectB,XLong,XLat)**2)
        RegTot   (GloMap(SelectA,XLong,XLat)) = RegTot   (GloMap(SelectA,XLong,XLat)) + &
        						 GloData(SelectB,XLong,XLat)
        BoxValid (GloMap(SelectA,XLong,XLat)) = BoxValid (GloMap(SelectA,XLong,XLat)) + 1.0
      else
        BoxMiss  (GloMap(SelectA,XLong,XLat)) = BoxMiss  (GloMap(SelectA,XLong,XLat)) + 1.0
      end if
    end if
  end do
end do

do XReg = 1, GloRegN(SelectA)
  RegThresh = (100.0 - MissAccept) * (BoxValid(XReg) + BoxMiss(XReg)) / 100.0
  
  if (BoxValid(XReg).GE.RegThresh.AND.BoxValid(XReg).GT.0) then
    OpMean = RegTot(XReg)/BoxValid(XReg)
    
    ThisGlo (XReg) = (RegSqTot(XReg)/BoxValid(XReg)) - (OpMean*OpMean)
    ThisGlo (XReg) = sqrt (ThisGlo (XReg))
    
    if (OpMean.GT.0) then
      ThisGlo (XReg) = ThisGlo (XReg) / OpMean
    else
      ThisGlo (XReg) = MissVal
    end if
  else
    ThisGlo (XReg) = MissVal
  end if

  write (99,"(i4,f12.4)"), XReg, ThisGlo (XReg)  		! #########################
end do

call SaveGlo (LongN,LatN,GloRegN(SelectA),GridFilePath,Blank,Blank,ThisGlo,ThisMap)

deallocate (ThisGlo,ThisMap,BoxValid,BoxMiss,RegTot,RegSqTot, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcCoeffVar: Deallocation failure #####"

end subroutine CalcCoeffVar

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

subroutine CompareAB

OpMiss  = 0.0
OpTot   = 0.0
OpEn    = 0.0
OpSqTot = 0.0

do XLong = 1, LongN
  do XLat = 1, LatN
    if (GloData(SelectA,XLong,XLat).NE.MissVal.AND.GloData(SelectB,XLong,XLat).NE.MissVal) then
      OpTot   = OpTot   +  (GloData(SelectA,XLong,XLat)-GloData(SelectB,XLong,XLat))
      OpSqTot = OpSqTot + ((GloData(SelectA,XLong,XLat)-GloData(SelectB,XLong,XLat)) ** 2)
      OpEn    = OpEn    +   1.0
    else
      OpMiss  = OpMiss  +   1.0
    end if
  end do
end do

print "(4e12.5)", OpMiss, OpEn, OpTot, OpSqTot

end subroutine CompareAB

! ******************************************************************************
! measures straight line distance from box value to A=B and prints global results

subroutine HypABEqual

OpMiss  = 0.0
OpTot   = 0.0
OpEn    = 0.0
OpMean  = MissVal

do XLong = 1, LongN
  do XLat = 1, LatN
    if (GloData(SelectA,XLong,XLat).NE.MissVal.AND.GloData(SelectB,XLong,XLat).NE.MissVal) then
      OpTot   = OpTot  + (sqrt(0.5) * abs(GloData(SelectA,XLong,XLat)-GloData(SelectB,XLong,XLat)))
      OpEn    = OpEn   + 1.0
    else
      OpMiss  = OpMiss + 1.0
    end if
  end do
end do

if (OpEn.GE.1) OpMean = OpTot / OpEn

print "(4e12.5)", OpMiss, OpEn, OpTot, OpMean

end subroutine HypABEqual

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

subroutine MakeRing

allocate (Ring (LongN,LatN,8,2), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeRing: Allocation failure #####"
Ring = MissVal

do XLong = 1, LongN			! identify the 8 grid boxes adjacent to each 
  do XLat = 1, LatN
    XAdjacent = 0

    do YLong = (XLong-1), (XLong+1)
      if      (YLong.LT.1)     then 
      		NearLong = LongN
      else if (YLong.GT.LongN) then
      		NearLong = 1
      else
      		NearLong = YLong
      end if

      do YLat = (XLat-1), (XLat+1)
        if (YLat.LT.1.OR.YLat.GT.LatN) then 
      		NearLat = MissVal
        else
      		NearLat = YLat
        end if
    
	if (NearLong.NE.XLong.OR.NearLat.NE.XLat) then
          if (NearLong.NE.MissVal.AND.NearLat.NE.MissVal) then
	          XAdjacent = XAdjacent + 1
	          Ring (XLong,XLat,XAdjacent,1) = NearLong
	          Ring (XLong,XLat,XAdjacent,2) = NearLat
	  end if
	end if
      end do
    end do
  end do
end do          

end subroutine MakeRing

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

subroutine ClosestNeighbour

call MakeRing

allocate (ThisMap  (LongN,LatN),       &
	  ThisGlo  (GloRegN(SelectB)), &
	  ScoreExp (GloRegN(SelectB)), &
	  ScoreAct (GloRegN(SelectB)), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClosestNeighbour: Allocation failure #####"
ScoreExp = 0.0
ScoreAct = 0.0
					! calc actual score for each region
do XLong = 1, LongN			! calc a priori E(score) for each region
 do XLat = 1, LatN
  if (GloMap(SelectB,XLong,XLat).NE.MissVal) then		! grid box of interest is part of a region
   if (GloData(SelectA,XLong,XLat).NE.MissVal) then		! grid box of interest has data

    OpTot  = 0.0						! no. of ring in same reg as box, + with data
    OpEn   = 0.0						! no. in ring with data
    
    OpMin    = 100000.0						! min diff between box and adjacent
    OpAdjMin = MissVal						! adjacent box with min diff
    
    do XAdjacent = 1, 8
     if (Ring(XLong,XLat,XAdjacent,1).NE.MissVal) then 		! adjacent box identified
      			 					! adjacent box has data
      if (GloData(SelectA,Ring(XLong,XLat,XAdjacent,1),Ring(XLong,XLat,XAdjacent,2)).NE.MissVal) then

        OpEn = OpEn + 1.0					! increment no. in ring with data
									 
        if (GloMap(SelectB,XLong,XLat).EQ.&			! adjacent box is in same region
        	GloMap(SelectB,Ring(XLong,XLat,XAdjacent,1),Ring(XLong,XLat,XAdjacent,2))) &
        	OpTot = OpTot + 1.0				! increment no. of ring in same reg
        
        OpDiff = abs (GloData(SelectA,XLong,XLat) - &		! diff between box and adjacent
        		GloData(SelectA,Ring(XLong,XLat,XAdjacent,1),Ring(XLong,XLat,XAdjacent,2)) )
        
        if (OpDiff.LT.OpMin) then				! diff is min so far
          OpMin    = OpDiff					! make min diff the current diff
          OpAdjMin = XAdjacent					! make min adjacent the current adjacent
        end if        
        
      end if
      
     end if
    end do
			
    if (OpEn.GT.0) then						! box(es) in ring have data

      if ((OpTot/OpEn).NE.1) then				! ignore box if all adjacent are in region
    								! increment a priori E(score)	
        ScoreExp(GloMap(SelectB,XLong,XLat)) = ScoreExp(GloMap(SelectB,XLong,XLat)) + (OpTot/OpEn)

        if (GloMap(SelectB,XLong,XLat).EQ.&			! min adjacent is of same region
      		GloMap(SelectB,Ring(XLong,XLat,OpAdjMin,1),Ring(XLong,XLat,OpAdjMin,2))) &    	
    		ScoreAct(GloMap(SelectB,XLong,XLat)) = ScoreAct(GloMap(SelectB,XLong,XLat)) + 1.0
								! increment actual score
      end if
      
    end if
    			
   end if
  end if
 end do
end do

do XLong = 1, LongN						! load ThisMap
  do XLat = 1, LatN
    ThisMap (XLong,XLat) = GloMap (SelectB,XLong,XLat)
  end do
end do

do XReg = 1, GloRegN(SelectB)					! load ThisReg
  if (ScoreExp(XReg).GT.0) then
    ThisGlo (XReg) = ScoreAct (XReg) / ScoreExp (XReg)
  end if
end do

call SaveGlo (LongN,LatN,GloRegN(SelectB),GridFilePath,Blank,Blank,ThisGlo,ThisMap)

deallocate (ScoreExp,ScoreAct,Ring,ThisMap,ThisGlo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClosestNeighbour: Deallocation failure #####"

end subroutine ClosestNeighbour

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

subroutine PrintBoxStats

print "(a4,5a12)", " Col", "         Sum", "  Sum Squares", "        Mean", "       StDev", "       AvMag" 

do XGlo = 1, (FirstFree-1)
  OpMiss   = 0.0
  OpTot    = 0.0
  OpMagTot = 0.0
  OpEn     = 0.0
  OpSqTot  = 0.0
  OpMean   = MissVal
  OpStDev  = MissVal
  
  do XLong = 1, LongN
   do XLat = 1, LatN
    if (GloData(XGlo,XLong,XLat).NE.MissVal) then
      OpTot    = OpTot    +  GloData(XGlo,XLong,XLat)
      OpSqTot  = OpSqTot  + (GloData(XGlo,XLong,XLat) ** 2)
      OpMagTot = OpMagTot +  abs (GloData(XGlo,XLong,XLat))
      OpEn     = OpEn     +  1.0
    else
      OpMiss  = OpMiss    +  1.0
    end if
   end do
  end do
  
  if (OpEn.GE.1) then
  	OpMean  =  OpTot/OpEn
  	OpStDev = (OpSqTot/OpEn) - ((OpTot/OpEn) ** 2)
  	OpStDev = sqrt (OpStDev)
  	OpMagAv =  OpMagTot/OpEn
  end if
  
  print "(i4,5e12.5)", XGlo, OpTot, OpSqTot, OpMean, OpStDev, OpMagAv
end do

end subroutine PrintBoxStats

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

subroutine PrintRegStats

print "(a4,5a12)", " Col", "         Sum", "  Sum Squares", "        Mean", "       StDev", "       AvMag" 

do XGlo = 1, (FirstFree-1)
  OpMiss   = 0.0
  OpTot    = 0.0
  OpMagTot = 0.0
  OpEn     = 0.0
  OpSqTot  = 0.0
  OpMean   = MissVal
  OpStDev  = MissVal
  
  do XReg = 1, GloRegN (XGlo)
   GotThisReg = 0
   
   do XLong = 1, LongN
    do XLat = 1, LatN
     if (GotThisReg.EQ.0) then
       if (GloMap(XGlo,XLong,XLat).EQ.XReg.AND.GloData(XGlo,XLong,XLat).NE.MissVal) then
      	OpTot    = OpTot    +  GloData(XGlo,XLong,XLat)
      	OpSqTot  = OpSqTot  + (GloData(XGlo,XLong,XLat) ** 2)
      	OpMagTot = OpMagTot +  abs (GloData(XGlo,XLong,XLat))
      	OpEn     = OpEn     +  1.0

      	GotThisReg = 1
       end if
     end if
    end do
   end do
  end do
  
  if (OpEn.GE.1) then
  	OpMean  =  OpTot/OpEn
  	OpStDev = (OpSqTot/OpEn) - ((OpTot/OpEn) ** 2)
  	OpStDev = sqrt (OpStDev)
  	OpMagAv =  OpMagTot/OpEn
  end if
  
  print "(i4,5e12.5)", XGlo, OpTot, OpSqTot, OpMean, OpStDev, OpMagAv
end do

end subroutine PrintRegStats

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

end program CompareReg